text:
using Test
using OceanParameterizations
using Statistics

@testset "OceanParameterizations" begin
    include("test_feature_scaling.jl")
    include("test_coarse_graining.jl")
end

ext: jl


text:
"""
      SeisPatchProcess(in,out;<keyword arguments>)
Read data from disk, split into multidimensional overlapping patches,
apply processes, merge patches with tapers, and write to disk.

Processing of patches can be done in parallel by running the code using
(for example): julia -p 2 script_name.jl

Important Notice: You must make declare global variables f and f_param on
every processor. You can do this in your main function by typing (for
example):
@everywhere global f_param = ["style"=>"mxmyhxhy", "Niter"=> 100, "alpha"=>1,"fmax"=>80]
@everywhere global f = [SeisPOCS]

f is an array of functions that have the following syntax:
d2, h2 = f(d1,h1,f_param), where
param is a dictionary (Dict) of parameters for the function.

note that param should contain parameters for the patching and unpatching operations.

to execute the code type:  julia -p 4 my_code.jl where 4 can be replaced
    with the number of processors you wish to use.
*Credits: A. Stanton*
"""
function SeisPatchProcess(in,out,param=Dict())

    patch_names = SeisPatch(in,out,param)
    a = pmap(MyProcess,patch_names)
    SeisUnPatch(patch_names,out,param)
    #for ipatch = 1 : length(patch_names)
    #	SeisRemove(patch_names[ipatch])
    #end

end

function MyProcess(filename)
    d1,h1 = SeisRead(filename)
    for ifunc = 1 : length(f)
	func = f[ifunc]
	d2,h2 = func(d1,h1,f_param)
	d1 = copy(d2)
	h1 = copy(h2)
    end
    SeisWrite(filename,d1,h1)
    return(1)
end


text:
# tools for dealing with compiler debug information

# generate a pseudo-backtrace from LLVM IR instruction debug information
#
# this works by looking up the debug information of the instruction, and inspecting the call
# sites of the containing function. if there's only one, repeat the process from that call.
# finally, the debug information is converted to a Julia stack trace.
function backtrace(inst::LLVM.Instruction, bt = StackTraces.StackFrame[])
    # look up the debug information from the current instruction
    if haskey(metadata(inst), LLVM.MD_dbg)
        loc = metadata(inst)[LLVM.MD_dbg]
        while loc !== nothing
            scope = LLVM.scope(loc)
            if scope !== nothing
                name = replace(LLVM.name(scope), r";$"=>"")
                file = LLVM.file(scope)
                path = joinpath(LLVM.directory(file), LLVM.filename(file))
                line = LLVM.line(loc)
                push!(bt, StackTraces.StackFrame(name, path, line))
            end
            loc = LLVM.inlined_at(loc)
        end
    end

    # move up the call chain
    f = LLVM.parent(LLVM.parent(inst))
    ## functions can be used as a *value* in eg. constant expressions, so filter those out
    callers = filter(val -> isa(user(val), LLVM.CallInst), collect(uses(f)))
    ## get rid of calls without debug info
    filter!(callers) do call
        md = metadata(user(call))
        haskey(md, LLVM.MD_dbg)
    end
    if !isempty(callers)
        # figure out the call sites of this instruction
        call_sites = unique(callers) do call
            # there could be multiple calls, originating from the same source location
            md = metadata(user(call))
            md[LLVM.MD_dbg]
        end

        if length(call_sites) > 1
            frame = StackTraces.StackFrame("multiple call sites", "unknown", 0)
            push!(bt, frame)
        elseif length(call_sites) == 1
            backtrace(user(first(call_sites)), bt)
        end
    end

    return bt
end


text:
using AlphaZero

gspec = Examples.games["tictactoe"]
mcts = MCTS.Env(gspec, MCTS.RolloutOracle(gspec))
computer = MctsPlayer(mcts, niters=1, timeout=1.0, τ=ConstSchedule(0.5))

# interactive!(gspec, computer, Human())
explore(computer, gspec)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Composition

Structuring of objects moving in 3D. Most important constructors
(dof are the `degrees-of-freedom`):

| Function                                    | dof | Description |
|:--------------------------------------------|:----:|:---------------------------------------------|
| [`Object3D`](@ref)`(...)`                  | 0,6 | Return a reference Object3D, if a parent is given fixed/moving w.r.t. `parent`|
| [`Modia3D.Revolute`](@ref)`(obj1,obj2;..)`  | 1 | Return a revolute joint |
| [`Modia3D.Prismatic`](@ref)`(obj1,obj2;..)` | 1 | Return a prismatic joint |

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Composition

export initialize, initAnalysis!, performAnalysis!, closeAnalysis!, visualize!, visualizeWorld!
export updatePosition!, update!, driveJoint!
export Object3D
export RotationVariables, RCardan123
export WStartVariables, WCartesian, WCardan123
export SimulationModel, Model
export printObject3DAndAllChildren, writeObject3DAndAllChildrenOnJsonFile
export Fix, FixTranslation
export Revolute, setAngle!, connect
export Prismatic, setDistance!
export FreeMotion

# export getModelResidues!, model_f1!, model_f2!
export multibodyResiduals

export distanceAndAngles, distance, planarRotationAngle
export measFrameRotation, measFramePosition, measFrameDistance
export measFrameRotVelocity, measFrameTransVelocity, measFrameDistVelocity
export measFrameRotAcceleration, measFrameTransAcceleration, measFrameDistAcceleration
export applyFrameTorque!, applyFrameForce!, applyFrameForceTorque!
export applyFrameTorquePair!, applyFrameForcePair!, applyFrameForceTorquePair!

export deleteMaterialLastContactDictContactEnd

export NoGravityField, UniformGravityField, PointGravityField, gravityAcceleration
export G, EarthMass, EarthRadius
export Scene, SceneOptions
export upwardsDirection, cameraPosition
export animationData, animationStep

export DLR_Visualization, DLR_Visualization_renderer

export ContactDetectionMPR, ContactDetectionMPR_handler
export initializeContactDetection!, selectContactPairsNoEvent!, selectContactPairsWithEvent!, getDistances!, setComputationFlag, closeContactDetection!

export responseCalculation, print_ModelVariables
export contactStart, contactEnd, regularize, resultantCoefficientOfRestitution, resultantDampingCoefficient
export ElasticContactPairResponseMaterial, elasticContactPairCoefficients, normalRelativeVelocityAtContact

export isFree, isNotFree, changeJointFromFreeMotionToFix!, changeJointFromFixToFreeMotion!

export fullName, instanceName

# export updateSimulationModel!

export InteractionBehavior, Gripper, Movable, Lockable, NoInteraction

export rot123fromR

export supportPoint, boundingBox!

@enum InteractionBehavior Gripper Movable Lockable NoInteraction

using StaticArrays
using OrderedCollections
using LinearAlgebra
using Unitful

import Base
import Modia3D
import Modia3D.Basics
import Modia3D.Frames
import Modia3D.Shapes
import JSON
import Printf
import ModiaLang
import TimerOutputs

include(joinpath("joints", "object3DMotion.jl"))
include("object3D.jl")
include("supportPoints.jl")
include("superObjects.jl")

# Joints must be included before scene, because joint memory is stored in scene
include(joinpath("joints", "FreeMotion.jl"))
include(joinpath("joints", "Fix.jl"))
include(joinpath("joints", "Revolute.jl"))
include(joinpath("joints", "Prismatic.jl"))

include("scene.jl") # must be included after superObjects.jl
include(joinpath("joints", "joints.jl"))
include(joinpath("joints", "changeJointState.jl"))

include("sensors.jl")
include("frameMeasurements.jl")
include("frameForceTorque.jl")
include("assignObjects.jl")
include(joinpath("responseCalculation", "othersCollisionResponse.jl"))
include(joinpath("responseCalculation", "elasticCollisionResponse.jl"))
include("massPropertiesComputation.jl")
include("printObject3D.jl")
include("handler.jl")
include("responseCalculation.jl")
include("dynamicCollision.jl")
include("dynamics.jl")

end


text:
"""
    get_horizon(c::Contract) :: Union{Missing, Date}

Returns the furthest date of choice or maturity for a given contract.

Horizon is not defined for spot positions. In this case,
this function returns `missing`.
"""
get_horizon(c::Contract) :: Union{Missing, Date} = error("Not implemented.")
get_horizon(::Union{Worthless, Unit}) :: Union{Missing, Date} = missing
get_horizon(c::Union{Give, Scale}) :: Union{Missing, Date} = get_horizon(c.c)
get_horizon(c::FixedIncomeEvent) :: Date = get_horizon(c.c)

function get_horizon(c::FixedIncomeContract) :: Union{Missing, Date}
    if isempty(c.events)
        return missing
    else
        result = get_horizon(c.events[1])
        if length(c.events) > 1
            for i in 2:length(c.events)
                result = max(result, get_horizon(c.events[i]))
            end
        end

        return result
    end
end

function get_horizon(c::SwapContract) :: Union{Missing, Date}
    return max(get_horizon(c.asset), get_horizon(c.liability))
end

function get_horizon(c::Union{Both, Either}) :: Union{Missing, Date}
    h1 = get_horizon(c.c1)
    h2 = get_horizon(c.c2)
    if !ismissing(h1) && !ismissing(h2)
        return max(h1, h2)
    elseif ismissing(h1) && !ismissing(h2)
        return h2
    elseif !ismissing(h1) && ismissing(h2)
        return h1
    elseif ismissing(h1) && ismissing(h2)
        return missing
    else
        error("Unexpected state.")
    end
end

function get_horizon(c::WhenAt) :: Union{Missing, Date}
    if c.expires_at_maturity
        return c.maturity
    else
        inner_horizon = get_horizon(c.c)
        if ismissing(inner_horizon)
            return c.maturity
        else
            return max(c.maturity, inner_horizon)
        end
    end
end

"""
    is_expired(c::Contract, pricing_date::Date) :: Union{Missing, Bool}

Returns wether this contract has expired.

# Possible Results

* before contract horizon: `false`

* at contract horizon: `false`

* after contract horizon: `true`

* horizon is not defined for this contract: `false`, meaning it never expires.

See also: [`is_expired_or_expires_today`](@ref).
"""
function is_expired(c::Contract, pricing_date::Date) :: Bool
    h = get_horizon(c)
    if ismissing(h)
        return false
    else
        return pricing_date > h
    end
end

"""
    is_expired_or_expires_today(c::Contract, pricing_date::Date) :: Union{Missing, Bool}

Similar to [`is_expired`](@ref), but with different behavior in the case
the contract expires today.

# Possible Results

* before contract horizon: `false`

* at contract horizon: `true`

* after contract horizon: `true`

* horizon is not defined for this contract: `false`, meaning it never expires.
"""
function is_expired_or_expires_today(c::Contract, pricing_date::Date) :: Bool
    h = get_horizon(c)
    if ismissing(h)
        return false
    else
        return pricing_date >= h
    end
end


text:
"""
    Fermi.Geometry

Module handling molecule and atoms data

# Objects

    Atom      Object storing information about an atom
    Molecule  Object storing information about a molecule (group of atoms)
"""
module Geometry

export Molecule
export Atom

using Fermi
using Fermi.Options
using Fermi.PhysicalConstants: atomic_number, bohr_to_angstrom

using LinearAlgebra
using Formatting

import Base: show

"""
    Fermi.Atom

Object storing information about an atom.

# Fields

| Name  | Description |
|:------|:-----------------------------------------------------------|
|`AtomicSymbol` |  Atomic symbol          |
|`Z`            |  Atomic number          |
|`xyz`          |  xyz array in Angstrom  |
"""
struct Atom
    AtomicSymbol::String
    Z::Int
    xyz::Tuple{Float64,Float64,Float64}
end

"""
    Fermi.Molecule

Object storing information about a molecule (group of atoms).

# Fields

    atoms         Array with Fermi.Atom objects
    charge        Charge of the molecule
    multiplicity  Multiplicity ``(2Mₛ + 1)``
    Nα            Number of Alpha electrons
    Nβ            Number of Beta electrons
    Vnuc          Nuclear repulsion energy

# Examples:

A Molecule object can be created by providing the desired keyworded arguments:

    molstring:      A string representing the XYZ of the molecule
    unit:           Unit used for the distance between atoms (Bohr or Angstrom)
    charge:         Molecular charge
    multiplicity    Molecular multiplicity

Any argument not given explicitly will be read from the Options.
```
julia> Fermi.Geometry.Molecule()
Molecule:

O    1.209153654800    1.766411818900   -0.017161397200
H    2.198480007500    1.797710062700    0.012116171900
H    0.919788188200    2.458018557000    0.629793883200

Charge: 0   Multiplicity: 1
Nuclear repulsion:    8.8880641737

julia> Fermi.Geometry.Molecule(charge=2, multiplicity=3)
Molecule:

O    1.209153654800    1.766411818900   -0.017161397200
H    2.198480007500    1.797710062700    0.012116171900
H    0.919788188200    2.458018557000    0.629793883200

Charge: 2   Multiplicity: 3
Nuclear repulsion:    8.8880641737
```
"""
struct Molecule
    atoms::Tuple
    charge::Int
    multiplicity::Int
    Nα::Int
    Nβ::Int
    Vnuc::Float64
end

function Molecule(;
    molstring = Options.get("molstring"),
    unit = Options.get("unit"),
    charge = Options.get("charge"),
    multiplicity = Options.get("multiplicity")
    )

    atoms = get_atoms(molstring, unit=unit)
    Molecule(atoms, charge, multiplicity)
end

function Molecule(atoms::Array{Atom,1}, charge::Int, multiplicity::Int)

    # Compute Nuclear repulsion
    Vnuc = 0.0
    for i in eachindex(atoms)
        for j in 1:(i-1)
            Vnuc += nuclear_repulsion(atoms[i], atoms[j])
        end
    end

    # Compute number of electrons
    nelec = -charge
    for i in eachindex(atoms)
        nelec += atoms[i].Z
    end

    # If the number of electrons turns out negative returns an error
    if nelec ≤ 0
        throw(FermiException("Invalid charge ($charge) for given molecule"))
    end

    # Mult = 2Ms + 1 thus the number of unpaired electrons (taken as α) is Mult-1 = 2Ms
    αexcess = multiplicity-1

    # The number of β electrons must be equal the number of doubly occupied orbitals (Ndo)
    # Ndo = (nelec - n_of_unpaired_elec)/2 this must be an integer
    if isodd(nelec) != isodd(αexcess)
        throw(FermiException("Incompatible charge $(charge) and multiplicity $(multiplicity)"))
    end

    Nβ = (nelec - αexcess)/2
    Nα = nelec - Nβ
    out =  Molecule(Tuple(atoms), charge, multiplicity, Nα, Nβ, Vnuc)
    return out
end

"""
    Fermi.get_atoms(molstring::String; unit::String="angstrom")

From a XYZ string, produces an array of Atom objects.
"""
function get_atoms(molstring::String; unit::String="angstrom")

    # Get a list of Atom objects from String
    if unit == "bohr"
        conv = bohr_to_angstrom
    elseif unit == "angstrom"
        conv = 1.0
    else
        throw(FermiException("unknown unit in molecule construction: $unit"))
    end

    atoms = Atom[]
    for line in split(strip(molstring), "\n")
        m = split(line)
        if length(m) != 4
            throw(FermiException("4 columns expected on XYZ string. Got $(length(m))"))
        end

        if m !== nothing
            # Convert SubString to String
            AtomicSymbol = String(m[1])
            Z = atomic_number(AtomicSymbol)

            # Convert Substring to String and then String to Float64 Array
            xyz = [0.0, 0.0, 0.0]
            try
                xyz .= parse.(Float64, m[2:4]).*conv
            catch ArgumentError
                throw(FermiException("Failed to process XYZ string: $(m[2:4])"))
            end

            push!(atoms, Atom(AtomicSymbol, Z, Tuple(xyz)))
        end
    end
    return atoms
end

"""
    Fermi.Geometry.nuclear_repulsion(A::Atom, B::Atom)

Returns the repulsion energy between atoms A and B.
"""
function nuclear_repulsion(A::Atom, B::Atom)
    return (A.Z*B.Z)/(√((A.xyz.-B.xyz)⋅(A.xyz.-B.xyz))/bohr_to_angstrom)
end

"""
    Fermi.get_xyz(M::Molecule)

Returns a XYZ string in angstrom for the given Molecule.
"""
function get_xyz(M::Molecule)

    molstring = ""

    for i in eachindex(M.atoms)

        A = M.atoms[i]
        molstring *= format("{}   {: 15.12f}   {: 15.12f}   {: 15.12f}\n", A.AtomicSymbol, A.xyz...)
    end
    return molstring
end

"""
    Fermi.string_repr(M::Molecule)

Returns a nicely formatted string with all the molecule's information
"""
function string_repr(M::Molecule)
    out = ""
    out = out*format("Molecule:\n\n")
    out = out*format(Fermi.Geometry.get_xyz(M))
    out = out*format("\n")
    out = out*format("\nCharge: {}   ", M.charge)
    out = out*format("Multiplicity: {}   \n", M.multiplicity)
    out = out*format("Nuclear repulsion: {:15.10f}", M.Vnuc)
    return out
end

"""
    Fermi.string_repr(A::Atom)

Returns a nicely formatted string for Atom object.
"""
function string_repr(A::Atom)
    out = "Symbol:           $(A.AtomicSymbol)\n"
    out = out*"Atomic Number:    $(A.Z)\n"
    out = out*"Position: $(format("{:2.5f}  {:2.5f}  {:2.5f}", A.xyz...))"
end

# Pretty printing
function show(io::IO, ::MIME"text/plain", X::T) where T<:Union{Molecule, Atom}
    print(io, string_repr(X))
end
end #Module


text:
module NumericallyIntegrateArrays
using OffsetArrays, LinearAlgebra

# package code goes here
export simps,trapz
# export sphericalIntegral,cylindricalIntegral

############################################################################################
## helper functions
############################################################################################

function check_last_value(a::AbstractArray,val=2π;var="ϕ",switch="ϕ2π")
	if a[end]≈val
		warn("$(var) appears to have $(val) as the last element, "*
			"but the keyword argument $(switch) is set to false.\n"*
			"To fix this, add $(switch)=true to your function call")
	end
end

scalar_to_0dim_array(z::Number)::Array{<:Number,0} = setindex!(Array{eltype(z)}(),z)
scalar_to_0dim_array(T::Type{<:Number},z::Number)::Array{<:Number,0} = setindex!(Array{T}(),z)

##############################################################################################

function leading_trailing_indices(y::AbstractArray{<:Number},dim::Integer)
	ax = axes(y)
	leading_axes = CartesianIndices(ax[1:dim-1])
	trailing_axes = CartesianIndices(ax[dim+1:end])
	inds_type = NTuple{ndims(y)-1,<:AbstractUnitRange}
	inds = (leading_axes.indices...,trailing_axes.indices...) :: inds_type
	return leading_axes,trailing_axes,inds
end

function trailing_indices(y::AbstractArray{<:Number})
	inds = Base.tail(axes(y))
	trailing_axes = CartesianIndices(inds)
	return trailing_axes,inds
end

include("trapz.jl")
include("simps.jl")

#####################################################################################

end # module


text:
using Images, ImageView

function main()
    img = Gray.(load("cat.png"))
    template = Gray.(load("template.png"))
    imshow(img)
    println(typeof(img))
    rows, cols = size(img)

    if check_size(img, template)
        match_template(img, template, "SSD")
    end

end

function check_size(img, template)
    r_img, c_img = size(img)
    r_tmp, c_tmp = size(template)
    if r_tmp > r_img || c_tmp > c_img
        return false
    end
    return true
end

function ssd(img, template)
    r_img, c_img = size(img)
    r_tmp, c_tmp = size(template)
    ssd_map = zeros((r_img - r_tmp + 1, c_img - c_tmp + 1))
    for y in 1:r_img - r_tmp + 1
        for x in 1:c_img - c_tmp + 1
            ssd_value = 0
            for yt in 1:c_tmp
                for xt in 1:r_tmp
                    ssd_value += (img[y+yt-1, x+xt-1] - template[yt, xt])^2
                end
            end
            ssd_map[y, x] = ssd_value
        end
    end
    return ssd_map
end

function match_template(img, template, mode)
    println("Template matching")
    ssd_img = ssd(img, template)
    imshow(ssd_img)
end

main()


text:
module IZCondMarkram
using ...AgentParts.IZNeuron
using ...AgentParts.ConductanceInjection
using ...AgentParts.MarkramTransmitter
using ...Types
using ...Network
using ...Signals

struct IZCondMrkrmStates{T <: AbstractFloat} # <: AgentStates{T}
    iz::IZStates{T}
    cond::ConductanceStates{T}
    mrkrm::MarkramStates{T}
end

struct IZCondMrkrmParams{T <: AbstractFloat}
    iz::IZParams
    cond::ConductanceParams
    mrkrm::Markramparams
    delta_v::AbstractFloat
end

struct IZCondMrkrmAgent{T <: AbstractFloat} <: Agent
    states::IZCondMrkrmStates{T}
    params::IZCondMrkrmParams{T}
    # donors_t_delta_v::Vector{DonorTimedDeltaV}
    # donors_t_delta_g::Vector{DonorTimedDeltaCond}
    acceptors_t_delta_v::Vector{Address}
    acceptors_mrkrm::Vector{Address}
    donors_t_delta_v::Vector{Address}
    stack_t_delta_v::Vector{TimedDeltaV}
    donors_t_exct_delta_g::Vector{Address}
    stack_t_exct_delta_g::Vector{TimedExctDeltaCond}
    donors_t_inhbt_delta_g::Vector{Address}
    stack_t_inhbt_delta_g::Vector{TimedInhbtDeltaCond}
end

# accept(acceptor) return [] of msgs
function act(agent::IZCondMrkrmAgent, t, dt, put_task)

    function cond_update(cond_states::ConductanceStates)
        agent.states.cond = cond_states
    end

    function inject_fn()
        agent.acceptors_t_delta_v.take(t)
        agent.acceptors_t_exct_delta_g.take(t)
        agent.acceptors_t_inhbt_delta_g.take(t)
        delta_v_signals = vcat(map(accptr -> take_due_signals(t, accptr),
                                   agent.acceptors_t_delta_v))
        delta_exct_signals = vcat(map(accptr -> take_due_signals(t, accptr),
                                      agent.acceptors_t_exct_delta_g))
        delta_inhbt_signals = vcat(map(accptr -> take_due_signals(t, accptr),
                                      agent.acceptors_t_inhbt_delta_g))

        if ! isnothing(agent.states.iz.idle_end)
            delta_v_signals = filter(s -> s.t >= agent.states.iz.idle_end, delta_v_signals)
            delta_exct_signals = filter(s -> s.t >= agent.states.iz.idle_end, delta_exct_signals)
            delta_inhbt_signals = filter(s -> s.t >= agent.states.iz.idle_end, delta_inhbt_signals)
        end

        delta_v = reduce((acc, x) -> acc + x.delta_v, delta_v_signals, 0.)
        delta_exct = reduce((acc, x) -> acc + x.delta_cond, delta_exct_signals, 0.)
        delta_inhbt = reduce((acc, x) -> acc + x.delta_cond, delta_inhbt_signals, 0.)

        i_syn = gen_syn_current(
            dt, delta_exct, delta_inhbt, agent.states.cond, agent.params.cond, agent.state.iz.v, cond_update
        )
        return (i_syn, delta_v)
    end

    function iz_update(iz_states::IZStates)
        agent.states.iz = iz_states
    end

    function mrkrm_update(mrkrm_states::MarkramStates)
        agent.states.mrkrm = mrkrm_states
    end

    function put_mrkrm_signal(signal::TimedMarkram)
        foreach(dnr -> dnr.put(signal), agent.acceptors_mrkrm)
    end

    function fire_fn(t, dt)
        fire(t, dt, agent.states.mrkrm, agent.params.mrkrm, mrkrm_update, put_mrkrm_signal)
        foreach(dnr -> put_task(t + dt, dnr.address), agent.donors_mrkrm)
        for dnr in agent.donors_simple
            dnr.put(TimedDelta(t, agent.params.delta_v))
            put_task(t + dt, dnr.address)
        end
    end

    evolve(t,
           dt,
           agent.states.iz,
           agent.params.iz,
           inject_fn,
           iz_update,
           fire_fn,
           t -> put_task(t, agent.Address))

end

end # Module end


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

mutable struct IoK8sApiAutoscalingV2beta1MetricSpec <: SwaggerModel
    external::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ExternalMetricSource } # spec name: external
    object::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ObjectMetricSource } # spec name: object
    pods::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1PodsMetricSource } # spec name: pods
    resource::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ResourceMetricSource } # spec name: resource
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAutoscalingV2beta1MetricSpec(;external=nothing, object=nothing, pods=nothing, resource=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1MetricSpec, Symbol("external"), external)
        setfield!(o, Symbol("external"), external)
        validate_property(IoK8sApiAutoscalingV2beta1MetricSpec, Symbol("object"), object)
        setfield!(o, Symbol("object"), object)
        validate_property(IoK8sApiAutoscalingV2beta1MetricSpec, Symbol("pods"), pods)
        setfield!(o, Symbol("pods"), pods)
        validate_property(IoK8sApiAutoscalingV2beta1MetricSpec, Symbol("resource"), resource)
        setfield!(o, Symbol("resource"), resource)
        validate_property(IoK8sApiAutoscalingV2beta1MetricSpec, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAutoscalingV2beta1MetricSpec

const _property_map_IoK8sApiAutoscalingV2beta1MetricSpec = Dict{Symbol,Symbol}(Symbol("external")=>Symbol("external"), Symbol("object")=>Symbol("object"), Symbol("pods")=>Symbol("pods"), Symbol("resource")=>Symbol("resource"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAutoscalingV2beta1MetricSpec = Dict{Symbol,String}(Symbol("external")=>"IoK8sApiAutoscalingV2beta1ExternalMetricSource", Symbol("object")=>"IoK8sApiAutoscalingV2beta1ObjectMetricSource", Symbol("pods")=>"IoK8sApiAutoscalingV2beta1PodsMetricSource", Symbol("resource")=>"IoK8sApiAutoscalingV2beta1ResourceMetricSource", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1MetricSpec }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1MetricSpec))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1MetricSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1MetricSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1MetricSpec }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1MetricSpec[property_name]

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

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1MetricSpec }, name::Symbol, val)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

@recipe function f(solution::EstimationSolution{<:PointSet}; variables=nothing)
  # retrieve domain geometry
  sdomain = domain(solution)
  dim = ndims(sdomain)

  # valid variables
  validvars = sort(collect(keys(solution.mean)))

  # plot all variables by default
  variables == nothing && (variables = validvars)
  @assert variables ⊆ validvars "invalid variable name"

  # plot layout: mean and variance for each variable
  layout := (length(variables), 2)

  for (i,var) in enumerate(variables)
    X  = hcat([coordinates(sdomain, i) for i in 1:npoints(sdomain)]...)
    μ  = solution.mean[var]
    σ² = solution.variance[var]

    if dim == 1
      @series begin
        subplot := 2i - 1
        seriestype := :scatter
        legend := false
        title --> string(var, " mean")
        X[1,:], μ
      end
      @series begin
        subplot := 2i
        seriestype := :scatter
        legend := false
        title --> string(var, " variance")
        X[1,:], σ²
      end
    elseif dim == 2
      @series begin
        subplot := 2i - 1
        seriestype := :scatter
        color := μ
        markercolor --> :bluesreds
        legend --> false
        title --> string(var, " mean")
        X[1,:], X[2,:]
      end
      @series begin
        subplot := 2i
        seriestype := :scatter
        color := σ²
        markercolor --> :bluesreds
        legend --> false
        title --> string(var, " variance")
        X[1,:], X[2,:]
      end
    elseif dim == 3
      @series begin
        subplot := 2i - 1
        seriestype := :scatter
        color := μ
        markercolor --> :bluesreds
        legend --> false
        title --> string(var, " mean")
        X[1,:], X[2,:], X[3,:]
      end
      @series begin
        subplot := 2i
        seriestype := :scatter
        color := σ²
        markercolor --> :bluesreds
        legend --> false
        title --> string(var, " variance")
        X[1,:], X[2,:], X[3,:]
      end
    else
      error("cannot plot solution in more than 3 dimensions")
    end
  end
end


text:
##################################################################
#                                             twalkexample.jl
#
#  Example for the twalk implementation in Julia,
#  Created by J Andres Christen, jac at cimat.mx .
#  Julia implementation by Nicolás Kuschinski
#
#  See http://www.cimat.mx/~jac/twalk/ for more details.
#  and for the Python, R, Matlab and C implementations.
#
##################################################################

### Differences between jtwalk and pytwalk versions:
## The Julia version does not perform output analysis
## Does not estimate the remaining time
## and does not include specialized plotting functions.

## Plotting is easy in julia! The following example
## includes sample plots using the Plots.jl plotting library

using JTwalk

####################################################################
######### Product of Independent Exponentilas #########

###   In the jtwalk object, we define the objective function with an
###   energy function U.
###   The support is defined in a separate function, and is defined in
###   jtwalk as Supp.

###   NOTE: The objective function is defined as -log of the target
###   density function

lambdas = [ 1., 2., 3., 4., 5.]

function ExpU(x)
    # -log of a product of exponentials
    return sum(x .* lambdas)
end

function ExpSupp(x)
	return all(0. .< x)
end

###   The following initializes the TWalk object
###   The dimension of the parameter space is n
Exp = jtwalk( n=5, U=ExpU, Supp=ExpSupp)

#### This runs the T-Walk

Run!(Exp, T=50000, x0=30*ones(5), xp0=40*ones(5))

##### The trajectory of the chain is stored into Exp.Output & Exp.Output_u
##### These are (n+1) by T matrices, with one row per iteration.
##### Each row has the output followed by the energy function evaluated
##### at the corresponding location. Both matrices are valid MCMC chains
##### but THEY ARE NOT INDEPENDENT.

### The following shows how to use the Plots.jl plotting library
### to make relevant plots. In julia even plotting complex thinning can
### be done with simple one line instructions.

using Plots

### The following produces a trace plot:

plot(Exp.Output[:,end])

### The following shows a scatter plot of the first and third variables with
### thinning of 1000 iteration burn-in and 250 iteration autocorrelation time:

scatter(Exp.Output[1000:250:end,1], Exp.Output[1000:250:end,3])

### The following shows a histogram of the 4th variable with the same thinning:

histogram(Exp.Output[1000:250:end,4])

max_stars_repo_name: tuerda/JTwalk

max_issues_repo_name: tuerda/JTwalk

max_forks_repo_name: tuerda/JTwalk


text:
module style

import Parameters: @with_kw

import Term: unspace_commas, NAMED_MODES, CODES, ANSICode

import ..markup: MarkupTag, extract_markup, has_markup, clean_nested_tags
import ..color:
    AbstractColor, NamedColor, is_color, is_background, get_color, is_hex_color, hex2rgb

export MarkupStyle, extract_style

is_mode(string::AbstractString) = string ∈ NAMED_MODES

# ---------------------------------------------------------------------------- #
#                                    STYLEX                                    #
# ---------------------------------------------------------------------------- #
"""
    MarkupStyle

Holds information about the style specification set out by a `MarkupTag`.
"""
@with_kw mutable struct MarkupStyle
    default::Bool = false
    bold::Bool = false
    dim::Bool = false
    italic::Bool = false
    underline::Bool = false
    blink::Bool = false
    inverse::Bool = false
    hidden::Bool = false
    striked::Bool = false

    color::Union{Nothing,AbstractColor} = nothing
    background::Union{Nothing,AbstractColor} = nothing

    tag::MarkupTag
end

"""
    MarkupStyle(tag::MarkupTag)

Builds a MarkupStyle definition from a MarkupTag.
"""
function MarkupStyle(tag::MarkupTag)
    codes = split(unspace_commas(tag.markup))

    style = MarkupStyle(; tag = tag)
    for code in codes
        if is_mode(code)
            setproperty!(style, Symbol(code), true)
        elseif is_color(code)
            setproperty!(style, :color, get_color(code))
        elseif is_background(code)
            setproperty!(style, :background, get_color(code; bg = true))
        elseif code != "nothing"
            @debug "Code type not recognized: $code" tag tag.markup typeof(code)
        end
    end
    return style
end

function toDict(style::MarkupStyle)
    return Dict(
        fieldnames(typeof(style)) .=> getfield.(Ref(style), fieldnames(typeof(style)))
    )
end

# ------------------------------ extract style ------------------------------ #
"""
    extract_style(text::AbstractString)

Extract style information from a `text`.
"""
function extract_style(text::AbstractString)
    tags = extract_markup(text)
    styles = [MarkupStyle(tag) for tag in tags]
    return styles
end

# -------------------------------- apply style ------------------------------- #

"""
    get_style_codes(style::MarkupStyle)

Get `ANSICode`s corresponding to a `MarkupStyle`.
"""
function get_style_codes(style::MarkupStyle)
    # start applying styles
    style_init, style_finish = "", ""
    for (attr, value) in toDict(style)
        # BACKGROUND
        if attr == :background
            if !isnothing(value)
                code = ANSICode(value; bg = true)
            else
                code = nothing
            end

            # COLOR
        elseif attr == :color
            if !isnothing(value)
                code = ANSICode(value; bg = false)
            else
                code = nothing
            end

            # MODES
        elseif attr != :tag && value == true
            code = CODES[attr]

        else
            if value != false && attr != :tag
                @debug "Attr/value not recognized or not set" attr value
            end
            continue
        end

        if !isnothing(code)
            style_init *= code.open
            style_finish *= code.close
        end
    end

    return style_init, style_finish
end

"""
    apply_style(text::AbstractString, style::MarkupStyle)::AbstractString

Apply a `style` to a `text`.
"""
function apply_style(text::AbstractString, style::MarkupStyle)::AbstractString
    style_init, style_finish = get_style_codes(style)

    return text = style_init * text * style_finish
end

"""
    apply_style(text::AbstractString, tag::MarkupTag)::AbstractString

Appliy the style of a markup tag and it's nested tags
"""
function apply_style(
    text::AbstractString, tag::MarkupTag; isinner::Bool = false
)::AbstractString
    style = MarkupStyle(tag)
    style_init, _ = get_style_codes(style)

    # if no inner tags, just style the text
    if length(tag.inner_tags) == 0
        return apply_style(tag.text, style)
    end

    # apply inner tags
    for inner in tag.inner_tags
        inner_text = apply_style(inner.text, inner; isinner = true)

        text = replace(
            isinner ? text : tag.text,
            "[$(inner.open.markup)]$(inner.text)[$(inner.close.markup)]" => "$(inner_text)$(style_init)",
        )
        # @info "\e[34minner style" inner.open.markup style_init text
    end

    # apply outer tag
    text = apply_style(text, style)
    # @info "\e[33mdone a style" tag.open.markup style_init text

    return text
end

"""
    apply_style(text::AbstractString)

Extract and apply all markup style in a string.
"""
function apply_style(text::AbstractString;)::AbstractString
    # @info "Applying style to " text
    text = clean_nested_tags(text)

    while has_markup(text)
        tag = extract_markup(text; firstonly = true)
        # @info "tag" tag tag.markup tag.open.start tag.close.stop

        pre = text[1:(tag.open.start - 1)]
        post = text[(tag.close.stop + 1):end]

        text = pre * apply_style(text, tag) * post
        # @info "     \e[31mdoing a style: " tag.markup pre post tag.open.start tag.close.stop

    end
    # @info "  \e[32mAfter styling" text has_markup(text)
    return text
end

end


text:
"""
    Vocabulary{T}(voc::Vector{T}, unk::T) where T

struct for holding the vocabulary list to encode/decode input tokens.
"""
struct Vocabulary{T}
    siz::Int
    list::Vector{T}
    unk::T
    unki::Int
    function Vocabulary(voc::Vector{T}, unk::T) where T
        if !(unk ∈ voc)
            pushfirst!(voc, unk)
            unki = 1
        else
            unki = findfirst(isequal(unk), voc)
        end

        new{T}(length(voc), voc, unk, unki)
    end
end

Base.length(vocab::Vocabulary) = vocab.siz
Base.getindex(vocab::Vocabulary, is) = decode(vocab, is)
Base.getindex(vocab::Vocabulary, i, is...) = (decode(vocab, i), map(i->decode(vocab, i), is)...)

"""
    encode(vocab::Vocabulary, x)

encode the given data to the index encoding.
"""
encode(vocab::Vocabulary{T}, i::Union{T,W}) where {T,W} = something(findfirst(isequal(i), vocab.list), vocab.unki)

encode(vocab::Vocabulary{T}, xs::Container{<:Union{T,W}}) where {T,W} = indices = map(x->encode(vocab, x), xs)

function encode(vocab::Vocabulary{T}, xs::Container{<:Container{<:Union{T,W}}}) where {T,W}
    lens = map(length, xs)
    indices = fill(vocab.unki, maximum(lens), length(xs))
    for (i, x) ∈ enumerate(xs)
        for (j, xi) ∈ enumerate(x)
            @inbounds indices[j, i] = encode(vocab, xi)
        end
    end
    indices
end

"""
    (vocab::Vocabulary)(x)

encode the given data to the index encoding.
"""
(vocab::Vocabulary)(x) = encode(vocab, x)
(vocab::Vocabulary)(x, xs...) = (encode(vocab, x), map(x->encode(vocab, x), xs)...)

decode(vocab::Vocabulary{T}, i::Int) where T = 0 <= i <= length(vocab) ? vocab.list[i] : vocab.unk

decode(vocab::Vocabulary{T}, is::Container{Int}) where T = map(i->decode(vocab, i), is)

function decode(vocab::Vocabulary{T}, is::Container{<:Container{Int}}) where T
    tokens = Vector{Vector{T}}(undef, length(is))
    for (idx, i) ∈ enumerate(is)
        token = decode(vocab, i)
        tokens[idx] = token
    end
    tokens
end

function decode(vocab::Vocabulary{T}, is::AbstractMatrix{Int}) where T
    ilen, olen = size(is)
    tokens = Vector{Vector{T}}(undef, olen)
    for idx ∈ 1:olen
        token = Vector{T}(undef, ilen)
        for idy ∈ 1:ilen
            token[idy] = decode(vocab, is[idy, idx])
        end
        tokens[idx] = token
    end
    tokens
end

Base.show(io::IO, v::Vocabulary) = print(io, "Vocabulary($(v.siz), unk=$(v.unk))")

Flux.onehot(v::Vocabulary{T}, x::Container{<:Union{T, Container{T}}}) where T = onehot(v, v(x))
Flux.onehot(v::Vocabulary, x::AbstractArray{Int}) = onehotarray(length(v), x)

Flux.onecold(v::Vocabulary{T}, p) where T = decode(v, _onecold(p))
function _onecold(p)
    p = collect(p)
    y = Array{Int}(undef, Base.tail(size(p)))
    for i = 1:length(y)
        ind = Tuple(CartesianIndices(y)[i])
        y[i] = argmax(@view(p[:, ind...]))
    end
    y
end

"""
    getmask(ls::Container{<:Container})

get the mask for batched data.
"""
function getmask(ls::Container{<:Container})
    lens = map(length, ls)
    m = zeros(Float32, maximum(lens), length(lens))

    for (i, l) ∈ enumerate(ls)
        selectdim(selectdim(m, 2, i), 1, 1:length(l)) .= 1
    end
    reshape(m, (1, size(m)...))
end

"""
    getmask(m1::A, m2::A) where A <: Abstract3DTensor

get the mask for the covariance matrix.
"""
getmask(m1::A, m2::A) where A <: Abstract3DTensor = permutedims(m1, [2,1,3]) .* m2


text:
using TightlyBound
using Test
using Suppressor

@testset "testing makescf" begin

    @suppress begin
        types=["Li"]
        pos=zeros((1,3))
        A=[ [10.0 0 0]; [0 10.0 0]; [ 0 0 10.0]]
        c=makecrys(A, pos, types, units="Bohr")

        tmpdir, prefix, inputfile_str = TightlyBound.QE.makeSCF(c)

        @test occursin("Li", inputfile_str)

    end
end


text:
module QNaN

export qnan, isqnan, isqnanPos, isqnanNeg, isjnan

#=
  A float64 quiet NaN is represented with these 2^52-2 UInt64 hexadecimal patterns:
  (positive) UnitRange(0x7ff8000000000000,0x7fffffffffffffff) has 2^51-1 realizations
  (negative) UnitRange(0xfff8000000000000,0xffffffffffffffff) has 2^51-1 realizations

  A float32 quiet NaN is represented with these 2^23-2 UInt32 hexadecimal patterns:
  (positive) UnitRange(0x7fc00000,0x7fffffff) has 2^22-1 realizations
  (negative) UnitRange(0xffc00000,0xffffffff) has 2^22-1 realizations

  A float16 quiet NaN is represented with these 2^10-2 UInt16 hexadecimal patterns:
  (positive) UnitRange(0x7e00,0x7fff) has 2^9-1 realizations
  (negative) UnitRange(0xfe00,0xffff) has 2^9-1 realizations

  Julia appears to assign as NaNs (generally denoting an Indeterminacy)
  0xfff8000000000000, 0xffc00000, 0xfe00, 0x7ff8000000000000, 0x7fc00000, 0x7e00
  Preferring not to give that more than one use, they excluded from settable QNaNs.
=#

for (FL, I, UI, UPos, UNeg) in [(:Float64, :Int64, :UInt64, :0x7ff8000000000000, :0xfff8000000000000),
                                (:Float32, :Int32, :UInt32, :0x7fc00000, :0xffc00000),
                                (:Float16, :Int16, :UInt16, :0x7e00, :0xfe00) ]
  @eval begin
      @inline isjnan(x::$(UI))    = (x == $(UNeg)) || (x == $(UPos))
      @inline isqnan(x::$(UI))    = (x & $(UPos)) == $(UPos)
      @inline isqnanPos(x::$(UI)) = (x & $(UNeg)) == $(UPos)
      @inline isqnanNeg(x::$(UI)) = (x & $(UNeg)) == $(UNeg)
      @inline isjnan(x::$(FL))    = isjnan(reinterpret($(UI),x))
      @inline isqnan(x::$(FL))    = isqnan(reinterpret($(UI),x))
      @inline isqnanPos(x::$(FL)) = isqnanPos(reinterpret($(UI),x))
      @inline isqnanNeg(x::$(FL)) = isqnanPos(reinterpret($(UI),x))

      function qnan(si::$(I))
          u = reinterpret($(UI), abs(si))
          if (u & ~$(UPos)) == 0
              throw(ArgumentError("The magnitude of n cannot be 0, that NaN is used by Julia."))
          elseif (u > ~$(UNeg)) # 2^51-1
              throw(ArgumentError("The magnitude of n $(si) exceeds available QNaN."))
          end
          u |= ((si > 0) ? $(UPos) : $(UNeg))
          reinterpret($(FL),u)
      end

      function qnan(fp::$(FL))
          u = reinterpret($(UI), fp)
          if !isqnan(u)
              throw(ArgumentError("The value $(fp) ($(u)) is not a QNaN."))
          end
          a = u & ~$(UNeg)
          b =  reinterpret($(I),a)
          (((u & ($(UPos) $ $(UNeg))) == 0) ? b : -b)
      end
  end
end

@doc """
Module QNaN

  qnan(si::{Int16,Int32,Int64}) generates a quiet NaN with a payload of `si`

  qnan(fp::FloatingPoint) recovers the Signed payload from `fp`

  isqnan(`fp`::FloatingPoint) true iff `fp` is any QNaN

  isjnan(`fp`::FloatingPoint) true iff `fp` is a QNaN that Julia uses

  isqnanPos(`fp`::FloatingPoint) true iff isqnan(`fp`) and !signbit(`fp`)

  isqnanNeg(`fp`::FloatingPoint) true iff isqnan(`fp`) and  signbit(`fp`)

""" -> QNaN

@doc """

qnan(`si`::{Int16|32|64}) generates a quiet NaN using `si` for the payload

qnan(`fp`::FloatingPoint) recovers the signed payload from a quiet NaN `fp`

""" -> qnan

@doc """

isqnan(`fp`::FloatingPoint) true iff `fp` is any QNaN

isjnan(`fp`::FloatingPoint) true iff `fp` is one QNaN that Julia uses

""" -> isqnan

@doc """

isjnan(`fp`::FloatingPoint) true iff `fp` is one QNaN that Julia uses

isqnan(`fp`::FloatingPoint) true iff `fp` is any QNaN

""" -> isjnan

@doc """

isqnanPos(`fp`::FloatingPoint) true iff isqnan(`fp`) and !signbit(`fp`)

""" -> isqnanPos

@doc """

isqnanNeg(`fp`::FloatingPoint) true iff isqnan(`fp`) and  signbit(`fp`)

""" -> isqnanNeg

end # module


text:
# Configuration for running spatial simulation
#=
Recommended command line to run:
>  julia -L ContVarEvolution.jl run_cv.jl examples/example1
or
>  julia -p 4 -L ContVarEvolution.jl run_cv.jl examples/example1
=#
export simtype
@everywhere simtype = 2
#@everywhere const N = 8        # Meta-population size
#@everywhere const N_list = [25,50,100,200]        # Meta-population size list
const num_trials = 1
@everywhere const N_list = [20]        # Meta-population size list
const num_attributes_list = [1]        # number attributes for quantitative representation
#const N_mut_list = [0.5, 1.0, 2.0, 4.0]
const N_mut_list = [2.0]
const num_subpops = 1
const ngens = 120      # Generations after burn-in
const burn_in= 0.0    # generations of burn_in as a multiple of N
const ideal = 1.0
const fit_slope = 1.0
const neutral=false


text:
@inline modv2(x, n) = rem(x - one(x), n) + one(x)

struct HighDegreeInterpPoly{N,V}
    c::NTuple{N,V}
    w::V
end

struct ExponentPairs{N} end
ExponentPairs(::Val{N}) where {N} = ExponentPairs{N}()
Base.@pure ExponentPairs(N) = ExponentPairs{N}()

"""
@inline cmod(n, d) = n & (d - one(d))
function Base.getindex(::ExponentPairs{N}, i::T) where {N, T}
    d, r = divrem(i-one(T), T(N))
    r += one(T)
    r, one(T) + cmod(cmod(d, (T(N)-one(T))) + r, T(N))
end

function Base.getindex(::ExponentPairs{N}, i::T) where {N, T}
    d, r = fldmod1(i, T(N))
    r,  mod1( mod1(d, T(N)-one(T)) + r, T(N) )
end
"""
function Base.getindex(::ExponentPairs{N}, i::T) where {N, T}
    d, r = divrem(i-one(T), T(N))
    r += one(T)
    r, one(T) + ((d % (T(N)-one(T))) + r) % T(N)
end

Base.Vector(::ExponentPairs{N}, n::T) where {T,N} = getindex.(ExponentPairs{N}(), one(T):n)
like_size(::Type{Float64}, x) = Int64(x)
like_size(::Type{Float32}, x) = Cuint(x)

function coefs(seq::A, ::Val{N}) where {N, T, A <: AbstractArray{T,2}}
    n = size(seq, 2) #n * N = length(seq)
    #Need to iterative over pairs.

    l = cld( 2n, N ) #+ 1

    #cycles
    num_cycles = 2cld(n, binomial(N, 2))
    base_polynomial = fill( one(Complex{T}), num_cycles, N)

    #F it, lets stick with a vector!
    #Hmm, we actually want to guarantee last set does not feature overlap.
    #
    k = Vector(ExponentPairs{N}(), like_size(T, n))
    for (i,j) ∈ k
        for w ∈ 1:num_cycles
            base_polynomial[:,i]
            base_polynomial[:,j]
        end
    end

    cycle1 = binomial(N, 2)
    cycles1, remainder1 = divrem(n, cycle1)

    for i ∈ x

    end
end

struct ComplexViewPolynomial{T,A,D,N} where {T,A<:AbstractArray{T},D,N}
    ω::NTuple{D,T}
    c::A
end

@generated gen_omega(::Val{D}, ::Type{T} = Float64) where {D,T} = :(@ntuple $D i -> Complex{T}(exp(-2π*im*(i-one(T))/D)))

@generated function gen_omega_r(::Val{D}, ::Type{T} = Float64) where {D,T}
    l = (D-1) ÷ 2 #Truncated 1st element ( = 1 ), and last non-symmetric element if D is even ( = -1 )
    :(@ntuple $l i -> T(cos(2π*i/D)))
end
@generated function gen_omega_i(::Val{D}, ::Type{T} = Float64) where {D,T}
    l = (D-1) ÷ 2 #Truncated 1st elment ( = 0 ) and last non-symmetric element if D is even ( = 0 )
    :(@ntuple $l i -> T(sin(2π*i/D)))
end

function complex_polynomial_kernel(state, A::AbstractArray{T}, ω) where T
    d,p = @cartesianidx A #poly coeficient, and the variable in the multi-dimensional polynomial.
    return
end

"""
D really ought to be
D = cld(2N, P)+1
"""
struct EvaluationKernel{D,P,N,A} end

@generated function split_roots(roots::A, ::EvaluationKernel{D,P,N})
    where {D,P,N,T,A <: AbstractArray{T,2}}
    complete_cycles, leftover = divrem( N, P )

    quote
        root_poly = Array{T,2}(uninitialized, 2*$complete_cycles, $P)
        ind = 0
        @nexprs $P j -> begin
            @nexprs $complete_cycles i -> begin
                ind += 1
                k = modv2( j + i + cld( i, $P - 1) - 1 , $P )
                root_poly[2(i-1)+1,j] = roots[j,ind]
                root_poly[2i,k] = roots[k,ind]
            end
        end
        root_poly
    end
end
@generated function split_roots_and_inds(roots::A, ::EvaluationKernel{D,P,N})
                where {D,P,N,T,A <: AbstractArray{T,2}}
    complete_cycles, leftover = divrem( N, P )

    quote
        root_poly = Array{T,2}(uninitialized, 2*$complete_cycles, $P)
        root_inds = Array{T,2}(uninitialized, $complete_cycles, $P)
        ind = 0
        @nexprs $P j -> begin
            @nexprs $complete_cycles i -> begin
                ind += 1
                k = modv2( j + i + cld( i, $P - 1) - 1 , $P )
                root_inds[2(i-1)+1,j] = ind
                root_inds[2i,k] = ind
                root_poly[2(i-1)+1,j] = roots[j,ind]
                root_poly[2i,k] = roots[k,ind]
            end
        end
        root_poly,root_inds
    end
end
@generated function split_root_inds(::EvaluationKernel{D,P,N}, ::Type{T} = Cuint)
                where {D,P,N,T}
    complete_cycles, leftover = divrem( N, P )

    quote
        root_inds = Array{T,2}(uninitialized, $complete_cycles, $P)
        ind = 0
        @nexprs $complete_cycles i -> begin
            @nexprs $P j -> begin
                ind += 1
                k = modv2( j + i + cld( i, $P - 1) - 1 , $P )
                root_inds[2(i-1)+1,j] = ind
                root_inds[2i,k] = ind
            end
        end
        root_inds
    end
end
@inline function update_poly!(base_poly, ωᵣ, ωᵢ, d, i, r)
    dᵣ = 2d
    dᵢ = dᵣ+1
    ωrτ = fma($ωᵣ[d], -r , one(eltype(base_poly)))
    bpᵣ = base_poly[dᵣ, i]
    bpᵢ = base_poly[dᵢ, i]
    base_poly[dᵣ, i] = bpᵣ*ωrτ + bpᵢ*r*ωᵢ[d]
    base_poly[dᵢ, i] = bpᵢ*ωrτ - bpᵣ*r*ωᵢ[d]
end
@inline function initial_write_poly!(base_poly, ωᵣ, ωᵢ, d, i, r)
    base_poly[2d, i] = fma($ωᵣ[d], -r , one(eltype(base_poly)))
    base_poly[2d+1, i] = - r*ωᵢ[d]
end
@inline function initial_write_poly!(base_poly, d, i)
    base_poly[2d, i] = one(eltype(base_poly))
    base_poly[2d+1, i] = zero(eltype(base_poly))
end

@generated function ComplexViewPolynomial(roots::AbstractArray{T,2}, ::EvaluationKernel{D,P,N,A})
                where {D,P,N,T,A <: AbstractArray}
    ωᵣ = gen_omega_r(Val{D}(), T)
    ωᵢ = gen_omega_i(Val{D}(), T)
    complete_cycles, leftover = divrem( N, P )
    offset = complete_cycles * P
    half_P, Podd = divrem(P, 2)
    first_set, repeats = divrem(leftover, half_P)
    if first_set == 0
        first_set = leftover
        finish_initialization = half_P - first_set
        repeats = 0
        Podd = 0
    else
        first_set = half_P
        finish_initialization = 0
        if repeats != 0
            repeats -= Podd ##Now we need Podd && repeats >= 0 to do the Podd branch.
        else
            Podd = 0
        end
    end
    #Now, first set can assume it is writing into the matrix for the first time.
    #Then, repeats may modify earlier entires.
    #If not, we will have to fill out the rest of the matrix with appropriate 1s (r) and 0s (i).

    D_clipped, Deven = divrem( D - 1, 2)#Deven asks whether we have a trailing real.
    quote
        base_poly = Array{T,2}(uninitialized, $D, $P)
        #first_set is set to however many need to be initialized.
        @nexprs $first_set i -> begin
            ind = i+$offset
            p₁ = mod1( 2ind - 1, $P )
            p₂ = mod1( 2ind, $P )
            if p₁ < p₂
                r₁ = roots[p₁, ind]
                r₂ = roots[p₂, ind]

                base_poly[1, p₁] = one(T) - r₁
                base_poly[1, p₂] = one(T) - r₂

                @nexprs $D_clipped d -> begin
                    #Originally, I expanded things here.
                    #But that was ugly and unclear.
                    #Inlining + constant propagation mean ugliness does not help performance.
                    #If ugliness does help, some argue we should avoid it anyway.
                    #The clean way is likely to catch up, be more maintainable, etc.
                    #On average, algorithmic advantages >>>> micro optimizations.
                    #Algorithmic advantages are enabled by code we understand.
                    initial_write_poly!(base_poly, ωᵣ, ωᵢ, d, p₁, r₁)
                    initial_write_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                end
                if $Deven == 1
                    base_poly[$D, p₁] = one(T) + r₁
                    base_poly[$D, p₂] = one(T) + r₂
                end
            else
                r₂ = roots[p₂, ind]
                r₁ = roots[p₁, ind]

                base_poly[1, p₂] = one(T) - r₂
                base_poly[1, p₁] = one(T) - r₁
                @nexprs $D_clipped d -> begin
                    initial_write_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                    initial_write_poly!(base_poly, ωᵣ, ωᵢ, d, p₁, r₁)
                end
                if $Deven == 1
                    base_poly[$D, p₂] = one(T) + r₂
                    base_poly[$D, p₁] = one(T) + r₁
                end
            end
        end
        # These are dimensions that have not been hit yet.
        @nexprs $finish_initialization i -> begin
            ind = i+$offset+$first_set
            p₁ = mod1( 2ind - 1, $P )
            p₂ = mod1( 2ind, $P )
            if p₁ < p₂#Need to set entries to one and zero, as appropriate.
                base_poly[1, p₁] = one(T)
                base_poly[1, p₂] = one(T)
                @nexprs $D_clipped d -> begin
                    initial_write_poly!(base_poly, d, p₁)
                    initial_write_poly!(base_poly, d, p₂)
                end
                if $Deven == 1
                    base_poly[$D, p₂] = one(T)
                    base_poly[$D, p₁] = one(T)
                end
            else
                base_poly[1, p₂] = one(T)
                base_poly[1, p₁] = one(T)
                @nexprs $D_clipped d -> begin
                    initial_write_poly!(base_poly, d, p₂)
                    initial_write_poly!(base_poly, d, p₁)
                end
                if $Deven == 1
                    base_poly[$D, p₂] = one(T)
                    base_poly[$D, p₁] = one(T)
                end
            end
        end

        if $Podd == 1
            ind = 1+$offset+$first_set
            p₁ = mod1( 2ind - 1, $P ) ###This has to be the one that was not hit before.
            p₂ = mod1( 2ind, $P )
            r₂ = roots[p₂, ind]
            if p₁ < p₂#Need to set entries to one and zero, as appropriate.
                @nexprs $D_clipped d -> begin
                    initial_write_poly!(base_poly, d, p₁)
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                end
                if $Deven == 1
                    base_poly[$D, p₁] = one(T)
                    base_poly[$D, p₂] = one(T) + r₂
                end
            else
                @nexprs $D_clipped d -> begin
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                    initial_write_poly!(base_poly, d, p₁)
                end
                if $Deven == 1
                    base_poly[$D, p₂] = one(T) + r₂
                    base_poly[$D, p₁] = one(T)
                end
            end

        end

        @nexprs $repeats i -> begin
            ind = i+$offset+$first_set+$Podd
            p₁ = mod1( 2ind - 1, $P )
            p₂ = mod1( 2ind, $P )
            if p₁ < p₂
                r₁ = roots[p₁, ind]
                r₂ = roots[p₂, ind]
                base_poly[1, p₁] *= one(T) - r₁
                base_poly[1, p₂] *= one(T) - r₂
                @nexprs $D_clipped d -> begin
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₁, r₁)
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                end
                if $Deven == 1
                    base_poly[$D, p₁] *= one(T) + r₁
                    base_poly[$D, p₂] *= one(T) + r₂
                end
            else
                r₂ = roots[p₂, ind]
                r₁ = roots[p₁, ind]
                base_poly[1, p₂] *= one(T) - r₂
                base_poly[1, p₁] *= one(T) - r₁
                @nexprs $D_clipped d -> begin
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₂, r₂)
                    update_poly!(base_poly, ωᵣ, ωᵢ, d, p₁, r₁)
                end
                if $Deven == 1
                    base_poly[$D, p₂] *= one(T) + r₂
                    base_poly[$D, p₁] *= one(T) + r₁
                end
            end
        end

        roots_sorted = split_roots(roots, EvaluationKernel{$D,$P,$N}())
        for p ∈ 1:$P, i ∈ 1:2*$complete_cycles
            r = roots_sorted[i,p]
            base_poly[1, p] *= one(T) - r
            if $Deven == 1
                base_poly[$D, p] *= one(T) + r
            end
        end

        base_poly_a = A(base_poly)
        nodes = A(roots)
        roots_sorted_a = A(roots_sorted)

        ComplexViewPolynomial($ω, p)
    end
end

function meta_info_gbk(N, P, D)
    complete_cycles = N ÷ P
    L, Deven = divrem(D-1, 2)
    D_two_factor = Deven
    if Deven == 1
        for i ∈ 1:trunc(Int,log2(D))
            if D >> D_two_factor <<  D_two_factor != D
                D_two_factor -= 1
                break
            else
                D_two_factor += 1
            end
        end
    end
    complete_cycles, L, Deven, D_two_factor
end

"""
roots are PxN
basis_poly is Dx2PxN
base_poly is DxP
"""
@generated function generate_basis_kernel(state, root_inds::GPUArray{Tuple{Cuint,Cuint}} roots::GPUArray, basis_poly::GPUArray{Float32,4}, base_poly::GPUArray, ωᵣ::NTuple{L,T}, ωᵢ::NTuple{L,T}, ::EvaluationKernel{D,P,N}) where {D,P,N, L, T}
    complete_cycles, L, Deven, D_two_factor = meta_info_gbk(N, P, D)
    quote
        n = @linearidx rootinds
        @nextract $P p i -> root_inds[Cuint(i),n] #We skip first two p_inds are factored out.

        @nexprs 2 j -> begin
            r = roots[p_j,n]
            basis_poly[one(Cuint),p_j,n] = base_poly[one(Cuint),p_j] / ( one(T) - r )
            @nexprs $L i -> begin
                realpart = one(T) - r * ωᵣ[Cuint(i)]
                imagpart = r * ωᵢ[Cuint(i)]
                bpreal = base_poly[Cuint(2)*Cuint(i),p_j]
                bpimag = base_poly[one(Cuint)+Cuint(2)*Cuint(i),p_j]

                denom = abs2(realpart) + abs2(imagpart)
                basis_poly[Cuint(2)*Cuint(i),p_j,n] = ( bpreal * realpart + bpimag * imagpart ) / denom
                basis_poly[one(Cuint)+Cuint(2)*Cuint(i),p_j,n] = ( bpimag * realpart - bpreal * imagpart ) / denom
            end
            if $Deven == 1
                basis_poly[$D,p_j,n] = base_poly[$D,p_j] / ( one(T) + r )
            end
        end

        @nexprs $P-2 j -> begin
            @nexprs $D d -> begin
                basis_poly[d,p_{j+2},n] = base_poly[d,p_{j+2}]
            end
        end
        @nexprs 2 j -> begin
            @nexprs $D d -> begin
                basis_poly[d,p_j+$P,n] = base_poly[d,p_j]
            end
        end
        @nexprs $P-2 j -> begin
            r = roots[p_{j+2},n]
            basis_poly[one(Cuint),$P+p_{j+2},n] = basis_poly[one(Cuint),p_{j+2},n] * ( one(T) - r )
            @nexprs $L i -> begin
                realpart = one(T) - r * ωᵣ[Cuint(i)]
                imagpart = r * ωᵢ[Cuint(i)]
                bpreal = basis_poly[Cuint(2)*Cuint(i),p_{j+2},n]
                bpimag = basis_poly[one(Cuint)+Cuint(2)*Cuint(i),p_{j+2},n]

                basis_poly[Cuint(2)*Cuint(i),$P+p_{j+2},n] = bpreal * realpart - bpimag * imagpart
                basis_poly[one(Cuint)+Cuint(2)*Cuint(i),$P+p_{j+2},n] = bpimag * realpart + bpreal * imagpart
            end
            if $Deven == 1
                basis_poly[$D,$P+p_{j+2},n] = basis_poly[$D,p_{j+2}] / ( one(T) + r )
            end
        end
        ### Now you have to ifft all of these.
        ### For now, this is just the naive algorithm.
        ### An FFT algorithm, like Cooley-Tukey, should be implemented later if this is actually a good idea.
        ### Should be easy to code an optimal one.
        ### May need to adjust D automatically via some heurstiics on number of calculations.
        @nexprs $P j -> begin
            @nexprs $D k -> begin
                @nexprs l -> begin

                    @nexprs 2 i -> begin

                    end
                end
            end
        end
        nothing
    end
end

@inline function base_poly_ind(n::T, d::T) where T
    t, r = divrem(n-one(T), d)
    t + one(T), r + one(T)
end
"""
Generates base polynomial.
L = (D-1) ÷ 2
"""
@generated function generate_base_poly_kernel(state, base_poly::GPUArray, root_split::GPUArray, ωᵣ::NTuple{L,T}, ωᵢ::NTuple{L,T}, ::EvaluationKernel{D,P,N}) where {D,P,N,L,T}
    complete_cycles = N ÷ P
    quote
        p, d = base_poly_ind(GPUArrays.linear_index(state), Cuint($L))
        ωr = ωᵣ[p]
        ωi = ωᵢ[p]
        pᵣ = T(2)*p
        pᵢ = pᵣ + one(T)
        cᵣ = base_poly[pᵣ, d]
        cᵢ = base_poly[pᵢ, d]
        @nexprs 2*$complete_cycles rind -> begin
            r = root_split[Cuint(rind), d]
            rᵣ = one(Cuint) - r*ωr
            rᵢ = r*ωi
            cᵣ, cᵢ = ( rᵣ*cᵣ - rᵢ*cᵢ, rᵣ*cᵢ + rᵢ*cᵣ )
        end
        base_poly[pᵣ, d] = cᵣ
        base_poly[pᵢ, d] = cᵢ
        nothing
    end
end


text:
export PrioritizedSweepingSampleModel

using DataStructures: PriorityQueue, dequeue!

import StatsBase: sample

"""
    PrioritizedSweepingSampleModel(θ::Float64=1e-4)

See more details at Section (8.4) on Page 168 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
mutable struct PrioritizedSweepingSampleModel <: AbstractEnvironmentModel
    experiences::Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}
    PQueue::PriorityQueue{Tuple{Any,Any},Float64}
    predecessors::Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}
    θ::Float64
    sample_count::Int
    PrioritizedSweepingSampleModel(θ::Float64 = 1e-4) =
        new(
            Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}(),
            PriorityQueue{Tuple{Any,Any},Float64}(Base.Order.Reverse),
            Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}(),
            θ,
            0,
        )
end

function extract_experience(
    t::AbstractTrajectory,
    model::PrioritizedSweepingSampleModel,
)
    if length(t) > 0
        get_trace(t, :state)[end],
        get_trace(t, :action)[end],
        get_trace(t, :reward)[end],
        get_trace(t, :terminal)[end],
        get_trace(t, :next_state)[end]
    else
        nothing
    end
end

function RLBase.update!(m::PrioritizedSweepingSampleModel, t::AbstractTrajectory, p::AbstractPolicy)
    experience = extract_experience(t, m)
    isnothing(experience) || update!(m, (experience..., get_priority(p, experience)))
end

function RLBase.update!(m::PrioritizedSweepingSampleModel, transition::Tuple)
    s, a, r, d, s′, P = transition
    m.experiences[(s, a)] = (r, d, s′)
    if P >= m.θ
        m.PQueue[(s, a)] = P
    end
    if !haskey(m.predecessors, s′)
        m.predecessors[s′] = Set{Tuple{Any,Any,Float64,Bool}}()
    end
    push!(m.predecessors[s′], (s, a, r, d))
end

function sample(m::PrioritizedSweepingSampleModel)
    if length(m.PQueue) > 0
        s, a = dequeue!(m.PQueue)
        r, d, s′ = m.experiences[(s, a)]
        m.sample_count += 1
        s, a, r, d, s′
    else
        nothing
    end
end


text:
"""
    linearRadius(initRadius::Float64, iteration::Int64, decay::String, epochs::Int64)

Return a neighbourhood radius. Use as the `radiusFun` parameter for `trainGigaSOM`.

# Arguments
- `initRadius`: Initial Radius
- `finalRadius`: Final Radius
- `iteration`: Training iteration
- `epochs`: Total number of epochs
"""
function linearRadius(
    initRadius::Float64,
    finalRadius::Float64,
    iteration::Int64,
    epochs::Int64,
)

    scaledTime = scaleEpochTime(iteration, epochs)
    return initRadius * (1 - scaledTime) + finalRadius * scaledTime
end

"""
    expRadius(steepness::Float64)

Return a function to be used as a `radiusFun` of `trainGigaSOM`, which causes
exponencial decay with the selected steepness.

Use: `trainGigaSOM(..., radiusFun = expRadius(0.5))`

# Arguments
- `steepness`: Steepness of exponential descent. Good values range
  from -100.0 (almost linear) to 100.0 (really quick decay).

"""
function expRadius(steepness::Float64 = 0.0)
    return (initRadius::Float64, finalRadius::Float64, iteration::Int64, epochs::Int64) ->
        begin

            scaledTime = scaleEpochTime(iteration, epochs)

            if steepness < -100.0
                # prevent floating point underflows
                error("Sanity check: steepness too low, use linearRadius instead.")
            end

            # steepness is simulated by moving both points closer to zero
            adjust = finalRadius * (1 - 1.1^(-steepness))

            if initRadius <= 0 || (initRadius - adjust) <= 0 || finalRadius <= 0
                error(
                    "Radii must be positive. (Possible alternative cause: steepness is too high.)",
                )
            end

            initRadius -= adjust
            finalRadius -= adjust

            return adjust + initRadius * ((finalRadius / initRadius)^scaledTime)
        end
end

"""
    gridRectangular(xdim, ydim)

Create coordinates of all neurons on a rectangular SOM.

The return-value is an array of size (Number-of-neurons, 2) with
x- and y- coordinates of the neurons in the first and second
column respectively.
The distance between neighbours is 1.0.
The point of origin is bottom-left.
The first neuron sits at (0,0).

# Arguments
- `xdim`: number of neurons in x-direction
- `ydim`: number of neurons in y-direction
"""
function gridRectangular(xdim, ydim)

    grid = zeros(Float64, (xdim * ydim, 2))
    for ix = 1:xdim
        for iy = 1:ydim
            grid[ix+(iy-1)*xdim, 1] = ix - 1
            grid[ix+(iy-1)*xdim, 2] = iy - 1
        end
    end
    return grid
end

"""
    gaussianKernel(x, r::Float64)

Return the value of normal distribution PDF (σ=`r`, μ=0) at `x`
"""
function gaussianKernel(x, r::Float64)

    return Distributions.pdf.(Distributions.Normal(0.0, r), x)
end

function bubbleKernelSqScalar(x::Float64, r::Float64)
    if x >= r
        return 0
    else
        return sqrt(1 - x / r)
    end
end

"""
    bubbleKernel(x, r::Float64)

Return a "bubble" (spherical) distribution kernel.

"""
function bubbleKernel(x, r::Float64)
    return bubbleKernelSqScalar.(x .^ 2, r^2)
end

"""
    thresholdKernel(x, r::Float64)

Simple FlowSOM-like hard-threshold kernel
"""
function thresholdKernel(x, r::Float64, maxRatio = 4 / 5, zero = 1e-6)
    if r >= maxRatio * maximum(x) #prevent smoothing everything to a single point
        r = maxRatio * maximum(x)
    end
    return zero .+ (x .<= r)
end

"""
    distMatrix(metric=Chebyshev())

Return a function that uses the `metric` (compatible with metrics from package `Distances`) calculates distance matrixes from normal row-wise data matrices, using the `metric`.

Use as a parameter of `trainGigaSOM`.
"""
function distMatrix(metric = Chebyshev())
    return (grid::Matrix{Float64}) -> begin
        n = size(grid, 1)
        dm = zeros(Float64, n, n)

        for i = 1:n
            for j = 1:n
                dm[i, j] = metric(grid[i, :], grid[j, :])
            end
        end

        return dm::Matrix{Float64}
    end
end


text:
struct DataDryad <: DataRepo
end

base_url(::DataDryad) = "https://datadryad.org/resource/doi:"

description(repo::DataDryad, mainpage) = replace(text_only(first(eachmatch(sel".article-abstract", mainpage.root))), "Abstract " => "")

function author(::DataDryad, mainpage)
    author = missing
    try
        author = text_only(first(eachmatch(sel".pub-authors a", mainpage.root)))
    catch
        author = string(split(text_only(first(eachmatch(sel".pub-authors", mainpage.root))), ", ")[1])
    end
    [author]
end

license(::DataDryad, mainpage) = getattr(first(eachmatch(sel".single-image-link", mainpage.root)), "href")

function published_date(::DataDryad, mainpage)
    dateelem = eachmatch(sel".publication-header p", mainpage.root)
    replace(text_only(dateelem[length(dateelem)-1]), "Date Published: " => "")
end

paper_cite(::DataDryad, mainpage) = text_only(first(eachmatch(sel".citation-sample", mainpage.root)))

dataset_cite(::DataDryad, mainpage) = replace(text_only(last(eachmatch(sel".publication-header p", mainpage.root))), "DOI: " => "")

function get_urls(repo::DataDryad, page)
    urls = []
    links = eachmatch(sel".package-file-description tbody tr td a", page.root)
    for link in links
        urlhref = getattr(link, "href")
        dryadurl = "https://datadryad.org" * string(urlhref)
        if occursin("Download" , string(link.parent.parent))
            push!(urls, dryadurl)
        end
    end
    urls
end

function get_checksums(repo::DataDryad, page)
    checksums = []
    links = eachmatch(sel"a", page.root)
    regex = r"\bresource\/doi:[0-9]*.[0-9]*\/dryad.[a-z, 0-9]*\/[0-9]+\b"
    for checksum_link in links
        if occursin(regex,checksum_link.attributes["href"])
            checksum = match(regex, checksum_link.attributes["href"])
            url = replace(checksum.match, "resource/doi:" => "https://datadryad.org/mn/checksum/doi:")
            md5 = (:md5, text_only(getpage(url).root))
            push!(checksums, md5)
        end
    end
    if length(checksums) > 0
        @info("The generated registration block uses md5 hash, " *
            "the MD5.jl package must be loaded to run the registration")
    end
    checksums
end

function data_fullname(::DataDryad, mainpage)
    # mainpage = replace(mainpage, "resource", "mn/object")
    text_only(first(eachmatch(sel".pub-title", mainpage.root)))
end


text:
module InterpretBig
using Base: Integer, Bool
using ..LambdaC: AExpr, @λc_str
export empty_env, Environment, interpret, init_interpret, std_env, si, init_sinterpret
import MLStyle
using SoftPredicates
using SoftPredicates: DualSoftBool

const Environment = NamedTuple
empty_env() = NamedTuple()
std_env() = empty_env()# (unif = interpret(λc"(λ i (λ ω (ω i)))", empty_env()),)

"Produces new environment Γ' s.t. `Γ(x) = v` (and everything else unchanged)"
update(Γ, x::Symbol, v) = merge(Γ, NamedTuple{(x,)}((v,)))

# unif(i, ω) = ω[(:unif, i)]

        #  (let unif (λ j (λ w (w j)))
        #  (let x (λ ω ((unif 1) ω))
        #  (let ω (λ i (if (== i 1) then 0.2 else 0.5))
        #    (x ω))))

prim_to_func = Dict(:+ => +,
                    :- => -,
                    :* => *,
                    :/ => /,
                    :& => &,
                    :! => !,
                    :| => |,
                    :> => >,
                    :>= => >=,
                    :< => <,
                    :<= => <=,
                    :(==) => ==,
                    :>ₛ => >ₛ)

isprim(f) = f in keys(prim_to_func)
primapl(f, x...) = prim_to_func[f](x...)

# For Lazy Evaluation
struct Closure
  env::Environment
  expr
end

clo(env, expr) = Closure(env, expr)

function arg(aex::AExpr)
  arr = [aex.head, aex.args...]
  MLStyle.@match arr begin
    [:fn, arg, body]                                       => arg
    _                                                      => error("Could find body of $arr")
  end
end

function body(aex::AExpr)
  arr = [aex.head, aex.args...]
  MLStyle.@match arr begin
    [:fn, arg, body]                                       => body
    _                                                      => error("Could find body of $arr")
  end
end

"substitute x for v withih expr"
function sub(aex::AExpr, (x, v))
  # print("SUb")
  # @show aex
  # @show x
  # @show v
  arr = [aex.head, aex.args...]
  # next(x) = interpret(x, Γ)
  isaexpr(x) = x isa AExpr
  MLStyle.@match arr begin
    [:fn, args, body]                                       => AExpr(:fn, args, sub(body, x => v))
    [:call, f, arg1] && if isprim(f) end                    => AExpr(:call, f, sub(arg1, x => v))
    [:call, f, arg1, arg2] && if isprim(f) end              => AExpr(:call, f, sub(arg1, x => v),  sub(arg2, x => v))
    [:if, c, t, e]                                          => AExpr(:if, sub(c, x => v), sub(t, x => v), sub(e, x => v))
    [:let, x, val, body] && if x ∉ keys(Γ) end              => AExpr(:let, x, val, sub(body, x => v))
    [:call, t1, t2]                                         => AExpr(:call, sub(t1, x => v),  sub(t2, x => v))
    _                                                       => error("Could not sub $arr")
  end
end

sub(aex::Symbol, (x, v)) = aex == x ? v : aex
sub(aex::Real, (x, v)) = aex

function interpret(aex::AExpr, Γ::Environment)
  arr = [aex.head, aex.args...]
  # println()
  # println("Env:")
  # display(Γ)
  # @show arr
  # next(x) = interpret(x, Γ)
  isaexpr(x) = x isa AExpr
  t = MLStyle.@match arr begin
    [:fn, args, body]                                       => clo(Γ, aex)
    [:call, f, arg1] && if isprim(f) end                    => primapl(f, interpret(arg1, Γ))
    [:call, f, arg1, arg2] && if isprim(f) end              => primapl(f, interpret(arg1, Γ), interpret(arg2, Γ))
    [:if, c, t, e] && if interpret(c, Γ) == true end        => interpret(t, Γ)
    [:if, c, t, e] && if interpret(c, Γ) == false end       => interpret(e, Γ)
    [:let, x, val, body] && if x ∉ keys(Γ) end              => interpret(body, update(Γ, x, clo(Γ, val)))
    [:call, t1, t2]                                         => let clo_ = interpret(t1, Γ),
                                                                   v1 = interpret(t2, Γ),
                                                                   t3_up = sub(body(clo_.expr), arg(clo_.expr) => v1)
                                                                  #  o = println("  --  type of ", typeof(t3_up))
                                                                 interpret(t3_up, clo_.env)
                                                                end
    _                                                       => error("Could not interpret $arr")
  end
  # println("FINSIH", arr)
  # @show(t)
  t
end

interpret(c::Closure, Γ) = c

canparse(t, x) = !isnothing(tryparse(t, x))

function interpret(x::Symbol, Γ::Environment)
  if x == Symbol("false")
    false
  elseif x == Symbol("true")
    true
  else
    MLStyle.@match x begin
      x::Symbol  && if x ∈ keys(Γ) end                        => interpret(Γ[x].expr, Γ[x].env)  # Var
      x && if canparse(Float64, string(x)) end                => parse(Float64, string(x))
      _                                                       => error("Could not interpret $x")
    end
  end
end

interpret(x::Integer, Γ) = x
interpret(x::AbstractFloat, Γ) = x
interpret(x::Bool, Γ) = x

init_interpret(aex) = interpret(aex, std_env())

## Soft Cosntraints

struct CtxVal
  val
  ctx
end

a ⊕ b = error()

function soft_interpret(aex, Γ, ctx)
  arr = [aex.head, aex.args...]
  # next(x) = interpret(x, Γ)
  MLStyle.@match arr begin
    [:fn, args, body]                                       => clo(Γ, aex)
    [:call, f, arg1] && if isprim(f) end                    => let cv1 = si(arg1, Γ, ctx)
                                                                 CtxVal(primapl(f, cv1.val), cv1.ctx)
                                                               end
    [:call, f, arg1, arg2] && if isprim(f) end              => let cv1 = si(arg1, Γ, ctx),
                                                                   cv2 = si(arg2, Γ, cv1.ctx)
                                                                   fx = primapl(f, cv1.val, cv2.val)
                                                                 CtxVal(fx, cv2.ctx) # zt : do we need this composition?
                                                                end
    [:if, c, t, e] && if si(c, Γ, ctx).val == true end      => si(t, Γ, si(c, Γ, ctx).ctx)
    [:if, c, t, e] && if si(c, Γ, ctx).val == false end     => si(e, Γ, si(c, Γ, ctx).ctx)
    [:let, x, val, body] && if x ∉ keys(Γ) end              => si(body, update(Γ, x, clo(Γ, val)), ctx)
    [:call, t1, t2]                                         => let clo_ = si(t1, Γ, ctx),
                                                                   v1 = si(t2, Γ, clo_.ctx),
                                                                   t3_up = sub(body(clo_.val.expr), arg(clo_.val.expr) => v1.val)
                                                                  #  o = println("  --  type of ", typeof(t3_up))
                                                                 si(t3_up, clo_.env, v1.ctx)
                                                                end
    [:cond, t1, t2]
    _                                                       => error("Could not interpret $arr")
  end
end

soft_interpret(x::Integer, Γ, ctx) = CtxVal(x, ctx)
soft_interpret(x::AbstractFloat, Γ, ctx) = CtxVal(x, ctx)
soft_interpret(x::Bool, Γ, ctx) = CtxVal(x, ctx)'

function soft_interpret(x::Symbol, Γ::Environment, ctx)
  if x == Symbol("false")
    CtxVal(false, ctx)
  elseif x == Symbol("true")
    CtxVal(true, ctx)
  else
    MLStyle.@match x begin
      x::Symbol  && if x ∈ keys(Γ) end                        => si(Γ[x].expr, Γ[x].env, ctx)  # Var
      x && if canparse(Float64, string(x)) end                => CtxVal(parse(Float64, string(x)), ctx)
      _                                                       => error("Could not interpret $x")
    end
  end
end
# end

const si = soft_interpret

init_sinterpret(t) = soft_interpret(t, empty_env(), one(DualSoftBool{Float64}))

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# surrogate_ground_truth_mnist_7_or_9_mala
#
# This script runs a long chain based on the
# Metropolis-adjusted Langevin algorithm (MALA). This is used as a gold
# standard as this procedure yields an asymptotically unbiased chain
# and thus is used to compare the efficiency of the biased chains produced
# in the experiment mnist_7_or_9_sgfs_sgld.

using StatsFuns: logistic
using SteinDistributions: SteinLogisticRegressionGaussianPrior
using SteinSamplers: runmala

include("experiment_utils.jl")
include("mnist/mnist_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7
# setup the data
data_set = "mnist"
# number of components in reduced matrix
L = 50
# labels of the y data to keep around
labels = [7,9]
# get eigenfaces
V = mnist_eigenfaces(L, labels)
# setup training
X, y, Xtest, ytest = mnist_reduce_dimension_pca(L, labels)
# get the training and test sizes
ntrain = size(X, 1); ntest = size(Xtest, 1)
# give X an intercept term
X = hcat(ones(ntrain) ./ sqrt(ntrain), X)
# give Xtest an intercept term (the sqrt(ntrain) is not a mistake!)
Xtest = hcat(ones(ntest) ./ sqrt(ntrain), Xtest)
# make y be -1 for 7 and 1 for 9
y = 2.0 * (y .== 9) - 1.0
# make ytest be 0 for 7 and 1 for 9
ytest = round(Int, ytest .== 9)
# Parameter dimension
d = L + 1
# name of prior + posterior
distname = "logisticgaussianprior"
# sigma2 term for prior
sigma2 = 1.0
# setup distribution
target = SteinLogisticRegressionGaussianPrior(X, y, zeros(d), sigma2 .* ones(d))
# the max number of samples to run in the solver
n = 500000
# number of burnin samples
b = 10000
# the epsilons to run the experiemnt
epsilons = [5e-2]
# initial guess for beta
beta0 = zeros(d)

for epsilon = epsilons
    # Draw MALA samples
    (betas_mala, numgrad_mala) =
        @setseed runmala(target, beta0;
                         epsilonfunc=t -> epsilon,
                         numiter=n+b)

    betas_mala = betas_mala[(b+1):(b+n),:]
    # beta estimates
    betahat = vec(mean(betas_mala, 1))
    sigmahat = cov(betas_mala)
    # Keep track of predictive probabilities
    predprobs = logistic(Xtest * betas_mala')
    yhat = vec(mean(predprobs, 2))
    # Package and save results
    summary_data = Dict{Any, Any}(
        "epsilon" => epsilon,
        "distname" => distname,
        "n" => n,
        "d" => d,
        "seed" => seed,
        "eigenfaces" => V,
        "betahat" => betahat,
        "sigmahat" => sigmahat,
        "mean.pred.probs" => yhat,
        "true.pred" => ytest,
        "epsilon" => epsilon,
        "numgrad" => numgrad_mala,
        "sampler" => "MALA-long-nodal",
    )
    save_json(
        summary_data;
        dir="mnist_7_or_9_mala",
        distname=distname,
        dataset=data_set,
        n=n,
        epsilon=epsilon,
        priorvariance=sigma2,
        d=d,
        seed=seed
    )
end


text:
export Scale

"""
`Scale(α::Number,A::AbstractOperator)`

Shorthand constructor:

`*(α::Number,A::AbstractOperator)`

Scale an `AbstractOperator` by a factor of `α`.

```julia
julia> A = FiniteDiff((10,2))
δx  ℝ^(10, 2) -> ℝ^(9, 2)

julia> S = Scale(10,A)
αδx  ℝ^(10, 2) -> ℝ^(9, 2)

julia> 10*A         #shorthand
αℱc  ℝ^10 -> ℝ^10

```

"""
struct Scale{T <: RealOrComplex, L <: AbstractOperator} <: AbstractOperator
  coeff::T
  coeff_conj::T
  A::L
end

# Constructors

function Scale(coeff::T, L::R) where {T <: RealOrComplex, R <: AbstractOperator}
    coeff_conj = conj(coeff)
    coeff, coeff_conj = promote(coeff, coeff_conj)
    cT = codomainType(L)
    isCodomainReal = typeof(cT) <: Tuple ? all([t <: Real for t in cT]) : cT <: Real
    if isCodomainReal && T <: Complex
        error("Cannot Scale AbstractOperator with real codomain with complex scalar. Use `DiagOp` instead.")
    end
    Scale{typeof(coeff), R}(coeff, coeff_conj, L)
end

# Special Constructors
# scale of scale
Scale(coeff::T2, L::S) where {T1 <: RealOrComplex,
			      T2 <: RealOrComplex,
			      R <: AbstractOperator,
			      S <: Scale{T1, R}}=
Scale(*(promote(coeff,L.coeff)...), L.A)
# scale of DiagOp
Scale(coeff::T,L::DiagOp) where {T<:RealOrComplex} = DiagOp(coeff*diag(L))

# Mappings

function mul!(y::C, L::Scale{T, A}, x::D) where {T, C <: AbstractArray, D, A <: AbstractOperator}
  mul!(y, L.A, x)
  y .*= L.coeff
end

function mul!(y::C, L::Scale{T, A}, x::D) where {T, C <: Tuple, D, A <: AbstractOperator}
  mul!(y, L.A, x)
  for k in eachindex(y)
    y[k] .*= L.coeff
  end
end

function mul!(y::D, S::AdjointOperator{Scale{T, A}}, x::C) where {T, C, D <: AbstractArray, A <: AbstractOperator}
    L = S.A
    mul!(y, L.A', x)
    y .*= L.coeff_conj
end

function mul!(y::D, S::AdjointOperator{Scale{T, A}}, x::C) where {T, C, D <: Tuple, A <: AbstractOperator}
    L = S.A
    mul!(y, L.A', x)
    for k in eachindex(y)
        y[k] .*= L.coeff_conj
    end
end

# Properties

size(L::Scale) = size(L.A)

domainType(L::Scale) = domainType(L.A)
codomainType(L::Scale) = codomainType(L.A)

is_linear(L::Scale) = is_linear(L.A)
is_null(L::Scale) = is_null(L.A)
is_eye(L::Scale) = is_diagonal(L.A)
is_diagonal(L::Scale) = is_diagonal(L.A)
is_invertible(L::Scale) = L.coeff == 0 ? false : is_invertible(L.A)
is_AcA_diagonal(L::Scale) = is_AcA_diagonal(L.A)
is_AAc_diagonal(L::Scale) = is_AAc_diagonal(L.A)
is_full_row_rank(L::Scale) = is_full_row_rank(L.A)
is_full_column_rank(L::Scale) = is_full_column_rank(L.A)

fun_name(L::Scale) = "α$(fun_name(L.A))"
fun_type(L::Scale) = fun_type(L.A)

diag(L::Scale) = L.coeff*diag(L.A)
diag_AcA(L::Scale) = (L.coeff)^2*diag_AcA(L.A)
diag_AAc(L::Scale) = (L.coeff)^2*diag_AAc(L.A)
remove_displacement(S::Scale) = Scale(S.coeff, S.coeff_conj, remove_displacement(S.A) )


text:
using IncrementalInference
using Test

@testset "saving to and loading from FileDFG" begin

  fg = generateCanonicalFG_Kaess()
  addVariable!(fg, :x4, ContinuousScalar)
  addFactor!(fg, [:x2;:x3;:x4], LinearRelative(Normal()), multihypo=[1.0;0.6;0.4])

  saveFolder = "/tmp/dfg_test"
  saveDFG(fg, saveFolder)
  # VERSION above 1.0.x hack required since Julia 1.0 does not seem to havfunction `splitpath`
  if v"1.1" <= VERSION
    retDFG = initfg()
    retDFG = loadDFG!(retDFG, saveFolder)
    Base.rm(saveFolder*".tar.gz")

    @test symdiff(ls(fg), ls(retDFG)) == []
    @test symdiff(lsf(fg), lsf(retDFG)) == []

    @show getFactor(fg, :x2x3x4f1).solverData.multihypo
    @show getFactor(retDFG, :x2x3x4f1).solverData.multihypo

    # check for match
    @test getFactor(fg, :x2x3x4f1).solverData.multihypo - getFactor(retDFG, :x2x3x4f1).solverData.multihypo |> norm < 1e-10
    @test getFactor(fg, :x2x3x4f1).solverData.certainhypo - getFactor(retDFG, :x2x3x4f1).solverData.certainhypo |> norm < 1e-10
  end

end

@testset "saving to and loading from FileDFG with nullhypo, eliminated, solveInProgress" begin
  fg = generateCanonicalFG_Kaess()
  addVariable!(fg, :x4, ContinuousScalar)
  addFactor!(fg, [:x2;:x3;:x4], LinearRelative(Normal()), multihypo=[1.0;0.6;0.4])
  addFactor!(fg, [:x1;], Prior(Normal(10,1)), nullhypo=0.5)

  solveTree!(fg)

  #manually change a few fields to test if they are preserved
  fa = getFactor(fg, :x2x3x4f1)
  fa.solverData.eliminated = true
  fa.solverData.solveInProgress = 1
  fa.solverData.nullhypo = 0.5

  saveFolder = "/tmp/dfg_test"
  saveDFG(fg, saveFolder)

  retDFG = initfg()
  loadDFG!(retDFG, saveFolder)
  Base.rm(saveFolder*".tar.gz")

  @test issetequal(ls(fg), ls(retDFG))
  @test issetequal(lsf(fg), lsf(retDFG))

  @show getFactor(fg, :x2x3x4f1).solverData.multihypo
  @show getFactor(retDFG, :x2x3x4f1).solverData.multihypo

  # check for match
  @test isapprox(getFactor(fg, :x2x3x4f1).solverData.multihypo, getFactor(retDFG, :x2x3x4f1).solverData.multihypo)
  @test isapprox(getFactor(fg, :x2x3x4f1).solverData.certainhypo, getFactor(retDFG, :x2x3x4f1).solverData.certainhypo)

  fb = getFactor(retDFG, :x2x3x4f1)
  @test fa == fb

  fa = getFactor(fg, :x1f2)
  fb = getFactor(retDFG, :x1f2)

  @test fa == fb
end


text:
# The @compat macro is used to implement compatibility rules that require
# syntax rewriting rather than simply new function/constant/module definitions.

export @compat

function _compat(ex::Expr)
    if ex.head === :quote && isa(ex.args[1], Symbol)
        # Passthrough
        return ex
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end

_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end

# Only the import-renaming usage results in multiple arguments.
# On versions where renaming is supported natively it parses as a
# single expression and thus takes the single-argument @compat
# code-path above
macro compat(ex, exprs...)
    @assert VERSION < v"1.6.0-DEV.1157"
    if length(exprs) == 2 && Meta.isexpr(ex, :import, 1) && Meta.isexpr(ex.args[1], :.) &&
        exprs[1] === :as && exprs[2] isa Symbol
        # Matches "import Module(.Submodule) as a"
        lhs = Symbol[ex.args[1].args...]
        alias = exprs[2]
        return _create_import_expression([lhs => alias])
    elseif Meta.isexpr(ex, [:import, :using], 1) &&
           Meta.isexpr(ex.args[1], :(:)) && length(exprs) % 2 == 0 &&
           all(x -> x === :as, exprs[1:2:end-1]) &&
           all(x -> Meta.isexpr(x, :tuple, 2), exprs[2:2:end-2]) &&
           exprs[end] isa Symbol
        # Matches "(import|using) Module(.Submodule): a as b(, c as d)
        syms = Symbol[ ex.args[1].args[2].args[1] ]
        foreach(x -> x isa Symbol ? push!(syms, x) : append!(syms, x.args), exprs[2:2:end])

        path_aliases = Pair{Vector{Symbol},Symbol}[]
        init = ex.args[1].args[1].args
        for i in 1:2:length(syms)
            push!(path_aliases, Symbol[init; syms[i]] => syms[i+1])
        end
        return _create_import_expression(path_aliases)
    else
        throw(ArgumentError("invalid use of @compat"))
    end
end

function _create_import_expression(path_aliases::Vector{Pair{Vector{Symbol}, Symbol}})
    # Create an gensymd baremodule to hide names in
    s = gensym()
    # Create all import/const exprs
    import_exprs = Expr[]
    const_exprs = Expr[]
    for (path, alias) in path_aliases
        import_expr = Expr(:import, Expr(:., path...))
        push!(import_exprs, import_expr)
        rhs_expr = Expr(:escape, Expr(:., s, QuoteNode(last(path))))
        const_expr = Expr(:const, Expr(:global, Expr(:(=), alias, rhs_expr)))
        push!(const_exprs, const_expr)
    end
    module_expr = Expr(:module, false, Expr(:escape, s), Expr(:block, import_exprs...))
    return_expr = Expr(:toplevel, module_expr, const_exprs..., nothing)
    return return_expr
end


text:
@recipe function f(p::AbstractVolumePrimitive)
    fs = surfaces(p)
    # In principle it would be better to directly get the edges of the primitive
    # Can be improved later...
    linecolor --> :black
    if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
        aspect_ratio --> 1.0
    end
    @series begin
        label --> "$(nameof(typeof(p)))"
        [fs...]
    end
end


text:
using MeuExemplo
using Test

@testset "MeuExemplo.jl" begin
    # Write your own tests here.
end


text:
# Bio.Intervals
# =============
#
# Module for genomic intervals.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Intervals

export
    Strand,
    Interval,
    IntervalCollection,
    IntervalStream,
    STRAND_NA,
    STRAND_POS,
    STRAND_NEG,
    STRAND_BOTH,
    coverage,
    isoverlapping,
    seqname,
    strand,
    BED,
    BEDMetadata,
    BEDInterval,
    BigBed,
    BigWig

import ..Ragel: tryread!
export tryread!

using
    BufferedStreams,
    Colors,
    Compat,
    IntervalTrees,
    Libz,
    Bio.Ragel,
    Bio.StringFields

import Bio.IO:
    FileFormat,
    AbstractParser,
    AbstractWriter

using Base.Collections:
    heappush!,
    heappop!

import Iterators

include("strand.jl")
include("interval.jl")
include("stream_buffer.jl")
include("intervalcollection.jl")
include("intervalstream.jl")

# Parsing file types
include("bed.jl")
include("bigbed.jl")

end # module Intervals


text:
using DataFrames
using CSV
using PGFPlotsX

push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\pgfplotsset{compat=1.7}")
df_gpu = DataFrame(CSV.File("benchmark-results/bench_diffusion_scaling_gpu.csv"))

op_intensity_shmem = df_gpu[1, :Intensity]
op_intensity_hide_comm = df_gpu[2, :Intensity]

ax = Axis(@pgf({ylabel=raw"[Flop/Bytes]",
                ylabel_style={align="center", font=raw"\Large"},
                title="Operational Intensity",
                title_style={font=raw"\large"},
                xmajorticks="false",
                ybar, ymin=0, bar_width=0.5, area_legend, legend_pos="outer north east"}),
          PlotInc(Table(:x => 0, :y => [op_intensity_shmem])),
          PlotInc(Table(:x => 0, :y => [op_intensity_hide_comm])),
          Legend(["shared memory", "hide communication"])
)

pgfsave("figs/diffusion_scaling_experiments_operational_intensity.png", ax)
pgfsave("figs/diffusion_scaling_experiments_operational_intensity.pdf", ax)


text:
function xavier_uniform(dims...)
    bound = sqrt(1 / dims[2])
    return Float32.(rand(Uniform(-bound, bound), dims...))
end

struct Linear{F,S,T}
    W::S
    b::T
    σ::F
end

Linear(W, b) = Linear(W, b, identity)

function Linear(in::Integer, out::Integer, σ = identity;
                initW = xavier_uniform, initb = nothing)
    if initb == nothing
        bias_bound = 1 / sqrt(in)
        initb = (out) -> Float32.(rand(Uniform(-bias_bound, bias_bound), out))
    end
    return Linear(param(initW(out, in)), param(initb(out)), σ)
end

Flux.@treelike Linear

function (a::Linear)(x::AbstractArray)
    W, b, σ = a.W, a.b, a.σ
    σ.(W*x .+ b)
end

function Base.show(io::IO, l::Linear)
    print(io, "Linear(", size(l.W, 2), ", ", size(l.W, 1))
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ")")
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module Mustache

using Compat

include("utils.jl")
include("tokens.jl")
include("scanner.jl")
include("context.jl")
include("writer.jl")
include("parse.jl")

export @mt_str, @mt_mstr, render, render_from_file

"""

Macro to comile simple parsing outside of loops
use as mt"{{a}} and {{b}}", say

"""
macro mt_str(s)
    parse(s)
end

macro mt_mstr(s)
    parse(s)
end

"""

Render a set of tokens with a view, using optional `io` object to print or store.

Arguments
---------

* `io::IO`: Optional `IO` object.

* `tokens`: Either Mustache tokens, or a string to parse into tokens

* `view`: A view provides a context to look up unresolved symbols
  demarcated by mustache braces. A view may be specified by a
  dictionary, a module, a composite type, a vector, or keyword
  arguments.

"""
function render(io::IO, tokens::MustacheTokens, view)
    _writer = Writer()
    render(io, _writer, tokens, view)
end
function render(io::IO, tokens::MustacheTokens; kwargs...)
    d = Dict(kwargs) # [k => v for (k,v) in kwargs]
    render(io, tokens, d)
end

render(tokens::MustacheTokens, view) = sprint(io -> render(io, tokens, view))
render(tokens::MustacheTokens; kwargs...) = sprint(io -> render(io, tokens; kwargs...))

## Exported call without first parsing tokens via mt"literal"
##
## @param template a string containing the template for expansion
## @param view a Dict, Module, CompositeType, DataFrame holding variables for expansion
function render(io::IO, template::AbstractString, view)
    _writer = Writer()
    render(io, _writer, parse(template), view)
end
function render(io::IO, template::AbstractString; kwargs...)
    _writer = Writer()
    view = Dict(kwargs) # [k => v for (k,v) in kwargs]
    render(io, _writer, parse(template), view)
end
render(template::AbstractString, view) = sprint(io -> render(io, template, view))
render(template::AbstractString; kwargs...) = sprint(io -> render(io, template; kwargs...))

## Dict for storing parsed templates
TEMPLATES = Dict{AbstractString, MustacheTokens}()

## Load template from file
function template_from_file(filepath)
    f = open(readstring, filepath)
    tpl = parse(f)
#    f = open(filepath)
#    tpl = parse(readstring(f))
#    close(f)
    return tpl
end

"""

Renders a template from `filepath` and `view`. If it has seen the file
before then it finds the compiled `MustacheTokens` in `TEMPLATES` rather
than calling `parse` a second time.

"""
function render_from_file(filepath, view)
    if haskey(TEMPLATES, filepath)
        render(TEMPLATES[filepath], view)
    else
        try
            tpl = template_from_file(filepath)
            TEMPLATES[filepath] = tpl
            render(tpl, view)
        catch
            nothing
        end
    end
end
function render_from_file(filepath::AbstractString; kwargs...)
    d = Dict(kwargs) # [k => v for (k,v) in kwargs]
    render_from_file(filepath, d)
end

end


text:
"""
    QuadraticFit{T} <: UnivariateAlgorithm

The quadratic fit algorithm constructs a quadratic function from three points at each
iteration as an approximation to the objective. The analytical solution to the quadratic
model is examined as the new minimizer candidate to update the three points.
"""
struct QuadraticFit{T} <: UnivariateAlgorithm
    ϵ::T

    function QuadraticFit{T}(ϵ) where {T}
        0 < ϵ < 1 ? new{T}(ϵ) : error("`ϵ` must be in the interval (0, 1)")
    end
end

"""
    QuadraticFit([ϵ=0.01])

Initialize quadratic fit search algorithm with optional argument `ϵ` in the interval
(0, 1) to choose a midpoint as `mid = ϵ * lower + (1 - ϵ) * upper`.
"""
QuadraticFit(ϵ=0.5) = QuadraticFit{typeof(ϵ)}(ϵ)

function _optimize(f, lower::T, upper::T, alg::QuadraticFit;
        reltol, abstol, maxiter) where {T}
    ϵ = alg.ϵ
    mid = T(ϵ * lower + (1 - ϵ) * upper)
    y_lower, y_upper, y_mid = f(lower), f(upper), f(mid)
    converged = false
    iter = 0
    x = T(NaN)
    yx = (typeof(y_lower))(NaN)
    while iter < maxiter
        iter += 1
        x = T(0.5 * (y_lower * (mid^2 - upper^2) + y_mid * (upper^2 - lower^2) +
            y_upper * (lower^2 - mid^2)) /
            (y_lower * (mid - upper) + y_mid * (upper - lower) + y_upper * (lower - mid)))
        yx = f(x)
        x_tol = reltol * abs(x) + abstol
        if max(upper - x, x - lower) < 2x_tol
            converged = true
            break
        end
        if x > mid
            if yx > y_mid
                upper, y_upper = x, yx
            else
                lower, y_lower, mid, y_mid = mid, y_mid, x, yx
            end
        elseif x < mid
            if yx > y_mid
                lower, y_lower = x, yx
            else
                upper, y_upper, mid, y_mid = mid, y_mid, x, yx
            end
        end
    end
    return Solution(converged, iter, x, yx)
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, NetCDF
using YAXArrays.Datasets: open_mfdataset
CairoMakie.activate!() # HIDE
let
    dataset = open_mfdataset("./data/GPPdata_Beer_etal_2010_Science.nc")
    lon = dataset.lon.values
    lat = dataset.lat.values
    data = dataset.gpp.data[:,:]

    fig = Figure(resolution = (1250,700), fontsize = 25)
    ax = Axis(fig, xlabel= "lon",ylabel = "lat", xticksize=15, yticksize=15,
        xgridstyle=:dash, ygridstyle=:dash,)
    pltobj=heatmap!(ax,lon,lat,replace(data,-9999=>NaN),colormap = :Spectral_11)
    cbar  =Colorbar(fig,pltobj,label = "gpp",ticklabelsize = 25,tickalign =1)
    ax.xticks = -180:60:180
    ax.yticks = -90:30:90
    ax.xtickformat = "{:d}ᵒ"
    ax.ytickformat = "{:d}ᵒ"
    hidespines!(ax, :t, :r)
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "gppBeer2010.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["CairoMakie","NetCDF", "YAXArrays"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module Infrastructure

import ..Counter: Domain, API

using .Domain, .API

using Serialization

export add_to_file, read_from_file, create_new_invoice_nbr

function read_from_file(file::String)
    io = open(file, "r")

    r = []
    while !eof(io)
        push!(r, deserialize(io))
    end

    close(io)

    return r
end # read_from_file

function add_to_file(file::String, data::Array{T, 1} where T <: Any)
    io = open(file, "a+")

    [serialize(io, r) for r in data]

    close(io)
end # add_to_file

function create_new_invoice_nbr(file::String)
    value = 1000
    try
        value = last(read_from_file(file))
    catch e
        add_to_file(file, [value])
    end

    invoice_nbr = value += 1
    add_to_file(file, [invoice_nbr])
    return invoice_nbr
end # create_new_invoice_nbr

end # module


text:
function  mae(predictions, actualdata )
	# mean absolute error
	if (size(predictions,1)>1 && size(predictions,2)>1) || (size(actualdata,1)>1 && size(actualdata,2)>1)
		out = mean(abs(predictions-actualdata),1)
	else
		out = mean(abs(predictions-actualdata))
	end
	return out
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer