text:
###############################################################################
#
# Temperature dependency of the photosynthetic parameters
#
###############################################################################
"""
    photo_TD_from_set(td_set::ArrheniusTD{FT}, T::FT) where {FT<:AbstractFloat}
    photo_TD_from_set(td_set::Q10TD{FT}, T::FT) where {FT<:AbstractFloat}

Make temperature correction from parameter set, given
- `td_set` [`ArrheniusTD`](@ref) type parameter set, which has a `VAL_25` field
- `T` Leaf temperature

Useful for Kc, Ko, Kpep, and ``Γ^{*}``.
"""
function photo_TD_from_set(
            td_set::ArrheniusTD{FT},
            T::FT
) where {FT<:AbstractFloat}
    return td_set.VAL_25 * temperature_correction(td_set, T)
end

function photo_TD_from_set(
            td_set::Q10TD{FT},
            T::FT
) where {FT<:AbstractFloat}
    return td_set.VAL_REF * temperature_correction(td_set, T)
end

"""
    photo_TD_from_val(
                td_set::AbstractTDParameterSet{FT},
                val::FT,
                T::FT
    ) where {FT<:AbstractFloat}

Make temperature correction from a given value, given
- `td_set` [`ArrheniusTD`](@ref) or [`ArrheniusPeakTD`](@ref) type struct
- `val` Uncorrected value at 298.15 K
- `T` Leaf temperature

Useful for Vcmax, Vomax, Vpmax, Jmax, and Respiration.
"""
function photo_TD_from_val(
            td_set::AbstractTDParameterSet{FT},
            val::FT,
            T::FT
) where {FT<:AbstractFloat}
    return val * temperature_correction(td_set, T)
end

###############################################################################
#
# Functions to update the TD individually
#
###############################################################################
"""
    leaf_jmax!(td_set::AbstractTDParameterSet{FT},
               leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update maximal electron transport rate at leaf temperature, given
- `td_set` [`AbstractTDParameterSet`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_jmax!(
            td_set::AbstractTDParameterSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Jmax = photo_TD_from_val(td_set, leaf.Jmax25, leaf.T);

    return nothing
end

"""
    leaf_kc!(td_set::ArrheniusTD{FT}, leaf::Leaf{FT}) where {FT<:AbstractFloat}

Update Kc at leaf temperature, given
- `td_set` [`ArrheniusTD`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_kc!(
            td_set::ArrheniusTD{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Kc = photo_TD_from_set(td_set, leaf.T);

    return nothing
end

"""
    leaf_km!(photo_set::C3ParaSet{FT},
             leaf::Leaf{FT},
             envir::AirLayer{FT}
    ) where {FT<:AbstractFloat}

Update Ko at leaf temperature, given
- `photo_set` [`C3ParaSet`](@ref) type photosynthesis parameter set
- `leaf` [`Leaf`](@ref) type struct
- `envir` [`AirLayer`](@ref) type struct
"""
function leaf_km!(
            photo_set::C3ParaSet{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    leaf.Km = leaf.Kc * (1 + envir.p_O₂/leaf.Ko);

    return nothing
end

"""
    leaf_ko!(td_set::ArrheniusTD{FT}, leaf::Leaf{FT}) where {FT<:AbstractFloat}

Update Ko at leaf temperature, given
- `td_set` [`ArrheniusTD`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_ko!(
            td_set::ArrheniusTD{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Ko = photo_TD_from_set(td_set, leaf.T);

    return nothing
end

"""
    leaf_kpep!(td_set::ArrheniusTD{FT},
               leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update Kpep at leaf temperature, given
- `td_set` [`ArrheniusTD`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_kpep!(
            td_set::ArrheniusTD{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Kpep = photo_TD_from_set(td_set, leaf.T);

    return nothing
end

"""
    leaf_rd!(td_set::AbstractTDParameterSet{FT},
             leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update leaf dark respiration rate at leaf temperature, given
- `td_set` [`AbstractTDParameterSet`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_rd!(
            td_set::AbstractTDParameterSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Rd = photo_TD_from_val(td_set, leaf.Rd25, leaf.T);

    return nothing
end

"""
    leaf_vcmax!(td_set::AbstractTDParameterSet{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update leaf maximal carboxylation rate at leaf temperature, given
- `td_set` [`AbstractTDParameterSet`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_vcmax!(
            td_set::AbstractTDParameterSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Vcmax = photo_TD_from_val(td_set, leaf.Vcmax25, leaf.T);

    return nothing
end

"""
    leaf_vpmax!(td_set::AbstractTDParameterSet{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update leaf maximal PEP carboxylation rate at leaf temperature, given
- `td_set` [`AbstractTDParameterSet`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_vpmax!(
            td_set::AbstractTDParameterSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Vpmax = photo_TD_from_val(td_set, leaf.Vpmax25, leaf.T);

    return nothing
end

"""
    leaf_Γstar!(td_set::ArrheniusTD{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update ``Γ^{*}`` at leaf temperature, given
- `td_set` [`ArrheniusTD`](@ref) type TD parameter set
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_Γstar!(
            td_set::ArrheniusTD{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Γ_star = photo_TD_from_set(td_set, leaf.T);

    return nothing
end

###############################################################################
#
# Calculate photosynthesis using Leaf
#
###############################################################################
"""
    leaf_temperature_dependence!(
                photo_set::AbstractPhotoModelParaSet{FT},
                leaf::Leaf{FT},
                envir::AirLayer{FT}
    ) where {FT<:AbstractFloat}
    leaf_temperature_dependence!(
                photo_set::AbstractPhotoModelParaSet{FT},
                leaf::Leaf{FT},
                envir::AirLayer{FT},
                T::FT
    ) where {FT<:AbstractFloat}

Update the temperature dependent photosynthesis only, given
- `photo_set` [`AbstractPhotoModelParaSet`](@ref) type parameter set
- `leaf` [`Leaf`](@ref) type struct
- `envir` [`AirLayer`](@ref) type struct
- `T` Given leaf temperature
"""
function leaf_temperature_dependence!(
            photo_set::C3ParaSet{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    if leaf.T_old != leaf.T
        leaf.T_old = leaf.T;
        leaf.p_sat = saturation_vapor_pressure(leaf.T);
        leaf_rd!(photo_set.ReT, leaf);
        leaf_vcmax!(photo_set.VcT, leaf);
        leaf_jmax!(photo_set.JT , leaf);
        leaf_kc!(photo_set.KcT, leaf);
        leaf_ko!(photo_set.KoT, leaf);
        leaf_km!(photo_set, leaf, envir);
        leaf_Γstar!(photo_set.ΓsT, leaf);
    end

    return nothing
end

function leaf_temperature_dependence!(
            photo_set::C4ParaSet{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    if leaf.T_old != leaf.T
        leaf.T_old = leaf.T
        leaf.p_sat = saturation_vapor_pressure(leaf.T);
        leaf_rd!(photo_set.ReT, leaf);
        leaf_vcmax!(photo_set.VcT, leaf);
        leaf_vpmax!(photo_set.VpT, leaf);
        leaf_kpep!(photo_set.KpT, leaf);
    end

    return nothing
end

function leaf_temperature_dependence!(
            photo_set::AbstractPhotoModelParaSet{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT},
            T::FT
) where {FT<:AbstractFloat}
    leaf.T = T;
    leaf_temperature_dependence!(photo_set, leaf, envir);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module BenchFindfirst

using BenchmarkTools
using Folds
using FoldsThreads

const SUITE = BenchmarkGroup()
const DATA_LENGTH = 2^20

for param in [
    (basesize = 2^10, needleloc = 2^18),
    (basesize = 2^10, needleloc = 2^19),
    (basesize = 2^12, needleloc = 2^20),
    (basesize = 2^14, needleloc = 2^20),
]
    basesize = param.basesize
    needleloc = param.needleloc
    xs = rand(DATA_LENGTH)
    xs[needleloc] = 2
    SUITE[param] =
        @benchmarkable(Folds.findfirst(>(1), $xs, WorkStealingEx(basesize = $basesize)))
end
SUITE

end  # module
BenchFindfirst.SUITE

ext: jl


text:
module PWscf

using AbInitioSoftwareBase.Inputs: asstring
using StructArrays: StructArray
using Test
using QuantumESPRESSOBase.Inputs.PWscf
using QuantumESPRESSOFormatter

@testset "Test constructing `AtomicSpeciesCard`" begin
    a = ["Al", "As"]
    m = [24590.7655930491, 68285.4024548272]
    pp = ["Al.pbe-n-kjpaw_psl.1.0.0.UPF", "As.pbe-n-kjpaw_psl.1.0.0.UPF"]
    card = AtomicSpeciesCard(StructArray{AtomicSpecies}((a, m, pp)))
    @test asstring(card) ==
          "ATOMIC_SPECIES\n      Al 24590.765593049 Al.pbe-n-kjpaw_psl.1.0.0.UPF\n      As 68285.402454827 As.pbe-n-kjpaw_psl.1.0.0.UPF"
end

@testset "Test constructing `AtomicPositionsCard`" begin
    # Data from https://github.com/QEF/q-e/blob/7be27df/PW/examples/gatefield/run_example#L129-L132.
    a = ["S", "Mo", "S"]
    pos = [
        [0.500000000, 0.288675130, 1.974192764],
        [0.000000000, 0.577350270, 2.462038339],
        [0.000000000, -0.577350270, 2.950837559],
    ]
    card = AtomicPositionsCard(
        StructArray{AtomicPosition}((a, pos, [[1, 1, 1], [1, 1, 1], [1, 1, 1]])),
        "alat",
    )
    @test asstring(card) ==
          "ATOMIC_POSITIONS { alat }\n       S    0.500000000    0.288675130    1.974192764\n      Mo    0.000000000    0.577350270    2.462038339\n       S    0.000000000   -0.577350270    2.950837559"
end

end


text:
let n
for n in [2, 4, 16]
    H = matrixdepot("hadamard", Int, n)
    @test H*H' == n * Matrix(1.0I, n, n)
end

@test_throws ArgumentError matrixdepot("hadamard", 0)
end
println("'hadamard' passed test...")


text:
# # Power Flow

# `PowerSystems.jl` provides the capability to run a power flow with the intention of providing a valid initial AC operating point to the system. It is not meant as an analytics tool; the main issue is to determine if the system has feasible AC data.

# The power flow solver uses [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) under the hood and takes any keyword argument accepted by NLsolve

# ```@autodocs
# Modules = [PowerSystems]
# Pages = ["power_flow.jl"]
# Private = false
# ```


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

# Skipping MacroDefinition: AFAPI __attribute__ ( ( visibility ( "default" ) ) )

#const SIZE_T_FRMT_SPECIFIER = "%zu"

# Skipping MacroDefinition: DEPRECATED ( msg ) __attribute__ ( ( deprecated ) )

#const AF_VERSION = "3.2.2"
#const AF_VERSION_MAJOR = 3
#const AF_VERSION_MINOR = 2
#const AF_VERSION_PATCH = 2
#const AF_API_VERSION = AF_API_VERSION_CURRENT
const AF_MAX_DIMS = 4

typealias dim_t Clonglong
typealias intl Clonglong
typealias uintl Culonglong

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const AF_SUCCESS = (UInt32)(0)
const AF_ERR_NO_MEM = (UInt32)(101)
const AF_ERR_DRIVER = (UInt32)(102)
const AF_ERR_RUNTIME = (UInt32)(103)
const AF_ERR_INVALID_ARRAY = (UInt32)(201)
const AF_ERR_ARG = (UInt32)(202)
const AF_ERR_SIZE = (UInt32)(203)
const AF_ERR_TYPE = (UInt32)(204)
const AF_ERR_DIFF_TYPE = (UInt32)(205)
const AF_ERR_BATCH = (UInt32)(207)
const AF_ERR_NOT_SUPPORTED = (UInt32)(301)
const AF_ERR_NOT_CONFIGURED = (UInt32)(302)
const AF_ERR_NONFREE = (UInt32)(303)
const AF_ERR_NO_DBL = (UInt32)(401)
const AF_ERR_NO_GFX = (UInt32)(402)
const AF_ERR_LOAD_LIB = (UInt32)(501)
const AF_ERR_LOAD_SYM = (UInt32)(502)
const AF_ERR_ARR_BKND_MISMATCH = (UInt32)(503)
const AF_ERR_INTERNAL = (UInt32)(998)
const AF_ERR_UNKNOWN = (UInt32)(999)
# end enum ANONYMOUS_1

# begin enum af_err
typealias af_err UInt32
const AF_SUCCESS = (UInt32)(0)
const AF_ERR_NO_MEM = (UInt32)(101)
const AF_ERR_DRIVER = (UInt32)(102)
const AF_ERR_RUNTIME = (UInt32)(103)
const AF_ERR_INVALID_ARRAY = (UInt32)(201)
const AF_ERR_ARG = (UInt32)(202)
const AF_ERR_SIZE = (UInt32)(203)
const AF_ERR_TYPE = (UInt32)(204)
const AF_ERR_DIFF_TYPE = (UInt32)(205)
const AF_ERR_BATCH = (UInt32)(207)
const AF_ERR_NOT_SUPPORTED = (UInt32)(301)
const AF_ERR_NOT_CONFIGURED = (UInt32)(302)
const AF_ERR_NONFREE = (UInt32)(303)
const AF_ERR_NO_DBL = (UInt32)(401)
const AF_ERR_NO_GFX = (UInt32)(402)
const AF_ERR_LOAD_LIB = (UInt32)(501)
const AF_ERR_LOAD_SYM = (UInt32)(502)
const AF_ERR_ARR_BKND_MISMATCH = (UInt32)(503)
const AF_ERR_INTERNAL = (UInt32)(998)
const AF_ERR_UNKNOWN = (UInt32)(999)
# end enum af_err

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const f32 = (UInt32)(0)
const c32 = (UInt32)(1)
const f64 = (UInt32)(2)
const c64 = (UInt32)(3)
const b8 = (UInt32)(4)
const s32 = (UInt32)(5)
const u32 = (UInt32)(6)
const u8 = (UInt32)(7)
const s64 = (UInt32)(8)
const u64 = (UInt32)(9)
const s16 = (UInt32)(10)
const u16 = (UInt32)(11)
# end enum ANONYMOUS_2

# begin enum af_dtype
typealias af_dtype UInt32
const f32 = (UInt32)(0)
const c32 = (UInt32)(1)
const f64 = (UInt32)(2)
const c64 = (UInt32)(3)
const b8 = (UInt32)(4)
const s32 = (UInt32)(5)
const u32 = (UInt32)(6)
const u8 = (UInt32)(7)
const s64 = (UInt32)(8)
const u64 = (UInt32)(9)
const s16 = (UInt32)(10)
const u16 = (UInt32)(11)
# end enum af_dtype

# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const afDevice = (UInt32)(0)
const afHost = (UInt32)(1)
# end enum ANONYMOUS_3

# begin enum af_source
typealias af_source UInt32
const afDevice = (UInt32)(0)
const afHost = (UInt32)(1)
# end enum af_source

typealias af_array Ptr{Void}

# begin enum ANONYMOUS_4
typealias ANONYMOUS_4 UInt32
const AF_INTERP_NEAREST = (UInt32)(0)
const AF_INTERP_LINEAR = (UInt32)(1)
const AF_INTERP_BILINEAR = (UInt32)(2)
const AF_INTERP_CUBIC = (UInt32)(3)
const AF_INTERP_LOWER = (UInt32)(4)
# end enum ANONYMOUS_4

# begin enum af_interp_type
typealias af_interp_type UInt32
const AF_INTERP_NEAREST = (UInt32)(0)
const AF_INTERP_LINEAR = (UInt32)(1)
const AF_INTERP_BILINEAR = (UInt32)(2)
const AF_INTERP_CUBIC = (UInt32)(3)
const AF_INTERP_LOWER = (UInt32)(4)
# end enum af_interp_type

# begin enum ANONYMOUS_5
typealias ANONYMOUS_5 UInt32
const AF_PAD_ZERO = (UInt32)(0)
const AF_PAD_SYM = (UInt32)(1)
# end enum ANONYMOUS_5

# begin enum af_border_type
typealias af_border_type UInt32
const AF_PAD_ZERO = (UInt32)(0)
const AF_PAD_SYM = (UInt32)(1)
# end enum af_border_type

# begin enum ANONYMOUS_6
typealias ANONYMOUS_6 UInt32
const AF_CONNECTIVITY_4 = (UInt32)(4)
const AF_CONNECTIVITY_8 = (UInt32)(8)
# end enum ANONYMOUS_6

# begin enum af_connectivity
typealias af_connectivity UInt32
const AF_CONNECTIVITY_4 = (UInt32)(4)
const AF_CONNECTIVITY_8 = (UInt32)(8)
# end enum af_connectivity

# begin enum ANONYMOUS_7
typealias ANONYMOUS_7 UInt32
const AF_CONV_DEFAULT = (UInt32)(0)
const AF_CONV_EXPAND = (UInt32)(1)
# end enum ANONYMOUS_7

# begin enum af_conv_mode
typealias af_conv_mode UInt32
const AF_CONV_DEFAULT = (UInt32)(0)
const AF_CONV_EXPAND = (UInt32)(1)
# end enum af_conv_mode

# begin enum ANONYMOUS_8
typealias ANONYMOUS_8 UInt32
const AF_CONV_AUTO = (UInt32)(0)
const AF_CONV_SPATIAL = (UInt32)(1)
const AF_CONV_FREQ = (UInt32)(2)
# end enum ANONYMOUS_8

# begin enum af_conv_domain
typealias af_conv_domain UInt32
const AF_CONV_AUTO = (UInt32)(0)
const AF_CONV_SPATIAL = (UInt32)(1)
const AF_CONV_FREQ = (UInt32)(2)
# end enum af_conv_domain

# begin enum ANONYMOUS_9
typealias ANONYMOUS_9 UInt32
const AF_SAD = (UInt32)(0)
const AF_ZSAD = (UInt32)(1)
const AF_LSAD = (UInt32)(2)
const AF_SSD = (UInt32)(3)
const AF_ZSSD = (UInt32)(4)
const AF_LSSD = (UInt32)(5)
const AF_NCC = (UInt32)(6)
const AF_ZNCC = (UInt32)(7)
const AF_SHD = (UInt32)(8)
# end enum ANONYMOUS_9

# begin enum af_match_type
typealias af_match_type UInt32
const AF_SAD = (UInt32)(0)
const AF_ZSAD = (UInt32)(1)
const AF_LSAD = (UInt32)(2)
const AF_SSD = (UInt32)(3)
const AF_ZSSD = (UInt32)(4)
const AF_LSSD = (UInt32)(5)
const AF_NCC = (UInt32)(6)
const AF_ZNCC = (UInt32)(7)
const AF_SHD = (UInt32)(8)
# end enum af_match_type

# begin enum ANONYMOUS_10
typealias ANONYMOUS_10 UInt32
const AF_YCC_601 = (UInt32)(601)
const AF_YCC_709 = (UInt32)(709)
const AF_YCC_2020 = (UInt32)(2020)
# end enum ANONYMOUS_10

# begin enum af_ycc_std
typealias af_ycc_std UInt32
const AF_YCC_601 = (UInt32)(601)
const AF_YCC_709 = (UInt32)(709)
const AF_YCC_2020 = (UInt32)(2020)
# end enum af_ycc_std

# begin enum ANONYMOUS_11
typealias ANONYMOUS_11 UInt32
const AF_GRAY = (UInt32)(0)
const AF_RGB = (UInt32)(1)
const AF_HSV = (UInt32)(2)
const AF_YCbCr = (UInt32)(3)
# end enum ANONYMOUS_11

# begin enum af_cspace_t
typealias af_cspace_t UInt32
const AF_GRAY = (UInt32)(0)
const AF_RGB = (UInt32)(1)
const AF_HSV = (UInt32)(2)
const AF_YCbCr = (UInt32)(3)
# end enum af_cspace_t

# begin enum ANONYMOUS_12
typealias ANONYMOUS_12 UInt32
const AF_MAT_NONE = (UInt32)(0)
const AF_MAT_TRANS = (UInt32)(1)
const AF_MAT_CTRANS = (UInt32)(2)
const AF_MAT_CONJ = (UInt32)(4)
const AF_MAT_UPPER = (UInt32)(32)
const AF_MAT_LOWER = (UInt32)(64)
const AF_MAT_DIAG_UNIT = (UInt32)(128)
const AF_MAT_SYM = (UInt32)(512)
const AF_MAT_POSDEF = (UInt32)(1024)
const AF_MAT_ORTHOG = (UInt32)(2048)
const AF_MAT_TRI_DIAG = (UInt32)(4096)
const AF_MAT_BLOCK_DIAG = (UInt32)(8192)
# end enum ANONYMOUS_12

# begin enum af_mat_prop
typealias af_mat_prop UInt32
const AF_MAT_NONE = (UInt32)(0)
const AF_MAT_TRANS = (UInt32)(1)
const AF_MAT_CTRANS = (UInt32)(2)
const AF_MAT_CONJ = (UInt32)(4)
const AF_MAT_UPPER = (UInt32)(32)
const AF_MAT_LOWER = (UInt32)(64)
const AF_MAT_DIAG_UNIT = (UInt32)(128)
const AF_MAT_SYM = (UInt32)(512)
const AF_MAT_POSDEF = (UInt32)(1024)
const AF_MAT_ORTHOG = (UInt32)(2048)
const AF_MAT_TRI_DIAG = (UInt32)(4096)
const AF_MAT_BLOCK_DIAG = (UInt32)(8192)
# end enum af_mat_prop

# begin enum ANONYMOUS_13
typealias ANONYMOUS_13 UInt32
const AF_NORM_VECTOR_1 = (UInt32)(0)
const AF_NORM_VECTOR_INF = (UInt32)(1)
const AF_NORM_VECTOR_2 = (UInt32)(2)
const AF_NORM_VECTOR_P = (UInt32)(3)
const AF_NORM_MATRIX_1 = (UInt32)(4)
const AF_NORM_MATRIX_INF = (UInt32)(5)
const AF_NORM_MATRIX_2 = (UInt32)(6)
const AF_NORM_MATRIX_L_PQ = (UInt32)(7)
const AF_NORM_EUCLID = (UInt32)(2)
# end enum ANONYMOUS_13

# begin enum af_norm_type
typealias af_norm_type UInt32
const AF_NORM_VECTOR_1 = (UInt32)(0)
const AF_NORM_VECTOR_INF = (UInt32)(1)
const AF_NORM_VECTOR_2 = (UInt32)(2)
const AF_NORM_VECTOR_P = (UInt32)(3)
const AF_NORM_MATRIX_1 = (UInt32)(4)
const AF_NORM_MATRIX_INF = (UInt32)(5)
const AF_NORM_MATRIX_2 = (UInt32)(6)
const AF_NORM_MATRIX_L_PQ = (UInt32)(7)
const AF_NORM_EUCLID = (UInt32)(2)
# end enum af_norm_type

# begin enum ANONYMOUS_14
typealias ANONYMOUS_14 UInt32
const AF_COLORMAP_DEFAULT = (UInt32)(0)
const AF_COLORMAP_SPECTRUM = (UInt32)(1)
const AF_COLORMAP_COLORS = (UInt32)(2)
const AF_COLORMAP_RED = (UInt32)(3)
const AF_COLORMAP_MOOD = (UInt32)(4)
const AF_COLORMAP_HEAT = (UInt32)(5)
const AF_COLORMAP_BLUE = (UInt32)(6)
# end enum ANONYMOUS_14

# begin enum af_colormap
typealias af_colormap UInt32
const AF_COLORMAP_DEFAULT = (UInt32)(0)
const AF_COLORMAP_SPECTRUM = (UInt32)(1)
const AF_COLORMAP_COLORS = (UInt32)(2)
const AF_COLORMAP_RED = (UInt32)(3)
const AF_COLORMAP_MOOD = (UInt32)(4)
const AF_COLORMAP_HEAT = (UInt32)(5)
const AF_COLORMAP_BLUE = (UInt32)(6)
# end enum af_colormap

# begin enum ANONYMOUS_15
typealias ANONYMOUS_15 UInt32
const AF_FIF_BMP = (UInt32)(0)
const AF_FIF_ICO = (UInt32)(1)
const AF_FIF_JPEG = (UInt32)(2)
const AF_FIF_JNG = (UInt32)(3)
const AF_FIF_PNG = (UInt32)(13)
const AF_FIF_PPM = (UInt32)(14)
const AF_FIF_PPMRAW = (UInt32)(15)
const AF_FIF_TIFF = (UInt32)(18)
const AF_FIF_PSD = (UInt32)(20)
const AF_FIF_HDR = (UInt32)(26)
const AF_FIF_EXR = (UInt32)(29)
const AF_FIF_JP2 = (UInt32)(31)
const AF_FIF_RAW = (UInt32)(34)
# end enum ANONYMOUS_15

# begin enum af_image_format
typealias af_image_format UInt32
const AF_FIF_BMP = (UInt32)(0)
const AF_FIF_ICO = (UInt32)(1)
const AF_FIF_JPEG = (UInt32)(2)
const AF_FIF_JNG = (UInt32)(3)
const AF_FIF_PNG = (UInt32)(13)
const AF_FIF_PPM = (UInt32)(14)
const AF_FIF_PPMRAW = (UInt32)(15)
const AF_FIF_TIFF = (UInt32)(18)
const AF_FIF_PSD = (UInt32)(20)
const AF_FIF_HDR = (UInt32)(26)
const AF_FIF_EXR = (UInt32)(29)
const AF_FIF_JP2 = (UInt32)(31)
const AF_FIF_RAW = (UInt32)(34)
# end enum af_image_format

# begin enum ANONYMOUS_16
typealias ANONYMOUS_16 UInt32
const AF_HOMOGRAPHY_RANSAC = (UInt32)(0)
const AF_HOMOGRAPHY_LMEDS = (UInt32)(1)
# end enum ANONYMOUS_16

# begin enum af_homography_type
typealias af_homography_type UInt32
const AF_HOMOGRAPHY_RANSAC = (UInt32)(0)
const AF_HOMOGRAPHY_LMEDS = (UInt32)(1)
# end enum af_homography_type

# begin enum ANONYMOUS_17
typealias ANONYMOUS_17 UInt32
const AF_BACKEND_DEFAULT = (UInt32)(0)
const AF_BACKEND_CPU = (UInt32)(1)
const AF_BACKEND_CUDA = (UInt32)(2)
const AF_BACKEND_OPENCL = (UInt32)(4)
# end enum ANONYMOUS_17

# begin enum af_backend
typealias af_backend UInt32
const AF_BACKEND_DEFAULT = (UInt32)(0)
const AF_BACKEND_CPU = (UInt32)(1)
const AF_BACKEND_CUDA = (UInt32)(2)
const AF_BACKEND_OPENCL = (UInt32)(4)
# end enum af_backend

# begin enum ANONYMOUS_18
typealias ANONYMOUS_18 UInt32
const AF_ID = (UInt32)(0)
# end enum ANONYMOUS_18

# begin enum af_someenum_t
typealias af_someenum_t UInt32
const AF_ID = (UInt32)(0)
# end enum af_someenum_t

immutable af_seq
    _begin::Cdouble
    _end::Cdouble
    step::Cdouble
end

immutable af_index_t
    idx::Void
    isSeq::Bool
    isBatch::Bool
end

immutable af_cfloat
    real::Cfloat
    imag::Cfloat
end

immutable af_cdouble
    real::Cdouble
    imag::Cdouble
end

typealias af_features Ptr{Void}
typealias af_window Culonglong

immutable af_cell
    row::Cint
    col::Cint
    title::Ptr{UInt8}
    cmap::af_colormap
end


text:
module HCE_ZBlock_Module

import ..ED
using ..SimLib
using ..SimLib: FArray
using LinearAlgebra: eigvals!, Hermitian, mul!, svdvals!
using SharedArrays: sdata
using SpinSymmetry
using XXZNumerics: entropy

export HalfChainEntropyZBlock, HCEDataDescriptor, HCEData, load_entropy

### Descriptor

struct HCEDataDescriptor <: ED.EDDerivedDataDescriptor
    L::Int
    symm::Bool
    derivedfrom::ED.EDDataDescriptor
end

HCEDataDescriptor(L, args...; kwargs...) = HCEDataDescriptor(L, true, EDDataDescriptor(args...; kwargs...))
HCEDataDescriptor(L, symm::Bool, args...; kwargs...) = HCEDataDescriptor(L, symm, EDDataDescriptor(args...; kwargs...))

### Data obj

struct HCEData{N} <: SimLib.AbstractSimpleData
    descriptor::HCEDataDescriptor
    data::FArray{N}
end

ED._default_folder(::HCEDataDescriptor) = "entropy"
ED._filename_addition(hcedd::HCEDataDescriptor) = "-l_$(hcedd.L)" * (hcedd.symm ? "_symm" : "")

"""
    load_entropy(L[, symm], edd)
    load_entropy(L[, symm], model[, diagtype][, location])
"""

load_entropy(args...; kwargs...) =  load(HCEDataDescriptor(args...; kwargs...))

## Functions

_zblock_inds(N, k) = SpinSymmetry._indices(zbasis(N, k))

struct SymmZBlockEntanglementEntropy
    Nfull::Int
    kfull::Int
    N1::Int
    indexLookup::Vector{Int}## ToDo: is that a sensible data structure?
    indsA::Vector{Vector{Int}} # note that indsA contains fewer indices -> should be used as column
    indsB::Vector{Vector{Int}}
    size::Int
    function SymmZBlockEntanglementEntropy(zblockbasis, N1)
        N = zblockbasis.N
        k = zblockbasis.k
        N1 = max(N1, N-N1)
        indexLookup = zeros(2^N)
        indexLookup[_zblock_inds(N, k)] = 1:binomial(N,k)
        krange = max(0,N1+k-N):min(N1, k)
        indsA = [2^N1 .* (_zblock_inds(N-N1, k-ki) .- 1) for ki in krange]
        indsB = [_zblock_inds(N1, ki) .- 1 for ki in krange]
        size = 2*N1 == N ? N1 : N
        new(N, k, N1, indexLookup, indsA, indsB, size)
    end
end

# scalar for broadcasting
Base.broadcastable(s::SymmZBlockEntanglementEntropy) = Ref(s)

entanglement_entropy(s::SymmZBlockEntanglementEntropy, ψ) = entanglement_entropy!(zeros(Float64, s.size),s,ψ)

function entanglement_entropy!(out, s::SymmZBlockEntanglementEntropy, ψ)
    NB = s.Nfull
    fill!(out, 0)
    for (indA, indB) in zip(s.indsA, s.indsB)
        ## No point in optimizing this allocation further.
        mat = Matrix{eltype(ψ)}(undef, length(indA), length(indB))
        for shift in 1:s.size
            mat .= getindex.(Ref(ψ), getindex.(Ref(s.indexLookup), SpinSymmetry._roll_bits.(NB, indA .+ indB', shift-1) .+ 1))
            out[shift] += entropy(svdvals!(mat) .^ 2) #  most allocations come from svdvals!
		end
	end
    out
end

### Task

mutable struct HalfChainEntropyTask{S} <: ED.EDTask
    L::Int
    symm::Bool
    entropy_strategy::S
    data
end

## ToDo: This always symmetrizes over the chain right now.

HalfChainEntropyZBlock(basis::SymmetrizedBasis, L=div(basis.basis.N,2)) = HalfChainEntropyTask(L, true, SymmZBlockEntanglementEntropy(basis.basis, L), nothing)
HalfChainEntropyZBlock(; basis, L=div(basis.basis.N,2)) = HalfChainEntropyZBlock(basis, L)

function ED.initialize!(task::HalfChainEntropyTask, arrayconstructor, spectral_size)
    task.data = arrayconstructor(Float64, task.entropy_strategy.size, spectral_size)
end

function ED.compute_task!(task::HalfChainEntropyTask, evals, evecs, inds...)
    n = min(size(task.data,1), size(evecs,2))
    for (i, ψ) in enumerate(eachcol(evecs))
        i <= n || break
        entanglement_entropy!(view(task.data, :, i, inds...), task.entropy_strategy, ψ)
    end
end

function ED.failed_task!(task::HalfChainEntropyTask, inds...)
    task.data[:, :, inds...] .= NaN64
end

function ED.assemble(task::HalfChainEntropyTask, edd)
    HCEData(HCEDataDescriptor(task.L, task.symm, edd), sdata(task.data))
end

Base.summary(task::HalfChainEntropyTask) = string(typeof(task)) * "(L=$(task.L), symm=$(task.symm))"

end # module


text:
# load ImpvolEquilibrium first so that methods are accessible
@everywhere include("../../../equilibrium.jl")
@everywhere using ImpvolEquilibrium, Logging
@everywhere Logging.configure(level=INFO)

@everywhere using FileIO
@everywhere parameters = load("../common_parameters.jld2")["parameters"]

# parameters that govern counterfactual
@everywhere include("change_parameters.jl")

@everywhere srand(7094)

T = 3
@time results = pmap(t -> (t, ImpvolEquilibrium.period_wrapper(parameters, t)), 1:T)
GDP = zeros(T,)
for t=1:T
	GDP[t] = sum(results[t][2][:real_GDP][1,end,:,1])
end
display(GDP)


text:
function fib(n)
  if n < 1
    0
  elseif n < 2
    1
  else
    fib(n - 1) + fib(n - 2)
  end
end

@time fib(big"33")

max_stars_repo_name: kmizu/mlcombinators


text:
using OrdinaryDiffEq, Flux, Plots
using DifferentialEquations
using ForwardDiff
using LinearAlgebra
using BSON: @save
using BSON: @load
LinearAlgebra.BLAS.set_num_threads(16)

pz = [4.0e1, 3.0e2, 1.0e2]
u0 = [1.0,0.0,0.0]
tspan = (0,1.0)

function rober(du,u,p,t)
    y₁,y₂,y₃ = u
    k₁,k₂,k₃ = p
    du[1] = -k₁*y₁+k₃*y₂*y₃
    du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
    du[3] =  k₂*y₂^2
    nothing
end

solver = Rosenbrock23(autodiff=false,diff_type=Val{:forward})
prob = ODEProblem(rober, u0, tspan, pz)
sol = solve(prob)
sol.t

#Create the NN
dudt2 = Chain(Dense(1,8, gelu),
             Dense(8, 8, gelu),
             Dense(8, 8, gelu),
             Dense(8, 3))

p,re = Flux.destructure(dudt2) # use this p as the initial condition!
dudt(u,p) = re(p)(u)

function grab_∇u(p)
    reduce(hcat,[ForwardDiff.jacobian(u-> dudt(u,p), [t]) for t in sol.t])
end
function grab_pred(p)
    reduce(hcat, [dudt([t],p) for t in sol.t])
end
#test grab_∇u
grab_∇u(p)
grab_pred(p)

function loss_PINNS(p)
   #loss for stiff-PINNS
   uₛ = grab_pred(p)
   u₁ = uₛ[1, :]
   u₂ = uₛ[2, :]
   u₃ = uₛ[3, :]
   k₁,k₂,k₃ = pz
   duₛ = grab_∇u(p)
   du₁ = duₛ[1, :]
   du₂ = duₛ[2, :]
   du₃ = duₛ[3, :]
   e1 = abs.(du₁ .- k₁*u₁ .+ k₃.*u₂.*u₃)
   e2 = abs.(du₂ .- k₁*u₁ .- k₂*u₂.^2 .- k₃.*u₂.*u₃)
   e3 = abs.(du₃ .- k₂*u₂.^2)
   #loss_phys = sum(abs, e1 + e2 + e3)
   loss_ic = sum(abs, uₛ[:, 1] - u0)
   #total_loss= loss_phys
   total_loss = loss_ic
   return total_loss
end
#test loss_PINNS
loss_PINNS(p*1e-4)

function grab_∇L(par)
 #Grabs gradients of NN loss w.r.t params
 ForwardDiff.gradient(loss_PINNS, par)
end
#test grab_∇L
@time grab_∇L(p)

function dθ_dt!(du,u,k,t)
 #ODE used to evolve parmaeters based on the negative gradient definition (evolving in the direction of steepest descent)
 du .= -1 .* grab_∇L(u)
end
#test dθ_dt!
dθ_dt!(zeros(length(p)), p, 0.2, 0.2)

tspan_θ = (0.0,1.0)
prob_θ = ODEProblem(dθ_dt!,p*1e-4,tspan_θ)
integrator = DifferentialEquations.init(prob_θ,Tsit5())
loss_arr = []
params_arr = []

for (u,t) in tuples(integrator)
    loss = loss_PINNS(u)
    push!(loss_arr, loss)
    push!(params_arr, u)
    println(loss)
end

plot(log.(loss_arr))

#savefig("Loss_PINN-ROBER.png")

test_p = params_arr[end]
grab_pred(test_p)
plot(sol.t, grab_pred(test_p)')
plot!(sol.t, Array(sol)')


text:
using ANTsRegistration, Images, TestImages, Rotations, CoordinateTransformations, MappedArrays, FileIO, NRRD, Glob
using Unitful
using Test, Statistics

diff2_0(f, m) = m == 0 ? 0.0 : (Float64(f) - Float64(m))^2
abs2f(x) = abs2(Float64(x))

@testset "Single stage" begin
    img = testimage("cameraman")
    tfm = Translation(125,250) ∘ LinearMap(RotMatrix(pi/50)) ∘ Translation(-125,-250)
    img_rotated = warp(img, tfm)
    fixed = img[50:300, 50:300]
    moving = img_rotated[50:300, 50:300]
    mreg = register(fixed, moving, Stage(fixed, Global("Rigid")))/255
    @test sum(mappedarray(diff2_0, fixed, mreg)) < 1e-3*sum(abs2f, fixed)

    # For a bigger rotation we need to customize
    tfm = Translation(125,250) ∘ LinearMap(RotMatrix(pi/12)) ∘ Translation(-125,-250)
    img_rotated = warp(img, tfm)
    fixed = img[50:300, 50:300]
    moving = img_rotated[50:300, 50:300]
    mreg = register(fixed, moving, Stage(fixed, Global("Rigid"), MeanSquares(), (8,4,2,1), (10,7,3,0)))/255
    @test sum(mappedarray(diff2_0, fixed, mreg)) < 1e-3*sum(abs2f, fixed)
end

@testset "Multistage" begin
    img = testimage("cameraman")
    tfm = Translation(125,250) ∘ LinearMap(RotMatrix(pi/50)) ∘ Translation(-125,-250)
    img_rotated = warp(img, tfm)
    fixed = img[50:300, 50:300]
    moving = img_rotated[50:300, 50:300]
    rigid = Stage(fixed, Global("Rigid"))
    syn = Stage(fixed, SyN())
    mreg = register(fixed, moving, [rigid,syn])/255
    @test sum(mappedarray(diff2_0, fixed, mreg)) < 1e-3*sum(abs2f, fixed)
end

@testset "Options" begin
    # Histogram matching
    fixed = [0 0 0 0 0 0 0;
             0 0 0 0 0 0 0;
             0 0 10 0 0 0 0;
             0 0 0 0 0 0 0;
             0 0 0 0 1 1 0;
             0 0 0 0 0 0 0]
    moving = [11 0 0 0 0 0 0;
              0 0 0 0 0 0 0;
              0 0 8 8 0 0 0;
              zeros(3, 7)]
    pfixed  = padarray(fixed, Fill(0, (3,3))).parent
    pmoving = padarray(moving, Fill(0, (3,3))).parent
    stage = Stage(pfixed, Global("Translation"), MeanSquares(), (1,1), (3,0), (1000,1000))
    imgw = register(pfixed, pmoving, stage)
    @test sum(mappedarray(diff2_0, pfixed, adjust_histogram(Matching(), imgw, pfixed, 256))) > 1e-3*sum(abs2f, fixed)
    imgw = register(pfixed, pmoving, stage; histmatch=true)
    @test sum(mappedarray(diff2_0, pfixed, adjust_histogram(Matching(), imgw, pfixed, 0:12))) < 1e-3*sum(abs2f, fixed)

    # Winsorizing
    fixed = rand(100,100)
    moving = copy(fixed)
    fixed[50,50] = 1000
    moving[50,51] = 1000
    imgw = register(fixed, moving, stage)
    @test median(abs.(imgw[:,1:end-1] - moving[:,2:end])) < 1e-3
    imgw = register(fixed, moving, stage; winsorize=(0.05,0.95))
    @test median(abs.(imgw - moving)) < 0.01
end

@testset "TimeSeries" begin
    img = testimage("cameraman")
    fixed = img[50:300, 50:300]
    ## Create a movie with a time dimension and progressive rigid transformation
    testdir = ANTsRegistration.userpath()
    rawfile = joinpath(testdir, "series.raw")
    θs = ((0:5)/5)*(π/50)
    open(rawfile, "w") do io
        for θ in θs
            tfm = Translation(125,250) ∘ LinearMap(RotMatrix(θ)) ∘ Translation(-125,-250)
            img_rotated = warp(img, tfm)
            moving = img_rotated[50:300, 50:300]
            write(io, moving)
        end
    end
    hdrfile = joinpath(testdir, "series.nhdr")
    header = NRRD.headerinfo(N0f8, (Axis{:y}(axes(fixed, 1)), Axis{:x}(axes(fixed, 2)), Axis{:time}(axes(θs, 1))))
    header["datafile"] = rawfile
    # Hack to overcome ITK's limitations in reading NRRD files
    delete!(header, "space dimension")
    delete!(header, "space origin")
    open(hdrfile, "w") do io
        write(io, magic(format"NRRD"))
        NRRD.write_header(io, "0004", header)
    end

    ## Register the sequence
    output = (joinpath(testdir, "seriesinfo"), joinpath(testdir, "series_warped.nrrd"))
    motioncorr(output, fixed, hdrfile, Stage(fixed, Global("Rigid"), CC(), (2,1), (1,0), (15,3)))
    imgw = load(output[2])
    for i = 1:size(imgw, 3)
        @test sum(mappedarray(diff2_0, fixed, imgw[:,:,i]/255)) < 1e-3*sum(abs2f, fixed)
    end
    for tfmfile in glob("seriesinfo*.csv", testdir)
        rm(tfmfile)
    end
    rm(output[2])

    # Also do it via register. Exercise the low-level call that uses only filenames.
    moving = load(hdrfile)
    fixedname = joinpath(testdir, "series_fixed.png")
    save(fixedname, fixed)
    movingtmp = joinpath(testdir, "series_frame.png")
    stageMS = Stage(fixed, Global("Rigid"), MeanSquares())
    stageCC = Stage(fixed, Global("Rigid"), CC())
    for i = 2:size(moving, 3)
        save(movingtmp, moving[:,:,i])
        output = (joinpath(testdir, "info$i"), joinpath(testdir, "series_warp.nrrd"))
        register(output, sdims(fixed), fixedname, movingtmp, [stageMS])
        imgw = load(output[2])'/255  # FIXME NRRD orientation
        @test sum(mappedarray(diff2_0, fixed, imgw)) < 1e-3*sum(abs2f∘Float64, fixed)
        # register(output, sdims(fixed), fixedname, movingtmp, [stageCC])
        # imgw = load(output[2])/255
        # @test sum(mappedarray(diff2_0, fixed, imgw)) < 1e-3*sum(abs2f, fixed)
    end
    for tfmfile in glob("info*.mat", testdir)
        rm(tfmfile)
    end
    rm(output[2])
    rm(rawfile)
    rm(hdrfile)
    rm(fixedname)
    rm(movingtmp)
end

@testset "Anisotropic spacing" begin
    img = testimage("cameraman")
    fixed0 = img[50:300, 50:300]
    tfm = Translation(250,250) ∘ LinearMap(RotMatrix(-pi/4)) ∘ Translation(-250,-250)
    img_rotated = warp(img, tfm)
    moving0 = img_rotated[75:320, 75:325]
    ps = (1u"mm", 2u"mm")
    # If you don't trim the edges, the mixing with zeros "anchors" the image in place.
    # see resolution of https://github.com/ANTsX/ANTs/issues/675
    fixed = AxisArray(restrict(fixed0, 2)[:,2:end-1], (:y, :x), ps)
    moving = AxisArray(restrict(moving0, 2)[:,2:end-1], (:y, :x), ps)
    stage = Stage(fixed, Global("Rigid"), MeanSquares(), (1,), (25u"mm",), (1000,))
    # Both of the following should register well: fixed0 and moving0 are isotropic, so
    # the rotation is a rotation; fixed and moving encode their pixel spacing, so even
    # though "squashed" horizontally they also differ by a rotation (once the pixel spacing is accounted for)
    imgw0 = register(fixed0, moving0, stage)/255
    imgw = AxisArray(register(fixed, moving, stage), (:y, :x), ps)
    @test 2*mean(mappedarray(diff2_0, fixed0, imgw0)) >= mean(mappedarray(diff2_0, fixed, imgw))
    # In contrast, if we take away the pixel spacing info, we should not register well
    imgws = register(fixed.data, moving.data, stage)
    @test mean(mappedarray(diff2_0, fixed.data, imgws)) > 3*mean(mappedarray(diff2_0, fixed, imgw))
end


text:
"""
  TCG is a trust region sub-problem solver that approximately solves the quadratic
  programming problem subject to the solution being in a Euclidean Ball. It is
  appropriate for all types of Hessians. It solves the problem by using a conjugate
  gradient approach to solving the unconstrained quadratic optimization problem. If
  an interior solution exists it will be found according to the usual CG theory. If the
  problem is unbounded (negative curvature) or the solution is outside of the trust region
  we step to the boundary.

  TCG accepts all Hessians and is therefore well-suited for Newton's method for
  any vexity. It is relatively cheap for (very) large problems.

  The implementation follows Algorithm 7.5.1 on p. 205 of [ConnGouldTointBook] and
  the algorithm on p. 628 in [Steihaug1983]
"""
struct TCG <: TRSPSolver
end
summary(::TCG) = "Steihaug-Toint Truncated CG"

function (ms::TCG)(∇f, H, Δ::T, s, scheme, λ0=0; abstol=1e-10, maxiter=min(5,length(s)), κeasy=T(1)/10, κhard=T(2)/10) where T
    M = I
    # We only know that TCG has its properties if we start with s = 0
    s .= T(0)
    # g is the gradient of the quadratic model with the step α*p
    if all(iszero, ∇f)
        return tr_return(;λ=NaN, ∇f=∇f, H=H, s=s, interior=true, solved=true, hard_case=false, Δ=Δ)
    end
    g = copy(∇f)
    v = M\g
    p = -v
    for i = 1:maxiter
        # Check for negative curvature
        κ = dot(p, H*p)

        c = dot(s, M*s)
        b′ = dot(s, M*p) # b/2
        a = dot(p, M*p)

        if κ <= 0
          	# This branch just catches that sigma can be nan if a is exactly 0.
        	  if !iszero(a)
	              σ = (-b′ + sqrt(b′^2 + a*(Δ^2 - c)))/a
                @. s = s + σ*p
            end
            return tr_return(;λ=NaN, ∇f=∇f, H=H, s=s, interior=false, solved=true, hard_case=false, Δ=Δ)
        end

        α = dot(g, v)/κ

        if a*α^2 + α*2*b′ + c >= Δ^2
            σ = (-b′ + sqrt(b′^2 + a*(Δ^2 - c)))/a
            @. s = s + σ*p
            return tr_return(;λ=NaN, ∇f=∇f, H=H, s=s, interior=false, solved=true, hard_case=false, Δ=Δ)
        end
        @. s = s + α*p
        den = dot(g, v)
        g .= g .+ α.*(H*p)
        v .= M\g
        β = dot(g, v)/den
        @. p = -v + β*p
    end
    tr_return(;λ=NaN, ∇f=∇f, H=H, s=s, interior=true, solved=true, hard_case=false, Δ=Δ)
end


text:
module TestDecisionTree

using Test
using AMLPipelineBase
using Random
using DataFrames: nrow, DataFrame

function generateXY()
    Random.seed!(123)
    iris = getiris()
    indx = Random.shuffle(1:nrow(iris))
    features=iris[indx,1:4]
    sp = iris[indx,5] |> Vector
    (features,sp)
end

const X,Y = generateXY()

function test_decisiontree()
    learners = Dict(:rf=>RandomForest(),:ada=>Adaboost(),:ptree=>PrunedTree())
    results = Dict(:rf=>50.0,:ada=>50.0,:ptree=>50.0)
    for (name,obj) in learners
        fit!(obj,X,Y)
        res = transform!(obj,X)
        @testset "$name: Full dataset mutate" begin
            @test sum(res .== Y)/length(Y)*100 |> floor > results[name]
        end
    end
    for (name,obj) in learners
        ob=fit(obj,X,Y)
        res = transform(ob,X)
        @testset "$name: Full dataset" begin
            @test sum(res .== Y)/length(Y)*100 |> floor > results[name]
        end
    end
    trndx = 1:80
    tstndx = 81:nrow(X)
    results = Dict(:rf=>50.0,:ada=>50.0,:ptree=>50.0)
    for (name,obj) in learners
        fit!(obj,X[trndx,:],Y[trndx])
        res = transform!(obj,X[tstndx,:])
        @testset "$name: partial dataset mutate" begin
            @test sum(res .== Y[tstndx])/length(Y[tstndx])*100 |> floor > results[name]
        end
    end
    for (name,obj) in learners
        ob=fit(obj,X[trndx,:],Y[trndx])
        res = transform(ob,X[tstndx,:])
        @testset "$name: partial dataset" begin
            @test sum(res .== Y[tstndx])/length(Y[tstndx])*100 |> floor > results[name]
        end
    end
end
@testset "DecisionTrees" begin
    Random.seed!(123)
    test_decisiontree()
end

end


text:
@testset "Robots: quadruped" begin
    q0 = [0.0; 0.0; 0.0; 0.5 * π * ones(8)]
    v0 = zeros(11)
    @test norm(lagrangian(quadruped, q0, v0)) < 1.0e-8

    # visualize
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, quadruped, [q0], Δt=0.1);
end


text:
@testset "Testing Slicing Submodule" begin

I = ones(1,3,3,3,3)

IA = makeAxisArray(I, [1.0,1.0,1.0]./256, -[0.5,0.5,0.0].+[1.0,1.0,0.0]./256, 1.0 )

@test sliceTimeDim(IA, 3) == ones(1,3,3,3)
@test sliceTimeDim(IA, "MIP") == ones(1,3,3,3)

IAL = sliceTimeDim(IA, 3)

x,y,z = threeSlices(IAL,"MIP")
@test x == ones(1,3,3)
@test y == ones(1,3,3)
@test z == ones(1,3,3)
x,y,z = threeSlices(IAL,[1,1,1])
@test x == ones(1,3,3)
@test y == ones(1,3,3)
@test z == ones(1,3,3)

I = ones(3,3,3)

IAL = makeAxisArray(I, [1.0,1.0,1.0]./256, -[0.5,0.5,0.0].+[1.0,1.0,0.0]./256 )

x,y,z = threeSlices(IAL,"MIP")
@test x == ones(3,3)
@test y == ones(3,3)
@test z == ones(3,3)
x,y,z = threeSlices(IAL,[1,1,1])
@test x == ones(3,3)
@test y == ones(3,3)
@test z == ones(3,3)

I = ones(1,3,3,3,3)

IA = makeAxisArray(I, [1.0,1.0,1.0]./256, -[0.5,0.5,0.0].+[1.0,1.0,0.0]./256, 1.0 )
IAL = sliceTimeDim(IA, 3)

params = Dict{Symbol,Any}()
params[:sliceX] = 1
params[:sliceY] = 1
params[:sliceZ] = 1
params[:spatialMIP] = true
params[:blendChannels] = true
params[:complexBlending] = false
params[:activeChannel] = 1

xx,yy,zz = getColoredSlices(IAL, nothing, [ColoringParams(0,1,"gray")], [0.0], [1.0], params)
@show typeof(xx)

xx,yy,zz = getColoredSlicesMovie(IA, nothing, [ColoringParams(0,1,"gray")], params)
@show typeof(xx)

end


text:
import Base: +, -, *, /

# TODO: more operators, like AbstractArray

# Unary ops
@inline -(a::StaticArray) = map(-, a)

# Binary ops
# Between arrays
@inline +(a::StaticArray, b::StaticArray) = map(+, a, b)
@inline +(a::AbstractArray, b::StaticArray) = map(+, a, b)
@inline +(a::StaticArray, b::AbstractArray) = map(+, a, b)

@inline -(a::StaticArray, b::StaticArray) = map(-, a, b)
@inline -(a::AbstractArray, b::StaticArray) = map(-, a, b)
@inline -(a::StaticArray, b::AbstractArray) = map(-, a, b)

# Scalar-array
@inline +(a::Number, b::StaticArray) = broadcast(+, a, b)
@inline +(a::StaticArray, b::Number) = broadcast(+, a, b)

@inline -(a::Number, b::StaticArray) = broadcast(-, a, b)
@inline -(a::StaticArray, b::Number) = broadcast(-, a, b)

@inline *(a::Number, b::StaticArray) = broadcast(*, a, b)
@inline *(a::StaticArray, b::Number) = broadcast(*, a, b)

@inline /(a::StaticArray, b::Number) = broadcast(/, a, b)
@inline \(a::Number, b::StaticArray) = broadcast(\, a, b)

# With UniformScaling
@generated function +(a::StaticMatrix, b::UniformScaling)
    n = size(a,1)
    newtype = similar_type(a, promote_type(eltype(a), eltype(b)))

    if size(a,2) != n
        error("Dimension mismatch")
    end

    exprs = [i == j ? :(a[$(sub2ind(size(a), i, j))] + b.λ) : :(a[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function +(a::UniformScaling, b::StaticMatrix)
    n = size(b,1)
    newtype = similar_type(b, promote_type(eltype(a), eltype(b)))

    if size(b,2) != n
        error("Dimension mismatch")
    end

    exprs = [i == j ? :(a.λ + b[$(sub2ind(size(a), i, j))]) : :(b[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function -(a::StaticMatrix, b::UniformScaling)
    n = size(a,1)
    newtype = similar_type(a, promote_type(eltype(a), eltype(b)))

    if size(a,2) != n
        error("Dimension mismatch")
    end

    exprs = [i == j ? :(a[$(sub2ind(size(a), i, j))] - b.λ) : :(a[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function -(a::UniformScaling, b::StaticMatrix)
    n = size(b,1)
    newtype = similar_type(b, promote_type(eltype(a), eltype(b)))

    if size(b,2) != n
        error("Dimension mismatch")
    end

    exprs = [i == j ? :(a.λ - b[$(sub2ind(size(a), i, j))]) : :(-b[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@inline *(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ * b
@inline *(a::StaticMatrix, b::UniformScaling) = a * b.λ
@inline \(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ \ b
@inline /(a::StaticMatrix, b::UniformScaling) = a / b.λ

# Transpose, conjugate, etc
@inline conj(a::StaticArray) = map(conj, a)

@generated function transpose(v::StaticVector)
    n = length(v)
    newtype = similar_type(v, (1,n))
    exprs = [:(v[$j]) for j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function ctranspose(v::StaticVector)
    n = length(v)
    newtype = similar_type(v, (1,n))
    exprs = [:(conj(v[$j])) for j = 1:n]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function transpose(m::StaticMatrix)
    (s1,s2) = size(m)
    if s1 == s2
        newtype = m
    else
        newtype = similar_type(m, (s2,s1))
    end

    exprs = [:(m[$j1, $j2]) for j2 = 1:s2, j1 = 1:s1]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function ctranspose(m::StaticMatrix)
    (s1,s2) = size(m)
    if s1 == s2
        newtype = m
    else
        newtype = similar_type(m, (s2,s1))
    end

    exprs = [:(conj(m[$j1, $j2])) for j2 = 1:s2, j1 = 1:s1]

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@inline vcat(a::Union{StaticVector,StaticMatrix}) = a
@generated function vcat(a::Union{StaticVector, StaticMatrix}, b::Union{StaticVector,StaticMatrix})
    if size(a,2) != size(b,2)
        error("Dimension mismatch")
    end

    if a <: StaticVector && b <: StaticVector
        newtype = similar_type(a, (length(a) + length(b),))
        exprs = vcat([:(a[$i]) for i = 1:length(a)],
                     [:(b[$i]) for i = 1:length(b)])
    else
        newtype = similar_type(a, (size(a,1) + size(b,1), size(a,2)))
        exprs = [((i <= size(a,1)) ? ((a <: StaticVector) ? :(a[$i]) : :(a[$i,$j]))
                                   : ((b <: StaticVector) ? :(b[$(i-size(a,1))]) : :(b[$(i-size(a,1)),$j])))
                                   for i = 1:(size(a,1)+size(b,1)), j = 1:size(a,2)]
    end

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end
# TODO make these more efficient
@inline vcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}, c::Union{StaticVector,StaticMatrix}) =
    vcat(vcat(a,b), c)
@inline vcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}, c::Union{StaticVector,StaticMatrix}...) =
    vcat(vcat(a,b), c...)

@generated function hcat(a::StaticVector)
    newtype = similar_type(a, (length(a),1))
    exprs = [:(a[$i]) for i = 1:length(a)]
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end
@inline hcat(a::StaticMatrix) = a
@generated function hcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix})
    if size(a,1) != size(b,1)
        error("Dimension mismatch")
    end

    exprs1 = [:(a[$i]) for i = 1:length(a)]
    exprs2 = [:(b[$i]) for i = 1:length(b)]

    newtype = similar_type(a, (size(a,1), size(a,2) + size(b,2)))

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs1..., exprs2...)))
    end
end
# TODO make these more efficient
@inline hcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}, c::Union{StaticVector,StaticMatrix}) =
    hcat(hcat(a,b), c)
@inline hcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}, c::Union{StaticVector,StaticMatrix}...) =
    hcat(hcat(a,b), c...)

@generated function eye{SA <: StaticArray}(::Type{SA})
    s = size(SA)
    T = eltype(SA)
    if T == Any
        T = Float64
    end
    if length(s) != 2
        error("Must call `eye` with a two-dimensional array. Got size $s")
    end
    e = eye(T, s...)
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, SA, Expr(:tuple, e...)))
    end
end

@generated function eye{SA <: StaticMatrix}(::SA)
    s = size(SA)
    T = eltype(SA)
    e = eye(T, s...)
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, SA, Expr(:tuple, e...)))
    end
end

@generated function cross(a::StaticVector, b::StaticVector)
    if length(a) === 3 && length(b) === 3
        return quote
            $(Expr(:meta, :inline))
            similar_type(a, promote_type(eltype(a), eltype(b)))((a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]))
        end
    else
        error("Cross product only defined for 3-vectors")
    end
end

@generated function dot(a::StaticVector, b::StaticVector)
    if length(a) == length(b)
        expr = :(a[1] * b[1])
        for j = 2:length(a)
            expr = :($expr + a[$j] * b[$j])
        end

        return quote
            $(Expr(:meta, :inline))
            @inbounds return $expr
        end
    else
        error("dot() expects vectors of same length. Got lengths $(length(a)) and $(length(b)).")
    end
end

@generated function vecdot(a::StaticArray, b::StaticArray)
    if length(a) == length(b)
        expr = :(a[1] * b[1])
        for j = 2:length(a)
            expr = :($expr + a[$j] * b[$j])
        end

        return quote
            $(Expr(:meta, :inline))
            @inbounds return $expr
        end
    else
        error("vecdot() expects arrays of same length. Got sizes $(size(a)) and $(size(b)).")
    end
end

@generated function vecnorm(a::StaticArray)
    if length(a) == 0
        return zero(eltype(a))
    end

    expr = :(real(conj(a[1]) * a[1]))
    for j = 2:length(a)
        expr = :($expr + real(conj(a[$j]) * a[$j]))
    end

    return quote
        $(Expr(:meta, :inline))
        @inbounds return sqrt($expr)
    end
end

@generated function trace(a::StaticMatrix)
    s = size(a)
    if s[1] != s[2]
        error("matrix is not square")
    end

    if s[1] == 0
        return zero(eltype(a))
    end

    exprs = [:(a[$(sub2ind(s, i, i))]) for i = 1:s[1]]
    total = reduce((ex1, ex2) -> :($ex1 + $ex2), exprs)

    return quote
        $(Expr(:meta, :inline))
        @inbounds return $total
    end
end

#=

abstract SVector{Size,T,D} <: StaticArray{T,D}
@pure size{Size}(::SVector{Size}) = (Size,)
@pure size{T<:SVector}(::Type{T}) = (T.parameters[1],)

#typealias SVector{Size,T,D} SArray{Size,T,1,D}
typealias MVector{Size,T,D} MArray{Size,T,1,D}
typealias StaticVector{Size,T,D} Union{SVector{Size,T,D},MVector{Size,T,D}}

typealias SMatrix{Sizes,T,D} SArray{Sizes,T,2,D}
typealias MMatrix{Sizes,T,D} MArray{Sizes,T,2,D}
typealias StaticMatrix{Size,T,D} Union{SMatrix{Size,T,D},MMatrix{Size,T,D}}

# Constructors not covered by SArray
@generated Base.call{N,T}(::Type{SVector}, d::NTuple{N,T}) = :($(SVector{(N,),T,d})(d))
@generated function Base.call{N}(::Type{SVector}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(SVector{(N,),T})(convert_ntuple($T, d)))
end
@generated Base.call{Size,N,T}(::Type{SVector{Size}}, d::NTuple{N,T}) = :($(SVector{Size,T,d})(d))
@generated function Base.call{Size,N}(::Type{SVector{Size}}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(SVector{Size,T})(convert_ntuple($T, d)))
end

@generated Base.call{Sizes,N,T}(::Type{SMatrix{Sizes}}, d::NTuple{N,T}) = :($(SMatrix{Sizes,T,d})(d))
@generated function Base.call{Sizes,N}(::Type{SMatrix{Sizes}}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(SMatrix{Sizes,T})(convert_ntuple($T, d)))
end

# Conversions to SVector from AbstractVector
@generated function Base.convert{Sizes,T}(::Type{SVector{Sizes}}, a::AbstractVector{T})
    M = prod(Sizes)
    NewType = SArray{Sizes,T,1,NTuple{M,T}}
    quote
        size(a) == Sizes || error("Input array must have size $Sizes, got $(size(a))")
        $(NewType)((a...))
    end
end

# Conversions to SMatrix from AbstractMatrix
@generated function Base.convert{Sizes,T}(::Type{SMatrix{Sizes}}, a::AbstractMatrix{T})
    M = prod(Sizes)
    NewType = SArray{Sizes,T,2,NTuple{M,T}}
    quote
        size(a) == Sizes || error("Input array must have size $Sizes, got $(size(a))")
        $(NewType)((a...))
    end
end

# Constructors not covered by MArray
@generated Base.call{N,T}(::Type{MVector}, d::NTuple{N,T}) = :($(MVector{(N,),T,d})(d))
@generated function Base.call{N}(::Type{MVector}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(MVector{(N,),T})(convert_ntuple($T, d)))
end
@generated Base.call{Size,N,T}(::Type{MVector{Size}}, d::NTuple{N,T}) = :($(MVector{Size,T,d})(d))
@generated function Base.call{Size,N}(::Type{MVector{Size}}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(MVector{Size,T})(convert_ntuple($T, d)))
end

@generated Base.call{Sizes,N,T}(::Type{MMatrix{Sizes}}, d::NTuple{N,T}) = :($(MMatrix{Sizes,T,d})(d))
@generated function Base.call{Sizes,N}(::Type{MMatrix{Sizes}}, d::NTuple{N})
    T = promote_tuple_eltype(d)
    return :($(MMatrix{Sizes,T})(convert_ntuple($T, d)))
end
# Conversions to MVector from AbstractVector
@generated function Base.convert{Sizes,T}(::Type{MVector{Sizes}}, a::AbstractVector{T})
    M = prod(Sizes)
    NewType = MArray{Sizes,T,1,NTuple{M,T}}
    quote
        size(a) == Sizes || error("Input array must have size $Sizes, got $(size(a))")
        $(NewType)((a...))
    end
end

# Conversions to MMatrix from AbstractMatrix
@generated function Base.convert{Sizes,T}(::Type{MMatrix{Sizes}}, a::AbstractMatrix{T})
    M = prod(Sizes)
    NewType = MArray{Sizes,T,2,NTuple{M,T}}
    quote
        size(a) == Sizes || error("Input array must have size $Sizes, got $(size(a))")
        $(NewType)((a...))
    end
end

# Conversions to Vector from StaticVector
function Base.convert{Sizes,T}(::Type{Vector},a::StaticVector{Sizes,T})
    out = Vector{T}(Sizes...)
    out[:] = a[:]
    return out
end

# Conversions to Matrix from StaticMatrix
function Base.convert{Sizes,T}(::Type{Matrix},a::StaticMatrix{Sizes,T})
    out = Matrix{T}(Sizes...)
    out[:] = a[:]
    return out
end

# Level-0 ops: Reductions *should* already work

# Level-1 ops: Mostly covered by map(), etc. Vector transpose... need operators for Julia 0.4.

# Level-2 ops: matrix-vector multiplication and matrix transpose

# Level-3 ops: matrix-matrix multiplication
=#


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    SmoothingProcess(σ=3)

A smoothing process (i.e. Gaussian filter) with `σ` parameter.
"""
struct SmoothingProcess <: TimelessProcess
  σ::Float64
end

SmoothingProcess() = SmoothingProcess(3.0)

evolve!(land::Matrix, proc::SmoothingProcess) =
  imfilter!(land, land, centered(Kernel.gaussian(proc.σ)), "replicate")


text:
include("FpOutput.jl")

@recipe function f(fpds::FpDataset)
    if fpds.dataset isa Matrix
        attr = attrib(fpds.fpoutput, fpds.varname)
        lons, lats = Flexpart.mesh(fpds.fpoutput)
        sel = selected(fpds)
        title --> "$(attr["name"]) - $(haskey(attr, "long_name") ? attr["long_name"] : nothing) [$(attr["units"])] time = $(sel[:time]) height = $(sel[:height])"
        (lons, lats, fpds.dataset' .|> log10)
    else
        error("The dataset must be 2D")
    end
end
# function plot_filtered_grid(file, timestep)
#     rellon, rellat = relloc(file)
#     lon, lat, conc = filtered_fields(file, timestep)

#     mg_lon = lon .* ones(length(lat))'
#     mg_lat = lat' .* ones(length(lon))

#     scatter(rellon, rellat, markershape=:circle)
#     scatter!(mg_lon, mg_lat, markershape=:x, legend=false)
# end

# function with_pyplot(f::Function)
#     f()
#     fig = gcf()
#     close(fig)
#     return fig
# end

# function plot_contour(file, timestep)
#     lon, lat, conc = fields(file, timestep)

#     lonf, latf, concf = filtered_fields(file, timestep)
#     z0max = maximum(concf)
#     z0min = minimum(concf)

#     with_pyplot() do
#         # p = PyPlot.pcolormesh(lon, lat, conc')
#         p = PyPlot.pcolormesh(lon, lat, conc', shading="nearest", norm=matplotlib[:colors][:LogNorm](vmin=z0min, vmax=z0max), cmap=:jet, vmax=z0max)
#         colorbar(p)
#         gca().set_aspect(1)
#     end
# end

# function heat(file, timestep)
#     lon, lat, conc = filtered_fields(file, timestep)
#     heatmap(lon, lat, conc')
# end


text:
function x2c!(c::Array{Complex{Float64},1},ψ::Array{Complex{Float64},1},tinfo::ProjectedGPE.Tinfo)
    c = tinfo.Tx'*(tinfo.W.*ψ)
end

function x2c!(c::Array{Complex{Float64},2},ψ::Array{Complex{Float64},2},tinfo::ProjectedGPE.Tinfo)
    c .= tinfo.Tx'*(tinfo.W.*ψ)*tinfo.Ty
end

function x2c!(c::Array{Complex{Float64},3},ψ::Array{Complex{Float64},3},tinfo::ProjectedGPE.Tinfo)
    sa = size(ψ);sx = size(tinfo.Tx');sy = size(tinfo.Ty');sz = size(tinfo.Tz')
    Ax = reshape(tinfo.W.*ψ,(sa[1],sa[2]*sa[3]))
    Ax = reshape(tinfo.Tx'*Ax,(sx[1],sa[2],sa[3]))
    Ax = permutedims(Ax,(2,3,1));sa = size(Ax)

    Ay = reshape(Ax,(sa[1],sa[2]*sa[3]))
    Ay = reshape(tinfo.Ty'*Ay,(sy[1],sa[2],sa[3]))
    Ay = permutedims(Ay,(2,3,1));sa = size(Ay)

    Az = reshape(Ay,(sa[1],sa[2]*sa[3]))
    Az = reshape(tinfo.Tz'*Az,(sz[1],sa[2],sa[3]))
    ψ = permutedims(Az,(2,3,1))

end


text:
export
MarginalRule,
@marginalRule

"""
`MarginalRule{factor_type}` specifies a joint marginal update rule with respect
to a node of type `factor_type`.
"""
abstract type MarginalRule{factor_type} <: MarginalUpdateRule end

"""
Construct a `MarginalEntry` for computing the marginal over `cluster`
through a node-specific joint marginal update rule.
"""
function MarginalEntry(target::Cluster, outbound_types::Dict{Interface, Type})
    inbound_types = collectInboundTypes(target, outbound_types)
    marginal_update_rule = inferMarginalRule(target, inbound_types)

    # Collect inbound interfaces
    inbound_interfaces = Interface[]
    for edge in target.edges
        if edge.a in target.node.interfaces
            push!(inbound_interfaces, edge.a.partner) # Partner is the required inbound interface
        else
            push!(inbound_interfaces, edge.b.partner)
        end
    end

    return MarginalEntry(target, inbound_interfaces, marginal_update_rule)
end

"""
Infer the rule that computes the joint marginal over `cluster`
"""
function inferMarginalRule(cluster::Cluster, inbound_types::Vector{<:Type})
    # Find applicable rule(s)
    applicable_rules = Type[]
    for rule in leaftypes(MarginalRule{typeof(cluster.node)})
        if isApplicable(rule, inbound_types)
            push!(applicable_rules, rule)
        end
    end

    # Select and set applicable rule
    if isempty(applicable_rules)
        error("No applicable marginal update rule for $(typeof(cluster.node)) node with inbound types: $(join(inbound_types, ", "))")
    elseif length(applicable_rules) > 1
        error("Multiple applicable marginal update rules for $(typeof(cluster.node)) node with inbound types: $(join(inbound_types, ", ")): $(join(applicable_rules, ", "))")
    else
        marginal_update_rule = first(applicable_rules)
    end

    return marginal_update_rule
end

"""
Construct the marginal computations table for a given posterior factor.
The marginal table defines which marginal posteriors are computed.
"""
function marginalTable(pf::PosteriorFactor)
    # Construct outbound types dictionary
    outbound_types = Dict{Interface, Type}()
    for entry in pf.schedule
        outbound_types[entry.interface] = outboundType(entry.message_update_rule)
    end

    variable_table = marginalTable(sort(collect(pf.target_variables)))
    cluster_table = [MarginalEntry(cluster, outbound_types) for cluster in sort(collect(pf.target_clusters))]

    return [variable_table; cluster_table]
end

"""
Find the inbound types that are required to compute a joint marginal over `target`.
Returns a vector with inbound types that correspond with required interfaces.
"""
function collectInboundTypes(cluster::Cluster, outbound_types::Dict{Interface, Type})
    inbound_types = Type[]
    cluster_pf = posteriorFactor(first(cluster.edges))
    encountered_external_regions = Set{Region}()
    for node_interface in cluster.node.interfaces
        current_region = region(cluster.node, node_interface.edge) # Returns a Variable if no cluster is assigned
        current_pf = posteriorFactor(node_interface.edge) # Returns an Edge if no posterior factor is assigned
        inbound_interface = ultimatePartner(node_interface)

        if (current_pf === cluster_pf) && (current_region === cluster)
            # Edge is internal and in cluster, accept message
            push!(inbound_types, outbound_types[inbound_interface])
        elseif (current_pf === cluster_pf)
            # Edge is internal but not in cluster, signal to rule signature that edge will be marginalized out by appending "Nothing"
            push!(inbound_types, Nothing)
        elseif !(current_region in encountered_external_regions)
            # Edge is external and cluster is not yet encountered, accept probability distribution
            push!(inbound_types, ProbabilityDistribution)
            push!(encountered_external_regions, current_region) # Register current region with encountered external regions
        end
    end

    return inbound_types
end

"""
`@marginalRule` registers a marginal update rule for a (joint) marginal
by defining the rule type and the corresponding methods for the `isApplicable` function.
If no name (type) for the new rule is passed, a unique name (type) will be generated.
Returns the rule type.
"""
macro marginalRule(fields...)
    # Init required fields in macro scope
    node_type = :unknown
    inbound_types = :unknown
    name = :auto # Triggers automatic naming unless overwritten

    # Loop over fields because order is unknown
    for arg in fields
        (arg.args[1] == :(=>)) || error("Invalid call to @marginalRule")

        if arg.args[2].value == :node_type
            node_type = arg.args[3]
        elseif arg.args[2].value == :inbound_types
            inbound_types = arg.args[3]
            (inbound_types.head == :tuple) || error("Inbound types should be passed as Tuple")
        elseif arg.args[2].value == :name
            name = arg.args[3]
        else
            error("Unrecognized field $(arg.args[2].value) in call to @marginalRule")
        end
    end

    # Assign unique name if not set already
    if name == :auto
        # Added hash ensures that the rule name is unique
        msg_types_hash = string(hash(vcat([outbound_type], inbound_types)))[1:6]
        name = Symbol("Marginal$(node_type)$(msg_types_hash)")
    end

    # Build validators for isApplicable
    input_type_validators = Expr[]

    push!(input_type_validators, :(length(input_types) == $(length(inbound_types.args))))
    for (i, i_type) in enumerate(inbound_types.args)
        if i_type != :Nothing
            push!(input_type_validators, :(ForneyLab.matches(input_types[$i], $i_type)))
        else # For a marginal rule, "Nothing" means that the input is marginalized out
            push!(input_type_validators, :(input_types[$i] == Nothing))
        end
    end

    expr = quote
        struct $name <: MarginalRule{$node_type} end
        ForneyLab.isApplicable(::Type{$name}, input_types::Vector{<:Type}) = begin
            $(reduce((current, item) -> :($current && $item), input_type_validators, init = :true))
        end
    end

    return esc(expr)
end

"""
Construct the inbound code that computes the marginal for `entry`. Allows for
overloading and for a user the define custom node-specific inbounds collection.
Returns a vector with inbounds that correspond with required interfaces.
"""
collectInbounds(entry::MarginalEntry) = collectMarginalNodeInbounds(entry.target.node, entry)

function collectMarginalNodeInbounds(::FactorNode, entry::MarginalEntry)
    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry
    target_to_marginal_entry = current_inference_algorithm.target_to_marginal_entry
    inbound_cluster = entry.target # Entry target is a cluster

    inbounds = Any[]
    entry_pf = posteriorFactor(first(entry.target.edges))
    encountered_external_regions = Region[]
    for node_interface in entry.target.node.interfaces
        current_region = region(inbound_cluster.node, node_interface.edge) # Note: edges that are not assigned to a posterior factor are assumed mean-field
        current_pf = posteriorFactor(node_interface.edge) # Returns an Edge if no posterior factor is assigned
        inbound_interface = ultimatePartner(node_interface)

        if isClamped(inbound_interface) && (current_pf === entry_pf)
            # Edge is clamped and internal, hard-code message of clamp node
            push!(inbounds, assembleClamp!(copy(inbound_interface.node), Message)) # Copy Clamp before assembly to prevent overwriting dist_or_msg field
        elseif isClamped(inbound_interface)
            # Edge is clamped and external, hard-code marginal of clamp node
            push!(inbounds, assembleClamp!(copy(inbound_interface.node), ProbabilityDistribution)) # Copy Clamp before assembly to prevent overwriting dist_or_msg field
        elseif (current_pf === entry_pf)
            # Edge is internal, collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        elseif !(current_region in encountered_external_regions)
            # Edge is external and region is not yet encountered, collect marginal from marginal dictionary
            push!(inbounds, target_to_marginal_entry[current_region])
            push!(encountered_external_regions, current_region) # Register current region with encountered external regions
        end
    end

    return inbounds
end


text:
using Test, Example

@test hello("Julia") == "Hello, Julia"
@test domath(2.0) ≈ 7.0

# julia-actions/julia-runtest
touch(joinpath(ENV["HOME"], "julia-runtest"))


text:
module TestUtils
    using Base.Test
    using DataFrames

    import DataFrames: identifier

    @test identifier("%_B*_\tC*") == :_B_C_
    @test identifier("2a") == :x2a
    @test identifier("!") == :x!
    @test identifier("\t_*") == :_
    @test identifier("begin") == :_begin
    @test identifier("end") == :_end

    @test DataFrames.make_unique([:x, :x, :x_1, :x2]) == [:x, :x_2, :x_1, :x2]

    # Check that reserved words are up to date
    f = "$JULIA_HOME/../../src/julia-parser.scm"
    if isfile(f)
        r1 = r"define reserved-words '\(([^)]+)"
        r2 = r"define \(parse-block s(?: \([^)]+\))?\)\s+\(parse-Nary s (?:parse-eq '\([^(]*|down '\([^)]+\) '[^']+ ')\(([^)]+)"
        body = readall(f)
        m1, m2 = match(r1, body), match(r2, body)
        if m1 == nothing || m2 == nothing
            error("Unable to extract keywords from 'julia-parser.scm'.")
        else
            rw = Set(split(m1.captures[1]*" "*m2.captures[1], r"\W+"))
            @test rw == DataFrames.RESERVED_WORDS
        end
    else
        warn("Unable to validate reserved words against parser. ",
             "Expected if Julia was not built from source.")
    end

    @test DataFrames.countna([1:3;]) == 0

    data = @data rand(20)
    @test DataFrames.countna(data) == 0
    data[sample(1:20, 11, replace=false)] = NA
    @test DataFrames.countna(data) == 11
    data[1:end] = NA
    @test DataFrames.countna(data) == 20

    pdata = @data sample(1:5, 20)
    @test DataFrames.countna(pdata) == 0
    pdata[sample(1:20, 11, replace=false)] = NA
    @test DataFrames.countna(pdata) == 11
    pdata[1:end] = NA
    @test DataFrames.countna(pdata) == 20

    funs = [mean, sum, var, x -> sum(x)]
    @test DataFrames._fnames(funs) == ["mean", "sum", "var", "λ1"]
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Host <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Host(; kwargs...)
        obj = new(meta(Host), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Host
const __meta_Host = Ref{ProtoMeta}()
function meta(::Type{Host})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Host)
            __meta_Host[] = target = ProtoMeta(Host)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :port => Int64, :load => Int64, :entities => Base.Vector{AbstractString}, :id => AbstractString]
            meta(target, Host, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Host[]
    end
end
function Base.getproperty(obj::Host, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :port
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :load
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :entities
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :id
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct PlacementTable_LoadMapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementTable_LoadMapEntry(; kwargs...)
        obj = new(meta(PlacementTable_LoadMapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementTable_LoadMapEntry (mapentry)
const __meta_PlacementTable_LoadMapEntry = Ref{ProtoMeta}()
function meta(::Type{PlacementTable_LoadMapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementTable_LoadMapEntry)
            __meta_PlacementTable_LoadMapEntry[] = target = ProtoMeta(PlacementTable_LoadMapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Host]
            meta(target, PlacementTable_LoadMapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementTable_LoadMapEntry[]
    end
end
function Base.getproperty(obj::PlacementTable_LoadMapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Host
    else
        getfield(obj, name)
    end
end

mutable struct PlacementTable_HostsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementTable_HostsEntry(; kwargs...)
        obj = new(meta(PlacementTable_HostsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementTable_HostsEntry (mapentry)
const __meta_PlacementTable_HostsEntry = Ref{ProtoMeta}()
function meta(::Type{PlacementTable_HostsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementTable_HostsEntry)
            __meta_PlacementTable_HostsEntry[] = target = ProtoMeta(PlacementTable_HostsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => UInt64, :value => AbstractString]
            meta(target, PlacementTable_HostsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementTable_HostsEntry[]
    end
end
function Base.getproperty(obj::PlacementTable_HostsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct PlacementTable <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementTable(; kwargs...)
        obj = new(meta(PlacementTable), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementTable
const __meta_PlacementTable = Ref{ProtoMeta}()
function meta(::Type{PlacementTable})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementTable)
            __meta_PlacementTable[] = target = ProtoMeta(PlacementTable)
            pack = Symbol[:sorted_set]
            allflds = Pair{Symbol,Union{Type,String}}[:hosts => Base.Dict{UInt64,AbstractString}, :sorted_set => Base.Vector{UInt64}, :load_map => Base.Dict{AbstractString,Host}, :total_load => Int64]
            meta(target, PlacementTable, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementTable[]
    end
end
function Base.getproperty(obj::PlacementTable, name::Symbol)
    if name === :hosts
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{UInt64,AbstractString}
    elseif name === :sorted_set
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{UInt64}
    elseif name === :load_map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Host}
    elseif name === :total_load
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct PlacementTables_EntriesEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementTables_EntriesEntry(; kwargs...)
        obj = new(meta(PlacementTables_EntriesEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementTables_EntriesEntry (mapentry)
const __meta_PlacementTables_EntriesEntry = Ref{ProtoMeta}()
function meta(::Type{PlacementTables_EntriesEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementTables_EntriesEntry)
            __meta_PlacementTables_EntriesEntry[] = target = ProtoMeta(PlacementTables_EntriesEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => PlacementTable]
            meta(target, PlacementTables_EntriesEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementTables_EntriesEntry[]
    end
end
function Base.getproperty(obj::PlacementTables_EntriesEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::PlacementTable
    else
        getfield(obj, name)
    end
end

mutable struct PlacementTables <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementTables(; kwargs...)
        obj = new(meta(PlacementTables), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementTables
const __meta_PlacementTables = Ref{ProtoMeta}()
function meta(::Type{PlacementTables})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementTables)
            __meta_PlacementTables[] = target = ProtoMeta(PlacementTables)
            allflds = Pair{Symbol,Union{Type,String}}[:entries => Base.Dict{AbstractString,PlacementTable}, :version => AbstractString]
            meta(target, PlacementTables, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementTables[]
    end
end
function Base.getproperty(obj::PlacementTables, name::Symbol)
    if name === :entries
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,PlacementTable}
    elseif name === :version
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct PlacementOrder <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function PlacementOrder(; kwargs...)
        obj = new(meta(PlacementOrder), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct PlacementOrder
const __meta_PlacementOrder = Ref{ProtoMeta}()
function meta(::Type{PlacementOrder})
    ProtoBuf.metalock() do
        if !isassigned(__meta_PlacementOrder)
            __meta_PlacementOrder[] = target = ProtoMeta(PlacementOrder)
            allflds = Pair{Symbol,Union{Type,String}}[:tables => PlacementTables, :operation => AbstractString]
            meta(target, PlacementOrder, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_PlacementOrder[]
    end
end
function Base.getproperty(obj::PlacementOrder, name::Symbol)
    if name === :tables
        return (obj.__protobuf_jl_internal_values[name])::PlacementTables
    elseif name === :operation
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

# service methods for Placement
const _Placement_methods = MethodDescriptor[
        MethodDescriptor("ReportDaprStatus", 1, Channel{Host}, Channel{PlacementOrder})
    ] # const _Placement_methods
const _Placement_desc = ServiceDescriptor("dapr.proto.placement.v1.Placement", 1, _Placement_methods)

Placement(impl::Module) = ProtoService(_Placement_desc, impl)

mutable struct PlacementStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    PlacementStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_Placement_desc, channel))
end # mutable struct PlacementStub

mutable struct PlacementBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    PlacementBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_Placement_desc, channel))
end # mutable struct PlacementBlockingStub

ReportDaprStatus(stub::PlacementStub, controller::ProtoRpcController, inp::Channel{Host}, done::Function) = call_method(stub.impl, _Placement_methods[1], controller, inp, done)
ReportDaprStatus(stub::PlacementBlockingStub, controller::ProtoRpcController, inp::Channel{Host}) = call_method(stub.impl, _Placement_methods[1], controller, inp)

export PlacementOrder, PlacementTables_EntriesEntry, PlacementTables, PlacementTable_HostsEntry, PlacementTable_LoadMapEntry, PlacementTable, Host, Placement, PlacementStub, PlacementBlockingStub, ReportDaprStatus
# mapentries: "PlacementTable_HostsEntry" => ("UInt64", "AbstractString"), "PlacementTables_EntriesEntry" => ("AbstractString", "PlacementTable"), "PlacementTable_LoadMapEntry" => ("AbstractString", "Host")


text:
using Documenter
using SVDD

makedocs(
    sitename = "SVDD Documentation",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "start.md",
        "Custom Solvers" => [
            "SMO" => "smo.md"
            ]
        ],
    format = :html,
    modules = [SVDD]
)
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#
deploydocs(
    repo = "github.com/holtri/SVDD.jl.git",
)


text:
"""
sample!([rng, ]a, k)

Sample `k` element from array `a` without repetition and eventually excluding elements in `exclude`.

### Optional Arguments
- `exclude=()`: elements in `a` to exclude from sampling.

### Implementation Notes
Changes the order of the elements in `a`. For a non-mutating version, see [`sample`](@ref).
"""
function sample!(rng::AbstractRNG, a::AbstractVector, k::Integer; exclude=())
    minsize = k + length(exclude)
    length(a) < minsize && throw(ArgumentError("vector must be at least size $minsize"))
    res = Vector{eltype(a)}()
    sizehint!(res, k)
    n = length(a)
    i = 1
    while length(res) < k
        r = rand(rng, 1:(n - i + 1))
        if !(a[r] in exclude)
            push!(res, a[r])
            a[r], a[n - i + 1] = a[n - i + 1], a[r]
            i += 1
        end
    end
    res
end

sample!(a::AbstractVector, k::Integer; exclude=()) = sample!(getRNG(), a, k; exclude=exclude)

"""
    sample([rng,] r, k)

Sample `k` element from unit range `r` without repetition and eventually excluding elements in `exclude`.

### Optional Arguments
- `exclude=()`: elements in `a` to exclude from sampling.

### Implementation Notes
Unlike [`sample!`](@ref), does not produce side effects.
"""
sample(a::AbstractVector, k::Integer; exclude=()) = sample!(getRNG(), collect(a), k; exclude=exclude)

getRNG(seed::Integer=-1) = seed >= 0 ? MersenneTwister(seed) : GLOBAL_RNG

"""
    insorted(item, collection)

Return true if `item` is in sorted collection `collection`.

### Implementation Notes
Does not verify that `collection` is sorted.
"""
function insorted(item, collection)
    index = searchsortedfirst(collection, item)
    @inbounds return (index <= length(collection) && collection[index] == item)
end

"""
    generate_min_set(g, gen_func, Reps)
Generate a vector `Reps` times using `gen_func(g)` and return the vector with the least elements.
"""
generate_min_set(g::AbstractGraph{T}, gen_func, Reps::Integer) where T<: Integer =
mapreduce(gen_func, (x, y)->length(x)<length(y) ? x : y, Iterators.repeated(g, Reps))

"""
    generate_max_set(g, gen_func, Reps)
Generate a vector `Reps` times using `gen_func(g)` and return the vector with the most elements.
"""
generate_max_set(g::AbstractGraph{T}, gen_func, Reps::Integer) where T<: Integer =
mapreduce(gen_func, (x, y)->length(x)>length(y) ? x : y, Iterators.repeated(g, Reps))

"""
    greedy_contiguous_partition(weight, required_partitions, num_items=length(weight))

Partition `1:num_items` into atmost `required_partitions` number of contiguous partitions with
the objective of minimising the largest partition.
The size of a partition is equal to the num of the weight of its elements.
`weight[i] > 0`.

### Performance
Time: O(num_items+required_partitions)
Requires only one iteration over `weight` but may not output the optimal partition.

### Implementation Notes
`Balance(wt, left, right, n_items, n_part) =
max(sum(wt[left:right])*(n_part-1), sum(wt[right+1:n_items]))`.
Find `right` that minimises `Balance(weight, 1, right, num_items, required_partitions)`.
Set the first partition as `1:right`.
Repeat on indices `right+1:num_items` and one less partition.
"""
function greedy_contiguous_partition(
    weight::Vector{<:Integer},
    required_partitions::Integer,
    num_items::U=length(weight)
    ) where U <: Integer

    suffix_sum = cumsum(reverse(weight))
    reverse!(suffix_sum)
    push!(suffix_sum, 0) #Eg. [2, 3, 1] => [6, 4, 1, 0]

    partitions = Vector{UnitRange{U}}()
    sizehint!(partitions, required_partitions)

    left = one(U)
    for partitions_remain in reverse(1:(required_partitions-1))

        left >= num_items && break

        partition_size = weight[left]*partitions_remain #At least one item in each partition
        right = left

        #Find right: sum(wt[left:right])*partitions_remain and sum(wt[(right+1):num_items]) is balanced
        while right+one(U) < num_items && partition_size < suffix_sum[right+one(U)]
            right += one(U)
            partition_size += weight[right]*partitions_remain
        end
        #max( sum(wt[left:right]), sum(wt[(right+1):num_items]) ) = partition_size
        #max( sum(wt[left:(right-1)]), sum(wt[right:num_items]) ) = suffix_sum[right]
        if left != right && partition_size > suffix_sum[right]
            right -= one(U)
        end

        push!(partitions, left:right)
        left = right + one(U)
    end
    push!(partitions, left:num_items)

    return partitions
end

"""
    optimal_contiguous_partition(weight, required_partitions, num_items=length(weight))

Partition `1:num_items` into atmost `required_partitions` number of contiguous partitions such
that the largest partition is minimised.
The size of a partition is equal to the sum of the weight of its elements.
`weight[i] > 0`.

### Performance
Time: O(num_items*log(sum(weight)))

### Implementation Notes
Binary Search for the partitioning over `[fld(sum(weight)-1, required_partitions), sum(weight)]`.
"""
function optimal_contiguous_partition(
    weight::Vector{<:Integer},
    required_partitions::Integer,
    num_items::U=length(weight)
    ) where U <: Integer

    item_it = Iterators.take(weight, num_items)

    up_bound = sum(item_it) # Smallest known possible value
    low_bound = fld(up_bound-1, required_partitions) # Largest known impossible value

    # Find optimal balance
    while up_bound > low_bound+1
        search_for = fld(up_bound+low_bound, 2)

        sum_part = 0
        remain_part = required_partitions

        possible = true
        for w in item_it
            sum_part += w
            if sum_part > search_for
                sum_part = w
                remain_part -= 1
                if remain_part == 0
                    possible = false
                    break
                end
            end
        end
        if possible
            up_bound = search_for
        else
            low_bound = search_for
        end
    end
    best_balance = up_bound

    # Find the partition with optimal balance
    partitions = Vector{UnitRange{U}}()
    sizehint!(partitions, required_partitions)
    sum_part = 0
    left = 1
    for (i, w) in enumerate(item_it)
        sum_part += w
        if sum_part > best_balance
            push!(partitions, left:(i-1))
            sum_part = w
            left = i
        end
    end
    push!(partitions, left:num_items)

    return partitions
end

"""
    isbounded(n)

Returns true if `typemax(n)` of an integer `n` exists.
"""
isbounded(n::Integer) = true
isbounded(n::BigInt) = false

"""
    isbounded(T)

Returns true if `typemax(T)` of a type `T <: Integer` exists.
"""
isbounded(::Type{T}) where {T <: Integer} = isconcretetype(T)
isbounded(::Type{BigInt}) = false


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

MOIU.@model(IpoptModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

# Without fixed_variable_treatment set, duals are not computed for variables
# that have lower_bound == upper_bound.
const ipopt_optimizer = Ipopt.Optimizer(print_level=0, fixed_variable_treatment="make_constraint")
const ipopt_config = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                               optimal_status=MOI.LOCALLY_SOLVED)

@testset "SolverName" begin
    @test MOI.get(ipopt_optimizer, MOI.SolverName()) == "Ipopt"
end

@testset "MOI Linear tests" begin
    exclude = ["linear8a", # Behavior in infeasible case doesn't match test.
               "linear12", # Same as above.
               "linear8b", # Behavior in unbounded case doesn't match test.
               "linear8c", # Same as above.
               "linear7",  # VectorAffineFunction not supported.
               "linear15", # VectorAffineFunction not supported.
               ]
    model_for_ipopt = MOIU.UniversalFallback(IpoptModelData{Float64}())
    linear_optimizer = MOI.Bridges.SplitInterval{Float64}(
                         MOIU.CachingOptimizer(model_for_ipopt, ipopt_optimizer))
    MOIT.contlineartest(linear_optimizer, ipopt_config, exclude)
end

MOI.empty!(ipopt_optimizer)

@testset "MOI QP/QCQP tests" begin
    qp_optimizer = MOIU.CachingOptimizer(IpoptModelData{Float64}(), ipopt_optimizer)
    MOIT.qptest(qp_optimizer, ipopt_config)
    exclude = ["qcp1", # VectorAffineFunction not supported.
              ]
    MOIT.qcptest(qp_optimizer, ipopt_config, exclude)
end

MOI.empty!(ipopt_optimizer)

@testset "MOI NLP tests" begin
    MOIT.nlptest(ipopt_optimizer, ipopt_config)
end


text:
###
# Tainting of values coming from random number generators
###

# Cassette.@primitive function rand(::Type{T}) where {T <: Integer, __CONTEXT__<:TraceCtx}
#     resevoir = __trace__.metadata.rands
#     if length(resevoir) >= 1
#         vv = pop!(resevoir)::T
#     else
#         vv = rand(T)
#     end
#     sv = Sym(:rand, typeof(vv))
#     return Cassette.Box(__trace__.context, vv, sv)
# end


text:
# Code for video 2_4:

Pkg.init()
# INFO: Initializing package repository C:\Users\CVO\.julia\v0.4
# INFO: Cloning METADATA from git://github.com/JuliaLang/METADATA.jl

Pkg.add("IJulia")
# INFO: Updating cache of ZMQ...
# INFO: Installing BinDeps v0.3.19
# INFO: Installing Compat v0.7.7
# INFO: Installing Conda v0.1.8
# INFO: Installing IJulia v1.1.8
# INFO: Installing JSON v0.5.0
# INFO: Installing LibExpat v0.1.0
# INFO: Installing Nettle v0.2.0
# INFO: Installing SHA v0.1.2
# INFO: Installing URIParser v0.1.1
# INFO: Installing WinRPM v0.1.13
# INFO: Installing ZMQ v0.3.1
# INFO: Installing Zlib v0.1.12
# INFO: Building WinRPM
# INFO: Recompiling stale cache file C:\Users\CVO\.julia\lib\v0.4\Compat.ji for mo
# dule Compat.
# INFO: Recompiling stale cache file C:\Users\CVO\.julia\lib\v0.4\URIParser.ji for
#  module URIParser.

Pkg.status()
# 1 required packages:
#  - IJulia                        1.1.8
# 11 additional packages:
#  - BinDeps                       0.3.19
#  - Compat                        0.7.7
#  - Conda                         0.1.8
#  - JSON                          0.5.0
#  - LibExpat                      0.1.0
#  - Nettle                        0.2.0
#  - SHA                           0.1.2
#  - URIParser                     0.1.1
#  - WinRPM                        0.1.13
#  - ZMQ                           0.3.1
#  - Zlib                          0.1.12

using IJulia
# INFO: Recompiling stale cache file C:\Users\CVO\.julia\lib\v0.4\IJulia.ji for mo
# dule IJulia.
# INFO: Recompiling stale cache file C:\Users\CVO\.julia\lib\v0.4\ZMQ.ji for modul
# e ZMQ.
# INFO: Recompiling stale cache file C:\Users\CVO\.julia\lib\v0.4\Nettle.ji for mo
# dule Nettle.

Pkg.update()
# INFO: Updating METADATA...
# INFO: Computing changes...
# INFO: No packages to install, update or remove

Pkg.pin("HDF5", v"0.4.3")

Pkg.clone("git@github.com:EricChiang/ANN.jl.git")

Pkg.add("Winston")
using Winston

# example 1:
plot(randn(100))

# example 2:
plot(x -> x * sin(3x) * exp(-0.3x), 0, 3pi)

# example 3:
x = linspace(0, 10, 100)
g = gamma(x)
p = FramedPlot(
         title="Gamma function",
         xlabel="x-coord",
         ylabel="y-coord")
add(p, Curve(x, g, color="green"))

max_stars_repo_path: Section 2/2_4.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/2_4.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/2_4.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video