text:
module Graphs

export
    # Types
    NodeID,
    SequenceDistanceGraphLink,
    SDGLink,
    SequenceDistanceGraph,
    SDG,
    SequenceDistanceGraphPath,

    # Nodes and sequences
    name,
    n_nodes,
    sequence,
    each_node_id,

    # Graph topology
    source,
    destination,
    distance,
    is_forwards_from,
    is_backwards_from,
    find_link,
    forward_links,
    backward_links,
    get_next_nodes,
    get_previous_nodes,
    find_tip_nodes,
    find_tip_nodes!,
    find_all_unitigs,
    find_all_unitigs!,

    # IO
    load_from_gfa1!,
    write_to_gfa1

using BioSequences, FASTX

###
### Node type for SequenceDistanceGraph
###

const NodeID = Int64

"""
The SDGNode type represents a node in a SequenceDistanceGraph.

At present it contains only two fields, first it holds an instance of a BioSequences
BioSequence type.

Secondly, it tracks a flag which indicates if the node has been
deleted or not.

!!! note
    The deleted flag allows us to mark nodes in the graph as deleted, which can be of
    help in some algorithms where the graph structure is being edited (merging nodes for example).

    Actually deleting the node would shift node IDs and require redoing all links in the graph
    and so on.

    So just marking a node as deleted and not using it anymore is a lazy but sometimes
    helpful choice.
"""
struct SDGNode{S<:BioSequence}
    seq::S
    deleted::Bool
end

# TODO: This is a hacked copy of the dna string literal macro from BioSequences,
# except it creates 2-bit based DNA sequences rather than 4 bit based ones.
# This ability to choose the bit encoding should make its way to BioSequences.jl
# in the future, but for now, it's here.
#
# I basically want this as it lets me create a single literal empty sequence, shared
# by all deleted SDG nodes. Rather than having each deleted SDG node create a new empty
# sequence.
macro dna2_str(seq, flag)
    if flag == "s"
        return LongSequence{DNAAlphabet{2}}(BioSequences.remove_newlines(seq))
    elseif flag == "d"
        return quote
            LongSequence{DNAAlphabet{2}}($(BioSequences.remove_newlines(seq)))
        end
    end
    error("Invalid DNA flag: '$(flag)'")
end

"""
    empty_seq(::Type{LongSequence{DNAAlphabet{2}}})

Get a reference to an empty 2-bit long dna sequence.

In this module, a single empty sequence is allocated, and every call to
`empty_seq` returns a reference to that same allocated sequence.

This saves on allocations if you are deleting many nodes in a
SequenceDistanceGraph.
"""
@inline empty_seq(::Type{LongSequence{DNAAlphabet{2}}}) = dna2""s

"""
empty_seq(::Type{LongSequence{DNAAlphabet{4}}})

Get a reference to an empty 4-bit long dna sequence.

In this module, a single empty sequence is allocated, and every call to
`empty_seq` returns a reference to that same allocated sequence.

This saves on allocations if you are deleting many nodes in a
SequenceDistanceGraph.
"""
@inline empty_seq(::Type{LongSequence{DNAAlphabet{4}}}) = dna""s

"Create a SequenceDistanceGraph node with a refence to an empty sequence."
function empty_node(::Type{S}) where {S<:BioSequence}
    return SDGNode{S}(empty_seq(S), true)
end

"Check if a SDGNode is a deleted."
@inline is_deleted(n::SDGNode{S}) where {S<:BioSequence} = n.deleted

"Get the length of an SDGNode's sequence."
@inline Base.length(n::SDGNode{S}) where {S<:BioSequence} = length(n.seq)

"""
    unsafe_sequence(n::SDGNode{S}) where {S<:BioSequence}

Get the **reference** to a node's underlying sequence object.

!!! warning
    You get a reference to the node's sequence object - not a copy.
    So doing transformation operations (reverse_complement, setindex, etc.) to
    it will probably screw up the graph!
"""
@inline unsafe_sequence(n::SDGNode{S}) where {S<:BioSequence} = n.seq

###
### Link type for SequenceDistanceGraph
###

"""
A SequenceDistanceGraphLink or SDGLink represents a single distance between
two sequences in a SequenceDistanceGraph.

Every link describes a connection between two node ends and contains a distance
(they take the form `([+, -]n1, [+, -]n2, [+, -]dist)`).

!!! note
    A link connects two node ends, and so the order of the signed nodes in the
    links does not change the link.

!!! note
    If the distance in a link is negative, this represents an overlap between
    two sequences. These overlaps must be **perfect overlaps**.

!!! note
    Two links are considered equal (`==` will return `true`) if the source and
    the destination fields of the two links are the same, the distance or dist
    field is not considered!
"""
struct SequenceDistanceGraphLink
    source::NodeID
    destination::NodeID
    dist::Int64
end

"Shorthand for SequenceDistanceGraphLink"
const SDGLink = SequenceDistanceGraphLink

SequenceDistanceGraphLink(src::NodeID, dst::NodeID) = SequenceDistanceGraphLink(src, dst, 0)

"Get the source node end of a link."
source(l::SequenceDistanceGraphLink) = l.source

"Get the destination node end of a link."
destination(l::SequenceDistanceGraphLink) = l.destination

"Get the distance between two node ends involved in a link."
distance(l::SequenceDistanceGraphLink) = l.dist

function Base.:(==)(x::SequenceDistanceGraphLink, y::SequenceDistanceGraphLink)
    return source(x) == source(y) && destination(x) == destination(y)
end

"Test if link `l` is a forward link leaving node `n`."
is_forwards_from(l::SequenceDistanceGraphLink, n::NodeID) = source(l) == -n

"Test if link `l` is a forward link leaving node `n`."
is_backwards_from(l::SequenceDistanceGraphLink, n::NodeID) = source(l) == n

###
### The SequenceDistanceGraph type
###

const LinksVecVec = Vector{Vector{SequenceDistanceGraphLink}}

"""
The SequenceDistanceGraph is a representation of a genome assembly.

The SequenceDistanceGraph stores sequences in a vector of
`SequenceDistanceGraphNode`s (or `SDGNode` for short).

A singe node represents a sequence *and* its reverse complement.

Every node also has an orientaton: Each node has a positive end (+),
and a negative end (-).
So when a node is accessed with (or traversed by entering) the positive end the
node yields the stored sequence.
Conversely, when a node is accessed with (or traversed by entering) the negative
end the node yelds the reverse complement of the stored sequence.
In this way the positive end can be thought of as the sequence start of the
node's canonical sequence, and the negative end can be thought of as its end.

Thankfully, you don't need to think about "positive node ends" and "negative
node ends" most of the time you interact with a `SequenceDistanceGraph`:
In a `SequenceDistanceGraph`, every has a correlative ID starting from 1.
For every node stored in the graph with an ID X, the negative ID -X is mapped to
the reverse complement of that node.
This mapping is virtual: Only one node is stored in the graph.
To see more about how these positive and negative node IDs work, refer to the
docs for methods such as `get_next_nodes` and `get_previous_nodes`. These high
level functions simplify coding graph traversal A LOT.

This `SequenceDistanceGraph` type is only intended to be interacted with
directly by developers and users who know what they are doing. Most tasks an end
user wants to do to manipulate or query a graph can be achieved through the
`WorkSpace` type and the NodeView interface.
"""
struct SequenceDistanceGraph{S<:BioSequence}
    nodes::Vector{SDGNode{S}}
    links::LinksVecVec
end

"Shorthand for SequenceDistanceGraph"
const SDG = SequenceDistanceGraph

"Construct an empty sequence distance graph."
function SequenceDistanceGraph{S}() where {S<:BioSequence}
    return SequenceDistanceGraph{S}(Vector{SDGNode{S}}(), LinksVecVec())
end

struct SequenceDistanceGraphPath{G<:SequenceDistanceGraph}
    sg::G
    nodes::Vector{NodeID}
end

Base.summary(io::IO, sdg::SequenceDistanceGraph) = print(io, "Sequence distance graph (", n_nodes(sdg), " nodes)")
Base.show(io::IO, sdg::SequenceDistanceGraph) = println(io, summary(sdg))

## Internal / not-nessecerily-safe

### Nodes and sequences

"""
    check_node_id(sg::SequenceDistanceGraph, i::NodeID)

The method responsible for checking that a node id used as input to a method
that queries or edits a `SequenceDistanceGraph` is a sensible value.
"""
@inline function check_node_id(sg::SequenceDistanceGraph, i::NodeID)
    if 0 < abs(i) ≤ n_nodes(sg)
        return true
    end
    @error "Sequence graph has no node with ID of $i"
end

"""
Get a **reference** to the vector of nodes in a graph `sg`.

!!! warning
    It is a bad idea to edit this vector yourself unless you know what you are
    doing.
"""
@inline nodes(sg::SequenceDistanceGraph) = sg.nodes

"""
    node_unsafe(sg::SequenceDistanceGraph, n::NodeID)

Get a reference specific node from a sequence distance graph `sg` using its
correlative node id `n`.

!!! note
    `node_unsafe` accepts a NodeID that can be positive or negative.
    E.g. providing either 5 or -5 both mean node 5 in a graph,
    and so you will get the node for node 5.

!!! warning
    This method returns a **reference** to a graph's underlying `SDGNode`.
    NOT a copy! Messing with this will screw up your graph. So this method is
    not recommended or exported - unless you really know what you're doing.

!!! warning
    This method is explicitly marked unsafe for a reason.
    It does *zero* checking of the value it is passed as the node id.
    It also makes use of the `@inbounds` macro.
    This makes it faster, but you need to be 100% sure that your code won't try
    to call this with a bad id. So it is not exported or recommended
    - unless you really know what you're doing.
"""
@inline node_unsafe(sg::SequenceDistanceGraph, n::NodeID) = @inbounds nodes(sg)[abs(n)]

"""
    node(sg::SequenceDistanceGraph, n::NodeID)

Get a **reference** to a specific node from a sequence distance graph `sg` using
its correlative node id `n`.

!!! note
    `node` accepts a NodeID that can be positive or negative.
    E.g. providing either 5 or -5 both mean node 5 in a graph,
    and so you will get the links for node 5.

!!! warning
    This method returns a **reference** to a graph's underlying `SDGNode`.
    NOT a copy! Messing with this will screw up your graph. So this method is
    not recommended or exported - unless you really know what you're doing.
"""
@inline function node(sg::SequenceDistanceGraph, n::NodeID)
    check_node_id(sg, n)
    return node_unsafe(sg, n)
end

"""
    unsafe_sequence(sg::SequenceDistanceGraph, n::NodeID)

Get the **reference** to a node's underlying sequence object.

!!! warning
    This method is unsafe, no checking of node id's occurs and you get a
    reference to the node's sequence object - not a copy, so doing transformation
    operations (reverse_complement, setindex, etc.) to it will probably screw up
    the graph!
"""
@inline unsafe_sequence(sg::SequenceDistanceGraph, n::NodeID) = unsafe_sequence(node_unsafe(sg, n))

### Graph topology

"""
Get a **reference** to the vector of vectors of links in a graph `sg`.

!!! warning
    It is a bad idea to edit this vector yourself unless you know what you are
    doing.
"""
@inline links(sg::SequenceDistanceGraph) = sg.links

"""
    linksof_unsafe(sg::SequenceDistanceGraph, n::NodeID)

Get a **reference** to a vector storing all the links of a node in a
`SequenceDistanceGraph`, the node is specified using its correlative node id `n`.

!!! note
    `links_unsafe` accepts a NodeID that can be positive or negative.
    E.g. providing either 5 or -5 both mean node 5 in a graph,
    and so you will get the links for node 5.

!!! warning
    This method returns a **reference** to an underlying links vector that the
    `SequenceDistanceGraph` owns - NOT a copy!
    Messing with this will screw up your graph. So this method is not
    recommended or exported - unless you really know what you're doing.

!!! warning
    This method is explicitly marked unsafe for a reason.
    It does *zero* checking of the value it is passed as the node id.
    It also makes use of the `@inbounds` macro.
    This makes it faster, but you need to be 100% sure that your code won't try
    to call this with a bad id. So it is not exported or recommended
    - unless you really know what you're doing.
"""
@inline linksof_unsafe(sg::SequenceDistanceGraph, n::NodeID) = @inbounds links(sg)[abs(n)]

"""
    linksof(sg::SequenceGraph, n::NodeID)

Get a **reference** to a vector storing all the links of a node in a
`SequenceDistanceGraph`, the node is specified using its correlative node id `n`.

!!! note
    `links` accepts a NodeID that can be positive or negative.
    E.g. providing either 5 or -5 both mean node 5 in a graph,
    and so you will get the links for node 5.

!!! warning
    This method returns a **reference** to an underlying links vector that the
    `SequenceDistanceGraph` owns - NOT a copy!
    Messing with this will screw up your graph. So this method is not
    recommended or exported - unless you really know what you're doing.
"""
@inline function linksof(sg::SequenceDistanceGraph, n::NodeID)
    check_node_id(sg, n)
    return linksof_unsafe(sg, n)
end

### Graph editing and manipulation

"""
    add_node!(sg::SequenceDistanceGraph{S}, n::SDGNode{S}) where {S<:BioSequence}

Add a node `n` to a graph `sg`.

Returns the node ID used to access the new node added in the graph.

!!! warning
    We don't enforce the sequence in the node is canonical here. We just trust
    that it is canonical.

!!! warning
    Adding a node to the graph does just that. After adding the node it still
    will not be linked to any other nodes.
"""
function add_node!(sg::SequenceDistanceGraph{S}, n::SDGNode{S}) where {S<:BioSequence}
    newlen = length(push!(nodes(sg), n))
    push!(links(sg), Vector{SequenceDistanceGraphLink}())
    return newlen
end

"""
   add_node!(sg::SequenceDistanceGraph{S}, seq::BioSequence) where {S<:BioSequence}

Add a sequence to a sequence distance graph as a node.

Returns the node ID used to access the new node added in the graph.

Can accept any sequence type and will attempt to coerce the input sequence to the
type required by the graph.

!!! warning
    We don't enforce the sequence in the node is canonical here. We just trust
    that it is canonical.

!!! warning
    Adding a node to the graph does just that. After adding the node it still
    will not be linked to any other nodes.
"""
function add_node!(sg::SequenceDistanceGraph{S}, seq::BioSequence) where {S<:BioSequence}
    return add_node!(sg, SDGNode{S}(convert(S, seq), false))
end

"""
    remove_node!(sg::SequenceDistanceGraph{S}, n::NodeID) where {S<:BioSequence}

Remove a node from a sequence distance graph.

!!! note
    This method can accepts a NodeID that can be positive or negative.
    E.g. providing either 5 or -5 both mean node 5 in a graph,
    and so you will end up deleting node 5.

!!! note
    Links involving this node will also be removed from the graph.
"""
function remove_node!(sg::SequenceDistanceGraph{S}, n::NodeID) where {S<:BioSequence}
    oldlinks = copy(linksof(sg, n))
    for oldlink in oldlinks
        remove_link!(sg, source(oldlink), destination(oldlink))
    end
    # TODO: This is a lazy solution to getting rid of the node.
    @inbounds nodes(sg)[abs(n)] = empty_node(S)
end

"""
    add_link!(sg::SequenceDistanceGraph, source::NodeID, dest::NodeID, dist::Int)

Construct a link between two node ends in a sequence Graph.
"""
function add_link!(sg::SequenceDistanceGraph, source::NodeID, dest::NodeID, dist::Int)
    push!(linksof(sg, source), SequenceDistanceGraphLink(source, dest, dist))
    push!(linksof(sg, dest), SequenceDistanceGraphLink(dest, source, dist))
end

"""
    remove_link!(sg::SequenceDistanceGraph, src::NodeID, dst::NodeID)

Remove a link between two nodes in a SequenceDistanceGraph.
Returns a boolean indicating whether the removal was successful.
Reasons this function would not return `true` include that the link
didn't exist in the graph, and so could not be removed.
"""
function remove_link!(sg::SequenceDistanceGraph, src::NodeID, dst::NodeID)
    slinks = linksof(sg, src)
    slinkslen = length(slinks)
    filter!(!isequal(SequenceDistanceGraphLink(src, dst, 0)), slinks)
    dlinks = linksof(sg, dst)
    dlinkslen = length(dlinks)
    filter!(!isequal(SequenceDistanceGraphLink(dst, src, 0)), dlinks)
    return slinkslen != length(slinks) || dlinkslen != length(dlinks)
end

remove_link!(sg::SequenceDistanceGraph, lnk::SequenceDistanceGraphLink) = remove_link!(sg, source(lnk), destination(lnk))

"""
Removes all the links in the collection from and to a given nodeID.
"""
function disconnect_node!(sg::SequenceDistanceGraph, n::NodeID)
    for flink in forward_links(sg, n)
        remove_link!(sg, flink)
    end
    for rlink in backward_links(sg, n)
        remove_link!(sg, rlink)
    end
end

"""
    collapse_all_unitigs!(unitigs::Vector{SequenceDistanceGraphPath{G}}, newnodes::Vector{NodeID}, sg::G, min_nodes::Integer, consume::Bool) where {G<:SequenceDistanceGraph}

Detects all of the trivial paths through the graph that define unitigs. Such
paths are defined as a chain of nodes with only one neighbour. Each such simple
paths through the graph can safely be collapsed into one larger node, which is
what this function does.

!!! note
    This method will edit `unitigs` so as after this method returns it contains
    All the paths that were detected as unitigs in the graph.
    This method also edits the `newnodes` vector, so as after this method returns
    it contains the IDs of all new nodes that were created as a result of
    collapsing the paths in `unitigs`.

!!! note
    Modifies the SequenceDistanceGraph `sg`.
"""
function collapse_all_unitigs!(unitigs::Vector{SequenceDistanceGraphPath{G}},
                               newnodes::Vector{NodeID},
                               sg::G,
                               min_nodes::Integer,
                               consume::Bool) where {G<:SequenceDistanceGraph}

    find_all_unitigs!(unitigs, sg, min_nodes)
    resize!(newnodes, length(unitigs))
    @inbounds for i in eachindex(unitigs)
        newnodes[i] = join_path!(unitigs[i], consume)
    end
end

"""
    collapse_all_unitigs!(sg::SequenceDistanceGraph, min_nodes::Integer, consume::Bool)

Detects all of the trivial paths through the graph that define unitigs. Such
paths are defined as a chain of nodes with only one neighbour. Each such simple
paths through the graph can safely be collapsed into one larger node, which is
what this function does.

!!! note
    Modifies the SequenceDistanceGraph `sg`.
"""
function collapse_all_unitigs!(sg::SequenceDistanceGraph, min_nodes::Integer, consume::Bool)
    unitigs = Vector{SequenceGraphPath{typeof(sg)}}()
    newnodes = Vector{NodeID}()
    return collapse_all_unitigs!(unitigs, newnodes, sg, min_nodes, consume)
end

## Public / safe

"Get the name of the graph. Defaults to the symbol :sdg."
@inline name(sg::SequenceDistanceGraph) = :sdg

### Nodes and sequences

"Get the number of nodes in the sequence distance graph `sg`."
@inline n_nodes(sg::SequenceDistanceGraph) = length(nodes(sg))

"Iterate over every node ID in the sequence distance graph `sg`."
@inline each_node_id(sg::SequenceDistanceGraph) = eachindex(nodes(sg))

"""
    sequence(sg::SequenceDistanceGraph, n::NodeID)

Get the full sequence of a node in a sequence distance graph using its
correlative node id `n`.

!!! note
    `sequence` accepts a NodeID that can be positive or negative.
    Nodes represent stretches of sequence in a canonical orientation, if you ask
    for for the sequence of say the third node, the positive node id 3
    (which denotes traversing the third node in the forward direction),
    gives you the canonical sequence. If you use the negative ID -3
    (which denotes traversing the third node in the reverse direction), you will
    get the reverse complement of the node's canonical (forward) sequence.

!!! note
    It is safe to modify the returned sequence without screwing up your graph,
    yet thanks to BioSequences.jl's copy on write system for LongSequences, data
    copying will only occur if nessecery. You get the best of both worlds.
"""
function sequence(sg::SequenceDistanceGraph, n::NodeID)
    check_node_id(sg, n)
    seqref = unsafe_sequence(sg, n)
    outseq = typeof(seqref)(seqref, 1:lastindex(seqref))
    if n < 0
        reverse_complement!(outseq)
    end
    return outseq
end

### Graph Topology

"""
    find_link(sg::SequenceDistanceGraph, src::NodeID, dst::NodeID)

Find and return the link that exists between a source node and a destination
node, using their correlative node ids.

In this instance, the IDs also denote the "ends" of a node, a link connects.

Recall that every node has an orientaton: Each node has a positive end (+),
and a negative end (-).
So when a node is accessed with (or traversed by entering) the positive end the
node yields the stored sequence.
Conversely, when a node is accessed with (or traversed by entering) the negative
end the node yelds the reverse complement of the stored sequence.

Thus `find_link(sg, -5, 1)` means you want to find a link in the graph that
allows you to exit the (-) end of node 5 (meaning you just traversed it in the
canonical + --> - orientation), and enter the (+) end of node 1 (meaning you will
also traverse node one in the canonical + --> - orientation).

By contrast `find_link(sg, -5, -1)` means you want to find a link in the graph
that allows you to exit the (-) end of node 5 - again having just traverse it in
the + --> - orientation, but want to enter node 1 through it's (-) end - and thus
traversing node 1 in the + <-- - orientation, yielding the reverse complement of
the canonical sequence of node 1.

If a link connecting the two node ends is present, you get it. If not, you get
`nothing`. So make sure to check the output.

!!! note
    This method is safe and public, but not reeeeaaallly intended for the end
    user, as they have to worry about links between "ends" of a node.
    In contrast to the rest of this framework, which rather has the user thinking
    in node centric terms.
    e.g. "I am at node 5 which means I have traversed it forwards or I am at
    node -5, which means I have traversed node 5 backwards. From here I can visit
    node 1, -3, and 2, so I may go through nodes 1 and 2 forwards, and through 3
    backwards. The methods `get_next_nodes` and `get_previous_nodes` are more
    user-friendly and follow this more node-centric mindset.
"""
function find_link(sg::SequenceDistanceGraph, src::NodeID, dst::NodeID)
    query = SequenceDistanceGraphLink(src, dst)
    for l in linksof(sg, src)
        if l == query
            return l
        end
    end
    return nothing
end

"""
    forward_links(sg::SequenceDistanceGraph, n::NodeID)

Get a vector of the links that are ahread of you, as you traverse node `n`,
continuing forward in you present direction of travel (`n` can be a positive or
negative node id).

The node id can be positive or negative. For example if you use a positive ID,
such as 5, this means you are traversing node 5 in the canonical orientation,
and so you will get the links that allow you to leave node 5 in the canonical
direction. If you used a negative ID such as -5, that means you are traversing
node 5 in the non-canonical orientation, and so you will get the links that
allow you to leave node 5 in the non-canonical direction.
"""
function forward_links(sg::SequenceDistanceGraph, n::NodeID)
    r = Vector{SequenceDistanceGraphLink}()
    nodelinks = linksof(sg, n)
    sizehint!(r, length(nodelinks))
    for link in nodelinks
        if is_forwards_from(link, n)
            push!(r, link)
        end
    end
    return r
end

"""
    backward_links(sg::SequenceDistanceGraph, n::NodeID)

Get a vector of the links that are behind of you, as you traverse node `n`,
continuing forward in you present direction of travel (`n` can be a positive or
negative node id).

The node id can be positive or negative. For example if you use a positive ID,
such as 5, this means you are traversing node 5 in the canonical orientation,
and so you will get the links that would allow you to enter node 5 in the
canonical direction. If you used a negative ID such as -5, that means you are
traversing node 5 in the non-canonical orientation, and so you will get the links
that would allow you to enter node 5 in the non-canonical direction.
"""
backward_links(sg::SequenceDistanceGraph, n::NodeID) = forward_links(sg, -n)

"""
    get_next_nodes(sg::SequenceDistanceGraph, n::NodeID)

Get the nodes you may visit next as you exit node `n`, maintaining your current
direction of travel (`n` can be a positive or negative node ID).

!!! note
    The node id can be positive or negative. For example if you use a positive ID,
    such as 5, this means you are traversing node 5 in the canonical orientation,
    and so you will get the nodes you may visit as you leave node 5 in the canonical
    direction. If you use a negative ID, such as -5, this means you are
    traversing node 5 in the non-canonical orientation, and so you will get the
    nodes you may visit as you leave node 5 in the non-canonical direction.

!!! note
    The list of nodes returned are also signed to denote direction: Say
    you got `[1, -2, 3]` as a result of `get_next_nodes(sg, 5)`. That means as
    you leave node 5 after traversing it in the canonical direction, you may
    proceed to travel through nodes 1 and 3 in the canonical direction, and node
    2 in the non-canonical direction.
"""
function get_next_nodes(sg::SequenceDistanceGraph, n::NodeID)
    r = Vector{NodeID}()
    nodelinks = linksof(sg, n)
    sizehint!(r, length(nodelinks))
    for link in nodelinks
        if is_forwards_from(link, n)
            push!(r, destination(link))
        end
    end
    return r
end

"""
    get_previous_nodes(sg::SequenceDistanceGraph, n::NodeID)

Get the nodes you may have previously been on before you enterd node `n`,
maintaining your current direction of travel (`n` can be a positive or negative
node ID).

!!! note
    The node id can be positive or negative. For example if you use a positive ID,
    such as 5, this means you are traversing node 5 in the canonical orientation,
    and so you will get the nodes you may have visited prior to entering node 5
    in the canonical direction. If you use a negative ID, such as -5, this means
    you are traversing node 5 in the non-canonical orientation, and so you will
    get the nodes you may have visited prior to entering node 5 in the
    non-canonical direction.

!!! note
    The list of nodes returned are also signed to denote direction: Say
    you got `[1, -2, 3]` as a result of `get_previous_nodes(sg, 5)`. That means
    prior to entering node 5 in the canonical direction, you may have traveled
    through nodes 1 and 3 in the canonical direction, and node 2 in the
    non-canonical direction.
"""
function get_previous_nodes(sg::SequenceDistanceGraph, n::NodeID)
    r = Vector{NodeID}()
    nodelinks = linksof(sg, n)
    sizehint!(r, length(nodelinks))
    for link in nodelinks
        if is_backwards_from(link, n)
            push!(r, -destination(link))
        end
    end
    return r
end

"""
    find_tip_nodes!(result::Set{NodeID}, sg::SequenceDistanceGraph, min_size::Integer)

Get a set of IDs of all the nodes in the graph `sg` that count as "tips".

Here a tip node is defined as a node that only has a single neighbouring node
at one of it's ends, and no neighbouring nodes out of one of its other end.

The node's sequence must also be larger than `min_size` base pairs in length.

!!! note
    This method modifys a `result` input. So if you want to find tips in a
    graph repeatedly, you can resuse a `Set`, saving you some allocation costs.
"""
function find_tip_nodes!(result::Set{NodeID}, sg::SequenceDistanceGraph, min_size::Integer)
    empty!(result)
    for n in each_node_id(sg)
        nd = node(sg, n)
        if is_deleted(nd) || length(nd) > min_size
            continue
        end
        fwl = forward_links(sg, n)
        bwl = backward_links(sg, n)
        if length(fwl) == 1 && length(bwl) == 0
            if length(backward_links(sg, destination(first(fwl)))) > 1
                push!(result, n)
            end
        end
        if length(fwl) == 0 && length(bwl) == 1
            if length(forward_links(sg, -destination(first(bwl)))) > 1
                push!(result, n)
            end
        end
        if isempty(fwl) && isempty(bwl)
            push!(result, n)
        end
    end
    return result
end

"""
    find_tip_nodes(sg::SequenceDistanceGraph, min_size::Integer)

Get a set of IDs of all the nodes in the graph `sg` that count as "tips".

Here a tip node is defined as a node that only has a single neighbouring node
at one of it's ends, and no neighbouring nodes out of one of its other end.

The node's sequence must also be larger than `min_size` base pairs in length.
"""
function find_tip_nodes(sg::SequenceDistanceGraph, min_size::Integer)
    return find_tip_nodes!(Set{NodeID}(), sg, min_size)
end

"""
    find_all_unitigs!(unitigs::Vector{SequenceDistanceGraphPath{G}}, sg::G, min_nodes::Integer) where {G<:SequenceDistanceGraph}

Find and return a vector of paths through the graph that represent all the
unitigs or transitive paths in the graphs. Such paths are defined as a chain
of nodes with only one neighbour. Such simple regions of the graph can safely be
collapsed into one larger node.

!!! note
    This vector modifies an input `unitigs` vector to contain the result.
    This is useful for situations where you want to repeatedly find unitigs
    in the graph to save on additional allocations.
"""
function find_all_unitigs!(unitigs::Vector{SequenceDistanceGraphPath{G}},
    sg::G, min_nodes::Integer) where {G<:SequenceDistanceGraph}
    empty!(unitigs)
    consumed = falses(n_nodes(sg))
    for n in each_node_id(sg)
        if consumed[n] || is_deleted(node(sg, n))
            continue
        end
        consumed[n] = true
        path = SequenceDistanceGraphPath(sg, [n])

        # Two passes, fw and bw, path is inverted twice, so still n is +
        for pass in 1:2
            fn = forward_links(sg, last(path))
            while length(fn) == 1
                dest = destination(first(fn))
                if (!consumed[abs(dest)]) && (length(backward_links(sg, dest)) == 1)
                    push!(path, dest)
                    consumed[abs(dest)] = true
                else
                    break
                end
                fn = forward_links(sg, last(path))
            end
            reverse!(path)
        end
        if n_nodes(path) >= min_nodes
            push!(unitigs, path)
        end
    end
    return unitigs
end

"""
    find_all_unitigs(sg::G, min_nodes::Integer) where {G<:SequenceDistanceGraph}

Find and return a vector of paths through the graph that represent all the
unitigs or transitive paths in the graphs. Such paths are defined as a chain
of nodes with only one neighbour. Such simple regions of the graph can safely be
collapsed into one larger node.
"""
function find_all_unitigs(sg::G, min_nodes::Integer) where {G<:SequenceDistanceGraph}
    return find_all_unitigs!(Vector{SequenceDistanceGraphPath{G}}(), sg, min_nodes)
end

### IO

"""
    write_to_gfa1(sg::SequenceDistanceGraph, filename::String)

Write a graph from a GFAv1 formatted file, and associated FASTA file of node sequences.

!!! note
    The GFA format permits storing sequences of the graph nodes in a seperate fasta
    file, instead of in the GFA file. This is so as the sequences of the graph
    nodes can be easily fed into other tools that typically accept FASTA files
    as input. Many assemblers also output a GFA + FASTA combo.
    Therefore, this method writes node sequences to a file called `[filename].fasta`,
    and the graph structure to a file called `[filename].gfa`.
"""
function write_to_gfa1(sg::SequenceDistanceGraph, filename::String)
    @info string("Saving graph to ", filename)
    fasta_filename = "$filename.fasta"
    gfa = open("$filename.gfa", "w")
    fasta = open(FASTA.Writer, fasta_filename)
    println(gfa, "H\tVN:Z:1.0")
    for nid in eachindex(nodes(sg))
        n = node(sg, nid)
        if n.deleted
            continue
        end
        println(gfa, "S\tseq", nid, "\t*\tLN:i:", length(n), "\tUR:Z:", fasta_filename)
        write(fasta, FASTA.Record(string("seq", nid), unsafe_sequence(n)))
    end
    close(fasta)
    for ls in links(sg)
        for l in ls
            if source(l) <= destination(l)
                print(gfa, "L\t")
                if source(l) > 0
                    print(gfa, "seq", source(l), "\t-\t")
                else
                    print(gfa, "seq", -source(l), "\t+\t")
                end
                if destination(l) > 0
                    print(gfa, "seq", destination(l), "\t+\t")
                else
                    print(gfa, "seq", -destination(l), "\t-\t")
                end
                println(gfa, distance(l) < 0 ? -distance(l) : 0, "M")
            end
        end
    end
    close(gfa)
end

"""
    load_from_gfa1!(sg::SequenceDistanceGraph{S}, gfafile::AbstractString, fafile::AbstractString) where {S<:BioSequence}

Load a graph from a GFAv1 formatted file, and associated FASTA file of node sequences.

!!! note
    The GFA format permits storing sequences of the graph nodes in a seperate fasta
    file, instead of in the GFA file. This is so as the sequences of the graph
    nodes can be easily fed into other tools that typically accept FASTA files
    as input. Many assemblers also output a GFA + FASTA combo.
    Therefore, this method asks for the filepath of a GFAv1 file, as well as a
    filepath to a FASTA formatted file. This method reads the node sequences from
    the FASTA file, before getting the links between nodes from the GFAv1
    file.
"""
function load_from_gfa1!(sg::SequenceDistanceGraph{S},
                         gfafile::AbstractString,
                         fafile::AbstractString) where {S<:BioSequence}
    @info "Loading graph"
    @info string("Graph FASTA filename: ", fafile)
    # Load all the sequences from FASTA file. If they are not canonical, flip them,
    # and remember that they are flipped.
    @info string("Loading sequences from ", fafile)
    fardr = open(FASTA.Reader, fafile)
    add_nodes!(sg, fardr)
    @info string("Loading links from ", gfafile)
    gfas = open(gfafile, "r")
    for line in eachline(gfas)
        error("Ooops this function is not finished! ;)")
    end
end

###
### SequenceDistanceGraph path
###

SequenceDistanceGraphPath(sg::G) where {G<:SequenceDistanceGraph} = SequenceDistanceGraphPath{G}(sg, Vector{NodeID}())

@inline nodes(p::SequenceDistanceGraphPath) = p.nodes
@inline n_nodes(p::SequenceDistanceGraphPath) = length(nodes(p))
@inline Base.push!(p::SequenceDistanceGraphPath, n::NodeID) = push!(nodes(p), n)
@inline graph(p::SequenceDistanceGraphPath) = p.sg
@inline Base.first(p::SequenceDistanceGraphPath) = first(nodes(p))
@inline Base.last(p::SequenceDistanceGraphPath) = last(nodes(p))

function Base.reverse!(p::SequenceDistanceGraphPath)
    nds = nodes(p)
    i = firstindex(nds)
    j = lastindex(nds)
    while i ≤ j
        @inbounds x = nds[i]
        @inbounds y = nds[j]
        @inbounds nds[i] = -y
        @inbounds nds[j] = -x
        i = i + 1
        j = j - 1
    end
    return p
end

# TODO: See about potential performance improvement using BioSequences.jl `mismatches`.
function check_overlap(a::BioSequence, b::BioSequence, by::Integer)
    i′ = lastindex(a) + 1
    i = i′ - by
    j = 1
    while i < i′
        if a[i] != b[j]
            return false
        end
        i = i + 1
        j = j + 1
    end
    return true
end

# TODO: This function currently does not check for deleted nodes in a graph.
# It probably should throw an error if a node is deleted and you're trying to
# use it in a path.
function sequence(p::SequenceDistanceGraphPath{SequenceDistanceGraph{S}}) where {S<:BioSequence}
    s = S()
    pnode = 0
    for n in nodes(p)
        nseq = sequence(graph(p), n)
        if pnode != 0
            # Find link between pnode output (+pnode) and n's sink (-n)
            lnk = find_link(graph(p), pnode, n)
            if !isnothing(lnk)
                dst = distance(lnk)
                if dst > 0

                else
                    # Check that the overlap is valid in terms of its sequence.
                    ovl = -dst
                    if !check_overlap(s, nseq, ovl)
                        error("Sequences do not overlap.")
                    end
                    # TODO: Investigate if this can be done more efficiently
                    # using resize! and copyto!.... Or perhaps something more
                    # hacky.
                    nseq = nseq[ovl + 1:end]
                end
            else
                error("This path is not a valid path.")
            end
        end
        append!(s, nseq)
        pnode = -n
    end
    return s
end

function join_path!(p::SequenceDistanceGraphPath, consume::Bool)
    pnodes = Set{NodeID}()
    for n in nodes(p)
        push!(pnodes, n)
        push!(pnodes, -n)
    end
    # If p is not canonical, flip it.
    ps = sequence(p)
    if !iscanonical(ps)
        reverse_complement!(ps)
        reverse!(p)
    end
    # Add a new node that is this sequence to the graph.
    newid = add_node!(graph(p), ps)
    # Connect the new node up.
    # TODO: Circular structures in the graph MIGHT screw this up, but I'm not sure.
    for l in backward_links(graph(p), first(p))
        add_link!(graph(p), newid, destination(l), distance(l))
    end
    for l in forward_links(graph(p), last(p))
        add_link!(graph(p), -newid, destination(l), distance(l))
    end

    if consume
        for n in pnodes
            # Check if the node has neighbours that are not involved in the path.
            ext_neigh = false
            if n != last(p)
                for l in forward_links(graph(p), n)
                    if destination(l) ∉ pnodes
                        ext_neigh = true
                    end
                end
            end
            if n != first(p)
                for l in backward_links(graph(p), n)
                    if destination(l) ∉ pnodes
                        ext_neigh = true
                    end
                end
            end
            if ext_neigh
                continue
            end
            remove_node!(graph(p), n)
        end
    end
    return newid
end

end # module

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
module MiniOmega

using Cassette
import Base:~
export sample, unif, normal, pointwise, <|, rt, Ω

const ID = NTuple{N, Int} where N

"Ω maps primitives to values"
struct Ω
  data::Dict{Any, Any}
end

unif(id, ω) = ω[(unif, @show(id))]

"Standard uniform"
unif(ω) = unif(1, ω)

normal(id, ω) = ω[(normal, id)]

"Standard normal"
normal(ω) = normal(1, ω)

Base.getindex(ω::Ω, prim::Tuple{typeof(normal), I}) where I =
  get!(ω.data, prim, randn())

Base.getindex(ω::Ω, prim::Tuple{typeof(unif), I}) where I =
  get!(ω.data, prim, rand())

## Logpdf
logpdf(prim::Tuple{typeof(unif), I}, v) where I = 0.0
# logpdf(prim::Tuple{typeof(unif), I}, v) = 0.0
logpdf(ω::Ω) = sum(logpdf(ωi, vi) for (ωi, vi) in ω.data)

## Sampling

"Sample a random ω ∈ Ω"
sample(::Type{Ω}) = Ω(Dict{Any, Any}())
sample(f) = f(sample(Ω))

# unif(i::ID) = ω -> ω[(i)]

# "Single primitive random variable"
# unif(ω::Ω) = ω[(1,)]

## (Conditional) independence

# Use cassette to augment enivonrment with extra state
Cassette.@context IIDCtx

"""Conditionally independent and identically distributed given `shared` parents

Returns the `id`th element of an (exchangeable) sequence of random variables
that are identically distributed with `f` but conditionally independent given
random variables in `shared`.

This is essentially constructs a plate where all variables in shared are shared,
and all other parents are not.
"""
function iid(f, id::Integer)
  let ctx = IIDCtx(metadata = (id = (id,),))
    ω -> withctx(ctx, f, ω)
  end
end

Base.:~(x, y) = iid(y, x)

withctx(ctx, f, ω) = Cassette.overdub(ctx, f, ω)
function Cassette.overdub(ctx::IIDCtx, ::typeof(withctx), ctxinner, f, ω)
  # Merge the context
  id = (ctxinner.metadata.id..., ctx.metadata.id...)
  Cassette.overdub(IIDCtx(metadata = (id = id,)), f, ω)
end

Primitives = Union{typeof(unif), typeof(normal)}

# Intercept `unif(id, ω)`
Cassette.overdub(ctx::IIDCtx, f::T, id, ω::Ω) where T <: Primitives =
  f(ctx.metadata.id, ω)

## Plates

## A plate is any function of the form `f(id, ω)``

"`i`th member of plate `f`"
struct PlateMember{F, S}
  f::F
  i::S
end

@inline <|(f, i) = PlateMember(f, i)
(m::PlateMember)(ω) = m.f(m.i, ω)

## Syntactic Sugar (to make model-building nicer)

"Random tuple"
rt(fs...) = ω -> map(f -> f(ω), fs)

end

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
import Pkg

modelproto(;kwargs...) = ModelProto(;
    ir_version=7,
    opset_import=[OperatorSetIdProto(version=14)],
    producer_name="ONNX.jl",
    producer_version=string(Pkg.Types.Context().env.project.version), # TODO: Ugh....
    kwargs...)

"""
    graphproto()
Return an [`ONNX.GraphProto`](@ref) with all fields initialized to empty arrays.
"""
graphproto(;kwargs...) = GraphProto(;
    node = NodeProto[],
    initializer = TensorProto[],
    input = ValueInfoProto[],
    output = ValueInfoProto[],
    value_info = ValueInfoProto[],
    kwargs...
)

add!(gp::GraphProto, np::NodeProto) = push!(gp.node, np)

add!(gp::GraphProto, tp::TensorProto) = push!(gp.initializer, tp)

##############################################################################
#                                 Utils                                    #
##############################################################################

# can we make it more robust?
iskwfunc(f) = endswith(string(f), "##kw")

function kwargs2dict(op::Ghost.Call)
    kw = iskwfunc(op.fn) ? op.args[1] : (;)
    return Dict(zip(keys(kw), values(kw)))
end

macro opconfig_kw(backend, fn)
    return quote
        $OpConfig{$backend, <:Union{typeof($fn), typeof(Core.kwfunc($fn))}}
    end
end

function NodeProto(op_type::String, op::Ghost.Call, attrs::Dict=Dict())
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    return NodeProto(
        input=[onnx_name(v) for v in args],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto.(keys(attrs), values(attrs)),
        op_type=op_type
    )
end

ValueInfoProto(op::Ghost.AbstractOp) = ValueInfoProto(
    onnx_name(op),
    # something down the road reverses the shape, so we don't do it here
    # try the following for example:
    #     TypeProto_Tensor((4, 3), Float64).shape.dim[1].dim_value
    # which gives 3 instead of 4
    size(op.val),
    eltype(op.val)
)

##############################################################################
#                                 Methods                                    #
##############################################################################

onnx_name(v::Variable) = "x$(v.id)"
onnx_name(op::Ghost.AbstractOp) = "x$(op.id)"

"""
    save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g::GraphProto, ::OpConfig{:Backend, Fn}, op::Ghost.Call)

Serialize a single operation from a tape to graph.
"""
function save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g, OpConfig{:ONNX, typeof(op.fn)}(), op)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(*)}, op::Ghost.Call)
    nd = NodeProto(
        input=[onnx_name(v) for v in reverse(op.args)],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto[],
        op_type="Gemm"
    )
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, onnx_gemm), op::Ghost.Call)
    kw_dict = kwargs2dict(op)
    attrs = rename_keys(kw_dict, Dict(
        :tA => :transA,
        :tB => :transB,
        :α => :alpha,
        :β => :beta
    ))
    nd = NodeProto("Gemm", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, conv), op::Ghost.Call)
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    w = args[2]._op.val
    # ONNXRuntime gives the following error for Float64:
    # NOT_IMPLEMENTED : Could not find an implementation for the node x3:Conv(11)')
    eltype(w) == Float64 && @warn "Not all ONNX runtimes support input & weights as Float64"
    attrs = from_nnlib_conv(kwargs2dict(op), ndims(w) - 2)
    nd = NodeProto("Conv", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::@opconfig_kw(:ONNX, maxpool), op::Ghost.Call)
    args = iskwfunc(op.fn) ? op.args[3:end] : op.args
    x = args[1]._op.val
    attrs = from_nnlib_conv(kwargs2dict(op), ndims(x) - 2)
    nd = NodeProto("MaxPool", op, attrs)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(add)}, op::Ghost.Call)
    nd = NodeProto("Add", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(mul)}, op::Ghost.Call)
    nd = NodeProto("Mul", op)
    push!(g.node, nd)
end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(relu)}, op::Ghost.Call)
    nd = NodeProto("Relu", op)
    push!(g.node, nd)
end

##############################################################################
#                                    API                                     #
##############################################################################

"""
    save(io::IO, tape::Ghost.Tape{ONNXCtx})
    save(filename::String, tape::Ghost.Tape{ONNXCtx})

Save tape as an ONNX model. The way a particular operation is serialized is
controlled by methods of [save_node!](@ref).

See also: [`load!`](@ref)
"""
function save(io::IO, tape::Tape{ONNXCtx})
    g = graphproto()
    g.name = "generated_model"
    for op in tape
        if op isa Ghost.Input
            # Ghost.Tape represents both - arguments and model parameters
            # as Input ops. In ONNX, parameters (and constants) must be added
            # to .initializer, so will we have to do when we understand how to
            # handle parameters (and complex structs) in general.
            # For reference, the following commented line shows how to add Input
            # to the .initializer:
            # add!(g, TensorProto(op.val, onnx_name(op)))
            push!(g.input, ValueInfoProto(op))
        elseif op isa Ghost.Call
            save_node!(g, op)
        else
            error("$(typeof(op)) is not yet supported in model export")
        end
    end
    push!(g.output, ValueInfoProto(tape[tape.result]))
    m = modelproto();
    m.graph = g;
    writeproto(io, m)
end

function save(filename::String, tape::Tape{ONNXCtx})
    open(filename, "w") do io
        save(io, tape)
    end
end


text:
module TestKernelOptimization
using Test
using LinearAlgebra
using GaussianDistributions
using GaussianDistributions: nlml, optimize, optimize!

@testset "kernel optimization" begin
    function k(θ)
        a, l = exp.(θ)
        (x, y) -> a*exp(-(x-y)^2/(2l^2))
    end
    f(x) = sin(x)
    n = 8
    x = randn(n)
    y = f.(x)
    θ = fill(-2., 2) # WARNING: without noise covariance, this has to lead to full rank system
    k0 = k(θ)
    optimize!(k, θ, x, y)
    @test nlml(k(θ), x, y) < nlml(k0, x, y)

    # using Plots
    # G = Gaussian(k(θ))
    # C = G | (x, y)
    # plot(C)
    # scatter!(x, y)
    # gui()

    # with noisy observations
    σ = 5e-2
    @. y += σ * randn()

    θ = fill(-2., 2) # WARNING: without noise covariance, this has to lead to full rank system
    k0 = k(θ)
    θ = optimize(k, θ, x, y, σ^2)
    @test nlml(k(θ), x, y, σ^2) < nlml(k0, x, y, σ^2)

    # G = Gaussian(k(θ))
    # C = G | (x, y, σ^2)
    # plot(C)
    # scatter!(x, y)
    # gui()
end

@testset "noise variance optimization" begin
    l = .5
    k(x, y) = exp(-(x-y)^2/(2l^2))
    f(x) = sin(2π*x)
    n = 128
    x = randn(n)
    y = f.(x)
    σ_truth = 5e-2     # with noisy observations
    @. y += σ_truth * randn()

    σ²_new = .1
    σ²_new = optimize(σ²_new, k, x, y) # optimizing noise variance
    @test nlml(k, x, y, σ²_new) < nlml(k, x, y, σ_truth^2)
    @test isapprox(σ²_new, σ_truth^2, rtol = .5)

    # using Plots
    # G = Gaussian(k)
    # C = G | (x, y, σ²_new)
    # plot(C)
    # scatter!(x, y)
    # gui()
end

end # TestKernelOptimization


text:
# Automatically generated using Clang.jl

const NVTX_VERSION = 3

# Skipping MacroDefinition: NVTX_INLINE_STATIC inline static

# Skipping MacroDefinition: NVTX_DECLSPEC

# Skipping MacroDefinition: NVTX_VERSIONED_IDENTIFIER_L3 ( NAME , VERSION ) NAME ## _v ## VERSION
# Skipping MacroDefinition: NVTX_VERSIONED_IDENTIFIER_L2 ( NAME , VERSION ) NVTX_VERSIONED_IDENTIFIER_L3 ( NAME , VERSION )
# Skipping MacroDefinition: NVTX_VERSIONED_IDENTIFIER ( NAME ) NVTX_VERSIONED_IDENTIFIER_L2 ( NAME , NVTX_VERSION )

const NVTX_SUCCESS = 0
const NVTX_FAIL = 1
const NVTX_ERR_INIT_LOAD_PROPERTY = 2
const NVTX_ERR_INIT_ACCESS_LIBRARY = 3
const NVTX_ERR_INIT_LOAD_LIBRARY = 4
const NVTX_ERR_INIT_MISSING_LIBRARY_ENTRY_POINT = 5
const NVTX_ERR_INIT_FAILED_LIBRARY_ENTRY_POINT = 6
const NVTX_ERR_NO_INJECTION_LIBRARY_AVAILABLE = 7

# Skipping MacroDefinition: NVTX_EVENT_ATTRIB_STRUCT_SIZE ( ( uint16_t ) ( sizeof ( nvtxEventAttributes_t ) ) )
# Skipping MacroDefinition: NVTX_NO_PUSH_POP_TRACKING ( ( int ) - 2 )
# Skipping MacroDefinition: NVTX_RESOURCE_MAKE_TYPE ( CLASS , INDEX ) ( ( ( ( uint32_t ) ( NVTX_RESOURCE_CLASS_ ## CLASS ) ) << 16 ) | ( ( uint32_t ) ( INDEX ) ) )

const NVTX_RESOURCE_CLASS_GENERIC = 1

# Skipping MacroDefinition: NVTX_RESOURCE_ATTRIB_STRUCT_SIZE ( ( uint16_t ) ( sizeof ( nvtxResourceAttributes_v0 ) ) )

const nvtxRangeId_t = UInt64
const nvtxDomainRegistration_st = Cvoid
const nvtxDomainRegistration = nvtxDomainRegistration_st
const nvtxDomainHandle_t = Ptr{nvtxDomainRegistration}
const nvtxStringRegistration_st = Cvoid
const nvtxStringRegistration = nvtxStringRegistration_st
const nvtxStringHandle_t = Ptr{nvtxStringRegistration}

@cenum nvtxColorType_t::UInt32 begin
    NVTX_COLOR_UNKNOWN = 0
    NVTX_COLOR_ARGB = 1
end

@cenum nvtxMessageType_t::UInt32 begin
    NVTX_MESSAGE_UNKNOWN = 0
    NVTX_MESSAGE_TYPE_ASCII = 1
    NVTX_MESSAGE_TYPE_UNICODE = 2
    NVTX_MESSAGE_TYPE_REGISTERED = 3
end

struct nvtxMessageValue_t
    ascii::Cstring
end

@cenum nvtxPayloadType_t::UInt32 begin
    NVTX_PAYLOAD_UNKNOWN = 0
    NVTX_PAYLOAD_TYPE_UNSIGNED_INT64 = 1
    NVTX_PAYLOAD_TYPE_INT64 = 2
    NVTX_PAYLOAD_TYPE_DOUBLE = 3
    NVTX_PAYLOAD_TYPE_UNSIGNED_INT32 = 4
    NVTX_PAYLOAD_TYPE_INT32 = 5
    NVTX_PAYLOAD_TYPE_FLOAT = 6
end

struct payload_t
    ullValue::UInt64
end

struct nvtxEventAttributes_v2
    version::UInt16
    size::UInt16
    category::UInt32
    colorType::Int32
    color::UInt32
    payloadType::Int32
    reserved0::Int32
    payload::payload_t
    messageType::Int32
    message::nvtxMessageValue_t
end

const nvtxEventAttributes_t = nvtxEventAttributes_v2

@cenum nvtxResourceGenericType_t::UInt32 begin
    NVTX_RESOURCE_TYPE_UNKNOWN = 0
    NVTX_RESOURCE_TYPE_GENERIC_POINTER = 65537
    NVTX_RESOURCE_TYPE_GENERIC_HANDLE = 65538
    NVTX_RESOURCE_TYPE_GENERIC_THREAD_NATIVE = 65539
    NVTX_RESOURCE_TYPE_GENERIC_THREAD_POSIX = 65540
end

struct identifier_t
    ullValue::UInt64
end

struct nvtxResourceAttributes_v0
    version::UInt16
    size::UInt16
    identifierType::Int32
    identifier::identifier_t
    messageType::Int32
    message::nvtxMessageValue_t
end

const nvtxResourceAttributes_t = nvtxResourceAttributes_v0
const nvtxResourceHandle = Cvoid
const nvtxResourceHandle_t = Ptr{nvtxResourceHandle}
const NVTX_RESOURCE_CLASS_CUDA = 4

@cenum nvtxResourceCUDAType_t::UInt32 begin
    NVTX_RESOURCE_TYPE_CUDA_DEVICE = 262145
    NVTX_RESOURCE_TYPE_CUDA_CONTEXT = 262146
    NVTX_RESOURCE_TYPE_CUDA_STREAM = 262147
    NVTX_RESOURCE_TYPE_CUDA_EVENT = 262148
end


text:
using OmegaGrav
using MatterPower

# %% Specify a redshift
redshift = 0

# %% Cosmological parameters
As, ns, kpivot = 2.097e-9, 0.9652, 0.05
Ωm, ΩΛ, Ωb = 0.315, 0.685, 0.049
Ωk = 1 - Ωm - ΩΛ
h0 = 0.674
ωm, fb = Ωm * h0^2, Ωb / Ωm

# %% Tabulate linear growth factor as a function of scale factor, a
sol = setup_growth(Ωm, ΩΛ)
a = 1 / (1 + redshift)
D1 = sol(a)[1]

# %% Define a function to return a linear matter power spectrum (in units of Mpc^3/h^3)
# as a function of the comoving wavenumber, kovh, in units of h/Mpc.
# Here is an example using Einstein & Hu's analytical transfer function
pk(kovh) =
   D1^2 *
   As *
   (kovh * h0 / kpivot)^(ns - 1) *
   (2 * kovh^2 * 2998^2 / 5 / Ωm)^2 *
   t_nowiggle(kovh * h0, ωm, fb)^2 *
   2 *
   π^2 / kovh^3

# %% Compute `Ωgrav = Ωm * W / 2` from halos, using M200m and Mvir, and tinker08MF and tinker10MF
z = redshift
println("redshift = ", z)
Ωgrav_halo = ograv_halo(pk, z, Ωm)
println("Ωgrav(Δm=200, T08) = ", Ωgrav_halo)
Ωgrav_halo = ograv_halo(pk, z, Ωm, virial = true)
println("Ωgrav(Δvir, T08) = ", Ωgrav_halo)
Ωgrav_halo = ograv_halo(pk, z, Ωm, t10MF = true)
println("Ωgrav(Δm=200, T10) = ", Ωgrav_halo)
Ωgrav_halo = ograv_halo(pk, z, Ωm, virial = true, t10MF = true)
println("Ωgrav(Δvir, T10) = ", Ωgrav_halo)


text:
using Base.Cartesian

import Base.getindex

Base.linearindexing{T<:AbstractInterpolation}(::Type{T}) = Base.LinearSlow()

define_indices{IT}(::Type{IT}, N, pad) = Expr(:block, Expr[define_indices_d(iextract(IT, d), d, padextract(pad, d)) for d = 1:N]...)

coefficients{IT}(::Type{IT}, N) = Expr(:block, Expr[coefficients(iextract(IT, d), N, d) for d = 1:N]...)

function gradient_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim)
    exs = Expr[d==dim ? gradient_coefficients(iextract(IT, dim), d) :
                        coefficients(iextract(IT, d), N, d) for d = 1:N]
    Expr(:block, exs...)
end
function hessian_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim1, dim2)
    exs = if dim1 == dim2
        Expr[d==dim1==dim2 ? hessian_coefficients(iextract(IT, dim), d) :
                             coefficients(iextract(IT, d), N, d) for d in 1:N]
    else
        Expr[d==dim1 || d==dim2 ? gradient_coefficients(iextract(IT, dim), d) :
                                  coefficients(iextract(IT, d), N, d) for d in 1:N]
    end
    Expr(:block, exs...)
end

index_gen{IT}(::Type{IT}, N::Integer, offsets...) = index_gen(iextract(IT, min(length(offsets)+1, N)), IT, N, offsets...)

function getindex_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
    meta = Expr(:meta, :inline)
    quote
        $meta
        @nexprs $N d->(x_d = xs[d])

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

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

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

@generated function getindex{T,N}(itp::BSplineInterpolation{T,N}, xs::Number...)
    getindex_impl(itp)
end

function gradient_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
    meta = Expr(:meta, :inline)
    # For each component of the gradient, alternately calculate
    # coefficients and set component
    n = count_interp_dims(IT, N)
    exs = Array(Expr, 2n)
    cntr = 0
    for d = 1:N
        if count_interp_dims(iextract(IT, d), 1) > 0
            cntr += 1
            exs[2cntr-1] = gradient_coefficients(IT, N, d)
            exs[2cntr] = :(@inbounds g[$cntr] = $(index_gen(IT, N)))
        end
    end
    gradient_exprs = Expr(:block, exs...)
    quote
        $meta
        length(g) == $n || throw(ArgumentError(string("The length of the provided gradient vector (", length(g), ") did not match the number of interpolating dimensions (", n, ")")))
        @nexprs $N d->(x_d = xs[d])

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

        $gradient_exprs

        g
    end
end

function getindex_return_type{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}, argtypes)
    Tret = eltype(TCoefs)
    for a in argtypes
        Tret = Base.promote_op(@functorize(*), Tret, a) # the macro is used to support julia 0.4
    end
    Tret
end

@generated function gradient!{T,N}(g::AbstractVector, itp::BSplineInterpolation{T,N}, xs::Number...)
    length(xs) == N || error("Can only be called with $N indexes")
    gradient_impl(itp)
end

@generated function gradient!{T,N}(g::AbstractVector, itp::BSplineInterpolation{T,N}, index::CartesianIndex{N})
    args = [:(index[$d]) for d = 1:N]
    :(gradient!(g, itp, $(args...)))
end

# @eval uglyness required for disambiguation with method in Base
for R in [:Real, :Any]
    @eval @generated function gradient{T,N}(itp::AbstractInterpolation{T,N}, xs::$R...)
        n = count_interp_dims(itp, N)
        Tg = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
        xargs = [:(xs[$d]) for d in 1:length(xs)]
        :(gradient!(Array($Tg,$n), itp, $(xargs...)))
    end
end

gradient1{T}(itp::AbstractInterpolation{T,1}, x) = gradient(itp, x)[1]

function hessian_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
    meta = Expr(:meta, :inline)
    # For each component of the hessian, alternately calculate
    # coefficients and set component
    n = count_interp_dims(IT, N)
    exs = Expr[]
    cntr = 0
    for d1 in 1:N, d2 in 1:N
        if count_interp_dims(iextract(IT,d1), 1) > 0 && count_interp_dims(iextract(IT,d2),1) > 0
            cntr += 1
            push!(exs, hessian_coefficients(IT, N, d1, d2))
            push!(exs, :(@inbounds H[$cntr] = $(index_gen(IT, N))))
        end
    end
    hessian_exprs = Expr(:block, exs...)

    quote
        $meta
        size(H) == ($n,$n) || throw(ArgumentError(string("The size of the provided Hessian matrix wasn't a square matrix of size ", size(H))))
        @nexprs $N d->(x_d = xs[d])

        $(define_indices(IT, N, Pad))

        $hessian_exprs

        H
    end
end

@generated function hessian!{T,N}(H::AbstractMatrix, itp::BSplineInterpolation{T,N}, xs::Number...)
    length(xs) == N || throw(ArgumentError("Can only be called with $N indexes"))
    hessian_impl(itp)
end

@generated function hessian!{T,N}(H::AbstractMatrix, itp::BSplineInterpolation{T,N}, index::CartesianIndex{N})
    args = [:(index[$d]) for d in 1:N]
    :(hessian!(H, itp, $(args...)))
end

@generated function hessian{T,N}(itp::AbstractInterpolation{T,N}, xs...)
    n = count_interp_dims(itp,N)
    TH = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
    xargs = [:(xs[$d]) for d in 1:length(xs)]
    :(hessian!(Array($TH,$n,$n), itp, $(xargs...)))
end

hessian1{T}(itp::AbstractInterpolation{T,1}, x) = hessian(itp, x)[1,1]

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


text:
using Revise
using AutomotivePOMDPs
using SARSOP
using QMDP
using POMDPs
using AutomotiveDrivingModels
using POMDPPolicies
using BeliefUpdaters
includet("scripts/old_policy_plot.jl")

params = CrosswalkParams(obstacles_visible = true)
env = CrosswalkEnv(params)
pomdp = SingleOCPOMDP(env = env,
                      no_ped_prob = 0.0,
                      p_birth = 0.3,
                      collision_cost=-1.5)

@show length(states(pomdp))
@show length(actions(pomdp))
@show length(observations(pomdp))

sarsop_solver = SARSOPSolver(precision = 1e-3,
                      randomization=true,
                      timeout = 8*3600.0,
                      pomdp_filename = "crosswalk.pomdpx",
                      policy_filename = "crosswalk_long.policy",
                      policy_interval = 3*3600.0)

# qmdp_solver = QMDPSolver(belres = 1e-4, verbose=true)

# qmdp_policy = solve(qmdp_solver, pomdp);

# pomdp.collision_cost = -1.5

sarsop_policy = solve(sarsop_solver, pomdp);

sarsop_policy = load_policy(pomdp, "crosswalk_long.policy");

# b0 = initialstate_distribution(pomdp)

# a = action(policy, b0)
# actionvalues(policy, b0)

using POMDPGifs, POMDPSimulators, Random

rng = MersenneTwister(1);
hr = HistoryRecorder(max_steps=40, rng=rng);
s0 = initialstate(pomdp, rng);
up = DiscreteUpdater(pomdp);
b0 = initialstate_distribution(pomdp);
hist = simulate(hr, pomdp, sarsop_policy, up, b0, s0);
makegif(pomdp, hist, filename=joinpath("./", "crosswalk.gif"), spec="s,a", fps=2);

# includet("legacy/policy_plot.jl")

function uncertainty_plots(pomdp, policy)
    plots = []
    for sig in [1e-3, 0.5, 1.0, 2.0]
        p = policy_plot(pomdp, policy, sig=sig, n_bins=5, n_pts=150)
        push!(plots, p)
    end
    return plots
end

qp = QMDPPolicy(qmdp_policy);
sp = SARSOPPolicy(sarsop_policy);
qmdp_plots = uncertainty_plots(pomdp, qp);
sarsop_plots = uncertainty_plots(pomdp, sp);
g = GroupPlot(4, 4, groupStyle="horizontal sep = 1.75cm, vertical sep = 1.5cm");
push!.(Ref(g), qmdp_plots);
push!.(Ref(g), sarsop_plots);
g;
save("qmdp_sarsop_policies.tex", g)
g

save("policy_plots2.pdf", g)

qmdp_plots[end]

sarsop_plots[1]

using Profile
using ProfileView
Profile.clear()

p = policy_plot(pomdp, sp, sig=1.0, n_bins=5, n_pts=100, v_ego=5.0)

ProfileView.view()

policy = sarsop_policy

v_ego = 5.0
v_ped = 1.0
X, Y = get_XY_grid(pomdp);
X = X[1:21]
grid = RectangleGrid(X, Y);
vals = get_grid_data(policy, grid, pomdp, v_ego, v_ped);
acts = actions(pomdp)
action_map2 = [acts[ci[1]].acc for ci in  vec(argmax(vals, dims=1))]
action_map2 = reshape(action_map2, (21, 11))
acts = action_map2
acts = acts'
acts = acts[end:-1:1,1:end]

acts2 = compute_acts3(pomdp, policy)

cmap = ColorMaps.RGBArrayMap(colormap("RdBu"))
PGFPlots.Image(acts, (X[1], X[end]), (Y[1], Y[end]), colormap=cmap)

value_plots = [];
for i=1:length(actions(pomdp))
    xyval = reshape(vals[i, :], (21, 11))
    xyval = xyval'
    xyval = xyval[end:-1:1,1:end]
    p = PGFPlots.Image(xyval, (X[1], X[end]), (Y[1], Y[end]), colormap=cmap)
    push!(value_plots, p)
end;
g = GroupPlot(4, 1);
push!.(Ref(g), value_plots);
g

xyval = reshape(maximum(vals, dims=1), (21, 11))
xyval = xyval'
xyval = xyval[end:-1:1,1:end]
p = PGFPlots.Image(xyval, (X[1], X[end]), (Y[1], Y[end]), colormap=cmap)

length(X)
length(Y)

v_ego = 5.0
v_ped = 1

PGFPlots.save("sarsop_cw.pdf", p)

# using POMDPTesting

# trans_prob_consistency_check(pomdp)

# obs_prob_consistency_check(pomdp)

# ego = AutomotivePOMDPs.xv_to_state(pomdp, 5.0, 0.0)
# ped = AutomotivePOMDPs.yv_to_state(pomdp, 1.0, 0.0)
# crash = false

# s = SingleOCState(crash, ego, ped)

# od = observation(pomdp, SingleOCAction(-4.0), s)

# ospace = observations(pomdp);
# p = 0.0
# for o in ospace
#     global p
#     p += pdf(od, o)
#     if !(pdf(od, o) ≈ 0.0)
#         println(o)
#     end
# end

# for o in od.it
#     println(o)
# end

# op = SingleOCObs(false, ego,AutomotivePOMDPs.yv_to_state(pomdp, 2.0, 0.0))

# pdf(od, op)

# findfirst(isequal(op), ospace)

# for (i, o) in enumerate(ospace)
#     if o.ego == ego && o.ped.posG.y ≈ 2.0
#         println(o)
#         println(i)
#     end
# end

# pdf(od, ospace[22])

# ospace[22].ped
# op.ped

# op.ped ≈ ospace[22].ped

# sum(o.p)

# o.it

# t = transition(pomdp, s, SingleOCAction(2.0))

# t.p

# sum(t.p)

# for sp in t.it
#     println(sp)
# end

# sspace = states(pomdp);

# p = 0.0
# for sp in sspace
#     global p
#     p += pdf(o, sp)
#     if pdf(o, sp) != 0.0
#         println(sp)
#     end
# end

# o.it

# pedp = AutomotivePOMDPs.yv_to_state(pomdp, 2.0, 0.0)
# sp = SingleOCState(false, ego, pedp)

# findfirst(isequal(sp), sspace)

# pdf(o, sp)

# pdf(o, sp)

# for (i, s) in enumerate(sspace)
#     if s.ped.posG.y ≈ 2.0 && s.ego.posG.x ≈ 5.0
#         println(s)
#         println(i)
#     end
# end

# a = findfirst(isequal(sp), sspace)

# sp = t.it[end]

# is_crash(pomdp, sp.ego, sp.ped)

# collision_checker(sp.ego, sp.ped, pomdp.ego_type, pomdp.ped_type)

# using AutoViz
# AutoViz.set_render_mode(:fancy)

# ego = Vehicle(sp.ego, pomdp.ego_type, EGO_ID)
# ped = Vehicle(sp.ped, pomdp.ped_type, PED_ID)

# render(Scene([ego, ped]), pomdp.env.roadway, cam=CarFollowCamera(1, 20.0))

# sum(t.p) ≈ 1.0

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# ------------------------------------------------------------------- #
# Useful helper functions
# ------------------------------------------------------------------- #
"Use the FOC wrt c1, c2, a1, a2 to comptue the exchange ratio p2/p1"
function get_exchangerate(m::AbstractModel, c1, c2, a1, a2)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    # Get numerator and denominator
    p2 = c1.^(1-σ1) .* (1-ω1) .* a1.^(σ1-1)
    p1 = c2.^(1-σ2) .* ω2 .* a2.^(σ2-1)

    # fxr = p2/p1 = num/den
    return p2/p1
end

"Use the FOC wrt c1, c2, b1, b2 to comptue the exchange ratio p2/p1"
function get_exchangerateb(m::AbstractModel, c1, c2, b1, b2)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    # Get numerator and denominator
    p2 = c1.^(1-σ1) .* ω1 .* b1.^(σ1-1)
    p1 = c2.^(1-σ2) .* (1-ω2) .* b2.^(σ2-1)

    # fxr = p2/p1 = num/den
    return p2/p1
end

"""
Given the model and a matrix `Nshocks × Ngrid` matrix of optimal policy choices,
construct an interpolant for lλ☆p as a function of (lλ☆, lzh, v, ϵ1, ϵ2, ϵ3)

(NOTE if `m` is a `ModelCV` then the variables `v` and `ϵ3` do not apply.)
"""
function get_lλ☆p_itp(m::AbstractModel, lλ☆p::Matrix{Float64})
    # Get the unique shocks for each shock
    eps_grids = Vector[unique(m.ϵ[:, i]) for i in 1:n_exog(m.exog)]

    # Build a gridded interpolant (so we can have uneven grids)
    linspace_grids = map(collect, m.m_ls)
    knots = (linspace_grids..., eps_grids...)
    Nknots = map(length, knots)
    lλ☆p_rs = reshape(lλ☆p', Nknots...)

    # Create interpolant
    lλ☆p_itp = interpolate(knots, lλ☆p_rs, Gridded(Linear()))

    return lλ☆p_itp
end

# ------------------------------------------------------------------- #
# Simulation functions
# ------------------------------------------------------------------- #
"""
Simulate the model given an interpolant for lλ☆p, a sequence of shocks for
`ϵ` and an initial level of lλ☆p.

Note that ϵ is assumed to be `Nshocks × T`, where `Nshocks` is the number of
shocks in the model (2 for CV, 3 for SV) and `T` is the desired simulation
length.

The function returns simulated paths of a1, b1, c1, a2, b2, c2, lλ☆, exogvals,
c2_c1_ratio, fxr,

where

- `exogvals` = `(lzh, lg)`` in the CV model and `(lzh, v, lg)` in the SV model
- `fxr = p2/p1`
"""
function simulate(m::AbstractModel, lλ☆p_itp::AbstractInterpolation,
                  ϵ::Matrix{Float64}; lλ0=0.0)

    # Simulate the exogenous parameters
    nshocks, capT = size(ϵ)
    initvals = init(m.exog)
    exogvals = simulate(m.exog, initvals, ϵ)

    # Allocate space for remaining info
    a1 = Array(Float64, capT-1)
    b1 = Array(Float64, capT-1)
    c1 = Array(Float64, capT-1)
    a2 = Array(Float64, capT-1)
    b2 = Array(Float64, capT-1)
    c2 = Array(Float64, capT-1)
    lλ☆ = Array(Float64, capT-1)
    fxr = Array(Float64, capT-1)

    # Initial condition and simulate
    lλ☆[1] = lλ0
    for t=1:capT-1
        # Get period t allocation
        lλ☆_t = lλ☆[t]

        # Get allocation values
        a1[t], b1[t], c1[t], a2[t], b2[t], c2[t] = get_allocation(m, lλ☆_t, exogvals[1, t])
        fxr[t] = get_exchangerate(m, c1[t], c2[t], a1[t], a2[t])

        if t<capT-1
            lλ☆[t+1] = lλ☆p_itp[lλ☆[t], exogvals[1:end-1, t]..., ϵ[:, t+1]...]
        end
    end
    c2_c1_ratio = c2 ./ c1

    return a1, b1, c1, a2, b2, c2, lλ☆, exogvals, c2_c1_ratio, fxr
end

# API method above
simulate(m::AbstractModel, lλ☆p_itp::AbstractInterpolation; capT=10_000, lλ0=0.0) =
    simulate(m, lλ☆p_itp, randn(n_exog(m), capT), lλ0=lλ0)

"""
Approximate the stochastic steady state of lλ☆ by simulating the stochastic
version of the model and computing the observed mean of lλ☆.
"""
function find_steadystate_lλ(m::AbstractModel, lλ☆p_itp::AbstractInterpolation;
                             capT=50_000)

    return mean(simulate(m, lλ☆p_itp; capT=capT)[7][end-10_000:end])
end

# ------------------------------------------------------------------- #
# Euler Errors
# ------------------------------------------------------------------- #
"""
Compute euler errors given the model, a level of lλ☆, the exogenous state,
and interpolands for the value functions and policy function. `exog_state` will
contain `(lzh, v)`
"""
function euler_err_t(m::ModelSV, lλ☆::Float64, exog_state::Array{Float64},
                     vf_itps::VFITP, lλ☆p_itp::AbstractInterpolation)

    # Unpack parameters
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ = length(m.Π)

    # Evaluate time t allocations
    lzh, v = exog_state[1], exog_state[2]
    a1, b1, c1, a2, b2, c2 = get_allocation(m, lλ☆, lzh)

    # Get all possible lλ☆p values for tomorrow
    lλ☆p = Float64[lλ☆p_itp[lλ☆, exog_state..., m.ϵ[i, :]...] for i in 1:Nϵ]

    # Get t+1 possible exogenous values
    exog_p = exog_step(m.exog, exog_state, m.ϵ)
    lgp = exog_p[end]
    gp = exp(lgp)

    # Get t+1 value functions
    ljp, lup = map(n->Array(Float64, n), fill(Nϵ, 2))
    for j=1:Nϵ
        curr_exogp_state = map(x->x[j], exog_p[1:end-1])
        ljp[j] = vf_itps[1][lλ☆p[j], curr_exogp_state...]
        lup[j] = vf_itps[2][lλ☆p[j], curr_exogp_state...]
    end

    # Evaluate μ values
    lμ1 = log(_CE(m.agent1, gp, exp(ljp), m.Π))
    lμ2 = log(_CE(m.agent2, gp, exp(lup), m.Π))

    # Evaluate current value functions
    lj, lu = eval_vfs(m, c1, c2, exp(lμ1), exp(lμ2))

    # Evaluate Euler Equation
    lk = lk_residual(m, lμ1, lμ2, 0.0)
    ee = exp(up_residual(m, ljp, lup, lλ☆p, lgp, lλ☆, lk)) - 1.0

    return ee
end

"Compute the Euler errors for the model along a random simulation path"
function euler_err_sim(m::AbstractModel, vf_itps::VFITP,
                       lλ☆p_itp::AbstractInterpolation; capT=10_000)

    # Allocate memory
    eul_err = Array(Float64, length(m.Π), capT)

    # Simulate economy
    simvals = simulate(m, lλ☆p_itp, randn(3, capT+1))
    lλ☆ = simvals[7]
    exogvals = simvals[8]

    for t=1:capT
        eul_err[:, t] = euler_err_t(m, lλ☆[t], exogvals[1:end-1, t], vf_itps, lλ☆p_itp)
    end

    return eul_err
end


text:
### test the mathematical part of the spline library ###

#+


text:
function stone_features(pos::GoPosition)
  features = zeros(UInt8, pos.env.N, pos.env.N, 16)

  num_deltas_avail = size(pos.board_deltas, 3)
  cumulative_deltas = cumsum(pos.board_deltas, 3)
  last_eight = repeat(pos.board, outer = [1, 1, 8])
  # apply deltas to compute previous board states
  last_eight[:, :, 2:num_deltas_avail + 1] .= last_eight[:, :, 2:num_deltas_avail + 1] .- cumulative_deltas
  # if no more deltas are available, just repeat oldest board.
  last_eight[:, :, num_deltas_avail + 2:end] .= last_eight[:, :, num_deltas_avail + 1]

  features[:, :, 1:2:end] .= last_eight .== pos.to_play
  features[:, :, 2:2:end] .= last_eight .== -pos.to_play
  return features + 0.
end

color_to_play_feature(pos::GoPosition) = pos.to_play * ones(UInt8, pos.env.N, pos.env.N, 1)

get_feats(player::MCTSPlayer) = get_feats(player.root.position)

get_feats(pos::GoPosition) = cat(3, stone_features(pos), color_to_play_feature(pos))


text:
"""
    Fermi.HartreeFock

Module for running Hartree--Fock computations in Fermi.
"""
module HartreeFock

using Fermi
using Fermi.Integrals: IntegralHelper
using Fermi.Geometry: Molecule
using Fermi.Output
using Lints
using LinearAlgebra

function print_header()
    @output repeat("=",80)*"\n"
    @output "|    {:<74}|\n" "Hartree Fock"
    @output "|    {:<74}|\n" "Hartree Fock"
    @output "|    {:<74}|\n" "Hartree Fock"
    @output "|        {:<70}|\n" "Module by M.M. Davis and G.J.R Aroeira"
    @output repeat("=",80)*"\n"
end

"""
    Fermi.HartreeFock.AbstractHFWavefunction

Abstract type common to all Hartree-Fock wave functions.

_struct tree:_

**AbstractHFWavefunction** <: AbstractReferenceWavefunction <: AbstractWavefunction
"""
abstract type AbstractHFWavefunction <: Fermi.AbstractReferenceWavefunction end

#Restricted HF
include("RHF/RHF.jl")

end #module


text:
function _can_be_restricted_to(trian_i,trian)
  if have_compatible_domains(trian_i,trian) ||
    have_compatible_domains(trian_i,get_background_triangulation(trian)) ||
    have_compatible_domains(trian_i,get_background_triangulation(get_background_triangulation(trian))) ||
    Geometry.is_included(trian,trian_i)
    true
  else
    false
  end
end

abstract type MultiFieldStyle end

struct ConsecutiveMultiFieldStyle <: MultiFieldStyle end

struct StridedMultiFieldStyle <: MultiFieldStyle end

"""
    struct MultiFieldFESpace{S<:MultiFieldStyle,B} <: FESpace
      spaces::Vector{<:SingleFieldFESpace}
      multi_field_style::S
      constraint_style::Val{B}
    end
"""
struct MultiFieldFESpace{MS<:MultiFieldStyle,CS<:ConstraintStyle,V} <: FESpace
  vector_type::Type{V}
  spaces::Vector{<:SingleFieldFESpace}
  multi_field_style::MS
  constraint_style::CS
  function MultiFieldFESpace(
    ::Type{V},
    spaces::Vector{<:SingleFieldFESpace},
    multi_field_style::MultiFieldStyle) where V
    @assert length(spaces) > 0

    MS = typeof(multi_field_style)
    if any( map(has_constraints,spaces) )
      constraint_style = Constrained()
    else
      constraint_style = UnConstrained()
    end
    CS = typeof(constraint_style)
    new{MS,CS,V}(V,spaces,multi_field_style,constraint_style)
  end
end

"""
    MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
"""
function MultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
  Ts = map(get_dof_value_type,spaces)
  T = typeof(*(map(zero,Ts)...))
  MultiFieldFESpace(Vector{T},spaces,ConsecutiveMultiFieldStyle())
end

function MultiFieldFESpace(::Type{V},spaces::Vector{<:SingleFieldFESpace}) where V
  MultiFieldFESpace(V,spaces,ConsecutiveMultiFieldStyle())
end

MultiFieldStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = S()
MultiFieldStyle(f::MultiFieldFESpace) = MultiFieldStyle(typeof(f))

# Implementation of FESpace

function FESpaces.get_triangulation(f::MultiFieldFESpace)
  s1 = first(f.spaces)
  trian = get_triangulation(s1)
  @check all(map(i->trian===get_triangulation(i),f.spaces))
  trian
end

function FESpaces.num_free_dofs(f::MultiFieldFESpace)
  n = 0
  for U in f.spaces
    n += num_free_dofs(U)
  end
  n
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace)
  get_free_dof_ids(f,MultiFieldStyle(f))
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace,::MultiFieldStyle)
  @abstractmethod
end

function FESpaces.get_free_dof_ids(f::MultiFieldFESpace,::ConsecutiveMultiFieldStyle)
  block_num_dofs = Int[]
  for U in f.spaces
    push!(block_num_dofs,num_free_dofs(U))
  end
  BlockArrays.blockedrange(block_num_dofs)
end

FESpaces.get_dof_value_type(f::MultiFieldFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V)

FESpaces.get_vector_type(f::MultiFieldFESpace) = f.vector_type

FESpaces.ConstraintStyle(::Type{MultiFieldFESpace{S,B,V}}) where {S,B,V} = B()

function FESpaces.get_fe_basis(f::MultiFieldFESpace)
  nfields = length(f.spaces)
  all_febases = SingleFieldFEBasis[]
  for field_i in 1:nfields
    dv_i = get_fe_basis(f.spaces[field_i])
    cell_basis = lazy_map(BlockMap(nfields,field_i),get_data(dv_i))
    trian = get_triangulation(dv_i)
    bs = BasisStyle(dv_i)
    @assert bs == TestBasis()
    ds = DomainStyle(dv_i)
    dv_i_b = SingleFieldFEBasis(cell_basis,trian,bs,ds)
    push!(all_febases,dv_i_b)
  end
  MultiFieldCellField(all_febases)
end

function FESpaces.get_trial_fe_basis(f::MultiFieldFESpace)
  nfields = length(f.spaces)
  all_febases = SingleFieldFEBasis[]
  for field_i in 1:nfields
    du_i = get_trial_fe_basis(f.spaces[field_i])
    cell_basis = lazy_map(BlockMap((1,nfields),field_i),get_data(du_i))
    trian = get_triangulation(du_i)
    bs = BasisStyle(du_i)
    @assert bs == TrialBasis()
    ds = DomainStyle(du_i)
    du_i_b = SingleFieldFEBasis(cell_basis,trian,bs,ds)
    push!(all_febases,du_i_b)
  end
  MultiFieldCellField(all_febases)
end

function FESpaces.FEFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    FEFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function FESpaces.EvaluationFunction(fe::MultiFieldFESpace, free_values)
  blocks = map(1:length(fe.spaces)) do i
    free_values_i = restrict_to_field(fe,free_values,i)
    EvaluationFunction(fe.spaces[i],free_values_i)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

function CellData.CellField(fe::MultiFieldFESpace,cell_values)
  single_fields = map(1:length(fe.spaces)) do i
    cell_values_field = lazy_map(a->a.array[i],cell_values)
    CellField(fe.spaces[i],cell_values_field)
  end
  MultiFieldCellField(single_fields)
end

"""
    restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
"""
function restrict_to_field(f::MultiFieldFESpace,free_values::AbstractVector,field::Integer)
  _restrict_to_field(f,MultiFieldStyle(f),free_values,field)
end

function  _restrict_to_field(f,::MultiFieldStyle,free_values,field)
  @notimplemented
end

function  _restrict_to_field(f,::ConsecutiveMultiFieldStyle,free_values,field)
  offsets = compute_field_offsets(f)
  U = f.spaces
  pini = offsets[field] + 1
  pend = offsets[field] + num_free_dofs(U[field])
  SubVector(free_values,pini,pend)
end

"""
    compute_field_offsets(f::MultiFieldFESpace)
"""
function compute_field_offsets(f::MultiFieldFESpace)
  @assert MultiFieldStyle(f) == ConsecutiveMultiFieldStyle()
  U = f.spaces
  n = length(U)
  offsets = zeros(Int,n)
  for i in 1:(n-1)
    Ui = U[i]
    offsets[i+1] = offsets[i] + num_free_dofs(Ui)
  end
  offsets
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrain flag for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_isconstrained(f,trian)
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace,trian::Triangulation)
  data = map(f.spaces) do space
    trian_i = get_triangulation(space)
    if _can_be_restricted_to(trian_i,trian)
      get_cell_isconstrained(space,trian)
    else
      Fill(false,num_cells(trian))
    end
  end
  lazy_map( (args...) -> +(args...)>0,  data...)
end

function FESpaces.get_cell_isconstrained(f::MultiFieldFESpace,trian::SkeletonTriangulation)
  plus = get_cell_isconstrained(f,trian.plus)
  minus = get_cell_isconstrained(f,trian.minus)
  lazy_map((l,r)-> l||r,plus,minus)
end

function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrain flag for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_is_dirichlet(f,trian)
end

function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace,trian::Triangulation)
  data = map(f.spaces) do space
    trian_i = get_triangulation(space)
    if _can_be_restricted_to(trian_i,trian)
      get_cell_is_dirichlet(space,trian)
    else
      Fill(false,num_cells(trian))
    end
  end
  lazy_map( (args...) -> +(args...)>0,  data...)
end

function FESpaces.get_cell_is_dirichlet(f::MultiFieldFESpace,trian::SkeletonTriangulation)
  plus = get_cell_is_dirichlet(f,trian.plus)
  minus = get_cell_is_dirichlet(f,trian.minus)
  lazy_map((l,r)-> l||r,plus,minus)
end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get
  the constrains for the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_constraints(f,trian)
end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace,trian::Triangulation)
  nfields = length(f.spaces)
  blockmask = [ _can_be_restricted_to(get_triangulation(Vi),trian) for Vi in f.spaces ]
  active_block_ids = findall(blockmask)
  active_block_data = Any[ get_cell_constraints(f.spaces[i],trian) for i in active_block_ids ]
  blockshape = (nfields,nfields)
  blockindices = [(i,i) for i in active_block_ids]
  lazy_map(BlockMap(blockshape,blockindices),active_block_data...)
end

function FESpaces.get_cell_constraints(f::MultiFieldFESpace,trian::SkeletonTriangulation)
  cell_values_plus = get_cell_constraints(f,trian.plus)
  cell_values_minus = get_cell_constraints(f,trian.minus)
  lazy_map(BlockMap((2,2),[(1,1),(2,2)]),cell_values_plus,cell_values_minus)
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get the DOF ids
  on top of the corresponding cells.
  """
  trians = map(get_triangulation,f.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_dof_ids(f,trian)
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::Triangulation)
  get_cell_dof_ids(f,trian,MultiFieldStyle(f))
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::SkeletonTriangulation)
  cell_values_plus = get_cell_dof_ids(f,trian.plus)
  cell_values_minus = get_cell_dof_ids(f,trian.minus)
  lazy_map(BlockMap(2,[1,2]),cell_values_plus,cell_values_minus)
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,::Triangulation,::MultiFieldStyle)
  @notimplemented
end

function FESpaces.get_cell_dof_ids(f::MultiFieldFESpace,trian::Triangulation,::ConsecutiveMultiFieldStyle)
  offsets = compute_field_offsets(f)
  nfields = length(f.spaces)
  blockmask = [ _can_be_restricted_to(get_triangulation(Vi),trian) for Vi in f.spaces ]
  active_block_ids = findall(blockmask)
  active_block_data = Any[]
  for i in active_block_ids
    cell_dofs_i = get_cell_dof_ids(f.spaces[i],trian)
    if i == 1
      push!(active_block_data,cell_dofs_i)
    else
      offset = Int32(offsets[i])
      o = Fill(offset,length(cell_dofs_i))
      cell_dofs_i_b = lazy_map(Broadcasting(_sum_if_first_positive),cell_dofs_i,o)
      push!(active_block_data,cell_dofs_i_b)
    end
  end
  lazy_map(BlockMap(nfields,active_block_ids),active_block_data...)
end

function _sum_if_first_positive(a,b)
  if a>0
    return a+b
  else
    return a
  end
end

# API for multi field case

"""
    num_fields(f::MultiFieldFESpace)
"""
function num_fields(f::MultiFieldFESpace)
  length(f.spaces)
end

Base.iterate(m::MultiFieldFESpace) = iterate(m.spaces)

Base.iterate(m::MultiFieldFESpace,state) = iterate(m.spaces,state)

Base.getindex(m::MultiFieldFESpace,field_id::Integer) = m.spaces[field_id]

# API for the ConsecutiveMultiFieldStyle
import Gridap.FESpaces: interpolate
import Gridap.FESpaces: interpolate_everywhere
import Gridap.FESpaces: interpolate_dirichlet

@deprecate(
  interpolate(fs::MultiFieldFESpace, object),
  interpolate(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_everywhere(fs::MultiFieldFESpace, object),
  interpolate_everywhere(object, fs::MultiFieldFESpace)
)

@deprecate(
  interpolate_dirichlet(fs::MultiFieldFESpace, object),
  interpolate_dirichlet(object, fs::MultiFieldFESpace)
)

"""
The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs
even in the case that the given cell field does not fulfill them)
"""
function FESpaces.interpolate(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    uhi = interpolate!(object, free_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
like interpolate, but also compute new degrees of freedom for the dirichlet component.
The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
"""
function FESpaces.interpolate_everywhere(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_everywhere!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end

"""
"""
function FESpaces.interpolate_dirichlet(objects, fe::MultiFieldFESpace)
  free_values = zero_free_values(fe)
  blocks = SingleFieldFEFunction[]
  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
    free_values_i = restrict_to_field(fe,free_values,field)
    dirichlet_values_i = zero_dirichlet_values(U)
    uhi = interpolate_dirichlet!(object, free_values_i,dirichlet_values_i,U)
    push!(blocks,uhi)
  end
  MultiFieldFEFunction(free_values,fe,blocks)
end


text:
@testset "helpers" begin
    r = renamer("a" => "A", "b" => "B", "c" => "C")
    @test r("a") == Sorted(1, "A")
    @test r("b") == Sorted(2, "B")
    @test r("c") == Sorted(3, "C")
    @test_throws KeyError r("d")
    @test string(r("a")) == "A"
    @test string(r("b")) == "B"
    @test string(r("c")) == "C"
    @test r("a") < r("b") < r("c")
    @test r("a") == r("a")
    @test r("a") != r("b")
    r̂ = renamer(["a" => "A", "b" => "B", "c" => "C"])
    @test r̂("a") == Sorted(1, "A")
    @test r̂("b") == Sorted(2, "B")
    @test r̂("c") == Sorted(3, "C")
    @test_throws KeyError r̂("d")
    @test string(r̂("a")) == "A"
    @test string(r̂("b")) == "B"
    @test string(r̂("c")) == "C"
    @test r̂("a") < r̂("b") < r̂("c")
    @test r̂("a") == r̂("a")
    @test r̂("a") != r̂("b")

    s = sorter("b", "c", "a")
    @test s("a") == Sorted(3, "a")
    @test s("b") == Sorted(1, "b")
    @test s("c") == Sorted(2, "c")
    @test_throws KeyError s("d")
    @test string(s("a")) == "a"
    @test string(s("b")) == "b"
    @test string(s("c")) == "c"
    @test s("b") < s("c") < s("a")
    @test s("a") == s("a")
    @test s("a") != s("b")
    ŝ = sorter(["b", "c", "a"])
    @test ŝ("a") == Sorted(3, "a")
    @test ŝ("b") == Sorted(1, "b")
    @test ŝ("c") == Sorted(2, "c")
    @test_throws KeyError ŝ("d")
    @test string(ŝ("a")) == "a"
    @test string(ŝ("b")) == "b"
    @test string(ŝ("c")) == "c"
    @test ŝ("b") < ŝ("c") < ŝ("a")
    @test ŝ("a") == ŝ("a")
    @test ŝ("a") != ŝ("b")

    a = Sorted(1, [1, 2])
    b = Sorted(1, [1, 2])
    c = Sorted(1, [1, 3])
    @test a == b
    @test hash(a) == hash(b)
    @test a != c
    @test hash(a) != hash(c)

    @test string(nonnumeric(1)) == "1"
    @test isless(nonnumeric(1), nonnumeric(2))
end


text:
using Documenter, FastTransforms

makedocs(modules=[FastTransforms],
			doctest = true,
			clean = true,
			format = :html,
			sitename = "FastTransforms.jl",
			authors = "Richard Mikael Slevinsky",
			pages = Any[
					"Home" => "index.md"
					]
			)

deploydocs(
    repo   = "github.com/JuliaApproximation/FastTransforms.jl.git",
    latest = "master",
    julia  = "0.6",
    osname = "linux",
    target = "build",
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make   = nothing
    )


text:
# collatz conjenture

#Jl Plots

function collatz_next(n)
    if n % 2 == 0
        return n ÷ 2
    else
        return 3n+1
    end
end

function collatz_path(n)
    steps = Vector{Int}()
    while n != 1
        push!(steps, n)
        n = collatz_next(n)
    end
    push!(steps, 1)
    return steps
end

using Plots

x = 1:10000
y = length.(collatz_path.(x))
Plots.scatter(x, y, label="Stopping time", legend=:topleft, size=(800,600))

x = 2 .^ (0:13)
y = length.(collatz_path.(x))
Plots.scatter!(x, y, color=:red, label="for powers of two")

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# Julia 1.1.0
s = "Let's build a better world, together."

println("\n\n---------------- firstindex:lastindex ----------------")
for i = firstindex(s):lastindex(s)
    try
        println("$(s[i : end]) <- s[$(i):end]")
    catch
        # ignore the index error
    end
end

println("\n\n---------------- for in ----------------")
for c in s
    println(c)
end

println("\n\n---------------- foreach ----------------")
foreach(println, s)

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire