text:
module PSMV

# package code goes here

using Iterators

#simple_partition(n::Int, k::Int) = push!(map(first, partition(1:n,floor(Int,n/k))),n+1)
function simple_partition(n::Int,k::Int)
  rval = ones(Int,k+1)
  m,r = divrem(n,k)
  fill!(rval, m)
  rval[1] = 1
  for i = 2:k+1
    rval[i] += r .> (i-2)
  end
  cumsum!(rval,rval)
end

import Base.LinAlg.At_mul_B, Base.LinAlg.At_mul_B!, Base.LinAlg.A_mul_B!

"""
A type to encapsulate the multithreaded matrix-vector product operation.

Usage
-----
~~~~
A = sprandn(10^5,10^5,100/10^5); # create a fairly dense sparse matrix.
x = randn(10^5);
M = MultithreadedMatVec(A, simple_partition(A.n, Threads.nthreads()));
y = M.'*x;
norm(y - A.'*x)
~~~~
"""
mutable struct MultithreadedMatVec{T,I} <: AbstractArray{T,2}
  A::SparseMatrixCSC{T,I}
  regions::Vector{Int}

  function MultithreadedMatVec(A::SparseMatrixCSC{T,I}, regions::Vector{Int}) where {T,I}
    new{T,I}(A, regions)
  end
end

MultithreadedMatVec(A::SparseMatrixCSC) = MultithreadedMatVec(A,Threads.nthreads())
MultithreadedMatVec(A::SparseMatrixCSC, k::Int) = MultithreadedMatVec(A, simple_partition(A.n,k))

import Base.size
size(M::MultithreadedMatVec,inputs...) = size(M.A,inputs...)

import Base.getindex
getindex(M::MultithreadedMatVec,inputs...) = getindex(M.A, inputs...)

import Base.setindex!
setindex!(M::MultithreadedMatVec,inputs...) = setindex!(M.A, inputs...)

# Julia's built in multiplication operations are called with
# A_mul_B!
# Ac_mul_B!
# At_mul_B!
# which take in
# α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat
# and compute
# βC += αA B
# βC += αA^* B
# βC += αA^T B
# respectively
# look in base/sparse/linalg.jl
# for their implementations

""" Run the internal loop """
function internal_loop_transmult(C,B,nzv,rv,colptr,i1,i2,α::Number)
  for k=1:size(C,2)
    @inbounds for col=i1:i2
      tmp = zero(eltype(C))
      for j=colptr[col]:(colptr[col+1]-1)
        tmp += nzv[j]*B[rv[j],k]
      end
      C[col,k] += α*tmp
    end
  end
  return
end

"""
we are going to make these work with MuilthreadedMatVec types
"""
function At_mul_B!(α::Number, M::MultithreadedMatVec, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)
  M.A.m == size(B,1) || throw(DimensionMismatch())
  M.A.n == size(C,1) || throw(DimensionMismatch())
  size(B,2) == size(C,2) || throw(DimensionMismatch())
  nzv = M.A.nzval
  rv = M.A.rowval
  colptr = M.A.colptr
  if β != 1
    β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
  end
  # this is the parallel construction
  Threads.@threads for t=1:(length(M.regions)-1)
    internal_loop_transmult(C,B,nzv,rv,colptr,M.regions[t],M.regions[t+1]-1,α)
  end
  C
end

function At_mul_B(M::MultithreadedMatVec{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx}
  T = promote_type(TA,Tx)
  At_mul_B!(one(T), M, x, zero(T), similar(x, T, M.A.n))
end

module MyAtomics

lt = "double" # llvmtype
ilt = "i64" # llvmtype of int representation

import Base.Sys: ARCH, WORD_SIZE

@inline function atomic_cas!(x::Array{Float64}, ind::Int, oldval::Float64, newval::Float64)
    xptr = Base.unsafe_convert(Ptr{Float64}, x)
    xptr += 8*(ind-1)
    # on new versions of Julia, this should be
    # %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
    Base.llvmcall( """%iptr = bitcast $lt* %0 to $ilt*
                      %icmp = bitcast $lt %1 to $ilt
                      %inew = bitcast $lt %2 to $ilt
                      %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                      %irv = extractvalue { $ilt, i1 } %irs, 0
                      %rv = bitcast $ilt %irv to $lt
                      ret $lt %rv
                """,
    Float64,  # return type
    Tuple{Ptr{Float64}, Float64, Float64},  # argument types
    xptr, oldval, newval # arguments
    )
end

@inline function atomic_add!(x::Array{Float64}, ind::Int, val::Float64)
  IT = Int64
  xptr = Base.unsafe_convert(Ptr{Float64}, x)
  xptr += 8*(ind-1)

  while true
    oldval = x[ind]
    newval = oldval + val

    #old2 = atomic_cas!(x, ind, oldval, newval)
    # inline this call and optimize out some stuff
    old2 = Base.llvmcall( """%iptr = bitcast $lt* %0 to $ilt*
                      %icmp = bitcast $lt %1 to $ilt
                      %inew = bitcast $lt %2 to $ilt
                      %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                      %irv = extractvalue { $ilt, i1 } %irs, 0
                      %rv = bitcast $ilt %irv to $lt
                      ret $lt %rv
                """,
    Float64,  # return type
    Tuple{Ptr{Float64}, Float64, Float64},  # argument types
    xptr, oldval, newval # arguments
    )

    if reinterpret(IT,oldval) == reinterpret(IT,old2)
      return newval
    end
  end
end
end

function internal_loop_mult(C,B,nzv,rv,colptr,i1,i2,α::Number)
  for k=1:size(C,2)
    koffset = size(C,1)*(k-1)
    @inbounds for col=i1:i2
      αxj = α*B[col,k]
      for j=colptr[col]:(colptr[col+1]-1)
        # need to be done atomically...
        #Threads.atomic_add!(C[rv[j],k], nzv[j]*αxj)
        MyAtomics.atomic_add!(C,koffset+rv[j],nzv[j]*αxj)
      end
    end
  end
  return
end

function A_mul_B!(α::Number, M::MultithreadedMatVec, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)
  M.A.n == size(B,1) || throw(DimensionMismatch())
  M.A.m == size(C,1) || throw(DimensionMismatch())
  size(B,2) == size(C,2) || throw(DimensionMismatch())
  nzv = M.A.nzval
  rv = M.A.rowval
  colptr = M.A.colptr
  if β != 1
    β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
  end
  # this is the parallel construction
  Threads.@threads for t=1:(length(M.regions)-1)
    internal_loop_mult(C,B,nzv,rv,colptr,M.regions[t],M.regions[t+1]-1,α)
  end
  C
end

import Base.eltype
eltype(M::MultithreadedMatVec{T,I}) where {T,I} = T

function A_mul_B!(C::StridedVecOrMat, M::MultithreadedMatVec, B::StridedVecOrMat)
  T = promote_type(eltype(M), eltype(B))
  A_mul_B!(one(T), M, B, zero(T), C)
end

"""

Multithreaded Sparse Matrix Vector Multiply using Transpose matrix for reduced communication.

"""

mutable struct MultithreadedTransMatVec{T,I} <: AbstractArray{T,2}
  A::SparseMatrixCSC{T,I}
  Q::SparseMatrixCSC{T,I}
  regions::Vector{Int}

  function MultithreadedTransMatVec(A::SparseMatrixCSC{T,I}, Q::SparseMatrixCSC, regions::Vector{Int}) where {T,I}
    new{T,I}(A, Q, regions)
  end

end

function MultithreadedTransMatVec(A::SparseMatrixCSC, Q::SparseMatrixCSC, k::Int)
  MultithreadedTransMatVec(A, Q, simple_partition(A.n,k))
end

function MultithreadedTransMatVec(A::SparseMatrixCSC, Q::SparseMatrixCSC)
  MultithreadedTransMatVec(A,Q,Threads.nthreads())
end

size(M::MultithreadedTransMatVec,inputs...) = size(M.A,inputs...)

import Base.getindex
getindex(C::MultithreadedTransMatVec, inputs...) = getindex(C.A, inputs...)

#import Base.linearindexing
#linearindexing(C::MultithreadedTransMatVec) = Base.LinearSlow()
import Base.setindex!
function setindex!(C::MultithreadedTransMatVec{Tv, Ti}, _v, _i::Int, _j::Int) where {Tv, Ti<:Integer}
  v = convert(Tv, _v)
  i = convert(Ti, _i)
  j = convert(Ti, _j)
  if !((1 <= i <= C.A.m) & (1 <= j <= C.A.n))
      throw(BoundsError(C.A, (i,j)))
  end

  setindex!(C.A, v, i, j)
  setindex!(C.Q, v, j, i)
end

function internal_loop_mult_opt(C,B,nzv,rv,colptr,i1,i2,α::Number)
  k=1
  #for k=1:size(C,1)
    @inbounds for col=i1:i2
      tmp = zero(eltype(C))
      for j=colptr[col]:(colptr[col+1]-1)
        tmp += nzv[j]*B[rv[j],k]
      end
      C[col,k] += α*tmp
    end
  #end
  return
end

function A_mul_B!(α::Number, M::MultithreadedTransMatVec, B::StridedVector, β::Number, C::StridedVector)
  M.A.n == size(B,1) || throw(DimensionMismatch())
  M.A.m == size(C,1) || throw(DimensionMismatch())
  size(B,2) == size(C,2) || throw(DimensionMismatch())
  nzv = M.Q.nzval
  rv = M.Q.rowval
  colptr = M.Q.colptr

  if β != 1
    β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
  end

  # this is the parallel construction
  Threads.@threads for t=1:(length(M.regions)-1)
    internal_loop_mult_opt(C,B,nzv,rv,colptr,M.regions[t],M.regions[t+1]-1,α)
  end
  C
end

import Base.eltype
eltype(M::MultithreadedTransMatVec{T,I}) where {T,I} = T

function A_mul_B!(C::StridedVector, M::MultithreadedTransMatVec, B::StridedVector)
  T = promote_type(eltype(M), eltype(B))
  A_mul_B!(one(T), M, B, zero(T), C)
end

function test_perf(n::Int)
  println("Constructing Sparse Matrix")
  @time A = sprandn(n,n,100/n); # create a fairly dense sparse matrix.
  println("Constructing Transpose")
  @time Q = A'

  x = randn(n);
  M = MultithreadedMatVec(A, simple_partition(n, 1));
  @show norm(M*x - A*x)
  println("Ax")
  @time A*x;
  @time A*x;
  @time A*x;

  println("Mx")
  @time M*x;
  @time M*x;
  @time M*x;

  M = MultithreadedMatVec(A, simple_partition(n, 2));
  println("Mx - 2")
  @show norm(M*x - A*x)
  @time M*x;
  @time M*x;
  @time M*x;

  M = MultithreadedMatVec(A, simple_partition(n, 4));
  println("Mx - 4")
  @show norm(M*x - A*x)
  @time M*x;
  @time M*x;
  @time M*x;

  M = MultithreadedMatVec(A, simple_partition(n, 6));
  println("Mx - 6")
  @show norm(M*x - A*x)
  @time M*x;
  @time M*x;
  @time M*x;

  M = MultithreadedMatVec(A, simple_partition(n, 12));
  println("Mx - 12")
  @show norm(M*x - A*x)
  @time M*x;
  @time M*x;
  @time M*x;

  println("My benchmarks")
  @time C = MultithreadedTransMatVec(A, Q, 1)
  @time @show norm(C*x - A*x)
  println("Ax")
  @time A*x
  @time A*x
  @time A*x

  println("Cx")
  @time C*x
  @time C*x
  @time C*x

  C = MultithreadedTransMatVec(A, Q, 2)
  println("Cx - 2")
  @time C*x
  @time C*x
  @time C*x

  C = MultithreadedTransMatVec(A, Q, 4)
  println("Cx - 4")
  @time C*x
  @time C*x
  @time C*x

  C = MultithreadedTransMatVec(A, Q, 6)
  println("Cx - 6")
  @time C*x
  @time C*x
  @time C*x

  C = MultithreadedTransMatVec(A, Q, 12)
  println("Cx - 12")
  @time C*x
  @time C*x
  @time C*x

  C = MultithreadedTransMatVec(A, Q)
  println("Cx - ?")
  @time C*x
  @time C*x
  @time C*x

  return 1
end

end # module

max_stars_repo_name: camruggles/PSVM

max_issues_repo_name: camruggles/PSVM

max_forks_repo_name: camruggles/PSVM


text:
# Run this in a separate process (with sudo privileges) to gather performance counter data.
function runserver()
    # Create a NamedPipe to listen for commands from the other process.
    server = listen(PIPEPATH)

    # Make readable and writeable by non-sudo processes
    chmod(PIPEPATH, 0x666)

    # Set up
    while true
        sock = accept(server)

        # Setup defaults
        sampletime = 1
        filepath = "test.jls"
        measurements = NamedTuple()
        params = NamedTuple()

        # Listen on the socket for instructions.
        while isopen(sock)
            cmd = readline(sock)

            # Shutdown command
            if cmd == "shutdown"
                return nothing

            # Set filepath
            elseif startswith(cmd, "filepath")
                filepath = last(split(cmd))
                println("filepath = $filepath")

            # Update `measurements` tuple.
            elseif cmd == "measurements"
                measurements = getmeasurements()
                println("measurements = $measurements")

            # Update `params` tuple
            elseif cmd == "params"
                params = getparams()
                println("params = $params")

            # Update `sampletime`.
            elseif startswith(cmd, "sampletime")
                sampletime = parse(Float64, last(split(cmd)))
                println("sampletime = $sampletime")

            # start sampling.
            elseif cmd == "start"
                GC.gc()
                sample(sock, sampletime, filepath, params, measurements)
                GC.gc()

            # who knows what happened.
            else
                println("Unknown Command: $cmd")
            end
        end
    end
end

# Deserialize the contents at `path`. It should be a function that can be called to get
# the requested measurements `NamedTuple`.
#
# The reason this has to be a function instead of just the NamedTuple itself is because the
# NamedTuple can contain some non-serializable structs.
function getmeasurements()
    x = deserialize(TRANSFERPATH)
    nt = x()::NamedTuple
    return nt
end

getparams() = deserialize(PARAMPATH)::NamedTuple

# Sampline routine.
function sample(sock, sampletime, filepath, params, measurements)
    # Sample time regularly.
    sampler = SystemSnoop.SmartSample(Second(sampletime))

    # Normally, SystemSnoop requires a PID for what it's snooping. I no longer like this API,
    # but I guess I kind of have to deal with it for now until I get around to changing it.
    #
    # The PID is not needed for this snooping since we're monitoring the entire system,
    # so just pass the current PID
    local data
    @sync begin
        # Spawn a task to sample the buffer and notify when a `stop` command is reached
        canexit = false
        @async begin
            while true
                ln = readline(sock)
                if ln == "stop"
                    canexit = true
                    break
                else
                    println("Unhandled command: $ln")
                end
            end
        end

        # Measuring loop.
        trace, data = SystemSnoop.snoop(measurements) do snooper
            while true
                # Sleep until it's time to sample.
                sleep(sampler)
                SystemSnoop.measure!(snooper)

                # See if we have any communication from the master process. If so, pull it
                # out and make sure it's a `stop` command.
                if canexit
                    println("Stopping Sampling")
                    return SystemSnoop.postprocess(snooper)
                end
            end
        end
    end

    # Serialize the data to the given filepath.
    # If the path is already valid, than first deserialize it, add our data to it, then
    # reserialize it.
    dir = dirname(filepath)
    !ispath(dir) && mkpath(dir)
    if ispath(filepath)
        x = deserialize(filepath)
    else
        x = Slab()
    end

    # Only save the post-processed data for now.
    dict = Dict(k => v for (k,v) in pairs(data))

    # Merge this data in with the rest
    x[params] = dict
    serialize(filepath, x)
    return true
end

ext: jl

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
struct FixedKernel{K<:Kernel} <: Kernel
    kernel::K
    free::Vector{Int} # vector of *free* parameters
end

@deprecate FixedKern FixedKernel

get_params(k::FixedKernel) = get_params(k.kernel)[k.free]
get_param_names(k::FixedKernel) = get_param_names(k.kernel)[k.free]
function set_params!(k::FixedKernel, hyp)
    p = get_params(k.kernel)
    p[k.free] = hyp
    set_params!(k.kernel, p)
end
num_params(k::FixedKernel) = length(k.free)

# convenience functions to fix a parameter
function fix(k::Kernel, par::Symbol)
    npars = num_params(k)
    free = collect(1:npars)
    names = get_param_names(k)
    tofix = findfirst(==(par), names)
    tofix == nothing || deleteat!(free, tofix)
    return FixedKernel(k, free)
end

function fix(k::FixedKernel, par::Symbol)
    free = k.free
    names = get_param_names(k)
    tofix = findfirst(==(par), names)
    tofix == nothing || deleteat!(free, tofix)
    return FixedKernel(k.kernel, free)
end
function fix(k::Kernel)
    return FixedKernel(k, Int[])
end

# convenience functions to free a parameter
function free(k::FixedKernel)
    return k.kernel
end
function free(k::FixedKernel, par::Symbol)
    all_names = get_param_names(k.kernel)
    ipar = findfirst(==(par), all_names)
    if ipar == nothing || ipar ∉ k.free
        free = k.free
    else
        free = sort(unique([k.free; ipar]))
    end

    return FixedKernel(k.kernel, free)
end

function grad_slice!(dK::MatF64, k::FixedKernel, X::MatF64, data::KernelData, p::Int)
    return grad_slice!(dK, k.kernel, X, data, k.free[p])
end
function grad_slice!(dK::MatF64, k::FixedKernel, X::MatF64, data::EmptyData, p::Int)
    return grad_slice!(dK, k.kernel, X, data, k.free[p])
end
@inline function dKij_dθp(fk::FixedKernel,X::MatF64,i::Int,j::Int,p::Int,dim::Int)
    return dKij_dθp(fk.kernel, X, i, j, fk.free[p], dim)
end
@inline function dKij_dθp(fk::FixedKernel,X::MatF64,data::KernelData,i::Int,j::Int,p::Int,dim::Int)
    return dKij_dθp(fk.kernel, X, data, i, j, fk.free[p], dim)
end

# delegate everything else to the wrapped kernel
# (is there a better way to do this?)
cov_ij(fk::FixedKernel, X::MatF64, i::Int, j::Int, dim::Int) = cov_ij(fk.kernel, X, i, j, dim)
cov_ij(fk::FixedKernel, X::MatF64, data::KernelData, i::Int, j::Int, dim::Int) = cov_ij(fk.kernel, X, data, i, j, dim)
cov_ij(fk::FixedKernel, X::MatF64, data::EmptyData, i::Int, j::Int, dim::Int) = cov_ij(fk, X, i, j, dim)
Statistics.cov(fk::FixedKernel, x::VecF64, y::VecF64) = cov(fk.kernel, x, y)
KernelData(fk::FixedKernel, args...) = KernelData(fk.kernel, args...)
KernelData(fk::FixedKernel, X::MatF64) = KernelData(fk.kernel, X)
kernel_data_key(fk::FixedKernel, args...) = kernel_data_key(fk.kernel, args...)
kernel_data_key(fk::FixedKernel, X::MatF64) = kernel_data_key(fk.kernel, X)

##########
# Priors #
##########

function get_priors(k::FixedKernel)
    free_priors = get_priors(k.kernel)
    if isempty(free_priors) return []
    else
        return free_priors[k.free]
    end
end

function set_priors!(k::FixedKernel, priors)
    p = get_priors(k.kernel)
    p[k.free] = priors
    set_priors!(k.kernel, p)
end

function prior_logpdf(k::FixedKernel)
    return 0.0
end

function prior_gradlogpdf(k::FixedKernel)
    return zeros(num_params(k))
end


text:
function setup_mplhep()
    @eval begin
        using Conda
        Conda.add_channel("conda-forge")
        Conda.add("mplhep")
    end
end

using PyCall, PyPlot

@pyimport mplhep as hep
hep.set_style(hep.style.CMS)
hep.set_style(Dict("font.sans-serif" => "DejaVu Sans"))
matplotlib.rc("text", usetex=true)
#matplotlib.rc("text", usetex=false)

@pyimport matplotlib.patheffects as path_effects

function annotate_cms(ax; sel_pos=(.05, .05), kw...)
    hep.cms.text("Private Work"; fontsize=20, ax)
    hep.cms.lumitext(L"41.5 $\mathrm{fb^{-1}}$"; fontsize=16, ax)
    ax.text(sel_pos..., L"1 Lepton, 6 Jets, $\ge$3 Btags";
        transform=ax.transAxes, fontsize=16, kw...)
end


text:
using DataFrames
using Test

m = test_page_model()
include("../src/components/SLRDamages.jl")

slrdamages = addslrdamages(m)

set_param!(m, :SLRDamages, :y_year, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
set_param!(m, :SLRDamages, :atl_adjustedtolerablelevelofsealevelrise, readpagedata(m, "test/validationdata/atl_adjustedtolerablelevelofsealevelrise.csv"))
set_param!(m, :SLRDamages, :imp_actualreductionSLR, readpagedata(m, "test/validationdata/imp_actualreductionSLR.csv"))
set_param!(m, :SLRDamages, :y_year_0, 2008.)
set_param!(m, :SLRDamages, :s_sealevel, readpagedata(m, "test/validationdata/s_sealevelrise.csv"))
set_param!(m, :SLRDamages, :cons_percap_consumption, readpagedata(m, "test/validationdata/cons_percap_consumption.csv"))
set_param!(m, :SLRDamages, :tct_per_cap_totalcostspercap, readpagedata(m, "test/validationdata/tct_per_cap_totalcostspercap.csv"))
set_param!(m, :SLRDamages, :act_percap_adaptationcosts, readpagedata(m, "test/validationdata/act_percap_adaptationcosts.csv"))
set_param!(m, :SLRDamages, :isatg_impactfxnsaturation, 28.333333333333336)
set_param!(m, :SLRDamages, :wincf_weightsfactor, readpagedata(m, "data/wincf_weightsfactor.csv"))

##running Model
run(m)

cons_percap_aftercosts = m[:SLRDamages, :cons_percap_aftercosts]
cons_percap_aftercosts_compare = readpagedata(m, "test/validationdata/cons_percap_aftercosts.csv")
@test ones(10, 8) ≈ cons_percap_aftercosts ./ cons_percap_aftercosts_compare rtol=.001

i_regionalimpactSLR = m[:SLRDamages, :i_regionalimpactSLR]
i_regionalimpactSLR_compare = readpagedata(m, "test/validationdata/i_regionalimpact_SLRise.csv")
@test i_regionalimpactSLR ≈ i_regionalimpactSLR_compare rtol=.001

iref_ImpactatReferenceGDPperCapSLR = m[:SLRDamages, :iref_ImpactatReferenceGDPperCapSLR]
iref_ImpactatReferenceGDPperCapSLR_compare = readpagedata(m, "test/validationdata/iref_ImpactatReferenceGDPperCap_sea.csv")
@test iref_ImpactatReferenceGDPperCapSLR ≈ iref_ImpactatReferenceGDPperCapSLR_compare rtol=.0001

rcons_per_cap_SLRRemainConsumption = m[:SLRDamages, :rcons_per_cap_SLRRemainConsumption]
rcons_per_cap_SLRRemainConsumption_compare = readpagedata(m, "test/validationdata/rcons_per_cap_SLRRemainConsumption.csv")
@test rcons_per_cap_SLRRemainConsumption ≈ rcons_per_cap_SLRRemainConsumption_compare rtol=.01

rgdp_per_cap_SLRRemainGDP = m[:SLRDamages, :rgdp_per_cap_SLRRemainGDP]
rgdp_per_cap_SLRRemainGDP_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_SLRRemainGDP.csv")
@test rgdp_per_cap_SLRRemainGDP ≈ rgdp_per_cap_SLRRemainGDP_compare rtol=.01

igdp_ImpactatActualGDPperCapSLR = m[:SLRDamages, :igdp_ImpactatActualGDPperCapSLR]
igdp_ImpactatActualGDPperCapSLR_compare = readpagedata(m, "test/validationdata/igdp_ImpactatActualGDPperCap_sea.csv")
@test igdp_ImpactatActualGDPperCapSLR ≈ igdp_ImpactatActualGDPperCapSLR_compare rtol=.01

isat_ImpactinclSaturationandAdaptationSLR = m[:SLRDamages, :isat_ImpactinclSaturationandAdaptationSLR]
isat_ImpactinclSaturationandAdaptationSLR_compare = readpagedata(m, "test/validationdata/isat_ImpactinclSaturationandAdaptation_SLRise.csv")
@test isat_ImpactinclSaturationandAdaptationSLR ≈ isat_ImpactinclSaturationandAdaptationSLR_compare rtol=.001


text:
function createdgnodes(p,t,f,curvedboundary,curvedboundaryexpr,porder)

nve,ne = size(t);
nd=size(p,1);

elemtype = 0;
if (nd==2) && (nve==4)
    elemtype=1;
end
if (nd==3) && (nve==8)
    elemtype=1;
end

plocal,~,~,~,perm = masternodes(porder,nd,elemtype);
npl = size(plocal,1);

philocal = zeros(Float64, npl, nve);

if nd==1
    xi  = plocal[:,1];
    philocal[:,1] = 1.0 .- xi;
    philocal[:,2] = xi;
elseif nd==2 && nve==3 # tri
    xi  = plocal[:,1];
    eta = plocal[:,2];
    philocal[:,1] = 1.0 .- xi .- eta;
    philocal[:,2] = xi;
    philocal[:,3] = eta;
elseif nd==2 && nve==4 # quad
    xi  = plocal[:,1];
    eta = plocal[:,2];
    philocal[:,1] = (1.0 .- xi).*(1.0 .- eta);
    philocal[:,2] = xi.*(1.0 .- eta);
    philocal[:,3] = xi.*eta;
    philocal[:,4] = (1.0 .- xi).*eta;
elseif nd==3 && nve==4 # tet
    xi   = plocal[:,1];
    eta  = plocal[:,2];
    zeta = plocal[:,3];
    philocal[:,1] = 1.0 .- xi .- eta .- zeta;
    philocal[:,2] = xi;
    philocal[:,3] = eta;
    philocal[:,4] = zeta;
elseif nd==3 && nve==8 # hex
    xi   = plocal[:,1];
    eta  = plocal[:,2];
    zeta = plocal[:,3];
    philocal[:,1] = (1.0 .-xi).*(1.0 .-eta).*(1.0 .-zeta);
    philocal[:,2] = xi.*(1.0 .-eta).*(1.0 .-zeta);
    philocal[:,3] = xi.*eta.*(1.0 .-zeta);
    philocal[:,4] = (1.0 .-xi).*eta.*(1.0 .-zeta);
    philocal[:,5] = (1.0 .-xi).*(1.0 .-eta).*(zeta);
    philocal[:,6] = xi.*(1.0 .-eta).*(zeta);
    philocal[:,7] = xi.*eta.*(zeta);
    philocal[:,8] = (1.0 .-xi).*eta.*(zeta);
end

# Allocate nodes
dgnodes = zeros(Float64,npl,ne,nd);
for dim=1:nd
  for node=1:nve
    dp = reshape(philocal[:,node],npl,1)*reshape(p[dim,t[node,:]],1,ne);
    dgnodes[:,:,dim] = dgnodes[:,:,dim] + reshape(dp,npl,ne);
  end
end

dgnodes = permutedims(dgnodes,[1, 3, 2]);

if length(curvedboundaryexpr)>0 && porder>1 && maximum(curvedboundary[:])>0
    print("Project dgnodes onto the curved boundaries...\n");
    fd = curvedboundaryexpr;
    nd = size(p,1);
    nfe,ne = size(f);
    if nd == 2
        for i = 1:ne
            for j = 1:nfe
                if f[j,i]!=0 # boundary element
                    k = abs(f[j,i]); # get boundary index
                    if curvedboundary[k]==1 # if this boundary is curved
                        p = dgnodes[perm[:,j],:,i];
                        deps = 1e-8*maximum(maximum(p,dims=1)-minimum(p,dims=1));
                        d = fd[k](p');
                        dgradx = (fd[k]([p[:,1].+deps p[:,2]]') - d)/deps;
                        dgrady = (fd[k]([p[:,1] p[:,2].+deps]') - d)/deps;
                        dgrad2 = dgradx.*dgradx + dgrady.*dgrady;
                        dgrad2[dgrad2.==0] .= 1;
                        p = p - [d.*dgradx./dgrad2 d.*dgrady./dgrad2];
                        dgnodes[perm[:,j],:,i] = p;
                    end
                end
            end
        end
    elseif nd==3
        for i = 1:ne
            for j = 1:nfe
                if f[j,i]!=0 # boundary element
                    k = abs(f[j,i]); # get boundary index
                    if curvedboundary[k]==1 # if this boundary is curved
                        p = dgnodes[perm[:,j],:,i];
                        deps = 1e-8*maximum(maximum(p,dim=1)-minimum(p,dim=1));
                        d = fd[k](p');
                        dgradx = (fd[k]([p[:,1].+deps p[:,2] p[:,3]]')-d)/deps;
                        dgrady = (fd[k]([p[:,1] p[:,2].+deps p[:,3]]')-d)/deps;
                        dgradz = (fd[k]([p[:,1] p[:,2] p[:,3].+deps]')-d)/deps;
                        dgrad2 = dgradx.*dgradx+dgrady.*dgrady+dgradz.*dgradz;
                        dgrad2[dgrad2.==0] .= 1;
                        p = p-[d.*dgradx./dgrad2 d.*dgrady./dgrad2 d.*dgradz./dgrad2];
                        dgnodes[perm[:,j],:,i] = p;
                    end
                end
            end
        end
    end
end

return dgnodes

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
## Get the optimal position of nodes given the list of nodes and plates

function getpos(N, plateList)

    nPlate = length(plateList)

    nNodes = 0
    for i in 1:length(N)
        nNodes += length(N[i])
    end

    nodeList = collect(1:nNodes)

    plateLevelList = []

    for i in 1: nPlate
        h = []
        for k in 1:length(plateList[i])
            for j in 1: length(N)
                if plateList[i][k] ∈ N[j]
                    push!(h,j)
                end
            end
        end
        push!(plateLevelList, h)
    end

    rangeList=[]
    for i in 1: nPlate
        h = []
        push!(h,minimum(plateLevelList[i]), maximum(plateLevelList[i]))
        push!(rangeList, h)
    end

    #model = Model(with_optimizer(Mosek.Optimizer, QUIET=false, INTPNT_CO_TOL_DFEAS=1e-7))
    model = Model(GLPK.Optimizer)

    nX = nNodes^2 - nNodes

    nA = 0
    for j in 1:nPlate
        for i in rangeList[j][1] : rangeList[j][2]
            nA += length(setdiff(N[i], plateList[j]))
        end
    end

    @variable(model, A[1:nA] >= 0)
    @variable(model, B[1:nA], Bin)
    @variable(model, X[1:nNodes,1:nNodes], Int)
    #@variable(model, Y[1:14,1:2], Bin)

    @objective(model, Min, sum(A[i] for i in 1:nA))

    w = 1
    for j in 1:nPlate
        for i in rangeList[j][1] : rangeList[j][2]
            for k in setdiff(N[i], plateList[j])
                @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                    A[w] + length(plateList[j]))
                @constraint(model, -sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                    A[w] + length(plateList[j]))
                @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) >=
                    length(plateList[j]) - A[w] - 1000*B[w])
                @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                    1000 - length(plateList[j]) + A[w] - 1000*B[w])
                w = w + 1
            end
        end
    end

    for i in 1:(nNodes-1)
        for j in (i+1):nNodes
            @constraint(model, X[i,j] + X[j,i] == 0)
        end
    end

    for i in 1:nNodes
            @constraint(model, X[i,i] == 0)
    end

    for i in 1:nNodes
        for j ∈ nodeList[setdiff(1:nNodes, i)]
            for k ∈ nodeList[setdiff(1:nNodes, [i,j])]
                @constraint(model, X[i,j] + X[j,k] - X[i,k]  <= 1)
            end
        end
    end

    @constraint(model, -1 .<= X .<= 1)

    optimize!(model)
    @show value.(X)
    #return(nNodes)
end


text:
typealias MPIDatatype Union{Char,
                            Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64,
                            UInt64,
                            Float32, Float64, Complex64, Complex128}
typealias MPIBuffertype{T} Union{Ptr{T}, Array{T}, Ref{T}}

if VERSION >= v"0.5.0-"
    # TODO: Use Compat for this
    fieldoffsets{T}(::Type{T}) = Int[fieldoffset(T, i) for i in 1:nfields(T)]
end

# Define a function mpitype(T) that returns the MPI datatype code for
# a given type T. In the case the the type does not exist, it is created and
# then returned. The dictonary is defined in __init__ so the module can be
# precompiled

# accessor and creation function for getting MPI datatypes
function mpitype{T}(::Type{T})

    if haskey(mpitype_dict, T)  # if the datatype already exists
      return mpitype_dict[T]
    end

    if !isbits(T)
        throw(ArgumentError("Type must be isbits()"))
    end

    # get the data from the type
    fieldtypes = T.types
    offsets = fieldoffsets(T)
    nfields = Cint(length(fieldtypes))

    # put data in MPI format
    blocklengths = ones(Cint, nfields)
    displacements = zeros(Cptrdiff_t, nfields)  # size_t == MPI_Aint ?
    types = zeros(Cint, nfields)
    for i=1:nfields
        displacements[i] = offsets[i]
        # create an MPI_Datatype for the current field if it does not exist yet
        types[i] = mpitype(fieldtypes[i])
    end

    # create the datatype
    newtype_ref = Ref{Cint}()
    flag = Ref{Cint}()
    ccall(MPI_TYPE_CREATE_STRUCT, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cptrdiff_t},
          Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), &nfields, blocklengths,
          displacements, types, newtype_ref, flag)

    if flag[] != 0
        throw(ErrorException("MPI_Type_create_struct returned non-zero exit status"))
    end

    # commit the datatatype
    flag2 = Ref{Cint}()
    ccall(MPI_TYPE_COMMIT, Void, (Ptr{Cint}, Ptr{Cint}), newtype_ref, flag2)

    if flag2[] != 0
        throw(ErrorException("MPI_Type_commit returned non-zero exit status"))
    end

    # add it to the dictonary of known types
    mpitype_dict[T] = newtype_ref[]
    mpitype_dict_inverse[newtype_ref[]] = T

    return mpitype_dict[T]
end

type Comm
    val::Cint
    Comm(val::Integer) = new(val)
end
const COMM_NULL  = Comm(MPI_COMM_NULL)
const COMM_SELF  = Comm(MPI_COMM_SELF)
const COMM_WORLD = Comm(MPI_COMM_WORLD)

type Op
    val::Cint
end
const OP_NULL = Op(MPI_OP_NULL)
const BAND    = Op(MPI_BAND)
const BOR     = Op(MPI_BOR)
const BXOR    = Op(MPI_BXOR)
const LAND    = Op(MPI_LAND)
const LOR     = Op(MPI_LOR)
const LXOR    = Op(MPI_LXOR)
const MAX     = Op(MPI_MAX)
const MIN     = Op(MPI_MIN)
const PROD    = Op(MPI_PROD)
const SUM     = Op(MPI_SUM)

type Request
    val::Cint
    buffer
end
const REQUEST_NULL = Request(MPI_REQUEST_NULL, nothing)

type Status
    val::Array{Cint,1}
    Status() = new(Array(Cint, MPI_STATUS_SIZE))
end
Get_error(stat::Status) = Int(stat.val[MPI_ERROR])
Get_source(stat::Status) = Int(stat.val[MPI_SOURCE])
Get_tag(stat::Status) = Int(stat.val[MPI_TAG])

const ANY_SOURCE = Int(MPI_ANY_SOURCE)
const ANY_TAG    = Int(MPI_ANY_TAG)
const TAG_UB     = Int(MPI_TAG_UB)
const UNDEFINED  = Int(MPI_UNDEFINED)

function serialize(x)
    s = IOBuffer()
    Base.serialize(s, x)
    Base.takebuf_array(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Base.deserialize(s)
end

# Administrative functions

function Init()
    ccall(MPI_INIT, Void, (Ptr{Cint},), &0)
end

function Finalize()
    ccall(MPI_FINALIZE, Void, (Ptr{Cint},), &0)
end

function Abort(comm::Comm, errcode::Integer)
    ccall(MPI_ABORT, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &comm.val, &errcode, &0)
end

function Initialized()
    flag = Ref{Cint}()
    ccall(MPI_INITIALIZED, Void, (Ptr{Cint}, Ptr{Cint}), flag, &0)
    flag[] != 0
end

function Finalized()
    flag = Ref{Cint}()
    ccall(MPI_FINALIZED, Void, (Ptr{Cint}, Ptr{Cint}), flag, &0)
    flag[] != 0
end

function Comm_dup(comm::Comm)
    newcomm = Ref{Cint}()
    ccall(MPI_COMM_DUP, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &comm.val, newcomm, &0)
    MPI.Comm(newcomm[])
end

function Comm_free(comm::Comm)
    ccall(MPI_COMM_FREE, Void, (Ptr{Cint}, Ptr{Cint}), &comm.val, &0)
end

function Comm_rank(comm::Comm)
    rank = Ref{Cint}()
    ccall(MPI_COMM_RANK, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &comm.val, rank, &0)
    Int(rank[])
end

function Comm_size(comm::Comm)
    size = Ref{Cint}()
    ccall(MPI_COMM_SIZE, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &comm.val, size, &0)
    Int(size[])
end

function type_create(T::DataType)
    if !isbits(T)
        throw(ArgumentError("Type must be isbits()"))
    end

    if haskey(mpitype_dict, T)  # if the datatype already exists
        return nothing
    end

    # get the data from the type
    fieldtypes = T.types
    nfields = Cint(length(fieldtypes))
    if VERSION < v"0.5.0-dev+2285"
        offsets = fieldoffsets(T)
    else
        offsets = map(idx->fieldoffset(T, idx), 1:nfields)
    end

    # put data in MPI format
    blocklengths = ones(Cint, nfields)
    displacements = zeros(Cptrdiff_t, nfields)  # size_t == MPI_Aint ?
    types = zeros(Cint, nfields)
    for i=1:nfields
        displacements[i] = offsets[i]
        # create an MPI_Datatype for the current field if it does not exist yet
        if !haskey(mpitype_dict, fieldtypes[i])
            type_create(fieldtypes[i])
        end
        types[i] = mpitype(fieldtypes[i])
    end

    # create the datatype
    newtype_ref = Ref{Cint}()
    flag = Ref{Cint}()
    ccall(MPI_TYPE_CREATE_STRUCT, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cptrdiff_t},
          Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), &nfields, blocklengths, displacements,
          types, newtype_ref, flag)

    if flag[] != 0
        throw(ErrorException("MPI_Type_create_struct returned non-zero exit status"))
    end

    # commit the datatatype
    flag2 = Ref{Cint}()
    ccall(MPI_TYPE_COMMIT, Void, (Ptr{Cint}, Ptr{Cint}), newtype_ref, flag2)

    if flag2[] != 0
        throw(ErrorException("MPI_Type_commit returned non-zero exit status"))
    end

    # add it to the dictonary of known types
    mpitype_dict[T] = newtype_ref[]

    return nothing
end

# Point-to-point communication

function Probe(src::Integer, tag::Integer, comm::Comm)
    stat = Status()
    ccall(MPI_PROBE, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &src, &tag, &comm.val, stat.val, &0)
    stat
end

function Iprobe(src::Integer, tag::Integer, comm::Comm)
    flag = Ref{Cint}()
    stat = Status()
    ccall(MPI_IPROBE, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &src, &tag, &comm.val, flag, stat.val, &0)
    if flag[] == 0
        return false, nothing
    end
    true, stat
end

function Get_count{T}(stat::Status, ::Type{T})
    count = Ref{Cint}()
    ccall(MPI_GET_COUNT, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          stat.val, &mpitype(T), count, &0)
    Int(count[])
end

function Send{T}(buf::MPIBuffertype{T}, count::Integer,
                              dest::Integer, tag::Integer, comm::Comm)
    ccall(MPI_SEND, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}),
          buf, &count, &mpitype(T), &dest, &tag, &comm.val, &0)
end

function Send{T}(buf::Array{T}, dest::Integer, tag::Integer,
                              comm::Comm)
    Send(buf, length(buf), dest, tag, comm)
end

function Send{T}(obj::T, dest::Integer, tag::Integer, comm::Comm)
    buf = [obj]
    Send(buf, dest, tag, comm)
end

function send(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Send(buf, dest, tag, comm)
end

function Isend{T}(buf::MPIBuffertype{T}, count::Integer,
                               dest::Integer, tag::Integer, comm::Comm)
    rval = Ref{Cint}()
    ccall(MPI_ISEND, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Ptr{Cint}),
          buf, &count, &mpitype(T), &dest, &tag, &comm.val, rval, &0)
    Request(rval[], buf)
end

function Isend{T}(buf::Array{T}, dest::Integer, tag::Integer,
                               comm::Comm)
    Isend(buf, length(buf), dest, tag, comm)
end

function Isend{T}(obj::T, dest::Integer, tag::Integer, comm::Comm)
    buf = [obj]
    Isend(buf, dest, tag, comm)
end

function isend(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Isend(buf, dest, tag, comm)
end

function Recv!{T}(buf::MPIBuffertype{T}, count::Integer,
                               src::Integer, tag::Integer, comm::Comm)
    stat = Status()
    ccall(MPI_RECV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Ptr{Cint}),
          buf, &count, &mpitype(T), &src, &tag, &comm.val, stat.val, &0)
    stat
end

function Recv!{T}(buf::Array{T}, src::Integer, tag::Integer,
                               comm::Comm)
    Recv!(buf, length(buf), src, tag, comm)
end

function Recv{T}(::Type{T}, src::Integer, tag::Integer, comm::Comm)
    buf = Ref{T}()
    stat = Recv!(buf, 1, src, tag, comm)
    (buf[], stat)
end

function recv(src::Integer, tag::Integer, comm::Comm)
    stat = Probe(src, tag, comm)
    count = Get_count(stat, UInt8)
    buf = Array(UInt8, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (MPI.deserialize(buf), stat)
end

function Irecv!{T}(buf::MPIBuffertype{T}, count::Integer,
                                src::Integer, tag::Integer, comm::Comm)
    val = Ref{Cint}()
    ccall(MPI_IRECV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Ptr{Cint}),
          buf, &count, &mpitype(T), &src, &tag, &comm.val, val, &0)
    Request(val[], buf)
end

function Irecv!{T}(buf::Array{T}, src::Integer, tag::Integer,
                                comm::Comm)
    Irecv!(buf, length(buf), src, tag, comm)
end

function irecv(src::Integer, tag::Integer, comm::Comm)
    (flag, stat) = Iprobe(src, tag, comm)
    if !flag
        return (false, nothing, nothing)
    end
    count = Get_count(stat, UInt8)
    buf = Array(UInt8, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (true, MPI.deserialize(buf), stat)
end

function Wait!(req::Request)
    stat = Status()
    ccall(MPI_WAIT, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &req.val, stat.val, &0)
    req.buffer = nothing
    stat
end

function Test!(req::Request)
    flag = Ref{Cint}()
    stat = Status()
    ccall(MPI_TEST, Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &req.val, flag, stat.val, &0)
    if flag[] == 0
        return (false, nothing)
    end
    req.buffer = nothing
    (true, stat)
end

function Waitall!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    statvals = Array(Cint, MPI_STATUS_SIZE, count)
    ccall(MPI_WAITALL, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, statvals, &0)
    stats = Array(Status, count)
    for i in 1:count
        reqs[i].val = reqvals[i]
        reqs[i].buffer = nothing
        stats[i] = Status()
        for v in 1:MPI_STATUS_SIZE
            stats[i].val[v] = statvals[v,i]
        end
    end
    stats
end

function Testall!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    flag = Ref{Cint}()
    statvals = Array(Cint, MPI_STATUS_SIZE, count)
    ccall(MPI_TESTALL, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, flag, statvals, &0)
    if flag[] == 0
        return (false, nothing)
    end
    stats = Array(Status, count)
    for i in 1:count
        reqs[i].val = reqvals[i]
        reqs[i].buffer = nothing
        stats[i] = Status()
        for v in 1:MPI_STATUS_SIZE
            stats[i].val[v] = statvals[v,i]
        end
    end
    (true, stats)
end

function Waitany!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    ind = Ref{Cint}()
    stat = Status()
    ccall(MPI_WAITANY, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, ind, stat.val, &0)
    index = Int(ind[])
    reqs[index].val = reqvals[index]
    reqs[index].buffer = nothing
    (index, stat)
end

function Testany!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    ind = Ref{Cint}()
    flag = Ref{Cint}()
    stat = Status()
    ccall(MPI_TESTANY, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, ind, flag, stat.val, &0)
    index = Int(ind[])
    if flag[] == 0 || index == MPI_UNDEFINED
        return (false, 0, nothing)
    end
    reqs[index].val = reqvals[index]
    reqs[index].buffer = nothing
    (true, index, stat)
end

function Waitsome!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    outcnt = Ref{Cint}()
    inds = Array(Cint, count)
    statvals = Array(Cint, MPI_STATUS_SIZE, count)
    ccall(MPI_WAITSOME, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, outcnt, inds, statvals, &0)
    outcount = Int(outcnt[])
    # This can happen if there were no valid requests
    if outcount == UNDEFINED
        outcount = 0
    end
    indices = Array(Int, outcount)
    stats = Array(Status, outcount)
    for i in 1:outcount
        ind = Int(inds[i])
        reqs[ind].val = reqvals[ind]
        reqs[ind].buffer = nothing
        indices[i] = inds[i]
        stats[i] = Status()
        for v in 1:MPI_STATUS_SIZE
            stats[i].val[v] = statvals[v,i]
        end
    end
    (indices, stats)
end

function Testsome!(reqs::Array{Request,1})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    outcnt = Ref{Cint}()
    inds = Array(Cint, count)
    statvals = Array(Cint, MPI_STATUS_SIZE, count)
    ccall(MPI_TESTSOME, Void,
          (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          &count, reqvals, outcnt, inds, statvals, &0)
    outcount = outcnt[]
    # This can happen if there were no valid requests
    if outcount == UNDEFINED
        outcount = 0
    end
    indices = Array(Int, outcount)
    stats = Array(Status, outcount)
    for i in 1:outcount
        ind = Int(inds[i])
        reqs[ind].val = reqvals[ind]
        reqs[ind].buffer = nothing
        indices[i] = inds[i]
        stats[i] = Status()
        for v in 1:MPI_STATUS_SIZE
            stats[i].val[v] = statvals[v,i]
        end
    end
    (indices, stats)
end

function Cancel!(req::Request)
    ccall(MPI_CANCEL, Void, (Ptr{Cint},Ptr{Cint}), &req.val, &0)
    req.buffer = nothing
    nothing
end

# Collective communication

function Barrier(comm::Comm)
    ccall(MPI_BARRIER, Void, (Ptr{Cint},Ptr{Cint}), &comm.val, &0)
end

function Bcast!{T}(buffer::MPIBuffertype{T}, count::Integer,
                   root::Integer, comm::Comm)
    ccall(MPI_BCAST, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          buffer, &count, &mpitype(T), &root, &comm.val, &0)
    buffer
end

function Bcast!{T}(buffer::Array{T}, root::Integer, comm::Comm)
    Bcast!(buffer, length(buffer), root, comm)
end

#=
function Bcast{T}(obj::T, root::Integer, comm::Comm)
    buf = [T]
    Bcast!(buf, root, comm)
    buf[1]
end
=#

function bcast(obj, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    count = Array(Cint, 1)
    if isroot
        buf = MPI.serialize(obj)
        count[1] = length(buf)
    end
    Bcast!(count, root, comm)
    if !isroot
        buf = Array(UInt8, count[1])
    end
    Bcast!(buf, root, comm)
    if !isroot
        obj = MPI.deserialize(buf)
    end
    obj
end

function Reduce{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                   op::Op, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    recvbuf = Array(T, isroot ? count : 0)
    ccall(MPI_REDUCE, Void,
          (Ptr{T}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Ptr{Cint}),
          sendbuf, recvbuf, &count, &mpitype(T), &op.val, &root, &comm.val,
          &0)
    isroot ? recvbuf : nothing
end

function Reduce{T}(sendbuf::Array{T}, op::Union{Op,Function}, root::Integer, comm::Comm)
    Reduce(sendbuf, length(sendbuf), op, root, comm)
end

function Reduce{T}(object::T, op::Union{Op,Function}, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    sendbuf = T[object]
    recvbuf = Reduce(sendbuf, op, root, comm)
    isroot ? recvbuf[1] : nothing
end

function Allreduce!{T}(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                      count::Integer, op::Op, comm::Comm)
    ccall(MPI_ALLREDUCE, Void, (Ptr{T}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
          Ptr{Cint}, Ptr{Cint}), sendbuf, recvbuf, &count, &mpitype(T),
          &op.val, &comm.val, &0)

    recvbuf
end

function Allreduce!{T}(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
                      op::Union{Op,Function}, comm::Comm)
    Allreduce!(sendbuf, recvbuf, length(recvbuf), op, comm)
end

function Allreduce{T}(obj::T, op::Union{Op,Function}, comm::Comm)
    objref = Ref(obj)
    outref = Ref{T}()
    Allreduce!(objref, outref, 1, op, comm)

    outref[]
end

# allocate receive buffer automatically
function allreduce{T}(sendbuf::MPIBuffertype{T}, op::Union{Op,Function}, comm::Comm)

  recvbuf = similar(sendbuf)
  Allreduce!(sendbuf, recvbuf, length(recvbuf), op, comm)
end

include("mpi-op.jl")

function Scatter{T}(sendbuf::MPIBuffertype{T},count::Integer, root::Integer,
                    comm::Comm)
    recvbuf = Array(T, count)
    ccall(MPI_SCATTER, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), sendbuf, &count, &mpitype(T),
           recvbuf, &count, &mpitype(T), &root, &comm.val, &0)
    recvbuf
end

function Scatterv{T}(sendbuf::MPIBuffertype{T},
                     counts::Vector{Cint}, root::Integer,
                     comm::Comm)
    recvbuf = Array(T, counts[Comm_rank(comm) + 1])
    recvcnt = counts[Comm_rank(comm) + 1]
    disps = cumsum(counts) - counts
    ccall(MPI_SCATTERV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, counts, disps, &mpitype(T), recvbuf, &recvcnt, &mpitype(T), &root, &comm.val, &0)
    recvbuf
end

function Gather{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                   root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    recvbuf = Array(T, isroot ? Comm_size(comm) * count : 0)
    ccall(MPI_GATHER, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, &count, &mpitype(T), recvbuf, &count, &mpitype(T), &root, &comm.val, &0)
    isroot ? recvbuf : nothing
end

function Gather{T}(sendbuf::Array{T}, root::Integer, comm::Comm)
    Gather(sendbuf, length(sendbuf), root, comm)
end

function Gather{T}(object::T, root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    sendbuf = T[object]
    recvbuf = Gather(sendbuf, root, comm)
    isroot ? recvbuf : nothing
end

function Allgather{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                      comm::Comm)
    recvbuf = Array(T, Comm_size(comm) * count)
    ccall(MPI_ALLGATHER, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, &count, &mpitype(T), recvbuf, &count, &mpitype(T), &comm.val, &0)
    recvbuf
end

function Allgather{T}(sendbuf::Array{T}, comm::Comm)
    Allgather(sendbuf, length(sendbuf), comm)
end

function Allgather{T}(object::T, comm::Comm)
    sendbuf = T[object]
    recvbuf = Allgather(sendbuf, comm)
    recvbuf
end

function Gatherv{T}(sendbuf::MPIBuffertype{T}, counts::Vector{Cint},
                    root::Integer, comm::Comm)
    isroot = Comm_rank(comm) == root
    displs = cumsum(counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    recvbuf = Array(T, isroot ? sum(counts) : 0)
    ccall(MPI_GATHERV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, &sendcnt, &mpitype(T), recvbuf, counts, displs, &mpitype(T), &root, &comm.val, &0)
    isroot ? recvbuf : nothing
end

function Allgatherv{T}(sendbuf::MPIBuffertype{T}, counts::Vector{Cint},
                       comm::Comm)
    displs = cumsum(counts) - counts
    sendcnt = counts[Comm_rank(comm) + 1]
    recvbuf = Array(T, sum(counts))
    ccall(MPI_ALLGATHERV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, &sendcnt, &mpitype(T), recvbuf, counts, displs, &mpitype(T), &comm.val, &0)
    recvbuf
end

function Alltoall{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                     comm::Comm)
    recvbuf = Array(T, Comm_size(comm)*count)
    ccall(MPI_ALLTOALL, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, &count, &mpitype(T), recvbuf, &count, &mpitype(T), &comm.val, &0)
    recvbuf
end

function Alltoallv{T}(sendbuf::MPIBuffertype{T}, scounts::Vector{Cint},
                      rcounts::Vector{Cint}, comm::Comm)
    recvbuf = Array(T, sum(rcounts))
    sdispls = cumsum(scounts) - scounts
    rdispls = cumsum(rcounts) - rcounts
    ccall(MPI_ALLTOALLV, Void,
          (Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, scounts, sdispls, &mpitype(T), recvbuf, rcounts, rdispls, &mpitype(T), &comm.val, &0)
    recvbuf
end

function Scan{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                 op::Op, comm::Comm)
    recvbuf = Array(T, count)
    ccall(MPI_SCAN, Void,
          (Ptr{T}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, recvbuf, &count, &mpitype(T), &op.val, &comm.val, &0)
    recvbuf
end

function Scan{T}(object::T, op::Op, comm::Comm)
    sendbuf = T[object]
    Scan(sendbuf,1,op,comm)
end

function Exscan{T}(sendbuf::MPIBuffertype{T}, count::Integer,
                   op::Op, comm::Comm)
    recvbuf = Array(T, count)
    ccall(MPI_EXSCAN, Void,
          (Ptr{T}, Ptr{T}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
          sendbuf, recvbuf, &count, &mpitype(T), &op.val, &comm.val, &0)
    recvbuf
end

function Exscan{T}(object::T, op::Op, comm::Comm)
    sendbuf = T[object]
    Exscan(sendbuf,1,op,comm)
end

# Conversion between C and Fortran Comm handles:
if HAVE_MPI_COMM_C2F
    # use MPI_Comm_f2c and MPI_Comm_c2f
    Base.convert(::Type{CComm}, comm::Comm) =
        ccall(:MPI_Comm_f2c, CComm, (Cint,), comm.val)
    Base.convert(::Type{Comm}, ccomm::CComm) =
        Comm(ccall(:MPI_Comm_c2f, Cint, (CComm,), ccomm))
elseif sizeof(CComm) == sizeof(Cint)
    # in MPICH, both C and Fortran use identical Cint comm handles
    # and MPI_Comm_c2f is not provided.
    Base.convert(::Type{CComm}, comm::Comm) = reinterpret(CComm, comm.val)
    Base.convert(::Type{Comm}, ccomm::CComm) = Comm(reinterpret(Cint, ccomm))
else
    warn("No MPI_Comm_c2f found - conversion to/from MPI.CComm will not work")
end


text:
using InteractiveUtils

function _strip_comments(out::IO, inp::IO)
    for ln in eachline(inp)
        if !startswith(ln, ';')
            println(out, ln)
        end
    end
end

function code_native_nocomment(
    io::IO, @nospecialize(f), @nospecialize(types=Tuple);
    syntax::Symbol = :intel)

    tmpio = IOBuffer()
    code_native(tmpio, f, types, syntax=syntax)
    seekstart(tmpio)
    _strip_comments(io, tmpio)
end

code_native_nocomment(
    @nospecialize(f), @nospecialize(types=Tuple); syntax::Symbol = :att) =
        code_native_nocomment(stdout, f, types, syntax = syntax)

code_native_nocomment(::IO, ::Any, ::Symbol) =
    error("illegal code_native call") # resolve ambiguous call

function code_llvm_nocomment(
    io::IO, @nospecialize(f), @nospecialize(types=Tuple),
    strip_ir_metadata=true, dump_module=false, optimize=true)

    tmpio = IOBuffer()
    code_llvm(tmpio, f, types, strip_ir_metadata, dump_module, optimize)
    seekstart(tmpio)
    _strip_comments(io, tmpio)
end

code_llvm_nocomment(
    @nospecialize(f), @nospecialize(types=Tuple);
    raw=false, dump_module=false, optimize=true) =
        code_llvm_nocomment(stdout, f, types, !raw, dump_module, optimize)

macro code_native_nocomment(ex0)
    InteractiveUtils.gen_call_with_extracted_types(
        __module__, :code_native_nocomment, ex0)
end

macro code_llvm_nocomment(ex0)
    InteractiveUtils.gen_call_with_extracted_types(
        __module__, :code_llvm_nocomment, ex0)
end

max_stars_repo_name: raphbacher/julia-intro

max_issues_repo_name: raphbacher/julia-intro

max_forks_repo_name: raphbacher/julia-intro


text:
using FastGaussQuadrature
import Base: size, show, minimum, maximum, eltype

struct Grid{T<:Number,U<:Integer}
    X::Matrix{T} # Quadrature roots
    W::Matrix{T} # Quadrature weights
    N::Matrix{T} # Inverted weights for matrix elements
    nel::U # Number of finite elements
    n::U # Polynomial order of basis functions
    bl::Symbol # Left boundary condition
    br::Symbol # Right boundary condition
end

"""
    Grid(r, n[, bl[, br]])

Construct FEDVR grid, where the boundaries of the finite elements are
given by `r`. `n` is the order of the basis functions of the DVR. `bl`
and `br` decide the boundary conditions for the left and right edge of
the grid, respectively (:dirichlet0 by default)."""
function Grid(r::AbstractVector, n::Integer,
              bl::Symbol=:dirichlet0,
              br::Symbol=:dirichlet0)
    xₘ,wₘ = gausslobatto(n)
    xₘ = (xₘ .+ 1)/2
    # https://github.com/JuliaLang/julia/pull/18777
    lerp(a::T,b::T,t) where T = T(fma(t, b, fma(-t, a, a)))
    nel = length(r)-1
    X = zeros(nel, n)
    W = zeros(nel, n)
    for i = 1:nel
        X[i,:] = lerp.(r[i], r[i+1], xₘ)
        W[i,:] = 0.5*(r[i+1]-r[i])*wₘ
    end
    # Precalculate inverse weights for matrix elements
    N = zeros(nel, n)
    for i = 1:nel
        N[i,1] = 1/√(W[i,1] + (i == 1 ? 0 : W[i-1,end]))
        for m = 2:n-1
            N[i,m] = 1/√(W[i,m])
        end
        N[i,end] = 1/√(W[i,end] + (i == nel ? 0 : W[i+1,1]))
    end
    Grid(X,W,N,nel,n,bl,br)
end

"""
    elcount(grid)

Return number of finite elements in grid.
"""
elcount(grid::Grid) = grid.nel

"""
    elems(grid)

Return range of finite element indices in grid.
"""
elems(grid::Grid) = 1:elcount(grid)

"""
    order(grid)

Return polynomial order of basis functions for the DVR.
"""
order(grid::Grid) = grid.n

"""
    basecount(grid)

Return total number of basis functions in grid.
"""
basecount(grid::Grid) = elcount(grid)*order(grid) - (elcount(grid)-1) - (grid.bl == :dirichlet0 ? 1 : 0) - (grid.br == :dirichlet0 ? 1 : 0)

"""
    bases(grid, i)

Return range of basis function indices in element `i`. For an interior
element, this is always `1:n`, where `n` is the order of the `i`th
element. For the first element, in the case of a `:dirichlet0`
boundary condition, the result is `2:n`, otherwise
`1:n`. Equivalently, for the last element, in case of `:dirichlet0`
boundary conditions, the result is `1:n-1`."""
function bases(grid::Grid, i::Integer)
    n = order(grid)
    if i == 1 && grid.bl == :dirichlet0
        2:n
    elseif i == grid.nel && grid.br == :dirichlet0
        1:n-1
    else
        1:n
    end
end

minimum(grid::Grid) = min(grid.X[1],grid.X[end])
maximum(grid::Grid) = max(grid.X[1],grid.X[end])

eltype(grid::Grid) = eltype(grid.X)

function boundary_sel(grid::Grid, v::AbstractVector)
    if grid.bl == :dirichlet0 && grid.br == :dirichlet0
        v[2:end-1]
    elseif grid.bl == :dirichlet0
        v[2:end]
    elseif grid.br == :dirichlet0
        v[1:end-1]
    else
        v
    end
end

"""
    locs(grid)

Return locations of Gauss–Lobatto quadrature points.
"""
locs(grid::Grid) = boundary_sel(grid, [grid.X[:,1:end-1]'[:]..., grid.X[end]]::Vector{eltype(grid.X)})

"""
    weights(grid)

Return weights of Gauss–Lobatto quadrature points.
"""
function weights(grid::Grid)
    bridges = grid.W[1:end-1,end] + grid.W[2:end,1]
    wl = grid.bl == :dirichlet0 ? zero(grid.W[1]) : grid.W[1]
    wr = grid.br == :dirichlet0 ? zero(grid.W[1]) : grid.W[end]
    boundary_sel(grid, [wl, [grid.W[:,2:end-1] [bridges;wr]]'[:]...])
end

"""
    find_interval(X, x, i, sel)

Find the interval of points in `x` covering the `i`:th row in
`X`. The search is started from the first element in `sel`."""
function find_interval(X, x, i, sel)
    lt(v) = x -> x < v
    a = something(findlast(lt(X[i,1]), x[first(sel):end]), 0) + first(sel)
    b = max(a,last(sel))
    b = something(findlast(lt(X[i,end]), x[b:end]), 0) + b - 1
    if i == size(X,1) && b < length(x)
        # If last finite element and *not* :dirichlet0 boundary
        # conditions, include finite-element end-point in interval.
        b += 1
    end
    a:b
end

show(io::IO, grid::Grid) =
    write(io, "FEDVR Grid{$(eltype(grid))} with $(elcount(grid)) elements of order $(order(grid)) [$(grid.bl),$(grid.br)]")

@recipe function plot(grid::Grid)
    markershape --> :circle
    label --> "FEDVR grid locations"
    xticks --> 1:order(grid)-1:basecount(grid)
    yticks --> [grid.X[:,1]..., grid.X[end,end]]
    xlabel --> "#"
    ylabel --> "x"
    locs(grid)
end

export elcount, elems, order, basecount, locs, weights


text:
"""
    ADVI(samples_per_step = 10, max_iters = 5000)

Automatic Differentiation Variational Inference (ADVI) for a given model.
"""
struct ADVI{AD} <: VariationalInference{AD}
    samples_per_step # number of samples used to estimate the ELBO in each optimization step
    max_iters        # maximum number of gradient steps used in optimization
end

ADVI(args...) = ADVI{ADBackend()}(args...)
ADVI() = ADVI(10, 5000)

alg_str(::ADVI) = "ADVI"

function vi(model::Model, alg::ADVI; optimizer = TruncatedADAGrad())
    # setup
    varinfo = Turing.VarInfo(model)
    num_params = sum([size(varinfo.metadata[sym].vals, 1) for sym ∈ keys(varinfo.metadata)])

    dists = vcat([varinfo.metadata[sym].dists for sym ∈ keys(varinfo.metadata)]...)

    num_ranges = sum([length(varinfo.metadata[sym].ranges) for sym ∈ keys(varinfo.metadata)])
    ranges = Vector{UnitRange{Int}}(undef, num_ranges)
    idx = 0
    range_idx = 1
    for sym ∈ keys(varinfo.metadata)
        for r ∈ varinfo.metadata[sym].ranges
            ranges[range_idx] = idx .+ r
            range_idx += 1
        end

        # append!(ranges, [idx .+ r for r ∈ varinfo.metadata[sym].ranges])
        idx += varinfo.metadata[sym].ranges[end][end]
    end

    q = Variational.MeanField(zeros(num_params), zeros(num_params), dists, ranges)

    # construct objective
    elbo = ELBO()

    Turing.DEBUG && @debug "Optimizing ADVI..."
    θ = optimize(elbo, alg, q, model; optimizer = optimizer)
    μ, ω = θ[1:length(q)], θ[length(q) + 1:end]

    return MeanField(μ, ω, dists, ranges)
end

# TODO: implement optimize like this?
# (advi::ADVI)(elbo::EBLO, q::MeanField, model::Model) = begin
# end

function optimize(elbo::ELBO, alg::ADVI, q::MeanField, model::Model; optimizer = TruncatedADAGrad())
    θ = randn(2 * length(q))
    optimize!(elbo, alg, q, model, θ; optimizer = optimizer)

    return θ
end

function (elbo::ELBO)(
    alg::ADVI,
    q::MeanField,
    model::Model,
    θ::AbstractVector{T},
    num_samples
) where T <: Real
    # setup
    varinfo = Turing.VarInfo(model)

    num_params = length(q)
    μ, ω = θ[1:num_params], θ[num_params + 1: end]

    elbo_acc = 0.0

    # TODO: instead use `rand(q, num_samples)` and iterate through?
    # Requires new interface for Bijectors.jl

    for i = 1:num_samples
        # iterate through priors, sample and update
        idx = 0
        z = zeros(T, num_params)

        for sym ∈ keys(varinfo.metadata)
            md = varinfo.metadata[sym]

            for i = 1:size(md.dists, 1)
                prior = md.dists[i]
                r = md.ranges[i] .+ idx

                # mean-field params for this set of model params
                μ_i = μ[r]
                ω_i = ω[r]

                # obtain samples from mean-field posterior approximation
                η = randn(length(μ_i))
                ζ = center_diag_gaussian_inv(η, μ_i, exp.(ω_i))

                # inverse-transform back to domain of original priro
                z[r] .= invlink(prior, ζ)

                # update
                # @info θ
                # z[md.ranges[i]] .= θ
                # @info md.vals

                # add the log-det-jacobian of inverse transform;
                # `logabsdet` returns `(log(abs(det(M))), sign(det(M)))` so return first entry
                # add `eps` to ensure SingularException does not occurr in `logabsdet`
                elbo_acc += logabsdet(jac_inv_transform(prior, ζ) .+ eps(T))[1] / num_samples
            end

            idx += md.ranges[end][end]
        end

        # compute log density
        varinfo = VarInfo(varinfo, SampleFromUniform(), z)
        model(varinfo)
        elbo_acc += varinfo.logp / num_samples
    end

    # add the term for the entropy of the variational posterior
    variational_posterior_entropy = sum(ω)
    elbo_acc += variational_posterior_entropy

    return elbo_acc
end

function (elbo::ELBO)(alg::ADVI, q::MeanField, model::Model, num_samples)
    # extract the mean-field Gaussian params
    θ = vcat(q.μ, q.ω)

    return elbo(alg, q, model, θ, num_samples)
end


text:
@testset "size(A, i)" begin
    tV = SimpleRDArray(V)
    @test_throws BoundsError size(tV, 0)
    @test size(tV, 0x1) == 3
    @test size(tV, 0x2) == 1
end

@testset "ArrayInterface.strides" begin
    tV = SimpleRDArray(V)
    @test ArrayInterface.strides(tV) == (StaticInt(1),)
    @test ArrayInterface.strides(tV, 1) == 1
    @test ArrayInterface.stride(tV, StaticInt(1)) isa StaticInt
    @test ArrayInterface.stride(tV, 1) == 1
    tM = SimpleRDArray(M)
    @test strides(tM) == (StaticInt(1), 3)
    @test ArrayInterface.strides(tM, 1) == 1
    @test ArrayInterface.strides(tM, 2) == 3
    @test ArrayInterface.stride(tM, StaticInt(1)) isa StaticInt
    @test ArrayInterface.stride(tM, 1) == 1
    @test ArrayInterface.stride(tM, 2) == 3
    tT = SimpleRDArray(T)
    @test strides(tT) == (StaticInt(1), 3, 9)
    @test ArrayInterface.strides(tT, 1) == 1
    @test ArrayInterface.strides(tT, 2) == 3
    @test ArrayInterface.strides(tT, 3) == 9
    @test ArrayInterface.stride(tT, StaticInt(1)) isa StaticInt
    @test ArrayInterface.stride(tT, 1) == 1
    @test ArrayInterface.stride(tT, 2) == 3
    @test ArrayInterface.stride(tT, 3) == 9
end

@testset "non-resizable arrays" begin
    @test_throws MethodUndefineError sizehint!(1:2, 3)
    @test_throws MethodUndefineError resize!(1:2, (3,))
    @test_throws MethodUndefineError resize!(1:2, 3, 1)
end

@testset "wrong index" begin
    tV = SimpleRDArray(V)
    @test_throws DimBoundsError resize!(tV, 1, Bool[0, 0, 1, 0])
    @test_throws DimBoundsError resize!(tV, 2, Bool[0, 0, 1, 0])
    @test_throws DimBoundsError resize!(tV, 1, 3:4)
    @test_throws DimBoundsError resize!(tV, 2, 1:2)
end

@testset "tailn" begin
    @test tailn(Val(1), 1, 2, 3, 4) == (4,)
    @test tailn(Val(2), 1, 2, 3, 4) == (3, 4)
    @test tailn(Val(3), 1, 2, 3, 4) == (2, 3, 4)
    @test tailn(Val(4), 1, 2, 3, 4) == (1, 2, 3, 4)
    @test tailn(Val(5), 1, 2, 3, 4) == (1, 2, 3, 4)
end

@testset "Size" begin
    tp = ntuple(UInt, Val(4))
    sz = Size(4, 4, 4, 4)
    @test length(sz) == 4
    set!(sz, tp)
    @test sz == tp
    @test setindex!(sz, 0, 1) == (0, 2, 3, 4)
end

# test for LoopVectorization
@testset "@turbo" begin
    tV = SimpleRDArray(V)
    tM = SimpleRDArray(M)
    @test @avx(tV .* tV) ≈ (tV .* tV)
    @test @avx(tV .* tM) ≈ (tV .* tM)
    @test @avx(tM .* tM) ≈ (tM .* tM)
    @test @avx(V .* tV) ≈ (V .* tV)
    @test @avx(tV .* V) ≈ (tV .* V)
    @test @avx(V .* tM) ≈ (V .* tM)
    @test @avx(tV .* M) ≈ (tV .* M)
    @test @avx(M .* tM) ≈ (M .* tM)
    @test @avx(tM .* M) ≈ (tM .* M)
    function A_mul_B(A, B) # this function is just for test
        C = Array{Base.promote_op(*, eltype(A), eltype(B))}(undef, size(A, 1), size(B, 2))
        @avx for n ∈ axes(B, 2), m ∈ axes(A, 1)
            Cmn = zero(eltype(C))
            for k ∈ axes(B, 1)
                Cmn += A[m,k] * B[k,n]
            end
            C[m,n] = Cmn
        end
        return C
    end
    @test A_mul_B(M, M) ≈ M * M
    @test A_mul_B(M, tM) ≈ M * tM
    @test A_mul_B(tM, M) ≈ tM * M
    @test A_mul_B(tM, tM) ≈ tM * tM
end


text:
using FindPeaks1D

using Test

@testset "FindPeaks1D" begin
    x = [13, 12, 14, 18, 19, 19, 19, 15, 11, 6, 4, 10, 8, 13, 8, 13, 3, 18, 7, 4]

    pkindices, _ = FindPeaks1D.localmaxima1d(x)
    @test pkindices == [6, 12, 14, 16, 18]

    pkindicesh1, _ = findpeaks1d(x; height=11)
    pkindicesh2, _ = findpeaks1d(x; height=(11, 16))
    @test pkindicesh1 == [6, 14, 16, 18]
    @test pkindicesh2 == [14, 16]

    pkindicesd1, _ = findpeaks1d(x, distance=1)
    pkindicesd2, _ = findpeaks1d(x, distance=3)
    @test pkindicesd1 == [6, 12, 14, 16, 18]
    @test pkindicesd2 == [6, 14, 18]

    prominences1, leftbases1, rightbases1 = peakprominences1d(x, pkindices)
    prominences2, leftbases2, rightbases2 = peakprominences1d(x, pkindices, 2)
    @test prominences1 == [7, 2, 9, 9, 14]
    @test leftbases1 == [2, 11, 11, 11, 17]
    @test rightbases1 == [17, 13, 17, 17, 20]
    @test prominences2 == [0, 2, 5, 5, 11]
    @test leftbases2 == [6, 11, 13, 15, 17]
    @test rightbases2 == [6, 13, 15, 17, 19]

    widths1, widthheight1, leftips1, rightips1 = peakwidths1d(x, pkindices, 1.0)
    @test all(isapprox.(widths1 ,[6.75, 1.3333, 5.900, 5.900, 2.9333], atol=0.001))
    @test widthheight1 == [12, 8, 4, 4, 4]
    @test all(isapprox.(leftips1 ,[2.0, 11.6667, 11.0, 11.0, 17.0667], atol=0.001))
    @test all(isapprox.(rightips1 ,[8.75, 13.0, 16.9, 16.9, 20.0], atol=0.001))

    pkindices1, _ = findpeaks1d(x, width=3.0, relheight=1.0)
    @test pkindices1 == [6,14,16]
    pkindices2, _ = findpeaks1d(x, width=2.0, relheight=0.5)
    @test pkindices2 == [6]
    pkindices3, _ = findpeaks1d(x, prominence=3)
    @test pkindices3 == [6,14,16,18]
    pkindices4, _ = findpeaks1d(x, prominence=100)
    @test pkindices4 == Int64[]
    # test for less strict type signature
    pkindices1, _ = findpeaks1d(x, width=3//1, relheight=1.0)
    @test pkindices1 == [6,14,16]
    pkindices2, _ = findpeaks1d(x, width=2, relheight=1//2)
    @test pkindices2 == [6]
    pkindices3, _ = findpeaks1d(x, prominence=3//1)
    @test pkindices3 == [6,14,16,18]
    pkindices4, _ = findpeaks1d(x, prominence=100//1)
    @test pkindices4 == Int64[]

    @test_throws ArgumentError pkindiceswlenerror, _ = findpeaks1d(x, prominence=0.1, wlen=0)

    x = [0, 1, 1, 1, 1, 1, 0]
    pkindices, _ = FindPeaks1D.localmaxima1d(x)
    @test length(pkindices) == 1
    @test pkindices[1] == 4

end


text:
using NNlib: gather, gather!

@testset "gather scalar index" begin
    T = Float32

    ## 1d src, 2d index of ints -> 2d output
    src = T[3, 4, 5, 6, 7]
    index = [1 2 3 4;
            4 2 1 3;
            3 5 5 3]
    output = T[3 4 5 6;
              6 4 3 5;
              5 7 7 5]

    y = gather(src, index)
    @test y isa Array{T,2}
    @test size(y) == size(index)
    @test y == output
    @test gather!(T.(zero(index)), src, index) == output
    @test_throws ArgumentError gather!(zeros(T, 3, 5), src, index)

    index2 = [1 2 3 4;
              4 2 1 3;
              3 6 5 3]
    @test_throws BoundsError gather!(T.(zero(index)), src, index2)

    ## 1d src, 3d index of ints -> 3d output
    src = T[3, 4, 5, 6, 7]
    index = [1 2 3 4;
            4 2 1 3;
            3 5 5 3][:,:,1:1]
    output = T[3 4 5 6;
              6 4 3 5;
              5 7 7 5][:,:,1:1]

    y = gather(src, index)
    @test y isa Array{T,3}
    @test size(y) == size(index)
    @test y == output

    ## 2d src, 2d index of ints -> 3d output
    src = T[3 5 7
            4 6 8]
    index = [1 2 3;
            2 2 1;
            3 1 3]

    output = zeros(T, 2, 3, 3)

    output[:,:,1] = [3 5 7
                    4 6 8]

    output[:,:,2] = [5 5 3
                    6 6 4]

    output[:,:,3] = [7 3 7
                    8 4 8]

    y = gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa Array{T,3}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
    @test y == output
end

@testset "gather tuple index" begin
    T = Float32

    ## 2d src, 1d index of 2-tuples -> 1d output
    src = T[3 5 7
            4 6 8]

    index = [(1,1), (1,2), (1,3), (2,1), (2,2), (2,3)]

    output = T[3, 5, 7, 4, 6, 8]

    y = gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa Array{T,1}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
    @test y == output

    ## 3d src, 2d index of 2-tuples -> 3d output
    n1, nsrc, nidx = 2, 3, 6
    src = rand(Float32, n1, nsrc, nsrc)
    index = [(rand(1:nsrc), rand(1:nsrc)) for i=1:nidx, j=1:nidx]

    y = gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa Array{T,3}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
end

@testset "gather cartesian index" begin
    T = Float32

    ## 2d src, 1d index of 2-tuples -> 1d output
    src = T[3 5 7
            4 6 8]

    index = CartesianIndex.([(1,1), (1,2), (1,3), (2,1), (2,2), (2,3)])

    output = T[3, 5, 7, 4, 6, 8]

    y = gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa Array{T,1}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
    @test y == output

    ## 3d src, 2d index of 2-tuples -> 3d output
    n1, nsrc, nidx = 2, 3, 6
    src = rand(Float32, n1, nsrc, nsrc)
    index = [CartesianIndex((rand(1:nsrc), rand(1:nsrc))) for i=1:nidx, j=1:nidx]

    y = gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa Array{T,3}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
end

@testset "gather gradient for scalar index" begin
    T = Float64
    src = T[3, 4, 5, 6, 7]
    index = [1 2 3 4;
            4 2 1 3;
            3 5 5 3]
    dst = T[3 4 5 6;
            6 4 3 5;
            5 7 7 5]

    gradtest(xs -> gather!(dst, xs, index), src)
    gradtest(xs -> gather(xs, index), src)
end

@testset "gather gradient for tuple index" begin
    T = Float64
    src = T[3 5 7
            4 6 8]
    index = [(1,1), (1,2), (1,3), (2,1), (2,2), (2,3)]
    dst = T[3, 5, 7, 4, 6, 8]

    gradtest(xs -> gather!(dst, xs, index), src)
    gradtest(xs -> gather(xs, index), src)
end


text:
function main()

    configData = getConfig()
    wallpts, Interiorpts, outerpts, shapepts = 0,0,0,0

    # file_name = string(ARGS[2])
    folder_name = string(ARGS[2])

    numPoints = findTotalPoints(folder_name)
    println(numPoints)
    # globaldata = Array{Point,1}(undef, numPoints)
    # table = Array{Int32,1}(undef, numPoints)
    # wallptsidx = Array{Int32,1}(undef, 0)
    # Interiorptsidx = Array{Int32,1}(undef, 0)
    # outerptsidx = Array{Int32,1}(undef, 0)
    # shapeptsidx = Array{Int32,1}(undef, 0)
    # print(splitdata[1:3])
    defprimal = getInitialPrimitive(configData)

    # local_points_holder = Array{Array{Point,1},1}(undef, nworkers())
    # res_old_holder = Array{Float64,1}(undef, nworkers())
    ghost_holder = Array{Dict{Int32,Point},1}(undef, nworkers())
    ghost_holder_mutable = Array{Dict{Int32, Array{Float64,1}}, 1}(undef, nworkers())
    # q_ghost_holder = Array{Dict{Int64,TempQ},1}(undef, nworkers())
    # dq_ghost_holder = Array{Dict{Int64,TempDQ},1}(undef, nworkers())
    # prim_ghost_holder = Array{Dict{Int64,TempPrim},1}(undef, nworkers())

    global_local_map_index = Dict{Tuple{Float64, Float64},Int64}()
    global_local_direct_index = Array{Int32, 1}(undef, numPoints)
    files_length = Array{Int64, 1}(undef, length(workers()))
    actual_files_length = Array{Int64, 1}(undef, length(workers()))
    println("Indexing")

    createGlobalLocalMapIndex(global_local_map_index, global_local_direct_index, folder_name::String, files_length, actual_files_length)
    # println(global_local_map_index)

    println("Start Read")
    # count = 0
    # readFile(file_name::String, globaldata, table, defprimal, wallptsidx, outerptsidx, Interiorptsidx, shapeptsidx,
        # wallpts, Interiorpts, outerpts, shapepts, numPoints)

    println("Read Local Files")
    # readDistribuedFile(folder_name::String, defprimal, 12, global_local_map_index)
    globaldata_parts = [@spawnat p readDistribuedFile(folder_name::String, defprimal, p, global_local_map_index) for p in workers()]
    globaldata_parts = reshape(globaldata_parts, (nworkers()))
    dist_globaldata = DArray(globaldata_parts)
    # globaldata_parts = nothing

    q_parts = [@spawnat p readDistribuedFileQ(folder_name::String, defprimal, p, global_local_map_index, files_length) for p in workers()]
    q_parts = reshape(q_parts, (nworkers()))
    dist_q = DArray(q_parts)
    # q_parts = nothing

    dq_parts = [@spawnat p readDistribuedFileDQ(folder_name::String, defprimal, p, global_local_map_index, files_length) for p in workers()]
    dq_parts = reshape(dq_parts, (nworkers()))
    dist_dq = DArray(dq_parts)
    # dq_parts = nothing

    globaldata_parts_mutable = [@spawnat p readDistribuedFileMutables(folder_name::String, p, actual_files_length) for p in workers()]
    globaldata_parts_mutable = reshape(globaldata_parts_mutable, (1,(nworkers())))
    dist_globaldata_mutable = DArray(globaldata_parts_mutable)
    # globaldata_parts_mutable = nothing

    println("Reading Ghost")
    readGhostFile(folder_name, ghost_holder, global_local_map_index, dist_globaldata)
    ghost_holder = distribute(ghost_holder, procs=workers(), dist=(length(workers()),))

    # readGhostFileMutables(folder_name, ghost_holder_mutable)
    # ghost_holder_mutable = distribute(ghost_holder_mutable, procs=workers(), dist=(length(workers()),))

    format = configData["format"]["type"]
    if format == 1
        interior = configData["point"]["interior"]
        wall = configData["point"]["wall"]
        outer = configData["point"]["outer"]
        @sync for pid in procs(dist_globaldata)
            @spawnat pid begin
                placeNormals(dist_globaldata[:L], dist_globaldata, ghost_holder[:L], interior, wall, outer)
            end
        end
    end

    println("Start Table Sorting")
    @sync for pid in procs(dist_globaldata)
        @spawnat pid begin
            calculateConnectivity(dist_globaldata[:L], dist_globaldata, ghost_holder[:L])
            # setConnectivity(dist_globaldata[idx], connectivity)
        end
    end
    # @showprogress 3 "Computing Table" for idx in table
    #     connectivity = calculateConnectivity(globaldata, idx)
    #     setConnectivity(globaldata[idx], connectivity)
    #     # smallest_dist(globaldata, idx)
    #     # if idx % (length(table) * 0.25) == 0
    #     #     println("Bump In Table")
    #     # end
    # end

    # points_holder = DArray(reshape([lph[1]..., lph[2]...,lph[3]...,lph[4]...], :))
    @sync for pid in procs(dist_globaldata)
        @spawnat pid begin
            println(localindices(dist_globaldata))
        end
    end

    println("! Transferring FixedData to GPU")
    @sync for pid in workers()
        @spawnat pid begin
            numPoints = files_length[pid-1]
            totalPoints = actual_files_length[pid-1]
            global gpuLocNumPoints = numPoints
            locDataFixedPoint = Array{FixedPoint,1}(undef, totalPoints)
            locDataConn = zeros(Int32, 55, numPoints)
            locGlobalData = dist_globaldata[:L]

            for idx in 1: numPoints
                convertToFixedArray(locDataFixedPoint, locGlobalData[idx], idx)
                convertToNeighbourArray(locDataConn, locGlobalData[idx], idx)
            end

            locGhostGlobalData = ghost_holder[:L][1]
            numPoints += 1
            for idx in numPoints:totalPoints
                convertToFixedArray(locDataFixedPoint, locGhostGlobalData[idx], idx)
            end

            # locGhostGlobalData = ghost_holder[:L][1]
            # localkeys = keys(locGhostGlobalData)
            # locGhostDataFixedPoint = Array{Dict{Int32,FixedPoint},1}(undef, 1)
            # locGhostDataFixedPoint[1] = Dict{Int32,FixedPoint}()
            # for iter in localkeys
            #     convertToFixedArray(locGhostDataFixedPoint[1], locGhostGlobalData[iter], iter)
            # end

            global gpuLocDataFixedPoint = CuArray(locDataFixedPoint)
            global gpuLocDataConn = CuArray(locDataConn)
            # global gpuLocGhostDataFixedPoint = CuArray(locGhostDataFixedPoint)
            # global gpuLocGhostGlobalDataMutable = CuArray(locGhostGlobalDataMutable)
            # @cuda changeToOne(cutest)
            # part_test[:L] = Array(cutest)

            global gpuConfigData = CuArray([
                            getConfig()["core"]["points"],#1
                            getConfig()["core"]["cfl"],
                            getConfig()["core"]["max_iters"],
                            getConfig()["core"]["mach"],
                            getConfig()["core"]["aoa"],#5
                            getConfig()["core"]["power"],
                            getConfig()["core"]["limiter_flag"],
                            getConfig()["core"]["vl_const"],
                            getConfig()["core"]["initial_conditions_flag"],
                            getConfig()["core"]["interior_points_normal_flag"],#10
                            getConfig()["core"]["shapes"],
                            getConfig()["core"]["rho_inf"],
                            getConfig()["core"]["pr_inf"],
                            getConfig()["core"]["threadsperblock"],
                            getConfig()["core"]["gamma"],#15
                            getConfig()["core"]["clcd_flag"],
                            getConfig()["point"]["wall"],
                            getConfig()["point"]["interior"],
                            getConfig()["point"]["outer"]
                        ])
        end
    end

    # println(dist_globaldata[3])

    # println(wallptsidx)

    # return

    # configData = distribute(configData)
    # println(globaldata[1])
    # println(globaldata[end])
    # println("Size is before: ",length(globaldata))
    # dist_globaldata = distribute(globaldata, procs=workers(), dist=(length(workers()),))
    # println("Size is: ", length(globaldata))

    println(Int(getConfig()["core"]["max_iters"]) + 1)
    function run_code(ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
        fpi_solver((Int(getConfig()["core"]["max_iters"])), ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
    end

    res_old = dzeros(nworkers())
    res_new = dzeros(nworkers())
    function test_code(ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
        println("! Starting warmup function")
        fpi_solver(1, ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
        # res_old[1] = 0.0
        # Profile.clear_malloc_data()
        # @trace(fpi_solver(1, globaldata, configData, wallptsidx, outerptsidx, Interiorptsidx, res_old), maxdepth = 3)
        # res_old[1] = 0.0
        # fpi_solver(1, globaldata, configData, wallptsidx, outerptsidx, Interiorptsidx, res_old)
        # @profile fpi_solver(1, globaldata, configData, wallptsidx, outerptsidx, Interiorptsidx, res_old)
        # Profile.print()
        # res_old[1] = 0.0
        println("! Starting Main Function")
        @timeit to "nest 1" begin
            run_code(ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
        end
    end

    test_code(ghost_holder, dist_globaldata, dist_q, dist_dq, dist_globaldata_mutable, configData, res_old, res_new, numPoints)
    println("! Work Completed")
    # # println(to)
    open("temp/timercuda" * string(numPoints) * "_" * string(getConfig()["core"]["max_iters"]) *
        "_" * string(length(workers())) *".txt", "w") do io
        print_timer(io, to)
    end

    println("! Free GPU")
    @sync for pid in workers()
        @spawnat pid begin
            Array(gpuLocDataFixedPoint)
            Array(gpuLocDataConn)
            # Array(gpuLocGhostDataFixedPoint)
            # Array(gpuLocGhostGlobalDataMutable)
            # @cuda changeToOne(cutest)
            # part_test[:L] = Array(cutest)
        end
    end

    # t = rmprocs(2, 3, waitfor=0)
    # wait(t)
    # rmprocs(length(workers()))
    # @spawnat 2 begin
    # println(dist_globaldata[1])
    # println(dist_globaldata[end])
    # globaldata = makelocal(dist_globaldata)
    # close(dist_globaldata)
    # end
    # compute_cl_cd_cm(globaldata, configData, shapeptsidx)

    # println(IOContext(stdout, :compact => false), globaldata[1].q)
    # println(IOContext(stdout, :compact => false), globaldata[1].dq)
    # println(IOContext(stdout, :compact => false), globaldata[100].q)
    # println(IOContext(stdout, :compact => false), globaldata[100].dq)
    # println(IOContext(stdout, :compact => false), globaldata[1000].q)
    # println(IOContext(stdout, :compact => false), globaldata[1000].dq)
    # println()
    # println(IOContext(stdout, :compact => false), globaldata[1].flux_res)
    # println(IOContext(stdout, :compact => false), globaldata[100].flux_res)
    # println(IOContext(stdout, :compact => false), globaldata[1000].flux_res)
    # println()
    # println(IOContext(stdout, :compact => false), globaldata[1].delta)
    # println(IOContext(stdout, :compact => false), globaldata[100].delta)
    # println(IOContext(stdout, :compact => false), globaldata[1000].delta)
    # println()
    # println(IOContext(stdout, :compact => false), globaldata[1].prim)
    # println(IOContext(stdout, :compact => false), globaldata[100].prim)
    # println(IOContext(stdout, :compact => false), globaldata[1000].prim)
    # println(IOContext(stdout, :compact => false), globaldata[100].ypos_conn)
    # println(IOContext(stdout, :compact => false), globaldata[100].yneg_conn)
    # println(globaldata[1])

    # file = open("results/primvals"* "_" * string(getConfig()["core"]["max_iters"]) *
    #  "_" * string(numPoints) * ".txt", "w")
    # @showprogress 1 "This takes time" for (idx, _) in enumerate(dist_globaldata)
    #     primtowrite = dist_globaldata[global_local_direct_index[idx]].prim
    #     for element in primtowrite
    #         @printf(file,"%0.17f", element)
    #         @printf(file, " ")
    #     end
    #     print(file, "\n")
    # end
    # close(file)

    println("! Write To Files")
    # @sync for ip in procs(dist_globaldata)
    #     @spawnat ip begin
    #         writeToFile(dist_globaldata[:L], numPoints)
    #     end
    # end

    println("! Close Distributed Arrays")
    # d_closeall()
    exit()
    close(ghost_holder)
    close(ghost_holder_mutable)
    close(dist_dq)
    close(dist_q)
    close(dist_globaldata)
    close(dist_globaldata_mutable)
end

function testraid()
    test = [@spawnat p returnArrayOne() for p in workers()]
    test_parts = reshape(test, (nworkers()))
    part_test = DArray(test_parts)
    println(part_test)
    @timeit to "nest 4" begin
        @sync for pid in workers()
            @spawnat pid begin
                global cutest = CuArray(part_test[:L])
                # @cuda changeToOne(cutest)
                # part_test[:L] = Array(cutest)
            end
        end

        @sync for pid in workers()
            @spawnat pid begin
                # println(part_test[:L])
                @cuda changeToOne(cutest)
                # part_test[:L] = Array(cutest)
            end
        end

        @sync for pid in workers()
            @spawnat pid begin
                # println(part_test[:L])
                # @cuda changeToOne(cutest)
                part_test[:L] = Array(cutest)
            end
        end
    end

    println(part_test)
    open("temp/timer_stop.txt", "w") do io
        print_timer(io, to)
    end
    close(part_test)
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
# Driver script for all unit tests

using Jive

runtests(@__DIR__, skip=["runtests.jl"])


text:
type CVTimer
    tfd::Float64
    tc::Float64
    ts::Float64
    tr::Float64
    tt::Float64
    function CVTimer()
        this = new();
        tfd = 0;
        tc = 0;
        ts = 0;
        tr = 0;
        tt = 0;
        return this
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
@testset "cutensor macro for elementary operations" begin
    t0 = time()
    A = randn(Float64, (3,5,4,6))
    p = (4,1,3,2)
    @cutensor C1[4,1,3,2] := A[1,2,3,4]
    C2 = permutedims(A, p)
    @test collect(C1) ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @cutensor C[1,2,3,4] := A[1,2,3]
    end
    @test_throws TensorContractionsXD.IndexError begin
        @cutensor C[1,2,3,4] := A[1,2,2,4]
    end
    println("tensorcopy: $(time()-t0) seconds")
    t0 = time()

    B = randn(Float64, (5,6,3,4))
    p = [3,1,4,2]
    α = randn(Float64)
    @cutensor C1[3,1,4,2] := A[3,1,4,2] + α * B[1,2,3,4] # mixed eltype add/permutation
    C2 = A + α * permutedims(B, p)
    @test collect(C1) ≈ C2
    @test_throws CUDA.CUTENSOR.CUTENSORError begin
        @cutensor C[1,2,3,4] := A[1,2,3,4] + B[1,2,3,4]
    end
    println("tensoradd: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (50,100,100))
    @cutensor C1[a] := A[a, b', b']
    C2=zeros(50)
    for i=1:50
        for j=1:100
            C2[i]+=A[i, j, j]
        end
    end
    @test collect(C1) ≈ C2
    A=randn(Float64, (3,20,5,3,20,4,5))
    @cutensor C1[e, a, d] := A[a, b, c, d, b, e, c]
    C2=zeros(4,3,3)
    for i1=1:4, i2=1:3, i3=1:3
        for j1=1:20, j2=1:5
            C2[i1, i2, i3]+=A[i2, j1, j2, i3, j1, i1, j2]
        end
    end
    @test collect(C1) ≈ C2
    println("tensortrace: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (3,20,5,3,4))
    B=randn(Float64, (5,6,20,3))
    @cutensor C1[a, g, e, d, f] := A[a, b, c, d, e] * B[c, f, b, g]
    C2=zeros(3,3,4,3,6)
    for a=1:3, b=1:20, c=1:5, d=1:3, e=1:4, f=1:6, g=1:3
        C2[a, g, e, d, f] += A[a, b, c, d, e] * B[c, f, b, g]
    end
    @test collect(C1) ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @cutensor A[a, b, c, d] * B[c, f, b, g]
    end
    println("tensorcontract 1: $(time()-t0) seconds")
    t0 = time()

    A=randn(Float64, (5,5,5,5))
    B=rand(ComplexF64, (5,5,5,5))
    # @cutensor C1[1,2,5,6,3,4,7,8] := A[1,2,3,4] * B[5,6,7,8]
    # C2=reshape(kron(reshape(B, (25,25)), reshape(A, (25,25))), (5,5,5,5,5,5,5,5))
    # @test C1 ≈ C2
    @test_throws TensorContractionsXD.IndexError begin
        @cutensor C[a, b, c, d, e, f, g, i] := A[a, b, c, d] * B[e, f, g, h]
    end
    println("tensorcontract 2: $(time()-t0) seconds")
    t0 = time()
end

@testset "cutensor macro for more complicated expressions" begin
    t0 = time()
    Da=10
    Db=15
    Dc=4
    Dd=8
    De=6
    Df=7
    Dg=3
    Dh=2
    A=rand(ComplexF64, (Dc, Da, Df, Da, De, Db, Db, Dg))
    B=rand(ComplexF64, (Dc, Dh, Dg, De, Dd))
    C=rand(ComplexF64, (Dd, Dh, Df))
    @cutensor D1[d, f, h] := A[c, a, f, a, e, b, b, g] * B[c, h, g, e, d]+0.5 * C[d, h, f]
    D2=zeros(ComplexF64, (Dd, Df, Dh))
    for d=1:Dd, f=1:Df, h=1:Dh
        D2[d, f, h] += 0.5 * C[d, h, f]
        for a=1:Da, b=1:Db, c=1:Dc, e=1:De, g=1:Dg
            D2[d, f, h] += A[c, a, f, a, e, b, b, g] * B[c, h, g, e, d]
        end
    end
    @test collect(D1) ≈ D2
    @test norm(vec(D1)) ≈ sqrt(abs((@cutensor scalar(D1[d, f, h] * conj(D1[d, f, h])))))
    println("tensorcontract 3: $(time()-t0) seconds")
    t0 = time()

    # Example from README.md
    α=randn()
    A=randn(5,5,5,5,5,5)
    B=randn(5,5,5)
    C=randn(5,5,5)
    D=zeros(5,5,5)
    @cutensor begin
        D[a, b, c] = A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
        E[a, b, c] := A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
    end
    @test D == collect(E)
    println("readme example: $(time()-t0) seconds")
    t0 = time()

    # Some tensor network examples
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        D1, D2, D3 = 30, 40, 20
        d1, d2 = 2, 3
        A1 = randn(T, D1, d1, D2)
        A2 = randn(T, D2, d2, D3)
        rhoL = randn(T, D1, D1)
        rhoR = randn(T, D3, D3)
        H = randn(T, d1, d2, d1, d2)
        A12 = reshape(reshape(A1, D1 * d1, D2) * reshape(A2, D2, d2 * D3), (D1, d1, d2, D3))
        rA12 = reshape(reshape(rhoL * reshape(A12, (D1, d1*d2*D3)), (D1*d1*d2, D3)) * rhoR, (D1, d1, d2, D3))
        HrA12 = permutedims(reshape(reshape(H, (d1 * d2, d1*d2)) * reshape(permutedims(rA12, (2,3,1,4)), (d1 * d2, D1 * D3)), (d1, d2, D1, D3)), (3,1,2,4))
        E = dot(A12, HrA12)
        @cutensor HrA12′[a, s1, s2, c] := rhoL[a, a'] * A1[a', t1, b] * A2[b, t2, c'] * rhoR[c', c] * H[s1, s2, t1, t2]
        @cutensor HrA12′′[:] := rhoL[-1, 1] * H[-2, -3, 4, 5] * A2[2, 5, 3] * rhoR[3, -4] * A1[1, 4, 2] # should be contracted in exactly same order
        @test collect(HrA12′) == collect(HrA12′′) # should be exactly equal
        @test HrA12 ≈ collect(HrA12′)
        @test HrA12 ≈ collect(HrA12′′)
        @test E ≈ @cutensor scalar(rhoL[a', a] * A1[a, s, b] * A2[b, s', c] * rhoR[c, c'] * H[t, t', s, s'] * conj(A1[a', t, b']) * conj(A2[b', t', c']))
    end
    println("tensor network examples: $(time()-t0) seconds")
    t0 = time()
end


text:
using DifferentialEquations, GeometricIntegratorsDiffEq, ParameterizedFunctions, Plots, LinearAlgebra, DelimitedFiles, LaTeXStrings

## 3.1 Determining A Benchmark Value for the Infinite-Time average

include("..\\NewtonCotes\\NewtonCotes11pt.jl")

# Set BigFloat precision. This gives about 425 digits, enough to capture all the precision from Kehlet & Logg
setprecision(1400)

# Extract the z-variable of the reference solution trajectory
# This proved troublesome for accomadating the BigFloat type.
# There may be a cleaner way to do this.
u_ref = readdlm("lorenz\\Kehlet&Logg_solution_u.txt", ' ', String, '\n')
z_ref = Array{BigFloat,1}(undef, size(test)[1])

for i = 1:size(u_ref)[1]
    z_ref[i] = parse(BigFloat, u_ref[i,3])
end

# Integrate and then average over t ∈ [20, 1000]
# Kehlet & Logg's data is provided at time intervals of 0.1
Z_true_Lorenz = NewtonCotes.NewtonCotes11pt(z_ref[201:10001], 0.1)/(1000-20)

## 3.2 Initial Condition Determination

# define the Lorenz system and parameters
lorenz = @ode_def begin
 dx = σ * (y - x)
 dy = x * (ρ - z) - y
 dz = x * y - β*z
end σ ρ β

p = [10.0,28.0,8/3]

# Plots from Fig. 1

probLorenz10a = ODEProblem(lorenz, [1, 0, 0], (0, 10), p)
probLorenz10b = ODEProblem(lorenz, [-1, -1, -1], (0, 10), p)
probLorenz50a = ODEProblem(lorenz, [1, 0, 0], (0, 50), p)
probLorenz50b = ODEProblem(lorenz, [-1, -1, -1], (0, 50), p)

solLorenz10a = solve(probLorenz10a, alg=GIGLRK(4), dt=0.001)
solLorenz10b = solve(probLorenz10b, alg=GIGLRK(4), dt=0.001)
solLorenz50a = solve(probLorenz10a, alg=GIGLRK(4), dt=0.001)
solLorenz50b = solve(probLorenz10b, alg=GIGLRK(4), dt=0.001)

# Figure 1a
plot(solLorenz10a, vars = (1, 2, 3), label=L"u_0 = [1, 0, 0]", linecolor=:blue)
plot!(solLorenz10b, vars = (1, 2, 3), label=L"u_0 = [-1, -1, -1]", linecolor=:orange)

# Figure 1b
plot(solLorenz50a, vars = (1, 2, 3), label=L"u_0 = [1, 0, 0]", linecolor=:blue)
plot!(solLorenz50b, vars = (1, 2, 3), label=L"u_0 = [-1, -1, -1]", linecolor=:orange)

# Figure 2 and new initial condition coordinates
probLorenz22 = ODEProblem(lorenz, [-1, -1, -1], (0, 22), p)
solLorenz22 = solve(probLorenz22, alg=GIGLRK(4), dt=0.001)
plot(solLorenz22, vars = (1, 2, 3), label=L"u_0 = [-1, -1, -1]", linecolor=:orange)

u0Lorenz = solLorenz22.u[end] # = [-5.276131687990789, -8.274832825254473, 16.36301350713599]

## 3.3 Step Size Determination

errors = [0.04 0;
          0.02 0;
          0.01 0;
          0.005 0;
          0.004 0]

probLorenzStep = ODEProblem(lorenz, u0Lorenz, (0, 1), p)

for i=1:size(errors)[1]
    local h = errors[i,1]
    sol1 = solve(probLorenzStep, alg=GIHeun(), dt=h)
    sol2 = solve(probLorenzStep, alg=GIGLRK(4), dt=h)
    errors[i,2] = abs(norm(sol1.u[end] - sol2.u[end])/norm(sol2.u[end]))*100
end

# The errors matrix now shows Table 1

## 3.4 Time-March Testing
# NOTE: At time of initial investigation, the programmer was more experienced with MATLAB
# than Julia. Thus a good portion of the statistical analysis code was written in MATLAB, and
# the workflow for gathering statistics looks like this:
#    1. Run simulations in Julia
#    2. Export simulation results to .txt file
#    3. Import results from .txt file into MATLAB
#    4. Perform statistical analysis and plotting within MATLAB
#
# The code may, in the future, be updated to use Julia entirely. However, the code below shows
# steps 1 and 2 of the above workflow; while step 3 is performed by importer.m, and step 4
# is performed by runner.m and stats_plotter.m
#

probLorenz = ODEProblem(lorenz, u0Lorenz, (0, 300), p)

h = [0.32, 0.16, 0.08, 0.04, 0.02, 0.01, 0.005, 0.0025]
for i = h
    # GLRK stable for all step sizes
    solGLRK = solve(probLorenz, alg=GIGLRK(4), dt=i)
    writedlm("lorenz\\Simulated Solutions\\GLRK(4)_u_"*string(i)*".txt", solGLRK.u)

    # Explicit RK4 stable for h < 0.32
    if i < 0.32
        solRK4 = solve(probLorenz, alg=GIERK4(), dt=i)
        writedlm("lorenz\\Simulated Solutions\\RK4_u_"*string(i)*".txt", solRK4.u)
    end

    # Heun's method stable for h < 0.08
    if i < 0.08
        solHeun = solve(probLorenz, alg=GIHeun(), dt=i)
        writedlm("lorenz\\Simulated Solutions\\Heun_u_"*string(i)*".txt", solHeun.u)
    end
end

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos