text:
"""
    Class

A type used for controlling which classes of `TimeZone` are considered valid. Instances of
`Class` can be combined using bitwise operators to generate a mask which allows multiple
classes to be considered valid at once.

Currently supported masks are:

- `Class(:FIXED)`: Class indicating the time zone name is parsable as a fixed UTC offset.
- `Class(:STANDARD)`: The time zone name is included in the primary IANA tz source files.
- `Class(:LEGACY)`: The time zone name is included in the deprecated IANA tz source files.
- `Class(:NONE)`: Mask that will match nothing.
- `Class(:DEFAULT)`: Default mask used by functions: `Class(:FIXED) | Class(:STANDARD)`
- `Class(:ALL)`: Mask allowing all supported classes.
"""
struct Class
    val::UInt8
end

function Class(name::Symbol)
    if name === :NONE
        Class(0x00)
    elseif name === :FIXED
        Class(0x01)
    elseif name === :STANDARD
        Class(0x02)
    elseif name === :LEGACY
        Class(0x04)
    elseif name === :DEFAULT
        Class(:FIXED) | Class(:STANDARD)
    elseif name === :ALL
        Class(:FIXED) | Class(:STANDARD) | Class(:LEGACY)
    else
        throw(ArgumentError("Unknown class name: $name"))
    end
end

function Class(str::AbstractString, regions::AbstractSet{<:AbstractString})
    class = Class(:NONE)
    occursin(FIXED_TIME_ZONE_REGEX, str) && (class |= Class(:FIXED))
    !isempty(intersect(regions, TZData.STANDARD_REGIONS)) && (class |= Class(:STANDARD))
    !isempty(intersect(regions, TZData.LEGACY_REGIONS)) && (class |= Class(:LEGACY))
    return class
end

Class(str::AbstractString, regions::AbstractVector) = Class(str, Set{String}(regions))

Base.:(|)(a::Class, b::Class) = Class(a.val | b.val)
Base.:(&)(a::Class, b::Class) = Class(a.val & b.val)
Base.:(==)(a::Class, b::Class) = a.val == b.val
Base.:(~)(a::Class) = Class(~a.val)

function Base.show(io::IO, mask::Class)
    C = repr(Class)
    mask == Class(:NONE) && return print(io, "$C(:NONE)")

    names = String[]
    mask & Class(:FIXED) == Class(:FIXED) && push!(names, "$C(:FIXED)")
    mask & Class(:STANDARD) == Class(:STANDARD) && push!(names, "$C(:STANDARD)")
    mask & Class(:LEGACY) == Class(:LEGACY) && push!(names, "$C(:LEGACY)")

    unused = mask & ~Class(:ALL)
    unused != Class(:NONE) && push!(names, "$C($(repr(unused.val)))")

    join(io, names, " | ")
end


text:
module Optimise

export train!,
	SGD, Descent, ADAM, Momentum, Nesterov, RMSProp,
	ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, ADAMW,
	InvDecay, ExpDecay, WeightDecay, stop, Optimiser

include("optimisers.jl")
include("train.jl")
include("deprecations.jl")

end

ext: jl


text:
using Mocking

create_cluster_1_patch = @patch function create_cluster(a...; b...)
    return Dict{String, Any}()
end

delete_cluster_1_patch = @patch function delete_cluster(a...; b...)
    return Dict{String, Any}()
end

describe_cluster_1_patch = @patch function describe_cluster(a...; b...)
    return Dict{String, Any}("cluster" => Dict{String, Any}("roleArn" => "arn:aws:iam::012345678910:role/eks-service-role-AWSServiceRoleForAmazonEKS-J7ONKE3BQ4PI", "name" => "devel", "certificateAuthority" => Dict{String, Any}("data" => "LS0tLS1CRUdJTiBDRVJUSUZJQ0FURS0tLS0tCk1JSUN5RENDQWJDZ0F3SUJBZ0lCQURBTkJna3Foa2lHOXcwQkFRc0ZBREFWTVJNd0VRWURWUVFERXdwcmRXSmwKY201bGRHVnpNQjRYRFRFNE1EVXpNVEl6TVRFek1Wb1hEVEk0TURVeU9ESXpNVEV6TVZvd0ZURVRNQkVHQTFVRQpBeE1LYTNWaVpYSnVaWFJsY3pDQ0FTSXdEUVlKS29aSWh2Y05BUUVCQlFBRGdnRVBBRENDQVFvQ2dnRUJBTTZWCjVUaG4rdFcySm9Xa2hQMzRlVUZMNitaRXJOZGIvWVdrTmtDdWNGS2RaaXl2TjlMVmdvUmV2MjlFVFZlN1ZGbSsKUTJ3ZURyRXJiQyt0dVlibkFuN1ZLYmE3ay9hb1BHekZMdmVnb0t6b0M1N2NUdGVwZzRIazRlK2tIWHNaME10MApyb3NzcjhFM1ROeExETnNJTThGL1cwdjhsTGNCbWRPcjQyV2VuTjFHZXJnaDNSZ2wzR3JIazBnNTU0SjFWenJZCm9hTi8zODFUczlOTFF2QTBXb0xIcjBFRlZpTFdSZEoyZ3lXaC9ybDVyOFNDOHZaQXg1YW1BU0hVd01aTFpWRC8KTDBpOW4wRVM0MkpVdzQyQmxHOEdpd3NhTkJWV3lUTHZKclNhRXlDSHFtVVZaUTFDZkFXUjl0L3JleVVOVXM3TApWV1FqM3BFbk9RMitMSWJrc0RzQ0F3RUFBYU1qTUNFd0RnWURWUjBQQVFIL0JBUURBZ0trTUE4R0ExVWRFd0VCCi93UUZNQU1CQWY4d0RRWUpLb1pJaHZjTkFRRUxCUUFEZ2dFQkFNZ3RsQ1dIQ2U2YzVHMXl2YlFTS0Q4K2hUalkKSm1NSG56L2EvRGt0WG9YUjFVQzIrZUgzT1BZWmVjRVZZZHVaSlZCckNNQ2VWR0ZkeWdBYlNLc1FxWDg0S2RXbAp1MU5QaERDSmEyRHliN2pVMUV6VThTQjFGZUZ5ZFE3a0hNS1E1blpBRVFQOTY4S01hSGUrSm0yQ2x1UFJWbEJVCjF4WlhTS1gzTVZ0K1Q0SU1EV2d6c3JRSjVuQkRjdEtLcUZtM3pKdVVubHo5ZEpVckdscEltMjVJWXJDckxYUFgKWkUwRUtRNWEzMHhkVWNrTHRGQkQrOEtBdFdqSS9yZUZPNzM1YnBMdVoyOTBaNm42QlF3elRrS0p4cnhVc3QvOAppNGsxcnlsaUdWMm5SSjBUYjNORkczNHgrYWdzYTRoSTFPbU90TFM0TmgvRXJxT3lIUXNDc2hEQUtKUT0KLS0tLS1FTkQgQ0VSVElGSUNBVEUtLS0tLQo="), "status" => "ACTIVE", "resourcesVpcConfig" => Dict{String, Any}("subnetIds" => Any["subnet-6782e71e", "subnet-e7e761ac"], "vpcId" => "vpc-950809ec", "securityGroupIds" => Any["sg-6979fe18"]), "endpoint" => "https://A0DCCD80A04F01705DD065655C30CC3D.yl4.us-west-2.eks.amazonaws.com", "version" => "1.10", "arn" => "arn:aws:eks:us-west-2:012345678910:cluster/devel", "createdAt" => 1.527807879988e9))
end

list_clusters_1_patch = @patch function list_clusters(a...; b...)
    return Dict{String, Any}("clusters" => Any["devel", "prod"])
end

list_tags_for_resource_1_patch = @patch function list_tags_for_resource(a...; b...)
    return Dict{String, Any}("tags" => Dict{String, Any}("aws:tag:domain" => "beta"))
end


text:
# example of 5 taxa tree in ipad notes
# used to test createHybrid!
# Claudia September 2014
#
# in julia: include("tree_example.jl")

ed11=Edge(11,1.5);
ed12=Edge(12,0.2);
ed3=Edge(3,0.9);
ed4=Edge(4,0.1);
ed5=Edge(5,0.2);
ed6=Edge(6,0.1);
ed7=Edge(7,0.1);
ed8=Edge(8,0.1);
ed9=Edge(9,0.1);
ed10=Edge(10,0.1);
ed13=Edge(13,0.1);

n1=Node(1,false,false,[ed11,ed5,ed6]);
n4=Node(4,true,false,[ed11]);
n5=Node(5,false,false,[ed12,ed5,ed9]);
n6=Node(6,true, false,[ed6]);
n7=Node(7,true, false,[ed12]);
n8=Node(8,true, false,[ed8]);
n9=Node(9,false,false,[ed8,ed9,ed10]);
n10=Node(10,true, false,[ed10]);

setNode!(ed5,[n5,n1]);
setNode!(ed6,[n1,n6]);
setNode!(ed8,[n8,n9]);
setNode!(ed9,[n5,n9]);
setNode!(ed10,[n9,n10]);
setNode!(ed11,[n1,n4]);
setNode!(ed12,[n5,n7]);

net=HybridNetwork([n1,n4,n5,n6,n7,n8,n9,n10],[ed5,ed6,ed8,ed9,ed10,ed11,ed12]);


text:
"""
```
DSGEVECM{T} <: AbstractDSGEVECMModel{T}
```
implements a simple interface for combining a given DSGE model
with a VECM (a VAR with cointegrating terms) to create a DSGE-VECM.
Confer with Del Negro and Schorfheide (2006)
and/or Del Negro, Schorfheide, Smets, and Wouters (2007) for details about DSGE-VECMs.
We recommend the first paper as an initial introduction to DSGE-VECMs.

The recommended constructor requires the user to provide
(1) an `AbstractDSGEModel` object, (2) which structural shocks from
the DSGE to use, and (3) the subspec (optional, defaults to "ss0").
If the subspec "ss0" is used, then the result is a `DSGEVECM`
whose VECM component is "empty" in that the
observables, cointegrating relationships, lags, and λ weight are not specified.
The reason why this constructor requires the user to specify which
structural shocks of DSGE to use is that this information is
DSGE-specific rather than information about the VECM.

However, we can also construct a `DSGEVECM` without having to
specify the structural shocks when calling the constructor,
although we still need to give an instance of an `AbstractDSGEModel`.

### Example
The code below instantiates an empty `DSGEVECM`
with `AnSchorfheide` as the underlying DSGE and then
calls `update!` on the empty DSGE-VECM
to add information about the desired DSGE-VECM spec.

```jldoctest; output = false
dsgevecm = DSGEVECM(AnSchorfheide())
DSGE.update!(dsgevecm, shocks = [:rm_sh, :z_sh, :g_sh],
    observables = [:obs_gdp, :obs_cpi, :obs_nominalrate],
    λ = 1., lags = 4)

# output

DSGE-VECM Model
observables:      3
data_vintage:     200310
DSGE model:       an_schorfheide
DSGE description: Julia implementation of model defined in 'Bayesian Estimation of DSGE Models' by Sungbae An and Frank Schorfheide: AnSchorfheide, ss0
```

### Fields

#### DSGE object
* `dsge::AbstractDSGEModel{T}`: underlying DSGE model object

#### DSGE-VECM Information
* `observables::OrderedDict{Symbol,Int}`: dictionary mapping observables
    of the VECM to their index in the matrices representing the DSGE-VECM.
* `cointegrating::OrderedDict{Symbol,Int}`: dictionary mapping cointegrating relationships
    of the VECM to their index in the matrices representing the DSGE-VECM.
    When creating the state space representation of a DSGE-VECM, the
    cointegrating relationships will come after the observables. Accordingly,
    the first `cointegrating` index will be after the last observable.
* `cointegrating_add::OrderedDict{Symbol,Int}`: dictionary mapping
    additional cointegrating relationships of the VECM to their index
    in the matrices representing the DSGE-VECM. These relationships are
    intended to be ones that only add to the `DD` matrix in DSGE.jl's
    representation of a state space model. Importantly, these relationships
    do not require changing the `ZZ` matrix. As a result, the indices for
    these relationships start at 1.
* `shocks::OrderedDict{Symbol,Int}`: dictionary mapping structural
    shocks in the DSGE to their index in the matrices representing the DSGE-VECM
* `lags::Int`: number of lags in the VECM
* `λ::T`: weight on the DSGE prior

#### Auxiliary Information
* `spec::String`: concatenates `dsgevar` with the spec of the DSGE, e.g.
    for `AnSchorfheide`, we have `dsgevar_an_schorfheide`.
* `subspec::String`: specifies the model subspecification.
    Cached here for filepath computation.
* `testing::Bool`: indicates whether the model is in testing mode.
    Currently, this setting has no uses in practice
"""
mutable struct DSGEVECM{T} <: AbstractDSGEVECMModel{T}
    dsge::AbstractDSGEModel{T}
    observables::OrderedDict{Symbol,Int}
    cointegrating::OrderedDict{Symbol,Int}
    cointegrating_add::OrderedDict{Symbol,Int}
    shocks::OrderedDict{Symbol,Int}
    lags::Int
    λ::T
    spec::String
    subspec::String
    testing::Bool
end

function Base.show(io::IO, m::DSGEVECM)
    @printf io "DSGE-VECM Model\n"
    @printf io "observables:      %i\n" n_observables(m)
    @printf io "data_vintage:     %s\n" data_vintage(m.dsge)
    @printf io "DSGE model:       %s\n" spec(get_dsge(m))
    @printf io "DSGE description: %s\n" description(get_dsge(m))
end

function DSGEVECM(dsge::AbstractDSGEModel{T}, shocks::Vector{Symbol}, subspec::String = "ss0";
                 custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                 copy_dsge::Bool = false, testing = false) where {T<:Real}

    # Initialize specs
    spec     = "dsgevecm_" * ModelConstructors.spec(dsge)
    subspec  = subspec

    # Initialize empty DSGEVECM
    if copy_dsge
        dsge = deepcopy(dsge)
    end
    m = DSGEVECM{T}(dsge, OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
                    OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), 0,
                    0., spec, subspec, testing)

    # Initialize shocks from DSGE
    update!(m; shocks = shocks, check_valid = false)

    # Initialize subspec
    init_subspec!(m)

    # Do checks of observables, shocks, and lags
    check_valid_dsgevecm(m)

    return m
end

# Empty constructor
function DSGEVECM(dsge::AbstractDSGEModel{T}, subspec::String = "ss0";
                 custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                 copy_dsge::Bool = false, testing = false) where {T<:Real}

    # Initialize specs
    spec     = "dsgevar_" * ModelConstructors.spec(dsge)
    subspec  = subspec

    # Initialize empty DSGEVECM
    if copy_dsge
        dsge = deepcopy(dsge)
    end
    m = DSGEVECM{T}(dsge, OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
                    OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), 0,
                    0., spec, subspec, testing)

    # Initialize shocks from DSGE
    update!(m; check_valid = false)

    # Initialize subspec
    init_subspec!(m)

    # Do checks of observables, shocks, and lags
    check_valid_dsgevecm(m)

    return m
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "renumber element nodes" begin
    mesh = Mesh()
    add_element!(mesh, 1, :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
    mapping = Dict{Symbol, Vector{Int}}(
        :Tet10 => [1, 2, 4, 3, 5, 6, 7, 8, 9, 10])
    reorder_element_connectivity!(mesh, mapping)
    @test mesh.elements[1] == [1, 2, 4, 3, 5, 6, 7, 8, 9, 10]
    invmapping = Dict{Symbol, Vector{Int}}()
    invmapping[:Tet10] = invperm(mapping[:Tet10])
    reorder_element_connectivity!(mesh, invmapping)
    @test mesh.elements[1] == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
end


text:
# Sequences to Matrix
# -------------------

"""
    seqmatrix(vseq::AbstractVector{BioSequence{A}}, major::Symbol) where {A<:Alphabet}

Construct a matrix of nucleotides or amino acids from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{BioSequences.BioSequence{BioSequences.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site)
4x3 Array{BioSequences.DNA,2}:
 DNA_A  DNA_A  DNA_A
 DNA_T  DNA_T  DNA_T
 DNA_C  DNA_C  DNA_C
 DNA_G  DNA_G  DNA_G

 julia> seqmatrix(seqs, :seq)
 3x4 Array{BioSequences.DNA,2}:
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
```
"""
function seqmatrix(vseq::AbstractVector{LongSequence{A}}, major::Symbol) where {A<:Alphabet}
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length"))
    end
    if major == :site
        mat = Matrix{eltype(A)}(undef, (nseqs, nsites))
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = vseq[seq][site]
        end
        return mat
    elseif major == :seq
        mat = Matrix{eltype(A)}(undef, (nsites, nseqs))
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = vseq[seq][site]
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

"""
    seqmatrix(::Type{T}, vseq::AbstractVector{BioSequence{A}}, major::Symbol) where {T,A<:Alphabet}

Construct a matrix of `T` from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{BioSequences.BioSequence{BioSequences.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site, UInt8)
4×3 Array{UInt8,2}:
 0x01  0x01  0x01
 0x08  0x08  0x08
 0x02  0x02  0x02
 0x04  0x04  0x04

julia> seqmatrix(seqs, :seq, UInt8)
3×4 Array{UInt8,2}:
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
```
"""
function seqmatrix(::Type{T}, vseq::AbstractVector{LongSequence{A}}, major::Symbol) where {T,A<:Alphabet}
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length."))
    end
    if major == :site
        mat = Matrix{T}(undef, (nseqs, nsites))
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = convert(T, vseq[seq][site])
        end
        return mat
    elseif major == :seq
        mat = Matrix{T}(undef, (nsites, nseqs))
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = convert(T, vseq[seq][site])
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

# Consensus
# ---------

"""
    majorityvote(seqs::AbstractVector{LongSequence{A}}) where {A<:NucleicAcidAlphabet}

Construct a sequence that is a consensus of a vector of sequences.

The consensus is established by a simple majority vote rule, where ambiguous
nucleotides cast an equal vote for each of their possible states.
For each site a winner(s) out of A, T(U), C, or G is determined, in the cases
of ties the ambiguity symbol that unifies all the winners is returned.
E.g if A and T tie, then W is inserted in the consensus. If all A, T, C, and G
tie at a site, then N is inserted in the consensus. Note this means that if a
nucletide e.g. 'C' and a gap '-' draw, the nucleotide will always win over the
gap, even though they tied.

# Examples
```julia
julia> seqs = [dna"CTCGATCGATCC", dna"CTCGAAAAATCA", dna"ATCGAAAAATCG", dna"ATCGGGGGATCG"]

4-element Array{BioSequences.BioSequence{BioSequences.DNAAlphabet{4}},1}:
 CTCGATCGATCC
 CTCGAAAAATCA
 ATCGAAAAATCG
 ATCGGGGGATCG

julia> majorityvote(seqs)
12nt DNA Sequence:
MTCGAAARATCG
```
"""
function majorityvote(seqs::AbstractVector{LongSequence{A}}) where {A<:NucleicAcidAlphabet}
    mat = seqmatrix(UInt8, seqs, :site)
    nsites = size(mat, 2)
    nseqs = size(mat, 1)
    result = BioSequence{A}(nsites)
    votes = Array{Int}(undef, 16)
    @inbounds for site in 1:nsites
        fill!(votes, 0)
        for seq in 1:nseqs
            nuc = mat[seq, site]
            votes[1] += nuc == 0x00
            votes[2] += (nuc & 0x01) != 0x00
            votes[3] += (nuc & 0x02) != 0x00
            votes[5] += (nuc & 0x04) != 0x00
            votes[9] += (nuc & 0x08) != 0x00
        end
        m = maximum(votes)
        merged = 0x00
        for i in 0x01:0x10
            merged |= ifelse(votes[i] == m, i - 0x01, 0x00)
        end
        result[site] = reinterpret(eltype(A), merged)
    end
    return result
end


text:
struct RPCond
    density
    pressure
    velocity
    sound_speed
    energy
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
# Literal value.

mutable struct LiteralNode <: AbstractSQLNode
    val

    LiteralNode(; val) =
        new(val)
end

LiteralNode(val) =
    LiteralNode(val = val)

"""
    Lit(; val)
    Lit(val)

A SQL literal.

In a scalar context, `missing`, numbers, strings and datetime values are
automatically converted to SQL literals.

# Examples

```jldoctest
julia> q = Select(:null => missing,
                  :boolean => true,
                  :integer => 42,
                  :text => "SQL is fun!",
                  :date => Date(2000));

julia> print(render(q))
SELECT NULL AS "null", TRUE AS "boolean", 42 AS "integer", 'SQL is fun!' AS "text", '2000-01-01' AS "date"
```
"""
Lit(args...; kws...) =
    LiteralNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Lit), pats::Vector{Any}) =
    dissect(scr, LiteralNode, pats)

Base.convert(::Type{AbstractSQLNode}, val::SQLLiteralType) =
    LiteralNode(val)

Base.convert(::Type{AbstractSQLNode}, ref::Base.RefValue) =
    LiteralNode(ref.x)

PrettyPrinting.quoteof(n::LiteralNode, ctx::QuoteContext) =
    Expr(:call, nameof(Lit), n.val)


text:
print("""
#include "giac/giac.h"
#include "giac_c.h"
#include<string>

using namespace std;
using namespace giac;

extern "C" {

    static _Gen* _error(const string &msg)
    {
        cerr << msg << endl;
        return giac_undef();
    }

""")

open("giac_funs.txt","r") do f
for line in eachline(f)
fun=strip(line)
if fun[1]=='#'   #comment
    continue
end

print("""

    _Gen* giac_$fun(_Gen* g, void *cp)
    {
        gen* r0;
        try{
            gen* g0 = reinterpret_cast<gen*>(g);
            context* cp0 = reinterpret_cast<context*>(cp);
            r0 = new gen();
            *r0 = _$fun(*g0, cp0);
        }
        catch(runtime_error &e) {return _error(e.what());}
        return reinterpret_cast<_Gen*>(r0);
    }
""")

end
end

println("""

}
""")


text:
const Vehicle1D = Entity{State1D, VehicleDef, Int}
const Scene1D = Frame{Vehicle1D}
Scene1D(n::Int=100) = Frame(Vehicle1D, n)
Scene1D(arr::Vector{Vehicle1D}) = Frame{Vehicle1D}(arr, length(arr))

get_center(veh::Vehicle1D) = veh.state.s
get_footpoint(veh::Vehicle1D) = veh.state.s
get_front(veh::Vehicle1D) = veh.state.s + veh.def.length/2
get_rear(veh::Vehicle1D) = veh.state.s - veh.def.length/2

function get_headway(veh_rear::Vehicle1D, veh_fore::Vehicle1D, roadway::StraightRoadway)
    return get_headway(get_front(veh_rear), get_rear(veh_fore), roadway)
end
function get_neighbor_fore(scene::Scene1D, vehicle_index::Int, roadway::StraightRoadway)
    ego = scene[vehicle_index]
    best_ind = nothing
    best_gap = Inf
    for (i,veh) in enumerate(scene)
        if i != vehicle_index
            Δs = get_headway(ego, veh, roadway)
            if Δs < best_gap
                best_gap, best_ind = Δs, i
            end
        end
    end
    return NeighborLongitudinalResult(best_ind, best_gap)
end
function get_neighbor_rear(scene::Scene1D, vehicle_index::Int, roadway::StraightRoadway)
    ego = scene[vehicle_index]
    best_ind = nothing
    best_gap = Inf
    for (i,veh) in enumerate(scene)
        if i != vehicle_index
            Δs = get_headway(veh, ego, roadway)
            if Δs < best_gap
                best_gap, best_ind = Δs, i
            end
        end
    end
    return NeighborLongitudinalResult(best_ind, best_gap)
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
struct ValueRewardTest <: TrainingTest
    idx::Integer
    state::AbstractState
    image::AbstractImage
    value::Float32
    reward::Float32
end

mutable struct ValueRewardTests <: TrainingTests
    tests::Vector{ValueRewardTest}
    images::AbstractArray{Float32, 4}
    actual_value::AbstractArray{Float32, 2}
    actual_reward::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function ValueRewardTests(tests; period=1, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            Flux.unsqueeze(get_value.(tests), 1),
            Flux.unsqueeze(get_reward.(tests), 1),
            period,
            warmup,
        )
    end
end

function gettestfn(tests::ValueRewardTests)
    function fn(n::Network)
        initial_output = initialinferencefn(tests)(n)
        value_mse = getvaluemse(tests, initial_output)
        reward_mse = getrewardmse(
            tests,
            recurrentinferencefn(tests, initial_output)(n)
        )
        value_mse, reward_mse
    end
end

struct ValueRewardTestsSVG
    tests::ValueRewardTests
end
svg(x::ValueRewardTests) = ValueRewardTestsSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::ValueRewardTestsSVG)
    num_tests = length(x.tests)
    range = 1:(min(10, num_tests))
    subset = 1:num_tests
    subset_indices = subset[range]
    svgbatch = [gridsvg(x.tests[i]) for i in subset_indices]
    Base.show(io, m, GridSVG(svgbatch, main_title="ValueReward Tests"))
end


text:
module MicroArchitectures

const augment = """
    using Base.BinaryPlatforms: arch, arch_march_isa_mapping, CPUID, HostPlatform, Platform

    function augment_microarchitecture!(platform::Platform)
        haskey(platform, "march") && return platform

        host_arch = arch(HostPlatform())
        host_isas = arch_march_isa_mapping[host_arch]
        idx = findlast(((name, isa),) -> isa <= CPUID.cpu_isa(), host_isas)
        platform["march"] = first(host_isas[idx])
        return platform
    end
    """

end

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Documenter, LibPQ, Memento

setlevel!(getlogger(LibPQ), "critical")

makedocs(;
    modules=[LibPQ],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Type Conversions" => "pages/type-conversions.md",
        "API" => "pages/api.md",
        "FAQ" => "pages/faq.md",
    ],
    repo="https://github.com/invenia/LibPQ.jl/blob/{commit}{path}#L{line}",
    sitename="LibPQ.jl",
    checkdocs=:exports,
    linkcheck=true,
    strict=true,
    authors="Eric Davies",
)

deploydocs(;
    repo="github.com/invenia/LibPQ.jl",
)


text:
using BenchmarkTools

# see this section of docs https://docs.julialang.org/en/v1/base/arrays/#Broadcast-and-vectorization-1

function G_single(t::T, params::Vector{T}) where T<:Real

   G₀, G₁, η₁, G₂, η₂ = params

   G = G₀ + G₁*exp(-t*G₁/η₁) + G₂*exp(-t*G₂/η₂)

end

function G_vec(t::Vector{T}, params::Vector{T}) where T<:Real

   G₀, G₁, η₁, G₂, η₂ = params

   G = G₀ .+ G₁*exp.(-t*G₁/η₁) .+ G₂*exp.(-t*G₂/η₂)

end

function G_loop(t::Vector{T}, params::Vector{T}) where T<:Real

    out = Array{T}(undef, length(t))

    @inbounds for (i,ti) in enumerate(t)
        out[i] = G_single(ti, params)
    end

    out
end

t = collect(0.0:0.1:10000)
params = [1.0, 1.0, 1.0, 1.0, 1.0]

# always worst by far, ignore for now
# println("test vectorised version of function")
# @btime G_vec($t, $params)
# println("\n")

# println("test broadcasting using . notation")
# @btime G_single.($t, $(params,))
# println("\n")

# println("test broadcast function explicitly (should be the same as . notation)")
# @btime broadcast($G_single, $t, $(params,))
# println("\n")

# println("test using array comprehension")
# @btime [G_single(ti, $params) for ti in $t]
# println("\n")

println("test using map")
@btime map(G_single, t, (params,))
println("\n")

# println("test using for loop")
# @btime G_loop($t, $params)
# println("\n")

println("Benchmarks complete")

max_stars_repo_path: 20190319_broadcast benchmarks.jl

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_path: 20190319_broadcast benchmarks.jl

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_path: 20190319_broadcast benchmarks.jl

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
export natural_adaptive_walk, random_adaptive_walk, greedy_adaptive_walk,
  reluctant_adaptive_walk, neutral_walk, neutral_or_fitter_walk,
  fitter_then_neutral_walk, fitness_range_walk

# TODO: Add a tabu walk that looks at its own history, possibly only as far back as when it got to a certain level of fitness, and doesn't revisit nodes.
# TODO: Add an optional limit on walk length.
# TODO: Another tabu strategy is to avoid mutating a given locus more than once every k iterations.
# TODO: Consider optionally doing a depth-first search with backtracking to ensure total exploration of a plateau before declaring it a local optimum.

@doc """A record of the path of an adaptive or random walk.

`strategy` - Can be one of :natural, :random, :greedy, :reluctant,
:fitness_range, :fitter_then_neutral for now.

`length` - The number of transitions from the original genotype to the ultimate
optimum. This will be one less than the width of the `history_list` matrix.

`history_list` - Vector of genotypes found along the path from the original to
the optimum, stored in the order they were observed. The first element is
always the starting genotype.

`history_set` - Dictionary of genotypes found along the path from the original
to the optimum. Used in plateau search to avoid repeated states
"""
type Walk
  strategy::Symbol
  length::Int64
  history_list::Vector{Genotype}
  history_set::Set{Genotype}
end

Walk(s::Symbol, g::Genotype) = Walk(s, 0, [g], Set{Genotype}([g]))

function add_step!(w::Walk, g::Genotype)
  w.length += 1
  w.history_list = [w.history_list; g]
  push!(w.history_set, g)
end

function natural_adaptive_walk(g::Genotype)
end

@doc """An adaptive walk in which a random, fitter neighbor is chosen at each
step.
"""
function random_adaptive_walk(g::Genotype)
  g0 = g
  w = Walk(:random, g0)
  while true
    nbrs = fitter_neighbors(g0)
    if length(nbrs) == 0
      break
    end
    g0 = nbrs[rand(1:end)]
    add_step!(w, g0)
  end
  return w
end

@doc """An adaptive walk in which the fittest neighbor is chosen at each step.
"""
function greedy_adaptive_walk(g::Genotype)
  g0 = g
  w = Walk(:greedy, g0)
  while true
    nbrs = fitter_neighbors(g0)
    if length(nbrs) == 0
      break
    end
    g0 = nbrs[end]

    add_step!(w, g0)
  end
  return w
end

@doc """An adaptive walk in which the least fit, fitter neighbor is chosen at
each step.
"""
function reluctant_adaptive_walk(g::Genotype)
  g0 = g
  w = Walk(:random, g0)
  while true
    nbrs = fitter_neighbors(g0)
    if length(nbrs) == 0
      break
    end
    g0 = nbrs[1]

    add_step!(w, g0)
  end
  return w
end

function generic_neutral_walk(g::Genotype, walktype::Symbol, nbrfunc::Function, maxsteps::Int64)
  g0 = g
  w = Walk(walktype, g0)

  step = 1
  while true
    nbrs = nbrfunc(g0)
    if length(nbrs) == 0
      break
    end

    # Choose a new genotype from neutral neighbors
    # or terminate if there isn't one we haven't yet
    # visited.
    i = findfirst(j -> !(nbrs[j] in w.history_set), 1:length(nbrs))
    if i == 0
      break
    end
    g0 = nbrs[i]
    add_step!(w, g0)

    step += 1
    if maxsteps > 0 && step > maxsteps
      break
    end
  end
  return w
end

@doc """
A purely neutral walk in which a random, equally fit neighbor is chosen at each
step.

Note that this walk is unlikely to generate a path of length greater than one
with the `NKLandscape` because it requires exact fitness equality.

To avoid cycles, this walk will not revisit genotypes. It uses the
`history_set` attribute of the `Walk` type to enforce this condition.

The resulting walk will have no more than `maxsteps` steps.
"""
function neutral_walk(g::Genotype, maxsteps::Int64=0)
  generic_neutral_walk(g, :neutral, neutral_neighbors, maxsteps)
end

@doc """
An approximatly neutral walk in which a random neighbor with fitness in the
closed interval [lb, ub] is chosen at each step.  This walk is suitable for
generating an approximately neutral walk with an `NKLandscape` because it
does not require exact fitness equality.

To avoid cycles, this walk will not revisit genotypes. It uses the
`history_set` attribute of the `Walk` type to enforce this condition.

The resulting walk will have no more than `maxsteps` steps.
"""
function fitness_range_walk(g::Genotype, lb::Float64, ub::Float64, maxsteps::Int64=0)
  nbrfunc = (g) -> fitness_range_neighbors(g, lb, ub)
  generic_neutral_walk(g, :fitness_range, nbrfunc, maxsteps)
end

@doc """
A walk in which a random neighbor is chosen from the set of fitter or equal
neighbors at each step. Note that this walk is equivalent to a random adaptive
walk for an `NKLandscape` because the likelihood of equal fitness is
vanishingly small for an `NKLandsape`.

To avoid cycles, this walk will not revisit genotypes. It uses the
`history_set` attribute of the `Walk` type to enforce this condition.

The resulting walk will have no more than `maxsteps` steps.
"""
function neutral_or_fitter_walk(g::Genotype, maxsteps::Int64=0)
  # TODO: Considering optimizing so history_set is cleared when fitness changes.
  nbrfunc = (g) -> fitter_or_equal_neighbors(g, sorted=false)
  generic_neutral_walk(g, :neutral, nbrfunc, maxsteps)
end

@doc """
A walk in which a random fitter neighbor is chosen if there is one. If there is
no fitter neighbor, then a random neutral neighbor is chosen. Note that this
walk is equivalent to a random adaptive walk for an `NKLandscape` because the
likelihood of equal fitness is vanishingly small for an NKLandsape.

To avoid cycles, this walk will not revisit genotypes at the same fitness
level. It uses the `history_set` attribute of the `Walk` type to enforce this
condition.

The resulting walk will have no more than `maxsteps` steps.
"""
function fitter_then_neutral_walk(g::Genotype, maxsteps::Int64=0)
  # TODO: Considering optimizing so history_set is cleared when fitness changes.
  function nbrfunc(g)
    nbrs = fitter_neighbors(g, sorted=false)
    if length(nbrs) > 0
      return nbrs
    end
    return neutral_neighbors(g)
  end
  generic_neutral_walk(g, :fitter_then_neutral, nbrfunc, maxsteps)
end


text:
function Grid=HexagonGrid(NumNodes,OrientFace)

Grid.Dim=3;
Grid.Type='Hexa';

Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
for k=1:NumNodes
  xx=cos(2*pi*k/NumNodes+pi/NumNodes);
  yy=sin(2*pi*k/NumNodes+pi/NumNodes);
  Nodes(NodeNumber)=Node([xx;yy;0],NodeNumber);
  NodeNumber=NodeNumber+1;
end
Grid.Nodes=Nodes;
Grid.NumNodes=NumNodes;

NumEdges=NumNodes;
Edges(1:NumEdges)=Edge([1 2],Grid,0,0);
N1=1;
N2=2;
EdgeNumber=1;
for i=1:NumEdges-1
  Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
  N1=N1+1;
  N2=N2+1;
end
Edges(EdgeNumber)=Edge([N1 1],Grid,EdgeNumber,EdgeNumber);
Grid.Edges=Edges;
Grid.NumEdges=NumEdges;

NumFaces=1;
Faces(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
Type='o';
[Faces(1),Grid]=Face(1:NumEdges,Grid,1,Type,OrientFace);
Grid.Faces=Faces;
Grid.NumFaces=NumFaces;
end


text:
using SparseArrays
using LinearAlgebra

@testset "matrixnetwork" begin
    A = load_matrix_network("dfs_example")
    M = MatrixNetwork(A)
    B = sparse(M)
    @test A == B
    C = sparse_transpose(M)
    @test A' == C
    @test size(M) == size(B)
    @test_throws DomainError size(M, 0)
    @test size(M, 4) == 1
    @test ndims(M) == 2

    load_matrix_network_all("minnesota")
    load_matrix_network_all("U3A")

    #@show issym(sparse([0 1; 0 0]'))
    @test is_undirected(MatrixNetwork(sparse([0 1; 1 0]))) == true
    @test is_undirected(MatrixNetwork(sparse([0 1; 0 1]))) == false
    @test is_undirected(MatrixNetwork(load_matrix_network("dfs_example"))) == false
    @test is_undirected(sparse([0 1; 1 0])) == true

    @test is_connected(havel_hakimi_graph([1,1,1,1])) == false
    @test is_connected(havel_hakimi_graph([2,2,2])) == true
    @test is_connected(sparse([0 1 0; 0 0 1; 1 0 0])) == true
    @test is_connected(spzeros(1,1)) == true
    @test is_connected(empty_graph(0)) == false
    @test is_connected(empty_graph(1)) == true
    @test is_connected(empty_graph(2)) == false

    @test is_empty(empty_graph(0)) == true
    @test is_empty(MatrixNetwork(Int[],Int[],0)) == true
    @test is_empty(MatrixNetwork(Int[1],Int[1])) == false
    @test is_empty(erdos_renyi_undirected(0,0)) == true
    @test is_empty(erdos_renyi_undirected(1,0)) == false

    @testset "matvec" begin
        G = erdos_renyi_directed(5,2)
        y = collect(1:5)
        A = sparse(G)
        @test norm(A*y - G*y,2) <= 10*eps(1.0)
        @test norm(G'*y - A'*y,2) <= 10*eps(1.0)
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
VERSION

# %%
using RCall

x = 1

f(var::String) = eval(Meta.parse("@rput ("*var*")"))

module M1
using RCall
import RCall: @rput

f(var::String) = eval(Meta.parse("@rput ("*var*")"))

end

# %%
f("x")

# %%
M1.f("x")

# %%
M1.f("Main.x")

# %%
using RCall

x = 1

f(var::String) = eval(:(@rput($(Symbol(var)))))

module M2

using RCall

x = 2

f(var::String) = eval(:(@rput($(Symbol(var)))))

g(var::String) = Main.eval(:(@rput($(Symbol(var)))))
g(m::Module, var::String) = Core.eval(m, :(@rput($(Symbol(var)))))

end

# %%
M2.f("x")

# %%
R"x"

# %%
M2.g("x")

# %%
R"x"

# %%
M2.g(M2, "x")

# %%
R"x"

# %% [markdown]
# ```julia
# using RCall
#
# x = 1
#
# f(var::String) = eval(:(@rput($(Symbol(var)))))
#
# module M2
#
# using RCall
#
# x = 2
#
# f(var::String) = eval(:(@rput($(Symbol(var)))))
#
# g(var::String) = Main.eval(:(@rput($(Symbol(var)))))
#
# g(m::Module, var::String) = Core.eval(m, :(@rput($(Symbol(var)))))
#
# end
# ```
#
# ```julia
# julia> M2.f("x")
# 2
#
# julia> R"x"
# RObject{IntSxp}
# [1] 2
# ```
#
# ```julia
# julia> M2.g("x")
# 1
#
# julia> R"x"
# RObject{IntSxp}
# [1] 1
# ```
#
# ```julia
# julia> M2.g(M2, "x") # equivalent to M2.f("x")
# 2
#
# julia> R"x"
# RObject{IntSxp}
# [1] 2
# ```
#

# %%
using RCall

x = 1

f(var::String) = eval(:(@rput($(Symbol(var)))))

module M2

using RCall

x = 2

f(var::String) = eval(:(@rput($(Symbol(var)))))

g(var::String) = Main.eval(:(@rput($(Symbol(var)))))

g(m::Module, var::String) = Core.eval(m, :(@rput($(Symbol(var)))))

end

# %%
using RCall

x = 1

f(var::String) = eval(:(@rput($(Symbol(var)))))

module M2

using RCall

x = 2

f(var::String) = eval(:(@rput($(Symbol(var)))))
g(var::String) = Main.eval(:(@rput($(Symbol(var)))))
g(m::Module, var::String) = Core.eval(m, :(@rput($(Symbol(var)))))

end

# %%

max_stars_repo_path: 0008/eval example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0008/eval example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0008/eval example.jl

max_forks_repo_name: genkuroki/public


text:
# Utilities for `TimeArray`.

# Import modules.

using Dates, TimeSeries

import Base.map

"""
Align two timeseries according to their times, padding values from earlier times.

# Arguments
- `x: TimeArray{X,1,T,Array{X,1}}`: the first timeseries
- `y: TimeArray{Y,1,T,Array{Y,1}}`: the first timeseries

# Example
```
x = TimeArray(
    [Time(0), Time(1), Time(3), Time(5)],
    [    0.1,     1.1,     3.1,     5.1]
)
y = TimeArray(
    [Time(0), Time(1), Time(2), Time(5), Time(6)],
    [    0.1,     1.1,     2.1,     5.1,     6.1]
)
z = align(x, y)
```
"""
function aligntimes(x :: TimeArray{X,1,T,Array{X,1}}, y :: TimeArray{Y,1,T,Array{Y,1}}) :: TimeArray{Tuple{X,Y},1,T,Array{Tuple{X,Y},1}} where X where Y where T <: TimeType
    nx = length(x)
    ny = length(y)
    ix = 1
    iy = 1
    xt = timestamp(x)
    yt = timestamp(y)
    zt = Array{T,1}()
    xv = values(x)
    yv = values(y)
    zv = Array{Tuple{X,Y},1}()
    if nx == 0 || ny == 0
        throw(ArgumentError("align: TimeArrays must not be empty."))
    elseif xt[1] != yt[1]
        throw(ArgumentError("align: TimeArrays must start at the same timestamp."))
    end
    while ix <= nx || iy <= ny
        if ix <= nx && (iy > ny || xt[ix] < yt[iy])
            push!(zt, xt[ix])
            push!(zv, (xv[ix], yv[iy-1]))
            ix = ix + 1
        elseif iy <= ny && (ix > nx || xt[ix] > yt[iy])
            push!(zt, yt[iy])
            push!(zv, (xv[ix-1], yv[iy]))
            iy = iy + 1
        else
            push!(zt, xt[ix])
            push!(zv, (xv[ix], yv[iy]))
            ix = ix + 1
            iy = iy + 1
        end
    end
    TimeArray(zt, zv)
end

"""
Compute the time difference corresponding to a fractional number of hours.

# Arguments
- `delta` : the fractional number of hours
"""
function fractionalhours(delta) :: Nanosecond
    onehour = Time(1) - Time(0)
    if false # Old implementation.
        resolution = 60 * 60 * 1000 # Resolve to the nearest millisecond.
        onehour * trunc(Int64, resolution * delta) / resolution
    else # New implementation.
        typeof(onehour)(trunc(Int64, delta * onehour.value))
    end
end

"""
Add a fractional number of hours to a time.

# Arguments
- `t :: Time`: the time
- `delta`    : the number of hours to add
"""
function addhours(t :: Time, delta) :: Time
    onehour = Time(1) - Time(0)
    t1 = t.instant / onehour + delta
    if t1 > 24
        Time(23, 59, 59, 999, 999, 999)
    else
        t + fractionalhours(delta)
    end
end

"""
Add a fractional number of hours to a time.

# Arguments
- `t :: Time`: the time
- `delta`    : the number of hours to add
"""
function addhours(t :: TimeType, delta) :: TimeType
    onehour = Time(1) - Time(0)
    t + delta * onehour
end

"""
Map the values in a time array.

# Arguments
- `f`                                : the function to apply to the values
- `x :: TimeArray{X,1,T,Array{x,1}}` : the time array
"""
function map(f, x :: TimeArray{X,1,T,Array{X,1}}) where X where T <: TimeType
    xt = timestamp(x)
    xv = values(x)
    TimeArray(xt, map(f, xv))
end


text:
function read_raster(path::String, T)
    # Check if file exists (ArchGDAL error is cryptic)
    check_path = endswith(path, ".gz") ? path[10:lastindex(path)] : path
    !isfile(check_path) && error("the file \"$(check_path)\" does not exist")

    raw = ArchGDAL.unsafe_read(path)
    transform = ArchGDAL.getgeotransform(raw)
    wkt = ArchGDAL.getproj(raw)

    # Extract 1st band (should only be one band anyway)
    # to get a 2D array instead of 3D
    band = ArchGDAL.getband(raw, 1)

    # Extract the array
    array_t = ArchGDAL.read(band)

    # This handles UInt tiff rasters that can still have negative NoData values
    # Need to convert the NoData value to Int64 in these cases
    if eltype(array_t) <: Integer
        ras_type = Int64
    else
        ras_type = eltype(array_t)
    end

    # Extract no data value, first converting it to the proper type (based on
    # the raster). Then, need to convert to T. Weird, yes,
    # but it's the only way I could get it to work for all raster types... -VL
    nodata_val = convert(T, convert(ras_type, ArchGDAL.getnodatavalue(band)))

    # Transpose the array -- ArchGDAL returns a x by y array, need y by x
    # Line to handle NaNs in datasets read from tifs
    array_t[isnan.(array_t)] .= nodata_val

    array = convert(Array{Union{Missing, T}}, permutedims(array_t, [2, 1]))

    array[array .== nodata_val] .= missing

    # Close connection to dataset
    ArchGDAL.destroy(raw)

    array, wkt, transform # wkt and transform are needed later for write_raster
end

# Write a single band raster, either in .tif or .asc format,
# inspired by GeoArrays.write()
# Intended to only write rasters with no values < 0, so NoData as -9999 will
# always be safe
function write_raster(fn_prefix::String,
                      array::Array{T, 2} where T <: Number,
                      wkt::String,
                      transform,
                      file_format::String)
    dtype = eltype(array)

    # transponse array back to columns by rows
    array_t = permutedims(array, [2, 1])

    width, height = size(array_t)

    # Define extension and driver based in file_format
    file_format == "tif" ? (ext = ".tif"; driver = "GTiff") :
            (ext = ".asc"; driver = "AAIGrid")

    file_format == "tif" ? (options = ["COMPRESS=LZW"]) :
                           (options = [])

    # Append file extention to filename
    fn = string(fn_prefix, ext)

    # Create raster in memory *NEEDED* because no create driver for .asc
    ArchGDAL.create(fn_prefix,
                    driver = ArchGDAL.getdriver("MEM"),
                    width = width,
                    height = height,
                    nbands = 1,
                    dtype = dtype,
                    options = options) do dataset
        band = ArchGDAL.getband(dataset, 1)
        # Write data to band
        ArchGDAL.write!(band, array_t)

        # Write nodata and projection info
        ArchGDAL.setnodatavalue!(band, -9999)
        ArchGDAL.setgeotransform!(dataset, transform)
        ArchGDAL.setproj!(dataset, wkt)

        # Copy memory object to disk (necessary because ArchGDAL.create
        # does not support creation of ASCII rasters)
        ArchGDAL.destroy(ArchGDAL.copy(dataset,
                                       filename = fn,
                                       driver = ArchGDAL.getdriver(driver),
                                       options = options))
    end

end


text:
# core mechanism of extending Taine Zhao's @thautwarm 's MatchCore pattern matching.

using MatchCore

## compile (quote) left and right hands of a rule
# escape symbols to create MLStyle compatible patterns

# compile left hand of rule
# if it's the first time seeing v, add it to the "seen symbols" Set
# insert a :& expression only if v has not been seen before
function c_left(v::Symbol, s)
    if Base.isbinaryoperator(v) return v end
    (v ∉ s ? (push!(s, v); v) : amp(v)) |> dollar
end
c_left(v::Expr, s) = v.head ∈ add_dollar ? dollar(v) : v
# support symbol literals in left hand
c_left(v::QuoteNode, s) = v.value isa Symbol ? dollar(v) : v
c_left(v, s) = v # ignore other types

c_right(v::Symbol) = Base.isbinaryoperator(v) ? v : dollar(v)
function c_right(v::Expr)
    v.head ∈ add_dollar ? dollar(v) : v
end
c_right(v::QuoteNode) = v.value isa Symbol ? v.value : v
c_right(v) = v #ignore other types

# add dollar in front of the expressions with those symbols as head
const add_dollar = [:(::), :(...)]
# don't walk down on these symbols
const skips = [:(::), :(...)]

# Compile rules from Metatheory format to MatchCore format
function compile_rule(rule::Rule)::Expr
    le = df_walk(c_left, rule.left, Vector{Symbol}(); skip=skips, skip_call=true) |> quot

    if rule.mode == :dynamic # regular pattern matching
        # right side not quoted! needed to evaluate expressions in right hand.
        re = rule.right
    elseif rule.mode == :rewrite || rule.mode == :equational
		# right side is quoted, symbolic replacement
        re = df_walk(c_right, rule.right; skip=skips, skip_call=true) |> quot
	else
        error(`rule "$e" is not in valid form.\n`)
    end

    return :($le => $re)
end

# catch-all for reductions
const identity_axiom = :($(quot(dollar(:i))) => i)

# TODO analyse theory before compiling. Identify associativity and commutativity
# and other loop-creating problems. Generate a pattern matching block with the
# correct rule order and expansions for associativity and distributivity
# import Iterators: flatten

function theory_block(t::Vector{Rule})
	tn = Vector{Expr}()

	for r ∈ t
		push!(tn, compile_rule(r))
		if r.mode == :equational
			mirrored = Rule(r.right, r.left, r.expr, r.mode, nothing)
			push!(tn, compile_rule(mirrored))
		end
	end

	block(tn..., identity_axiom)
end

"""
Compile a theory to a closure that does the pattern matching job
Returns a RuntimeGeneratedFunction, which does not use eval and
is as fast as a regular Julia anonymous function 🔥
"""
function compile_theory(theory::Vector{Rule}, mod::Module; __source__=LineNumberNode(0))
    # generate an unique parameter name
    parameter = Meta.gensym(:reducing_expression)
    block = theory_block(theory)

    matching = MatchCore.gen_match(parameter, block, __source__, mod)
    matching = MatchCore.AbstractPatterns.init_cfg(matching)

    ex = :(($parameter) -> $matching)
    closure_generator(mod, ex)
end

"""
Compile a theory at runtime to a closure that does the pattern matching job
"""
macro compile_theory(theory)
    gettheory(theory, __module__)
end

# Retrieve a theory from a module at compile time. Not exported
function gettheory(var, mod; compile=true)
	t = nothing
    if Meta.isexpr(var, :block) # @matcher begine rules... end
		t = rmlines(macroexpand(mod, var)).args .|> Rule
	else
		t = mod.eval(var)
	end

	if compile && !(t isa Function)
		t = compile_theory(t, mod)
	end

	return t
end


text:
module ParametrisedConvexApproximators

using Flux
using UnPack
using Convex
using SCS  # 211207; currently not compatible with apple silicon
using Optim
using Transducers
using Random

# approximators
export AbstractApproximator, FNN
export ConvexApproximator, MA, LSE
export ParametrisedConvexApproximator, PMA, PLSE
export PICNN, project_nonnegative!
export optimise, number_of_parameters
# data manipulation
export partitionTrainTest

# approximators
include("approximators/approximators.jl")
include("optimise.jl")
include("data_manipulation/data_manipulation.jl")

end  # module


text:
using PyCall
win32 = pyimport_conda("win32com.client", "pywin32", "conda-forge")
genc = win32.gencache
igor = genc.EnsureDispatch("IgorPro.Application")
oo = igor._oleobj_

igor.

pycom = genc.pythoncom
using Markdown
h = genc.GetModuleForProgID("IgorPro.Application") .|> propertynames
using ElectronDisplay
ElectronDisplay.niceh = Mark

max_stars_repo_name: pacific-zitter/GingerGUI

max_issues_repo_name: pacific-zitter/GingerGUI

max_forks_repo_name: pacific-zitter/GingerGUI


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry
using DiskTracer.trace

using LaTeXStrings
using Plots
pyplot()

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

# Plot the vertical profile of Upsilon for a bunch of radii.

# Set up the grid interpolators
nr=128
rmax = 500 * AU
nz=128
zmax = 200 * AU

# choose some radii we want to sample
zs = linspace(0., 4 * AU, nz)
Upsilons = Float64[DiskTracer.model.Upsilon_nu(1. * AU, z, mol.nu_0, pars, mol) for z in zs]
p1 = plot(zs./AU, Upsilons, ylabel=L"$\Upsilon$", xlabel=L"$r_\mathrm{cyl}$")

zs = linspace(0., 10 * AU, nz)
Upsilons = Float64[DiskTracer.model.Upsilon_nu(10. * AU, z, mol.nu_0, pars, mol) for z in zs]
p2 = plot(zs./AU, Upsilons, ylabel=L"$\Upsilon$", xlabel=L"$r_\mathrm{cyl}$")

zs = linspace(0., 40 * AU, nz)
Upsilons = Float64[DiskTracer.model.Upsilon_nu(100. * AU, z, mol.nu_0, pars, mol) for z in zs]
p3 = plot(zs./AU, Upsilons, ylabel=L"$\Upsilon$", xlabel=L"$r_\mathrm{cyl}$")

zs = linspace(0., 150 * AU, nz)
Upsilons = Float64[DiskTracer.model.Upsilon_nu(300. * AU, z, mol.nu_0, pars, mol) for z in zs]
p4 = plot(zs./AU, Upsilons, ylabel=L"$\Upsilon$", xlabel=L"$r_\mathrm{cyl}$")

p = plot(p1,p2,p3,p4,layout=(4,1), legend=false, size=(800,1200))

savefig(p, "Upsilon-traces.png")


text:
using Revise
st0 = time_ns()

import CSV
import DataFrames
import DataFrames: DataFrame
import Random
import OnlineStats
import StatsBase

using Distributed # only for Tests?
using Statistics
using Pkg # only for Tests?
using Test # only for Tests?
using DecisionTrees

# set directories
global pkgdir = ""
global testdir = joinpath(dirname(@__FILE__), "") # joinpath(@__DIR__,"..","test",filename)
# the following if clause is for "direct execution" (e.g. in VScode where @__FILE__ does not seem to work)
if testdir == ""
    try
        global pkgdir = Pkg.dir("DecisionTrees") # how long will this remain supported?
        global testdir = joinpath(pkgdir, "test")
    catch
    end
else
    global pkgdir = normpath(joinpath(testdir, ".."))
end
datadir = joinpath(pkgdir, "data")

Random.seed!(1239946)

t0runtests = time_ns()
println("DTM: Running tests...") # Time = $(t0runtests)")

global eltypesData1 = [Int,	Float64,	Float64,	Float64,	Float64,	Float64,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	String,	String,	String,	String,	String,	String,	Int,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	String,	Int,	String,	Int,	String,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	Int,	Int,	Int,	String,	String,	String,	String,	String,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Float64,	Int,	Float64,	Float64,	Float64,	Float64,	Int,	Int,	Int,	Int,	Int,	Int]

import DataFrames
import Random: rand,randstring

##################################################
# Read the data
##################################################

thisfile = joinpath(datadir, "GermanMotorPremiums", "data1Small.csv")
@test isfile(thisfile)
@time df_tmp = CSV.read(thisfile, allowmissing=:none, types=eltypesData1, categorical=false, rows_for_type_detect=10000);
df_tmp_orig = deepcopy(df_tmp)

selected_explanatory_vars = ["PLZ_WOHNORT","ART_DES_WOHNEIGENTUM","GEBURTSDATUM","FAMILIENSTAND","NATIONALITAET","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]

# keep 10 largest PLZ only
vorig = deepcopy(df_tmp[:PLZ_WOHNORT])
counts, freqs, vals, keep = getCounts(vorig, threshold=10)
vnew = df_tmp[:PLZ_WOHNORT]
mapToOther!(vnew,keep,9999999)

##################################################
# Run Mini Tree
##################################################
dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);
sett.minWeight = -.2
# strs,resm=dtm(dtmtableMini,sett)
# boosting
settB = deepcopy(sett)
settB.model_type = "boosted_tree"
strs, resm = dtm(dtmtableMini, settB)

# selected_explanatory_vars=[    "VORSCHAEDEN_ANZAHL",    "MALLORCA_POLICE",	"SCHUTZBRIEF_INKL",	"FREIE_WERKSTATTWAHL",	"AUTOMOBILCLUB_MITGLIED_SEIT",	"BAHNCARD",	"ZAHLUNGSWEISE",	"JAHRESKARTE_OEPNV",	"MOTORRAD_BESITZER",	"AUTOMOBILCLUB",	"SFKLASSE_VOLLKASKO",	"SFKLASSE_HAFTPFLICHT",	"STELLPLATZ_ABSCHLIESSBAR",	"NAECHTLICHER_STELLPLATZ",	"NUTZUNGSWEISE",	"JAEHRLICHE_FAHRLEISTUNG",	"TSN",	"ERSTZULASSUNG",	"HSN",	"FINANZIERUNGSART",	"ZULASSUNG_AUF_VERSICHERUNGSNEHM",	"STADT",	"KENNZEICHEN",	"PLZ_DES_HALTER",	"SELBSTGENUTZTES_WOHNEIGENTUM",	"ART_DES_WOHNEIGENTUM",	"GEBURTSDATUM",	"FAMILIENSTAND",	"NATIONALITAET",	"GESCHLECHT",	"FUEHRERSCHEIN_ERWORBEN_AM",	"VORSCHAEDEN0_typeKH",	"VORSCHAEDEN0_typetk",	"VORSCHAEDEN0_month",	"VORSCHAEDEN0_year",	"VORSCHAEDEN1_typetk",	"VORSCHAEDEN1_month",	"VORSCHAEDEN1_year",	"VORSCHAEDEN2_typevk",	"VORSCHAEDEN2_month",	"VORSCHAEDEN2_year",	"adacid",	"name",	"marke",	"modell",	"preis",	"getriebeart",	"antriebsart",	"Fahrzeugklasse",	"co2klasse",	"kw",	"ps",	"tueranzahl",	"Motorart",	"Kraftstoffart",	"Motorbauart",	"Schadstoffklasse",	"Karosserie",	"Sitzanzahl",	"typklasseh_num",	"typklassetk_num",	"typklassevk_num",	"hubraum2",	"drehmoment2",	"breite2",	"radstand2",	"laenge2",	"hoehe2",	"leergewicht2",	"gesamtgewicht2",	"zuladung2",	"kofferraumvolumen_num",	"hoechstgeschwindigkeit2",	"verbrauchgesamt2",	"verbrauchausserorts2",	"verbrauchinnerorts2",	"beschleunigung2",	"tank2",	"kfzsteuer2",	"anzahlgaenge2",	"anzahlzylinder2",	"co2_wert",	"modellstart_y"]
##################################################
# run tree
##################################################
dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);
# strs,resm=dtm(dtmtable,sett)

##################################################
# run boosting
##################################################

sett.iterations = 10
sett.model_type = "boosted_tree"
# strs,resm2=dtm(dtmtable,sett)

# reduce iterations for the remaining tests
sett.iterations = 3

sett.statsByVariables = Int[1,2]
# strs,resm2=dtm(dtmtable,sett)
sett.statsByVariables = Int[]

sett.model_type = "boosted_tree"
cvsampler = CVOptions(-3, 0.0, true)
@info("DTM:testing cvsampler")

try
    statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)
catch e
    println("a")
    @show e
    @assert false
end

cvsampler = CVOptions(3, 0.5, false)
try
    statsdf, settsdf, cvModels = dtm(dtmtableMini, sett, cvsampler)
catch e
    println("b")
    @show e
    @assert false
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    HalfEdge(head, cell, prev, next, half)

Stores the indices of the `head` vertex, the `prev`
and `next` edges in the left `cell`, and the opposite
`half`-edge. For some half-edges the `cell` may be
`nothing`, e.g. border edges of the mesh.

See [`HalfEdgeStructure`](@ref) for more details.
"""
mutable struct HalfEdge
  head::Int
  cell::Union{Int,Nothing}
  prev::HalfEdge
  next::HalfEdge
  half::HalfEdge
  HalfEdge(head, cell) = new(head, cell)
end

function Base.show(io::IO, e::HalfEdge)
  print(io, "HalfEdge($(e.head), $(e.cell))")
end

"""
    HalfEdgeStructure(halfedges, edgeoncell, edgeonvertex)

A data structure for orientable 2-manifolds based
on half-edges.

Two types of half-edges exist (Kettner 1999). This
implementation is the most common type that splits
the incident cells.

A vector of `halfedges` together with a vector of
`edgeoncell` and a vector of `edgeonvertex` can be
used to retrieve topolological relations in optimal
time. In this case, `edgeonvertex[i]` returns the
index of the half-edge in `halfedges` with head equal
to `i`. Similarly, `edgeoncell[i]` returns the index
of a half-edge in `halfedges` that is in the cell `i`.

Such data structure is usually constructed from another
data structure such as [`ExplicitStructure`](@ref) via
`convert` methods:

```julia
he = convert(HalfEdgeStructure, structure)
```

See also [`TopologicalStructure`](@ref).

## References

* Kettner, L. (1999). [Using generic programming for
  designing a data structure for polyhedral surfaces]
  (https://www.sciencedirect.com/science/article/pii/S0925772199000073)
"""
struct HalfEdgeStructure <: TopologicalStructure
  halfedges::Vector{HalfEdge}
  edgeoncell::Vector{Int}
  edgeonvertex::Vector{Int}
end

"""
    halfedges(s)

Return the half-edges of the half-edge structure `s`.
"""
halfedges(s::HalfEdgeStructure) = s.halfedges

"""
    edgeoncell(s, c)

Return a half-edge of the half-edge structure `s` on the `c`-th cell.
"""
edgeoncell(s::HalfEdgeStructure, c) = s.halfedges[s.edgeoncell[c]]

"""
    edgeonvertex(s, v)

Return the half-edge of the half-edge structure `s` for which the
head is the `v`-th index.
"""
edgeonvertex(s::HalfEdgeStructure, v) = s.halfedges[s.edgeonvertex[v]]

"""
    ncells(s)

Return the number of cells in the half-edge structure `s`.
"""
ncells(s::HalfEdgeStructure) = length(s.edgeoncell)


text:
using Logging

@testset "Graph traversal" begin
    g = StorageGraph()
    @test walkdep(g, (a=1,)=>(b=1,)) == ((a=1,), Set{Int}())
    StorageGraphs.add_vertex!(g, (a=1,))
    @test nextid(g, (a=1,)=>(b=1,)) == 1
    add_nodes!(g, (a=1,)=>(b=1,))
    @test walkdep(g, (a=1,)=>(b=1,)) == ((b=1,), Set(1))

    dep = (a=1,)=>(b=1,)=>(c=1,)
    @test walkdep(g, dep) == ((b=1,), Set(1))
    @test nextid(g, dep) == 1
    add_nodes!(g, dep)
    @test walkdep(g, dep) == ((c=1,), Set(1))
    @test walkdep(g, dep, stopcond=(g,v)->has_prop(g,g[v],:b)) == ((b=1,), Set(1))

    add_nodes!(g, (a=1,)=>(b=1,)=>(c=2,))
    add_nodes!(g, (a=1,)=>(b=2,)=>(c=1,))
    v = walkpath(g, 1, g[(a=1,)])
    @test g[v] == (c=1,)
    @test length(v) == 1
    v = walkpath(g, 1, g[(a=1,)], stopcond=(g,v)->has_prop(g,v,:b))
    @test g[v] == (b=1,)
    a(g,v,n) = @debug v
    @test @test_logs((:debug, 1), (:debug, 2), (:debug, 3),
        min_level=Logging.Debug, match_mode=:all,
        StorageGraphs.walkpath!(g, 1, g[(a=1,)], outneighbors, a)) == 3

    g = StorageGraph()
    dep = ((a=1,),(b=1,),(c=1,),(d=1,))
    add_bulk!(g, foldr(=>, dep), (e=[1,2],))
    add_nodes!(g, (a=1,)=>(b=1,)=>(c=2,)=>(d=1,))
    f_dep = foldr(=>, (dep..., (e=3,)))
    @test walkdep(g, f_dep) == ((d=1,), Set([1,2]))
    @test nextid(g, f_dep) == 4
end


text:
using GraphTools
using Base.Test

g = @graph A -> B
@test 1 == g["A"].index
@test "A" == g["A"].label

g = @graph "A" -> "B"
@test 1 == g["A"].index
@test "A" == g["A"].label


text:
#=
This file is run several times
* with grad=Base vs grad=Dual
* with Tracker, Zygote
* using KernelAbstractions, LoopVectorization, TensorOperations
=#

using Tullio, Test, ForwardDiff, Random
# using Tracker; _gradient(x...) = Tracker.gradient(x...); GRAD = :Tracker; macro printline() end

function gradtest(f, dims)
    x = randn(dims...)
    grad = ForwardDiff.gradient(x -> sum(sin, f(x)), x)
    grad ≈ _gradient(x -> sum(sin, f(x)), x)[1]
end

@testset "simple" begin

    @test _gradient(x -> sum(@tullio y[i] := 2*x[i]), rand(3))[1] == [2,2,2]
    @test _gradient(x -> sum(@tullio y[i] := 2*x[i] + i), rand(3))[1] == [2,2,2]

    # two contributions
    g2(x) = @tullio y[i, j] := 1 * x[i] + 1000 * x[j]  avx=false
    mat = [1 1 3; 1 1 5; 7 7 7]
    g_fd = ForwardDiff.gradient(x -> sum(mat .* g2(x)), rand(3))
    @test g_fd ≈ _gradient(x -> sum(mat .* g2(x)), rand(3))[1]

    # larger array, no shared indices -- https://github.com/mcabbott/Tullio.jl/issues/14
    r100 = randn(100)
    g_fd = ForwardDiff.gradient(x -> sum(sin, g2(x)), r100)
    @test g_fd ≈ _gradient(x -> sum(sin, g2(x)), r100)[1]

    r100 = randn(100)

    # scalar output
    s2(x) = @tullio s := exp(x[i]) / x[j]
    @test _gradient(s2, r100)[1] ≈ ForwardDiff.gradient(s2, r100)

    # two arrays, and a sum
    h2(x,y) = @tullio z[i] := x[i,j] + y[j,i]
    @test _gradient(sum∘h2, rand(2,3), rand(3,2)) == (ones(2,3), ones(3,2))

    # nontrivial function
    flog(x,y) = @tullio z[i] := log(x[i,j]) / y[j,i]  avx=false  # new failure LoopVectorization v0.12.14? only on CI?
    r_x, r_y = rand(2,3), rand(3,2)
    fx = ForwardDiff.gradient(x -> sum(flog(x, r_y)), r_x)
    fy = ForwardDiff.gradient(y -> sum(flog(r_x, y)), r_y)
    @test fx ≈ _gradient(sum∘flog, r_x, r_y)[1]
    @test fy ≈ _gradient(sum∘flog, r_x, r_y)[2]

    # classic
    mm(x,y) = @tullio z[i,j] := 2 * x[i,k] * y[k,j]  avx=false # new?
    x1 = rand(3,4);
    y1 = rand(4,5);
    z1 = x1 * y1
    dx, dy = _gradient(sum∘mm, x1, y1)
    @test dx ≈ 2 * ones(3,5) * y1'
    @test dy ≈ 2 * x1' * ones(3,5)

    # abs, abs2
    va = [1,-2,3,-4,5]
    abs_grad = ForwardDiff.gradient(v -> sum(abs, 1 .+ v.^2), va)
    @test abs_grad ≈ _gradient(v -> (@tullio s := abs(1 + v[i]^2)), va)[1]
    abs2_grad = ForwardDiff.gradient(v -> sum(abs2, 1 .+ v.^2), va)
    @test abs2_grad ≈ _gradient(v -> (@tullio s := abs2(1 + v[i]^2)), va)[1]

end

@printline

@testset "zero-arrays" begin

    # Using zero-dim arrays fails on ReverseDiff & Tracker
    # Tracker.gradient(x -> x[], fill(1.0))
    # ReverseDiff.gradient(x -> x[], fill(1.0)) # is ambiguous
    if GRAD in [:Tracker, :ReverseDiff]
        @test_skip _gradient(x -> sum(@tullio y[] := log(x[i])), collect(1:3.0))[1] == 1 ./ (1:3)
    else
        @test _gradient(x -> sum(@tullio y[] := log(x[i])), collect(1:3.0))[1] == 1 ./ (1:3)
    end
    # one-element vectors are fine:
    @test _gradient(x -> sum(@tullio y[1] := log(x[i]) avx=false), collect(1:3.0))[1] == 1 ./ (1:3)  # new failure LoopVectorization v0.12.14? only on CI?
    # which is what's now used for this:
    @test _gradient(x -> (@tullio y := log(x[i])), collect(1:3.0))[1] == 1 ./ (1:3)

end
@testset "gather/scatter" begin

    inds = vcat(1:3, 1:2)
    @test _gradient(x -> sum(@tullio y[i] := x[inds[i]]), rand(3))[1] == [2,2,1]

    _gradient(x -> sum(@tullio y[inds[i]] := x[i]), rand(5))[1] == [1,1,1,1,1]
    ForwardDiff.gradient(x -> sum(@tullio y[inds[i]] := x[i]), rand(5)) == [0,0,1,1,1]
    # This difference may be another edge case like multiple maxima?

    ind2 = rand(1:10, 1024) # many repeats
    dx2 = ForwardDiff.gradient(x -> sum(@tullio y[i] := x[ind2[i]] + x[i]), rand(1024))
    @test dx2 ≈ _gradient(x -> sum(@tullio y[i] := x[ind2[i]] + x[i]), rand(1024))[1]

    ind3 = vcat(unique(rand(2:1024, 10)), 1) # many missing, no repeats, but always includes 1
    g3 = ForwardDiff.gradient(x -> sum(@tullio y[ind3[i]] := i^2 * x[i]), ones(size(ind3)))
    @test g3 ≈ _gradient(x -> sum(@tullio y[ind3[i]] := i^2 * x[i]), ones(size(ind3)))[1]
    # You get weird errors here if indices of y don't start at 1.

    # 1.6 failure on CI, with rand(1:1024, 10)
    # [1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 0.0, 64.0, 81.0, 100.0, 121.0] ≈ [1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0, 121.0]

end

@printline

@testset "shifts, etc" begin

    c1(N,K) = @tullio M[x,y,c] := N[x+i-1, y+j-1,c] * K[i,j]
    m1 = rand(10,10,2)
    k1 = rand(3,3)
    g_m = ForwardDiff.gradient(N -> sum(sin, c1(N, k1)), m1)
    g_k = ForwardDiff.gradient(K -> sum(sin, c1(m1, K)), k1)
    @test_skip g_m ≈ _gradient(N -> sum(sin, c1(N, k1)), m1)[1]  atol=0.01 # works at repl, fails in tests
    @test g_k ≈ _gradient(K -> sum(sin, c1(m1, K)), k1)[1]  atol=0.01

    c2(mat, kern) = @tullio out[x,y,n] := begin
            i = mod(x+a, axes(mat,1))
            j = mod(y+b, axes(mat,2))
            @inbounds mat[i,j,n] * abs(kern[a,b])
        end (x in axes(mat,1), y in axes(mat,2)) grad=Dual

    if Tullio._GRAD[] == :Dual
        g_m = ForwardDiff.gradient(N -> sum(sin, c2(N, k1)), m1)
        g_k = ForwardDiff.gradient(K -> sum(sin, c2(m1, K)), k1)
        @test g_m ≈ _gradient(N -> sum(sin, c2(N, k1)), m1)[1]  atol=0.01
        @test g_k ≈ _gradient(K -> sum(sin, c2(m1, K)), k1)[1]  atol=0.01
    end

end
@testset "mod, clamp, pad" begin

    fmod(x) = @tullio y[i] := x[mod(i)]  i in 1:5  avx=false # fails on 1.4, LV 0.8
    fclamp(x) = @tullio y[i] := x[clamp(i)]  i in 1:5  avx=false
    fpad(x) = @tullio y[i] := x[pad(i-2,2)]
    @test _gradient(sum∘fmod, ones(3))[1] == [2,2,1]
    @test _gradient(sum∘fclamp, ones(3))[1] == [1,1,3]
    @test _gradient(sum∘fpad, ones(3))[1] == [1,1,1]

end
@testset "@inferred" begin

    h2(x,y) = @tullio z[i] := x[i,j] + y[j,i]  # as above
    flog(x,y) = @tullio z[i] := log(x[i,j]) / y[j,i]  avx=false  # new failure LoopVectorization v0.12.14? only on CI?

    mat = rand(3,3)
    @test @inferred(h2(mat, mat)) ≈ vec(sum(mat .+ mat', dims=2))
    @test @inferred(flog(mat, mat)) isa Vector

    if GRAD == :Zygote
        @test_broken @inferred(_gradient(sum∘h2, rand(2,3), rand(3,2))) isa Tuple
        @test_broken @inferred(_gradient(sum∘flog, mat, mat)) isa Tuple
    else
        @test @inferred(_gradient(sum∘h2, rand(2,3), rand(3,2))) isa Tuple
        @test @inferred(_gradient(sum∘flog, mat, mat)) isa Tuple
    end

end

@printline

@testset "from TensorTrace" begin
    # These can all be handled using TensorOperations

    triv1(x) = @tullio A[i,j] := 2 * x[i,j]
    @test gradtest(triv1, (2,3))

    r32 = randn(3,2);
    r312 = randn(3,1,2);

    ## trace!
    tr1(x) = @tullio T[k] := 22 * x[i,i,k]
    @test gradtest(tr1, (3,3,4))

    tr2(x) = @tullio T[k] := 22 * x[i,i,k,j,j]
    @test gradtest(tr2, (3,3,4,7,7))

    ## contract! A
    con1(x) = @tullio C[i,j] := 5 * x[i,k] * r32[k,j]
    @test gradtest(con1, (2,3))

    r22 = rand(2,2);

    con3(x) = @tullio C[i,j,m,n] := x[i,j,k] * r312[k,m,n]  avx=false # I think leading size-1 dims are the problem
    @test gradtest(con3, (1,2,3))

    con4(x) = @tullio C[i,m] := x[i,kk,k] * r312[k,m,kk]  avx=false
    @test gradtest(con4, (1,2,3))

    con5(x) = @tullio C[j,i,n,m] := 44 * x[i,j,k] * r312[k,m,n]  avx=false
    @test gradtest(con5, (1,2,3))

    r392 = randn(3,9,2);
    con6(x) = @tullio C[n,i,m,j] := x[i,j,k] * r392[k,m,n]
    @test gradtest(con6, (9,2,3))

    con7(x) = @tullio C[m,n,j,i] := 44 * x[i,j,k] * r392[k,m,n]
    @test gradtest(con7, (9,2,3))

    @printline

    ## contract! B
    con8b(x) = @tullio K[i,j] := 5 * r32[i,k] * x[k,j]  avx=false
    @test gradtest(con8b, (2,3))

    con9b(x) = @tullio K[i,j,m,n] := r312[i,j,k] * x[m,k,n]
    @test gradtest(con9b, (1,2,3))

    con10b(x) = @tullio K[n,j,m,i] := r392[i,j,k] * x[m,k,n]  avx=false
    @test gradtest(con10b, (9,2,3))

    r3399 = randn(3,3,9,9);

    con13(x) = @tullio K[i,j] := r3399[s,s,j,k] * x[t,t,k,i]
    @test gradtest(con13, (3,3,9,9))

    r33 = rand(3,3);
    con14(x) = @tullio K[i,j] := r3399[a,b,j,k] * x[b,c,k,i] * r33[a,c]  avx=false
    @test gradtest(con14, (3,3,9,9))

    @printline

    ## scalar -- one with :=, one without
    sc1(x) = @tullio s = r22[b,β] * x[a,b,c] * r312[c,a,β]  avx=false verbose=true
    @test gradtest(sc1, (1,2,3)) # UndefVarError: ####op#798_0 not defined

    @printline

    sc2(x) = @tullio s := x[γ,c] * r3399[c,γ,i,i]  avx=false verbose=true
    @test gradtest(sc2, (3,3))

end

@printline

if Tullio._GRAD[] != :Dual
#=
    @testset "products" begin

        p1(x) = @tullio (*) z = x[i]
        @test _gradient(p1, 1:4)[1] == ForwardDiff.gradient(p1, 1:4)
        @test _gradient(p1, -1:3)[1] == ForwardDiff.gradient(p1, -1:3) # one zero
        @test _gradient(p1, [1,0,2,0])[1] == ForwardDiff.gradient(p1, [1,0,2,0])

        p2(m,v) = @tullio (*) y[i] := (m[i,j] + 3*v[j])^2 # / sqrt(v[i])
        p2(m,v) = @tullio (*) y[i] := m[i,j] * v[j]
        m1 = rand(4,4) .+ 1
        v1 = rand(4) .+ 1
        dm = ForwardDiff.gradient(m -> sum(p2(m,v1)), m1)
        @test dm ≈ _gradient(sum∘p2, m1, v1)[1]
        dv = ForwardDiff.gradient(v -> sum(p2(m1,v)), v1)
        @test_broken dv ≈ _gradient(sum∘p2, m1, v1)[2]

        m1[2,3] = 0
        p3(m) = @tullio (*) y[i] := 4 * m[i,j]
        @test _gradient(sum∘p3, m1)[1] ≈ ForwardDiff.gradient(sum∘p3, m1)
        m1[3,4] = -1
        p4(m) = @tullio (*) y[i] := sin(1 + m[i,j])
        @test _gradient(sum∘p4, m1)[1] ≈ ForwardDiff.gradient(sum∘p4, m1)

    end
=#
    @testset "min/max" begin

        f1(x) = @tullio (max) z = x[i]
        f2(x) = @tullio (min) z = x[i] # avx=false

        @test _gradient(f1, 1:4)[1] == ForwardDiff.gradient(f1, 1:4)
        @test _gradient(f2, 1:4)[1] == ForwardDiff.gradient(f2, 1:4)

        @test _gradient(f1, [2,2,3,3])[1] in ([0,0,1,0], [0,0,0,1]) # changes with @avx
        ForwardDiff.gradient(f1, [2,2,3,3]) == [0,0,0,1] # different sub-gradient, OK
        @test _gradient(f2, [2,2,3,3])[1] == [1,0,0,0]

        m4 = reshape(shuffle(1:3*4*5*2), 3,4,5,2);
        m2 = reshape(shuffle(1:16), 4,4);
        v2 = shuffle(1:4)

        f3(x) = @tullio (max) y[i,k,l] := x[i,j,k,l]

        @test all(==(1), sum(_gradient(sum∘f3, m4)[1], dims=2))
        @test _gradient(sum∘f3, m4)[1] ≈ ForwardDiff.gradient(sum∘f3, m4)

        f4(x) = @tullio (min) y[j] := x[i,j,k,l]

        @test all(==(1), sum(_gradient(sum∘f4, m4)[1], dims=(1,3,4)))
        @test _gradient(sum∘f4, m4)[1] ≈ ForwardDiff.gradient(sum∘f4, m4)

        f5(x,y) = @tullio (max) z[i] := x[i,j] + 0.01*y[i]

        dm = ForwardDiff.gradient(m -> sum(f5(m,v2)), m2)
        @test dm ≈_gradient(sum∘f5, m2, v2)[1]
        dv = ForwardDiff.gradient(v -> sum(f5(m2,v)), v2)
        @test dv ≈_gradient(sum∘f5, m2, v2)[2]

        f6(x,y) = @tullio (max) z[i] := x[i,j] + 0.01*y[j] # max is now along y, not perp

        dm = ForwardDiff.gradient(m -> sum(f6(m,v2)), m2)
        @test dm ≈ _gradient(sum∘f6, m2, v2)[1]
        dv = ForwardDiff.gradient(v -> sum(f6(m2,v)), v2)
        @test dv ≈ _gradient(sum∘f6, m2, v2)[2]

        f7(x,y) = @tullio (max) z[i] := x[i,j]^2 / sqrt(y[i]) + exp(y[j])  avx=false

        dm = ForwardDiff.gradient(m -> sum(f7(m,v2)), m2)
        @test dm ≈ _gradient(sum∘f7, m2, v2)[1]  # avx: broken in tests, Julia 1.4
        dm .- _gradient(sum∘f7, m2, v2)[1]
        dv = ForwardDiff.gradient(v -> sum(f7(m2,v)), v2)
        @test dv ≈ _gradient(sum∘f7, m2, v2)[2]

        f8(x,y) = @tullio (max) z[i,l] := log(x[i,j,k,l]) / y[j]^1/3  avx=false
        f9(x,y) = @tullio (min) z[i,j] := log(x[i,j,k,l]) / y[j]^1/3  avx=false
        @tullio z89[i,j,k,l] := log(m4[i,j,k,l]) / v2[j]^1/3
        length(z89), length(unique(z89))

        dm = ForwardDiff.gradient(m -> sum(f8(m,v2)), m4)
        @test dm ≈ _gradient(sum∘f8, m4, v2)[1]  # avx: OK with 0.8, broken with 0.9
        dm .- _gradient(sum∘f8, m4, v2)[1]       # at exactly one element
        dv = ForwardDiff.gradient(v -> sum(f8(m4,v)), v2)
        @test dv ≈ _gradient(sum∘f8, m4, v2)[2]

        dm = ForwardDiff.gradient(m -> sum(f9(m,v2)), m4)
        @test dm ≈_gradient(sum∘f9, m4, v2)[1]  # avx: broken with 0.8 and 0.9
        dm .- _gradient(sum∘f9, m4, v2)[1]
        dv = ForwardDiff.gradient(v -> sum(f9(m4,v)), v2)
        @test dv ≈ _gradient(sum∘f9, m4, v2)[2]  # avx: broken with 0.8 and 0.9
        dv .- _gradient(sum∘f9, m4, v2)[2]       # but broken in different elements
        # I suspect that @avx is re-ordering loops, which makes onlyone() incorrect.

    end

    @printline

    @testset "finalisers" begin

        norm2(m) = @tullio n[i] := m[i,j]^2 |> sqrt

        gradtest(norm2, (3,4))
        mat = rand(3,3)
        @test _gradient(sum∘norm2, mat)[1] ≈ ForwardDiff.gradient(sum∘norm2, mat)
        @test gradtest(norm2, (3,4))

        layer(x) = @tullio y[i,k] := mat[i,j] * x[j,k] |> tanh  avx=false # this takes 15 mins +?
        @test gradtest(layer, (3,4))

        @printline

        lse1(mat) = @tullio lse[j] := log <| exp(mat[i,j])
        @test gradtest(lse1, (3,4))

        # relu(x) = max(x, zero(x))
        # lay2(x) = @tullio y[i,k] := mat[i,j] * x[j,k] |> relu

        @printline

        mx3(x) = @tullio (max) r[i] := x[i,j]^3 |> cbrt  avx=false # sometimes gets stuck here?
        mx3(mat) # hmmm what is this?
        _gradient(sum∘mx3, mat)[1] # zero

    end
end

@printline

if GRAD == :Zygote
    @testset "nograd keyword" begin

        f2(x,y) = @tullio out[i,j] := x[i] + y[j]  nograd=y threads=false
        @test _gradient(sum∘f2, rand(2), rand(2)) == ([2,2], nothing)

        f3(x,y,z) = @tullio out[i,j] := x[i] + y[j] * z[k]  nograd=(x,z) threads=false
        @test _gradient(sum∘f3, rand(2), rand(2), ones(2)) == (nothing, [4,4], nothing)

        f0(x,y) = @tullio out[i,j] := x[i]/y[j]  nograd=(y,x) threads=false
        @test _gradient(sum∘f0, rand(2), rand(2)) == (nothing, nothing)

    end
end

@printline


text:
#EasyPlotMPL AbstractBuilder interface implementation
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
#HACK: showable works on Figure objects (not types)
#HACK: hardcoding supported mimes (might not be true for different backends)
const SUPPORTED_MIMES = Set(["image/svg+xml", "image/png"])
const SUPPORTED_BACKENDS = Set(["tk", "gtk3", "gtk", "qt", "wx"])

#==Types
===============================================================================#
mutable struct Builder <: EasyPlot.AbstractBuilder
	guimode::Bool
	backend::Symbol
	postproc::Optional{Function}
	args::Tuple
	kwargs::Base.Iterators.Pairs
end
EasyPlot.AbstractBuilder(::DS{:PyPlot}) = Builder #Register builder

#==Constructors
===============================================================================#
MPLState(b::Builder) = MPLState(false, b.backend, b.guimode)

#==Constructor interface
===============================================================================#
#Use same builder irrespective of target application:
function EasyPlot.getbuilder(::Target{T}, ::Type{Builder}, args...;
		postproc=nothing, backend=nothing, guimode::Bool=true, kwargs...) where T
	if isnothing(backend); backend = defaults.backend; end
	return Builder(guimode, backend, postproc, args, kwargs)
end

#==I/O interface
===============================================================================#
function EasyPlot._show(io::IO, mime::MIME, opt::EasyPlot.ShowOptions, fig::PyPlot.Figure)
	if !EasyPlot.isauto(opt.dim)
		@warn("Cannot currently specify plot size. Using defaults.")
	end
	show(io, mime, fig)
end

EasyPlot._showable(mime::MIME{T}, b::Builder) where T = in(string(T), SUPPORTED_MIMES)
	#hasmethod(show, (IO, typeof(mime), PyPlot.Figure)) #Apparently not enough

#==Display interface
===============================================================================#
EasyPlot.displaygui(f::PyPlot.Figure) = f.show()

#==Plot building functions
===============================================================================#
function _build(ax, eplot::EasyPlot.Plot, theme::EasyPlot.Theme)
	ax.set_title(eplot.title)
	fold = isa(eplot.xaxis, EasyPlot.FoldedAxis) ? eplot.xaxis : nothing

	builder = WfrmBuilder(ax, theme, fold)
	for (i, wfrm) in enumerate(eplot.wfrmlist)
		EasyPlot.addwfrm(builder, wfrm, i)
	end

	#x-axis properties:
	xscale = Symbol(eplot.xaxis)
	xmin = eplot.xext.min; xmax = eplot.xext.max

	#y-axis properties:
	ylabel = ""
	yscale = :lin
	ymin, ymax = (NaN, NaN)
	if length(eplot.ystriplist) > 0
		strip = eplot.ystriplist[1]
		ylabel = strip.axislabel
		yscale = strip.scale
		ymin = strip.ext.min; ymax = strip.ext.max
	end

	#Apply x/y labels:
	ax.set_xlabel(eplot.xlabel)
	ax.set_ylabel(ylabel)

	#Apply x/y scales:
	ax.set_xscale(scalemap[xscale])
	ax.set_yscale(scalemap[yscale])

	#Apply x/y extents:
	(c_xmin, c_xmax) = ax.set_xlim() #Read in current limits
		isnan(xmin) && (xmin = c_xmin); isnan(xmax) && (xmax = c_xmax)
		ax.set_xlim(xmin, xmax)
	(c_ymin, c_ymax) = ax.set_ylim() #Read in current limits
		isnan(ymin) && (ymin = c_ymin); isnan(ymax) && (ymax = c_ymax)
		ax.set_ylim(ymin, ymax)

	return ax
end

function _build(fig::PyPlot.Figure, ecoll::EasyPlot.PlotCollection)
	ecoll = EasyPlot.condxfrm_multistrip(ecoll, "EasyPlotMPL") #Emulate multi-strip plots
	ncols = ecoll.ncolumns
	fig.suptitle(ecoll.title)
	nrows = div(length(ecoll.plotlist)-1, ncols)+1
	iplot = 0

	for plot in ecoll.plotlist
#		row = div(iplot, ncols) + 1
#		col = mod(iplot, ncols) + 1
		ax = fig.add_subplot(nrows, ncols, iplot+1)
		_build(ax, plot, ecoll.theme)
		if plot.legend; ax.legend(); end
		iplot += 1
	end

	fig.canvas.draw()
	return fig
end

function EasyPlot.build(b::Builder, ecoll::EasyPlot.PlotCollection)
	local fig=nothing
	origstate = _getstate()
	newstate = MPLState(b)
	try
		_applystate(newstate)
		fig = PyPlot.figure(b.args...; b.kwargs...)
		_build(fig, ecoll)
	finally
		#Do not restore guimode... PyPlot will not display properly
		origstate.guimode = newstate.guimode
		_applystate(origstate)
	end
	return fig
end

#Last line


text:
using Statistics

@testset "Config parser" begin
    cfg_filename = string(@__DIR__, "/dualcgp_test.yaml")
    game_name = "enduro"
    seed=0
    game = Game(game_name, seed)
    n_in = 3  # RGB images
    n_out = length(getMinimalActionSet(game.ale))  # One output per legal action
    rgb = get_rgb(game)
    img_size = size(rgb[1])
    close!(game)

    encoder_cfg, controller_cfg, reducer = IICGP.dualcgp_config(
        cfg_filename, game_name, seed=seed
    )

    # I/O dimensionality
    @test encoder_cfg.n_in == n_in
    @test encoder_cfg.n_out == 2
    @test controller_cfg.n_in == 50
    @test controller_cfg.n_out == n_out

    @test encoder_cfg.img_size == img_size
    @test encoder_cfg.id == controller_cfg.id
    @test encoder_cfg.n_parameters == controller_cfg.n_parameters == 1

    @test encoder_cfg.seed == controller_cfg.seed == 0
    @test encoder_cfg.d_fitness == controller_cfg.d_fitness == 1
    @test encoder_cfg.n_gen == controller_cfg.n_gen == 1
    @test encoder_cfg.log_gen == controller_cfg.log_gen == 10
    @test encoder_cfg.save_gen == controller_cfg.save_gen == 10

    @test reducer.parameters["pooling_function"] == Statistics.mean
    @test reducer.parameters["size"] == 5

    # Construction and processing test
    enco = IICGP.IPCGPInd(encoder_cfg)
    cont = CartesianGeneticProgramming.CGPInd(controller_cfg)
    features, out = IICGP.process_f(enco, reducer, cont, rgb)
    @test typeof(features) == Array{Array{Float64,2},1}
    @test length(features) == encoder_cfg.n_out
    @test typeof(out) == Array{Float64,1}
    @test length(out) == n_out
end