text:
using LinearAlgebra, MixedModels, StableRNGs, Test, SparseArrays

include("modelcache.jl")

@testset "femat" begin
    trm = MixedModels.FeMat(hcat(ones(30), repeat(0:9, outer = 3)), ["(Intercept)", "U"])
    piv = trm.piv
    ipiv = invperm(piv)
    @test size(trm) == (30, 2)
    @test length(trm) == 60
    @test size(trm') == (2, 30)
    @test eltype(trm) == Float64
    @test trm.x === trm.wtx
    prd = trm'trm
    @test typeof(prd) == Matrix{Float64}
    @test prd == [30.0 135.0; 135.0 855.0][piv, piv]
    wts = rand(StableRNG(123454321), 30)
    MixedModels.reweight!(trm, wts)
    @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)

    # empty fixed effects
    trm = MixedModels.FeMat(ones(10,0), String[])
    @test size(trm) == (10, 0)
    @test length(trm) == 0
    @test size(trm') == (0, 10)
    @test eltype(trm) == Float64
    @test trm.rank == 0
end

@testset "fematSparse" begin

    @testset "sparse and dense yield same fit" begin
        # deepcopy because we're going to modify
        m = deepcopy(last(models(:insteval)))
        # this is kinda sparse:
        # julia> mean(first(m.feterms).x)
        # 0.10040140325753434

        fe = first(m.feterms)
        X =  MixedModels.FeMat(SparseMatrixCSC(fe.x), fe.cnames)
        @test typeof(X.x) <: SparseMatrixCSC
        @test typeof(X.wtx) <: SparseMatrixCSC
        @test X.rank == 28
        @test X.cnames == fe.cnames

        dense_θ = copy(m.θ)
        dense_feval = m.optsum.feval
        m.feterms[1] = X
        refit!(m)

        @test dense_θ ≈ m.θ
    end

    @testset "rank defiency in sparse FeMat" begin
        trm = MixedModels.FeMat(SparseMatrixCSC(hcat(ones(30),
                                                     repeat(0:9, outer = 3),
                                                     2repeat(0:9, outer = 3))),
                                ["(Intercept)", "U", "V"])
        piv = trm.piv
        ipiv = invperm(piv)
        @test_broken rank(trm) == 2
        @test size(trm) == (30, 3)
        @test length(trm) == 90
        @test size(trm') == (3, 30)
        @test eltype(trm) == Float64
        @test trm.x === trm.wtx
        prd = trm'trm
        @test typeof(prd) == typeof(prd)
        @test prd == [30.0 135.0 270.0; 135.0 855.0 1710.0; 270.0 1710.0 3420][piv, piv]
        wts = rand(StableRNG(123454321), 30)
        MixedModels.reweight!(trm, wts)
        @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
    end

end


text:
using ModelingToolkit, OrdinaryDiffEq, Test

@variables t
sts = @variables x1(t) x2(t) x3(t) x4(t)
params = @parameters u1(t) u2(t) u3(t) u4(t)
D = Differential(t)
eqs = [x1 + x2 + u1 ~ 0
       x1 + x2 + x3 + u2 ~ 0
       x1 + D(x3) + x4 + u3 ~ 0
       2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + u4 ~ 0]
@named sys = ODESystem(eqs, t)

let dd = dummy_derivative(sys)
    has_dx1 = has_dx2 = false
    for eq in equations(dd)
        vars = ModelingToolkit.vars(eq)
        has_dx1 |= D(x1) in vars || D(D(x1)) in vars
        has_dx2 |= D(x2) in vars || D(D(x2)) in vars
    end
    @test has_dx1 ⊻ has_dx2 # only one of x1 and x2 can be a dummy derivative
    @test length(states(dd)) == length(equations(dd)) == 9
    @test length(states(structural_simplify(dd))) < 9
end

let pss = partial_state_selection(sys)
    @test length(equations(pss)) == 1
    @test length(states(pss)) == 2
    @test length(equations(ode_order_lowering(pss))) == 2
end

@parameters t σ ρ β
@variables x(t) y(t) z(t) a(t) u(t) F(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x)
       D(y) ~ x * (ρ - z) - y + β
       0 ~ z - x + y
       0 ~ a + z
       u ~ z + a]

lorenz1 = ODESystem(eqs, t, name = :lorenz1)
let al1 = alias_elimination(lorenz1)
    let lss = partial_state_selection(al1)
        @test length(equations(lss)) == 2
    end
end

# 1516
let
    @variables t
    D = Differential(t)

    @connector function Fluid_port(; name, p = 101325.0, m = 0.0, T = 293.15)
        sts = @variables p(t)=p m(t)=m [connect = Flow] T(t)=T [connect = Stream]
        ODESystem(Equation[], t, sts, []; name = name)
    end

    #this one is for latter
    @connector function Heat_port(; name, Q = 0.0, T = 293.15)
        sts = @variables T(t)=T Q(t)=Q [connect = Flow]
        ODESystem(Equation[], t, sts, []; name = name)
    end

    # like ground but for fluid systems (fluid_port.m is expected to be zero in closed loop)
    function Compensator(; name, p = 101325.0, T_back = 273.15)
        @named fluid_port = Fluid_port()
        ps = @parameters p=p T_back=T_back
        eqs = [fluid_port.p ~ p
               fluid_port.T ~ T_back]
        compose(ODESystem(eqs, t, [], ps; name = name), fluid_port)
    end

    function Source(; name, delta_p = 100, T_feed = 293.15)
        @named supply_port = Fluid_port() # expected to feed connected pipe -> m<0
        @named return_port = Fluid_port() # expected to receive from connected pipe -> m>0
        ps = @parameters delta_p=delta_p T_feed=T_feed
        eqs = [supply_port.m ~ -return_port.m
               supply_port.p ~ return_port.p + delta_p
               supply_port.T ~ instream(supply_port.T)
               return_port.T ~ T_feed]
        compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port])
    end

    function Substation(; name, T_return = 343.15)
        @named supply_port = Fluid_port() # expected to receive from connected pipe -> m>0
        @named return_port = Fluid_port() # expected to feed connected pipe -> m<0
        ps = @parameters T_return = T_return
        eqs = [supply_port.m ~ -return_port.m
               supply_port.p ~ return_port.p # zero pressure loss for now
               supply_port.T ~ instream(supply_port.T)
               return_port.T ~ T_return]
        compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port])
    end

    function Pipe(; name, L = 1000, d = 0.1, N = 100, rho = 1000, f = 1)
        @named fluid_port_a = Fluid_port()
        @named fluid_port_b = Fluid_port()
        ps = @parameters L=L d=d rho=rho f=f N=N
        sts = @variables v(t)=0.0 dp_z(t)=0.0
        eqs = [fluid_port_a.m ~ -fluid_port_b.m
               fluid_port_a.T ~ instream(fluid_port_a.T)
               fluid_port_b.T ~ fluid_port_a.T
               v * pi * d^2 / 4 * rho ~ fluid_port_a.m
               dp_z ~ abs(v) * v * 0.5 * rho * L / d * f  # pressure loss
               D(v) * rho * L ~ (fluid_port_a.p - fluid_port_b.p - dp_z)]
        compose(ODESystem(eqs, t, sts, ps; name = name), [fluid_port_a, fluid_port_b])
    end
    function System(; name, L = 10.0)
        @named compensator = Compensator()
        @named source = Source()
        @named substation = Substation()
        @named supply_pipe = Pipe(L = L)
        @named return_pipe = Pipe(L = L)
        subs = [compensator, source, substation, supply_pipe, return_pipe]
        ps = @parameters L = L
        eqs = [connect(compensator.fluid_port, source.supply_port)
               connect(source.supply_port, supply_pipe.fluid_port_a)
               connect(supply_pipe.fluid_port_b, substation.supply_port)
               connect(substation.return_port, return_pipe.fluid_port_b)
               connect(return_pipe.fluid_port_a, source.return_port)]
        compose(ODESystem(eqs, t, [], ps; name = name), subs)
    end

    @named system = System(L = 10)
    @unpack supply_pipe = system
    sys = structural_simplify(system)
    u0 = [system.supply_pipe.v => 0.1, system.return_pipe.v => 0.1, D(supply_pipe.v) => 0.0]
    # This is actually an implicit DAE system
    @test_throws Any ODEProblem(sys, u0, (0.0, 10.0), [])
    @test_throws Any ODAEProblem(sys, u0, (0.0, 10.0), [])
    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, u0, (0.0, 10.0), [])
    @test solve(prob, DFBDF()).retcode == :Success
end

# 1537
let
    @variables t
    @variables begin
        p_1(t)
        p_2(t)
        rho_1(t)
        rho_2(t)
        rho_3(t)
        u_1(t)
        u_2(t)
        u_3(t)
        mo_1(t)
        mo_2(t)
        mo_3(t)
        Ek_1(t)
        Ek_2(t)
        Ek_3(t)
    end

    @parameters dx=100 f=0.3 pipe_D=0.4

    D = Differential(t)

    eqs = [p_1 ~ 1.2e5
           p_2 ~ 1e5
           u_1 ~ 10
           mo_1 ~ u_1 * rho_1
           mo_2 ~ u_2 * rho_2
           mo_3 ~ u_3 * rho_3
           Ek_1 ~ rho_1 * u_1 * u_1
           Ek_2 ~ rho_2 * u_2 * u_2
           Ek_3 ~ rho_3 * u_3 * u_3
           rho_1 ~ p_1 / 273.11 / 300
           rho_2 ~ (p_1 + p_2) * 0.5 / 273.11 / 300
           rho_3 ~ p_2 / 273.11 / 300
           D(rho_2) ~ (mo_1 - mo_3) / dx
           D(mo_2) ~ (Ek_1 - Ek_3 + p_1 - p_2) / dx - f / 2 / pipe_D * u_2 * u_2]

    @named trans = ODESystem(eqs, t)

    sys = structural_simplify(trans)

    n = 3
    u = 0 * ones(n)
    rho = 1.2 * ones(n)

    u0 = [p_1 => 1.2e5
          p_2 => 1e5
          u_1 => 0
          u_2 => 0.1
          u_3 => 0.2
          rho_1 => 1.1
          rho_2 => 1.2
          rho_3 => 1.3
          mo_1 => 0
          mo_2 => 1
          mo_3 => 2]
    prob = ODAEProblem(sys, u0, (0.0, 0.1))
    @test solve(prob, FBDF()).retcode == :Success
end

ext: jl


text:
module Alt

using ProximalBase: shrink
using CoordinateDescent: CDOptions
using SparseArrays

function differencePrecisionNaive(Σx, Σy, λ, Ups, options = CDOptions())
  maxIter = options.maxIter
  optTol = options.optTol
  p = size(Σx, 1)

  Δ = zeros(p, p)
  A = zeros(p, p)

  iter = 1;
  while iter < maxIter
    fDone = true
    for a=1:p
      for b=a:p
        if a==b
          # diagonal elements
          x0 = Σx[a,a]*Σy[a,a] / 2.
          x1 = A[a,a] - Σy[a,a] + Σx[a,a]
          tmp = shrink(-x1/x0/2. + Δ[a,b], λ*Ups[a,b] / 2. / x0)
        else
          # off-diagonal elements
          x0 = (Σx[a,a]*Σy[b,b] + Σx[b,b]*Σy[a,a])/2. + Σx[a,b]*Σy[a,b]
          x1 = A[a,b] + A[b,a] - 2. * (Σy[a,b] - Σx[a,b])
          tmp = shrink(-x1/(2. * x0) + Δ[a,b], λ*Ups[a,b] / x0)
        end

        h = tmp - Δ[a,b]
        Δ[a,b] = tmp
        Δ[b,a] = tmp
        if abs(h) > optTol
          fDone = false
        end
        for j=1:p
          for k=1:p
            if a == b
              A[j,k] = A[j,k] + h * Σx[j,a]*Σy[a,k]
            else
              A[j,k] = A[j,k] + h * (Σx[j,a]*Σy[b,k] + Σx[j,b]*Σy[a,k])
            end
          end
        end
      end
    end

    iter = iter + 1;
    if fDone
      break
    end
  end
  Δ
end

function updateDelta!(
  Δ::SparseMatrixCSC{T},
  A::StridedMatrix{T}, Σx::StridedMatrix{T}, Σy::StridedMatrix{T},
  λ::T, Ups::StridedMatrix{T},
  options::CDOptions = CDOptions()) where {T<:AbstractFloat}
  # Delta is a sparse matrix stored in CSC format
  # only diagonal and lower triangular elements are used

  # extract fields from sparse matrix Delta
  rows = rowvals(Δ)
  vals = nonzeros(Δ)

  maxInnerIter = options.maxIter
  maxChangeTol = options.optTol

  p = size(Δ, 1)

  for iter=1:maxInnerIter
    fDone = true
    for colInd=1:p
      for j in nzrange(Δ, colInd)
        rowInd = rows[j]

        if rowInd==colInd
          # diagonal elements
          ix0 = one(T) / (Σx[rowInd,rowInd] * Σy[rowInd,rowInd])
          x1 = A[rowInd,rowInd] - Σy[rowInd,rowInd] + Σx[rowInd,rowInd]
          tmp = shrink(-x1*ix0 + vals[j], λ*Ups[rowInd,colInd]*ix0)
        else
          # off-diagonal elements
          ix0 = one(T) / ((Σx[rowInd,rowInd]*Σy[colInd,colInd] + Σx[colInd,colInd]*Σy[rowInd,rowInd])/2.
                            + Σx[rowInd,colInd]*Σy[rowInd,colInd])
          x1 = A[rowInd,colInd] + A[colInd,rowInd] - 2. * (Σy[rowInd,colInd] - Σx[rowInd,colInd])
          tmp = shrink(-x1*ix0/2. + vals[j], λ*Ups[rowInd,colInd] * ix0)
        end

        # size of update
        h = tmp - vals[j]
        # update is too big -- not done
        if abs(h) > maxChangeTol
          fDone = false
        end
        # update Δ
        vals[j] = tmp

        # update A -- only active set
        for ci=1:p
          for k in nzrange(Δ, ci)
            ri = rows[k]

            if rowInd == colInd
              A[ri,ci] = A[ri,ci] + h * Σx[ri,rowInd]*Σy[rowInd,ci]
              if ri != ci
                A[ci,ri] = A[ci,ri] + h * Σx[ci,rowInd]*Σy[rowInd,ri]
              end
            else
              A[ri,ci] = A[ri,ci] + h * (Σx[ri,colInd]*Σy[rowInd,ci] + Σx[ri,rowInd]*Σy[colInd,ci])
              if ri != ci
                A[ci,ri] = A[ci,ri] + h * (Σx[ci,colInd]*Σy[rowInd,ri] + Σx[ci,rowInd]*Σy[colInd,ri])
              end
            end
          end
        end   # end update A
      end
    end

    if fDone
      break
    end
  end  # while

  sparse(Δ)
end

function findViolator!(
  Δ::SparseMatrixCSC{T},
  A::StridedMatrix{T}, Σx::StridedMatrix{T}, Σy::StridedMatrix{T},
  λ::T, Ups::StridedMatrix{T},
  options::CDOptions = CDOptions()) where {T<:AbstractFloat}

  p = size(Σx, 1)

  vmax = zero(T)
  im = 0
  for c=1:p
    for r=c:p
      t = abs((A[r,c] + A[c,r]) / 2. - (Σy[r,c] - Σx[r,c])) - λ * Ups[r,c]
      if t > vmax
        vmax = t
        im = (c-1)*p+r
      end
    end
  end
  if vmax > options.optTol
    Δ[im] = eps(T)
  else
    im = 0
  end
  im
end

function updateA!(
  A::StridedMatrix{T},
  Δ::SparseMatrixCSC{T},
  Σx::StridedMatrix{T}, Σy::StridedMatrix{T}) where {T<:AbstractFloat}

  p = size(A, 1)

  rows = rowvals(Δ)
  vals = nonzeros(Δ)

  for ac=1:p
    for ar=1:p
      v = zero(T)
      #
      for ci=1:p
        for j in nzrange(Δ, ci)
          ri = rows[j]
          if ri == ci
            v += Σx[ri, ar] * Σy[ci, ac] * vals[j]
          else
            v += (Σx[ri, ar] * Σy[ci, ac] + Σx[ci, ar] * Σy[ri, ac]) * vals[j]
          end
        end
      end
      #
      A[ar,ac] = v
    end
  end
  A
end

differencePrecisionActiveShooting(
  Σx::StridedMatrix{T}, Σy::StridedMatrix{T},
  λ::T, Ups::StridedMatrix{T},
  options::CDOptions = CDOptions()) where {T<:AbstractFloat} =
    differencePrecisionActiveShooting!(spzeros(size(Σx)...), Σx, Σy, λ, Ups, options)

function differencePrecisionActiveShooting!(
  Δ::SparseMatrixCSC{T},
  Σx::StridedMatrix{T}, Σy::StridedMatrix{T},
  λ::T, Ups::StridedMatrix{T},
  options::CDOptions = CDOptions()) where {T<:AbstractFloat}

  maxIter = options.maxIter
  p = size(Σx, 1)
  A = zeros(p,p)
  tril!(Δ)

  # lDelta = tril(Δ, -1)
  # dDelta = spdiagm(diag(Δ))
  # A = Σx * (lDelta + lDelta' + dDelta) * Σy
  updateA!(A, Δ, Σx, Σy)
  if iszero(Δ)
    ind = findViolator!(Δ, A, Σx, Σy, λ, Ups, options)
  end

  for iter=1:maxIter
    updateDelta!(Δ, A, Σx, Σy, λ, Ups, options)
    updateA!(A, Δ, Σx, Σy)

    if findViolator!(Δ, A, Σx, Σy, λ, Ups, options) == 0
      break
    end
  end

  lDelta = tril(Δ, -1)
  dDelta = spdiagm(diag(Δ))
  Δ = (lDelta + lDelta') + dDelta
end

end


text:
using Plot3D: read_multi_block, Block, find_block_connectivity
using JLD2
using JSON

## Cell 1 Reads the block data
if !isfile("blocks.jld")
    blocks = read_multi_block("../../../testfiles/finalmesh.xyz")
    @save "blocks.jld" blocks
end
##
@load "blocks.jld" blocks

## Cell 2 Read the Block data
print("Block read, Finding matching faces\n")
block_match_indices, block_match_corners = find_block_connectivity(blocks[1],blocks[2])
json_string = JSON.json(block_match_indices)
open("block_match_indices.json","w") do f
    JSON.print(f, json_string, 4)
end

json_string = JSON.json(block_match_corners)
open("block_match_corners.json","w") do f
    JSON.print(f, json_string)
end

# @load "block_match_indices.jld" block_match_indices
# @load "block_match_corners.jld" block_match_corners

print("check")
##


text:
@recipe function f(::Val{:InvertedCoax}, dim::Symbol;
                    r=missing,
                    φ=missing,
                    z=missing,
                    half=true)
    T = typeof(d.crystal_radius)

    cross_section::Symbol, value::T = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 0
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, φ
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, r
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, z
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end

    if cross_section == :φ
        points_to_connect::AbstractVector = [SVector{2,T}(0.0,0.0)]
        add_point!(points_to_connect, d.groove_rInner,T(0.0))
        add_point!(points_to_connect, d.groove_rInner,d.groove_depth)
        add_point!(points_to_connect, d.groove_rInner+d.groove_width,d.groove_depth)
        add_point!(points_to_connect, d.groove_rInner+d.groove_width,T(0.0))
        add_point!(points_to_connect, d.crystal_radius,T(0.0))
        add_point!(points_to_connect, d.crystal_radius,d.crystal_length-d.taper_outer_length)
        add_point!(points_to_connect, d.taper_outer_rInner,d.crystal_length)
        add_point!(points_to_connect, d.borehole_radius,d.crystal_length)
        add_point!(points_to_connect, d.borehole_radius,d.crystal_length-d.borehole_length)
        add_point!(points_to_connect, T(0.0),d.crystal_length-d.borehole_length)
        if half == false
            for ipoint in points_to_connect[end-1:-1:1 ]
                add_point!(points_to_connect,-1*ipoint[1],ipoint[2])
            end
        end
        @series begin
            lw --> 1
            label --> ""
            linestyle --> :dash
            color --> :red
            points_to_connect
        end
    end
end

@recipe function f(::Val{:BEGe}, dim::Symbol;
                    r=missing,
                    φ=missing,
                    z=missing,
                    half=true)
    T = typeof(d.crystal_radius)

    cross_section::Symbol, value::T = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 0
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, φ
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, r
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, z
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end

    if cross_section == :φ
        points_to_connect::AbstractVector = [SVector{2,T}(0.0,0.0)]
        add_point!(points_to_connect, d.groove_rInner,T(0.0))
        if d.groove_endplate == "bot"
            add_point!(points_to_connect, d.groove_rInner,T(0.0))
            add_point!(points_to_connect, d.groove_rInner,T(0.0)+d.groove_depth)
            add_point!(points_to_connect, d.groove_rInner+d.groove_width,T(0.0)+d.groove_depth)
            add_point!(points_to_connect, d.groove_rInner+d.groove_width,T(0.0))
        end
        add_point!(points_to_connect, d.taper_bot_rInner,T(0.0))
        add_point!(points_to_connect, d.crystal_radius,d.taper_bot_length)
        add_point!(points_to_connect, d.crystal_radius,d.crystal_length-d.taper_top_length)
        add_point!(points_to_connect, d.taper_top_rInner,d.crystal_length)
        if d.groove_endplate == "top"
            add_point!(points_to_connect, d.groove_rInner+d.groove_width,d.crystal_length)
            add_point!(points_to_connect, d.groove_rInner+d.groove_width,d.crystal_length-d.groove_depth)
            add_point!(points_to_connect, d.groove_rInner,d.crystal_length-d.groove_depth)
            add_point!(points_to_connect, d.groove_rInner,d.crystal_length)
        end
        add_point!(points_to_connect, T(0.0),d.crystal_length)
        if half == false
            for ipoint in points_to_connect[end-1:-1:1 ]
                add_point!(points_to_connect,-1*ipoint[1],ipoint[2])
            end
        end
        @series begin
            lw --> 1
            label --> ""
            linestyle --> :dash
            color --> :red
            points_to_connect
        end
    end
end

function add_point!(points_to_connect::AbstractVector{<:SVector{2,T}},x::T,y::T) where T <: Real
    push!(points_to_connect,SVector{2,T}(x,y))
end
@recipe function f(::Val{:Coax})
    size:=(800,800)
    linewidth:=2
    linecolor:=:black
    label-->""
    # aspect_ratio := equal
    f=1/c.geometry_unit_factor
    @series begin
        partialcircle_3d(f*c.crystal_radius,0,2π,[0,0,0])
    end
    @series begin
        partialcircle_3d(f*c.crystal_radius,0,2π,[0,0,f*c.crystal_length])
    end
    @series begin
        partialcircle_3d(f*c.taper_inner_top_rOuter,0,2π,[0,0,f*c.crystal_length])
    end
    linestyle:=:dash
    @series begin
        partialcircle_3d(f*c.borehole_radius,0,2π,[0,0,f*c.crystal_length-f*c.taper_inner_top_length])
    end
    @series begin
        partialcircle_3d(f*c.taper_inner_bot_rOuter,0,2π,[0,0,0])
    end
    @series begin
        partialcircle_3d(f*c.borehole_radius,0,2π,[0,0,f*c.taper_inner_bot_length])
    end

    for i_rep_seg in 2:size(c.segmentation_r_ranges,1)
        @series begin
            partialcircle_3d(f*c.crystal_radius,c.segmentation_phi_ranges[i_rep_seg][1],c.segmentation_phi_ranges[i_rep_seg][2],[0,0,c.segmentation_z_ranges[i_rep_seg][1]*f])
        end
        @series begin
            partialcircle_3d(f*c.crystal_radius,c.segmentation_phi_ranges[i_rep_seg][1],c.segmentation_phi_ranges[i_rep_seg][2],[0,0,c.segmentation_z_ranges[i_rep_seg][2]*f])
        end
        @series begin
            x::Vector{AbstractFloat}=[f*c.crystal_radius*cos.(c.segmentation_phi_ranges[i_rep_seg][1]),f*c.crystal_radius*cos.(c.segmentation_phi_ranges[i_rep_seg][1])]
            y::Vector{AbstractFloat}=[f*c.crystal_radius*sin.(c.segmentation_phi_ranges[i_rep_seg][1]),f*c.crystal_radius*sin.(c.segmentation_phi_ranges[i_rep_seg][1])]
            z::Vector{AbstractFloat}=[c.segmentation_z_ranges[i_rep_seg][1]*f,c.segmentation_z_ranges[i_rep_seg][2]*f]
            x,y,z
        end
        @series begin
            x=[f*c.crystal_radius*cos.(c.segmentation_phi_ranges[i_rep_seg][2]),f*c.crystal_radius*cos.(c.segmentation_phi_ranges[i_rep_seg][2])]
            y=[f*c.crystal_radius*sin.(c.segmentation_phi_ranges[i_rep_seg][2]),f*c.crystal_radius*sin.(c.segmentation_phi_ranges[i_rep_seg][2])]
            z=[c.segmentation_z_ranges[i_rep_seg][1]*f,c.segmentation_z_ranges[i_rep_seg][2]*f]
            x,y,z
        end
    end
end

struct outer_taper{T <: SSDFloat}
    rStart::T
    rStop::T
    φStart::T
    φStop::T
    zStart::T
    zStop::T
    orientation::String
    function outer_taper(rStart::T, rStop::T, φStart::T, φStop::T, zStart::T, zStop::T, orientation) where T <: SSDFloat
        return new{T}(rStart, rStop, φStart, φStop, zStart, zStop, orientation)
    end
end

function outer_taper(rStart, rStop, φStart, φStop, zStart, zStop, orientation)
    return outer_taper{typeof(rStart)}(rStart, rStop, φStart, φStop, zStart, zStop, orientation)
end

@recipe function f(::Val{:InvertedCoax}; coloring=[], labeling=[])
    if b.name =="ExampleInvertedCoax"
        coloring = [:blue, :orange, :orange, :orange, :orange, :orange, :orange, :orange]
        labeling = ["Core", "Mantle", "", "", "", "", "", ""]
    end
    # aspect_ratio := equal
    f=1/b.geometry_unit_factor
    for i_part in 1:size(b.segmentation_r_ranges,1)
        rStart = f*b.segmentation_r_ranges[i_part][1]
        rStop = f*b.segmentation_r_ranges[i_part][2]
        ΘStart = b.segmentation_phi_ranges[i_part][1]
        ΘStop = b.segmentation_phi_ranges[i_part][2]
        zStart = f*b.segmentation_z_ranges[i_part][1]
        zStop = f*b.segmentation_z_ranges[i_part][2]
        if !isempty(coloring)
            color := coloring[i_part]
        else
            color := :black
        end
        if b.segmentation_types[i_part]=="Tubs"
            if !isempty(labeling)
                label := labeling[i_part]
            else
                label := ""
            end
            if b.borehole_modulation == true && i_part in [b.borehole_segment_idx,b.borehole_bot_segment_idx]
                @series begin
                    phirange = [i for i in 0:0.05:2π+0.5]
                    xs = [(b.borehole_ModulationFunction(phirange[i])+b.borehole_radius) * cos(phirange[i]) for i in eachindex(phirange)] .* f
                    ys = [(b.borehole_ModulationFunction(phirange[i])+b.borehole_radius)* sin(phirange[i]) for i in eachindex(phirange)] .* f
                    zs = [zStart for i in eachindex(phirange)]
                    xs, ys, zs
                end
            elseif b.borehole_modulation == true && i_part == b.borehole_top_segment_idx
                @series begin
                    phirange = [i for i in 0:0.05:2π+0.5]
                    xs = [(b.borehole_ModulationFunction(phirange[i])+b.borehole_radius) * cos(phirange[i]) for i in eachindex(phirange)] .* f
                    ys = [(b.borehole_ModulationFunction(phirange[i])+b.borehole_radius)* sin(phirange[i]) for i in eachindex(phirange)] .* f
                    zs = [zStart for i in eachindex(phirange)]
                    xs, ys, zs
                end
                @series begin
                    phirange = [i for i in 0:0.05:2π+0.5]
                    xs = [rStop * cos(phirange[i]) for i in eachindex(phirange)]
                    ys = [rStop * sin(phirange[i]) for i in eachindex(phirange)]
                    zs = [zStart for i in eachindex(phirange)]
                    xs, ys, zs
                end
            else
                @series begin
                    Tubs("",1,1.,"",0.0,rStart,rStop,ΘStart,ΘStop,zStart,zStop)
                end
            end
        else
            @series begin
                if !isempty(labeling)
                    label := labeling[i_part]
                else
                    label := ""
                end
                # Taper("",1,1.,"",0.0,rStart,rStop,ΘStart,ΘStop,zStart,zStop,"bl")
                outer_taper(rStart,rStop,ΘStart,ΘStop,zStart,zStop,"c//")
            end
        end
    end
end

@recipe function f(::Val{:BEGe}; coloring=[], labeling=[])
    if b.name =="ExampleSegmentedBEGe"
        coloring = [:blue, :red, :purple, :orange, :purple, :grey, :purple,
                            :red, :purple, :orange, :purple, :grey, :purple,
                            :red, :purple, :orange, :purple, :grey, :purple, :purple]
        labeling = ["Core", "Seg. 1", "", "Seg. 2", "", "Seg.3", "Seg. 4",
                            "","","","","","",
                            "","","","","","",""]
        xlims --> (-45,45)
        ylims --> (-45,45)
        zlims --> (-25,65)
    end
    # aspect_ratio := equal
    f=1/b.geometry_unit_factor
    for i_part in 1:size(b.segmentation_r_ranges,1)
        rStart = f*b.segmentation_r_ranges[i_part][1]
        rStop = f*b.segmentation_r_ranges[i_part][2]
        φStart = b.segmentation_phi_ranges[i_part][1]
        φStop = b.segmentation_phi_ranges[i_part][2]
        zStart = f*b.segmentation_z_ranges[i_part][1]
        zStop = f*b.segmentation_z_ranges[i_part][2]
        if !isempty(coloring)
            color := coloring[i_part]
        else
            color := :black
        end
        if b.segmentation_types[i_part]=="Tubs"
            @series begin
                if !isempty(labeling)
                    label := labeling[i_part]
                else
                    label := ""
                end
                Tubs("",1,1.,"",0.0,rStart,rStop,φStart,φStop,zStart,zStop)
            end
        else
            @series begin
                if !isempty(labeling)
                    label := labeling[i_part]
                else
                    label := ""
                end
                # Taper("",1,1.,"",0.0,rStart,rStop,φStart,φStop,zStart,zStop,"bl")
                outer_taper(rStart,rStop,φStart,φStop,zStart,zStop,"c//")
            end
        end
    end
end

function mylinspace(Start,Stop,nSteps)
    return collect(Start:(Stop-Start)/nSteps:Stop)
end

function partialcircle(radius,phiStart,phiStop,Translate::Vector=[0.0,0.0,0.0];nSteps=400)
    phirange = mylinspace(phiStart,phiStop,nSteps)
    x::Vector{AbstractFloat}=map(x->radius*cos.(x),phirange)
    y::Vector{AbstractFloat}=map(x->radius*sin.(x),phirange)
    x.+=Translate[1]
    y.+=Translate[2]
    return x,y
end

function line_3d(r1::T, r2::T, phi1::T, phi2::T, z1::T, z2::T; translate::AbstractVector = SVector{3,T}(0.0,0.0,0.0)) where T <:Real
    x1::SVector{3,T} = [r1*cos(phi1),r1*sin(phi1),z1]+translate
    x2::SVector{3,T} = [r2*cos(phi2),r2*sin(phi2),z2]+translate
    return [x1[1],x2[1]] , [x1[2],x2[2]], [x1[3],x2[3]]
end

@recipe function f(contact::AbstractContact{T}) where T
    c-->:orange
    for (i,g) in enumerate(vcat(contact.geometry_positive, contact.geometry_negative))
        @series begin
            i==1 ? label --> "$(contact.id)" : label := ""
            g
        end
    end
end

@recipe function f(c::SolidStateDetector{T}; coloring = [], labeling = []) where T
    legendfont --> 16
    # aspect_ratio := 1
    tickfontsize --> 9
    guidefontsize --> 13
    lw --> 2
    # xlabel --> "\n x / m"
    # ylabel --> "\n y / m"
    # zlabel --> "\n z / m"
    camera --> (15,25)
    coloring = coloring
    labeling = labeling
        if c.name == "Public Segmented BEGe"
        @series begin
            # zlims --> (-0.023,0.063)
            # aspect_ratio --> 0.5
            c, Val(:segBEGe)
        end
    elseif c.name == "Public Inverted Coax"
        @series begin
            c, Val(:ivc)
        end
    elseif c.name == "Public Coax"
        @series begin
            c, Val(:coax)
        end
    else
        for (ic,contact) in enumerate(c.contacts)
            @series begin
                !isempty(coloring) ? c --> coloring[ic] : nothing
                !isempty(labeling) ? label --> labeling[ic] : nothing
                contact
            end
        end
    end
end

@recipe function f(c,a::Val{:ivc})
    labeling = ["Core","Mantle"]
    coloring = [:blue,:orange]
    for (ic, contact) in enumerate(c.contacts)
        @series begin
            color --> coloring[ic]
            label --> labeling[ic]
            contact
        end
    end
end
@recipe function f(c,a::Val{:coax})
    labeling = ["" for i in 1:19]
    coloring = vcat([:blue], [:orange for i in 1:18])
    for (ic, contact) in enumerate(c.contacts)
        @series begin
            color --> coloring[ic]
            label --> labeling[ic]
            contact
        end
    end
end

@recipe function f(c,a::Val{:segBEGe}, coloring = [])
    isempty(coloring) ? coloring = [:blue, :red, :orange, :grey, :purple] : nothing
    labeling = [ "Core", "Seg. 1", "Seg. 2", "Seg. 3", "Seg. 4"]
    for (ic, contact) in enumerate(c.contacts)
        @series begin
            color --> coloring[ic]
            label --> labeling[ic]
            contact
        end
    end
end

@recipe function f(geometry::AbstractGeometry)
    pos, neg = get_decomposed_volumes(geometry)
    for g in pos
        @series begin
            g
        end
    end
end
@recipe function f(object::AbstractObject)
    @series begin
        object.geometry
    end
end

## 2D from Felix

function line_2d(r1,r2,z1,z2)
   x1::Vector = [r1,z1]
   x2::Vector = [r2,z2]
   return [x1[1],x2[1]],[x1[2],x2[2]]
end

function polar_circle_2d(r::T, phiStart::T, phiStop::T;nSteps=360) where T <: SSDFloat
    if r == 0; return [],[]; end
    phiRange = collect(phiStart:(phiStop-phiStart)/nSteps:phiStop)
    rRange = [r for a in 1:length(phiRange)]
    return phiRange, rRange
end

@recipe function f(d::SolidStateDetector{T}, dim::Symbol; φ = missing, z = missing) where{T <: SSDFloat}
    if ismissing(z) && ismissing(φ)
        print("Please specify φ or z.")
    elseif ismissing(z)
        for c in d.contacts
            for g in c.geometry_positive
                @series begin
                    if d.name == "Public Inverted Coax"
                        if typeof(c) == SSD.Contact{T,:N}; color --> :orange
                        elseif typeof(c) == SSD.Contact{T,:P}; color --> :blue; end
                    else color --> :black
                    end
                    lw --> 2
                    label = ""
                    g, :φ, T(deg2rad(mod(φ,2π)))
                end
            end
        end
    elseif ismissing(φ)
        proj --> :polar
        for c in d.contacts
            for g in c.geometry_positive
                @series begin
                    if d.name == "Public Inverted Coax"
                        if typeof(c) == SSD.Contact{T,:N}; color --> :orange
                        elseif typeof(c) == SSD.Contact{T,:P}; color --> :blue; end
                    else color --> :black
                    end
                    lw --> 2
                    label = ""
                    g, :z, T(z)
                end
            end
        end
    else
        print("Please specify only one of either φ or z.")
    end
end

@recipe function f(Vol::SSD.Tube{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
    if dim == :φ
        if parameter in Vol.φ_interval
           rStart = Vol.r_interval.left
           rStop = Vol.r_interval.right
           zStart = Vol.z_interval.left
           zStop = Vol.z_interval.right

           @series begin
                   label --> ""
                   line_2d(rStop,rStop,zStart,zStop)
               end

            @series begin
                   label --> ""
                   line_2d(rStart,rStop,zStart,zStart)
            end

           if rStart != rStop
                @series begin
                   label --> ""
                   line_2d(rStart,rStop,zStop,zStop)
               end
            end

            if zStart != zStop && rStart != 0
                @series begin
                   label --> ""
                   line_2d(rStart,rStart,zStart,zStop)
               end
            end
        end

    elseif dim == :z
        if parameter in Vol.z_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            phiStart = Vol.φ_interval.left
            phiStop = Vol.φ_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right

            if rStart != 0
                @series begin
                    label --> ""
                    polar_circle_2d(rStart,phiStart,phiStop)
                end
            end
            @series begin
                label --> ""
                polar_circle_2d(rStop,phiStart,phiStop)
            end

            if !isapprox(mod(phiStart,2π),mod(phiStop,2π),atol = 1e-5)
                @series begin
                    label --> ""
                    line_2d(phiStart,phiStart,rStart,rStop)
                end
                @series begin
                    label --> ""
                    line_2d(phiStop,phiStop,rStart,rStop)
                end
            end
        end
    end
end

# @recipe function f(Vol::SSD.ConeMantle{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
#     newVol = Vol.cone
#     @series begin
#         newVol, dim, parameter
#     end
# end

@recipe function f(Vol::SSD.Cone{T}, dim::Symbol, parameter::T) where{T <: SSDFloat}
    if dim == :φ
        if parameter in Vol.φ_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right
            orientation = Vol.orientation

            if orientation in [:top_right, :bottom_left]
                @series begin
                    label --> ""
                   line_2d(rStop,rStart,zStart,zStop)
                end
            else
                @series begin
                    label --> ""
                   line_2d(rStart,rStop,zStart,zStop)
                end
            end
        end

    elseif dim == :z
        if parameter in Vol.z_interval
            rStart = Vol.r_interval.left
            rStop = Vol.r_interval.right
            phiStart = Vol.φ_interval.left
            phiStop = Vol.φ_interval.right
            zStart = Vol.z_interval.left
            zStop = Vol.z_interval.right
            orientation = Vol.orientation

            if orientation in [:top_right, :bottom_left]
                @series begin
                    label --> ""
                    r = (rStop-rStart)*(parameter-zStop)/(zStart-zStop) + rStart
                    polar_circle_2d(r,phiStart,phiStop)
                end
            else
                @series begin
                    label --> ""
                    r = (rStop-rStart)*(parameter-zStart)/(zStop-zStart) + rStart
                    polar_circle_2d(r,phiStart,phiStop)
                end
            end
        end

    end
end


text:
module CholeskyExperiment

using LightGraphs
using MetaGraphs
using Mongoc
using CSV
using UUIDs
import JSON
import MPI

include("DecompositionDB/DecompositionDB.jl")
include("Generate/Generate.jl")
include("Solve/Solve.jl")
include("ReadFeatures/ReadFeatures.jl")
include("VisualizeFeatures/VisualizeFeatures.jl")
#include("utils/chordalextension.jl")
include("utils/constructgraph.jl")
#include("utils/misc.jl")
include("utils/maximalcliques.jl")
include("utils/merge.jl")

import .DecompositionDB; const db = DecompositionDB;
import .Generate; const gen = Generate;
import .Solve; const slv = Solve;
import .ReadFeatures; const readf = ReadFeatures;
import .VisualizeFeatures; const viz = VisualizeFeatures;

struct ExperimentManager
    instances::Mongoc.Collection
    decompositions::Mongoc.Collection
end

function ExperimentManager(;port::String="mongodb://localhost:27017", db="choleskyexp")
    client = Mongoc.Client(port)
    try
        Mongoc.ping(client)
    catch err
        if isa(err, MethodError)
            error("You need to have a MongoDB daemon running.")
        end
    end
    instances = DecompositionDB.getcollection(db, "instances", client)
    decompositions = DecompositionDB.getcollection(db, "decompositions", client)
    return ExperimentManager(instances, decompositions)
end

include("./readfeatures.jl")
include("./loadinstance.jl")
include("./generate.jl")
include("./solve.jl")
include("./mpi.jl")
include("./visualizefeatures.jl")

end # module

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
@testset "Duality" begin
    @test dual(A) == A ⋅ inv(I)
    @test dual(A) ⋅ I == A
    for X ∈ (Zs..., NHZs...)
        @test dual(dual(X)) == I^2 * X
    end
end


text:
"""
    This example shows how to use panels in Term.jl

A `Panel` prints a piece of content surrounded by a border.
It's one of the best and easiest way to style your terminal output,
especially if you nest multiple panels or use Term's layout functionality
to create a structured terminal output.
"""

import Term: Panel

# Creating a panel is very simple

print(Panel("This is my first panel!"))

"""
By default Panel expands to fit its content.
You can use different options to change this beahavior.
"""

print("\n\n")
print(
    Panel("this panel has fixed width", width = 44, justify = :right),
    Panel("this one too, but the text is at the center!", width = 66, justify = :center),
    Panel("this one fits its content", width = :fit),
    Panel("this one fits the terminal window size!", width = nothing),
)

"""
As you can see you can also specify how the panel's content should be justified.
Of course you can use markup styled text in your panel.
"""

print("\n\n")
print(
    Panel(
        "[red]This is the panel's first line.[/red]",
        "[bold green]and this is another, panel just stacks all inputs into one piece of content[/bold green]",
        width = :fit,
    ),
)

"""
When it comes to style, panel offers a lot of options.
You can style the box itself and have titles and subtitles.
"""

print("\n\n")

print(
    Panel("content "^10, title = "My Panel", title_style = "bold red", width = 44),
    Panel(
        "content "^10,
        subtitle = "another panel",
        subtitle_style = "dim underline",
        subtitle_justify = :right,
        width = 44,
    ),
    Panel("content "^10, box = :ASCII_DOUBLE_HEAD, style = "red", width = 44),
    Panel("content "^10, box = :DOUBLE, style = "blue", width = 44),
)

"""
Finally, you can layout panels to create structured content
"""

print("\n\n")
pleft = Panel("content "^30, box = :DOUBLE, style = "blue", width = 66)
pright = Panel(
    "content [red]with style[/red] "^26,
    title = "My Panel",
    title_style = "bold red",
    width = 44,
)

print(
    Panel(
        pleft / pright,
        style = "green dim",
        title_style = "green",
        title = "vertically stacked!",
        title_justify = :center,
        subtitle = "styled by Term.jl",
        subtitle_justify = :right,
        subtitle_style = "dim",
        justify = :center,
    ),
)


text:
push!(LOAD_PATH, "Modules")
using Colors
using JuliRay

φ=MathConstants.φ

## Tetrahedron
vertices=[[1.0,1.0,1.0],[1.0,-1.0,-1.0],[-1.0,1.0,-1.0],[-1.0,-1.0,1.0]]
edges=[[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]]
faces=[[2,3,4],[1,3,4],[1,2,4],[1,2,3]]

V=rgbColor(csgUnion((p->Sphere(p,0.05)).(vertices)),RGB(0.1,0.1,0.1))
E=rgbColor(csgUnion([Cylinder(vertices[I...],vertices[J...],0.025) for (I,J) in edges]),RGB(0.2,0.2,0.2))
F=rgbftColor(csgUnion([Polygon((i->vertices[i]).(I)) for I in faces]),RGB(1,0,0),FT(0.2,0.2))
render(csgUnion(V,E,F),camera=LngLatCamera(lng=π/5,lat=π/5,pers=0.2,zoom=0.25),name="PlatonicSolid_4")

## Cube
vertices=[[[i,j,k] for i ∈ [1,-1], j ∈ [1,-1], k ∈ [1,-1]]...]
edges=[[[1+i+j,2+i+j] for i ∈ [0,2], j ∈ [0,4]]...]
edges=vcat(edges,[[[1+i+j,3+i+j] for i ∈ [0,1], j ∈ [0,4]]...])
edges=vcat(edges,[[[1+i+j,5+i+j] for i ∈ [0,1], j ∈ [0,2]]...])
faces=[[1,2,4,3],[1,2,4,3].+4,[1,3,7,5],[1,3,7,5].+1,[1,5,6,2],[1,5,6,2].+2]

V=rgbColor(csgUnion((p->Sphere(p,0.05)).(vertices)),RGB(0.1,0.1,0.1))
E=rgbColor(csgUnion([Cylinder(vertices[I...],vertices[J...],0.025) for (I,J) in edges]),RGB(0.2,0.2,0.2))
F=rgbftColor(csgUnion([Polygon((i->vertices[i]).(I)) for I in faces]),RGB(1,0,0),FT(0.2,0.2))
render(csgUnion(V,E,F),camera=LngLatCamera(lng=π/5,lat=π/5,pers=0.2,zoom=0.25),name="PlatonicSolid_6")

## Octahedron
vertices=[[1,0,0],[-1,0,0],[0,1,0],[0,-1,0],[0,0,1],[0,0,-1]]*sqrt(3)
edges=[[1,3],[3,2],[2,4],[4,1],[3,5],[5,4],[4,6],[6,3],[5,1],[1,6],[6,2],[2,5]]
faces=[[1,3,5],[2,3,5],[2,4,5],[1,4,5],[1,3,6],[2,3,6],[2,4,6],[1,4,6]]

V=rgbColor(csgUnion((p->Sphere(p,0.05)).(vertices)),RGB(0.1,0.1,0.1))
E=rgbColor(csgUnion([Cylinder(vertices[I...],vertices[J...],0.025) for (I,J) in edges]),RGB(0.2,0.2,0.2))
F=rgbftColor(csgUnion([Polygon((i->vertices[i]).(I)) for I in faces]),RGB(1,0,0),FT(0.2,0.2))
render(csgUnion(V,E,F),camera=LngLatCamera(lng=π/5,lat=π/5,pers=0.2,zoom=0.25),name="PlatonicSolid_8")

## Dodecahedron
vertices0=[[[i,j,k] for i ∈ [1,-1], j ∈ [1,-1], k ∈ [1,-1]]...]
vertices1=[[[0,i/φ,j*φ] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices2=[[[j*φ,0,i/φ] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices3=[[[i/φ,j*φ,0] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices=vcat(vertices0,vertices1,vertices2,vertices3)
edges1=[[1,9],[2,9],[3,10],[4,10],[9,10]]
edges2=[[5,11],[6,11],[7,12],[8,12],[11,12]]
edges3=[[1,13],[3,13],[5,14],[7,14],[13,14]]
edges4=[[2,15],[4,15],[6,16],[8,16],[15,16]]
edges5=[[1,17],[5,17],[2,18],[6,18],[17,18]]
edges6=[[3,19],[7,19],[4,20],[8,20],[19,20]]
edges=vcat(edges1,edges2,edges3,edges4,edges5,edges6)
faces=[[1,9,10,3,13],[2,9,10,4,15],[5,11,12,7,14],[6,11,12,8,16],[1,13,14,5,17],[3,13,14,7,19],[2,15,16,6,18],[4,15,16,8,20],[1,17,18,2,9],[5,17,18,6,11],[3,19,20,4,10],[7,19,20,8,12]]

V=rgbColor(csgUnion((p->Sphere(p,0.05)).(vertices)),RGB(0.1,0.1,0.1))
E=rgbColor(csgUnion([Cylinder(vertices[I...],vertices[J...],0.025) for (I,J) in edges]),RGB(0.2,0.2,0.2))
F=rgbftColor(csgUnion([Polygon((i->vertices[i]).(I)) for I in faces]),RGB(1,0,0),FT(0.2,0.2))
render(csgUnion(V,E,F),camera=LngLatCamera(lng=π/5,lat=π/5,pers=0.2,zoom=0.25),name="PlatonicSolid_12")

## Icosahedron
vertices1=[[[0,i,j*φ] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices2=[[[j*φ,0,i] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices3=[[[i,j*φ,0] for i ∈ [1,-1], j ∈ [1,-1]]...]
vertices=vcat(vertices1,vertices2,vertices3)*sqrt(3)/√(φ+2)
edges=[[2i-1,2i] for i ∈ 1:6]
edges=[edges...,[1,5],[5,9],[9,1],[2,5],[5,11],[11,2],[3,6],[6,9],[9,3],[4,6],[6,11],[11,4],[1,7],[7,10],[10,1],[2,7],[7,12],[12,2],[3,8],[8,10],[10,3],[4,8],[8,12],[12,4]]
faces=[[1,5,9],[2,5,11],[3,6,9],[4,6,11],[1,7,10],[2,7,12],[3,8,10],[4,8,12],[1,2,5],[1,2,7],[3,4,6],[3,4,8],[5,6,9],[5,6,11],[7,8,10],[7,8,12],[9,10,1],[9,10,3],[11,12,2],[11,12,4]]

V=rgbColor(csgUnion((p->Sphere(p,0.05)).(vertices)),RGB(0.1,0.1,0.1))
E=rgbColor(csgUnion([Cylinder(vertices[I...],vertices[J...],0.025) for (I,J) in edges]),RGB(0.2,0.2,0.2))
F=rgbftColor(csgUnion([Polygon((i->vertices[i]).(I)) for I in faces]),RGB(1,0,0),FT(0.2,0.2))
render(csgUnion(V,E,F),camera=LngLatCamera(lng=π/5,lat=π/5,pers=0.2,zoom=0.25),name="PlatonicSolid_20")

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
using SpecialFunctions, NearestNeighbors, Distances, Distributions, Random

"""
    n_ball(n::Number)
Computes the volume of a n-dimensional unit sphere.
"""
function n_ball(n::Number)
    return π^(n/2.) / gamma(n/2. + 1.)
end

"""
    kl_entropy(data::Array{Float64, 2}; k=5)
Compute the nearest-neighbor estimate of the differential entropy of data.

data is a 2d array, with every column representing one data point.
For further information, see

"A class of Rényi information estimators for multidimensional densities"
Nikolai Leonenko, Luc Pronzato, and Vippal Savani
The Annals of Statistics, 2008
https://projecteuclid.org/euclid.aos/1223908088

keyword arguments:
k=5: number of nearest neighbors
"""
function kl_entropy(data::Array{Float64, 2}; k=5)
    d, N = size(data)
    kdtree = KDTree(data)
    _, dist = knn(kdtree, data, k+1, true)
    H = log(N) - digamma(k) + log(n_ball(d)) + d/N * sum(map(l->log(l[end]), dist))
end

"""
    kl_mutual_information(x, y; k=5, bias_correction=true)
compute the nearest-neighbor 'KGS' estimate of the mutual information between x and y.

x and y are 2d arrays, with every column representing one data point.
For further information, see

"Estimating Mutual Information"
Alexander Kraskov, Harald Stoegbauer, and Peter Grassberger
Physical Review E
https://arxiv.org/pdf/cond-mat/0305641.pdf

"Demystifying Fixed k-Nearest Neighbor Information Estimators"
Weihao Gao, Sewoong Oh, Pramod Viswanath
EEE International Symposium on Information Theory - Proceedings
https://arxiv.org/pdf/1604.03006.pdf

keyword arguments:
k=5: number of nearest neighbors
bias_correction=true: flag to apply Gao's bias correction
"""
function kl_mutual_information(x, y; k=5, bias_correction=true)
    dist = bias_correction ? Euclidean() : Chebyshev()

    d_x, N = size(x)
    d_y, _ = size(y)

    xy = vcat(x,y)

    kdtree_x = KDTree(x, dist)
    kdtree_y = KDTree(y, dist)
    kdtree_xy = KDTree(xy, dist)

    I = 0.

    for i in 1:N
        _, dist = knn(kdtree_xy, xy[:,i], k+1, true)
        n_x = length(inrange(kdtree_x, x[:,i], dist[end], false))-1
        n_y = length(inrange(kdtree_y, y[:,i], dist[end], false))-1
        I += digamma(n_x+1) + digamma(n_y+1)
    end

    I = -I/N + log(N) + digamma(k)

    if bias_correction
        I += log( n_ball(d_x) * n_ball(d_y)/n_ball(d_x+d_y) )
    end

    return I
end

"""
    kl_renyi(data::Array{Float64, 2}, q; k=5)
Compute the nearest-neighbor estimate of the Renyi-alpha entropy of data.

data is a 2d array, with every column representing one data point.
For further information, see

"A class of Rényi information estimators for multidimensional densities"
Nikolai Leonenko, Luc Pronzato, and Vippal Savani
The Annals of Statistics, 2008
https://projecteuclid.org/euclid.aos/1223908088

keyword arguments:
k=5: number of nearest neighbors
"""
function kl_renyi(data, q, k=5)
    d, N = size(data)
    Vd = π^(d/2.) / gamma(d/2. + 1.)
    Ck = ( gamma(k)/gamma(k+1-q) )^(1/(1-q))
    kdtree = KDTree(data)
    _, dist = knn(kdtree, data, k+1, true)
    Iq = 1/N * sum(map(l->( (N-1)*Ck*Vd*l[end]^d )^(1-q), dist))

    return log(Iq)/(1-q)
end

"""
    kl_cond_mi(x, y, z; k=5, bias_correction=true)
compute the nearest-neighbor 'KGS' estimate of the conditional mutual information between x and y given z.

x, y, and z are 2d arrays, with every column representing one data point.
keyword arguments:
k=5: number of nearest neighbors
bias_correction=true: flag to apply Gao's bias correction
"""
function kl_cond_mi(x, y, z; k=5, bias_correction=true)
    dist = bias_correction ? Euclidean() : Chebyshev()

    xz = vcat(x,z)
    yz = vcat(y,z)
    xyz = vcat(x,y,z)

    d_x, N = size(x)
    d_y, _ = size(y)
    d_z, _ = size(z)

    kdtree_z = KDTree(z, dist)
    kdtree_xz = KDTree(xz, dist)
    kdtree_yz = KDTree(yz, dist)
    kdtree_xyz = KDTree(xyz, dist)

    CMI = 0.

    for i in 1:N
        _, dist = knn(kdtree_xyz, xyz[:,i], k+1, true)
        n_xz = length(inrange(kdtree_xz, xz[:,i], dist[end], false))-1
        n_yz = length(inrange(kdtree_yz, yz[:,i], dist[end], false))-1
        n_z = length(inrange(kdtree_z, z[:,i], dist[end], false))-1
        CMI = CMI + digamma(n_xz) + digamma(n_yz) - digamma(n_z)
    end

    CMI = digamma(k) - CMI/N

    if bias_correction
        CMI += log( (n_ball(d_x + d_z) * n_ball(d_y+d_z))/(n_ball(d_x+d_y+d_z)*n_ball(d_z)))
    end

    return CMI
end

"""
    kl_perm_mi_test(x, y; k=5, B=100, bias_correction=true)
compute permutation test of independence of x and y.

keyword arguments:
k=5: number of nearest neighbors to use for mutual information estimate
B=100: number of permutations
bias_correction=true: flag to apply Gao's bias correction
"""
function kl_perm_mi_test(x, y; k=5, B=100, bias_correction=true)
    MI = kl_mutual_information(x, y, k=k, bias_correction=bias_correction)
    samples = Float64[]

    for i in 1:B
        push!(samples, kl_mutual_information(x, y[:, shuffle(1:end)], k=k, bias_correction=bias_correction))
    end

    p = length(filter(d->MI<d, samples))/B
    return p
end

"""
    kl_perm_cond_mi_test(x, y, z; k=5, B=100, kp=5, bias_correction=true)
compute permutation test of conditional independence of x and y given z.

For further information, see:
"Conditional independence testing based on a nearest-neighbor estimator of conditional mutual information"
Jakob Runge
Proceedings of the 21st International Conference on Artificial Intelligence and Statistics (AISTATS) 2018, Lanzarote, Spain.
http://proceedings.mlr.press/v84/runge18a/runge18a.pdf

keyword arguments:
k=5: number of nearest neighbors to use for mutual information estimate
B=100: number of permutations
bias_correction=true: flag to apply Gao's bias correction

"""
function kl_perm_cond_mi_test(x, y, z; k=5, B=100, kp=5, bias_correction=true)
    d, N = size(z)
    CMI = kl_cond_mi(x,y,z,k=k,bias_correction=bias_correction)
    samples = Float64[]

    kdtree_z = KDTree(z, Chebyshev())
    z_knn, _ = knn(kdtree_z, z, kp+1)

    for b in 1:B
        # create permutation for independence test
        U = Int64[]
        P = collect(1:N)
        j = 0
        Ns = map(shuffle, z_knn)

        for i in shuffle(collect(1:N))
            j = Ns[i][1]
            m = 1
            while (j ∈ U) && (m < kp)
                m += 1
                j = Ns[i][m]
            end
            P[i]=j
            push!(U,j)
        end
        push!(samples, kl_cond_mi(x[:,P],y,z,k=k,bias_correction=bias_correction))
    end

    p = length(filter(d->CMI<d, samples))/B

    return p
end


text:
"""
    LogisticSoftMaxLikelihood(num_class::Int)

## Arguments
- `num_class::Int` : Total number of classes

---
The multiclass likelihood with a logistic-softmax mapping: :
```math
p(y=i|{fₖ}₁ᴷ) = σ(fᵢ)/∑ₖ σ(fₖ)
```
where `σ` is the logistic function.
This likelihood has the same properties as [softmax](https://en.wikipedia.org/wiki/Softmax_function).
---

For the analytical version, the likelihood is augmented multiple times.
More details can be found in the paper [Multi-Class Gaussian Process Classification Made Conjugate: Efficient Inference via Data Augmentation](https://arxiv.org/abs/1905.09670)
"""
mutable struct LogisticSoftMaxLikelihood{T<:Real,A<:AbstractVector{T}} <:
               MultiClassLikelihood{T}
    nClasses::Int
    class_mapping::Vector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    Y::Vector{BitVector} #Mapping from instances to classes (one hot encoding)
    y_class::Vector{Int} # GP Index for each sample
    c::Vector{A} # Second moment of fₖ
    α::A # First variational parameter of Gamma distribution
    β::A # Second variational parameter of Gamma distribution
    θ::Vector{A} # Variational parameter of Polya-Gamma distribution
    γ::Vector{A} # Variational parameter of Poisson distribution
    function LogisticSoftMaxLikelihood{T}(nClasses::Int) where {T<:Real}
        return new{T,Vector{T}}(nClasses)
    end
    function LogisticSoftMaxLikelihood{T}(
        nClasses::Int, labels::AbstractVector, ind_mapping::Dict
    ) where {T<:Real}
        return new{T,Vector{T}}(nClasses, labels, ind_mapping)
    end
    function LogisticSoftMaxLikelihood{T}(
        nClasses::Int,
        Y::AbstractVector{<:BitVector},
        class_mapping::AbstractVector,
        ind_mapping::Dict{<:Any,<:Int},
        y_class::AbstractVector{<:Int},
    ) where {T<:Real}
        return new{T,Vector{T}}(nClasses, class_mapping, ind_mapping, Y, y_class)
    end
    function LogisticSoftMaxLikelihood{T}(
        nClasses,
        Y::AbstractVector{<:BitVector},
        class_mapping::AbstractVector,
        ind_mapping::Dict{<:Any,<:Int},
        y_class::AbstractVector{<:Int},
        c::AbstractVector{A},
        α::A,
        β::A,
        θ::AbstractVector{A},
        γ::AbstractVector{A},
    ) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(nClasses, class_mapping, ind_mapping, Y, y_class, c, α, β, θ, γ)
    end
end

LogisticSoftMaxLikelihood(nClasses::Int) = LogisticSoftMaxLikelihood{Float64}(nClasses)
function LogisticSoftMaxLikelihood(ylabels::AbstractVector)
    return LogisticSoftMaxLikelihood{Float64}(
        length(ylabels), ylabels, Dict(value => key for (key, value) in enumerate(ylabels))
    )
end

function implemented(
    ::LogisticSoftMaxLikelihood, ::Union{<:AnalyticVI,<:MCIntegrationVI,<:GibbsSampling}
)
    return true
end

function logisticsoftmax(f::AbstractVector{<:Real})
    return normalize(logistic.(f), 1)
end

function logisticsoftmax(f::AbstractVector{<:Real}, i::Integer)
    return logisticsoftmax(f)[i]
end

function (::LogisticSoftMaxLikelihood)(f::AbstractVector)
    return logisticsoftmax(f)
end

function (::LogisticSoftMaxLikelihood)(y::Integer, f::AbstractVector)
    return logisticsoftmax(f)[y]
end

function Base.show(io::IO, model::LogisticSoftMaxLikelihood{T}) where {T}
    return print(io, "Logistic-Softmax Likelihood")
end

function init_likelihood(
    l::LogisticSoftMaxLikelihood{T},
    inference::AbstractInference{T},
    nLatent::Integer,
    nSamplesUsed::Integer,
) where {T}
    if inference isa AnalyticVI || inference isa GibbsSampling
        c = [ones(T, nSamplesUsed) for i in 1:nLatent]
        α = nLatent * ones(T, nSamplesUsed)
        β = nLatent * ones(T, nSamplesUsed)
        θ = [abs.(rand(T, nSamplesUsed)) * 2 for i in 1:nLatent]
        γ = [abs.(rand(T, nSamplesUsed)) for i in 1:nLatent]
        LogisticSoftMaxLikelihood{T}(
            num_class(l), l.Y, l.class_mapping, l.ind_mapping, l.y_class, c, α, β, θ, γ
        )
    else
        return l
    end
end

## Local Updates ##
function local_updates!(
    l::LogisticSoftMaxLikelihood,
    y,
    μ::NTuple{N,<:AbstractVector},
    Σ::NTuple{N,<:AbstractVector},
) where {N}
    @. l.c = broadcast((Σ, μ) -> sqrt.(Σ + abs2.(μ)), Σ, μ)
    for _ in 1:2
        broadcast!(
            (β, c, μ, ψα) -> 0.5 / β * exp.(ψα) .* safe_expcosh.(-0.5 * μ, 0.5 * c),
            l.γ,
            Ref(l.β),
            l.c,
            μ,
            Ref(digamma.(l.α)),
        )
        l.α .= 1.0 .+ (l.γ...)
    end
    broadcast!(
        (y, γ, c) -> 0.5 * (y + γ) ./ c .* tanh.(0.5 .* c),
        l.θ, # target
        y, # argument 1
        l.γ, # argument 2
        l.c, # argument 3
    )
    return nothing
end

function sample_local!(l::LogisticSoftMaxLikelihood{T}, y::AbstractVector, f) where {T}
    broadcast!(f -> rand.(Poisson.(0.5 * l.α .* safe_expcosh.(-0.5 * f, 0.5 * f))), l.γ, f)
    l.α .= rand.(Gamma.(one(T) .+ (l.γ...), 1.0 ./ l.β))
    set_ω!(l, broadcast((y, γ, f) -> rand.(PolyaGamma.(y .+ Int.(γ), abs.(f))), y, l.γ, f))
    return nothing
end

## Global Gradient Section ##

@inline function ∇E_μ(l::LogisticSoftMaxLikelihood, ::AOptimizer, y::AbstractVector)
    return 0.5 .* (y .- l.γ)
end
@inline ∇E_Σ(l::LogisticSoftMaxLikelihood, ::AOptimizer, ::AbstractVector) = 0.5 .* l.θ

## ELBO Section ##
function expec_loglikelihood(
    l::LogisticSoftMaxLikelihood{T}, ::AnalyticVI, y, μ, Σ
) where {T}
    tot = -length(y) * logtwo
    tot += -sum(sum(l.γ .+ y)) * logtwo
    tot += 0.5 * sum(zip(l.θ, l.γ, y, μ, Σ)) do (θ, γ, y, μ, Σ)
        dot(μ, (y - γ)) - dot(θ, abs2.(μ)) - dot(θ, Σ)
    end
    return tot
end

function AugmentedKL(l::LogisticSoftMaxLikelihood, y::AbstractVector)
    return PolyaGammaKL(l, y) + PoissonKL(l) + GammaEntropy(l)
end

function PolyaGammaKL(l::LogisticSoftMaxLikelihood, y)
    return sum(broadcast(PolyaGammaKL, y .+ l.γ, l.c, l.θ))
end

function PoissonKL(l::LogisticSoftMaxLikelihood)
    return sum(broadcast(PoissonKL, l.γ, Ref(l.α ./ l.β), Ref(digamma.(l.α) .- log.(l.β))))
end

##  Compute the equivalent of KL divergence between an improper prior p(λ) (``1_{[0,\\infty]}``) and a variational Gamma distribution ##
function GammaEntropy(l::LogisticSoftMaxLikelihood)
    return -sum(l.α) + sum(log, first(l.β)) - sum(x -> first(logabsgamma(x)), l.α) -
           dot(1.0 .- l.α, digamma.(l.α))
end

## Numerical Gradient Section ##

function grad_samples(
    model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:NumericalVI},
    samples::AbstractMatrix{T},
    index::Int,
) where {T}
    class = model.likelihood.y_class[index]::Int
    grad_μ = zeros(T, nLatent(model))
    grad_Σ = zeros(T, nLatent(model))
    g_μ = similar(grad_μ)
    nSamples = size(samples, 1)
    @views @inbounds for i in 1:nSamples
        σ = logistic.(samples[i, :])
        samples[i, :] .= logisticsoftmax(samples[i, :])
        s = samples[i, class]
        g_μ .= grad_logisticsoftmax(samples[i, :], σ, class) / s
        grad_μ += g_μ
        grad_Σ += diaghessian_logisticsoftmax(samples[i, :], σ, class) / s - abs2.(g_μ)
    end
    for k in 1:nLatent(model)
        get_opt(inference(model), k).ν[index] = -grad_μ[k] / nSamples
        get_opt(inference(model), k).λ[index] = grad_Σ[k] / nSamples
    end
end

function log_like_samples(
    model::AbstractGP{T,<:LogisticSoftMaxLikelihood},
    samples::AbstractMatrix,
    index::Integer,
) where {T}
    class = model.likelihood.y_class[index]
    nSamples = size(samples, 1)
    loglike = zero(T)
    for i in 1:nSamples
        σ = logistic.(samples[i, :])
        loglike += log(σ[class]) - log(sum(σ))
    end
    return loglike / nSamples
end

function grad_logisticsoftmax(
    s::AbstractVector{T}, σ::AbstractVector{T}, i::Integer
) where {T<:Real}
    return s[i] * (δ.(T, i, eachindex(σ)) .- s) .* (1.0 .- σ)
end

function diaghessian_logisticsoftmax(
    s::AbstractVector{T}, σ::AbstractVector{T}, i::Integer
) where {T<:Real}
    return s[i] * (1.0 .- σ) .* (
        abs2.(δ.(T, i, eachindex(σ)) - s) .* (1.0 .- σ) - s .* (1.0 .- s) .* (1.0 .- σ) -
        σ .* (δ.(T, i, eachindex(σ)) - s)
    )
end

function hessian_logisticsoftmax(
    s::AbstractVector{T}, σ::AbstractVector{T}, i::Integer
) where {T<:Real}
    m = length(s)
    hessian = zeros(T, m, m)
    @inbounds for j in 1:m
        for k in 1:m
            hessian[j, k] =
                (1 - σ[j]) *
                s[i] *
                (
                    (δ(T, i, k) - s[k]) * (1.0 - σ[k]) * (δ(T, i, j) - s[j]) -
                    s[j] * (δ(T, j, k) - s[k]) * (1.0 - σ[k]) -
                    δ(T, k, j) * σ[j] * (δ(T, i, j) - s[j])
                )
        end
    end
    return hessian
end


text:
using Base.Test
using UsingOverride

data = rand(100)

#load M1 for the first time
module M1

    export T1, x, do_thing

    immutable T1
        x::Int
    end
    x = T1(1)

    do_thing(x::Vector) = sum(x.^2)
    do_thing(x::T1) = 1
end

@override M1

old_T1 = T1(1)

@assert x == T1(1)
@assert do_thing(x) == 1
@assert do_thing(data) == sum(data.^2)

# Make some changs to M1
module M1

    export T1, x, do_thing

    immutable T1
        x::Int
        y::Float64
    end
    global const x = T1(2,1.0)

    do_thing(x::Vector) = sin(x.^2)
    do_thing(x::T1) = 2
end

@override M1

@assert x == T1(2,1.0)
@assert do_thing(x) == 2
@assert do_thing(data) == sin(data.^2)

@test_throws MethodError do_thing(old_T1)


text:
# deterministic optimization models
function detBuild(Δt, T, fData, bData, dData, solveOpt = true)
    # deterministic formulation
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = round(fData.g[k]/(fData.g[k]^2 + fData.b[k]^2), digits = 6);
        Xdict[k] = -round(fData.b[k]/(fData.g[k]^2 + fData.b[k]^2), digits = 6);
    end

    # construct the first stage without disruption occurring
    mp = Model(solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,Threads = 1));

    # set up the variables
    @variable(mp, fData.Pmin[i] <= sp[i in fData.genIDList, t in 1:T] <= fData.Pmax[i]);
    @variable(mp, fData.Qmin[i] <= sq[i in fData.genIDList, t in 1:T] <= fData.Qmax[i]);
    sphatsum = Dict();
    for t in 1:T
        for i in fData.IDList
            sphatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sphatsum[i,t] += sp[j,t];
                end
            end
        end
    end
    sqhatsum = Dict();
    for t in 1:T
        for i in fData.IDList
            sqhatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sqhatsum[i,t] += sq[j,t];
                end
            end
        end
    end

    @variable(mp, p[k in fData.brList, t in 1:T]);
    @variable(mp, q[k in fData.brList, t in 1:T]);
    @variable(mp, fData.Vmin[i]^2 <= v[i in fData.IDList, t in 1:T] <= fData.Vmax[i]^2);
    @variable(mp, 0 <= w[i in bData.IDList, t in 0:T] <= bData.cap[i]);
    @variable(mp, y[i in bData.IDList, t in 1:T]);
    @variable(mp, zp[i in bData.IDList, t in 1:T]);
    @variable(mp, zq[i in bData.IDList, t in 1:T]);
    @variable(mp, lpp[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lqp[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lpm[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lqm[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, 0 <= u[i in bData.IDList] <= bData.uCap[i]);

    # set up the constraints
    @constraint(mp, pBalance[i in fData.IDList, t in 1:T], sum(zp[b,t] for b in bData.IDList if bData.Loc[b] == i) + lpp[i,t] - lpm[i,t] +
        sphatsum[i,t] - dData.pd[i][t] == sum(p[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, qBalance[i in fData.IDList, t in 1:T], sum(zq[b,t] for b in bData.IDList if bData.Loc[b] == i) + lqp[i,t] - lqm[i,t]  +
        sqhatsum[i,t] - dData.qd[i][t] == sum(q[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, pequal[k in fData.brList, t in 1:T], p[k,t] == -p[(k[2],k[1],k[3]),t]);
    @constraint(mp, qequal[k in fData.brList, t in 1:T], q[k,t] == -q[(k[2],k[1],k[3]),t]);
    @constraint(mp, lineThermal[k in fData.brList, t in 1:T], norm([p[k,t], q[k,t]]) <= fData.rateA[k]);
    @constraint(mp, powerflow[k in fData.brList, t in 1:T], v[k[2],t] == v[k[1],t] - 2*(Rdict[k]*p[k,t] + Xdict[k]*q[k,t]));
    @constraint(mp, rampUp[i in fData.genIDList, t in 2:T], sp[i,t] - sp[i,t - 1] <= fData.RU[i]);
    @constraint(mp, rampDown[i in fData.genIDList, t in 2:T], sp[i,t] - sp[i,t - 1] >= fData.RD[i]);
    @constraint(mp, bInv[i in bData.IDList, t in 1:T], w[i,t] == w[i,t-1] - y[i,t]*Δt);
    @constraint(mp, bThermal[i in bData.IDList, t in 1:T], norm([zp[i,t], zq[i,t]]) <= u[i]);
    @constraint(mp, bEfficient[i in bData.IDList, l in 1:length(bData.ηα[i]), t in 1:T], zp[i,t] <= bData.ηα[i][l]*y[i,t] + bData.ηβ[i][l]);
    @constraint(mp, bInvmax[i in bData.IDList, t in 1:T], w[i,t] <= bData.cap[i]);
    @constraint(mp, bInvIni[i in bData.IDList], w[i,0] == bData.bInv[i]);

    # set up the objective function
    @variable(mp,fs[i in fData.genIDList, t in 1:T] >= 0);
    @variable(mp,tAux1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3] >= 0);
    @variable(mp,tAux2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3]);
    @variable(mp,tAux3[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3]);

    @constraint(mp,gcAux1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux1[i,t] == fs[i,t] +
        (fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    @constraint(mp,gcAux2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux2[i,t] == fs[i,t] +
        (-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    @constraint(mp,gcAux3[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux3[i,t] == sqrt(fData.cp[i].params[1])*sp[i,t] +
        fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1])));
    @constraint(mp, genCost1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3],
        norm([tAux2[i,t],tAux3[i,t]]) <= tAux1[i,t]);
    @constraint(mp, genCost2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 2],
        fs[i,t] == fData.cp[i].params[1]*sp[i,t]);

    # deterministic objective function
    @objective(mp, Min, sum(bData.cost[i]*u[i] for i in bData.IDList) +
        sum(fData.cz*sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList) +
        sum(fs[i,t] for i in fData.genIDList) for t in 1:T));

    if solveOpt
        # solve the problem
        statusMp = solve(mp);
        mpObj = getobjectivevalue(mp);
        println("First stage, solving status $(statusMp)");

        # obtain the solutions
        solSp = Dict();
        solSq = Dict();
        solw = Dict();
        solu = Dict();
        solLp = Dict();
        solLq = Dict();
        solzp = Dict();
        for i in fData.genIDList
            for t in 1:T
                if abs(getvalue(sp[i,t])) > 1e-5
                    solSp[i,t] = getvalue(sp[i,t]);
                else
                    solSp[i,t] = 0;
                end
                if abs(getvalue(sq[i,t])) > 1e-5
                    solSq[i,t] = getvalue(sq[i,t]);
                else
                    solSq[i,t] = 0;
                end
            end
        end
        for i in bData.IDList
            if abs(getvalue(u[i])) > 1e-5
                solu[i] = getvalue(u[i]);
                for t in 1:T
                    if getvalue(w[i,t]) > 1e-5
                        solw[i,t] = getvalue(w[i,t]);
                    else
                        solw[i,t] = 0;
                    end
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            else
                solu[i] = 0;
                for t in 1:T
                    solw[i,t] = 0;
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            end
        end
        for i in fData.IDList
            for t in 1:T
                if (abs(getvalue(lpp[i,t])) > 1e-8)|(abs(getvalue(lpm[i,t])) > 1e-8)
                    solLp[i,t] = getvalue(lpp[i,t]) - getvalue(lpm[i,t]);
                else
                    solLp[i,t] = 0;
                end
                if (abs(getvalue(lqp[i,t])) > 1e-8)|(abs(getvalue(lqm[i,t])) > 1e-8)
                    solLq[i,t] = getvalue(lqp[i,t]) - getvalue(lqm[i,t]);
                else
                    solLq[i,t] = 0;
                end
            end
        end

        sol = solData(solSp,solSq,solw,solu,solLp,solLq,solzp);
        return sol,mpObj;
    else
        return mp;
    end
end

function fDetBuild(td, ωd, currentSol, τ, Δt, T, fData, bData, dData, solveOpt = true)
    # precalculate data
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = round(fData.g[k]/(fData.g[k]^2 + fData.b[k]^2), digits = 6);
        Xdict[k] = -round(fData.b[k]/(fData.g[k]^2 + fData.b[k]^2), digits = 6);
    end
    Bparams = Dict();
    for t in td:T
        # create B parameters
        for k in fData.brList
            # if the line is disrupted and it is within disruption time
            if (((k[1],k[2]) in ωd)|((k[2],k[1]) in ωd))&(t <= td + τ)
                if ((k[1],k[2]) in hardened) | ((k[2],k[1]) in hardened)
                    Bparams[k,t] = 1;
                else
                    Bparams[k,t] = 0;
                end
            else
                Bparams[k,t] = 1;
            end
        end
        for i in fData.genIDList
            if (i in ωd)&(t <= td + τ)
                if i in hardened
                    Bparams[i,t] = 1;
                else
                    Bparams[i,t] = 0;
                end
            else
                Bparams[i,t] = 1;
            end
        end
    end
    for i in fData.genIDList
        Bparams[i,td - 1] = 1;
    end

    mp = Model(solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,Threads = 1));

    # set up the variables
    @variable(mp, fData.Pmin[i]*Bparams[i,t] <= sp[i in fData.genIDList,t in (td - 1):T] <= fData.Pmax[i]*Bparams[i,t]);
    @variable(mp, fData.Qmin[i]*Bparams[i,t] <= sq[i in fData.genIDList,t in td:T] <= fData.Qmax[i]*Bparams[i,t]);
    sphatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sphatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sphatsum[i,t] += sp[j,t];
                end
            end
        end
    end
    sqhatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sqhatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sqhatsum[i,t] += sq[j,t];
                end
            end
        end
    end

    @variable(mp, p[k in fData.brList, t in td:T]);
    @variable(mp, q[k in fData.brList, t in td:T]);
    @variable(mp, fData.Vmin[i]^2 <= v[i in fData.IDList, t in td:T] <= fData.Vmax[i]^2);
    @variable(mp, 0 <= w[i in bData.IDList, t in (td - 1):T] <= bData.cap[i]);
    @variable(mp, y[i in bData.IDList, t in td:T]);
    @variable(mp, zp[i in bData.IDList, t in td:T]);
    @variable(mp, zq[i in bData.IDList, t in td:T]);
    @variable(mp, lpp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lpm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, 0 <= u[i in bData.IDList] <= bData.uCap[i]);

    # set up the constraints
    bigM = 1000;
    @constraint(mp, pBalance[i in fData.IDList, t in td:T], sum(zp[b,t] for b in bData.IDList if bData.Loc[b] == i) + lpp[i,t] - lpm[i,t] +
        sphatsum[i,t] - dData.pd[i][t] == sum(p[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, qBalance[i in fData.IDList, t in td:T], sum(zq[b,t] for b in bData.IDList if bData.Loc[b] == i) + lqp[i,t] - lqm[i,t] +
        sqhatsum[i,t] - dData.qd[i][t] == sum(q[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, pequal[k in fData.brList, t in td:T], p[k,t] == -p[(k[2],k[1],k[3]),t]);
    @constraint(mp, qequal[k in fData.brList, t in td:T], q[k,t] == -q[(k[2],k[1],k[3]),t]);
    @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], norm([p[k,t],q[k,t]]) <= fData.rateA[k]);
    @constraint(mp, lineThermal2[k in fData.brList, t in td:T;Bparams[k,t] == 0], p[k,t] == 0);
    @constraint(mp, lineThermal3[k in fData.brList, t in td:T;Bparams[k,t] == 0], q[k,t] == 0);
    @constraint(mp, powerflow1[k in fData.brList, t in td:T;Bparams[k,t] == 1], v[k[2],t] == v[k[1],t] - 2*(Rdict[k]*p[k,t] + Xdict[k]*q[k,t]));
    @constraint(mp, rampUp[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] <= fData.RU[i]);
    @constraint(mp, rampDown[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] >= fData.RD[i]);
    @constraint(mp, bInv[i in bData.IDList, t in td:T], w[i,t] == w[i,t-1] - y[i,t]*Δt);
    #@constraint(mp, bThermal[i in bData.IDList, t in td:T], norm([zp[i,t],zq[i,t]]) <= u[i]);
    @constraint(mp, bThermal1[i in bData.IDList, t in td:T], zp[i,t] <= u[i]);
    @constraint(mp, bThermal2[i in bData.IDList, t in td:T], zq[i,t] <= u[i]);
    @constraint(mp, bEfficient[i in bData.IDList, l in 1:length(bData.ηα[i]), t in td:T], zp[i,t] <= bData.ηα[i][l]*y[i,t] + bData.ηβ[i][l]);
    @constraint(mp, bInvmax[i in bData.IDList, t in td:T], w[i,t] <= bData.cap[i]);
    @constraint(mp, bInvIni[i in bData.IDList], w[i,td - 1] == currentSol.w[i,td - 1]);
    @constraint(mp, spIni[i in fData.genIDList], sp[i,td - 1] == currentSol.sp[i,td - 1]);
    @constraint(mp, uIni[i in bData.IDList], u[i] == currentSol.u[i]);

    # set up the objective function
    # set up the objective function
    @variable(mp,fs[i in fData.genIDList, t in td:T]);
    @variable(mp,tAux1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3]);
    @variable(mp,tAux2[i in fData.genIDList, t in td:T; fData.cp[i].n == 3]);
    @variable(mp,tAux3[i in fData.genIDList, t in td:T; fData.cp[i].n == 3]);

    @constraint(mp,gcAux1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux1[i,t] == fs[i,t] +
        (fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    @constraint(mp,gcAux2[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux2[i,t] == fs[i,t] +
        (-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
    @constraint(mp,gcAux3[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux3[i,t] == sqrt(fData.cp[i].params[1])*sp[i,t] +
        fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1])));
    @constraint(mp, genCost1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3],
        norm([tAux2[i,t],tAux3[i,t]]) <= tAux1[i,t]);
    @constraint(mp, genCost2[i in fData.genIDList, t in td:T; fData.cp[i].n == 2],
        fs[i,t] == fData.cp[i].params[1]*sp[i,t]);

    objExpr = @expression(mp, sum(fData.cz*sum(lpp[i,t] + lqp[i,t] + lpm[i,t] + lqm[i,t] for i in fData.IDList) +
        sum(fs[i,t] for i in fData.genIDList) for t in td:T));

    @objective(mp, Min, objExpr);

    if solveOpt
        # solve the problem
        statusMp = solve(mp);
        # optimize!(mp, with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0,
        #     QCPDual = 1, NumericFocus = 3, BarQCPConvTol = 1e-9, FeasibilityTol = 1e-9));
        println("Disruption time $(td), scenario $(ωd), solving status $(statusMp)");

        mpObj = getobjectivevalue(mp);
        # obtain the solutions
        solSp = Dict();
        solSq = Dict();
        solw = Dict();
        solu = Dict();
        solLp = Dict();
        solLq = Dict();
        solzp = Dict();
        for i in fData.genIDList
            for t in td:T
                if abs(getvalue(sp[i,t])) > 1e-5
                    solSp[i,t] = getvalue(sp[i,t]);
                else
                    solSp[i,t] = 0;
                end
                if abs(getvalue(sq[i,t])) > 1e-5
                    solSq[i,t] = getvalue(sq[i,t]);
                else
                    solSq[i,t] = 0;
                end
            end
        end
        for i in bData.IDList
            if abs(getvalue(u[i])) > 1e-5
                solu[i] = getvalue(u[i]);
                for t in td:T
                    if abs(getvalue(w[i,t])) > 1e-5
                        solw[i,t] = getvalue(w[i,t]);
                    else
                        solw[i,t] = 0;
                    end
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            else
                solu[i] = 0;
                for t in td:T
                    solw[i,t] = 0;
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            end
        end
        for i in fData.IDList
            for t in td:T
                if (abs(getvalue(lpp[i,t])) > 1e-8)|(abs(getvalue(lpm[i,t])) > 1e-8)
                    solLp[i,t] = getvalue(lpp[i,t]) - getvalue(lpm[i,t]);
                else
                    solLp[i,t] = 0;
                end
                if (abs(getvalue(lqp[i,t])) > 1e-8)|(abs(getvalue(lqm[i,t])) > 1e-8)
                    solLq[i,t] = getvalue(lqp[i,t]) - getvalue(lqm[i,t]);
                else
                    solLq[i,t] = 0;
                end
            end
        end

        sol = solData(solSp,solSq,solw,solu,solLp,solLq,solzp);
        return sol,mpObj;
    else
        return mp;
    end
end

function constructDetM(td, ωd, sol, τ, Δt, T, fData, bData, dData)
    # construct the math program given the state variables and current stage
    if td == 1
        # if it is the no-disruption problem
        sol,objV = detBuild(Δt, T, fData, bData, dData);
    else
        # if it is f_{ht}^ω
        sol,objV = fDetBuild(td, ωd, sol, τ, Δt, T, fData, bData, dData);
    end
    return sol,objV;
end

function buildPathDet(T, Δt, fData, bData, dData, pDistr, pathList = [])
    disT = 1;
    ωd = 0;
    costn = 0;
    iter = 1;
    solHist = [];
    τω = 0;
    currentSol = solData(Dict(),Dict(),Dict(),Dict(),Dict(),Dict(),Dict());
    while disT <= T
        # solve the current stage problem, state variables are passed
        nowT = disT;
        currentSol,objV = constructDetM(disT, ωd, currentSol, τω, Δt, T, fData, bData, dData);
        push!(solHist,(currentSol,nowT,ωd));

        # generate disruption
        if pathList == []
            tp,ωd,τω = genScenario(pDistr);
        else
            tp,ωd,τω = pathList[iter];
        end
        iter += 1;
        if nowT == 1
            disT += tp;
            disT = min(disT, T + 1);
            # calculate the cost of the solution until the next disruption time
            costn += sum(sum(fData.cz*(abs(currentSol.lp[i,t]) + abs(currentSol.lq[i,t])) for i in fData.IDList) for t in nowT:(disT - 1)) +
                sum(currentSol.u[i]*bData.cost[i] for i in bData.IDList);
            costn = calCostF(costn, currentSol, T, fData, nowT, disT);
        else
            disT += tp + τω;
            disT = min(disT, T + 1);
            # calculate the cost of the solution until the next disruption time
            costn += sum(sum(fData.cz*(abs(currentSol.lp[i,t]) + abs(currentSol.lq[i,t])) for i in fData.IDList) for t in nowT:(disT - 1));
            costn = calCostF(costn, currentSol, T, fData, nowT, disT);
        end
    end
    println("Det Path Built!");
    return [solHist,costn];
end

function exeDet(T, Δt, fData, bData, dData, pDistr, N, pathDict = Dict())
    # execution of forward pass
    # input: N: the number of trial points;
    # output: solList: a list of solution paths
    solDict = Dict();
    costDict = Dict();
    if pathDict == Dict()
        for i in 1:N
            pathDict[i] = [];
        end
    end
    returnData = pmap(i -> buildPathDet(T, Δt, fData, bData, dData, pDistr, pathDict[i]), 1:N);
    for n in 1:N
        solDict[n] = returnData[n][1];
        costDict[n] = returnData[n][2];
    end
    return solDict, costDict;
end


text:
mutable struct PardisoSolver <: AbstractPardisoSolver
    pt::Vector{Int}
    iparm::Vector{Int32}
    dparm::Vector{Float64}
    mtype::MatrixType
    solver::Solver
    phase::Phase
    msglvl::MessageLevel
    maxfct::Int32
    mnum::Int32
    perm::Vector{Int32}
end

function PardisoSolver()
    if !PARDISO_LOADED[]
      error("pardiso library was not loaded")
    end

    pt = zeros(Int, 64)
    iparm = zeros(Int32, 64)
    dparm = zeros(Float64, 64)
    mtype = REAL_NONSYM
    solver = DIRECT_SOLVER
    phase = ANALYSIS_NUM_FACT_SOLVE_REFINE
    msglvl = MESSAGE_LEVEL_OFF
    # Set numper of processors to CPU_CORES unless "OMP_NUM_THREADS" is set
    if haskey(ENV, "OMP_NUM_THREADS")
        iparm[3] = parse(Int, ENV["OMP_NUM_THREADS"])
    else
        # Assume 2 threads per core
        iparm[3] = max(div(Sys.CPU_THREADS, 2), 1)
    end

    mnum = 1
    maxfct = 1
    perm = Int32[]

    ps = PardisoSolver(pt, iparm, dparm, mtype, solver,
                  phase, msglvl, maxfct, mnum, perm)

    return ps
end

show(io::IO, ps::PardisoSolver) = print(io, string("$PardisoSolver:\n",
                                  "\tSolver: $(SOLVER_STRING[get_solver(ps)])\n",
                                  "\tMatrix type: $(MATRIX_STRING[get_matrixtype(ps)])\n",
                                  "\tPhase: $(PHASE_STRING[get_phase(ps)])\n",
                                  "\tNum processors: $(get_nprocs(ps))"))

phases(ps::PardisoSolver) = PHASES

set_transposed(ps::PardisoSolver, t::Bool) = t ? set_iparm(ps, 12, 1) : set_iparm(ps, 12, 0)

get_dparm(ps::PardisoSolver, i::Integer) = ps.dparm[i]
get_dparms(ps::PardisoSolver) = ps.dparm
set_dparm!(ps::PardisoSolver, i::Integer, v::AbstractFloat) = ps.dparm[i] = v
get_nprocs(ps::PardisoSolver) = ps.iparm[3]

set_solver!(ps::PardisoSolver, v::Int) = set_solver!(ps, Solver(v))
function set_solver!(ps::PardisoSolver, v::Solver)
    ps.solver = v
end
get_solver(ps::PardisoSolver) = ps.solver

@inline function ccall_pardisoinit(ps::PardisoSolver)
    ERR = Ref{Int32}(0)
    ccall(init[], Cvoid,
          (Ptr{Int}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Float64}, Ptr{Int32}),
          ps.pt, Ref(Int32(ps.mtype)), Ref(Int32(ps.solver)), ps.iparm, ps.dparm, ERR)
     check_error(ps, ERR[])
end

@inline function ccall_pardiso(ps::PardisoSolver, N::Int32, AA::Vector{Tv},
                                   IA, JA, NRHS::Int32, B::StridedVecOrMat{Tv}, X::StridedVecOrMat{Tv}) where {Tv}
    ERR = Ref{Int32}(0)
    ccall(pardiso_f[], Cvoid,
          (Ptr{Int}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Tv}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Tv},
           Ptr{Int32}, Ptr{Float64}),
          ps.pt, Ref(ps.maxfct), Ref(Int32(ps.mnum)), Ref(Int32(ps.mtype)), Ref(Int32(ps.phase)),
          Ref(N), AA, IA, JA, ps.perm,
          Ref(NRHS), ps.iparm, Ref(Int32(ps.msglvl)), B, X,
          ERR, ps.dparm)
    check_error(ps, ERR[])
end

@inline function ccall_pardiso_get_schur(ps::PardisoSolver, S::Vector{Tv},
                                   IS::Vector{Int32}, JS::Vector{Int32}) where Tv
    ccall(pardiso_get_schur_f[], Cvoid,
          (Ptr{Int}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Tv},
           Ptr{Int32}, Ptr{Int32}),
          ps.pt, Ref(ps.maxfct), Ref(Int32(ps.mnum)), Ref(Int32(ps.mtype)),S,
          IS, JS)
end

# Different checks
function printstats(ps::PardisoSolver, A::SparseMatrixCSC{Tv, Ti},
                    B::StridedVecOrMat{Tv}) where {Ti,Tv <: PardisoNumTypes}
    N = Int32(size(A, 2))
    AA = A.nzval
    IA = convert(Vector{Int32}, A.colptr)
    JA = convert(Vector{Int32}, A.rowval)
    NRHS = Int32(size(B, 2))
    ERR = Ref{Int32}(0)
    if Tv <: Complex
        f = pardiso_printstats_z[]
    else
        f = pardiso_printstats[]
    end
    ccall(f, Cvoid,
          (Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Int32},
           Ptr{Int32}, Ptr{Int32}, Ptr{Tv},
           Ptr{Int32}),
          Ref(Int32(ps.mtype)), Ref(N), AA, IA, JA, Ref(NRHS), B, ERR)

    check_error(ps, ERR[])
    return
end

function checkmatrix(ps::PardisoSolver, A::SparseMatrixCSC{Tv, Ti}) where {Ti,Tv <: PardisoNumTypes}
    N = Int32(size(A, 1))
    AA = A.nzval
    IA = convert(Vector{Int32}, A.colptr)
    JA = convert(Vector{Int32}, A.rowval)
    ERR = Ref{Int32}(0)

    if Tv <: Complex
        f = pardiso_chkmatrix_z[]
    else
        f = pardiso_chkmatrix[]
    end

    ccall(f, Cvoid,
          (Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Int32},
           Ptr{Int32}, Ptr{Int32}),
          Ref(Int32(ps.mtype)), Ref(N), AA, IA,
          JA, ERR)

    check_error(ps, ERR[])
    return
end

function checkvec(ps, B::StridedVecOrMat{Tv}) where {Tv <: PardisoNumTypes}
    N = Int32(size(B, 1))
    NRHS = Int32(size(B, 2))
    ERR = Int32[0]

    if Tv <: Complex
        f = pardiso_chkvec_z[]
    else
        f = pardiso_chkvec[]
    end
    ccall(f, Cvoid,
          (Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Int32}),
          Ref(N), Ref(NRHS), B, ERR)

    check_error(ps, ERR[])
    return
end

function check_error(ps::PardisoSolver, err::Integer)
    err != -1  || throw(PardisoException("Input inconsistent."))
    err != -2  || throw(PardisoException("Not enough memory."))
    err != -3  || throw(PardisoException("Reordering problem."))
    err != -4  || throw(PardisoPosDefException("Zero pivot, numerical fact. or iterative refinement problem."))
    err != -5  || throw(PardisoException("Unclassified (internal) error."))
    err != -6  || throw(PardisoException("Preordering failed (matrix types 11, 13 only)."))
    err != -7  || throw(PardisoException("Diagonal matrix problem."))
    err != -8  || throw(PardisoException("32-bit integer overflow problem."))
    err != -10 || throw(PardisoException("No license file pardiso.lic found."))
    err != -11 || throw(PardisoException("License is expired."))
    err != -12 || throw(PardisoException("Wrong username or hostname."))
    err != -100|| throw(PardisoException("Reached maximum number of Krylov-subspace iteration in iterative solver."))
    err != -101|| throw(PardisoException("No sufficient convergence in Krylov-subspace iteration within 25 iterations."))
    err != -102|| throw(PardisoException("Error in Krylov-subspace iteration."))
    err != -103|| throw(PardisoException("Break-Down in Krylov-subspace iteration."))
    return
end


text:
#######################################################################################################################################################################################################
#
# Changes to this type
# General
#     2022-Feb-01: add abstract type for vulnerability curve
#     2022-Feb-01: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Hierachy of AbstractVulnreabilityCurve:
- [`LogisticVC`](@ref)
- [`PowerVC`](@ref)
- [`WeibullVC`](@ref)
- [`ComplexVC`](@ref)
"""
abstract type AbstractVulnreabilityCurve{FT<:AbstractFloat} end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Feb-01: add LogisticVC (use mutable so as to fit)
#     2022-Feb-01: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Modified logistic function for vulnerability curve

# Fields

$(TYPEDFIELDS)

"""
mutable struct LogisticVC{FT<:AbstractFloat} <: AbstractVulnreabilityCurve{FT}
    "Multiplier to exponential component"
    A::FT
    "Multiplier to pressure `[MPa⁻¹]`"
    B::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Feb-01: add PowerVC (use mutable so as to fit)
#     2022-Feb-01: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Power function for vulnerability curve

# Fields

$(TYPEDFIELDS)

"""
mutable struct PowerVC{FT<:AbstractFloat} <: AbstractVulnreabilityCurve{FT}
    "Multiplier to power component `[MPa⁻ᵇ]`"
    A::FT
    "Power to pressure"
    B::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Feb-01: add WeibullVC (use mutable so as to fit)
#     2022-Feb-01: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Weibull cumulative distribution function for vulnerability curve

# Fields

$(TYPEDFIELDS)

"""
mutable struct WeibullVC{FT<:AbstractFloat} <: AbstractVulnreabilityCurve{FT}
    "Numerator in the exponential component `[MPa]`"
    B::FT
    "Power to pressure component"
    C::FT
end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Feb-01: add ComplexVC (use mutable so as to fit)
#     2022-Feb-01: add documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

A complex struct for segmented vulnerability curve such as dual Weibull function

# Fields

$(TYPEDFIELDS)

"""
mutable struct ComplexVC{FT} <: AbstractVulnreabilityCurve{FT}
    "Percentages of each VC component"
    PS::Vector{FT}
    "Vector of vulnerability curve components"
    VCS::Union{Vector{LogisticVC{FT}}, Vector{PowerVC{FT}}, Vector{WeibullVC{FT}}, Vector{AbstractVulnreabilityCurve{FT}}}
end


text:
module Stochy

using Base.Collections
import Base: factor, show

import Distributions
import Distributions: Distribution, support, rand, params
const score = Distributions.logpdf

export @pp, sample, factor, score, observe, mem, cache
export Discrete, Dir, flip, randominteger, hellingerdistance, kl

const primitives = [:!,:+,:*,:-,:/,:sqrt,:√,
                    :tuple, :mem, :println, :cons, :list, :tail, :cat,
                    :reverse, :.., :first, :second, :third, :fourth]

debugflag = false
debug(b::Bool) = global debugflag = b
debug() = debugflag

partial(f,arg) = (args...) -> f(arg,args...)

macro pp(expr)
    tx = storetransform(cps(addressing(desugar(expr)), :(Stochy.store_exit)))
    debug() && println(striplineinfo(tx))
    esc(tx)
end

abstract Ctx
type Prior <: Ctx end
ctx = Prior()

include("cps.jl")
include("store.jl")
include("addressing.jl")
include("erp.jl")
include("rand.jl")
include("enumerate.jl")
include("mh.jl")
include("pmcmc.jl")
include("variational.jl")
include("dp.jl")
include("plotting/gadfly.jl")
include("plotting/pyplot.jl")

# Dispatch based on current context.
sample(s::Store, k::Function, a, e::ERP) = sample(s,k,a,e,ctx)
factor(s::Store, k::Function, a, score) = factor(s,k,a,score,ctx)

sample(s::Store, k::Function, a, e::ERP, ::Prior) = k(s,rand(e))

# @pp
score(s::Store, k::Function, a, e::ERP, x) = k(s,score(e,x))

observe(s::Store, k::Function, a, erp::ERP, x) = factor(s, k, a, score(erp,x))
observe(s::Store, k::Function, a, erp::ERP, xs...) = factor(s, k, a, sum([score(erp,x) for x in xs]))

# cache and mem are similar. The primary difference is that mem uses a
# cache local to the "thread", cache uses a global cache. Another
# difference is that the cache used by mem is local to a @pp block,
# whereas the cache shares its cache between @pp blocks.

# @pp
function mem(f::Function)
    key = gensym()
    (store::Store, k::Function, address, args...) -> begin
        cachekey = (key,args)
        if haskey(store,cachekey)
            k(store,store[cachekey])
        else
            partial(f,store)(address, args...) do s,val
                s = copy(s)
                s[cachekey] = val
                k(s,val)
            end
        end
    end
end

const CACHE = Dict()

# @pp
function cache(s::Store, k::Function, address, f::Function)
    partial(k,s)() do s2, k2, a2, args...
        # This is the body of the memoized function.
        if haskey(CACHE, args)
            k2(s2, CACHE[args])
        else
            partial(f,s2)(a2, args...) do s3, value
                CACHE[args] = value
                k2(s3,value)
            end
        end
    end
end

import Base.==, Base.hash, Base.first, Base.isless
export .., first, second, third, fourth, repeat

using DataStructures: Cons, Nil, head, tail, cons, list, cat
export cons, list

const .. = cons

==(x::Cons,y::Cons) = head(x) == head(y) && tail(x) == tail(y)
==(x::Nil,y::Cons) = false
==(x::Cons,y::Nil) = false
==(x::Nil,y::Nil) = true

hash(x::Cons,h::Uint64) = hash(tail(x), hash(head(x), h))
hash(_::Nil,h::Uint64) = hash(object_id(list()), h)

function isless(xs::Cons,ys::Cons)
    x, y = head(xs), head(ys)
    x==y ? isless(tail(xs), tail(ys)) : isless(x,y)
end

isless(xs::Nil,ys::Cons) = true
isless(xs::Cons,ys::Nil) = false
isless(xs::Nil,ys::Nil) = true

first(l::Cons)  = head(l)
second(l::Cons) = head(tail(l))
third(l::Cons)  = head(tail(tail(l)))
fourth(l::Cons) = head(tail(tail(tail(l))))

# TODO: Change the base case to n==1 so that the created list is
# tightly typed. Throw error for n<1.
@pp function repeat(f::Function, n::Int64)
    n < 1 ? list() : cons(f(), repeat(f, n-1))
end

end