text:
mutable struct Generator
    next_node_id::NODE_ID
    next_scenario_id::SCENARIO_ID
end

function Generator()::Generator
    return Generator(0,0)
end

function _generate_node_id(gen::Generator)::NodeID
    id = NodeID(gen.next_node_id)
    gen.next_node_id += 1
    return id
end

function _generate_scenario_id(gen::Generator)::ScenarioID
    id = ScenarioID(gen.next_scenario_id)
    gen.next_scenario_id += 1
    return id
end

"""
Struct representing a node in a scenario tree.
"""
struct ScenarioNode
    id::NodeID # id of this node
    stage::StageID # stage of this node
    scenario_bundle::Set{ScenarioID} # scenarios that are indistiguishable
    parent::Union{Nothing,ScenarioNode}
    children::Set{ScenarioNode}
end

function stringify(node::ScenarioNode)
    node_str = "$(value(node.id)):{"
    for s in sort!(collect(node.scenario_bundle))
        node_str *= "$(value(s)),"
    end
    rstrip(node_str, ',')
    node_str *= "}"
    return node_str
end

Base.show(io::IO, sn::ScenarioNode) = print(io, "ScenarioNode($(sn.id), $(sn.stage), $(sn.scenario_bundle))")

function _add_child(parent::ScenarioNode,
                    child::ScenarioNode,
                    )::Nothing
    push!(parent.children, child)
    return
end

function _create_node(gen::Generator,
                      parent::Union{Nothing,ScenarioNode}
                      )::ScenarioNode
    nid = _generate_node_id(gen)
    stage = (parent==nothing ? StageID(1) : _increment(parent.stage))
    sn = ScenarioNode(nid, stage,
                      Set{ScenarioID}(),
                      parent,
                      Set{ScenarioNode}())
    if parent != nothing
        _add_child(parent, sn)
    end
    return sn
end

function id(node::ScenarioNode)::NodeID
    return node.id
end

function stage(node::ScenarioNode)::StageID
    return node.stage
end

"""
Struct representing the scenario structure of a stochastic program.

Can be built up by the user using the functions `add_node` and `add_leaf`.

**Constructor**

ScenarioTree()

Default constructor generates the root node of the tree. Can get the root node with `root`.
"""
struct ScenarioTree
    root::ScenarioNode
    tree_map::Dict{NodeID, ScenarioNode} # map from NodeID to tree node
    prob_map::Dict{ScenarioID, Float64}
    id_gen::Generator
end

function ScenarioTree(root_node::ScenarioNode,
                      gen::Generator,
                      )::ScenarioTree

    tree_map = Dict{NodeID, ScenarioNode}()
    prob_map = Dict{ScenarioID, Float64}()

    st = ScenarioTree(root_node,
                      tree_map,
                      prob_map,
                      gen)

    _add_node(st, root_node)
    return st
end

function ScenarioTree()::ScenarioTree
    gen = Generator()
    rn = _create_node(gen, nothing)
    st = ScenarioTree(rn, gen)
    return st
end

"""
    root(tree::ScenarioTree)

Return the root node of the given ScenarioTree
"""
function root(tree::ScenarioTree)::ScenarioNode
    return tree.root
end

function last_stage(tree::ScenarioTree)::StageID
    return maximum(getfield.(values(tree.tree_map), :stage))
end

function is_leaf(tree::ScenarioTree, vid::VariableID)::Bool
    return is_leaf(node(tree, vid.scenario, vid.stage))
end

function is_leaf(node::ScenarioNode)::Bool
    return length(node.children) == 0 ||
        (length(node.children) == 1 && is_leaf(first(node.children)))
end

function is_leaf(tree::ScenarioTree, nid::NodeID)::Bool
    return is_leaf(tree.tree_map[nid])
end

function is_leaf(tree::ScenarioTree, node::ScenarioNode)::Bool
    return is_leaf(tree.tree_map[node.id])
end

function _add_node(tree::ScenarioTree, node::ScenarioNode)::Nothing
    tree.tree_map[node.id] = node
    return
end

"""
    add_node(tree::ScenarioTree, parent::ScenarioNode)

Add a node to the ScenarioTree `tree` with parent node `parent`. Return the added node. If the node to add is a leaf, use `add_leaf` instead.
"""
function add_node(tree::ScenarioTree,
                  parent::ScenarioNode,
                  )::ScenarioNode
    new_node = _create_node(tree.id_gen, parent)
    _add_node(tree, new_node)
    return new_node
end

function _add_scenario_to_bundle(tree::ScenarioTree,
                                 nid::NodeID,
                                 scid::ScenarioID,
                                 )::Nothing
    push!(tree.tree_map[nid].scenario_bundle, scid)
    return
end

"""
    add_leaf(tree::ScenarioTree, parent::ScenarioNode, probability<:Real)

Add a leaf to the ScenarioTree `tree` with parent node `parent`. The probability of this scenario occuring is given by `probability`. Returns the ScenarioID representing the scenario.
"""
function add_leaf(tree::ScenarioTree,
                  parent::ScenarioNode,
                  probability::R,
                  )::ScenarioID where R <: Real

    if probability < 0.0 || probability > 1.0
        error("Invalid probability value: $probability")
    end

    leaf = add_node(tree, parent)
    scid = _assign_scenario_id(tree)
    tree.prob_map[scid] = probability

    node = leaf
    while node != nothing
        id = node.id
        _add_scenario_to_bundle(tree, id, scid)
        node = node.parent
    end

    return scid
end

_assign_scenario_id(tree::ScenarioTree)::ScenarioID = _generate_scenario_id(tree.id_gen)

function node(tree::ScenarioTree,
              nid::NodeID
              )::ScenarioNode
    return tree.tree_map[nid]
end

function node(tree::ScenarioTree,
              scid::ScenarioID,
              stid::StageID
              )::Union{Nothing,ScenarioNode}

    ret_node = nothing

    for node in values(tree.tree_map)

        if stage(node) == stid && scid in scenario_bundle(node)
            ret_node = node
            break
        end

    end

    return ret_node

end

function scenario_bundle(node::ScenarioNode)::Set{ScenarioID}
    return node.scenario_bundle
end

function scenario_bundle(tree::ScenarioTree, nid::NodeID)::Set{ScenarioID}
    return scenario_bundle(tree.tree_map[nid])
end

scenarios(tree::ScenarioTree)::Set{ScenarioID} = tree.root.scenario_bundle

function scenario_nodes(tree::ScenarioTree)
    return collect(values(tree_map))
end

function scenario_nodes(tree::ScenarioTree, scid::ScenarioID)
    nodes = Vector{ScenarioNode}()
    for node in values(tree.tree_map)
        if scid in scenario_bundle(node)
            push!(nodes, node)
        end
    end
    return nodes
end


text:
module DBInterface

"Database packages should subtype `DBInterface.Connection` which represents a connection to a database"
abstract type Connection end

"""
    DBInterface.connect(DB, args...; kw...) => DBInterface.Connection

Database packages should overload `DBInterface.connect` for a specific `DB` `DBInterface.Connection` subtype
that returns a valid, live database connection that can be queried against.
"""
function connect end

# Different `close!` signatures have their own docstrings.
function close! end

"""
    DBInterface.close!(conn::DBInterface.Connection)

Immediately closes a database connection so further queries cannot be processed.
"""
close!(conn::Connection)

"Database packages should provide a `DBInterface.Statement` subtype which represents a valid, prepared SQL statement that can be executed repeatedly"
abstract type Statement end

"""
    DBInterface.prepare(conn::DBInterface.Connection, sql::AbstractString) => DBInterface.Statement
    DBInterface.prepare(f::Function, sql::AbstractString) => DBInterface.Statement

Database packages should overload `DBInterface.prepare` for a specific `DBInterface.Connection` subtype, that validates and prepares
a SQL statement given as an `AbstractString` `sql` argument, and returns a `DBInterface.Statement` subtype. It is expected
that `DBInterface.Statement`s are only valid for the lifetime of the `DBInterface.Connection` object against which they are prepared.
For convenience, users may call `DBInterface.prepare(f::Function, sql)` which first calls `f()` to retrieve a valid `DBInterface.Connection`
before calling `DBInterface.prepare(conn, sql)`; this allows deferring connection retrieval and thus statement preparation until runtime,
which is often convenient when building applications.
"""
function prepare end

prepare(f::Function, sql::AbstractString) = prepare(f(), sql)

const PREPARED_STMTS = Dict{Symbol, Statement}()

"""
    DBInterface.@prepare f sql

Takes a `DBInterface.Connection`-retrieval function `f` and SQL statement `sql` and will return a prepared statement, via usage of `DBInterface.prepare`.
If the statement has already been prepared, it will be re-used (prepared statements are cached).
"""
macro prepare(getDB, sql)
    key = gensym()
    return quote
        get!(DBInterface.PREPARED_STMTS, $(QuoteNode(key))) do
            DBInterface.prepare($(esc(getDB)), $sql)
        end
    end
end

"Any object that iterates \"rows\", which are objects that are property-accessible and indexable. See `DBInterface.execute` for more details on fetching query results."
abstract type Cursor end

"""
    DBInterface.execute(conn::DBInterface.Connection, sql::AbstractString, [params]) => DBInterface.Cursor
    DBInterface.execute(stmt::DBInterface.Statement, [params]) => DBInterface.Cursor

Database packages should overload `DBInterface.execute` for a valid, prepared `DBInterface.Statement` subtype (the first method
signature is defined in DBInterface.jl using `DBInterface.prepare`), which takes an optional `params` argument, which should be
an indexable collection (`Vector` or `Tuple`) for positional parameters, or a `NamedTuple` for named parameters.
`DBInterface.execute` should return a valid `DBInterface.Cursor` object, which is any iterator of "rows",
which themselves must be property-accessible (i.e. implement `propertynames` and `getproperty` for value access by name),
and indexable (i.e. implement `length` and `getindex` for value access by index). These "result" objects do not need
to subtype `DBInterface.Cursor` explicitly as long as they satisfy the interface. For DDL/DML SQL statements, which typically
do not return results, an iterator is still expected to be returned that just iterates `nothing`, i.e. an "empty" iterator.

Note that `DBInterface.execute` returns ***a single*** `DBInterface.Cursor`, which represents a single resultset from the database.
For use-cases involving multiple resultsets from a single query, see `DBInterface.executemultiple`.
"""
function execute end

execute(conn::Connection, sql::AbstractString, params=()) = execute(prepare(conn, sql), params)

struct LazyIndex{T} <: AbstractVector{Any}
    x::T
    i::Int
end

Base.IndexStyle(::Type{<:LazyIndex}) = Base.IndexLinear()
Base.IteratorSize(::Type{<:LazyIndex}) = Base.HasLength()
Base.size(x::LazyIndex) = (length(x.x),)
Base.getindex(x::LazyIndex, i::Int) = x.x[i][x.i]

"""
    DBInterface.executemany(conn::DBInterface.Connection, sql::AbstractString, [params]) => Nothing
    DBInterface.executemany(stmt::DBInterface.Statement, [params]) => Nothing

Similar in usage to `DBInterface.execute`, but allows passing multiple sets of parameters to be executed in sequence.
`params`, like for `DBInterface.execute`, should be an indexable collection (`Vector` or `Tuple`) or `NamedTuple`, but instead
of a single scalar value per parameter, an indexable collection should be passed for each parameter. By default, each set of
parameters will be looped over and `DBInterface.execute` will be called for each. Note that no result sets or cursors are returned
for any execution, so the usage is mainly intended for bulk INSERT statements.
"""
function executemany(stmt::Statement, params=())
    if !isempty(params)
        param = params[1]
        len = length(param)
        all(x -> length(x) == len, params) || throw(ParameterError("parameters provided to `DBInterface.executemany!` do not all have the same number of parameters"))
        for i = 1:len
            xargs = LazyIndex(params, i)
            execute(stmt, xargs)
        end
    else
        execute(stmt)
    end
    return
end

executemany(conn::Connection, sql::AbstractString, params=()) = executemany(prepare(conn, sql), params)

"""
    DBInterface.executemultiple(conn::DBInterface.Connection, sql::AbstractString, [params]) => Cursor-iterator
    DBInterface.executemultiple(stmt::DBInterface.Statement, [params]) => Cursor-iterator

Some databases allow returning multiple resultsets from a "single" query (typically semi-colon (`;`) separated statements, or from calling stored procedures).
This function takes the exact same arguments as `DBInterface.execute`, but instead of returning a single `Cursor`, it returns an iterator of `Cursor`s.
This function defines a generic fallback that just returns `(DBInterface.execute(stmt, params),)`, a length-1 tuple for a single `Cursor` resultset.
"""
function executemultiple end

executemultiple(stmt::Statement, params=()) = (DBInterface.execute(stmt, params),)
executemultiple(conn::Connection, sql::AbstractString, params=()) = executemultiple(prepare(conn, sql), params)

"""
    DBInterface.close!(stmt::DBInterface.Statement)

Close a prepared statement so further queries cannot be executed.
"""
close!(stmt::Statement)

"""
    DBInterface.lastrowid(x::Cursor) => Int

If supported by the specific database cursor, returns the last inserted row id after executing an INSERT statement.
"""
function lastrowid end

"""
    DBInterface.close!(x::Cursor) => Nothing

Immediately close a resultset cursor. Database packages should overload for the provided resultset `Cursor` object.
"""
close!(x::Cursor)

# exception handling
"Error for signaling that parameters are used inconsistently or incorrectly."
struct ParameterError <: Exception
    msg::String
end

"Standard warning object for various database operations"
struct Warning
    msg::String
end

"Fallback, generic error object for database operations"
struct Error <: Exception
    msg::String
end

end # module

ext: jl


text:
__precompile__()

module RandomFerns

using FunctionalData

export RandFern, getcell
export RandFerns, predict, predict!, votes!

type RandFern
    base::Array{Float32,2}
    splits::Array{Float32,2}
    cells::Array{Array{Float32,2},1}
    cellinds::Dict{Int,Int}
    C::Int
    proj::Array{Float32,2}
    zeros::Array{Float32,2}
end

function RandFern(data, targets; regression = false, nsplits = 10)
    assert(size(data,2)==size(targets,2))
    base = randn(nsplits,size(data,1))
    proj = base*data
    mi = minimum(proj,2)
    ma = maximum(proj,2)
    splits = col(mi + (ma - mi).*rand(nsplits,1))
    binary = proj .< vec(splits)
    cells = Any[]
    cellinds = Dict{Int,Int}()
    counts = Dict{Int,Int}()
    if regression
        C = size(targets,1)
    else
        C = round(Int,maximum(targets))
        labelcounts = [sum(targets.==i) for i in 1:C]
    end
    for i = 1:size(binary,2)
        ind = bin2ind(binary[:,i])
        if !haskey(cellinds,ind)
            push!(cells, zeros(Float32,C,1))
            cellinds[ind] = length(cells)
        end
        if regression
            cells[cellinds[ind]] += targets[:,i]
            counts[ind] = get(counts, ind, 0) + 1
        else
            c = round(Int,targets[i])
            cells[cellinds[ind]][c] += 1/labelcounts[c]
        end
    end
    if regression
        for i = keys(cellinds)
            cells[cellinds[i]] /= counts[i]
        end
    end
    RandFern(base, convert(Array{Float32,2},splits), cells, cellinds, C, zeros(Float32,nsplits,1), zeros(Float32,C,1))
end

function getcell(a::RandFern, data::Array{Float32,2})
    A_mul_B!(a.proj, a.base, data)
    ind = 1
    n = length(a.splits)
    for i = 1:n
        if a.proj[i] < a.splits[i]
            ind += 2^(n-i)
        end
    end
    haskey(a.cellinds,ind) ? (a.cells[a.cellinds[ind]]::Array{Float32,2},true) : (a.zeros, false)
end

bin2ind(a) = (local r = 1; for i = 1:length(a) r += 2^(length(a)-i)*a[i] end; r)

type RandFerns
    ferns::Array{RandFern}
    regression::Bool
    buf::Array{Float32,2}
end

function RandFerns(data, targets; nferns = 32, nsplits = 10, regression = false)
    ferns = [RandFern(data, targets; nsplits = nsplits, regression = regression) for i in 1:nferns]
    RandFerns(ferns, regression, zeros(Float32,ferns[1].C,1))
end

predict{T<:Real}(a::RandFerns, data::Array{T,2}) = @p map convert(Array{Float32,2},data) x->predict!(a,x)

function predict!(a::RandFerns, data::Array{Float32,2})
    assert(size(data,2) == 1)
    a.regression ? votes!(a,data) : indmax(votes!(a, data))
end

function votes!(a::RandFerns, data)
    nfound = 0
    a.buf[:] = 0
    for i = 1:length(a.ferns)
        temp, found = getcell(a.ferns[i], data)
        nfound += found
        for j = 1:length(a.buf)
            a.buf[j] += temp[j]
        end
    end
    if a.regression
        for i = 1:length(a.buf)
            a.buf[i] /= nfound
        end
    end
    a.buf
end

end


text:
@testset "floydwarshall" begin
    A = load_matrix_network("all_shortest_paths_example")
    #A.nzval += abs(minimum(A)) + 1 # remove negative edges
    nzvals = nonzeros(A)
    val = abs(minimum(A)) + 1
    for ii=1:length(nzvals)
        nzvals[ii] += val
    end
    m = size(A,1)
    D2 = zeros(Float64,m,m)
    P2 = zeros(Float64,m,m)
    for ii=1:m
        (d,p) = dijkstra(A,ii)
        D2[ii,:] = d
        P2[ii,:] = p
    end

    (D,P) = floydwarshall(A)

    @test D == D2
    @test P == P2

    (D,P) = floydwarshall(MatrixNetwork(A))

    @test D == D2
    @test P == P2
end


text:
"""
    build_julia_model_object(path_to_model_file::String, path_to_output_dir::String;
        defaults_file_name::String="Defaults.toml", model_type::Symbol=:static)::VLResult

Build a static or dynamic data object in Julia language. If a TOML file with default parameter values is not provided by the user, a `Defaults.toml` file is generated automatically. Additionally, the user can edit the generated `Defaults.toml` file with their own values. After generating code, if a directory already exists at the user specified location, it can be deleted or backed-up before new code is written based on user input.

Input arguments:
`path_to_model_file::String` - user-specified path to read the Network.vff file from as well as the Defaults.toml file (if it already exists)
`path_to_output_dir::String` - path to where generated code will be written
`defaults_file_name::String` - name of the TOML file with default parameter values (optional).
`model_type::Symbol` - type of model to generate (static/ dynamic) (optional).

Output arguments:
`model_object::VLResult{VLJuliaModelObject}` - abstract Julia data object holding information for generating model code.

"""
function build_julia_model_object(path_to_model_file::String, path_to_output_dir::String;
    defaults_file_name::String="Defaults.toml", model_type::Symbol=:static)::VLResult

    # Check: is path_to_model_file legit?
    check_result = is_file_path_ok(path_to_model_file)
    if (isnothing(check_result.value) == false)
        return check_result
    end

    # Check: do we have a trailing slash on the path_to_output_dir?
    last_char = path_to_output_dir[end]
    if (last_char == '/')
        path_to_output_dir = path_to_output_dir[1:(end-1)]
    end

    # Check: is path_to_output_dir legit?
    # before we get too far along, we need to check if the user already has code in the location that we want to generate at -
    # if they do, then move it -
    if (isdir(path_to_output_dir) == true)

        # let the user know that we found code -
        @info "Ooops! We found some code where you wanted to generate your model code."

        # ask the user: should we move or delete?
        user_input_result = _request_user_input("Overwrite the existing project files [Y]/N ? ")
        if (lowercase(user_input_result) == "no" || lowercase(user_input_result) == "n")

            @info "Got it. Backing up the existing files ..."

            # ok, looks like we may have a conflict - mv the offending code
            if (move_existing_project_at_path(path_to_output_dir) == false)

                # Something happend ... the world is ending ...
                error = ArgumentError("automatic directory conflict resolution failed. Unable to move existing directory $(path_to_output_dir)")
                return VLResult(error)
            end
        else
            @info "Ok! Overwriting the existing files ..."
        end
    end

    # Check: do we have a defaults file?
    path_to_defaults_file = joinpath(dirname(path_to_model_file), defaults_file_name)
    if (isfile(path_to_defaults_file) == false)

        # ok, we do not have a defaults file, so lets build one, and let the user know -
        generate_defaults_result = generate_default_project_file(path_to_defaults_file)
        if (isnothing(generate_defaults_result.value) == false)
            return generate_defaults_result
        end

        # let the user know that we have generated a new Defaults.toml file -
        @info "We generated Defaults.toml using the generate_default_project_file function, see the help system in the REPL"
    end

    # build the model wrapper -
    model_object = VLJuliaModelObject(path_to_model_file, path_to_output_dir;
        defaults_file_path=path_to_defaults_file, model_type=model_type)

    # return -
    return VLResult{VLJuliaModelObject}(model_object)
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
using Flux, DiffEqFlux, OrdinaryDiffEq, StatsBase, RecursiveArrayTools
using Distributions

const tspan = (0.0,25)
const RANGE = (Float32(-8.),Float32(8.))
const BS = 200

target(u) = u.^3

function gen_data(batchsize,target_fun)
    x = Float32.(rand(Uniform(RANGE...,),batchsize))'|>collect
    return x,target(x)
end
data = Iterators.repeated(gen_data(BS,target), 10000)
u0_test,y0_test= first(data)

dudt = Chain(Dense(1,20,tanh),Dense(20,1))

function n_ode(u0)
    neural_ode(dudt,u0,tspan,Tsit5(),
               save_start=false,
               save_everystep=false,
               reltol=1e-5,abstol=1e-12)
end
ps = Flux.params(dudt)
n_ode(u0_test)

loss(x,y) = mean(abs.(n_ode(x).-y))

opt = ADAM(0.1)
cb = function () #callback function to observe training
    tx,ty = gen_data(BS,target)
    display(loss(tx,ty))
end

n_ode(u0_test)

g = Tracker.gradient(()->loss(u0_test,y0_test),ps)

g[dudt[1].W]

[print(g) for g in Tracker.gradient(()->loss(u0_test,y0_test),ps)] # no gradients

# Display the ODE with the initial parameter values.
cb()

Flux.train!(loss, ps, data, opt, cb = cb)

max_stars_repo_name: jessebett/model-zoo

max_issues_repo_name: jessebett/model-zoo

max_forks_repo_name: jessebett/model-zoo


text:
# ---
# title: 830. Positions of Large Groups
# id: problem830
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/positions-of-large-groups/description/>
# hidden: true
# ---
#
# In a string `s` of lowercase letters, these letters form consecutive groups of
# the same character.
#
# For example, a string like `s = "abbxxxxzyy"` has the groups `"a"`, `"bb"`,
# `"xxxx"`, `"z"`, and `"yy"`.
#
# A group is identified by an interval `[start, end]`, where `start` and `end`
# denote the start and end indices (inclusive) of the group. In the above
# example, `"xxxx"` has the interval `[3,6]`.
#
# A group is considered  **large**  if it has 3 or more characters.
#
# Return  _the intervals of every **large** group sorted in  **increasing order
# by start index**_.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abbxxxxzzy"
#     Output: [[3,6]]
#     **Explanation** : "xxxx" is the only large group with start index 3 and end index 6.
#
#
# **Example 2:**
#
#
#
#     Input: s = "abc"
#     Output: []
#     **Explanation** : We have groups "a", "b", and "c", none of which are large groups.
#
#
# **Example 3:**
#
#
#
#     Input: s = "abcdddeeeeaabbbcd"
#     Output: [[3,5],[6,9],[12,14]]
#     **Explanation** : The large groups are "ddd", "eeee", and "bbb".
#
#
# **Example 4:**
#
#
#
#     Input: s = "aba"
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 1000`
#   * `s` contains lower-case English letters only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using TensorNetworkEvolve
using Test

@testset "peps.jl" begin
    include("peps.jl")
end

@testset "tebd.jl" begin
    include("tebd.jl")
end

@testset "TensorAD.jl" begin
    include("TensorAD/TensorAD.jl")
end

@testset "timeevolve.jl" begin
    include("timeevolve.jl")
end


text:
"""
    solve!(calc::Calc)

Solves the object `sim::Calc`.

See also: [`Calc`](@ref)
"""
function solve!(calc::Calc)
    N = length(calc.λ)
    ψₜ = Array{Complex{Float64}}(undef, calc.box.Nₜ, calc.box.Nₓ, N+1)
    if calc.seed == "0"
        ψₜ[:,:,1] .= zeros(calc.box.Nₜ, calc.box.Nₓ)
    elseif calc.seed == "exp"
        ψₜ[:,:,1] .= exp.(im*calc.box.x').*ones(calc.box.Nₜ, calc.box.Nₓ)
    elseif calc.seed == "dn"
        ψₜ[:,:,1] .= exp.(im*calc.box.x'.*(1-calc.m/2)).*dn.(calc.box.t,calc.m)
    elseif calc.seed == "cn"
        ψₜ[:,:,1] .= exp.(im*calc.box.x'.*(calc.m - 1/2)).*sqrt(calc.m).*cn.(calc.box.t,calc.m)
    end
    calc_rs(calc, N, 1, ψₜ)
    calc.ψ .= ψₜ[:,:,N+1]
    calc.ψ̃ .= fftshift(fft(calc.ψ, 1),1)./calc.box.Nₜ

    return nothing
end

"""
    calc_rs(c:Calc, n, p, ψₜ)

Calculate the lax pair generating functions recursively. See the paper for more details.

See also: [`solve!`](@ref), [`Calc`](@ref)
"""
@memoize function calc_rs(c::Calc, n, p, ψₜ)
    @info "Calculating Lax pair generating functions rₙₚ(x,t) and sₙₚ(x,t) for (n,p) = ($n,$p)"
    # Base case
    if n == 1
        if c.seed == "0"
            # Set up x and t in the proper way
            t = c.box.t .- c.tₛ[p]
            x = c.box.x' .- c.xₛ[p]
            # Compute r and s directly
            α = c.f[:α]
            γ = c.f[:γ]
            δ = c.f[:δ]
            d = (1 - 4*α*c.λ[p] - 8*γ*c.λ[p]^2 + 16*δ*c.λ[p]^3)
            rf = exp.(+im.*(c.λ[p] .* t .+ c.λ[p]^2 * d .* x .- π/4))
            sf = exp.(-im.*(c.λ[p] .* t .+ c.λ[p]^2 * d .* x .- π/4))
        elseif c.seed == "exp"
            # Set up x and t in the proper way
            t = c.box.t .- c.tₛ[p]
            x = c.box.x' .- c.xₛ[p]
            # Compute r and s directly
            α = c.f[:α]
            γ = c.f[:γ]
            δ = c.f[:δ]
            d = 2*(α+3*δ) + (1+4*γ)*c.λ[p] - 4*(α+2*δ)*c.λ[p]^2 - 8*γ*c.λ[p]^3 + 16*δ*c.λ[p]^4
            A = +c.χ[p].+0.5.*(c.Ω[p].*t .+ c.Ω[p]*d.*x) .- π/4
            B = -c.χ[p].+0.5.*(c.Ω[p].*t .+ c.Ω[p]*d.*x) .- π/4
            rf = 2im.*exp.(-im.*c.box.x'./2) .* sin.(A)
            sf = 2  .*exp.(+im.*c.box.x'./2) .* cos.(B)
        elseif c.seed == "dn"
            if c.f[:α] != 0.0 || c.f[:γ] !=0 || c.f[:δ] != 0
                @error "Computing DT with dn seed for extended NLS isnot currently supported"
            end
            # Allocate empty arrays to hold lax pair generating functions
            rf = Array{Complex{Float64}}(undef, c.box.Nₜ, c.box.Nₓ)
            sf = Array{Complex{Float64}}(undef, c.box.Nₜ, c.box.Nₓ)
            # Set up x in the proper way
            # We currently do not suppport t shifts for this seed
            x = c.box.x' .- c.xₛ[p]

            # Set up a funnction for the ODE
            function dn_ab!(du,u,p,t)
                du[1] = +im*p[1].*u[1] .+ im.*u[2].*dn.(t,real(p[2]))
                du[2] = -im*p[1].*u[2] .+ im.*u[1].*dn.(t,real(p[2]))
            end
            # And the initial condition
            function ab_dn_t0(x)
                A = +c.χ[p] .+ 0.5.*c.Ω[p]*c.λ[p].*(x - c.xₛ[p]) .- π/4
                B = -c.χ[p] .+ 0.5.*c.Ω[p]*c.λ[p].*(x - c.xₛ[p]) .- π/4
                return [2*im*sin(A); 2*cos(B)]
            end

            # Prepare the ODE intnegrator
            tspan = (0.0,abs(minimum(c.box.t)))
            u0 = ab_dn_t0(c.box.x[1])
            prob = ODEProblem(dn_ab!, u0, tspan, [c.λ[p], c.m])
            dt = c.box.t[2]-c.box.t[1]
            integrator = init(prob, Tsit5(), saveat=abs.(c.box.t[c.box.Nₜ÷2+1:-1:1]), dt = dt, adaptive=false)

            # Loop over x and solve the ODE for each initial condition
            # Should do this using the ensemble interface in the future
            for i in eachindex(c.box.x)
                # Set up the new initial condition
                u0 = ab_dn_t0(c.box.x[i])
                reinit!(integrator, u0)
                # Solve
                DiffEqBase.solve!(integrator)
                # Save results
                rf[c.box.Nₜ÷2+1:-1:1, i] .= integrator.sol[1, :].*exp.(+im*(c.box.x[i])./4*(c.m-2))
                sf[c.box.Nₜ÷2+1:-1:1, i] .= integrator.sol[2, :].*exp.(-im*(c.box.x[i])./4*(c.m-2))
                # Mirror x -> x
                rf[c.box.Nₜ÷2+2:end, i] .= rf[c.box.Nₜ÷2:-1:2, i]
                sf[c.box.Nₜ÷2+2:end, i] .= sf[c.box.Nₜ÷2:-1:2, i]
            end
        elseif c.seed == "cn"
            if c.f[:α] != 0.0 || c.f[:γ] !=0 || c.f[:δ] != 0
                @error "Computing DT with cn seed for extended NLS isnot currently supported"
            end
            # Allocate empty arrays to hold lax pair generating functions
            rf = Array{Complex{Float64}}(undef, c.box.Nₜ, c.box.Nₓ)
            sf = Array{Complex{Float64}}(undef, c.box.Nₜ, c.box.Nₓ)
            # Set up x in the proper way
            # We currently do not suppport t shifts for this seed
            x = c.box.x' .- c.xₛ[p]

            # Set up a funnction for the ODE
            function cn_ab!(du,u,p,t)
                du[1] = +im*p[1].*u[1] .+ im.*u[2]*sqrt(p[2]).*cn.(t,real(p[2]))
                du[2] = -im*p[1].*u[2] .+ im.*u[1]*sqrt(p[2]).*cn.(t,real(p[2]))
            end
            # And the initial condition
            function ab_cn_t0(x)
                A = +c.χ[p] .+ 0.5.*c.Ω[p]*c.λ[p].*(x - c.xₛ[p]) .- π/4
                B = -c.χ[p] .+ 0.5.*c.Ω[p]*c.λ[p].*(x - c.xₛ[p]) .- π/4
                return [2*im*sin(A); 2*cos(B)]
            end

            # Prepare the ODE intnegrator
            tspan = (0.0,abs(minimum(c.box.t)))
            u0 = ab_cn_t0(c.box.x[1])
            prob = ODEProblem(cn_ab!, u0, tspan, [c.λ[p], c.m])
            dt = c.box.t[2] - c.box.t[1]
            integrator = init(prob, Tsit5(), saveat=abs.(c.box.t[c.box.Nₜ÷2+1:-1:1]), dt=dt, adaptive=false)

            # Loop over x and solve the ODE for each initial condition
            # Should do this using the ensemble interface in the future
            for i in eachindex(c.box.x)
                # Set up the new initial condition
                u0 = ab_cn_t0(c.box.x[i])
                reinit!(integrator, u0)
                # Solve
                DiffEqBase.solve!(integrator)
                # Save results
                rf[c.box.Nₜ÷2+1:-1:1, i] .= integrator.sol[1, :].*exp.(+im*(c.box.x[i])./4*(2*c.m-1))
                sf[c.box.Nₜ÷2+1:-1:1, i] .= integrator.sol[2, :].*exp.(-im*(c.box.x[i])./4*(2*c.m-1))
                # Mirror x -> x
                rf[c.box.Nₜ÷2+2:end, i] .= rf[c.box.Nₜ÷2:-1:2, i]
                sf[c.box.Nₜ÷2+2:end, i] .= sf[c.box.Nₜ÷2:-1:2, i]
            end
        end

        # Compute ψ₁ when p = 1
        if  p == 1
            ψₜ[:,:,n+1] .= ψₜ[:,:,1] .+ (2*(c.λ[n]' - c.λ[n]).*sf.*conj.(rf)) ./
                              (abs2.(rf) + abs2.(sf))
        end
    else
        # Recursion when n != 1
        # Calculate r_{n_1, 1} and r_{n-1, p} and similarly for s
        r1, s1 = calc_rs(c, n-1, 1, ψₜ)
        r2, s2 = calc_rs(c, n-1, p+1, ψₜ)

        # Apply the DT equations
        rf = ((c.λ[n-1]'  - c.λ[n-1] ) .* conj.(s1) .* r1 .* s2 +
              (c.λ[p+n-1] - c.λ[n-1] ) .* abs2.(r1) .* r2       +
              (c.λ[p+n-1] - c.λ[n-1]') .* abs2.(s1) .* r2        ) ./
              (abs.(r1).^2 + abs.(s1).^2);
        sf = ((c.λ[n-1]'  - c.λ[n-1] ) .*       s1   .*conj(r1) .* r2 +
              (c.λ[p+n-1] - c.λ[n-1] ) .*  abs2.(s1) .*     s2        +
              (c.λ[p+n-1] - c.λ[n-1]') .*  abs2.(r1) .*     s2         ) ./
              (abs.(r1).^2 + abs.(s1).^2);
        if p == 1
            # Compute ψₙ when p = 1
            ψₜ[:,:,n+1] .= ψₜ[:,:,n] + (2*(c.λ[n]' - c.λ[n]).*sf.*conj.(rf)) ./
                                      (abs2.(rf) + abs2.(sf))
        end
    end

    return rf, sf
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_native_jl_rho0_015(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_native_jl_rho0_015!(A_copy)
end

@inline function exp_native_jl_rho0_015!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=4
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 Ua U V Z X P
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing Ua = x*I+x*A2
    coeff1=60.0
    coeff2=1.0
    memslots3 .= coeff2.*memslots2
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing U with operation: mult
    mul!(memslots4,memslots3,memslots1)
    # Deallocating Ua in slot 3
    # Deallocating A in slot 1
    # Computing V = x*I+x*A2
    coeff1=120.0
    coeff2=12.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots4
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots4
    # Deallocating U in slot 4
    # Computing P with operation: ldiv
    memslots3 .=memslots1\memslots2
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    return memslots3 # Returning P
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module SpringCollab2020StrawberryIgnoringLighting

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Strawberry (Ignore Lighting) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        Maple.Strawberry,
        "point",
        Dict{String, Any}(
            "SpringCollab2020_ignoreLighting" => true
        )
    )
)

end


text:
"""
    ExponentialLikelihood(l=exp)

Exponential likelihood with scale given by `l(f)`.

```math
    p(y|f) = \\operatorname{Exponential}(y | l(f))
```
"""
struct ExponentialLikelihood{Tl<:AbstractLink} <: AbstractLikelihood
    invlink::Tl
end

ExponentialLikelihood(l=exp) = ExponentialLikelihood(link(l))

(l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f))

(l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))

function expected_loglikelihood(
    ::AnalyticExpectation,
    ::ExponentialLikelihood{ExpLink},
    q_f::AbstractVector{<:Normal},
    y::AbstractVector{<:Real},
)
    f_μ = mean.(q_f)
    return sum(-f_μ - y .* exp.((var.(q_f) / 2) .- f_μ))
end

default_expectation_method(::ExponentialLikelihood{ExpLink}) = AnalyticExpectation()


text:
struct Keyboard
  pressed::Array{UInt32, 1}
end

Keyboard() = Keyboard(Array{Any, 1}())

function Base.getproperty(k::Keyboard, s::Symbol)
  s = Symbol(uppercase(string(s)))
  if isdefined(GameZero.Keys, s)
    return UInt32(getfield(GameZero.Keys, s)) in getfield(k, :pressed)
  end
  if isdefined(GameZero.Keymods, s)
    return UInt32(getfield(GameZero.Keymods, s)) in getfield(k, :pressed)
  end
  return false
end

Base.push!(k::Keyboard, item) = push!(getfield(k, :pressed), item)
function Base.delete!(k::Keyboard, item)
  a = getfield(k, :pressed)
  deleteat!(a, a.==item)
end


text:
export to_device, to_host

using CUDAapi
using Flux
using Knet
using Adapt
import Adapt: adapt, adapt_storage

to_host(x) = to_device(Val(:cpu), x)
to_device(::Val{:cpu}, x) = x  # cpu(x) is not very efficient! So by default we do nothing here.

if has_cuda()
    using CuArrays

    to_device(::Val{:Zygote_gpu}, x) = Flux.fmap(a -> adapt(CuArray{Float32}, a), x)
    to_device(::Val{:cpu}, x::Union{KnetArray, CuArray}) = adapt(Array, x)

    to_device(::Val{:Knet_gpu}, x) = adapt(KnetArray, x)
    to_device(::Val{:Knet_gpu}, x::SubArray) = KnetArray(x)
end

adapt_storage(T::Type{<:KnetArray}, x::AbstractArray{<:Real}) = T(x)
adapt_storage(T::Type{<:KnetArray}, x::Knet.Param{<:AbstractArray}) = Knet.Param(T(value(x)))


text:
export Matsuno_fast,Matsuno

"""
	Matsuno_fast(param;dealias,label,verbose)

Same as `Matsuno`, but faster.
"""
mutable struct Matsuno_fast <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function

    function Matsuno_fast(param::NamedTuple;
							dealias=false, label = "Matsuno", verbose=true )

			# Set up
			ϵ 	= param.ϵ

			if verbose  # Print information
				info = "Build the Matsuno model.\n"
				info *= "Steepness parameter ϵ=$ϵ (infinite depth case).\n"
				if dealias == true || dealias == 1
					info *= "Dealiasing with Orszag’s 3/2 rule. "
				else
					info *= "No dealiasing. "
				end
				info *= "\nShut me up with keyword argument `verbose = false`."
				@info info
				@warn "The velocity is consistent with the \
				derivative of the trace of the velocity potential \
				used for water waves only when they are null."
			end

		# Pre-allocate useful data
		ϵ        = param.ϵ
        mesh     = Mesh(param)
        x        = mesh.x
		k 		 = mesh.k
        Γ        = abs.(k)
        Dx       =  1im * k        	# Differentiation
        H        = -1im * sign.(k) 	# Hilbert transform
		if dealias == true || dealias == 1
			Π⅔    = Γ .< mesh.kmax * 2/3    # Dealiasing low-pass filter
		else
			Π⅔    = zero(Γ) .+ 1     		# No dealiasing (Π⅔=Id)
		end

        hnew = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₀ = zeros(Complex{Float64}, mesh.N)
        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

        Px  = plan_fft(hnew; flags = FFTW.MEASURE)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		    for i in eachindex(hnew)
		        hnew[i] = Γ[i] * U[i,1]
		    end

		    ldiv!(unew, Px, hnew )

		    for i in eachindex(hnew)
		        hnew[i] = Dx[i] * U[i,1]
		    end

		    ldiv!(I₁, Px, hnew)

		    unew  .*= I₁

		    mul!(I₁, Px, unew)

		    I₁  .*= ϵ .* Π⅔
		    I₁  .-= hnew

		    ldiv!(hnew, Px, view(U,:,1))
		    ldiv!(unew, Px, view(U,:,2))

		    I₂    .= hnew .* unew

		    mul!(I₃, Px, I₂)

		    I₃    .*= Dx

		    for i in eachindex(H)
		        U[i,1]  = H[i] * U[i,2]
		        I₀[i] = Γ[i] * U[i,2]
		    end

		    ldiv!(I₂, Px, I₀)

		    I₂    .*= hnew

		    mul!(hnew, Px, I₂)

		    hnew  .*= H
		    I₃    .+= hnew
		    I₃    .*= ϵ .* Π⅔

		    for i in eachindex(I₃)
		        U[i,1] -= I₃[i]
		    end

		    I₃    .=  unew.^2

		    mul!(unew, Px, I₃)

		    unew  .*= Dx
		    unew  .*= ϵ/2 .* Π⅔
		    I₁    .-= unew

		    for i in eachindex(I₁)
		        U[i,2] =  I₁[i]
		    end

		end

		# Build raw data from physical data.
		# Discrete Fourier transform with, possibly, dealiasing.
		# This function is correct only when the initial data for v is zero.
		function mapto(data::InitialData)

		    [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]

		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the velocity variable.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
		    real(ifft(view(U,:,1))),real(ifft(view(U,:,2)))
		end

		new(label, f!, mapto, mapfro )
    end
end

"""
	Matsuno(param;dealias,label,verbose)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the quadratic deep-water model proposed by [Matsuno](https://doi.org/10.1103/PhysRevLett.69.609).

# Arguments
`param` is of type `NamedTuple` and must contain
- the dimensionless parameters `ϵ` (nonlinearity);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`.

## Optional keyword arguments
- `dealias`: dealiasing with `1/3` Orlicz rule if `true` or no dealiasing if `false` (by default);
- `label`: a label for future references (default is `"Matsuno"`);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `DeepQuadratic.f!` to be called in explicit time-integration solvers;
2. a function `DeepQuadratic.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `DeepQuadratic.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is a velocity variable which is *not* the derivative of the trace of the velocity potential (if not null).

"""
mutable struct Matsuno <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function

    function Matsuno(param::NamedTuple ;
						dealias=false, label = "Matsuno", verbose=true )

		# Set up
		ϵ 	= param.ϵ

		if verbose  # Print information
			info = "Build the Matsuno model.\n"
			info *= "Steepness parameter ϵ=$ϵ (infinite depth case).\n"
			if dealias == true || dealias == 1
				info *= "Dealiasing with Orszag’s 3/2 rule. "
			else
				info *= "No dealiasing. "
			end
			info *= "\nShut me up with keyword argument `verbose = false`."
			@info info
			@warn "The velocity is consistent with the \
			derivative of the trace of the velocity potential \
			used for water waves only when they are null."
		end

		# Pre-allocate useful data
		mesh = Mesh(param)
		x   = mesh.x
		k 	= mesh.k
        Γ 	= abs.(k)
    	∂ₓ	=  1im * k            # Differentiation
        H 	= -1im * sign.(k)     # Hilbert transform
		if dealias == true || dealias == 1
			Π⅔    = Γ .< mesh.kmax * 2/3     # Dealiasing low-pass filter
		else
			Π⅔    = zero(Γ) .+ 1     # Dealiasing low-pass filter
		end

        hnew = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		   hnew .= ifft(U[:,1])
		   unew .= ifft(U[:,2])
		   I₃ .= fft(ifft((∂ₓ).*U[:,1]).*ifft((Γ).*U[:,1]))
		   I₁ .= H.*U[:,2].-ϵ*Π⅔.*(H.*fft(hnew.*ifft(Γ.*U[:,2])).+∂ₓ.*fft(hnew.*unew))
		   I₂ .= -(∂ₓ.*U[:,1])+ϵ*Π⅔.*(I₃-∂ₓ.*fft(unew.^2)/2)
		   #
		   U[:,1] .= I₁
		   U[:,2] .= I₂

		end

		# Build raw data from physical data.
		# Discrete Fourier transform with, possibly, dealiasing.
		# This function is correct only when the initial data for v is zero.
		function mapto(data::InitialData)

		    [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]

		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the velocity variable.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
		    real(ifft(view(U,:,1))),real(ifft(view(U,:,2)))
		end

		new(label, f!, mapto, mapfro )
    end
end


text:
# Implementation of sensitivities for unary linalg optimisations.
_ϵ, lb, ub = 3e-2, -3.0, 3.0
unary_linalg_optimisations = [
    (:-,          ∇Array,  ∇Array,  :(-Ȳ),                               (lb, ub)),
    (:tr,         ∇Array,  ∇Scalar, :(Diagonal(fill!(similar(X), Ȳ))),   (lb, ub)),
    (:inv,        ∇Array,  ∇Array,  :(-transpose(Y) * Ȳ * transpose(Y)), (lb, ub)),
    (:det,        ∇Array,  ∇Scalar, :(Y * Ȳ * transpose(inv(X))),        (_ϵ, ub)),
    (:logdet,     ∇Array,  ∇Scalar, :(Ȳ * transpose(inv(X))),            (_ϵ, ub)),
    (:transpose,  ∇Array,  ∇Array,  :(transpose(Ȳ)),                     (lb, ub)),
    (:adjoint,    ∇Scalar, ∇Scalar, :(adjoint(Ȳ)),                       (_ϵ, ub)),
    (:adjoint,    ∇Array,  ∇Array,  :(adjoint(Ȳ)),                       (lb, ub)),
    (:norm,       ∇Array,  ∇Scalar, :(Ȳ ./ Y .* abs2.(X) ./ X),          (lb, ub)),
    (:norm,       ∇Scalar, ∇Scalar, :(Ȳ * sign(X)),                      (lb, ub))
]
for (f, T_In, T_Out, X̄, bounds) in unary_linalg_optimisations
    if f === :-
        @eval import Base: -
    else
        @eval import LinearAlgebra: $f
    end
    @eval begin
        @explicit_intercepts $f Tuple{$T_In}
        ∇(::typeof($f), ::Type{Arg{1}}, p, Y::$T_Out, Ȳ::$T_Out, X::$T_In) = $X̄
    end
end

# Implementation of sensitivities for binary linalg optimisations.
const A = ∇Array
const S = ∇Scalar
const AS = Union{∇Scalar, ∇Array}
δ = 1e-5
binary_linalg_optimisations = [
    (:*, A, A, AS,
        :(Ȳ * B'),
        :(A' * Ȳ)),
    (:/, A, A, AS,
        :(Ȳ / transpose(B)),
        :(-transpose(Y) * (Ȳ / transpose(B)))),
    (:\, A, A, AS,
        :(-(transpose(A) \ Ȳ) * transpose(Y)),
        :(transpose(A) \ Ȳ)),
    (:norm, A, S, S,
        :(Ȳ .* Y^(1 - B) .* abs.(A).^B ./ A),
        :(Ȳ * (Y^(1 - B) * sum(abs.(A).^B .* log.(abs.(A))) - Y * log(Y)) / B)),
    (:norm, S, S, S,
        :(Ȳ * sign(A)),
        :(0)),
]
import Base: *, /, \
import LinearAlgebra: norm
for (f, T_A, T_B, T_Y, Ā, B̄) in binary_linalg_optimisations
    @eval begin
        @explicit_intercepts $f Tuple{$T_A, $T_B}
        ∇(::typeof($f), ::Type{Arg{1}}, p, Y::$T_Y, Ȳ::$T_Y, A::$T_A, B::$T_B) = $Ā
        ∇(::typeof($f), ::Type{Arg{2}}, p, Y::$T_Y, Ȳ::$T_Y, A::$T_A, B::$T_B) = $B̄
    end
end

# Sensitivities for the Kronecker product:
import LinearAlgebra: kron
@explicit_intercepts kron Tuple{A, A}

# The allocating versions simply allocate and then call the in-place versions.
∇(::typeof(kron), ::Type{Arg{1}}, p, Y::A, Ȳ::A, A::A, B::A) =
    ∇(zeroslike(A), kron, Arg{1}, p, Y, Ȳ, A, B)
∇(::typeof(kron), ::Type{Arg{2}}, p, Y::A, Ȳ::A, A::A, B::A) =
    ∇(zeroslike(B), kron, Arg{2}, p, Y, Ȳ, A, B)

function ∇(Ā::A, ::typeof(kron), ::Type{Arg{1}}, p, Y::A, Ȳ::A, A::A, B::A)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        aij, āij = A[i, j], Ā[i, j]
        for k = reverse(1:K)
            āij += Ȳ[m] * B[k, l]
            m -= 1
        end
        Ā[i, j] = āij
    end
    return Ā
end
function ∇(B̄::A, ::typeof(kron), ::Type{Arg{2}}, p, Y::A, Ȳ::A, A::A, B::A)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        aij = A[i, j]
        for k = reverse(1:K)
            B̄[k, l] += Ȳ[m] * aij
            m -= 1
        end
    end
    return B̄
end

@explicit_intercepts Base.:+ Tuple{A, UniformScaling}
∇(::typeof(+), ::Type{Arg{1}}, p, Y::∇Array, Ȳ::∇Array, A::∇Array, B::UniformScaling) = Ȳ

@explicit_intercepts Base.:+ Tuple{UniformScaling, A}
∇(::typeof(+), ::Type{Arg{2}}, p, Y::∇Array, Ȳ::∇Array, A::UniformScaling, B::∇Array) = Ȳ

# Short-form `dot`.
@explicit_intercepts LinearAlgebra.dot Tuple{∇Array, ∇Array}
∇(::typeof(LinearAlgebra.dot), ::Type{Arg{1}}, p, z, z̄, x::A, y::A) = z̄ .* y
∇(::typeof(LinearAlgebra.dot), ::Type{Arg{2}}, p, z, z̄, x::A, y::A) = z̄ .* x
∇(x̄, ::typeof(LinearAlgebra.dot), ::Type{Arg{1}}, p, z, z̄, x::A, y::A) = (x̄ .= x̄ .+ z̄ .* y)
∇(ȳ, ::typeof(LinearAlgebra.dot), ::Type{Arg{2}}, p, z, z̄, x::A, y::A) = (ȳ .= ȳ .+ z̄ .* x)

# `copy` materializes `Adjoint` and `Transpose` wrappers but can be called on anything
import Base: copy
@explicit_intercepts copy Tuple{Any}
∇(::typeof(copy), ::Type{Arg{1}}, p, Y, Ȳ, A) = copy(Ȳ)

# Matrix exponential
# Ported from Theano, see https://github.com/Theano/Theano/blob/3b8a5b342b30c7ffd2f89f0...
# e9efef601b7492411/theano/tensor/slinalg.py#L518-L553
# Implementation there is based on Kalbfleisch and Lawless, 1985, The Analysis of Panel
# Data Under a Markov Assumption.
import Base: exp
@explicit_intercepts exp Tuple{AbstractMatrix{<:∇Scalar}}
function ∇(::typeof(exp), ::Type{Arg{1}}, p, Y, Ȳ, X::AbstractMatrix)
    # TODO: Make this work for asymmetric matrices
    issymmetric(X) || throw(ArgumentError("input is not symmetric; eigenvalues are complex"))
    n = LinearAlgebra.checksquare(X)
    λ, U = eigen(X)
    eλ = exp.(λ)
    Z = @inbounds begin
        eltype(eλ)[i == j ? eλ[i] : (eλ[i] - eλ[j]) / (λ[i] - λ[j]) for i = 1:n, j = 1:n]
    end
    Uᵀ = transpose(U)
    F = factorize(Uᵀ)
    return real(F \ (Uᵀ * Ȳ / F .* Z) * Uᵀ)
end


text:
function affordable_real(
    irreducible_characters,
    multiplicities=fill(1, length(irreducible_characters)),
)
    irr_real = similar(irreducible_characters, 0)
    mls_real = similar(multiplicities, 0)
    for (i, χ) in pairs(irreducible_characters)
        ι = Characters.frobenius_schur(χ)
        if abs(ι) == 1 # real or quaternionic
            @debug "real/quaternionic:" χ
            push!(irr_real, χ)
            push!(mls_real, multiplicities[i])
        else # complex one...
            cχ = conj(χ)
            k = findfirst(==(cχ), irreducible_characters)
            @assert k !== nothing
            @debug "complex" χ conj(χ)=irreducible_characters[k]
            if k > i # ... we haven't already observed a conjugate of
                @assert multiplicities[i] == multiplicities[k]
                push!(irr_real, χ + cχ)
                push!(mls_real, multiplicities[i])
            end
        end
    end

    return irr_real, mls_real
end

function symmetry_adapted_basis(
    G::PermutationGroups.AbstractPermutationGroup,
    S::Type = Rational{Int};
    semisimple::Bool = false,
)
    tbl = CharacterTable(S, G)
    return symmetry_adapted_basis(eltype(tbl), tbl, semisimple = semisimple)
end

"""
    symmetry_adapted_basis([T::Type,] G::AbstractPermutationGroup[, S=Rational{Int};
        semisimple=false])
Compute a decomposition of `V=𝕂ⁿ` (with `n = degree(G)`) into
subspaces invariant under the natural permutation action of `G`.

Consider the decomposition of `V` into irreducible (simple) subspaces
    `V ≅ m₁V₁ ⊕ ⋯ ⊕ mᵣVᵣ`
and write `Wᵢ = mᵢVᵢ` (the `mᵢ`-fold direct sum of `Vᵢ`). The decomposition is
returned as a vector of `DirectSummand{T}`s (blocks) corresponding
to the distinct irreducible characters of `G` (i.e. types of action, here from `1`
to `r`). Each block contains a basis for a `G`-invariant subspace in `V` (`Vᵢ` or `Wᵢ`). The blocks are guaranteed to be orthogonal to **each other**, however vectors within a single block may *not* be orthogonal.

If `T<:LinearAlgebra.BlasFloat` BLAS routines will be used to orthogonalize
vectors within each block.

!!! note:
Each returned block is invariant (as a subspace) under the action of `G`, which
means that the action may still e.g. permute (row) vectors, but only *within* each block.

Arguments:
* `T` controls the type of coefficients of the returned basis. Unless specified,
the coefficients will be computed exactly in the field of cyclotomic numbers.
If you know that the group has rational characters only (which happens e.g. for
the full symmetric groups) You may specify `Rational{Int}` here. For a group
with complex characters specifying `T<:Real` will result in the computation of the realified basis.
* `S` controls the types of `Cyclotomic`s used in the computation of
character table. Exact type are preferred. For larger groups `G` (or if
overflow occurs during the computation of characters) specifying
`Rational{BigInt}` might be necessary.
* `semisimple::Bool` controls the nature of the the returned basis.
  - `semisimple=true`: the returned basis consists of orthogonal blocks
  (`DirectSummand`s) which define an isomorphism
  `V ≅ Wₖ ⊕ ⋯ ⊕ Wₖ`.
  associated to isotypical components `Wₖ`, which are (in general) semi-simple.
  I.e. each direct summand `ds` may further decopose into `G`-invariant
  subspaces `Wₖ ≅ mₖVₖ`, all of the same _type_. Multiplicity `mₖ` can be
  obtained by calling `multiplicity(@ref)`.
  - `semisimple=false`: (the default) In addition to finding blocks `Wₖ`, an
  effort to find _minimal projection system_ is made, i.e all, some (or none!)
  of the returned blocks corresponds to a **projection** `V → Wₖ ≅ mₖVₖ → Vₖ` for a single irreducible subspace `Vₖ`. This means that some blocks can not be further decomposed into nontrivial `G`-invariant subspaces.
"""
symmetry_adapted_basis(
    T::Type,
    G::PermutationGroups.AbstractPermutationGroup,
    S::Type = Rational{Int};
    semisimple::Bool = false,
) = symmetry_adapted_basis(T, CharacterTable(S, G), semisimple = semisimple)

function symmetry_adapted_basis(
    T::Type,
    tbl::CharacterTable;
    semisimple::Bool = false,
)
    irr, multips = _constituents_decomposition(
        action_character(T, conjugacy_classes(tbl), tbl),
        tbl,
    )

    if T <: Real
        irr, multips = affordable_real(irr, multips)
    end

    if semisimple || all(isone ∘ degree, irr)
        return _symmetry_adapted_basis(T, irr, multips)
    else
        RG = _group_algebra(parent(tbl))
        return _symmetry_adapted_basis(T, irr, multips, RG)
    end
end

"""
    symmetry_adapted_basis([T::Type,] G::Group, action, basis[, S=Rational{Int}];
        semisimple=false)
Compute a decomposition of `V=⟨basis⟩` into subspaces invariant under the
given action of `G`.

It is assumed that `G` acts on a subset of basis and the action needs to be
extended to the whole `basis`. If `G` is a permutation group already acting on
the whole `basis`, a call to `symmetry_adapted_basis(G)` is preferred.

* For inducing the action `basis` needs to be indexable and iterable
(e.g. in the form of an `AbstractVector`).

"""
function symmetry_adapted_basis(
    G::Group,
    action::Action,
    basis,
    S::Type = Rational{Int};
    semisimple=false,
)
    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true)
    return symmetry_adapted_basis(eltype(tbl), tbl, ehom, semisimple=semisimple)
end

function symmetry_adapted_basis(
    T::Type,
    G::Group,
    action::Action,
    basis,
    S::Type = Rational{Int};
    semisimple=false,
)
    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true)
    return symmetry_adapted_basis(T, tbl, ehom, semisimple=semisimple)
end

function symmetry_adapted_basis(
    T::Type,
    tbl::CharacterTable,
    ehom::InducedActionHomomorphism;
    semisimple=false,
)
    ψ = action_character(ehom, tbl)

    irr, multips = _constituents_decomposition(ψ, tbl)
    if T <: Real
        irr, multips = affordable_real(irr, multips)
        @debug "Decomposition into real character spaces:
        degrees:        $(join([lpad(d, 6) for d in degree.(irr)], ""))
        multiplicities: $(join([lpad(m, 6) for m in multips], ""))"

    end

    if semisimple || all(isone ∘ degree, irr)
        return _symmetry_adapted_basis(T, irr, multips, ehom)
    else
        RG = _group_algebra(parent(tbl))
        return _symmetry_adapted_basis(T, irr, multips, RG, ehom)
    end
end

function _constituents_decomposition(ψ::Character, tbl::CharacterTable)
    irr = irreducible_characters(tbl)
    degrees = degree.(irr)
    multiplicities = constituents(ψ)

    @debug "Decomposition into character spaces:
    degrees:        $(join([lpad(d, 6) for d in degrees], ""))
    multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))"

    @assert dot(multiplicities, degrees) == degree(ψ)
    "Something went wrong: characters do not constitute a complete basis for action:
    $(dot(multiplicities, degrees)) ≠ $(degree(ψ))"

    present_irreps = [i for (i, m) in enumerate(multiplicities) if m ≠ 0]
    return irr[present_irreps], multiplicities[present_irreps]
end

function _group_algebra(G::Group)
    @assert isfinite(G)
    b = StarAlgebras.Basis{UInt16}(vec(collect(G)))
    RG = if order(Int, G) <= (typemax(UInt16)>>2)
        StarAlgebra(G, b, (length(b), length(b)), precompute=true)
        # cache is about ~ 1Gb
    else
        StarAlgebra(G, b)
    end
    return RG
end

function _symmetry_adapted_basis(
    T::Type,
    irr::AbstractVector{<:Character},
    multiplicities::AbstractVector{<:Integer},
    hom=nothing
)
    res = map(zip(irr, multiplicities)) do (µ, m)
        Threads.@spawn begin
            µT = eltype(µ) == T ? µ : Character{T}(µ)
            image = hom === nothing ? image_basis(µT) : image_basis(hom, µT)
            simple = size(image, 1) == m
            deg = degree(µ)
            @assert size(image, 1) == (simple ? m : m*deg) "incompatible projection dimension: $(size(image, 1)) ≠ $(simple ? m : m*deg)"
            if deg == 1
                @assert simple "Central projection associated to character is not simple unless its degree == 1"
            end
            DirectSummand(image, m, deg, simple)
        end
    end
    return fetch.(res)
end

function _symmetry_adapted_basis(
    T::Type,
    irr::AbstractVector{<:Character},
    multiplicities::AbstractVector{<:Integer},
    RG::StarAlgebra{<:Group},
    hom=nothing,
)
    mps, simples = minimal_projection_system(irr, RG)
    degrees = degree.(irr)
    res = map(zip(mps, multiplicities, degrees, simples)) do (µ, m, deg, simple)
        Threads.@spawn begin
            µT = eltype(µ) == T ? µ : AlgebraElement{T}(µ)
            image = hom === nothing ? image_basis(µT) : image_basis(hom, µT)
            @assert size(image, 1) == (simple ? m : m*deg) "incompatible projection dimension: $(size(image, 1)) ≠ $(simple ? m : m*deg)"
            DirectSummand(image, m, deg, simple)
        end
    end
    direct_summands = fetch.(res)

    for (χ, ds) in zip(irr, direct_summands)
        if issimple(ds) && (d = size(ds, 1)) != (e = multiplicity(ds)*sum(constituents(χ).>0))
            throw("The dimension of the projection doesn't match with simple summand multiplicity: $d ≠ $e")
        end
    end

    return direct_summands
end


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

function imfilter_gpu(img::AbstractArray{T,N},
                      kern::AbstractArray{K,N}) where {T<:Real,K<:Real,N}
  # GPU metadata
  ctx = GPU.ctx; queue = GPU.queue
  mult_kernel = GPU.mult_kernel

  # operations with complex type
  img  = map(Complex64, img)
  kern = map(Complex64, kern)

  # kernel may require padding
  prepad  = ntuple(d->(size(kern,d)-1) ÷ 2, N)
  postpad = ntuple(d->(size(kern,d)  ) ÷ 2, N)

  # OpenCL FFT expects powers of 2, 3, 5, 7, 11 or 13
  clpad = clfftpad(img)
  A = padarray(img, Pad(:symmetric, zeros(Int, ndims(img)), clpad))
  A = parent(A)

  krn = zeros(Complex64, size(A))
  indexesK = ntuple(d->[size(A,d)-prepad[d]+1:size(A,d);1:size(kern,d)-prepad[d]], N)
  krn[indexesK...] = reflect(kern)

  # plan FFT
  p = clfft.Plan(Complex64, ctx, size(A))
  clfft.set_layout!(p, :interleaved, :interleaved)
  clfft.set_result!(p, :inplace)
  clfft.bake!(p, queue)

  # populate GPU memory
  bufA   = cl.Buffer(Complex64, ctx, :copy, hostbuf=A)
  bufkrn = cl.Buffer(Complex64, ctx, :copy, hostbuf=krn)
  bufRES = cl.Buffer(Complex64, ctx, length(A))

  # compute ifft(fft(A).*fft(kern))
  clfft.enqueue_transform(p, :forward, [queue], bufA, nothing)
  clfft.enqueue_transform(p, :forward, [queue], bufkrn, nothing)
  queue(mult_kernel, length(A), nothing, bufA, bufkrn, bufRES)
  clfft.enqueue_transform(p, :backward, [queue], bufRES, nothing)

  # get result back
  AF = reshape(cl.read(queue, bufRES), size(A))

  # undo OpenCL FFT paddings
  AF = view(AF, ntuple(d->1:size(AF,d)-clpad[d], N)...)

  out = Array{realtype(eltype(AF))}(ntuple(d->size(img,d) - prepad[d] - postpad[d], N))
  indexesA = ntuple(d->postpad[d]+1:size(img,d)-prepad[d], N)
  copyreal!(out, AF, indexesA)

  out
end

@generated function reflect{T,N}(A::AbstractArray{T,N})
    quote
        B = Array{T}(size(A))
        @nexprs $N d->(n_d = size(A, d)+1)
        @nloops $N i A d->(j_d = n_d - i_d) begin
            @nref($N, B, j) = @nref($N, A, i)
        end
        B
    end
end

for N = 1:5
    @eval begin
        function copyreal!{T<:Real}(dst::Array{T,$N}, src, I::Tuple{Vararg{UnitRange{Int}}})
            @nexprs $N d->(I_d = I[d])
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = real(@nref $N src j)
            end
            dst
        end
        function copyreal!{T<:Complex}(dst::Array{T,$N}, src, I::Tuple{Vararg{UnitRange{Int}}})
            @nexprs $N d->I_d = I[d]
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = @nref $N src j
            end
            dst
        end
    end
end

realtype{R<:Real}(::Type{R}) = R
realtype{R<:Real}(::Type{Complex{R}}) = R


text:
#= Testing data access patterns

The following methods will be evaluated:
 	_variable_
	varsubset
 	datasubset
	getx/getx!
	getindex
	setindex
=#

function t_data_access_patterns()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)
AV = [A01, A02, A03, A04, A05, A06]
AVv = [A01, A02, A03]
AVm = [A04, A05, A06]

# Data selection
# --------------

# Select observations
#println("Checking data selection...")
#print("\t|> Observations: datasubset - getindex - getx/getx!...")
for idx in [1,1:1,1:2,[1,2]]
	# Test for vector data
	for A in AVv
		Test.@test all(+j4pr.datasubset(A,idx) .== +A[idx] .== j4pr.getx(A)[idx] .== j4pr.getx!(A)[idx])
	end

	# Test for matrix data
	for A in AVm
		Test.@test all(+j4pr.datasubset(A,idx) .== +A[:,idx] .== j4pr.getx(A)[:,idx] .== j4pr.getx!(A)[:,idx])
	end
end
#println("PASSED")

# Select variables
#print("\t|> Variables: varsubset - _variable_ - getindex - getx/getx!...")
for idx in [1,1:1,1:2,[1,2]]
	# Test for vector data
	for A in AVv
		if idx isa Int 	# idx has to be 1 otherwise most of the variable access methods fail, as they should;
				# Note that there is no getindex method for A
			Test.@test all(+j4pr.varsubset(A,idx) .== j4pr._variable_(A,idx) .== +A .== j4pr.getx(A)[:,idx] .== j4pr.getx!(A)[:,idx])
		end
	end

	# Test for matrix data
	for A in AVm
		if idx isa Int
			# data cell indexing returns a matrix data if single row selected, result must be transposed as all the other methods use vectors
			Test.@test all(+j4pr.varsubset(A,idx) .== j4pr._variable_(A,idx) .== (+A[idx,:])' .== j4pr.getx(A)[idx,:] .== j4pr.getx!(A)[idx,:])
		else
			Test.@test all(+j4pr.varsubset(A,idx) .== j4pr._variable_(A,idx) .== +A[idx,:] .== j4pr.getx(A)[idx,:] .== j4pr.getx!(A)[idx,:])
		end
	end
end
#println("PASSED")

# Selecting observations + variables is done by applying successively one of the methods for variable selection and observation selection

# Data modification
# -----------------

# Modify observations
#println("Checking data modification...")
#print("\t|> Observations: datasubset - settindex - getx! and NOT getx...")
for idx in [1,1:1,1:2,[1,2]]
	# Test for vector data
	for A in AVv
		# datasubset approach
		Ac = deepcopy(A)
		newvals = idx isa Int ? rand() : rand(length(idx))
		j4pr.datasubset(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr.datasubset(Ac,idx) .== newvals)

		# setindex! approach
		Ac = deepcopy(A)
		newvals = idx isa Int ? rand() : rand(length(idx))
		Ac[idx] = newvals;
		Test.@test all(+Ac[idx] .== newvals)

		# getx! approach
		Ac = deepcopy(A)
		newvals = idx isa Int ? rand() : rand(length(idx))
		j4pr.getx!(Ac)[idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[idx] .== newvals)

		# getx approach (MUST NOT MODIDY)
		Ac = deepcopy(A)
		newvals = idx isa Int ? rand() : rand(length(idx))
		j4pr.getx(Ac)[idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[idx] .!= newvals)
		Test.@test all(j4pr.getx!(Ac)[idx] .== j4pr.getx!(A)[idx])
	end

	# Test for matrix data
	for A in AVm
		# datasubset approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nvars(Ac), length(idx))
		j4pr.datasubset(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr.datasubset(Ac,idx) .== newvals)

		# setindex! approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nvars(Ac), length(idx))
		Ac[:,idx] = newvals;
		Test.@test all(+Ac[:,idx] .== newvals)

		# getx! approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nvars(Ac), length(idx))
		j4pr.getx!(Ac)[:,idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[:,idx] .== newvals)

		# getx approach (MUST NOT MODIDY)
		Ac = deepcopy(A)
		newvals = rand(j4pr.nvars(Ac), length(idx))
		j4pr.getx(Ac)[:,idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[:,idx] .!= newvals)
		Test.@test all(j4pr.getx!(Ac)[:,idx] .== j4pr.getx!(A)[:,idx])
	end
end
#println("PASSED")

# Modify variables
#print("\t|> Variables: varsubset - _variable_ - settindex - getx! and NOT getx...")
for idx in [1,1:1,1:2,[1,2]]
	# Test for vector data
	for A in AVv
		if idx isa Int
		# varsubset approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nobs(A))
		j4pr.varsubset(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr.varsubset(Ac,idx) .== newvals)

		# _variable_ approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nobs(A))
		j4pr._variable_(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr._variable_(Ac,idx) .== newvals)

		# setindex! approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nobs(A))
		Ac[:] = newvals;
		Test.@test all(+Ac .== newvals)

		# getx! approach
		Ac = deepcopy(A)
		newvals = rand(j4pr.nobs(A))
		j4pr.getx!(Ac)[:,idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[:,idx] .== newvals)

		# getx approach (MUST NOT MODIDY)
		Ac = deepcopy(A)
		newvals = rand(j4pr.nobs(A))
		j4pr.getx(Ac)[:,idx] = newvals;
		Test.@test all(j4pr.getx!(Ac)[:,idx] .!= newvals)
		Test.@test all(j4pr.getx!(Ac)[:,idx] .== j4pr.getx!(A)[:,idx])
		end
	end
	# Test for matrix data
	for A in AVm
		# varsubset approach
		Ac = deepcopy(A)
		newvals = rand(size(j4pr.varsubset(Ac,idx)[:]))
		j4pr.varsubset(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr.varsubset(Ac,idx) .== newvals)

		# _variable_ approach
		Ac = deepcopy(A)
		newvals = rand(size(j4pr._variable_(Ac,idx)))
		j4pr._variable_(Ac,idx)[:] = newvals;
		Test.@test all(+j4pr._variable_(Ac,idx) .== newvals)

		# setindex! approach
		Ac = deepcopy(A)
		newvals = rand(size(Ac[idx,:]))
		Ac[idx,:] = newvals;
		Test.@test all(+Ac[idx,:] .== newvals)

		# getx! approach
		Ac = deepcopy(A)
		newvals = rand(size(j4pr.getx!(Ac)[idx,:]))
		j4pr.getx!(Ac)[idx,:] = newvals;
		Test.@test all(j4pr.getx!(Ac)[idx,:] .== newvals)

		# getx approach (MUST NOT MODIDY)
		Ac = deepcopy(A)
		newvals = rand(size(j4pr.getx!(Ac)[idx,:]))
		j4pr.getx(Ac)[idx,:] = newvals;
		Test.@test all(j4pr.getx!(Ac)[idx,:] .!= newvals)
		Test.@test all(j4pr.getx!(Ac)[idx,:] .== j4pr.getx!(A)[idx,:])
	end
end
#println("PASSED")

end


text:
## Optimization setup

# Derivative stuff for learning. Consider OpenMDAO to use MAUD generalizations.
solve_direct(x, p, ∂R∂x, ∂R∂p) = @views reshape(reduce(hcat, ∂R∂p \ -(∂R∂x)[:,:,i][:] for i in eachindex(x)), (size(p)..., length(x)))

solve_adjoint(f, p, ∂R∂p, dfdp) = @views reshape(reduce(hcat, ∂R∂p' \ -(dfdp)'[i,:] for i in eachindex(f)), (length(f), size(p)...)

total_derivative_direct(∂f∂x, ψ, ∂f∂p) = @views ∂f∂x + [ sum(∂f∂p * ψ[n]) for n in eachindex(∂f∂x) ]

total_derivative_adjoint(∂f∂x, φ, ∂R∂x) = @views ∂f∂x + [ sum(permutedims(φ) * reshape(∂R∂x, (length(R[:]), length(∂f∂x)))[:,n]) for n in eachindex(∂f∂x) ]

max_stars_repo_name: GodotMisogi/AdjointCFD

max_issues_repo_name: GodotMisogi/AdjointCFD

max_forks_repo_name: GodotMisogi/AdjointCFD


text:
using Gadfly
using LinearAlgebra
using VoronoiDelaunay
using GeometricalPredicates

function distmesh(fd, fh, h0, bbox, pfix) end

# Co jest ten broadcasting aka .+ w 12

function fh(x::T, y::T) where {T}
    return 0.2
end

function generate(fd, fh, bbox, h0, pfix)
    ttol = 0.1
    geps = 0.001 * h0
    h1 = h0 * (sqrt(3) / 2)
    v1 = bbox[1, 1]:h0:bbox[2, 1]
    v2 = bbox[1, 2]:h1:bbox[2, 2]
    x, y = meshgrid(v1, v2)
    x[2:2:end, :] = x[2:2:end, :] .+ h0 / 2
    p = [x[:] y[:]]
    z = [vcat(a) for a in eachrow(p) if fd(a) < geps]
    f = transpose(reshape(vcat(z...), 2, length(z)))
    r0 = 1 ./ fh(f[:, 1], f[:, 2]) .^ 2
    r0_max = r0 ./ max(r0)
    l = [vcat(b) for b in eachrow(f) if (rand(Float64, size(f)))[1] < r0_max]
    g = transpose(reshape(vcat(z...), 2, length(z)))
    g
    N = size(g, 1)
    N
    println(N)

    pold = Inf
    scale = abs(1 / (bbox[1] - bbox[2]))
    transx = 1 + (0 - bbox[1] * scale)
    transy = 1 + (0 - bbox[1, 2] * scale)

    maximum(sqrt.(sum((g .- pold) .^ 2, dims = 2) / h0)) > ttol
    tess = DelaunayTessellation()
    aa = map(el -> Point(el[1] * scale + transx, el[2] * scale + transy), eachrow(g))
    aa
    println(aa)

    push!(tess, aa)
    zda = Array{Point2D,1}()
    for tri in tess
        push!(zda, centroid(Primitive(geta(tri), getb(tri), getc(tri))))
    end
    zza = filter(cent -> dc([getx(cent) gety(cent)]) > -geps, zda)
    zdda = Array{Float64,1}()
    barvec = Array{Array{Float64,1},1}()
    for edge in delaunayedges(tess)
        push!(zdda, sqrt(length2(Line(geta(edge), getb(edge)))) / scale)
        push!(
            barvec,
            [
                (getx(getb(edge)) - getx(geta(edge))) / scale,
                ((gety(getb(edge)) - gety(geta(edge)))) / scale,
            ],
        )
    end
    hbars = 0.2
    ss = 0.2 * 1.2 * sqrt(sum(zdda .^ 2) / (0.2^2))
    ttt = Array{Float64,1}()
    for edge in delaunayedges(tess)
        push!(ttt, ss - sqrt(length2(Line(geta(edge), getb(edge)))) / scale)
    end
    rttt = reshape(ttt, :, 1)
    mzdda = reshape(zdda, :, 1)
    init_Fvec = rttt ./ mzdda
    #    Fvec = Fvec.*barvec
    l
    #    println(sort(with_f[:]))

    #mzzda
    #    width = max_coord - min_coord
    #    a = Point2D[Point(min_coord + rand() * width, min_coord + rand() * width) for i in 1:100]
    #    push!(tess, a)
    #    x, y = getplotxy(delaunayedges(tess))
    #    Gadfly.plot(x=x, y=y, Geom.path)
end

function evaluate(fd, p, geps)
    for x in eachrow(p)
        if (fd(x) < geps)

        end
    end
end

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)

    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)

    return gx, gy
end

function feval(f, x...) end

function generateNodes()
    x = [1, 2, 3, 4]
    y = [1, 1, 2, 3]
    plot(x, y)
end

function dc(p)
    return sqrt(sum(p .^ 2)) - 1
end

function dcircle(p, xc, yc, r)
    return sqrt.((p[:, 1] .- xc) .^ 2 .+ (p[:, 2] .- yc) .^ 2) .- r
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
using Oceananigans.Grids: AbstractGrid

""" Store source terms for `u`, `v`, and `w`. """
@kernel function store_field_tendencies!(G⁻, grid::AbstractGrid{FT}, G⁰) where FT
    i, j, k = @index(Global, NTuple)
    G⁻[i, j, k] = G⁰[i, j, k]
end

""" Store previous source terms before updating them. """
function store_tendencies!(model)

    barrier = Event(device(model.architecture))

    workgroup, worksize = work_layout(model.grid, :xyz)

    model_fields = fields(model)

    store_field_tendencies_kernel! = store_field_tendencies!(device(model.architecture), workgroup, worksize)

    events = []

    for (i, field) in enumerate(model_fields)

        field_event = store_field_tendencies_kernel!(model.timestepper.G⁻[i],
                                                     model.grid,
                                                     model.timestepper.Gⁿ[i],
                                                     dependencies=barrier)

        push!(events, field_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end


text:
include( "../src/AgnosticBayesEnsemble.jl" )
include( "../src/dirichletPosteriorEstimation.jl" )
using DataFrames
using Random
using Statistics
using Test
using Optim
using MultivariateStats

println( "running dirichletian algorithm unit tests" );
predMatTraining, predMatEval, tTraining, tEval, errMatTraining = makeupPredictions();
posterior = dirichletPosteriorEstimation( errMatTraining, 1000, 2.34 );

@test all( posterior .>= 0.0 )
@test all( posterior .<= 1.0 )
@test sum( posterior ) ≈ 1.0

posterior = dirichletPosteriorEstimationV2( errMatTraining, 100, 2.34, 10 );

@test all( posterior .>= 0.0 )
@test all( posterior .<= 1.0 )
@test sum( posterior ) ≈ 1.0

posterior = zeros( Float64, size( errMatTraining, 2 ) );
dirichletPosteriorEstimation!( errMatTraining, 1000, 2.34, posterior )

@test all( posterior .>= 0.0 )
@test all( posterior .<= 1.0 )
@test sum( posterior ) ≈ 1.0

αSequence, performance = metaParamSearchValidationDirichlet( predMatTraining, tTraining, 1000, 0.05, 10.0, 20, 0.25, lossFunctions.MSE );
best_index            = argmaxUniProb.argmaxProb( performance );

posterior_tuned = dirichletPosteriorEstimation( errMatTraining, 5000, αSequence[best_index] );
@test all( posterior_tuned .>= 0.0 )
@test all( posterior_tuned .<= 1.0 )
@test sum( posterior_tuned ) ≈ 1.0

posteriorUni = zeros( Float64, size( errMatTraining, 2 ) )
fill!( posteriorUni, 1 / size( errMatTraining, 2 ) );

yUni = round.( AgnosticBayesEnsemble.predictEnsemble( predMatEval, posteriorUni ) );
yDir = round.( AgnosticBayesEnsemble.predictEnsemble( predMatEval, posterior_tuned ) );


text:
module MockFieldsTests

using Gridap.Fields
using Gridap.Fields: MockField, MockBasis, OtherMockBasis
using Gridap.TensorValues

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
f = MockField{d}(v)
fx = fill(v,np)
test_field(f,x,fx)

∇fx = fill(VectorValue(v,0.0),np)
∇f = gradient(f)

test_field(∇f,x,∇fx)
test_field(f,x,fx,grad=∇fx)

ndof = 8
b = MockBasis{d}(v,ndof)
bx = fill(v,np,ndof)
∇bx = fill(VectorValue(v,0.0),np,ndof)
test_field(b,x,bx,grad=∇bx)

b = OtherMockBasis{d}(ndof)
bx = fill(2*p,np,ndof)
∇bx = fill(TensorValue(2.0,0.0,0.0,2.0),np,ndof)
test_field(b,x,bx,grad=∇bx)

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Conversions related to the orbit elements.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Schwarz, R (2014). Memorandum No. 2: Cartesian State Vectors to
#       Keplerian Orbit Elements. Available at www.rene-schwarz.com
#
#       https://downloads.rene-schwarz.com/dc/category/18
#       (Accessed on 2017-08-09).
#
#   [2] Vallado, D. A., McClain, W. D (2013). Fundamentals of Astrodynamics
#       and Applications. Microcosm Press.
#
#   [3] Kuga, H. K., Carrara, V., Rao, K. R (2005). Introdução à Mecânica
#       Orbital. 2ª ed. Instituto Nacional de Pesquisas Espaciais.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export change_oe_frame
export kepler_to_rv
export rv_to_kepler
export kepler_to_sv, sv_to_kepler

"""
    change_oe_frame(a::Number, e::Number, i::Number, Ω::Number, ω::Number, f::Number, conv_args...)
    change_oe_frame(oe::Orbit, conv_args...)

Change the reference frame of orbit elements. The orbit elements can be
specified by `a`, `e`, `i`, `Ω`, `ω`, and `f`, or the structure `oe` (see
`Orbit`).

The conversion arguments `conv_args` are **the same** arguments that one should
pass to the function `rECItoECI` to convert between the desired frames. For more
information, see the documentation of the function `rECItoECI`.

# Args

* `a`: Semi-major axis [m].
* `e`: Excentricity.
* `i`: Inclination [rad].
* `Ω`: Right-ascension of the ascending node [rad].
* `ω`: Argument of perigee [rad].
* `f`: True anomaly [rad].
* `conv_args...`: Conversion arguments, which are the same arguments that one
                  would pass to the function `rECItoECI` to convert between the
                  desired frames.

* `oe`: An instance of the structure `Orbit` with the orbit elements that will
        be converted [SI units].

# Returns

An instance of the structure `Orbit` with the Keplerian elements [SI units]
converted to the new frame.

# Examples

```julia-repl
julia> eop = get_iers_eop(:IAU1980);

julia> teme_epoch = DatetoJD(2016,6,1,11,0,0);

julia> tod_epoch  = DatetoJD(2016,1,1,0,0,0);

julia> oe_teme    = Orbit(0,
                          7130.982e3,
                          0.001111,
                          98.405*pi/180,
                          227.336*pi/180,
                          90*pi/180,
                          320*pi/180)

                 Orbit
  =====================================
                  t =        0.0
    Semi-major axis =     7130.9820 km
       Eccentricity =        0.001111
        Inclination =       98.4050 ˚
               RAAN =      227.3360 ˚
    Arg. of Perigee =       90.0000 ˚
       True Anomaly =      320.0000 ˚

julia> oe_j2000 = change_oe_frame(oe_teme, TEME(), J2000(), teme_epoch, eop)

                 Orbit
  ======================================
                  t =        0.0
    Semi-major axis =     7130.9820 km
       Eccentricity =        0.001111
        Inclination =       98.3365 ˚
               RAAN =      227.1345 ˚
    Arg. of Perigee =       90.0604 ˚
       True Anomaly =      320.0000 ˚

julia> oe_tod   = change_oe_frame(oe_teme, TEME(), teme_epoch, TOD(), tod_epoch, eop)

                 Orbit
  ======================================
                  t =        0.0
    Semi-major axis =     7130.9820 km
       Eccentricity =        0.001111
        Inclination =       98.4037 ˚
               RAAN =      227.3306 ˚
    Arg. of Perigee =       90.0014 ˚
       True Anomaly =      320.0000 ˚
```
"""
function change_oe_frame(a::Number,
                         e::Number,
                         i::Number,
                         Ω::Number,
                         ω::Number,
                         f::Number,
                         conv_args...)

    # The approach is to transform the orbit elements to Cartesian
    # representation, convert the frame, and then convert back to orbit
    # elements.
    #
    # NOTE: In my benchmarks, the operation with DCMs are faster than with
    # quaternions after the DCM representation was changed to SMatrix.

    r_o, v_o    = kepler_to_rv(a, e, i, Ω, ω, f)
    D_ECIf_ECIo = rECItoECI(DCM, conv_args...)
    r_f         = D_ECIf_ECIo*r_o
    v_f         = D_ECIf_ECIo*v_o

    rv_to_kepler(r_f, v_f)
end

function change_oe_frame(oe::Orbit, conv_args...)
    change_oe_frame(oe.a, oe.e, oe.i, oe.Ω, oe.ω, oe.f, conv_args...)
end

"""
    kepler_to_rv(a::Number, e::Number, i::Number, Ω::Number, ω::Number, f::Number)
    kepler_to_rv(o::Orbit)

Convert the Keplerian elements (`a`, `e`, `i`, `Ω`, `ω`, and `f`) to a Cartesian
representation (position vector `r` and velocity vector `v`). The Keplerian
elements can also be passed inside an instance of the `Orbit` structure.

# Args

* `a`: Semi-major axis [m].
* `e`: Excentricity.
* `i`: Inclination [rad].
* `Ω`: Right ascension of the ascending node [rad].
* `ω`: Argument of perigee [rad].
* `f`: True anomaly [rad].

# Returns

* The position vector represented in the inertial reference frame [m].
* The velocity vector represented in the inertial reference frame [m].

# References

This algorithm was adapted from [1] and [3, p. 37-38].

"""
function kepler_to_rv(a::Number,
                      e::Number,
                      i::Number,
                      Ω::Number,
                      ω::Number,
                      f::Number)
    # Check eccentricity.
    if !(0 <= e < 1)
        throw(ArgumentError("Eccentricity must be in the interval [0,1)."))
    end

    # Auxiliary variables.
    sin_f, cos_f = sincos(f)

    # Compute the geocentric distance.
    r = a*(1-e^2)/(1+e*cos_f)

    # Compute the position vector in the orbit plane, defined as:
    #   - The X axis points towards the perigee;
    #   - The Z axis is perpendicular to the orbital plane (right-hand);
    #   - The Y axis completes a right-hand coordinate system.
    r_o = SVector{3}(r*cos_f, r*sin_f, 0)

    # Compute the velocity vector in the orbit plane without perturbations.
    n = angvel(a, e, i, :J0)
    v_o = n*a/sqrt(1-e^2)*SVector{3}(-sin_f, e+cos_f, 0)

    # Compute the matrix that rotates the orbit reference frame into the
    # inertial reference frame.
    Dio = angle_to_dcm(-ω, -i, -Ω, :ZXZ)

    # Compute the position and velocity represented in the inertial frame.
    r_i = Dio*r_o
    v_i = Dio*v_o

    (r_i, v_i)
end

kepler_to_rv(o::Orbit) = kepler_to_rv(o.a, o.e, o.i, o.Ω, o.ω, o.f)

"""
    rv_to_kepler(r_i::AbstractVector, v_i::AbstractVector, t::Number = 0)

Convert a Cartesian representation (position vector `r` [m] and velocity vector
`v` [m/s²]) to the Keplerian elements. Optionally, the user can specify the
epoch of the returned elements using the parameter `t`. It it is omitted, then
it default to 0.

# Returns

An instance of the structure `Orbit` with the Keplerian elements [SI units].

# Remarks

The special cases are treated as follows:

* **Circular and equatorial**: the right ascension of the ascending node and the
  argument of perigee are set to 0. Hence, the true anomaly is equal to the true
  longitude.
* **Elliptical and equatorial**: the right ascension of the ascending node is
  set to 0. Hence, the argument of perigee is equal to the longitude of
  periapsis.
* **Circular and inclined**: the argument of perigee is set to 0. Hence, the
  true anomaly is equal to the argument of latitude.

# References

The algorithm was adapted from [1].

"""
function rv_to_kepler(r_i::AbstractVector, v_i::AbstractVector, t::Number = 0)
    # Check inputs.
    length(r_i) != 3 && error("The vector r_i must have 3 dimensions.")
    length(v_i) != 3 && error("The vector v_i must have 3 dimensions.")

    @inbounds begin

        # Position and velocity vector norms and auxiliary dot products.
        r2 = dot(r_i,r_i)
        v2 = dot(v_i,v_i)

        r  = sqrt(r2)
        v  = sqrt(v2)

        rv = dot(r_i,v_i)

        # The type of `r` will be the type of the orbit elements.
        T   = typeof(r)
        Tm0 = T(m0)

        # Angular momentum vector.
        h_i = cross( r_i, v_i )
        h   = norm(h_i)

        # Vector that points to the right ascension of the ascending node (RAAN).
        n_i = SVector{3}(0,0,1) × h_i
        n   = norm(n_i)

        # Eccentricity vector.
        e_i = ( (v2 - Tm0/r)*r_i - rv*v_i )/Tm0

        # Orbit energy.
        ξ = v2/2 - Tm0/r

        # Eccentricity
        # ============

        ecc = norm(e_i)

        # Semi-major axis
        # ===============

        if abs(ecc) <= 1.0-1e-6
            a = -Tm0/(2ξ)
        else
            error("Could not convert the provided Cartesian values to Kepler elements.\n" *
                  "The computed eccentricity was not between 0 and 1");
        end

        # Inclination
        # ===========

        cos_i = h_i[3]/h
        cos_i = abs(cos_i) > 1 ? sign(cos_i) : cos_i
        i     = acos(cos_i)

        # Check the type of the orbit to account for special cases
        # ======================================================================

        # Equatorial
        # ----------------------------------------------------------------------

        if abs(n) <= 1e-6

            # Right Ascension of the Ascending Node.
            # ======================================

            Ω = T(0)

            # Equatorial and elliptical
            # ------------------------------------------------------------------

            if abs(ecc) > 1e-6

                # Argument of Perigee
                # ===================

                cos_ω = e_i[1]/ecc
                cos_ω = abs(cos_ω) > 1 ? sign(cos_ω) : cos_ω
                ω     = acos(cos_ω)

                (e_i[3] < 0) && (ω = T(2π) - ω)

                # True anomaly
                # ============

                cos_v = dot(e_i,r_i)/(ecc*r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (rv < 0) && (v = T(2π) - v)

            # Equatorial and circular
            # ------------------------------------------------------------------

            else
                # Argument of Perigee
                # ===================

                ω = T(0)

                # True anomaly
                # ============

                cos_v = r_i[1]/r
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (r_i[2] < 0) && (v = T(2π) - v)
            end

        # Inclined
        # ----------------------------------------------------------------------
        else

            # Right Ascension of the Ascending Node.
            # ======================================

            cos_Ω = n_i[1]/n
            cos_Ω = abs(cos_Ω) > 1 ? sign(cos_Ω) : cos_Ω
            Ω     = acos(cos_Ω)

            (n_i[2] < 0) && (Ω = T(2π) - Ω)

            # Circular and inclined
            # ------------------------------------------------------------------

            if abs(ecc) < 1e-6

                # Argument of Perigee
                # ===================

                ω = T(0)

                # True anomaly
                # ============

                cos_v = dot(n_i,r_i)/(n*r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (r_i[3] < 0) && (v = T(2π) - v)
            else

                # Argument of Perigee
                # ===================

                cos_ω = dot(n_i,e_i)/(n*ecc)
                cos_ω = abs(cos_ω) > 1 ? sign(cos_ω) : cos_ω
                ω     = acos(cos_ω)

                (e_i[3] < 0) && (ω = T(2π) - ω)

                # True anomaly
                # ============

                cos_v = dot(e_i,r_i)/(ecc*r)
                cos_v = abs(cos_v) > 1 ? sign(cos_v) : cos_v
                v     = acos(cos_v)

                (rv < 0) && (v = T(2π) - v)
            end
        end
    end

    # Return the Keplerian elements.
    # ==============================

    Orbit(t,a,ecc,i,Ω,ω,v)
end

"""
    rv_to_kepler(x::Number, y::Number, z::Number, vx::Number, vy::Number, vz::Number)

Convert a Cartesian representation (position vector `[x;y;z]` [m] and velocity
vector `[vx;vy;vz]` [m/s²]) to the Keplerian elements.

# Returns

An instance of the structure `Orbit` with the Keplerian elements [SI units].

"""
function rv_to_kepler(x::Number,  y::Number,  z::Number,
                      vx::Number, vy::Number, vz::Number)
    # Create the position and velocity vectors.
    r_i = SVector{3}( x, y, z)
    v_i = SVector{3}(vx,vy,vz)

    # Compute the Keplerian orbit elements.
    rv_to_kepler(r_i,v_i)
end

"""
    kepler_to_sv(orb::Orbit)

Convert the Keplerian elements in the structure `orb` to a state vector.

"""
function kepler_to_sv(orb::Orbit)
    r_i, v_i = kepler_to_rv(orb)
    satsv(orb.t, r_i, v_i)
end

"""
    sv_to_kepler(sv::SatelliteStateVector)

Convert the state vector `sv` to Keplerian elements represented by an instance
of the structure `Orbit`.

"""
sv_to_kepler(sv::SatelliteStateVector) = rv_to_kepler(sv.r, sv.v, sv.t)


text:
module Lattices

using ..LatticeVars
import ..LatticeVars.init!, ..LatticeVars.get_num_sites_visited
using ..Points
using Distributions
using JDistributions: BernoulliBool

export AbstractLattice, Lattice, AbstractLatticeTypes, Traps, TrapsLattice
export SitePercolationLattice, SitePercolation

abstract type AbstractLattice end

abstract type AbstractLatticeTypes end
struct Traps <: AbstractLatticeTypes
end

"""
    Lattice{N, T} <: AbstractLattice

Type representing variables on the sites of an `N`-dimensional lattice.
`T` is a `Tuple` type of variables. These variables are meant to be a mixture of
types `Distribution` (annealed variables) and `LatticeVars` (quenched variables).
"""
struct Lattice{N, T, Opts} <: AbstractLattice
    vars::T
    Lattice{N}(vars::Tuple; lattice_type = Traps) where N = new{N, typeof(vars), lattice_type}(vars)
    function Lattice{N}(vars...; lattice_type = Traps) where N
        vartup = (vars...,)
        return new{N, typeof(vartup), lattice_type}(vartup)
    end
end

Lattice(args...; kwargs...) = Lattice{1}(args...; kwargs...)

const TrapsLattice = Lattice{N, T, Traps} where {N, T}

struct SitePercolation <: AbstractLatticeTypes
end

const SitePercolationLattice = Lattice{N, T, SitePercolation} where {N, T}

function SitePercolationLattice{N}(p) where {N}
    return Lattice{N}(LatticeVar{N}(BernoulliBool(p)); lattice_type = SitePercolation)
end

SitePercolationLattice(p) = SitePercolationLattice{1}(p)

_getindex(dist::Distribution, inds...) = rand(dist)
_getindex(lattice_var::AbstractLatticeVar, inds...) = lattice_var[inds...]
Base.getindex(lattice::Lattice, i::Integer) = lattice.vars[i]
Base.getindex(lattice::Lattice, i::Integer, inds::Integer...) = _getindex(lattice[i], inds)
Base.getindex(lattice::Lattice, i::Integer, p::Point) = _getindex(lattice[i], get_coords(p))
init!(lattice::Lattice) = foreach(x -> init!(x), lattice.vars)
get_vars(lattice::Lattice) = lattice.vars
get_num_sites_visited(lattice::Lattice) = get_num_sites_visited(get_vars(lattice)[1])

function Base.show(io::IO, lattice::Lattice{N}) where N
    println(io, N, "-dimensional Lattice")
    for v in get_vars(lattice)
        println(io, v)
    end
end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API for the orbit propagators in SatelliteToolbox.jl.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export get_epoch, get_mean_elements, init_orbit_propagator, propagate!,
       propagate_to_epoch!, step!

"""
    init_orbit_propagator(T, args...; kwargs...)

Initialize the orbit propagator of type `T`. The arguments `args` and keywords
`kwargs` depends of the propagator type.

"""
init_orbit_propagator

"""
    get_epoch(orbp)

Return the epoch of the propagator `orbp` [JD].

"""
get_epoch

"""
    get_mean_elements(orbp)

Return the mean elements of the latest propagation performed by `orbp`. This is
an optinal function in the API. It will return `nothing` if the propagator does
not support it.

"""
get_mean_elements(orbp::OrbitPropagator) = nothing

"""
    propagate!(orbp::OrbitPropagator{T}, t::Number) where T
    propagate!(orbp::OrbitPropagator{T}, t::AbstractVector) where T

If `t` is a number, then propagate `orbp` by `t` [s] from the orbit epoch.
Otherwise, if `t` is an array, then propagate the orbit in `orbp` using the time
instants defined in the vector `t` [s].

In both cases, the orbit propagator algorithm is the one related to the
structure `orbp`.

# Returns

* The position vector represented in inertial frame in each time instant [m].
* The velocity vector represented in inertial frame in each time instant [m].

If `t` is an array, then those values will be an array containing the
information related to each epoch in `t`.

"""
propagate!

function propagate!(orbp::OrbitPropagator{T}, t::AbstractVector) where T
    # Output.
    num_elems  = length(t)
    result_orb = Vector{KeplerianElements{T}}(undef, num_elems)
    result_r   = Vector{SVector{3,T}}(undef, num_elems)
    result_v   = Vector{SVector{3,T}}(undef, num_elems)

    @inbounds for k = 1:num_elems
        # Propagate the orbit.
        r_i_k, v_i_k = propagate!(orbp, t[k])
        result_r[k]  = r_i_k
        result_v[k]  = v_i_k
    end

    return result_r, result_v
end

"""
    propagate_to_epoch!(orbp::OrbitPropagator{T}, JD::Number) where T
    propagate_to_epoch!(orbp::OrbitPropagator{T}, JD::AbstractVector) where T

If `t` is a number, then propagate `orbp` until the epoch `JD` [Julian Day].
Otherwise, if `JD` is an array, then propagate the orbit in `orbp` using the
epochs defined in the vector `t` [Julian Day].

In both cases, the orbit propagator algorithm is the one related to the
structure `orbp`.

# Returns

* The position vector represented in inertial frame in each time instant [m].
* The velocity vector represented in inertial frame in each time instant [m].

If `JD` is an array, then those values will be an array containing the
information related to each epoch in `JD`.

"""
propagate_to_epoch!(orbp::OrbitPropagator,
                    JD::Union{Number, AbstractVector}) where T =
    propagate!(orbp, (JD .- get_epoch(orbp))*86400)

"""
    step!(orbp::OrbitPropagator{T}, Δt::Number)

Propagate the orbit in `orbp` by `Δt` [s] using the algorithm of `orbp`. The new
parameters will be written in `orbp`.

# Returns

* The position vector represented in the inertial frame after the step [m].
* The velocity vector represented in the inertial frame after the step [m].

"""
step!

################################################################################
#                              Iterator interface
################################################################################

# There functions allow broadcast when using the orbit propagators.
Base.iterate(orbp::OrbitPropagator) = (orbp, nothing)
Base.iterate(orbp::OrbitPropagator, ::Nothing) = nothing
Base.length(orbp::OrbitPropagator) = 1
Base.eltype(orbp::T) where T<:OrbitPropagator = T


text:
using Distributions

interaction_dict = Dict()

interaction_dict["binary"] = [-1.0,1.0]
interaction_dict["normal"] = Normal()
interaction_dict["uniform"] = Uniform(-1.0,1.0)


text:
###########################################################################
# Configuration
###########################################################################
D = 2     # choice size to decide from
final = 3 # approximately how many ideas to remain in the end

###########################################################################
# Actual Program
###########################################################################
include("tools.jl")
using Pkg
Pkg.activate(@__DIR__)
using BSON, DelimitedFiles, Random
using REPL.TerminalMenus

function main(D, final)

    ideas = sort!(vec(readdlm("ideas.txt", '\n', String)))
    N = length(ideas)
    rounds = round(Int, log(D, N/final))
    remaining = copy(ideas)

    for r in 1:rounds
        println("Round n. $r. Choose your ideas...")
        RN = length(remaining)
        rp = randperm(RN)
        todel = Int[]
        c = 0
        while c < RN
            j = min(RN - c, D)
            if j > 1
                selection = [remaining[rp[i+c]] for i in 1:j]
                si = request(RadioMenu(selection))
                if si == -1
                    error("Selection cancelled.")
                end
                selected = selection[si]
                append!(todel, [rp[c+sj] for sj in 1:j if sj ≠ si])
            else
                println("Only option remaining is: $(remaining[rp[j+c]]). "*
                "Accepted automatically.")
            end
            c += j
            move_cursor_up_while_clearing_lines(stdout, D)
        end
        deleteat!(remaining, sort!(todel))
    end

    println("Your final selections are:")
    println(remaining)
    result_idxs = [findfirst(isequal(res), ideas) for res in remaining]

    println("\nUpdating statistics now...")
    data = update_data(result_idxs, ideas)

    println("\nCurrent frequencies of your ideas:")
    sp = sortperm(data[:counts], rev = true)
    maxl = maximum(length(i) for i in ideas)
    for s ∈ sp
        idea = data[:ideas][s]
        cou = data[:counts][s]
        println(rpad(idea*":", maxl + 3), cou)
    end
    return data
end

main(D, final)

max_stars_repo_name: JuliaTagBot/WhatDoYouLike

max_issues_repo_name: JuliaTagBot/WhatDoYouLike

max_forks_repo_name: JuliaTagBot/WhatDoYouLike


text:
"""
Authors : Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020
"""

@doc """
Module contenant les types relatifs aux instances et aux solutions du problème de lot sizing avec contraintes environnementales "Green Lot Sizing", que l'on note G_ULS.
Il contient les types et les fonctions suivants:

- Instance
- Solution

Utilisez using .G_ULS pour utiliser le module et son contenu dans votre code.

Utiliser ?NOM dans l'interface en ligne de commande Julia pour avoir de la documentation sur NOM.
Par exemple ?Instance affiche la documentation sur le type Instance

""" ->
module G_ULS

  export Instance, Solution

  @doc """
  Type représentant une entrée du problème G_ULS.

  Il n'est pas nécessaire de construire vous même les entrées.
  Utilisez le fichier generate.jl pour cela.

  """
  mutable struct Instance
    T::Int # horizon
    M::Int # nombre de modes
    R::Int # taille de l'intervalle glissant
    Emax::Array{Int} # impact environnementale max
    f::Array{Array{Int}} # cout d'approvisionnement fixe
    e::Array{Array{Int}} # impact environnementale
    h::Array{Int} # cout de stockage unitaire
    p::Array{Array{Int}} # cout d'approvisionnement unitaire
    d::Array{Int} # demande
  end
  Instance() = Instance(0, 0, 0, [], [], [], [], [], []) # Constructeur par défaut

  @doc """
  Type représentant une solution du problème G_ULS.

  Vous n'avez pas à construire une solution à la main. Vous pouvez utiliser :
  - Solution(inst::Instance)

  """
  mutable struct Solution
    inst
    x
    y
    s
  end
  Solution(inst) = Solution(inst,[],[],[]) # Constructeur par défaut

end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
#=
n! means n × (n − 1) × ... × 3 × 2 × 1

For example, 10! = 10 × 9 × ... × 3 × 2 × 1 = 3628800,
and the sum of the digits in the number 10! is 3 + 6 + 2 + 8 + 8 + 0 + 0 = 27.

Find the sum of the digits in the number 100!
=#
@time println(sum(digits(factorial(BigInt(100)))))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler