text:
export @class

function allocclass(name, super)
  ptr = ccall(:objc_allocateClassPair, Ptr{Nothing}, (Ptr{Nothing}, Ptr{Cchar}, Csize_t),
              super, name, 0)
  ptr == C_NULL && error("Couldn't allocate class $name")
  return Class(ptr)
end

function register(class::Class)
  ccall(:objc_registerClassPair, Nothing, (Ptr{Nothing},),
        class)
  return class
end

createclass(name, super) = allocclass(name, super) |> register

getmethod(class::Class, sel::Selector) =
  ccall(:class_getInstanceMethod, Ptr{Nothing}, (Ptr{Nothing}, Ptr{Nothing}),
        class, sel)

methodtypeenc(method::Ptr) =
  ccall(:method_getTypeEncoding, Ptr{Cchar}, (Ptr{Nothing},),
        method) |> unsafe_string

methodtypeenc(class::Class, sel::Selector) = methodtypeenc(getmethod(class, sel))

methodtype(args...) = methodtypeenc(args...) |> parseencoding

replacemethod(class::Class, sel::Selector, imp::Ptr{Nothing}, types::String) =
  ccall(:class_replaceMethod, Bool, (Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Cchar}),
        class, sel, imp, types)

function setmethod(class::Class, sel::Selector, imp::Ptr{Nothing}, types::String)
  meth = getmethod(class, sel)
  meth ≠ C_NULL && methodtype(meth) != parseencoding(types) &&
    error("New method $(name(sel)) of $class must match $(methodtype(meth))")
  replacemethod(class, sel, imp, types)
end

# Syntax

function createdef(f, args, Ts, body, ret)
  def = :(function $f($(map((x,T)->:($x::$(ctype(T))), args, Ts)...)) end)
  for (arg, T) in zip(args, Ts)
    T in (Object, Selector) && push!(def.args[2].args, :($arg = $T($arg)))
  end
  push!(def.args[2].args, ret == Nothing ? # Ugly workaround for 0.3
                            :($body; nothing) :
                            :(convert($(ctype(ret)), $body)))
  return def
end

function createmethod(__module__, class, ex)
  instance = ex.args[1] == Symbol("@-")
  ret = Base.eval(__module__, :($(ex.args[3])))
  params = ex.args[4:end-1]
  sel = @as y params map(x->x.args[2], y) join(y, ":") string(y, ":") Selector
  Ts = @>> params map(x->x.args[3]) map(x->isa(x, Symbol) ? Object : x.args[2]) map(T->Base.eval(__module__, T))
  args = @>> params map(x->x.args[3]) map(x->isa(x, Symbol) ? x : x.args[3])
  body = ex.args[end]

  Ts = [Object, Selector, Ts...]
  args = [:self, :_cmd, args...]
  typ = encodetype(ret, Ts...)

  f = "$(class)_$(name(sel))" |> Symbol
  return quote
    $(createdef(f, args, Ts, body, ret))
    cfun = @cfunction($f, $(ctype(ret)), $(Expr(:tuple, map(ctype, Ts)...)))
    ObjectiveC.setmethod($(instance ? class : :(class(class))), $sel, cfun, $typ)
  end
end

macro class(def)
  name = namify(def.args[2])
  sym = Expr(:quote, name)
  super = isexpr(def.args[2], :(<:)) ? def.args[2].args[2] : :NSObject
  expr = quote
    ObjectiveC.classexists($sym) || ObjectiveC.createclass($sym, Class($(Expr(:quote, super))))
    isdefined($__module__, $sym) || const $name = Class($sym)
  end

  for ex in def.args[3].args
    if isexpr(ex, :macrocall) && ex.args[1] in (Symbol("@+"), Symbol("@-"))
      ex = createmethod(__module__, name, ex)
    end
    push!(expr.args, ex)
  end

  :(@objc $expr; $name)
end

# quote
#   @class type Foo
#     @- (Cdouble) multiply:(Cdouble)x by:(Cdouble)y begin
#       x*y
#     end
#   end
# end |> macroexpand

# @class type Foo
#   @- (Cdouble) multiply:(Cdouble)x by:(Cdouble)y begin
#     x*y
#   end
# end

# @objc begin
#   foo = [Foo new]
#   x = [foo multiply:28 by:26]
#   [foo release]
#   x
# end


text:
function run_mc_opf_bf_dvorkin(data::Union{Dict{String,<:Any},String}, model_type::Type, solver; kwargs...)
    return _PMD.run_mc_model(data, model_type, solver, build_mc_opf_bf_dvorkin; kwargs...)
end

"constructor for branch flow opf"
function build_mc_opf_bf_dvorkin(pm::_PMD.AbstractUBFModels)
    # Variables
    _PMD.variable_mc_bus_voltage(pm)
    _PMD.variable_mc_branch_current(pm)
    _PMD.variable_mc_branch_power(pm)
    # TODO: revert to bounded in v0.10
    _PMD.variable_mc_transformer_power(pm)
    _PMD.variable_mc_gen_power_setpoint(pm)
    _PMD.variable_mc_load_setpoint(pm)
    _PMD.variable_mc_storage_power(pm)

    # Constraints
    _PMD.constraint_mc_model_current(pm)

    for i in _IM.ids(pm, :ref_buses)
        _PMD.constraint_mc_theta_ref(pm, i)
    end

    # gens should be constrained before KCL, or Pd/Qd undefined
    for id in _IM.ids(pm, :gen)
        _PMD.constraint_mc_gen_setpoint(pm, id)
    end

    # loads should be constrained before KCL, or Pd/Qd undefined
    for id in _IM.ids(pm, :load)
        _PMD.constraint_mc_load_setpoint(pm, id)
    end

    for i in _IM.ids(pm, :bus)
        _PMD.constraint_mc_load_power_balance(pm, i)
    end

    for i in _IM.ids(pm, :storage)
        _PM.constraint_storage_state(pm, i)
        _PM.constraint_storage_complementarity_nl(pm, i)
        _PMD.constraint_mc_storage_losses(pm, i)
        _PMD.constraint_mc_storage_thermal_limit(pm, i)
    end

    for i in _IM.ids(pm, :branch)
        _PMD.constraint_mc_power_losses(pm, i)
        _PMD.constraint_mc_model_voltage_magnitude_difference(pm, i)

        _PMD.constraint_mc_voltage_angle_difference(pm, i)

        _PMD.constraint_mc_thermal_limit_from(pm, i)
        _PMD.constraint_mc_thermal_limit_to(pm, i)
    end

    for i in _IM.ids(pm, :transformer)
        _PMD.constraint_mc_transformer_power(pm, i)
    end

    # Objective
    _PM.objective_min_fuel_cost(pm)
end

ext: jl


text:
module Mie

using MPI
using OffsetArrays
using Parameters
using ..Constants
using ..Data
using ..SpecialFunctions

export MieData

const ci = 1.0im

@with_kw mutable struct MieData
    number_eqns::Int64 = 0
    max_order::Int64 = 0
    number_spheres::Int64 = 0
    order::Array{Int64,1} = zeros(Int64, 0)
    offset::Array{Int64,1} = zeros(Int64, 0)
    block::Array{Int64,1} = zeros(Int64, 0)
    block_offset::Array{Int64,1} = zeros(Int64, 0)
    number_field_expansions::Array{Int64,1} = zeros(Int64, 0)
    tm_offset::Array{Int64,1} = zeros(Int64, 0)
    is_optically_active::Array{Bool,1} = fill(false, 0)
    tm_on_file::Array{Bool,1} = fill(false, 0)
    qext::Array{Float64,1} = zeros(0)
    qabs::Array{Float64,1} = zeros(0)
    an::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    cn::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    un::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    vn::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    dn::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    an_inv::Array{ComplexF64,1} = zeros(ComplexF64, 0)
    stored_tm::Array{ComplexF64,1} = zeros(ComplexF64, 0)
end

function miecoefcalc(
    run_params::MSTMParameters,
    xsp::Array{Float64,1},
    ri::OffsetArray{ComplexF64,2},
    hostsphere::Array{Int64,1},
    numberfieldexp::Array{Int64,1},
    qeps::Float64;
    tmatrix::Union{TMaxtrix,Nothing} = nothing,
)::MieData
    nsphere = length(xsp)
    @assert size(ri) == (2, nsphere + 1)
    @assert length(hostsphere) == nsphere
    @assert length(numberfieldexp) == nsphere

    mie = MieData()
    mie.max_order = 0
    mie.number_spheres = nsphere
    mie.order = zeros(Int64, nsphere)
    mie.offset = zeros(Int64, nsphere + 1)
    mie.block = zeros(Int64, nsphere)
    mie.block_offset = zeros(Int64, nsphere + 1)
    mie.number_field_expansions = numberfieldexp
    mie.tm_offset = zeros(Int64, nsphere)
    mie.is_optically_active = fill(false, nsphere)
    mie.tm_on_file = fill(false, nsphere)
    mie.qext = zeros(nsphere)
    mie.qabs = zeros(nsphere)

    # TODO: Need to consider communicators other than COMM_WORLD
    rank = MPI.Comm_rank(MPI.COMM_WORLD)

    ntermstot = 0
    nblktot = 0
    sizetmstore = 0
    rihost = zeros(ComplexF64, 2)

    # TODO: Implement usage of T-Matrix files
    if !isnothing(tmatrix)
    end

    for i in 1:nsphere
        if !mie.tm_on_file[i]
            rihost = ri[:, hostsphere[i]]
            _, _, _, nodrn, _, _, _, _, _, _ = mieoa(xsp[i], parent(ri[:, i]), 0, qeps, ri_medium = parent(rihost))
            mie.order[i] = nodrn
        end
    end

    # Mie order of host sphere should be no lower than that of the inner sphere
    for i in 1:nsphere
        j = hostsphere[i]
        while j != 0
            mie.order[j] = max(mie.order[j], mie.order[i])
            j = hostsphere[j]
        end
    end

    for i in 1:nsphere
        rihost = ri[:, hostsphere[i]]
        mie.is_optically_active[i] = ri[1, i] != ri[2, i]

        # TODO: Implement usage of T-Matrix files
        if mie.tm_on_file[i]
        else
            qext, qsca, qabs, nodrn, _, _, _, _, _, _ =
                mieoa(xsp[i], parent(ri[:, i]), mie.order[i], 0.0, ri_medium = parent(rihost))
            nterms = 4nodrn
            mie.offset[i] = ntermstot
            ntermstot += nterms
            mie.qext[i] = qext
            mie.qabs[i] = qabs
            mie.block_offset[i] = nblktot
            mie.order[i] = nodrn
            mie.block[i] = nodrn * (nodrn + 2) * 2
        end

        nblktot += mie.block[i] * numberfieldexp[i]
        mie.max_order = max(mie.max_order, mie.order[i])
    end

    mie.offset[nsphere + 1] = ntermstot
    mie.block_offset[nsphere + 1] = nblktot
    mie.number_eqns = nblktot

    mie.an = zeros(ComplexF64, ntermstot)
    mie.cn = zeros(ComplexF64, ntermstot)
    mie.un = zeros(ComplexF64, ntermstot)
    mie.vn = zeros(ComplexF64, ntermstot)
    mie.dn = zeros(ComplexF64, ntermstot)
    mie.an_inv = zeros(ComplexF64, ntermstot)

    # TODO: mieoa() is called 3 times for each sphere, can there be any improvements?
    for i in 1:nsphere
        if !mie.tm_on_file[i]
            rihost = ri[:, hostsphere[i]]
            qext, qsca, qabs, nodrn, anp, cnp, unp, vnp, dnp, anpinv =
                mieoa(xsp[i], parent(ri[:, i]), mie.order[i], 0.0, ri_medium = parent(rihost))
            nterms = 4nodrn
            n1 = mie.offset[i] + 1
            n2 = mie.offset[i] + nterms
            mie.an[n1:n2] = reshape(anp, nterms)
            mie.cn[n1:n2] = reshape(cnp, nterms)
            mie.dn[n1:n2] = reshape(dnp, nterms)
            mie.un[n1:n2] = reshape(unp, nterms)
            mie.vn[n1:n2] = reshape(vnp, nterms)
            mie.an_inv[n1:n2] = reshape(anpinv, nterms)
        end
    end

    return mie
end

"""
    mieoa(x, ri, nodr0, qeps, ri_medium)

Optically active Lorenz/Mie coefficients
"""
function mieoa(
    x::Float64,
    ri::Array{ComplexF64,1},
    nodr0::Int64,
    qeps::Float64;
    ri_medium::Union{Array{ComplexF64,1},Nothing} = nothing,
)
    @assert length(ri) == 2
    @assert isnothing(ri_medium) || length(ri_medium) == 2

    rii = zeros(ComplexF64, 2, 2)
    rii[:, 1] = isnothing(ri_medium) ? [1.0 + 0.0im, 1.0 + 0.0im] : ri_medium
    rii[:, 2] = ri
    ribulk = 2.0 ./ (1.0 ./ rii[1, :] .+ 1.0 ./ rii[2, :])
    xri = x * rii
    nstop = if iszero(qeps)
        nodr0
    elseif qeps > 0.0
        floor(Int, x + 4x^(1 / 3)) + 5
    else
        ceil(Int, -qeps)
    end

    Ψ = OffsetArray(zeros(ComplexF64, nstop + 2, 2, 2), 0:(nstop + 1), 1:2, 1:2)
    Ξ = OffsetArray(zeros(ComplexF64, nstop + 2, 2, 2), 0:(nstop + 1), 1:2, 1:2)
    qext1 = zeros(nstop)
    for i in 1:2
        for p in 1:2
            Ψ[:, p, i] = cricbessel(nstop + 1, xri[p, i])
            Ξ[:, p, i] = crichankel(nstop + 1, xri[p, i])
        end
    end

    qext = 0.0
    qsca = 0.0

    Ψ_p = zeros(ComplexF64, 2, 2)
    Ξ_p = zeros(ComplexF64, 2, 2)
    amat = zeros(ComplexF64, 2, 2)
    bmat = zeros(ComplexF64, 2, 2)
    dmat = zeros(ComplexF64, 2, 2)
    gmat = zeros(ComplexF64, 2, 2)
    umat = zeros(ComplexF64, 2, 2)
    vmat = zeros(ComplexF64, 2, 2)
    anp = zeros(ComplexF64, 2, 2, nstop)
    cnp = zeros(ComplexF64, 2, 2, nstop)
    unp = zeros(ComplexF64, 2, 2, nstop)
    vnp = zeros(ComplexF64, 2, 2, nstop)
    dnp = zeros(ComplexF64, 2, 2, nstop)
    anpinv = zeros(ComplexF64, 2, 2, nstop)

    for n in 1:nstop
        for s in 1:2
            for t in 1:2
                Ψ_p[s, t] = Ψ[n - 1, s, t] - Ψ[n, s, t] / xri[s, t] * n
                Ξ_p[s, t] = Ξ[n - 1, s, t] - Ξ[n, s, t] / xri[s, t] * n
            end
        end

        for s in 1:2
            ss = (-1)^s
            for t in 1:2
                st = (-1)^t
                gmat[s, t] =
                    (ss * ribulk[1] + st * ribulk[2]) / (rii[s, 2] * rii[t, 1]) *
                    (Ψ_p[s, 2] * Ξ[n, t, 1] - ss * st * Ψ[n, s, 2] * Ξ_p[t, 1])
                amat[s, t] =
                    (ss * ribulk[1] + st * ribulk[2]) / (rii[s, 2] * rii[t, 1]) *
                    (Ψ_p[s, 2] * Ψ[n, t, 1] - ss * st * Ψ[n, s, 2] * Ψ_p[t, 1])
                umat[s, t] = (ss * ribulk[2] + st * ribulk[2]) / (rii[s, 2] * rii[t, 2]) * ci
                bmat[s, t] =
                    (ss * ribulk[2] + st * ribulk[1]) / (rii[s, 1] * rii[t, 2]) *
                    (Ξ_p[s, 1] * Ψ[n, t, 2] - ss * st * Ξ[n, s, 1] * Ψ_p[t, 2])
                dmat[s, t] = -(ss * ribulk[1] + st * ribulk[1]) / (rii[s, 1] * rii[t, 1]) * ci
                vmat[s, t] =
                    (ss * ribulk[2] + st * ribulk[1]) / (rii[s, 1] * rii[t, 2]) *
                    (Ξ_p[s, 1] * Ξ[n, t, 2] - ss * st * Ξ[n, s, 1] * Ξ_p[t, 2])
            end
        end

        gmatinv = twobytwoinverse(gmat)
        bmatinv = twobytwoinverse(bmat)

        amat = -gmatinv * amat
        umat = -gmatinv * umat
        dmat = -bmatinv * dmat
        vmat = -bmatinv * vmat

        anp[:, :, n] = amat
        dnp[:, :, n] = dmat
        unp[:, :, n] = umat
        vnp[:, :, n] = vmat
        amatinv = twobytwoinverse(amat)
        cnp[:, :, n] = dmat * amatinv
        amatinv = twobytwoinverse(amat)
        anpinv[:, :, n] = amatinv

        fn1 = 2n + 1
        for p in 1:2
            for q in 1:2
                qsca += fn1 * abs(amat[p, q]) * abs(amat[p, q])
            end
            qext1[n] -= fn1 * real(amat[p, p])
        end
        qext += qext1[n]
    end

    nodr = nstop
    if qeps > 0.0
        qextt = qext
        qext = 0.0
        for n in 1:nstop
            qext += qext1[n]
            err = abs(1.0 - qext / qextt)
            if err < qeps || n == nstop
                nodr = n
                break
            end
        end
    end

    qsca *= 2.0 / x / x
    qext *= 2.0 / x / x
    qabs = qext - qsca

    return qext,
    qsca,
    qabs,
    nodr,
    anp[:, :, 1:nodr],
    cnp[:, :, 1:nodr],
    unp[:, :, 1:nodr],
    vnp[:, :, 1:nodr],
    dnp[:, :, 1:nodr],
    anpinv[:, :, 1:nodr]
end

function onemiecoeffmult(mie::MieData, i::Int64, cx::OffsetArray{ComplexF64,3}, mie_coefficent::Char = 'a')
    @assert length(mie.order) >= i
    nodr = mie.order[i]

    @assert size(cx) == (nodr + 2, nodr, 2)

    nterms = 4nodr
    n1 = mie.offset[i] + 1
    n2 = mie.offset[i] + nterms

    an1 = if mie_coefficent == 'a'
        reshape(mie.an[n1:n2], 2, 2, nodr)
    elseif mie_coefficent == 'c'
        reshape(mie.cn[n1:n2], 2, 2, nodr)
    elseif mie_coefficent == 'd'
        reshape(mie.dn[n1:n2], 2, 2, nodr)
    elseif mie_coefficent == 'u'
        reshape(mie.un[n1:n2], 2, 2, nodr)
    else
        reshape(mie.vn[n1:n2], 2, 2, nodr)
    end

    cy = OffsetArray(zeros(ComplexF64, nodr + 2, nodr, 2), 0:(nodr + 1), 1:nodr, 1:2)

    for n in 1:nodr
        for p in 1:2
            for i in n:-1:1
                cy[n + 1, i, p] = an1[p, 1, n] * cx[n + 1, i, 1] + an1[p, 2, n] * cx[n + 1, i, 2]
            end
            for i in 0:n
                cy[i, n, p] = an1[p, 1, n] * cx[i, n, 1] + an1[p, 2, n] * cx[i, n, 2]
            end
        end
    end

    return cy
end

"""
    multmiecoeffmult(mie, nrhs, idir, ain, rhs_list)

Generalized Mie coefficient multiplication.

If idir != 1, we do the transpose.
"""
function multmiecoeffmult(
    mie::MieData,
    nrhs::Int64,
    idir::Int64,
    ain::Array{ComplexF64,1},
    rhs_list::Union{Array{Bool,1},Nothing} = nothing,
)
    neqns = mie.number_eqns
    @assert length(ain) == nrhs * neqns
    @assert isnothing(rhs_list) || length(rhs_list) == nrhs

    if isnothing(rhs_list)
        rhs_list = fill(true, nrhs)
    end

    noffi = 0
    aout = zeros(ComplexF64, nrhs * neqns)
    for i in 1:(mie.number_spheres)
        nodri = mie.order[i]
        nblki = 2nodri * (nodri + 2)
        b11 = noffi + 1

        # TODO: Implement usage of T-Matrix files
        if mie.tm_on_file[i]
        else
            aout_t = OffsetArray(zeros(ComplexF64, nodri + 2, nodri, 2, nrhs), 0:(nodri + 1), 1:nodri, 1:2, 1:nrhs)
            b12 = b11 + nblki * nrhs - 1
            gin_t = OffsetArray(reshape(ain[b11:b12], nodri + 2, nodri, 2, nrhs), 0:(nodri + 1), 1:nodri, 1:2, 1:nrhs)
            nterms = 4nodri
            n1 = mie.offset[i] + 1
            n2 = mie.offset[i] + nterms
            an1 = reshape(mie.an[n1:n2], 2, 2, nodri)
            if mie.number_field_expansions[i] == 1
                for j in 1:nrhs
                    if rhs_list[j]
                        for n in 1:nodri
                            for p in 1:2
                                for q in 1:2
                                    coeff = idir == 1 ? an1[p, q, n] : an1[q, p, n]
                                    for k in n:-1:1
                                        aout_t[n + 1, k, p, j] += gin_t[n + 1, k, q, j] * coeff
                                    end
                                    for k in 0:n
                                        aout_t[k, n, p, j] += gin_t[k, n, q, j] * coeff
                                    end
                                end
                            end
                        end
                    end
                end
                aout[b11:b12] = reshape(aout_t, nblki * nrhs)
            else
                # TODO: It seems that number of field expansions cannot exceed 2. Need to confirm this.

                b21 = noffi + nblki * nrhs + 1
                b22 = b21 + nblki * nrhs - 1
                bin_t =
                    OffsetArray(reshape(ain[b21:b22], nodri + 2, nodri, 2, nrhs), 0:(nodri + 1), 1:nodri, 1:2, 1:nrhs)
                dn1 = reshape(mie.dn[n1:n2], 2, 2, nodri)
                un1 = reshape(mie.un[n1:n2], 2, 2, nodri)
                vn1 = reshape(mie.vn[n1:n2], 2, 2, nodri)
                fout_t = OffsetArray(zeros(ComplexF64, nodri + 2, nodri, 2, nrhs), 0:(nodri + 1), 1:nodri, 1:2, 1:nrhs)
                for j in 1:nrhs
                    if rhs_list[j]
                        for n in 1:nodri
                            for p in 1:2
                                for q in 1:2
                                    coeff1, coeff2, coeff3, coeff4 =
                                        idir == 1 ? (an1[p, q, n], un1[p, q, n], dn1[p, q, n], vn1[p, q, n]) :
                                        (an1[q, p, n], dn1[q, p, n], un1[q, p, n], vn1[q, p, n])
                                    for k in n:-1:1
                                        aout_t[n + 1, k, p, j] +=
                                            gin_t[n + 1, k, q, j] * coeff1 + bin_t[n + 1, k, q, j] * coeff2
                                        fout_t[n + 1, k, p, j] +=
                                            gin_t[n + 1, k, q, j] * coeff3 + bin_t[n + 1, k, q, j] * coeff4
                                    end
                                    for k in 0:n
                                        aout_t[k, n, p, j] += gin_t[k, n, q, j] * coeff1 + bin_t[k, n, q, j] * coeff2
                                        fout_t[k, n, p, j] += gin_t[k, n, q, j] * coeff3 + bin_t[k, n, q, j] * coeff4
                                    end
                                end
                            end
                        end
                    end
                end
                aout[b11:b12] = reshape(aout_t, nblki * nrhs)
                aout[b21:b22] = reshape(fout_t, nblki * nrhs)
            end
        end
        noffi += nblki * nrhs * mie.number_field_expansions[i]
    end

    return aout
end

end # module Mie


text:
"""
`CALCULATE_AFFINE` - Compute left-hand affine transform from two Nxd point sets

```
tform = calculate_affine(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being homogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being homogeneous coords of a point, of the
    reference point set.
* tform: affine transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_affine(moving_pts, fixed_pts)
  return moving_pts \ fixed_pts
end

"""
`CALCULATE_RIGID` - Compute right-hand rigid transform from two Nxd point sets

```
tform = calculate_rigid(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being homogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being homogeneous coords of a point, of the
    reference point set.
* tform: rigid transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_rigid(moving_pts, fixed_pts)
  moving_pts = moving_pts[:,1:end-1]
  fixed_pts = fixed_pts[:,1:end-1]
  n, dim = size(moving_pts)
  moving_bar = mean(moving_pts, 1)
  fixed_bar = mean(fixed_pts, 1)
  moving_centered = moving_pts .- moving_bar
  fixed_centered = fixed_pts .- fixed_bar
  C = fixed_centered.' * moving_centered / n
  U,s,V = svd(C)
  # Rotation R in least squares sense:
  # moving_pts - moving_bar = (fixed_pts - fixed_bar)*R
  R = (U * diagm(vcat(ones(dim-1), det(U*V.'))) * V.' ).'
  t = fixed_bar - moving_bar*R
  return [R zeros(dim); t 1]
end

"""
`CALCULATE_TRANSLATION` - Compute left-hand translation from two Nxd point sets

```
tform = calculate_translation(moving_pts, fixed_pts)
```

* moving_pts: Nxd array, each row being nonhomogeneous coords of a point, of the
    point set that is moving to fit the fixed point set.
* fixed_pts: Nxd array, each row being nonhomogeneous coords of a point, of the
    reference point set.
* tform: translation transformation that produces a mapping of the moving_pts to
    fixed_pts.

```
fixed_pts = moving_pts * tform
```
"""
function calculate_translation(moving_pts, fixed_pts)
  n, dim = size(moving_pts)
  moving_bar = mean(moving_pts, 1)
  fixed_bar = mean(fixed_pts, 1)
  t = fixed_bar - moving_bar
  return [eye(2) zeros(dim); t 1]
end

function calculate_affine(matches::Matches)
  return calculate_affine(matches.src_points, matches.dst_points)
end

function calculate_rigid(matches::Matches)
  return calculate_rigid(matches.src_points, matches.dst_points)
end

function calculate_translation(matches::Matches)
  return calculate_translation(matches.src_points, matches.dst_points)
end

function calculate_affine(mesh::Mesh)
  return calculate_affine(mesh.src_nodes, mesh.dst_nodes)
end

function calculate_rigid(mesh::Mesh)
  return calculate_rigid(mesh.src_nodes, mesh.dst_nodes)
end

function calculate_translation(mesh::Mesh)
  return calculate_translation(mesh.src_nodes, mesh.dst_nodes)
end

"""
Matrix for rotation in degrees (clockwise rotation with left-handed Cartesian)
"""
function make_rotation_matrix(theta, image_size = nothing)
  angle = deg2rad(theta)
  tform = [cos(angle) sin(angle) 0; -sin(angle) cos(angle) 0; 0 0 1]
  tform[abs(tform) .< eps] = 0
  if image_size != nothing
    	rotated_bb = snap_bb(tform_bb(sz_to_bb([image_size[1], image_size[2]]), tform))
    	rotation_offset = [rotated_bb.i, rotated_bb.j]
	translation_matrix = make_translation_matrix(-rotation_offset);
	tform = tform * translation_matrix
  end
  return tform
end

function make_translation_matrix(offset)
  return [1 0 0; 0 1 0; offset[1] offset[2] 1]
end

function make_scale_matrix(scale)
  return [scale 0 0; 0 scale 0; 0 0 1]
end


text:
struct dynamic_problem
    DDEdynProblem::DDEProblem #e.g. =DDEProblem(....)
    alg::MethodOfSteps # e.g.: alg = MethodOfSteps(Tsit5())
    maxdelay::Float64
    StateSmaplingTime::Vector
    eigN::Integer # number of eigen vectors
    zerofixpont::Bool
    SolutionSet::Vector{ODESolution}
    StateCombinations::Matrix
    #EigenVectors::Matrix
    #EigenValues::Matrix
end

function dynamic_problem(prob,alg,maxdelay;Historyresolution=20,eigN=8,zerofixpont=true,reltol=1e-3,abstol=1e-3);

    #StateSmaplingTime=-maxdelay  :maxdelay/(Historyresolution-1) :0;
    StateSmaplingTime=LinRange(-maxdelay,0.0,Historyresolution);
    #first solution for random initial functions
    #TODO: the true random thing, leads to too small stap size at the initial simulation
   # solset= [
   #     solve(
   #         remake(prob;u0=CompRand(u0), h=(p,t)->u0+CompRand(u0))
   #         ,alg,reltol=1e-3,abstol=1e-3) for k in 1:eigN];

   #StateSmaplingTimeSMALL=LinRange(-maxdelay,0.0,100);
#TODO: melyik kezdő függvény a jó?!?!?!?! - lehet sima lépcsős is?

    interpolationfunctions=[
             LinearInterpolation(StateSmaplingTime,
             [[sin(k*2pi*t/maxdelay),sin(k*2pi*t/maxdelay)] for t in StateSmaplingTime]
             ) for k in 1:eigN
             ]

    solset= [
        solve(
            remake(prob;u0=CompRand(u0),
            h=(p,t)->(1.0+1.0im) .* interpolationfunctions[k](t))
            ,alg,reltol=1e-3,abstol=1e-3) for k in 1:eigN];
    #solset= [
    #    solve(
     #       remake(prob;u0=CompRand(u0),
     #       h=(p,t)->LinearInterpolation(StateSmaplingTimeSMALL,[[sin(k*t/maxdelay),sin(k*t/maxdelay)] for t in StateSmaplingTimeSMALL])(t))
     #       ,alg,reltol=1e-3,abstol=1e-3) for k in 1:eigN];
    StateCombinations=diagm(0 => ones(typeof(CompRand(prob.u0[1])),eigN))## at the initaliatio there is no combination, a random inital values are used

    dynamic_problem(prob,alg,maxdelay,StateSmaplingTime,eigN,zerofixpont,solset,StateCombinations);

end

function iterate!(dp::dynamic_problem);
    #TODO: az S-ek alapján csinálni az új history függvényt, bár ennek nem kellene, hogy értelme legyen.
    dp_prev=deepcopy(dp);
    for k in 1:dp.eigN
        hloc(p,t)=histopryremapp(t+T,dp_prev.StateCombinations,dp_prev.SolutionSet,k)
        dp.SolutionSet[k]= solve(
            remake(prob;u0= hloc(p,0),
            h=(p,t)-> hloc(p,t))
            ,alg,reltol=1e-3,abstol=1e-3)
    end
end

function CompRand(x::Vector)
    return rand(size(x)...) .- 0.5 + 1.0im*(rand(size(x)...) .- 0.5)
end
function CompRand(x::Real)
    return rand() .- 0.5 + 1.0im*(rand() .- 0.5)
end

function compute_eig!(dp::dynamic_problem)
    #TODO: az előző végén lévő adatokból meg lehet határozni: a Si=Vi-1*dp.StateCombinations
    Si=[getvalues(dp.SolutionSet[solind],t) for t in dp.StateSmaplingTime, solind in 1:dp.eigN];
    Vi=[getvalues(dp.SolutionSet[solind],(t+dp.DDEdynProblem.tspan[end])) for t in dp.StateSmaplingTime, solind in 1:dp.eigN];#initialization of the Starting Vector

    SV=Si'*Vi;
    SS=Si'*Si;

    Snorm=[norm(Si[:,k]) for k in 1:dp.eigN ];

    Hi=SS\SV .+ 0im;

    μs,Ai=eigen(Hi);
    μs,Ai=sorteigen(μs,Ai);
    dp.StateCombinations[:] .= (Ai ./ μs./ Snorm)[:];#length(dp.StateSmaplingTime)*  ./ Snorm   ./ μs
    return  μs
end

function SVi1real(dp::dynamic_problem,idx)
    Si=[(getvalues(dp.SolutionSet[solind],t))[idx] for t in dp.StateSmaplingTime, solind in 1:dp.eigN];
    Vi=[(getvalues(dp.SolutionSet[solind],(t+dp.DDEdynProblem.tspan[end])))[idx] for t in dp.StateSmaplingTime, solind in 1:dp.eigN];#initialization of the Starting Vector
 return  Si,Vi
end

function histopryremapp(t,A::Matrix,SolutionSetV,solindex)
    #sum([SolutionSetV[i](t) * A[i,solindex] for i in 1:length(SolutionSetV)])
    sum([getvalues(SolutionSetV[i],t) * A[i,solindex]  for i in 1:length(SolutionSetV)])
end

function getvalues(sol::ODESolution,t::Real)
    if t<0
        sol.prob.h(sol.prob.p,t);
    elseif t==0
        sol.prob.u0;
    else
        sol(t)
    end
end

function getvalues(sol::ODESolution,t::Vector)
    [getvalues(sol,ti)[1] for ti in t]
end

function sorteigen(evals::Vector{T},evecs::Matrix{T}) where {T}
    p = sortperm(abs.(evals))
    evals[p[end:-1:1]], evecs[:, p[end:-1:1]]
end


text:
import Base: show

"""
    Clock{T<:Number}

    Clock{T}(time, iteration)

Keeps track of the current `time` and `iteration` number.
"""
mutable struct Clock{T<:Number}
       time :: T
  iteration :: Int
end

Base.show(io::IO, c::Clock{FT}) where FT =
    println(io, "Clock{$FT}: time = ", prettytime(c.time), ", iteration = ", c.iteration)


text:
import Random
rng = Random.MersenneTwister(1234)
include("../utils/miscFunctions.jl")

"""
function modify(array)
    array[1]=20
    5
end

a = [1,2,3,4,5]
"""

"""
mutable struct myStruct
    a
end
function prueba(a)
    a.a=5
end

b = myStruct(9)
prueba(b)

println(b)
"""

"""

function argum(p=3;varargs...)
    return varargs
end

a = argum(x=5, y=6, b=10, c=15, j=20, l=16, p=10)

println(a[:c])
"""
a = Dict()

a["x"] = 5
a["y"] = 10
a["c"] = 20
a["p"] = 50
a["h"] = 10

a = Array{UInt8}(undef,1000)
b = Array{Int}(undef,1000)

c = rand(UInt8, 1000)
d = rand(Int, 1000)

a[1] = 0x5
b[1] = 5

local a1
local b1
function probar1(iteraciones)

    for i=1:iteraciones
        a::UInt8 = 0x8
        #randomIndexSelection2(829,243,rng)
    end

end

function probar2(iteraciones)

    for i=1:iteraciones
        a::UInt8 = 8
        #randomIndexSelection(699,81,rng)
        #randomIndexSelection2(100,75)
    end

end

@time probar1(10000)
@time probar2(10000)

"""
a = randn(10000)

@noinline function probar1(iteraciones)

    for i = 1:iteraciones
        len = length(a)
        b = Array{Float64}(undef, len)
        for i=1:len
            i^253+3451/12312.23251+ 953.2235*sqrt(exp(i))
            b[i] = a[i]
        end
        #randomIndexSelection2(829,243,rng)
    end

end

@noinline function probar2(iteraciones)

    for i = 1:iteraciones
        b = similar(Array{Float64}, axes(a))
        for i in eachindex(b)
            i^253+3451/12312.23251+ 953.2235*sqrt(exp(i))
            b[i] = a[i]
        end
        #randomIndexSelection2(829,243,rng)
    end

end

@time probar1(10000)
@time probar2(10000)

"""

"""
rng = Random.MersenneTwister(2)
rng2 = Random.MersenneTwister(2)

function prueba(rng::Random.MersenneTwister)
    rand(rng)
end

function prueba1(iteraciones)
    for i=1:iteraciones
        rng3=copy(rng)
    end
end

function prueba2(iteraciones)
    for i=1:iteraciones
        a=Random.MersenneTwister(5)
    end
end

@time prueba1(10000)
@time prueba2(10000)

rng = Random.MersenneTwister(2)
rng2 = Random.MersenneTwister(2)

rng2 = copy(rng)

println(rand(rng))
rng3 = copy(rng)
"""

"""
println("===================")

for i=1:15
    println(rand(rng2))
end

ncores = 4

println()

rngs=Array{Random.MersenneTwister}(undef, 4)

for i=1:4
    rngs[i] = copy(rng)
end

randoms = Threads.@threads for i=1:ncores
    println(prueba(rngs[i]))
end

#println(randoms)

Threads.@threads for N in 1:5:20
    println("The number of this iteration is N")
end

job = Distributed.@spawn prueba(rng)

println(job)
println(fetch(job))
"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
include("abstract.jl")

"""
Analysis data of determinant quantum Monte Carlo (DQMC) simulation
"""
@with_kw mutable struct DQMCAnalysis
    acc_rate::Float64 = 0.
    prop_local::Int = 0
    acc_local::Int = 0
    acc_rate_global::Float64 = 0.
    prop_global::Int = 0
    acc_global::Int = 0
end

"""
Parameters of determinant quantum Monte Carlo (DQMC)
"""
@with_kw struct DQMCParameters
    global_moves::Bool = false
    global_rate::Int = 5
    thermalization::Int = 100 # number of thermalization sweeps
    sweeps::Int = 100 # number of sweeps (after thermalization)

    all_checks::Bool = true # e.g. check if propagation is stable/instable
    safe_mult::Int = 10

    delta_tau::Float64 = 0.1
    beta::Float64
    slices::Int = beta / delta_tau
    @assert isinteger(beta / delta_tau) string("beta/delta_tau", "
        (= number of imaginary time slices) must be an integer but is",
        beta / delta_tau, ".")

    measure_rate::Int = 10
end

"""
Determinant quantum Monte Carlo (DQMC) simulation
"""
mutable struct DQMC{M<:Model, CB<:Checkerboard, ConfType<:Any,
        Stack<:AbstractDQMCStack} <: MonteCarloFlavor
    model::M
    conf::ConfType
    s::Stack

    p::DQMCParameters
    a::DQMCAnalysis
    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}

    DQMC{M, CB, ConfType, Stack}() where {M<:Model, CB<:Checkerboard,
        ConfType<:Any, Stack<:AbstractDQMCStack} = new()
end

include("stack.jl")
include("slice_matrices.jl")

"""
    DQMC(m::M; kwargs...) where M<:Model

Create a determinant quantum Monte Carlo simulation for model `m` with
keyword parameters `kwargs`.
"""
function DQMC(m::M;
        seed::Int=-1,
        checkerboard::Bool=false,
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default,
        kwargs...
    ) where M<:Model
    # default params
    # paramskwargs = filter(kw->kw[1] in fieldnames(DQMCParameters), kwargs)
    p = DQMCParameters(; kwargs...)

    geltype = greenseltype(DQMC, m)
    heltype = hoppingeltype(DQMC, m)
    conf = rand(DQMC, m, p.slices)
    mc = DQMC{M, checkerboard ? CheckerboardTrue : CheckerboardFalse,
        typeof(conf), DQMCStack{geltype, heltype}}()
    mc.model = m
    mc.p = p
    mc.s = DQMCStack{geltype, heltype}()

    init!(
        mc, seed = seed, conf = conf,
        thermalization_measurements = thermalization_measurements,
        measurements = measurements
    )
    return mc
end

"""
    DQMC(m::M, params::Dict)
    DQMC(m::M, params::NamedTuple)

Create a determinant quantum Monte Carlo simulation for model `m` with
(keyword) parameters as specified in the dictionary/named tuple `params`.
"""
DQMC(m::Model, params::Dict{Symbol, T}) where T = DQMC(m; params...)
DQMC(m::Model, params::NamedTuple) = DQMC(m; params...)

# convenience
@inline beta(mc::DQMC) = mc.p.beta
@inline nslices(mc::DQMC) = mc.p.slices
@inline model(mc::DQMC) = mc.model
@inline conf(mc::DQMC) = mc.conf
@inline current_slice(mc::DQMC) = mc.s.current_slice

# cosmetics
import Base.summary
import Base.show
Base.summary(mc::DQMC) = "DQMC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::DQMC)
    print(io, "Determinant quantum Monte Carlo simulation\n")
    print(io, "Model: ", mc.model, "\n")
    print(io, "Beta: ", mc.p.beta, " (T ≈ $(round(1/mc.p.beta, sigdigits=3)))\n")
    N_th_meas = length(mc.thermalization_measurements)
    N_me_meas = length(mc.measurements)
    print(io, "Measurements: ", N_th_meas + N_me_meas, " ($N_th_meas + $N_me_meas)")
end
Base.show(io::IO, m::MIME"text/plain", mc::DQMC) = print(io, mc)

"""
    init!(mc::DQMC[; seed::Real=-1])

Initialize the determinant quantum Monte Carlo simulation `mc`.
If `seed !=- 1` the random generator will be initialized with `Random.seed!(seed)`.
"""
function init!(mc::DQMC;
        seed::Real = -1,
        conf = rand(DQMC,model(mc),nslices(mc)),
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default
    )
    seed == -1 || Random.seed!(seed)

    mc.conf = conf

    init_hopping_matrices(mc, mc.model)
    initialize_stack(mc)

    mc.a = DQMCAnalysis()

    mc.thermalization_measurements = thermalization_measurements
    if measurements isa Dict{Symbol, AbstractMeasurement}
        mc.measurements = measurements
    elseif measurements == :default
        mc.measurements = default_measurements(mc, mc.model)
    else
        @warn(
            "`measurements` should be of type Dict{Symbol, AbstractMeasurement}, but is " *
            "$(typeof(measurements)). No measurements have been set."
        )
        mc.measurements = Dict{Symbol, AbstractMeasurement}()
    end

    nothing
end

"""
    run!(mc::DQMC[; verbose::Bool=true, sweeps::Int, thermalization::Int])

Runs the given Monte Carlo simulation `mc`.
Progress will be printed to `stdout` if `verbose=true` (default).
"""
function run!(mc::DQMC; verbose::Bool=true, sweeps::Int=mc.p.sweeps,
        thermalization=mc.p.thermalization)

    do_th_measurements = !isempty(mc.thermalization_measurements)
    do_me_measurements = !isempty(mc.measurements)
    !do_me_measurements && @warn(
        "There are no measurements set up for this simulation!"
    )
    total_sweeps = sweeps + thermalization

    start_time = now()
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # fresh stack
    verbose && println("Preparing Green's function stack")
    initialize_stack(mc) # redundant ?!
    build_stack(mc)
    propagate(mc)

    _time = time()
    verbose && println("\n\nThermalization stage - ", thermalization)
    do_th_measurements && prepare!(mc.thermalization_measurements, mc, mc.model)
    for i in 1:total_sweeps
        verbose && (i == thermalization + 1) &&println("\n\nMeasurement stage - ", sweeps)
        for u in 1:2 * nslices(mc)
            update(mc, i)

            # For optimal performance whatever is most likely to fail should be
            # checked first.
            if current_slice(mc) == nslices(mc) && i <= thermalization && mc.s.direction == -1 &&
                    iszero(mod(i, mc.p.measure_rate)) && do_th_measurements
                measure!(mc.thermalization_measurements, mc, mc.model, i)
            end
            if (i == thermalization+1)
                do_th_measurements && finish!(mc.thermalization_measurements, mc, mc.model)
                do_me_measurements && prepare!(mc.measurements, mc, mc.model)
            end
            if current_slice(mc) == nslices(mc) && mc.s.direction == -1 && i > thermalization &&
                    iszero(mod(i, mc.p.measure_rate)) && do_me_measurements
                measure!(mc.measurements, mc, mc.model, i)
            end

        end

        if mod(i, 10) == 0
            mc.a.acc_rate = mc.a.acc_rate / (10 * 2 * nslices(mc))
            mc.a.acc_rate_global = mc.a.acc_rate_global / (10 / mc.p.global_rate)
            sweep_dur = (time() - _time)/10
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
                @printf("\t\tacc rate (local) : %.1f%%\n", mc.a.acc_rate*100)
                if mc.p.global_moves
                  @printf("\t\tacc rate (global): %.1f%%\n", mc.a.acc_rate_global*100)
                  @printf("\t\tacc rate (global, overall): %.1f%%\n",
                    mc.a.acc_global/mc.a.prop_global*100)
                end
            end

            mc.a.acc_rate = 0.0
            mc.a.acc_rate_global = 0.0
            flush(stdout)
            _time = time()
        end
    end
    do_me_measurements && finish!(mc.measurements, mc, mc.model)

    mc.a.acc_rate = mc.a.acc_local / mc.a.prop_local
    mc.a.acc_rate_global = mc.a.acc_global / mc.a.prop_global

    end_time = now()
    verbose && println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
    verbose && @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)

    nothing
end

"""
    update(mc::DQMC, i::Int)

Propagates the Green's function and performs local and global updates at
current imaginary time slice.
"""
function update(mc::DQMC, i::Int)
    propagate(mc)

    # global move
    # if mc.p.global_moves && (current_slice(mc) == mc.p.slices &&
    #        mc.s.direction == -1 && iszero(mod(i, mc.p.global_rate)))
    #     mc.a.prop_global += 1
    #     b = global_move(mc, mc.model, mc.conf) # not yet in DQMC_optional, i.e. unsupported
    #     mc.a.acc_global += b
    # end

    # local moves
    sweep_spatial(mc)

    nothing
end

"""
    sweep_spatial(mc::DQMC)

Performs a sweep of local moves along spatial dimension at current
imaginary time slice.
"""
function sweep_spatial(mc::DQMC)
    m = model(mc)
    N = nsites(m)

    @inbounds for i in 1:N
        detratio, ΔE_boson, Δ = propose_local(mc, m, i, current_slice(mc), conf(mc))
        mc.a.prop_local += 1

        if abs(imag(detratio)) > 1e-6
            println("Did you expect a sign problem? imag. detratio: ",
                abs(imag(detratio)))
            @printf "%.10e" abs(imag(detratio))
        end
        p = real(exp(- ΔE_boson) * detratio)

        # Metropolis
        if p > 1 || rand() < p
            accept_local!(mc, m, i, current_slice(mc), conf(mc), Δ, detratio,
                ΔE_boson)
            mc.a.acc_rate += 1/N
            mc.a.acc_local += 1
        end
    end
    nothing
end

"""
    greens(mc::DQMC)

Obtain the current equal-time Green's function.

Internally, `mc.s.greens` is an effective Green's function. This method transforms
this effective one to the actual Green's function by multiplying hopping matrix
exponentials from left and right.
"""
function greens(mc::DQMC_CBFalse)
    eThalfminus = mc.s.hopping_matrix_exp
    eThalfplus = mc.s.hopping_matrix_exp_inv

    greens = copy(mc.s.greens)
    greens .= greens * eThalfminus
    greens .= eThalfplus * greens
    return greens
end
function greens(mc::DQMC_CBTrue)
    chkr_hop_half_minus = mc.s.chkr_hop_half
    chkr_hop_half_plus = mc.s.chkr_hop_half_inv

    greens = copy(mc.s.greens)

    @inbounds @views begin
        for i in reverse(1:mc.s.n_groups)
          greens .= greens * chkr_hop_half_minus[i]
        end
        for i in reverse(1:mc.s.n_groups)
          greens .= chkr_hop_half_plus[i] * greens
        end
    end
    return greens
end

include("DQMC_mandatory.jl")
include("DQMC_optional.jl")
include("measurements.jl")


text:
module ConnectionPools

import Base.delete!

export ConnectionPool,
    # Getters
       get_connection!,
       get_n_connections,
       get_n_unoccupied,
       get_n_occupied,
       get_target,
       get_peak,
       get_wait,
       get_n_tries,
    # Setters
       set_target!,
       set_peak!,
       set_wait!,
       set_n_tries!,
    # Cleaning up
       release!,
       delete!

include("constructors.jl")
include("getters.jl")
include("setters.jl")
include("cleanup.jl")

function new_connection end
function disconnect end

end # module


text:
using LorentzVectors
using Test
using Random

Random.seed!(8372946187652352328)

@testset "All tests" begin
    include("basics.jl")
    include("algebra.jl")
    include("random.jl")
    include("geometry.jl")
end;


text:
valsdomain_type_promote{T<:Complex}(::Type{T},::Type{T})=T,T
valsdomain_type_promote{T<:Real}(::Type{T},::Type{T})=T,T
valsdomain_type_promote(::Type{Int},::Type{Int})=Float64,Int
valsdomain_type_promote{T<:Real,V<:Real}(::Type{T},::Type{Complex{V}})=promote_type(T,V),Complex{promote_type(T,V)}
valsdomain_type_promote{T<:Real,V<:Real}(::Type{Complex{T}},::Type{V})=Complex{promote_type(T,V)},promote_type(T,V)
valsdomain_type_promote{T<:Integer}(::Type{T},::Type{Int})=Float64,Int
valsdomain_type_promote{T<:Real}(::Type{T},::Type{Int})=T,Int
valsdomain_type_promote{T<:Complex}(::Type{T},::Type{Int})=T,Int
valsdomain_type_promote{T<:Integer,V<:Real}(::Type{T},::Type{V})=valsdomain_type_promote(Float64,V)
valsdomain_type_promote{T<:Integer,V<:Complex}(::Type{T},::Type{V})=valsdomain_type_promote(Float64,V)
valsdomain_type_promote{T<:Real}(::Type{T},::Type{Vector{T}})=T,Vector{T}
valsdomain_type_promote{T,V}(::Type{T},::Type{V})=promote_type(T,V),promote_type(T,V)

function choosefuneltype(ftype,Td)
    if !( ftype<: Number || ( (ftype <: Array) && (ftype.parameters[1] <: Number) ) )
        warn("Function outputs type $(ftype), which is not a Number")
    end

    Tprom = ftype

    if ftype <: Number #TODO should also work for array-valued functions
        Tprom,Tpromd=valsdomain_type_promote(ftype,Td)

        if ftype != Int && Tprom != ftype
                warn("Promoting function output type from $(ftype) to $(Tprom)")
        end
        if Tpromd != Td
                warn("Space domain number type $(Td) is not compatible with coefficient type $(Tprom)")
                #TODO should construct a new Space that contains a domain where the numbers have been promoted
                #and call constructor with this Space.
        end
    end

    Tprom
end

function defaultFun{ReComp}(f,d::Space{ReComp},n::Integer)
    pts=points(d, n)
    if !hasnumargs(f,1)  # Splat out Vec
        return Fun(xy->f(xy...),d,n)
    end

    f1=f(pts[1])

    if isa(f1,Array) && !isa(d,ArraySpace)
        return Fun(f,ArraySpace(d,size(f1)...),n)
    end

    # we need 3 eltype calls for the case Interval(Point([1.,1.]))
    Tprom=choosefuneltype(typeof(f1),eltype(eltype(eltype(domain(d)))))

    vals=Tprom[f(x) for x in pts]
    Fun(transform(d,vals),d)
end

Fun{ReComp}(f,d::Space{ReComp},n::Integer)=defaultFun(f,d,n)

# the following is to avoid ambiguity
# Fun(f::Fun,d) should be equivalent to Fun(x->f(x),d)
#TODO: fall back to Fun(x->f(x),d) if conversion not implemented?
Fun(f::Fun,d::Space)=Fun(coefficients(f,d),d)
Fun{T<:Space}(f::Fun,::Type{T})=Fun(f,T(domain(f)))

Fun(f::AbstractVector,T::Type)=Fun(f,T())

Fun(f,T::Type)=Fun(f,T())
Fun(f,T::Type,n::Integer)=Fun(f,T(),n)

Fun(f::AbstractVector,d::Domain)=Fun(f,Space(d))

Fun(f,d::Domain,n)=Fun(f,Space(d),n)

# We do zero special since zero exists even when one doesn't
Fun{T<:Space}(c::Number,::Type{T})=c==0?zeros(T(AnyDomain())):c*ones(T(AnyDomain()))
Fun(c::Number,d::Domain)=c==0?c*zeros(d):c*ones(d)
Fun(c::Number,d::Space)=c==0?c*zeros(eltype(d),d):c*ones(eltype(d),d)

## List constructor

Fun{T<:Domain}(c::Number,dl::AbstractVector{T})=Fun(c,UnionDomain(dl))
Fun{T<:Domain}(f,dl::AbstractVector{T})=Fun(f,UnionDomain(dl))
Fun{T<:Domain}(f,dl::AbstractVector{T},n::Integer)=Fun(f,UnionDomain(dl),n)

## Adaptive constructors

function randomFun(f,d::IntervalDomain)
    @assert d == Interval()

    #TODO: implement other domains

    Fun(chebyshevtransform(randomadaptivebary(f)),d)
end

# function veczerocfsFun(f,d::IntervalDomain)
#     #reuse function values
#
#     tol = 200*eps()
#
#     for logn = 4:20
#         cf = Fun(f, d, 2^logn + 1)
#         cfs=coefficients(cf)  ##TODO: general domain
#
#         if norm(cfs[:,end-8:end],Inf) < tol*norm(cfs[:,1:8],Inf)
#             nrm=norm(cfs,Inf)
#             return map!(g->chop!(g,10eps()*nrm),cf)
#         end
#     end
#
#     warn("Maximum number of coefficients reached")
#
#     Fun(f,d,2^21 + 1)
# end

function zerocfsFun(f, d::Space)
    #TODO: reuse function values?
    T = eltype(domain(d))
    if T <: Complex
        T = T.parameters[1] #get underlying real representation
    end
    r=checkpoints(d)
    f0=f(first(r))

    if !isa(d,ArraySpace) && isa(f0,Array)
        return zerocfsFun(f,ArraySpace(d,size(f0)...))
    end

    tol =T==Any?200eps():200eps(T)

    fr=map(f,r)
    maxabsfr=norm(fr,Inf)

    for logn = 4:20
        #cf = Fun(f, d, 2^logn + 1)
        cf = defaultFun(f, d, 2^logn)
        maxabsc=maxabs(cf.coefficients)
        if maxabsc==0 && maxabsfr==0
            return(zeros(d))
        end

        # we allow for transformed coefficients being a different size
        ##TODO: how to do scaling for unnormalized bases like Jacobi?
        if ncoefficients(cf) > 8 && maxabs(cf.coefficients[end-8:end]) < tol*maxabsc &&
                all(k->norm(cf(r[k])-fr[k],1)<1E-4,1:length(r))
            return chop!(cf,tol*maxabsc/10)
        end
    end

    warn("Maximum number of coefficients "*string(2^20+1)*" reached")

    Fun(f,d,2^21)
end

function abszerocfsFun(f,d::Space)
    #reuse function values
    T = eltype(domain(d))
    if T <: Complex
        T = T.parameters[1] #get underlying real representation
    end

    tol = 200eps(T)

    for logn = 4:20
        #cf = Fun(f, d, 2^logn + 1)
        cf = Fun(f, d, 2^logn)

        if maxabs(cf.coefficients[end-8:end]) < tol
            return chop!(cf,10eps(T))
        end
    end

    warn("Maximum number of coefficients "*string(2^20+1)*" reached")

    Fun(f,d,2^21)
end

function Fun(f, d::Space; method="zerocoefficients")
    T = eltype(domain(d))

    if f==identity
        identity_fun(d)
    elseif f==zero # zero is always defined
        zeros(T,d)
    elseif f==one
        ones(T,d)
    elseif !hasnumargs(f,1)  # Splat out Vec
        Fun(xy->f(xy...),d;method=method)
    elseif !isinf(dimension(d))
        Fun(f,d,dimension(d))  # use exactly dimension number of sample points
    elseif method == "zerocoefficients"
        zerocfsFun(f,d)
    elseif method == "abszerocoefficients"
        abszerocfsFun(f,d)
    else
        randomFun(f,d)
    end
end
Fun(f,d::Domain;opts...)=Fun(f,Space(d);opts...)

# this supports expanding a Fun to a larger or smaller domain.
# we take the union and then intersection to get at any singularities
# TODO: singularities in space(f)
Fun(f::Fun,d::Domain;opts...)=Fun(f,Space((d ∪ domain(f)) ∩ d);opts...)

## Aliases

Fun(f,n::Integer)=Fun(f,Interval(),n)
Fun{T<:Number}(f,d::AbstractVector{T},n::Integer)=Fun(f,Domain(d),n)
Fun{T<:Number,M<:Number}(cfs::AbstractVector{M},d::AbstractVector{T})=Fun(1.0*cfs,Domain(d))
Fun{T<:Number}(f,d::AbstractVector{T})=Fun(f,Domain(d))
Fun{T<:Number}(f::Number,d::AbstractVector{T})=Fun(f,Domain(d))
Fun(f::AbstractVector)=Fun(Domain(f))

function Fun(cfs::AbstractVector{Any},s::Space)
    @assert isempty(cfs)
    Fun(Float64[],s)
end

Fun(f::Fun{SequenceSpace},s::Space) = Fun(f.coefficients,s)


text:
module TestConstraintSplitInterval

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_split_basic()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    T = Float64
    bridged_mock = MOI.Bridges.Constraint.SplitInterval{T}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_$(F)_$(S)" for F in [
                "VariableIndex",
                "ScalarAffineFunction",
                "ScalarQuadraticFunction",
            ] for S in ["Interval", "EqualTo"]
        ],
    )
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_$(F)_Zeros" for F in [
                "VectorOfVariables",
                "VectorAffineFunction",
                "VectorQuadraticFunction",
            ]
        ],
    )
    return
end

function test_linear_integration_Interval()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.SplitInterval{Float64}(mock)
    config = MOI.Test.Config()
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [5.0, 5.0],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.NONBASIC],
            ],
            variable_basis_status = [MOI.BASIC, MOI.BASIC],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [2.5, 2.5],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
            ],
            variable_basis_status = [MOI.BASIC, MOI.BASIC],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
            ],
            variable_basis_status = [MOI.BASIC, MOI.BASIC],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [6.0, 6.0],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.NONBASIC],
            ],
            variable_basis_status = [MOI.BASIC, MOI.BASIC],
        ),
    )
    MOI.Test.test_linear_integration_Interval(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [0.0, 0.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
            constraint_basis_status = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
            ],
            variable_basis_status = [
                MOI.NONBASIC_AT_LOWER,
                MOI.NONBASIC_AT_LOWER,
            ],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
        ),
    )
    MOI.Test.test_linear_Interval_inactive(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.Interval{Float64},
            }(),
        ),
    )
    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    newf =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], vis), 0.0)
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, newf)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ newf

    MOI.modify(bridged_mock, ci, MOI.ScalarCoefficientChange(vis[2], 1.0))
    modified_f =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(2), vis), 0.0)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ modified_f
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) == 2.0
        bridge = MOI.Bridges.bridge(bridged_mock, ci)
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 0.0
        end
        MOI.set(bridged_mock, attr, ci, -2.0)
        @test MOI.get(bridged_mock, attr, ci) == -2.0
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == -2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 0.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
        ),
    )
    return
end

function _test_linear_FEASIBILITY_SENSE(T)
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.SplitInterval{T}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [one(T), one(T)],
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}) =>
                zeros(T, 2),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => zeros(T, 1),
        ),
    )
    MOI.Test.test_linear_FEASIBILITY_SENSE(bridged_mock, MOI.Test.Config(T))
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{T},
                MOI.EqualTo{T},
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}, 1),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}, 0),
        ),
    )
    return
end

function test_linear_FEASIBILITY_SENSE()
    _test_linear_FEASIBILITY_SENSE(Float64)
    _test_linear_FEASIBILITY_SENSE(Int)
    return
end

function test_conic_linear_VectorOfVariables()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.SplitInterval{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) =>
                [[-3, -1]],
        )
    MOI.Test.test_conic_linear_VectorOfVariables(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives, 0),
        ),
    )
    return
end

function _test_interval(
    mock,
    bridged_mock,
    set::MOI.Interval{T},
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{T}},
) where {T}
    haslb = set.lower != typemin(T)
    hasub = set.upper != typemax(T)
    @test MOI.Bridges.is_bridged(bridged_mock, ci)
    bridge = MOI.Bridges.bridge(bridged_mock, ci)
    if haslb
        @test bridge.lower !== nothing
        MOI.set(mock, MOI.ConstraintBasisStatus(), bridge.lower, MOI.BASIC)
    else
        @test bridge.lower === nothing
    end
    if hasub
        @test bridge.upper !== nothing
        MOI.set(mock, MOI.ConstraintBasisStatus(), bridge.upper, MOI.BASIC)
    else
        @test bridge.upper === nothing
    end
    @test set == MOI.get(bridged_mock, MOI.ConstraintSet(), ci)
    attr = MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{T}}()
    @test 0 == MOI.get(bridged_mock, attr)
    @test (haslb ? 1 : 0) == MOI.get(mock, attr)
    attr = MOI.NumberOfConstraints{MOI.VariableIndex,MOI.LessThan{T}}()
    @test 0 == MOI.get(bridged_mock, attr)
    @test (hasub ? 1 : 0) == MOI.get(mock, attr)
    attr = MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}()
    @test isempty(MOI.get(bridged_mock, attr))
    @test (haslb ? 1 : 0) == length(MOI.get(mock, attr))
    attr = MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}()
    @test isempty(MOI.get(bridged_mock, attr))
    @test (hasub ? 1 : 0) == length(MOI.get(mock, attr))
    if haslb || hasub
        @test one(T) == MOI.get(bridged_mock, MOI.ConstraintPrimal(), ci)
        @test MOI.BASIC ==
              MOI.get(bridged_mock, MOI.ConstraintBasisStatus(), ci)
        for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
            MOI.set(bridged_mock, attr, ci, T(2))
            @test T(2) == MOI.get(bridged_mock, attr, ci)
        end
    else
        for attr in [
            MOI.ConstraintPrimal(),
            MOI.ConstraintPrimalStart(),
            MOI.ConstraintBasisStatus(),
        ]
            err = ErrorException(
                "Cannot get `$attr` for a constraint " *
                "in the interval `[-Inf, Inf]`.",
            )
            @test_throws err MOI.get(bridged_mock, attr, ci)
        end
        @test zero(T) == MOI.get(bridged_mock, MOI.ConstraintDualStart(), ci)
    end
    @test zero(T) == MOI.get(bridged_mock, MOI.ConstraintDual(), ci)
end

function test_infinite_bounds(::Type{T} = Float64) where {T<:AbstractFloat}
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.SplitInterval{T}(mock)
    x = MOI.add_variable(bridged_mock)
    MOI.set(mock, MOI.VariablePrimal(), x, one(T))
    set = MOI.Interval(typemin(T), one(T))
    ci = MOI.add_constraint(bridged_mock, x, set)
    _test_interval(mock, bridged_mock, set, ci)
    set = MOI.Interval(zero(T), typemax(T))
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, set)
    _test_interval(mock, bridged_mock, set, ci)
    set = MOI.Interval(typemin(T), typemax(T))
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, set)
    _test_interval(mock, bridged_mock, set, ci)
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        (
            (MOI.VariableIndex, MOI.GreaterThan{T}, 0),
            (MOI.VariableIndex, MOI.LessThan{T}, 0),
        ),
    )
    ci = MOI.add_constraint(bridged_mock, x, set)
    _test_interval(mock, bridged_mock, set, ci)
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        (
            (MOI.VariableIndex, MOI.GreaterThan{T}, 0),
            (MOI.VariableIndex, MOI.LessThan{T}, 0),
        ),
    )
    set = MOI.Interval(zero(T), typemax(T))
    ci = MOI.add_constraint(bridged_mock, x, set)
    return _test_interval(mock, bridged_mock, set, ci)
end

end  # module

TestConstraintSplitInterval.runtests()


text:
using IntegralHistograms
using Test

@testset "IntegralHistograms.jl" begin
    include("integral_image.jl")
end


text:
# plotting script for  PMCMC, MCWM, DA-GP-MCMC, and ADA-GP-MCMC

# load packages
using DataFrames
using PyPlot

# fix for removing missing values
remove_missing_values(x) = reshape(collect(skipmissing(x)),3,:)

################################################################################
###   Plot data                                                              ###
################################################################################

# load data
y = Array(readtable("Ricker model/y_data_set_2.csv"))[:,1]

text_size = 25
label_size = 20

PyPlot.figure(figsize=(20,10))
ax = axes()
PyPlot.plot(y)
#PyPlot.xlabel("time", fontsize=text_size)
#PyPlot.ylabel("y", fontsize=text_size)
ax[:tick_params]("both",labelsize=label_size)

################################################################################
###   Plot resutls                                                              ###
################################################################################

# load plotting function
include(pwd()*"/Ricker model/Results/plotting.jl")

################################################################################
###   Results for PMCMC, MCWM, DA-GP-MCMC, and ADA-GP-MCMC                          ###
################################################################################

# set job name
jobname = "_dagpmcmc_lunarc"
jobname = "_adagpmcmc_lunarc"
jobname = "mcwm"
jobname = "pmcmc"

savepath = "Ricker model/Results/" # we are currently in the Results folder

# load data
Theta = remove_missing_values(Array(readtable(savepath*"Theta"*jobname*".csv")))
loklik_avec_priorv = Array(readtable(savepath*"loglik_avec_priorvec"*jobname*".csv"))
algorithm_parameters = Array(readtable(savepath*"algorithm_parameters"*jobname*".csv"))

# set values
loglik = loklik_avec_priorv[1,:]
accept_vec = loklik_avec_priorv[2,:]
prior_vec = loklik_avec_priorv[3,:]
burn_in = Int64(algorithm_parameters[1,1])
theta_true = algorithm_parameters[2:4,1]
theta_0 = algorithm_parameters[5:7,1]
prior_parameters = algorithm_parameters[8:end,:]

# calc results and plot results
analyse_results(Theta, loglik, accept_vec, prior_vec,theta_true, burn_in, prior_parameters)

################################################################################
###    Compare results for MCWM, DA-GP-MCMC, and ADA-GP-MCMC                                            ###
################################################################################

jobname_mcwm = "mcwm"
jobname_da = "_dagpmcmc_lunarc"
jobname_ada = "_adagpmcmc_lunarc"

savepath = "Ricker model/Results/"

Theta_mcwm = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_mcwm*".csv")))
loklik_avec_priorv_mcwm = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_mcwm*".csv"))
algorithm_parameters_mcwm = Array(readtable(savepath*"algorithm_parameters"*jobname_mcwm*".csv"))

# set values:
loglik_mcwm = loklik_avec_priorv_mcwm[1,:]
accept_vec_mcwm = loklik_avec_priorv_mcwm[2,:]
prior_vec_mcwm = loklik_avec_priorv_mcwm[3,:]
burn_in_mcwm = Int64(algorithm_parameters_mcwm[1,1])
theta_true = algorithm_parameters_mcwm[2:4,1]
theta_0_mcwm = algorithm_parameters_mcwm[5:7,1]
prior_parameters = algorithm_parameters_mcwm[8:end,:]

Theta_da = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_da*".csv")))
loklik_avec_priorv_da = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_da*".csv"))
algorithm_parameters_da = Array(readtable(savepath*"algorithm_parameters"*jobname_ada*".csv"))

# set values:
loglik_da = loklik_avec_priorv_da[1,:]
accept_vec_da = loklik_avec_priorv_da[2,:]
prior_vec_da = loklik_avec_priorv_da[3,:]
burn_in_da = Int64(algorithm_parameters_da[1,1])
theta_true_da = algorithm_parameters_da[2:4,1]
theta_0_da = algorithm_parameters_da[5:7,1]

Theta_ada = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_ada*".csv")))
loklik_avec_priorv_ada = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_ada*".csv"))
algorithm_parameters_ada = Array(readtable(savepath*"algorithm_parameters"*jobname_ada*".csv"))

# set values:
loglik_ada = loklik_avec_priorv_ada[1,:]
accept_vec_ada = loklik_avec_priorv_ada[2,:]
prior_vec_ada = loklik_avec_priorv_ada[3,:]
burn_in_ada = Int64(algorithm_parameters_ada[1,1])
theta_true_ada = algorithm_parameters_ada[2:4,1]
theta_0_ada = algorithm_parameters_ada[5:7,1]

# Posterior
x_c1 = prior_parameters[1,1]-0.5:0.01:prior_parameters[1,2]+0.5
x_c2 = prior_parameters[2,1]-0.5:0.01:prior_parameters[2,2]+0.5
x_c3 = prior_parameters[3,1]-0.5:0.01:prior_parameters[3,2]+0.5

priordens_c1 = pdf.(Uniform(prior_parameters[1,1], prior_parameters[1,2]), x_c1)
priordens_c2 = pdf.(Uniform(prior_parameters[2,1], prior_parameters[2,2]), x_c2)
priordens_c3 = pdf.(Uniform(prior_parameters[3,1], prior_parameters[3,2]), x_c3)

h1_mcwm = kde(Theta_mcwm[1,burn_in_mcwm:end])
h2_mcwm = kde(Theta_mcwm[2,burn_in_mcwm:end])
h3_mcwm = kde(Theta_mcwm[3,burn_in_mcwm:end])

h1_da = kde(Theta_da[1,burn_in_da:end])
h2_da = kde(Theta_da[2,burn_in_da:end])
h3_da = kde(Theta_da[3,burn_in_da:end])

h1_ada = kde(Theta_ada[1,burn_in_ada:end])
h2_ada = kde(Theta_ada[2,burn_in_ada:end])
h3_ada = kde(Theta_ada[3,burn_in_ada:end])

text_size = 15
label_size = 15

PyPlot.figure()
PyPlot.plot(h1_mcwm.x,h1_mcwm.density, "b")
PyPlot.hold(true)
PyPlot.plot(h1_da.x,h1_da.density, "r")
PyPlot.plot(h1_ada.x,h1_ada.density, "r--")
PyPlot.plot(x_c1,priordens_c1, "g")
PyPlot.plot((theta_true[1], theta_true[1]), (0, maximum([maximum(h1_mcwm.density);maximum(h1_da.density); maximum(h1_ada.density)])), "k")
PyPlot.xlabel(L"log $r$",fontsize=text_size)
PyPlot.ylabel(L"Density",fontsize=text_size)
PyPlot.figure()
PyPlot.plot(h2_mcwm.x,h2_mcwm.density, "b")
PyPlot.hold(true)
PyPlot.plot(h2_da.x,h2_da.density, "r")
PyPlot.plot(h2_ada.x,h2_ada.density, "r--")
PyPlot.plot(x_c2,priordens_c2, "g")
PyPlot.plot((theta_true[2], theta_true[2]), (0, maximum([maximum(h2_mcwm.density);maximum(h2_da.density); maximum(h2_ada.density)])), "k")
PyPlot.xlabel(L"log $\phi$",fontsize=text_size)
PyPlot.ylabel(L"Density",fontsize=text_size)
PyPlot.figure()
PyPlot.plot(h3_mcwm.x,h3_mcwm.density, "b")
PyPlot.hold(true)
PyPlot.plot(h3_da.x,h3_da.density, "r")
PyPlot.plot(h3_ada.x,h3_ada.density, "r--")
PyPlot.plot(x_c3,priordens_c3, "g")
PyPlot.plot((theta_true[3], theta_true[3]), (0, maximum([maximum(h3_mcwm.density);maximum(h3_da.density); maximum(h3_ada.density)])), "k")
PyPlot.xlabel(L"log $\sigma$",fontsize=text_size)
PyPlot.ylabel(L"Density",fontsize=text_size)

################################################################################
###    Compare results for PMCMC, MCWM, DA-GP-MCMC, and ADA-GP-MCMC                                            ###
################################################################################

jobname_pmcmc = "pmcmc"
jobname_mcwm = "mcwm"
jobname_da = "_dagpmcmc_lunarc"
jobname_ada = "_adagpmcmc_lunarc"

savepath = "Ricker model/Results/"

Theta_pmcmc = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_pmcmc*".csv")))
loklik_avec_priorv_pmcmc = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_pmcmc*".csv"))
algorithm_parameters_pmcmc = Array(readtable(savepath*"algorithm_parameters"*jobname_pmcmc*".csv"))

# set values:
loglik_pmcmc = loklik_avec_priorv_pmcmc[1,:]
accept_vec_pmcmc = loklik_avec_priorv_pmcmc[2,:]
prior_vec_pmcmc = loklik_avec_priorv_pmcmc[3,:]
burn_in_pmcmc = Int64(algorithm_parameters_pmcmc[1,1])
theta_true_pmcmc = algorithm_parameters_pmcmc[2:4,1]
theta_0_pmcmc = algorithm_parameters_pmcmc[5:7,1]
prior_parameters_pmcmc = algorithm_parameters_pmcmc[8:end,:]

Theta_mcwm = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_mcwm*".csv")))
loklik_avec_priorv_mcwm = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_mcwm*".csv"))
algorithm_parameters_mcwm = Array(readtable(savepath*"algorithm_parameters"*jobname_mcwm*".csv"))

# set values:
loglik_mcwm = loklik_avec_priorv_mcwm[1,:]
accept_vec_mcwm = loklik_avec_priorv_mcwm[2,:]
prior_vec_mcwm = loklik_avec_priorv_mcwm[3,:]
burn_in_mcwm = Int64(algorithm_parameters_mcwm[1,1])
theta_true = algorithm_parameters_mcwm[2:4,1]
theta_0_mcwm = algorithm_parameters_mcwm[5:7,1]
prior_parameters = algorithm_parameters_mcwm[8:end,:]

Theta_da = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_da*".csv")))
loklik_avec_priorv_da = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_da*".csv"))
algorithm_parameters_da = Array(readtable(savepath*"algorithm_parameters"*jobname_ada*".csv"))

# set values:
loglik_da = loklik_avec_priorv_da[1,:]
accept_vec_da = loklik_avec_priorv_da[2,:]
prior_vec_da = loklik_avec_priorv_da[3,:]
burn_in_da = Int64(algorithm_parameters_da[1,1])
theta_true_da = algorithm_parameters_da[2:4,1]
theta_0_da = algorithm_parameters_da[5:7,1]

Theta_ada = remove_missing_values(Array(readtable(savepath*"Theta"*jobname_ada*".csv")))
loklik_avec_priorv_ada = Array(readtable(savepath*"loglik_avec_priorvec"*jobname_ada*".csv"))
algorithm_parameters_ada = Array(readtable(savepath*"algorithm_parameters"*jobname_ada*".csv"))

# set values:
loglik_ada = loklik_avec_priorv_ada[1,:]
accept_vec_ada = loklik_avec_priorv_ada[2,:]
prior_vec_ada = loklik_avec_priorv_ada[3,:]
burn_in_ada = Int64(algorithm_parameters_ada[1,1])
theta_true_ada = algorithm_parameters_ada[2:4,1]
theta_0_ada = algorithm_parameters_ada[5:7,1]

# kernel density for posteriors
h1_pmcmc = kde(Theta_pmcmc[1,burn_in_pmcmc:end])
h2_pmcmc = kde(Theta_pmcmc[2,burn_in_pmcmc:end])
h3_pmcmc = kde(Theta_pmcmc[3,burn_in_pmcmc:end])

h1_mcwm = kde(Theta_mcwm[1,burn_in_mcwm:end])
h2_mcwm = kde(Theta_mcwm[2,burn_in_mcwm:end])
h3_mcwm = kde(Theta_mcwm[3,burn_in_mcwm:end])

h1_da = kde(Theta_da[1,burn_in_da:end])
h2_da = kde(Theta_da[2,burn_in_da:end])
h3_da = kde(Theta_da[3,burn_in_da:end])

h1_ada = kde(Theta_ada[1,burn_in_ada:end])
h2_ada = kde(Theta_ada[2,burn_in_ada:end])
h3_ada = kde(Theta_ada[3,burn_in_ada:end])

# prior
x_c1 = prior_parameters[1,1]-0.5:0.01:prior_parameters[1,2]+0.5
x_c2 = prior_parameters[2,1]-0.5:0.01:prior_parameters[2,2]+0.5
x_c3 = prior_parameters[3,1]-0.5:0.01:prior_parameters[3,2]+0.5

priordens_c1 = pdf.(Uniform(prior_parameters[1,1], prior_parameters[1,2]), x_c1)
priordens_c2 = pdf.(Uniform(prior_parameters[2,1], prior_parameters[2,2]), x_c2)
priordens_c3 = pdf.(Uniform(prior_parameters[3,1], prior_parameters[3,2]), x_c3)

text_size = 25
label_size = 20

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plot(h1_pmcmc.x,h1_pmcmc.density, "b")
PyPlot.plot(h1_mcwm.x,h1_mcwm.density, "b--")
PyPlot.hold(true)
PyPlot.plot(h1_da.x,h1_da.density, "r")
PyPlot.plot(h1_ada.x,h1_ada.density, "r--")
PyPlot.plot(x_c1,priordens_c1, "g")
PyPlot.plot((theta_true[1], theta_true[1]), (0, maximum([maximum(h1_pmcmc.density); maximum(h1_mcwm.density);maximum(h1_da.density); maximum(h1_ada.density)]) ), "k")
#PyPlot.xlabel(L"log $r$",fontsize=text_size)
#PyPlot.ylabel(L"Density",fontsize=text_size)
PyPlot.xlim((minimum(h1_pmcmc.x)-0.2, maximum(h1_pmcmc.x)+0.2))
ax[:tick_params]("both",labelsize=label_size)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plot(h2_pmcmc.x,h2_pmcmc.density, "b")
PyPlot.plot(h2_mcwm.x,h2_mcwm.density, "b--")
PyPlot.hold(true)
PyPlot.plot(h2_da.x,h2_da.density, "r")
PyPlot.plot(h2_ada.x,h2_ada.density, "r--")
PyPlot.plot(x_c2,priordens_c2, "g")
PyPlot.plot((theta_true[2], theta_true[2]), (0, maximum([maximum(h2_pmcmc.density); maximum(h2_mcwm.density);maximum(h2_da.density); maximum(h2_ada.density)]) ), "k")
#PyPlot.xlabel(L"log $\phi$",fontsize=text_size)
#PyPlot.ylabel(L"Density",fontsize=text_size)
PyPlot.xlim((minimum(h2_pmcmc.x)-0.07, maximum(h2_pmcmc.x)+0.07))
ax[:tick_params]("both",labelsize=label_size)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plot(h3_pmcmc.x,h3_pmcmc.density, "b")
PyPlot.plot(h3_mcwm.x,h3_mcwm.density, "b--")
PyPlot.hold(true)
PyPlot.plot(h3_da.x,h3_da.density, "r")
PyPlot.plot(h3_ada.x,h3_ada.density, "r--")
PyPlot.plot(x_c3,priordens_c3, "g")
PyPlot.plot((theta_true[3], theta_true[3]), (0, maximum([maximum(h3_pmcmc.density); maximum(h3_mcwm.density);maximum(h3_da.density); maximum(h3_ada.density)]) ), "k")
#PyPlot.xlabel(L"log $\sigma$",fontsize=text_size)
#PyPlot.ylabel(L"Density",fontsize=text_size)
PyPlot.xlim((minimum(h3_pmcmc.x)-0.25, maximum(h3_pmcmc.x)+0.25))
ax[:tick_params]("both",labelsize=label_size)

max_stars_repo_path: Ricker model/Results/plotting_script.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: Ricker model/Results/plotting_script.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: Ricker model/Results/plotting_script.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text: INVALID_PRELIMINARYSELECTION = "selected preliminary selection not found"


text:
using JuMP
const MOIT = MOI.Test

MOIU.@model(NoFreeVariable,
            (), (MOI.EqualTo, MOI.LessThan, MOI.GreaterThan), (MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros, MOI.RotatedSecondOrderCone, MOI.PositiveSemidefiniteConeTriangle), (),
            (), (MOI.ScalarAffineFunction,), (MOI.VectorOfVariables,), (MOI.VectorAffineFunction,))
# No free variables to make sure variable bridges are used to increase coverage
MOI.supports_constraint(::NoFreeVariable, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals}) = false

function bridged_mock(mock_optimize!::Function...;
                      model = MOI.Utilities.Model{Float64}())
    mock = MOI.Utilities.MockOptimizer(model)
    bridged = MOI.Bridges.full_bridge_optimizer(mock, Float64)
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!...)
    return bridged
end

function cached_mock(
    args...; cache = () -> MOIU.UniversalFallback(MOIU.Model{Float64}()), kws...)
    # We want the MOI backend of the JuMP model to be in `MOIU.EMPTY_OPTIMIZER`
    # mode so that it's copied at `optimize!` to test that `copy_to` works.
    # If we just return `cached`, it will be emptied in `_model` and the state
    # will be `MOIU.ATTACHED_OPTIMIZER` which is not what we want. For this
    # reason we return a `JuMP.OptimizerFactory` which returns `cached` instead.
    with_optimizer(() -> begin
        cached = MOIU.CachingOptimizer(cache(), MOIU.AUTOMATIC)
        optimizer = bridged_mock(args...; kws...)
        MOIU.reset_optimizer(cached, optimizer)
        return cached
    end)
end

function mocks(args...; kws...)
    return [bridged_mock(args...; kws...), cached_mock(args...; kws...)]
end


text:
export visstatic, visdynamic, show_strain_stress, prepare_strain_stress_data1D, prepare_strain_stress_data2D,
prepare_sequence_strain_stress_data2D, VisualizeStress2D,visσ, VisualizeStrainStressSurface, visstate
import PyPlot:scatter3D
function visstatic(domain::Domain, vmin=nothing, vmax=nothing; scaling = 1.0)
    u,v = domain.state[1:domain.nnodes], domain.state[domain.nnodes+1:end]
    nodes = domain.nodes
    fig,ax = subplots()
    temp = nodes + [u v]
    x1, x2 = minimum(temp[:,1]), maximum(temp[:,1])
    y1, y2 = minimum(temp[:,2]), maximum(temp[:,2])
    σ =[]
    for e in domain.elements
        σs = e.stress
        if σs[1]==undef
            σs = ones(length(e.stress))
        end
        push!(σ,mean([sqrt(σ[1]^2-σ[1]*σ[2]+σ[2]^2+3*σ[3]^2) for σ in σs]))
    end
    vmin = vmin==nothing ? minimum(σ) : vmin
    vmax = vmax==nothing ? maximum(σ) : vmax
    #@show vmin, vmax
    cNorm  = colors.Normalize(
            vmin=vmin,
            vmax=vmax)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
    for (k,e) in enumerate(domain.elements)
        n_ = nodes[getNodes(e),:] + scaling*[u[getNodes(e),:] v[getNodes(e),:]]
        p = plt.Polygon(n_, facecolor = scalarMap.to_rgba(σ[k]), fill=true, alpha=0.5)
        ax.add_patch(p)
    end
    xlim(x1 .-0.1,x2 .+0.1)
    ylim(y1 .-0.1,y2 .+0.1)
end

function visdynamic(domain::Domain, name::String)

    # Set up formatting for the movie files
    Writer = animation.writers.avail["html"]
    writer = Writer(fps=15, bitrate=1800)

    close("all")
    fig = figure()
    # visualization
    scat0 = scatter(domain.nodes[:,1], domain.nodes[:,2], color="grey")
    grid(true)
    ims = Any[(scat0,)]

    N = size(domain.nodes,1)
    for k = 1:length(domain.history["state"])
        u1 = domain.history["state"][k][1:N] + domain.nodes[:,1]
        u2 = domain.history["state"][k][N+1:end] + domain.nodes[:,2]

        scat = scatter(u1, u2, color="orange")
        grid(true)
        tt = gca().text(.5, 1.05,"$k")
        # s2 = scatter(nodes[div(n+1,2)*n,1], nodes[div(n+1,2)*n,2], marker="x", color="red")
        # s3 = scatter(u1[div(n+1,2)*n], u2[div(n+1,2)*n], marker="*", color="red")
        push!(ims, (scat0,scat,tt))
    end

    im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000,
                                    blit=true)
    im_ani.save("$name.html", writer=writer)

end

function show_strain_stress(domain::Domain)
    strain = domain.history["strain"]
    stress = domain.history["stress"]
    strain = hcat(strain...)
    stress = hcat(stress...)
    close("all")
    for i = 1:size(strain,1)
        plot(strain[i,:], stress[i,:])
    end
    xlabel("strain")
    ylabel("stress")
end

function prepare_strain_stress_data1D(domain::Domain)
    strain = domain.history["strain"]
    stress = domain.history["stress"]
    strain = hcat(strain...)
    stress = hcat(stress...)
    ngp = size(strain,1)
    nt = size(strain,2)
    X = zeros(nt*ngp,3)
    y = zeros(nt*ngp)
    strain = [zeros(ngp, 1) strain]
    stress = [zeros(ngp,1) stress]
    nt += 1
    k = 1
    for i = 1:ngp
        for j = 2:nt
            X[k,:] = [strain[i,j] strain[i,j-1] stress[i,j-1]] #ε, ε0, σ0
            y[k] = stress[i,j]
            k = k + 1
        end
    end
    X,y
end

function prepare_strain_stress_data1D(strain::Array{Any, 1}, stress::Array{Any, 1},)
    strain = hcat(strain...)
    stress = hcat(stress...)
    ngp = size(strain,1)
    nt = size(strain,2)
    X = zeros(nt*ngp,3)
    y = zeros(nt*ngp)
    strain = [zeros(ngp, 1) strain]
    stress = [zeros(ngp,1) stress]
    nt += 1
    k = 1
    for i = 1:ngp
        for j = 2:nt
            X[k,:] = [strain[i,j] strain[i,j-1] stress[i,j-1]] #ε, ε0, σ0
            y[k] = stress[i,j]
            k = k + 1
        end
    end
    X,y
end

function prepare_strain_stress_data2D(domain::Domain, scale::Float64=1.0)
    strain = domain.history["strain"]
    stress = domain.history["stress"]
    ngp = size(strain[1],1)
    nt = length(strain)
    X = zeros(nt*ngp,9)
    y = zeros(nt*ngp,3)
    pushfirst!(strain, zero(strain[1]))
    pushfirst!(stress, zero(stress[1]))

    k = 1
    nt += 1
    for j = 2:nt
        for i = 1:ngp
            X[k,:] = [strain[j][i,:]; strain[j-1][i,:]; stress[j-1][i,:]/scale]#ε, ε0, σ0
            y[k,:] = stress[j][i,:]/scale
            k = k + 1
        end
    end
    X,y
end

function prepare_strain_stress_data2D(strain::Array{Any, 1}, stress::Array{Any, 1}, scale::Float64=1.0)
    ngp = size(strain[1],1)
    nt = length(strain)
    X = zeros(nt*ngp,9)
    y = zeros(nt*ngp,3)
    pushfirst!(strain, zero(strain[1]))
    pushfirst!(stress, zero(stress[1]))

    k = 1
    nt += 1
    for j = 2:nt
        for i = 1:ngp
            X[k,:] = [strain[j][i,:]; strain[j-1][i,:]; stress[j-1][i,:]/scale]#ε, ε0, σ0
            y[k,:] = stress[j][i,:]/scale
            k = k + 1
        end
    end
    X,y
end

# Return strain_seq[ngp, NT+1, 3], stress_seq[ngp, NT+1, 3]
function prepare_sequence_strain_stress_data2D(domain::Domain, scale::Float64=1.0)
    strain = domain.history["strain"]
    stress = domain.history["stress"]
    ngp = size(strain[1],1)
    nt = length(strain)
    strain_seq = zeros(nt+1, ngp, 3)
    stress_seq = zeros(nt+1, ngp, 3)

    for j = 2:nt+1
            #@show nt, ngp, size(strain[j-1])
            strain_seq[j,:,:] = strain[j-1]
            stress_seq[j,:,:] = stress[j-1]
    end
    strain_seq, stress_seq
end

# https://risa.com/risahelp/risa3d/Content/3D_2D_Only_Topics/Plates%20-%20Results.htm
function postprocess_stress(stress::Array{Float64}, name::String)
    if name == "vonMises"
        σx = stress[1]; σy = stress[2]; τxy = stress[3]
        σ1 = (σx+σy)/2+sqrt((σx-σy)^2/4 + τxy^2)
        σ2 = (σx+σy)/2-sqrt((σx-σy)^2/4 + τxy^2)
        sqrt(σ1^2-σ1*σ2+σ2^2)

        #σ[1]^2-σ[1]*σ[2]+σ[2]^2+3*σ[3]^2
    elseif name == "principal"
        σx = stress[1]; σy = stress[2]; τxy = stress[3]
        σ1 = (σx+σy)/2+sqrt((σx-σy)^2/4 + τxy^2)
        σ2 = (σx+σy)/2-sqrt((σx-σy)^2/4 + τxy^2)
        [σ1, σ2]
    else
        @error("postprocess_stress does not recognize ", name)
    end
end

function VisualizeStress2D(domain::Domain)
    strain = domain.history["strain"]
    stress = domain.history["stress"]
    NT = length(stress)
    ngp = size(stress[1], 1)
    V = zeros(NT, ngp, 2)
    for i = 1:NT
        for j = 1:ngp
            V[i,j,:] = postprocess_stress(stress[i][j,:], "principal")
        end
    end
    close("all")
    for i = rand(1:ngp, 10)
        x = V[:,i,1][:]/1e8
        y = V[:,i,2][:]/1e8
        plot(x, y, ".-")
    end
    V
end

function VisualizeStress2D(σ_ref::Array{Float64}, σ_comp::Array{Float64}, NT::Int64, Nlines::Int64=1)
    ngp = Int64(size(σ_ref,1)/NT)
    V_ref = zeros(NT, ngp, 3)
    V_comp = zeros(NT, ngp, 3)
    # 2
    for i = 1:NT
        for j = 1:ngp
            V_ref[i,j,1] = postprocess_stress(σ_ref[(i-1)*ngp + j,:], "vonMises")
            V_ref[i,j,2:3] = postprocess_stress(σ_ref[(i-1)*ngp + j,:], "principal")
            V_comp[i,j,1] = postprocess_stress(σ_comp[(i-1)*ngp + j,:], "vonMises")
            V_comp[i,j,2:3] = postprocess_stress(σ_comp[(i-1)*ngp + j,:], "principal")
        end
    end
    close("all")
    col = "r"
    col2 = "g"

    k = 0
    for i = rand(1:ngp, Nlines)
        k += 1
        x = V_ref[:,i,1][:]
        y = V_ref[:,i,2][:]
        plot(x, y, ".-"*col[1])

        x = V_comp[:,i,1][:]
        y = V_comp[:,i,2][:]
        plot(x, y, ".--"*col2[1])
    end
end

function VisualizeStrainStressSurface(X::Array{Float64}, Y::Array{Float64}, seed::Int64=233)
    n = size(X,1)
    Random.seed!(seed)
    idx = rand(1:n, 2000)
    scatter3D(X[idx,1], X[idx,2], Y[idx,1], marker=".")
end

function visσ(domain::Domain, vmin=nothing, vmax=nothing; scaling = 1.0)
    u,v = domain.state[1:domain.nnodes], domain.state[domain.nnodes+1:end]
    nodes = domain.nodes
    fig,ax = subplots()
    temp = nodes + [u v]
    x1, x2 = minimum(temp[:,1]), maximum(temp[:,1])
    y1, y2 = minimum(temp[:,2]), maximum(temp[:,2])
    σ =[]
    for e in domain.elements
        σs = e.stress
        # σs = [ones(3) for i = 1:length(e.stress)] # ! remove me
        push!(σ,mean([postprocess_stress(s, "vonMises")[1] for s in σs]))
    end
    vmin = vmin==nothing ? minimum(σ) : vmin
    vmax = vmax==nothing ? maximum(σ) : vmax
    #@show vmin, vmax
    cNorm  = colors.Normalize(
            vmin=vmin,
            vmax=vmax)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
    for (k,e) in enumerate(domain.elements)
        N = getNodes(e)
        if length(N)==9
            N = N[[1;5;2;6;3;7;4;8]]
        end
        n_ = nodes[N,:] + scaling*[u[N,:] v[N,:]]
        p = plt.Polygon(n_, facecolor = scalarMap.to_rgba(σ[k]), fill=true)
        ax.add_patch(p)
    end
    scalarMap.set_array(σ)
    colorbar(scalarMap)
    xlim(x1 .-0.1,x2 .+0.1)
    ylim(y1 .-0.1,y2 .+0.1)
end

function visstate(domain::Domain, state::Array{Float64}; kwargs...)
    u,v = state[1:domain.nnodes], state[domain.nnodes+1:end]
    nodes = domain.nodes
    fig,ax = subplots()
    temp = nodes + [u v]
    x1, x2 = minimum(temp[:,1]), maximum(temp[:,1])
    y1, y2 = minimum(temp[:,2]), maximum(temp[:,2])
    for (k,e) in enumerate(domain.elements)
        N = getNodes(e)
        if length(N)==9
            N = N[[1;5;2;6;3;7;4;8]]
        end
        n_ = nodes[N,:] + [u[N,:] v[N,:]]
        p = plt.Polygon(n_; kwargs...)
        ax.add_patch(p)
    end
    xlim(x1 .-0.1,x2 .+0.1)
    ylim(y1 .-0.1,y2 .+0.1)
end

function visstate(domain::Domain; kwargs...)
    visstate(domain, domain.state; kwargs...)
end

# need 1)empty domain for data structure
#      2)data in stresses and disps
function visσ(domain::Domain, nx::Int64, ny::Int64,  stress::Array{Float64, 2}, disps::Array{Float64, 1},
    vmin=nothing, vmax=nothing; scaling = [1.0, 1.0])

    L_scale, sigma_scale = scaling

    u,v = disps[1:domain.nnodes], disps[domain.nnodes+1:end]
    nodes = domain.nodes
    ngpt = length(domain.elements[1].weights)
    ngp = Int64(sqrt(ngpt))

    ### Compute X[nx*ngp+1, ny*ngp+1], Y[nx*ngp+1, ny*ngp+1], C[nx*ngp, ny*ngp] for pcolormesh
    X = zeros(Float64, nx*ngp+1, ny*ngp+1)
    Y = zeros(Float64, nx*ngp+1, ny*ngp+1)
    C = zeros(Float64, nx*ngp, ny*ngp)

    Gp_order_1D = zeros(Int64, ngp)
    if ngp == 2
        Gp_order_1D = [1, 2]
    elseif ngp == 3
        Gp_order_1D = [2, 1, 3]
    elseif ngp == 4
        Gp_order_1D = [4, 2, 1, 3]
    else
        error(" ngp == ", ngp, " has not implemented")
    end

    Gp_order = zeros(Int64, ngp, ngp)
    for gy = 1:ngp
        for gx = 1:ngp
        Gp_order[gx, gy] = (Gp_order_1D[gx] - 1)*ngp  +  Gp_order_1D[gy]
        end
    end

    for iy = 1:ny
        for ix = 1:nx
            eid = ix + (iy-1)*nx
            ele = domain.elements[eid]
            # (iy-1)*ngp + ngp
            # (iy-1)*ngp + 1
            # (iy-1)*ngp + 0       (ix-1)*ngp + 0, (ix-1)*ngp + 1, (ix-1)*ngp + ngp

            N = getNodes(ele)
            N = N[[1;2;3;4]] #take these nodes on the four corners
            xy = L_scale*(nodes[N,:] + [u[N,:] v[N,:]])

            for gy = 0:ngp
                for gx = 0:ngp
                    w = [(1.0 - gx/ngp) * (1.0 - gy/ngp), gx/ngp * (1.0 - gy/ngp),
                          gx/ngp * gy/ngp, (1.0 - gx/ngp) * gy/ngp]

                    X[1 + (ix-1)*ngp + gx, 1 + (iy-1)*ngp + gy] = w[1]*xy[1,1] +w[2]*xy[2,1] +w[3]*xy[3,1] +w[4]*xy[4,1]
                    Y[1 + (ix-1)*ngp + gx, 1 + (iy-1)*ngp + gy] = w[1]*xy[1,2] +w[2]*xy[2,2] +w[3]*xy[3,2] +w[4]*xy[4,2]
                end
            end

            for gy = 1:ngp
                for gx = 1:ngp
                    σ = stress[(eid-1)*ngpt + Gp_order[gx ,gy], :]
                    σvm = sigma_scale * postprocess_stress(σ, "vonMises")
                    C[(ix-1)*ngp + gx, (iy-1)*ngp + gy] = σvm

                end
            end

        end
    end

    vmin = vmin==nothing ? minimum(C) : vmin
    vmax = vmax==nothing ? maximum(C) : vmax

    @show minimum(C), maximum(C) , "!!!!!!"

    fig,ax = subplots()
    x1, x2 = minimum(X), maximum(X)
    y1, y2 = minimum(Y), maximum(Y)

    p = plt.pcolormesh(X, Y, C, cmap="rainbow")

    clim(vmin,vmax)
    # cNorm  = colors.Normalize(
    #         vmin=0,
    #         vmax=1.0)
    # scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)

    # #scalarMap.set_array(σ)
    # colorbar(scalarMap)
    colorbar()
    dx, dy = (x2 - x1)*0.1, (y2 - y1)*0.1
    xlim(x1 .- dx,x2 .+ dx)
    ylim(y1 .- dy,y2 .+ dy)
    axis("equal")

    return vmin, vmax
end


text:
using Plots
using Random
using Pkg; Pkg.add("LaTeXStrings")
using LaTeXStrings

# Set working directory
cd("/Users/elvis/Documents/DSI/courses/content/stochastic-approximations/julia-code/")

# Set seed
Random.seed!(123)

# Generate MC points
n = 4000
points_x = 3 .* rand(n)
points_y = -3 .+ 7 .* rand(n)

# Define functions
f1(x) = 3*abs(cos(x)) .+ 2*sin(x)
f2(x) = -3*abs(cos(x)) .+ 2*sin(x)

# Label points in circle
labels = [f2(x) <= y <= f1(x) for (x,y) in zip(points_x, points_y)]

# Create curves
x_vals = 0:.01:3
f1_x = 3*abs.(cos.(x_vals)) .+ 2*sin.(x_vals)
f2_x = -3*abs.(cos.(x_vals)) .+ 2*sin.(x_vals)

# Create plot of MC points
plot(points_x, points_y, seriestype = :scatter,
     group = .!labels,
     legend = false,
     ms = 3,
     bg = RGB(247/255, 236/255, 226/255))

# Add curves to plot
plot!(x_vals, f1_x, seriestype= :line, label = L"f_1(x)", lc = :blue,
    legend = false)
plot!(x_vals, f2_x, seriestype= :line, label = L"f_2(x)", lc = :violet)
plot!(x_vals, f1_x, fillrange = [f2_x, f1_x], fillcolor = :blue, alpha = .12)

# Area estimate
area = 21*sum(labels)/n

savefig("../images/area-ex1-monte-carlo.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using Oceananigans, Printf

#####
##### Model set-up
#####

Nx, Ny, Nz = 16, 16, 16       # No. of grid points in x, y, and z, respectively.
Lx, Ly, Lz = 100, 100, 100    # Length of the domain in x, y, and z, respectively (m).
tf = 5000                     # Length of the simulation (s)

model = Model(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
              closure=ConstantIsotropicDiffusivity())

# Add a cube-shaped warm temperature anomaly that takes up the middle 50%
# of the domain volume.
i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

#####
##### Set up output
#####

write_grid(model)

outputs =      Dict("u" => model.velocities.u,
                    "v" => model.velocities.v,
                    "w" => model.velocities.w,
                    "T" => model.tracers.T,
                    "S" => model.tracers.S)

outputattrib = Dict("u" => ["longname" => "Velocity in the x-direction", "units" => "m/s"],
                    "v" => ["longname" => "Velocity in the y-direction", "units" => "m/s"],
                    "w" => ["longname" => "Velocity in the z-direction", "units" => "m/s"],
                    "T" => ["longname" => "Temperature", "units" => "K"],
                    "S" => ["longname" => "Salinity", "units" => "g/kg"])

globalattrib = Dict("f" => 1e-4, "name" => "Thermal bubble expt 1")

# The following writer saves a yz plane at xC[5] for all fields that
# have xC as their dimension and at xF[6] for all fields that have xF
# as their dimension. Ranges also can be specified (e.g. xC=2:10)
subsetwriter = NetCDFOutputWriter(model, outputs;
                                  interval=10, filename="dump_subset.nc",
                                  output_attributes=outputattrib,
                                  global_attributes=globalattrib,
                                  xC=5, xF=6)
push!(model.output_writers, subsetwriter)

# The following writer saves a data from the entire domain
globalwriter = NetCDFOutputWriter(model, outputs, interval=10,
                                  filename="dump_global.nc")
push!(model.output_writers, globalwriter)

#####
##### Run the simulation
#####

function terse_message(model, walltime, Δt)
    cfl = Δt / Oceananigans.cell_advection_timescale(model)
    return @sprintf("i: %d, t: %.4f s, Δt: %.1f s, cfl: %.3f, wall time: %s\n",
                    model.clock.iteration, model.clock.time, Δt, cfl, prettytime(walltime))
end

# A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.2, Δt=10.0, max_change=1.1, max_Δt=50.0)

# Run the model
while model.clock.time < tf
    update_Δt!(wizard, model)
    walltime = @elapsed time_step!(model, 10, wizard.Δt)
    @printf "%s" terse_message(model, walltime, wizard.Δt)
end

# Close the NetCDFOutputWriter
close(subsetwriter)
close(globalwriter)


text:
#=
Some functions from *Universal resilience patterns in complex networks*
by Gao et al. (2016)
=#

using Statistics: mean, std

"""
    s(N::AbstractUnipartiteNetwork; dims::Union{Nothing,Integer}=nothing)

Computes the average weighted degree. This is proportional to the (weighted)
density of interactions.

If dims is provided, the incoming (`dims=1`) or outgoing (`dims=2`) is computed.
"""
function s(N::AbstractUnipartiteNetwork; dims::Union{Nothing,Integer}=nothing)
    dims == 1 && return sum(N.A, dims=2)
    dims == 2 && return sum(N.A, dims=1)'
    if dims == nothing
        return sum(N.A) / size(N)[1]
    end
end
#=
"""
    s(N::AbstractUnipartiteNetwork; dims)

Computes the vector of incoming weighted degrees of an unipartite network.
"""
s_in(N::AbstractUnipartiteNetwork) = sum(N.A, dims=2)

"""
    s_out(N::AbstractUnipartiteNetwork)

Computes the vector of outgoing weighted degrees of an unipartite network.
"""
s_out(N::AbstractUnipartiteNetwork) = sum(N.A, dims=1)'

=#

"""
    σ_in(N::AbstractUnipartiteNetwork)

Computes the standard deviation of the ingoing weighted degree of an unipartite
network.
"""
σ_in(N::AbstractUnipartiteNetwork) = std(s(N, dims=1), corrected=false)

"""
    σ_out(N::AbstractUnipartiteNetwork)

Computes the standard deviation of the outgoing weighted degree of an unipartite
network.
"""
σ_out(N::AbstractUnipartiteNetwork) = std(s(N, dims=2), corrected=false)

"""
    symmetry(N::AbstractUnipartiteNetwork)

Computes the symmetry between s^in and s^out (the in- and outgoing weighted
degree of an unipartite network). This is computed as the Pearson correlation
between the s^in and s^out. It is hence a value between -1 and 1, where high
positive values indicate that species with many outgoing degrees tend to have
many ingoing degrees and negative values mean the opposite. An undirected network
is perfectly symmetric but, for example, a food web where predators are less likely
to be prey would have a negative symmetry.

> Goa, J., Barzael, B. and Barabási 2016. Universal resilience patterns in complex networks.
> Nature 530(7590), 307-312. doi:10.1038/nature16948

"""
symmetry(N::AbstractUnipartiteNetwork) = (mean(s(N, dims=1) .* s(N, dims=2)) - mean(s(N, dims=1)) * mean(s(N, dims=2))) / (σ_in(N) * σ_out(N))

"""
    heterogeneity(N::AbstractUnipartiteNetwork)

Computes the heterogeneity for an unipartite network, a topological characteristic
which quantifies the difference in in- and outgoing degrees between species. It
is computed as σ_in * σ_out / s_mean. A value of 0 indicates that all species
have the same (weighted) in- and outdegrees.

> Goa, J., Barzael, B. and Barabási 2016. Universal resilience patterns in complex networks.
> Nature 530(7590), 307-312. doi:10.1038/nature16948

"""
heterogeneity(N::AbstractUnipartiteNetwork) = (σ_in(N) * σ_out(N)) / s(N)

"""
    resilience(N::AbstractUnipartiteNetwork)

A resilience parameters described by Gao et al. (2016). It is a global parameters
describing the dynamics of an unipartite network as an effective 1D equation of
the form

f(xeff) = F(xeff) + βeff G(xeff, xeff)

i.e. describing a second-order term representing the effect of the network on the
dynamics of the 'effective state' xeff of the system.

> Goa, J., Barzael, B. and Barabási 2016. Universal resilience patterns in complex networks.
> Nature 530(7590), 307-312. doi:10.1038/nature16948

"""
resilience(N::AbstractUnipartiteNetwork) = dot(s(N, dims=1), s(N, dims=2)) / sum(N)