text:
using Revise # reduce recompilation time
using Plots
using LinearAlgebra
using StaticArrays
using SparseArrays
using BenchmarkTools
using UnPack
using ForwardDiff

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh
using SetupDG
using ExplicitJacobians

push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

"Approximation parameters"
N = 1 # The order of approximation
K1D = 8
CFL = 10
T = 1 # endtime

"Mesh related variables"
VX, VY, EToV = uniform_tri_mesh(K1D)
# VX = @. VX - .3*sin(pi*VX)

# initialize ref element and mesh
rd = init_reference_tri(N)
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,FToF,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps!(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB,FToF)
mapP[mapB] = mapPB

## construct hybridized SBP operators
@unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd
Qr = Pq'*M*Dr*Pq
Qs = Pq'*M*Ds*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Bs = diagm(wf.*nsJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
        -Br*Ef  Br]
Qsh = .5*[Qs-Qs' Ef'*Bs;
        -Bs*Ef  Bs]

Vh = [Vq;Vf]
Ph = M\transpose(Vh)
VhP = Vh*Pq

# make sparse skew symmetric versions of the operators"
Qrhskew = .5*(Qrh-transpose(Qrh))
Qshskew = .5*(Qsh-transpose(Qsh))

# interpolate geofacs to both vol/surf nodes
@unpack rxJ, sxJ, ryJ, syJ = md
rxJ, sxJ, ryJ, syJ = (x->Vh*x).((rxJ, sxJ, ryJ, syJ)) # interp to hybridized points

## global matrices

Ax,Ay,Bx,By,B = assemble_global_SBP_matrices_2D(rd,md,Qrhskew,Qshskew)

# add off-diagonal couplings
Ax += Bx
Ay += By

Ax *= 2 # for flux differencing
Ay *= 2

AxTr = sparse(transpose(Ax))
AyTr = sparse(transpose(Ay))
Bx   = abs.(Bx) # create LF penalization term

# globalize operators and nodes
@unpack x,y,J = md
Vq   = kron(speye(K),sparse(Vq))
Pq   = kron(speye(K),sparse(Pq))
VhTr = kron(speye(K),sparse(transpose(Vh)))
M    = kron(spdiagm(0 => J[1,:]),sparse(M))
Vh   = kron(speye(K),sparse(Vh))
Ph   = kron(spdiagm(0 => 1 ./ J[1,:]), sparse(Ph))
VhP  = Vh*Pq
x,y = (a->a[:]).((x,y))
wJq = vec(diagm(rd.wq)*(rd.Vq*J))

println("Done building global ops")

## define Euler fluxes
function F(UL,UR)
    # convert to flux variables
    function UtoQ(U)
        rho,rhou,rhov,E = U
        beta = betafun(U...)
        return (rho,rhou./rho,rhov./rho,beta),(log.(rho),log.(beta))
    end
    QL,QlogL = UtoQ(UL)
    QR,QlogR = UtoQ(UR)
    Fx,Fy = euler_fluxes(QL...,QR...,QlogL...,QlogR...)
    return SVector{length(Fx)}(Fx...),SVector{length(Fy)}(Fy...)
end

# extract coordinate fluxes
Fx = (uL,uR)->F(uL,uR)[1]
Fy = (uL,uR)->F(uL,uR)[2]

# AD for jacobians
dFx(uL,uR) = ForwardDiff.jacobian(uR->F(uL,uR)[1],uR)
dFy(uL,uR) = ForwardDiff.jacobian(uR->F(uL,uR)[2],uR)

## dissipative flux

function LF(uL,uR,nxL,nyL,nxR,nyR)
        rhoL,rhouL,rhovL,EL = uL
        rhoR,rhouR,rhovR,ER = uR

        nx = nxL #.5*(nxL+nxR)
        ny = nyL #.5*(nyL+nyR)
        rhoUnL = @. rhouL*nx + rhovL*ny
        rhoUnR = @. rhouR*nx + rhovR*ny
        cL   = @. wavespeed(rhoL,rhoUnL,EL)
        cR   = @. wavespeed(rhoR,rhoUnR,ER)
        lam  = @. max(abs(cL),abs(cR))
        return (@. lam*(uL-uR))
end
dLF(uL,uR,args...) = ForwardDiff.jacobian(uR->LF(uL,uR,args...),uR)

## mappings between conservative and entropy variables and vice vera
dVdU_fun(U) = ForwardDiff.jacobian(U->SVector(v_ufun(U...)...),U)
dUdV_fun(V) = ForwardDiff.jacobian(V->SVector(u_vfun(V...)...),V)

## nonlinear solver setup - uses closure to initialize arrays and other things
function init_newton_fxn(Q,ops,rd::RefElemData,md::MeshData)

        # set up normals for use in penalty term
        @unpack Vq,Vf = rd
        @unpack nxJ,nyJ,sJ = md
        Nq,Np = size(Vq)
        Nf    = size(Vf,1)
        Nh    = Nq + Nf
        fids  = Nq + 1:Nh
        nxh,nyh = ntuple(x->zeros(Nh,K),2)
        nxh[fids,:] = @. nxJ/sJ
        nyh[fids,:] = @. nyJ/sJ
        nxh,nyh = vec.((nxh,nyh))

        # get lengths of arrays
        Ax,Ay,AxTr,AyTr,Bx,Vh,Vq,Pq = ops

        Nfields = length(Q)
        Id_fld  = speye(Nfields) # for Kronecker expansion to large matrices - fix later with lazy evals
        Vq_fld  = droptol!(kron(Id_fld,Vq),1e-12)
        VhP     = Vh*Pq
        VhP_fld = droptol!(kron(Id_fld,VhP),1e-12)
        Vh_fld  = droptol!(kron(Id_fld,Vh),1e-12)
        M_fld   = droptol!(kron(Id_fld,M),1e-12)
        Ph_fld  = droptol!(kron(Id_fld,Ph),1e-12)

        # init jacobian matrix (no need for entropy projection since we'll zero it out later)
        dFdU_h = repeat(I+Ax+Ay,Nfields,Nfields)

        function midpt_newton_iter!(Qnew, Qprev) # for Burgers' eqn specifically

                # perform entropy projection
                Uq    = (x->Vq*x).(SVector{Nfields}(Qnew))
                VUh   = SVector{Nfields}((x->VhP*x).(v_ufun(Uq...)))
                Qh    = SVector{Nfields}(u_vfun(VUh...))

                ftmp  = hadamard_sum(AxTr,Fx,Qh) + hadamard_sum(AyTr,Fy,Qh) + hadamard_sum(B,LF,Qh,nxh,nyh)
                f     = Ph_fld*vcat(ftmp...)
                res   = vcat(Qnew...) + .5*dt*f - vcat(Qprev...)

                fill!(dFdU_h.nzval,0.0)
                accum_hadamard_jacobian!(dFdU_h, Ax, dFx, Qh)
                accum_hadamard_jacobian!(dFdU_h, Ay, dFy, Qh)
                accum_hadamard_jacobian!(dFdU_h, B,  dLF, Qh, nxh,nyh) # flux term involving normals
                dVdU_h = banded_matrix_function(dVdU_fun, Uq)
                dUdV_h = banded_matrix_function(dUdV_fun, VUh)
                dFdU   = droptol!(transpose(Vh_fld)*(dFdU_h*dUdV_h*VhP_fld*dVdU_h*Vq_fld),1e-12)

                # solve and update
                dQ   = (M_fld + .5*dt*dFdU)\(M_fld*res)
                # TODO: compute damping factor
                Qnew = vcat(Qnew...) - dQ                            # convert Qnew to column vector for update
                Qnew = columnize(reshape(Qnew,length(Q[1]),Nfields)) # convert back to array of arrays

                return Qnew,norm(dQ)
        end
        return midpt_newton_iter!
end

# pack inputs together
ops = (Ax,Ay,copy(transpose(Ax)),copy(transpose(Ay)),Bx,Vh,Vq,Pq)

## init condition, rhs

rho = @. 2 + (abs(x)<.5)
rho = @. 2 + exp(-10*(x^2+y^2)) #ones(size(x))
rhou = zeros(size(x))
rhov = zeros(size(x))
E = rho.^1.4
Q = @SVector [rho,rhou,rhov,E]
# rho,u,v,p = vortex(x,y,0)
# Q = SVector{4}(primitive_to_conservative(rho,u,v,p))

# set time-stepping constants
CN = (N+1)*(N+2)/2  # estimated trace constant
h = minimum(J)
dt = CFL * 2 * h / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

# initialize jacobian
midpt_newton_iter! = init_newton_fxn(Q,ops,rd,md)

## newton time iteration

it_count = zeros(Nsteps)
energy   = zeros(Nsteps)
for i = 1:Nsteps
        global Q

        Qnew = copy(Q)  # copy / over-write at each timestep
        iter = 0
        dQnorm = 1
        while dQnorm > 1e-12
                Qnew,dQnorm = midpt_newton_iter!(Qnew,Q)
                iter += 1
                if iter > 15
                        println("iter = $iter")
                end
        end
        it_count[i] = iter
        Q = @. 2*Qnew - Q # implicit midpoint rule

        Qq = (x->Vq*x).(Q)
        energy[i] = sum(wJq.*Sfun(Qq...))

        if i%5==0 || i==Nsteps
                println("Number of time steps $i out of $Nsteps")
                # display(scatter(x,Q[1]))
        end
end

@unpack VDM = rd
rp, sp = Basis2DTri.equi_nodes_2D(25)
Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM

gr(aspect_ratio=1, legend=false,
   markerstrokewidth=0, markersize=2)
xp,yp,up = (x->Vp*reshape(x,size(Vp,2),K)).((x,y,Q[1]))
# display(scatter(xp,yp,up,zcolor=up,cam=(3,25),axis=false))
scatter(xp,yp,up,zcolor=up,cam=(0,90),border=:none,axis=false)
# png("sol_unif_mesh.png")

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
import PerronFrobenius

pts = rand(3, 10000)
ϵF = [0.4, 0.4, 0.4]
ϵs = [ϵF, ϵF]

induced_measures = [inducedrectangularinvariantmeasure(pts, ϵF, ϵᵢ) for ϵᵢ in [ϵF, ϵF]]
avg_measure = average_measure(induced_measures)

measure_ϵF = invariantmeasure(pts, RectangularBinning(ϵF))

@testset "Average measure constructor" begin
    ϵF = [0.4, 0.4, 0.4]
    ϵs = [ϵF, ϵF]
    pts = rand(3, 5000)

    ϵF = RectangularBinning(ϵF)
    ϵs = [RectangularBinning(ϵi) for (i, ϵi) in enumerate(ϵs)]

    @testset "Computing average measure directly from points" begin
        # Construct a AverageRectangularInvariantMeasure directly from data, final
        # resolution and resolutions from which we induce the measure.
        a2 = averagerectangularinvariantmeasure(pts, ϵF, ϵs)
        @test a2 isa AverageRectangularInvariantMeasure
    end
end

# PerronFrobenius.AverageRectangularInvariantMeasure(pts, ϵF,ϵs)

ext: jl


text:
function test_filt_base!(f::Function, order, delay, args...)
    srand(98765)
    dummy_input = rand(1024)
    dummy_sp = abs(fft(dummy_input))
    dummy_logsp = log(dummy_sp)
    dummy_mgceps = rand(order + 1)

    @assert applicable(f, first(dummy_input), dummy_mgceps, args..., delay)
    for x in dummy_input
        @test isfinite(f(x, dummy_mgceps, args..., delay))
    end
end

println("test waveform generation filters")

println("-- test_poledf")
for order in [20,  25, 30]
    println(" where order = $order")
    delay = poledf_delay(order)
    test_filt_base!(poledf, order, delay)
end

let
    @test_throws DimensionMismatch poledf(0.0, ones(10), ones(1))
end

println("-- test_lmadf")
for order in [20,  25, 30]
    for pd in 4:5
        println(" where order = $order, pd = $pd")
        delay = lmadf_delay(order, pd)
        test_filt_base!(lmadf, order, delay, pd)
    end
end

let
    # invalid delay length
    @test_throws DimensionMismatch lmadf(0.0, ones(10), 5, ones(1))
    # invalid pade
    @test_throws ArgumentError lmadf(0.0, ones(10), 3, lmadf_delay(9, 3))
end

println("-- test_lspdf")
for order in [20, 21, 25, 26, 30, 31]
    println(" where order = $order")
    delay = lspdf_delay(order)
    test_filt_base!(lspdf, order, delay)
end

let
    @test_throws DimensionMismatch lspdf(0.0, ones(10), ones(1))
    @test_throws DimensionMismatch lspdf(0.0, ones(9), ones(1))
end

println("-- test_ltcdf")
for order in [20, 21, 25, 26, 30, 31]
    println(" where order = $order")
    delay = ltcdf_delay(order)
    test_filt_base!(ltcdf, order, delay)
end

let
    @test_throws DimensionMismatch ltcdf(0.0, ones(10), ones(1))
end

println("-- test_glsadf")
for order in [20, 21, 25, 26, 30, 31]
    for stage in 1:10
        println(" where order = $order, stage = $stage")
        delay = glsadf_delay(order, stage)
        test_filt_base!(glsadf, order, delay, stage)
    end
end

let
    # invalid delay length
    @test_throws DimensionMismatch glsadf(0.0, ones(10), 1, ones(1))
    # invalid stage
    @test_throws ArgumentError glsadf(0.0, ones(10), 0, glsadf_delay(9, 0))
end

println("-- test_mlsadf")
for order in [20,  25, 30]
    for α in [0.0, 0.35, 0.41, 0.5]
        for pd in 4:5
            println(" where order = $order, α = $α, pd = $pd")
            delay = mlsadf_delay(order, pd)
            test_filt_base!(mlsadf, order, delay, α, pd)
        end
    end
end

let
    # invalid delay length
    @test_throws DimensionMismatch mlsadf(0.0, ones(10), 0.41, 5, ones(1))
    # invalid pade
    @test_throws ArgumentError mlsadf(0.0, ones(10), 0.41, 3, mlsadf_delay(9, 3))
end

println("-- test_mglsadf")
for order in [20,  25, 30]
    for α in [0.0, 0.35, 0.41, 0.5]
        for stage in 1:10
            println(" where order = $order, α = $α, stage = $stage")
            delay = mglsadf_delay(order, stage)
            test_filt_base!(mglsadf, order, delay, α, stage)
        end
    end
end

let
    # invalid delay length
    @test_throws DimensionMismatch mglsadf(0.0, ones(10), 0.41, 1, ones(1))
    # invalid stage
    @test_throws ArgumentError mglsadf(0.0, ones(10), 0.41, 0, glsadf_delay(9, 0))
end

let
    @test applicable(SPTK.poledf_delay_length, 10)
    @test applicable(SPTK.lmadf_delay_length, 10, 5)
    @test applicable(SPTK.lspdf_delay_length, 10)
    @test applicable(SPTK.ltcdf_delay_length, 10)
    @test applicable(SPTK.glsadf_delay_length, 10, 1)
    @test applicable(SPTK.mlsadf_delay_length, 10, 5)
    @test applicable(SPTK.mglsadf_delay_length, 10, 1)
end


text:
function mis2str(x::Missing) return "NA" end
function mis2str(x::Char) return x end
function mis2str(x) return string(x) end

function jags(
  model::Jagsmodel,
  data = Dict{String, Any}(),
  init = Dict{String, Any}(),
  ProjDir=pwd();
  updatedatafile::Bool=true,
  updateinitfiles::Bool=true
  )

  old = pwd()

  try
    cd(ProjDir)

    if updatedatafile
      if length(keys(data)) > 0
        print("\nCreating data file $(model.name)-data.R - ")
        @time update_R_file(joinpath(model.tmpdir, "$(model.name)-data.R"), data)
      end
    else
      println("\nData file not updated.")
    end

    if updateinitfiles
      if size(init, 1) > 0
        update_init_files(model, init)
      end
    else
      println("\nInit files not updated.")
    end

    println()
    curdir = pwd()
    cd(model.tmpdir)
    println("Executing $(model.ncommands) command(s), each with $(model.nchains) chain(s) took:")
    @time run(pipeline(par(model.command), stdout="$(model.name)-run.log"))
    sim = mchain(model)
    cd(old)
    return(sim)
  catch e
    println(e)
    cd(old)
  end
end

#### Update data and init files

function update_init_files(model::Jagsmodel, init)
  println()
  k = length(init)
  m = max(model.nchains, model.ncommands)
  indx = filter(x -> x!=0, [%(i, (k+1)) for i in 1:2m])
  for i in 1:m
    print("Creating init file $(model.name)-inits$(i).R - ")
    @time update_R_file(joinpath(model.tmpdir, "$(model.name)-inits$(i).R"), init[indx[i]])
  end
end

function update_R_file(file::String, dct; replaceNaNs::Bool=true)
  isfile(file) && rm(file)
  strmout = open(file, "w")

  str = ""
  for entry in dct

    str = "\""*entry[1]*"\""*" <- "
    val = entry[2]
    if replaceNaNs
      if typeof(entry[2]) == Array{Float64, 1}
        if true in isnan.(entry[2])
          val = convert(Array{Union{Float64,Missing}},entry[2])
          for i in 1:length(val)
            if isnan(val[i])
              val[i] = NA
            end
          end
        end
      end
      if typeof(entry[2]) == Array{Float64, 2}
        if true in isnan.(entry[2])
          val = convert(Array{Union{Float64,Missing}},entry[2])
          k,l = size(val)
          for i in 1:k
            for j in 1:l
              if isnan(val[i, j])
                val[i, j] = missing
              end
            end
          end
        end
      end
     end

    if typeof(val) <: String
      str = str*"\"$(val)\"\n"
    elseif length(val)==1 && length(size(val))==0
      # Scalar
      str = str * mis2str(val) * "\n"
    elseif length(val)>=1 && length(size(val))==1
      # Vector
      str = str*"structure(c("
      for i in 1:length(val)
        str = str*"$(mis2str(val[i]))"
        if i < length(val)
          str = str*", "
        end
      end
      str = str*"), .Dim=c($(length(val))))\n"
    elseif length(val)>1 && length(size(val))>1
      # Array
      str = str*"structure(c("
      for i in 1:length(val)
        str = str*"$(mis2str(val[i]))"
        if i < length(val)
          str = str*", "
        end
      end
      dimstr = "c"*string(size(val))
      str = str*"), .Dim=$(dimstr))\n"
    end
    write(strmout, str)
  end
  close(strmout)
end

#### use readdlm to read in all chains and create a Dict

function read_jagsfiles(model::Jagsmodel)
  index = readdlm(joinpath(model.tmpdir, "$(model.name)-cmd1-index.txt"), header=false)
  idxdct = Dict{String, Any}()
  for row in 1:size(index)[1]
    idxdct[index[row, 1]]=[Int(index[row, 2]), Int(index[row, 3])]
  end

  ## Collect the results of a chain in an array ##

  chainarray = Dict{String, Any}[]

  ## Each chain dictionary can contain up to 4 types of results ##

  result_type_files = ["samples"]
  rtdict = Dict{String, Any}()
  res_type = result_type_files[1]

  ## tdict contains the arrays of values ##
  tdict = Dict{String, Any}()

  println()
  for i in 1:model.ncommands
    tdict = Dict{String, Any}()
    for j in 1:model.nchains
      if isfile(joinpath(model.tmpdir, "$(model.name)-cmd$(i)-chain$(j).txt"))
        println("Reading $(model.name)-cmd$(i)-chain$(j).txt")
        res = readdlmcsv(joinpath(model.tmpdir, "$(model.name)-cmd$(i)-chain$(j).txt"))
        for key in index[:, 1]
          indx1 = idxdct[key][1]
          indx2 = idxdct[key][2]
          if length(keys(tdict)) == 0
            tdict = Dict(key => res[indx1:indx2, 2])
          else
            tdict = merge(tdict, Dict(key => res[indx1:indx2, 2]))
          end
        end
        ## End of processing result type file ##
        ## If any keys were found, merge it in the rtdict ##

        if length(keys(tdict)) > 0
          #println("Merging $(convert(Symbol, res_type)) with keys $(keys(tdict))")
          rtdict[res_type]=tdict
          tdict = Dict{String, Any}()
        end
      end

      ## If rtdict has keys, push it to the chain array ##

      if length(keys(rtdict)) > 0
        #println("Pushing the rtdict with keys $(keys(rtdict))")
        push!(chainarray, rtdict)
        rtdict = Dict{String, Any}()
      end
    end
  end
  (idxdct, chainarray)
end

#### Create a MCMCChains::Chains result

function mchain(model::Jagsmodel)
  println()
  local totalnchains, curchain
  index = readdlm(joinpath(model.tmpdir, "$(model.name)-cmd1-index.txt"),
    header=false)

  # Correct model.adapt for jagsthin != 1
  # if jagsthin != 1, adaptation samples are not included.

  if model.jagsthin != 1
    model.adapt = 1
  end

  cnames = String[]
  for i in 1:size(index)[1]
    append!(cnames, [index[i]])
  end

  totalnchains = model.nchains * model.ncommands
  a3d = fill(0.0, Int(index[1, 3]),
    size(index, 1), totalnchains);
  for i in 1:model.ncommands
    for j in 1:model.nchains
      if isfile(joinpath(model.tmpdir, "$(model.name)-cmd$(i)-chain$(j).txt"))
        println("Reading $(model.name)-cmd$(i)-chain$(j).txt")
        res = readdlmcsv(joinpath(model.tmpdir, "$(model.name)-cmd$(i)-chain$(j).txt"))
        curchain = (i-1)*model.nchains + j
        #println(curchain)
        k = 0
        for key in cnames
          k += 1
          a3d[:, k, curchain] = res[index[k, 2]:index[k, 3], 2]
        end
      end
    end
  end
  println()
  MCMCChains.Chains(a3d,cnames)
end

#### Read DIC related results

function read_pDfile(model::Jagsmodel)
  index = readdlm(joinpath(model.tmpdir, "$(model.name)-cmd1-index0.txt"), header=false);
  idxdct = Dict{String, Any}()
  for row in 1:size(index)[1]
    idxdct[index[row, 1]]=[Int(index[row, 2]), Int(index[row, 3])]
  end

  ## Collect the results of a chain in an array ##

  chainarray = Dict{String, Any}[]

  ## Each chain dictionary can contain up to 4 types of results ##

  result_type_files = ["samples"]
  rtdict = Dict{String, Any}()
  res_type = result_type_files[1]

  ## tdict contains the arrays of values ##
  tdict = Dict{String, Any}()

  for i in 0:0
    tdict = Dict{String, Any}()
    if isfile(joinpath(model.tmpdir, "$(model.name)-cmd1-chain$(i).txt"))
      println("Reading $(model.name)-cmd1-chain$(i).txt")
      res = readdlmcsv(joinpath(model.tmpdir, "$(model.name)-cmd1-chain$(i).txt"))
      for key in index[:, 1]
        indx1 = idxdct[key][1]
        indx2 = idxdct[key][2]
        if length(keys(tdict)) == 0
          tdict = Dict(key => res[indx1:indx2, 2])
        else
          tdict = merge(tdict, Dict(key => res[indx1:indx2, 2]))
        end
      end
      ## End of processing result type file ##
      ## If any keys were found, merge it in the rtdict ##

      if length(keys(tdict)) > 0
        #println("Merging $(convert(Symbol, res_type)) with keys $(keys(tdict))")
        rtdict = merge(rtdict, Dict(res_type => tdict))
        tdict = Dict{String, Any}()
      end
    end

    ## If rtdict has keys, push it to the chain array ##

    if length(keys(rtdict)) > 0
      #println("Pushing the rtdict with keys $(keys(rtdict))")
      push!(chainarray, rtdict)
      rtdict = Dict{String, Any}()
    end
  end
  (idxdct, chainarray)
end

function read_table_file(model::Jagsmodel, len::Int)
  pdpopt = Dict{String, Any}[]
  if model.dic && model.popt
      numrows = 2len
  else
      numrows = len
  end
  res = readdlmcsv(joinpath(model.tmpdir, "$(model.name)-cmd1-table0.txt"))
  if model.dic && model.popt
    pdpopt = Dict("pD.mean" => res[1:len, 2])
    pdpopt = merge(pdpopt, Dict("popt" => res[len+1:2len, 2]))
  else
    if model.dic
      pdpopt = Dict("pD.mean" => res[1:len, 2])
    else
      pdpopt = Dict("popt" => res[1:len, 2])
    end
  end
  pdpopt
end


text:
# ---------------------------------
#      Configuration functions
# ---------------------------------

"""
    set_config(key::AbstractString, val::AbstractString)

Set configuration string.

# Arguments
* `key::AbstractString`: The key to configure, following table shows possible `key` values, its default setting and its usage.
* `val::AbstractString`: The value to set to the key

Key                     |Default|Description
:-----------------------|:-----:|:---------------------------------------------------------------------------------
ALTERNATIVE_TABLES_DIRECTORY | "" | If provided, this path will be the root directory for the tabular data.  Otherwise, \${HOME}/.CoolProp/Tables is used
ALTERNATIVE_REFPROP_PATH | "" | An alternative path to be provided to the directory that contains REFPROP's fluids and mixtures directories.  If provided, the SETPATH function will be called with this directory prior to calling any REFPROP functions.
ALTERNATIVE_REFPROP_HMX_BNC_PATH | "" | An alternative path to the HMX.BNC file.  If provided, it will be passed into REFPROP's SETUP or SETMIX routines
VTPR_UNIFAQ_PATH | "" | The path to the directory containing the UNIFAQ JSON files.  Should be slash terminated
"""
function set_config(key::AbstractString, val::AbstractString)
    ccall( (:set_config_string, "CoolProp"), Void, (Cstring, Cstring), key, val)
    return get_global_param_string("errstring")
end

"""
    set_config(key::AbstractString, val::Real)

Set configuration numerical value as double.

# Arguments
* `key::AbstractString`: The key to configure, following table shows possible `key` values, its default setting and its usage.
* `val::Real`: The value to set to the key

Key                     |Default|Description
:-----------------------|:-----:|:---------------------------------------------------------------------------------
MAXIMUM_TABLE_DIRECTORY_SIZE_IN_GB |  1.0 | The maximum allowed size of the directory that is used to store tabular data
PHASE_ENVELOPE_STARTING_PRESSURE_PA |  100.0 | Starting pressure [Pa] for phase envelope construction
R_U_CODATA |  8.3144598 | The value for the ideal gas constant in J/mol/K according to CODATA 2014.  This value is used to harmonize all the ideal gas constants. This is especially important in the critical region.
SPINODAL_MINIMUM_DELTA |  0.5 | The minimal delta to be used in tracing out the spinodal; make sure that the EOS has a spinodal at this value of delta=rho/rho_r
"""
function set_config(key::AbstractString, val::Real)
    ccall( (:set_config_double, "CoolProp"), Void, (Cstring, Cdouble), key, val)
    return get_global_param_string("errstring")
end

"""
    set_config(key::AbstractString, val::Bool)

Set configuration value as a boolean.

# Arguments
* `key::AbstractString`: The key to configure, following table shows possible `key` values, its default setting and its usage.
* `val::Bool`: The value to set to the key

Key                     |Default|Description
:-----------------------|:-----:|:---------------------------------------------------------------------------------
NORMALIZE_GAS_CONSTANTS |  true | If true, for mixtures, the molar gas constant (R) will be set to the CODATA value
CRITICAL_WITHIN_1UK |  true | If true, any temperature within 1 uK of the critical temperature will be considered to be AT the critical point
CRITICAL_SPLINES_ENABLED |  true | If true, the critical splines will be used in the near-vicinity of the critical point
SAVE_RAW_TABLES |  false | If true, the raw, uncompressed tables will also be written to file
REFPROP_DONT_ESTIMATE_INTERACTION_PARAMETERS |  false | If true, if the binary interaction parameters in REFPROP are estimated, throw an error rather than silently continuing
REFPROP_IGNORE_ERROR_ESTIMATED_INTERACTION_PARAMETERS |  false | If true, if the binary interaction parameters in REFPROP are unable to be estimated, silently continue rather than failing
REFPROP_USE_GERG |  false | If true, rather than using the highly-accurate pure fluid equations of state, use the pure-fluid EOS from GERG-2008
REFPROP_USE_PENGROBINSON |  false | If true, rather than using the highly-accurate pure fluid equations of state, use the Peng-Robinson EOS
DONT_CHECK_PROPERTY_LIMITS |  false | If true, when possible, CoolProp will skip checking whether values are inside the property limits
HENRYS_LAW_TO_GENERATE_VLE_GUESSES |  false | If true, when doing water-based mixture dewpoint calculations, use Henry's Law to generate guesses for liquid-phase composition
OVERWRITE_FLUIDS |  false | If true, and a fluid is added to the fluids library that is already there, rather than not adding the fluid (and probably throwing an exception), overwrite it
OVERWRITE_DEPARTURE_FUNCTION |  false | If true, and a departure function to be added is already there, rather than not adding the departure function (and probably throwing an exception), overwrite it
OVERWRITE_BINARY_INTERACTION |  false | If true, and a pair of binary interaction pairs to be added is already there, rather than not adding the binary interaction pair (and probably throwing an exception), overwrite it
"""
function set_config(key::AbstractString, val::Bool)
    ccall( (:set_config_bool, "CoolProp"), Void, (Cstring, UInt8), key, val)
    return get_global_param_string("errstring")
end


text:
abstract type Shape{dim} end

getdim(::Shape{dim}) where {dim} = dim

#########################
# Triangles   #
#########################
struct Simplex{dim} <: Shape{dim} end
const Triangle = Simplex{2}
@inline get_num_faces(::Type{Triangle}) = 3
@inline get_num_vertices(::Type{Triangle}) = 3

function reference_coordinates(::Type{Triangle})
    [Tensors.Vec{2, Float64}((0.0, 0.0)),
    Tensors.Vec{2, Float64}((1.0, 0.0)),
    Tensors.Vec{2, Float64}((0.0, 1.0))]
end
function reference_edges(::Type{Triangle})
    [[Tensors.Vec{2, Float64}((1.0, 0.0)),Tensors.Vec{2, Float64}((0.0, 1.0))],
     [Tensors.Vec{2, Float64}((0.0, 1.0)),Tensors.Vec{2, Float64}((0.0, 0.0))],
     [Tensors.Vec{2, Float64}((0.0, 0.0)),Tensors.Vec{2, Float64}((1.0, 0.0))]]
end
@inline reference_edge_nodes(::Type{Triangle}) = ((1,2),(2,3),(3,1))

function gettopology(::Type{Triangle})
    return Dict(0=>3,1=>3,2=>1)
end

##########################
# Tetrahedron
######################

const Tetrahedron = Simplex{3}

get_num_faces(::Type{Tetrahedron}) = 3
get_num_vertices(::Type{Tetrahedron}) = 3

function reference_coordinates(::Type{Tetrahedron})
    [Tensors.Vec{3, Float64}((0.0, 0.0,0.0)),
     Tensors.Vec{3, Float64}((1.0, 0.0,0.0)),
     Tensors.Vec{3, Float64}((0.0,1.0, 0.0)),
     Tensors.Vec{3, Float64}((0.0,0.0, 1.0))]
end
function reference_edges(::Type{Tetrahedron})
    coords = reference_coordinates(Tetrahedron)
    [[coords[i[1]], coords[i[2]]] for i in reference_edge_nodes(Tetrahedron)]
end
function reference_faces(::Type{Tetrahedron})
    coords = reference_coordinates(Tetrahedron)
    [[coords[i[1]], coords[i[2]], coords[i[3]]] for i in reference_face_nodes(Tetrahedron)]
end
reference_edge_nodes(::Type{Tetrahedron}) = ((1,2),(2,3),(3,4),(4,1),(1,3),(2,4))
reference_face_nodes(::Type{Tetrahedron}) = ((1,2,3),(2,3,4),(1,3,4),(1,2,4))

function gettopology(::Type{Tetrahedron})
    return Dict(0=>4,1=>6,2=>4,3=>1)
end

"""
get points for a nodal basis of order `order` on a `dim`
    dimensional simplex
"""
function get_nodal_points(::Type{Simplex{dim}}, order) where {dim}
    points = Vector{Tensors.Vec{dim,Float64}}()
    vertices = reference_coordinates(Simplex{dim})
    topology = Dict{Int, Int}()
    append!(points, vertices)
    push!(topology, 0=>length(points))
    [append!(points, _interior_points(verts, order)) for verts in reference_edges(Simplex{dim})]
    push!(topology, 1=>length(points)-topology[0])
    if dim == 2
        append!(points, _interior_points(vertices, order))
        push!(topology, 2=>length(points)-topology[0]-topology[1])
    else
        [append!(points, _interior_points(verts, order)) for verts in reference_faces(Simplex{dim})]
        push!(topology, 2=>length(points)-topology[0]-topology[1])
        append!(points, _interior_points(vertices, order))
        push!(topology, 3=>length(points)-topology[0]-topology[1]-topology[2])
    end
    points, topology
end

function _interior_points(verts, order)
    n = length(verts)
    ls = [(verts[i] - verts[1])/order for i in 2:n]
    m = length(ls)
    grid_indices =  []
    if m == 1
        grid_indices = [[i] for i in 1:order-1]
    elseif m == 2 && order > 2
        grid_indices = [[i,j] for i in 1:order-1 for j in 1:order-i-1]
    elseif m == 3 && order > 3
        grid_indices = [[i,j,k] for i in 1:order-3 for j in 1:order-i-2 for k in 1:order-1-j-i]
    end
    pts = Vector{typeof(verts[1])}()
    for indices in grid_indices
        res = verts[1]
        for (i,ii) in enumerate(indices)
            res += (ii) * ls[m - i+1]
        end
        push!(pts,res)
    end
    pts
end

#########################
# Rectangles  #
#########################
struct HyperCube{dim} <: Shape{dim} end
const Rectangle = HyperCube{2}
@inline get_num_faces(::Type{Rectangle}) = 4
@inline get_num_vertices(::Type{Rectangle}) = 4

function reference_coordinates(::Type{Rectangle})
    [Tensors.Vec{2, Float64}((0.0, 0.0)),
    Tensors.Vec{2, Float64}((1.0, 0.0)),
    Tensors.Vec{2, Float64}((0.0, 1.0)),
    Tensors.Vec{2, Float64}((1.0, 1.0))]
end
function reference_edges(::Type{Rectangle})
    [[Tensors.Vec{2, Float64}((1.0, 0.0)),Tensors.Vec{2, Float64}((1.0, 1.0))],
     [Tensors.Vec{2, Float64}((1.0, 1.0)),Tensors.Vec{2, Float64}((0.0, 1.0))],
     [Tensors.Vec{2, Float64}((0.0, 1.0)),Tensors.Vec{2, Float64}((0.0, 0.0))],
     [Tensors.Vec{2, Float64}((0.0, 0.0)),Tensors.Vec{2, Float64}((1.0, 0.0))]]
end
@inline reference_edge_nodes(::Type{Rectangle}) = ((1,2),(2,3),(3,4),(4,1))

function gettopology(::Type{Rectangle})
    return Dict(0=>4,1=>4,2=>1)
end

##########################
# Hexahedron
######################

const Hexahedron = HyperCube{3}

get_num_faces(::Type{Hexahedron}) = 6
get_num_vertices(::Type{Hexahedron}) = 8

function reference_coordinates(::Type{Hexahedron})
    [Tensors.Vec{3}([0.0, 0.0, 0.0]),
    Tensors.Vec{3}([1.0, 0.0, 0.0]),
    Tensors.Vec{3}([1.0, 1.0, 0.0]),
    Tensors.Vec{3}([0.0, 1.0, 0.0]),
    Tensors.Vec{3}([0.0, 0.0, 1.0]),
    Tensors.Vec{3}([1.0, 0.0, 1.0]),
    Tensors.Vec{3}([1.0, 1.0, 1.0]),
    Tensors.Vec{3}([0.0, 1.0, 1.0])]
end

function reference_edges(::Type{Hexahedron} )
    coords = reference_coordinates(Hexahedron)
    [[coords[1],coords[2]],[coords[2],coords[3]],
    [coords[3],coords[4]],[coords[4],coords[1]],
    [coords[5],coords[6]],[coords[6],coords[7]],
    [coords[7],coords[8]],[coords[8],coords[5]],
    [coords[1],coords[4]],[coords[2],coords[3]],
    [coords[6],coords[7]],[coords[5],coords[8]]]
end

function reference_faces(::Type{Hexahedron})
    coords = reference_coordinates(Hexahedron)
    [[coords[1],coords[2],coords[3],coords[4]],
    [coords[5],coords[6],coords[7],coords[8]],
    [coords[2],coords[3],coords[7],coords[6]],
    [coords[1],coords[5],coords[8],coords[4]],
    [coords[1],coords[2],coords[5],coords[6]],
    [coords[3],coords[4],coords[8],coords[7]]]
end

reference_edge_nodes(::Type{Hexahedron}) = ((1,2),(2,3),(3,4),(4,1),(5,6),(6,7),(7,8),(8,5),(1,4),(2,3),(6,7),(5,8))
reference_face_nodes(::Type{Hexahedron}) = ((1,2,3,4),(6,5,8,7),(2,6,7,3),(1,4,8,5),(3,7,8,4),(1,5,6,2))

function gettopology(::Type{Hexahedron})
    return Dict(0=>8,1=>12,2=>6,3=>1)
end

"""
get points for a nodal basis of order `order` on a `dim`
    dimensional hypercube
"""
function get_nodal_points(::Type{HyperCube{dim}}, order) where {dim}
    points = Vector{Tensors.Vec{dim,Float64}}()
    vertices = reference_coordinates(HyperCube{dim})
    topology = Dict{Int, Int}()
    append!(points, vertices)
    push!(topology, 0=>length(points))
    [append!(points, _interior_points(verts, order)) for verts in reference_edges(HyperCube{dim})]
    push!(topology, 1=>length(points)-topology[0])
    if dim == 2
        append!(points, _cube_interior_points(vertices, order, dim))
        push!(topology, 2=>length(points)-topology[0]-topology[1])
    else
        [append!(points, _cube_interior_points(verts, order, dim)) for verts in reference_faces(HyperCube{dim})]
        push!(topology, 2=>length(points)-topology[0]-topology[1])
        append!(points, _cube_interior_points(vertices, order, dim))
        push!(topology, 3=>length(points)-topology[0]-topology[1]-topology[2])
    end
    points, topology
end

function _cube_interior_points(verts, order, dim)
    n = length(verts)
    ls = [(verts[i] - verts[1])/order for i in 2:n]
    m = length(ls)
    lx = ls[1]

    if dim == 2
        for point in ls
            if !any([x ≈ 0 for x in point])
                lx = point
            end
        end
        return [eltype(verts)((verts[1][1]+i*lx[1],verts[1][2]+j*lx[2])) for i in 1:order-1 for j in 1:order-1]
    elseif dim == 3
        if m == 3
            for point in ls
                if sum(point .≈ 0.0) == 1
                    lx = point
                end
            end
            if lx[1] ≈ 0.0
                return [eltype(verts)((verts[1][1],verts[1][2]+j*lx[2], verts[1][3]+k*lx[3])) for j in 1:order-1 for k in 1:order-1]
            elseif lx[2] ≈ 0.0
                return [eltype(verts)((verts[1][1]+i*lx[1],verts[1][2], verts[1][3]+k*lx[3])) for i in 1:order-1 for k in 1:order-1]
            elseif lx[3] ≈ 0.0
                return [eltype(verts)((verts[1][1]+i*lx[1],verts[1][2]+j*lx[2], verts[1][3])) for i in 1:order-1 for j in 1:order-1]
            end
        elseif m == 7 && order > 1
            for point in ls
                if !any([x ≈ 0 for x in point])
                    lx = point
                end
            end
            return [eltype(verts)((verts[1][1]+i*lx[1],verts[1][2]+j*lx[2], verts[1][3]+k*lx[3])) for i in 1:order-1 for j in 1:order-1 for k in 1:order-1]
        end
    else
        error("interior nodes for hypercube of dimension $dim not supported")
    end
    return []
end

############
# Segment
###########
const Segment = Simplex{1}
@inline get_num_faces(::Type{Segment}) = 1
@inline get_num_vertices(::Type{Segment}) = 2

function reference_coordinates(::Type{Segment})
    return [Tensors.Vec{1, Float64}((0.0,)),Tensors.Vec{1, Float64}((1.0,))]
end

function get_nodal_points(::Type{Segment}, order)
    points = Vector{Tensors.Vec{1,Float64}}()
    vertices = reference_coordinates(Segment)
    topology = Dict{Int, Int}()
    append!(points, vertices)
    push!(topology, 0=>length(points))
    append!(points, _interior_points(vertices, order))
    push!(topology, 1=>length(points)-topology[0])
    points, topology
end

### Utils
function map_shape_symbols(symbol)
    if symbol == :Triangle
        Triangle()
    elseif symbol == :Segment
        Segment()
    elseif symbol == :Hexahedron
        Hexahedron()
    elseif symbol == :Quad
        Rectangle()
    elseif symbol == :Tetrahedron
        Tetrahedron()
    else
        error("Shape not available")
    end
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
using Test
using Bigsimr
using LinearAlgebra
using Distributions

@testset "Correlation Functions" begin

    r_negdef = [
        1.00 0.82 0.56 0.44
        0.82 1.00 0.28 0.85
        0.56 0.28 1.00 0.22
        0.44 0.85 0.22 1.00
    ]

    @testset "Nearest positive definite correlation matrix" begin
        r = cor_nearPD(r_negdef)
        @test Bigsimr.iscorrelation(r)

        # Must respect input eltype
        for T in [Float64, Float32, Float64]
            @test eltype(cor_nearPD(Matrix{T}(r_negdef))) === T
        end
    end

    @testset "Nearest positive semi-definite correlation matrix" begin
        r = cor_nearPD(r_negdef, 0.0)
        λ = eigvals(r)
        @test issymmetric(r)
        @test all(λ .≥ 0)
        @test all(diag(r) .== 1.0)
        @test all(-1.0 .≤ r .≤ 1.0)
    end

    @testset "Fast near positive definite correlation matrix" begin
        r = cor_fastPD(r_negdef)
        @test Bigsimr.iscorrelation(r)

        # Must respect input eltype
        test_types = [Float64, Float32, Float16]
        for T in test_types
            @test eltype(cor_fastPD(Matrix{T}(r_negdef))) === T
        end
    end
end

@testset "Random Correlation Generation" begin

    @testset "Random postive definite correlation matrix" begin
        r = cor_randPD(100)
        @test Bigsimr.iscorrelation(r)

        # The element type must be respected
        test_types = [Float64, Float32, Float16]
        for T in test_types
            @test eltype(cor_randPD(T, 4)) === T
        end

        # Must work for numbers with integer representations
        test_types = [Float64, Float32, Float16, Rational, Int64, Int32, Int16]
        for T in test_types
            @test_nowarn cor_randPD(T(4))
            @test_nowarn cor_randPD(Float64, T(4))
            @test_nowarn cor_randPD(Float32, T(4))
            @test_nowarn cor_randPD(Float16, T(4))

            for S in test_types
                @test_nowarn cor_randPD(T(4), S(3))
                @test_nowarn cor_randPD(Float64, T(4), S(3))
                @test_nowarn cor_randPD(Float32, T(4), S(3))
                @test_nowarn cor_randPD(Float16, T(4), S(3))
            end
        end

        # `d` must not be less than 1
        @test_throws AssertionError cor_randPD(-1)
        @test_throws AssertionError cor_randPD(Float64, -1)
        @test_throws AssertionError cor_randPD(Float32, -1)
        @test_throws AssertionError cor_randPD(Float16, -1)

        # `k` must not be larger than `d`
        @test_throws AssertionError cor_randPD(4, 5)
        @test_throws AssertionError cor_randPD(Float64, 4, 5)
        @test_throws AssertionError cor_randPD(Float32, 4, 5)
        @test_throws AssertionError cor_randPD(Float16, 4, 5)

        # `k` must not be less than 1
        @test_throws AssertionError cor_randPD(4, 0)
        @test_throws AssertionError cor_randPD(Float64, 4, 0)
        @test_throws AssertionError cor_randPD(Float32, 4, 0)
        @test_throws AssertionError cor_randPD(Float16, 4, 0)
    end

    @testset "Random positive semi-definite correlation matrix" begin
        r = cor_randPSD(100)
        λ = eigvals(r)
        @test issymmetric(r)
        @test all(λ .≥ 0)
        @test all(diag(r) .== 1.0)
        @test all(-1.0 .≤ r .≤ 1.0)

        # The element type must be respected
        test_types = [Float64, Float32, Float16]
        for T in test_types
            @test eltype(cor_randPSD(T, 4)) === T
        end

        # Must work for numbers with integer representations
        test_types = [Float64, Float32, Float16, Rational, Int64, Int32, Int16]
        for T in test_types
            @test_nowarn cor_randPSD(T(4))
            @test_nowarn cor_randPSD(Float64, T(4))
            @test_nowarn cor_randPSD(Float32, T(4))
            @test_nowarn cor_randPSD(Float16, T(4))

            for S in test_types
                @test_nowarn cor_randPSD(T(4), S(3))
                @test_nowarn cor_randPSD(Float64, T(4), S(3))
                @test_nowarn cor_randPSD(Float32, T(4), S(3))
                @test_nowarn cor_randPSD(Float16, T(4), S(3))
            end
        end

        # `d` must not be less than 1
        @test_throws AssertionError cor_randPSD(-1)
        @test_throws AssertionError cor_randPSD(Float64, -1)
        @test_throws AssertionError cor_randPSD(Float32, -1)
        @test_throws AssertionError cor_randPSD(Float16, -1)

        # `k` must not be larger than `d`
        @test_throws AssertionError cor_randPSD(4, 5)
        @test_throws AssertionError cor_randPSD(Float64, 4, 5)
        @test_throws AssertionError cor_randPSD(Float32, 4, 5)
        @test_throws AssertionError cor_randPSD(Float16, 4, 5)

        # `k` must not be less than 1
        @test_throws AssertionError cor_randPSD(4, 0)
        @test_throws AssertionError cor_randPSD(Float64, 4, 0)
        @test_throws AssertionError cor_randPSD(Float32, 4, 0)
        @test_throws AssertionError cor_randPSD(Float16, 4, 0)
    end

end

@testset "Correlation Utilities" begin

    cor_types = (Pearson, Spearman, Kendall)

    @testset "Correlation calculation" begin
        # Must work for any floating point matrix or vector pair
        test_types = [Float64, Float32, Float16]
        for T in test_types
            A = rand(T, 200, 4)
            x, y = rand(T, 100), rand(T, 100)

            for C in cor_types
                @test_nowarn cor(A, C)
                @test_nowarn cor(x, y, C)

                @test_nowarn cor_fast(A, C)
            end
        end

    end

    @testset "Correlation to correlation conversion" begin

        # Converting type A -> A must result in the same matrix
        for C in cor_types
            r = cor_randPD(4)
            @test r == cor_convert(r, C, C)
        end

        # Must map (-1, 0, 1) onto itself within numerical error
        for C1 in cor_types
            for C2 in cor_types
                @test cor_convert( 0.0, C1, C2)  == 0.0
                @test cor_convert( 1.0, C1, C2)   ≤ 1.0
                @test cor_convert( 1.0, C1, C2)   ≈ 1.0
                @test cor_convert(-1.0, C1, C2)  ≥ -1.0
                @test cor_convert(-1.0, C1, C2)  ≈ -1.0
            end
        end

        # Must work for each type
        test_types = [Float64, Float32, Float16, Rational]
        for T in test_types
            for C1 in cor_types
                for C2 in cor_types
                    @test_nowarn cor_convert(T(0.5), C1,  C2)
                end
            end
        end

        # Must respect these input eltypes
        test_types = [Float64, Float32, Float16]
        for T in test_types
            for C1 in cor_types
                for C2 in cor_types
                    @test eltype(cor_convert(T(0.5), C1,  C2))  === T
                end
            end
        end
    end

    @testset "Constrain to Correlation" begin
        test_types = [Float64, Float32, Float16]
        for T in test_types
            A = rand(T, 10, 10)
            @test_nowarn cor_constrain(A)
            @test eltype(cor_constrain(A)) === T
        end
    end

    @testset "Covariance to Correlation" begin
        test_types = [Float64, Float32, Float16]
        for T in test_types
            R = cor_randPSD(T, 4, 3)
            @test_nowarn cov2cor(R)
            @test eltype(cov2cor(R)) === T
        end
    end

    @testset "Correlation Bounds" begin
        A, B = NegativeBinomial(20, 0.2), LogNormal(3, 1)
        dist_types = (A, B)

        # Must work for any univariate distribution and correlation type
        for D1 in dist_types
            for D2 in dist_types
                for C in cor_types
                    @test_nowarn cor_bounds(D1, D2, C)
                end
            end
        end

        # Must work for any number with an integer representation
        test_types = [Float64, Float32, Float16, Int64, Int32, Int16]
        for T in test_types
            @test_nowarn cor_bounds(A, B, n_samples=T(10_000))
        end

        test_types = [Float64, Float32]
        for T in test_types
            @test_throws InexactError cor_bounds(A, B, n_samples=T(10_000.5))
        end
    end

end


text:
module Cthulhu

using InteractiveUtils

using Core: MethodInstance
const Compiler = Core.Compiler

include("callsite.jl")
include("reflection.jl")
include("ui.jl")
include("codeview.jl")

export descend, @descend, descend_code_typed, descend_code_warntype, @descend_code_typed, @descend_code_warntype

"""
    @descend_code_typed

Evaluates the arguments to the function or macro call, determines their
types, and calls `code_typed` on the resulting expression.
"""
macro descend_code_typed(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_typed, ex0)
end

"""
    @descend_code_warntype

Evaluates the arguments to the function or macro call, determines their
types, and calls `code_warntype` on the resulting expression.
"""
macro descend_code_warntype(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_warntype, ex0)
end

"""
    @descend

Shortcut for [`@descend_code_typed`](@ref).
"""
macro descend(ex0...)
    esc(:(@descend_code_typed($(ex0...))))
end

"""
    descend_code_typed(f, tt; kwargs...)

Given a function and a tuple-type, interactively explore the output of
`code_typed` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
function foo()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

descend_code_typed(foo, Tuple{})
```
"""
descend_code_typed(f, @nospecialize(tt); kwargs...) =
    _descend_with_error_handling(f, tt; iswarn=false, kwargs...)

"""
    descend_code_warntype(f, tt)

Given a function and a tuple-type, interactively explore the output of
`code_warntype` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
function foo()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

descend_code_warntype(foo, Tuple{})
```
"""
descend_code_warntype(f, @nospecialize(tt)) =
    _descend_with_error_handling(f, tt; iswarn=true)

function _descend_with_error_handling(f, @nospecialize(tt); kwargs...)
    try
        _descend(f, tt; kwargs...)
    catch x
        if x isa InterruptException
            return nothing
        else
            rethrow(x)
        end
    end
    return nothing
end

"""
    descend

Shortcut for [`descend_code_typed`](@ref).
"""
const descend = descend_code_typed

##
# _descend is the main driver function.
# src/reflection.jl has the tools to discover methods
# src/ui.jl provides the user facing interface to which _descend responds
##
function _descend(mi::MethodInstance; iswarn::Bool, params=current_params(), optimize::Bool=true, kwargs...)
    display_CI = true
    debuginfo = true
    if :debuginfo in keys(kwargs)
        selected = kwargs[:debuginfo]
        # TODO: respect default
        debuginfo = selected == :source
    end

    while true
        (CI, rt, slottypes) = do_typeinf_slottypes(mi, optimize, params)
        preprocess_ci!(CI, mi, optimize)
        callsites = find_callsites(CI, mi, slottypes; params=params, kwargs...)

        debuginfo_key = debuginfo ? :source : :none
        if display_CI
            println()
            println("│ ─ $(string(Callsite(-1, MICallInfo(mi, rt))))")

            if iswarn
                cthulhu_warntype(stdout, CI, rt, debuginfo_key)
            elseif VERSION >= v"1.1.0-DEV.762"
                show(stdout, CI, debuginfo = debuginfo_key)
            else
                display(CI=>rt)
            end
            println()
        end
        display_CI = true

        TerminalMenus.config(cursor = '•', scroll = :wrap)
        menu = CthulhuMenu(callsites)
        cid = request(menu)
        toggle = menu.toggle

        if toggle === nothing
            if cid == length(callsites) + 1
                break
            end
            if cid == -1
                throw(InterruptException())
            end
            callsite = callsites[cid]

            if callsite.info isa MultiCallInfo
                sub_callsites = map(ci->Callsite(callsite.id, ci), callsite.info.callinfos)
                menu = CthulhuMenu(sub_callsites, sub_menu=true)
                cid = request(menu)
                if cid == length(sub_callsites) + 1
                    continue
                end
                if cid == -1
                    throw(InterruptException())
                end
                callsite = sub_callsites[cid]
            end

            if callsite.info isa GeneratedCallInfo || callsite.info isa FailedCallInfo
                @error "Calliste %$(callsite.id) failed to be extracted" callsite
            end

            # recurse
            next_mi = get_mi(callsite)
            if next_mi === nothing
                continue
            end
            _descend(next_mi; params=params, optimize=optimize,
                     iswarn=iswarn, debuginfo=debuginfo_key, kwargs...)
        elseif toggle === :warn
            iswarn ⊻= true
        elseif toggle === :optimize
            optimize ⊻= true
        elseif toggle === :debuginfo
            debuginfo ⊻= true
        elseif toggle === :llvm
            cthulhu_llvm(stdout, mi, optimize, debuginfo, params)
            display_CI = false
        elseif toggle === :native
            cthulhu_native(stdout, mi, optimize, debuginfo, params)
            display_CI = false
        elseif toggle === :dump_params
            @info "Dumping inference cache"
            Core.show(map(((i, x),) -> (i, x.result, x.linfo), enumerate(params.cache)))
            display_CI = false
        else
            error("Unknown option $toggle")
        end
    end
end

function _descend(@nospecialize(F), @nospecialize(TT); params=current_params(), kwargs...)
    mi = first_method_instance(F, TT; params=params)
    _descend(mi; params=params, kwargs...)
end

end


text:
## USER FRIENDLY INSPECTION OF COMPOSITE MACHINES

try_scalarize(v) = length(v) == 1 ? v[1] : v

# get the component models of a composite model `mach` (versions at
# last fit!!) and the corresponding field names (filtering out fields
# that are not models):
function models_and_names(mach::Machine{M}) where  M<:Composite
    raw_names = fieldnames(M)
    _values = Tuple(getproperty(mach.old_model, name) for name in raw_names)
    raw_meta = zip(_values, raw_names) |> collect
    meta = filter(raw_meta) do (value, name)
        value isa Model
    end
    return tuple(zip(meta...)...)
end

# given a dictionary keyed on machines, with "items" as values
# (use-cases are reports or fitted_params) return a named-tuple, keyed
# on model names, whose values are corresponding items or vectors of
# items, but restricted to the specified `models`.
function tuple_keyed_on_model_names(item_given_machine, _models, names)
    machs = keys(item_given_machine) |> collect
    named_tuple_values = map(_models) do model
        sub_machs =
            filter(mach -> mach.old_model == model, machs)
        [item_given_machine[mach] for mach in sub_machs] |> try_scalarize
    end
    return NamedTuple{names}(Tuple(named_tuple_values))
end

function report(mach::Machine{<:Composite})
    machines = mach.report.machines
    dict = mach.report.report_given_machine
    _models, names =  models_and_names(mach)
    return merge(tuple_keyed_on_model_names(dict, _models, names),
                 (machines=machines, report_given_machine=dict,))
end

function fitted_params(mach::Machine{<:Composite})
    fp = fitted_params(mach.model, mach.fitresult)
    machines = fp.machines
    dict = fp.fitted_params_given_machine
    _models, names =  models_and_names(mach)
    return merge(tuple_keyed_on_model_names(dict, _models, names),
                 (machines=machines, fitted_params_given_machine=dict,))
end


text:
# This file includes all the plot functions for the MG
function plot_operation(mg::Microgrid ; y=2, s=1)
    # Seaborn configuration
    Seaborn.set(context="notebook", style="ticks", palette="muted", font="serif", font_scale=1.5)

    # Parameters
    nh = mg.parameters.nh
    Δh = mg.parameters.Δh
    hours = range(1, length = nh, step = Δh) / Δh

    # Plots
    # Powers
    f = figure("Powers")
    for (i, type) in enumerate([typeof(Electricity()), typeof(Heat()), typeof(Hydrogen())])
        i == 1 ? subplot(3, 1, i) : subplot(3, 1, i, sharex = f.axes[1])
        # Demands
        for (k, a) in enumerate(mg.demands)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Demand ", k))
            end
        end
        # Generations
        for (k, a) in enumerate(mg.generations)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Generation ", k))
            end
        end
        # Storages
        for (k, a) in enumerate(mg.storages)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Storage ", k))
            end
        end
        # Converters
        for (k, a) in enumerate(mg.converters)
            for c in a.carrier
                if c isa type
                    plot(hours, c.power[:,y,s], label = string("Converter ", k))
                end
            end
        end
        for (k, a) in enumerate(mg.grids)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Grids ", k))
            end
        end
        legend()
    end
    # State of charge
    figure("State-of-charge")
    for (k, a) in enumerate(mg.storages)
        k == 1 ? subplot(length(mg.storages), 1, k) : subplot(length(mg.storages), 1, k, sharex = f.axes[1])
        plot(hours, a.soc[1:end-1, y, s], label = string("Storage ", k))
        legend()
    end
end
function plot_operation(mg::Microgrid, controller::AbstractController; y=2, s=1)
    # Seaborn configuration
    Seaborn.set(context="notebook", style="ticks", palette="muted", font="serif", font_scale=1.5)

    # Parameters
    nh = mg.parameters.nh
    Δh = mg.parameters.Δh
    hours = range(1, length = nh, step = Δh) / Δh

    # Plots
    # Powers
    f = figure("Powers")
    for (i, type) in enumerate([typeof(Electricity()), typeof(Heat()), typeof(Hydrogen())])
        i == 1 ? subplot(3, 1, i) : subplot(3, 1, i, sharex = f.axes[1])
        # Demands
        for (k, a) in enumerate(mg.demands)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Demand ", k))
            end
        end
        # Generations
        for (k, a) in enumerate(mg.generations)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Generation ", k))
            end
        end
        # Storages
        for (k, a) in enumerate(mg.storages)
            if a.carrier isa type
                plot(hours, controller.decisions.storages[k][:,y,s], label = string("Storage ", k))
            end
        end
        # Converters
        for (k, a) in enumerate(mg.converters)
            for c in a.carrier
                if c isa type
                    plot(hours, controller.decisions.converters[k][:,y,s], label = string("Converter ", k))
                end
            end
        end
        for (k, a) in enumerate(mg.grids)
            if a.carrier isa type
                plot(hours, a.carrier.power[:,y,s], label = string("Grids ", k))
            end
        end
        legend()
    end
end
# Statistics
function plot_metrics(metrics::Metrics; type = "hist")
    # Seaborn configuration
    Seaborn.set(context="notebook", style="ticks", palette="muted", font="serif", font_scale = 1.5)

    if type == "hist"
        figure("Renewable share")
        hist(reshape(metrics.renewable_share[2:end, :], :) * 100, color="steelblue")
        ylabel("Counts", size = "large"), yticks(size = "medium")
        xlabel("Renewable share (%)", size = "large"), xticks(size = "medium")
        figure("Cost")
        hist(reshape(metrics.eac.total[:, :], :) / 1000, color="steelblue")
        ylabel("Counts", size = "large"), yticks(size = "medium")
        xlabel("Annual cost (k€/y)", size = "large"), xticks(size = "medium")
        if !isa(metrics.lpsp.elec, Nothing)
            figure("LPSP elec")
            hist(reshape(metrics.lpsp.elec[2:end, :], :) * 100, color="steelblue")
            ylabel("Counts", size = "large"), yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
        if !isa(metrics.lpsp.heat, Nothing)
            figure("LPSP heat")
            hist(reshape(metrics.lpsp.heat[2:end, :], :) * 100, color="steelblue")
            ylabel("Counts", size = "large"), yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
        if !isa(metrics.lpsp.hydrogen, Nothing)
            figure("LPSP hydrogen")
            hist(reshape(metrics.lpsp.hydrogen[2:end, :], :) * 100, color="steelblue")
            ylabel("Counts", size = "large"), yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
    elseif type == "violin"
        figure("Renewable share")
        violinplot(reshape(metrics.renewable_share[2:end, :], :) * 100, color="steelblue")
        yticks(size = "medium")
        xlabel("Renewable share (%)", size = "large"), xticks(size = "medium")
        figure("Cost")
        violinplot(reshape(metrics.eac.total[:, :], :) / 1000, color="steelblue")
        yticks(size = "medium")
        xlabel("Annual cost (k€/y)", size = "large"), xticks(size = "medium")
        if !isa(metrics.lpsp.elec, Nothing)
            figure("LPSP elec")
            violinplot(reshape(metrics.lpsp.elec[2:end, :], :) * 100, color="steelblue")
            yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
        if !isa(metrics.lpsp.heat, Nothing)
            figure("LPSP heat")
            violinplot(reshape(metrics.lpsp.heat[2:end, :], :) * 100, color="steelblue")
            yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
        if !isa(metrics.lpsp.hydrogen, Nothing)
            figure("LPSP hydrogen")
            violinplot(reshape(metrics.lpsp.hydrogen[2:end, :], :) * 100, color="steelblue")
            yticks(size = "medium")
            xlabel("LPSP (%)", size = "large"), xticks(size = "medium")
        end
    else
        println("Only 'hist' or 'violin' type accepted")
    end
end


text:
function TetgenIO(
    points::Vector{Point{3,T}};
    pointattributes = SVector{0,T}[],
    pointmtrs = SVector{0,T}[],
    pointmarkers = Cint[],
    tetrahedrons = SimplexFace{4,Cint}[],
    tetrahedronattributes = T[],
    tetrahedronvolumes = T[],
    neighbors = Cint[],
    facets = JLFacet{T,Vector{Cint}}[],
    facetmarkers = Cint[],
    holes = Point{3,T}[],
    regions = Region{T}[],
    facetconstraints = FacetConstraint{inttype(T),T}[],
    segmentconstraints = SegmentationConstraint{inttype(T),T}[],
    trifaces = TriangleFace{Cint}[],
    trifacemarkers = Cint[],
    edges = LineFace{Cint}[],
    edgemarkers = Cint[],
) where T
    TetgenIO(
        points,
        pointattributes,
        pointmtrs,
        pointmarkers,
        tetrahedrons,
        tetrahedronattributes,
        tetrahedronvolumes,
        neighbors,
        facets,
        facetmarkers,
        holes,
        regions,
        facetconstraints,
        segmentconstraints,
        trifaces,
        trifacemarkers,
        edges,
        edgemarkers,
    )

end


text:
using DrWatson

@quickactivate

using Pluto

function real_main()
    Pluto.run(
        notebook="Mage.jl",
        workspace_use_distributed=false,
        port=8001
    )
end

function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

julia_main()


text:
# Each pixel is given a unique ID, missings are automatically considered no data
# and elements with value equal to no_data_val are also ignored (not considered valid nodes)
"""
    construct_nodemap(A::Matrix{T} where T <: Real;
                      no_data_val = nothing)
    construct_nodemap(A::GeoData.GeoArray)

Construct a nodemap from `A` and with dims equal to `size(A)` that can be used
to provide spatial references for each vertex in a graph representation of
`A`. Each element in the nodemap is given a unique integer ID. Elements in the
nodemap that correspond to NoData in `A` (`A.missingval` if `A` is a GeoArray,
or `no_data_val` if `A` is a matrix) are given a value of 0.
"""
function construct_nodemap(A::Matrix{T} where T<: Real;
                           no_data_val = nothing)
    dims = size(A)

    # Make an resistance of unique node identifiers
    nodemap = zeros(Int64, dims)
    is_node = (A .!= no_data_val) .&
        ((!).(isnan.(A)))
    nodemap[is_node] = 1:sum(is_node)

    nodemap
end

function construct_nodemap(A::GeoData.GeoArray)
    # Handle dimensions, get them so the order matches the input
    lat_lon_dims = get_lat_lon_dims(A)
    A_band1 = A[Band(Between(1, 1))]

    # Make an array of unique node identifiers
    nodemap = zeros(Int64, size(A_band1.data))
    is_node = (A_band1.data .!= A.missingval) .&
        ((!).(isnan.(A_band1.data)))

    nodemap[is_node] = 1:sum(is_node)

    nodemap = GeoData.GeoArray(nodemap, dims = lat_lon_dims)

    nodemap
end

"""
    construct_weighted_graph(cost_surface::Matrix{T} where T <: Real,
                             nodemap::Matrix{Int};
                             no_data_val = nothing,
                             cost_layer_is_conductance = false,
                             connect_four_neighbors_only = false,
                             connect_using_avg_cost = true)

    construct_weighted_graph(cost_surface::GeoData.GeoArray,
                             nodemap::GeoData.GeoArray;
                             cost_layer_is_conductance = false,
                             connect_four_neighbors_only = false,
                             connect_using_avg_cost = true)

Construct a impleWeightedGraph from a `cost_surface` representing traversal costs`
and a `nodemap` representing the locations in space (either geographic or
cartesian) of the vertices to be added to the graph. The resulting graph will
contain edges connecting neighboring elements in `nodemap` with weights
calculated using the corresponding elements in `cost_surface` for each element
in the nodemap. For cardinal neighbors with cartesian coordinates [i, j] and
[i+1, j], the corresponding edge connecting node `nodemap[i, j]` to node
`nodemap[i + 1, j]`. It's weight (traversal cost) equals
`(cost_surface[i, j] + cost_surface[i+1, j]) / 2`. For diagonal neighbors,
[i, j] and [i+1, j+1], the edge weight equals
`(2/√2)*(cost_surface[i, j] + cost_surface[i+1, j+1]) / 2`. The average of the
two costs is multiplied by 2/√2 to account for the increased distance between
diagonal neighbors.

## Keyword Arguments
- `no_data_val`: Real. The value corresponding to "no data" in `cost_surface`.
Elements equal to `no_data_val` correspond to pixels that cannot be traversed.
Defaults to `nothing`, or, if cost_surface is a GeoData.GeoArray,
cost_surface.missingval is used.

- `cost_layer_is_conductace`: Boolean. Does the `cost_surface` represent
permeability/conductance instead of cost/resistance? If `true`, cost is
calculated as `1 / cost_surface`. Defaults to `false`.

- `connect_four_neighbors_only`: Boolean. Connect only cardinal neighbors in
`nodemap`? Defaults to `false`.

- `connect_using_avg_cost`: `Boolean`. This is intended to offer methods that
complement those used in Circuitscape.jl and Omniscape.jl. In this context,
cost is in units of electrical resistance. If `false`, the cost between two
nodes with resistances R1 and R2 is calculated by converting resistance to
conductances, taking the average, then taking the inverse of the result to
convert back to resistance: `1 / ((1/R1 + 1/R2) / 2)`.
`connect_using_avg_cost = false` correspondes to the default settings in
Circuitscape. Defaults to `true'.
"""
function construct_weighted_graph(cost_surface::Matrix{T} where T <: Real,
                                  nodemap::Matrix{Int};
                                  no_data_val = nothing,
                                  cost_layer_is_conductance::Bool = false,
                                  connect_four_neighbors_only::Bool = false,
                                  connect_using_avg_cost::Bool = true)
    cost = float.(deepcopy(cost_surface))

    # Which averaging function to use
    card_avg = connect_using_avg_cost ? res_cardinal_avg : cond_cardinal_avg
    diag_avg = connect_using_avg_cost ? res_diagonal_avg : cond_diagonal_avg
    dims = size(cost)
    not_no_data = cost .!= no_data_val

    if sum((cost .<= 0) .& not_no_data) != 0
        @error("cost surface contains 0 or negative values (aside from the provided no_data_val), which is not supported")
    end

    if cost_layer_is_conductance
        cost[not_no_data] = 1 ./ cost[not_no_data]
    end

    # Construct graph
    #res_graph = SimpleWeightedGraph(sum(is_node))
    sources = Vector{Int64}()
    destinations = Vector{Int64}()
    node_weights = Vector{Float64}()

    # Add the edges
    # only need to do neighbors down or to the right because
    # the graph is undirected and edge additions will be redundant
    for row in 1:dims[1]
        for column in 1:dims[2]
            if nodemap[row, column] == 0 || cost[row, column] == no_data_val
                continue
            else
                ## Add cardinal neighbors
                # East
                if column != dims[2] && nodemap[row, column + 1] != 0 && cost[row, column + 1] != no_data_val
                    res = card_avg(cost[row, column],
                                   cost[row, column + 1])
                    push!(sources, nodemap[row, column])
                    push!(destinations, nodemap[row, column + 1])
                    push!(node_weights, res)
                end
                # South
                if row != dims[1] && nodemap[row + 1, column] != 0  && cost[row + 1, column] != no_data_val
                    res = card_avg(cost[row, column],
                                   cost[row + 1, column])
                    push!(sources, nodemap[row, column])
                    push!(destinations, nodemap[row + 1, column])
                    push!(node_weights, res)
                end

                ## Add diagonal neighbors if needed
                if !connect_four_neighbors_only
                    # Northeast
                    if column != dims[2] && row != 1 && nodemap[row - 1, column + 1] != 0 && cost[row - 1, column + 1] != no_data_val
                        res = diag_avg(cost[row, column],
                                       cost[row - 1, column + 1])
                        push!(sources, nodemap[row, column])
                        push!(destinations, nodemap[row - 1, column + 1])
                        push!(node_weights, res)
                    end
                    # Southeast
                    if row != dims[1] && column != dims[2] && nodemap[row + 1, column + 1] != 0  && cost[row + 1, column + 1] != no_data_val
                        res = diag_avg(cost[row, column],
                                       cost[row + 1, column + 1])
                        push!(sources, nodemap[row, column])
                        push!(destinations, nodemap[row + 1, column + 1])
                        push!(node_weights, res)
                    end
                end
            end
        end
    end

    # push!'ing to vectors then combining is way faster than add_edge!
    # Sometimes, depending on the nodemap, there may be multiple edges with
    # different weights connecting the same two vertices (only when vertices
    # span multiple "pixels"). Default to using the min of all duplicates via
    # combine = min
    resistance_graph = SimpleWeightedGraph(sources, destinations, node_weights, combine = min)

    return resistance_graph
end

function construct_weighted_graph(cost_surface::GeoData.GeoArray,
                                  nodemap::GeoData.GeoArray;
                                  cost_layer_is_conductance::Bool = false,
                                  connect_four_neighbors_only::Bool = false,
                                  connect_using_avg_cost::Bool = true)
    construct_weighted_graph(cost_surface.data[:, :, 1],
                             nodemap.data[:, :, 1],
                             no_data_val = cost_surface.missingval,
                             cost_layer_is_conductance = cost_layer_is_conductance,
                             connect_four_neighbors_only = connect_four_neighbors_only,
                             connect_using_avg_cost = connect_using_avg_cost)
end

function get_cartesian_indices(node_array::Matrix{Int})
    # Get Cartesian Index of all nodes once, this method ensures that the index
    # in all_coords (created below) matches the node value, e.g.
    # all_coords[10] is the cartesian index for the node with value 10
    ###########################################
    # THIS WILL ONLY WORK WITH NODEMAPS THAT ##
    # WERE CREATED WITH construct_nodemap()  ##
    ###########################################

    # Throw an error if node is formatted the right way
    nodes_bool = node_array .!= 0
    if node_array[nodes_bool] != collect(1:sum(nodes_bool))
        @error("This function does not support nodemaps that were not created with construct_nodemap")
    end

    all_coords = Vector{CartesianIndex{2}}(undef, maximum(node_array))
    size_dims = size(node_array)
    idx = 1
    # Need to do cols then rows to match how nodemap is created
    # this is WAY faster than a bunch of findall's
    for col in 1:size_dims[2]
        for row in 1:size_dims[1]
            if node_array[row, col] == 0
                continue
            end
            all_coords[idx] = CartesianIndex(row, col)
            idx+=1
        end
    end

    return all_coords
end

function get_cartesian_indices(nodemap::GeoArray)
    get_cartesian_indices(nodemap.data[:, :, 1])
end

max_stars_repo_name: Circuitscape/SpatialGraphs-OLD

max_issues_repo_name: Circuitscape/SpatialGraphs-OLD

max_forks_repo_name: Circuitscape/SpatialGraphs-OLD


text:
using DataFrames
using LinearAlgebra
using Statistics: std, var, mean, cov
using HighFrequencyCovariance
using Random

brownian_corr_matrix = Hermitian([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
assets = [:BARC, :HSBC, :VODL, :RYAL]
twister = MersenneTwister(1)

ts1, true_covar, micro_noise, update_rates = generate_random_path(4, 2000; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister))
ts1.df.value = exp.(ts1.df.Value)

# Testing making a ItoSet
ito = ItoSet(true_covar)

iscloser(a,b) = (a.Correlation_error + a.Volatility_error < b.Correlation_error + b.Volatility_error)

function subset_dict(dic::Dict, keys::Vector)
    new_dict = Dict()
    for k in keys
        new_dict[k] = dic[k]
    end
    return new_dict
end

# Testing Covariation matrix functions.
covar_functions =  [simple_covariance, bnhls_covariance, spectral_covariance, preaveraged_covariance, two_scales_covariance]
function test_covariance(settings = Dict())
    success = true
    for f in covar_functions
        k = Symbol(f)
        args = reduce(union, Base.kwarg_decl.(methods(f)))
        releventsettings = intersect(args, keys(settings) )
        sets = (length(releventsettings) > 0) ? subset_dict(settings, releventsettings)  : Dict()
        estimate   = f(ts1, assets; sets... )
        estimate_c = estimate_covariance(ts1, assets, k; sets...)
        cor_err = calculate_mean_abs_distance(estimate, estimate_c).Correlation_error
        t1 = (cor_err .< 10*eps()) | isnan(cor_err)
        if (t1==false) println("In convenience function covariance test there is a covariance error for ", k) end
        vol_err = calculate_mean_abs_distance(estimate, estimate_c).Volatility_error
        t2 = ( vol_err  .< 10*eps()) | isnan(vol_err)
        if (t2==false) println("In convenience function covariance test there is a volatility error for ", k) end
        success = success & (t1 & t2)
    end
    return success
end
test_covariance()
settings1 = Dict([:kernel , :num_grids] .=>
                 [fejer   , 15])
test_covariance(settings1)
settings2 = Dict([:theta , :num_grids, :only_regulise_if_not_PSD  , :regularisation          , :equalweight] .=>
                 [0.5    , 20        , true                       , :identity_regularisation , true        ])
test_covariance(settings2)
settings3 = Dict{Symbol,Any}([:numJ , :num_blocks, :block_width  ] .=>
                 [Int(110)   , 30        ,  20           ])
test_covariance(settings3)

# Volatility
vol_functions = [simple_volatility, two_scales_volatility]
function test_volatility(settings = Dict())
    success = true
    for f in vol_functions
        k = Symbol(f)
        args = reduce(union, Base.kwarg_decl.(methods(f)))
        releventsettings = intersect(args, keys(settings) )
        sets = (length(releventsettings) > 0) ? subset_dict(settings, releventsettings)  : Dict()
        estimate   = f(ts1, assets; sets... )
        estimate = (Symbol(f) == :two_scales_volatility) ? estimate[1] : estimate
        estimate_c = estimate_volatility(ts1, assets, k; sets...)
        cor_err = calculate_mean_abs_distance(estimate, estimate_c)
        t1 = all( (cor_err .< 10*eps()) .| isnan.(cor_err))
        if (t1==false) println("In convenience function volatility test there is a error for ", k) end
        success = success & (t1)
    end
    return success
end
test_volatility()
settings1 = Dict([:rough_guess_number_of_intervals , :num_grids] .=>
                 [7   ,  2*default_num_grids(ts1)  ])
test_volatility(settings1)

# Testing Regularisation
covar = two_scales_covariance(ts1; regularisation = missing)

reg_functions = [identity_regularisation, eigenvalue_clean, nearest_correlation_matrix, nearest_psd_matrix]
function test_regularisation(covar, settings = Dict())
    success = true
    for f in reg_functions
        k = Symbol(f)
        args = reduce(union, Base.kwarg_decl.(methods(f)))
        releventsettings = intersect(args, keys(settings) )
        sets = (length(releventsettings) > 0) ? subset_dict(settings, releventsettings)  : Dict()
        estimate   = f(covar, ts1; sets... )
        estimate_c = regularise(covar, ts1, k; sets...)
        cor_err = calculate_mean_abs_distance(estimate, estimate_c).Correlation_error
        t1 = (cor_err .< 10*eps()) | isnan(cor_err)
        if (t1==false) println("In convenience function covariance test there is a covariance error for ", k) end
        vol_err = calculate_mean_abs_distance(estimate, estimate_c).Volatility_error
        t2 = ( vol_err  .< 10*eps()) | isnan(vol_err)
        if (t2==false) println("In convenience function covariance test there is a volatility error for ", k) end
        success = success & (t1 & t2)
    end
    return success
end
test_regularisation(covar)
settings1 = Dict{Symbol,Any}([:apply_to_covariance , :identity_weight, :doDykstra] .=>
                 Any[false   ,  0.5,  false ])
test_regularisation(covar, settings1)


text:
using Pkg
Pkg.add("NetworkInference")
using NetworkInference

function run_pidc(Expr, filename = "links.txt")

        nodes = get_nodes(Expr)
        inferred_network = InferredNetwork(PIDCNetworkInference(), nodes)
        write_network_file(filename, inferred_network)

        end

run_pidc("Expr.txt")

max_stars_repo_name: MingBit/SCNetEnrich

max_issues_repo_name: MingBit/SCNetEnrich

max_forks_repo_name: MingBit/SCNetEnrich


text:
"""
    volfunc(when::AbstractVector, ts::TS, γ::Real, δ::Real = 0)::AbstractVector{<:Real}

Compute the volatility function that is the term in the rate of a SEPP that corresponds to the self-excitement, denoted by ν in Li2020.
"""
function volfunc(when::AbstractVector, ts::TS, γ::Real, δ::Real = 0, impact_func::Function = (d -> 1 + tanh(d)))::AbstractVector{<:Real}

    isempty(when) && (return when)
    isempty(ts.times) && return(zero(when))

    γ < 0 && error("γ must be positive or zero")

    (ts isa TimeSeries && δ != 0) && (@warn "no marks but δ non zero")

    first(when) isa TimeType && (when = Dates.value.(DateTime.(when)) ./ (1000*3600*24))
    times = first(ts.times) isa TimeType ? Dates.value.(DateTime.(ts.times)) ./ (1000*3600*24) : copy(ts.times)
    marks = ts isa MarkedTimeSeries ? copy(ts.marks) : fill(0, size(times))

    mts = hcat(times, marks)

    function self_ex(t)

        mts_to_t = mts[times .< t, :]

        term(markedpoint) = impact_func(δ * markedpoint[2]) * exp(-γ * (t - markedpoint[1]))

        isempty(mts_to_t) && return 0

        return sum(term.(eachrow(mts_to_t)))
    end

    return self_ex.(when)
end

include("ParametersEstimation/DiscreteProcessEstimation.jl")
include("ParametersEstimation/ContinuousProcessEstimation.jl")


text:
# some extensions to ITensors.jl

getindex(T::ITensor, pairs::Tuple{Index, Int}...) = getindex(T, [IndexVal(pair...) for pair ∈ pairs]...)
getindex(T::ITensor, pairs::Vector{Tuple{Index, Int}}) = getindex(T, pairs...)

function contract(A::ITensor, B::ITensor, i::Index, j::Index)
	@assert hasindex(IndexSet(A),i)
	@assert hasindex(IndexSet(B),j)
	@assert dim(i) == dim(j)
	pos_i = 1
	pos_j = 1
	while IndexSet(A)[pos_i]!= i
		pos_i+=1
	end
	while IndexSet(B)[pos_j]!= j
		pos_j+=1
	end
	A.inds[pos_i] = j
	C = A*B
	A.inds[pos_i] = i
	return C
end

function clamp(A::ITensor, ind::Index, j::Int)
	v = zeros(dim(ind))
	v[j] = 1.0
	projector = ITensor(v, ind)
	return A*projector
end


text:
module Subspaces

using DocStringExtensions

import Base.hcat, Base.vcat, Base.hvcat, Base.cat, Base.+, Base.*, Base.kron, Base.show, Base.iterate, Base.==, Base.in, Base.adjoint
import Base.|, Base.&, Base.~, Base./, Base.⊆, Base.⊇
import Base.size
using LinearAlgebra
using Convex
using Compat

export Subspace, dim, each_basis_element
export random_subspace, random_hermitian_subspace, empty_subspace, full_subspace
export tobasis, frombasis, random_element, projection, perp
export hermitian_basis
export variable_in_space
export ⟂

default_tol = 1e-6

"""
    Subspace(basis::AbstractArray{AbstractArray{N}, 1}; tol::Real=$(default_tol))

Create a subspace from the given basis, expressed as a list of basis vectors.

    Subspace(basis::AbstractArray{N+1}; tol::Real=$(default_tol))

Create a subspace from a basis given as a multi-dimensional array.  The last array index
enumerates the basis elements.  E.g., if given a matrix this constructor will create a
subspace representing the column span of that matrix.

The `tol` parameter sets the tolerance for determining whether vectors are linearly
dependent.

$(TYPEDFIELDS)

```jldoctest
julia> Subspace([[1, 2, 3], [4, 5, 6]]) == Subspace([ 1 4; 2 5; 3 6])
true
```
"""
struct Subspace{T, N}
    """An orthonormal basis for this subspace.  The final index of this array indexes the
    basis vectors."""
    basis::Array # FIXME can it somehow be Array{T, N+1}?
    """An orthonormal basis for the perpendicular subspace."""
    perp::Array # FIXME can it somehow be Array{T, N+1}?
    """The tolerance for determining whether vectors are linearly dependent."""
    tol::AbstractFloat

    function Subspace(basis::AbstractArray{<:Number}; tol::Real=default_tol)
        shape = size(basis)[1:end-1]
        d = size(basis)[end]
        mat = reshape(basis, prod(shape), d)
        if true
            q = qr(mat, Val(true))
            @compat take = [norm(x) >= tol for x in eachrow(q.R)]
            # This is the most efficient way to get the Q matrix.
            # Multiply by I is slow: https://github.com/JuliaLang/julia/issues/38972
            # Cast to Array changes size: https://github.com/JuliaLang/julia/issues/37102
            qQ = q.Q * Matrix(1.0*I, (size(q.Q)[2], size(q.Q)[2]))
            resize!(take, size(qQ)[2])
            good = qQ[:,take]
            perp = qQ[:,.!take]
        else
            s = svd(mat; full=true)
            lastgood = findlast(s.S .>= tol)
            if typeof(lastgood) == Nothing
                lastgood = 0
            end
            good = s.U[:,1:lastgood]
            perp = s.U[:,lastgood+1:end]
        end
        good = reshape(good, shape..., size(good)[2])
        perp = reshape(perp, shape..., size(perp)[2])
        return new{eltype(good), length(shape)}(good, perp, tol)
    end

    function Subspace(basis::AbstractArray{<:AbstractArray{<:Number}, 1}; tol::Real=default_tol)
        shape = size(basis[1])
        return Subspace(cat(basis...; dims=length(shape)+1); tol=tol)
    end
end

function show(io::IO, S::Subspace)
    print(io, "Subspace{$(eltype(S.basis))} size $(size(S)) dim $(dim(S))")
end

"""
$(TYPEDSIGNATURES)

Returns the size of the elements of a subspace.

```jldoctest
julia> size(Subspace([ [1,2,3], [4,5,6] ]))
(3,)
```
"""
size(S::Subspace)::Tuple{Vararg{<:Integer}} = size(S.basis)[1:end-1]

"""
$(TYPEDSIGNATURES)

Returns the linear dimension of this subspace.

```jldoctest
julia> dim(Subspace([ [1,2,3], [4,5,6] ]))
2
```
"""
dim(S::Subspace)::Integer = size(S.basis)[end]

"""
$(TYPEDSIGNATURES)

Returns the orthogonal subspace.  Can also be written as `~S`.

```jldoctest
julia> S = Subspace([ [1,2,3], [4,5,6] ])
Subspace{Float64} size (3,) dim 2

julia> perp(S)
Subspace{Float64} size (3,) dim 1

julia> perp(S) == ~S
true

julia> perp(S) ⟂ S
true
```
"""
perp(S::Subspace)::Subspace = Subspace(S.perp, tol=S.tol)

"""
$(TYPEDSIGNATURES)

Returns the perpendicular subspace.  Can also be written as `perp(S)`.
"""
(~)(S::Subspace) = perp(S)

# It'd be nice for these to all take SubspaceOrArray but then our overloads seem to be selected even
# when all args are Array.
#const SubspaceOrArray{T, N} = Union{Subspace{T, N}, AbstractArray{T, N}}

"""
$(TYPEDSIGNATURES)

Create a generator that iterates over orthonormal basis vectors of a subspace.
"""
each_basis_element(S::Subspace) = eachslice(S.basis; dims=length(size(S.basis)))

each_basis_element(arr::AbstractArray) = [arr]

function each_basis_element_or_zero(S::Subspace{T, N}) where {T, N}
    if dim(S) == 0
        return [ zeros(T, size(S)) ]
    else
        return each_basis_element(S)
    end
end

"""
    +(a::Subspace, b::Subspace)
    +(a::AbstractArray, b::Subspace)
    +(a::Subspace, b::AbstractArray)
    +(a::UniformScaling, b::Subspace)
    +(a::Subspace, b::UniformScaling)

Linear span of two subspaces, or of a subspace an and array.  Equivalent to |(a, b).
"""
function +(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N}
    if size(a) != size(b)
        throw(DimensionMismatch("Array size mismatch: $(size(a)) vs $(size(b))"))
    end
    return Subspace(cat(a.basis, b.basis; dims=N+1), tol=max(a.tol, b.tol))
end

(+)(a::Subspace{T, N}, b::AbstractArray{U, N}) where {T,U,N} = a + Subspace([b], tol=a.tol)
(+)(a::AbstractArray{T, N}, b::Subspace{U, N}) where {T,U,N} = Subspace([a], tol=b.tol) + b
(+)(S::Subspace{T, 2}, x::UniformScaling) where T = S + Subspace([ Array{T}(I, size(S)) ], tol=S.tol)
(+)(x::UniformScaling, S::Subspace{T, 2}) where T = S + I

"""
    |(a::Subspace, b::Subspace)
    |(a::AbstractArray, b::Subspace)
    |(a::Subspace, b::AbstractArray)
    |(a::UniformScaling, b::Subspace)
    |(a::Subspace, b::UniformScaling)

Linear span of two subspaces, or of a subspace an and array.  Equivalent to |(a, b).
"""
(|)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = a + b
(|)(a::Subspace{T, N}, b::AbstractArray{U, N}) where {T,U,N} = a + b
(|)(a::AbstractArray{T, N}, b::Subspace{U, N}) where {T,U,N} = a + b
(|)(S::Subspace{T, 2}, x::UniformScaling) where T = S + I
(|)(x::UniformScaling, S::Subspace{T, 2}) where T = S + I

"""
    *(a::Subspace, b::Subspace)
    *(a::AbstractArray, b::Subspace)
    *(a::Subspace, b::AbstractArray)

Linear span of products of elements of space `a` with elements of space `b`.
"""
function *(a::Subspace, b::Subspace)
    #if dim(a) == 0 || dim(b) == 0
    #    return Subspace([ zeros(T, size(a)) * zeros(U, size(b)) ], tol=max(a.tol, b.tol))
    #else
        return Subspace([ x*y for x in each_basis_element_or_zero(a) for y in each_basis_element_or_zero(b) ], tol=max(a.tol, b.tol))
    #end
end

(*)(a::Subspace, b::AbstractArray) = a * Subspace([b], tol=a.tol)
(*)(a::AbstractArray, b::Subspace) = Subspace([a], tol=b.tol) * b

"""
$(TYPEDSIGNATURES)

Linear span of Kronecker products of elements of space `a` with elements of space `b`.
"""
function kron(a::Subspace{T,N}, b::Subspace{U,N}) where {T,U,N}
    return Subspace([
        kron(x, y)
        for x in each_basis_element_or_zero(a)
        for y in each_basis_element_or_zero(b)
    ], tol=max(a.tol, b.tol))
end

kron(a::Subspace, b::AbstractArray) = kron(a, Subspace([b], tol=a.tol))
kron(a::AbstractArray, b::Subspace) = kron(Subspace([a], tol=b.tol), b)

"""
$(TYPEDSIGNATURES)

Vertical concatenations of vector subspaces (direct sum).
"""
vcat(S::Subspace...) = cat(S...; dims=1)

"""
$(TYPEDSIGNATURES)

Horizontal concatenations of vector subspaces (direct sum).
"""
hcat(S::Subspace...) = cat(S...; dims=2)

"""
$(TYPEDSIGNATURES)

Concatenations of vector subspaces (direct sum).
"""
function cat(S::Subspace...; dims)
    n = length(S)
    # FIXME doesn't work well with heterogenous types
    #T = promote_type(map((x)->eltype(x.basis), [S...])...)
    Subspace([
        cat([ i==j ? x : zeros(size(S[i])) for i in 1:n ]...; dims=dims)
        for j in 1:n
        for x in each_basis_element_or_zero(S[j])
    ], tol=maximum([ s.tol for s in S ]))
end

"""
$(TYPEDSIGNATURES)

Concatenations of vector subspaces (direct sum).
"""
function hvcat(rows::Tuple{Vararg{Int}}, S::Subspace{T, N}...) where {T, N}
    n = length(S)
    basis = Array{Array{T, N}, 1}()
    for j in 1:n
        for x in each_basis_element_or_zero(S[j])
            push!(basis, hvcat(rows, [ i==j ? x : zeros(T, size(S[i])) for i in 1:n ]...))
        end
    end
    if isempty(basis)
        push!(basis, hvcat(rows, [ zeros(T, size(S[i])) for i in 1:n ]...))
    end
    return Subspace(basis, tol=maximum([ s.tol for s in S ]))
end

"""
$(TYPEDSIGNATURES)

Adjoint of vector subspace (linear span of adjoints of members of a space).
"""
adjoint(S::Subspace) =
    Subspace([ x' for x in each_basis_element_or_zero(S) ], tol=S.tol)

"""
$(TYPEDSIGNATURES)

Check for membership of a vector in a vector subspace.
"""
function in(x::AbstractArray{<:Number, N}, S::Subspace{<:Number, N}) where N
    return norm(x - projection(S, x)) <= S.tol
end

function in(x::UniformScaling, S::Subspace{T, 2}) where T
    return Matrix{T}(I, size(S)) in S
end

"""
$(TYPEDSIGNATURES)

Check whether `a` is a subspace of `b`.
"""
function in(a::Subspace{<:Number, N}, b::Subspace{<:Number, N}) where N
    shp = size(a)
    if dim(a) > dim(b)
        return false
    end
    Ma = reshape(a.basis, prod(shp), dim(a))
    Mb = reshape(b.basis, prod(shp), dim(b))
    s = svdvals(Mb' * Ma)
    tol = max(a.tol, b.tol)
    return all(s .> (1.0 - tol))
end

"""
$(TYPEDSIGNATURES)

Check whether two vector subspaces are equal.
"""
function ==(a::Subspace{<:Number, N}, b::Subspace{<:Number, N}) where N
    return dim(a) == dim(b) && a in b
end

"""
$(TYPEDSIGNATURES)

Intersection of vector subspaces.
"""
(&)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = perp(perp(a) + perp(b))

"""
$(TYPEDSIGNATURES)

Intersection of vector subspaces.
"""
(∩)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = a & b

"""
$(TYPEDSIGNATURES)

Quotient of vector subspaces.  Throws error unless `b ⊆ a`.

```jldoctest
julia> a = random_subspace(ComplexF64, 2, 10)
Subspace{Complex{Float64}} size (10,) dim 2

julia> b = random_subspace(ComplexF64, 3, 10)
Subspace{Complex{Float64}} size (10,) dim 3

julia> a / b
ERROR: ArgumentError: divisor must be a subspace of dividend for subspace quotient

julia> (a+b) / a
Subspace{Complex{Float64}} size (10,) dim 3

julia> (a+b) / a == Subspace([ projection(~a, x) for x in each_basis_element(b) ])
true

julia> ((a+b) / b) ⟂ b
true

```
"""
function (/)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N}
    if !(b ⊆ a)
        throw(ArgumentError("divisor must be a subspace of dividend for subspace quotient"))
    end
    return perp(perp(a) + b)
end

(/)(a::Subspace{T, N}, b::AbstractArray{U, N}) where {T,U,N} = a / Subspace([b], tol=a.tol)
(/)(a::Subspace{T, 2}, b::UniformScaling) where T = a / Array{T}(I, size(a))

"""
$(TYPEDSIGNATURES)

Check whether `a` is a subspace of `b`.
"""
(⊆)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = a in b

"""
$(TYPEDSIGNATURES)

Check whether `b` is a subspace of `a`.
"""
(⊇)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = b in a

"""
$(TYPEDSIGNATURES)

Check whether `a` is orthogonal to `b`.
"""
(⟂)(a::Subspace{T, N}, b::Subspace{U, N}) where {T,U,N} = a ⊆ perp(b)

################
### Math
################

"""
$(TYPEDSIGNATURES)

Returns the basis components of `x` in the basis `S.basis`.

See also: [`frombasis`](@ref).

```jldoctest
julia> S = random_subspace(ComplexF64, 2, 10)
Subspace{Complex{Float64}} size (10,) dim 2

julia> x = randn(10);

julia> size(tobasis(S, x))
(2,)

julia> norm( frombasis(S, tobasis(S, x)) - projection(S, x) ) < 1e-8
true
```
"""
function tobasis(S::Subspace{<:Number, N}, x::AbstractArray{<:Number, N}) where N
    shp = size(S)
    basis_mat = reshape(S.basis, prod(shp), dim(S))
    return basis_mat' * vec(x)
end

"""
$(TYPEDSIGNATURES)

Returns a vector from the given basis components in the basis `S.basis`.

See also: [`tobasis`](@ref).

```jldoctest
julia> S = random_subspace(ComplexF64, 2, 10)
Subspace{Complex{Float64}} size (10,) dim 2

julia> x = randn(2);

julia> size(frombasis(S, x))
(10,)

julia> norm( tobasis(S, frombasis(S, x)) - x ) < 1e-8
true
```
"""
function frombasis(S::Subspace, x::AbstractArray{<:Number, 1})
    shp = size(S)
    basis_mat = reshape(S.basis, prod(shp), dim(S))
    return reshape(basis_mat * x, shp)
end

"""
$(TYPEDSIGNATURES)

Projection of vector `x` onto suspace `S`.
"""
function projection(S::Subspace{<:Number, N}, x::AbstractArray{<:Number, N}) where N
    return frombasis(S, tobasis(S, x))
end

"""
$(TYPEDSIGNATURES)

Returns a random element of a subspace.
"""
random_element(S::Subspace{T}) where T <: Number = frombasis(S, randn(T, dim(S)))

################
### Constructors
################

"""
$(TYPEDSIGNATURES)

Random dimension-`d` subspace of dimension-`siz` vector space, on base field `T`.

```jldoctest
julia> S = random_subspace(ComplexF64, 2, (3, 4))
Subspace{Complex{Float64}} size (3, 4) dim 2
```
"""
function random_subspace(T::Type, d::Int, siz, tol=default_tol)
    if d < 0
        throw(ArgumentError("subspace dimension was negative: $d"))
    elseif d == 0
        return empty_subspace(T, siz, tol)
    else
        b = [ randn(T, siz) for i in 1:d ]
        return Subspace(b, tol=tol)
    end
end

"""
$(TYPEDSIGNATURES)

Random dimension-`d` subspace of `n × n` matrices on base field `T`, satisfying
`x ∈ S => x' ∈ S`.

```jldoctest
julia> S = random_hermitian_subspace(ComplexF64, 2, 3)
Subspace{Complex{Float64}} size (3, 3) dim 2

julia> S == S'
true
```
"""
function random_hermitian_subspace(T::Type, d::Int, n::Int, tol=default_tol)
    if d < 0
        throw(ArgumentError("subspace dimension was negative: $d"))
    elseif d == 0
        return empty_subspace(T, (n, n), tol)
    else
        b = [ randn(T, n, n) for i in 1:d ]
        b = [ x + x' for x in b ]
        return Subspace(b, tol=tol)
    end
end

"""
$(TYPEDSIGNATURES)

Create an empty subspace of dimension-`siz` vector space, on base field `T`.

```jldoctest
julia> S = empty_subspace(ComplexF64, (3, 4))
Subspace{Complex{Float64}} size (3, 4) dim 0
```
"""
empty_subspace(T::Type, siz::Tuple, tol=default_tol) = Subspace([zeros(T, siz)], tol=tol)

"""
$(TYPEDSIGNATURES)

Create an full subspace of dimension-`siz` vector space, on base field `T`.

```jldoctest
julia> S = full_subspace(ComplexF64, (3, 4))
Subspace{Complex{Float64}} size (3, 4) dim 12
```
"""
full_subspace(T::Type, siz::Tuple) = perp(empty_subspace(T, siz))

#############
### Hermitian
#############

function hermit_to_vec(M::AbstractArray{Complex{T}, 2}) where T
    n = size(M)[1]
    size(M)[2] == n || throw(ArgumentError("matrix was not square: $(size(M))"))
    v = zeros(T, n*n)
    k = 0
    for i in 1:n
        v[k+1:k+i] = real(M[1:i, i])
        k += i
        v[k+1:k+i-1] = imag(M[1:i-1, i])
        k += i-1
    end
    return v
end

function vec_to_hermit(v::AbstractArray{T, 1}, n::Integer) where T
    @assert n*n == size(v)[1]
    M = zeros(Complex{T}, n, n)
    k = 0
    for i in 1:n
        M[1:i, i] = v[k+1:k+i]
        M[i, 1:i] = v[k+1:k+i]
        k += i
        M[1:i-1, i] += v[k+1:k+i-1] * 1im
        M[i, 1:i-1] -= v[k+1:k+i-1] * 1im
        k += i-1
    end
    return Hermitian(M)
end

"""
$(TYPEDSIGNATURES)

Given a subspace satisfying `S == S'`, returns a basis in which each basis element is Hermitian.
"""
function hermitian_basis(S::Subspace{Complex{T}})::Array{Hermitian{Complex{T},Array{Complex{T},2}},1} where T
    if dim(S) == 0
        return []
    end
    n = size(S)[1]
    size(S)[2] == n || throw(ArgumentError("subspace size was not square: $(size(S))"))
    M = hcat(
        [ hermit_to_vec( x     + x'     ) for x in each_basis_element(S) ]...,
        [ hermit_to_vec((x*1im)+(x*1im)') for x in each_basis_element(S) ]...
    )
    hb = [ vec_to_hermit(x, n) for x in each_basis_element(Subspace(M)) ]
    @assert (Subspace(hb, tol=S.tol) == S) "Hermitian basis didn't equal original space"
    return hb
end

#########################
### Support for Convex.jl
#########################

"""
$(TYPEDSIGNATURES)

Returns the basis components of `x` in the basis `S.basis`.
"""
function tobasis(S::Subspace{<:Number, N}, x::Convex.AbstractExpr) where N
    shp = size(S)
    basis_mat = reshape(S.basis, prod(shp), dim(S))
    return basis_mat' * vec(x)
end

"""
$(TYPEDSIGNATURES)

Returns a vector from the given basis components in the basis `S.basis`.
"""
function frombasis(S::Subspace{<:Number, 1}, x::Convex.AbstractExpr)
    shp = size(S)
    basis_mat = reshape(S.basis, prod(shp), dim(S))
    return basis_mat * x
end

"""
$(TYPEDSIGNATURES)

Returns a vector from the given basis components in the basis `S.basis`.
"""
function frombasis(S::Subspace{<:Number, 2}, x::Convex.AbstractExpr)
    shp = size(S)
    basis_mat = reshape(S.basis, prod(shp), dim(S))
    return reshape(basis_mat * x, shp...)
end

"""
$(TYPEDSIGNATURES)

Projection of vector `x` onto suspace `S`.
"""
function projection(S::Subspace{<:Number, N}, x::Convex.AbstractExpr) where N
    return frombasis(S, tobasis(S, x))
end

"""
$(TYPEDSIGNATURES)

Constrains the variable `x` to be in the subspace `S`.
"""
function in(x::Convex.AbstractExpr, S::Subspace{<:Number, N}) where N
    return tobasis(perp(S), x) == 0
end

"""
$(TYPEDSIGNATURES)

Creates a `Convex.jl` variable ranging over the given subspace.
"""
function variable_in_space(S::Subspace{<:Complex{<:Real}, N}) where N
    x = Convex.ComplexVariable(dim(S))
    return frombasis(S, x)
end

"""
$(TYPEDSIGNATURES)

Creates a `Convex.jl` variable ranging over the given subspace.
"""
function variable_in_space(S::Subspace{<:Real, N}) where N
    x = Convex.Variable(dim(S))
    return frombasis(S, x)
end

end # module


text:
Base.@kwdef struct RepoSpec
    reponame::String # e.g "Lyceum/LyceumBase.jl"
    username::String
    token::String
    masterbranch::String = MASTER_BRANCH
    compatbranch::String = COMPAT_BRANCH
    gitconfig::Dict{String,String} = Dict{String,String}()
end

MaybeStr = Union{String,Nothing}
Base.@kwdef mutable struct CompatResult
    project_file::String
    # (name, new_compat)
    new::Vector{Tuple{String,String}} = Vector{Tuple{String,String}}()
    # (name, old_compat, new_compat)
    updated::Vector{Tuple{String,String,MaybeStr}} = Vector{Tuple{String,String,MaybeStr}}()
    # (name, old_compat)
    unchanged::Vector{Tuple{String,MaybeStr}} = Vector{Tuple{String,MaybeStr}}()
    julia_compat::MaybeStr = nothing
    manifest_updated::Bool = false
end

function project_updated(r::CompatResult)
    !isempty(r.new) || !isempty(r.updated) || r.julia_compat !== nothing
end
manifest_updated(r::CompatResult) = r.manifest_updated
updated(r::CompatResult) = project_updated(r) || manifest_updated(r)

function format_message(pkg_dir::AbstractString, rs::AbstractVector{CompatResult})
    msg = []
    for (i, r) in enumerate(rs)
        header = "Summary ($(relpath(r.project_file, pkg_dir)))"
        push!(msg, Header(header, 3))

        misc = List()
        r.manifest_updated && push!(misc.items, Paragraph("Updated Manifest"))
        if r.julia_compat !== nothing
            push!(misc.items, Paragraph("Updated compat entry for Julia: v$(r.julia_compat)"))
        end
        push!(misc.items, Paragraph("$(length(r.new)) new compat entries"))
        push!(misc.items, Paragraph("$(length(r.updated)) updated compat entries"))
        push!(misc.items, Paragraph("$(length(r.unchanged)) unchanged compat entries"))
        push!(msg, misc)

        if !isempty(r.new)
            push!(msg, Header("New Compat Entries", 4))
            titles = map(Bold, ["Package", "New Compat"])
            table = Markdown.Table([titles], [:l, :c])
            for (name, new) in r.new
                push!(table.rows, [name, new])
            end
            push!(msg, table)
        end

        if !isempty(r.updated)
            push!(msg, Header("Updated Compat Entries", 4))
            titles = map(Bold, ["Package", "Old Compat", "New Compat"])
            table = Markdown.Table([titles], [:l, :c, :c])
            for (name, old, new) in r.updated
                push!(table.rows, [name, old, new])
            end
            push!(msg, table)
        end
        i != length(rs) && push!(msg, Paragraph()) # line break
    end

    # magic bytes for identifying if an existing PR came from LyceumDevTools.Compat
    push!(msg, HorizontalRule())
    push!(msg, Paragraph("Last updated: $(now(UTC)) UTC"))
    push!(msg, Paragraph("Magic Bytes: $COMPAT_UUID"))

    Markdown.MD(msg...)
end