text:
# Strongly-local code for minimizing the HCL objective
# Implemented with the thresholded linear hyperedge splitting penalty.

include("Helper_Functions.jl")
include("maxflow.jl")
include("../../src/SparseCard.jl")

"""
SparseCardHyperLocal: minimizes Hypergraph Local Conductance with a cardinality-based
                hypergraph cut function. This version of the code, unlike the
                original HyperLocal, uses a general graph expansion technique
                that allows you to define an approximate sparse reduced graph.

H:          Binary indicence matrix for hypergraph
Ht:         Transpose of binary incidence matrix
EdgesW:     For each hyperedge, a set of cardinality-based splitting penalties.
                This assumes symmetric splitting penalties. The first penalty should
                be w[1] = w_0 = 0, (zero penalty when zero nodes on small side of cut).
                Then one penalty is given for each possible # nodes on small side.
sparsityeps: The parameter controling how good the reduction has to be.
                0 leads to a perfect but denser reduction.
                A larger parameter leads to sparser graph.
order:      Order (number of nodes) in each hyperedge
d:          Degree vector, d[v] = number of hyperedges a node is in
R:          Set of nodes in seed/reference set
epsilon:    Locality parameter, must exceed vol(R)/vol(bar{R})
Rs_local:   Within the R vector, which are strong seed nodes? (I.e., need to be in output)
"""
function SparseCardHyperLocal(H::SparseMatrixCSC{Float64,Int64},Ht::SparseMatrixCSC{Float64,Int64},
    EdgesW::Vector{Vector{Float64}},Edges::Vector{Vector{Int64}},sparsityeps::Float64,
    order::Vector{Int64},d::Vector{Float64},R::Vector{Int64},
    epsilon::Float64,Rs_local::Vector{Int64},localflag::Bool=true)

    m,n = size(H)
    # Edges = incidence2elist(H)
    volA = sum(d)
    volR = sum(d[R])
    # @assert(volR <= volA/2)
    Rstrong = R[Rs_local]
    # Check Locality Parameter
    fR = volR/(volA - volR)
    # @show fR, volR, volA
    if epsilon < fR
        println("Locality parameter epsilon was set too small.
        Setting it to lower bound of $fR. Computations will not be local.")
        epsilon = fR
        localflag = false
    end
    A = 0; N = 0;
    if localflag

        if volA*epsilon/volR < 10
            println("Note that vol(R)/epsilon = O(vol(G)).
            For these parameters \nit may be faster to run the algorithm
            without the locality setting.")
        end

    else
        # A = tl_expansion_inc(H,order,delta)
        A = SymmetricCard_reduction(Edges,EdgesW,n,epsilon,false)
        N = round(Int64,size(A,1))
    end

    # Store useful sets
    # Rn = hyper_neighborhood(H,R)    # get the immediate neighbors of R...
    # Rn = setdiff(Rn,R)                # ...but we exclude R itself
    Rn = get_immediate_neighbors(H,Ht,R)
    Rc = setdiff(1:n,R)               # Complement set of R
    nR = length(R)

    # condR,volR, cutR = tl_cond(H,R,d,delta,volA,order)
    cutR = gen_hypergraph_cut(Edges,EdgesW,R,n)
    condR = cutR/(min(volR,volA-volR))

    println("\nRunning HyperLocal")
    println("----------------------------------------")
    println("Epsilon = $epsilon \t")
    println("|R| = $nR, cond(R) = $condR")
    println("-------------------------------------------------------------------------")

    S_best = R
    a_best = condR
    a_old = condR
    still_improving = true
    Iteration = 1
    while still_improving

        still_improving = false

        stepstart = time()
        if localflag
            S_new = HyperLocal_Step_General(H,Ht,EdgesW,sparsityeps,order,R,Rn,a_best,epsilon,d,Rs_local)
        else
            S_new = HLC_Step_General(A,R,Rc,a_best,epsilon,N,d,n,Rs_local)
        end
        stime = round(time()-stepstart,digits=1)

        # a_new = hlc_tl(H,order,R,S_new,d,volA,epsilon,delta)

        a_new = hlc_general(Edges,EdgesW,R,S_new,d,epsilon)
        if a_new < a_old
            still_improving = true
            S_best = S_new
            nS = length(S_best)
            a_old = a_new
            a_best = a_new
            println("Iter $Iteration: |S| = $nS, lcond(S) = $a_new, min-cut took $stime seconds")
        else
            println("Iter $Iteration: Algorithm converged. Last min-cut took $stime sec")
            println("-------------------------------------------------------------------------")
        end
        Iteration += 1
    end

    return S_best, a_best
end

function hlc_general(Edges::Vector{Vector{Int64}},EdgesW::Vector{Vector{Float64}},
    R::Vector{Int64},S::Vector{Int64},d::Vector{Float64},epsilon::Float64)
    """
    For a set S in a hypergraph, return the hypergraph local conductance
    score with general splitting penalty.
    """
    volS = sum(d[S])
    RnS = intersect(R,S)
    volRnS = sum(d[RnS])
    # cut = tl_cut(H,S,delta,order)
    n = length(d)
    cut = gen_hypergraph_cut(Edges,EdgesW,S,n)

    lcond = cut/((1+epsilon)*volRnS - epsilon*volS)

    return lcond

end

# A non-local version of the min-cut code that works by calling the same
# subroutine, but on the entire graph all at once
function HLC_Step_General(A::SparseMatrixCSC{Float64,Int64},R::Vector{Int64},Rbar::Vector{Int64},
    alpha::Float64, epsilon::Float64, N::Int64, d::Vector{Float64},n::Int64,Rs_local::Vector{Int64})

        Rstrong = R[Rs_local]
        # Directly set up the flow matrix
        sVec = zeros(N)
        tVec = zeros(N)
        sVec[R] .= alpha*d[R]
        sVec[Rstrong] .= N^2
        tVec[Rbar] .= alpha*epsilon*d[Rbar]
        F = maxflow(A,sVec,tVec,0)
        Src = source_nodes_min(F)[2:end].-1
        S = intersect(1:n,Src)

        return S
end

# Strongly-local subroutine for computing a minimum s-t cut
# This uses the thresholded linear splitting function for each hyperegde
function HyperLocal_Step_General(H::SparseMatrixCSC{Float64,Int64},Ht::SparseMatrixCSC{Float64,Int64},
    EdgesW::Vector{Vector{Float64}},sparsityeps::Float64,
    order::Vector{Int64}, R::Vector{Int64},Rn::Vector{Int64},alpha::Float64,
    epsilon::Float64,d::Vector{Float64},Rs_local::Vector{Int64})

    # Map from local node indices to global node indices
    Local2Global = [R; Rn]

    n = length(d)

    # Keep track of which nodes are in the local hypergraph L
    inL = zeros(Bool,n)
    inL[Local2Global] .= true

    # Number of nodes in the local graph
    Lsize = length(Local2Global)

    # Complete nodes = nodes whose hyperedge set in the local hypergraph
    #                   is the same as their global hyperedge set
    # Incomplete nodes = everything else in the local hypergraph
    #                   (must be a neighbor of a complete node)
    #
    # Initialize the complete set to be R
    # Incomplete set is R-complement
    C_global = R
    I_global = Rn

    # Indices, in the local graph, of complete and incomplete nodes
    C_local = collect(1:length(R))
    I_local = collect(length(R)+1:Lsize)
    R_local = collect(1:length(R))
    Rstrong_local = R_local[Rs_local]

    # Get the set of hyperedges to expand around.
    # At first this is every hyperedge that touches
    # a node from R.
    Hc = H[:,C_global]
    rp_c = Hc.rowval
    # ci_c = Hc.colptr
    L_edges = unique(rp_c)

    # Binary indicence matrix for the local hypergraph (without terminal edges)
    HL = H[L_edges,Local2Global]
    order_L = order[L_edges]

    # Expand into a directed graph
    # A_L = tl_expansion_inc(HL,order_L,delta)
    Edges_L = incidence2elist(HL)
    EdgesW_L = EdgesW[L_edges]
    A_L = SymmetricCard_reduction(Edges_L,EdgesW_L,Lsize,sparsityeps,false)
    N_L = size(A_L,1)           # includes auxiliary nodes
    n_L = length(Local2Global)  # number of non-auxiliary nodes in A_L

    # Find the first mincut, which can be done by calling HLC_Step
    # with localized objects
    S_local = HLC_Step_General(A_L,C_local,I_local,alpha,epsilon,N_L,d[Local2Global],n_L,Rstrong_local)

    # Find nodes to "expand" around:
    #   any nodes in the cut set tha are "incomplete" still
    E_local = intersect(S_local,I_local)
    E_global = Local2Global[E_local]

    # ne = length(E_global)
    # println("There are $ne new nodes to expand on")

    # As long as we have new nodes to expand around, we haven't yet found
    # the global minimum s-t cut, so we continue.
    while length(E_local) > 0

        # Update which nodes are complete and which are incomplete
        C_local = [C_local; E_local]
        C_global = Local2Global[C_local]

        # Take these away from I_local
        I_local = setdiff(I_local,E_local)

        # This is better
        Nbs_of_E = get_immediate_neighbors(H,Ht,E_global)
        Lnew = setdiff(Nbs_of_E,Local2Global)
        numNew = length(Lnew)
        # Update the set of indices in L
        Local2Global = [Local2Global; Lnew]

        # Store local indices for new nodes added to L
        Lnew_local = collect((Lsize+1):(Lsize+numNew))
        Lsize = length(Local2Global)

        # These are going to be "incomplete" nodes
        I_local = [I_local; Lnew_local]
        I_global = Local2Global[I_local]

        # Now we have a new set of complete and incomplete edges,
        # we do the same thing over again to find a localize min-cut
        Hc = H[:,C_global]
        rp_c = Hc.rowval
        # ci_c = Hc.colptr
        L_edges = unique(rp_c)

        # Binary indicence matrix for the local hypergraph (without terminal edges)
        HL = H[L_edges,Local2Global]
        order_L = order[L_edges]

        # Expand into a directed graph
        # A_L = tl_expansion_inc(HL,order_L,delta)
        Edges_L = incidence2elist(HL)
        EdgesW_L = EdgesW[L_edges]
        A_L = SymmetricCard_reduction(Edges_L,EdgesW_L,Lsize,sparsityeps,false)
        N_L = size(A_L,1)           # includes auxiliary nodes
        n_L = length(Local2Global)  # number of non-auxiliary nodes in A_L

        # Find the first mincut, which can be done by calling HLC_Step
        # with localized objects
        R_bar_l = setdiff(1:n_L,R_local)
        S_local = HLC_Step_General(A_L,R_local,R_bar_l,alpha,epsilon,N_L,d[Local2Global],n_L,Rstrong_local)

        # Find nodes to "expand" around:
        #   any nodes in the cut set tha are "incomplete" still
        E_local = intersect(S_local,I_local)
        E_global = Local2Global[E_local]
        # ne = length(E_global)
        # println("There are $ne new nodes to expand on")
    end

    return Local2Global[S_local]
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
@testset "Slope Limiter" begin

using ConservationLawsDiffEq

for limiter in (GeneralizedMinmodLimiter(), MinmodLimiter())
    printstyled("Limiter: ", string(typeof(limiter)) ,"\n"; color = :blue)
    @test limiter(1.0,2.0) == 1.0
    @test limiter(-1.0,2.0) == 0.0
    @test limiter(-1.0,-2.0) == -1.0
end

limiter  = OsherLimiter()
    printstyled("Limiter: ", string(typeof(limiter)) ,"\n"; color = :blue)
    @test limiter(1.0,2.0) == 1.0
    @test limiter(-1.0,2.0) == 0.0
    @test limiter(-1.0,-2.0) == 0.0

limiter = SuperbeeLimiter()
    printstyled("Limiter: ", string(typeof(limiter)) ,"\n"; color = :blue)
    @test limiter(1.0,2.0) == 2.0
    @test limiter(-1.0,2.0) == 0.0
    @test limiter(-1.0,-2.0) == 0.0
end

ext: jl


text:
using BinaryBuilder

# Collection of sources required to build FastTransforms
name = "FastTransforms"
version = v"0.3.3"
sources = [
    ArchiveSource("https://github.com/MikaelSlevinsky/FastTransforms/archive/v$(version).tar.gz",
                  "4566fd59d29f4bff4d68814334d8bc603056932d462e40bb22acaacbf353763a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/FastTransforms-*
if [[ ${target} == x86_64-* ]] || [[ ${target} == i686-* ]]; then
    export MSSE=-msse
    export MSSE2=-msse2
    export MAVX=-mavx
    export MFMA=-mfma
    if [[ ${target} == *-w64-mingw32 ]]; then
        export CFLAGS="-O3 -mavx -fno-asynchronous-unwind-tables "
    else
        export CFLAGS="-O3 -mavx "
        export MAVX512F=-mavx512f
    fi
else
    export CFLAGS="-O3 "
fi
if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
    SYMBOL_DEFS=()
    SYMBOLS=(dgemm dtrmm dtrmv dtrsm sgemm strmm strsm ztrmm)
    for sym in ${SYMBOLS[@]}; do
        SYMBOL_DEFS+=("-Dcblas_${sym}=cblas_${sym}64_")
    done
    CFLAGS+=${SYMBOL_DEFS[@]}
    BLAS=openblas64_
else
    BLAS=openblas
fi
if [[ ${target} == *apple* ]]; then
    export FT_OPENMP="-fopenmp=libgomp "
fi
make assembly
make lib FT_PREFIX=${prefix} FT_BLAS=${BLAS} FT_FFTW_WITH_COMBINED_THREADS=1
mv -f libfasttransforms.${dlext} ${libdir}
"""

platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libfasttransforms", :libfasttransforms),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FFTW_jll"),
    Dependency("MPFR_jll", v"4.0.2"),
    Dependency("OpenBLAS_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
module FieldArraysTests

using Test
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField, MockBasis
using InplaceArrays.Fields: Valued
using FillArrays
using InplaceArrays.TensorValues

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
f = MockField{d}(v)
fx = fill(v,np)
∇fx = fill(VectorValue(v,0.0),np)

l = 10
af = Fill(f,l)
ax = fill(x,l)
afx = fill(fx,l)
a∇fx = fill(∇fx,l)
test_array_of_fields(af,ax,afx,grad=a∇fx)

ag = apply_to_field_array(bcast(+),af,af)
gx = fill(v+v,np)
∇gx = fill(VectorValue(v+v,0.0),np)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

w = 2.0
aw = fill(w,l)
ag = apply_to_field_array(bcast(+),af,aw)
gx = fill(v+w,np)
∇gx = fill(VectorValue(v,0.0),np)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

l = 10
af = Fill(f,l)
ax = Fill(x,l)
ag = apply_to_field_array(bcast(+),af,af)
r1 = evaluate(ag,ax)
@test isa(r1,Fill)

np = 4
p = Point(1,2)
x = fill(p,np)
v = 2.0
d = 2
ndof = 8
wi = 3.0
w = fill(wi,ndof)
r = fill(v+wi,np,ndof)
f = MockBasis{d}(v,ndof)
∇fx = evaluate(∇(f),x)
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = apply_to_field_array(bcast(+),af,aw)
agx = fill(r,l)
a∇gx = fill(∇fx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

v = 2.0
d = 2
wi = 3.0
w = fill(wi,ndof)
r = fill(v+wi,np,ndof)
f = MockField{d}(v)
∇r = fill(VectorValue(v,0.0),np,ndof)
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = apply_to_field_array(bcast(+),af,aw)
agx = fill(r,l)
a∇gx = fill(∇r,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ContainerServiceMasterProfile <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    dnsPrefix::Any # spec type: Union{ Nothing, String } # spec name: dnsPrefix
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn

    function ContainerServiceMasterProfile(;count=1, dnsPrefix=nothing, fqdn=nothing)
        o = new()
        validate_property(ContainerServiceMasterProfile, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(ContainerServiceMasterProfile, Symbol("dnsPrefix"), dnsPrefix)
        setfield!(o, Symbol("dnsPrefix"), dnsPrefix)
        validate_property(ContainerServiceMasterProfile, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        o
    end
end # type ContainerServiceMasterProfile

const _property_map_ContainerServiceMasterProfile = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"), Symbol("dnsPrefix")=>Symbol("dnsPrefix"), Symbol("fqdn")=>Symbol("fqdn"))
const _property_types_ContainerServiceMasterProfile = Dict{Symbol,String}(Symbol("count")=>"Int32", Symbol("dnsPrefix")=>"String", Symbol("fqdn")=>"String")
Base.propertynames(::Type{ ContainerServiceMasterProfile }) = collect(keys(_property_map_ContainerServiceMasterProfile))
Swagger.property_type(::Type{ ContainerServiceMasterProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceMasterProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceMasterProfile }, property_name::Symbol) =  _property_map_ContainerServiceMasterProfile[property_name]

const _allowed_ContainerServiceMasterProfile_count = ["1", "3", "5"]

function check_required(o::ContainerServiceMasterProfile)
    (getproperty(o, Symbol("dnsPrefix")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ContainerServiceMasterProfile }, name::Symbol, val)
    if name === Symbol("count")
        Swagger.validate_param(name, "ContainerServiceMasterProfile", :enum, val, _allowed_ContainerServiceMasterProfile_count)
    end
end


text:
const BLOCK_INFO_FIELD_STOP = UInt64(0)
const BLOCK_INFO_FIELD_OVERFLOWS = UInt64(1)
const BLOCK_INFO_FIELD_BUCKET_NUM = UInt64(2)

struct BlockInfo
    is_overflows::Bool
    bucket_num::Int32

    BlockInfo() = new(false, -1)
    BlockInfo(is_overflows, bucket_num) = new(is_overflows, bucket_num)
end

function chread(sock::ClickHouseSock, ::Type{BlockInfo})::BlockInfo
    is_overflows = false
    bucket_num = -1

    while (field_num = UInt64(chread(sock, VarUInt))) != BLOCK_INFO_FIELD_STOP
        if field_num == BLOCK_INFO_FIELD_OVERFLOWS
            is_overflows = chread(sock, Bool)
        elseif field_num == BLOCK_INFO_FIELD_BUCKET_NUM
            bucket_num = chread(sock, Int32)
        else
            throw("Unknown block info field")
        end
    end

    BlockInfo(is_overflows, bucket_num)
end

function chwrite(sock::ClickHouseSock, x::BlockInfo)
    # This mirrors what the C++ client does.
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_OVERFLOWS))
    chwrite(sock, x.is_overflows)
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_BUCKET_NUM))
    chwrite(sock, x.bucket_num)
    chwrite(sock, VarUInt(BLOCK_INFO_FIELD_STOP))
end

struct Column
    name::String
    type::String
    data::Any
end

Base.:(==)(a::Column, b::Column) =
    a.name == b.name && a.type == b.type && a.data == b.data

# We can't just use chread here because we need the size to be passed
# in from the `Block` decoder that holds the row count.
function read_col(sock::ClickHouseSock, num_rows::VarUInt)::Column
    name = chread(sock, String)
    type_name = chread(sock, String)

    type = parse_typestring(type_name)
    data = if UInt64(num_rows) == 0
        result_type(type)(undef, 0)
    else
        try
            read_state_prefix(sock, type)
            read_col_data(sock, num_rows, type)
        catch e
            if e isa ArgumentError
                error("Error while reading col $(name) ($(type)): $(e.msg)")
            else
                rethrow(e)
            end
        end
    end
    Column(name, type_name, data)
end

function chwrite(sock::ClickHouseSock, x::Column)
    chwrite(sock, x.name)
    chwrite(sock, x.type)

    try
        type = parse_typestring(x.type)
        write_state_prefix(sock, type)
        write_col_data(sock, x.data, type)
    catch e
        if e isa ArgumentError
            error("Error while writing col $(x.name) ($(x.type)): $(e.msg)")
        else
            rethrow(e)
        end
    end
end

struct Block
    temp_table::String
    block_info::BlockInfo
    num_columns::VarUInt
    num_rows::VarUInt
    columns::Array{Column}
end

function chread(sock::ClickHouseSock, ::Type{Block})::Block
    temp_table = chread(sock, String)
    block_info = chread(sock, BlockInfo)
    num_columns = chread(sock, VarUInt)
    num_rows = chread(sock, VarUInt)
    columns = [read_col(sock, num_rows) for _ ∈ 1:UInt64(num_columns)]
    Block(temp_table, block_info, num_columns, num_rows, columns)
end

function chwrite(sock::ClickHouseSock, x::Block)
    chwrite(sock, x.temp_table)
    chwrite(sock, x.block_info)
    chwrite(sock, x.num_columns)
    chwrite(sock, x.num_rows)
    for x ∈ x.columns
        chwrite(sock, x)
    end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 1968f0a2-fde4-11ea-3a5e-c1fc5176627f
begin
	using Plots
	using IntervalRootFinding
	using IntervalArithmetic

	include("../AcousticPropagation.jl")
end

# ╔═╡ 34c9dff0-fdd8-11ea-383b-b3c3d5e3bbe8
md"""
# Beam Tracing
```math
\newcommand{\Par}[1]{\left({#1}\right)}
\newcommand{\Brace}[1]{\left\{{#1}\right\}}
\newcommand{\Brack}[1]{\left[{#1}\right]}
```
"""

# ╔═╡ 3e7e5a20-fdd9-11ea-3a3b-31fe9237e71e
md"""
In this document, we take the step forward from having solved the eikonal, first transport equation, and dynamic ray equations which gives us the following instantaneous ray parameters:
* Length $s$
* Range $r(s)$
* Depth $z(s)$
* Tangent range component $\xi = \dfrac{\cos\Par{\theta(s)}}{c(s)}$
* Tangent depth component $\zeta = \dfrac{\sin(\theta)}{c(s)}$
* Time $\tau\Par{s} = \dfrac{s}{c(s)}$
* Pace $p(s) = \dfrac{1}{c(s)}$
* Spreading $q(s)$
* Angle $\theta(s)$
* Speed $c(s)$
"""

# ╔═╡ 6edbdfb0-fddb-11ea-1b8d-3dd87aa01386
md"""
These variables allow us to approximate the ray series
```math
p\Par{\vec{x}} = e^{i\omega\tau\Par{\vec{x}}} \sum_{j = 0}^\infty \frac{A_j\Par{\vec{x}}}{\Par{i\omega}^j}.
```
Along the length of the ray,
```math
p(s) = e^{i\omega\tau(s)} \sum_{j = 0}^\infty \frac{A_j(s)}{\Par{i\omega}^j}
```
"""

# ╔═╡ d2c2a7f0-fde2-11ea-16ad-ab1114259f56
md"""
The initial conditions for $p(s)$ and $q(s)$ allow us to implement different beam tracing methods.
"""

# ╔═╡ fdf32220-fde1-11ea-25e6-9db949c15ec2
md"""
## Gaussian Beams
The initial conditions for the pace $p(s)$ and spreading $q(s)$ for calculating Gaussian beams are
```math
\begin{align*}
p(0) &= 1 \\
q(0) &= \frac{i\omega W(0)^2}{2}
\end{align*}
```
"""

# ╔═╡ 75c3a120-fde3-11ea-1682-87e3dc7aee71
md"""
The equation for the beam is given by
```math
p^\textrm{beam}(s, n) = A \sqrt{\frac{c(s)}{rq(s)}} \exp\Brace{-i\omega\Brack{\tau(s) + \frac{p(s)}{q(s)}\frac{n^2}{2}}}
```
"""

# ╔═╡ 73d6a8c0-fde4-11ea-237d-15e08ebdaa1c
begin
	z_c = [0., 300., 1200., 2e3, 5000.]
	c = [1520, 1500, 1515, 1495, 1545.]
	R = 250e3
	Z = z_c[end]
	ocn = AcousticPropagation.Medium(z_c, c, R)

	r₀ = 0.0
	z₀ = 20.0
	src = AcousticPropagation.Entity(r₀, z₀)

	bty = AcousticPropagation.Boundary(Z)

	θ_crit = acos(ocn.c(r₀, z₀)/ocn.c(r₀, Z))
	θ₀ = θ_crit*rand((0.5, 1.0))

	nothing
end

# ╔═╡ 57c41320-fde4-11ea-3fd8-351192415917
ray = AcousticPropagation.Ray(θ₀, src, ocn, bty);

# ╔═╡ b1c28750-fde7-11ea-24f0-23050e9d0d0e
plot(ray.Sol, vars = (1, 2),
	xaxis = "Range (m)",
	yaxis = ("Depth (m)", :flip))

# ╔═╡ 22d78b60-fde9-11ea-2a88-3118125cf30b
begin
	f = 200
	ω = 2π*f
end

# ╔═╡ ecf09052-fde8-11ea-27cf-c397d368a09c
begin
	A = 1
	p_beam(s) = A * sqrt(ray.c(s)/ray.r(s)/ray.q(s)) * exp(-im*ω)
	TL_beam(s) = -20log10(abs(p_beam(s)))
end

# ╔═╡ 54b58be0-fdea-11ea-356e-b9f30d0f15be
s = range(0, ray.S, length = 1000)

# ╔═╡ 4fcd84b0-fe27-11ea-29f2-07ff5129b488
begin
	r = R/100
	z = Z/2
end

# ╔═╡ 9e3552fe-fdee-11ea-1664-97194950ce10
fzero(s) = 2(ray.r(s) - r)*ray.∂r_∂s(s) + 2(ray.z(s) - z)*ray.∂z_∂s(s)

# ╔═╡ d5a98de0-fe27-11ea-0348-0f3eca4ef8ea
rts = roots(fzero, 0..10)

# ╔═╡ Cell order:
# ╟─34c9dff0-fdd8-11ea-383b-b3c3d5e3bbe8
# ╠═1968f0a2-fde4-11ea-3a5e-c1fc5176627f
# ╟─3e7e5a20-fdd9-11ea-3a3b-31fe9237e71e
# ╟─6edbdfb0-fddb-11ea-1b8d-3dd87aa01386
# ╟─d2c2a7f0-fde2-11ea-16ad-ab1114259f56
# ╟─fdf32220-fde1-11ea-25e6-9db949c15ec2
# ╟─75c3a120-fde3-11ea-1682-87e3dc7aee71
# ╠═73d6a8c0-fde4-11ea-237d-15e08ebdaa1c
# ╠═57c41320-fde4-11ea-3fd8-351192415917
# ╠═b1c28750-fde7-11ea-24f0-23050e9d0d0e
# ╠═22d78b60-fde9-11ea-2a88-3118125cf30b
# ╠═ecf09052-fde8-11ea-27cf-c397d368a09c
# ╠═54b58be0-fdea-11ea-356e-b9f30d0f15be
# ╠═9e3552fe-fdee-11ea-1664-97194950ce10
# ╠═4fcd84b0-fe27-11ea-29f2-07ff5129b488
# ╠═d5a98de0-fe27-11ea-0348-0f3eca4ef8ea


text:
struct Header
    sdf_magic::String
    endianness::Int32
    file_version::Int32
    file_revision::Int32
    code_name::String
    first_block_location::Int64
    summary_location::Int64
    summary_size::Int32
    nblocks::Int32
    block_header_length::Int32
    step::Int32
    time::Float64
    jobid1::Int32
    jobid2::Int32
    string_length::Int32
    code_io_version::Int32
end

# T Data type
# N Dimensionality of data
# B Dispatch puropeses (in read for choosing block type)
struct BlockHeader{T,N,B}
    next_block_location::Int64
    data_location::Int64
    id::String
    data_length::Int64
    name::String
end

function BlockHeader(f::IOStream, start, string_length, header_length)
    seek(f, start)
    next_block_location = read(f, Int64)
    data_location = read(f, Int64)
    id = simple_str(read(f, ID_LENGTH))
    data_length = read(f, Int64)
    block_type = read(f, Int32)
    data_type = read(f, Int32)
    n_dims = read(f, Int32)
    name = simple_str(read(f, string_length))

    d_type = typemap(Val(data_type))
    seek(f, start + header_length)

    BlockHeader{d_type, Int(n_dims), block_type}(
        next_block_location,
        data_location,
        id,
        data_length,
        name,
    )
end

abstract type AbstractBlockHeader{T,N} end

Base.ndims(::AbstractBlockHeader{T,N}) where {T,N} = N
Base.eltype(::AbstractBlockHeader{T}) where T = T
Base.nameof(block::BlockHeader) = block.name

struct ConstantBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}
    val::T
end

struct CPUSplitBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}
end

@enum Geometry begin
    Null        = 0
    Cartesian   = 1
    Cylindrical = 2
    Spherical   = 3
end

@doc """
    PlainMeshBlockHeader{T,N}

A mesh defines the locations at which variables are defined.
Since the geometry of a problem is fixed and most variables will be defined
at positions relative to a fixed grid, it makes sense to write this
position data once in its own block.
Each variable will then refer to one of these mesh blocks to provide their location data.

The `PlainMeshBlockHeader` is used for representing the positions at which
scalar field discretizations are defined.
The block header contains the `base_header` (a `BlockHeader`) and the following metadata
- `mults`: The normalisation factor applied to the grid data
in each direction.
- `labels`: The axis labels for this grid in each direction.
- `units`: The units for this grid in each direction after the
normalisation factors have been applied.
- `geometry`: The geometry of the block.
- `minval`: The minimum coordinate values in each direction.
- `maxval`: The maximum coordinate values in each direction.

The geometry of the block can take the following values
- `Null`: Unspecified geometry. This is an error.
- `Cartesian`: Cartesian geometry.
- `Cylindrical`: Cylindrical geometry.
- `Spherical`: Spherical geometry.

The last item in the header is `dims`, is the number of grid points in
each dimension.

The data written is the locations of node points for the mesh in each
of the simulation dimensions.
Therefore for a 3d simulation of resolution ``(nx; ny; nz)``, the data will consist of
a 1d array of X positions with ``(nx + 1)`` elements followed by
a 1d array of Y positions with ``(ny + 1)`` elements and finally
a 1d array of Z positions with ``(nz + 1)`` elements.
Here the resolution specifies the number of simulation cells and therefore
the nodal values have one extra element. In a 1d or 2d simulation, you would
write only the X or X and Y arrays respectively.
"""
struct PlainMeshBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T}

    mults::NTuple{N,Float64}
    labels::NTuple{N,String}
    units::NTuple{N,String}
    geometry::Geometry
    minval::Array{Float64,1}
    maxval::Array{Float64,1}
    dims::Array{Int32,1}
end

@doc """
    PointMeshBlockHeader{T,N}

A mesh defines the locations at which variables are defined.
Since the geometry of a problem is fixed and most variables will be defined
at positions relative to a fixed grid, it makes sense to write this
position data once in its own block.
Each variable will then refer to one of these mesh blocks to provide their location data.

The `PointMeshBlockHeader` is used for representing the positions at which
vector field discretizations are defined.
The block header contains the `base_header` (a `BlockHeader`) and the following metadata
- `mults`: The normalisation factor applied to the grid data
in each direction.
- `labels`: The axis labels for this grid in each direction.
- `units`: The units for this grid in each direction after the
normalisation factors have been applied.
- `geometry`: The geometry of the block.
- `minval`: The minimum coordinate values in each direction.
- `maxval`: The maximum coordinate values in each direction.
- `np`: The number of points.

The geometry of the block can take the following values
- `Null`: Unspecified geometry. This is an error.
- `Cartesian`: Cartesian geometry.
- `Cylindrical`: Cylindrical geometry.
- `Spherical`: Spherical geometry.

The data written is the locations of each point in the first direction followed by
the locations in the second direction and so on.
Thus, for a 3d simulation, if we define the first point as having coordinates
``(x_1; y_1; x_1)`` and the second point as ``(x_2; y_2; z_2)``, etc.
then the data written to file is a 1d array with
 elements ``(x_1; x_2; \\dots; x_{np})``, followed by
the array ``(y_1; y_2; \\dots; y_{np})`` and finally
the array ``(z_1; z_2; \\dots; z_{np})`` where ``np`` corresponds to the
number of points in the mesh. For a 1d simulation, only the x array is written
and for a 2d simulation only the x and y arrays are written.
"""
struct PointMeshBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}

    mults::NTuple{N,Float64}
    labels::NTuple{N,String}
    units::NTuple{N,String}
    geometry::Geometry
    minval::Array{Float64,1}
    maxval::Array{Float64,1}
    np::Int64
end

@enum Stagger begin
    CellCentre = 0
    FaceX      = 1
    FaceY      = 2
    FaceZ      = 4
    EdgeX      = 6
    EdgeY      = 5
    EdgeZ      = 3
    Vertex     = 7
end

@doc """
    PlainVariableBlockHeader{T,N}

The `PlainVariableBlockHeader` is used to describe a variable which is
located relative to the points given in a mesh block.

The block header contains the `base_header` (a `BlockHeader`) and the following metadata
- `mult`: The normalisation factor applied to the variable data.
- `units`: The units for this variable after the normalisation factor has
been applied.
- `mesh_id`: The name(`id`) of the mesh relative to which this block's data is defined.
- `dims`: The number of grid points in each dimension.
- `stagger`: The location of the variable relative to its associated mesh.

The mesh associated with a variable is always node-centred, i.e. the values written as
mesh data specify the nodal values of a grid. Variables may be defined at points
which are offset from this grid due to grid staggering in the code.
The `stagger` entry specifies where the variable is defined relative to the mesh.
Since we have already defined the number of points that the associated mesh contains,
this determines how many points are required to display the variable.

The `stagger` entry can take one of the following values
- `CellCentre`: Cell centred. At the midpoint between nodes. Implies an
``(nx; ny; nz)`` grid.
- `FaceX`: Face centred in X. Located at the midpoint between nodes on
the Y-Z plane. Implies an ``(nx + 1; ny; nz)`` grid.
- `FaceY`: Face centred in Y. Located at the midpoint between nodes on
the X-Z plane. Implies an ``(nx; ny + 1; nz)`` grid.
- `FaceZ`: Face centred in Z. Located at the midpoint between nodes on
the X-Y plane. Implies an ``(nx; ny; nz + 1)`` grid.
- `EdgeX`: Edge centred along X. Located at the midpoint between nodes
along the X-axis. Implies an ``(nx; ny + 1; nz + 1)`` grid.
- `EdgeY`: Edge centred along Y. Located at the midpoint between nodes
along the Y-axis. Implies an ``(nx + 1; ny; nz + 1)`` grid.
- `EdgeZ`: Edge centred along Z. Located at the midpoint between nodes
along the Z-axis. Implies an ``(nx + 1; ny + 1; nz)`` grid.
- `Vertex`: Node centred. At the same place as the mesh. Implies an ``(nx+
1; ny + 1; nz + 1)`` grid.

For a grid based variable, the data written contains the values of
the given variable at each point on the mesh.
This is in the form of a 1d, 2d or 3d array depending on the dimensions
of the simulation. The size of the array depends on the size of
the associated mesh and the grid staggering as indicated above. It
corresponds to the values written into the `dims` array written for
this block.
"""
struct PlainVariableBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}

    mult::Float64
    units::String
    mesh_id::String
    dims::NTuple{N,Int32}
    stagger::Stagger
end

@doc """
    PointVariableBlockHeader{T,N}

The `PointVariableBlockHeader` is used to describe a variable which is
located relative to the points given in a mesh block.

The block header contains the `base_header` (a `BlockHeader`) and the following metadata
- `mult`: The normalisation factor applied to the variable data.
- `units`: The units for this variable after the normalisation factor has
been applied.
- `mesh_id`: The name(`id`) of the mesh relative to which this block's data is defined.
- `np`: The number of points.

Similarly to the grid based variable, the data written contains the values
of the given variable at each point on the mesh. Since each the location of
each point in space is known fully, there is no need for a stagger variable.
The data is in the form of a 1d array with `np` elements.
"""
struct PointVariableBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}

    mult::Float64
    units::String
    mesh_id::String
    np::Int64
end

struct RunInfoBlockHeader{T,N} <: AbstractBlockHeader{T,N}
    base_header::BlockHeader{T,N}
end


text:
function process_momentum_runs(mechanism::Mechanism, storage::Vector{Storage{T,N}}, tcompute::Vector) where {T,N}
	H = mechanism.timestep * N
	Nsim = length(tcompute)
	speed = 0.0
	plin = 0.0
	pang = 0.0
	for i = 1:Nsim
		p = momentum(mechanism, storage[i], N)
		speed += H / tcompute[i]
		plin += norm(p[1:3])
		pang += norm(p[4:6])
	end
	speed /= Nsim
	plin /= Nsim
	pang /= Nsim
	return speed, plin, pang
end

function process_energy_runs(mechanism::Mechanism, storage::Vector{Storage{T,N}},
		tcompute::Vector, tctrl::T) where {T,N}
	H = mechanism.timestep * N
	Nsim = length(tcompute)
	speed = 0.0
	energy = 0.0
	for i = 1:Nsim
		speed += H / tcompute[i]
		initial = Int(floor(tctrl/mechanism.timestep+2))
		final = N
		energy += kinetic_energy(mechanism, storage[i], final) - kinetic_energy(mechanism, storage[i], initial)
	end
	speed /= Nsim
	energy /= Nsim
	return speed, energy
end


text:
using ProgressMeter
using FFTW, LinearAlgebra
using Plots
using JLD

abstract type AbstractModel end
abstract type TimeSolver    end
abstract type InitialData   end

struct Mesh

    N    :: Int64
    xmin :: Float64
    xmax :: Float64
    dx   :: Float64
    x    :: Vector{Float64}
    kmin :: Float64
    kmax :: Float64
    dk   :: Float64
    k    :: Vector{Float64}

    function Mesh(param :: NamedTuple)

        xmin = - Float64(param.L)
        xmax =   Float64(param.L)
        N    =   param.N

        dx   = (xmax-xmin)/N
        x    = zeros(Float64, N)
        x   .= range(xmin, stop=xmax, length=N+1)[1:end-1]
        dk   = 2π/(N*dx)
        kmin = -N/2*dk
        kmax = (N/2-1)*dk
        k    = zeros(Float64, N)
        k   .= dk .* vcat(0:N÷2-1, -N÷2:-1)

        new( N, xmin, xmax, dx, x, kmin, kmax, dk, k)

    end

end

struct Times

    Nt   :: Int
    Nr   :: Int
    nr   :: Int
    tfin :: Float64
    dt   :: Float64
    t    :: Vector{Float64}
    tr   :: Vector{Float64}

    function Times( dt , tfin)
        t  = 0:dt:tfin
        Nt = length(t)
        Nr = Nt
        nr = 1
        tr = t
        new( Nt, Nr, nr, tfin, dt, t, tr)
    end

end

"""
    RK4(params)

Runge-Kutta fourth order solver.

"""
mutable struct RK4 <: TimeSolver

    Uhat :: Array{ComplexF64,2}
    dU   :: Array{ComplexF64,2}

    function RK4( param::NamedTuple, model::AbstractModel )

        n = param.N

        Uhat = zeros(ComplexF64, (n,2))
        dU   = zeros(ComplexF64, (n,2))

        new( Uhat, dU)

    end

end

function step!(s  :: RK4,
               f! :: AbstractModel,
               U  :: Array{ComplexF64,2},
               dt :: Float64)

    @inbounds for i in eachindex(U)
        s.Uhat[i] = U[i]
    end

    f!( s.Uhat )

    @inbounds for i in eachindex(U)
        s.dU[i]   = s.Uhat[i]
        s.Uhat[i] = U[i] + dt/2 * s.Uhat[i]
    end

    f!( s.Uhat )

    @inbounds for i in eachindex(U)
        s.dU[i]   += 2 * s.Uhat[i]
        s.Uhat[i]  = U[i] + dt/2 * s.Uhat[i]
    end

    f!( s.Uhat )

    @inbounds for i in eachindex(U)
        s.dU[i]   += 2 * s.Uhat[i]
        s.Uhat[i]  = U[i] + dt * s.Uhat[i]
    end

    f!( s.Uhat )

    @inbounds for i in eachindex(U)
        s.dU[i] += s.Uhat[i]
        U[i]    += dt/6 * s.dU[i]
    end

end

struct BellCurve <: InitialData

    h :: Vector{ComplexF64}
    u :: Vector{ComplexF64}

    function BellCurve(p :: NamedTuple,theta :: Real)

        mesh  = Mesh(p)
        h     = zeros(ComplexF64, mesh.N)
        h    .= exp.(.-((abs.(mesh.x)).^theta).*log(2))
        u     = zeros(ComplexF64, mesh.N)

        new( h, u )

    end

end

"""
    Matsuno(params)

"""
mutable struct Matsuno <: AbstractModel

    label    :: String
    datasize :: Int
    Γ        :: Array{Float64,1}
    Dx       :: Array{ComplexF64,1}
    H        :: Array{ComplexF64,1}
    Π⅔       :: BitArray{1}
    ϵ        :: Float64
    hnew     :: Vector{ComplexF64}
    unew     :: Vector{ComplexF64}
    I₀       :: Vector{ComplexF64}
    I₁       :: Vector{ComplexF64}
    I₂       :: Vector{ComplexF64}
    I₃       :: Vector{ComplexF64}

    Px       :: FFTW.FFTWPlan

    function Matsuno(param::NamedTuple)

        label    = "Matsuno"
        datasize = 2
        ϵ        = param.ϵ
        mesh     = Mesh(param)
        Γ        = abs.(mesh.k)
        Dx       =  1im * mesh.k        # Differentiation
        H        = -1im * sign.(mesh.k) # Hilbert transform
        Π⅔       = Γ .< mesh.kmax * 2/3 # Dealiasing low-pass filter

        hnew = zeros(ComplexF64, mesh.N)
        unew = zeros(ComplexF64, mesh.N)

        I₀ = zeros(ComplexF64, mesh.N)
        I₁ = zeros(ComplexF64, mesh.N)
        I₂ = zeros(ComplexF64, mesh.N)
        I₃ = zeros(ComplexF64, mesh.N)

        Px  = plan_fft(hnew; flags = FFTW.MEASURE)

        new(label, datasize, Γ, Dx, H, Π⅔, ϵ,
            hnew, unew, I₀, I₁, I₂, I₃, Px)
    end
end

function (m::Matsuno)(U::Array{ComplexF64,2})

    @inbounds for i in eachindex(m.hnew)
        m.hnew[i] = m.Γ[i] * U[i,1]
    end

    ldiv!(m.unew, m.Px, m.hnew )

    @inbounds for i in eachindex(m.hnew)
        m.hnew[i] = m.Dx[i] * U[i,1]
    end

    ldiv!(m.I₁, m.Px, m.hnew)

    @inbounds for i in eachindex(m.unew)
        m.unew[i] *= m.I₁[i]
    end

    mul!(m.I₁, m.Px, m.unew)

    @inbounds for i in eachindex(m.hnew)
        m.I₁[i] = m.I₁[i] * m.ϵ * m.Π⅔[i] - m.hnew[i]
    end

    ldiv!(m.hnew, m.Px, view(U,:,1))
    ldiv!(m.unew, m.Px, view(U,:,2))

    @inbounds for i in eachindex(m.hnew)
        m.I₂[i] = m.hnew[i] * m.unew[i]
    end

    mul!(m.I₃, m.Px, m.I₂)

    @inbounds for i in eachindex(m.H)
        U[i,1]  = m.H[i] * U[i,2]
        m.I₀[i] = m.Γ[i] * U[i,2]
    end

    ldiv!(m.I₂, m.Px, m.I₀)

    @inbounds for i in eachindex(m.hnew)
        m.I₂[i] *= m.hnew[i]
    end

    mul!(m.hnew, m.Px, m.I₂)

    @inbounds for i in eachindex(m.unew)
        U[i,1] -= (m.I₃[i] * m.Dx[i] + m.hnew[i] * m.H[i]) * m.ϵ * m.Π⅔[i]
        m.I₃[i]  = m.unew[i]^2
    end

    mul!(m.unew, m.Px, m.I₃)

    @inbounds for i in eachindex(m.unew)
        U[i,2]  =  m.I₁[i] - m.unew[i] * m.Dx[i] * m.ϵ/2 * m.Π⅔[i]
    end

end

function create_animation( mesh, times, data )

    prog = Progress(times.Nr,1)

    hr = zeros(Float64,mesh.N)
    ur = zeros(Float64,mesh.N)

    anim = @animate for l in 1:size(data)[end]

        pl = plot(layout=(2,1))

        hr .= real(ifft(view(data,:,1,l)))
        ur .= real(ifft(view(data,:,2,l)))

        plot!(pl[1,1], mesh.x, hr;
              ylims=(-0.6,1),
              title="physical space",
              label=p.model.label)

        plot!(pl[2,1], fftshift(mesh.k),
              log10.(1e-18.+abs.(fftshift(fft(hr))));
              title="frequency",
              label=model.label)

        next!(prog)

    end when mod(l, 100) == 0

    gif(anim, "anim.gif", fps=15); nothing

end

function main()

    param = ( ϵ  = 1/2,
              N  = 2^12,
              L  = 10.,
              T  = 5.,
              dt = 0.001 )

    mesh    = Mesh(param)
    times   = Times(param.dt, param.T)
    init    = BellCurve(param,2.5)
    model   = Matsuno(param)
    solver  = RK4(param,model)

    U       = zeros(ComplexF64,(mesh.N,2))
    U[:,1] .= model.Π⅔ .* fft(init.h)
    U[:,2] .= model.Π⅔ .* fft(init.u)

    dt = times.dt
    nr = times.nr

    J = 1:times.Nt÷nr
    L = 1:nr
    data = zeros(ComplexF64,(mesh.N,2,length(J)))

    @showprogress 1 for j in J
        for l in L
            step!(solver, model, U, dt)
        end
        data[:,:,j] = U
    end

    print("\n")

    data

    #Uref =  problem.data[end]
    #save("reference.jld", "Uref", Uref)

    #Uref = load("reference.jld", "Uref")

    #println(norm(Uref .- data[end]))

    # @time create_animation( problem, data )

end

@time main()


text:
using Documenter:
    DocMeta,
    HTML,
    asset,
    deploydocs,
    makedocs
using PowerAnalyses

DocMeta.setdocmeta!(
    PowerAnalyses,
    :DocTestSetup,
    :(using PowerAnalyses);
    recursive=true
)

sitename = "PowerAnalyses.jl"
pages = [
    "PowerAnalyses" => "index.md"
]

prettyurls = get(ENV, "CI", nothing) == "true"
format = HTML(; prettyurls)
modules = [PowerAnalyses]
strict = true
checkdocs = :none
makedocs(; sitename, pages, format, modules, strict, checkdocs)

deploydocs(;
    branch="docs",
    devbranch="main",
    repo="github.com/rikhuijzer/PowerAnalyses.jl.git",
    push_preview=false
)


text:
# # Helper functions

Point2{T}(a::Vector{T}) where T = Point2(a[1], a[2])

"""
    toPointVecs(poly::Polygon)

Creates a vector of [`Point`](@ref)s from the given polygon.
Returns a `Vector{Vector{Point}}`, with one element.
"""
toPointVecs(poly::T) where T <: GeoInterface.AbstractPolygon = [Point2.(cs) for cs in GeoInterface.coordinates(poly)]

"""
    toPointVecs(mp::MultiPolygon)

Creates a vector of vectors of [`Point`](@ref)s from the given MultiPolygon.
Returns a `Vector{Vector{Point}}`.

!!! warning
    There is no support for holes in this function at present; that will need to
    wait for GeometryBasics.jl to become complete, possibly replacing GeoInterface.
"""
function toPointVecs(mp::T) where T <: GeoInterface.AbstractMultiPolygon

    polys = GeoInterface.coordinates(mp)

    return to_nansep_vec(identity, toPointVecs.(polys))

end

"""
    nonzero(pt::Point{2, T}) where T

Returns `true` if both of the elements are nonzero, and `false` otherwise.
"""
nonzero(pt::Point{2, T}) where {T} = pt[1] != zero(T) && pt[2] != zero(T)

"""
    imprecise(arr)

Returns `true` if the polygon is invalid when converted to Float32,
and `false` otherwise.
"""
function imprecise(arr)
    dif = Point2{Float32}.(arr) |> diff
    return !any(nonzero.(dif))
end

"""
    toMeshes(mp::Vector{Vector{Vector{Point2}}})::Vector{GLNormalUVMesh}
    toMeshes(p::Vector{Vector{Point2}})::Vector{GLNormalUVMesh}

Takes in a polygon, or a multipolygon, and returns a vector of meshes
(which may only contain one mesh).
"""
function toMeshes(mp::Vector{Vector{Vector{Point2{T}}}}) where T

    meshes = GLNormalMesh[]
    # give a hint for the size of the vector so we don't dynamically allocate
    sizehint!(meshes, length(mp))

    for pol in mp
        # this **should** account for holes...
        triangle_faces = EarCut.triangulate(pol)

        v = map(x-> Point3f0(x[1], x[2], 0), vcat(pol...))

        push!(meshes, GeometryBasics.Mesh(v, triangle_faces))
    end

    return meshes
end

function toMeshes(mp::Vector{Vector{Point2{T}}}) where T

    triangle_faces = EarCut.triangulate(mp)

    v = map(x-> Point3{T}(x[1], x[2], 0), vcat(mp...))

    return [GeometryBasics.Mesh(v, triangle_faces)]
end

function toMeshes(poly::T) where T <: GeoInterface.AbstractPolygon

    mp = rec_point(GeoInterface.coordinates(poly))

    triangle_faces = EarCut.triangulate(mp)

    v = vcat(mp...)

    return [Mesh(v, triangle_faces)]
end

function toMeshes(mps::Vector{<: GeoInterface.AbstractMultiPolygon})
    meshes = GeometryBasics.Mesh[]
    sizehint!(meshes, length(mps))
    for mp in mps
        polys = GeoInterface.coordinates(mp) |> rec_point
        push!(meshes, merge(vcat(toMeshes.(polys)...)))
    end

    return meshes
end

function toMeshes(mp::T where T <: GeoInterface.AbstractMultiPolygon)
    meshes = GeometryBasics.Mesh[]

    for p in GeoInterface.coordinates(mp)
        poly = rec_point.(p)
        append!(meshes, toMeshes(poly))
    end

    return meshes
end
# # Argument conversions

# ## Polygons

convert_arguments(::Type{<: Poly}, poly::T) where T <: GeoInterface.AbstractPolygon = (toPointVecs(poly)[1],)

convert_arguments(::Type{<: Poly}, mp::T) where T <: GeoInterface.AbstractMultiPolygon = (toPointVecs(mp),)

# Only converts polygons and multipolygons
function convert_arguments(::Type{<: AbstractPlotting.Poly}, fc::GeoInterface.AbstractFeatureCollection)
    features = GeoInterface.features(fc)

    cs = Vector{Point2f0}[]

    imp = false

    for feature in features
        for coord in toPointVecs(feature.geometry)
            if imprecise(coord) || length(coord) < 4
                @warn(
                    """
                    Imprecise coordinates!  You may want to consider reprojecting to a different coordinate system.
                    """,
                    coord
                )
                imp = true
            else
                push!(cs, coord)
            end
        end

    end

    imp && @warn(
                """Coordinates have been removed, because they were
                too ambiguous.  Number of coordinates is now $(length(cs)).
                """
            )

    (cs,)

end

function AbstractPlotting.convert_arguments(::Type{<: Mesh}, fc::GeoInterface.FeatureCollection)
    return (geo |> GeoInterface.features .|> GeoInterface.geometry .|> toMeshes,) # return a Vector of meshes
end

convert_arguments(::Type{<: Mesh}, polys::Vector{<:GeoInterface.AbstractMultiPolygon}) = (toMeshes(polys),)

# set the default plot type for Vectors of polygons,
# so that they are plotted using the most efficient method!
plottype(::Vector{<: GeoInterface.AbstractMultiPolygon}) = Mesh
plottype(::Vector{<: GeoInterface.AbstractPolygon}) = Mesh

function Proj4.transform(src, dest, pt::Point2{T}) where T
    if isnan(pt[1]) && isnan(pt[2])
        return Point2{T}(NaN)
    end
    return Point2{T}(transform(src, dest, Vector(pt)))
end
Base.Broadcast.broadcastable(p::Projection) = Ref(p)


text:
module SModels

  # dependencies
  using DecisionTree
  using ScikitLearn
  using ScikitLearn.GridSearch: GridSearchCV
  using ScikitLearn.CrossValidation: cross_val_predict
  using PyCall
  using PyCallJLD
  using JLD
  using DistributedArrays
  using Distributions

  # imports from ScikitLearn
  @sk_import neural_network: MLPRegressor  #to train the regressor
  @sk_import neural_network: MLPClassifier #to train the classifier (convergence or not of the model)
  @sk_import preprocessing: StandardScaler #to scale observations
  @sk_import metrics: accuracy_score       #accuracy of classifier

  # package code goes here
  include("types.jl")
  include("generic.jl")
  include("draws.jl")
  include("train.jl")

  # types.jl
  #---------
  export SModelsProblem, SModelsOptions, default_function

  #generic.jl
  #-----------
  export split_train_test, convert_to_array_of_arrays, convert_to_array
  export predict_scaled_robust
  export convert_vector_to_array, convert_vector_to_array, predict_one_obs_scaled
  export predict_one_obs_unscaled, predict_one_obs_unscaled, set_model_function!
  export calculate_mean_per_error, calculate_median_abs_per_error
  export calculate_maximum_abs_per_error
  export set_bounds!
  export date_now

  # draws.jl
  #---------
  export create_grid_stochastic, generate_std
  # train.jl
  #---------
  export train_surrogate_model
  export evaluateModel!

end # module


text:
function apply_dyson!(H, psi, dt::Real, dysonOrder::Integer)
    kq = copy(psi)
    tmp = similar(kq)
    for k = 1 : dysonOrder
        mul!(tmp, H, kq, -1im * dt/k, 0)
        psi .+= tmp
        copyto!(kq, tmp)
    end
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
export PetscMat, MatSetType, SetUp, MatSetValues, MatAssemblyBegin, MatAssemblyEnd, MatSetSizes, MatGetSize, MatGetValues, MatGetOwnershipRange, MatXAIJSetPreallocation, MatMPIAIJSetPreallocation, SetFromOptions, MatGetInfo, MatMatMult, MatNorm, MatZeroEntries, MatSetValuesBlocked, MatSetOption, MatCreateShell, MatShellSetOperation, MatShellGetContext, MatShellSetContext, MatGetType, MatCreateTranspose, MatTranspose

"""
  PetscMat type.  Currently a subtype of `AbstractArray`, although that may
  change.
"""
type PetscMat  <: AbstractArray{PetscScalar, 2}
  pobj::Ptr{Void}

  """
    Constructor

    **Input**

     * comm: MPI communicator
  """
  function PetscMat(comm::MPI_Comm)
#    comm = PETSC_COMM_SELF();
    vec = Array{Ptr{Void}}(1)
    err = ccall( (:MatCreate,  libpetsclocation), PetscErrorCode, (comm_type, Ptr{Ptr{Void}}),comm,vec);
    vec = new(vec[1])
#    finalizer(vec,PetscDestroy)
    # does not seem to be called immediately when vec is no longer visible, is it called later during garbage collection?
    return vec
  end

  """
    Constructs a PetscMat object from an existing pointer to a Petsc matrix.

    **Inputs**

     * pobj: a void pointer to a Petsc matrix
  """
  function PetscMat(pobj::Ptr{Void})  # default constructor
    return new(pobj)
  end
end

"""
  MatCreateShell
"""
function MatCreateShell(arg1::MPI_Comm,arg2::Integer,arg3::Integer,arg4::Integer,arg5::Integer, arg6::Ptr{Void})
  # arg6 is the user provided context
    arg7 = Ref{Ptr{Void}}()
    ccall((:MatCreateShell,petsc),PetscErrorCode,(comm_type,PetscInt,PetscInt,PetscInt,PetscInt,Ptr{Void}, Ref{Ptr{Void}}),arg1,arg2,arg3,arg4,arg5,arg6,arg7)
    return PetscMat(arg7[])
end

"""
  MatShellSetOperation
"""
function MatShellSetOperation(arg1::PetscMat,arg2::MatOperation,arg3::Ptr{Void})
# arg3 is a function pointer, and must have the signature:
# void fname(Mat, vec, vec) for MATOP_MULT

    ccall((:MatShellSetOperation,petsc),PetscErrorCode,(Ptr{Void},MatOperation,Ptr{Void}),arg1.pobj,arg2,arg3)
end

"""
  MatShellGetContext
"""
function MatShellGetContext(arg1::PetscMat)
# get the user provided context for the matrix shell
    arg2 = Ref{Ptr{Void}}()
    ccall((:MatShellGetContext,petsc),PetscErrorCode,(Ptr{Void},Ref{Ptr{Void}}),arg1.pobj,arg2)
    return arg2[]  # turn it into a julia object here?
end

function MatShellSetContext(arg1::PetscMat, arg2::Ptr{Void})
# get the user provided context for the matrix shell
    ccall((:MatShellSetContext,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void}),arg1.pobj,arg2)
    return nothing  # turn it into a julia object here?
end

"""
  Frees a Petsc object.  Safe to call multiple times
"""
function PetscDestroy(vec::PetscMat)
  if (vec.pobj != 0)
    err = ccall( (:MatDestroy,  libpetsclocation), PetscErrorCode, (Ptr{Ptr{Void}},), &vec.pobj);
  end
  vec.pobj = 0

#    println("Petsc Mat Destroy called")
#    sleep(5)
  return 0
end

"""
  Equivalent to Petsc's MatInfo struct
"""
immutable MatInfo
    block_size::PetscLogDouble
    nz_allocated::PetscLogDouble
    nz_used::PetscLogDouble
    nz_unneeded::PetscLogDouble
    memory::PetscLogDouble
    assemblies::PetscLogDouble
    mallocs::PetscLogDouble
    fill_ratio_given::PetscLogDouble
    fill_ratio_needed::PetscLogDouble
    factor_mallocs::PetscLogDouble

    function MatInfo()  # incomplete initialization
      return new()
    end
end

import Base.show
"""
  show() for a MatInfo
"""
function show(io::IO, obj::MatInfo)
# print the fields of PetscMatInfo
#  println("PetscMatInfo:")
  println(io, "  block_size : ", obj.block_size)
  println(io, "  nz_allocated : ", obj.nz_allocated)
  println(io, "  nz_used : ", obj.nz_used)
  println(io, "  nz_unneeded : ", obj.nz_unneeded)
  println(io, "  memory : ", obj.memory)
  println(io, "  assemblies : ", obj.assemblies)
  println(io, "  mallocs : ", obj.mallocs)
  println(io, "  fill_ratio_given : ", obj.fill_ratio_given)
  println(io, "  fill_ratio_needed : ", obj.fill_ratio_needed)
  println(io, "  factor_mallocs : ", obj.factor_mallocs)
end

"""
  MatSetFromOptions
"""
function SetFromOptions(mat::PetscMat)
  ccall((:MatSetFromOptions,petsc),PetscErrorCode,(Ptr{Void},), mat.pobj)

end

"""
  MatSetType
"""
function MatSetType(vec::PetscMat,name)
  err = ccall( (:MatSetType,  libpetsclocation), PetscErrorCode,(Ptr{Void}, Cstring), vec.pobj,name);

end

"""
  MatSetUp
"""
function SetUp(vec::PetscMat)
  err = ccall( ( :MatSetUp,  libpetsclocation), PetscErrorCode, (Ptr{Void},), vec.pobj);

end

"""
  MatGetType
"""
function MatGetType(arg1::PetscMat)
    arg2 = Ref{Ptr{UInt8}}()
    ccall((:MatGetType,petsc),PetscErrorCode,(Ptr{Void}, Ref{Ptr{UInt8}}),arg1.pobj,arg2)
    return unsafe_string(arg2[])
end

#=
  MAT_FLUSH_ASSEMBLY = 1;
  MAT_FINAL_ASSEMBLY = 0
=#

"""
  MatSetValues.  The `idx` and `idy` arrays must have PetscInt elements, and the
  `vals` array must have PetscScalar elements.
"""
function MatSetValues(vec::PetscMat,idi::Array{PetscInt},idj::Array{PetscInt},array::Array{PetscScalar},flag::Integer)
  # remember, only matrices can be inserted into a Petsc matrix
  # if array is a 3 by 3, then idi and idj are vectors of length 3
#    idi = idi
#    idj = idj

  @assert length(idi)*length(idj) == length(array)

  # do check here to ensure array is the right shape (remember tranpose)
  err = ccall( ( :MatSetValues,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar},Int32), vec.pobj,length(idi), idi, length(idj), idj,array,flag);
#    idi = idi
#    idj = idj
  return err
end

"""
  MatSetValues blocked.
"""
function MatSetValuesBlocked(mat::PetscMat, idi::Array{PetscInt}, idj::Array{PetscInt}, v::Array{PetscScalar}, flag::Integer)

    err = ccall((:MatSetValuesBlocked, libpetsclocation), PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}, Int32),mat.pobj, length(idi), idi, length(idj), idj, v,flag)

    return err
end

"""
  MatSetOption
"""
function MatSetOption(mat::PetscMat,arg2::MatOption,arg3::Union{Integer, Bool})
    ierr = ccall((:MatSetOption, libpetsclocation),PetscErrorCode, (Ptr{Void},MatOption,PetscBool), mat.pobj, arg2, arg3)

    if ierr != 0
      println(STDERR, "Error: MatSetOption returned non zero exit status")
    end
end

"""
  MatAssemblyBegin
"""
function MatAssemblyBegin(obj::PetscMat,flg::Integer)
  err = ccall( ( :MatAssemblyBegin,  libpetsclocation), PetscErrorCode,(Ptr{Void},Int32), obj.pobj,flg);
end

"""
  MatAssemblyBegin, impicitly using `MAT_FINAL_ASSEMBLY`
"""
function MatAssemblyBegin(obj::PetscMat)
  return MatAssemblyBegin(obj,MAT_FINAL_ASSEMBLY);
end

"""
  MatAssemblyEnd
"""
function MatAssemblyEnd(obj::PetscMat,flg::Integer)
  err = ccall( ( :MatAssemblyEnd,  libpetsclocation), PetscErrorCode,(Ptr{Void},Int32), obj.pobj,flg);
end

"""
  MatAssemblyEnd, implicitly using `MAT_FINAL_ASSEMBLY`
"""
function MatAssemblyEnd(obj::PetscMat)
  return MatAssemblyEnd(obj,MAT_FINAL_ASSEMBLY);
end

"""
  MatSetSizes
"""
function MatSetSizes(vec::PetscMat,m::Integer, n::Integer, M::Integer, N::Integer)
  err = ccall( ( :MatSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void},PetscInt, PetscInt, PetscInt, PetscInt), vec.pobj,m,n,M,N);

  return err
end

"""
  Constructs the transpose of a given matrix.  Can be in place or out of place.

  **Inputs**

   * A: matrix to take the transpose of

  **Keywords**

   * inplace: whether or not to do the transpose in place

  **Outputs**

   * A PetscMat object, either the original object A if the transpose was done
     in place, or a new matrix object if the transpose was done out of place
"""
function MatTranspose(A::PetscMat; inplace::Bool=false, mat_initialized=false)
  if inplace
    B = Ref{Ptr{Void}}(A.pobj)
    reuse = MAT_INPLACE_MATRIX
  else
    B = Ref{Ptr{Void}}(C_NULL)
    reuse = MAT_INITIAL_MATRIX
  end

    ccall((:MatTranspose,petsc),PetscErrorCode,(Ptr{Void},MatReuse,Ptr{Ptr{Void}}),A.pobj, reuse, B)

  if inplace
    return A
  else
    return PetscMat(B[])
  end
end

"""
  MatCreateTranspose
"""
function MatCreateTranspose(arg1::PetscMat)
  arg2 = Ref{Ptr{Void}}()

  ccall((:MatCreateTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Ptr{Void}}),arg1.pobj, arg2)

  return PetscMat(arg2[])
end

"""
  PetscViewer
"""
function PetscView(obj::PetscMat,viewer::PetscViewer=C_NULL)
  err = ccall( (:MatView,  libpetsclocation), PetscErrorCode, (Ptr{Void}, Ptr{Void}),obj.pobj,viewer);
end

"""
  MatGetSize

  **Inputs**

   * obj: a PetscMat

  **Outputs**

   * m: first dimension size
   * n: second dimension size
"""
function MatGetSize(obj::PetscMat)
  m = Array{PetscInt}(1)
  n = Array{PetscInt}(1)
  err = ccall(Libdl.dlsym(libpetsc, :MatGetSize), PetscErrorCode,(Ptr{Void}, Ptr{PetscInt},Ptr{PetscInt}), obj.pobj,m,n);
  return m[1],n[1]
end

export MatGetLocalSize

"""
  MatGetLocalsize.  Same interface as `MatGetSize`
"""
function MatGetLocalSize(mat::PetscMat)
    m = Array{PetscInt}(1)
    n = Array{PetscInt}(1)
    ccall((:MatGetLocalSize,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),mat.pobj, m, n)
    return m[1], n[1]
end

# Petsc populates v row major, so Julia will see the returned array as being transposed
"""
  MatGetValues

  **Inputs**

   * obj: the PetscMat
   * idxm: row indices
   * idxn: column indices
   * vals: logically 2 dimensional array of values (although it could be a
           vector too)
"""
function MatGetValues(obj::PetscMat, idxm::Array{PetscInt, 1}, idxn::Array{PetscInt, 1}, v::Array{PetscScalar, 2})
    # do check here to ensure v is the right shape

    ccall((:MatGetValues,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}), obj.pobj, length(idxm), idxm, length(idxn), idxn, v)

end

"""
  MatGetOwnershipRange

  **Inputs**

   * mat: a PetscMat

  **Outputs**

   * low: lowest (zero-based) index that is owned
   * high: highest + 1 (zero-based) index that is owned
"""
function MatGetOwnershipRange(mat::PetscMat)
    low = Array{PetscInt}(1)
    high = Array{PetscInt}(1)
    ccall((:MatGetOwnershipRange,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),mat.pobj, low, high)

    return low[1], high[1]
end

"""
  MatZeroEntries
"""
function MatZeroEntries(mat::PetscMat)
    ccall((:MatZeroEntries,petsc),PetscErrorCode,(Ptr{Void},),mat.pobj)
end

### new function ###

"""
  MatNorm
"""
function MatNorm(mat::PetscMat, ntype::NormType)
    nrm = Array{PetscReal}(1)
    ccall((:MatNorm,petsc),PetscErrorCode,(Ptr{Void}, NormType,Ptr{PetscReal}), mat.pobj, ntype, nrm)
    return nrm[1]
end

export MatAXPY, MatAYPX, MatScale, MatShift

"""
  MatAxPY
"""
function MatAXPY(Y::PetscMat, a::PetscScalar, X::PetscMat, str::PetscMatStructure)
    ccall((:MatAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void},PetscMatStructure), Y.pobj, a, X.pobj, str)
end

"""
  MatAYPX
"""
function MatAYPX( Y::PetscMat, a::PetscScalar, X::PetscMat, str::PetscMatStructure)
    ccall((:MatAYPX,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}, PetscMatStructure), Y.pobj, a, X.pobj, str)
end

"""
  MatScale
"""
function MatScale(mat::PetscMat, a::Number)
    ccall((:MatScale,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), mat.pobj, PetscScalar(a))
end

"""
  MatShift
"""
function MatShift(mat::PetscMat, a::Number)
    ccall((:MatShift,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), mat.pobj, a)
end

export MatMult, MatMultAdd, MatMultTranspose, MatMultHermitianTranspose

"""
  MatMult
"""
function MatMult(mat::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMult,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), mat.pobj, x.pobj, y.pobj)
end

"""
  MatMultAdd
"""
function MatMultAdd(A::PetscMat, x::PetscVec, y::PetscVec, z::PetscVec)
    ccall((:MatMultAdd,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}), A.pobj, x.pobj, y.pobj, z.pobj)
end

"""
  MatMultTranspose
"""
function MatMultTranspose(A::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMultTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), A.pobj, x.pobj, y.pobj)
end

"""
  MatMultHermitianTranspose
"""
function MatMultHermitianTranspose(mat::PetscMat, x::PetscVec, y::PetscVec)
    ccall((:MatMultHermitianTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}),mat.pobj, x.pobj, y.pobj)
end

"""
  MatMatMult
"""
function MatMatMult(A::PetscMat, B::PetscMat, scall::MatReuse, fill::PetscReal, C::PetscMat)
    arr = [C.pobj]
    ccall((:MatMatMult,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, MatReuse, PetscReal, Ptr{Ptr{Void}}), A.pobj, B.pobj, scall, fill, arr)
    if C.pobj != arr[1]
      PetscDestroy(C)
      C.pobj = arr[1]
    end

end

"""
  MatXAIJSetPreallocation
"""
function MatXAIJSetPreallocation(mat::PetscMat, bs::Integer, dnnz::AbstractArray{PetscInt, 1}, onnz::AbstractArray{PetscInt,1}, dnnzu::AbstractArray{PetscInt, 1}, onnzu::AbstractArray{PetscInt, 1})
# this is a unified interface for matrix preallocation for the Petsc built in
# matrix types: Aij, Bij, and their respective symmetric forms SAij, SBij
# for a non symmetric format matrix, dnnzu and onnzu are not required
# and vice versa for a non symmetric format matrix

    ccall((:MatXAIJSetPreallocation,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},Ptr{PetscInt},Ptr{PetscInt},Ptr{PetscInt}), mat.pobj, bs, dnnz, onnz, dnnzu, onnzu)
end

# Matrix preallocation
"""
  MatMPIAIJSetPreallocation
"""
function MatMPIAIJSetPreallocation(mat::PetscMat, d_nz::Integer, d_nnz::PetscInt_arr_or_null, o_nz::Integer, o_nnz::PetscInt_arr_or_null)

    ccall((:MatMPIAIJSetPreallocation,petsc),PetscErrorCode,(Ptr{Void}, PetscInt,Ptr{PetscInt},PetscInt,Ptr{PetscInt}), mat.pobj, d_nz, d_nnz, o_nz, o_nnz)
end

"""
  MatGetInfo
"""
function MatGetInfo(mat::PetscMat, info_type::Int32)
    matinfo = MatInfo()  # create uninitialized struct
    ref_matinfo = Ref{MatInfo}(matinfo)
    ccall((:MatGetInfo,petsc),PetscErrorCode,(Ptr{Void}, Int32,Ref{MatInfo}), mat.pobj ,info_type, ref_matinfo)
    return ref_matinfo[]
end

# some functions to make a PetscMat look like a Julia array
# these only work after a matrix has been assembled
# some of them only work for uniprocess matrices
import Base.size
import Base.getindex

function size(A::PetscMat)

  return MatGetLocalSize(A)
end

function size(A::PetscMat, dim::Integer)
  dims = MatGetLocalSize(A)
  return dims[dim]
end

function getindex(A::PetscMat, i::Integer, j::Integer)
  # not efficient

  i_ = [PetscInt(i - 1)]
  j_ = [PetscInt(j - 1)]
  val = Array{PetscScalar}(1, 1)

  MatGetValues(A, i_, j_, val)

  return val[1]
end

function getindex(A::PetscMat, idx::Integer)
  # linear indexing
  m,n = MatGetLocalSize(A)
  i = idx % m
  j = div(idx, m)

  # zero based indexing
  i -= 1
  j -= 1

  i_ = [PetscInt(i)]
  j_ = [PetscInt(j)]
  val = Array{PetscScalar}(1, 1)

  MatGetValues(A, i_, j_, val)
  return val[1]
end

include("mat_interface.jl")


text:
function f(n)
    for n = 1:n
        println(n)
    end
end

f(5)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
    pass3(project)

 - Remove `pass3/`
 - Copy files from `pass2/` into `pass3`
 - Translating scope syntax
"""
function pass3(project)
    println("Performing pass 3...")
    files = prep_pass(project, 3)

    for f in files

        println("   -----------------------")
        println("   Removing comments for file $(basename(file))")
        contents = open(f->read(f, String), file)

        println("   Single-line comments...")
        contents = join([rstrip(split(line, "#")[1]) for line in split(contents, "\n")], "\n")
        # Remove lines with only whitespaces:
        println("   Cleaning white spaces...")
        contents = clean_whitespaces_and_newlines(contents)

        contents = join([line for line in split(contents, "\n")], "\n")

        flag_set = FlagSet([
        FlagPair{GreedyType}(
        StartFlag("'''", ["\n", " ", ""], ["\n", " ", ""]),
        StopFlag("'''", ["\n", " ", ""], ["\n", " ", ""])
        ),
        FlagPair{GreedyType}(
        StartFlag("\"\"\"", ["\n", " ", ""], ["\n", " ", ""]),
        StopFlag("\"\"\"", ["\n", " ", ""], ["\n", " ", ""])
        )])
        println("   Creating token stream...")
        ts = TokenStream(contents, flag_set)
        println("   Finished creating token stream...")

        fs = ("'''-'''", "\"\"\"-\"\"\"")
        i_code = [i for i in eachindex(ts.code) if cond(ts, i, fs[1],fs[2])]
        contents = join([SubString(ts.code, i, i) for i in i_code])

        println("   Cleaning white spaces...")
        contents = clean_whitespaces_and_newlines(contents)

        open(file, "w") do io
            print(io, contents)
        end

    end

end


text:
# MKL now supports 3D transforms, so...
# ENV["JULIA_FFTW_PROVIDER"]="FFTW"
ENV["JULIA_FFTW_PROVIDER"]="MKL"