text:
"""
Lazy representation of a scaled map `λ * A = A * λ` with real or complex map
`A <: LinearMap{RealOrComplex}` and real or complex scaling factor
`λ <: RealOrComplex`.
"""
struct ScaledMap{T, S<:RealOrComplex, L<:LinearMap} <: LinearMap{T}
    λ::S
    lmap::L
    function ScaledMap{T}(λ::S, A::L) where {T, S <: RealOrComplex, L <: LinearMap{<:RealOrComplex}}
        @assert Base.promote_op(*, S, eltype(A)) == T "target type $T cannot hold products of $S and $(eltype(A)) objects"
        new{T,S,L}(λ, A)
    end
end

# constructor
ScaledMap(λ::RealOrComplex, lmap::LinearMap{<:RealOrComplex}) =
    ScaledMap{Base.promote_op(*, typeof(λ), eltype(lmap))}(λ, lmap)

# basic methods
Base.size(A::ScaledMap) = size(A.lmap)
Base.axes(A::ScaledMap) = axes(A.lmap)
Base.isreal(A::ScaledMap) = isreal(A.λ) && isreal(A.lmap)
LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap)
LinearAlgebra.ishermitian(A::ScaledMap) = isreal(A.λ) && ishermitian(A.lmap)
LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap)

Base.transpose(A::ScaledMap) = A.λ * transpose(A.lmap)
Base.adjoint(A::ScaledMap) = conj(A.λ) * adjoint(A.lmap)

# comparison (sufficient, not necessary)
Base.:(==)(A::ScaledMap, B::ScaledMap) =
    eltype(A) == eltype(B) && A.lmap == B.lmap && A.λ == B.λ

# scalar multiplication and division
Base.:(*)(α::RealOrComplex, A::LinearMap{<:RealOrComplex}) = ScaledMap(α, A)
Base.:(*)(A::LinearMap{<:RealOrComplex}, α::RealOrComplex) = ScaledMap(α, A)

Base.:(*)(α::Number, A::ScaledMap) = (α * A.λ) * A.lmap
Base.:(*)(A::ScaledMap, α::Number) = A.lmap * (A.λ * α)
# needed for disambiguation
Base.:(*)(α::RealOrComplex, A::ScaledMap) = (α * A.λ) * A.lmap
Base.:(*)(A::ScaledMap, α::RealOrComplex) = (A.λ * α) * A.lmap
Base.:(-)(A::LinearMap) = -1 * A

# composition (not essential, but might save multiple scaling operations)
Base.:(*)(A::ScaledMap, B::ScaledMap) = (A.λ * B.λ) * (A.lmap * B.lmap)
Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ

# multiplication with vectors/matrices
for (In, Out) in ((AbstractVector, AbstractVecOrMat),
                  (AbstractMatrix, AbstractMatrix))
    @eval begin
        # commutative case
        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In{<:RealOrComplex})
            return _unsafe_mul!(y, A.lmap, x, A.λ, false)
        end
        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In{<:RealOrComplex}, α::Number, β::Number)
            return _unsafe_mul!(y, A.lmap, x, A.λ * α, β)
        end
        # non-commutative case
        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In)
            return lmul!(A.λ, _unsafe_mul!(y, A.lmap, x))
        end
    end
end


text:
"""
Tableau of one-stage, 1st order explicit (forward) Euler method

```julia
TableauExplicitEuler(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Leonhard Euler.
    Institutiones calculi differentialis cum eius vsu in analysi finitorum ac doctrina serierum.
    Imp. Acad. Imper. Scient. Petropolitanae, Opera Omnia, Vol.X, [I.6], 1755.
    In: Opera Omnia, 1st Series, Volume 11, Institutiones Calculi Integralis. Teubner, Leipzig, Pages 424-434, 1913.
    Sectio secunda. Caput VII. De integratione aequationum differentialium per approximationem. Problema 85.

"""
function TableauExplicitEuler(::Type{T}=Float64) where {T}
    a = zeros(BigFloat, 1, 1)
    b = ones(BigFloat, 1)
    c = zeros(BigFloat, 1)
    o = 1

    Tableau{T}(:explicit_euler, o, a, b, c)
end

"Alias for [`TableauExplicitEuler`](@ref)"
TableauForwardEuler = TableauExplicitEuler

"""
Tableau of explicit two-stage, 2nd order midpoint method

```julia
TableauExplicitMidpoint(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Carl Runge.
    Über die numerische Auflösung von Differentialgleichungen.
    Mathematische Annalen, Volume 46, Pages 167-178, 1895.
    doi: 10.1007/BF01446807.
    Equation (2)

"""
function TableauExplicitMidpoint(::Type{T}=Float64) where {T}
    a = @big [[ 0     0    ]
              [ 1//2  0    ]]
    b = @big  [ 0,    1    ]
    c = @big  [ 0,    1//2 ]
    o = 2

    Tableau{T}(:explicit_midpoint, o, a, b, c)
end

"""
Tableau of Heun's two-stage, 2nd order method

```julia
TableauHeun2(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Karl Heun.
    Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen.
    Zeitschrift für Mathematik und Physik, Volume 45, Pages 23-38, 1900.
    Algorithm II.

"""
function TableauHeun2(::Type{T}=Float64) where {T}
    a = @big [[ 0      0    ]
              [ 1      0    ]]
    b = @big  [ 1//2,  1//2 ]
    c = @big  [ 0,     1    ]
    o = 2

    Tableau{T}(:heun2, o, a, b, c)
end

"""
Tableau of Heun's three-stage, 3rd order method

```julia
TableauHeun3(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Karl Heun.
    Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen.
    Zeitschrift für Mathematik und Physik, Volume 45, Pages 23-38, 1900.
    Algorithm VI.

"""
function TableauHeun3(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0    ]
              [ 1//3   0      0    ]
              [ 0      2//3   0    ]]
    b = @big  [ 1//4,  0,     3//4 ]
    c = @big  [ 0,     1//3,  2//3 ]
    o = 3

    Tableau{T}(:heun3, o, a, b, c)
end

"""
Tableau of Ralston's two-stage, 2nd order method

```julia
TableauRalston2(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Anthony Ralston.
    Runge-Kutta Methods with Minimum Error Bounds.
    Mathematics of Computation, Volume 16, Pages 431-437, 1962.
    doi: 10.1090/S0025-5718-1962-0150954-0.
    Equation (3.5)

"""
function TableauRalston2(::Type{T}=Float64) where {T}
    a = @big [[ 0      0    ]
              [ 2//3   0    ]]
    b = @big  [ 1//4,  3//4 ]
    c = @big  [ 0,     2//3 ]
    o = 2

    Tableau{T}(:ralston2, o, a, b, c)
end

"""
Tableau of Ralston's three-stage, 3rd order method

```julia
TableauRalston3(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Anthony Ralston.
    Runge-Kutta Methods with Minimum Error Bounds.
    Mathematics of Computation, Volume 16, Pages 431-437, 1962.
    doi: 10.1090/S0025-5718-1962-0150954-0.
    Equation (4.10)

"""
function TableauRalston3(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0    ]
              [ 1//2   0      0    ]
              [ 0      3//4   0    ]]
    b = @big  [ 2//9,  3//9,  4//9 ]
    c = @big  [ 0,     1//2,  3//4 ]
    o = 3

    Tableau{T}(:ralston3, o, a, b, c)
end

"""
Tableau of Runge's two-stage, 2nd order method

```julia
TableauRunge(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Carl Runge
    Über die numerische Auflösung von Differentialgleichungen.
    Mathematische Annalen, Volume 46, Pages 167-178, 1895.
    doi: 10.1007/BF01446807.
    Equation (3)

"""
function TableauRunge(::Type{T}=Float64) where {T}
    a = @big [[ 0      0    ]
              [ 1      0    ]]
    b = @big  [ 1//2,  1//2 ]
    c = @big  [ 0,     1    ]
    o = 2

    Tableau{T}(:runge, o, a, b, c)
end

"Alias for [`TableauRunge`](@ref)"
TableauRunge2 = TableauRunge

"""
Tableau of Kutta's three-stage, 3rd order method

```julia
TableauKutta(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Wilhelm Kutta
    Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
    Zeitschrift für Mathematik und Physik, Volume 46, Pages 435–453, 1901.
    Page 440

"""
function TableauKutta(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0    ]
              [ 1//2   0      0    ]
              [-1      2      0    ]]
    b = @big  [ 1//6,  4//6,  1//6 ]
    c = @big  [ 0,     1//2,  1    ]
    o = 3

    Tableau{T}(:kutta, o, a, b, c)
end

"Alias for [`TableauKutta`](@ref)"
TableauKutta3 = TableauKutta

"""
Tableau of explicit Runge-Kutta method of order four (1/6 rule)

```julia
TableauRK416(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Wilhelm Kutta
    Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
    Zeitschrift für Mathematik und Physik, Volume 46, Pages 435–453, 1901.
    Page 443

"""
function TableauRK416(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0      0    ]
              [ 1//2   0      0      0    ]
              [ 0      1//2   0      0    ]
              [ 0      0      1      0    ]]
    b = @big  [ 1//6,  1//3,  1//3,  1//6 ]
    c = @big  [ 0,     1//2,  1//2,  1    ]
    o = 4

    Tableau{T}(:erk4, o, a, b, c)
end

"Alias for [`TableauRK416`](@ref)"
TableauRK4 = TableauRK416

"""
Tableau of explicit Runge-Kutta method of order four (3/8 rule)

```julia
TableauRK438(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Wilhelm Kutta
    Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
    Zeitschrift für Mathematik und Physik, Volume 46, Pages 435–453, 1901.
    Page 441

"""
function TableauRK438(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0      0    ]
              [ 1//3   0      0      0    ]
              [-1//3   1      0      0    ]
              [ 1     -1      1      0    ]]
    b = @big  [ 1//8,  3//8,  3//8,  1//8 ]
    c = @big  [ 0,     1//3,  2//3,  1    ]
    o = 4

    Tableau{T}(:erk438, o, a, b, c)
end

"""
Tableau of 2rd order Strong Stability Preserving method with two stages and CFL ≤ 1

Alias for [`TableauHeun2`](@ref)

```julia
TableauSSPRK2(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Chi-Wang Shu, Stanley Osher.
    Efficient implementation of essentially non-oscillatory shock-capturing schemes.
    Journal of Computational Physics, Volume 77, Issue 2, Pages 439-471, 1988.
    doi: 10.1016/0021-9991(88)90177-5.
    Equation (2.16)

"""
TableauSSPRK2 = TableauHeun2

"""
Tableau of 3rd order Strong Stability Preserving method with three stages and CFL ≤ 1

```julia
TableauSSPRK3(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    Chi-Wang Shu, Stanley Osher.
    Efficient implementation of essentially non-oscillatory shock-capturing schemes.
    Journal of Computational Physics, Volume 77, Issue 2, Pages 439-471, 1988.
    doi: 10.1016/0021-9991(88)90177-5.
    Equation (2.18)

"""
function TableauSSPRK3(::Type{T}=Float64) where {T}
    a = @big [[ 0      0      0    ]
              [ 1      0      0    ]
              [ 1//4   1//4   0    ]]
    b = @big  [ 1//6,  1//6,  4//6 ]
    c = @big  [ 0,     1,     1//2 ]
    o = 3

    Tableau{T}(:ssprk3, o, a, b, c)
end

ext: jl


text:
export TFI

struct TFI{D} <: AbstractModel{D}
    size::NTuple{D, Int}
    h::Float64
    periodic::Bool
    TFI(size::Int...; h::Real, periodic::Bool) = new{length(size)}(size, Float64(h), periodic)
end

function get_bonds(model::TFI{1})
    nbit, = model.size
    [(i, i%nbit+1, 1.0) for i in 1:(model.periodic ? nbit : nbit-1)]
end

Base.size(model::TFI) = model.size

function hamiltonian(model::TFI{1})
    model.periodic && throw()
    nbit = nspin(model)
    sum(repeat(nbit, Z, (i,i+1)) for i=1:nbit-1)*0.25 +
    sum(put(nbit, i=>X) for i=1:nbit)*0.5model.h
end

function energy(chem::QuantumMPS, model::TFI)
   energy(chem, Z, model) +
    energy(chem, X, model)
end

function energy(chem::QuantumMPS, pauli::ZGate, model::TFI)
    res = gensample(chem, pauli)
    local eng = 0.0
    for bond in ((i, i+1, 1.0) for i=1:nspin(model)-1)
        eng += bond[3]*mean(res[:,bond[1]].*res[:,bond[2]])
    end
    eng/4
end

function energy(chem::QuantumMPS, pauli::XGate, model::TFI)
    res = gensample(chem, pauli)
    eng = mean(sum(res, dims=2))
    model.h*eng/2
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
function add_constraint_dual!(
    container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{T, D},
) where {T <: PSY.Component, D <: AbstractDeviceFormulation}
    if !isempty(get_duals(model))
        devices = get_available_components(T, sys)
        for constraint_type in get_duals(model)
            assign_dual_variable!(container, constraint_type, devices, D)
        end
    end
    return
end

function add_constraint_dual!(
    container::OptimizationContainer,
    sys::PSY.System,
    model::NetworkModel{T},
) where {T <: PM.AbstractPowerModel}
    if !isempty(get_duals(model))
        devices = PSY.get_components(PSY.Bus, sys)
        for constraint_type in get_duals(model)
            assign_dual_variable!(container, constraint_type, devices, T)
        end
    end
    return
end

function add_constraint_dual!(
    container::OptimizationContainer,
    sys::PSY.System,
    model::NetworkModel{T},
) where {T <: Union{CopperPlatePowerModel, StandardPTDFModel}}
    if !isempty(get_duals(model))
        for constraint_type in get_duals(model)
            assign_dual_variable!(container, constraint_type, sys, T)
        end
    end
    return
end

function add_constraint_dual!(
    container::OptimizationContainer,
    sys::PSY.System,
    model::ServiceModel{T, D},
) where {T <: PSY.Service, D <: AbstractServiceFormulation}
    return
end

function assign_dual_variable!(
    container::OptimizationContainer,
    constraint_type::Type{<:ConstraintType},
    devices::U,
    formulation,
) where {U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}} where {D <: PSY.Component}
    @assert !isempty(devices)
    time_steps = get_time_steps(container)
    add_dual_container!(
        container,
        constraint_type,
        D,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )
    return
end

function assign_dual_variable!(
    container::OptimizationContainer,
    constraint_type::Type{<:ConstraintType},
    sys::U,
    formulation,
) where {U <: PSY.System}
    time_steps = get_time_steps(container)
    add_dual_container!(container, constraint_type, U, time_steps)
    return
end


text:
"""
    GaborKernel(; ell::Real=1.0, p::Real=1.0)

Gabor kernel with lengthscale `ell` and period `p`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the Gabor kernel with lengthscale ``l_i > 0``
and period ``p_i > 0`` is defined as
```math
k(x, x'; l, p) = \\exp\\bigg(- \\cos\\bigg(\\pi\\sum_{i=1}^d \\frac{x_i - x'_i}{p_i}\\bigg)
                             \\sum_{i=1}^d \\frac{(x_i - x'_i)^2}{l_i^2}\\bigg).
```
"""
struct GaborKernel{K<:Kernel} <: Kernel
    kernel::K
    function GaborKernel(; ell=nothing, p=nothing)
        k = _gabor(; ell=ell, p=p)
        return new{typeof(k)}(k)
    end
end

@functor GaborKernel

(κ::GaborKernel)(x, y) = κ.kernel(x, y)

function _gabor(; ell=nothing, p=nothing)
    if ell === nothing
        if p === nothing
            return SqExponentialKernel() * CosineKernel()
        else
            return SqExponentialKernel() * transform(CosineKernel(), 1 ./ p)
        end
    elseif p === nothing
        return transform(SqExponentialKernel(), 1 ./ ell) * CosineKernel()
    else
        return transform(SqExponentialKernel(), 1 ./ ell) *
               transform(CosineKernel(), 1 ./ p)
    end
end

function Base.getproperty(k::GaborKernel, v::Symbol)
    if v == :kernel
        return getfield(k, v)
    elseif v == :ell
        kernel1 = k.kernel.kernels[1]
        if kernel1 isa TransformedKernel
            return 1 ./ kernel1.transform.s[1]
        else
            return 1.0
        end
    elseif v == :p
        kernel2 = k.kernel.kernels[2]
        if kernel2 isa TransformedKernel
            return 1 ./ kernel2.transform.s[1]
        else
            return 1.0
        end
    else
        error("Invalid Property")
    end
end

function Base.show(io::IO, κ::GaborKernel)
    return print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")
end

kernelmatrix(κ::GaborKernel, x::AbstractVector) = kernelmatrix(κ.kernel, x)

function kernelmatrix(κ::GaborKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix(κ.kernel, x, y)
end

kernelmatrix_diag(κ::GaborKernel, x::AbstractVector) = kernelmatrix_diag(κ.kernel, x)


text:
function inf_geometric_sum(r::Real, θ::Real)
    x = 1.0 - r * cos(θ)
    y = -r * sin(θ)
    return inv(apy2(x, y)), -atan(y, x)
end

function sqr2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return abs2(w)
    return abs2(w) * (1.0 + abs2(z / w))
end

function apy2(x::Real, y::Real)
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return w
    return w * sqrt(1.0 + abs2(z / w))
end


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

module TypesTest

using Test
using Dates

# Date internal algorithms
@testset "totaldays" begin
    @test Dates.totaldays(0, 2, 28) == -307
    @test Dates.totaldays(0, 2, 29) == -306
    @test Dates.totaldays(0, 3, 1) == -305
    @test Dates.totaldays(0, 12, 31) == 0
    # Rata Die Days # start from 0001-01-01
    @test Dates.totaldays(1, 1, 1) == 1
    @test Dates.totaldays(1, 1, 2) == 2
    @test Dates.totaldays(2013, 1, 1) == 734869
end
@testset "daysinmonth" begin
    @test Dates.daysinmonth(2000, 1) == 31
    @test Dates.daysinmonth(2000, 2) == 29
    @test Dates.daysinmonth(2000, 3) == 31
    @test Dates.daysinmonth(2000, 4) == 30
    @test Dates.daysinmonth(2000, 5) == 31
    @test Dates.daysinmonth(2000, 6) == 30
    @test Dates.daysinmonth(2000, 7) == 31
    @test Dates.daysinmonth(2000, 8) == 31
    @test Dates.daysinmonth(2000, 9) == 30
    @test Dates.daysinmonth(2000, 10) == 31
    @test Dates.daysinmonth(2000, 11) == 30
    @test Dates.daysinmonth(2000, 12) == 31
    @test Dates.daysinmonth(2001, 2) == 28
end
@testset "isleapyear" begin
    @test Dates.isleapyear(1900) == false
    @test Dates.isleapyear(2000) == true
    @test Dates.isleapyear(2004) == true
    @test Dates.isleapyear(2008) == true
    @test Dates.isleapyear(0) == true
    @test Dates.isleapyear(1) == false
    @test Dates.isleapyear(-1) == false
    @test Dates.isleapyear(4) == true
    @test Dates.isleapyear(-4) == true
end
# Create "test" check manually
y = Dates.Year(1)
m = Dates.Month(1)
w = Dates.Week(1)
d = Dates.Day(1)
h = Dates.Hour(1)
mi = Dates.Minute(1)
s = Dates.Second(1)
ms = Dates.Millisecond(1)
@testset "DateTime construction by parts" begin
    test = Dates.DateTime(Dates.UTM(63492681600000))
    @test Dates.DateTime(2013) == test
    @test Dates.DateTime(2013, 1) == test
    @test Dates.DateTime(2013, 1, 1) == test
    @test Dates.DateTime(2013, 1, 1, 0) == test
    @test Dates.DateTime(2013, 1, 1, 0, 0) == test
    @test Dates.DateTime(2013, 1, 1, 0, 0, 0) == test
    @test Dates.DateTime(2013, 1, 1, 0, 0, 0, 0) == test

    @test Dates.DateTime(y) == Dates.DateTime(1)
    @test Dates.DateTime(y, m) == Dates.DateTime(1, 1)
    @test Dates.DateTime(y, m, d) == Dates.DateTime(1, 1, 1)
    @test Dates.DateTime(y, m, d, h) == Dates.DateTime(1, 1, 1, 1)
    @test Dates.DateTime(y, m, d, h, mi) == Dates.DateTime(1, 1, 1, 1, 1)
    @test Dates.DateTime(y, m, d, h, mi, s) == Dates.DateTime(1, 1, 1, 1, 1, 1)
    @test Dates.DateTime(y, m, d, h, mi, s, ms) == Dates.DateTime(1, 1, 1, 1, 1, 1, 1)
    @test Dates.DateTime(Dates.Day(10), Dates.Month(2), y) == Dates.DateTime(1, 2, 10)
    @test Dates.DateTime(Dates.Second(10), Dates.Month(2), y, Dates.Hour(4)) == Dates.DateTime(1, 2, 1, 4, 0, 10)
    @test Dates.DateTime(Dates.Year(1), Dates.Month(2), Dates.Day(1),
                         Dates.Hour(4), Dates.Second(10)) == Dates.DateTime(1, 2, 1, 4, 0, 10)
end

@testset "Date construction by parts" begin
    test = Dates.Date(Dates.UTD(734869))
    @test Dates.Date(2013) == test
    @test Dates.Date(2013, 1) == test
    @test Dates.Date(2013, 1, 1) == test
    @test Dates.Date(y) == Dates.Date(1)
    @test Dates.Date(y, m) == Dates.Date(1, 1)
    @test Dates.Date(y, m, d) == Dates.Date(1, 1, 1)
    @test Dates.Date(m) == Dates.Date(1, 1, 1)
    @test Dates.Date(d, y) == Dates.Date(1, 1, 1)
    @test Dates.Date(d, m) == Dates.Date(1, 1, 1)
    @test Dates.Date(m, y) == Dates.Date(1, 1, 1)
    @test Dates.Date(Dates.Day(10), Dates.Month(2), y) == Dates.Date(1, 2, 10)
end

@testset "Time construction by parts" begin
    t = Dates.Time(Dates.Nanosecond(82800000000000))
    @test Dates.Time(23) == t
    @test Dates.Time(23, 0) == t
    @test Dates.Time(23, 0, 0) == t
    @test Dates.Time(23, 0, 0, 0) == t
    @test Dates.Time(23, 0, 0, 0, 0) == t
    @test Dates.Time(23, 0, 0, 0, 0, 0) == t
    us = Dates.Microsecond(1)
    ns = Dates.Nanosecond(1)
    @test Dates.Time(h) == Dates.Time(1)
    @test Dates.Time(h, mi) == Dates.Time(1, 1)
    @test Dates.Time(h, mi, s) == Dates.Time(1, 1, 1)
    @test Dates.Time(h, mi, s, ms) == Dates.Time(1, 1, 1, 1)
    @test Dates.Time(h, mi, s, ms, us) == Dates.Time(1, 1, 1, 1, 1)
    @test Dates.Time(h, mi, s, ms, us, ns) == Dates.Time(1, 1, 1, 1, 1, 1)
    @test Dates.Time(us, h, s, ns, mi, ms) == Dates.Time(1, 1, 1, 1, 1, 1)
    @test Dates.Time(Dates.Second(10), Dates.Minute(2), us, Dates.Hour(4)) == Dates.Time(4, 2, 10, 0, 1)
    @test Dates.Time(Dates.Hour(4), Dates.Second(10), Dates.Millisecond(15),
                     Dates.Microsecond(20), Dates.Nanosecond(25)) == Dates.Time(4, 0, 10, 15, 20, 25)
end

@testset "various input types for Date/DateTime" begin
    test = Dates.Date(1, 1, 1)
    @test Dates.Date(Int8(1), Int8(1), Int8(1)) == test
    @test Dates.Date(UInt8(1), UInt8(1), UInt8(1)) == test
    @test Dates.Date(Int16(1), Int16(1), Int16(1)) == test
    @test Dates.Date(UInt8(1), UInt8(1), UInt8(1)) == test
    @test Dates.Date(Int32(1), Int32(1), Int32(1)) == test
    @test Dates.Date(UInt32(1), UInt32(1), UInt32(1)) == test
    @test Dates.Date(Int64(1), Int64(1), Int64(1)) == test
    @test Dates.Date('\x01', '\x01', '\x01') == test
    @test Dates.Date(true, true, true) == test
    @test_throws ArgumentError Dates.Date(false, true, false)
    @test Dates.Date(false, true, true) == test - Dates.Year(1)
    @test_throws ArgumentError Dates.Date(true, true, false)
    @test Dates.Date(UInt64(1), UInt64(1), UInt64(1)) == test
    @test Dates.Date(-1, UInt64(1), UInt64(1)) == test - Dates.Year(2)
    @test Dates.Date(Int128(1), Int128(1), Int128(1)) == test
    @test_throws InexactError Dates.Date(170141183460469231731687303715884105727, Int128(1), Int128(1))
    @test Dates.Date(UInt128(1), UInt128(1), UInt128(1)) == test
    @test Dates.Date(big(1), big(1), big(1)) == test
    @test Dates.Date(big(1), big(1), big(1)) == test
    # Potentially won't work if can't losslessly convert to Int64
    @test Dates.Date(BigFloat(1), BigFloat(1), BigFloat(1)) == test
    @test Dates.Date(complex(1), complex(1), complex(1)) == test
    @test Dates.Date(Float64(1), Float64(1), Float64(1)) == test
    @test Dates.Date(Float32(1), Float32(1), Float32(1)) == test
    @test Dates.Date(Float16(1), Float16(1), Float16(1)) == test
    @test Dates.Date(Rational(1), Rational(1), Rational(1)) == test
    @test_throws InexactError Dates.Date(BigFloat(1.2), BigFloat(1), BigFloat(1))
    @test_throws InexactError Dates.Date(1 + im, complex(1), complex(1))
    @test_throws InexactError Dates.Date(1.2, 1.0, 1.0)
    @test_throws InexactError Dates.Date(1.2f0, 1.f0, 1.f0)
    @test_throws InexactError Dates.Date(3//4, Rational(1), Rational(1)) == test

    # Months, days, hours, minutes, seconds, and milliseconds must be in range
    @test_throws ArgumentError Dates.Date(2013, 0, 1)
    @test_throws ArgumentError Dates.Date(2013, 13, 1)
    @test_throws ArgumentError Dates.Date(2013, 1, 0)
    @test_throws ArgumentError Dates.Date(2013, 1, 32)
    @test_throws ArgumentError Dates.DateTime(2013, 0, 1)
    @test_throws ArgumentError Dates.DateTime(2013, 13, 1)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 0)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 32)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 25)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, -1)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, -1)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, 60)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, 0, -1)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, 0, 60)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, 0, 0, -1)
    @test_throws ArgumentError Dates.DateTime(2013, 1, 1, 0, 0, 0, 1000)
    @test_throws ArgumentError Dates.Time(24)
    @test_throws ArgumentError Dates.Time(-1)
    @test_throws ArgumentError Dates.Time(0, -1)
    @test_throws ArgumentError Dates.Time(0, 60)
    @test_throws ArgumentError Dates.Time(0, 0, -1)
    @test_throws ArgumentError Dates.Time(0, 0, 60)
    @test_throws ArgumentError Dates.Time(0, 0, 0, -1)
    @test_throws ArgumentError Dates.Time(0, 0, 0, 1000)
    @test_throws ArgumentError Dates.Time(0, 0, 0, 0, -1)
    @test_throws ArgumentError Dates.Time(0, 0, 0, 0, 1000)
    @test_throws ArgumentError Dates.Time(0, 0, 0, 0, 0, -1)
    @test_throws ArgumentError Dates.Time(0, 0, 0, 0, 0, 1000)
end
a = Dates.DateTime(2000)
b = Dates.Date(2000)
c = Dates.Time(0)
@testset "DateTime traits" begin
    @test Dates.calendar(a) == Dates.ISOCalendar
    @test Dates.calendar(b) == Dates.ISOCalendar
    @test eps(a) == Dates.Millisecond(1)
    @test eps(b) == Dates.Day(1)
    @test eps(c) == Dates.Nanosecond(1)
    @test string(typemax(Dates.DateTime)) == "146138512-12-31T23:59:59"
    @test string(typemin(Dates.DateTime)) == "-146138511-01-01T00:00:00"
    @test typemax(Dates.DateTime) - typemin(Dates.DateTime) == Dates.Millisecond(9223372017043199000)
    @test string(typemax(Dates.Date)) == "252522163911149-12-31"
    @test string(typemin(Dates.Date)) == "-252522163911150-01-01"
    @test string(typemax(Dates.Time)) == "23:59:59.999999999"
    @test string(typemin(Dates.Time)) == "00:00:00"
    @test isfinite(Dates.Date)
    @test isfinite(Dates.DateTime)
    @test isfinite(Dates.Time)
    @test c == c
    @test c == (c + Dates.Hour(24))
    @test hash(c) == hash(c + Dates.Hour(24))
    @test hash(c + Dates.Nanosecond(10)) == hash(c + Dates.Hour(24) + Dates.Nanosecond(10))
end
@testset "Date-DateTime conversion/promotion" begin
    global a, b, c, d
    @test Dates.DateTime(a) == a
    @test Dates.Date(a) == b
    @test Dates.DateTime(b) == a
    @test Dates.Date(b) == b
    @test a == b
    @test a == a
    @test b == a
    @test b == b
    @test !(a < b)
    @test !(b < a)
    c = Dates.DateTime(2000)
    d = Dates.Date(2000)
    @test ==(a, c)
    @test ==(c, a)
    @test ==(d, b)
    @test ==(b, d)
    @test ==(a, d)
    @test ==(d, a)
    @test ==(b, c)
    @test ==(c, b)
    b = Dates.Date(2001)
    @test b > a
    @test a < b
    @test a != b
    @test Dates.Date(Dates.DateTime(Dates.Date(2012, 7, 1))) == Dates.Date(2012, 7, 1)
end

@testset "min and max" begin
    for (a, b) in [(Dates.Date(2000), Dates.Date(2001)),
                    (Dates.Time(10), Dates.Time(11)),
                    (Dates.DateTime(3000), Dates.DateTime(3001)),
                    (Dates.Week(42), Dates.Week(1972))]
        @test min(a, b) == a
        @test min(b, a) == a
        @test min(a) == a
        @test max(a, b) == b
        @test max(b, a) == b
        @test max(b) == b
        @test minmax(a, b) == (a, b)
        @test minmax(b, a) == (a, b)
        @test minmax(a) == (a, a)
    end
end

@testset "issue #31524" begin
    dt1 = Libc.strptime("%Y-%M-%dT%H:%M:%SZ", "2018-11-16T10:26:14Z")
    dt2 = Base.Libc.TmStruct(14, 30, 5, 10, 1, 99, 3, 40, 0)

    time = Time(dt1)
    @test typeof(time) == Time
    @test time == Dates.Time(10, 26, 14)

    date = Date(dt2)
    @test typeof(date) == Date
    @test date == Dates.Date(1999, 2, 10)

    datetime = DateTime(dt2)
    @test typeof(datetime) == DateTime
    @test datetime == Dates.DateTime(1999, 2, 10, 5, 30, 14)

end

end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
# # Reading TIFFs

# Loading most TIFFs should just work, see [Writing TIFFs](@ref) for more
# advanced manipulation of TIFF objects. But we'll quickly run through a common
# use cases.

#md # ```@contents
#md # Pages = ["reading.md"]
#md # Depth = 5
#md # ```

get_example(name) = download("https://github.com/tlnagy/exampletiffs/blob/master/$name?raw=true") #hide
filepath = get_example("spring.tif")                                                              #hide

# ### Basic loading

# At it's most basic, we can just point `TiffImages.jl` to the filepath of an
# image and it will attempt to load it. Here, we're loading `spring.tif` from
# the [`tlnagy/exampletiffs`](https://github.com/tlnagy/exampletiffs) repo

using TiffImages
img = TiffImages.load(filepath)

# If you're a graphical environment, you can load the `Images.jl` repo to get
# a nice graphical representation of your image. If you're in the REPL, I highly
# recommend the `ImageInTerminal.jl` package for some visual feedback.

# Continuing on, `img` here behaves exactly as you would expect a Julian array
# to despite the funky type signature

typeof(img)

# Everything should behave as expected

eltype(img)
#---------

# Accessing and setting data should work as expected
img[160:180, 50]
#---------
img[160:180, 50] .= 1.0
img

# ### Memory-mapped files

# `TiffImages.jl` also supports memory-mapping so that you can load images that
# are larger than the available memory.

img = TiffImages.load(filepath, mmap=true)
img[:, :, 1]

# Naturally, this only really makes sense for large images, not the single pane
# image here. Currently, memory-mapped files are readonly so trying to set a
# value will throw an error

# ```julia
# img[1, 1, 1] = 0 # this won't work
# ```


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Edward J. Xu
# 2019.5.26
# include("$(homedir())/Documents/Github/DantzigWolfeDecomposition/Test/test.jl")
# push!(LOAD_PATH, "$(homedir())/Documents/Github/DantzigWolfeDecomposition")
# cd("$(homedir())/Documents/Github/DantzigWolfeDecomposition")

include("FuncData.jl")
using DantzigWolfeDecomposition

function main(strFileName::String, series::Int64)
    println("#### 1/2,  Prepare Data ########################################################") ########################
    (mat_a, vec_b, vec_c, vecSenseAll, indexMas, blocks, indexSub, numXPerSub) = setGeneralAssignProbDate(
        strFileName, series)
    println("#### 2/2,  Begin Optim #########################################################") ########################
    timeStart = time()
    ## Set parameter and start DW-Decomposition
    dualPen = 10
    dualPenMult = 0.1
    dualPenThreshold = 0.01 - 1e-5
    epsilon = 0.00001                # [maximum difference between two bounds]
    whePrint = false                 # [whether to print detailed info]
    whiSolver = 1                    # 1: Gurobi, 2: CPLEX, 3: GLPK
    doDWDecomp(
        mat_a, vec_b, vec_c,                                  # Data in LP Problem
        vecSenseAll, indexMas, blocks, indexSub, numXPerSub,  # Data for DW-Decomp
        dualPen, dualPenMult, dualPenThreshold,               # Para for Stable
        epsilon, whePrint, whiSolver                          # Control Para
        )
    println("Elapsed time is $(time() - timeStart) seconds.")
    println("#### End #######################################################################") ########################
end

main("Data/gapa.txt", 1)


text:
# we could just store intervals as a rational number
# but there's more than one way to represent a given rational
# so we don't want to lose user information
mutable struct Interval
    numerator::Int
    denominator::Int
    octave::Int
end

function Interval(; numerator = 1, denominator = 1, octave = 0)
    Interval(numerator, denominator, octave)
end

function Rational(interval::Interval)
    interval.numerator // interval.denominator * (2 // 1)^interval.octave
end

const INTERVAL_REGEX = r"(?<numerator>[^/o]+)(?:/(?<denominator>[^o]+))?(?:o(?<octave>.+))?"

# TODO: reuse from AudioSchedules
function parse(::Type{Interval}, text::AbstractString; line_number)
    a_match = match(INTERVAL_REGEX, text)
    if a_match === nothing
        throw_parse_error(text, "interval", line_number)
    else
        numerator_string = a_match["numerator"]
        numerator = tryparse(Int, numerator_string)
        denominator_string = a_match["denominator"]
        octave_string = a_match["octave"]
        Interval(
            if numerator === nothing
                throw_parse_error(numerator_string, "numerator", line_number)
            else
                numerator
            end,
            if denominator_string === nothing
                1
            else
                denominator = tryparse(Int, denominator_string)
                if denominator === nothing
                    throw_parse_error(denominator_string, "denominator", line_number)
                else
                    denominator
                end
            end,
            if octave_string === nothing
                0
            else
                octave = tryparse(Int, octave_string)
                if octave === nothing
                    throw_parse_error(octave_string, "octave", line_number)
                else
                    octave
                end
            end,
        )
    end
end

function print(io::IO, interval::Interval)
    denominator = interval.denominator
    octave = interval.octave
    print(io, interval.numerator)
    if denominator != 1
        print(io, '/')
        print(io, denominator)
    end
    if octave != 0
        print(io, 'o')
        print(io, octave)
    end
end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
using PyCall
using JLD

# @pyimport sympy.physics.quantum.cg as cg
@pyimport sympy.physics.wigner as wigner

const AngMom = Union{Int, Float64}
const AngMomVector = Union{Vector{Int}, Vector{Float64}} # unfortunately this seems to be necessary
isangmom(q::AngMom) = isinteger(2q)
isangmom(qs::AngMomVector) = all(isinteger, 2qs)

m(alpha::Int, ell::AngMom)::Float64 = (alpha - 1) - ell

Aplus(ell::AngMom, m::AngMom)::Float64 = sqrt( ell * (ell + 1) - m * (m + 1) )
Aminus(ell::AngMom, m::AngMom)::Float64 = sqrt( ell * (ell + 1) - m * (m - 1) )

# CG(j1, m1, j2, m2, j3, m3) = (cg.CG(j1, m1, j2, m2, j3, m3)[:doit]() |> float)
CG(j1, m1, j2, m2, j3, m3) = (wigner.clebsch_gordan(j1, j2, j3, m1, m2, m3)[:doit]() |> float)
Wigner3j(j1, j2, j3, m1, m2, m3) = (wigner.wigner_3j(j1, j2, j3, m1, m2, m3)[:doit]() |> float)
Wigner6j(j1, j2, j3, j4, j5, j6) = (wigner.wigner_6j(j1, j2, j3, j4, j5, j6)[:doit]() |> float)

Lqn(v) = 0.5 * (-1 + sqrt.(1 + 4v)) # v = L * (L + 1)

function build_CG_table(ell1::AngMom, ell2::AngMom)
    @assert isangmom(ell1)
    @assert isangmom(ell2)

    dim1 = Int(2ell1) + 1
    dim2 = Int(2ell2) + 1

    Ls = abs(ell1-ell2):ell1+ell2
    dim3 = length(Ls)

    println("Building CG table for ell1 = $(ell1), ell2 = $(ell2):")
    tic()
    CG_table = Array{Vector{Float64}}(dim1, dim2, dim3)
    for alpha1 in 1:dim1
        println("Working on alpha1 = $(alpha1): $(alpha1)/$(dim1)")
        for alpha2 in 1:dim2
            # println("Working on alpha2 = $(alpha2): $(alpha2)/$(dim2)")
            for alpha3 in 1:dim3
                tmp = Float64[]
                for M in -Ls[alpha3]:Ls[alpha3]
                    push!(tmp,  CG(float(ell1), m(alpha1, ell1), float(ell2), m(alpha2, ell2), Ls[alpha3], M))
                end
                CG_table[alpha1, alpha2, alpha3] = tmp
            end
        end
    end
    toc()
    println()

    return CG_table
end

build_CG_table(ell::AngMom) = build_CG_table(ell, ell)

function write_CG_table(ell1::AngMom, ell2::AngMom, dir::String)
    @assert isangmom(ell1)
    @assert isangmom(ell2)

    filename = "CG_$(float(ell1))_$(float(ell2)).jld"
    fullpath = joinpath(dir, filename)

    println("Writing CG table for ell1 = $(ell1), ell2 = $(ell2) to $(fullpath) ...")
    table = build_CG_table(ell1, ell2)
    save(fullpath, "table", table)
    println("Write successful")
end
write_CG_table(ell::AngMom, dir::String) = write_CG_table(ell, ell, dir)

function load_CG_table(ell1::AngMom, ell2::AngMom, dir::String)
    @assert isangmom(ell1)
    @assert isangmom(ell2)

    filename = "CG_$(float(ell1))_$(float(ell2)).jld"
    fullpath = joinpath(dir, filename)

    println("Loading CG table for ell1 = $(ell1), ell2 = $(ell2) from $(fullpath) ...")
    table = load(fullpath, "table")
    println("Load successful")
    return table
end
load_CG_table(ell::AngMom, dir::String) = load_CG_table(ell, ell, dir)

function compile_CG_tables(ell1vec::AngMomVector, ell2vec::AngMomVector, dir::String)
    for ell1 in ell1vec
        for ell2 in ell2vec
            write_CG_table(ell1, ell2, dir)
        end
    end
end

function compile_CG_tables(ellvec::AngMomVector, dir::String)
    for ell in ellvec
        write_CG_table(ell, dir)
    end
end

u(θ::Float64, ϕ::Float64)::Complex128 = cos(θ/2) * exp(im * ϕ/2)
v(θ::Float64, ϕ::Float64)::Complex128 = sin(θ/2) * exp(-im * ϕ/2)

function monopole_harmonic(Q::AngMom, ell::AngMom, m::AngMom, θ::Float64, ϕ::Float64)::Complex128
    @assert isangmom(Q)
    @assert isangmom(ell)
    @assert isangmom(m)

    @assert Q >= 0 # for simplicity
    @assert ell >= Q
    @assert abs(m) <= ell
    # could do more thorough consistency checks regarding integers versus half integers etc., but these will get caught below

    if m < Q
        N = prod([(ell - m - i)/(ell + Q - i) for i in 0:Q-m-1])
    elseif m > Q
        N = prod([(ell + m - i)/(ell - Q - i) for i in 0:m-Q-1])
    else
        N = 1.
    end
    N = sqrt((2ell + 1)/(4pi) * N)

    Y = sum([(-1.)^s * binomial(Int(ell-Q), s) * binomial(Int(ell+Q), Int(ell-m-s)) * sin(θ/2)^(2ell-2Q-2s) * cos(θ/2)^(2s) for s in 0:Int(ell-Q)])
    Y *= N * (-1.)^(ell - m) * v(θ, ϕ)^(Q - m) * u(θ, ϕ)^(Q + m)

    return Y
end

max_stars_repo_name: mishmash/QHED

max_issues_repo_name: mishmash/QHED

max_forks_repo_name: mishmash/QHED


text:
module Interpolations

using Base.Cartesian
using WoodburyMatrices

import Base:
    eltype,
    gradient,
    getindex,
    ndims,
    size

export
    Interpolation,
    Constant,
    Linear,
    Quadratic,
    ExtrapError,
    ExtrapNaN,
    ExtrapConstant,
    ExtrapLinear,
    ExtrapReflect,
    ExtrapPeriodic,
    OnCell,
    OnGrid,
    Flat,
    Line,
    Free,
    Periodic,
    Reflect,

    degree,
    boundarycondition,
    gridrepresentation,

    gradient,
    gradient!

abstract Degree{N}

abstract GridRepresentation
immutable OnGrid <: GridRepresentation end
immutable OnCell <: GridRepresentation end

abstract BoundaryCondition
immutable None <: BoundaryCondition end
immutable Flat <: BoundaryCondition end
immutable Line <: BoundaryCondition end
typealias Natural Line
immutable Free <: BoundaryCondition end
immutable Periodic <: BoundaryCondition end
immutable Reflect <: BoundaryCondition end

abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}

include("extrapolation.jl")

abstract AbstractInterpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractArray{T,N}
type Interpolation{T,N,TCoefs,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,IT,EB}
    coefs::Array{TCoefs,N}
end
function Interpolation{N,TCoefs,TWeights<:Real,IT<:InterpolationType,EB<:ExtrapolationBehavior}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, it::IT, ::EB)
    isleaftype(IT) || error("The interpolation type must be a leaf type (was $IT)")

    isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type T (eltype(A) == $(eltype(A)))")

    c = one(TWeights)
    for _ in 2:N
        c *= c
    end
    T = typeof(c * one(TCoefs))

    Interpolation{T,N,TCoefs,IT,EB}(prefilter(TWeights,A,it))
end
Interpolation(A::AbstractArray, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float64, A, it, eb)
Interpolation(A::AbstractArray{Float32}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float32, A, it, eb)
Interpolation(A::AbstractArray{Rational{Int}}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Rational{Int}, A, it, eb)

# Unless otherwise specified, use coefficients as they are, i.e. without prefiltering
# However, all prefilters copy the array, so do that here as well
prefilter{TWeights,T,N,IT<:InterpolationType}(::Type{TWeights}, A::AbstractArray{T,N}, ::IT) = copy(A)

size{T,N}(itp::Interpolation{T,N}, d::Integer) = d > N ? 1 : size(itp.coefs, d) - 2*padding(interpolationtype(itp))
size(itp::AbstractInterpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
ndims(itp::Interpolation) = ndims(itp.coefs)

offsetsym(off, d) = off == -1 ? symbol("ixm_", d) :
                    off ==  0 ? symbol("ix_", d) :
                    off ==  1 ? symbol("ixp_", d) :
                    off ==  2 ? symbol("ixpp_", d) : error("offset $off not recognized")

interpolationtype{T,N,TCoefs,IT<:InterpolationType}(itp::Interpolation{T,N,TCoefs,IT}) = IT()

boundarycondition{D,BC<:BoundaryCondition}(::InterpolationType{D,BC}) = BC()
gridrepresentation{D,BC,GR<:GridRepresentation}(::InterpolationType{D,BC,GR}) = GR()
degree{D<:Degree,BC,GR}(::InterpolationType{D,BC,GR}) = D()

# If nothing else is specified, don't pad at all
padding(::InterpolationType) = 0

for op in (:boundarycondition, :gridrepresentation, :degree, :padding)
    eval(:($(op)(itp::Interpolation) = $(op)(interpolationtype(itp))))
end

include("constant.jl")
include("linear.jl")
include("quadratic.jl")

function padded_index(sz::Tuple, pad)
    sz = Int[s+2pad for s in sz]
    N = length(sz)
    ind = cell(N)
    for i in 1:N
        ind[i] = 1+pad:sz[i]-pad
    end
    ind,sz
end
function copy_with_padding(A, it::InterpolationType)
    pad = padding(it)
    ind,sz = padded_index(size(A), pad)
    coefs = zeros(eltype(A), sz...)
    coefs[ind...] = A
    coefs, pad
end

supported(::Type{Constant{OnGrid}}) = true
supported(::Type{Constant{OnCell}}) = true
supported(::Type{Linear{OnGrid}}) = true
supported(::Type{Linear{OnCell}}) = true
supported(::Type{Quadratic{Flat,OnCell}}) = true
supported(::Type{Quadratic{Flat,OnGrid}}) = true
supported(::Type{Quadratic{Reflect,OnCell}}) = true
supported(::Type{Quadratic{Reflect,OnGrid}}) = true
supported(::Type{Quadratic{Line,OnGrid}}) = true
supported(::Type{Quadratic{Line,OnCell}}) = true
supported(::Type{Quadratic{Free,OnGrid}}) = true
supported(::Type{Quadratic{Free,OnCell}}) = true
supported(::Type{Quadratic{Periodic,OnGrid}}) = true
supported(::Type{Quadratic{Periodic,OnCell}}) = true
supported(::Any) = false

function getindex_impl(N, IT::Type, EB::Type)
    it = IT()
    eb = EB()
    gr = gridrepresentation(it)
    quote
        @nexprs $N d->(x_d = x[d])

        # Handle extrapolation, by either throwing an error,
        # returning a value, or shifting x to somewhere inside
        # the domain
        $(extrap_transform_x(gr,eb,N))

        # Calculate the indices of all coefficients that will be used
        # and define fx = x - xi in each dimension
        $(define_indices(it, N))

        # Calculate coefficient weights based on fx
        $(coefficients(it, N))

        # Generate the indexing expression
        @inbounds ret = $(index_gen(degree(it), N))
        ret
    end
end

# Resolve ambiguity with linear indexing
getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::Real) =
    error("Linear indexing is not supported for interpolation objects")
stagedfunction getindex{T,TCoefs,IT,EB}(itp::Interpolation{T,1,TCoefs,IT,EB}, x::Real)
    getindex_impl(1, IT, EB)
end

# Resolve ambiguity with real multilinear indexing
stagedfunction getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::Real...)
    n = length(x)
    n == N || return :(error("Must index ", $N, "-dimensional interpolation objects with ", $(nindexes(N))))
    getindex_impl(N, IT, EB)
end

# Resolve ambiguity with general array indexing,
#   getindex{T,N}(A::AbstractArray{T,N}, I::AbstractArray{T,N})
function getindex{T,N,TCoefs,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, I::AbstractArray{T})
    error("Array indexing is not defined for interpolation objects.")
end

# Resolve ambiguity with colon indexing for 1D interpolations
#   getindex{T}(A::AbstractArray{T,1}, C::Colon)
function getindex{T,TCoefs,IT,EB}(itp::Interpolation{T,1,TCoefs,IT,EB}, c::Colon)
    error("Colon indexing is not supported for interpolation objects")
end

# Allow multilinear indexing with any types
stagedfunction getindex{T,N,TCoefs,TIndex,IT,EB}(itp::Interpolation{T,N,TCoefs,IT,EB}, x::TIndex...)
    n = length(x)
    n == N || return :(error("Must index ", $N, "-dimensional interpolation objects with ", $(nindexes(N))))
    getindex_impl(N, IT, EB)
end

# Define in-place gradient calculation
#   gradient!(g::Vector, itp::Interpolation, NTuple{N}...)
stagedfunction gradient!{T,N,TCoefs,TOut,IT,EB}(g::Vector{TOut}, itp::Interpolation{T,N,TCoefs,IT,EB}, x...)
    n = length(x)
    n == N || return :(error("Must index ", $N, "-dimensional interpolation objects with ", $(nindexes(N))))
    it = IT()
    eb = EB()
    gr = gridrepresentation(it)
    quote
        length(g) == $N || error("g must be an array with exactly N elements (length(g) == "*string(length(g))*", N == "*string($N)*")")
        @nexprs $N d->(x_d = x[d])
        $(extrap_transform_x(gr,eb,N))
        $(define_indices(it,N))
        @nexprs $N dim->begin
            @nexprs $N d->begin
                (d==dim ? $(gradient_coefficients(it,N,:d))
                        : $(coefficients(it,N,:d)))
            end

            @inbounds g[dim] = $(index_gen(degree(it),N))
        end
        g
    end
end

gradient{T,N}(itp::Interpolation{T,N}, x...) = gradient!(Array(T,N), itp, x...)

# This creates prefilter specializations for all interpolation types that need them
stagedfunction prefilter{TWeights,TCoefs,N,IT<:Quadratic}(::Type{TWeights}, A::Array{TCoefs,N}, it::IT)
    quote
        ret, pad = copy_with_padding(A, it)

        szs = collect(size(ret))
        strds = collect(strides(ret))

        for dim in 1:$N
            n = szs[dim]
            szs[dim] = 1

            M, b = prefiltering_system(TWeights, TCoefs, n, it)

            @nloops $N i d->1:szs[d] begin
                cc = @ntuple $N i

                strt_diff = sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
                strt = 1 + strt_diff
                rng = range(strt, strds[dim], n)

                bdiff = ret[rng]
                b += bdiff
                ret[rng] = M \ b
                b -= bdiff
            end
            szs[dim] = n
        end
        ret
    end
end

nindexes(N::Int) = N == 1 ? "1 index" : "$N indexes"

end # module


text:
dv = dataones(3)
@assert isequal(dv, padNA(dv, 0, 0))
@assert length(padNA(dv, 2, 0)) == length(dv) + 2
@assert length(padNA(dv, 0, 2)) == length(dv) + 2
@assert length(padNA(dv, 2, 2)) == length(dv) + 4


text:
"""
    wshenk(object::Union{Pca, Plsr}, X; nlv = nothing)
Compute the Shenk et al. (1997) PLSR weights
* `object` : The fitted model.
* `X` : X-data on which the weights are computed.
* `nlv` : Nb. latent variables (LVs) to consider. If nothing,
    it is the maximum nb. of components.

For each observation (row) of `X`, the weights are returned
for the models with 1, ..., nlv LVs.

## References

Shenk, J., Westerhaus, M., Berzaghi, P., 1997. Investigation of a LOCAL calibration
procedure for near infrared instruments.
Journal of Near Infrared Spectroscopy 5, 223. https://doi.org/10.1255/jnirs.115

"""
function wshenk(object::Union{Plsr}, X; nlv = nothing, norm = false)
    X = ensure_mat(X)
    m = size(X, 1)
    p = size(X, 2)
    q = size(object.C, 1)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    E = similar(X)
    rms_r = similar(X, m, a)
    zB = similar(X, p + 1, q)
    rms_b = similar(X, a)
    W = similar(X, m, a)
    for j = 1:a
        E .= xresid(object, X, nlv = j)
        rms_r[:, j] .= vec(sqrt.(sum(E.^2, dims = 2)))
        z = coef(object; nlv = j)
        zB[1, :] .= vec(z.int)
        zB[2:end, :] .= z.B
        rms_b[j] = sqrt(mean(mean(zB.^2, dims = 2)))
        W[:, j] .= 1 ./ (vcol(rms_r, j) * rms_b[j])
    end
    for i = 1:m
        W[i, :] ./= sum(vrow(W, i))
    end
    (W = W, rms_r = rms_r, rms_b = rms_b)
end


text:
# compare-hyperparameters-multivariatetpseudohuber
#
# This script compares the quality of samples drawn
# from a Bayesian regession for the AIS dataset.

include("experiment_utils.jl")
include("ais/ais_utils.jl")

# random number generator seed
# this controls the seed used for the subsampling
@parseintcli seed "s" "seed" 7
# set the sampler
@parsestringcli sampler "q" "sampler" "runsgrld"
# set the delta param for the Huber prior
@parsestringcli clinu "n" "nu" "10.0"
nu = float(clinu)
# set the nu param degrees of freedom for the d-distribution
@parsestringcli clidelta "d" "delta" "0.1"
delta = float(clidelta)

using DataFrames

using SteinDiscrepancy: gsd
using SteinDistributions:
    SteinMultivariateStudentTRegressionPseudoHuberPrior,
    gradlogdensity
using SteinSamplers: runsgrld

# Select a solver for Stein discrepancy optimization problem
#solver = "clp"
solver = "gurobi"
# distname is multivariatetpseudohuber
distname = "multivariatetpseudohuber"
###            Pull in data here             ###
### ais dataset is available in sn R package ###
data_set = "ais"
yfeature = "LBM"
dat = readtable("src/experiments/ais/ais.tsv")
X = convert(Array{Float64,2}, dat[[:RCC, :WCC, :Fe]])
# standardize columns of X
X = X .- mean(X,1)
X = [ones(size(X,1)) X]
(N, d) = size(X)
column_norms = zeros(1, d)
for ii in 1:d
    column_norms[1,ii] = norm(X[:,ii])
end
X = X ./ column_norms
y = convert(Array{Float64,1}, dat[:LBM])
y = (y - mean(y)) / std(y)
# create posterior distribution
sigma = Diagonal(ones(size(X,1)))
target = SteinMultivariateStudentTRegressionPseudoHuberPrior(X, y, nu, sigma, zeros(d), delta)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# starting beta0
beta0 = getapproxmode(target)
# batchsize for SGRLD algorithm
batchsize = 30
# the SGRLD epsilons to run the experiemnt
sgrld_epsilons = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
# define the covariance terms for the Riemannian-Langevin diffusion
function volatility_covariance(x::Array{Float64,1})
    p = length(x)
    sqrt(1 + sum(x.^2) / delta^2) * eye(p)
end
function grad_volatility_covariance(x::Array{Float64,1})
    delta^(-2) * x / sqrt(1 + sum(x.^2) / delta^2)
end

# Define sampler
n = 2000
thinby = 100
numsweeps = round(Int, ceil(n * thinby / floor(N/batchsize)))
epsilons = sgrld_epsilons
runsampler(x0,epsilon) = runsgrld(target,
                                  x0,
                                  volatility_covariance,
                                  grad_volatility_covariance;
                                  epsilonfunc=t -> epsilon,
                                  numsweeps=numsweeps,
                                  batchsize=batchsize)

for epsilon in epsilons
    println("Generating samples for sampler=$(sampler),epsilon=$(epsilon)")

    betas = []; numeval = NaN
    @trycatchcontinue(
        begin
            (betas, numeval) = @setseed seed runsampler(beta0,epsilon);
            if any(isnan(betas))
                error("[$(distname):seed=$(seed),epsilon=$(epsilon)] NaNs found. Skipping.");
            end;
            betas = betas[thinby:thinby:(thinby*n),:]
        end,
        println("[$(distname):seed=$(seed),epsilon=$(epsilon)]:")
    )

    # Find stein discrepancy
    @printf("[Computing the stein discrepancy]\n")
    result = @setseed seed gsd(points=betas,
                               solver=solver,
                               gradlogdensity=gradlogp,
                               operator="riemannian-langevin",
                               volatility_covariance=volatility_covariance,
                               grad_volatility_covariance=grad_volatility_covariance)

    instance_data = Dict{Any, Any}(
        "distname" => distname,
        "n" => size(betas,1),
        "d" => d,
        "seed" => seed,
        "thinby" => thinby,
        "epsilon" => epsilon,
        "betas" => betas,
        "objectivevalue" => vec(result.objectivevalue),
        "numeval" => numeval,
        "sampler" => sampler,
        "nu" => nu,
        "delta" => delta,
        "batchsize" => batchsize,
        "nprocs" => nprocs(),
        "edgetime" => result.edgetime,
        "solvertime" => result.solvetime
    )

    save_json(
        instance_data;
        dir="compare-hyperparameters-multivariatetpseudohuber",
        distname=distname,
        n=size(betas,1),
        thinby=thinby,
        sampler=sampler,
        epsilon=epsilon,
        nu=nu,
        delta=delta,
        batchsize=batchsize,
        d=d,
        seed=seed,
    )
end


text:
@testset "Test Finite Elements" begin

using HDiscontinuousGalerkin
using Tensors
HDG = HDiscontinuousGalerkin

# Test ContinuousLagrange Element
finiteElement = ContinuousLagrange{2,RefTetrahedron,2}()
@test HDG.gettopology(finiteElement) == Dict(0=>3,1=>3,2=>0)
@test HDG.getgeomtopology(finiteElement) == Dict(0=>3,1=>0,2=>0)
x = [Vec{2}((1.0,1.0)),Vec{2}((2.0,1.0)),Vec{2}((2.0,2.0))]
y = [Vec{2}((1.0,1.0)),Vec{2}((2.0,1.0)),Vec{2}((2.0,2.0)),Vec{2}((1.5,1.5)),Vec{2}((1.5,1.0)),Vec{2}((1.0,1.5))]
for i in 1:3
    @test HDG.spatial_nodal_coordinate(finiteElement, i, x) == y[i]
end

finiteElement = ContinuousLagrange{2,RefTetrahedron,1}()
@test finiteElement.nodal_base_coefs ≈ [0.23570226039551578 0.2357022603955159 0.2357022603955158;
                                 -0.14433756729740646 0.14433756729740646 0.0; -1/12 -1/12 1/6]
@test HDG.gettopology(finiteElement) == Dict(0=>3,1=>0,2=>0)
for i in 1:3
    @test HDG.spatial_nodal_coordinate(finiteElement, i, x) == x[i]
end

# Test Lagrange Continuous Elements
function ref_value(i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 1 && return 1. - ξ_x - ξ_y
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

ξ_ref = Vec{2}((0.5,0.5))
for j = 1:3
    @test abs(ref_value(j,ξ_ref) - value(finiteElement,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value(j, ξ), ξ_ref) - gradient_value(finiteElement,j,ξ_ref), Inf) < eps(Float64)
end

function ref_value1(i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return 1-ξ_x
    i == 2 && return ξ_x
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end
ξ_ref = Vec{1}((0.5,))
finiteElement = ContinuousLagrange{1,RefTetrahedron,1}()
@test HDG.gettopology(finiteElement) == Dict(0=>2,1=>0)
for j = 1:2
    @test abs(ref_value1(j,ξ_ref) - value(finiteElement,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value1(j, ξ), ξ_ref) - gradient_value(finiteElement,j,ξ_ref), Inf) < eps(Float64)
end

end


text:
"""
Created on Fri May 28 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Algorithms for clustering data.
"""

"""
    build_conn(labels)

Build the connectivity matrix C of a list of labels of size N. C has size N x N and C_ij is 1 if i and j are have the
same label, 0 otherwise.
"""
function build_conn(labels::Array{Int64,1})
    indices_1 = repeat(transpose(labels), length(labels), 1)
    indices_2 = repeat(labels, 1, length(labels))
    convert(Matrix{Float64}, indices_1 .== indices_2)
end

"""
    _kmeans_cost(X, centroids, labels, dist_func)

Auxiliary function for kmeans_init.
"""
function _kmeans_cost(X::Matrix{T}, centroids::Matrix{T}, labels::Vector{Int64}, dist_func::SemiMetric) where T <: Real
    list_intra_dist = Array{Float64, 1}()

    for label in unique(labels)
        intra_dist = pairwise(dist_func, hcat(centroids[:,label], X[:, labels .== label]), dims=2)[1, 2:end]
        push!(list_intra_dist, mean(intra_dist))
    end

    mean(list_intra_dist)
end

"""
    kmeans_init(X, init, max_iter, dist_func, stopping_crit, stopping_tol, stopping_iter)

Run kmeans algorithm for a specified set of initial centroids. If `stopping_crit` is :dist, convergenced is reached when
the maximum distance between the new centroids and the ref centroids remains below `stopping_tol` for `stopping_iter`
consecutive iterations. If `stopping_crit` is :cost, convergence is reached  when  the relative change in the cost
remains below `stopping_tol` for `stopping_iter` consecutive iterations.
"""
function kmeans_init(;X::Matrix{T}, init::Vector{Int64}, max_iter::Integer, dist_func::SemiMetric, stopping_crit::Symbol, stopping_tol::Real, stopping_iter::Real, verbose::Bool=false) where T <: Real
    n_clus        = length(init)
    n_tot         = size(X, 2)
    centroids_ref = rand(T, size(X,1), n_clus)
    labels_ref    = zeros(Int64, n_tot)
    cost_ref      = 1e9
    n_iter_tol    = 0

    #### initial centroids
    centroids     = X[:, init]

    for i_iter = 1:max_iter
        #### compute distance of each data point to each centroid
        dist_to_centroids = pairwise(dist_func, hcat(centroids, X), dims=2)
        dist_to_centroids = transpose(dist_to_centroids[1:n_clus, (n_clus+1):(n_clus+n_tot)])

        #### update cluster labels using the specific structure of the matrix X
        #### The matrix X is an assembly of submatrices, each with n_clus columns, and
        #### labels are assigned so that for each submatrix, each column is assigned
        #### to one distinct label.
        labels = zeros(Int64, n_tot)

        for j in randperm(n_clus)
            for i in 1:n_clus:n_tot
                i_range = i:(i+n_clus-1)
                i_index = i + argmin(dist_to_centroids[i_range, j]) - 1
                labels[i_index] = j
                dist_to_centroids[i_index, :] .= Inf
            end
        end

        # labels = vec((x -> x[2]).(argmin(dist_to_centroids, dims=2)))

        #### update cluster centroids
        for label = 1:n_clus
            centroids[:, label] = mean(X[:, labels .== label], dims=2)
        end

        #### compute max distance between new and ref centroids
        dist_max = maximum(diag(pairwise(dist_func, centroids_ref, centroids, dims=2)))

        #### compute new cost
        cost = _kmeans_cost(X, centroids, labels, dist_func)

        #### check convergence
        if stopping_crit == :cost
            if abs((cost_ref-cost)/cost_ref) < stopping_tol
                n_iter_tol += 1
            else
                n_iter_tol = 0
                centroids_ref .= centroids
                labels_ref .= labels
                cost_ref = cost
            end
        elseif stopping_crit == :dist
            if dist_max < stopping_tol
                n_iter_tol += 1
            else
                n_iter_tol = 0
                centroids_ref .= centroids
                labels_ref .= labels
                cost_ref = cost
            end
        else
            throw(ArgumentError("Invalid `stopping_crit`"))
        end

        if n_iter_tol == stopping_iter
            if verbose
                @info "convergence at iter" i_iter
            end
            break
        end

        if i_iter == max_iter
            if verbose
                @info "maximum number of iterations reached" i_iter
            end
        end
    end

    centroids_ref, labels_ref, cost_ref
end

"""
    kmeans_cons(X, k_min, k_max, n_init, best_init, max_iter, dist_func, stopping_crit, stopping_tol, stopping_iter)

Build a consensus matrix from application of kmeans algorithm for a varying number of clusters. For each value of k,
kmeans is run for multiple initializations. If `best_init` is true, only the initialization resulting in the lowest cost
is used for building the consensus matrix for each value of k. Otherwise every run of kmeans is used for building the
consensus matrix. The X matrix should be in d × n format with d n the number of observations of a d-dimension vector of
characteristics.
"""
function kmeans_cons(;X::Matrix{T}, k_min::Integer, k_max::Integer, n_init::Integer, best_init::Bool, max_iter::Integer, dist_func::SemiMetric, stopping_crit::Symbol, stopping_tol::Real, stopping_iter::Real) where T <: Real
    list_conn = Array{Matrix{Float64}, 1}()

    for k = k_min:k_max
        list_cost = Array{Float64,1}()
        list_labels = Array{Array{Int64, 1}, 1}()

        for i_init = 1:n_init
            init = sample(1:size(X,2), k, replace=false)
            centroids, labels, cost = kmeans_init(
                X             = X,
                init          = init,
                max_iter      = max_iter,
                dist_func     = dist_func,
                stopping_crit = stopping_crit,
                stopping_tol  = stopping_tol,
                stopping_iter = stopping_iter
            )

            push!(list_cost, cost)
            push!(list_labels, labels)
        end

        if best_init
            labels_best = list_labels[argmin(list_cost)]
            conn = build_conn(labels_best)
            push!(list_conn, conn)
        else
            for labels in list_labels
                conn = build_conn(labels)
                push!(list_conn, conn)
            end
        end
    end

    mean(list_conn)
end

# using PyPlot
# using PyCall
# const plt = PyPlot
# const cm  = pyimport("matplotlib.cm")
#
# py"""
# from sklearn import datasets
#
# def noisy_circles(n_samples):
#     return datasets.make_circles(n_samples=n_samples, factor=.5, noise=.05)
#
# def noisy_moons(n_samples):
#     return datasets.make_moons(n_samples=n_samples, noise=.05)
# """
#
# X_moons, y_moons = py"noisy_moons"(100)
# X_circs, y_circs = py"noisy_circles"(100)
#
# fig, ax = plt.subplots(figsize=(14,10))
# ax.scatter(X_circs[y_circs .== 0, 1], X_circs[y_circs .== 0, 2], color="purple", label="true 0")
# ax.scatter(X_circs[y_circs .== 1, 1], X_circs[y_circs .== 1, 2], color="orange", label="true 0")
#
# X = convert(Array{Float64, 2}, transpose(X_circs))
#
# cons = kmeans_cons(
#     X             = X,
#     k_min         = 2,
#     k_max         = 15,
#     n_init        = 20,
#     best_init     = false,
#     max_iter      = 100,
#     dist_func     = Euclidean(),
#     stopping_crit = :dist,
#     stopping_tol  = 0.05,
#     stopping_iter = 100
# )
#
# hc = hclust(1 .- cons, linkage=:average, branchorder=:r)
# labels = cutree(hc, k=2)
#
# fig, ax = plt.subplots(figsize=(14,10))
# cmap = cm.get_cmap("tab20")
# for label in unique(labels)
#     ax.scatter(X[1, labels .== label], X[2, labels .== label], color=cmap(label))
# end
#
# plt.close()


text:
```Write a function that returns the largest element in a list.```

list = [1,2,3,4,5,6,7,46,9]

function max_n(x::Vector{Int64})
    maximum(x)
end

max_n(list)
maximum(list)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# Methods to update FDTD structs in several ways

"""
update the perturbation vector using the perturbed medium
in this case, medium will be treated as the background medium
* `δmod` is [δKI, δrhoI]
"""
function update_δmods!(pac::P_common, δmod::Vector{Float64})
    nx = pac.ic[:nx]
    nz = pac.ic[:nz]
    nznxd = prod(length.(pac.medium.mgrid))
    copyto!(pac.δmodall, δmod)
    fill!(pac.δmod[:KI], 0.0)
    δmodKI = view(pac.δmod[:KI], npml+1:nz-npml, npml+1:nx-npml)
    for i = 1:nznxd
        # put perturbation due to KI as it is
        δmodKI[i] = δmod[i]
    end
    fill!(pac.δmod[:rhoI], 0.0)
    δmodrr = view(pac.δmod[:rhoI], npml+1:nz-npml, npml+1:nx-npml)
    for i = 1:nznxd
        # put perturbation due to rhoI here
        δmodrr[i] = δmod[nznxd+i]
    end
    # project δmodrr onto the vz and vx grids
    get_rhovxI!(pac.δmod[:rhovxI], pac.δmod[:rhoI])
    get_rhovzI!(pac.δmod[:rhovzI], pac.δmod[:rhoI])
    return nothing
end

"""
This method should be executed only after the updating the main medium.
Update the `δmods` when a perturbed `medium_pert` is input.
The medium through which the waves are propagating
is assumed to be the background medium.
"""
function update_δmods!(pac::P_common, medium_pert::Medium)
    nznxd = prod(length.(pac.medium.mgrid))
    fill!(pac.δmodall, 0.0)
    copyto!(pac.δmodall, medium_pert, [:KI, :rhoI])

    for i = 1:nznxd
        pac.δmodall[i] -= pac.mod[:KI][i] # subtracting the background medium
        pac.δmodall[nznxd+i] -= pac.mod[:rhoI][i] # subtracting the background medium
    end
    update_δmods!(pac, pac.δmodall)
    return nothing
end

"""
```julia
update!(pa,medium_new)
```
Update `pa` with a new bundle of medium parameters `medium_new`, without additional memory allocation.
This routine is used during inversion, where medium parameters are iteratively updated.
The ability to iteratively run the forward mediuming task (with no additional memory allocation) on
various subsurface mediums is necessary while implementing inversion
algorithms.
"""
function update!(pa::PFdtd, medium::Medium)
    return update!(pa.c, medium)
end
function update!(pac::T, medium::Medium) where {T<:P_common{FdtdAcou}}
    copyto!(pac.medium, medium)
    padarray!(pac.exmedium, pac.medium, npml, pac.pml_edges)
    copyto!(pac.mod[:K], pac.exmedium, [:K])
    copyto!(pac.mod[:KI], pac.exmedium, [:KI])
    copyto!(pac.mod[:rhoI], pac.exmedium, [:rhoI])
    # get_rhovxI!(pac.mod[:rhovxI], pac.mod[:rhoI])
    # get_rhovzI!(pac.mod[:rhovzI], pac.mod[:rhoI])
    return nothing
end
function update!(pac::T, medium::Medium) where {T<:P_common{FdtdElastic}}
    copyto!(pac.medium, medium)
    padarray!(pac.exmedium, pac.medium, npml, pac.pml_edges)
    copyto!(pac.mod[:mu], pac.exmedium, [:mu])
    copyto!(pac.mod[:lambda], pac.exmedium, [:lambda])
    copyto!(pac.mod[:M], pac.exmedium, [:M])
    copyto!(pac.mod[:rho], pac.exmedium, [:rho])
    return nothing
end

"""
```julia
update!(pa,srcwav_new,sflags)
```
Update `pa` with a new bundle of source wavelets `srcwav_new`, without additional memory allocation.
Optionally, `sflags` can be changed.
"""
function update!(pa::PFdtd, srcwav::SrcWav, sflags::Any = nothing; verbose = false)
    update!(pa, [srcwav], sflags)
end
function update!(pa::PFdtd, srcwav::Vector{SrcWav}, sflags = nothing; verbose = false)
    # update srcwav in pa.c
    (length(srcwav) ≠ pa.c.ic[:npw]) && error("cannot update")
    for i = 1:length(srcwav)
        copyto!(pa.c.srcwav[i], srcwav[i])
    end
    if (!(sflags === nothing))
        copyto!(pa.c.sflags, sflags)
    end
    for ipw = 1:pa.c.ic[:npw]
        freqmin = 0.0
        freqmax = Inf
        freqpeaks = []
        # fill_wavelets for each supersource
        @sync begin
            for (ip, p) in enumerate(procs(pa.p))
                @async remotecall_wait(p) do
                    pap = localpart(pa.p)[ipw]
                    for is = 1:length(pap.ss)
                        iss = pap.ss[is].iss
                        wavelets = pap.ss[is].wavelets
                        broadcast(x -> fill!.(x, 0.0), wavelets)
                        fmin, fmax, fpeak =
                            fill_wavelets!(ipw, iss, wavelets, pa.c.srcwav, pa.c.sflags)
                        freqmin = max(freqmin, fmin)
                        freqmax = min(freqmax, fmax)
                        push!(freqpeaks, fpeak)
                    end
                end
            end
        end
        freqpeak = Statistics.mean(freqpeaks)
        # store frequency (in Hz) bounds for first propagating wavefield
        if (ipw == 1)
            for (f, nf) in
                zip([freqmin, freqmax, freqpeak], [:freqmin, :freqmax, :freqpeak])
                if (nf ∈ names(pa.c.fc)[1])
                    pa.c.fc[nf] = f
                else
                    pa.c.fc = vcat(pa.c.fc, NamedArray([f], [nf]))
                end
            end
        end
        if (verbose)
            freqmin = @sprintf("%0.2e", freqmin)
            freqmax = @sprintf("%0.2e", freqmax)
            freqpeak = @sprintf("%0.2e", freqpeak)
            @info "frequency bounds for propagating wavefield $ipw are: [$freqmin, $freqmax], with peak at: $freqpeak"
        end
    end
end

function update!(pass::P_x_worker_x_pw_x_ss, ipw, iss, ageomss::AGeomss, pac)
    @assert ageomss.ns == pac.ageom[ipw][iss].ns
    @assert ageomss.nr == pac.ageom[ipw][iss].nr

    ssprayw = pass.ssprayw
    rinterpolatew = pass.rinterpolatew

    for sfield in names(pass.sindices)[1]
        sindices=pass.sindices[sfield]
    for is = 1:ageomss.ns
        w = ssprayw[sfield][is]
        sindices[is], ww = get_indices_weights(eval(sfield)(),
            pac.exmedium.mgrid...,
            [s[is] for s in ageomss.s],
        )
        copyto!(w, ww)
    end
end
    for rfield in names(pass.rindices)[1]
    rindices = pass.rindices[rfield]
    for ir = 1:ageomss.nr
        w = rinterpolatew[rfield][ir]
        rindices[ir], ww = get_indices_weights(eval(rfield)(),
            pac.exmedium.mgrid...,
            [r[ir] for r in ageomss.r],
        )
        copyto!(w, ww)
    end
end

end

# if just one propagating field
update!(pa::PFdtd, ageom::AGeom) = update!(pa, [ageom])

function update!(pa::PFdtd, ageom::Vector{AGeom})
    for ipw = 1:pa.c.ic[:npw]
        copyto!(pa.c.ageom[ipw], ageom[ipw])
        @sync begin
            for (ip, p) in enumerate(procs(pa.p))
                @async remotecall_wait(p) do
                    pap = localpart(pa.p)[ipw]
                    for is = 1:length(pap.ss)
                        iss = pap.ss[is].iss
                        update!(localpart(pap).ss[is], ipw, iss, ageom[ipw][iss], pa.c)
                    end
                end
            end
        end
    end
end

"""
```julia
update!(pa)
```
In-place method to perform the experiment and update `pa` after wave propagation. After update, see
`pa[:data]` and `pa[:snaps]`.
"""
@fastmath function update!(pa::PFdtd)

    global to

    reset_timer!(to)

    # zero out all the results stored in pa.c
    initialize!(pa.c)

    # zero out results stored per worker
    @sync begin
        for (ip, p) in enumerate(procs(pa.p))
            @async remotecall_wait(p) do
                initialize!(localpart(pa.p))
            end
        end
    end

    # @timeit to "mod_x_proc!" begin
        # all localparts of DArray are input to this method
        # parallelization over shots
    @sync begin
        for (ip, p) in enumerate(procs(pa.p))
            @async remotecall_wait(p) do
                mod_x_proc!(pa.c, localpart(pa.p))
            end
        end
    end
    # end

    # stack gradients and illum over sources
    @sync begin
        for (ip, p) in enumerate(procs(pa.p))
            @sync remotecall_wait(p) do
                (pa.c.gmodel_flag) && stack_grads!(pa.c, localpart(pa.p))
                (pa.c.illum_flag) && stack_illums!(pa.c, localpart(pa.p))
            end
        end
    end

    # update gradient medium using grad_modKI_stack, grad_modrr_stack
    (pa.c.gmodel_flag) && update_gradient!(pa.c)

    for ipw in pa.c.activepw
        if (pa.c.rflags[ipw] ≠ 0) # record only if rflags is non-zero
            for rfield in pa.c.rfields

                fill!(pa.c.datamat, 0.0)
                @sync begin
                    for (ip, p) in enumerate(procs(pa.p))
                        @sync remotecall_wait(p) do
                            update_datamat!(rfield, ipw, pa.c, localpart(pa.p))
                        end
                    end
                end
                update_data!(rfield, ipw, pa.c)
            end
        end
    end
    if (pa.c.verbose)
        println("timer output of methods inside fdtd time loop")
        show(to)
        println("  ")
    end
    return nothing
end

# modelling for each worker
function mod_x_proc!(pac::P_common, pap::Vector{P_x_worker_x_pw{N,B}}) where {N,B}

    # source_loop
    for issp = 1:length(pap[1].ss) # note, all npw have same sources
        reset_w2!(pap)

        iss = pap[1].ss[issp].iss # same note as above

        # only for first propagating wavefield, i.e., pap[1]
        if (pac.backprop_flag == -1)
            "initial conditions from boundary for first propagating field only"
            boundary_force_snap_p!(issp, pac, pap)
            boundary_force_snap_vxvz!(issp, pac, pap)
        end

        prog = Progress(
            pac.ic[:nt],
            dt = 1,
            desc = "\tmodeling supershot $iss/$(length(pac.ageom[1])) ",
            color = :blue,
        )
        # time_loop
        """
        * don't use shared arrays inside this time loop, for speed when using multiple procs
        """
        for it = 1:pac.ic[:nt]

            pac.verbose && next!(prog, :white)

            @timeit to "advance time step" begin
                advance!(pac, pap)
            end

            # force p[1] on boundaries, only for ipw=1
            (pac.backprop_flag == -1) && boundary_force!(it, issp, pac, pap)

            @timeit to "add source" begin
                add_source!(it, issp, iss, pac, pap)
            end

            # no born flag for adjoint modelling
            if (!pac.gmodel_flag)
                (typeof(pac.attrib_mod) == FdtdAcouBorn) &&
                    add_born_sources!(issp, pac, pap)
            end

            # record boundaries after time reversal already
            (pac.backprop_flag == 1) && boundary_save!(pac.ic[:nt] - it + 1, issp, pac, pap)

            @timeit to "record at receivers" begin
                record!(it, issp, iss, pac, pap)
            end

            if(pac.gmodel_flag)
                @timeit to "compute gradient" begin
                    compute_gradient!(issp, pac, pap)
                end
            end

            (pac.illum_flag) && compute_illum!(issp, pap)

            if (eval(pac.snaps_field) <: Fields)
                iitsnaps = findall(x -> ==(x, it), pac.itsnaps)
                for itsnap in iitsnaps
                    snaps_save!(itsnap, issp, pac, pap)
                end
            end

        end # time_loop
        "now pressure is at [nt], velocities are at [nt-1/2]"

        "one more propagating step to save pressure at [nt+1] -- for time revarsal"
        advance!(pac, pap)

        "save last snap of pressure field"
        (pac.backprop_flag == 1) && boundary_save_snap_p!(issp, pac, pap)

        "one more propagating step to save velocities at [nt+3/2] -- for time reversal"
        advance!(pac, pap)

        "save last snap of velocity fields with opposite sign for adjoint propagation"
        (pac.backprop_flag == 1) && boundary_save_snap_vxvz!(issp, pac, pap)

        "scale gradients for each issp"
        (pac.gmodel_flag) && scale_gradient!(
            issp,
            pap,
            step(pac.medium.mgrid[2]) * step(pac.medium.mgrid[1]),
        )

    end # source_loop
end # mod_x_shot

function update_datamat!(
    rfield,
    ipw,
    pac::P_common,
    pap::Vector{P_x_worker_x_pw{N,B}},
) where {N,B}
    datamat = pac.datamat
    pass = pap[ipw].ss
    for issp = 1:length(pass)
        iss = pass[issp].iss
        records = pass[issp].records
        for ir = 1:pac.ageom[ipw][iss].nr
            for it = 1:pac.ic[:nt]
                datamat[it, ir, iss] = records[rfield][it, ir]
            end
        end
    end
end

function update_data!(rfield, ipw, pac::P_common)
    datamat = pac.datamat
    for iss = 1:length(pac.ageom[1])
        data = pac.data[ipw][iss].d[rfield]
        for ir = 1:pac.ageom[ipw][iss].nr
            for it = 1:pac.ic[:nt]
                data[it, ir] = datamat[it, ir, iss]
            end
        end
    end
end


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

## shell-like command parsing ##

function shell_parse(raw::AbstractString, interp::Bool)
    s = lstrip(raw)
    #Strips the end but respects the space when the string endswith "\\ "
    r = RevString(s)
    i = start(r)
    c_old = nothing
    while !done(r,i)
        c, j = next(r,i)
        if c == '\\' && c_old == ' '
            i -= 1
            break
        elseif !(c in _default_delims)
            break
        end
        i = j
        c_old = c
    end
    s = s[1:end-i+1]

    last_parse = 0:-1
    isempty(s) && return interp ? (Expr(:tuple,:()),last_parse) : ([],last_parse)

    in_single_quotes = false
    in_double_quotes = false

    args::Vector{Any} = []
    arg::Vector{Any} = []
    i = start(s)
    j = i

    function update_arg(x)
        if !isa(x,AbstractString) || !isempty(x)
            push!(arg, x)
        end
    end
    function append_arg()
        if isempty(arg); arg = Any["",]; end
        push!(args, arg)
        arg = []
    end

    while !done(s,j)
        c, k = next(s,j)
        if !in_single_quotes && !in_double_quotes && isspace(c)
            update_arg(s[i:j-1])
            append_arg()
            j = k
            while !done(s,j)
                c, k = next(s,j)
                if !isspace(c)
                    i = j
                    break
                end
                j = k
            end
        elseif interp && !in_single_quotes && c == '$'
            update_arg(s[i:j-1]); i = k; j = k
            if done(s,k)
                error("\$ right before end of command")
            end
            if isspace(s[k])
                error("space not allowed right after \$")
            end
            stpos = j
            ex, j = parse(s,j,greedy=false)
            last_parse = stpos:j
            update_arg(esc(ex)); i = j
        else
            if !in_double_quotes && c == '\''
                in_single_quotes = !in_single_quotes
                update_arg(s[i:j-1]); i = k
            elseif !in_single_quotes && c == '"'
                in_double_quotes = !in_double_quotes
                update_arg(s[i:j-1]); i = k
            elseif c == '\\'
                if in_double_quotes
                    if done(s,k)
                        error("unterminated double quote")
                    end
                    if s[k] == '"' || s[k] == '$' || s[k] == '\\'
                        update_arg(s[i:j-1]); i = k
                        c, k = next(s,k)
                    end
                elseif !in_single_quotes
                    if done(s,k)
                        error("dangling backslash")
                    end
                    update_arg(s[i:j-1]); i = k
                    c, k = next(s,k)
                end
            end
            j = k
        end
    end

    if in_single_quotes; error("unterminated single quote"); end
    if in_double_quotes; error("unterminated double quote"); end

    update_arg(s[i:end])
    append_arg()

    if !interp
        return (args,last_parse)
    end

    # construct an expression
    ex = Expr(:tuple)
    for arg in args
        push!(ex.args, Expr(:tuple, arg...))
    end
    (ex,last_parse)
end
shell_parse(s::AbstractString) = shell_parse(s,true)

function shell_split(s::AbstractString)
    parsed = shell_parse(s,false)[1]
    args = AbstractString[]
    for arg in parsed
       push!(args, string(arg...))
    end
    args
end

function print_shell_word(io::IO, word::AbstractString)
    if isempty(word)
        print(io, "''")
    end
    has_single = false
    has_special = false
    for c in word
        if isspace(c) || c=='\\' || c=='\'' || c=='"' || c=='$'
            has_special = true
            if c == '\''
                has_single = true
            end
        end
    end
    if !has_special
        print(io, word)
    elseif !has_single
        print(io, '\'', word, '\'')
    else
        print(io, '"')
        for c in word
            if c == '"' || c == '$'
                print(io, '\\')
            end
            print(io, c)
        end
        print(io, '"')
    end
end

function print_shell_escaped(io::IO, cmd::AbstractString, args::AbstractString...)
    print_shell_word(io, cmd)
    for arg in args
        print(io, ' ')
        print_shell_word(io, arg)
    end
end
print_shell_escaped(io::IO) = nothing

shell_escape(args::AbstractString...) = sprint(print_shell_escaped, args...)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Test, StatsAPI

@testset "StatsAPI" begin

# Nothing to test currently apart from checking that package can be loaded

end # @testset "StatsAPI"


text:
using ONNXML, Flux, ProtoBuf
using Base.Test

include("ops_tests.jl")

@testset "ONNXML" begin

include("conversions.jl")
include("constant.jl")
include("logical_ops.jl")
include("pooling.jl")
include("conv.jl")
include("reshape.jl")
include("arithmetic_ops.jl")
include("lstm.jl")

end