text:
"""
Package: Forecast

    ccf(x1::{AbstractVector,DataFrame},
        x2::{AbstractVector,DataFrame};
        type = "cor",
        lag = Integer(ceil(10*log10(length(x1)))),
        alpha = (0.95,0.99))

Compute the cross-correlation or cros-covariance of two univariate series.

The results are normalized to preserve homoscedasticity. The distribution used to normalize the data is an approximation of a Fisher Transformation via a Normal Distribution. There is a plot recipe for the returned object, if the type is `cor` the plot will also show confidence intervals for the given alpha values.

If, for a given integer `k`, `x2` repeats `x1` values such that x1[t] = x2[t+k] for all `i` then high correlation value will be placed *at the right from the center* in the results. That is, this convention will be represented in the plots as `x1_t = x2_{t+k} -> _____0__k__` meaning x2 behavior can be predicted by x1 in k units.

# Arguments
- `x1`: Vector or uni-dimensional DataFrame of data.
- `x2`: Vector or uni-dimensional DataFrame of data.
- `type`: Valid values are "cor" for correlation (default) and "cov" for convariance.
- `lag`: Maximum number of lags.
- `alpha`: A tuple with two thresholds (t1,t2) with t1 <= t2 to plot confidence intervals. The default values are 0.95 and 0.99.

# Returns
A CCF object.

# Examples
```julia-repl
julia> x1 = rand(100);
x2 = circshift(x1,6);
res = ccf(x1, x2; type="cor");
plot(res)
```
"""
function ccf(df1::DataFrame,
             df2::DataFrame;
             type::String = "cor",
             lag::Integer = Integer(ceil(10*log10(length(x1)))),
             alpha::Tuple{Real,Real} = (0.95,0.99))

    x1 = Array(df1[:,eltype.(eachcol(df1)) .<: Real])
    x2 = Array(df1[:,eltype.(eachcol(df1)) .<: Real])

    ccf(x1, x2; type = type, lag = lag, alpha = alpha)

end

function ccf(x1::AbstractVector{<:Real},
             x2::AbstractVector{<:Real};
             type::String = "cor",
             lag::Integer = Integer(ceil(10*log10(length(x1)))),
             alpha::Tuple{Real,Real} = (0.95,0.99))

    N = length(x1)
    @assert N == length(x2) "Vectors should be of equal size"
    @assert N >= 4 "Vectors should have at least four values"
    @assert isnothing(findfirst(isnan,x1)) "No missing values allowed"
    @assert isnothing(findfirst(isnan,x2)) "No missing values allowed"
    @assert isnothing(findfirst(ismissing,x1)) "No missing values allowed"
    @assert isnothing(findfirst(ismissing,x2)) "No missing values allowed"
    @assert type in  ["cor","cov"] "The options for `type` are `cor` and `cov`"
    @assert 1 <= lag <= N-4
    @assert length(alpha) == 2
    @assert 0.0 < alpha[1] < alpha[2] < 1.0

    # auto-ccf
    auto = x1 == x2

    ft = (type == "cor" ? Statistics.cor : Statistics.cov)

    x = []

    for i in 1:lag
        push!(x,ft(x1[i+1:N],x2[1:(N-i)]))
    end

    if auto
        kx = collect(N-1:-1:N-lag)
    else
        x = vcat(ft(x1,x2),x) # adding central value (lag=0)
        x = reverse(x)
        for i in 1:lag
            push!(x,ft(x2[i+1:N],x1[1:(N-i)]))
        end
        kx = vcat(collect(N-lag:N-1), N, collect(N-1:-1:N-lag))
    end

    # ccf normalized with the standard error for the
    # Normal distsribution of an approximation for
    # the Fisher transformation
    function fse(v::Real)::AbstractFloat
        1/sqrt(v-3)
    end
    k = fse.(kx)/fse(N)
    ccf_res = x ./ k

    call = "ccf("* (auto ? "x" : "x1, x2") *
        "; type=\""*type*
        "\", lag="*string(lag)*
        ", alpha="*string(alpha)*")"

    # Confidence Intervals
    a1 = alpha[1]
    a2 = alpha[2]
    z1 = Distributions.quantile(Normal(), a1 + (1-a1)/2)
    ci1 = z1*fse(N)
    z2 = Distributions.quantile(Normal(), a2 + (1-a2)/2)
    ci2 = z2*fse(N)
    ci = (ci1,ci2)

    CCF(ccf_res, N, type, lag, alpha, ci, auto, call)

end


text:
using Documenter
using Plots

push!(LOAD_PATH,"../src/")
using UnderwaterAcoustics

makedocs(
  sitename = "UnderwaterAcoustics.jl",
  format = Documenter.HTML(prettyurls = false),
  linkcheck = !("skiplinks" in ARGS),
  pages = Any[
    "Home" => "index.md",
    "Manual" => Any[
      "uw_basic.md",
      "pm_basic.md",
      "pm_envref.md",
      "pm_api.md"
    ],
    "Propagation models" => Any[
      "pm_pekeris.md",
      "pm_rays.md",
      "pm_bellhop.md"
    ],
    "Tutorials" => Any[
      "tut_turing.md",
      "tut_autodiff.md"
    ]
  ]
)

deploydocs(
  repo = "github.com/org-arl/UnderwaterAcoustics.jl.git",
  branch = "gh-pages",
  devbranch = "master",
  devurl = "dev",
  versions = ["stable" => "v^", "v#.#", "dev" => "dev"]
)

ext: jl


text:
facts("unbiased functions") do
    context("unbiased_indmax samples over ties") do
        @fact unbiased_indmax([1,2,3,4]) --> 4
        @fact unbiased_indmax([3,2,1]) --> 1
        @fact Set([unbiased_indmax([3,4,4,3]) for _ in 1:100]) --> Set([2, 3])
    end

    context("unbiased_indmax samples over ties") do
        @fact unbiased_indmin([1,2,3,4]) --> 1
        @fact unbiased_indmin([3,2,1]) --> 3
        @fact Set([unbiased_indmin([3,4,4,3]) for _ in 1:100]) --> Set([1, 4])
    end
end


text:
using Plots

numberslist=[[1],[1,1]]

for step in 3:100
    R=[1]
    for n in 2:(step-1)
        number=numberslist[step-1][n-1]+numberslist[step-1][n]
        push!(R, number)
    end
    push!(R, 1)
    push!(numberslist,R)

end
#println(numberslist)
odd_num_pos=[]
even_num_pos=[]
y=0.0
x=0.0
for row in numberslist
    x=y/2
    for number in row
        (number%2==0) ? (push!(even_num_pos, [x;y])) : (push!(odd_num_pos, [x;y]))
        x+=1
    end
    y-=1
end

scatter(hcat(odd_num_pos...)[1,:], hcat(odd_num_pos...)[2,:], color=:darkgreen, xlabel="x", ylabel="y",markersize=0.7,legend=false)
scatter!(hcat(even_num_pos...)[1,:], hcat(even_num_pos...)[2,:], color=:lightsalmon, xlabel="x", ylabel="y", markersize=0.7,legend=false)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter2\\report\\KhayyamTriangle.png")


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BallNeighborhood{N}(radius, metric=Euclidean())

A `N`-dimensional ball neighborhood with `radius` and `metric`.
"""
struct BallNeighborhood{T,N,M} <: AbstractNeighborhood{T,N}
  radius::T
  metric::M
end

function BallNeighborhood{N}(radius::T, metric::M=Euclidean()) where {M,N,T}
  @assert radius > 0 "radius must be positive"
  @assert N > 0 "number of dimensions must be positive"
  BallNeighborhood{T,N,M}(radius, metric)
end

"""
    radius(ball)

Return the radius of the `ball`.
"""
radius(ball::BallNeighborhood) = ball.radius

"""
    metric(ball)

Return the metric of the `ball`.
"""
metric(ball::BallNeighborhood) = ball.metric

isneighbor(ball::BallNeighborhood, xₒ::AbstractVector, x::AbstractVector) =
  evaluate(ball.metric, xₒ, x) ≤ ball.radius

function volume(ball::BallNeighborhood{T,N,Euclidean}) where {N,T}
  # https://en.wikipedia.org/wiki/Volume_of_an_n-ball
  R = ball.radius
  (π^(N/2) / gamma(N/2 + 1)) * R^N
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, ball::BallNeighborhood{T,N,M}) where {M,N,T}
  r = ball.radius
  d = ball.metric
  print(io, "$(N)D BallNeighborhood($r, $d)")
end

function Base.show(io::IO, ::MIME"text/plain",
                   ball::BallNeighborhood{T,N,M}) where {M,N,T}
  println(io, "$(N)D BallNeighborhood")
  println(io, "  radius: ", ball.radius)
  print(  io, "  metric: ", ball.metric)
end


text:
module GnssWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractAnalyticalPropagator = @jimport org.orekit.propagation.analytical.AbstractAnalyticalPropagator
const AbstractGNSSPropagator = @jimport org.orekit.propagation.analytical.gnss.AbstractGNSSPropagator
const AbstractPropagator = @jimport org.orekit.propagation.AbstractPropagator
const AdditionalStateProvider = @jimport org.orekit.propagation.AdditionalStateProvider
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BeidouOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.BeidouOrbitalElements
const BeidouPropagator = @jimport org.orekit.propagation.analytical.gnss.BeidouPropagator
const BoundedPropagator = @jimport org.orekit.propagation.BoundedPropagator
const Class = @jimport java.lang.Class
const ClockCorrectionsProvider = @jimport org.orekit.propagation.analytical.gnss.ClockCorrectionsProvider
const Collection = @jimport java.util.Collection
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const Frame = @jimport org.orekit.frames.Frame
const GLONASSAnalyticalPropagator = @jimport org.orekit.propagation.analytical.gnss.GLONASSAnalyticalPropagator
const GLONASSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GLONASSOrbitalElements
const GNSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GNSSOrbitalElements
const GPSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GPSOrbitalElements
const GPSPropagator = @jimport org.orekit.propagation.analytical.gnss.GPSPropagator
const GalileoOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GalileoOrbitalElements
const GalileoPropagator = @jimport org.orekit.propagation.analytical.gnss.GalileoPropagator
const IRNSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.IRNSSOrbitalElements
const IRNSSPropagator = @jimport org.orekit.propagation.analytical.gnss.IRNSSPropagator
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const PVCoordinates = @jimport org.orekit.utils.PVCoordinates
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const QZSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.QZSSOrbitalElements
const QZSSPropagator = @jimport org.orekit.propagation.analytical.gnss.QZSSPropagator
const SBASOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.SBASOrbitalElements
const SBASPropagator = @jimport org.orekit.propagation.analytical.gnss.SBASPropagator
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeStamped = @jimport org.orekit.time.TimeStamped
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates

export AbstractGNSSPropagator
export BeidouOrbitalElements
export BeidouPropagator
export ClockCorrectionsProvider
export GLONASSAnalyticalPropagator
export GLONASSOrbitalElements
export GNSSOrbitalElements
export GPSOrbitalElements
export GPSPropagator
export GalileoOrbitalElements
export GalileoPropagator
export IRNSSOrbitalElements
export IRNSSPropagator
export QZSSOrbitalElements
export QZSSPropagator
export SBASOrbitalElements
export SBASPropagator
export add_additional_state_provider
export add_event_detector
export clear_events_detectors
export equals
export get_a_gf0
export get_a_gf1
export get_additional_state
export get_additional_state_providers
export get_af0
export get_af1
export get_af2
export get_aodc
export get_aode
export get_attitude_provider
export get_beidou_orbital_elements
export get_bgd
export get_bgde1e5a
export get_bgde5b_e1
export get_cic
export get_cis
export get_class
export get_crc
export get_crs
export get_cuc
export get_cus
export get_date
export get_delta_i
export get_delta_t
export get_delta_t_dot
export get_e
export get_ecef
export get_eci
export get_events_detectors
export get_frame
export get_galileo_orbital_elements
export get_gamma_n
export get_generated_ephemeris
export get_glonass_orbital_elements
export get_gps_orbital_elements
export get_i0
export get_i_dot
export get_initial_state
export get_iod
export get_iod_nav
export get_iodc
export get_iode
export get_iodec
export get_iodn
export get_irnss_orbital_elements
export get_lambda
export get_m0
export get_managed_additional_states
export get_mean_motion
export get_mode
export get_mu
export get_n4
export get_na
export get_name
export get_omega0
export get_omega_dot
export get_pa
export get_prn
export get_pv_coordinates
export get_pv_provider
export get_qzss_orbital_elements
export get_sbas_orbital_elements
export get_sma
export get_tgd
export get_tgd1
export get_tgd2
export get_time
export get_tn
export get_toc
export get_week
export get_x
export get_x_dot
export get_x_dot_dot
export get_y
export get_y_dot
export get_y_dot_dot
export get_z
export get_z_dot
export get_z_dot_dot
export hash_code
export is_additional_state_managed
export notify
export notify_all
export propagate
export propagate_in_ecef
export reset_initial_state
export set_attitude_provider
export set_ephemeris_mode
export set_master_mode
export set_slave_mode
export to_string
export wait

include("abstract_gnss_propagator.jl")
include("beidou_orbital_elements.jl")
include("beidou_propagator.jl")
include("clock_corrections_provider.jl")
include("galileo_orbital_elements.jl")
include("galileo_propagator.jl")
include("glonass_analytical_propagator.jl")
include("glonass_orbital_elements.jl")
include("gnss_orbital_elements.jl")
include("gps_orbital_elements.jl")
include("gps_propagator.jl")
include("irnss_orbital_elements.jl")
include("irnss_propagator.jl")
include("qzss_orbital_elements.jl")
include("qzss_propagator.jl")
include("sbas_orbital_elements.jl")
include("sbas_propagator.jl")

end


text:
using SparseArrays
using LinearAlgebra
using SparseArrays
using Printf
using DelimitedFiles

include("sim_funs.jl")
include("square_circle.jl")

function compute_energy(q, rhsops, Np)
    nelems = length(rhsops)
    energy = 0
    Nq = round(Int, sqrt(Np))
    @assert Np == Nq^2
    lenq0 = 2 * Np + 2 * Nq + 2 * Nq

    # Loop over blocks are compute the energy in the blocks at this time
    for e in 1:nelems
        qe = @view q[(lenq0 * (e - 1) + 1):(e * lenq0)]

        qe_u = qe[1:Np]
        qe_v = qe[Np .+ (1:Np)]
        qe_û1 = qe[(2Np + 1):(2Np + Nq)]
        qe_û2 = qe[(2Np + 1 + Nq):(2(Np + Nq))]
        qe_û3 = qe[(2(Np + Nq) + 1):(2(Np + Nq) + Nq)]
        qe_û4 = qe[(2(Np + Nq) + Nq + 1):(2(Np + Nq + Nq))]

        τ̂ = (
            rhsops[e].nCB[1] * qe_u + 1 * rhsops[e].nCnΓ[1] * qe_û1 -
            1 * rhsops[e].nCnΓL[1] * qe_u,
            rhsops[e].nCB[2] * qe_u + 1 * rhsops[e].nCnΓ[2] * qe_û2 -
            1 * rhsops[e].nCnΓL[2] * qe_u,
            rhsops[e].nCB[3] * qe_u + 1 * rhsops[e].nCnΓ[3] * qe_û3 -
            1 * rhsops[e].nCnΓL[3] * qe_u,
            rhsops[e].nCB[4] * qe_u + 1 * rhsops[e].nCnΓ[4] * qe_û4 -
            1 * rhsops[e].nCnΓL[4] * qe_u,
        )

        Mu = rhsops[e].Ã
        Mv = rhsops[e].JH

        Mû1 = rhsops[e].H[1]
        Mû2 = rhsops[e].H[2]
        Mû3 = rhsops[e].H[3]
        Mû4 = rhsops[e].H[4]

        Mτ̂1 = rhsops[e].X[1] * rhsops[e].H[1]
        Mτ̂2 = rhsops[e].X[2] * rhsops[e].H[2]
        Mτ̂3 = rhsops[e].X[3] * rhsops[e].H[3]
        Mτ̂4 = rhsops[e].X[4] * rhsops[e].H[4]

        Mu1 =
            (rhsops[e].nCB[1])' *
            rhsops[e].X[1] *
            rhsops[e].H[1] *
            rhsops[e].nCB[1]
        Mu2 =
            (rhsops[e].nCB[2])' *
            rhsops[e].X[2] *
            rhsops[e].H[2] *
            rhsops[e].nCB[2]
        Mu3 =
            (rhsops[e].nCB[3])' *
            rhsops[e].X[3] *
            rhsops[e].H[3] *
            rhsops[e].nCB[3]
        Mu4 =
            (rhsops[e].nCB[4])' *
            rhsops[e].X[4] *
            rhsops[e].H[4] *
            rhsops[e].nCB[4]

        energy +=
            0.5 * qe_v' * Mv * qe_v +
            0.5 * qe_u' * Mu * qe_u +
            0.5 * (
                1 * τ̂[1]' * Mτ̂1 * τ̂[1] - 1 * qe_u' * Mu1 * qe_u +
                1 * τ̂[2]' * Mτ̂2 * τ̂[2] - 1 * qe_u' * Mu2 * qe_u +
                1 * τ̂[3]' * Mτ̂3 * τ̂[3] - 1 * qe_u' * Mu3 * qe_u +
                1 * τ̂[4]' * Mτ̂4 * τ̂[4] - 1 * qe_u' * Mu4 * qe_u
            )
    end
    return sqrt(energy)
end

let
    friction(V) = asinh(V)

    sbp_order = 6
    cfl = 1 / 2

    data = Dict()

    # This is the base mesh size in each dimension on each element.
    Nvals = 17 * 2 .^ (3:-1:1)
    for N in Nvals

        # solver times (output done at each break point)
        ts = 0:0.01:3

        bc_map = [
            BC_DIRICHLET,
            BC_DIRICHLET,
            BC_NEUMANN,
            BC_NEUMANN,
            BC_JUMP_INTERFACE,
        ]

        # Set dirichlet data at left/right boundaries
        μ1 = 0.1
        μ2 = 0.2
        σ1 = 0.005
        σ2 = 0.005
        ue(x, y, _) = exp.(-(x .- μ1) .^ 2 ./ σ1 .- (y .- μ2) .^ 2 ./ 2σ2)
        ue_t(_...) = 0
        gDfun(_...) = 0
        gDdotfun(_...) = 0
        gNfun(_...) = 0
        body_force = nothing

        Nq = N + 1

        Np = Nq * Nq

        λ1(x, y) = 1
        λ2(x, y) = 1 / 2
        θ(x, y) = π * (2 - x) * (2 - y) / 4
        cxx(x, y) = cos(θ(x, y))^2 * λ1(x, y) + sin(θ(x, y))^2 * λ2(x, y)
        cxy(x, y) = cos(θ(x, y)) * sin(θ(x, y)) * (λ2(x, y) - λ1(x, y))
        cyy(x, y) = sin(θ(x, y))^2 * λ1(x, y) + cos(θ(x, y))^2 * λ2(x, y)

        (metrics, rhsops, EToDomain, EToF, EToO, EToS, FToB, FToE, FToLF) =
            build_square_circle(
                sbp_order,
                Nq,
                Nq,
                gDfun,
                gDdotfun,
                gNfun,
                body_force,
                bc_map,
                true;
                cxx = cxx,
                cxy = cxy,
                cyy = cyy,
            )

        FToB[findall(FToB .== BC_JUMP_INTERFACE)] .= BC_LOCKED_INTERFACE

        nelems = length(rhsops)

        # initial conditions
        lenq0 = 2 * Np + 2 * Nq + 2 * Nq
        q = Array{Float64, 1}(undef, nelems * lenq0)
        for e in 1:nelems
            (xf1, xf2, xf3, xf4) = metrics[e].facecoord[1]
            (yf1, yf2, yf3, yf4) = metrics[e].facecoord[2]

            u0 = ue.(metrics[e].coord[1][:], metrics[e].coord[2][:], e)
            v0 = ue_t.(metrics[e].coord[1][:], metrics[e].coord[2][:], e)
            û10 = ue.(xf1[:], yf1[:], e)
            û20 = ue.(xf2[:], yf2[:], e)
            û30 = ue.(xf3[:], yf3[:], e)
            û40 = ue.(xf4[:], yf4[:], e)

            q[(lenq0 * (e - 1) + 1):(e * lenq0)] =
                [u0; v0; û10; û20; û30; û40]
        end

        if N == maximum(Nvals)
            mkpath("output")
            write_vtk(
                @sprintf(
                    "output/locked_N_blocks_sim_step_2p_%d_N0_%04d_%04d",
                    sbp_order,
                    N,
                    0
                ),
                metrics,
                q;
                cxx = cxx,
                cxy = cxy,
                cyy = cyy,
            )
        end

        hmin = mapreduce(m -> m.hmin, min, values(metrics))
        dt = cfl * hmin

        energy = zeros(length(ts))
        energy[1] = compute_energy(q, rhsops, Np)

        # Parameters to pass to the ODE solver
        params = (
            Nqr = Nq,
            Nqs = Nq,
            rhsops = rhsops,
            EToF = EToF,
            EToO = EToO,
            EToS = EToS,
            FToB = FToB,
            FToE = FToE,
            FToLF = FToLF,
            friction = friction,
        )

        # Loop over times and advance simulation (saving at given times before
        # continuing)
        for step in 1:(length(ts) - 1)
            tspan = (ts[step], ts[step + 1])
            @show tspan
            timestep!(q, waveprop!, params, dt, tspan)
            if N == maximum(Nvals)
                write_vtk(
                    @sprintf(
                        "output/locked_N_blocks_sim_step_2p_%d_N0_%04d_%04d",
                        sbp_order,
                        N,
                        step
                    ),
                    metrics,
                    q;
                    cxx = cxx,
                    cxy = cxy,
                    cyy = cyy,
                )
            end
            energy[step + 1] = compute_energy(q, rhsops, Np)
        end
        open(
            @sprintf("locked_energy_2p_%d_N0_%04d.csv", sbp_order, N),
            "w",
        ) do io
            writedlm(
                io,
                [Array(ts) sqrt.(energy / energy[1]) 1 .-
                                                     sqrt.(energy / energy[1])],
            )
        end

        data[N] = (q = q, rhsops = rhsops)
    end

    nelems = length(data[Nvals[1]].rhsops)
    ϵ = zeros(length(Nvals) - 1)

    for lvl in 1:(length(Nvals) - 1)
        N1 = Nvals[lvl + 1]
        Nq1 = N1 + 1
        Np1 = Nq1^2

        N2 = Nvals[lvl]
        Nq2 = N2 + 1
        Np2 = Nq2^2

        q1 = reshape(data[N1].q, Np1 + 2Nq1, 2, nelems)
        q2 = reshape(data[N2].q, Np2 + 2Nq2, 2, nelems)

        rhsops = data[N1].rhsops

        for e in 1:nelems
            u1 = reshape(q1[1:Np1, 1, e], Nq1, Nq1)
            u2 = reshape(q2[1:Np2, 1, e], Nq2, Nq2)
            Δu = u1 - u2[1:2:end, 1:2:end]

            Mv = rhsops[e].JH

            ϵ[lvl] += Δu[:]' * Mv * Δu[:]
        end
        ϵ[lvl] = sqrt(ϵ[lvl])
    end
    @show ϵ
    rate = (log.(ϵ[2:end]) - log.(ϵ[1:(end - 1)])) / log(2)
    @show rate
end


text:
function loglikeam(data::Data, s::Dict{Symbol, Any}, beta::Bool)
  state = StateAM(s[:gammaC], s[:gammaT],
                  s[:etaC], s[:etaT], s[:mu], s[:sigma],
                  s[:nu], s[:phi])
  return loglike(state, data, beta)
end

function loglike(s::StateAM, data::Data, beta::Bool)
  return loglike_C(s, data) + loglike_T(s, data, beta)
end

function loglike_C(s::StateAM, data::Data)
  f = Util.mixskewtlogpdf(s.mu', s.sigma', s.nu', s.phi', s.etaC',
                          data.yC_finite, dims=2)
  return data.ZC * log(s.gammaC) + data.NC_finite * log1p(-s.gammaC) + sum(f)
end

function loglike_T(s::StateAM, data::Data, beta::Bool)
  f = Util.skewtlogpdf.(s.mu', s.sigma', s.nu', s.phi', data.yT_finite)

  if beta
    gammaTstar = s.gammaT
    etaTstar = s.etaT
  else
    gammaTstar = s.gammaC
    etaTstar = s.etaC
  end

  f = Util.mixskewtlogpdf(s.mu', s.sigma', s.nu', s.phi', etaTstar',
                          data.yT_finite, dims=2)

  ll_inf = data.ZT * log(gammaTstar)
  ll_finite = data.NT_finite * log1p(-gammaTstar) + sum(f)

  return ll_inf + ll_finite
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# From https://github.com/microsoft/Microsoft-MPI/blob/v10.0/src/include/mpi.h

const MPI_Aint   = Int
const MPI_Offset = Int64
const MPI_Count  = Int64

for T in [:MPI_Comm, :MPI_Info, :MPI_Win, :MPI_Request, :MPI_Op, :MPI_Datatype, :MPI_File]
    @eval begin
        primitive type $T 32 end
        $T(c::Cint) = reinterpret($T, c)
    end
end

const MPI_COMM_NULL  = reinterpret(Cint, 0x04000000)
const MPI_COMM_SELF  = reinterpret(Cint, 0x44000001)
const MPI_COMM_WORLD = reinterpret(Cint, 0x44000000)

const MPI_INFO_NULL = reinterpret(Cint, 0x1c000000)

const MPI_WIN_NULL = reinterpret(Cint, 0x20000000)

const MPI_OP_NULL = reinterpret(Cint, 0x18000000)
const MPI_MAX     = reinterpret(Cint, 0x58000001)
const MPI_MIN     = reinterpret(Cint, 0x58000002)
const MPI_SUM     = reinterpret(Cint, 0x58000003)
const MPI_PROD    = reinterpret(Cint, 0x58000004)
const MPI_LAND    = reinterpret(Cint, 0x58000005)
const MPI_BAND    = reinterpret(Cint, 0x58000006)
const MPI_LOR     = reinterpret(Cint, 0x58000007)
const MPI_BOR     = reinterpret(Cint, 0x58000008)
const MPI_LXOR    = reinterpret(Cint, 0x58000009)
const MPI_BXOR    = reinterpret(Cint, 0x5800000a)
const MPI_MINLOC  = reinterpret(Cint, 0x5800000b)
const MPI_MAXLOC  = reinterpret(Cint, 0x5800000c)
const MPI_REPLACE = reinterpret(Cint, 0x5800000d)
const MPI_NO_OP   = reinterpret(Cint, 0x5800000e)

const MPI_REQUEST_NULL = reinterpret(Cint, 0x2c000000)

const MPI_DATATYPE_NULL         = reinterpret(Cint, 0x0c000000)
const MPI_PACKED                = reinterpret(Cint, 0x4c00010f)
const MPI_CHAR                  = reinterpret(Cint, 0x4c000101)
const MPI_SIGNED_CHAR           = reinterpret(Cint, 0x4c000118)
const MPI_UNSIGNED_CHAR         = reinterpret(Cint, 0x4c000102)
const MPI_WCHAR                 = reinterpret(Cint, 0x4c00020e)
const MPI_BYTE                  = reinterpret(Cint, 0x4c00010d)
const MPI_SHORT                 = reinterpret(Cint, 0x4c000203)
const MPI_UNSIGNED_SHORT        = reinterpret(Cint, 0x4c000204)
const MPI_INT                   = reinterpret(Cint, 0x4c000405)
const MPI_UNSIGNED              = reinterpret(Cint, 0x4c000406)
const MPI_LONG                  = reinterpret(Cint, 0x4c000407)
const MPI_UNSIGNED_LONG         = reinterpret(Cint, 0x4c000408)
const MPI_FLOAT                 = reinterpret(Cint, 0x4c00040a)
const MPI_DOUBLE                = reinterpret(Cint, 0x4c00080b)

const MPI_INT8_T                = reinterpret(Cint, 0x4c000133)
const MPI_INT16_T               = reinterpret(Cint, 0x4c000234)
const MPI_INT32_T               = reinterpret(Cint, 0x4c000435)
const MPI_INT64_T               = reinterpret(Cint, 0x4c000836)
const MPI_UINT8_T               = reinterpret(Cint, 0x4c000137)
const MPI_UINT16_T              = reinterpret(Cint, 0x4c000238)
const MPI_UINT32_T              = reinterpret(Cint, 0x4c000439)
const MPI_UINT64_T              = reinterpret(Cint, 0x4c00083a)
const MPI_C_FLOAT_COMPLEX       = reinterpret(Cint, 0x4c000813)
const MPI_C_DOUBLE_COMPLEX      = reinterpret(Cint, 0x4c001014)

const MPI_FILE_NULL = Cint(0)

const MPI_PROC_NULL = Cint(-1)
const MPI_ANY_SOURCE = Cint(-2)
const MPI_ANY_TAG = Cint(-1)
const MPI_SUCCESS = Cint(0)
const MPI_UNDEFINED = Cint(-32766)
const MPI_LOCK_EXCLUSIVE = Cint(234)
const MPI_LOCK_SHARED = Cint(235)
const MPI_MAX_INFO_KEY = Cint(255)
const MPI_MAX_INFO_VAL = Cint(1024)
const MPI_TAG_UB = reinterpret(Cint, 0x64400001)
const MPI_COMM_TYPE_SHARED = Cint(1)
const MPI_ORDER_C = Cint(56)
const MPI_ORDER_FORTRAN = Cint(57)
const MPI_UNIVERSE_SIZE = reinterpret(Cint, 0x64400009)
const MPI_MAX_ERROR_STRING = Cint(512)
const MPI_SUCCESS = Cint(0)
const MPI_MODE_RDONLY = reinterpret(Cint, 0x00000002)
const MPI_MODE_RDWR = reinterpret(Cint, 0x00000008)
const MPI_MODE_WRONLY = reinterpret(Cint, 0x00000004)
const MPI_MODE_CREATE = reinterpret(Cint, 0x00000001)
const MPI_MODE_EXCL = reinterpret(Cint, 0x00000040)
const MPI_MODE_DELETE_ON_CLOSE = reinterpret(Cint, 0x00000010)
const MPI_MODE_UNIQUE_OPEN = reinterpret(Cint, 0x00000020)
const MPI_MODE_SEQUENTIAL = reinterpret(Cint, 0x00000100)
const MPI_MODE_APPEND = reinterpret(Cint, 0x00000080)

const MPI_BOTTOM = reinterpret(SentinelPtr, 0)
const MPI_IN_PLACE = reinterpret(SentinelPtr, -1)

const MPI_STATUS_IGNORE = reinterpret(SentinelPtr, 1)
const MPI_STATUSES_IGNORE = reinterpret(SentinelPtr, 1)

struct Status
    _pad1::Cint
    _pad2::Cint
    source::Cint
    tag::Cint
    error::Cint
end
const STATUS_EMPTY = Status(zero(Cint), zero(Cint), MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_SUCCESS)


text:
"""
Some statistical functions, mainly related to (weighted) median computations.
"""
module Stats

import ..GraphIdx: Const, Vec, Ones

"""
    permcumsum(w, pi)

Return the cumulative sums over `w` in order `pi`.
"""
function permcumsum(w::Array{W}, pi::Vector{I}) where {W, I<:Integer}
    local acc::W = 0
    for j in pi
        acc += w[j]
        w[j] = acc
    end
    return w
end

permcumsum(w::Const{W}, pi::Vector{I}) where {W, I<:Integer} =
    pi

permcumsum(w::Ones{W}, pi::Vector{I}) where {W, I<:Integer} =
    pi

"""
    weighted_median(x, [w, [∂ ]])

Compute a generalized weighted order statistics.
Hereby, we minimize
```math
\\min_{t} \\sum_{i} w_i |t - x_i| + f(t)
```
whereby `∂` is the derivative of the convex function ``f``
(by default constant `0` for ordinary weighted median).

Return the interval ``[a, b]`` containing the minimizer (if unique ``a = b``).

!!! todo

    Compare the sorting to binary search and median strategies.
    Maybe use [SIMD for median](https://gist.github.com/jw3126/3a3c65009e96af2c4fcc96f701bd4913).

"""
weighted_median(x::Vector{X}) where {X} =
    weighted_median(x, Ones{X}())

function weighted_median(
    x::Vector{X},
    ww::W,
    ∂::F = _ -> 0,
    pi::Vector{Int} = Int[],
)::Tuple{X,X} where {X, W, F <: Function}
    if ww isa Array
        @assert length(x) == length(ww)
    end
    local n = length(x)
    resize!(pi, n)
    sortperm!(pi, x)

    # Precompute  wsum(i) = sum(-w[j] for j < i) + sum(+w[j] for j > i) + ∂(x[i])
    local w = permcumsum(ww, pi)

    local wsum0 = -w[pi[end]]
    w0(i::Int) = i == 1 ? w[pi[1]] : w[pi[i]] - w[pi[i-1]]
    wsum(i::Int) = wsum0 + ∂(x[pi[i]]) + w[pi[i]] + (i == 1 ? 0 : w[pi[i-1]])

    # sum of abs "derivative" if t < x[first] (and x is sorted)
    local first::Int, last::Int = 1, length(x)
    while last - first > 1
        local mid::Int = (first + last) ÷ 2
        if wsum(mid) < 0
            first = mid
        else
            last = mid
        end
    end
    # @show first, last, wsum(first), wsum(last)
    if wsum(first) + w0(first) > 0
        return x[pi[first]], x[pi[first]]
    elseif wsum(last) - w0(last) < 0
        return x[pi[last]], x[pi[last]]
    else
        return x[pi[first]], x[pi[last]]
    end
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#


text:
""" Serialize abstract wiring diagrams as GraphML.

Serialization of box, port, and wire values can be overloaded by data type
(see `convert_to_graphml_data` and `convert_from_graphml_data`).

GraphML is the closest thing to a de jure and de facto standard in the space of
graph data formats, supported by a variety of graph applications and libraries.
We depart mildly from the GraphML spec by allowing JSON data attributes for
GraphML nodes, ports, and edges.

References:

- GraphML Primer: http://graphml.graphdrawing.org/primer/graphml-primer.html
- GraphML DTD: http://graphml.graphdrawing.org/specification/dtd.html
"""
module GraphMLWiringDiagrams
export read_graphml, parse_graphml, write_graphml, generate_graphml,
  convert_from_graphml_data, convert_to_graphml_data

using Compat
using DataStructures: OrderedDict
import JSON
using LightXML

using ...CategoricalAlgebra.Graphs
using ..DirectedWiringDiagrams, ..WiringDiagramSerialization
import ..DirectedWiringDiagrams: PortData

# Data types
############

struct GraphMLKey
  id::String
  attr_name::String
  attr_type::String
  scope::String
  default::Union{Some,Nothing}
end
GraphMLKey(id::String, attr_name::String, attr_type::String, scope::String) =
  GraphMLKey(id, attr_name, attr_type, scope, nothing)

struct WriteState
  keys::OrderedDict{Tuple{String,String},GraphMLKey}
  WriteState() = new(OrderedDict{Tuple{String,String},GraphMLKey}())
end

struct ReadState
  keys::OrderedDict{String,GraphMLKey}
  BoxValue::Type
  PortValue::Type
  WireValue::Type
end

# Serialization
###############

""" Write a wiring diagram to a file as GraphML.
"""
function write_graphml(diagram::WiringDiagram, filename::String)
  LightXML.save_file(generate_graphml(diagram), filename)
end

""" Generate GraphML representing a wiring diagram.
"""
function generate_graphml(diagram::WiringDiagram)::XMLDocument
  # Create XML document.
  xdoc = XMLDocument()
  finalizer(free, xdoc) # Destroy all children when document is GC-ed.
  xroot = create_root(xdoc, "graphml")
  set_attributes(xroot, Pair[
    "xmlns" => "http://graphml.graphdrawing.org/xmlns",
    "xmlns:xsi" => "http://www.w3.org/2001/XMLSchema-instance",
    "xsi:schemaLocation" => "http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd"
  ])

  # Create top-level graph element.
  xgraph = new_element("graph")
  set_attribute(xgraph, "edgedefault", "directed")

  # Recursively create nodes.
  state = WriteState()
  generate_graphml_node(state, xgraph, diagram, Int[])

  # Add attribute keys (data declarations). The keys are collected while
  # writing the nodes and are stored in the state object.
  for key in values(state.keys)
    generate_graphml_key(xroot, key)
  end

  add_child(xroot, xgraph)
  return xdoc
end

function generate_graphml_node(
    state::WriteState, xgraph::XMLElement, diagram::WiringDiagram, path::Vector{Int})
  # Create node element for wiring diagram and graph subelement to contain
  # boxes and wires.
  xnode = new_child(xgraph, "node")
  set_attribute(xnode, "id", box_id(path))
  generate_graphml_ports(state, xnode, diagram)

  xsubgraph = new_child(xnode, "graph")
  set_attribute(xsubgraph, "id", string(box_id(path), ":graph"))
  generate_graphml_data(state, xsubgraph, "graph", diagram.value)

  # Add node elements for boxes.
  for v in box_ids(diagram)
    generate_graphml_node(state, xsubgraph, box(diagram, v), [path; v])
  end

  # Add edge elements for wires.
  for (i, wire) in enumerate(wires(diagram))
    xedge = new_child(xsubgraph, "edge")
    set_attributes(xedge, Pair[
      "id"         => wire_id(path, i),
      "source"     => box_id(diagram, [path; wire.source.box]),
      "sourceport" => port_name(diagram, wire.source),
      "target"     => box_id(diagram, [path; wire.target.box]),
      "targetport" => port_name(diagram, wire.target),
    ])
    generate_graphml_data(state, xedge, "edge", wire.value)
  end
end

function generate_graphml_node(
    state::WriteState, xgraph::XMLElement, box::Box, path::Vector{Int})
  xnode = new_child(xgraph, "node")
  set_attribute(xnode, "id", box_id(path))
  generate_graphml_data(state, xnode, "node", box.value)
  generate_graphml_ports(state, xnode, box)
end

function generate_graphml_ports(state::WriteState, xnode::XMLElement, box::AbstractBox)
  # Write input ports.
  for (i, port) in enumerate(input_ports(box))
    xport = new_child(xnode, "port")
    set_attribute(xport, "name", port_name(InputPort, i))
    generate_graphml_data(state, xport, "port", Dict("portkind" => "input"))
    generate_graphml_data(state, xport, "port", port)
  end
  # Write output ports.
  for (i, port) in enumerate(output_ports(box))
    xport = new_child(xnode, "port")
    set_attribute(xport, "name", port_name(OutputPort, i))
    generate_graphml_data(state, xport, "port", Dict("portkind" => "output"))
    generate_graphml_data(state, xport, "port", port)
  end
end

function generate_graphml_key(xroot::XMLElement, key::GraphMLKey)
  xkey = new_child(xroot, "key")
  set_attributes(xkey, Pair[
    "id" => key.id,
    "for" => key.scope,
    "attr.name" => key.attr_name,
    "attr.type" => key.attr_type,
  ])
  if !isnothing(key.default)
    xdefault = new_child(xkey, "default")
    set_content(xdefault, generate_graphml_data_value(something(key.default)))
  end
end

function generate_graphml_data(state::WriteState, xelem::XMLElement, scope::String, value)
  data = convert_to_graphml_data(value)
  for (attr_name, attr_value) in data
    # Retrieve or create key from state object.
    key = get!(state.keys, (attr_name, scope)) do
      nkeys = length(state.keys)
      id = "d$(nkeys+1)"
      attr_type = generate_graphml_data_type(typeof(attr_value))
      GraphMLKey(id, attr_name, attr_type, scope)
    end

    # Write attribute data to <key> element.
    xdata = new_child(xelem, "data")
    set_attribute(xdata, "key", key.id)
    xvalue = generate_graphml_data_value(attr_value)
    (xvalue isa AbstractString ? set_content : add_child)(xdata, xvalue)
  end
end

generate_graphml_data_type(::Type{Bool}) = "boolean"
generate_graphml_data_type(::Type{<:Integer}) = "int"
generate_graphml_data_type(::Type{<:Real}) = "double"
generate_graphml_data_type(::Type{String}) = "string"
generate_graphml_data_type(::Type{Symbol}) = "string"
generate_graphml_data_type(::Type{Dict{String,T}}) where T = "json"
generate_graphml_data_type(::Type{Vector{T}}) where T = "json"

generate_graphml_data_value(x::Number) = string(x)
generate_graphml_data_value(x::String) = x
generate_graphml_data_value(x::Symbol) = string(x)
generate_graphml_data_value(x::Dict) = JSON.json(x)
generate_graphml_data_value(x::Vector) = JSON.json(x)

convert_to_graphml_data(x) = convert_to_graph_data(x)

# Deserialization
#################

""" Read a wiring diagram from a GraphML file.
"""
function read_graphml(
    BoxValue::Type, PortValue::Type, WireValue::Type, filename::String)
  parse_graphml(BoxValue, PortValue, WireValue, LightXML.parse_file(filename))
end

""" Parse a wiring diagram from a GraphML string or XML document.
"""
function parse_graphml(
    BoxValue::Type, PortValue::Type, WireValue::Type, s::AbstractString)
  parse_graphml(BoxValue, PortValue, WireValue, LightXML.parse_string(s))
end
function parse_graphml(::Type{BoxValue}, ::Type{PortValue}, ::Type{WireValue},
    xdoc::XMLDocument)::WiringDiagram where {BoxValue, PortValue, WireValue}
  xroot = root(xdoc)
  @assert name(xroot) == "graphml" "Root element of GraphML document must be <graphml>"
  xgraphs = xroot["graph"]
  @assert length(xgraphs) == 1 "Root element of GraphML document must contain exactly one <graph>"
  xgraph = xgraphs[1]
  xnodes = xgraph["node"]
  @assert length(xnodes) == 1 "Root graph of GraphML document must contain exactly one <node>"
  xnode = xnodes[1]

  keys = parse_graphml_keys(xroot)
  state = ReadState(keys, BoxValue, PortValue, WireValue)
  diagram, ports = parse_graphml_node(state, xnode)
  diagram
end

function parse_graphml_node(state::ReadState, xnode::XMLElement)
  # Parse all port elements.
  ports, input_ports, output_ports = parse_graphml_ports(state, xnode)

  # Handle special cases: atomic boxes and malformed elements.
  xgraphs = xnode["graph"]
  if length(xgraphs) > 1
    error("Node element can contain at most one <graph> (subgraph element)")
  elseif isempty(xgraphs)
    data = parse_graphml_data(state, xnode)
    value = convert_from_graphml_data(state.BoxValue, data)
    return (Box(value, input_ports, output_ports), ports)
  end
  xgraph = xgraphs[1]

  # If we get here, we're reading a wiring diagram.
  # FIXME: We should not assume that diagram data has same type as box data.
  data = parse_graphml_data(state, xgraph)
  value = isempty(data) ? nothing : convert_from_graph_data(state.BoxValue, data)
  diagram = WiringDiagram(value, input_ports, output_ports)
  all_ports = Dict{Tuple{String,String},Port}()
  for (key, port_data) in ports
    all_ports[key] = port_data.kind == InputPort ?
      Port(input_id(diagram), OutputPort, port_data.port) :
      Port(output_id(diagram), InputPort, port_data.port)
  end

  # Read the node elements.
  for xsubnode in xgraph["node"]
    box, subports = parse_graphml_node(state, xsubnode)
    v = add_box!(diagram, box)
    for (key, port_data) in subports
      all_ports[key] = Port(v, port_data.kind, port_data.port)
    end
  end

  # Read the edge elements.
  for xedge in xgraph["edge"]
    data = parse_graphml_data(state, xedge)
    value = convert_from_graphml_data(state.WireValue, data)
    xsource = attribute(xedge, "source", required=true)
    xtarget = attribute(xedge, "target", required=true)
    xsourceport = attribute(xedge, "sourceport", required=true)
    xtargetport = attribute(xedge, "targetport", required=true)
    source = all_ports[(xsource, xsourceport)]
    target = all_ports[(xtarget, xtargetport)]
    add_wire!(diagram, Wire(value, source, target))
  end

  (diagram, ports)
end

function parse_graphml_ports(state::ReadState, xnode::XMLElement)
  ports = Dict{Tuple{String,String},PortData}()
  input_ports, output_ports = state.PortValue[], state.PortValue[]
  xnode_id = attribute(xnode, "id", required=true)
  xports = xnode["port"]
  for xport in xports
    xport_name = attribute(xport, "name", required=true)
    data = parse_graphml_data(state, xport)
    port_kind = pop!(data, "portkind")
    value = convert_from_graphml_data(state.PortValue, data)
    if port_kind == "input"
      push!(input_ports, value)
      ports[(xnode_id, xport_name)] = PortData(InputPort, length(input_ports))
    elseif port_kind == "output"
      push!(output_ports, value)
      ports[(xnode_id, xport_name)] = PortData(OutputPort, length(output_ports))
    else
      error("Invalid port kind: $portkind")
    end
  end
  (ports, input_ports, output_ports)
end

function parse_graphml_keys(xroot::XMLElement)
  keys = OrderedDict{String,GraphMLKey}()
  for xkey in xroot["key"]
    # Read attribute ID, name, type, and scope. Both name and type are optional.
    attrs = attributes_dict(xkey)
    id = attrs["id"]
    attr_name = get(attrs, "attr.name", id)
    attr_type = get(attrs, "attr.type", "string")
    scope = get(attrs, "for", "all")

    # Read attribute default value.
    xdefaults = xkey["default"]
    default = if isempty(xdefaults)
      nothing
    else
      @assert length(xdefaults) == 1 "GraphML key can have at most one <default>"
      xdefault = xdefaults[1]
      Some(parse_graphml_data_value(Val{Symbol(attr_type)}, content(xdefault)))
    end

    keys[id] = GraphMLKey(id, attr_name, attr_type, scope, default)
  end
  keys
end

function parse_graphml_data(keys::AbstractDict, xelem::XMLElement)
  # FIXME: We are not using the default values for the keys.
  data = Dict{String,Any}()
  for xdata in xelem["data"]
    xkey = attribute(xdata, "key", required=true)
    key = keys[xkey]
    data[key.attr_name] = parse_graphml_data_value(
      Val{Symbol(key.attr_type)}, xdata)
  end
  data
end
parse_graphml_data(state::ReadState, xelem::XMLElement) =
  parse_graphml_data(state.keys, xelem)

parse_graphml_data_value(type::Type, xdata::XMLElement) =
  parse_graphml_data_value(type, content(xdata))
parse_graphml_data_value(::Type{Val{:boolean}}, s::String) = parse(Bool, lowercase(s))
parse_graphml_data_value(::Type{Val{:int}}, s::String) = parse(Int, s)
parse_graphml_data_value(::Type{Val{:long}}, s::String) = parse(Int, s)
parse_graphml_data_value(::Type{Val{:float}}, s::String) = parse(Float32, s)
parse_graphml_data_value(::Type{Val{:double}}, s::String) = parse(Float64, s)
parse_graphml_data_value(::Type{Val{:string}}, s::String) = s
parse_graphml_data_value(::Type{Val{:json}}, s::String) = JSON.parse(s)

convert_from_graphml_data(T::Type, data) = convert_from_graph_data(T, data)

# Graphs
########

""" Read a property graph from a GraphML file.

The equivalent of NetworkX's `read_graphml` function.
"""
function read_graphml_metagraph(filename::String; kw...)
  parse_graphml_metagraph(LightXML.parse_file(filename); kw...)
end

""" Parse a property graph from a GraphML string or XML document.

The equivalent of NetworkX's `parse_graphml` function.
"""
function parse_graphml_metagraph(s::AbstractString; kw...)
  parse_graphml_metagraph(LightXML.parse_string(s); kw...)
end
function parse_graphml_metagraph(xdoc::XMLDocument;
                                 directed::Bool=false)::AbstractPropertyGraph
  xroot = root(xdoc)
  @assert name(xroot) == "graphml" "Root element of GraphML document must be <graphml>"
  xgraphs = xroot["graph"]
  @assert length(xgraphs) == 1 "Root element of GraphML document must contain exactly one <graph>"
  xgraph = xgraphs[1]

  # Read the GraphML keys.
  keys = parse_graphml_keys(xroot)
  parse_props = xnode::XMLElement -> Dict{Symbol,Any}(
    Symbol(k) => v for (k,v) in parse_graphml_data(keys, xnode))

  # Create attributed graph.
  graph = if directed || attribute(xgraph, "edgedefault") == "directed"
    PropertyGraph{Any}()
  else
    SymmetricPropertyGraph{Any}()
  end
  set_gprops!(graph, parse_props(xgraph))

  # Read the node elements.
  vertices = Dict{String,Int}()
  for (i, xnode) in enumerate(xgraph["node"])
    node_id = attribute(xnode, "id", required=true)
    vertices[node_id] = i
    add_vertex!(graph, parse_props(xnode))
  end

  # Read the edge elements.
  for xedge in xgraph["edge"]
    source_id = attribute(xedge, "source", required=true)
    target_id = attribute(xedge, "target", required=true)
    source, target = vertices[source_id], vertices[target_id]
    add_edge!(graph, source, target, parse_props(xedge))
  end

  graph
end

end


text:
export AkersNicholls_fast,AkersNicholls

"""
    AkersNicholls_fast(param;kwargs)

Same as `AkersNicholls`, but faster.
"""
mutable struct AkersNicholls_fast <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info    :: String

    function AkersNicholls_fast( param::NamedTuple;
									mesh = Mesh(param),
									IL	    = false,
									dealias=false,
									label="Akers-Nicholls" )

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used thereafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "Akers-Nicholls model.\n"
		if IL == true
			info *= "├─Steepness parameter ϵ=$ϵ (infinite depth case).\n"
		else
			info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
					scaling parameter ν=$nu.\n"
		end
		if dealias == true || dealias == 1
			info *= "└─Dealiasing with Orszag’s 3/2 rule. "
		else
			info *= "└─No dealiasing. "
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data
		x = mesh.x
		k = mesh.k
		if IL == true
			Tμ 	= -1im * sign.(k)
			Lμ 	= abs.(k);Lμ[1]=1;
		else
			Tμ 	= -1im* sign.(k) .* tanh.(sqrt(μ)*abs.(k))
			Lμ 	= sqrt(μ)*abs.(k)./tanh.(sqrt(μ)*abs.(k));Lμ[1]=1;
		end
		∂ₓ	=  1im * sqrt(μ)* k            # Differentiation
		Tμ∂ₓ= Tμ .* ∂ₓ
		if dealias == true || dealias == 1
			Π⅔    = abs.(k) .< (mesh.kmax-mesh.kmin)/3 	# Dealiasing low-pass filter
		else
			Π⅔    = zero(k) .+ 1     		# No dealiasing (Π⅔=Id)
		end

        ζ = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

        Px  = plan_fft(ζ; flags = FFTW.MEASURE)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		    ldiv!(ζ, Px , view(U,:,1))

		    I₁  .= view(U,:,2) .* Lμ
		    ldiv!(unew, Px , I₁)
		    unew .^= 2
		    mul!(I₁, Px , unew)
		    I₁ .*= Tμ

		    I₂  .= view(U,:,1) .* ∂ₓ
		    ldiv!(unew, Px , I₂)
		    unew .*= ζ
		    mul!(I₂, Px , unew)

		    I₃  .= view(U,:,1) .* Tμ∂ₓ
		    ldiv!(unew, Px, I₃)
		    unew .*= ζ
		    mul!(I₃ , Px , unew)
		    I₃ .*= Tμ

		    ζ  .= .- view(U,:,2) .* ∂ₓ

		    I₁ .*= ν
			I₁ .-= I₂
		    I₁ .-= I₃
		    I₁ .*= Π⅔
		    I₁ .*= ϵ

		    U[:,2] .= (view(U,:,1) .* Tμ .+ I₁)/sqrt(μ)/ν
		    U[:,1] .= ζ/sqrt(μ)

		end

		# Build raw data from physical data.
		function mapto(data::InitialData)
			ζ.=data.η(x);v=data.v(x);
			fftv=fft(v);
			[Π⅔ .* fft(ζ)  Π⅔ .* (fftv./Lμ+ϵ*fft(ζ.*v)+ ϵ*Tμ.*fft(ζ.*ifft(Tμ.*fftv)))/ν]
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U;n=10)
			ζ.=ifft(view(U,:,1));I₁.=ν*view(U,:,2);I₂.=Lμ.*I₁;
			for j=1:n
				I₂.=Lμ.*(I₁-ϵ*Π⅔ .* ( fft(ζ.*ifft(I₂)) + Tμ.*fft(ζ.*ifft(Tμ.*I₂))) )
			end
			real(ζ),real(ifft(I₂)),mesh.x
		end

		new(label, f!, mapto, mapfro, info )

    end
end

"""
    AkersNicholls(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the quadratic deep-water model proposed by [Akers and Nicholls](https://doi.org/10.1137/090771351)
and [Cheng, Granero-Belinchón, Shkoller and Milewski](https://doi.org/10.1007/s42286-019-00005-w)

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- optionally, `ν` the shallow/deep water scaling factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0`, or `μ=Inf`.
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `IL`: Set the infinite-layer case if `IL=true` (or `μ=Inf`, or `ν=0`), in which case `ϵ` is the steepness parameter. Default is `false`.
- `dealias`: dealiasing with `1/3` Orlicz rule if `true` or no dealiasing if `false` (by default);
- `label`: a label for future references (default is `"deep quadratic"`);

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `AkersNicholls.f!` to be called in explicit time-integration solvers;
2. a function `AkersNicholls.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `AkersNicholls.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct AkersNicholls <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info 	:: String

    function AkersNicholls( param::NamedTuple;
							mesh = Mesh(param),
							IL	    = false,
							dealias = false,
							label="Akers-Nicholls")

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used thereafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "Akers-Nicholls model.\n"
		if IL == true
			info *= "├─Steepness parameter ϵ=$ϵ (infinite depth case).\n"
		else
			info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
					scaling parameter ν=$nu.\n"
		end
		if dealias == true || dealias == 1
			info *= "└─Dealiasing with Orszag’s 3/2 rule. "
		else
			info *= "└─No dealiasing. "
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data
		x = mesh.x
		k = mesh.k
		if IL == true
			Tμ 	= -1im * sign.(k)
			Lμ 	= abs.(k);Lμ[1]=1;
		else
			Tμ 	= -1im* sign.(k) .* tanh.(sqrt(μ)*abs.(k))
			Lμ 	= sqrt(μ)*abs.(k)./tanh.(sqrt(μ)*abs.(k));Lμ[1]=1;
		end
		∂ₓ	=  1im * sqrt(μ)* k            # Differentiation
		Tμ∂ₓ= Tμ .* ∂ₓ

		if dealias == true || dealias == 1
			Π⅔    = abs.(k) .< (mesh.kmax-mesh.kmin)/3 	# Dealiasing low-pass filter
		else
			Π⅔    = zero(k) .+ 1     		# No dealiasing (Π⅔=Id)
		end

        ζ = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

			ζ .=ifft(U[:,1]);
			I₁ .=Tμ.*fft(ifft(Lμ.*U[:,2]).^2);
			# Tricomi identity : the above equals the two commented lines below
			# I₁ .=Tμ.*fft(ifft(Lμ.*U[:,2]).^2 .+ ifft(∂ₓ.*U[:,2]).^2)/2;
			# I₁ .-=fft(ifft(Lμ.*U[:,2]).*ifft(∂ₓ.*U[:,2]));

			I₂ .=fft(ζ.*ifft(∂ₓ.*U[:,1])) ;
			I₃ .=Tμ.*fft(ζ.*ifft(Tμ∂ₓ.*U[:,1]));
			ζ .= U[:,1] ;
			U[:,1] .= -∂ₓ.*U[:,2]/sqrt(μ) ;
			U[:,2] .= (Tμ.*ζ+ϵ*Π⅔.*(ν*I₁-I₂-I₃))/sqrt(μ)/ν ;

		end

		# Build raw data from physical data.
		function mapto(data::InitialData)
			ζ.=data.η(x);v=data.v(x);
			fftv=fft(v);
			[Π⅔ .* fft(ζ)  Π⅔ .* (fftv./Lμ+ϵ*fft(ζ.*v)+ ϵ*Tμ.*fft(ζ.*ifft(Tμ.*fftv)))/ν]
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U;n=10)
			ζ.=ifft(U[:,1]);I₁.=ν*U[:,2];I₂.=Lμ.*I₁;
			for j=1:n
				I₂.=Lμ.*(I₁-ϵ*Π⅔ .* ( fft(ζ.*ifft(I₂)) + Tμ.*fft(ζ.*ifft(Tμ.*I₂))) )
			end
			real(ζ),real(ifft(I₂)),mesh.x
		end

		new(label, f!, mapto, mapfro, info )

    end
end


text:
module SplitExplicit

using StaticArrays

using ..Ocean
using ..HydrostaticBoussinesq
using ..ShallowWater

using ...VariableTemplates
using ...MPIStateArrays
using ...Mesh.Geometry
using ...DGMethods
using ...BalanceLaws

import ...BalanceLaws:
    initialize_states!,
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

HBModel = HydrostaticBoussinesqModel
SWModel = ShallowWaterModel

function initialize_states!(
    ::HBModel{C},
    ::SWModel{C},
    _...,
) where {C <: Uncoupled}
    return nothing
end
function tendency_from_slow_to_fast!(
    ::HBModel{C},
    ::SWModel{C},
    _...,
) where {C <: Uncoupled}
    return nothing
end
function cummulate_fast_solution!(
    ::HBModel{C},
    ::SWModel{C},
    _...,
) where {C <: Uncoupled}
    return nothing
end
function reconcile_from_fast_to_slow!(
    ::HBModel{C},
    ::SWModel{C},
    _...,
) where {C <: Uncoupled}
    return nothing
end

include("VerticalIntegralModel.jl")
include("Communication.jl")
include("ShallowWaterCoupling.jl")
include("HydrostaticBoussinesqCoupling.jl")

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module ParserCombinator

using AutoHashEquals
import Base: iterate, getindex, isless, size, hash,
axes, lastindex, firstindex
import Base: ==, ~, +, &, |, >=, >, |>, !
using Printf: @printf

export Matcher,
diagnostic, forwards, LineSource, LineIter,
Config, Cache, NoCache, make, make_all, make_one, once,
parse_one, parse_one_cache, parse_one_nocache,
parse_all, parse_all_cache, parse_all_nocache,
parse_lines, parse_lines_cache,
Debug, Trace,
parse_dbg, parse_one_dbg, parse_one_cache_dbg, parse_one_nocache_dbg,
parse_all_dbg, parse_all_cache_dbg, parse_all_nocache_dbg,
parse_lines_dbg, parse_lines_cache_dbg,
Success, EMPTY, Failure, FAILURE, Execute,
State, Clean, CLEAN, Dirty, DIRTY,
ParserException, Value, Empty, EMPTY, Delegate, DelegateState,
Epsilon, Insert, Dot, Fail, Drop, Equal,
Repeat, Depth, Breadth, Depth!, Breadth!, ALL,
Series, Seq, And, Seq!, And!, Alt, Alt!, Lookahead, Not, Pattern, Delayed, Eos,
ParserError, Error,
Transform, App, Appl, ITransform, IApp, IAppl,
@p_str, @P_str, @e_str, @E_str, Opt, Opt!,
Parse, PUInt, PUInt8, PUInt16, PUInt32, PUInt64,
PInt, PInt8, PInt16, PInt32, PInt64, PFloat32, PFloat64,
Word, Space,
Star, Plus, Star!, Plus!, StarList, StarList!, PlusList, PlusList!,
@with_names, set_name,
@with_pre, @with_post, set_fix,
TrySource, Try, parse_try, parse_try_dbg, parse_try_cache, parse_try_cache_dbg,
Parsers, axes, lastindex

include("core/types.jl")
include("core/sources.jl")
include("core/names.jl")
include("core/matchers.jl")
include("core/parsers.jl")
include("core/debug.jl")
include("core/transforms.jl")
include("core/sugar.jl")
include("core/extras.jl")
include("core/print.jl")
include("core/fix.jl")
include("core/try.jl")

include("Parsers.jl")

end


text:
# used by data.jl and multinetwork.jl
function compare_dict(d1, d2)
    for (k1,v1) in d1
        if !haskey(d2, k1)
            #println(k1)
            return false
        end
        v2 = d2[k1]

        if isa(v1, Number)
            if isnan(v1)
                #println("1.1")
                if !isnan(v2)
                    #println(v1, " ", v2)
                    return false
                end
            else
                #println("1.2")
                if !isapprox(v1, v2)
                    #println(v1, " ", v2)
                    return false
                end
            end
        elseif isa(v1, Dict)
            if !compare_dict(v1, v2)
                #println(v1, " ", v2)
                return false
            end
        else
            #println("2")
            if v1 != v2
                #println(v1, " ", v2)
                return false
            end
        end
    end
    return true
end


text:
macro cxx_str(exprs...) return clang_str(__module__, __source__, :cxx, exprs...) end
macro cxx_cmd(exprs...) return clang_cmd(__module__, __source__, :cxx, exprs...) end

Base.String(::Type{Context{:cxx}}) = "C++"

# TODO: implement me!


text:
import CausalityToolsBase:
    BinningScheme,
    RectangularBinning

import ..BinningSchemes:
    InducedRectangularBinning,
    AverageRectangularBinning

"""
    invariantmeasure(data, binning_scheme::BinningScheme; kwargs...)

Partition `data` according to the the given `binning_scheme` and compute the invariant
measure over the partition elements.
"""
function invariantmeasure(data, binning_scheme::BinningScheme; kwargs...)

end

######################################
# Triangulation binnign schemes
######################################

"""
    invariantmeasure(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ExactIntersection) -> TriangulationExactInvariantMeasure

Estimate the invariant measure over the state space defined by `pts` using a triangulation
of the phase space as the partition, using exact simplex intersections to compute transition
probabilities between the states (simplices).

# Example

Assume we have sufficiently few points that a triangulation approach involving simplex
intersections is computationally feasible to compute the invariant measure. Then
a transfer operator can be computed as follows.

```julia
pts = [rand(3) for i = 1:30]
invariantmeasure(pts, TriangulationBinning(), ExactIntersection())
```
"""
function invariantmeasure(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ExactIntersection)
    triangulationexactinvariantmeasure(pts)
end

"""
    invariantmeasure(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ApproximateIntersection;
        n::Int = 100, sample_randomly::Bool = false) -> TriangulationApproxInvariantMeasure

Estimate the invariant measure over the state space defined by `pts` using a triangulation
of the phase space as the partition, using exact simplex intersections to compute transition
probabilities between the states (simplices).

# Example

Assume we have sufficiently few points that a triangulation approach involving simplex
intersections is computationally feasible to compute the invariant measure. Then
a transfer operator can be computed as follows.

```julia
pts = [rand(3) for i = 1:30]
invariantmeasure(pts, TriangulationBinning(), ApproximateIntersection())
```
"""
function invariantmeasure(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ApproximateIntersection;
        n::Int = 100, sample_randomly::Bool = false)
    triangulationapproxinvariantmeasure(pts, n_sample_pts = n, sample_randomly = sample_randomly)
end

######################################
# Providing customized binning schemes
######################################
# Regular rectnagular invariant measure

"""
    invariantmeasure(data, binning_scheme::RectangularBinning; kwargs...)

Partition `data` according to the the given `binning_scheme` and compute the invariant
measure over the partition elements.

# Example

Assume we have enough points that a rectangular partition yields a good estimate of the
invariant measure. Then the measure over the partition can be computed as follows.

```julia
pts = [rand(3) for i = 1:2000]

# Use rectangular boxes constructed by subdividing each coordinate
# axis into 10 subintervals of equal length.
binning_scheme = RectangularBinning(10)
invariantmeasure(pts, binning_scheme)
```
"""
function invariantmeasure(data, binning_scheme::RectangularBinning; kwargs...)
    rectangularinvariantmeasure(
        data,
        binning_scheme,
        kwargs...)
end

# Induced measures
function invariantmeasure(data, binning_scheme::InducedRectangularBinning; kwargs...)
    inducedrectangularinvariantmeasure(
        data,
        binning_scheme.ϵ_target_partition,
        binning_scheme.ϵ_source_partition,
        kwargs...)
end

# If you don't want to construct an induced invariant measure.
function invariantmeasure(data,
    ϵ_target_partition::RectangularBinning,
    ϵ_source_partition::RectangularBinning;
    kwargs...)

inducedrectangularinvariantmeasure(data, ϵ_target_partition.ϵ, ϵ_source_partition.ϵ,
    kwargs...)
end

function invariantmeasure(data, binning_scheme::AverageRectangularBinning; kwargs...)
    averagerectangularinvariantmeasure(
        data,
        binning_scheme.ϵ_target_partition,
        binning_scheme.ϵs_source_partitions,
        kwargs...)
end

function invariantmeasure(data,
    ϵ_target_partition::RectangularBinning,
    ϵs_source_partitions::Vector{RectangularBinning};
    kwargs...)

    averagerectangularinvariantmeasure(
        data,
        ϵ_target_partition.ϵ,
        [ϵs_source_partitions[i].ϵ for i in 1:length(ϵs_source_partitions)],
        kwargs...)
end

export invariantmeasure


text:
#=

MCMC for biogeographic competition model

Ignacio Quintero Mächler

t(-_-t)

April 27 2017

=#

"""
    tribe_mcmc(...)

Run MCMC for trait and range interspecitific biogeographic evolution (tribe)
**under the prior**.
"""
function tribe_mcmc(out_file::String;
                    niter   ::Int64             = 50_000,
                    nthin   ::Int64             = 100,
                    nburn   ::Int64             = 50_000,
                    ωxprior ::NTuple{2,Float64} = (0.,10.),
                    ω1prior ::NTuple{2,Float64} = (1.,1.),
                    ω0prior ::Float64           = 1e-1,
                    σ²prior ::Float64           = 1e-1,
                    λprior  ::Float64           = 1e-1,
                    weight  ::NTuple{4,Float64} = (0.15,0.05,0.02,0.02),
                    σ²i     ::Float64           = 1.,
                    ωxi     ::Float64           = 0.,
                    ω1i     ::Float64           = 0.01,
                    ω0i     ::Float64           = 0.01,
                    λ1i     ::Float64           = 1.0,
                    λ0i     ::Float64           = 0.2,
                    fix_ωx  ::Bool              = false,
                    fix_ω1  ::Bool              = false,
                    fix_ω0  ::Bool              = false)

  parvec = [1:6...]
  np = length(parvec)

  # add to parameter update vector according to weights
  append!(parvec, fill(1,ceil(Int64,np*weight[1])))
  append!(parvec, fill(2,ceil(Int64,np*weight[2])))
  append!(parvec, repeat(3:4,  inner = ceil(Int64,np*weight[3])))
  append!(parvec, repeat(5:6,  inner = ceil(Int64,np*weight[4])))

  # if fixed ωx, ω1 or ω0, then remove
  fix_ωx && filter!(x -> x ≠ 2, parvec)
  fix_ω1 && filter!(x -> x ≠ 3, parvec)
  fix_ω0 && filter!(x -> x ≠ 4, parvec)

  # burning phase
  llc, prc, σ²c, ωxc, ω1c, ω0c, λ1c, λ0c, ptn =
    burn_tribe(σ²i, ωxi, ω1i, ω0i, λ1i, λ0i,
               λprior, ωxprior, ω1prior, ω0prior, σ²prior,
               np, parvec, nburn)

  # log for nthin
  lthin = 0

  # progress bar
  p = Progress(niter, dt=5, desc="mcmc...", barlen=20, color=:green)

  # write to file
  open(out_file*".log", "w")

  open(out_file*".log", "a") do f

    print(f, "iteration\tlikelihood\tprior\tomega_x\tomega_1\tomega_0\tsigma2\tlambda_1\tlambda_0\tcollision_probability\n")

    #=
    start MCMC
    =#

    for it = Base.OneTo(niter)

      # update vector order
      # Update vector
      shuffle!(parvec)

      for up = parvec

        # update λ1
        if up == 5

          # update λ
          λ1p = mulupt(λ1c, ptn[5])::Float64

          llr = 0.0

          prr = llrdexp_x(λ1p, λ1c, λprior)

          if -randexp()()() < (llr + prr + log(λ1p/λ1c))
            llc    += llr::Float64
            prc    += prr::Float64
            λ1c     = λ1p::Float64
          end

        # update λ0
        elseif up == 6

          # update λ
          λ0p = mulupt(λ0c, ptn[6])::Float64

          llr = 0.0

          prr = llrdexp_x(λ0p, λ0c, λprior)

          if -randexp()()() < (llr + prr + log(λ0p/λ0c))
            llc     += llr::Float64
            prc     += prr::Float64
            λ0c      = λ0p::Float64
          end

        # if σ² is updated
        elseif up == 1

          σ²p = mulupt(σ²c, ptn[1])::Float64

          #likelihood ratio
          if ωxc >= 0.0
            llr = 0.0
          else
            llr = 0.0
          end

          # prior ratio
          prr = llrdexp_x(σ²p, σ²c, σ²prior)

          if -randexp()()() < (llr + prr + log(σ²p/σ²c))
            llc += llr::Float64
            prc += prr::Float64
            σ²c  = σ²p::Float64
          end

        #update ωx
        elseif up == 2

          ωxp = addupt(ωxc, ptn[2])::Float64

          #likelihood ratio

          llr = 0.0

          # prior ratio
          prr = llrdnorm_x(ωxp, ωxc, ωxprior[1], ωxprior[2])

          if -randexp()()() < (llr + prr)
            llc += llr::Float64
            prc += prr::Float64
            ωxc  = ωxp::Float64
          end

        #update ω1
        elseif up == 3

          ω1p = addupt(ω1c, ptn[3])::Float64

          # proposal likelihood and prior
          llr = 0.0

          # prior ratio
          prr = llrdbeta_x(ω1p, ω1c, ω1prior[1], ω1prior[2])

          if -randexp()()() < (llr + prr)
            llc += llr::Float64
            prc += prr::Float64
            ω1c  = ω1p::Float64
          end

        # update ω0
        else

          ω0p = mulupt(ω0c, ptn[4])

          llr = 0.0

          # prior ratio
          prr = llrdexp_x(ω0p, ω0c, ω0prior)

          if -randexp()()() < (llr + prr + log(ω0p/ω0c))
            llc += llr
            prc += prr
            ω0c  = ω0p
          end
        end

      end # end parameter loop

      # log parameters
      lthin += 1
      if lthin == nthin
        pci = Pc(f_λ(λ1c,ω1c,1.0), f_λ(λ0c,ω0c,1.0), 0.01)
        print(f, it,"\t", llc, "\t", prc,"\t",ωxc,"\t",ω1c,"\t",ω0c,"\t",
             σ²c,"\t",λ1c,"\t",λ0c,"\t",pci,"\n")
        lthin = 0
      end

      next!(p)
    end # end MCMC
  end # end print loop

  return llc, prc, ωxc, ω1c, ω0c, σ²c, λ1c, λ0c
end

"""
    burn_compete(...)

Burning & adaptive phase for MCMC **under the prior**.
"""
function burn_tribe(σ²c     ::Float64,
                    ωxc     ::Float64,
                    ω1c     ::Float64,
                    ω0c     ::Float64,
                    λ1c     ::Float64,
                    λ0c     ::Float64,
                    λprior  ::Float64,
                    ωxprior ::NTuple{2,Float64},
                    ω1prior ::NTuple{2,Float64},
                    ω0prior ::Float64,
                    σ²prior ::Float64,
                    np      ::Int64,
                    parvec  ::Array{Int64,1},
                    nburn   ::Int64;
                    obj_ar  ::Float64 = 0.234,
                    tune_int::Int64   = 100)

  llc = 0.0

  prc = allλpr(  λ1c, λ0c, λprior)            +
        logdexp( σ²c, σ²prior)                +
        logdnorm(ωxc, ωxprior[1], ωxprior[2]) +
        logdbeta(ω1c, ω1prior[1], ω1prior[2]) +
        logdexp( ω0c, ω0prior)

  # make scaling function
  scalef = makescalef(obj_ar)

  # rest of tuning parameters
  ptn = fill(.1, np)

  # initialize acceptance log
  ltn = zeros(Int64, np)
  lup = zeros(Int64, np)
  lac = zeros(Int64, np)

  # progress bar
  p = Progress(nburn, dt=5, desc="burn...", barlen=20, color=:green)

  # print number of parameters
  print_with_color(:green,
    "\nωx updates per iter = ", lastindex(filter(x -> x == 2,parvec)),
    "\nω1 updates per iter = ", lastindex(filter(x -> x == 3,parvec)),
    "\nω0 updates per iter = ", lastindex(filter(x -> x == 4,parvec)),
    "\nσ² updates per iter = ", lastindex(filter(x -> x == 1,parvec)),
    "\nλ1 updates per iter = ", lastindex(filter(x -> x == 5,parvec)),
    "\nλ0 updates per iter = ", lastindex(filter(x -> x == 6,parvec)),
    "\nθ  updates per iter = ", length(parvec), "\n")

  #start burnin
  for it = Base.OneTo(nburn)

    # Update vector
    shuffle!(parvec)

    for up = parvec

      # update λ1
      if up == 5

        # update λ
        λ1p = mulupt(λ1c, ptn[5])::Float64

        llr = 0.0

        prr = llrdexp_x(λ1p, λ1c, λprior)

        if -randexp()()() < (llr + prr + log(λ1p/λ1c))
          llc    += llr::Float64
          prc    += prr::Float64
          λ1c     = λ1p::Float64
          lac[5] += 1
        end

      # update λ0
      elseif up == 6

        # update λ
        λ0p = mulupt(λ0c, ptn[6])::Float64

        llr = 0.0

        prr = llrdexp_x(λ0p, λ0c, λprior)

        if -randexp()()() < (llr + prr + log(λ0p/λ0c))
          llc     += llr::Float64
          prc     += prr::Float64
          λ0c      = λ0p::Float64
          lac[6] += 1
        end

      # if σ² is updated
      elseif up == 1

        σ²p = mulupt(σ²c, ptn[1])::Float64

        #likelihood ratio
        if ωxc >= 0.0
          llr = 0.0
        else
          llr = 0.0
        end

        # prior ratio
        prr = llrdexp_x(σ²p, σ²c, σ²prior)

        if -randexp()()() < (llr + prr + log(σ²p/σ²c))
          llc += llr::Float64
          prc += prr::Float64
          σ²c  = σ²p::Float64
          lac[1] += 1
        end

      #update ωx
      elseif up == 2

        ωxp = addupt(ωxc, ptn[2])::Float64

        #likelihood ratio

        llr = 0.0

        # prior ratio
        prr = llrdnorm_x(ωxp, ωxc, ωxprior[1], ωxprior[2])

        if -randexp()()() < (llr + prr)
          llc += llr::Float64
          prc += prr::Float64
          ωxc  = ωxp::Float64
          lac[2] += 1
        end

      #update ω1
      elseif up == 3

        ω1p = addupt(ω1c, ptn[3])::Float64

        # proposal likelihood and prior
        llr = 0.0

        # prior ratio
        prr = llrdbeta_x(ω1p, ω1c, ω1prior[1], ω1prior[2])

        if -randexp()()() < (llr + prr)
          llc += llr::Float64
          prc += prr::Float64
          ω1c  = ω1p::Float64
          lac[3] += 1
        end

      # update ω0
      else

        ω0p = mulupt(ω0c, ptn[4])

        llr = 0.0

        # prior ratio
        prr = llrdexp_x(ω0p, ω0c, ω0prior)

        if -randexp()()() < (llr + prr + log(ω0p/ω0c))
          llc += llr
          prc += prr
          ω0c  = ω0p
          lac[4] += 1
        end
      end

      # log number of updates
      ltn[up] += 1
      lup[up] += 1

      if (in(tune_int,ltn))
        wts = find(map(x -> x == tune_int,ltn))      # which to scale
        for j = wts
          ar     = lac[j]/lup[j]
          ptn[j] = scalef(ptn[j],ar)
          ltn[j] = 0
        end
      end
    end

    next!(p)
  end

  return llc, prc, σ²c, ωxc, ω1c, ω0c, λ1c, λ0c, ptn
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
"""
`init_subspec!(m::Model990)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model990)
    if subspec(m) == "ss2"
        return
    elseif subspec(m) == "ss3"
        ss3!(m)
    elseif subspec(m) == "ss5"
        ss5!(m)
    elseif subspec(m) == "ss6"
        # Normalize s.t. sum of variances of anticipated shocks equals the
        # variance of the contemporaneous shock
        # See measurement equation
        return
    else
        error("This subspec is not defined.")
    end
end

"""
```
ss3!(m::Model990)
```

Initialize subspec 3 of `Model990`. This is the same as subspec 2, but with
Iskander's changes (including `betabar` defined correctly with `σ_c`, not
`σ_ω_star`).
"""
function ss3!(m::Model990)
    return
end

"""
`ss5(m::Model990)`

Initializes subspecification 5 for Model990. Specifically, fixes ι_w
and ι_p to 0 (so that intermediate goods producers who do not readjust
prices and wages in a given period do not index to inflation.)
"""
function ss5!(m::Model990)

    m <= parameter(:ι_p, 0.0, fixed=true,
                   description= "ι_p: The persistence of last period's inflation in
                   the equation that describes the intertemporal
                   change in prices for intermediate goods producers
                   who cannot adjust prices. The change in prices is a
                   geometric average of steady-state inflation
                   (π_star, with weight (1-ι_p)) and last period's
                   inflation (π_{t-1})).",
                   tex_label="\\iota_p")

    m <= parameter(:ι_w,   0.0, fixed=true,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")
end


text:
# elment Types
const CPE3 = "CPE3"
const CPE4 = "CPE4"
const CPS3 = "CPS3"
const CPS4 = "CPS4"
const C3D4 = "C3D4"
const C3D6 = "C3D6"
const C3D8 = "C3D8"

function apply_mesh_periodicity_constraints(
    uc::AbstractUnitCell
)
    dim::Int = dimension(uc)
    parents::Tuple = dim == 2 ? (:XLB, :YLB,) : (:XLB, :YLB, :ZLB)
    uc_side_lengths = side_lengths(uc)
    for a_parent in parents
        #
        affine_matrix = [
        1.0, 0.0, 0.0, 0.0,
        0.0, 1.0, 0.0, 0.0,
        0.0, 0.0, 1.0, 0.0,
        0.0, 0.0, 0.0, 1.0,
        ]
        #
        if a_parent == :XLB
            affine_matrix[4] = uc_side_lengths[1]
        elseif a_parent == :YLB
            affine_matrix[8] = uc_side_lengths[2]
        elseif a_parent == :ZLB
            affine_matrix[12] = uc_side_lengths[3]
        end
        #
        if isa(uc, UDC2D)
            pxl, pyl, pxu, pyu = buffer_bbox(uc, a_parent)
            pzl, pzu = 0.0, 0.0
        else
            pxl, pyl, pzl, pxu, pyu, pzu = buffer_bbox(uc, a_parent)
        end
        #
        parent_entities = gmsh.model.get_entities_in_bounding_box(pxl, pyl, pzl, pxu, pyu, pzu, dim-1)
        for (ap_dim, ap_tag) in parent_entities
            apxl, apyl, apzl, apxu, apyu, apzu = gmsh.model.get_bounding_box(ap_dim, ap_tag)
            apxl += (affine_matrix[4, ] - ϵ)
            apyl += (affine_matrix[8, ] - ϵ)
            apzl += (affine_matrix[12,] - ϵ)
            apxu += (affine_matrix[4, ] + ϵ)
            apyu += (affine_matrix[8, ] + ϵ)
            apzu += (affine_matrix[12,] + ϵ)
            chidren_entities = gmsh.model.get_entities_in_bounding_box(
                apxl, apyl, apzl, apxu, apyu, apzu, dim-1
            )
            for (ac_dim, ac_tag) in chidren_entities
                (acxl, acyl, aczl, acxu, acyu, aczu) = gmsh.model.get_bounding_box(ac_dim, ac_tag)
                if (
                    abs(acxl-apxl) < ϵ &&
                    abs(acyl-apyl) < ϵ &&
                    abs(aczl-apzl) < ϵ &&
                    abs(acxu-apxu) < ϵ &&
                    abs(acyu-apyu) < ϵ &&
                    abs(aczu-apzu) < ϵ
                )
                    gmsh.model.mesh.set_periodic(dim, [ac_tag,], [ap_tag,], affine_matrix)
                end
            end
        end
    end
end

function check_generated_ele_types(
    ele_types::Tuple,
    dim::Int64
)
    num_tet     = convert(Int64, gmsh.option.get_number("Mesh.NbTetrahedra"))  # C3D4
    num_hex     = convert(Int64, gmsh.option.get_number("Mesh.NbHexahedra"))  # C3D8
    num_prism   = convert(Int64, gmsh.option.get_number("Mesh.NbPrisms"))  # C3D6
    num_pyramids= convert(Int64, gmsh.option.get_number("Mesh.NbPyramids"))  #C3D5
    num_trihedra= convert(Int64, gmsh.option.get_number("Mesh.NbTrihedra"))
    num_tri  = convert(Int64, gmsh.option.get_number("Mesh.NbTriangles"))  # CPS3 || CPE3
    num_quad = convert(Int64, gmsh.option.get_number("Mesh.NbQuadrangles"))  # CPS4 || CPE4
    if dim==3
        if !(:C3D4 in ele_types) && (num_tet > 0)
            @warn "C3D4 / Tet elements are generated which are not requested...!"
        end
        if !(:C3D6 in ele_types) && (num_prism > 0)
            @warn "C3D6 / Prism elements are generated which are not requested...!"
        end
        if !(:C3D8 in ele_types) && (num_hex > 0)
            @warn "C3D8 / Brick elements are generated which are not requested...!"
        end
        if num_pyramids >0
            @warn "Pyramidal elements are generated which are not requested...!"
        end
        if num_trihedra >0
            @warn "Trihedra elements are generated which are not requested...!"
        end
    else
        if !(:CPS3 in ele_types || CPE3 in ele_types) && (num_tri > 0)
            @warn "Triangular elements are generated which are not requested...!"
        end
        if !(:CPS4 in ele_types || CPE4 in ele_types) && (num_quad > 0)
            @warn "Quad elements are generated which are not requested...!"
        end
    end
end

function get_mesh_statistics(dim::Int64=-1)::Dict{String, Int64}
    mesh_stats = Dict{String, Int64}()
    if (dim == -1 ) || (dim == 3)
        mesh_stats["num_tet"]     = convert(Int64, gmsh.option.get_number("Mesh.NbTetrahedra"))
        mesh_stats["num_hex"]     = convert(Int64, gmsh.option.get_number("Mesh.NbHexahedra"))
        mesh_stats["num_prism"]   = convert(Int64, gmsh.option.get_number("Mesh.NbPrisms"))
        mesh_stats["num_pyramids"]= convert(Int64, gmsh.option.get_number("Mesh.NbPyramids"))
        mesh_stats["num_trihedra"]= convert(Int64, gmsh.option.get_number("Mesh.NbTrihedra"))
    end
    if (dim == -1 ) || (dim == 2)
        mesh_stats["num_tri"] = convert(Int64, gmsh.option.get_number("Mesh.NbTriangles"))
        mesh_stats["num_quad"] = convert(Int64, gmsh.option.get_number("Mesh.NbQuadrangles"))
    end
    #
    mesh_stats["total_nodes"] = convert(Int64, gmsh.option.get_number("Mesh.NbNodes"))
    #
    mesh_stats["num_3D_elements"] = (
        mesh_stats["num_tet"] +
        mesh_stats["num_hex"] +
        mesh_stats["num_prism"] +
        mesh_stats["num_pyramids"] +
        mesh_stats["num_trihedra"]
    )
    mesh_stats["num_2D_elements"] = (
        mesh_stats["num_tri"] +
        mesh_stats["num_quad"]
    )
    #
    return mesh_stats
end

function write_mesh_statistics(dim::Int64=-1)
    #
    mesh_stats = get_mesh_statistics(dim)
    #
    println(repeat("=", 50))
    println("\t Mesh Statistics")
    println(repeat("=", 50))
    #
    println("--- NODES ---")
    println("Number of nodes ::                 \t", mesh_stats["total_nodes"])
    println(repeat("-", 10))
    println("--- 3D ELEMENTS ---")
    println("Number of Tetrahedra elements ::   \t", mesh_stats["num_tet"])
    println("Number of Hexahedral elements ::   \t", mesh_stats["num_hex"])
    println("Number of Prsim elements ::        \t", mesh_stats["num_prism"])
    println("Number of Pyramidal elements ::    \t", mesh_stats["num_pyramids"])
    println("Number of Trihedral elements ::    \t", mesh_stats["num_trihedra"])
    println(repeat("-", 20))
    println("TOTAL ::                           \t", mesh_stats["num_3D_elements"])
    println(repeat("-", 20))
    #

    println("--- 2D ELEMENTS ---")
    println("Number of Triangles ::             \t", mesh_stats["num_tri"])
    println("Number of Quadrangles ::           \t", mesh_stats["num_quad"])
    println(repeat("-", 20))
    println("TOTAL ::                           \t", mesh_stats["num_2D_elements"])
    println(repeat("-", 20))
    println(repeat("=", 50))
end

#  , , , mesh_min_size, mesh_max_size, mesh_opt_algorithm, show_mesh_stats

function generate_mesh(
    uc::AbstractUnitCell,
    mesh_periodicity::Bool,
    element_types::Tuple{Vararg{Symbol}},
    extr_dir_num_ele::Vector{Int64},
    min_ele_size::Float64,
    max_ele_size::Float64,
    mesh_opt_algorithm::String,
    show_mesh_stats::Bool,
)
    #
    if mesh_periodicity
        apply_mesh_periodicity_constraints(uc)
    end
    #
    element_types = begin
        if (element_types[1:end] == (:DEFAULT, ) )
            if isa(uc, UDC2D)
                (CPS3,)
            elseif isa(uc, UDC3D)
                isempty(extr_dir_num_ele) ? (C3D4,) : (C3D6, C3D8,)
            elseif isa(uc, PRC)
                (C3D4,)
            end
        else
            element_types
        end
    end
    #
    # Setting constraints to get the mesh of desired element types
    if any([i in (:CPS4, :CPE4, :C3D8) for i in element_types])
        gmsh.option.set_number("Mesh.RecombineAll", 1)
    end
    if (length(element_types)==1 && (element_types[1] == :CPE4 || element_types[1] == :CPS4 || element_types[1] == :C3D8))
        # Purely quad elements in the mesh
        gmsh.option.set_number("Mesh.Algorithm", 8)
        gmsh.option.set_number("Mesh.RecombinationAlgorithm", 3)
    end
    #
    gmsh.option.set_number("Mesh.MeshSizeMin", min_ele_size)
    gmsh.option.set_number("Mesh.MeshSizeMax", max_ele_size)
    #
    gmsh.model.mesh.generate(dimension(uc))
    gmsh.model.mesh.optimize(mesh_opt_algorithm, true)
    gmsh.model.mesh.remove_duplicate_nodes()
    gmsh.model.mesh.renumber_nodes()
    #
    check_generated_ele_types(element_types, dimension(uc))
    #
    if show_mesh_stats
        write_mesh_statistics()
    end
    #
end

"""
    get_mesh_data()

Returns a dictionary of mesh data with the following key-values pairs

- `"all_node_tags"` => `Vector{Int}`
- `"all_node_coordinates"` => `Matrix{Float64}`
- `"matrix_element_connectivity"` => `Matrix{Int}`
- `"inclusions_element_connectivity"` => `Matrix{Int}`
- `"total_nodes"` => Int
- `"num_2D_elements"` => Int
- `"num_3D_elements"` => Int

"""
function get_mesh_data(
    geom_dim::Int,
    geometry_tags::Vector{Int},
)::Dict{String, Any}
    mesh_data = Dict{String, Any}()
    # -----------------------------------
    #   Nodal data collection
    # -----------------------------------
    nt, nc, _ = gmsh.model.mesh.get_nodes()
    mesh_data["all_node_tags"] = convert(Vector{Int64}, nt)  # Vector{Int64}
    mesh_data["all_node_coordinates"] = reshape(nc, 3, :)  # Matrix{Float64}
    #
    # -----------------------------------
    #   Collecting element connectivity
    # -----------------------------------
    element_connectivity::Dict{Int64, Matrix{Int64}} = Dict()
    ele_tags::Vector{Vector{Int}} = Vector{Int}[]
    ele_node_tags::Vector{Vector{Int}} = Vector{Int}[]
    ele_types::Vector{Int} = Int[]
    elt_num_nodes::Int = 0

    function _element_connectivity(gtags::Vector{Int64})::Dict{Int64, Matrix{Int64}}
        element_connectivity = Dict{Int64, Matrix{Int64}}()
        for ag_tag in gtags
            ele_types, ele_tags, ele_node_tags = gmsh.model.mesh.get_elements(geom_dim, ag_tag)
            #
            for (etk, ael_type) in enumerate(ele_types)
                _, _, _, elt_num_nodes, _, _ = gmsh.model.mesh.get_element_properties(ael_type)
                elt_ele_connectivity = [
                    reshape(ele_tags[etk], 1, :);
                    reshape(ele_node_tags[etk], elt_num_nodes, :)
                ]
                # if element type (ael_type,) connectivity data is exisiting then append new data.
                if ael_type in keys(element_connectivity)
                    elt_ele_connectivity = hcat(element_connectivity[ael_type], elt_ele_connectivity)
                end
                #
                element_connectivity[ael_type] = elt_ele_connectivity
            end
        end
        return element_connectivity
    end
    mesh_data["matrix_element_connectivity"] = _element_connectivity(geometry_tags[1:1])
    mesh_data["inclusions_element_connectivity"] = _element_connectivity(geometry_tags[2:end])
    #
    # Adding mesh Statistics like number of nodes, elements....etc.
    merge!(mesh_data, get_mesh_statistics())
    #
    return mesh_data
end


text:
module OneDimensional

using DataFrames, Gadfly

export goldenSectionSearch, fibSearch

function ternarySearch(f, left, right, ɛ)
  callCnt, ls, rs, rels = 0, [left], [right], [1.0]

  while right - left > ɛ
    a = (left * 2 + right) / 3
    b = (left + right * 2) / 3
    if f(a) < f(b)
      right = b
    else
      left = a
    end
    push!(rels, (right - left) / (rs[end] - ls[end]))
    push!(ls, left)
    push!(rs, right)
    callCnt += 2
  end
  (left + right) / 2.0, callCnt, ls, rs, rels
end

function dichotomySearch(f, left, right, ɛ)
  callCnt, ls, rs, rels = 0, [left], [right], [1.0]

  δ = ɛ / 4
  while right - left > ɛ
    a = (left + right - δ) / 2
    b = (left + right + δ) / 2
    if f(a) < f(b)
      right = b
    else
      left = a
    end
    push!(rels, (right - left) / (rs[end] - ls[end]))
    push!(ls, left)
    push!(rs, right)
    callCnt += 2
  end
  (left + right) / 2.0, callCnt, ls, rs, rels
end

function goldenSectionSearch(f, left, right, ɛ)
  callCnt, ls, rs, rels = 2, [left], [right], [1.0]

  resphi = 2 - φ
  a = left + resphi * (right - left)
  b = right - resphi * (right - left)
  fa, fb = f(a), f(b)

  while right - left > ɛ
    if fa < fb
      right = b
      b = a
      fb = fa
      a = left + resphi * (right - left)
      fa = f(a)
    else
      left = a
      a = b
      fa = fb
      b = right - resphi * (right - left)
      fb = f(b)
    end
    push!(rels, (right - left) / (rs[end] - ls[end]))
    push!(ls, left)
    push!(rs, right)
    callCnt += 1
  end
  (left + right) / 2, callCnt, ls, rs, rels
end

function fibSearch(f, left, right, ɛ)
  callCnt, ls, rs, rels = 2, [left], [right], [1.0]

  fibs = [1, 2]
  while (right - left) / fibs[end] > ɛ
    push!(fibs, fibs[end] + fibs[end - 1])
  end

  x1 = left + (right - left) * fibs[end - 2] / fibs[end]
  #x2 = left + (right - left) * fibs[end - 1] / fibs[end]
  x2 = left + right - x1
  f1, f2 = f(x1), f(x2)

  for k=1:length(fibs)-3
    if f1 < f2
      right = x2
      x2 = x1
      #x1 = left + (right - left) * fibs[end - k - 2] / fibs[end - k]
      x1 = left + right - x2
      f2 = f1
      f1 = f(x1)
    else
      left = x1
      x1 = x2
      #x2 = left + (right - left) * fibs[end - k - 1] / fibs[end - k]
      x2 = left + right - x1
      f1 = f2
      f2 = f(x2)
    end
    push!(rels, (right - left) / (rs[end] - ls[end]))
    push!(ls, left)
    push!(rs, right)
    callCnt += 1
  end
  (left + right) / 2, callCnt, ls, rs, rels
end

function uniformSearch(f, left::Float64, right::Float64, n::Int32)
  fs = [f(left + i * (right - left) / (n + 1)) for i in 0:n+1]
  minval, j, = fs[1], 0
  for i in 2:length(fs)
    if minval > fs[i]
      minval, j = fs[i], i - 1
    end
  end
  left + j * (right - left) / (n + 1), minval, n + 1
end

function uniformSearch(f, left::Float64, right::Float64, ɛ::Float64)
  n = convert(Int32, (right - left) / ɛ - 1.0)
  return uniformSearch(f, left, right, n)
end

function uniformSequenceSearch(f, left, right, n::Int64, ɛ)
  callCnt, ls, rs, rels = 0, [left], [right], [1.0]

  while right - left > ɛ
    fs = [f(left + i * (right - left) / n) for i in 0:n]
    minval, j = fs[1], 1
    for i in 2:length(fs)
      if minval > fs[i]
        minval, j = fs[i], i - 1
      end
    end

    if j == 0
      right = left + (right - left) / n
    elseif j == length(fs) - 1
      left = right - (right - left) / n
    else
      a = left + (j - 1) * (right - left) / n
      b = left + (j + 1) * (right - left) / n
      left, right = a, b
    end

    push!(rels, (right - left) / (rs[end] - ls[end]))
    push!(ls, left)
    push!(rs, right)
    callCnt += n + 1
  end
  (left + right) / 2, callCnt, ls, rs, rels
end

type Line
  a :: Float64
  b :: Float64
  c :: Float64
end

function polylineSearch(f, L, left, right, ɛ)
  lineLipshitz(L, fx0, x0) = Line(-L, 1, L * x0 - fx0)

  function segmentLine(x1, y1, x2, y2)
    a = y1 - y2
    b = x2 - x1
    c = -a * x1 - b * y1
    Line(a, b, c)
  end

  function lineIntersection(l1::Line, l2::Line)
    Δ = det([[l1.a, l1.b] [l2.a, l2.b]])
    det([[-l1.c, l1.b] [-l2.c, l2.b]]) / Δ, det([[l1.a, -l1.c] [l2.a, -l2.c]]) / Δ
  end

  function minDiff(yc, fs)
    miny, j = yc[1], 1
    for i in 2:length(yc)
      if miny > yc[i]
        miny, j = yc[i], i
      end
    end
    fs[j] - yc[j]
  end

  fl, fr = f(left), f(right)
  pnt = lineIntersection(lineLipshitz(-L, fl, left), lineLipshitz(L, fr, right))
  xc = [left, pnt[1], right]
  yc = [fl, pnt[2], fr]
  fs = [fl, f(xc[2]), fr]
  callCnt = 3

  while minDiff(yc, fs) > ɛ
    # finding minimum value through polyline
    minval, j = yc[1], 1
    for i = 2:length(yc)
      if minval > yc[i]
        minval = yc[i]
        j = i
      end
    end

    # create two tangents
    newPike = f(xc[j])
    l1, l2 = lineLipshitz(-L, newPike, xc[j]), lineLipshitz(L, newPike, xc[j])
    callCnt += 1

    # cnage polyline insertint two new segments
    pnt1 = lineIntersection(l2, segmentLine(xc[j], yc[j], xc[j - 1], yc[j - 1]))
    pnt2 = lineIntersection(l1, segmentLine(xc[j], yc[j], xc[j + 1], yc[j + 1]))

    xc = [xc[1:j - 1]; pnt1[1]; xc[j]; pnt2[1]; xc[j + 1:end]]
    yc = [yc[1:j - 1]; pnt1[2]; newPike; pnt2[2]; yc[j + 1:end]]
    fs = [fs[1:j - 1]; f(pnt1[1]); newPike; f(pnt2[1]); xc[j + 1:end]]
    callCnt += 2
  end
  xc, yc, callCnt
end

minx, callNum, ls, rs, rels = ternarySearch(x -> x^2 + 2x, -5.0, 5.0, 0.001)
#df = DataFrame(L=ls, R=rs, Rel=rels)
#writetable("ternary.txt", df)
minx, callNum, ls, rs, rels = dichotomySearch(x -> x^2 + 2x, -5.0, 5.0, 0.001)
#df = DataFrame(L=ls, R=rs, Rel=rels)
#writetable("dichotomy.txt", df)
minx, callNum, ls, rs, rels = goldenSectionSearch(x -> x^2 + 2x, -5.0, 5.0, 0.001)
#df = DataFrame(L=ls, R=rs, Rel=rels)
#writetable("goldenSection.txt", df)
minx, callNum, ls, rs, rels = fibSearch(x -> x^2 + 2x, -5.0, 5.0, 0.001)
#df = DataFrame(L=ls, R=rs, Rel=rels)
#writetable("fib.txt", df)
#minf, fval = uniformSearch(x -> x^2 + 2x, -5.0, 5.0, 0.1)
minx, callNum, ls, rs, rels = uniformSequenceSearch(x -> x^2 + 2x, -5.0, 5.0, 5, 0.001)
#df = DataFrame(L=ls, R=rs, Rel=rels)
#writetable("uniform.txt", df)

eps = [0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000000001]
ds, gs, fs, us = Int32[], Int32[], Int32[], Int32[]
f(x) = x^2 + 2x
for ɛ in eps
  _, dscnt, _, _, _ = dichotomySearch(f, -5.0, 5.0, ɛ)
  _, gscnt, _, _, _ = goldenSectionSearch(f, -5.0, 5.0, ɛ)
  _, fscnt, _, _, _ = fibSearch(f, -5.0, 5.0, ɛ)
  _, uscnt, _, _, _ = uniformSequenceSearch(f, -5.0, 5.0, 5, ɛ)
  push!(ds, dscnt)
  push!(gs, gscnt)
  push!(fs, fscnt)
  push!(us, uscnt)
end

dsit = [calls / 2.0 for calls in ds]
gsit = [calls - 2 for calls in gs]
fsit = [calls - 2 for calls in fs]
usit = [calls / 6 for calls in fs]

plot(layer(x=dsit, y=ds, Theme(default_color=color("purple")), Geom.point, Geom.line),
     layer(x=gsit, y=gs, Theme(default_color=color("orange")), Geom.point, Geom.line),
     layer(x=fsit, y=fs, Theme(default_color=color("green")), Geom.point, Geom.line),
     layer(x=usit, y=us, Theme(default_color=color("red")), Geom.point, Geom.line),
     Guide.xticks(ticks=[0:5:60]), Guide.yticks(ticks=[0:10:150]),
     Guide.title("Связь числа итераций и числа вычислений функции"),
     Guide.xlabel("Число итераций"), Guide.ylabel("Число вызовов функций"))

epsLogs = [-log(ɛ) for ɛ in eps]

plot(layer(x=epsLogs, y=ds, Theme(default_color=color("purple")), Geom.point, Geom.line),
     layer(x=epsLogs, y=gs, Theme(default_color=color("orange")), Geom.point, Geom.line),
     layer(x=epsLogs, y=fs, Theme(default_color=color("green")), Geom.point, Geom.line),
     layer(x=epsLogs, y=us, Theme(default_color=color("red")), Geom.point, Geom.line),
     Guide.title("Зависимость числа вызовов функции от минус логарифма точности"),
     Guide.xlabel("Точность"), Guide.ylabel("Число вызовов функций"))

xlines, ylines, cntCall = polylineSearch(x -> 0.5 * x^2 * sin(2x), 6, -3.0, 3.0, 0.5)

Gadfly.set_default_plot_size(12cm, 12cm)
plot(layer(x -> 0.5 * x^2 * sin(2x), -3, 3),
     layer(x=xlines, y=ylines, Theme(default_color=color("orange")), Geom.point, Geom.line),
     Guide.xticks(ticks=[-3:1:3]),
     Guide.title("Метод ломанныx функции y=1/2 ⋅ x^2 ⋅ sin(2x)"),
     Theme(default_color=color("purple")))

end

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
"""
    within_interval(x, interval)

Return the indices of the elements of `x` that lie within the given
`interval`.
"""
within_interval(x::AbstractRange, interval::Interval) = findall(in(interval), x)


text:
using SpikingNN
using Plots

# SRM0 params
η₀ = 5.0
τᵣ = 1.0
vth = 0.5

# Input spike train params
rate = 0.01
T = 15
∂t = 0.01
n = convert(Int, ceil(T / ∂t))

srm = Neuron(QueuedSynapse(Synapse.Alpha()), SRM0(η₀, τᵣ), Threshold.Ideal(vth))
input = ConstantRate(rate)
spikes = excite!(srm, input, n)

# callback to record voltages
voltages = Float64[]
record = function ()
    push!(voltages, getvoltage(srm))
end

# simulate
@time output = simulate!(srm, n; dt = ∂t, cb = record, dense = true)

# plot raster plot
raster_plot = rasterplot(∂t .* spikes, ∂t .* output, label = ["Input", "Output"], xlabel = "Time (sec)",
                         title = "Raster Plot (α response)")

# plot dense voltage recording
t = ∂t:∂t:(∂t * n)
voltage_plot = plot(t, voltages, title = "SRM Membrane Potential with Varying Presynaptic Responses",
                                 xlabel = "Time (sec)", ylabel = "Potential (V)", label = "α response")

# resimulate using presynaptic response
voltages = Float64[]
srm = Neuron(QueuedSynapse(Synapse.EPSP(ϵ₀ = 2, τm = 0.5, τs = 1)), SRM0(η₀, τᵣ), Threshold.Ideal(vth))
excite!(srm, spikes)
@time simulate!(srm, n; dt = ∂t, cb = record, dense = true)

# plot voltages with response function
voltage_plot = plot!(t, voltages, label = "EPSP response")

plot(raster_plot, voltage_plot, layout = grid(2, 1))


text:
using Documenter, GeneralizedSDistributions

makedocs(;
    modules=[GeneralizedSDistributions],
    authors="Alex Knudson <aknudson@nevada.unr.edu>",
    sitename="GeneralizedSDistributions.jl",
    doctest=false,
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/adknudson/GeneralizedSDistributions.jl",
    devbranch="main",
    versions = ["stable" => "v^", "v#.#", "dev" => "main"]
)


text:
# For exporting names in imported packages to scripts using this module
# TODO: maybe don't reexport anything (except maybe distances?)
using Reexport

@reexport using Levenshtein
@reexport using BioSequences
@reexport using Distributions
@reexport using StatsBase
@reexport using Distances
@reexport using JuMP
@reexport using Clp

include("smdfuncs.jl")


text:
import FileIO

include("model.jl")

function generate_prior_samples()
    for i=1:100
        trace = simulate(model, ())
        img = get_call_record(trace).retval
        img = max.(zeros(img), img)
        img = min.(ones(img), img)
        id = @sprintf("%03d", i)
        FileIO.save("prior/img-$id.png", colorview(Gray, img))
    end
end
println("generating prior samples..")
generate_prior_samples()

const num_train = 100 #10000

function generate_training_data()
    traces = Vector{Any}(num_train)
    for i=1:num_train
        traces[i] = get_assignment(simulate(model, ()))
        if i % 100 == 0
            println("$i of $num_train")
        end
    end
    traces
end

mutable struct ADAMOptimizerState
    m::Dict{Symbol,Any}
    v::Dict{Symbol,Any}
    t::Float64
end

function ADAMOptimizerState()
    m = Dict{Symbol,Any}()
    v = Dict{Symbol,Any}()
    for param in [:W1, :b1, :W2, :b2, :W3, :b3]
        m[param] = zero(proposal.params[param])
        v[param] = zero(proposal.params[param])
    end
    ADAMOptimizerState(m, v, 1.)
end

function update_params!(state::ADAMOptimizerState)
    beta1 = 0.9
    beta2 = 0.999
    eps = 1e-8
    alpha = 0.00001
    for param in [:W1, :b1, :W2, :b2, :W3, :b3]
        grad = -get_param_grad(proposal, param)
        m_new = beta1 .* state.m[param] + (1-beta1).*grad
        v_new = beta2 .* state.v[param] + (1-beta2).*(grad.^2)
        m_hat = m_new ./ (1. - beta1.^state.t)
        v_hat = v_new ./ (1. - beta2.^state.t)
        current = get_param(proposal, param)
        new = current - alpha * m_hat ./ (sqrt.(v_hat) + eps)
        set_param!(proposal, param, new)
        zero_param_grad!(proposal, param)
        state.m[param] = m_new
        state.v[param] = v_new
    end
    state.t += 1
end

function save_params(gen_func, filename)
    data = Dict(String(name) => param for (name, param) in gen_func.params)
    FileIO.save(filename, data)
end

function train_inference_network(all_traces, num_iter)

    # initialize parameters
    init_param!(proposal, :W1, randn(num_hidden1, num_input) * sqrt(2. / num_input))
    init_param!(proposal, :b1, zeros(num_hidden1))
    init_param!(proposal, :W2, randn(num_hidden2, num_hidden1) * sqrt(2. / num_hidden1))
    init_param!(proposal, :b2, zeros(num_hidden2))
    init_param!(proposal, :W3, randn(num_output, num_hidden2) * sqrt(2. / num_hidden2))
    init_param!(proposal, :b3, zeros(num_output))

    # initialize optimizer
    state = ADAMOptimizerState()

    # do training
    minibatch_size = 100
    tic()
    for iter=1:num_iter
        minibatch = randperm(num_train)[1:minibatch_size]
        traces = all_traces[minibatch]
        @assert length(traces) == minibatch_size
        scores = Vector{Float64}(minibatch_size)
        for (i, model_choices) in enumerate(traces)
            (proposal_trace, _) = project(proposal, (model_choices,), model_choices)
            backprop_params(proposal, proposal_trace, nothing)
            scores[i] = get_call_record(proposal_trace).score
        end
        update_params!(state)
        println("iter: $iter, score: $(mean(scores))")
        if iter % 10 == 0
            save_params(proposal, "params.jld2")
        end
    end
    toc()
end

println("generating training data...")
tic()
const traces = generate_training_data()
toc()

println("training...")
tic()
train_inference_network(traces, 100)# 10000)
toc()

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen