text:
using Test, GenericLinearAlgebra, LinearAlgebra

@testset "The General eigenvalue problem" begin

cplxord = t -> (real(t), imag(t))

@testset "General eigen problem with n=$n and element type=$T" for
    n in (10, 23, 100),
    T in (Float64, Complex{Float64})

    A = randn(T, n, n)
    vGLA    = GenericLinearAlgebra._eigvals!(copy(A))
    vLAPACK = eigvals(A)
    vBig    = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra
    @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord)
    @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord)
    if VERSION > v"1.2.0-DEV.0"
        @test issorted(vBig, by = cplxord)
    end
end

@testset "make sure that solver doesn't hang" begin
    for i in 1:1000
        A = randn(8, 8)
        sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A)))
    end
end

@testset "Convergence in corner cases. Issue 29." begin
    function H(n::Int)
        H = zeros(2n, 2n)
        for i = 1 : 2 : 2n
            H[i, i+1] = 1
            H[i+1, i] = 1
        end
        H
    end

    function E(n::Int)
        E = zeros(2n, 2n)
        for i = 1 : (n - 1)
            E[2i + 1, 2i] = 1
        end
        E[1, 2n] = 1
        E
    end

    my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n)

    A = my_matrix(4, 1e-3);
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
end

@testset "Convergence in with 0s Issue 58." begin
    A = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 0.0]
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
    B = [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 -1.0 0.0]
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(B))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(B), by = t -> (real(t), imag(t)))
end

@testset "Extract Schur factor" begin
    A = randn(5, 5)
    @test sum(eigvals(schur(A).T))     ≈ sum(eigvals(Float64.(schur(big.(A)).T)))
    @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur)))
end

@testset "Issue 63" begin
    A = [1 2 1 1 1 1
         0 1 0 1 0 1
         1 1 2 0 1 0
         0 1 0 2 1 1
         1 1 1 0 2 0
         0 1 1 1 0 2]

    z1 = [complex(1, sqrt(big(3))), complex(1, -sqrt(big(3)))]
    z2 = [
        4*2^(big(2)/3)/complex(5, 3*sqrt(big(111)))^(big(1)/3),
        complex(5, 3*sqrt(big(111)))^(big(1)/3)/2^(big(2)/3)
        ]

    truevals = real.([
        (7 - transpose(z1)*z2),
        3,
        3,
        3,
        (7 - adjoint(z1)*z2),
        (7 + [2, 2]'*z2)])/3

    if VERSION > v"1.2.0-DEV.0"
        @test eigvals(big.(A)) ≈ truevals
    else
        @test sort(eigvals(big.(A)), by=real) ≈ truevals
    end
end

Demmel(η) = [0  1 0 0
             1  0 η 0
             0 -η 0 1
             0  0 1 0]

@testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8)
    # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it"
    A = Demmel(t)
    vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter=35))
    @test abs.(vals) ≈ ones(4)
end

function Hevil2(θ, κ, α, γ)
    # Eq (13) and (14)
    β = ω = 0.0
    ν = cos(θ)*cos(2γ) + cos(α + β + ω)*sin(2γ)*κ/2
    σ = 1 + κ*sin(2γ)*cos(α + β + ω - θ) + κ^2*sin(γ)^2
    μ = -sin(θ)*cos(2γ) - sin(α + β + ω)*sin(2γ)*κ/2
    ρ = sqrt(σ - ν^2)

    return [ν    (cos(2θ) - ν^2)/ρ              μ/ρ*(cos(2θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)      (-2*μ*ν - sin(2*θ))/sqrt(ρ^2 - μ^2)
            ρ   -ν - sin(2*θ)*μ/ρ^2            -μ/ρ^2*(μ*sin(2*θ) + 2*ν*ρ^2)/sqrt(ρ^2 - μ^2)    -μ/ρ*(cos(2*θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)
            0    sin(2*θ)*sqrt(ρ^2 - μ^2)/ρ^2   ν + sin(2*θ)*μ/ρ^2                              (cos(2*θ) - ν^2)/ρ
            0    0                              ρ                                               -ν]
end

@testset "Complicate matrix from Sandia technical report" begin
    H = Hevil2(0.111866322512629152, 1.08867072154101741, 0.338146383137297168, -0.313987810419091240)

    @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4)
end

@testset "Issue 67" for (A, λs) in (
    ([1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]    ,
        [0.0, 0.0, (3 - √big(3)*im)/2, (3 + √big(3)*im)/2, (3 - √big(3)im)/2, (3 + √big(3)im)/2]),
    ([1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]      ,
        zeros(6)),
    ([1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0],
        [(1 - √big(3)*im)/2, (1 + √big(3)*im)/2, (1 - √big(3)im)/2, (1 + √big(3)im)/2, 2, 2]),
    ([2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]    ,
        ones(6)),
    ([1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1] ,
        [-1, -1, 0, 0, 0, 0]))

    @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin
        GenericLinearAlgebra._schur!(float(A))
    end

    # For BigFloats, many iterations are required for convergence
    # Improving this is probably a hard problem
    vals = eigvals(big.(A), maxiter=1500)

    # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately
    @test sort(real(vals)) ≈ sort(real(λs)) atol=1e-25
    @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25
end

end


text:
struct Trajectory{T}
    q::Vector{Vector{T}} # generalized coordinates
    v::Vector{Vector{T}} # generalized velocities
    u::Vector{Vector{T}} # control inputs
    γ::Vector{Vector{T}} # impact impulses
    b::Vector{Vector{T}} # friction impulses
    w::Vector{Vector{T}} # disturbances
end

function Trajectory(model, T; nv=model.nq, nc=model.nc, nb=model.nc)
    nq = model.nq
    nu = model.nu
    nw = model.nw

    q = [zeros(nq) for t = 1:T+2]
    v = [zeros(nv) for t = 1:T+1] # midpoint velocities; TODO: knotpoint velocities
    u = [zeros(nu) for t = 1:T]
    γ = [zeros(nc) for t = 1:T]
    b = [zeros(nb) for t = 1:T]
    w = [zeros(nw) for t = 1:T]
    Trajectory(q, v, u, γ, b, w)
end

function reset!(traj::Trajectory)
    T = length(traj.u)
    for t = 1:T
        fill!(traj.q[t], 0.0)
        fill!(traj.v[t], 0.0)
        fill!(traj.u[t], 0.0)
        fill!(traj.γ[t], 0.0)
        fill!(traj.b[t], 0.0)
        fill!(traj.w[t], 0.0)
    end
    fill!(traj.q[T+1], 0.0)
    fill!(traj.q[T+2], 0.0)
    fill!(traj.v[T+1], 0.0)
    return nothing
end

struct GradientTrajectory{T}
    ∂q3∂q1::Vector{Matrix{T}}
    ∂q3∂q2::Vector{Matrix{T}}
    ∂q3∂u1::Vector{Matrix{T}}
    ∂γ1∂q1::Vector{Matrix{T}}
    ∂γ1∂q2::Vector{Matrix{T}}
    ∂γ1∂u1::Vector{Matrix{T}}
    ∂b1∂q1::Vector{Matrix{T}}
    ∂b1∂q2::Vector{Matrix{T}}
    ∂b1∂u1::Vector{Matrix{T}}
end

function GradientTrajectory(model, T; nc=model.nc, nb=model.nc)
    nq = model.nq
    nu = model.nu

    ∂q3∂q1 = [zeros(nq, nq) for t = 1:T]
    ∂q3∂q2 = [zeros(nq, nq) for t = 1:T]
    ∂q3∂u1 = [zeros(nq, nu) for t = 1:T]
    ∂γ1∂q1 = [zeros(nc, nq) for t = 1:T]
    ∂γ1∂q2 = [zeros(nc, nq) for t = 1:T]
    ∂γ1∂u1 = [zeros(nc, nu) for t = 1:T]
    ∂b1∂q1 = [zeros(nb, nq) for t = 1:T]
    ∂b1∂q2 = [zeros(nb, nq) for t = 1:T]
    ∂b1∂u1 = [zeros(nb, nu) for t = 1:T]

    GradientTrajectory(∂q3∂q1, ∂q3∂q2, ∂q3∂u1,
                       ∂γ1∂q1, ∂γ1∂q2, ∂γ1∂u1,
                       ∂b1∂q1, ∂b1∂q2, ∂b1∂u1)
end

function reset!(traj::GradientTrajectory)
    T = length(traj.∂q3∂q1)
    for t = 1:T
        fill!(traj.∂q3∂q1[t], 0.0)
        fill!(traj.∂q3∂q2[t], 0.0)
        fill!(traj.∂q3∂u1[t], 0.0)
        fill!(traj.∂γ1∂q1[t], 0.0)
        fill!(traj.∂γ1∂q2[t], 0.0)
        fill!(traj.∂γ1∂u1[t], 0.0)
        fill!(traj.∂b1∂q1[t], 0.0)
        fill!(traj.∂b1∂q2[t], 0.0)
        fill!(traj.∂b1∂u1[t], 0.0)
    end
    return nothing
end

ext: jl


text:
module Grid

function grid_xp(n::Int64; pad=false::Bool)
        if pad
                x = zeros(n + 1)
        else
                x = zeros(n)
        end
        h = 1.0 / (n - 1)
        for i in 1:n
                x[i] = (i - 1) * h
        end

        return x, h
end

function grid_xm(n::Int64)
        x = zeros(n)
        h = 1.0 / (n - 2)
        for i in 1:n
                x[i] = (i - 1 - 0.5) * h
        end
        x[1] = 0.0
        x[end] = 1.0

        return x, h
end

"""
Construct a 2D grid vector from a 1D grid vector `x1`. The resulting vector
repeats all x-values: x1[1], x1[1], ... , x1[1], x1[2], x1[2] ... x1[2] `ny`
times.
"""
function grid_2d_x(x1::AbstractArray, ny::Int64)
        nx = size(x1, 1)
        x = zeros(nx * ny)
        for i=1:nx
                for j=1:ny
                        @inbounds x[j + (i - 1) * ny] = x1[i]
                end
        end
        return x
end

"""
Construct a 2D grid vector from a 1D grid vector `y1`. The resulting vector
repeats all y-values: y1[1], y1[1], ... , y1[ny], y1[1], y1[2] ... y1[ny] `nx`
times.
"""
function grid_2d_y(y1::AbstractArray, nx::Int64)
        ny = size(y1, 1)
        y = zeros(nx * ny)
        for i=1:nx
                for j=1:ny
                        @inbounds y[j + (i - 1) * ny] = y1[j]
                end
        end
        return y
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMICore

# check float size (32 or 64 bits)
juliaArch = Sys.WORD_SIZE
@assert (juliaArch == 64 || juliaArch == 32) "FMICore: Unknown Julia Architecture with $(juliaArch)-bit, must be 64- or 32-bit."
Creal = Cdouble
if juliaArch == 32
    Creal = Cfloat
end

"""
The mutable struct representing an abstract (version unknown) FMU.
"""
abstract type FMU end
abstract type FMUSolution end

include("FMI2_c.jl")
include("FMI2.jl")

include("FMI3_c.jl")
include("FMI3.jl")

### EXPORTING LISTS START ###

export FMU
export FMU2ExecutionConfiguration, FMU_EXECUTION_CONFIGURATION_RESET, FMU_EXECUTION_CONFIGURATION_NO_RESET, FMU_EXECUTION_CONFIGURATION_NO_FREEING

# FMI2.jl
export FMU2, FMU2Component, FMU2ComponentEnvironment, FMU2Solution, FMU2Event
export fmi2StatusToString
export fmi2CausalityToString, fmi2StringToCausality
export fmi2VariabilityToString, fmi2StringToVariability
export fmi2InitialToString, fmi2StringToInitial
export fmi2DependencyKindToString, fmi2StringToDependencyKind

# FMI2_c.jl

# data types
export fmi2ValueReference, fmi2Component
export fmi2Char, fmi2String, fmi2Boolean, fmi2Real, fmi2Integer, fmi2Byte, fmi2FMUstate, fmi2ComponentEnvironment, fmi2Enum
export fmi2Status, fmi2StatusOK, fmi2StatusWarning, fmi2StatusDiscard, fmi2StatusError, fmi2StatusFatal, fmi2StatusPending
export fmi2CallbackFunctions
export fmi2Annotation
export fmi2Causality, fmi2CausalityParameter, fmi2CausalityCalculatedParameter, fmi2CausalityInput, fmi2CausalityOutput, fmi2CausalityLocal, fmi2CausalityIndependent
export fmi2Variability, fmi2VariabilityConstant, fmi2VariabilityFixed, fmi2VariabilityTunable, fmi2VariabilityDiscrete, fmi2VariabilityContinuous
export fmi2Initial, fmi2InitialExact, fmi2InitialApprox, fmi2InitialCalculated
export fmi2True, fmi2False
export fmi2Type, fmi2TypeModelExchange, fmi2TypeCoSimulation
export fmi2DatatypeVariable
export fmi2ScalarVariable
export fmi2StatusKind, fmi2StatusKindDoStepStatus, fmi2StatusKindPendingStatus, fmi2StatusKindLastSuccessfulTime, fmi2StatusKindTerminated
export fmi2EventInfo
export fmi2Unit, fmi2SimpleType
export fmi2DependencyKind, fmi2DependencyKindDependent, fmi2DependencyKindConstant, fmi2DependencyKindFixed, fmi2DependencyKindTunable, fmi2DependencyKindUnknown
export fmi2VariableNamingConvention, fmi2VariableNamingConventionFlat, fmi2VariableNamingConventionStructured
export fmi2VariableDependency
export fmi2DefaultExperiment, fmi2Unknown, fmi2ModelStructure
export fmi2ModelDescription
export fmi2ComponentState, fmi2ComponentStateInstantiated, fmi2ComponentStateInitializationMode, fmi2ComponentStateEventMode, fmi2ComponentStateContinuousTimeMode, fmi2ComponentStateTerminated, fmi2ComponentStateError, fmi2ComponentStateFatal

# functions
export fmi2CallbackLogger, fmi2CallbackAllocateMemory, fmi2CallbackFreeMemory, fmi2CallbackStepFinished
export fmi2Instantiate, fmi2FreeInstance!, fmi2GetTypesPlatform, fmi2GetVersion
export fmi2SetDebugLogging, fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
export fmi2GetFMUstate!, fmi2SetFMUstate, fmi2FreeFMUstate!, fmi2SerializedFMUstateSize!, fmi2SerializeFMUstate!, fmi2DeSerializeFMUstate!
export fmi2GetDirectionalDerivative!, fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives!
export fmi2DoStep, fmi2CancelStep, fmi2GetStatus!, fmi2GetRealStatus!, fmi2GetIntegerStatus!, fmi2GetBooleanStatus!, fmi2GetStringStatus!
export fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscreteStates!, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!
export fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!

# FMI3.jl
export FMU3, FMU3Component
export fmi3StatusToString

# FMI3_c.jl

# data types
export fmi3ValueReference, fmi3ModelDescription, fmi3Instance, fmi3ModelVariable, fmi3DatatypeVariable
export fmi3False, fmi3True
export fmi3Status, fmi3StatusOK, fmi3StatusWarning, fmi3StatusDiscard, fmi3StatusError, fmi3StatusFatal
export fmi3Causality, fmi3CausalityParameter, fmi3CausalityCalculatedParameter, fmi3CausalityInput, fmi3CausalityOutput, fmi3CausalityLocal, fmi3CausalityIndependent, fmi3CausalityStructuralParameter
export fmi3Initial, fmi3InitialExact, fmi3InitialApprox, fmi3InitialCalculated
export fmi3Variability, fmi3VariabilityConstant, fmi3VariabilityFixed, fmi3VariabilityTunable, fmi3VariabilityDiscrete, fmi3VariabilityContinuous
export fmi3Type, fmi3TypeModelExchange, fmi3TypeCoSimulation, fmi3TypeScheduledExecution
export fmi3Clock
export fmi3IntervalQualifier, fmi3IntervalQualifierIntervalChanged, fmi3IntervalQualifierIntervalNotYetKnown, fmi3IntervalQualifierIntervalUnchanged
export fmi3Char, fmi3String, fmi3Boolean, fmi3Binary, fmi3Float32, fmi3Float64, fmi3Int8, fmi3UInt8, fmi3Int16, fmi3UInt16, fmi3Int32, fmi3UInt32, fmi3Int64, fmi3UInt64, fmi3Byte, fmi3FMUState, fmi3InstanceEnvironment, fmi3Enum
export fmi3DependencyKind, fmi3DependencyKindDependent, fmi3DependencyKindIndependent, fmi3DependencyKindConstant, fmi3DependencyKindTunable, fmi3DependencyKindDiscrete, fmi3DependencyKindFixed, fmi3DependencyKindGetVariableDependencies

# functions
# ToDo

### EXPORTING LISTS END ###

end # module


text:
"""
Gives a reward of +1 if sp is the target lane
"""
type TargetLaneReward <: AbstractMLRewardModel
    target_lane::Int
end
function reward{D<:AbstractMLDynamicsModel}(mdp::MLMDP{MLState, MLAction, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState)
    return isnull(s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

function reward{D<:AbstractMLDynamicsModel}(mdp::MLPOMDP{MLState, MLAction, MLPhysicalState, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState)
    return isnull(s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
using WaterLily
using BenchmarkTools
using LinearAlgebra: norm2
using Profile

radius = 8; Re = 250
body = AutoBody((x,t)-> norm2(x .- 2radius) - radius)
body_sim = ()->Simulation((6radius+2,4radius+2),[1.,0.],radius; body, ν=radius/Re)

L = 2^5
function uλ(i,vx)
    x,y,z = @. (vx-1.5)*π/L              # scaled coordinates
    i==1 && return -sin(x)*cos(y)*cos(z) # u_x
    i==2 && return  cos(x)*sin(y)*cos(z) # u_y
    return 0.                            # u_z
end
vort_sim = ()->Simulation((L+2,L+2,L+2),zeros(3),L;uλ,ν=L/Re,U=1.)

function benchmark_step()
    @benchmark mom_step!($sim.flow,$sim.pois)
end

benchmark_vort() = @benchmark sim_step!(sim,10) setup=(sim=vort_sim()) seconds=30
benchmark_body() = @benchmark sim_step!(sim,10) setup=(sim=body_sim())


text:
# ----- Number -----
typealias AbstractUniform Union{DiscreteUniform, Uniform}

immutable NumberStrategy{D <: AbstractUniform, T} <: AbstractStrategy
    dist::D
    specialvalues::Vector{T}
    p::Bernoulli
    shrinker::NoShrinker
end
Base.minimum(strategy::NumberStrategy) = dist.params(strategy.dist)[1]
shrinker(strategy::NumberStrategy) = strategy.shrinker
returntype{D, T}(::Type{NumberStrategy{D, T}}) = T

# ----- Float -----
function _specialvalues{T <: AbstractFloat}(::Type{T}, min::T, max::T)
    specialvalues = [convert(T, NaN)]
    append!(specialvalues, [min, max])
    if (min < zero(T)) && (max > zero(T))
        push!(specialvalues, zero(T))
    end
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    if (min < -one(T)) && (max > -one(T))
        push!(specialvalues, -one(T))
    end
    return specialvalues
end

"""
Generate arbitrary float values between a desired range
"""
function floats{T <: AbstractFloat}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(Uniform(min, max), specialvalues, p, shrinker)
end

# ----- Int -----
function _specialvalues{T <: Signed}(::Type{T}, min::T, max::T)
    specialvalues = [min, max]
    if (min < zero(T)) && (max > zero(T))
        push!(specialvalues, zero(T))
    end
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    if (min < -one(T)) && (max > -one(T))
        push!(specialvalues, -one(T))
    end
    return specialvalues
end
"""
Generate arbitrary signed values between a desired range
"""
function signeds{T <: Signed}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(DiscreteUniform(min, max), specialvalues, p, shrinker)
end

function _specialvalues{T <: Unsigned}(::Type{T}, min::T, max::T)
    specialvalues = [min, max]
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    return specialvalues
end
"""
Generate arbitrary unsigned values between a desired range
"""
function unsigneds{T <: Unsigned}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(DiscreteUniform(min, max), specialvalues, p, shrinker)
end

function arbitrary{D, T}(strategy::NumberStrategy{D, T})
    if Bool(rand(strategy.p))
        if length(strategy.specialvalues) > 0
            return convert(T, rand(strategy.specialvalues))
        end
    end
    return convert(T, rand(strategy.dist))
end
shrink{T <: Number}(x::T, strategy::NumberStrategy) = Task(() -> nothing)

# ----- Booleans -----
immutable BoolStrategy <: AbstractStrategy
    bernoulli::Bernoulli
end
BoolStrategy(p::Float64) = BoolStrategy(Bernoulli(p))

bools(; p=0.5) = BoolStrategy(p)

arbitrary(strategy::BoolStrategy) = rand(strategy.bernoulli)
shrink(x::Bool, strategy::BoolStrategy) = Bool[]

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
using Distributed
if nprocs() == 1
    addprocs(3)
end

mat = Dict()
println("----------------alskdjflaskfjasklf laksjflaskj-=-------")
mat["5x5x5"] = [1,2,3]
mat["15x15x15"] = [5,6]

matice = @spawn mat;
# mat[1] = 7
# matice = @spawn mat;

@everywhere function used_in_smallfunc(par1, par2)
  println(par1, par2, " ")
  return par1 + par2
end

function bigfunc(iterator)
  localvar = 2.5
  pars = (1, 15)

  smallfunc = function (elm)
    # hypothetical (expensive) operation with elm
    # (in my code this is solving a linear system)
    vysledek = localvar * elm
    global matice

    # println("$localvar * $elm = $vysledek ")
    # matice
    # print("pred fetch $matice")
    # fetch(matice);
    # print("po fetch $matice")
    local_matrix_storage = fetch(matice)
    print(local_matrix_storage)
    if local_matrix_storage["5x5x5"][1] == 1
      local_matrix_storage["5x5x5"][1] = 2
      matice = @spawn local_matrix_storage
    else
      print(" else...")
    end
    # used_in_smallfunc(pars...)  # LoadError
    used_in_smallfunc(pars[1], pars[2])
  end

  pmap(smallfunc, iterator)
end

bigfunc(1:10)


text:
include("../src/DataFile.jl")
include("../src/SolveBalances.jl")
using Debug

# =============================================================================== #
# Solve.jl
# Function library to solve the EMT model for different experimental/test cases
# Author: J.Varner
#
# Version: 1.0
# Date: 12/17/15
# =============================================================================== #

# addition indices
addition_start,addition_end = 85,100

# mRNA / transcriptional indices
mdeg_start,mdeg_end = 101,123
script_start,script_end = 124,146

# protein / translational indices
pdeg_start,pdeg_end = 39,78
lation_start,lation_end = 147,169

# Transcriptional mult
tscript_mult = 0.2 # 0.033

#@debug
# Run the model to steady state -
 function EstimateSteadyState(current_data_dictionary)
# WARNING UNLESS DEEPCOPY IS USED, CURRENT_DATA_DICTIONARY IS THE SAME OBJECT AS DATA_DICIONARY!
  # Do we have a data_dictionary (if not load a new one)
  if (isempty(current_data_dictionary))
    # DB 02/27/2016
    TSTART = 0.0
    Ts=1.0
    TSTOP=100
    # DB
    current_data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  #test = data_dictionary["INITIAL_CONDITION_ARRAY"]
  #println(data_dictionary)
  include("../src/edit_file_global.jl")
  UpdateArray(current_data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(current_data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

  rate_constant_array = copy(current_data_dictionary["RATE_CONSTANT_ARRAY"]);
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  current_data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array

# 137 - 139
# rate_constant_array[137:139]=30*rate_constant_array[137:139]*degradation_constant_mRNA

  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)
  # Setup loop -
  #TSIM_all = []
  #X_all    = Array(Float64,length(data_dictionary["INITIAL_CONDITION_ARRAY"]))
  # tmpI = 0 # save
  EPSILON = 1e-3;
  TSTART = 0.0;
  Ts = 1.0;
  TSTOP = 100;
  did_reach_steady_state = false
  #@show current_data_dictionary["RATE_CONSTANT_ARRAY"][141] #here

  while (!did_reach_steady_state)
    # println("ITERATION",)
    # solve the balances -
  #  @show current_data_dictionary["RATE_CONSTANT_ARRAY"][141]

    (TSIM,X1) = SolveBalances(TSTART,TSTOP,Ts,current_data_dictionary);
  #  @show current_data_dictionary["RATE_CONSTANT_ARRAY"][141]
    ###save
    # if tmpI == 0
    #   X_all = X1
    # else
    # X_all=vcat(X_all,X1) # Save
    # end
    # tmpI = tmpI+1
    # Take a few additional steps -
    TNEXT_START = TSTOP+Ts;
    TNEXT_STOP = TNEXT_START+1.0;
    Ts = 0.1;

    # solve the balances again 0

    initial_condition_array = vec(X1[end,:])
    current_data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

    (TSIM,X2) = SolveBalances(TNEXT_START,TNEXT_STOP,Ts,current_data_dictionary);
    # X_all=vcat(X_all,X2) # Save
    # view simulation trajectory

    # Find the difference -
    DIFF = norm((X2[end,:] - X1[end,:]));
    # Take the end, the last step is the one we compare
    # Should we stop -or- go around again?

# println(size(X1))
# println(size(X2))
#@bp

    if (DIFF<EPSILON)
      did_reach_steady_state = true;
      # savepath = "../output/all_XSS.dat"
      # println("Save",size(X_all))
      # writedlm(savepath,X_all)# save
      return (X2[end,:]);

    else

      # No, we did *not* reach steady state ....
      TSTART = TSTOP+Ts
      TSTOP = 1000 + TSTART;
      Ts = 100.0;

      initial_condition_array = vec(X2[end,:])
      current_data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
      # println("else")
    end
  end
  # return
  # println("final return")
  return XSS;

end

@debug function SolveModelGeneric_ram(TSTART,TSTOP,Ts,data_dictionary_active,sim_RCA)

  include("../src/edit_file_global.jl")
  UpdateArray(data_dictionary_active,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)

  ### THIS WAS MISSING BECAUSE IT WAS REPLACED BY SIM_RCA... BUT NEEDS TO COME BEFORE IT
  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  ### IT HAD BEEN COMPESATED FOR BY THE MEMORY LEAK FROM THE STEADYSTATE FUNCTION ..

  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",sim_RCA)
  # @show sim_RCA
  # @show data_dictionary_active["RATE_CONSTANT_ARRAY"]
  # @bp
  # Run the simulation -
  rate_constant_array = data_dictionary_active["RATE_CONSTANT_ARRAY"];
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  #@show rate_constant_array # inspect scaled array
  data_dictionary_active["RATE_CONSTANT_ARRAY"] = rate_constant_array

# 137 - 139
# rate_constant_array[137:139]=30*rate_constant_array[137:139]*degradation_constant_mRNA

#

  #@bp
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary_active);
  #@show data_dictionary_active["RATE_CONSTANT_ARRAY"][141]
  #  @show data_dictionary_active["RATE_CONSTANT_ARRAY"][182]

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

function SolveModelGenericPOETS(TSTART,TSTOP,Ts,local_data_dictionary,simulation)

  include(simulation)
  # !!!! I think th initial condition array update is not required
  #UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(local_data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

  # Run the simulation -
  # println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,local_data_dictionary);

  # return tuple -
  return (TSIM,X)

end

#  #
#  #
#  #
#  #

function SolveModelGenericPOETS_ram(TSTART,TSTOP,Ts,local_data_dictionary,sim_RCA)
  # Set simulation rate constants
  UpdateArray(local_data_dictionary,"RATE_CONSTANT_ARRAY",sim_RCA)
  # Scale parameters
  rate_constant_array = local_data_dictionary["RATE_CONSTANT_ARRAY"];
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  local_data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array
  # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,local_data_dictionary);

  # return tuple -
  return (TSIM,X)

end

# Simulate the addition of TGFB -
function SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant
  include("../src/edit_file_global.jl")

  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  #println(data_dictionary)
  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add TGFB1/2 -
  # RESUME TGFB signaling for simulation
  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_ss.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  # Run the simulation -
  #println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

# Simulate the addition of VEGF -
function SolveModelAddVEGF(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  # Modify rate constant array

  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_global.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add VEGF -
  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_vegf.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

    # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

function SolveModelGeneric(TSTART,TSTOP,Ts,data_dictionary_active,simulation)

  # Do we have a data_dictionary (if not load a new one)
  # if (isempty(data_dictionary_active))
  #   data_dictionary_active = DataFile(TSTART,TSTOP,Ts)
  # end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  include("../src/edit_file_global.jl")
  UpdateArray(data_dictionary_active,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  #UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  #println(data_dictionary)
  # Run the model to steady-state *before* we add the TGFB -
  #XSS = EstimateSteadyState(data_dictionary);
#  println(data_dictionary["INITIAL_CONDITION_ARRAY"])
  # ok, we established a steady state -
  # reset the initial conditions -
  #initial_condition_array = vec(XSS);
  #data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add TGFB1/2 -
  # RESUME TGFB signaling for simulation
  include(simulation)
  #UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #println(data_dictionary["INITIAL_CONDITION_ARRAY"]) ###
  # Run the simulation -
  # println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary_active);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

#


text:
function check_data_sizes(data::AbstractMatrix, raxis, (l,u), (λ,μ))
    if blocksize(data,1) ≠ l + u + 1 && !(blocksize(data,1) == 0 && (-l > u || -λ > μ))
        throw(ArgumentError("Data matrix must have number of row blocks equal to number of block bands"))
    end
    for K = blockaxes(data,1)
        if length(axes(data,1)[K]) ≠ λ + μ + 1 && !(-λ > μ)
            throw(ArgumentError("Data matrix must have row block sizes equal to number of subblock bands"))
        end
    end
end

function _BandedBlockBandedMatrix end

# Represents a block banded matrix with banded blocks
#   similar to BandedMatrix{BandedMatrix{T,Matrix{T}}}
# Here the data is stored by blocks, in a way that is consistent with
# BandedMatrix
#

struct BandedBlockBandedMatrix{T, BLOCKS, RAXIS<:AbstractUnitRange{Int}} <: AbstractBlockBandedMatrix{T}
    data::BLOCKS
    raxis::RAXIS

    l::Int  # block lower bandwidth
    u::Int  # block upper bandwidth
    λ::Int  # sub lower bandwidth
    μ::Int  # sub upper bandwidth

    global function _BandedBlockBandedMatrix(data::AbstractMatrix,
                                             raxis::AbstractUnitRange{Int},
                                             lu::NTuple{2,Int}, λμ::NTuple{2,Int})
      @boundscheck check_data_sizes(data, raxis, lu, λμ)
      new{eltype(data), typeof(data), typeof(raxis)}(data, raxis, lu..., λμ...)
    end
end

const DefaultBandedBlockBandedMatrix{T} = BandedBlockBandedMatrix{T, PseudoBlockMatrix{T, Matrix{T}, NTuple{2,DefaultBlockAxis}}, DefaultBlockAxis}

@inline _BandedBlockBandedMatrix(data::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    _BandedBlockBandedMatrix(PseudoBlockArray(data,(blockedrange(Fill(sum(λμ)+1,sum(lu)+1)),axes[2])), axes[1], lu, λμ)

_BandedBlockBandedMatrix(data::AbstractMatrix,rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    _BandedBlockBandedMatrix(data, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

_blocklengths2blocklasts(b::Fill) = iszero(FillArrays.getindex_value(b)) ? (1:1:0) : cumsum(b)
_bbb_data_axes(caxes, lu, λμ) = (blockedrange(Fill(max(0,sum(λμ)+1),max(0,sum(lu)+1))),caxes)

BandedBlockBandedMatrix{T,B,R}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
    _BandedBlockBandedMatrix(B(undef, _bbb_data_axes(axes[2],lu,λμ)), axes[1], lu, λμ)

BandedBlockBandedMatrix{T,B,R}(::UndefInitializer, rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
    BandedBlockBandedMatrix{T,B,R}(undef, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

BandedBlockBandedMatrix{T,B}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
    BandedBlockBandedMatrix{T,B,typeof(axes[1])}(undef,axes,lu,λμ)

BandedBlockBandedMatrix{T,B}(::UndefInitializer, rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
    BandedBlockBandedMatrix{T,B}(undef, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

"""
    BandedBlockBandedMatrix{T}(M::Union{UndefInitializer,UniformScaling,AbstractMatrix},
                               rows, cols, (l, u), (λ, μ))

returns a `sum(rows)`×`sum(cols)` banded-block-banded matrix `A` of
type `T` with block-bandwidths `(l,u)` and where `A[Block(K,J)]` is a
`BandedMatrix{T}` of size `rows[K]`×`cols[J]` with bandwidths `(λ,μ)`.

# Examples

```jldoctest
julia> BandedBlockBandedMatrix(I, [3,4,3], [3,4,3], (1,1), (1,1))
3×3-blocked 10×10 BandedBlockBandedMatrix{Bool,BlockArrays.PseudoBlockArray{Bool,2,Array{Bool,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
 1  0  ⋅  │  0  0  ⋅  ⋅  │  ⋅  ⋅  ⋅
 0  1  0  │  0  0  0  ⋅  │  ⋅  ⋅  ⋅
 ⋅  0  1  │  ⋅  0  0  0  │  ⋅  ⋅  ⋅
 ─────────┼──────────────┼─────────
 0  0  ⋅  │  1  0  ⋅  ⋅  │  0  0  ⋅
 0  0  0  │  0  1  0  ⋅  │  0  0  0
 ⋅  0  0  │  ⋅  0  1  0  │  ⋅  0  0
 ⋅  ⋅  0  │  ⋅  ⋅  0  1  │  ⋅  ⋅  0
 ─────────┼──────────────┼─────────
 ⋅  ⋅  ⋅  │  0  0  ⋅  ⋅  │  1  0  ⋅
 ⋅  ⋅  ⋅  │  0  0  0  ⋅  │  0  1  0
 ⋅  ⋅  ⋅  │  ⋅  0  0  0  │  ⋅  0  1

julia> BandedBlockBandedMatrix(Ones{Int}(10,13), [3,4,3], [4,5,4], (1,1), (1,1))
3×3-blocked 10×13 BandedBlockBandedMatrix{Int64,BlockArrays.PseudoBlockArray{Int64,2,Array{Int64,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
 1  1  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 1  1  1  ⋅  │  1  1  1  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  1  1  1  │  ⋅  1  1  1  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ────────────┼─────────────────┼────────────
 1  1  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅
 1  1  1  ⋅  │  1  1  1  ⋅  ⋅  │  1  1  1  ⋅
 ⋅  1  1  1  │  ⋅  1  1  1  ⋅  │  ⋅  1  1  1
 ⋅  ⋅  1  1  │  ⋅  ⋅  1  1  1  │  ⋅  ⋅  1  1
 ────────────┼─────────────────┼────────────
 ⋅  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  │  1  1  1  ⋅  ⋅  │  1  1  1  ⋅
 ⋅  ⋅  ⋅  ⋅  │  ⋅  1  1  1  ⋅  │  ⋅  1  1  1
```
"""
BandedBlockBandedMatrix

BandedBlockBandedMatrix{T}(::UndefInitializer, args...) where T =
    DefaultBandedBlockBandedMatrix{T}(undef, args...)

# Auxiliary outer constructors

function convert(::Type{<:BandedBlockBandedMatrix}, B::BandedMatrix)
    if isdiag(B) # TODO really should be one block
        _BandedBlockBandedMatrix(copy(B.data), Fill(1,size(B,1)), Fill(1,size(B,2)), (0,0), (0,0))
    else
        _BandedBlockBandedMatrix(copy(B.data), [size(B,1)], [size(B,2)], (0,0), (B.l,B.u))
    end
end

convert(::Type{BandedBlockBandedMatrix{T,BLOCKS,RAXIS}}, A::BandedBlockBandedMatrix) where {T,BLOCKS,RAXIS} =
    _BandedBlockBandedMatrix(convert(BLOCKS, A.data), convert(RAXIS, A.raxis), (A.l, A.u), (A.λ, A.μ))

function BandedBlockBandedMatrix{T,B,R}(Z::Zeros, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
   if size(Z) ≠ map(length,axes)
       throw(DimensionMismatch())
   end
   blocks = fill!(B(undef, _bbb_data_axes(axes[2],lu,λμ)), zero(T))
   _BandedBlockBandedMatrix(blocks, convert(R, axes[1]), lu, λμ)
end

function BandedBlockBandedMatrix{T,B,R}(E::Eye, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
    if size(E) ≠ map(length,axes)
        throw(DimensionMismatch())
    end
    ret = BandedBlockBandedMatrix{T,B,R}(Zeros{T}(E), axes, lu, λμ)
    ret[diagind(ret)] .= one(T)
    ret
end

function BandedBlockBandedMatrix{T,B,R}(A::UniformScaling, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
    ret = BandedBlockBandedMatrix{T,B,R}(Zeros{T}(map(length,axes)), axes, lu, λμ)
    ret[diagind(ret)] .= convert(T, A.λ)
    ret
end

BandedBlockBandedMatrix{T,B}(m::Union{AbstractMatrix, UniformScaling},
                            axes::NTuple{2,AbstractUnitRange{Int}},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
  BandedBlockBandedMatrix{T,B,typeof(axes[1])}(m, axes, lu, λμ)

BandedBlockBandedMatrix{T}(m::Union{AbstractMatrix, UniformScaling},
                            axes::NTuple{2,AbstractUnitRange{Int}},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  DefaultBandedBlockBandedMatrix{T}(m, axes, lu, λμ)

BandedBlockBandedMatrix{T,B,R}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
  BandedBlockBandedMatrix{T,B,R}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix{T,B}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
  BandedBlockBandedMatrix{T,B}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix{T}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  BandedBlockBandedMatrix{T}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        axes::NTuple{2, AbstractUnitRange{Int}}, lu::NTuple{2,Int},
                        λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix{eltype(A)}(A, axes, lu, λμ)
BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                        lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix(A, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

function BandedBlockBandedMatrix{T,Blocks,RR}(A::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,Blocks,RR<:AbstractUnitRange{Int}}
    ret = BandedBlockBandedMatrix{T,Blocks,RR}(Zeros{T}(size(A)), axes, lu, λμ)
    L,M = λμ
    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
        kr, jr = axes[1][K], axes[2][J]

        # We have the correct block - now we need to only add the entries from
        # the correct bands
        B = view(A, kr, jr)
        R = view(ret, K, J)
        rows, cols = size(B)
        for λ = 1:L, j = 1:min(rows-λ, cols)
            view(R, j+λ, j) .= view(B, j+λ, j)
        end
        for i = 1:min(rows, cols)
            view(R, i, i) .= view(B, i, i)
        end
        for μ = 1:M, k = 1:min(rows, cols-μ)
            view(R, k, k+μ) .= view(B, k, k+μ)
        end
    end
    ret
end

BandedBlockBandedMatrix{T}(A::AbstractMatrix) where T =
    copyto!(BandedBlockBandedMatrix{T}(undef, axes(A), blockbandwidths(A), subblockbandwidths(A)), A)

BandedBlockBandedMatrix(A::AbstractMatrix{T}) where T = BandedBlockBandedMatrix{T}(A)

BandedBlockBandedMatrix(A::AbstractMatrix, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix(A, (blockedrange(rdims), blockedrange(cdims)), lu, λμ)

similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))

@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{BlockedUnitRange,AbstractUnitRange{Int}}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))
@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))
@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{BlockedUnitRange,BlockedUnitRange}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))

similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
    similar(Matrix{T}, map(length,axes)...)

axes(A::BandedBlockBandedMatrix) = (A.raxis, axes(A.data,2))

function sparse(A::BandedBlockBandedMatrix)
    i = Vector{Int}()
    j = Vector{Int}()
    z = Vector{eltype(A)}()
    for J = blockaxes(A,2), K = blockcolsupport(A, J)
        B = view(A, K, J)
        ĩ = _banded_rowval(B)
        j̃ = _banded_colval(B)
        z̃ = _banded_nzval(B)
        ĩ .+= first(axes(A,1)[K])-1
        j̃ .+= first(axes(A,2)[J])-1
        append!(i, ĩ)
        append!(j, j̃)
        append!(z, z̃)
    end
    sparse(i, j, z, size(A)...)
end

################################
# BandedBlockBandedMatrix Interface #
################################

bandedblockbandedcolumns(L::AbstractColumnMajor) = BandedBlockBandedColumnMajor()
bandedblockbandedcolumns(_) = BandedBlockBandedColumns{UnknownLayout}()

MemoryLayout(::Type{<:BandedBlockBandedMatrix{<:Any,BLOCKS}}) where BLOCKS =
    bandedblockbandedcolumns(MemoryLayout(BLOCKS))
bandedblockbandedbroadcaststyle(_) = BandedBlockBandedStyle()
BroadcastStyle(::Type{<:BandedBlockBandedMatrix{<:Any,BLOCKS}}) where BLOCKS =
    bandedblockbandedbroadcaststyle(BroadcastStyle(BLOCKS))

isbandedblockbanded(_) = false
isbandedblockbanded(::BandedBlockBandedMatrix) = true

blockbandwidths(A::BandedBlockBandedMatrix) = (A.l, A.u)

"""
    subblockbandwidths(A)

returns the sub-block bandwidths of `A`, where `A` is a banded-block-banded
matrix. In other words, `A[Block(K,J)]` will return a `BandedMatrix` with
bandwidths given by `subblockbandwidths(A)`.
"""
subblockbandwidths(A::BandedBlockBandedMatrix) = (A.λ, A.μ)

# default is to use bandwidth
subblockbandwidths(A::AbstractMatrix) = bandwidths(A)

"""
    subblockbandwidth(A, i)

returns the sub-block lower (`i == 1`) or upper (`i == 2`) bandwidth of `A`,
where `A` is a banded-block-banded matrix. In other words, `A[Block(K,J)]` will
return a `BandedMatrix` with the returned lower/upper bandwidth.
"""
subblockbandwidth(A::AbstractArray, k::Integer) = subblockbandwidths(A)[k]

isdiag(A::BandedBlockBandedMatrix) = A.λ == A.μ == A.l == A.u

################################
# AbstractBlockArray Interface #
################################

zeroblock(A::BandedBlockBandedMatrix, K::Int, J::Int) =
    BandedMatrix(Zeros{eltype(A)}(length.(getindex.(axes(A),(Block(K),Block(J))))), (A.λ, A.μ))

@inline function getblock(A::BandedBlockBandedMatrix, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)
    if -A.l ≤ J - K ≤ A.u
        convert(BandedMatrix, view(A, Block(K, J)))
    else
        zeroblock(A, K, J)
    end
end

# @inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
#     @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
#     @inbounds block = block_arr.blocks[blockindex.I...]
#     @boundscheck checkbounds(block, blockindex.α...)
#     @inbounds v = block[blockindex.α...]
#     return v
# end

###########################
# AbstractArray Interface #
###########################

# @inline function Base.similar(block_array::BandedBlockBandedMatrix{T}, ::Type{T2}) where {T,N,T2}
#     BlockArray(similar(block_array.blocks, Array{T2, N}), copy(block_array.block_sizes))
# end

Base.size(arr::BandedBlockBandedMatrix) =
    @inbounds return map(length,axes(arr))

@inline function getindex(A::BandedBlockBandedMatrix{T}, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = findblockindex.(axes(A), (i,j))
    @inbounds v = view(A, block.(bi)...)[blockindex.(bi)...]
    return v
end

@inline function setindex!(A::BandedBlockBandedMatrix{T}, v, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    BI,BJ = findblockindex.(axes(A), (i,j))
    if -A.l ≤ Int(block(BJ)-block(BI)) ≤ A.u
        V = view(A, block(BI),block(BJ))
         @inbounds V[blockindex(BI),blockindex(BJ)] = convert(T, v)::T
    elseif !iszero(v)
        throw(BandError(A))
    end
    return v
end

## structured matrix methods ##
function layout_replace_in_print_matrix(::AbstractBandedBlockBandedLayout, A, i, j, s)
    bi = findblockindex.(axes(A), (i,j))
    I,J = block.(bi)
    i,j = blockindex.(bi)
    l,u = blockbandwidths(A)
    λ,μ = subblockbandwidths(A)
    -l ≤ Int(J-I) ≤ u && -λ ≤ j-i ≤ μ ? s : Base.replace_with_centered_mark(s)
end

############
# Indexing #
############

# function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int}) where {T,N}
#     for i in 1:N
#         if size(v, i) != blocksize(block_arr.block_sizes, (i, block[i]))
#             throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
#         end
#     end
# end
#
#
@inline function setblock!(A::BandedBlockBandedMatrix, v, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)

    @boundscheck (bandwidth(v, 1) > A.λ || bandwidth(v, 2) > A.μ) && throw(BandError())
    V = view(A, Block(K), Block(J))
    V .= v
    return A
end
#
# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
#     getblock(block_array, block_index.I...)[block_index.α...] = v
# end

########
# Misc #
########

# @generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
#     # TODO: This will fail for empty block array
#     return quote
#         block_sizes = block_array.block_sizes
#         arr = similar(block_array.blocks[1], size(block_array)...)
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             arr[indices...] = getblock(block_array, block_index...)
#         end
#
#         return arr
#     end
# end
#
# @generated function Base.copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
#     return quote
#         block_sizes = block_array.block_sizes
#
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             copyto!(getblock(block_array, block_index...), arr[indices...])
#         end
#
#         return block_array
#     end
# end
#
# function Base.fill!(block_array::BlockArray, v)
#     for block in block_array.blocks
#         fill!(block, v)
#     end
# end

##################
# BandedBlockBandedBlock
#
#   views of the blocks satisfy the BandedMatrix interface, and are memory-compatible
#   with BLASBandedMatrix.
##################

const SubBandedBlockBandedMatrix{T,R1,R2} =
    SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}

sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()

sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedLayout()

sub_materialize(::AbstractBandedBlockBandedLayout, V, _) = BandedBlockBandedMatrix(V)

isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()

subblockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{BlockRange1}}}) =
    subblockbandwidths(parent(V))

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{Block1}}})
    A = parent(V)
    KR = parentindices(V)[1].block.indices[1]
    J = parentindices(V)[2].block
    shift = Int(KR[1])-Int(J)
    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:Block1},<:BlockSlice{<:BlockRange1}}})
    A = parent(V)
    K = parentindices(V)[1].block
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(K)-Int(JR[1])
    l,u = blockbandwidths(A)
    l - shift, u + shift
end

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1},<:BlockSlice{<:BlockRange1}}})
    A = parent(V)

    KR = parentindices(V)[1].block.indices[1]
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(first(KR))-Int(first(JR))

    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

const BandedBlockBandedBlock{T, BLOCKS, RAXIS} = SubArray{T,2,BandedBlockBandedMatrix{T, BLOCKS, RAXIS},<:Tuple{<:BlockSlice{Block1},<:BlockSlice{Block1}},false}

BroadcastStyle(::Type{<: BandedBlockBandedBlock}) = BandedStyle()

function inblockbands(V::SubArray{<:Any,2,<:AbstractMatrix,<:Tuple{<:BlockSlice{Block1},<:BlockSlice{Block1}},false})
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    l,u = blockbandwidths(A)
    -l ≤ Int(J-K) ≤ u
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block.block, J_sl.block.block)
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block, J_sl.block.block)
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block.block, J_sl.block)
end

# gives the columns of parent(V).data that encode the block
parentblocks2Int(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
                                                          Int(last(parentindices(V)).block)

######################################
# BandedMatrix interface  for Blocks #
######################################
@inline function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) where T
    inblockbands(V) && return subblockbandwidths(parent(V))
    (-720,-720)
end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
   B = parentblock(V)
   K_sl, J_sl = parentindices(V)
   bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],J_sl.block.indices[1])
end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    B = parentblock(V)
    K_sl, J_sl = parentindices(V)
    bandwidths(B) .+ (-1,1) .* bandshift(Base.OneTo(1),J_sl.block.indices[1])
 end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    B = parentblock(V)
    K_sl, J_sl = parentindices(V)
    bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],Base.OneTo(1))
 end

function bandeddata(V::BandedBlockBandedBlock{T}) where T
    inblockbands(V) || return Array{T}(undef, 0, size(V,2))
    A = parent(V)
    u = A.u
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    view(A.data, u + K - J + 1, J)
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    bandeddata(parentblock(V))
end

@inline function inbands_getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    A = parent(V)
    banded_getindex(bandeddata(V), A.λ, A.μ, k, j)
end

@inline function inbands_setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    A = parent(V)
    banded_setindex!(bandeddata(V), A.λ, A.μ, v, k, j)
end

# we need getindex as we use block-wise get index above
@propagate_inbounds function getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = parentblocks2Int(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_getindex(V, k, j)
    else
        zero(eltype(V))
    end
end

@propagate_inbounds function setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = parentblocks2Int(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_setindex!(V, v, k, j)
    elseif iszero(v) # allow setindex for 0 datya
        v
    else
        throw(BandError(V, J-K))
    end
end

BLAS.axpy!(a::T, A::BandedBlockBandedMatrix{T}, B::BandedBlockBandedMatrix{T}) where T =
    B .= a .* A .+ B

lmul!(x::Number, A::BandedBlockBandedMatrix) = (lmul!(x, A.data); A)
rmul!(A::BandedBlockBandedMatrix, x::Number) = (rmul!(A.data, x); A)
*(x::Number, A::BandedBlockBandedMatrix) = _BandedBlockBandedMatrix(x*A.data, axes(A,1), blockbandwidths(A), subblockbandwidths(A))
*(A::BandedBlockBandedMatrix, x::Number) = _BandedBlockBandedMatrix(A.data*x, axes(A,1), blockbandwidths(A), subblockbandwidths(A))


text:
using HTTP

export start_runtime

Base.@kwdef struct Invocation
  body::Any
  aws_request_id::String
  deadline_ms::Int
  invoked_function_arn::String
  trace_id::String
end

function get_endpoint(host::String)::String
  "http://$host/2018-06-01/runtime/invocation/"
end

function lambda_respond(response::String, endpoint::String, aws_request_id::String)
  HTTP.request(
    "POST",
    "$(endpoint)$(aws_request_id)/response",
    [],
    response,
  )
end

function lambda_error(error::String, endpoint::String, aws_request_id::String)
  HTTP.request(
    "POST",
    "$(endpoint)$(aws_request_id)/error",
    [("Lambda-Runtime-Function-Error-Type", "Unhandled")],
    JSON3.write(error),
  )
end

function start_runtime(host::String, func_name::String, param_type::String; single_shot=false)
  param_type = eval(Meta.parse(param_type))
  start_runtime(host, eval(Meta.parse(func_name)), param_type)
end

function start_runtime(host::String, react_function::Function, ::Type{T}; single_shot=false) where {T}
  endpoint = get_endpoint(host)
  println("Starting runtime at $endpoint")

  while true
    http = HTTP.request("GET", "$(endpoint)next"; verbose=3)
    body_raw = String(http.body)
    request_id = string(HTTP.header(http, "Lambda-Runtime-Aws-Request-Id"))

    body = try
      JSON3.read(body_raw, T)
    catch e
      body_sample = length(body_raw) > 50 ? "$(body_raw[start:50])..." : body_raw
      lambda_error("Unable to parse input JSON $body_sample", endpoint, request_id)
      continue
    end

    reaction = try
      react_function(body)
    catch e
      err(msg) = lambda_error(msg, endpoint, request_id)
      if isa(e, MethodError) && e.f == String(Symbol(react_function))
        err("react function is not a valid method for parsed JSON type")
      else
        err("$e")
      end
      continue
    end

    reaction_json = try
      JSON3.write(reaction)
    catch e
      err("Unable to parse function return value $(reaction) to JSON")
      continue
    end

    lambda_respond(reaction_json, endpoint, request_id)
    single_shot && break
  end
end


text:
using Distributions
using Base.Test
import Distributions: PoissonADSampler, PoissonCountSampler

function test_poissonsampler(s, μ::Float64, ns::Int, tol::Float64)
    ub = iceil(μ * 3) + 1
    pv = Distributions.poissonpvec(μ, ub)
    @assert length(pv) == ub + 1

    cnts = zeros(Int, ub+1)
    for i = 1:ns
        x = rand(s)
        @assert x >= 0
        if x <= ub
            cnts[x + 1] += 1
        end
    end

    q = cnts ./ ns
    @test_approx_eq_eps q pv tol
end

# test cases

for μ in [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0]
    test_poissonsampler(PoissonCountSampler(μ), μ, 10^5, 0.02)
end

for μ in [5.0, 10.0, 15.0, 20.0, 30.0]
    test_poissonsampler(PoissonADSampler(μ), μ, 10^5, 0.02)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using HTTP
using JSON
using SHA

# =========================================================================================
"""
    download_file(
        fromUrl::AbstractString,
        toFile::AbstractString
    )::String

Download the the URL and safe it to the file.  Returns the GitHub SHA1 value for the file
that was downloaded.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | fromUrl  |         | The URL to download.
| R | toFile   |         | The file to save the contents to.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| FileError    | There was an issue writing to the file.
| NetworkError | There was an issue downloading the file.

__Details:__

If the file exists it is overwritten.  The GitHub SHA1 is slightly different that a sha1sum
calculated locally.  First "blob $size\0" is hashed followed by the contents of the file.
"""
function download_file(
            fromUrl::AbstractString,
            toFile::AbstractString
        )::String
    local sha1 = SHA1_CTX()

    try
        # ---------------------------------------------------------------------------------
        # Open the file for writing, and make the HTTP request.
        open(toFile; create = true, write = true) do output

            HTTP.open("GET", fromUrl; redirect = false) do incoming

                # -------------------------------------------------------------------------
                # Start reading the data.
                startread(incoming)

                # -------------------------------------------------------------------------
                # Check for an error.
                if incoming.message.status != 200
                    @warn "Failed to download $fromUrl, HTTP error: $(incoming.message.status)"
                    throw(NetworkError(fromUrl))
                end

                # -------------------------------------------------------------------------
                # Start the hash.
                local size = parse(Int, HTTP.header(incoming, "Content-Length"))
                update!(sha1, Vector{UInt8}("blob $size\0"))

                # -------------------------------------------------------------------------
                # Download the file.
                local buffer = Vector{UInt8}(undef, 8192)

                while eof(incoming) == false
                    local bytes = readbytes!(incoming, buffer)
                    write(output, view(buffer, 1:bytes))
                    update!(sha1, buffer, bytes)
                end
            end
        end
    catch ex
        # ---------------------------------------------------------------------------------
        # We had an error nuke the file.
        if isfile(toFile)
            rm(toFile)
        end

        # ---------------------------------------------------------------------------------
        # SystemError means we had an file error.
        if isa(ex, SystemError)
            throw(FileError(toFile))
        end

        # ---------------------------------------------------------------------------------
        # Unknown error, rethrow it.
        rethrow(ex)
    end

    # -------------------------------------------------------------------------------------
    # Return the DataFile object.
    return bytes2hex(digest!(sha1))
end

# =========================================================================================
"""
    download_json(
        url::AbstractString,
    )::String

Download the the URL and parse the response as JSON.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | url  |         | The URL to download.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| NetworkError | There was an issue downloading the JSON object.
"""
function download_json(
            url::AbstractString
        )::Union{AbstractDict, AbstractArray}
    local response = HTTP.request("GET",  url)

    if response.status != 200
        @warn "Failed to download $url, HTTP error: $(response.status)"
        throw(NetworkError(url))
    end

    local raw = String(response.body)
    local json

    try
        json = JSON.parse(raw)
    catch ex
        if isa(ex, ErrorException)
            @warn "Unexpected response from the server: $raw"
            throw(NetworkError(url))
        end
        rethrow(ex)
    end

    return json
end


text:
# pick a nice default x range given a distribution
function default_range(dist::Distribution, alpha = 0.0001)
    minval = isfinite(minimum(dist)) ? minimum(dist) : quantile(dist, alpha)
    maxval = isfinite(maximum(dist)) ? maximum(dist) : quantile(dist, 1-alpha)
    minval, maxval
end

function default_range(m::Distributions.MixtureModel, alpha = 0.0001)
    minval = maxval = 0.0
    for c in m.components
        thismin = isfinite(minimum(c)) ? minimum(c) : quantile(c, alpha)
        thismax = isfinite(maximum(c)) ? maximum(c) : quantile(c, 1-alpha)
        if thismin < minval
            minval = thismin
        end
        if thismax > maxval
            maxval = thismax
        end
    end
    minval, maxval
end

yz_args(dist) = default_range(dist)
yz_args(dist::Distribution{N, T}) where N where T<:Discrete = (UnitRange(default_range(dist)...),)

# this "user recipe" adds a default x vector based on the distribution's μ and σ
@recipe function f(dist::Distribution)
    if dist isa Distribution{Univariate,Discrete}
        seriestype --> :scatterpath
    end
    (dist, yz_args(dist)...)
end

@recipe function f(m::Distributions.MixtureModel; components = true)
    if components
        for c in m.components
            @series begin
                (c, yz_args(c)...)
            end
        end
    else
        (m, yz_args(m)...)
    end
end

@recipe function f(distvec::AbstractArray{<:Distribution}, yz...)
    for di in distvec
        @series begin
            seriesargs = isempty(yz) ? yz_args(di) : yz
            if di isa Distribution{Univariate,Discrete}
                seriestype --> :scatterpath
            end
            (di, seriesargs...)
        end
    end
end

# this "type recipe" replaces any instance of a distribution with a function mapping xi to yi
@recipe function f(::Type{T}, dist::T; func = pdf) where T<:Distribution
    xi -> func(dist, xi)
end

#-----------------------------------------------------------------------------
# qqplots

@recipe function f(h::QQPair; qqline = :identity)
    if qqline in (:fit, :quantile, :identity, :R)
        xs = [extrema(h.qx)...]
        if qqline == :identity
            ys = xs
        elseif qqline == :fit
            itc, slp = hcat(fill!(similar(h.qx), 1), h.qx) \ h.qy
            ys = slp .* xs .+ itc
        else # if qqline == :quantile || qqline == :R
            quantx, quanty = quantile(h.qx, [0.25, 0.75]), quantile(h.qy, [0.25, 0.75])
            slp = diff(quanty) ./ diff(quantx)
            ys = quanty .+ slp .* (xs .- quantx)
        end

        @series begin
            primary := false
            seriestype := :path
            xs, ys
        end
    end

    seriestype --> :scatter
    legend --> false
    h.qx, h.qy
end

loc(D::Type{T}, x) where T<:Distribution = fit(D, x), x
loc(D, x) = D, x

@userplot QQPlot
recipetype(::Val{:qqplot}, args...) = QQPlot(args)
@recipe f(h::QQPlot) = qqbuild(loc(h.args[1], h.args[2])...)

@userplot QQNorm
recipetype(::Val{:qqnorm}, args...) = QQNorm(args)
@recipe f(h::QQNorm) = QQPlot((Normal, h.args[1]))


text:
abstract AbstractFilter{T<:Number,D<:AbstractDomain}

immutable Asymmetric1DFilter{T<:Number} <: AbstractFilter{T,SpatialDomain{1}}
    neg::Vector{T}
    pos::Vector{T}
    zero::T
end

immutable Symmetric1DFilter{T<:Number} <: AbstractFilter{T,SpatialDomain{1}}
    leg::Vector{T}
    zero::T
end

# product with scalar number is commutative
Base.:*(b::Number, ψ::AbstractFilter) = ψ * b

# element-wise multiplication operator ".*" with scalar falls back to "*"
Base.:.*(ψ::AbstractFilter, b::Number) = ψ * b

function spin!{T,D}(ψs::Array{AbstractFilter{T,D},3})
    (nΘs, nΧs, nJs) = size(ψs)
    for χ in 0:(nΧs-1), j in 0:(nJs-1)
        ψs[2:end, 1+χ, 1+j] = spin(ψs[1, 1+χ, 1+j])
    end
end


text:
function get_commandline_applescript(exe_dir, executable_path)
    """
    set RootPath to POSIX path of (path to me)
    tell application id "com.apple.terminal"
      do script ("exec '" & RootPath & "$exe_dir/$executable_path'")
      activate
    end tell
    """
end

# Build a wrapper app that opens a terminal and runs the provided binary.
# Returns a new applet name only if binary_name is already "applet" (to prevent collision).
function build_commandline_app_bundle(builddir, binary_name, appname, verbose)
    println("~~~~~~ Creating commandline-app wrapper applet. ~~~~~~~")

    mkpath(builddir)  # Create builddir if it doesn't already exist.
    applescript_file = "$builddir/build_$appname.applescript"
    exe_dir = "Contents/MacOS"  # Put the binaries next to the applet in MacOS.
    verbose && println("    Creating Applescript wrapper file: $applescript_file")
    write(applescript_file, get_commandline_applescript(exe_dir, binary_name))

    app_path = "$builddir/$appname.app"
    run(`osacompile -o $app_path $applescript_file`)
    applet_name = "applet"
    if binary_name == applet_name  # Prevent collisions.
        applet_name = "applet_wrapper"
        Compat.mv("$app_path/Contents/MacOS/applet", "$app_path/Contents/MacOS/$applet_name", force=true)
    end
    # Remove unneeded applet files it creates.
    rm("$app_path/Contents/Resources/applet.icns")

    # Cleaning up temporary file:
    verbose && println("    Cleaning up temporary Applescript file: $applescript_file")
    rm(applescript_file)

    return applet_name
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test
import Libdl

# these could fail on an embedded installation
# but for now, we don't handle that case
dlls = Libdl.dllist()
@test !isempty(dlls)
@test length(dlls) > 3 # at a bare minimum, probably have some version of libstdc, libgcc, libjulia, ...
if !Sys.iswindows() || Sys.windows_version() >= Sys.WINDOWS_VISTA_VER
    for dl in dlls
        if isfile(dl) && (Libdl.dlopen(dl; throw_error=false) !== nothing)
            @test Base.samefile(Libdl.dlpath(dl), dl)
        end
    end
end
@test length(filter(dlls) do dl
      if Base.DARWIN_FRAMEWORK
          return occursin(Regex("^$(Base.DARWIN_FRAMEWORK_NAME)(?:_debug)?\$"), basename(dl))
      else
          return occursin(Regex("^libjulia(?:.*)\\.$(Libdl.dlext)(?:\\..+)?\$"), basename(dl))
      end
    end) == 1 # look for something libjulia-like (but only one)

# library handle pointer must not be NULL
@test_throws ArgumentError Libdl.dlsym(C_NULL, :foo)
@test_throws ArgumentError Libdl.dlsym_e(C_NULL, :foo)

cd(@__DIR__) do

# Find the library directory by finding the path of libjulia (or libjulia-debug, as the case may be)
# and then adding on /julia to that directory path to get the private library directory, if we need
# to (where "need to" is defined as private_libdir/julia/libccalltest.dlext exists
private_libdir = if Base.DARWIN_FRAMEWORK
    if ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME * "_debug")))
    else
        joinpath(dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME))),"Frameworks")
    end
elseif ccall(:jl_is_debugbuild, Cint, ()) != 0
    dirname(abspath(Libdl.dlpath("libjulia-debug")))
else
    dirname(abspath(Libdl.dlpath("libjulia")))
end

if isfile(joinpath(private_libdir,"julia","libccalltest."*Libdl.dlext))
    private_libdir = joinpath(private_libdir, "julia")
end

@test !isempty(Libdl.find_library(["libccalltest"], [private_libdir]))
@test !isempty(Libdl.find_library("libccalltest", [private_libdir]))
@test !isempty(Libdl.find_library(:libccalltest, [private_libdir]))

# dlopen should be able to handle absolute and relative paths, with and without dlext
let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest.$(Libdl.dlext)")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(relpath(joinpath(private_libdir, "libccalltest")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(relpath(joinpath(private_libdir, "libccalltest.$(Libdl.dlext)")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen("./foo"; throw_error=false)
        @test dl === nothing
    finally
        Libdl.dlclose(dl)
    end
end

# unqualified names present in DL_LOAD_PATH
let dl = C_NULL
    try
        dl = Libdl.dlopen("libccalltest"; throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(string("libccalltest",".",Libdl.dlext); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

# path with dlopen-able file first in load path
#=
let dl = C_NULL,
    tmpdir = mktempdir(),
    fpath = joinpath(tmpdir,"libccalltest")
    try
        write(open(fpath,"w"))
        push!(Libdl.DL_LOAD_PATH, @__DIR__)
        push!(Libdl.DL_LOAD_PATH, dirname(fpath))
        dl = Libdl.dlopen_e("libccalltest")
        @test dl != C_NULL
    finally
        pop!(Libdl.DL_LOAD_PATH)
        pop!(Libdl.DL_LOAD_PATH)
        rm(tmpdir, recursive=true)
    end
end
=#

# path with dlopen-able file second in load path
#=
let dl = C_NULL,
    tmpdir = mktempdir(),
    fpath = joinpath(tmpdir,"libccalltest")
    try
        write(open(fpath,"w"))
        push!(Libdl.DL_LOAD_PATH, dirname(fpath))
        push!(Libdl.DL_LOAD_PATH, @__DIR__)
        dl = Libdl.dlopen_e("libccalltest")
        @test dl != C_NULL
    finally
        pop!(Libdl.DL_LOAD_PATH)
        pop!(Libdl.DL_LOAD_PATH)
        rm(tmpdir, recursive=true)
    end
end
=#

# test dlpath
let dl = C_NULL
    try
        path = abspath(joinpath(private_libdir, "libccalltest"))
        dl = Libdl.dlopen(path; throw_error=false)
        @test dl !== nothing
        @test Base.samefile(abspath(Libdl.dlpath(dl)),
                            abspath(Libdl.dlpath(path)))
        @test Base.samefile(abspath(Libdl.dlpath(dl)),
                            string(path,".",Libdl.dlext))
    finally
        Libdl.dlclose(dl)
    end
end

# opening a library that does not exist throws an ErrorException
@test_throws ErrorException Libdl.dlopen("./foo")

# opening a versioned library that does not exist does not result in adding extension twice
err = @test_throws ErrorException Libdl.dlopen("./foo.$(Libdl.dlext).0")
@test !occursin("foo.$(Libdl.dlext).0.$(Libdl.dlext)", err.value.msg)
err = @test_throws ErrorException Libdl.dlopen("./foo.$(Libdl.dlext).0.22.1")
@test !occursin("foo.$(Libdl.dlext).0.22.1.$(Libdl.dlext)", err.value.msg)

# test dlsym
let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest")))
        fptr = Libdl.dlsym(dl, :set_verbose)
        @test fptr !== nothing
        @test_throws ErrorException Libdl.dlsym(dl, :foo)

        fptr = Libdl.dlsym_e(dl, :set_verbose)
        @test fptr != C_NULL
        fptr = Libdl.dlsym_e(dl, :foo)
        @test fptr == C_NULL
    finally
        Libdl.dlclose(dl)
    end
end

# test dlclose
# If dl is NULL, jl_dlclose should return -1 and dlclose should return false
# dlclose should return true on success and false on failure
let dl = C_NULL
    @test -1 == ccall(:jl_dlclose, Cint, (Ptr{Cvoid},), dl)
    @test !Libdl.dlclose(dl)

    dl = Libdl.dlopen("libccalltest"; throw_error=false)
    @test dl !== nothing

    @test Libdl.dlclose(dl)
    @test_skip !Libdl.dlclose(dl)   # Syscall doesn't fail on Win32
end

end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
function g(x) x = 42 end
function f(x) x = 4; g(x); println(x) end
f(3)

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:
#=
    Julia is a bit even more snappy when indicies, bitshifts and other register-centric ops
        are performed with Int16 or Int32 values.  Int64 values are very fast, too.
     The functions I use to optimize the type of the indexing value are needed
        when they are being leaned upon to provide time-minimized access.
=#

@inline bitsof(::Type{T}) where T = sizeof(T) * 8
bitsof(x::T) where T = bitsof(T)

function bits_required(posint::T)::T where T<:Integer
    posint > 0 || throw(ErrorException("posint ($(posint)) must be > 0"))
    return bitsof(T) - leading_zeros(posint)%T
end

const TYPES_FOR_INDEXING = [ Int8, Int16, Int32, Int64, Int128 ]

function type_for_indexing(posint::T) where T<:Union{Int8, Int16, Int32, Int64, Int128}
    posint < 0 && throw(ErrorException("posint ($(posint)) must be => 0"))
    posint = max(127%T, posint)
    type_index = nextpow2(bits_required(posint) + 1) >>> 2
    type_index = trailing_zeros(type_index)
    type_index = max(1, min(5, type_index))
    return TYPES_FOR_INDEXING[ type_index ]
end


text:
using PyPlot
include("style_colorscheme.jl")
using3D()

DrWatson._wsave(s, fig::Figure) = fig.savefig(s, dpi = 600, transparent = false)

"""
    rdspl(x, n = 2)
Round `x` to `n` sigdigits for display purposes.
"""
rdspl(x::Real, n = 2) = round(x, digits=n)
rdspl(x::AbstractVector, n = 2) = Tuple((round.(Float64.(x); sigdigits=n)))

PyPlot.rc("lines", lw = 3.0)
PyPlot.rc("errorbar", capsize = 6)
PyPlot.rc("axes", grid = true)
PyPlot.rc("grid", color = "0.75", alpha = 0.75)
PyPlot.rc("axes", axisbelow=true) # makes grid behind other plotting elements

PyPlot.rc("font", size = 24) # set default fontsize
PyPlot.rc("xtick", labelsize = 24)
PyPlot.rc("ytick", labelsize = 24)
PyPlot.rc("axes", labelsize = 28)
PyPlot.rc("legend";
    fontsize = 24, handletextpad = 0.5, handlelength = 2,
    labelspacing = 0.25, columnspacing = 1.0,
)
# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

for z in ("x", "y")
    PyPlot.rc("$(z)tick.major", size = 7, width = 1.2)
    PyPlot.rc("$(z)tick.minor", size = 3, visible = false)
end

figx = 12 # default width correspoding to full text width
figy = 9  # default height corresponding to 1 row of plots
PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("savefig", dpi = 600, transparent = true, format = "png")

# set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORSCHEME))

if false # test color scheme
    fig = figure(figsize = (20, 15)) # show colors
    ax1 = subplot(231)
    ax2 = subplot(232)
    ax3 = subplot(233)
    ax4 = subplot(223)
    ax5 = subplot(224)
    lw = 60
    L = length(COLORSCHEME)
    for (i, c) in enumerate(COLORS)
        chsv = matplotlib.colors.rgb_to_hsv(matplotlib.colors.to_rgb(c))
        ax1.plot([0, 1], [0, 0] .+ i, color = c, lw = lw)
        ax1.set_title("color")
        ax2.plot([0, 1], [0, 0] .+ i, color = string(chsv[3]), lw = lw)
        ax2.set_title("brightness")
        ax3.plot([0, 1], [0, 0] .+ i, color = string(chsv[2]), lw = lw)
        ax3.set_title("saturation")
        x = 0:0.05:5π
        ax4.plot(x, cos.(x .+ i/2) .+ rand(length(x))/2; color=c, lw = 2)
        ax5.bar(collect(1:4) .+ (i-1)/L, 0.5rand(4) .+ 0.5, 1/L; color=c)
    end
    fig = tight_layout()
end

bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)

"`add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc = 0.985, yloc = 0.975)`"
function add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc = 0.975, yloc = 1)
    bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)
    for (i, ax) in enumerate(axs)
        l = ('a':'z')[i]
        try
            ax.text(xloc, yloc, "$(l)"; transform = ax.transAxes,
            bbox = bbox, zorder = 99, va = "top")
        catch e
            @warn string(e)
        end
    end
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
#=
https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/04.%20Loops.ipynb
=#

#=
let's use for loops to create some addition tables,
where the value of every entry is the sum of its row and column indices.
=#

function create_table1(m,n)
    A = fill(0, (m, n))
    for i in 1:m
        for j in 1:n
            A[i,j] = i+j
        end
    end
    return A
end

#=
Here is some syntactic sugar for the same nested `for` loop
=#

function create_table2(m,n)
    B = fill(0, (m, n))
    for i in 1:m, j in 1:n
        B[i,j] = i+j
    end
    return B
end

#=
The more "Julia" way to create this addition table would have been
with an array comprehension.
=#

function create_table3(m,n)
    C = [i+j for i in 1:m, j in 1:n]
    return C
end

@show create_table1(4,5)
@show create_table2(4,5)
@show create_table3(4,5)

using InteractiveUtils
versioninfo()

using BenchmarkTools
m,n=1000,1000
@btime create_table1(m,n)
@btime create_table2(m,n)
@btime create_table3(m,n)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
cd("lesson-packages")
using Pkg
Pkg.activate("TestPackage1")

#%%
using TestPackage1

TestPackage1.greet()
#%%
mySpecialFunction(3)
#%%
using Pkg
Pkg.instantiate()

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
struct ∇{T}
  operand::T
  string::String
end

struct ∇o{T}
  operand::T
  string::String
end

struct ∇x{T}
  operand::T
  string::String
end

String(x::∇o) = "∇o" * x.string
String(x::∇x) = "∇×" * x.string

(g::∇)(x) = CovariantVector(x->grad(g.operand, x), true)
(c::∇x)(x) = ContravariantVector(x->curl(c.operand, x), true)
(d::∇o)(v::VectorField) = Scalar(dot(d.operand, v))
(d::∇o)(c::∇x) = 0
(c::∇x)(g::∇) = zeros(3)

Base.show(io::IO, x::∇o) = print(io, String(x))
Base.show(io::IO, x::∇x) = print(io, String(x))

max_stars_repo_name: jwscook/DifferentialGeometries

max_issues_repo_name: jwscook/DifferentialGeometries

max_forks_repo_name: jwscook/DifferentialGeometries


text:
function save_features_instance!(db::SQLite.DB, name, scenario, source_type, source_path)
    println("Saving instance features: ($name, $scenario)")
    network = PowerFlowNetwork(source_path, source_type)
    merge_duplicate_branch!(network)
    features = get_features_instance(network)
    query = "UPDATE instances SET "
    for (feature_name, feature_value) in features
        if feature_value == Inf
            query *= "$feature_name = '+Infinity', "
        elseif feature_value == -Inf
            query *= "$feature_name = '-Infinity', "
        elseif isnan(feature_value)
            query *= "$feature_name = 'NaN', "
        else
            query *= "$feature_name = $feature_value, "
        end
    end
    query = query[begin:end - 2]
    query *= " WHERE name = '$name' AND scenario = $scenario"
    execute_query(db, query)
end

function save_features_instance_dfrow!(db::SQLite.DB, row)
    save_features_instance!(db, row[:name], row[:scenario], row[:source_type], row[:source_path])
end

function save_features_instances!(db::SQLite.DB; min_nv=typemin(Int), max_nv=typemax(Int), recompute=false, subset=nothing)
    query = "SELECT name, scenario, source_type, source_path FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !recompute
        query *= " AND nb_edge IS NULL"
    end
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Saving instance features config: min_nv=$min_nv, max_nv=$max_nv, recompute=$recompute")
    save_function!(row) = save_features_instance_dfrow!(db, row)
    save_function!.(eachrow(results[!, [:name, :scenario, :source_type, :source_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using CUDA
using OrderedCollections: OrderedDict

using Oceananigans.Advection
using Oceananigans: AbstractOutputWriter, AbstractDiagnostic, TimeStepper
using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Buoyancy: validate_buoyancy
using Oceananigans.TurbulenceClosures: ν₀, κ₀, with_tracers

mutable struct IncompressibleModel{TS, E, A<:AbstractArchitecture, G, T, B, R, SW, U, C, Φ, F,
                                   V, S, K} <: AbstractModel
       architecture :: A         # Computer `Architecture` on which `Model` is run
               grid :: G         # Grid of physical points on which `Model` is solved
              clock :: Clock{T}  # Tracks iteration number and simulation time of `Model`
          advection :: V         # Advection scheme for velocities _and_ tracers
           buoyancy :: B         # Set of parameters for buoyancy model
           coriolis :: R         # Set of parameters for the background rotation rate of `Model`
      surface_waves :: SW        # Set of parameters for surfaces waves via the Craik-Leibovich approximation
            forcing :: F         # Container for forcing functions defined by the user
            closure :: E         # Diffusive 'turbulence closure' for all model fields
         velocities :: U         # Container for velocity fields `u`, `v`, and `w`
            tracers :: C         # Container for tracer fields
          pressures :: Φ         # Container for hydrostatic and nonhydrostatic pressure
      diffusivities :: K         # Container for turbulent diffusivities
        timestepper :: TS        # Object containing timestepper fields and parameters
    pressure_solver :: S         # Pressure/Poisson solver
end

"""
    IncompressibleModel(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                  clock = Clock{float_type}(0, 0, 1),
              advection = CenteredSecondOrder(),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
                closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
    boundary_conditions = (u=UVelocityBoundaryConditions(grid),
                           v=VVelocityBoundaryConditions(grid),
                           w=WVelocityBoundaryConditions(grid)),
                tracers = (:T, :S),
            timestepper = :QuasiAdamsBashforth2,
             velocities = VelocityFields(architecture, grid, boundary_conditions),
              pressures = PressureFields(architecture, grid, boundary_conditions),
          diffusivities = DiffusivityFields(architecture, grid, tracernames(tracers), boundary_conditions, closure),
        pressure_solver = PressureSolver(architecture, grid, PressureBoundaryConditions(grid))
    )

Construct an incompressible `Oceananigans.jl` model on `grid`.

Keyword arguments
=================

    - `grid`: (required) The resolution and discrete geometry on which `model` is solved.
    - `architecture`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
    - `float_type`: `Float32` or `Float64`. The floating point type used for `model` data.
    - `advection`: The scheme that advects velocities and tracers. See `Oceananigans.Advection`.
    - `buoyancy`: The buoyancy model. See `Oceananigans.Buoyancy`.
    - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
    - `coriolis`: Parameters for the background rotation rate of the model.
    - `forcing`: User-defined forcing functions that contribute to solution tendencies.
    - `boundary_conditions`: `NamedTuple` containing field boundary conditions.
    - `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
                 preallocated `CellField`s.
    - `timestepper`: A symbol that species the time-stepping method. Either `:QuasiAdamsBashforth2` or
                     `:RungeKutta3`.
"""
function IncompressibleModel(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                  clock = Clock{float_type}(0, 0, 1),
              advection = CenteredSecondOrder(),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
                closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
    boundary_conditions = (u=UVelocityBoundaryConditions(grid),
                           v=VVelocityBoundaryConditions(grid),
                           w=WVelocityBoundaryConditions(grid)),
                tracers = (:T, :S),
            timestepper = :QuasiAdamsBashforth2,
             velocities = VelocityFields(architecture, grid, boundary_conditions),
              pressures = PressureFields(architecture, grid, boundary_conditions),
          diffusivities = DiffusivityFields(architecture, grid, tracernames(tracers), boundary_conditions, closure),
        pressure_solver = PressureSolver(architecture, grid, PressureBoundaryConditions(grid))
    )

    if architecture == GPU() && !has_cuda()
         throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
    end

    validate_buoyancy(buoyancy, tracernames(tracers))

    # Regularize forcing and closure for given tracer fields.
    forcing = ModelForcing(tracernames(tracers), forcing)
    closure = with_tracers(tracernames(tracers), closure)

    # Instantiate tracer fields if not already instantiated
    tracer_fields = TracerFields(architecture, grid, tracers, boundary_conditions)

    # Instantiate timestepper if not already instantiated
    timestepper = TimeStepper(timestepper, architecture, grid, tracernames(tracers))

    return IncompressibleModel(architecture, grid, clock, advection, buoyancy, coriolis, surface_waves,
                               forcing, closure, velocities, tracer_fields, pressures, diffusivities,
                               timestepper, pressure_solver)
end


text:
#---------------------------------------------------------
# # [SPECTrecon PSF](@id 3-psf)
#---------------------------------------------------------

# This page explains the PSF portion of the Julia package
# [`SPECTrecon.jl`](https://github.com/JeffFessler/SPECTrecon.jl).

# ### Setup

# Packages needed here.

using SPECTrecon
using MIRTjim: jim, prompt
using Plots: scatter, scatter!, plot!, default
default(markerstrokecolor=:auto, markersize=3)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# After rotating the image and the attenuation map,
# second step in SPECT image forward projection
# is to apply depth-dependent point spread function (PSF).
# Each (rotated) image plane
# is a certain distance from the SPECT detector
# and must be convolved with the 2D PSF appropriate
# for that plane.

# Because SPECT has relatively poor spatial resolution,
# the PSF is usually fairly wide,
# so convolution using FFT operations
# is typically more efficient
# than direct spatial convolution.

# Following other libraries like
# [FFTW.jl](https://github.com/JuliaMath/FFTW.jl),
# the PSF operations herein start with a `plan`
# where work arrays are preallocated
# for subsequent use.
# The `plan` is a `Vector` of `PlanPSF` objects:
# one for each thread.
# (Parallelism is across planes for a 3D image volume.)
# The number of threads defaults to `Threads.nthreads()`.

# ### Example

# Start with a 3D image volume.

T = Float32 # work with single precision to save memory
nx = 32
nz = 30
image = zeros(T, nx, nx, nz) # ny = nx required
image[1nx÷4, 1nx÷4, 3nz÷4] = 1
image[2nx÷4, 2nx÷4, 2nz÷4] = 1
image[3nx÷4, 3nx÷4, 1nz÷4] = 1
jim(image, "Original image")

# Create a synthetic gaussian depth-dependent PSF for a single view

px = 11
nview = 1 # for simplicity in this illustration
psf = repeat(psf_gauss( ; ny, px), 1, 1, 1, nview)
jim(psf, "PSF for each of $nx planes")

# Now plan the PSF modeling
# by specifying
# * the image size (must be square)
# * the PSF size: must be `px × pz × ny × nview`
# * the `DataType` used for the work arrays.

plan = plan_psf( ; nx, nz, px, T)

# Here are the internals for the plan for the first thread:

plan[1]

# With this `plan` pre-allocated, now we can apply the depth-dependent PSF
# to the image volume (assumed already rotated here).

result = similar(image) # allocate memory for the result
fft_conv!(result, image, psf[:,:,:,1], plan) # mutates the first argument
jim(result, "After applying PSF")

# ### Adjoint

# To ensure adjoint consistency between SPECT forward- and back-projection,
# there is also an adjoint routine:

adj = similar(result)
fft_conv_adj!(adj, result, psf[:,:,:,1], plan)
jim(adj, "Adjoint of PSF modeling")

# The adjoint is *not* the same as the inverse
# so one does not expect the output here to match the original image!

# ### LinearMap

# One can form a linear map corresponding to PSF modeling using `LinearMapAA`.
# Perhaps the main purpose is simply for verifying adjoint correctness.

using LinearMapsAA: LinearMapAA

nx, nz, px = 10, 7, 5 # small size for illustration
psf3 = psf_gauss( ; ny, px)
plan = plan_psf( ; nx, nz, px, T)
idim = (nx,nx,nz)
odim = (nx,nx,nz)
forw! = (y,x) -> fft_conv!(y, x, psf3, plan)
back! = (x,y) -> fft_conv_adj!(x, y, psf3, plan)
A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)

Afull = Matrix(A)
Aadj = Matrix(A')
jim(cat(dims=3, Afull, Aadj'), "Linear map for PSF modeling and its adjoint")

# The following check verifies adjoint consistency:

@assert Afull ≈ Aadj'