text:
include("utils.jl")

mode = "data"

# Global parameters
K_CONST =  9.869232667160130e-16 * 86400 * 1e3
ALPHA = 1.0
GRAV_CONST = 9.8    # gravity constant
SRC_CONST = 86400.0 #
n = 15
m = 30
h = 30.0 # meter
NT  = 50
Δt = 1000/NT
# Δt = 20.0 # day

ρw = 501.9
ρo = 1053.0
μw = 0.1
μo = 1.0
K_init = 20.0 .* ones(n,m) # initial guess of permeability
g = GRAV_CONST
φ0 = 0.25 .* ones(n,m)
qw = zeros(NT+1, n, m)
qw[:,9,3] .= 0.005 * (1/h^2)/10.0 * SRC_CONST
qo = zeros(NT+1, n, m)
qo[:,9,28] .= -0.005 * (1/h^2)/10.0 * SRC_CONST
sw0 = zeros(n, m)
K = 20.0 .* ones(n,m) # millidarcy
K[8:10,:] .= 120.0

# E = 6e9
E = 6.e9
ν = 0.35
D = E/(1+ν)/(1-2ν)*[1-ν ν 0;ν 1-ν 0;0 0 (1-2ν)/2]

Z = zeros(n, m)
for j = 1:n
    for i = 1:m
        Z[j, i] = (j-0.5)*h
    end
end

bdnode = bcnode("left | right", m, n, h)

InteractionM = constant(compute_interaction_matrix(m, n, h)')
StiffM = compute_fem_stiffness_matrix(D, m, n, h)/E
StiffM, _ = fem_impose_Dirichlet_boundary_condition_experimental(StiffM, bdnode, m, n, h)
StiffM = StiffM*E

# pl = placeholder([1.0])
# StiffM = StiffM*pl[1]

qo, qw = constant(qo), constant(qw)
function porosity(u)
    ε = compute_fvm_mechanics_term(u, m, n, h)/h^2
    ε = reshape(ε, (n, m))
    out = 1 - constant(1 .- φ0) .* exp(-ε)
end

################### solid equations
invη = constant(1e-11*ones(4*m*n))
E = 6.e9
ν = 0.35

μ = E/(2(1+ν))
λ = E*ν/(1+ν)/(1-2ν)

fn_G = invη->begin
  G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
    -μ/3*invη 1/Δt+2/3*μ*invη 0.0
    0.0 0.0 1/Δt+μ*invη])
  invG = inv(G)
end
S = map(fn_G, invη)/Δt
Temp = tensor([2μ+λ λ 0.0
    λ 2μ+λ 0.0
    0.0 0.0 μ])
H = S*Temp
Hmat = compute_fem_stiffness_matrix(H, m, n, h)

Hmat_adj, _ = fem_impose_Dirichlet_boundary_condition_experimental(Hmat/E, bdnode, m, n, h)
Hmat_adj = Hmat_adj*E

function solid(Ψ2, σ, uold)

    p = Ψ2 #+ ρw*g*Z
    # p = p - ave_normal(p)
    rhs = InteractionM*reshape(p, (-1,))

    Sσ = batch_matmul(S, σ)
    rhs = rhs + Hmat * uold - compute_strain_energy_term(Sσ, m, n, h)

    mask = ones(2*(m+1)*(n+1))
    mask[[bdnode;bdnode .+ (m+1)*(n+1)]] .= 0.0
    rhs = rhs .* mask

    u = Hmat_adj\rhs

    εold = eval_strain_on_gauss_pts(uold, m, n, h)
    ε = eval_strain_on_gauss_pts(u, m, n, h)
    σ = Sσ + batch_matmul(H, ε - εold)
    # σ2 = batch_matmul(H, ε)
    # op = tf.print(norm(σ-σ2))
    # u = bind(u, op)

    return u, σ
end
###################

################### fluid  equations
function Krw(Sw)
    return Sw ^ 2
end

function Kro(So)
    return So ^ 2
end

function ave_normal(quantity)
    aa = sum(quantity)
    return aa/(m*n)
end

function fluid(i, u, uold, sw, p)
    # step 0: compute porosity and its rate
    φ = porosity(u)
    φold = porosity(uold)
    dotφ = (φ-φold)/Δt

    # φ = constant(φ0)
    # φold = constant(φ0)
    # φ = φ * pl[1]
    # φold = φold * pl[1]
    # dotφ = tf.zeros_like(dotφ)
    # dotφ = constant(zeros(n, m))

    # step 1: update p
    # λw = Krw(sw)/μw
    # λo = Kro(1-sw)/μo
    λw = sw.*sw/μw
    λo = (1-sw).*(1-sw)/μo
    λ = λw + λo
    q = qw[i] + qo[i] + λw/(λo+1e-16).*qo[i]
    # q = qw + qo
    potential_c = (ρw - ρo)*g .* Z

    # Step 2: implicit potential
    Θ = upwlap_op(K * K_CONST, λo, potential_c, h, constant(0.0))

    load_normal = (Θ+q/ALPHA+dotφ) - ave_normal(Θ+q/ALPHA+dotφ)

    # p = poisson_op(λ.*K* K_CONST, load_normal, h, constant(0.0), constant(1))
    p = upwps_op(K * K_CONST, λ, load_normal, tf.zeros_like(p), h, constant(0.0), constant(0))
    # potential p = pw - ρw*g*h

    # step 3: implicit transport
    # sw = sat_op2(sw, dotφ, p, K * K_CONST, φ , qw[i], qo[i], μw, μo, sw, Δt, h)

    # ε = compute_fvm_mechanics_term(u, m, n, h)/h^2
    # ε = reshape(ε,(m, n))'
    # φ1 = 1 - constant(1 .- φ0) .* exp(-ε)

    sw1 = sat_op2(sw, dotφ, p, K * K_CONST , φ, qw[i], qo[i], μw, μo, sw, Δt, h)
    # sw2 = sat_op(sw, p, K * K_CONST , φ, qw[i], qo[i], μw, μo, sw, Δt, h)

    # op = tf.print(i, norm(gradients(sum(sw1),sw)-gradients(sum(sw2), sw)))
    # p = bind(p, op)

    # op = tf.print(i, "*", norm(sw1-sw2))#-gradients(sum(sw1), sw)))
    # p = bind(p, op)

    # op = tf.print(sum(sw2))
    # p = bind(p, op)
    return sw1, p, φ
end
###################

function simulate()
    ta_u, ta_sw, ta_p, ta_φ, ta_σ =
        TensorArray(NT+1), TensorArray(NT+1), TensorArray(NT+1), TensorArray(NT+1), TensorArray(NT+1)
    ta_u = write(ta_u, 1, constant(zeros(2*(m+1)*(n+1))))
    ta_u = write(ta_u, 2, constant(zeros(2*(m+1)*(n+1))))
    ta_sw = write(ta_sw, 1, constant(sw0))
    ta_sw = write(ta_sw, 2, constant(sw0))
    ta_p = write(ta_p, 1, constant(zeros(n, m)))
    ta_p = write(ta_p, 2, constant(zeros(n, m)))
    ta_φ = write(ta_φ, 1, constant(φ0))
    ta_φ = write(ta_φ, 2, constant(φ0))
    ta_σ = write(ta_σ, 1, constant(zeros(4*m*n, 3)))
    ta_σ = write(ta_σ, 2, constant(zeros(4*m*n, 3)))
    i = constant(2, dtype=Int32)
    function condition(i, tas...)
        i <= NT
    end
    function body(i, tas...)
        ta_u, ta_sw, ta_p, ta_φ, ta_σ = tas
        u = read(ta_u, i)
        uold = read(ta_u, i-1)
        σ = read(ta_σ, i)
        sw, p = read(ta_sw, i), read(ta_p, i)
        sw, p, φ = fluid(i, u, uold, sw, p)
        unew, σnew = solid(p, σ, u)
        ta_sw = write(ta_sw, i+1, sw)
        ta_p = write(ta_p, i+1, p)
        ta_u = write(ta_u, i+1, unew)
        ta_φ = write(ta_φ, i+1, φ)
        ta_σ = write(ta_σ, i+1, σnew)
        i+1, ta_u, ta_sw, ta_p, ta_φ, ta_σ
    end
    _, ta_u, ta_sw, ta_p, ta_φ,  ta_σ = while_loop(condition, body, [i, ta_u, ta_sw, ta_p, ta_φ, ta_σ])
    out_u, out_sw, out_p, out_φ, out_σ = stack(ta_u), stack(ta_sw), stack(ta_p), stack(ta_φ), stack(ta_σ)
    return set_shape(out_u, NT+1, 2*(m+1)*(n+1)), set_shape(out_sw, NT+1, n, m),
            set_shape(out_p, NT+1, n, m), set_shape(out_φ, NT+1, n, m),
            set_shape(out_σ, NT+1, 4*m*n, 3)
end

u, S2, Ψ2, φ, σ = simulate()
uobs = u[:, 1:m+1]

if mode!="data"
    uobs_ = readdlm("obs.txt")
    global loss = mean((uobs-uobs_)^2)
end

sess = Session(); init(sess)
if mode=="data"
    writedlm("obs.txt", run(sess, uobs))

    global o1, o2, o3, o4, o5 = run(sess, [u, S2, Ψ2, φ, σ])
    error("Stop")
end

BFGS!(sess, loss)

# # step 1
# @show run(sess, loss)
# # gd = gradients(loss, pl)
# # @show run(sess, gd)
# # step 2
# lineview(sess, pl, loss, [1.0;1.0], [5.0;1.5])
# gradview(sess, pl, loss, [5.0;1.5])

# #
# # visualize_potential(o4)
visualize_displacement(10*o1[end,:], m, n, h); savefig("two_u.pdf")
visualize_saturation(o2[1,:,:], m, n, h); savefig("two_s0.pdf")
visualize_saturation(o2[25,:,:], m, n, h); savefig("two_s25.pdf")
visualize_saturation(o2[51,:,:], m, n, h); savefig("two_s51.pdf")
# # visualize_potential(o3)

plot(o1[:,1:(m+1)])

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
#==============================================================================================#
#                              runtests.jl -- TypeTree test file                               #
#==============================================================================================#

using TypeTree
using Test

#==============================================================================================#
#                                          Test Data                                           #
#==============================================================================================#

module M

#----------------------------------------------------------------------------------------------#
#                                        Abstract types                                        #
#----------------------------------------------------------------------------------------------#

export Ab

abstract type tyTree <: Any end
abstract type A <: tyTree end
abstract type Aa <: A end
abstract type Ab <: A end
abstract type B <: tyTree end
abstract type Ba <: B end

#----------------------------------------------------------------------------------------------#
#                                        Concrete types                                        #
#----------------------------------------------------------------------------------------------#

export Aα

struct Aα <: A end
struct Aaα <: Aa end
struct Aaβ <: Aa end
struct Abα <: Ab end

end # module M

#==============================================================================================#
#                                            Tests                                             #
#==============================================================================================#

using Main.M

@testset "quiz.type.abstract:                                                     " begin
    # All types
    @test length(tt(M.A)) == 7
    @test length(tt(M.Aa)) == 3
    @test length(tt(Ab)) == 2
    @test length(tt(M.B)) == 2
    @test length(tt(M.Ba)) == 1
    @test length(tt(M.tyTree)) == 10
    # Abstract types only
    @test length(tt(M.A, concrete=false)) == 3
    @test length(tt(M.Aa, concrete=false)) == 1
    @test length(tt(Ab, concrete=false)) == 1
    @test length(tt(M.B, concrete=false)) == 2
    @test length(tt(M.Ba, concrete=false)) == 1
    @test length(tt(M.tyTree, concrete=false)) == 6

    @test contains(join(tt(M.tyTree)), " Main.M.Aa\n")
    @test contains(join(tt(M.tyTree)), " Ab\n")
    @test contains(join(tt(M.tyTree)), " Aα\n")

    @test !contains(join(tt(M.tyTree)), " Aa\n")
    @test !contains(join(tt(M.tyTree)), " Main.M.Ab\n")
    @test !contains(join(tt(M.tyTree)), " Main.M.Aα\n")

    @test contains(join(tt(M.tyTree; mod=Main.M)), " Aa\n")
    @test contains(join(tt(M.tyTree; mod=Main.M)), " Ab\n")
    @test contains(join(tt(M.tyTree; mod=Main.M)), " Aα\n")

    @test !contains(join(tt(M.tyTree; mod=Main.M)), " Main.M.Aa\n")
    @test !contains(join(tt(M.tyTree; mod=Main.M)), " Main.M.Ab\n")
    @test !contains(join(tt(M.tyTree; mod=Main.M)), " Main.M.Aα\n")
end

ext: jl


text:
using BenchmarkTools

module test
  using Dates
  using Future: randjump
  import Random
  const N = 10^4

  function f0()
    s = 0.
    for i in 1:N
      s += rand()
    end
    return s
  end

  function f1()
    RNG1 = Random.MersenneTwister()
    s = 0.
    for i in 1:N
      s += rand(RNG1)
    end
    return s
  end

  const RNG = Random.MersenneTwister()
  function f2()
    s = 0.
    for i in 1:N
      s += rand(RNG)
    end
    return s
  end

  const RNGS = [randjump(Random.MersenneTwister(),big(10)^20)]
  init_random() = foreach(_ -> push!(RNGS, randjump(last(RNGS),big(10)^20)), 2:Threads.nthreads())
  random() = rand(RNGS[Threads.threadid()])
  random(arg) = rand(RNGs[Threads.threadid()],arg)

  function g3()
    return random()
  end
  function f3()
    s = 0.
    for i in 1:N
      s += g3()
    end
    return s
  end

  function t(func)
    nthreads = Threads.nthreads()
    dates = Vector{DateTime}(undef,nthreads)
    task = Vector{Task}(undef,nthreads)
    for i in 1:nthreads
      task[i] = Threads.@spawn func()
      dates[i] = Dates.now()
    end
    s = 0.
    for i in 1:nthreads
      s += fetch(task[i])
    end
    #println(s)
  end

end

test.init_random()

println("f0: rand()")
@btime test.t(test.f0)

println("f1: RNG inside f")
@btime test.t(test.f1)

println("f2: using const RNG = ...")
@btime test.t(test.f2)

println("f3: using const RNGS = ...")
@btime test.t(test.f3)

println("init_random():")
@btime test.init_random()


text:
# import Base: getindex, size, *

struct SFMDistribution{T<:SFFM.Mesh} <: AbstractArray{Float64,2}
    coeffs::Array{Float64,2}
    model::SFFM.Model
    mesh::T
    Fil::SFFM.IndexDict
    function SFMDistribution{T}(
        coeffs::Array{Float64,2},
        model::SFFM.Model,
        mesh::T,
        Fil::SFFM.IndexDict=SFFM.MakeFil(model,mesh.Nodes),
        ) where T<:SFFM.Mesh
        return (size(coeffs,1)==1) ? new(coeffs,model,mesh,Fil) : throw(DimensionMismatch("coeffs must be a row-vector"))
    end
end

size(d::SFMDistribution) = size(d.coeffs)
getindex(d::SFMDistribution,i::Int,j::Int) = d.coeffs[i,j]
setindex!(d::SFMDistribution,x,i::Int,j::Int) = throw(DomainError("inserted value(s) must be Float64/Int"))
setindex!(d::SFMDistribution,x::Float64,i::Int,j::Int) = (d.coeffs[i,j]=x)
setindex!(d::SFMDistribution,x::Int,i::Int,j::Int) = (d.coeffs[i,j]=convert(Float64,x))
*(u::SFMDistribution,B::AbstractArray{Float64,2}) = *(u.coeffs,B)
*(B::AbstractArray{Float64,2},u::SFMDistribution) =
    (size(B,2)==1) ? *(B*u.coeffs) : throw(DimensionMismatch("u is a row-vector and B has more than 1 column"))

function _error_on_nothing(idx)
    (idx===nothing) && throw(DomainError("x is not in the support of the mesh"))
end

function _get_nodes_coeffs_from_index(cell_idx::Int,i::Int,mesh::Mesh,model::Model)
    cell_nodes = CellNodes(mesh)[:,cell_idx]
    N₋ = sum(model.C.<=0)
    coeff_idx = (N₋ + (i-1)*TotalNBases(mesh) + (cell_idx-1)*NBases(mesh)) .+ (1:NBases(mesh))
    return cell_nodes, coeff_idx
end

function _get_coeff_index_pos(x::Float64,i::Int,mesh::Mesh,model::Model)
    cell_idx = findlast(x.>=mesh.Nodes)
    _error_on_nothing(cell_idx)
    cell_nodes, coeff_idx = _get_nodes_coeffs_from_index(cell_idx,i,mesh,model)
    return cell_idx, cell_nodes, coeff_idx
end
function _get_coeff_index_neg(x::Float64,i::Int,mesh::Mesh,model::Model)
    cell_idx = findfirst(x.<=mesh.Nodes) - 1
    _error_on_nothing(cell_idx)
    cell_nodes, coeff_idx = _get_nodes_coeffs_from_index(cell_idx,i,mesh,model)
    return cell_idx, cell_nodes, coeff_idx
end

function _get_point_mass_data_pos(i::Int,mesh::Mesh,model::Model)
    cell_nodes = mesh.Nodes[end]
    N₋ = sum(model.C.<=0)
    coeff_idx = N₋ + TotalNBases(mesh)*NPhases(model) + sum(model.C[1:i].>=0)
    return cell_nodes, coeff_idx
end
function _get_point_mass_data_neg(i::Int,mesh::Mesh,model::Model)
    cell_nodes = mesh.Nodes[1]
    coeff_idx = sum(model.C[1:i].<=0)
    return cell_nodes, coeff_idx
end
_is_left_point_mass(x::Float64,i::Int,mesh::Mesh,model::Model) =
    (x==mesh.Nodes[1])&&(model.C[i]<=0)
_is_right_point_mass(x::Float64,i::Int,mesh::Mesh,model::Model) =
    (x==mesh.Nodes[end])&&(model.C[i]>=0)

function _get_coeffs_index(x::Float64,i::Int,model::Model,mesh::Mesh)
    !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))

    # find which inteval Dₗ,ᵢ x is in
    if _is_left_point_mass(x,i,mesh,model)
        cell_idx = "point mass"
        cell_nodes, coeff_idx = _get_point_mass_data_neg(i,mesh,model)
    elseif _is_right_point_mass(x,i,mesh,model)
        cell_idx = "point mass"
        cell_nodes, coeff_idx = _get_point_mass_data_pos(i,mesh,model)
    else # not a point mass
        if model.C[i]>=0
            cell_idx, cell_nodes, coeff_idx = _get_coeff_index_pos(x,i,mesh,model)
        elseif model.C[i]<0
            cell_idx, cell_nodes, coeff_idx = _get_coeff_index_neg(x,i,mesh,model)
        end
    end

    return cell_idx, cell_nodes, coeff_idx
end

function legendre_to_lagrange(coeffs)
    order = length(coeffs)
    V = vandermonde(NBases(mesh))
    return V.V*coeffs
end

function pdf(d::SFMDistribution{T},model::Model) where T<:Mesh
    throw(DomainError("unknown SFMDistribution{<:Mesh}"))
end

function pdf(d::SFMDistribution{DGMesh},model::Model)
    function f(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.mesh
        if ((x<mesh.Nodes[1])||(x>mesh.Nodes[end]))
            fxi = 0.0
        else
            cell_idx, cell_nodes, coeff_idx = _get_coeffs_index(x,i,model,mesh)
            coeffs = d.coeffs[coeff_idx]
            # if not a point mass, then reconstruct solution
            if !(cell_idx=="point mass")
                if Basis(mesh) == "legendre"
                    coeffs = legendre_to_lagrange(coeffs)
                else
                    V = vandermonde(NBases(mesh))
                    coeffs = (2/(Δ(mesh)[cell_idx]))*(1.0./V.w).*coeffs
                end
                basis_values = lagrange_poly_basis(cell_nodes, x)
                fxi = LinearAlgebra.dot(basis_values,coeffs)
            else
                fxi = coeffs
            end
        end
        return fxi
    end
    return f
end
pdf(d::SFMDistribution{T},model::Model,x,i) where T<:Mesh =
    throw(DomainError("x must be Float64/Int/Array{Float64/Int,1}, i must be Int/Array{Int,1}"))
pdf(d::SFMDistribution{T},model::Model,x::Float64,i::Int) where T<:Mesh = pdf(d,model)(x,i)
pdf(d::SFMDistribution{T},model::Model,x::Int,i::Int) where T<:Mesh = pdf(d,model)(convert(Float64,x),i)
# pdf(d::SFMDistribution{T},model::Model,x::Array{Float64,1},i::Int) where T<:Mesh = pdf(d,model).(x,i)
# pdf(d::SFMDistribution{T},model::Model,x::Array{Float64,1},i::Array{Int,1}) where T<:Mesh = pdf(d,model).(x,i)

function pdf(d::SFMDistribution{FRAPMesh},model::Model)
    function f(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.mesh
        if ((x<mesh.Nodes[1])||(x>mesh.Nodes[end]))
            fxi = 0.0
        else
            cell_idx, cell_nodes, coeff_idx = _get_coeffs_index(x,i,model,mesh)
            coeffs = d.coeffs[coeff_idx]
            # if not a point mass, then reconstruct solution
            if !(cell_idx=="point mass")
                if model.C[i]>0
                    yₖ₊₁ = mesh.Nodes[cell_idx+1]
                    to_go = yₖ₊₁-x
                elseif model.C[i]<0
                    yₖ = mesh.Nodes[cell_idx]
                    to_go = x-yₖ
                end
                me = MakeME(CMEParams[NBases(mesh)], mean = Δ(mesh)[cell_idx])
                fxi = (pdf(Array(coeffs'),me,to_go) + pdf(Array(coeffs'),me,2*Δ(mesh)[cell_idx]-to_go))./cdf(Array(coeffs'),me,2*Δ(mesh)[cell_idx])
            else
                fxi = coeffs
            end
        end
        return fxi
    end
    return f
end

function pdf(d::SFMDistribution{FVMesh},model::Model)
    function f(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.mesh
        if ((x<mesh.Nodes[1])||(x>mesh.Nodes[end]))
            fxi = 0.0
        else
            cell_idx, cell_nodes, coeff_idx = _get_coeffs_index(x,i,model,mesh)
            # if not a point mass, then reconstruct solution
            if !(cell_idx=="point mass")
                ptsLHS = Int(ceil(Order(mesh)/2))
                if cell_idx-ptsLHS < 0
                    nodesIdx = 1:Order(mesh)
                    nodes = CellNodes(mesh)[nodesIdx]
                    poly_vals = lagrange_poly_basis(nodes,x)
                elseif cell_idx-ptsLHS+Order(mesh) > TotalNBases(mesh)
                    nodesIdx = (TotalNBases(mesh)-Order(mesh)+1):TotalNBases(mesh)
                    nodes = CellNodes(mesh)[nodesIdx]
                    poly_vals = lagrange_poly_basis(nodes,x)
                else
                    nodesIdx =  (cell_idx-ptsLHS) .+ (1:Order(mesh))
                    poly_vals = lagrange_poly_basis(CellNodes(mesh)[nodesIdx],x)
                end
                coeff_idx = (sum(model.C.<=0) + (i-1)*TotalNBases(mesh)) .+ nodesIdx
                coeffs = d.coeffs[coeff_idx]#./Δ(mesh)[cell_idx]
                fxi = LinearAlgebra.dot(poly_vals,coeffs)
            else
                coeffs = d.coeffs[coeff_idx]
                fxi = coeffs
            end
        end
        return fxi
    end
    return f
end

############
### CDFs ###
############

function cdf(d::SFMDistribution{T},model::Model) where T<:Mesh
    throw(DomainError("unknown SFMDistribution{<:Mesh}"))
end

"""

    _sum_cells_left(d::SFMDistribution, i::Int, cell_idx::Int, mesh::Mesh, model::Model)

Add up all the probability mass in phase `i` in the cells to the left of `cell_idx`.
"""
function _sum_cells_left(d::SFMDistribution, i::Int, cell_idx::Int, mesh::Mesh, model::Model)
    c = 0.0
    if Basis(mesh) == "legendre"
        for cell in 1:(cell_idx-1)
            # first legendre basis function =1 & has all the mass
            idx = (sum(model.C.<=0) + (i-1)*TotalNBases(mesh) + (cell-1)*NBases(mesh)) .+ 1
            c += d.coeffs[idx]
        end
    else
        for cell in 1:(cell_idx-1)
            idx = (sum(model.C.<=0) + (i-1)*TotalNBases(mesh) + (cell-1)*NBases(mesh)) .+ (1:NBases(mesh))
            c += sum(d.coeffs[idx])
        end
    end
    return c
end

function cdf(d::SFMDistribution{DGMesh},model::Model)
    function F(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        mesh = d.mesh
        Fxi = 0.0
        if (x<mesh.Nodes[1])
            # Fxi = 0.0
        else
            # Fxi = 0.0
            # left pm
            if (x>=mesh.Nodes[1])&&(model.C[i]<=0)
                ~, left_pm_idx = _get_point_mass_data_neg(i,mesh,model)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            (x.>=mesh.Nodes[end]) ? (xd=mesh.Nodes[end]-sqrt(eps())) : xd = x
            cell_idx, ~, ~ = _get_coeffs_index(xd,i,model,mesh)
            if !(cell_idx=="point mass")
                # add all mass from cells to the left
                Fxi += _sum_cells_left(d, i, cell_idx, mesh, model)

                # integrate up to x in the cell which contains x
                temp_pdf(y) = pdf(d,model)(y,i)
                quad = gauss_lobatto_quadrature(temp_pdf,mesh.Nodes[cell_idx]+sqrt(eps()),xd,NBases(mesh))
                Fxi += quad
            end
            # add the RH point mass if  required
            if (x>=mesh.Nodes[end])&&(model.C[i]>=0)
                ~, ~, right_pm_idx = _get_coeffs_index(mesh.Nodes[end],i,model,mesh)
                right_pm = d.coeffs[right_pm_idx]
                Fxi += right_pm
            end
        end
        return Fxi
    end
    return F
end
cdf(d::SFMDistribution{T},model::Model,x,i) where T<:Mesh =
    throw(DomainError("x must be Float64/Int/Array{Float64/Int,1}, i must be Int/Array{Int,1}"))
cdf(d::SFMDistribution{T},model::Model,x::Float64,i::Int) where T<:Mesh = cdf(d,model)(x,i)
cdf(d::SFMDistribution{T},model::Model,x::Int,i::Int) where T<:Mesh = cdf(d,model)(convert(Float64,x),i)

function cdf(d::SFMDistribution{FRAPMesh},model::Model)
    function F(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.mesh
        if (x<mesh.Nodes[1])
            Fxi = 0.0
        else
            Fxi = 0.0
            # left pm
            if (x>=mesh.Nodes[1])&&(model.C[i]<=0)
                ~, ~, left_pm_idx = _get_coeffs_index(mesh.Nodes[1],i,model,mesh)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            (x.>=mesh.Nodes[end]) ? (xd=mesh.Nodes[end]-sqrt(eps())) : xd = x
            cell_idx, cell_nodes, coeff_idx = _get_coeffs_index(xd,i,model,mesh)
            coeffs = d.coeffs[coeff_idx]

            if !(cell_idx=="point mass")
                # add all mass from cells to the left
                Fxi += _sum_cells_left(d, i, cell_idx, mesh, model)

                # integrate up to x in the cell which contains x
                me = MakeME(CMEParams[NBases(mesh)], mean = Δ(mesh)[cell_idx])
                a = Array(coeffs')
                if model.C[i]>=0
                    yₖ₊₁ = mesh.Nodes[cell_idx+1]
                    Fxi += (ccdf(a,me,yₖ₊₁-x) - ccdf(a,me,2*Δ(mesh)[cell_idx]-(yₖ₊₁-x)))/(cdf(a,me,2*Δ(mesh)[cell_idx]))
                elseif model.C[i]<0
                    yₖ = mesh.Nodes[cell_idx]
                    Fxi += sum(a) - (ccdf(a,me,x-yₖ) - ccdf(a,me,2*Δ(mesh)[cell_idx]-(x-yₖ)))/(cdf(a,me,2*Δ(mesh)[cell_idx]))
                end
            end
            if (x>=mesh.Nodes[end])&&(model.C[i]>=0)
                ~, ~, right_pm_idx = _get_coeffs_index(mesh.Nodes[end],i,model,mesh)
                right_pm = d.coeffs[right_pm_idx]
                Fxi += right_pm
            end
        end
        return Fxi
    end
    return F
end

function _sum_cells_left_fv(d, i, cell_idx, mesh, model)
    c = 0
    for cell in 1:(cell_idx-1)
        # first legendre basis function =1 & has all the mass
        idx = sum(model.C.<=0) + (i-1)*TotalNBases(mesh) + cell
        c += d.coeffs[idx]*Δ(mesh)[cell]
    end
    return c
end

function cdf(d::SFMDistribution{FVMesh},model::Model)
    function F(x::Float64,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(model)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.mesh
        Fxi = 0.0
        if (x<mesh.Nodes[1])
            # Fxi = 0.0
        else
            # Fxi = 0.0
            # left pm
            if (x>=mesh.Nodes[1])&&(model.C[i]<=0)
                ~, ~, left_pm_idx = _get_coeffs_index(mesh.Nodes[1],i,model,mesh)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            (x.>=mesh.Nodes[end]) ? (xd=mesh.Nodes[end]-sqrt(eps())) : xd = x
            cell_idx, cell_nodes, coeff_idx = _get_coeffs_index(xd,i,model,mesh)
            # if not a point mass, then reconstruct solution
            if !(cell_idx=="point mass")
                # add all mass from cells to the left
                Fxi += _sum_cells_left_fv(d, i, cell_idx, mesh, model)

                # integrate up to x in the cell which contains x
                temp_pdf(y) = pdf(d,model)(y,i)
                quad = gauss_lobatto_quadrature(temp_pdf,mesh.Nodes[cell_idx]+sqrt(eps()),xd,Order(mesh))
                Fxi += quad
            end
            if (x>=mesh.Nodes[end])&&(model.C[i]>=0)
                ~, ~, right_pm_idx = _get_coeffs_index(mesh.Nodes[end],i,model,mesh)
                right_pm = d.coeffs[right_pm_idx]
                Fxi += right_pm
            end
        end
        return Fxi
    end
    return F
end

abstract type SFFMDistribution end

"""

    SFFMDistribution(
        pm::Array{<:Real},
        distribution::Array{<:Real,3},
        x::Array{<:Real},
        type::String,
    )

- `pm::Array{Float64}`: a vector containing the point masses, the first
    `sum(model.C.<=0)` entries are the left hand point masses and the last
    `sum(model.C.>=0)` are the right-hand point masses.
- `distribution::Array{Float64,3}`: "probability" or "density"`
- `x::Array{Float64,2}`:
    - if `type="probability"` is a `1×NIntervals×NPhases` array
        containing the cell centers.
    - if `type="density"` is a `NBases×NIntervals×NPhases` array
        containing the cell nodes at which the denisty is evaluated.
- `type::String`: either `"probability"` or `"density"`. `"cumulative"` is
    not possible.
"""
struct SFFMDensity <: SFFMDistribution
    pm::Array{<:Real}
    distribution::Array{<:Real,3}
    x::Array{<:Real}
end
struct SFFMProbability <: SFFMDistribution
    pm::Array{<:Real}
    distribution::Array{<:Real,3}
    x::Array{<:Real}
end
struct SFFMCDF <: SFFMDistribution
    pm::Array{<:Real}
    distribution::Array{<:Real,3}
    x::Array{<:Real}
end

"""
Convert from a vector of coefficients for the DG system to a distribution.

    Coeffs2Dist(
        model::SFFM.Model,
        mesh::SFFM.Mesh,
        Coeffs;
        type::String = "probability",
    )

# Arguments
- `model`: a Model object
- `mesh`: a Mesh object as output from MakeMesh
- `Coeffs::Array`: a vector of coefficients from the DG method
- `type::String`: an (optional) declaration of what type of distribution you
    want to convert to. Options are `"probability"` to return the probabilities
    ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``is the kth cell, `"cumulative"` to
    return the CDF evaluated at cell edges, or `"density"` to return an
    approximation to the density ar at the SFFM.CellNodes(mesh).

# Output
- a tuple with keys
(pm=pm, distribution=yvals, x=xvals, type=type)
    - `pm::Array{Float64}`: a vector containing the point masses, the first
        `sum(model.C.<=0)` entries are the left hand point masses and the last
        `sum(model.C.>=0)` are the right-hand point masses.
    - `distribution::Array{Float64,3}`:
        - if `type="cumulative"` returns a `2×NIntervals×NPhases` array
            containing the CDF evaluated at the cell edges as contained in
            `x` below. i.e. `distribution[1,:,i]` returns the cdf at the
            left-hand edges of the cells in phase `i` and `distribution[2,:,i]`
            at the right hand edges.
        - if `type="probability"` returns a `1×NIntervals×NPhases` array
            containing the probabilities ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``
            is the kth cell.
        - if `type="density"` returns a `NBases×NIntervals×NPhases` array
            containing the density function evaluated at the cell nodes as
            contained in `x` below.
    - `x::Array{Float64,2}`:
        - if `type="cumulative"` returns a `2×NIntervals×NPhases` array
            containing the cell edges as contained. i.e. `x[1,:]`
            returns the left-hand edges of the cells and `x[2,:]` at the
            right-hand edges.
        - if `type="probability"` returns a `1×NIntervals×NPhases` array
            containing the cell centers.
        - if `type="density"` returns a `NBases×NIntervals×NPhases` array
            containing the cell nodes.
    - `type`: as input in arguments.
"""
function Coeffs2Dist(
    model::SFFM.Model,
    mesh::SFFM.DGMesh,
    Coeffs::AbstractArray,
    type::Type{T} = SFFMProbability,
    v::Bool = false,
) where {T<:SFFMDistribution}

    V = SFFM.vandermonde(NBases(mesh))
    N₋ = sum(model.C .<= 0)
    N₊ = sum(model.C .>= 0)

    if type == SFFMDensity
        xvals = SFFM.CellNodes(mesh)
        if Basis(mesh) == "legendre"
            yvals = reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model))
            for i in 1:NPhases(model)
                yvals[:,:,i] = V.V * yvals[:,:,i]
            end
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            yvals =
                Coeffs[N₋+1:end-N₊] .* repeat(1.0 ./ V.w, NIntervals(mesh) * NPhases(model)) .*
                (repeat(2.0 ./ Δ(mesh), 1, NBases(mesh) * NPhases(model))'[:])
            yvals = reshape(yvals, NBases(mesh), NIntervals(mesh), NPhases(model))
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
        if NBases(mesh) == 1
            yvals = [1;1].*yvals
            xvals = [SFFM.CellNodes(mesh)-Δ(mesh)'/2;SFFM.CellNodes(mesh)+Δ(mesh)'/2]
        end
    elseif type == SFFMProbability
        if NBases(mesh) > 1
            xvals = SFFM.CellNodes(mesh)[1, :] + (Δ(mesh) ./ 2)
        else
            xvals = SFFM.CellNodes(mesh)
        end
        if Basis(mesh) == "legendre"
            yvals = (reshape(Coeffs[N₋+1:NBases(mesh):end-N₊], 1, NIntervals(mesh), NPhases(model)).*Δ(mesh)')./sqrt(2)
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            yvals = sum(
                reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
                dims = 1,
            )
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
    elseif type == SFFMCDF
        if NBases(mesh) > 1
            xvals = SFFM.CellNodes(mesh)[[1;end], :]
        else
            xvals = [SFFM.CellNodes(mesh)-Δ(mesh)'/2;SFFM.CellNodes(mesh)+Δ(mesh)'/2]
        end
        if Basis(mesh) == "legendre"
            tempDist = (reshape(Coeffs[N₋+1:NBases(mesh):end-N₊], 1, NIntervals(mesh), NPhases(model)).*Δ(mesh)')./sqrt(2)
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            tempDist = sum(
                reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
                dims = 1,
            )
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
        tempDist = cumsum(tempDist,dims=2)
        temppm = zeros(Float64,1,2,NPhases(model))
        temppm[:,1,model.C.<=0] = pm[1:N₋]
        temppm[:,2,model.C.>=0] = pm[N₊+1:end]
        yvals = zeros(Float64,2,NIntervals(mesh),NPhases(model))
        yvals[1,2:end,:] = tempDist[1,1:end-1,:]
        yvals[2,:,:] = tempDist
        yvals = yvals .+ reshape(temppm[1,1,:],1,1,NPhases(model))
        pm[N₋+1:end] = pm[N₋+1:end] + yvals[end,end,model.C.>=0]
    end

    out = type(pm, yvals, xvals)
    v && println("UPDATE: distribution object created with keys ", fieldnames(type))
    return out
end
function Coeffs2Dist(
    model::SFFM.Model,
    mesh::Union{FRAPMesh, FVMesh},
    Coeffs::AbstractArray,
    type::Type{T} = SFFMProbability,
    v::Bool = false,
) where {T<:SFFMDistribution}

    if type != SFFMProbability
        args = [
            model;
            mesh;
            Coeffs;
            type;
        ]
        error("Input Error: no functionality other than 'probability' implemented, yet...")
    end

    N₋ = sum(model.C .<= 0)
    N₊ = sum(model.C .>= 0)

    xvals = SFFM.CellNodes(mesh)

    yvals = sum(
        reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
        dims = 1,
    )
    pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]

    out = type(pm, yvals, xvals)
    v && println("UPDATE: distribution object created with keys ", fieldnames(type))
    return out
end

"""
Converts a distribution as output from `Coeffs2Dist()` to a vector of DG
coefficients.

    Dist2Coeffs(
        model::SFFM.Model,
        mesh::SFFM.Mesh,
        Distn::SFFMDistribution,
    )

# Arguments
- `model`: a Model object
- `mesh`: a Mesh object as output from MakeMesh
- `Distn::SFFMDistribution
    - if `type="probability"` is a `1×NIntervals×NPhases` array containing
        the probabilities ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``
        is the kth cell.
    - if `type="density"` is a `NBases×NIntervals×NPhases` array containing
        either the density function evaluated at the cell nodes which are in
        `x` below, or, the inner product of the density function against the
        lagrange polynomials.

# Output
- `coeffs` a row vector of coefficient values of length
    `TotalNBases*NPhases + N₋ + N₊` ordered according to LH point masses, RH
    point masses, interior basis functions according to basis function, cell,
    phase. Used to premultiply operators such as B from `MakeB()`
"""
function Dist2Coeffs(
    model::SFFM.Model,
    mesh::SFFM.DGMesh,
    Distn::SFFMDensity,
)
    V = SFFM.vandermonde(NBases(mesh))
    theDistribution =
        zeros(Float64, NBases(mesh), NIntervals(mesh), NPhases(model))
    if Basis(mesh) == "legendre"
        theDistribution = Distn.distribution
        for i = 1:NPhases(model)
            theDistribution[:, :, i] = V.inv * theDistribution[:, :, i]
        end
    elseif Basis(mesh) == "lagrange"
        theDistribution .= Distn.distribution
        # convert to probability coefficients by multiplying by the
        # weights in V.w/2 and cell widths Δ
        theDistribution = ((V.w .* theDistribution).*(Δ(mesh) / 2)')[:]
    end
    # also put the point masses on the ends
    coeffs = [
        Distn.pm[1:sum(model.C .<= 0)]
        theDistribution[:]
        Distn.pm[sum(model.C .<= 0)+1:end]
    ]
    coeffs = Matrix(coeffs[:]')
    return coeffs
end

function Dist2Coeffs(
    model::SFFM.Model,
    mesh::SFFM.DGMesh,
    Distn::SFFMProbability,
)
    V = SFFM.vandermonde(NBases(mesh))
    theDistribution =
        zeros(Float64, NBases(mesh), NIntervals(mesh), NPhases(model))
    if Basis(mesh) == "legendre"
        # for the legendre basis the first basis function is ϕ(x)=Δ√2 and
        # all other basis functions are orthogonal to this. Hence, we map
        # the cell probabilities to the first basis function only.
        theDistribution[1, :, :] = Distn.distribution./Δ(mesh)'.*sqrt(2)
    elseif Basis(mesh) == "lagrange"
        theDistribution .= Distn.distribution
        # convert to probability coefficients by multiplying by the
        # weights in V.w/2
        theDistribution = (V.w .* theDistribution / 2)[:]
    end
    # also put the point masses on the ends
    coeffs = [
        Distn.pm[1:sum(model.C .<= 0)]
        theDistribution[:]
        Distn.pm[sum(model.C .<= 0)+1:end]
    ]
    coeffs = Matrix(coeffs[:]')
    return coeffs
end
function Dist2Coeffs(
    model::SFFM.Model,
    mesh::Union{SFFM.FRAPMesh,SFFM.FVMesh},
    Distn::SFFMDistribution
)

    # also put the point masses on the ends
    coeffs = [
        Distn.pm[1:sum(model.C .<= 0)]
        Distn.distribution[:]
        Distn.pm[sum(model.C .<= 0)+1:end]
    ]

    coeffs = Matrix(coeffs[:]')
    return coeffs
end

"""
Computes the error between distributions.

    starSeminorm(d1::SFFMProbability, d2::SFFMProbability)

# Arguments
- `d1`: a distribution object as output from `Coeffs2Dist`
- `d2`: a distribution object as output from `Coeffs2Dist`
"""
function starSeminorm(d1::SFFMProbability, d2::SFFMProbability)
    e = sum(abs.(d1.pm-d2.pm)) + sum(abs.(d1.distribution-d2.distribution))
    return e
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
const Compiler = Core.Compiler
using .Compiler: ReturnNode, argextype, SSAValue, ⊑, widenconst, userefs, LineInfoNode, GotoNode,
    IRCode, GotoIfNot, PiNode, PhiNode, foreachssa, block_for_inst, insert_node!,
    dominates, dominated, DominatedBlocks
using InteractiveUtils
using ..XLA: HloParameter, HloOp, XLAComputationConfig, HloModuleProto, HloProto,
    HloSnapshot, XLAComputation, HloMap, HloGetTupleElement, XLAScalar,
    HloReduceWindow, HloReduce, HloTuple

Base.iterate(d::Compiler.DominatedBlocks, state...) = Compiler.iterate(d, state...)

function grab_ir(f, argtypes)
    method = @which f(W, x, b)
    sig = typeof((f, W, x, b))
    mi = Compiler.code_for_method(method, sig, Core.svec(), params.world)
    sv = Compiler.OptimizationState(mi, params)
    ir = Compiler.inflate_ir(ci, mi)
end

macro check_ir(cond, args...)
    @assert length(args) <= 1
    reason = length(args) == 0 ? nothing : esc(args[1])
    quote
        if !$(esc(cond))
            throw(NotOffloadableError($(esc(:ir)), $(esc(:sv)), $(reason)))
        end
    end
end

include("outlining.jl")

function process_function_argument(argvals, sig, idx, ir, stmt, sparams)
    sz = Compiler.isconstType(sig.parameters[idx]) ? 0 : sizeof(sig.parameters[idx])
    if sz != 0
        # If the function being called has non-zero information, we need to make sure it's a constant
        # In the future, XLA will allow this.
        fty = argextype(stmt.args[2+idx], ir, sparams)
        if !isa(fty, Const)
            error("At the moment applied computations must not capture non-constant values")
        end
        argvals = argvals === nothing ? Vector{Any}(undef, 1) : argvals
        argvals[1] = fty.val
    end
    argvals
end

function is_known_call(e::Expr, @nospecialize(func), src, spvals::Core.SimpleVector, slottypes::Vector{Any} = Compiler.empty_slottypes)
    if e.head !== :call
        return false
    end
    f = Compiler.argextype(e.args[1], src, spvals, slottypes)
    return (isa(f, Const) && f.val === func) ||
           (isa(f, DataType) && isdefined(f, :instance) && f.instance === func)
end

function _compile_to_xla!(computations, comp, ir, sv)
    ir = outline_control_flow!(ir, sv)
    arg_instrs = Vector{HloInstructionProto}(undef, length(ir.argtypes))
    xla_args = Type[]
    sparams = sv === nothing ? Core.svec() : sv.sp
    for i = 1:length(ir.argtypes)
        # TODO: We could represent this as a kConstant
        isa(ir.argtypes[i], Const) && continue
        if sv !== nothing && sv.linfo.def.isva && i == length(ir.argtypes)
            # To XLA, we express varargs as separate parameters, so gather all of
            # them and then put them in a tuple
            vargs = HloInstructionProto[]
            for T in ir.argtypes[i].parameters
                instr = HloInstructionProto(comp, "parameter")
                instr.shape = dtype_to_shape(T)
                instr.parameter_number = length(xla_args)
                push!(vargs, instr)
                push!(xla_args, T)
            end
            arg_instrs[i] = HloInstructionProto(comp, HloTuple(), vargs...)
        elseif ir.argtypes[i] ⊑ XRTArray
            AT = widenconst(ir.argtypes[i])
            eltype = AT.parameters[1]
            dims = AT.parameters[2]
            arg_instrs[i] = HloInstructionProto(comp, HloParameter(eltype, dims, length(xla_args)))
            push!(xla_args, ir.argtypes[i])
        elseif representable(ir.argtypes[i])
            T = widenconst(ir.argtypes[i])
            instr = HloInstructionProto(comp, "parameter")
            instr.shape = dtype_to_shape(T)
            instr.parameter_number = length(xla_args)
            arg_instrs[i] = instr
            push!(xla_args, T)
        end
    end
    ssa_vals = Vector{HloInstructionProto}(undef, length(ir.stmts))
    function hlo_eval(arg)
        if isa(arg, Compiler.Argument)
            return arg_instrs[arg.n]
        elseif isa(arg, SSAValue)
            @check_ir isassigned(ssa_vals, arg.id) "SSA value $(arg) not mapped to HLO"
            return ssa_vals[arg.id]
        else
            @check_ir false "Tried to evaluate $(arg) as an HLO argument"
        end
    end
    for (stmt_idx, stmt) in pairs(ir.stmts)
        isa(stmt, Nothing) && continue
        isa(stmt, ReturnNode) && continue
        @check_ir !isa(stmt, PhiNode) "Unsupported control flow found in function"
        if isa(stmt, PiNode)
            if !(widenconst(argextype(stmt.val, ir, sparams)) <: XLAScalar)
                ssa_vals[stmt_idx] = hlo_eval(stmt.val)
            end
            continue
        end
        if isexpr(stmt, :new) || (isexpr(stmt, :call) && is_known_call(stmt, tuple, ir, sparams))
            args = Any[]
            for arg in stmt.args[2:end]
                ty = argextype(arg, ir, sparams)
                representable(arg) || continue
                push!(args, hlo_eval(arg))
            end
            ssa_vals[stmt_idx] = HloInstructionProto(comp, HloTuple(), args...)
            continue
        end
        if isexpr(stmt, :call) && is_known_call(stmt, getfield, ir, sparams)
            structT = widenconst(argextype(stmt.args[2], ir, sparams))
            elt = argextype(stmt.args[3], ir, sparams)
            isa(elt, Const) || error("non-constant structure indexing (1)")
            idx = Compiler.try_compute_fieldidx(structT, elt.val)
            idx !== nothing || error("non-constant structure indexing (2)")
            if !representable(fieldtype(structT, idx))
                continue
            end
            tuple_idx = count(i->representable(fieldtype(structT, i)), 1:idx-1) + 1
            proto = HloInstructionProto(comp, HloGetTupleElement(tuple_idx - 1),
                hlo_eval(stmt.args[2]))
            ssa_vals[stmt_idx] = proto
            continue
        end
        if isa(stmt, GlobalRef) && isa(argextype(stmt, ir, sparams), Const)
            continue
        end
        if !isexpr(stmt, :invoke) && !(isexpr(stmt, :call) && (isa(stmt.args[1], _HloIf) || isa(stmt.args[1], _HloWhile)))
            @Base.show stmt
            Base.display(ir)
            error("Unrecognized expr")
        end
        hlo_inst = isexpr(stmt, :invoke) ? stmt.args[2] : stmt.args[1]
        hlo_inst = argextype(hlo_inst, ir, sparams)
        hlo_inst = hlo_inst.val
        if isa(hlo_inst, HloOp)
            if isa(hlo_inst, HloMap) || isa(hlo_inst, HloReduceWindow) || isa(hlo_inst, HloReduce) || isa(hlo_inst, HloSort)
                args = map(hlo_eval, stmt.args[4:end])
                proto = HloInstructionProto(comp, hlo_inst, args...)
                sig = Tuple{typeof(hlo_inst).parameters[1], (XRTArray{eltype(argextype(stmt.args[i], ir, sparams)), (), 0} for i = 4:length(stmt.args))...}
                argvals = process_function_argument(nothing, sig, 1, ir, stmt, sparams)
                res = code_typed_xla(sig; argvals=argvals)
                if res === nothing
                    throw(NotOffloadableError(ir, sv, sprint(io->showerror(io, MethodError(
                        Compiler.isconstType(sig.parameters[1]) ? sig.parameters[1].parameters[1] : sig.parameters[1].instance, Tuple{sig.parameters[2:end]...}
                    )))))
                end
                ir′, sv′ = res
                if !Compiler.isconstType(sig.parameters[1]) && sizeof(sig.parameters[1]) != 0
                    ir′.argtypes[1] = Const(argvals[1])
                end
                comp′ = HloComputationProto(
                    name = "comp$(length(computations))",
                    instructions = HloInstructionProto[ ],
                    id = length(computations)
                )
                pushfirst!(computations, comp′)
                _compile_to_xla!(computations, comp′, ir′, sv′)
                proto.called_computation_ids = [comp′.id]
            elseif isa(hlo_inst, Union{HloScatter, HloCrossReplicaSum})
                args = map(hlo_eval, stmt.args[4:end])
                shape = dtype_to_shape(infer_rt(hlo_inst, map(arg->widenconst(argextype(arg, ir, sparams)), stmt.args[3:end])...))
                proto = HloInstructionProto(comp, hlo_inst, args...; shape=shape)
                sig = Tuple{typeof(hlo_inst).parameters[1], (XRTArray{eltype(argextype(stmt.args[4], ir, sparams)), (), 0} for i = 1:2)...}
                argvals = process_function_argument(nothing, sig, 1, ir, stmt, sparams)
                ir′, sv′ = code_typed_xla(sig; argvals=argvals)
                if sizeof(sig.parameters[1]) != 0
                    ir′.argtypes[1] = Const(argvals[1])
                end
                comp′ = HloComputationProto(
                    name = "comp$(length(computations))",
                    instructions = HloInstructionProto[ ],
                    id = length(computations)
                )
                pushfirst!(computations, comp′)
                _compile_to_xla!(computations, comp′, ir′, sv′)
                proto.called_computation_ids = [comp′.id]
            elseif isa(hlo_inst, HloSelectAndScatter)
                args = map(hlo_eval, stmt.args[5:end])
                proto = HloInstructionProto(comp, hlo_inst, args...)
                proto.called_computation_ids = Int64[]
                # Compile comparison function
                let
                    select_sig = Tuple{typeof(hlo_inst).parameters[1], (XRTArray{eltype(argextype(stmt.args[5], ir, sparams)), (), 0} for i = 1:2)...}
                    select_argvals = process_function_argument(nothing, select_sig, 1, ir, stmt, sparams)
                    ir′, sv′ = code_typed_xla(select_sig; argvals=select_argvals)
                    if sizeof(select_sig.parameters[1]) != 0
                        ir′.argtypes[1] = Const(select_argvals[1])
                    end
                    comp′ = HloComputationProto(
                        name = "comp$(length(computations))",
                        instructions = HloInstructionProto[ ],
                        id = length(computations)
                    )
                    pushfirst!(computations, comp′)
                    _compile_to_xla!(computations, comp′, ir′, sv′)
                    push!(proto.called_computation_ids, comp′.id)
                end
                # Compile Scatter function
                let
                    scatter_sig = Tuple{typeof(hlo_inst).parameters[2], (XRTArray{eltype(argextype(stmt.args[5], ir, sparams)), (), 0} for i = 1:2)...}
                    scatter_argvals = process_function_argument(nothing, scatter_sig, 2, ir, stmt, sparams)
                    ir′, sv′ = code_typed_xla(scatter_sig; argvals=scatter_argvals)
                    if sizeof(scatter_sig.parameters[1]) != 0
                        ir′.argtypes[1] = Const(scatter_argvals[1])
                    end
                    comp′ = HloComputationProto(
                        name = "comp$(length(computations))",
                        instructions = HloInstructionProto[ ],
                        id = length(computations)
                    )
                    pushfirst!(computations, comp′)
                    _compile_to_xla!(computations, comp′, ir′, sv′)
                    push!(proto.called_computation_ids, comp′.id)
                end
            elseif isa(hlo_inst, _HloIf)
                @assert isexpr(stmt, :call)
                args = map(hlo_eval, stmt.args[2:end])
                if_type = ir.types[stmt_idx]
                proto = HloInstructionProto(comp,
                    GenericHloOp{:conditional}(if_type.parameters[1], if_type.parameters[2]), args...)
                proto.called_computation_ids = Int64[]
                for ir′ in (hlo_inst.if_func, hlo_inst.else_func)
                    comp′ = HloComputationProto(
                        name = "comp$(length(computations))",
                        instructions = HloInstructionProto[ ],
                        id = length(computations)
                    )
                    # Having a paremeter instruction is required, so if we don't
                    # otherwise have it, add it manually
                    if !representable(widenconst(ir′.argtypes[2]))
                        instr = HloInstructionProto(comp′, "parameter")
                        instr.shape = dtype_to_shape(Tuple{})
                        instr.parameter_number = 0
                    end
                    pushfirst!(computations, comp′)
                    _compile_to_xla!(computations, comp′, ir′, nothing)
                    push!(proto.called_computation_ids, comp′.id)
                end
            elseif isa(hlo_inst, _HloWhile)
                initial = hlo_eval(stmt.args[2])
                while_type = ir.types[stmt_idx]
                tuple_T = hlo_inst.condition_func.argtypes[2]
                shape = dtype_to_shape(tuple_T)
                proto = HloInstructionProto(comp,
                    GenericHloOp{:while}(Float32, ()), initial,
                    shape = shape)
                proto.called_computation_ids = Int64[]
                for ir′ in (hlo_inst.body_func, hlo_inst.condition_func)
                    comp′ = HloComputationProto(
                        name = "comp$(length(computations))",
                        instructions = HloInstructionProto[ ],
                        id = length(computations)
                    )
                    if !representable(widenconst(ir′.argtypes[2]))
                        instr = HloInstructionProto(comp′, "parameter")
                        instr.shape = dtype_to_shape(Tuple{})
                        instr.parameter_number = 0
                    end
                    pushfirst!(computations, comp′)
                    _compile_to_xla!(computations, comp′, ir′, nothing)
                    push!(proto.called_computation_ids, comp′.id)
                end
            elseif isa(hlo_inst, Union{HloInfeed, HloAfterAll})
                args = map(hlo_eval, stmt.args[3:end])
                shape = dtype_to_shape(infer_rt(hlo_inst, map(arg->argextype(arg, ir, sparams), stmt.args[3:end])...); tensorflow_order=isa(hlo_inst, HloInfeed))
                proto = HloInstructionProto(comp,
                    hlo_inst, args...,
                    shape = shape)
            elseif isa(hlo_inst, HloOutfeed)
                args = map(hlo_eval, stmt.args[3:end])
                shape = dtype_to_shape(infer_rt(hlo_inst, map(arg->argextype(arg, ir, sparams), stmt.args[3:end])...))
                proto = HloInstructionProto(comp,
                    hlo_inst, args...,
                    shape = shape)
                proto.outfeed_shape = dtype_to_shape(widenconst(argextype(stmt.args[3], ir, sparams)); tensorflow_order=true)
            else
                args = map(hlo_eval, stmt.args[3:end])
                shape = Shape(shape_infer(hlo_inst, map(arg->argextype(arg, ir, sparams), stmt.args[3:end])...)...)
                proto = HloInstructionProto(comp, hlo_inst, args...; shape=shape)
            end
            ssa_vals[stmt_idx] = proto
        elseif isa(hlo_inst, Type) && hlo_inst <: XRTArray
            ty = argextype(stmt.args[3], ir, sparams)
            if isa(ty, Const) && isa(ty.val, XLAScalar)
                ssa_vals[stmt_idx] = HloInstructionProto(comp, HloConstant(ty.val))
            else
                # Allow treating an XRTArray as a scalar inside another XRTArray
                # This is a no-op at the HLO level
                ty = widenconst(ty)
                @assert (hlo_inst <: XRTArray{<:Any, (), 0})
                ssa_vals[stmt_idx] = hlo_eval(stmt.args[3])
            end
        elseif hlo_inst == Base.convert
            ty = argextype(stmt.args[3], ir, sparams)
            arg = argextype(stmt.args[4], ir, sparams)
            if arg <: XRTArray{<:Any, (), 0} && isa(ty, Const) && ty.val == eltype(arg)
                # Allow conversions from a scalar array to its element
                # type and codegen it as a no-op.
                ssa_vals[stmt_idx] = hlo_eval(stmt.args[4])
            else
                @show stmt
                error("Unrecognized convert")
            end
        else
            @show stmt
            error("Unrecognized")
        end
    end
    ret = ir.stmts[end]
    @assert isa(ret, ReturnNode)
    comp.root_id = hlo_eval(ret.val).id
    xla_args
end

const ComputationDevice = DeviceAssignment_ComputationDevice
const DeviceMeshCoordinates = DeviceAssignment_ComputationDevice_DeviceMeshCoordinates
function compile_to_xla(ir, sv; replica_device_coords = nothing)
    # TODO: Since all HLO operations are essentially effect free, we could just have
    # a compilation result that throws this error.
    @check_ir isa(ir.stmts[end], ReturnNode) && isdefined(ir.stmts[end], :val) analyze_unreachable
    rt = widenconst(argextype(ir.stmts[end].val, ir, sv.sp))
    computations = HloComputationProto[]
    comp = HloComputationProto(
        name = "comp",
        instructions = HloInstructionProto[ ],
        id = 0
    )
    push!(computations, comp)
    xla_args = _compile_to_xla!(computations, comp, ir, sv)

    pshape = ProgramShape(
        parameters = collect(map(dtype_to_shape, xla_args)),
        result = dtype_to_shape(rt)
    )
    if replica_device_coords !== nothing
        coordinates = vec([ DeviceMeshCoordinates(value=[dev.x, dev.y, dev.core]) for dev in replica_device_coords ])
        cd = ComputationDevice(
            replica_devices = coordinates
        )
        da = DeviceAssignment(
            computation_devices = [cd]
        )
        config = XLAComputationConfig(
            program_shape = pshape,
            device_assignment = da,
            num_replicas = length(coordinates)
        )
    else
        config = XLAComputationConfig(
            program_shape = pshape
        )
    end
    hlo_module = HloModuleProto(
        name = "test",
        computations = computations,
        entry_computation_name = "comp",
        entry_computation_id = 0,
        id = 0,
        host_program_shape = pshape,
    )
    hlo = HloProto(
        hlo_module = hlo_module
    )
    hlo_snap = HloSnapshot(
        hlo = hlo
    )
    xlac = XLAComputation(
        config=config,
        hlo_snapshot = hlo_snap
    )
    xlac, rt
end

function representable(T)
    return sizeof(T) != 0
end

function dtype_to_shape(T::Type{<:XRTArray}; tensorflow_order=false)
    shp = convert(Shape, T)
    if tensorflow_order
        shp.layout.minor_to_major = collect((length(shp.layout.minor_to_major) - 1):-1:0)
    end
    shp
end
dtype_to_shape(T::Type{<:XLA.XLAScalar}; tensorflow_order=false) = dtype_to_shape(XRTArray{T, (), 0})
function dtype_to_shape(T::Type{HloToken}; tensorflow_order=false)
    Shape(
        element_type = xla.PrimitiveType.TOKEN
    )
end
function dtype_to_shape(T::DataType; tensorflow_order=false)
    Shape(
        element_type = xla.PrimitiveType.TUPLE,
        tuple_shapes = collect(dtype_to_shape(fieldtype(T, i); tensorflow_order=tensorflow_order) for i = 1:fieldcount(T) if representable(fieldtype(T, i)))
    )
end


text:
module TestTableTraits
using Test, TableTraits, DataFrames

struct ColumnSource
end

TableTraits.supports_get_columns_copy_using_missing(::ColumnSource) = true

function TableTraits.get_columns_copy_using_missing(x::ColumnSource)
    return (a=[1,2,3], b=[4.,5.,6.], c=["A", "B", "C"])
end

@testset "TableTraits" begin

    df = DataFrame(ColumnSource())

    @test size(df)==(3,3)
    @test df.a==[1,2,3]
    @test df.b==[4.,5.,6.]
    @test df.c==["A", "B", "C"]

end

end


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
using LaTeXStrings
theme(:wong)
gr(aspect_ratio=1,legend=nothing,axis=nothing,border=nothing,ticks=nothing,
   markerstrokewidth=0,markersize=1)
r1D, w1D = gauss_lobatto_quad(0,0,4)

# plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
# plot!([-1;1],[1;1],seriescolor=:black,lw=3,label = nothing)
# plot!([1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
# plot!([1;-1],[-1;-1],seriescolor=:black,lw=3,label = nothing)
plot([-1;1],[1;1],seriescolor=:black,lw=3,label = nothing)
plot!([r1D[2];0],[1;1],seriescolor=:grey,linestyle=:dash,lw=3,label = nothing)
scatter!(r1D,ones(size(r1D)),markersize=7,markercolor=:royalblue1,markerstrokecolor=:black,markerstrokewidth=2)
scatter!(r1D[2:4],ones(3),markersize=7,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2)
#savefig("LGL.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using MultivariateStats
using LinearAlgebra
using Test
using StableRNGs
import Statistics: mean, var, cov
import SparseArrays

@testset "PCA" begin

    rng = StableRNG(34568)

    ## PCA with zero mean

    X = randn(rng, 5, 10)
    Y = randn(rng, 3, 10)

    P = qr(randn(rng, 5, 5)).Q[:, 1:3]
    pvars = [5., 4., 3.]
    l = [-0.236399  -0.706558  -1.21106
          1.27426   -0.047799  -0.212282
          0.990598   1.22718   -0.882179
         -1.51861    0.725771  -0.631665
          0.181386  -1.21069   -0.557706]
    M = PCA(Float64[], P, pvars, 15.0)

    @test size(M) == (5, 3)
    @test mean(M) == zeros(5)
    @test projection(M) == P
    @test principalvars(M) == pvars
    @test principalvar(M, 2) == pvars[2]
    @test var(M) == 15.0
    @test tprincipalvar(M) == 12.0
    @test tresidualvar(M) == 3.0
    @test principalratio(M) == 0.8
    @test isapprox(loadings(M), l, atol = 0.001)

    @test predict(M, X[:,1]) ≈ P'X[:,1]
    @test predict(M, X) ≈ P'X

    @test reconstruct(M, Y[:,1]) ≈ P * Y[:,1]
    @test reconstruct(M, Y) ≈ P * Y

    ## PCA with non-zero mean

    mval = rand(rng, 5)
    M = PCA(mval, P, pvars, 15.0)

    @test size(M) == (5,3)
    @test mean(M) == mval
    @test projection(M) == P
    @test principalvars(M) == pvars
    @test principalvar(M, 2) == pvars[2]
    @test var(M) == 15.0
    @test tprincipalvar(M) == 12.0
    @test tresidualvar(M) == 3.0
    @test principalratio(M) == 0.8

    @test predict(M, X[:,1]) ≈ P' * (X[:,1] .- mval)
    @test predict(M, X) ≈ P' * (X .- mval)

    @test reconstruct(M, Y[:,1]) ≈ P * Y[:,1] .+ mval
    @test reconstruct(M, Y) ≈ P * Y .+ mval

    ## prepare training data

    d = 5
    n = 1000

    R = collect(qr(randn(rng, d, d)).Q)
    @test R'R ≈ Matrix(I, 5, 5)
    rmul!(R, Diagonal(sqrt.([0.5, 0.3, 0.1, 0.05, 0.05])))

    X = R'randn(rng, 5, n) .+ randn(rng, 5)
    mval = vec(mean(X, dims=2))
    Z = X .- mval

    C = cov(X, dims=2)
    pvs0 = sort(eigvals(C); rev=true)
    tv = sum(diag(C))

    ## pcacov (default using cov when d < n)

    M = fit(PCA, X)
    P = projection(M)
    pvs = principalvars(M)

    @test size(M) == (5,5)
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 5, 5)
    @test C*P ≈ P*Diagonal(pvs)
    @test issorted(pvs; rev=true)
    @test pvs ≈ pvs0
    @test var(M) ≈ tv
    @test sum(pvs) ≈ var(M)
    @test reconstruct(M, predict(M, X)) ≈ X

    M = fit(PCA, X; mean=mval)
    @test projection(M) ≈ P

    M = fit(PCA, Z; mean=0)
    @test projection(M) ≈ P

    M = fit(PCA, X; maxoutdim=3)
    P = projection(M)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)
    @test issorted(pvs; rev=true)

    M = fit(PCA, X; pratio=0.85)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)
    @test issorted(pvs; rev=true)

    ## pcastd

    M = fit(PCA, X; method=:svd)
    P = projection(M)
    pvs = principalvars(M)

    @test size(M) == (5,5)
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 5, 5)
    @test isapprox(C*P, P*Diagonal(pvs), atol=1.0e-3)
    @test issorted(pvs; rev=true)
    @test isapprox(pvs, pvs0, atol=1.0e-3)
    @test isapprox(var(M), tv, atol=1.0e-2)
    @test sum(pvs) ≈ var(M)
    @test reconstruct(M, predict(M, X)) ≈ X

    M = fit(PCA, X; method=:svd, mean=mval)
    @test projection(M) ≈ P

    M = fit(PCA, Z; method=:svd, mean=0)
    @test projection(M) ≈ P

    M = fit(PCA, X; method=:svd, maxoutdim=3)
    P = projection(M)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)
    @test issorted(pvs; rev=true)

    M = fit(PCA, X; method=:svd, pratio=0.85)

    @test size(M) == (5,3)
    @test P'P ≈ Matrix(I, 3, 3)
    @test issorted(pvs; rev=true)

    # Different data types
    # --------------------

    XX = convert.(Float32, X)
    YY = convert.(Float32, Y)
    p = 0.085
    pp = convert(Float32, p)

    MM = fit(PCA, XX; maxoutdim=3)

    # mix types
    fit(PCA, X ; pratio=pp)
    fit(PCA, XX ; pratio=p)
    fit(PCA, XX ; pratio=pp)
    predict(M, XX)
    predict(MM, X)
    reconstruct(M, YY)
    reconstruct(MM, Y)

    # type consistency
    for func in (mean, projection, principalvars, tprincipalvar, tresidualvar, var, principalratio)
        @test eltype(func(M)) == Float64
        @test eltype(func(MM)) == Float32
    end

    # views
    M = fit(PCA, view(X, :, 1:500), pratio=0.85)

    # sparse
    @test_throws AssertionError fit(PCA, SparseArrays.sprandn(rng, 100d, n, 0.6))

end


text:
function ob(r,Side)
    x = 0
    Flg = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide((Side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        Flg = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide(Side))
    end
    for i = 1:x
        move!(r,HorizonSide((Side+1) % 4))
    end
    return Flg
end

function ob2(r,Side)
    a = 0
    b = 0
    x = 0
    Flg = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide((Side + 1) % 4))
        x += 1
        a += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        Flg = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide(Side))
        b += 1
    end
    for i = 1:x
        move!(r,HorizonSide((Side+3) % 4))
    end
    return Flg, (a)*(b)
end

function corner(r)
    Flg = true
    y = true
    while !isborder(r,Sud) || Flg
        if !isborder(r,Sud)
            move!(r,Sud)
        end
        if isborder(r,Sud)
            Flg =  ob(r,2)
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
end

function scan(r)
    Flg = true
    sum = 0
    Z = 0
    if !isborder(r,Nord)
        Flg = true
        move!(r,Nord)
        while !isborder(r,Ost) || Flg
            if !isborder(r,Ost)
                move!(r,Ost)
            else
                Flg, Z = ob2(r,3)
                if sum < Z
                    sum = Z
                end
            end
        end
    end
    Flg = true
    while !isborder(r, West) || Flg
        if !isborder(r, West)
            move!(r, West)
        else
            Flg = ob(r,1)
        end
    end
    return sum
end

function task_24(r)
    sum = 0
    Z = 0
    corner(r)
    while !isborder(r,Nord)
        Z = scan(r)
        if sum < Z
            sum = Z
        end
    end
    return sum
end

task_24(r)


text:
function calculate_initial_conditions!(sys::PSY.System, initial_guess::Vector{Float64})

    res = PSY.solve_powerflow!(sys)
    if !res
        error("Power Flow fail")
    end
    var_count = get_variable_count(sys)

    #Index Setup
    bus_size = get_bus_count(sys)
    bus_vars_count = 2 * bus_size
    bus_range = 1:bus_vars_count
    injection_start = get_injection_pointer(sys)
    injection_count = 1
    branches_start = get_branches_pointer(sys)
    branches_count = 1

    #Update Voltage guess
    for bus in PSY.get_components(PSY.Bus, sys)
        #Write voltage initial guess
        bus_n = PSY.get_number(bus)
        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
        initial_guess[bus_ix] = PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        initial_guess[bus_ix + bus_size] = PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
    end

    #Update Source internal voltages
    sources = PSY.get_components(PSY.Source, sys)
    if !isempty(sources)
        for s in sources
            initialize_device(s)
        end
    end

    #Update Dynamic Injection internal references and guesses
    for d in PSY.get_components(PSY.DynamicInjection, sys)
        bus = PSY.get_bus(d)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
        n_states = PSY.get_n_states(d)
        ix_range = range(injection_start, length = n_states)
        injection_start = injection_start + n_states
        x0_device = initialize_device(d)
        @assert length(x0_device) == n_states
        initial_guess[ix_range] = x0_device
    end

    #Update Dynamic Branch guess
    dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)
    if !isempty(dyn_branches)
        for br in dyn_branches
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = PSY.get_ext(sys)[LOOKUP][from_bus_number]
            bus_ix_to = PSY.get_ext(sys)[LOOKUP][to_bus_number]
            ix_range = range(branches_start, length = n_states)
            branches_start = branches_start + n_states
            x0_branch = initialize_device(br)
            @assert length(x0_branch) == n_states
            initial_guess[ix_range] = x0_branch
        end
    end

    dx0 = zeros(var_count) #Define a vector of zeros for the derivative
    inif! = (out, x) -> system!(
        out,    #output of the function
        dx0,    #derivatives equal to zero
        x,      #states
        sys,    #Parameters
        0.0,    #time equals to zero.
    )
    #Refine initial solution
    sys_solve = NLsolve.nlsolve(
        inif!,
        initial_guess,
        xtol = :1e-9,
        ftol = :1e-9,
        method = :trust_region,
    ) #Solve using initial guess x0
    if !NLsolve.converged(sys_solve)
        @warn("Initialization failed, initial conditions do not meet conditions for an stable equilibrium.\nTrying to solve again reducing numeric tolerance:")
        sys_solve = NLsolve.nlsolve(
            inif!,
            initial_guess,
            xtol = :1e-6,
            ftol = :1e-6,
            method = :trust_region,
        ) #Solve using initial guess x0
        if NLsolve.converged(sys_solve)
            @info("Initialization succeeded with less tolerance. Saving solution")
        else
            @warn("Initialization failed again. Initial conditions do not meet conditions for an stable equilibrium\nSaving best result.")
        end
    end
    initial_guess .= sys_solve.zero
    return NLsolve.converged(sys_solve)

end

"""
Returns a Dictionary with the resulting initial conditions of the simulation
"""
function get_initial_conditions(sim::Simulation)
    bus_size = get_bus_count(sim.system)
    V_R = Dict{Int64, Float64}()
    V_I = Dict{Int64, Float64}()
    Vm = Dict{Int64, Float64}()
    θ = Dict{Int64, Float64}()
    for bus in PSY.get_components(PSY.Bus, sim.system)
        bus_n = PSY.get_number(bus)
        bus_ix = PSY.get_ext(sim.system)[LOOKUP][bus_n]
        V_R[bus_n] = sim.x0_init[bus_ix]
        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
        θ[bus_n] = angle(sim.x0_init[bus_ix] + sim.x0_init[bus_ix + bus_size] * 1im)
    end
    results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in PSY.get_components(PSY.DynamicInjection, sim.system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = PSY.get_ext(sim.system)[GLOBAL_INDEX][name]
        x0_device = Dict{Symbol, Float64}()
        for (i, s) in enumerate(states)
            x0_device[s] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, sim.system)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = PSY.get_ext(sim.system)[GLOBAL_INDEX][name]
            x0_br = Dict{Symbol, Float64}()
            for (i, s) in enumerate(states)
                x0_br[s] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end


text:
# test cross-validation, when the simulation includes home-ground advantage
using RatPack
using Random
using Distributions
using PlotlyJS

# generate large simulated dataset
include("test_sim_3.jl")

# test prediction/scoring
Random.seed!(2)
n2 = 1000
sim2 = SimRoundRobinF(;n=n2, factor_scale=real_f)
(test_competitions, winner) = simulate( real_ratings, sim2, perf_model )
srule = ScoreBrier(;normalise=true)
s0 = score_ratings( srule, test_competitions, irule, initial_ratings )
s1 = score_ratings( srule, test_competitions, irule, simple_ratings )
s2 = score_ratings( srule, test_competitions, irule, real_ratings )

println("score with real ratings = $s2, score with simulated ratings = $s1")

# test cross-validation
k = 10
Random.seed!(1)
S1 = cross_validate(irule, input_competitions;
                    k = k,
                    n_samples = 1000,
                    batch_size = 1,
                    srule = ScoreBrier(;normalise=true) )

println(" S1bar = $(sum(S1)/k)")

# see what happens when prediction rule mismatch with rule used to create the ratings
s3 = score_ratings( srule, test_competitions, UpdateColley(), real_ratings )

f_s = 0.0 : 20.0 : 600.0
score_f = zeros(size(f_s))
for (i,f) in enumerate( f_s )
    tmp_rule = UpdateIterate( rule=UpdateEloF(;r0 = r0, K = 30.0, θ=θ, factor_scale=f), batch_size=batch_size )
    new_ratings = update_ratings(tmp_rule, initial_ratings, input_competitions)
    score_f[i] = score_ratings( srule, test_competitions, tmp_rule, new_ratings )
end

data1 = scatter(;x=f_s, y=score_f, mode="lines+markers", name="misfit scale factor")
# data2 = scatter(;x=[θ, θ], y=[s2, s2], mode="markers", name="optimal")
# data3 = scatter(;x=[0.0, 0.0], y=[s3, s3], mode="markers", name="deterministic")
# data = [data1, data2, data3]
data = [data1]
layout = Layout(;
                xaxis=attr(title="f_s" ),
                yaxis=attr(title="(normalised) Brier score" ),
                )
PlotlyJS.plot( data, layout )


text:
function MCTS.node_tag(s::MLState)
    if s.crashed
        return "CRASH"
    else
        desc = "("
        for c in s.cars
            desc = string(desc, @sprintf("[%.1f,%.1f]",c.x,c.y))
        end
        return string(desc,")")
    end
end

function MCTS.tooltip_tag(s::MLState)
    if s.crashed
        return "CRASH"
    else
        desc = "( "
        for c in s.cars
            desc = string(desc,
            @sprintf("[%.1f,%.1f,v:%.1f,l:%.1f] ", c.x, c.y, c.vel, c.lane_change))
        end
        return string(desc,")")
    end
end

function MCTS.node_tag(a::MLAction)
    return @sprintf("[%.1f,%.1f]", a.acc, a.lane_change)
end

function MCTS.tooltip_tag(a::MLAction)
    return @sprintf("[a:%.1f,l:%.1f]", a.acc, a.lane_change)
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
######## Multiple fBm ########
"""
Multiple fBm defined on [0, Tmax]

# Members
- hurstfunc: the time varying Hurst function
- hurstvec: the time varying Hurst vector
- Tmax: the ending time
"""
struct MultipleFractionalBrownianMotion <: ContinuousTimeStochasticProcess
    hurstfunc::Function  # hurst(t) is the hurst exponent at time t
    hurstvec::AbstractVector{<:Real}
    Tmax::Real

    function FractionalBrownianMotion(hurstfunc::Function)
        new(hurstfunc, Float64[], +Inf)
    end

    function FractionalBrownianMotion(hurstvec::AbstractVector{<:Real}, Tmax::Real=1.)
        all(0. .< hurstvec .< 1.) || error("Hurst exponent must be bounded in 0 and 1.")
        @assert Tmax > 0.

        N = length(hurstvec)
        hurstfunc = t -> (t<0 || t>Tmax) ? NaN : hurst[max(1, round(Int, (t/Tmax)*N))]
        new(hurstfunc, hurstvec, Tmax)
    end
end

# Auxiliary functions for the covariance of the  multifractional field
_D(h1, h2) = (h1==h2) ? 1. : sqrt(gamma(2*h1+1)*sin(pi*h1)*gamma(2*h2+1)*sin(pi*h2)) / (gamma(h1+h2+1)*sin(pi*(h1+h2)/2))
# or equivalently
# _F(h) = lgamma(2h+1) + log(sin(π*h))
# _D(h1,h2) = exp((_F(h1)+_F(h2))/2 - _F((h1+h2)/2))

# h actually corresponds to 2h here:
_gn(t, s, h) = (abs(t)^(h) + abs(s)^(h) - abs(t-s)^(h))/2  # non-stationary
_gs(t, h) = (abs(t+1)^(h) + abs(t-1)^(h) - 2*abs(t)^(h))/2  # stationary

_mfBm_autocov(t1, t2, h1, h2) = _D(h1, h2) * _gn(t1, t2, h1+h2)

function autocov(P::MultipleFractionalBrownianMotion, t::Real, s::Real)
    return _mfBm_autocov(t, s, P.hurst(t), P.hurst(s))
end

function autocov(P::MultipleFractionalBrownianMotion, t::Int, s::Int)
    N = length(P.hurstvec)
    return _mfBm_autocov(t/N * P.Tmax, s/N * P.Tmax, P.hurstvec[t], P.hurstvec[s])
end

struct MultipleFractionalGaussianNoise <: IncrementProcess{MultipleFractionalBrownianMotion}
    parent_process::MultipleFractionalBrownianMotion
    step::Real

    hurstfunc::Function  # hurst(t) is the hurst exponent at time t
    hurstvec::AbstractVector{<:Real}
    Tmax::Real

    function MultipleFractionalGaussianNoise(hurstfunc::Function, step::Real=1.)
        0. < step || error("Invalid value for step.")
        new(MultipleFractionalBrownianMotion(hurstfunc), step)
    end

    function MultipleFractionalGaussianNoise(hurstvec::Real, Tmax::Real, step::Real=1.)
        0. < step < Tmax || error("Invalid value for step.")
        new(MultipleFractionalBrownianMotion(hurstvec, Tmax), step)
    end
end

"""
(Approximate) Covariance of multifractional Gaussian noise.
"""
_mfGn_autocov(t1, t2, h1, h2) = _D(h1, h2) * _gs(t1-t2, h1+h2)

function autocov(P::MultipleFractionalGaussianNoise, t::Real, s::Real)
    return _mfGn_autocov(t, s, P.hurst(t), P.hurst(s))
end

function autocov(P::MultipleFractionalGaussianNoise, t::Int, s::Int)
    N = length(P.hurstvec)
    return _mfGn_autocov(t/N * P.Tmax, s/N * P.Tmax, P.hurstvec[t], P.hurstvec[s])
end

function mGn_covmat(G1::AbstractVector{<:Real}, H1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, H2::AbstractVector{<:Real})
    @assert length(G1) == length(H1)
    @assert all(0 .< H1 .< 1)
    @assert length(G2) == length(H2)
    @assert all(0 .< H2 .< 1)

    N1, N2 = length(G1), length(G2)
    Σ = zeros(N1,N2)

    for c=1:N2, r=1:N1
        Σ[r,c] = mGn_cov(G1[r], G2[c], H1[r], H2[c])
    end
    return Σ
end

mGn_covmat(G::AbstractVector{<:Real}, H::AbstractVector{<:Real}) = mGn_covmat(G, H, G, H)

mGn_covmat(G::AbstractVector{<:Real}, H::Real) = mGn_covmat(G, H)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
# Autogenerated wrapper script for SLEEF_jll for armv6l-linux-gnueabihf
export libsleef, libsleefquad

using CompilerSupportLibraries_jll
using MPFR_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("SLEEF")
JLLWrappers.@declare_library_product(libsleef, "libsleef.so.3")
JLLWrappers.@declare_library_product(libsleefquad, "libsleefquad.so.3")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, MPFR_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libsleef,
        "lib/libsleef.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsleefquad,
        "lib/libsleefquad.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# CUDA-specific operations on pointers with address spaces

## adrspace aliases

export AS

module AS

const Generic  = 0
const Global   = 1
const Shared   = 3
const Constant = 4
const Local    = 5

end

## ldg

const LDGTypes = (UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64,
                  Float32, Float64)

# TODO: this functionality should throw <sm_32
# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet
for T in LDGTypes
    class = if T <: Integer
        :i
    elseif T <: AbstractFloat
        :f
    end
    # TODO: p class
    width = sizeof(T)*8 # in bits
    typ = Symbol(class, width)

    intr = "llvm.nvvm.ldg.global.$class.$typ.p1$typ"
    @eval @inline function pointerref_ldg(base_ptr::LLVMPtr{$T,AS.Global}, i::Integer,
                                          ::Val{align}) where align
        offset = i-one(i) # in elements
        ptr = base_ptr + offset*sizeof($T)
        @typed_ccall($intr, llvmcall, $T, (LLVMPtr{$T,AS.Global}, Int32), ptr, align)
    end
end

# interface

export unsafe_cached_load

unsafe_cached_load(p::LLVMPtr{<:Union{LDGTypes...},AS.Global}, i::Integer=1, align::Val=Val(1)) =
    pointerref_ldg(p, i, align)
# NOTE: fall back to normal unsafe_load for unsupported types. we could be smarter here,
#       e.g. destruct/load/reconstruct, but that's too complicated for what it's worth.
unsafe_cached_load(p::LLVMPtr, i::Integer=1, align::Val=Val(1)) =
    unsafe_load(p, i, align)


text:
function nlp_cvx_105_011(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    model = Model(optimizer)

    @variable(model, x, start = 0.1)
    @variable(model, y)

    @objective(model, Min, x + y)
    @NLconstraint(model, exp(x - 2.0) - 0.5 <= y)
    @NLconstraint(model, log(x) + 0.5 >= y)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 0.16878271368156372, tol = objective_tol)
    return check_solution(
        [x, y],
        [0.45538805755556067, -0.28660534387399694],
        tol = primal_tol,
    )
end


text:
module Sesi

using JuMP
using Gurobi
using CSV
using DataFrames
using JSON
using IterTools

include("design.jl")
include("forecast.jl")
include("model.jl")
include("utils.jl")

end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
function stacked_spectrum!(dobs; fs=1.0)
	D=rfft(dobs,1) # perform rfft
	window=zeros(size(D,1))
	nr=size(dobs,2)
	# stack the spectrum of data
	for j in 1:nr
		for i in eachindex(window)
			window[i] += (abs2(D[i,j]))
		end
	end
	rmul!(window,inv(maximum(abs,window)))
	window = 10. * log10.(window)
	return Float64.(DSP.Util.rfftfreq((size(dobs,1)),fs)), window
end


text:
const _all_defaults = KW[_series_defaults, _plot_defaults, _subplot_defaults]

const _initial_defaults = deepcopy(_all_defaults)
const _initial_axis_defaults = deepcopy(_axis_defaults)

# add defaults for the letter versions
const _axis_defaults_byletter = KW()

function reset_axis_defaults_byletter!()
    for letter in (:x, :y, :z)
        _axis_defaults_byletter[letter] = KW()
        for (k, v) in _axis_defaults
            _axis_defaults_byletter[letter][k] = v
        end
    end
end
reset_axis_defaults_byletter!()

# to be able to reset font sizes to initial values
const _initial_plt_fontsizes =
    Dict(:plot_titlefontsize => _plot_defaults[:plot_titlefontsize])

const _initial_sp_fontsizes = Dict(
    :titlefontsize => _subplot_defaults[:titlefontsize],
    :legend_font_pointsize => _subplot_defaults[:legend_font_pointsize],
    :legend_title_font_pointsize => _subplot_defaults[:legend_title_font_pointsize],
    :annotationfontsize => _subplot_defaults[:annotationfontsize],
    :colorbar_tickfontsize => _subplot_defaults[:colorbar_tickfontsize],
    :colorbar_titlefontsize => _subplot_defaults[:colorbar_titlefontsize],
)

const _initial_ax_fontsizes = Dict(
    :tickfontsize  => _axis_defaults[:tickfontsize],
    :guidefontsize => _axis_defaults[:guidefontsize],
)

const _initial_fontsizes =
    merge(_initial_plt_fontsizes, _initial_sp_fontsizes, _initial_ax_fontsizes)

const _internal_args =
    [:plot_object, :series_plotindex, :markershape_to_add, :letter, :idxfilter]

const _axis_args = sort(union(collect(keys(_axis_defaults))))
const _series_args = sort(union(collect(keys(_series_defaults))))
const _subplot_args = sort(union(collect(keys(_subplot_defaults))))
const _plot_args = sort(union(collect(keys(_plot_defaults))))

const _magic_axis_args = [:axis, :tickfont, :guidefont, :grid, :minorgrid]
const _magic_subplot_args =
    [:title_font, :legend_font, :legend_title_font, :plot_title_font, :colorbar_titlefont]
const _magic_series_args = [:line, :marker, :fill]

const _all_axis_args = sort(union([_axis_args; _magic_axis_args]))
const _all_subplot_args = sort(union([_subplot_args; _magic_subplot_args]))
const _all_series_args = sort(union([_series_args; _magic_series_args]))
const _all_plot_args = _plot_args

const _all_args = sort(union([
    _all_axis_args
    _all_subplot_args
    _all_series_args
    _all_plot_args
]))

# add all pluralized forms to the _keyAliases dict
for arg in _all_args
    add_aliases(arg, makeplural(arg))
end

# fill symbol cache
for letter in (:x, :y, :z)
    _attrsymbolcache[letter] = Dict{Symbol,Symbol}()
    for k in _axis_args
        # populate attribute cache
        lk = Symbol(letter, k)
        _attrsymbolcache[letter][k] = lk
        # allow the underscore version too: xguide or x_guide
        add_aliases(lk, Symbol(letter, "_", k))
    end
    for k in (_magic_axis_args..., :(_discrete_indices))
        _attrsymbolcache[letter][k] = Symbol(letter, k)
    end
end

# add all non_underscored forms to the _keyAliases
add_non_underscore_aliases!(_keyAliases)


text:
using HTTP
using EzXML
using Dates

# QiitaのJuliaタグの文書のタイトルとURLを取得
function get_qiita()
    qiita_feed = "https://qiita.com/tags/julia/feed"
    # ページの取得
    response = HTTP.request("GET", qiita_feed)
    if response.status != 200
        return []
    end

    code = String(response.body)
    doc = parsexml(code)
    # feed>entry にある url と title を抜き出す
    feed = root(doc)
    list = []
    for tag in eachelement(feed)
        if tag.name == "entry"
            url = nothing
            title = nothing
            for item in eachelement(tag)
                if item.name == "url"
                    url = item.content
                elseif item.name == "title"
                    title = item.content
                end
            end
            if !isnothing(url)
                pair = (title, url)
                push!(list, pair)
            end
        end
    end
    list
end

function main()
    # 日付付きのファイル名を作成
    filename = "qiita" * Dates.format(now(), "yyyymmdd") * ".txt"
    open(filename, "w") do f
        list = get_qiita()
        for (title, url) in list
            println(f, title * "\t" * url)
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: julialangjp/ScrapingQiitaTagByJulia

max_issues_repo_name: julialangjp/ScrapingQiitaTagByJulia

max_forks_repo_name: julialangjp/ScrapingQiitaTagByJulia