text:
module CheqFacial

function BubbleSort(v::AbstractVector, lo::Int, hi::Int)
    while true
        clean = true
        for i = lo:hi-1
            if ( v[i+1]<v[i])
                v[i+1], v[i] =v[i], v[i+1]
                clean = false
            end
        end
        clean && break
    end
    return v
end

#println("start")
v=rand(5)*10
println(v)
println(BubbleSort(v,1,5))

end

ext: jl


text:
using C2Julia
using Clang
using Test

@testset "unary operator" begin
    tu = parse_header(joinpath(@__DIR__, "c", "operator.h"), includes=[LLVM_INCLUDE])
    GC.@preserve tu begin
        root_cursor = getcursor(tu)
        child_cursors = children(root_cursor)
        op_func = search(root_cursor, "operators")[1]
        body = children(op_func)[1]
        op = children(body)[4]
        @test string(translate(op)) == "+5"
        op = children(body)[5]
        @test string(translate(op)) == "-5"
        op = children(body)[6]
        @test translate(op).expr == :(!z)
        op = children(body)[7]
        @test translate(op).expr == :(Cint(z))
        op = children(body)[8]
        @test string(translate(op)) == "@++ y"
        op = children(body)[9]
        @test translate(op).expr == :(y+=1)
        op = children(body)[10]
        @test string(translate(op)) == "@- x"
        op = children(body)[11]
        @test translate(op).expr == :(x-=1)
        op = children(body)[12]
        @test translate(op)[].expr == :(xp = Ref(x))
        op = children(body)[13]
        @test string(translate(op)) == "xp[] = 1"
        op = children(body)[14]
        @test translate(op).expr == :(Bool(z+1))
        op = children(body)[15]
        @test string(translate(op).expr) == "x == 0 && y < 1 || z == 3"
        op = children(body)[16]
        @test translate(op)[].expr == :(f = foo(x))
        op = children(body)[17]
        @test string(translate(op).expr) == "f.x += 1"
        op = children(body)[18]
        string(translate(op).expr) == string(:(y = z ? 3 : 30))
    end
end


text:
using REPL
using REPL.Terminals
using REPL.TerminalMenus
using ObjectFile.ELF
using HTTP
import PkgLicenses
using MbedTLS
using JLD2
using Pkg.BinaryPlatforms
using Dates

# It's Magic (TM)!
export run_wizard

include("wizard/state.jl")
include("wizard/github.jl")
include("wizard/yggdrasil.jl")
include("wizard/utils.jl")
include("wizard/obtain_source.jl")
include("wizard/interactive_build.jl")
include("wizard/deploy.jl")

function save_last_wizard_state(state::WizardState)
    create_and_bind_mutable_artifact!("wizard_state") do dir
        jldopen(joinpath(dir, "wizard.state"), "w") do f
            serialize(f, state)
        end
    end
    return state
end

function load_last_wizard_state()
    wizard_state_dir = get_mutable_artifact_path("wizard_state")

    # If no state dir exists, early-exit
    if wizard_state_dir === nothing
        return WizardState()
    end

    try
        state = jldopen(joinpath(wizard_state_dir, "wizard.state"), "r") do f
            return unserialize(f)
        end

        # Looks like we had an incomplete build; ask the user if they want to continue
        if !(state.step in (:done, :step1))
            terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
            choice = request(terminal,
                "Would you like to resume the previous incomplete wizard run?",
                RadioMenu([
                    "Resume previous run",
                    "Start from scratch",
                ]),
            )

            if choice == 1
                return state
            else
                return WizardState()
            end
        end
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @error(e)
    end

    # Either something went wrong, or there was nothing interesting stored.
    # Either way, just return a blank slate.
    return WizardState()
end

function run_wizard(state::Union{Nothing,WizardState} = nothing)
    global last_wizard_state

    if state === nothing
        # If we weren't given a state, check to see if we'd like to resume a
        # previous run or start from scratch again.
        state = load_last_wizard_state()
    end

    try
        while state.step != :done
            if state.step == :step1
                step1(state)
                state.step = :step2
            elseif state.step == :step2
                step2(state)
                state.step = :step3
            elseif state.step == :step3
                step34(state)
            elseif state.step == :step3_retry
                step3_retry(state)
            elseif state.step == :step5a
                step5a(state)
            elseif state.step == :step5b
                step5b(state)
            elseif state.step == :step5c
                step5c(state)
                state.step = :step6
            elseif state.step == :step6
                step6(state)
            elseif state.step == :step7
                step7(state)
                state.step = :done
            end

            # Save it every step along the way
            save_last_wizard_state(state)
        end
    catch err
        # If anything goes wrong, immediately save the current wizard state
        save_last_wizard_state(state)
        if isa(err, InterruptException)
            msg = "\n\nWizard stopped, use run_wizard() to resume.\n\n"
            printstyled(state.outs, msg, bold=true, color=:red)
        else
            bt = catch_backtrace()
            Base.showerror(stderr, err, bt)
            println(state.outs, "\n")
        end
        return state
    end

    # We did it!
    save_last_wizard_state(state)

    println(state.outs, "\nWizard Complete. Press any key to exit...")
    read(state.ins, Char)

    state
end


text:
@symbol_func function P_H(cur_reactor::AbstractReactor)
  cur_P = P_F(cur_reactor)

  cur_P /= cur_reactor.Q

  cur_P
end


text:
#=
Methods to simplify the plotting functionality for TS objects.
=#

const DEFAULT_COLORS = [:black, :red, :green, :blue, :cyan, :magenta, :orange, :pink]

import RecipesBase: @recipe

function beautify_string(s::String;
                         preserve_acronyms::Bool=true,
                         space_before_nums::Bool=true,
                         rm_underscores::Bool=true)
    idx_upper = findall(isuppercase, s)
    out_arr = String[string(s[1])]
    @inbounds for i in 2:length(s)
        if isuppercase(s[i]) && preserve_acronyms && !isuppercase(s[i-1])
            push!(out_arr, " $(s[i])")
        elseif isnumeric(s[i]) && space_before_nums && !isnumeric(s[i-1])
            push!(out_arr, " $(s[i])")
        else
            push!(out_arr, string(s[i]))
        end
    end
    return rm_underscores ? replace(join(out_arr), "_" => "") : join(out_arr)
end

function tslab(flds::Vector{Symbol}; beautify=true, args...)
    arr::Vector{String} = beautify ? beautify_string.(string.(flds); args...) : string.(flds)
    out::Matrix{String} = repeat([""], 1, length(flds))
    @inbounds for i in 1:length(flds)
        out[1,i] = arr[i]
    end
    return out
end

function get_xticks(t::Vector{T};
                    nticks::Int=5,
                    fmt_str::String=(T==Date ? "yyyy-mm-dd" : "yyyy-mm-ddTHH:MM:SS")) where T<:Dates.TimeType
    xticks_idx::Vector{Int} = round.(Int, range(1, stop=size(t,1), length=nticks))
    xticks_lab::Vector{String} = Dates.format.(t[xticks_idx], fmt_str)
    return xticks_idx, xticks_lab
end

@recipe function f(X::TS)
    #x = 1:size(X,1)
    xticks --> get_xticks(X.index)
    lab --> tslab(X.fields)
    X.values
end


text:
export futureValue, Dividend, CapitalizedDividend

struct Dividend
    amount::Number
    exDate::Number
    payDate::Number
    isProportional::Bool
    isKnown::Bool
end

struct CapitalizedDividend
    dividend::Dividend
    capitalizationFactor::Number #growth factor from the dividend ex date up to a maturity τ
end

function futureValue(cd::CapitalizedDividend)::Number
    return cd.dividend.amount*cd.capitalizationFactor
end


text:
using ProjectEuler
using Base.Test, BenchmarkTools

files = readdir(@__DIR__)
sort!(files)

@testset "Testing all existing problems" begin
    @testset "Testing problem $i" for i in 1:(length(files)-1)
        include(files[i])
    end
end

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler


text:
using PSIS
using CairoMakie
using Test

@testset "Makie.jl recipes" begin
    PSIS.plotting_backend!(:Makie)

    @testset "paretoshapeplot/plot" begin
        @testset "$f(::$T; showlines=$showlines) sz=$sz" for f in
                                                             (paretoshapeplot, Makie.plot),
            T in (PSISResult, Vector),
            showlines in (true, false),
            sz in (100, (10, 100))

            f === Makie.plot && T <: Vector && continue
            result = psis(randn(sz...))
            values = T <: PSISResult ? result : result.pareto_shape
            plt = f(values; showlines=showlines, linewidth=5)
            @test plt.plot isa PSIS.MakieRecipe.ParetoShapePlot
            plots = plt.plot.plots
            @test plt.axis.xlabel[] == "Parameter index"
            @test plt.axis.ylabel[] == "Pareto shape"
            @test length(plots) == 5
            @test plots[5] isa Scatter
            @test plots[5].converted[1][] ==
                Point2f0.(Base.OneTo(result.nparams), PSIS.as_array(result.pareto_shape))
            linestyles = [:dot, :dashdot, :dash, :solid]
            yvals = [0, 0.5, 0.7, 1]
            line_xmin, line_xmax = sz isa Tuple ? (0.91, 10.09) : (1, 1)
            for i in 1:4
                @test plots[i] isa LineSegments
                @test plots[i].converted[1][] ==
                    [Point2f0(line_xmin, yvals[i]), Point2f0(line_xmax, yvals[i])]
                @test plots[i].attributes[:linestyle][] == linestyles[i]
                @test plots[i].attributes[:linewidth][] == 5
                @test plots[i].attributes[:visible][] == showlines
            end
        end
    end

    @testset "paretoshapeplot!/plot!" begin
        @testset "$f!(::$T; showlines=$showlines) sz=$sz" for f! in (
                paretoshapeplot!, Makie.plot!
            ),
            T in (PSISResult, Vector),
            showlines in (true, false),
            sz in (100, (10, 100))

            f! === Makie.plot! && T <: Vector && continue
            result = psis(randn(sz...))
            values = T <: PSISResult ? result : result.pareto_shape
            fig = Figure()
            ax = Axis(fig[1, 1])
            plt = f!(values; showlines=showlines, linewidth=5)
            axis = fig.content[1]
            @test last(axis.scene.plots) isa PSIS.MakieRecipe.ParetoShapePlot
            plots = last(axis.scene.plots).plots
            @test axis.xlabel[] == "Parameter index"
            @test axis.ylabel[] == "Pareto shape"
            @test length(plots) == 5
            @test plots[5] isa Scatter
            @test plots[5].converted[1][] ==
                Point2f0.(Base.OneTo(result.nparams), PSIS.as_array(result.pareto_shape))
            linestyles = [:dot, :dashdot, :dash, :solid]
            yvals = [0, 0.5, 0.7, 1]
            line_xmin, line_xmax = sz isa Tuple ? (0.91, 10.09) : (1, 1)
            for i in 1:4
                @test plots[i] isa LineSegments
                @test plots[i].converted[1][] ==
                    [Point2f0(line_xmin, yvals[i]), Point2f0(line_xmax, yvals[i])]
                @test plots[i].attributes[:linestyle][] == linestyles[i]
                @test plots[i].attributes[:linewidth][] == 5
                @test plots[i].attributes[:visible][] == showlines
            end
        end
    end
end


text:
"""
    AbstractSlices{S,N} <: AbstractArray{S,N}

Supertype for arrays of slices into a parent array over some dimension(s),
returning views that select all the data from the other dimensions.

`parent` will return the parent array.
"""
abstract type AbstractSlices{T,N} <: AbstractArray{T,N} end

"""
    Slices{P,SM,AX,S,N} <: AbstractSlices{S,N}

An `AbstractArray` of slices into a parent array over specified dimension(s),
returning views that select all the data from the other dimension(s).

These should typically be constructed by [`eachslice`](@ref), [`eachcol`](@ref) or
[`eachrow`](@ref).

[`parent(s::Slices)`](@ref) will return the parent array.
"""
struct Slices{P,SM,AX,S,N} <: AbstractSlices{S,N}
    """
    Parent array
    """
    parent::P
    """
    A tuple of length `ndims(parent)`, denoting how each dimension should be handled:
      - an integer `i`: this is the `i`th dimension of the outer `Slices` object.
      - `:`: an "inner" dimension
    """
    slicemap::SM
    """
    A tuple of length `N` containing the [`axes`](@ref) of the `Slices` object.
    """
    axes::AX
end

unitaxis(::AbstractArray) = Base.OneTo(1)

function Slices(A::P, slicemap::SM, ax::AX) where {P,SM,AX}
    N = length(ax)
    S = Base._return_type(view, Tuple{P, map((a,l) -> l === (:) ? Colon : eltype(a), axes(A), slicemap)...})
    Slices{P,SM,AX,S,N}(A, slicemap, ax)
end

_slice_check_dims(N) = nothing
function _slice_check_dims(N, dim, dims...)
    1 <= dim <= N || throw(DimensionMismatch("Invalid dimension $dim"))
    dim in dims && throw(DimensionMismatch("Dimensions $dims are not unique"))
    _slice_check_dims(N,dims...)
end

@inline function _eachslice(A::AbstractArray{T,N}, dims::NTuple{M,Integer}, drop::Bool) where {T,N,M}
    _slice_check_dims(N,dims...)
    if drop
        # if N = 4, dims = (3,1) then
        # axes = (axes(A,3), axes(A,1))
        # slicemap = (2, :, 1, :)
        ax = map(dim -> axes(A,dim), dims)
        slicemap = ntuple(dim -> something(findfirst(isequal(dim), dims), (:)), N)
        return Slices(A, slicemap, ax)
    else
        # if N = 4, dims = (3,1) then
        # axes = (axes(A,1), OneTo(1), axes(A,3), OneTo(1))
        # slicemap = (1, :, 3, :)
        ax = ntuple(dim -> dim in dims ? axes(A,dim) : unitaxis(A), N)
        slicemap = ntuple(dim -> dim in dims ? dim : (:), N)
        return Slices(A, slicemap, ax)
    end
end
@inline function _eachslice(A::AbstractArray, dim::Integer, drop::Bool)
    _eachslice(A, (dim,), drop)
end

"""
    eachslice(A::AbstractArray; dims, drop=true)

Create a [`Slices`](@ref) object that is an array of slices over dimensions `dims` of `A`, returning
views that select all the data from the other dimensions in `A`. `dims` can either by an
integer or a tuple of integers.

If `drop = true` (the default), the outer `Slices` will drop the inner dimensions, and
the ordering of the dimensions will match those in `dims`. If `drop = false`, then the
`Slices` will have the same dimensionality as the underlying array, with inner
dimensions having size 1.

See also [`eachrow`](@ref), [`eachcol`](@ref), [`mapslices`](@ref) and [`selectdim`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

!!! compat "Julia 1.9"
     Prior to Julia 1.9, this returned an iterator, and only a single dimension `dims` was supported.

# Example

```jldoctest
julia> m = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> s = eachslice(m, dims=1)
3-element RowSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}}:
 [1, 2, 3]
 [4, 5, 6]
 [7, 8, 9]

julia> s[1]
3-element view(::Matrix{Int64}, 1, :) with eltype Int64:
 1
 2
 3

julia> eachslice(m, dims=1, drop=false)
3×1 Slices{Matrix{Int64}, Tuple{Int64, Colon}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}, 2}:
 [1, 2, 3]
 [4, 5, 6]
 [7, 8, 9]
```
"""
@inline function eachslice(A; dims, drop=true)
    _eachslice(A, dims, drop)
end

"""
    eachrow(A::AbstractVecOrMat) <: AbstractVector

Create a [`RowSlices`](@ref) object that is a vector of rows of matrix or vector `A`.
Row slices are returned as `AbstractVector` views of `A`.

See also [`eachcol`](@ref), [`eachslice`](@ref) and [`mapslices`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

!!! compat "Julia 1.9"
     Prior to Julia 1.9, this returned an iterator.

# Example

```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> s = eachrow(a)
2-element RowSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}}:
 [1, 2]
 [3, 4]

julia> s[1]
2-element view(::Matrix{Int64}, 1, :) with eltype Int64:
 1
 2
```
"""
eachrow(A::AbstractMatrix) = _eachslice(A, (1,), true)
eachrow(A::AbstractVector) = eachrow(reshape(A, size(A,1), 1))

"""
    eachcol(A::AbstractVecOrMat) <: AbstractVector

Create a [`ColumnSlices`](@ref) object that is a vector of columns of matrix or vector `A`.
Column slices are returned as `AbstractVector` views of `A`.

See also [`eachrow`](@ref), [`eachslice`](@ref) and [`mapslices`](@ref).

!!! compat "Julia 1.1"
     This function requires at least Julia 1.1.

!!! compat "Julia 1.9"
     Prior to Julia 1.9, this returned an iterator.

# Example

```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> s = eachcol(a)
2-element ColumnSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [1, 3]
 [2, 4]

julia> s[1]
2-element view(::Matrix{Int64}, :, 1) with eltype Int64:
 1
 3
```
"""
eachcol(A::AbstractMatrix) = _eachslice(A, (2,), true)
eachcol(A::AbstractVector) = eachcol(reshape(A, size(A, 1), 1))

"""
    RowSlices{M,AX,S}

A special case of [`Slices`](@ref) that is a vector of row slices of a matrix, as
constructed by [`eachrow`](@ref).

[`parent`](@ref) can be used to get the underlying matrix.
"""
const RowSlices{P<:AbstractMatrix,AX,S<:AbstractVector} = Slices{P,Tuple{Int,Colon},AX,S,1}

"""
    ColumnSlices{M,AX,S}

A special case of [`Slices`](@ref) that is a vector of column slices of a matrix, as
constructed by [`eachcol`](@ref).

[`parent`](@ref) can be used to get the underlying matrix.
"""
const ColumnSlices{P<:AbstractMatrix,AX,S<:AbstractVector} = Slices{P,Tuple{Colon,Int},AX,S,1}

IteratorSize(::Type{Slices{P,SM,AX,S,N}}) where {P,SM,AX,S,N} = HasShape{N}()
axes(s::Slices) = s.axes
size(s::Slices) = map(length, s.axes)

@inline function _slice_index(s::Slices, c...)
    return map(l -> l === (:) ? (:) : c[l], s.slicemap)
end

Base.@propagate_inbounds getindex(s::Slices{P,SM,AX,S,N}, I::Vararg{Int,N}) where {P,SM,AX,S,N} =
    view(s.parent, _slice_index(s, I...)...)
Base.@propagate_inbounds setindex!(s::Slices{P,SM,AX,S,N}, val, I::Vararg{Int,N}) where {P,SM,AX,S,N} =
    s.parent[_slice_index(s, I...)...] = val

parent(s::Slices) = s.parent

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
isautodifferentiable(alg::DEAlgorithm) = false
isadaptive(alg::DEAlgorithm) = true # Default to assuming adaptive, safer error("Adaptivity algorithm trait not set.")
isdiscrete(alg::DEAlgorithm) = false


text:
module Poisson

using CairoMakie
using DDF
using DifferentialForms
using IterativeSolvers
using LinearAlgebra
using SparseArrays
using SpecialFunctions
using StaticArrays
using WriteVTK

const CSI = "\e["

function vsplit(A::Array, sizes...)
    size(A, 1) == sum(sizes) ||
        throw(BoundsError("Excpected: size(A,1) == sum(sizes)"))
    @assert ndims(A) == 1       # TODO
    limits = (0, cumsum(sizes)...)
    # return Tuple((@view A[(limits[i] + 1):limits[i + 1]]) for i in sizes)
    return Tuple(A[(limits[i] + 1):limits[i + 1]] for i in 1:length(sizes))
end

################################################################################

function poisson(::Val{D}, levels::Int) where {D}
    D::Int
    S = Float64
    T = Float64

    println("Create manifold...")

    if D ≤ 2
        mfd = simplex_manifold(Val(D), S; optimize_mesh=false)
        for level in 1:levels
            mfd = refined_manifold(mfd; optimize_mesh=false)
        end
    else
        n = 2^levels
        nelts = Tuple(n for d in 1:D)
        mfd = large_hypercube_manifold(Val(D), S; nelts=nelts,
                                       optimize_mesh=true)
    end

    ############################################################################

    # DA:
    #     ⟨σ|τ⟩ - ⟨u|dτ⟩ = 0
    #     ⟨dσ|v⟩ + ⟨du|dv⟩ + ⟨p|v⟩ = ⟨f|v⟩
    #     ⟨u|q⟩ = 0
    #
    #     ⟨σ|τ⟩ - [u|τ] + ⟨δu|τ⟩ = 0
    #     ⟨dσ|v⟩ + [du|v] - ⟨δdu|v⟩ + ⟨p|v⟩ = ⟨f|v⟩
    #     ⟨u|q⟩ = 0
    #
    #     σ - [u] + δu = 0
    #     dσ + [du] - δdu + p = f
    #     u = 0
    #
    #     σ + δu = 0
    #     dσ + p = f

    # BH:
    #     -⋆f + dp = 0
    #     df       = 0

    # With von Neumann boundary conditions:
    #     L u == ρ
    #     B d u == 0
    #
    #     d u - f == 0
    #     δ f == ρ
    #     B f == 0

    # With Dirichlet boundary conditions:
    #     L u == ρ
    #     B u == 0
    #
    #     d u - f == 0
    #     δ f == ρ
    #     B u == 0

    println("Sample RHS...")

    x₀ = zero(Form{D,1,S}) .+ S(0.1)
    σ = S(0.2)
    function ρ⁼(x)
        local r = norm(x - x₀)
        local r′ = r / (sqrt(S(2)) * σ)
        local ρ = 1 / sqrt(2 * S(π) * σ^2)^D * exp(-r′^2)
        return Form{D,0}((ρ,))
    end
    function u⁼(x)
        local r = norm(x - x₀)
        local r′ = r / (sqrt(S(2)) * σ)
        if D == 1
            local u = σ / sqrt(2 * S(π)) * (exp(-r′^2) - 1) + r / 2 * erf(r′)
        elseif D == 2
            if r^8 ≤ eps(S)
                local u = 1 / (4 * S(π)) * r′^2 - 1 / (16 * S(π)) * r′^4 +
                          1 / (72 * S(π)) * r′^6
            else
                local u = -1 / (4 * S(π)) * (-MathConstants.eulergamma +
                           expinti(-r′^2) +
                           log(1 / r′^2))
            end
        elseif D == 3
            local u = -1 / (4 * S(π) * r) * erf(r′)
        elseif D == 4
            if r^8 ≤ eps(S)
                local u = -1 / (8 * S(π)^2 * σ^2) +
                          1 / (16 * S(π)^2 * σ^2) * r′^2 -
                          1 / (48 * S(π)^2 * σ^2) * r′^4 +
                          1 / (192 * S(π)^2 * σ^2) * r′^6
            else
                local u = 1 / (4 * S(π) * r^2) * (exp(-r′^2) - 1)
            end
        else
            @assert false
        end
        return Form{D,0}((u,))
    end

    ρ = sample(Fun{D,Pr,0,D,S,T}, ρ⁼, mfd)
    u₀ = sample(Fun{D,Pr,0,D,S,T}, u⁼, mfd)

    ############################################################################

    println("Determine boundaries...")

    # Boundary operator for vertices
    B0 = isboundary(Val(Pr), Val(0), mfd)

    N0 = zero(Op{D,Pr,0,Pr,0}, mfd)
    N1 = zero(Op{D,Pr,1,Pr,1}, mfd)
    E0 = one(Op{D,Pr,0,Pr,0}, mfd)
    E1 = one(Op{D,Pr,1,Pr,1}, mfd)
    n0 = zero(Fun{D,Pr,0,D,S}, mfd)
    n1 = zero(Fun{D,Pr,1,D,S}, mfd)

    N01 = zero(Op{D,Pr,0,Pr,1}, mfd)
    N10 = zero(Op{D,Pr,1,Pr,0}, mfd)

    E = [
        E0.values N01.values
        N10.values E1.values
    ]
    @assert issparse(E)
    @assert E * E == E

    ############################################################################

    println("Define operators...")

    # We want to solve this:
    #     [0   δ] [u] = [ρ]
    #     [-d  1] [f] = [0]
    #     A x == b
    # but need to handle boundary conditions. The boundary conditions are:
    #     [B  0] [u] = [u₀]
    #     [0  0] [f] = [0 ]
    #     B x == B c
    # Projecting these out of the linear system:
    #     (1-B) A x == (1-B) b
    #     A′ x = b′

    d = deriv(Val(Pr), Val(0), mfd)
    δ = coderiv(Val(Pr), Val(1), mfd)

    # The call to `sparse` is necessary for efficiency
    # <https://github.com/JuliaLang/julia/issues/38209>
    A = [
        N0.values δ.values
        -d.values sparse(E1.values)
    ]
    @assert issparse(A)
    b = [ρ.values.vec; n1.values.vec]

    B = [
        B0.values N01.values
        N10.values N1.values
    ]
    @assert issparse(B)
    @assert B * B == B
    c = [u₀.values.vec; n1.values.vec]

    A′ = (E - B) * A + B
    @assert issparse(A′)
    b′ = (E - B) * b + B * c

    ############################################################################

    println("Solve...")

    if D ≤ 2
        x = A′ \ b′

    else
        t0 = time_ns()

        Pl = Identity()

        x = zeros(T, size(A′, 1))
        solve = bicgstabl_iterator!(x, A′, b′; Pl=Pl)
        lastlog = t0
        iter = 0
        for _ in solve
            iter += 1
            t1 = time_ns()
            tlog = (t1 - lastlog) / 1.0e+9
            if tlog ≥ 1
                lastlog = t1
                tsol = (t1 - t0) / 1.0e+9
                print("\r    iter=$iter time=$(round(tsol; digits=3)) ",
                      "res=$(round(solve.residual; sigdigits=3))$(CSI)K")
            end
        end
        x = solve.x

        t1 = time_ns()
        tsol = (t1 - t0) / 1.0e+9
        print("\r    iter=$iter time=$(round(tsol; digits=3)) ",
              "res=$(round(solve.residual; sigdigits=3))$(CSI)K")
        println()
    end

    ############################################################################

    println("Analyse solution...")

    uv, fv = vsplit(x, nsimplices(mfd, 0), nsimplices(mfd, 1))
    u = Fun{D,Pr,0,D,S,T}(mfd, IDVector{0}(uv))
    f = Fun{D,Pr,1,D,S,T}(mfd, IDVector{1}(fv))

    # nint = norm((E0 - B0) * (laplace(u) - ρ), Inf)
    # nbnd = norm(B0 * (u - u₀), Inf)
    # println("    ‖residual[int]‖∞=$nint")
    # println("    ‖residual[bnd]‖∞=$nbnd")

    r = (E0 - B0) * (laplace(u) - ρ) + B0 * (u - u₀)
    nr = sqrt(norm(r)^2 / length(r))
    println("    ‖residual‖₂=$nr")

    e = u - u₀
    ne = norm(e, 1) / length(e)
    println("    ‖e‖₁=$ne")
    ne = sqrt(norm(e)^2 / length(e))
    println("    ‖e‖₂=$ne")
    ne = norm(e, Inf)
    println("    ‖e‖∞=$ne")

    ############################################################################

    if D ≤ 3
        if nsimplices(mfd, 0) ≤ 10000
            println("Plot result...")
            mkpath("poisson")
            # Can plot ρ, u, u₀, r, e
            plot_function(ρ, "poisson/poisson$(D)d-ρ.png")
            plot_function(u₀, "poisson/poisson$(D)d-u₀.png")
            plot_function(u, "poisson/poisson$(D)d-u.png")
            plot_function(r, "poisson/poisson$(D)d-r.png")
            plot_function(e, "poisson/poisson$(D)d-e.png")
        end
    end

    ############################################################################

    if D ≤ 3
        println("Write result to file...")

        points = [get_coords(mfd)[i][d]
                  for d in 1:D, i in axes(get_simplices(mfd, D), 1)]
        celltype = D == 1 ? VTKCellTypes.VTK_LINE :
                   D == 2 ? VTKCellTypes.VTK_TRIANGLE :
                   D == 3 ? VTKCellTypes.VTK_TETRA : nothing
        cells = [MeshCell(celltype,
                          SVector{D + 1}(Int(i)
                                         for i in
                                             sparse_column_rows(get_simplices(mfd,
                                                                              D),
                                                                j)))
                 for j in axes(get_simplices(mfd, D), 2)]

        mkpath("poisson")
        # append=false, ascii=true
        vtk_grid("poisson/poisson$(D)d", points, cells) do vtkfile
            vtkfile["ρ", VTKPointData()] = ρ.values.vec
            vtkfile["u₀", VTKPointData()] = u₀.values.vec
            vtkfile["u", VTKPointData()] = u.values.vec
            vtkfile["r", VTKPointData()] = r.values.vec
            vtkfile["e", VTKPointData()] = e.values.vec
            return nothing
        end
    end

    ############################################################################

    println("Done.")

    return nothing
end

poisson(::Val{1}) = poisson(Val(1), 5)
poisson(::Val{2}) = poisson(Val(2), 5)
poisson(::Val{3}) = poisson(Val(3), 6)
poisson(::Val{4}) = poisson(Val(4), 4)

poisson(D::Int) = poisson(Val(D))
poisson(D::Int, levels::Int) = poisson(Val(D), levels)

end


text:
using Plots, Colors, Compose, LazySets, Polyhedra, ConcaveHull
# using PyCall
include("adjacency_generator.jl")
gr()

# Used to convert obstacles to grid
function set2grid(set::Hyperrectangle)
    l = Tuple(Int.(low(set) .- 1))
    u = Tuple(Int.(high(set)))

    set_points = base_grid(u .- l) .+ [CartesianIndex(l)]
end

"""
    segment_blocked(::LineSegment, ::LazySet)

Returns true if the line and set intersect in a non-trivial way, i.e. at more than one point.
Naturally, this does not work on the intersection of two lines (since it is always at a point)
"""
function segment_blocked(line, obs)
    # any(iszero, radius_hyperrectangle(obs)) && error("The obstacle is infinitesimally thin, which will
    #                                                  break this function. obs radius: $(radius_hyperrectangle(obs))")

    if line.q ∈ obs || line.p ∈ obs
        return true
    end
    inter = intersection(line, obs)
    if isempty(inter) || length(vertices_list(inter)) <= 1
        return false
    else
        return true
    end
end

get_pos_tup(G::MetaGraph, i) = Tuple(G[i, :pos])
get_pos_vec(G::MetaGraph, i) = collect(Tuple(G[i, :pos]))

"""
    inf_visible(G,obs)

Takes in MetaGraph and array of obstacles. obs needs to have array of the Hyperrectangles
for the obstacles
"""
function inf_visible(G, obstacles)
    N = nv(G)
    # initialize a visibility vector for each point in the graph
    visibility = [CartesianIndex{2}[] for i in 1:N]
    for i in 1:N
        i_pos = Float64.(get_pos_vec(G, i))
        for p in (i+1):N
            p_pos = Float64.(get_pos_vec(G, p))

            vect = LineSegment(i_pos, p_pos)
            # @show vect
            isobstructed = any(segment_blocked(vect, o) for o in obstacles)
            if !isobstructed
                # If there is no obstruction,
                # then i can see p and p can see i
                push!(visibility[i], G[p, :pos])
                push!(visibility[p], G[i, :pos])
            end
        end
    end
    return visibility
end

"""
TOMERRRRRRRR

With a lot of nodes this takes a while to run (2.5 minutes when doing a 20x20)
Need to figure out how to give inputted nodes and output a vector with lots of
empty values and just the visible cells I care about in their correct positions.

Note: this would only be run once per env/obs pairing so it wouldn't be that big
of a deal.
"""
# Used to get all of the visible cells from all nodes assuming limited views
function limited_visible(G, obs, R, V_x)
    L_V_x = []
    for i in 1:nv(G)
        L_V_x_i = limited_visibility_node(G, obs, i, R, V_x)
        push!(L_V_x, collect(L_V_x_i))
    end
    return L_V_x
end

# Used to get all of the visible cells from a single node assuming limited views
function limited_visibility_node(G, obs, node, R, V_x)
    L_V_x_i = []
    if typeof(node) != CartesianIndex{2}
        V_x_i = V_x[node]
        node = get_prop(G, node, :pos)
    else
        i = G[node, :pos]
        V_x_i = V_x[node]
    end

    [push!(L_V_x_i, pos) for pos in V_x_i if (norm(Tuple(pos - node)) <= R) ]
    return L_V_x_i
end

# Takes in visible nodes (infinite or limited) and outputs all of the visible cells
# from given subset of nodes
function total_coverage(G, V_x, nodes, obs)
    T_V_x = []

    if typeof(nodes[1]) == CartesianIndex{2}
        node_idx = [G[node, :pos] for node in nodes]
    elseif typeof(nodes[1]) == Int64
        node_idx = [node for node in nodes]
    end

    for idx in node_idx
        push!(T_V_x, vcat(V_x[idx]))
    end
    T_V_x = unique(vcat(T_V_x...))

    grid_obs = []
    for ob in obs
        push!(grid_obs, vcat(set2grid(ob)...))
    end

    grid_obs = unique(vcat(grid_obs...))
    T_env = nv(G) - length(grid_obs)
    return T_V_x, T_env
end

# Plots visible cells given a Total Coverage visible set (works for one node too)
function plot_visible(G, T_V_x, nodes, obs, R = [])

    finalplot = graph2gridplot(G,obs)
    xlims0 = xlims(finalplot)
    ylims0 = ylims(finalplot)

    if typeof(nodes[1]) == CartesianIndex{2}
        nodes = [G[node, :pos] for node in nodes]
    end

    for pos in T_V_x
        scatter!(finalplot, Tuple(pos), color="lightgreen", markersize=5)
    end
    for node in nodes
        scatter!(finalplot, Tuple(G[node,:pos]), color="blue", markersize=5)
        if !isempty(R)
            plot!(finalplot, Ball2(Float64.(collect(Tuple(G[node,:pos]))), Float64.(R)),
                    color=:lightgreen, 1e-3, aspectratio=1, opacity=0.5)
        end
    end
    plot!(finalplot, xlims=xlims0, ylims=ylims0)
    png(finalplot, "test.png")

    return finalplot
end

Base.CartesianIndex(s::SVector) = CartesianIndex(Tuple(s))

function voronoi_cells(G, obstacles, nodes)
    N = length(nodes)
    V_c = [CartesianIndex{2}[] for i in 1:N]
    geo_dist = [Float64[] for i in 1:N]
    if typeof(nodes[1]) == CartesianIndex{2}
        nodeidx = [G[n, :pos] for n in nodes]
    elseif eltype(nodes) <: Vector{<:Number}
        nodeidx = nodes
    elseif eltype(nodes) <: SArray
        nodeidx = [G[CartesianIndex(c), :pos] for c in nodes]
    end

    for i in 1:N
        ds = dijkstra_shortest_paths(G, nodeidx[i])
        push!(geo_dist[i], ds.dists...)
    end
    for j in 1:nv(G)
        isobstacle = any(G[j, :pos] ∈ set2grid(o) for o in obstacles)
        if !isobstacle
            dists = []
            for k in 1:N
                push!(dists, geo_dist[k][j])
            end
            closest_node = findall(x-> x==minimum(dists), dists)
            for i in 1:length(closest_node)
                push!(V_c[closest_node[i]], G[j, :pos])
            end
        end
    end

    return V_c
end

function plot_voronoi(G, V_c, nodes, obs)
    N = length(nodes)
    finalplot = graph2gridplot(G,obs)
    xlims0 = xlims(finalplot)
    ylims0 = ylims(finalplot)
    colorchoices = distinguishable_colors(N)

    if typeof(nodes[1]) == CartesianIndex{2}
        nodeidx = [G[nodes[i], :pos] for i in 1:N]
    else
        nodeidx = nodes
    end

    for i in 1:N
        V_tup_conc = [collect(Float64.(Tuple(j) .+ (randn()/1000000))) for j in V_c[i]]
        V_tup = [collect(Float64.(Tuple(j))) for j in V_c[i]]
        P = concave_hull(vec(V_tup_conc))
        plot!(P, color=colorchoices[i], fill=true, opacity=0.5)
        for pos in V_c[i]
            scatter!(finalplot, Tuple(pos), color=colorchoices[i], markersize=5)
        end
        scatter!(finalplot, Tuple(G[nodeidx[i],:pos]), color="blue", markersize=5)
    end
    plot!(finalplot, xlims=xlims0, ylims=ylims0)
    png(finalplot, "test.png")
end

# Example workflow. From graph and obstacles to plot. THIS MIGHT TAKE A WHILE
# RUN DON'T SAY I DIDN'T WARN YOU.
function main()
    R = 4
    obs = [Hyperrectangle(low=[3, 3], high=[6, 8]),Hyperrectangle(low=[4, 9], high=[8, 10]),
            Hyperrectangle(low=[15,15], high=[16,19]),Hyperrectangle(low=[16,3],high=[18,10]),
            Hyperrectangle(low=[3,16],high=[10,18])]
    @time G = create_graph(20, 20)
    @time G = remove_obstacles(G,obs)
    nodes = [CartesianIndex((2, 3)),CartesianIndex((9, 2)),CartesianIndex((5,13)),
            CartesianIndex((14,11)),CartesianIndex((13,7))]
    # @time V_x = inf_visible(G, obs)
    @time L_V_x = limited_visible(G, obs, R, V_x)
    T_L_V_x, n_view = total_coverage(G, L_V_x, nodes, obs)
    fig = plot_visible(G, T_L_V_x, nodes, obs, R)
end


text:
# Baseline methods.
module BaselineMethods

importall Orchestra.Types
importall Orchestra.Util
import StatsBase: mode

export Baseline,
       Identity,
       fit!,
       transform!

# Baseline learner that by default assigns the most frequent label.
type Baseline <: Learner
  model
  options

  function Baseline(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      # Label assignment strategy.
      # Function that takes a label vector and returns the required output.
      :strategy => mode
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(bl::Baseline, instances::Matrix, labels::Vector)
  bl.model = bl.options[:strategy](labels)
end

function transform!(bl::Baseline, instances::Matrix)
  return fill(bl.model, size(instances, 1))
end

# Identity transformer passes the instances as is.
type Identity <: Transformer
  model
  options

  function Identity(options=Dict())
    default_options = Dict{Symbol, Any}()
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(id::Identity, instances::Matrix, labels::Vector)
  nothing
end

function transform!(id::Identity, instances::Matrix)
  return instances
end

end # module


text:
"""
    MoistThermodynamics

Moist thermodynamic functions, e.g., for air pressure (atmosphere equation
of state), latent heats of phase transitions, saturation vapor pressures, and
saturation specific humidities.

## AbstractParameterSet's

Many functions defined in this module rely on CLIMAParameters.jl.
CLIMAParameters.jl defines several functions (e.g., many planet
parameters). For example, to compute the mole-mass ratio:

```julia
using CLIMAParameters.Planet: molmass_ratio
using CLIMAParameters: AbstractEarthParameterSet
struct EarthParameterSet <: AbstractEarthParameterSet end
param_set = EarthParameterSet()
_molmass_ratio = molmass_ratio(param_set)
```

Because these parameters are widely used throughout this module,
`param_set` is an argument for many MoistThermodynamics functions.
"""
module MoistThermodynamics

using DocStringExtensions
using RootSolvers
using RootSolvers: AbstractTolerance

using CLIMAParameters: AbstractParameterSet
using CLIMAParameters.Planet
const APS = AbstractParameterSet

@inline q_pt_0(::Type{FT}) where {FT} = PhasePartition{FT}(FT(0), FT(0), FT(0))

include("states.jl")
include("relations.jl")
include("isentropic.jl")

end #module MoistThermodynamics.jl


text:
using PEPSAD
using ITensors

@testset "test local hamiltonian builder" begin
    Nx = 2
    Ny = 3
    sites = siteinds("S=1/2", Nx * Ny)
    sites = reshape(sites, Ny, Nx)
    H = mpo(Model("tfim"), sites; h = 1.0)
    H_local = localham(Model("tfim"), sites; h = 1.0)
    @test checklocalham(H_local, H, sites)
end


text:
__precompile__()

module CatViews

using Base.Cartesian
import Base.ReshapedArray
using Base.Iterators: repeated, flatten

export CatView, splitview, vecidx

include("catview.jl")
include("splitview.jl")

# end module
end


text:
function prettypercent(v, pad=false)
    if isnan(v)
        str = pad ? "--- " : "---"
    else
        i = round(Int, v)
        str = (i == 0 ? "0" : i > 999 ? ">999" : i < -999 ? "<-999" : @sprintf("%+d", i)) * "%"
    end
    return pad ? lpad(str, 5, " ") : str
end

const TIME_UNITS = [(1e3, 1, "ns"), (1e3, 1e3, "μs"), (1e3, 1e6, "ms"), (Inf, 1e9, "s")]

prettytime(v, pad=false) = pretty(TIME_UNITS, v, pad)

const MEMORY_UNITS = [(1e3, 1, "B"), (1e3, 1024, "KB"), (1e3, 1024^2, "MB"), (Inf, 1024^3, "GB")]

prettymemory(v, pad=false) = pretty(MEMORY_UNITS, v, pad)

function pretty(units, v, pad)
    _, scale, unit = units[findfirst(x -> round(Int, v/x[2]) < x[1], units)]
    str = three_significant_digits(v / scale) * " $unit"
    return pad ? lpad(str, 7, " ") : str
end

three_significant_digits(v) =
    round(v) >= 100 ? @sprintf("%.0f", v) :
    round(v * 10) >= 100 ? @sprintf("%.1f", v) : @sprintf("%.2f", v)

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
# Midiagonal, Ridiagonal matrices
# ============================================================

"""
R <: Ridiagonal ≡ R₂ ⋯ Rₙ₋₁Rₙ is block bidiagonal on subdiagonal
M <: Midiagonal ≡ block_diag(M₁, ⋯ ,Mₙ) is block diagonal

with R.data == [R₂, …, Rₙ]
with M.data == [M₁, ⋯ ,Mₙ]
"""

struct Ridiagonal{T,M<:AbstractMatrix{T}} <: VecchiaFactor{T}
    data::Vector{M} # sub diagonal matrices stored along the last dimension
end

struct Midiagonal{T,M<:AbstractMatrix{T}} <: VecchiaFactor{T}
    data::Vector{M} # diagonal matrices stored along the last dimension
end

const MiRi{T,M}    = Union{Midiagonal{T,M}, Ridiagonal{T,M}}

# Interface for VecchiaFactors
# ===================================

# bypass adjoint and inv
inv(M::Midiagonal)      = Midiagonal(map(inv, M.data))
adjoint(M::Midiagonal)  = Midiagonal(map(adjoint, M.data))

# merge products of Midiagonals
*(M1::Midiagonal, M2::Midiagonal) = Midiagonal(map(*, M1.data, M2.data))

# ldiv!
# ---------------------------------

# inv(M) * w
function ldiv!(M::A, w::AbstractVector)  where {A<:Midiagonal}
    rwB, wB = block_array(M, rw, w)
    for i = 1:length(wB)
        copyto!(rwB[i], M.data[i] \ wB[i])
    end
    return rw
end

# inv(M') * w
function ldiv!(Mᴴ::A, w::AbstractVector) where {A<:Adjoint{<:Any,<:Midiagonal}}
    M = parent(Mᴴ)
    rwB, wB = block_array(M, rw, w)
    for i = 1:length(wB)
        copyto!(rwB[i], M.data[i]' \ wB[i])
    end
    return rw
end

# inv(R) * w ≡ (R₂ ⋯ RₙRₙ₋₁)⁻¹ * w
function ldiv!(R::Ridiagonal, w::AbstractVector{T}) where T
    wbB = block_array(R, w)
    for i in 1:length(wbB)-1
        mul!(wbB[i+1], R.data[i], wbB[i], -1, true)
    end
    return w
end

# inv(R)' * w ≡ (Rₙ'Rₙ₋₁' ⋯ R₂')⁻¹ * w
function ldiv!(Rᴴ::Adjoint{<:Any,<:Ridiagonal}, w::AbstractVector{T}) where T
    R   = parent(Rᴴ)
    wbB = block_array(R, w)
    for i in length(wbB)-1:-1:1
        mul!(wbB[i], R.data[i]', wbB[i+1], -1, true)
    end
    return w
end

# lmul!
# ---------------------------------

# R * w ≡ (R₂ ⋯ RₙRₙ₋₁) * w
function lmul!(R::Ridiagonal, w::AbstractVector{T}) where T
    wbB = block_array(R, w)
    for i in length(wbB)-1:-1:1
        mul!(wbB[i+1], R.data[i], wbB[i], true, true)
    end
    return w
end

# R' * w ≡ (Rₙ'Rₙ₋₁' ⋯ R₂')*w
function lmul!(Rᴴ::Adjoint{<:Any, <:Ridiagonal}, w::AbstractVector{T}) where T
    R   = parent(Rᴴ)
    wbB = block_array(R, w)
    for i in 1:length(wbB)-1
        mul!(wbB[i], R.data[i]', wbB[i+1], true, true)
    end
    return w
end

# 3 and 5 arg mul!
# ---------------------------------

function mul!(rw::AbstractVector, R::A, w::AbstractVector) where {A<:Union{Ridiagonal, Adjoint{<:Any,<:Ridiagonal}}}
    copyto!(rw,w)
    lmul!(R, rw)
end

function mul!(rw::AbstractVector, M::A, w::AbstractVector) where {A<:Union{Midiagonal, Adjoint{<:Any,<:Midiagonal}}}
    mul!(rw, M, w, true, false)
end

function mul!(rw::AbstractVector, M::A, w::AbstractVector, α::Number, β::Number) where {A<:Midiagonal}
    rwB, wB = block_array(M, rw, w)
    for i = 1:length(wB)
        mul!(rwB[i], M.data[i], wB[i], α, β)
    end
    return rw
end

function mul!(rw::AbstractVector, Mᴴ::A, w::AbstractVector, α::Number, β::Number) where {A<:Adjoint{<:Any,<:Midiagonal}}
    M = parent(Mᴴ)
    rwB, wB = block_array(M, rw, w)
    for i = 1:length(wB)
        mul!(rwB[i], M.data[i]', wB[i], α, β)
    end
    return rw
end

# Non-interface methods
# ================================================

# TODO: where are these used ???
Ridiagonal(A::Ridiagonal) = A
Midiagonal(A::Midiagonal) = A
Ridiagonal{T}(A::Ridiagonal{T}) where {T} = A
Midiagonal{T}(A::Midiagonal{T}) where {T} = A
Ridiagonal{T}(A::Ridiagonal) where {T} = Ridiagonal{T}(A.data)
Midiagonal{T}(A::Midiagonal) where {T} = Midiagonal{T}(A.data)

nblocks(M::Midiagonal) = length(M.data)
nblocks(R::Ridiagonal) = length(R.data)+1

block_size(M::Midiagonal, d) = 0 < d::Integer <= 2 ? size.(M.data, 1) : 1
function block_size(R::Ridiagonal, d)
    if 0 < d::Integer <= 2
        bs = size.(R.data, 2)
        append!(bs, size(R.data[end],1))
    else
        return 1
    end
end

# used to prep incoming arrays for lmul! or ldiv!
function block_array(MR::MiRi, w::AbstractVector)
    bs  = block_size(MR,2)
    blocks(PseudoBlockArray(w, bs))
end
function block_array(MR::MiRi, w::AbstractVector, v::AbstractVector)
    bs  = block_size(MR,2)
    blocks(PseudoBlockArray(w, bs)), blocks(PseudoBlockArray(v, bs))
end

# appears to be useful for generating random Midiagonals or Ridiagonals
sizes_from_blocksides(::Type{<:Midiagonal}, bs::Vector{Int}) = [(bs[i],bs[i]) for i=1:length(bs)]
sizes_from_blocksides(::Type{<:Ridiagonal}, bs::Vector{Int}) = [(bs[i+1],bs[i]) for i=1:length(bs)-1]

# AbstractMatrix methods
# =================================================

size(MR::MiRi)    = (n = sum(block_size(MR,1)); (n,n))
size(MR::MiRi, d) = d::Integer <= 2 ? size(MR)[d] : 1

function getindex(R::Ridiagonal{T}, i::Integer, j::Integer) where T
    row_or_col_Ix = blockedrange(block_size(R,1))
    fbi = findblockindex(row_or_col_Ix, i)
    fbj = findblockindex(row_or_col_Ix, j)
    Block4i = fbi.I[1]
    Block4j = fbj.I[1]
    if Block4i  == Block4j
        wbi = blockindex(fbi)
        wbj = blockindex(fbj)
        wbi == wbj ? one(T) : zero(T)
    elseif  Block4i == Block4j + 1
        wbi = blockindex(fbi)
        wbj = blockindex(fbj)
        R.data[Block4j][wbi, wbj]
    else
        return zero(T)
    end
end

function getindex(M::Midiagonal{T}, i::Integer, j::Integer) where T
    row_or_col_Ix = blockedrange(block_size(M,1))
    fbi = findblockindex(row_or_col_Ix, i)
    fbj = findblockindex(row_or_col_Ix, j)
    Block4i = fbi.I[1]
    Block4j = fbj.I[1]
    if Block4i == Block4j
        wbi = blockindex(fbi)
        wbj = blockindex(fbj)
        M.data[Block4i][wbi, wbj]
    else
        return zero(T)
    end
end

# Convienence (not sure if we want these ?)
# =================================================

rand(::Type{A}, bs::Vector{Int}) where {T,A<:Midiagonal{T}} = Midiagonal(map(sz -> rand(T,sz), sizes_from_blocksides(A, bs)))
rand(::Type{A}, bs::Vector{Int}) where {T,A<:Ridiagonal{T}} = Ridiagonal(map(sz -> rand(T,sz), sizes_from_blocksides(A, bs)))

randn(::Type{A}, bs::Vector{Int}) where {T,A<:Midiagonal{T}} = Midiagonal(map(sz -> rand(T,sz), sizes_from_blocksides(A, bs)))
randn(::Type{A}, bs::Vector{Int}) where {T,A<:Ridiagonal{T}} = Ridiagonal(map(sz -> rand(T,sz), sizes_from_blocksides(A, bs)))

# an alternative to these is to simply overload map to reconsititue a Midiagonal?
sqrt(M::Midiagonal)      = Midiagonal(map(sqrt, M.data))
Hermitian(M::Midiagonal) = Midiagonal(map(x->Hermitian(x), M.data))
Symmetric(M::Midiagonal) = Midiagonal(map(x->Symmetric(x), M.data))
cholesky(M::Midiagonal)  = Midiagonal(map(x->cholesky(x).L, M.data))


text:
module SmartAsserts

export @smart_assert

include("smart_asserts.jl")

end


text:
# Declare&Document interface functions of PNML.jl
# Any method defined in this file should operate on `Any`.

"""
    pid(x) -> Symbol

Return pnml id symbol.
"""
function pid end

"""
    tag(x) -> Symbol

Return tag symbol.
"""
function tag end

"""
    name(x) -> String

Return name String. Default to empty string.
"""
function name end

"""
    has_xml(x) -> Bool

Return `true` if has XML attached. Defaults to `false`.
"""
function has_xml end
has_xml(x::Any) = hasproperty(x, :xml)

"""
    xmlnode(x) -> XMLNode

Return attached xml node.
"""
function xmlnode end

"""
    has_labels(x) -> Bool

Does x have any labels.
"""
function has_labels end

"""
    has_label(x, tag::Symbol) -> Bool

Does any label have a matching `tagvalue`.
"""
function has_label end

"""
    get_label(x, tag::Symbol) -> PnmlLabel

Return first label with a matching `tagvalue`.
"""
function get_label end

#--------------------------------------------
"""
$(TYPEDSIGNATURES)

Return the [`PnmlType`](@ref) subtype representing the flavor (or pntd) of this kind of
Petri Net Graph.

See also [`pnmltype`](@ref)
"""
function nettype end

"""
$(TYPEDSIGNATURES)

Return vector of pages.
"""
function pages end

"""
$(TYPEDSIGNATURES)

Return vector of places.
"""
function places end

"""
$(TYPEDSIGNATURES)

Return vector of place IDs.
"""
function place_ids end

"""
$(TYPEDSIGNATURES)
Return vector of transitions.
"""
function transitions end

"""
$(TYPEDSIGNATURES)

Return vector of arcs.
"""
function arcs end

"""
$(TYPEDSIGNATURES)
Return vector of reference places.
"""
function refplaces end

"""
$(TYPEDSIGNATURES)

Return vector of reference transitions.
"""
function reftransitions end
"""
$(TYPEDSIGNATURES)

Return `true` if any `arc` has `id`.
"""
function has_arc end

"""
$(TYPEDSIGNATURES)
Return arc with `id` if found, otherwise `nothing`.
"""
function arc end

"""
$(TYPEDSIGNATURES)

Return vector of arc ids.
"""
function arc_ids end

"""
$(TYPEDSIGNATURES)
Return vector of arcs that have a source or target of transition `id`.

See also [`src_arcs`](@ref), [`tgt_arcs`](@ref).
"""
function all_arcs end

"""
$(TYPEDSIGNATURES)

Return vector of arcs that have a source of transition `id`.

See also [`all_arcs`](@ref), [`tgt_arcs`](@ref).
"""
function src_arcs end

"""
$(TYPEDSIGNATURES)

Return vector of arcs that have a target of transition `id`.

See also [`all_arcs`](@ref), [`src_arcs`](@ref).
"""
function tgt_arcs end

"""
$(TYPEDSIGNATURES)
Return incription value of `arc`.
"""
function inscription end

"""
$(TYPEDSIGNATURES)
"""
function has_refP end

"""
$(TYPEDSIGNATURES)
"""
function has_refT end

"""
$(TYPEDSIGNATURES)
"""
function refplace_ids end

"""
$(TYPEDSIGNATURES)
"""
function reftransition_ids end

"""
Return reference place matching `id`.
$(TYPEDSIGNATURES)
"""
function refplace end

"""
Return reference transition matching `id`.
$(TYPEDSIGNATURES)
"""
function reftransition end

"""
$(TYPEDSIGNATURES)

Is there any place with `id`?
"""
function has_place end

"""
$(TYPEDSIGNATURES)

Return the place with `id`.
"""
function place end

"""
$(TYPEDSIGNATURES)

Return vector of place ids.
"""

"""
$(TYPEDSIGNATURES)

Return marking value of a place `p`.
"""
function marking end

"""
$(TYPEDSIGNATURES)

Return a labelled vector with key of place id and value of marking.
"""
function initialMarking end

"""
$(TYPEDSIGNATURES)

Is there a transition with `id`?
"""
function has_transition end

"""
$(TYPEDSIGNATURES)
"""
function transition end

"""
$(TYPEDSIGNATURES)
"""
function transition_ids end

"""
$(TYPEDSIGNATURES)

Return a labelled vector of condition values for net `s`. Key is transition id.
"""
function conditions end

"""
$(TYPEDSIGNATURES)

Return condition value of `transition`.
"""
function condition end


text:
function readphase(fn; rescale=true, keyargs...)
    phase = niread(fn; keyargs...)
    if phase.header.scl_slope == 0 || rescale
        minp, maxp = Float32.(approxextrema(phase))
        if isapprox(maxp - minp, 2π; atol=0.1) # no rescaling required
            return phase
        end
        minp, maxp = Float32.(approxextrema(phase.raw))
        if isapprox(maxp - minp, 2π; atol=0.1) # no rescaling required, but header wrong
            phase.header.scl_slope = 1
            phase.header.scl_inter = 0
        else # rescaling
            phase.header.scl_slope = 2pi / (maxp - minp)
            phase.header.scl_inter = -pi - minp * phase.header.scl_slope
        end
    end
    return phase
end

function readmag(fn; rescale=false, keyargs...)
    mag = niread(fn; keyargs...)
    if mag.header.scl_slope == 0 || rescale
        mini, maxi = Float32.(approxextrema(mag.raw))
        mag.header.scl_slope = 1 / (maxi - mini)
        mag.header.scl_inter = - mini * mag.header.scl_slope
    end
    return mag
end

Base.copy(x::NIfTI.NIfTI1Header) = NIfTI.NIfTI1Header([getfield(x, k) for k ∈ fieldnames(NIfTI.NIfTI1Header)]...)

function Base.similar(header::NIfTI.NIfTI1Header)
    hdr = copy(header)
    hdr.scl_inter = 0
    hdr.scl_slope = 1
    return hdr
end

header(v::NIVolume) = similar(v.header)

function savenii(image, name, writedir, header=nothing)
    if isnothing(writedir) return end
    if splitext(name)[2] != ".nii"
        name = name * ".nii"
    end
    savenii(image, joinpath(writedir, name); header)
end
"""
    savenii(image, filepath; header=nothing)
save the image at the path
Warning: MRIcro can only open images with types Int32, Int64, Float32, Float64
"""
function savenii(image::AbstractArray, filepath; header=nothing)
    vol = NIVolume([h for h in [header] if h !== nothing]..., image)
    niwrite(filepath, vol)
    return filepath
end

ConvertTypes = Union{BitArray, AbstractArray{UInt8}} #TODO debug NIfTI
MriResearchTools.savenii(image::ConvertTypes, args...;kwargs...) = savenii(Float32.(image), args...;kwargs...)

function write_emptynii(sz, path; datatype=Float64, header=NIVolume(zeros(datatype, 1)).header)
    header = copy(header)
    header.dim = Int16.((length(sz), sz..., ones(8-1-length(sz))...))
    header.datatype = NIfTI.nidatatype(datatype)
    header.bitpix = NIfTI.nibitpix(datatype)

    if isfile(path) rm(path) end
    file = open(path, "w")
    write(file, header)
    close(file)
    GC.gc()
    return niread(path; mmap=true, mode="r+")
end

mmtovoxel(sizemm, nii::NIVolume) = mmtovoxel(sizemm, nii.header)
mmtovoxel(sizemm, header::NIfTI.NIfTI1Header) = mmtovoxel(sizemm, header.pixdim)
mmtovoxel(sizemm, pixdim) = sizemm ./ pixdim

getcomplex(mag::NIVolume, phase::NIVolume) = getcomplex(mag.raw, phase.raw)


text:
# project.jl

export project, project!

"""
    project!(view, plan, image, viewidx)
SPECT projection of `image` into a single `view` with index `viewidx`.
The `view` must be pre-allocated but need not be initialized to zero.
"""
function project!(
    view::AbstractMatrix{<:RealU},
    image::AbstractArray{<:RealU, 3},
    plan::SPECTplan,
    viewidx::Int,
)
    # rotate image and mumap using multiple processors

    Threads.@threads for z = 1:plan.imgsize[3] # 1:nz
        thid = Threads.threadid() # thread id
        # rotate image in plan.imgr
        imrotate!(
            (@view plan.imgr[:, :, z]),
            (@view image[:, :, z]),
            plan.viewangle[viewidx],
            plan.planrot[thid],
        )
        # rotate mumap and store in plan.mumapr
        imrotate!(
            (@view plan.mumapr[:, :, z]),
            (@view plan.mumap[:, :, z]),
            plan.viewangle[viewidx],
            plan.planrot[thid],
        )
    end

    Threads.@threads for y = 1:plan.imgsize[2] # 1:ny
        thid = Threads.threadid() # thread id
        # account for half of the final slice thickness
        scale3dj!(plan.exp_mumapr[thid], plan.mumapr, y, -0.5)

        for j = 1:y
            plus3dj!(plan.exp_mumapr[thid], plan.mumapr, j)
        end

        broadcast!(*, plan.exp_mumapr[thid], plan.exp_mumapr[thid], - plan.dy)

        broadcast!(exp, plan.exp_mumapr[thid], plan.exp_mumapr[thid])

        # apply depth-dependent attenuation
        mul3dj!(plan.imgr, plan.exp_mumapr[thid], y)

        fft_conv!(
            (@view plan.add_img[:, y, :]),
            (@view plan.imgr[:, y, :]),
            (@view plan.psfs[:, :, y, viewidx]),
            plan.planpsf[thid],
        )
    end

    copy3dj!(view, plan.add_img, 1) # initialize accumulator
    for y in 2:plan.imgsize[2] # accumulate to get total view
        plus3dj!(view, plan.add_img, y)
    end

    # plan.add_img[2] # why does julia allocate (on heap!?) here?

    return view
end

"""
    project!(view, plan, image, thid, viewidx)
SPECT projection of `image` into a single `view` with index `viewidx`.
The `view` must be pre-allocated but need not be initialized to zero.
"""
function project!(
    view::AbstractMatrix{<:RealU},
    image::AbstractArray{<:RealU, 3},
    plan::SPECTplan,
    thid::Int,
    viewidx::Int,
)
    # rotate image and mumap using multiple processors

    for z = 1:plan.imgsize[3] # 1:nz
        # rotate image in plan.imgr
        imrotate!(
            (@view plan.imgr[thid][:, :, z]),
            (@view image[:, :, z]),
            plan.viewangle[viewidx],
            plan.planrot[thid],
        )
        # rotate mumap and store in plan.mumapr
        imrotate!(
            (@view plan.mumapr[thid][:, :, z]),
            (@view plan.mumap[:, :, z]),
            plan.viewangle[viewidx],
            plan.planrot[thid],
        )
    end

    for y = 1:plan.imgsize[2] # 1:ny
        # account for half of the final slice thickness
        scale3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], y, -0.5)

        for j = 1:y
            plus3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], j)
        end

        broadcast!(*, plan.exp_mumapr[thid], plan.exp_mumapr[thid], - plan.dy)

        broadcast!(exp, plan.exp_mumapr[thid], plan.exp_mumapr[thid])

        # apply depth-dependent attenuation
        mul3dj!(plan.imgr[thid], plan.exp_mumapr[thid], y)

        fft_conv!(
            (@view plan.add_img[thid][:, y, :]),
            (@view plan.imgr[thid][:, y, :]),
            (@view plan.psfs[:, :, y, viewidx]),
            plan.planpsf[thid],
        )

    end

    copy3dj!(view, plan.add_img[thid], 1) # initialize accumulator
    for y in 2:plan.imgsize[2] # accumulate to get total view
        plus3dj!(view, plan.add_img[thid], y)
    end

    # plan.add_img[2] # why does julia allocate (on heap!?) here?

    return view
end

"""
    project!(views, image, plan; index)
Project `image` into multiple `views` with indexes `index` (defaults to `1:nview`).
The 3D `views` array must be pre-allocated, but need not be initialized.
"""
function project!(
    views::AbstractArray{<:RealU,3},
    image::AbstractArray{<:RealU,3},
    plan::SPECTplan;
    index::AbstractVector{<:Int} = 1:plan.nview, # all views
)

    # loop over each view index
    if plan.mode === :fast
        Threads.@threads for i in index
            thid = Threads.threadid()
            project!((@view views[:,:,i]), image, plan, thid, i)
        end
    else
        for i in index
            project!((@view views[:,:,i]), image, plan, i)
        end
    end

    return views
end

"""
    views = project(image, plan ; kwargs...)
Convenience method for SPECT forward projector that allocates and returns views.
"""
function project(
    image::AbstractArray{<:RealU,3},
    plan::SPECTplan;
    kwargs...,
)
    views = Array{plan.T}(undef, plan.imgsize[1], plan.imgsize[3], plan.nview)
    project!(views, image, plan; kwargs...)
    return views
end

"""
    views = project(image, mumap, psfs, dy; interpmeth, kwargs...)
Convenience method for SPECT forward projector that does all allocation
including initializing `plan`.

In
* `image` : 3D array `(nx,ny,nz)`
* `mumap` : `(nx,ny,nz)` 3D attenuation map, possibly zeros()
* `psfs` : 4D PSF array
* `dy::RealU` : pixel size
Option
* `interpmeth` : `:one` or `:two`
"""
function project(
    image::AbstractArray{<:RealU, 3},
    mumap::AbstractArray{<:RealU, 3}, # (nx,ny,nz) 3D attenuation map
    psfs::AbstractArray{<:RealU, 4}, # (px,pz,ny,nview)
    dy::RealU;
    interpmeth::Symbol = :two,
    mode::Symbol = :fast,
#   nthread::Int = Threads.nthreads(), # todo: option for plan
    kwargs...,
)
    size(mumap) == size(image) || throw(DimensionMismatch("image/mumap size"))
    size(image,2) == size(psfs,3) || throw(DimensionMismatch("image/psfs size"))
    plan = SPECTplan(mumap, psfs, dy; interpmeth, mode, kwargs...)
    return project(image, plan; kwargs...)
end


text:
using Weave
using Test

# Make a dummy codehunk with figure
chunk = Weave.CodeChunk("plot(x)", 1, 1, "", Dict())
options = merge(Weave.rcParams[:chunk_defaults], chunk.options)
merge!(chunk.options, options)
chunk.figures = ["figs/figures_plot1.png"]

@test Weave.formatfigures(chunk, Weave.md2tex) == "\\includegraphics{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.tex) == "\\includegraphics[]{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.texminted) == "\\includegraphics[]{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![](figs/figures_plot1.png)\\ \n\n"
@test Weave.formatfigures(chunk, Weave.github) == "![](figs/figures_plot1.png)\n"
@test Weave.formatfigures(chunk, Weave.hugo) == "{{
}}"
@test Weave.formatfigures(chunk, Weave.multimarkdown) == "![][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png \n"
@test Weave.formatfigures(chunk, Weave.md2html) == "\n"

chunk.options[:out_width] = "100%"
@test Weave.formatfigures(chunk, Weave.adoc) == "image::figs/figures_plot1.png[width=100%]\n"
@test Weave.formatfigures(chunk, Weave.rst) == ".. image:: figs/figures_plot1.png\n   :width: 100%\n\n"

chunk.options[:fig_cap] = "Nice plot"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=100%]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![Nice plot](figs/figures_plot1.png){width=100%}\n"
@test Weave.formatfigures(chunk, Weave.md2tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=100%]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.md2html) == "
\n\n
Nice plot
\n
\n"
@test Weave.formatfigures(chunk, Weave.rst) == ".. figure:: figs/figures_plot1.png\n   :width: 100%\n\n   Nice plot\n\n"
@test Weave.formatfigures(chunk, Weave.multimarkdown) == "![Nice plot][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png width=100%\n"
@test Weave.formatfigures(chunk, Weave.adoc) == "image::figs/figures_plot1.png[width=100%,title=\"Nice plot\"]"

chunk.options[:label] = "somefig"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![Nice plot](figs/figures_plot1.png){width=100% #fig:somefig}\n"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=100%]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\label{fig:somefig}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.tex) == Weave.formatfigures(chunk, Weave.md2tex)


text:
@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end

function MKT_Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function MKT_AbstractTransistor(;name)
    @named d = Pin()
    @named g = Pin()
    @named s = Pin()
    sts = @variables vds(t)=1.0 ids(t)=1.0 vgs(t)=1.0
    eqs = [
           vds ~ d.v - s.v
           vgs ~ g.v - s.v
           0 ~ d.i + s.i + g.i
           g.i ~ 0
           ids ~ d.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), d, g,s)
end

function MKT_OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function MKT_Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function MKT_Transistor(;name, W=1.0,L=1.0,m=QuadraticModel())
    @named tr = AbstractTransistor()
    @unpack vds, ids, vgs = tr
    ps = @parameters K=K Vt=Vt L=L
    eqs = [
           ids ~ K*(vgs-Vt)^2*(1+vds*L)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), tr)
end

function MKT_Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function MKT_AcVoltage(;name, V = 1.0, Vac = 0.1, f=1e3)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           v ~ V + Vac*sin(2*π*f*t)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function MKT_ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end


text:
using Test
include("ExtendableSparseTest.jl")

ExtendableSparseTest.benchmark()

@time begin

    @test ExtendableSparseTest.check(m=10,n=10,nnz=5)
    @test ExtendableSparseTest.check(m=100,n=100,nnz=500,nsplice=2)
    @test ExtendableSparseTest.check(m=1000,n=1000,nnz=5000,nsplice=3)

    @test ExtendableSparseTest.check(m=20,n=10,nnz=5)
    @test ExtendableSparseTest.check(m=200,n=100,nnz=500,nsplice=2)
    @test ExtendableSparseTest.check(m=2000,n=1000,nnz=5000,nsplice=3)

    @test ExtendableSparseTest.check(m=10,n=20,nnz=5)
    @test ExtendableSparseTest.check(m=100,n=200,nnz=500,nsplice=2)
    @test ExtendableSparseTest.check(m=1000,n=2000,nnz=5000,nsplice=3)

    for irun=1:10
        m=rand((1:10000))
        n=rand((1:10000))
        nnz=rand((1:10000))
        nsplice=rand((1:5))
        @test ExtendableSparseTest.check(m=m,n=n,nnz=nnz,nsplice=nsplice)
    end

end


text:
function step_sim!(θ, p_θ, dt)
    p_θ[] = p_θ[] - 9.8*sin(θ[])*dt
    θ[]   = θ[]   + p_θ[]*dt

    if (θ[] > π)
        θ[] = θ[] - 2*π
    elseif (θ[] < -π)
        θ[] = θ[] + 2*π
    end
end

function advance_sim!(θ, p_θ, dt, Δt)

  Δt_partial = zero(dt)

  while Δt_partial <= Δt
    step_sim!(θ, p_θ, dt)
    Δt_partial += dt
  end

  return Δt_partial - Δt
end

function sim!(θ, p_θ; dt=1/24/4, fps = 24, scene)

    frame_duration = 1/fps
    t_overshoot = zero(dt)
    limiter = Limiter(frame_duration)

    while scene.events.window_open[]
        t_overshoot = advance_sim!(θ, p_θ, dt, frame_duration - t_overshoot)
        limiter()
    end
end

max_stars_repo_name: tlqmj/phasespace

max_issues_repo_name: tlqmj/phasespace

max_forks_repo_name: tlqmj/phasespace


text:
function padding(lambda, den)
    tran_lambda = cast(lambda, Float64)
    tran_den = cast(den, Float64)
    lambda_pad = tf.pad(tran_lambda, [nPml (nPml+nPad); nPml nPml], constant_values=1.0/3.0*3500.0^2*2200.0/1e6)
    den_pad = tf.pad(tran_den, [nPml (nPml+nPad); nPml nPml], constant_values=2200.0)
    return lambda_pad, den_pad
end

# global loss

function fwi_sep(iSur)
  if !isdir("./$(args["version"])/FWI_stage$iSur/")
      mkdir("./$(args["version"])/FWI_stage$iSur/")
  end
  if iSur <= 2
    lambda_init = den .* (cp_nopad.^2 .- 2.0 .* cp_nopad.^2 ./3.0) /1e6
    den_init = den
  else
    lp = readdlm("./$(args["version"])/FWI_stage$(iSur-1)/loss.txt")
    Lp = Int64((lp[end,1]))
    lambda_init = readdlm("./$(args["version"])/FWI_stage$(iSur-1)/Lambda$Lp.txt")
    den_init = readdlm("./$(args["version"])/FWI_stage$(iSur-1)/Den$Lp.txt")
    # den_init = den
  end

  tf_lambda_inv = Variable(lambda_init, dtype=Float64)
  tf_den_inv = Variable(den_init, dtype=Float64)
  tf_lambda_inv_pad, tf_den_inv_pad = padding(tf_lambda_inv, tf_den_inv)

  shot_id_points = Int32.(trunc.(collect(LinRange(0, length(x_src)-1, nGpus+1))))
  shot_jump = 1
  # NOTE Compute FWI loss
  para_fname = "./$(args["version"])/para_file$iSur.json"

  loss = constant(0.0)
  for i = 1:nGpus
      # global loss
      tf_shot_ids = constant(collect(shot_id_points[i]:shot_jump:shot_id_points[i+1]), dtype=Int32)
      loss += fwi_op(tf_lambda_inv_pad, tf_shear_pad, tf_den_inv_pad, tf_stf, tf_gpu_id_array[i], tf_shot_ids0, para_fname)
  end
  gradLambda = gradients(loss, tf_lambda_inv)
  gradDen= gradients(loss, tf_den_inv)

  # loss = fwi_op(tf_lambda_inv_pad, tf_shear_pad, tf_den_inv_pad, tf_stf, tf_gpu_id_array[1], tf_shot_ids0, para_fname)

  # Optimization

  __cnt = 0
  __l = 0
  __Lambda = zeros(nz,nx)
  __gradLambda = zeros(nz,nx)
  __Den = zeros(nz,nx)
  __gradDen = zeros(nz,nx)
  # invK = zeros(m,n)
  function print_loss(l, Lambda, Den, gradLambda, gradDen)
      # global __l, __Lambda, __gradLambda, __Den, __gradDen
      if mod(__cnt,1)==0
          println("\niter=$__iter, eval=$__cnt, current loss=",l)
          # println("a=$a, b1=$b1, b2=$b2")
      end
      __cnt += 1
      __l = l
      __Lambda = Lambda
      __gradLambda = gradLambda
      __Den = Den
      __gradDen = gradDen
  end

  __iter = 0
  function print_iter(rk)
      # global __iter
      if mod(__iter,1)==0
          println("\n************* ITER=$__iter *************\n")
      end
      __iter += 1
      open("./$(args["version"])/FWI_stage$iSur/loss.txt", "a") do io
          println("\n outer_iter=$__iter, current loss=", __l)
          writedlm(io, Any[__iter __l])
      end
      open("./$(args["version"])/FWI_stage$iSur/Lambda$__iter.txt", "w") do io
          writedlm(io, __Lambda)
      end
      open("./$(args["version"])/FWI_stage$iSur/gradLambda$__iter.txt", "w") do io
          writedlm(io, __gradLambda)
      end
      open("./$(args["version"])/FWI_stage$iSur/Den$__iter.txt", "w") do io
          writedlm(io, __Den)
      end
      open("./$(args["version"])/FWI_stage$iSur/gradDen$__iter.txt", "w") do io
          writedlm(io, __gradDen)
      end
  end

  config = tf.ConfigProto()
  config.intra_op_parallelism_threads = 24
  config.inter_op_parallelism_threads = 24
  sess = Session(config=config); init(sess);

  lambda_lb = 5800.0
  lambda_ub = 9000.0
  opt = ScipyOptimizerInterface(loss, var_list=[tf_lambda_inv, tf_den_inv], var_to_bounds=Dict(tf_lambda_inv=> (lambda_lb, lambda_ub),tf_den_inv=> (2100.0, 2200.0)), method="L-BFGS-B", options=Dict("maxiter"=> 100, "ftol"=>1e-12, "gtol"=>1e-12))

# #   lambda_lb = (2500.0^2 - 2.0 * 3000.0^2/3.0) * 2500.0 /1e6
# #   lambda_ub = (4000.0^2 - 2.0 * 3000.0^2/3.0) * 2500.0 /1e6
#   lambda_lb = 5800.0
#   lambda_ub = 9000.0
#   opt = ScipyOptimizerInterface(loss, var_list=[tf_lambda_inv], var_to_bounds=Dict(tf_lambda_inv=> (lambda_lb, lambda_ub)), method="L-BFGS-B", options=Dict("maxiter"=> 100, "ftol"=>1e-12, "gtol"=>1e-12))
  @info "Optimization Starts..."
  ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=print_iter, fetches=[loss,tf_lambda_inv,tf_den_inv,gradLambda,gradDen])
end


text:
# In case more abstraction is needed
abstract type AbstractSLM{T} end

"""
```
SLM
```

is a port of Shape Language Modeling (SLM) by John D'Errico,
who implements least squares spline modeling for curve fitting.

See https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling for details about the SLM toolbox.

Note that only the features described in the Keywords section have been ported for now.
Other features may be added on as-needed basis.

The main constructor has the form

```
SLM(x::AbstractVector{T}, y::AbstractVector{T}; calculate_stats::Bool = false,
    verbose::Symbol = :low, kwargs...) where {T <: Real}
```

### Inputs
- The vectors `x` and `y` are the points with which the least-squares spline will be fit.

### Keywords
- `degree::Int = 3`: Degree of the spline. A degree 3 spline is a piecewise cubic Hermite spline.
- `scaling::Bool = true`: If true, scale the `y` values and coefficients to avoid numerical issues
- `knots::Int = 6`: Number of knots for the spline
- `C2::Bool = true`: If true, enforce that the spline satisfies C2 continuity
- `λ::T = 1e-4`: Regularization parameter for how smooth the spline is (smaller λ ⇒ smoother spline)
- `left_value::T = NaN`: Controls the spline's value at the left end point. If it is a `NaN`, then the
    spline's value at the left end point is not controlled, i.e. it can be any number.
- `right_value::T = NaN`: Controls the spline's value at the right end point. If it is a `NaN`, then the
    spline's value at the right end point is not controlled, i.e. it can be any number.
- `min_value::T = NaN`: Sets the spline's global minimum value. If it is a `NaN`, then
    no minimum value is imposed.
- `max_value::T = NaN`: Set the spline's global maximum value. If it is a `NaN`, then the
    no maximum value is imposed.
- `min_max_sample_points::AbstractVector{T} =
    [.017037, .066987, .14645, .25, .37059, .5, .62941, .75, .85355, .93301, .98296]`:
    Sample points used to enforce global min/max values. The default points are Chebyshev nodes.
- `increasing::Bool = false`: If true, the spline will monotonically increase over the interpolation domain.
- `decreasing::Bool = false`: If true, the spline will monotonically decrease over the interpolation domain.
- `increasing_intervals::AbstractMatrix{T} = []`: If nonempty, each row of the n × 2 matrix `increasing_intervals` specifies
    the left and right end points of an interval over which the spline monotonically increases.
- `decreasing_intervals::AbstractMatrix{T} = []`: If nonempty, each row of the n × 2 matrix `decreasing_intervals` specifies
    the left and right end points of an interval over which the spline monotonically decreases.
- `concave_up_intervals::AbstractMatrix{T} = []`: If nonempty, each row of the n × 2 matrix `concave_up_intervals` specifies
    the left and right end points of an interval over which the spline is concave up (i.e. convex in the normal sense).
- `concave_down_intervals::AbstractMatrix{T} = []`: If nonempty, each row of the n × 2 matrix `concave_down_intervals` specifies
    the left and right end points of an interval over which the spline is concave down (i.e. concave in the normal sense).
- `init::AbstractVector{T} = []`: If nonempty, `init` is just as an initial guess for the coefficients of the spline.
- `calculate_stats::Bool = false`: If true, calculate statistics about the spline regression (e.g. R²). Currently,
    no statistics are calculated, so setting `calculate_stats` to true does nothing.
- `verbose::Symbol`: Verbosity of information printed during construction of an SLM object.
"""
mutable struct SLM{T} <: AbstractSLM{T}
    stats::NamedTuple
    type::Symbol
    x::AbstractVector{T}
    y::AbstractVector{T}
    knots::AbstractVector{T}
    coef::AbstractArray{T}
    extrapolation::Symbol
end

function Base.show(io::IO, slm::AbstractSLM{T}) where {T <: Real}
    @printf io "SLM with element type %s\n" string(T)
    @printf io "degree: %i\n" get_stats(slm)[:degree]
    @printf io "knots:  %i\n" length(get_stats(slm)[:knots])
end

# Access functions
get_stats(slm::AbstractSLM) = slm.stats
get_knots(slm::AbstractSLM) = slm.knots
get_type(slm::AbstractSLM) = slm.type
get_x(slm::AbstractSLM) = slm.x
get_y(slm::AbstractSLM) = slm.y
get_coef(slm::AbstractSLM) = slm.coef
get_extrapolation(slm::AbstractSLM) = slm.extrapolation
eltype(slm::AbstractSLM) = slm.stats

function getindex(slm::AbstractSLM, x::Symbol)
    if x == :stats
        get_stats(slm)
    elseif x == :knots
        get_knots(slm)
    elseif x == :x
        get_x(slm)
    elseif x == :y
        get_y(slm)
    elseif x == :coef
        get_coef(slm)
    elseif x == :type
        get_type(slm)
    elseif x == :extrapolation || x == :extrap
        get_extrapolation(slm)
    else
        error("type " * typeof(slm) * " has no field " * string(x))
    end
end

# Main user interface for constructing an SLM object
function SLM(x::AbstractVector{T}, y::AbstractVector{T}; calculate_stats::Bool = false,
             verbose::Symbol = :low, kwargs...) where {T <: Real}

    @assert length(x) == length(y) "x and y must be the same size"
    if verbose == :high
        calculate_stats = true # Statistics will be calculated if verbose is high
    end
    kwargs = Dict{Symbol, Any}(kwargs)

    if haskey(kwargs, :increasing_intervals)
        @assert size(kwargs[:increasing_intervals], 2) == 2 "The matrix for the keyword increasing_intervals must be n × 2."
    end
    if haskey(kwargs, :decreasing_intervals)
        @assert size(kwargs[:decreasing_intervals], 2) == 2 "The matrix for the keyword decreasing_intervals must be n × 2."
    end
    if haskey(kwargs, :concave_up_intervals)
        @assert size(kwargs[:concave_up_intervals], 2) == 2 "The matrix for the keyword concave_up_intervals must be n × 2."
    end
    if haskey(kwargs, :concave_down_intervals)
        @assert size(kwargs[:concave_down_intervals], 2) == 2 "The matrix for the keyword concave_down_intervals must be n × 2."
    end

    # Remove nans
    to_remove = isnan.(x) .| isnan.(y)
    if any(to_remove)
        x = x[to_remove]
        y = y[to_remove]

        if haskey(kwargs, :weights)
            error("Weights are not implemented currently.")
            kwargs[:weights] = kwargs[:weights][to_remove]
        end
    end

    # Additional checks
    if haskey(kwargs, :weights)
        @assert length(kwargs[:weights]) == length(x)
    end

    # Add default keyword arguments
    default_slm_kwargs!(kwargs, T)

    # Scale y. This updates the kwargs
    ŷ = scale_problem!(x, y, kwargs)

    slm = if kwargs[:degree] == 0
        error("degree 0 has not been implemented")
    elseif kwargs[:degree] == 1
        error("degree 1 has not been implemented")
    elseif kwargs[:degree] == 3
        SLM_cubic(x, ŷ, kwargs[:y_scale], kwargs[:y_shift]; init = kwargs[:init],
                  nk = kwargs[:knots], C2 = kwargs[:C2], λ = kwargs[:λ], increasing = kwargs[:increasing],
                  decreasing = kwargs[:decreasing], increasing_intervals = kwargs[:increasing_intervals],
                  decreasing_intervals = kwargs[:decreasing_intervals],
                  concave_up = kwargs[:concave_up], concave_down = kwargs[:concave_down],
                  concave_up_intervals = kwargs[:concave_up_intervals],
                  concave_down_intervals = kwargs[:concave_down_intervals],
                  left_value = kwargs[:left_value], right_value = kwargs[:right_value],
                  min_value = kwargs[:min_value], max_value = kwargs[:max_value],
                  min_max_sample_points = kwargs[:min_max_sample_points],
                  use_sparse = kwargs[:use_sparse], use_lls = kwargs[:use_lls])
    else
        error("degree $(kwargs[:degree]) has not been implemented")
    end

    if calculate_stats
        # Need to write a function that calculates stats like total df or R²
    end

    if kwargs[:scaling]
        coef = get_coef(slm)
        coef[:, 1] .-= kwargs[:y_shift]
        coef       ./= kwargs[:y_scale]
    end

    if verbose == :high
        @info "Model Statistics Report"
        println("Number of data points:      $(length(y))")
        println("Scale factor applied to y   $(y_scale)")
        println("Shift applied to y          $(y_shift)")
        # println("Total degrees of freedom:   $(get_stats(slm)[:total_df])")
        # println("Net degrees of freedom:     $(get_stats(slm)[:net_df])")
        # println("R-squared:                  $(get_stats(slm)[:R2])")
        # println("Adjusted R-squared:         $(get_stats(slm)[:R2_adj])")
        # println("RMSE:                       $(get_stats(slm)[:RMSE])")
        # println("Range of prediction errors: $(get_stats(slm)[:error_range])")
        # println("Error quartiles (25%, 75%): $(get_stats(slm)[:quartiles])")
    end

    return slm
end

function SLM_cubic(x::AbstractVector{T}, y::AbstractVector{T}, y_scale::T, y_shift::T;
                   init::AbstractVector{T} = Vector{T}(undef, 0),
                   nk::Int = 6, C2::Bool = true, λ::T = 1e-4, increasing::Bool = false, decreasing::Bool = false,
                   increasing_intervals::AbstractMatrix{T} = Matrix{T}(undef, 0, 0),
                   decreasing_intervals::AbstractMatrix{T} = Matrix{T}(undef, 0, 0),
                   concave_up::Bool = false, concave_down::Bool = false,
                   concave_up_intervals::AbstractMatrix{T} = Matrix{T}(undef, 0, 0),
                   concave_down_intervals::AbstractMatrix{T} = Matrix{T}(undef, 0, 0),
                   left_value::T = NaN, right_value::T = NaN,
                   min_value::T = NaN, max_value::T = NaN,
                   min_max_sample_points::AbstractVector{T} = [.017037, .066987, .1465, .25, .37059,
                                                               .5, .62941, .75, .85355, .93301, .98296],
                   use_sparse::Bool = false, use_lls::Bool = true) where {T <: Real}

    nₓ = length(x)

    # Choose knots
    knots = choose_knots(nk, x)
    dknots = diff(knots)
    if any(dknots .== 0.)
        error("Knots must be distinct.")
    end

    ### Calculate coefficients

    ## Set up
    nc = 2 * nk
    Mineq = zeros(0, nc)
    rhsineq = Vector{T}(undef, 0)
    Meq = zeros(0, nc)
    rhseq = Vector{T}(undef, 0)

    ## Build design matrix

    # Bin data so that xbin has nₓ and xbin specifies into which bin each x value falls
    xbin = bin_sort(x, knots)

    # design matrix
    Mdes = construct_design_matrix(x, knots, dknots, xbin, nₓ, nk, nc; use_sparse = use_sparse)
    rhs = y

    ## Regularizer
    Mreg = construct_regularizer(dknots, nk; use_sparse = use_sparse)
    rhsreg = zeros(T, nk)

    ## C2 continuity across knots
    if C2
        Meq, rhseq = C2_matrix(nk, nc, dknots, Meq, rhseq; use_sparse = use_sparse)
    end

    ## Left and right values
    if !isnan(left_value)
        Meq, rhseq = set_left_value(left_value, nc, Meq, rhseq; use_sparse = use_sparse)
    end

    if !isnan(right_value)
        Meq, rhseq = set_right_value(right_value, nc, nk, Meq, rhseq; use_sparse = use_sparse)
    end

    # Global minimum and maximum values
    if !isnan(min_value)
        Mineq, rhsineq = set_min_value(min_value, nk, nc, dknots, Mineq, rhsineq; sample_points = min_max_sample_points,
                                       use_sparse = use_sparse)
    end

    if !isnan(max_value)
        Mineq, rhsineq = set_max_value(max_value, nk, nc, dknots, Mineq, rhsineq; sample_points = min_max_sample_points,
                                       use_sparse = use_sparse)
    end

    ## Monotonicity restrictions
    @assert !(increasing && decreasing) "Only one of increasing and decreasing can be true"

    monotone_settings = Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}}(undef, 0)
    total_monotone_intervals = 0

    # Monotone increasing
    if increasing
        @assert isempty(increasing_intervals) && isempty(decreasing_intervals) "The spline cannot be monotone increasing and " *
            "have nonempty entries for the keywords increasing_intervals and/or decreasing_intervals"
        total_monotone_intervals += monotone_increasing!(monotone_settings, nk)
    end

    # Increasing intervals: nx2 array where each row is an interval over which we have increasing
    if !isempty(increasing_intervals)
        total_monotone_intervals += increasing_intervals_info!(monotone_settings, knots, increasing_intervals, nk)
    end

    # Monotone decreasing
    if decreasing
        @assert isempty(increasing_intervals) && isempty(decreasing_intervals) "The spline cannot be monotone decreasing and " *
            "have nonempty entries for the keywords increasing_intervals and/or decreasing_intervals"
        total_monotone_intervals += monotone_decreasing!(monotone_settings, nk)
    end

    # Decreasing intervals
    if !isempty(decreasing_intervals)
        total_monotone_intervals += decreasing_intervals_info!(monotone_settings, knots, decreasing_intervals, nk)
    end

    # Add inequalities enforcing monotonicity
    if !isempty(monotone_settings)
        Mineq, rhsineq = construct_monotonicity_matrix(monotone_settings, nc, nk, dknots, total_monotone_intervals, Mineq, rhsineq;
                                                       use_sparse = use_sparse)
    end

    ## Concavity
    @assert !(concave_up && concave_down) "Only one of concave_up and concave_down can be true"

    curvature_settings = Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}}(undef, 0)

    if concave_up
        @assert isempty(concave_up_intervals) && isempty(concave_down_intervals) "The spline cannot be concave up and " *
            "have nonempty entries for the keywords concave_up_intervals and/or concave_down_intervals"
        concave_up_info!(curvature_settings)
    end

    if concave_down
        @assert isempty(concave_up_intervals) && isempty(concave_down_intervals) "The spline cannot be concave down and " *
            "have nonempty entries for the keywords concave_up_intervals and/or concave_down_intervals"
        concave_down_info!(curvature_settings)
    end

    if !isempty(concave_up_intervals)
        concave_up_intervals_info!(curvature_settings, concave_up_intervals)
    end

    if !isempty(concave_down_intervals)
        concave_down_intervals_info!(curvature_settings, concave_down_intervals)
    end

    # Add inequalities enforcing curvature
    if !isempty(curvature_settings)
        Mineq, rhsineq = construct_curvature_matrix(curvature_settings, nc, nk, knots, dknots, Mineq, rhsineq; use_sparse = use_sparse)
    end

    # scale equalities for unit absolute row sum
    if !isempty(Meq)
        rs = Diagonal(vec(1. ./ sum(abs.(Meq), dims = 2)))
        Meq = rs * Meq
        rhseq = rs * rhseq
    end

    # scale equalities for unit absolute row sum
    if !isempty(Mineq)
        rs = Diagonal(vec(1. ./ sum(abs.(Mineq), dims = 2)))
        Mineq = rs * Mineq
        rhsineq = rs * rhsineq
    end

    # Currently, we only implement the standard regularizer parameter, while
    # SLM allows matching a specified RMSE and cross-validiation, too.
    coef = solve_slm_system(Mdes, rhs, Mreg, rhsreg, λ,
                            Meq, rhseq, Mineq, rhsineq; init = init,
                            use_lls = use_lls)
    coef = reshape(coef, nk, 2)

    ## Calculate model statistics
    # Currently, we just add degree and knots to the statistics NamedTuple
    stats = (degree = 3, knots = knots,
             y_scale = y_scale, y_shift = y_shift)

    # Unpack coefficients into the result structure
    return SLM{T}(stats, :cubic, x, y, knots, coef, :none)
end


text:
using Turing

# https://turing.ml/dev/tutorials/7-poissonregression/
@model poisson_regression(traps, complaints) = begin
    α ~ Normal(log(4), .1)
    β ~ Normal(-0.25, .1)
    for i in eachindex(traps)
        λ = exp(α + β * traps[i])
        complaints[i] ~ Poisson(λ)
    end
end