text:
# DataTrack
#

export DataTrack, plotTracks, pdf

"""
    DataTrack(;params...)

Create a data track consisting of one or multiple layers.

# Layer parameters

  - `:x`
  - `:y`
  - `:color`
  - `:fill`
  - `:size`
  - `:linetype`
  - `:alpha`
  - `:shape`
    - `:asterisk`
    - `:circle`
    - `:cross`
    - `:diamond`
    - `:dot`
    - `:plus`
    - `:square`
    - `:triangle`
    - `:down-triangle`
    - `:right-triangle`
    - `:left-triangle`

  - `:geom`
    - `:point`
    - `:line`
    - `:col`

If a single layer is used, the parameters are specified directly.
If multiple layers are used, the parameters for each layer are
specified in a separate dictionary, and a vector of layer
dictionaries can be supplied as parameter `:layers`.

# Track parameters

The following parameters affect the track as a whole

- `:title` track title
- `:title_color`
- `:title_font`

- `:ymin` minimum y-value to plot, or `nothing`
- `:ymax` maximum y-value to plot, or `nothing`

- `:axis_line_width`
- `:axis_color`
- `:axis_font`
- `:axis_tick_length`

- `:base_line_y` y-coordinates for horizontal line, or `nothing`
- `:base_line_color`
- `:base_line_width`

- `:width` width of track
- `:height` height of track

- `:margin_top`
- `:margin_bottom`
- `:margin_left`
- `:margin_right`

"""
struct DataTrack <: AbstractRefFrame
    # AbstractRefFrame
    transform::AbstractTransform
    children::Vector{AbstractRefFrame}
    params::Dict

    function DataTrack(;params...)
        new(Translation(0., 0.), AbstractRefFrame[], Dict(params))
    end
end

"""
    pdf(ref::AbstractRefFrame, path, xmin, xmax; absolute = true)

"""
function pdf(ref::AbstractRefFrame, path, xmin, xmax; absolute = true)
    s = Cairo.CairoRecordingSurface()
    cc = Cairo.CairoContext(s)
    bbox = draw(cc, AbstractRefFrame[ref], ref, xmin, xmax)
    if absolute
        bbox = union(Rect(0,0,0,0), bbox)
    end
    pdfs = Cairo.CairoPDFSurface(path, bbox.w, bbox.h)
    pdfcc = CairoContext(pdfs)
    Cairo.set_source_surface(pdfcc, s, -bbox.x, -bbox.y)
    Cairo.rectangle(pdfcc, 0, 0, bbox.w, bbox.h)
    Cairo.fill(pdfcc)
    Cairo.finish(pdfs)
    Cairo.destroy(pdfcc)
    Cairo.destroy(cc)
end

"""
    plotTracks(tracks::Vector{T}, path, xmin, xmax;
                    track_spacing = 5mm, params...)

Plot `tracks` into PDF file at `path`.
"""
function plotTracks(tracks::Vector{T}, path, xmin, xmax;
                    track_spacing = 5mm, params...) where T <: AbstractRefFrame

    s = Cairo.CairoRecordingSurface()
    cc = Cairo.CairoContext(s)

    # keep track of union of bounding boxes
    ubbox = nothing

    offset = 0.0
    for track in tracks
        ref = RefFrame(transform = Translation(0.0, offset),
                       children = [track]; params...)
        bbox = draw(cc, [ref], ref, xmin, xmax)

        offset = offset + bbox.h + track_spacing
        if ubbox == nothing ubbox = bbox end
        ubbox = union(ubbox, bbox)
    end

    pdfs = Cairo.CairoPDFSurface(path, ubbox.w, ubbox.h)
    pdfcc = CairoContext(pdfs)
    Cairo.set_source_surface(pdfcc, s, ubbox.x, ubbox.y)
    Cairo.rectangle(pdfcc, 0, 0, ubbox.w, ubbox.h)
    Cairo.fill(pdfcc)
    Cairo.destroy(pdfcc)
    Cairo.destroy(cc)
    Cairo.finish(pdfs)

end

function draw(cc::Cairo.CairoContext, tc::Vector{T}, ref::AbstractRefFrame, xmin, xmax) where T <: AbstractRefFrame
    ubbox = nothing
    for c::AbstractRefFrame in ref.children
        bbox = draw(cc, vcat(tc, c), c, xmin, xmax)
        if ubbox == nothing ubbox = bbox end
        ubbox = union(ubbox, bbox)

        #Cairo.rectangle(cc, bbox.x, bbox.y, bbox.w, bbox.h)
        #Cairo.stroke(cc)
    end
    user2dev(ref, ubbox)
end

# draw y-axis at xpos using ticklabels at tickypos
# tc should include mapping to user coordinates
function yaxis(cc::Cairo.CairoContext, tc::Vector{T}, xpos, tickypos, ticklabels;
               offset = 2mm, label = "") where T <: AbstractRefFrame

    devx, devy = user2dev(tc, repmat([xpos], length(tickypos)), tickypos)
    devx .-= offset

    axis_line_width = getparam(tc, :axis_line_width)
    Cairo.set_source_rgba(cc, rgba(getparam(tc, :axis_color))...)
    Cairo.set_font_face(cc, getparam(tc, :axis_font))
    Cairo.set_line_width(cc, axis_line_width)

    Cairo.move_to(cc, devx[1], devy[1] + axis_line_width/2)
    Cairo.line_to(cc, devx[end], devy[end] - axis_line_width/2)

    for i in 1:length(tickypos)
        Cairo.move_to(cc, devx[i] - getparam(tc, :axis_tick_length), devy[i])
        Cairo.line_to(cc, devx[i], devy[i])
        Cairo.text(cc, devx[i] - 2*getparam(tc, :axis_tick_length), devy[i], string(ticklabels[i]),
                   halign="right", valign="center")
    end

    Cairo.stroke(cc)
end

# draw baseline at position at env var :base_line_y
function baseline(cc::Cairo.CairoContext, tc::Vector{T}) where T <: AbstractRefFrame
    if getparam(tc, :base_line_y) != nothing
        Cairo.set_source_rgba(cc, rgba(getparam(tc, :base_line_color))...)
        Cairo.set_line_width(cc, getparam(tc, :base_line_width))

        base_line_y = getparam(tc, :base_line_y)
        s = transform(tc[end])
        xmin = s.to.x ##tc[end].units.x
        xmax = s.to.x + s.to.w ##tc[end].units.x + tc[end].units.w
        bl_x, bl_y = user2dev(tc, [xmin, xmax], [base_line_y, base_line_y])

        Cairo.move_to(cc, bl_x[1], bl_y[1])
        Cairo.line_to(cc, bl_x[2], bl_y[2])

        Cairo.stroke(cc)
    end
end

# draw DataTrack
function draw(cc::Cairo.CairoContext, tc::Vector{T}, dt::DataTrack, xmin, xmax) where T <: AbstractRefFrame
    # tc already contains dt as last element

    layers = getparam(tc, :layers, params(dt))

    if ~ (typeof(layers) <: AbstractArray)
        layers = [layers]
    end

    # scale
    ymin = getparam(tc, :ymin, nothing)
    ymax = getparam(tc, :ymax, nothing)

    ymin == nothing && (ymin = minimum([minimum(getparam(vcat(tc, RefFrame(;l...)), :y)) for l in layers]))
    ymax == nothing && (ymax = maximum([maximum(getparam(vcat(tc, RefFrame(;l...)), :y)) for l in layers]))

    xspan = xmax - xmin
    yspan = ymax - ymin

    width = getparam(tc, :width)
    height = getparam(tc, :height)

    units = Rect(xmin, ymax, xspan, -yspan)
    bbox = Rect(getparam(tc, :margin_left),
                getparam(tc, :margin_top),
                width - getparam(tc, :margin_left) - getparam(tc, :margin_right),
                height - getparam(tc, :margin_top) - getparam(tc, :margin_bottom))

    datactx = RefFrame(transform = Scaling(bbox, units))

    # Cairo stuff

    # draw baseline
    baseline(cc, vcat(tc, datactx))

    # draw yaxis
    yticks = _ticks_default_linear([ymin, ymax])
    yaxis(cc, vcat(tc, datactx), xmin, yticks, _format_ticklabel.(yticks); offset = 2mm, label = getparam(tc, :ylab, ""))

    # title
    title = getparam(tc, :title, "")
    if title != ""
        Cairo.set_source_rgba(cc, rgba(getparam(tc, :title_color))...)
        Cairo.set_font_face(cc, getparam(tc, :title_font))
        xp, yp = user2dev(tc, 0.0, bbox.h/2)
        Cairo.text(cc, xp, yp, title, halign="center", valign="top", angle = 90.0)
        Cairo.stroke(cc)
    end

    # plot data using CairoRecordingSurface
    s = Cairo.CairoRecordingSurface()
    fcc = Cairo.CairoContext(s)

    for l in layers
        ctxs = vcat(tc, RefFrame(transform = Scaling(bbox, units); l...))

        # dispatch on geometry
        geometry(fcc, ctxs)

    end

    # copy relevant bbox from CairoRecordingSurface
    dbbox = user2dev(tc, bbox)
    Cairo.set_source_surface(cc, s, 0, 0)
    Cairo.rectangle(cc, dbbox.x, dbbox.y, dbbox.w, dbbox.h)
    Cairo.fill(cc)
    Cairo.destroy(fcc)
    Cairo.finish(s)

    user2dev(dt, Rect(0.0, 0.0, width, height))
end


text:
# All Geometries
"""Returns the geometry type, such as `Polygon` or `Point`."""
function geomtype(geom)
    throw(ErrorException(string("Unknown Geometry type. ",
    "Define GeoInterface.geomtype(::$(typeof(geom))) to return the desired type.")))
end

# All types
"""
    ncoord(geom) -> Integer

Return the number of coordinate dimensions (such as XYZ) for the geometry.
"""
ncoord(geom) = ncoord(geomtype(geom), geom)

"""Return `true` when the geometry is empty."""
isempty(geom) = ngeom(geom) == 0
"""Return `true` when the geometry doesn't cross itself."""
issimple(geom) = issimple(geomtype(geom), geom)

# Point
getcoord(geom, i::Integer) = getcoord(geomtype(geom), geom, i)

# Curve, LineString, MultiPoint
npoint(geom) = npoint(geomtype(geom), geom)
getpoint(geom, i::Integer) = getpoint(geomtype(geom), geom, i)

# Curve
startpoint(geom) = getpoint(geom, 1)
isclosed(geom) = isclosed(geomtype(geom), geom)
isring(geom) = isclosed(geom) && issimple(geom)
length(geom) = length(geomtype(geom), geom)
endpoint(geom) = getcoord(geom, length(geom))

# Surface
area(geom) = area(geomtype(geom), geom)
centroid(geom) = centroid(geomtype(geom), geom)
pointonsurface(geom) = pointonsurface(geomtype(geom), geom)
boundary(geom) = boundary(geomtype(geom), geom)

# Polygon/Triangle
nring(geom) = nring(geomtype(geom), geom) # TODO If this is more than one, it has interior rings (holes)
getring(geom) = getring(geomtype(geom), geom)

"""Returns the exterior ring of this Polygon as a `LineString`."""
getexterior(geom) = getexterior(geomtype(geom), geom)  # getring(geom, 1)
"""Returns the number of interior rings in this Polygon."""
nhole(geom)::Integer = nhole(geomtype(geom), geom)  # nrings - 1
"""Returns the Nth interior ring for this Polygon as a `LineString`."""
gethole(geom, i::Integer) = gethole(geomtype(geom), geom, i)  # getring + 1

# PolyHedralSurface
"""Returns the number of including polygons."""
npatch(geom)::Integer = npatch(geomtype(geom), geom)
"""Returns a polygon in this surface, the order is arbitrary."""
getpatch(geom, i::Integer) = getpatch(geomtype(geom), geom, i)
"""Returns the collection of polygons in this surface that bounds the given polygon “p” for any polygon “p” in the surface."""
boundingpolygons()

# GeometryCollection
ngeom(geom) = ngeom(geomtype(geom), geom)
getgeom(geom, i::Integer) = getgeom(geomtype(geom), geom, i)

# MultiLineString
nlinestring(geom) = nlinestring(geomtype(geom), geom)
getlinestring(geom, i::Integer) = getlinestring(geomtype(geom), geom, i)

# MultiPolygon
npolygon(geom) = npolygon(geomtype(geom), geom)
getpolygon(geom, i::Integer) = getpolygon(geomtype(geom), geom, i)

# Other methods
crs(geom) = missing  # or conforming to <:CoordinateReferenceSystemFormat in GeoFormatTypes

# DE-9IM, see https://en.wikipedia.org/wiki/DE-9IM
equals(a, b)::Bool = equals(geomtype(a), geomtype(b), a, b)
disjoint(a, b)::Bool = disjoint(geomtype(a), geomtype(b), a, b)
touches(a, b)::Bool = touches(geomtype(a), geomtype(b), a, b)
within(a, b)::Bool = within(geomtype(a), geomtype(b), a, b)
overlaps(a, b)::Bool = overlaps(geomtype(a), geomtype(b), a, b)
crosses(a, b)::Bool = crosses(geomtype(a), geomtype(b), a, b)
intersects(a, b)::Bool = intersects(geomtype(a), geomtype(b), a, b)
contains(a, b)::Bool = contains(geomtype(a), geomtype(b), a, b)
relate(a, b)::Bool = relate(geomtype(a), geomtype(b), a, b)

# Set theory
symdifference(a, b) = difference(geomtype(a), geomtype(b), a, b)
difference(a, b, inverse=false) = difference(geomtype(a), geomtype(b), a, b, inverse)
intersection(a, b, inverse=false) = intersection(geomtype(a), geomtype(b), a, b, inverse)
union(a, b) = union(geomtype(a), geomtype(b), a, b)

ext: jl


text:
using BenchmarkTools
using Suppressor

println("day01/day01_1.jl")
@btime include("day01/day01_1.jl")
println("day01/day01_2.jl")
@btime include("day01/day01_2.jl")

println("day02/day02_1.jl")
@btime include("day02/day02_1.jl")
println("day02/day02_2.jl")
@btime include("day02/day02_2.jl")

println("day03/day03_1.jl")
@btime include("day03/day03_1.jl")
println("day03/day03_2.jl")
@btime include("day03/day03_2.jl")

println("day04/day04_1.jl")
@btime include("day04/day04_1.jl")
println("day04/day04_2.jl")
@btime include("day04/day04_2.jl")

println("day05/day05_1.jl")
@btime include("day05/day05_1.jl")
println("day05/day05_2.jl")
@btime include("day05/day05_2.jl")

println("day06/day06_1.jl")
@btime include("day06/day06_1.jl")
println("day06/day06_2.jl")
@btime include("day06/day06_2.jl")

println("day07/day07_1.jl")
@btime include("day07/day07_1.jl")
println("day07/day07_2.jl")
@btime include("day07/day07_2.jl")

println("day08/day08_1.jl")
@btime include("day08/day08_1.jl")
println("day08/day08_2.jl")
@btime include("day08/day08_2.jl")

println("day09/day09_1.jl")
@btime include("day09/day09_1.jl")
println("day09/day09_2.jl")
@btime include("day09/day09_2.jl")

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot([sin, cos], 0, 25)


text:
module Spectra

export power_spectrum

using MPI
using ClimateMachine
using ..BalanceLaws
using ..ConfigTypes
using ..DGMethods
using ..Mesh.Interpolation
using ..Mesh.Grids
using ..MPIStateArrays
using ..VariableTemplates
using ..Writers

include("power_spectrum.jl")

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
const s = ArgParseSettings()

@add_arg_table s begin
    "make"
        help = "Make a system backup."
        action = :command
    "delete"
        help = "Delete a previous backup."
        action = :command
    "list"
        help = "List all backups."
        action = :command
    "config"
        help = "Edit configuration."
        action = :command
end

@add_arg_table! s["make"] begin
    "target"
        help = "Directory to back up."
        required = true
    "-d", "--dst"
        help = "Directory where the backup will be stored."
        default = nothing
    "-n", "--no-confirm"
        help = "Do not prompt for confirmation."
        action = :store_true
end

@add_arg_table s["delete"] begin
    "id"
        help = "Backup number to delete."
        required = true
        arg_type = Int
end

@add_arg_table s["list"] begin
    "-d", "--dst"
        help = "Directory containing the backups."
        default = nothing
end

@add_arg_table s["config"] begin
    "-d", "--dst"
        help = "Set the directory in which backups will be managed."
end

function await_confirmation(default; msg="Confirm?")
    prompt = string(msg * " (", crayon"green", default ? 'Y' : 'y', RESET, '/', crayon"red", !default ? 'N' : 'n', RESET, ')')
    println(prompt)
    try
        while true
            answer = strip(readline(), [' ', '\n'])
            isempty(answer) && return default
            if lowercase(answer) ∉ ["y", "n"]
                println("Answer not understood. ", prompt)
            elseif answer == "y"
                return true
            else
                @info "Operation aborted."
                return false
            end
        end
    catch e
        if e isa InterruptException
            @info "Operation aborted."
            exit(0)
        else
            rethrow(e)
        end
    end
end


text:
@testset "Nystrom Preconditioner" begin
    Random.seed!(0)
    n, r_true = 200, 20
    A = randn(n, r_true)
    A = A*A'

    r = round(Int, r_true * 1.2)
    k = round(Int, 0.9*r)
    μ = 1e-3
    Anys = RP.NystromSketch(A, k, r)

    D, V = eigen(A; sortby=x->-x)
    Vk = V[:,1:k]
    Dk = D[1:k]

    Ptrue = Symmetric(1/(D[k+1] + μ) * Vk*(Diagonal(Dk) + μ*I)*Vk' + (I - Vk*Vk'))
    Ptrue_inv = Symmetric((D[k+1] + μ) * Vk * Diagonal(1 ./ (Dk .+ μ)) * Vk' + (I - Vk*Vk'))

    @testset "Standard" begin
        P = RP.NystromPreconditioner(Anys, μ)
        @test eltype(P) == Float64
        @test opnorm(Ptrue - Matrix(P)) < μ
        for _ in 1:10
            x = randn(n)
            y = P \ x
            @test Ptrue \ x ≈ y
            ldiv!(P, x)
            @test y ≈ x
        end

    end

    @testset "Inverse" begin
        Pinv = RP.NystromPreconditionerInverse(Anys, μ)
        @test eltype(Pinv) == Float64
        @test opnorm(Ptrue_inv - Matrix(Pinv)) < μ
        for _ in 1:10
            x = randn(n)
            @test Ptrue_inv*x ≈ Pinv*x
        end
    end

end


text:
import Base: copy, show, print
using Base.Enums

@enum Player begin
    void = 0
    black = 1
    white = -1
end
other(player::Player)::Player = if (player == black) white elseif (player == white) black else void end
tosymbol(player::Player) = if (player == black) "x" elseif (player==white) "o" else "." end

# ----------------------
Point = NamedTuple{(:row, :col), Tuple{Int, Int}}

up(p::Point)::Point = (row = p.row - 1, col = p.col)
down(p::Point)::Point = (row = p.row + 1, col = p.col)
left(p::Point)::Point = (row = p.row, col = p.col - 1)
right(p::Point)::Point = (row = p.row, col = p.col + 1)
nbrs(p::Point)::Vector{Point} = [up(p), down(p), left(p), right(p)]

# ----------------------
struct Move
    point::Point
    is_play::Bool
    is_pass::Bool
    is_resign::Bool
end

play(point::Point)::Move = Move(point, true, false, false)
pass_turn()::Move = Move((row=0, col=0), false, true, false)
resign()::Move = Move((row=0, col=0), false, false, true)

# ----------------------
IdType = UInt16

struct GoString
    color::Player
    stones::Set{Point}
    liberties::Set{Point}
    id::IdType
end

remove_liberty!(s::GoString, p::Point)::Nothing = (delete!(s.liberties, p); nothing)
add_liberty!(s::GoString, p::Point)::Nothing = (push!(s.liberties, p); nothing)
num_liberties(s::GoString)::Int = length(s.liberties)

function merge!(receiver::GoString, sender::GoString)::Nothing
    @assert receiver.color === sender.color
    union!(receiver.stones, sender.stones)
    union!(receiver.liberties, sender.liberties)
    setdiff!(receiver.liberties, receiver.stones)
    nothing
end

const EMPTY_STRING = GoString(void, Set{Point}(), Set{Point}(), IdType(0))

isvoid(s::GoString) = (s.color == void)
copy(s::GoString) = isvoid(s) ? EMPTY_STRING : GoString(s.color, copy(s.stones), copy(s.liberties), s.id)
show(io::IO, s::GoString) = print(io, tosymbol(s.color))
print(io::IO, s::GoString) = show(io, s)

mutable struct GoStringBuffer
    dict::Dict{IdType, GoString}
    lastid::IdType
end


text:
# Basic Calculator in Julia:
# Includes : Addition, Subtraction, Multiplication, Division

function add(x,y)
    return x+y
end

function sub(x,y)
    return x-y
end

function mul(x,y)
    return x*y
end

function division(x,y)
    return x/y
end

println("1. Addition \n2. Subtraction \n3. Multiplication \n4. Division ")

choice = parse(Int64,readline())

println("Enter parameter 1: X= ")

x = parse(Int64,readline())

println("Enter parameter 2: Y= ")

y = parse(Int64,readline())

println("Result= ")
if choice==1
    println("addition")
    print("$x + $y = ",add(x,y))
elseif choice == 2
    println("subtraction")
    print("$x - $y = ",sub(x,y))
elseif choice == 3
    println("Multiplication")
    print("$x * $y = ",mul(x,y))
elseif choice == 4
    println("Division")
    print("$x / $y = ",division(x,y))
end

max_stars_repo_name: chartews/Simple-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects


text:
using Documenter, RAD
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/MurrellGroup/RAD.jl.git",
    julia  = "nightly",
    osname = "osx"
makedocs()


text:
abstract type MultiMethod end
struct WeightedSum <: MultiMethod end
struct EpsilonCons <: MultiMethod end

struct NBI <: MultiMethod
    inequality::Bool
end
NBI() = NBI(false)

abstract type ProblemLinearity end
struct LinearProblem <: ProblemLinearity end
struct NonLinearProblem <: ProblemLinearity end

"""
    multisolve(m::Model, data::MultiData, met::MultiMethod, lin::ProblemLinearity)

Solves the multiobjective problem of type `lin` with method `met`.
"""
function multisolve end

function islinear(m::MultiData)
    m.linear ? LinearProblem() : NonLinearProblem()
end

function solvehook(m::Model; method::MultiMethod = WeightedSum(),
                   kwargs...)
    data = get_multidata(m)
    linearity = islinear(data)
    status = multisolve(m, data, method, linearity)
    return status
end

"""
Extract the Utopia point from the MultiData struct.
"""
function getutopia(m::MultiData)
    Fmin = diag(m.Phi)
    # Convert values depending on whether
    # we maximize or minimize them.
    multiplier = (obj.sense == MOI.MAX_SENSE ? -1. : 1. for obj in m.objectives)
    return [mult * f for (mult, f) in zip(multiplier, Fmin)]
end

"""
Extract the Nadir point from the MultiData struct.
"""
function getnadir(m::MultiData)
    Fmax = (f for f in maximum(m.Phi, dims=2))
    # Convert values depending on whether
    # we maximize or minimize them.
    multiplier = (obj.sense == MOI.MAX_SENSE ? -1. : 1. for obj in m.objectives)
    return [mult * f for (mult, f) in zip(multiplier, Fmax)]
end


text:
module Foo
    macro foo(x)
        esc(:(bar() = $x))
    end
    macro foo_y(x)
        esc(:(bar(y) = y + $x))
    end
    macro multifoo(x)
        _multifoo_impl(x)
    end
    _multifoo_impl(x) = esc(quote
        multibar() = $x
        multibar(y) = y
        multibar(z::Int, args...) = z + 1
    end)
end


text:
# This file includes the Discrete Systems

import DifferentialEquations: FunctionMap, DiscreteProblem
import UUIDs: uuid4

"""
    @def_discrete_system ex

where `ex` is the expression to define to define a new AbstractDiscreteSystem component type. The usage is as follows:
```julia
@def_discrete_system mutable struct MyDiscreteSystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractDiscreteSystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    righthandside::RH = righthandside_function  # mandatory field
    readout::RO = readout_function              # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MyDiscreteSystem` has `N` parameters. `MyDiscreteSystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the state, input port and output port of `MyDiscreteSystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside(dx, x, u, t, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New discrete system must be a subtype of `AbstractDiscreteSystem` to function properly.

# Example
```jldoctest
julia> @def_discrete_system mutable struct MyDiscreteSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
       α::Float64 = 1.
       β::Float64 = 2.
       righthandside::RH = (dx, x, u, t, α=α) -> (dx[1] = α * x[1] + u[1](t))
       state::Vector{Float64} = [1.]
       readout::RO = (x, u, t) -> x
       input::IP = Inport(1)
       output::OP = Outport(1)
       end

julia> ds = MyDiscreteSystem();

julia> ds.input
1-element Inport{Inpin{Float64}}:
 Inpin(eltype:Float64, isbound:false)
```
"""
macro def_discrete_system(ex)
    #= NOTE:
        Generate the parameter names TR, HS, CB, ID, MA, MK, SA, SK, AL, IT  of types of trigger, handshake, callbacks,
        id, modelargs, modelkwargs, solverargs, solverkwargs, alg, integrator by using `gensym` to avoid duplicate type
        parameter names so that the users can parametrizde their types as

        @def_discrete_system mutable struct MyDiscreteSystem{RH, RO, ST, IP, TR} <: AbstractDiscreteSystem
            righthandside::RH
            readout::RO
            state::ST
            input::IP
            output::TR
        end

        Note that the parameter name of output is TR. But, since the parameter name TR of trigger of the compnent
        is generated by `gensym`, we get rid of duplicate parameter names.
    =#
    checksyntax(ex, :AbstractDiscreteSystem)
    # TR, HS, CB, ID, MA, MK, SA, SK, AL, IT = [gensym() for i in 1 : 10]
    fields = quote
        trigger::$(TRIGGER_TYPE_SYMBOL) = Inpin()
        handshake::$(HANDSHAKE_TYPE_SYMBOL) = Outpin{Bool}()
        callbacks::$(CALLBACKS_TYPE_SYMBOL) = nothing
        name::Symbol = Symbol()
        id::$(ID_TYPE_SYMBOL)= Causal.uuid4()
        t::Float64 = 0.
        modelargs::$(MODEL_ARGS_TYPE_SYMBOL) = ()
        modelkwargs::$(MODEL_KWARGS_TYPE_SYMBOL) = NamedTuple()
        solverargs::$(SOLVER_ARGS_TYPE_SYMBOL) = ()
        solverkwargs::$(SOLVER_KWARGS_TYPE_SYMBOL) = NamedTuple()
        alg::$(ALG_TYPE_SYMBOL) = Causal.FunctionMap()
        integrator::$(INTEGRATOR_TYPE_SYMBOL) = Causal.construct_integrator(Causal.DiscreteProblem, input, righthandside, state, t,
            modelargs, solverargs; alg=alg, modelkwargs=modelkwargs, solverkwargs=solverkwargs, numtaps=3)
    end, [TRIGGER_TYPE_SYMBOL, HANDSHAKE_TYPE_SYMBOL, CALLBACKS_TYPE_SYMBOL, ID_TYPE_SYMBOL, MODEL_ARGS_TYPE_SYMBOL, MODEL_KWARGS_TYPE_SYMBOL, SOLVER_ARGS_TYPE_SYMBOL, SOLVER_KWARGS_TYPE_SYMBOL, ALG_TYPE_SYMBOL, INTEGRATOR_TYPE_SYMBOL]
    _append_common_fields!(ex, fields...)
    deftype(ex)
end

##### Define Discrete system library

"""
    DiscreteSystem(; righthandside, readout, state, input, output)

Constructs a generic discrete system

# Example
```jldoctest
julia> sfuncdiscrete(dx,x,u,t) = (dx .= 0.5x);

julia> ofuncdiscrete(x, u, t) = x;

julia> DiscreteSystem(righthandside=sfuncdiscrete, readout=ofuncdiscrete, state=[1.], input=nothing, output=Outport())
DiscreteSystem(righthandside:sfuncdiscrete, readout:ofuncdiscrete, state:[1.0], t:0.0, input:nothing, output:Outport(numpins:1, eltype:Outpin{Float64}))

```
"""
@def_discrete_system mutable struct DiscreteSystem{RH, RO, ST, IP, OP} <: AbstractDiscreteSystem
    righthandside::RH
    readout::RO
    state::ST
    input::IP
    output::OP
end

@doc raw"""
    DiscreteLinearSystem(input, output, modelargs=(), solverargs=();
        A=fill(-1, 1, 1), B=fill(0, 1, 1), C=fill(1, 1, 1), D=fill(0, 1, 1), state=rand(size(A,1)), t=0.,
        alg=ODEAlg, modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `DiscreteLinearSystem` with `input` and `output`. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

The `DiscreteLinearSystem` is represented by the following state and output equations.
```math
\begin{array}{l}
    \dot{x} = A x + B u \\[0.25cm]
    y = C x + D u
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation.
"""
@def_discrete_system mutable struct DiscreteLinearSystem{IP, OP, RH, RO} <: AbstractDiscreteSystem
    A::Matrix{Float64} = fill(-1., 1, 1)
    B::Matrix{Float64} = fill(0., 1, 1)
    C::Matrix{Float64} = fill(1., 1, 1)
    D::Matrix{Float64} = fill(-1., 1, 1)
    input::IP = Inport(1)
    output::OP = nothing
    state::Vector{Float64} = rand(size(A, 1))
    righthandside::RH = input === nothing ? (dx, x, u, t) -> (dx .= A * x) :
        (dx, x, u, t) -> (dx .= A * x + B * map(ui -> ui(t), u.itp))
    readout::RO = input === nothing ? (x, u, t) -> (C * x) :
           ( (C === nothing || D === nothing) ? nothing : (x, u, t) -> (C * x + D * map(ui -> ui(t), u)) )
end

@doc raw"""
    Henon()

Constructs a `Henon` system evolving with the dynamics
```math
\begin{array}{l}
    \dot{x}_1 = 1 - \alpha (x_1^2) + x_2 \\[0.25cm]
    \dot{x}_2 = \beta x_1
\end{array}
```
"""
@def_discrete_system mutable struct HenonSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
    α::Float64 = 1.4
    β::Float64 = 0.3
    γ::Float64 = 1.
    righthandside::RH = function henonrhs(dx, x, u, t, α=α, β=β, γ=γ)
        dx[1] = 1 - α * x[1]^2 + x[2]
        dx[2] = β * x[1]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(2)
end

@doc raw"""
    LoziSystem()

Constructs a `Lozi` system evolving with the dynamics
```math
\begin{array}{l}
    \dot{x}_1 = 1 - \alpha |x_1| + x_2 \\[0.25cm]
    \dot{x}_2 = \beta x_1
\end{array}
```
"""
@def_discrete_system mutable struct LoziSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
    α::Float64 = 1.4
    β::Float64 = 0.3
    γ::Float64 = 1.
    righthandside::RH = function lozirhs(dx, x, u, t, α=α, β=β, γ=γ)
        dx[1] = 1 - α * abs(x[1]) + x[2]
        dx[2] = β * x[1]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(2)
end

@doc raw"""
    BogdanovSystem()

Constructs a Bogdanov system with equations
```math
\begin{array}{l}
    \dot{x}_1 = x_1 + \dot{x}_2 \\[0.25cm]
    \dot{x}_2 = x_2 + \epsilon + x_2 + k x_1 (x_1 - 1) + \mu  x_1 x_2
\end{array}
```
"""
@def_discrete_system mutable struct BogdanovSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
    ε::Float64 = 0.
    μ::Float64 = 0.
    k::Float64 = 1.2
    γ::Float64 = 1.
    righthandside::RH = function bogdanovrhs(dx, x, u, t, ε=ε, μ=μ, k=k, γ=γ)
        dx[2]= x[2] + ε * x[2] + k * x[1] * (x[1] - 1) + μ * x[1] * x[2]
        dx[1] = x[1] + dx[2]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(2)
end

@doc raw"""
    GingerbreadmanSystem()

Constructs a GingerbreadmanSystem with the dynamics
```math
\begin{array}{l}
    \dot{x}_1 = 1 - x_2 + |x_1|\\[0.25cm]
    \dot{x}_2 = x_1
\end{array}
```
"""
@def_discrete_system mutable struct GingerbreadmanSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
    γ::Float64 = 1.
    righthandside::RH = function gingerbreadmanrhs(dx, x, u, t, γ=γ)
        dx[1] = 1 - x[2] + abs(x[1])
        dx[2] = x[1]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(2)
end

@doc raw"""
    LogisticSystem()

Constructs a LogisticSystem with the dynamics
```math
\begin{array}{l}
    \dot{x} = r x (1 - x)
\end{array}
```
"""
@def_discrete_system mutable struct LogisticSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
    r::Float64 = 1.
    γ::Float64 = 1.
    righthandside::RH = function logisticrhs(dx, x, u, t, r = r, γ=γ)
        dx[1] = r * x[1] * (1 - x[1])
        dx[1] *= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(1)
    input::IP = nothing
    output::OP = Outport(1)
end

##### Pretty-printting

show(io::IO, ds::DiscreteSystem) = print(io,
    "DiscreteSystem(righthandside:$(ds.righthandside), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::DiscreteLinearSystem) = print(io,
    "DiscreteLinearystem(A:$(ds.A), B:$(ds.B), C:$(ds.C), D:$(ds.D), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::HenonSystem) = print(io,
    "HenonSystem(α:$(ds.α), β:$(ds.β), γ:$(ds.γ),state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::LoziSystem) = print(io,
    "LoziSystem(α:$(ds.α), β:$(ds.β), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::BogdanovSystem) = print(io,
    "BogdanovSystem(ε:$(ds.ε), μ:$(ds.μ), k:$(ds.k), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::GingerbreadmanSystem) = print(io,
    "GingerbreadmanSystem(γ:$(ds.γ), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::LogisticSystem) = print(io,
    "LogisticSystem(r:$(ds.r), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
function DiffEqBase.__solve(prob::DiffEqBase.AbstractRODEProblem,
                            alg::Union{AbstractRODEAlgorithm,AbstractSDEAlgorithm},
                            timeseries=[],ts=[],ks=nothing, # needed for variable rate
                            recompile::Type{Val{recompile_flag}}=Val{true};
                            kwargs...) where recompile_flag
  integrator = DiffEqBase.__init(prob,alg,timeseries,ts,recompile;kwargs...)
  solve!(integrator)
  integrator.sol
end

# Make it easy to grab the RODEProblem/SDEProblem/DiscreteProblem from the keyword arguments
concrete_prob(prob) = prob
concrete_prob(prob::JumpProblem) = prob.prob

function DiffEqBase.__init(
  _prob::Union{DiffEqBase.AbstractRODEProblem,JumpProblem},
  alg::Union{AbstractRODEAlgorithm,AbstractSDEAlgorithm},timeseries_init=typeof(_prob.u0)[],
  ts_init=eltype(concrete_prob(_prob).tspan)[],
  ks_init=nothing,
  recompile::Type{Val{recompile_flag}}=Val{true};
  saveat = (),
  tstops = (),
  d_discontinuities= (),
  save_idxs = nothing,
  save_everystep = isempty(saveat),
  save_noise = save_everystep && (typeof(concrete_prob(_prob).f) <: Tuple ?
               DiffEqBase.has_analytic(concrete_prob(_prob).f[1]) : DiffEqBase.has_analytic(concrete_prob(_prob).f)),
  save_on = true,
  save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : concrete_prob(_prob).tspan[1] in saveat,
  save_end = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : concrete_prob(_prob).tspan[2] in saveat,
  callback=nothing,
  dense = false, # save_everystep && isempty(saveat),
  calck = false, #(!isempty(setdiff(saveat,tstops)) || dense),
  dt = eltype(concrete_prob(_prob).tspan)(0),
  adaptive = isadaptive(_prob,alg),
  gamma= isadaptive(alg) ? 9//10 : 0,
  abstol=nothing,
  reltol=nothing,
  qmax=qmax_default(alg),qmin=qmin_default(alg),
  qoldinit= isadaptive(alg) ? 1//10^4 : 0, fullnormalize=true,
  failfactor = 2,
  beta2=beta2_default(alg),
  beta1=beta1_default(alg,beta2),
  delta=delta_default(alg),
  maxiters = adaptive ? 1000000 : typemax(Int),
  dtmax=eltype(concrete_prob(_prob).tspan)((concrete_prob(_prob).tspan[end]-concrete_prob(_prob).tspan[1])),
  dtmin = typeof(one(eltype(concrete_prob(_prob).tspan))) <: AbstractFloat ? eps(eltype(concrete_prob(_prob).tspan)) :
          typeof(one(eltype(concrete_prob(_prob).tspan))) <: Integer ? 0 :
          eltype(concrete_prob(_prob).tspan)(1//10^(10)),
  internalnorm = ODE_DEFAULT_NORM,
  isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
  unstable_check = ODE_DEFAULT_UNSTABLE_CHECK,
  verbose = true,force_dtmin = false,
  timeseries_errors = true, dense_errors=false,
  advance_to_tstop = false,stop_at_next_tstop=false,
  initialize_save = true,
  progress=false,progress_steps=1000,progress_name="SDE",
  progress_message = ODE_DEFAULT_PROG_MESSAGE,
  userdata=nothing,
  initialize_integrator=true,
  seed = UInt64(0), alias_u0=false, alias_jumps = Threads.threadid()==1,
  kwargs...) where recompile_flag

  prob = concrete_prob(_prob)

  _seed = if iszero(seed)
    if (!(typeof(prob) <: DiffEqBase.AbstractRODEProblem) || iszero(prob.seed))
      seed_multiplier()*rand(UInt64)
    else
      prob.seed
    end
  else
    seed
  end

  if typeof(_prob) <: JumpProblem
    if alias_jumps
      _prob = DiffEqJump.resetted_jump_problem(_prob, _seed)
    elseif _seed !== 0
      DiffEqJump.reset_jump_problem!(_prob, _seed)
    end
  end

  # Grab the deepcopied version for caching reasons.
  prob = concrete_prob(_prob)

  if typeof(prob.f)<:Tuple
    if any(mm != I for mm in prob.f.mass_matrix)
      error("This solver is not able to use mass matrices.")
    end
  elseif typeof(prob) <: DiffEqBase.AbstractRODEProblem && prob.f.mass_matrix != I && !alg_mass_matrix_compatible(alg)
    error("This solver is not able to use mass matrices.")
  end

  if !isempty(saveat) && dense
    @warn("Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.")
  end

  if typeof(prob) <: DiffEqBase.AbstractRODEProblem && typeof(prob.noise)<:NoiseProcess && prob.noise.bridge === nothing && adaptive
    error("Bridge function must be given for adaptivity. Either declare this function in noise process or set adaptive=false")
  end

  if !alg_compatible(_prob,alg)
    error("The algorithm is not compatible with the chosen noise type. Please see the documentation on the solver methods")
  end

  progress && @logmsg(-1,progress_name,_id=_id = :StochasticDiffEq,progress=0)

  tType = eltype(prob.tspan)
  noise = prob isa DiffEqBase.AbstractRODEProblem ? prob.noise : nothing
  tspan = prob.tspan
  tdir = sign(tspan[end]-tspan[1])

  t = tspan[1]

  if !adaptive && iszero(dt) && isempty(tstops)
      error("Fixed timestep methods require a choice of dt or choosing the tstops")
  end

  f = prob.f
  p = prob.p
  g = prob isa DiffEqBase.AbstractSDEProblem ? prob.g : nothing

  if typeof(prob.u0) <: Tuple
    u = ArrayPartition(prob.u0,Val{true})
  else
    if alias_u0
      u = prob.u0
    else
      u = recursivecopy(prob.u0)
    end
  end

  uType = typeof(u)
  uBottomEltype = recursive_bottom_eltype(u)
  uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

  ks = ()

  uEltypeNoUnits = recursive_unitless_eltype(u)
  tTypeNoUnits   = typeof(one(tType))

  if abstol === nothing
    if uBottomEltypeNoUnits === uBottomEltype && !(uBottomEltype <: Integer)
      abstol_internal = real(convert(uBottomEltype,oneunit(uBottomEltype)*1//10^2))
    elseif uBottomEltype <: Integer
      abstol_internal = real(oneunit(uBottomEltype)*1//10^2)
    else
      abstol_internal = 0
    end
  else
    abstol_internal = real.(abstol)
  end

  if reltol === nothing
    if uBottomEltypeNoUnits === uBottomEltype && !(uBottomEltype <: Integer)
      reltol_internal = real(convert(uBottomEltype,oneunit(uBottomEltype)*1//10^2))
    elseif uBottomEltype <: Integer
      reltol_internal = real(oneunit(uBottomEltype)*1//10^2)
    else
      reltol_internal = 0
    end
  else
    reltol_internal = real.(reltol)
  end

  dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert
  # dtmin is all abs => does not care about sign already.

  if isinplace(prob) && typeof(u) <: AbstractArray && eltype(u) <: Number && uBottomEltypeNoUnits == uBottomEltype # Could this be more efficient for other arrays?
    if !(typeof(u) <: ArrayPartition)
      rate_prototype = recursivecopy(u)
    else
      rate_prototype = similar(u, typeof.(oneunit.(recursive_bottom_eltype.(u.x))./oneunit(tType))...)
    end
  else
    if uBottomEltypeNoUnits == uBottomEltype
      rate_prototype = u
    else # has units!
      rate_prototype = u/oneunit(tType)
    end
  end
  rateType = typeof(rate_prototype) ## Can be different if united

  if is_diagonal_noise(prob)
    noise_rate_prototype = rate_prototype
  elseif prob isa DiffEqBase.AbstractRODEProblem
    if prob isa DiffEqBase.AbstractSDEProblem
      noise_rate_prototype = copy(prob.noise_rate_prototype)
    else
      noise_rate_prototype = copy(prob.rand_prototype)
    end
  else
    noise_rate_prototype = nothing
  end

  tstops_internal, saveat_internal, d_discontinuities_internal =
    tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan)

  ### Algorithm-specific defaults ###
  # if save_idxs === nothing
  #   ksEltype = Vector{rateType}
  # else
  #   ks_prototype = rate_prototype[save_idxs]
  #   ksEltype = Vector{typeof(ks_prototype)}
  # end

  # Have to convert incase passed in wrong.
  if save_idxs === nothing
    timeseries = convert(Vector{uType},timeseries_init)
  else
    u_initial = u[save_idxs]
    timeseries = convert(Vector{typeof(u_initial)},timeseries_init)
  end
  ts = convert(Vector{tType},ts_init)

  alg_choice = typeof(alg) <: StochasticDiffEqCompositeAlgorithm ? Int[] : ()

  if !adaptive && save_everystep && tspan[2]-tspan[1] != Inf
    iszero(dt) ? steps = length(tstops) : steps = ceil(Int,internalnorm((tspan[2]-tspan[1])/dt,tspan[1]))
    sizehint!(timeseries,steps+1)
    sizehint!(ts,steps+1)
  elseif save_everystep
    sizehint!(timeseries,50)
    sizehint!(ts,50)
  elseif !isempty(saveat_internal)
    sizehint!(timeseries,length(saveat_internal)+1)
    sizehint!(ts,length(saveat_internal)+1)
  else
    sizehint!(timeseries,2)
    sizehint!(ts,2)
  end

  if save_start
    saveiter = 1 # Starts at 1 so first save is at 2
    copyat_or_push!(ts,1,t)
    if save_idxs === nothing
      copyat_or_push!(timeseries,1,u)
    else
      copyat_or_push!(timeseries,1,u_initial,Val{false})
    end
    if typeof(alg) <: StochasticDiffEqCompositeAlgorithm
      copyat_or_push!(alg_choice,1,1)
    end
  else
    saveiter = 0
  end

  QT = tTypeNoUnits <: Integer ? typeof(qmin) : tTypeNoUnits

  uprev = recursivecopy(u)

  if !(uType <: AbstractArray)
    rand_prototype = zero(u/u) # Strip units and type info
    randType = typeof(rand_prototype)
  else
    randElType = uBottomEltypeNoUnits # Strip units and type info
    if is_diagonal_noise(prob)
      if typeof(u) <: SArray
        rand_prototype = zero(u) # TODO: Array{randElType} for units
      else
        rand_prototype = (u .- u)./sqrt(oneunit(t))
      end
    elseif prob isa DiffEqBase.AbstractSDEProblem
      rand_prototype = false .* noise_rate_prototype[1,:]
    elseif prob isa DiffEqBase.AbstractRODEProblem
      rand_prototype = copy(prob.rand_prototype)
    else
      rand_prototype = nothing
    end
    randType = typeof(rand_prototype) # Strip units and type info
  end

  if typeof(_prob) <: JumpProblem
    callbacks_internal = CallbackSet(callback,_prob.jump_callback)
  else
    callbacks_internal = CallbackSet(callback)
  end

  max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
  if max_len_cb isa DiffEqBase.VectorContinuousCallback
    callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltype,uBottomEltype)
  else
    callback_cache = nothing
  end

  if typeof(prob) <: DiffEqBase.AbstractRODEProblem && prob.noise === nothing
    rswm = isadaptive(alg) ? RSWM(adaptivealg=:RSwM3) : RSWM(adaptivealg=:RSwM1)
    if isinplace(prob)
      #if isadaptive(alg) || callback !== nothing
        if alg_needs_extra_process(alg)
          W = WienerProcess!(t,rand_prototype,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        else
          W = WienerProcess!(t,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        end
      #=
      else
        if alg_needs_extra_process(alg)
          W = SimpleWienerProcess!(t,rand_prototype,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        else
          W = SimpleWienerProcess!(t,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        end
      end
      =#
    else
      #if isadaptive(alg) || callback !== nothing
        if alg_needs_extra_process(alg)
          W = WienerProcess(t,rand_prototype,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        else
          W = WienerProcess(t,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        end
      #=
      else
        if alg_needs_extra_process(alg)
          W = SimpleWienerProcess(t,rand_prototype,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        else
          W = SimpleWienerProcess(t,rand_prototype,
                             save_everystep=save_noise,
                             rng = Xorshifts.Xoroshiro128Plus(_seed))
        end
      end
      =#
    end
  elseif typeof(prob) <: DiffEqBase.AbstractRODEProblem
    W = prob.noise
    if hasfield(typeof(W),:t) && W.reset
      if W.t[end] != t
        reinit!(W,t)
      end
      # Reseed
      if typeof(W) <: NoiseProcess && W.reseed
        Random.seed!(W.rng,_seed)
      end
    elseif hasfield(typeof(W),:t) && W.t[end] != t
      error("Starting time in the noise process is not the starting time of the simulation. The noise process should be re-initialized for repeated use")
    end
  else # Only a jump problem
    @assert typeof(_prob) <: JumpProblem
    W = nothing
  end

  if typeof(_prob) <: JumpProblem && _prob.regular_jump !== nothing

    if !isnothing(_prob.regular_jump.mark_dist) == nothing # https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/250
      error("Mark distributions are currently not supported in SimpleTauLeaping")
    end

    jump_prototype = zeros(_prob.regular_jump.numjumps)
    c = _prob.regular_jump.c

    if isinplace(_prob.regular_jump)
      rate_constants = zeros(_prob.regular_jump.numjumps)
      _prob.regular_jump.rate(rate_constants,u./u,prob.p,tspan[1])
      P = CompoundPoissonProcess!(_prob.regular_jump.rate,t,jump_prototype,
                                  computerates = !alg_control_rate(alg) || !adaptive,
                                  save_everystep=save_noise,
                                  rng = Xorshifts.Xoroshiro128Plus(_seed))
      alg_control_rate(alg) && adaptive && P.cache.rate(P.cache.currate,u,p,tspan[1])
    else
      rate_constants = _prob.regular_jump.rate(u./u,prob.p,tspan[1])
      P = CompoundPoissonProcess(_prob.regular_jump.rate,t,jump_prototype,
                                 save_everystep=save_noise,
                                 computerates = !alg_control_rate(alg) || !adaptive,
                                 rng = Xorshifts.Xoroshiro128Plus(_seed))
      alg_control_rate(alg) && adaptive && (P.cache.currate = P.cache.rate(u,p,tspan[1]))
    end

  else
    jump_prototype = nothing
    c = nothing
    P = nothing
    rate_constants = nothing
  end

  dW,dZ = isnothing(W) ? (nothing,nothing) : (W.dW,W.dZ)

  cache = alg_cache(alg,prob,u,dW,dZ,p,rate_prototype,noise_rate_prototype,jump_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,Val{isinplace(_prob)})

  if typeof(_prob) <: JumpProblem && typeof(prob) <: DiscreteProblem && typeof(prob) <: Integer
    id = DiffEqBase.ConstantInterpolation(ts,timeseries)
  else
    id = LinearInterpolationData(timeseries,ts)
  end

  opts = SDEOptions(maxiters,save_everystep,
                    adaptive,abstol_internal,
                    reltol_internal,QT(gamma),
                    QT(qmax),QT(qmin),
                    QT(failfactor),
                    tType(dtmax),tType(dtmin),internalnorm,save_idxs,
                    tstops_internal,saveat_internal,
                    d_discontinuities_internal,
                    tstops,saveat,d_discontinuities,
                    userdata,
                    progress,progress_steps,
                    progress_name,progress_message,
                    timeseries_errors,dense_errors,
                    QT(beta1),QT(beta2),
                    convert.(uBottomEltypeNoUnits,delta),
                    QT(qoldinit),
                    dense,save_on,save_start,save_end,save_noise,
                    callbacks_internal,isoutofdomain,unstable_check,
                    verbose,calck,force_dtmin,
                    advance_to_tstop,stop_at_next_tstop)

  destats = DiffEqBase.DEStats(0)
  if typeof(alg) <: Union{StochasticDiffEqCompositeAlgorithm,
                          StochasticDiffEqRODECompositeAlgorithm}
    sol = DiffEqBase.build_solution(prob,alg,ts,timeseries,W=W,
                                    destats = destats,
                                    calculate_error = false, alg_choice=alg_choice,
                                    interp = id, dense = dense, seed = _seed)
  else
    sol = DiffEqBase.build_solution(prob,alg,ts,timeseries,W=W,
                                    destats = destats,
                                    calculate_error = false,
                                    interp = id, dense = dense, seed = _seed)
  end

  if recompile_flag == true
    FType = typeof(f)
    GType = typeof(g)
    CType = typeof(c)
    SolType = typeof(sol)
    cacheType = typeof(cache)
  else
    FType = Function
    GType = Function
    CType = Function
    SolType = DiffEqBase.AbstractRODESolution
    cacheType = StochasticDiffEqCache
  end

  tprev = t
  dtcache = tType(dt)
  iter = 0
  u_modified = false
  eigen_est = inv(one(tType)) # rate/state = (state/time)/state = 1/t units
  EEst = tTypeNoUnits(1)
  just_hit_tstop = false
  isout = false
  accept_step = false
  force_stepfail = false
  last_stepfail = false
  event_last_time = 0
  vector_event_last_time = 1
  last_event_error = zero(uBottomEltypeNoUnits)
  dtchangeable = true
  q11 = tTypeNoUnits(1)
  success_iter = 0
  q = tTypeNoUnits(1)

  integrator =    SDEIntegrator{typeof(alg),isinplace(prob),uType,
                  uBottomEltype,tType,typeof(tdir),typeof(p),
                  typeof(eigen_est),QT,
                  uEltypeNoUnits,typeof(W),typeof(P),rateType,typeof(sol),typeof(cache),
                  FType,GType,CType,typeof(opts),typeof(noise),typeof(last_event_error),typeof(callback_cache),typeof(rate_constants)}(
                  f,g,c,noise,uprev,tprev,t,u,p,tType(dt),tType(dt),tType(dt),dtcache,tspan[2],tdir,
                  just_hit_tstop,isout,event_last_time,vector_event_last_time,last_event_error,accept_step,
                  last_stepfail,force_stepfail,
                  dtchangeable,u_modified,
                  saveiter,
                  alg,sol,
                  cache,callback_cache,tType(dt),W,P,rate_constants,
                  opts,iter,success_iter,eigen_est,EEst,q,
                  QT(qoldinit),q11,destats)

  if initialize_integrator
    initialize_callbacks!(integrator, initialize_save)
    initialize!(integrator,integrator.cache)
    save_start && typeof(alg) <: Union{StochasticDiffEqCompositeAlgorithm,
                                       StochasticDiffEqRODECompositeAlgorithm} && copyat_or_push!(alg_choice,1,integrator.cache.current)
  end

  handle_dt!(integrator)

  ## Modify the first dt for tstops
  modify_dt_for_tstops!(integrator)
  ### Needs to be done before first rand
  integrator.sqdt = integrator.tdir*sqrt(abs(integrator.dt))

  !isnothing(integrator.W) && (integrator.W.dt = integrator.dt)
  !isnothing(integrator.P) && (integrator.P.dt = integrator.dt)

  DiffEqNoiseProcess.setup_next_step!(integrator::SDEIntegrator)

  integrator
end

function DiffEqBase.solve!(integrator::SDEIntegrator)
  @inbounds while !isempty(integrator.opts.tstops)
    while integrator.tdir*integrator.t < top(integrator.opts.tstops)
      loopheader!(integrator)
      if check_error!(integrator) != :Success
        return integrator.sol
      end
      perform_step!(integrator,integrator.cache)
      loopfooter!(integrator)
      if isempty(integrator.opts.tstops)
        break
      end
    end
    handle_tstop!(integrator)
  end
  postamble!(integrator)

  f = typeof(integrator.sol.prob.f) <: Tuple ? integrator.sol.prob.f[1] : integrator.sol.prob.f

  if DiffEqBase.has_analytic(f)
    DiffEqBase.calculate_solution_errors!(integrator.sol;timeseries_errors=integrator.opts.timeseries_errors,dense_errors=integrator.opts.dense_errors)
  end
  if integrator.sol.retcode != :Default
    return integrator.sol
  end
  integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol,:Success)
end

# Helpers

function handle_dt!(integrator)
  if iszero(integrator.dt) && integrator.opts.adaptive
    auto_dt_reset!(integrator)
    if sign(integrator.dt)!=integrator.tdir && !iszero(integrator.dt) && !isnan(integrator.dt)
      error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
    end
    if isnan(integrator.dt)
      if integrator.opts.verbose
        @warn("Automatic dt set the starting dt as NaN, causing instability.")
      end
    end
  elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) && integrator.tdir < 0
    integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
  end
end

function tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan)
  t0, tf = tspan
  tType = eltype(tspan)
  tdir = sign(tf - t0)

  tdir_t0 = tdir * t0
  tdir_tf = tdir * tf

  # time stops
  tstops_internal = BinaryMinHeap{tType}()
  if isempty(d_discontinuities) && isempty(tstops) # TODO: Specialize more
    push!(tstops_internal, tdir_tf)
  else
    for t in tstops
      tdir_t = tdir * t
      tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
    end

    for t in d_discontinuities
      tdir_t = tdir * t
      tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
    end

    push!(tstops_internal, tdir_tf)
  end

  # saving time points
  saveat_internal = BinaryMinHeap{tType}()
  if typeof(saveat) <: Number
    directional_saveat = tdir * abs(saveat)
    for t in (t0 + directional_saveat):directional_saveat:tf
      push!(saveat_internal, tdir * t)
    end
  elseif !isempty(saveat)
    for t in saveat
      tdir_t = tdir * t
      tdir_t0 < tdir_t < tdir_tf && push!(saveat_internal, tdir_t)
    end
  end

  # discontinuities
  d_discontinuities_internal = BinaryMinHeap{tType}()
  sizehint!(d_discontinuities_internal.valtree, length(d_discontinuities))
  for t in d_discontinuities
    push!(d_discontinuities_internal, tdir * t)
  end

  tstops_internal, saveat_internal, d_discontinuities_internal
end

function initialize_callbacks!(integrator, initialize_save = true)
  t = integrator.t
  u = integrator.u
  callbacks = integrator.opts.callback
  integrator.u_modified = true

  u_modified = initialize!(callbacks,u,t,integrator)

  # if the user modifies u, we need to fix previous values before initializing
  # FSAL in order for the starting derivatives to be correct
  if u_modified

    if isinplace(integrator.sol.prob)
      recursivecopy!(integrator.uprev,integrator.u)
    else
      integrator.uprev = integrator.u
    end

    if initialize_save &&
      (any((c)->c.save_positions[2],callbacks.discrete_callbacks) ||
      any((c)->c.save_positions[2],callbacks.continuous_callbacks))
      savevalues!(integrator,true)
    end
  end

  # reset this as it is now handled so the integrators should proceed as normal
  integrator.u_modified = false
end


text:
# This file contains some auxiliary code that is needed to process the DCE-MRI data through Julia

# This function will do non-linear fitting for the constrained reference region model
function fitCNRRM{N}(Ct::Array{Float64,N}, Crr::Vector{Float64}, t::Vector{Float64}; estKepRR=0.0, doTime=false)
  # Options:
  # estKepRR - Estimate for kepRR to use with CNRRM. If estKepRR<=0, then it will be estimated by first using NRRM
  # doTime - Records how long it takes to do the fitting
  DCEMRI.@dprint "fitting DCE data"
  (nT,nV) = size(Ct)
  @assert nT == length(t)
  mask = ones(nV) .> 0
  idxs = find(mask)
  nidxs = length(idxs)
  didNRRM = false
  if (estKepRR<=0.0)
    # First do NRRM, then get kepRR estimate, then do CNRRM
    didNRRM = true
    DCEMRI.@dprint "attempting Unconstrained Reference Region Tofts-Kety model"

    paramsNRRM = zeros(4, nV)
    residNRRM = Inf*ones(nV)
    runtimeNRRM = 0.0
    (paramsNRRM, residNRRM, runtimeNRRM) = fitNRRM(Ct,Crr,t,doTime=true)
    maskNegative = paramsNRRM .<= 0
    maskNegative = sum(maskNegative,1)
    maskPositive = vec(maskNegative .== 0)
    # Apply the mask, and obtain the median from the voxel-wise kepRR estimates
    paramsTmp = paramsNRRM[:,maskPositive]
    estKepRR = median(vec(paramsTmp[4,:]))
  end

  DCEMRI.@dprint "attempting Constrained Reference Region Tofts-Kety model"
  if doTime
    tic()
  end
  paramsCNRRM = zeros(3,nV)
  residCNRRM = Inf*ones(nV)
  p0 = [0.01, 0.01]
  fr2(x,p) = CRRM(x, p, Crr) # Cp is reference region curve. NOT AIF.
  p, r, dof = DCEMRI.nlsfit(fr2, Ct, idxs, [t; [estKepRR]], p0)
  r = squeeze(sumabs2(r, 1), 1) / dof
  for k in idxs
    residCNRRM[k] = r[k]
    paramsCNRRM[1,k] = p[1,k] # Relative KTrans
    paramsCNRRM[2,k] = estKepRR*p[1,k]/p[2,k] # Relative ve
    paramsCNRRM[3,k] = p[2,k] # kep
	end
  if doTime
    runtimeCNRRM = toc()
  end
  if didNRRM
    if doTime
      return (paramsCNRRM, residCNRRM, estKepRR, paramsNRRM, residNRRM, runtimeNRRM, runtimeCNRRM)
    else
      return (paramsCNRRM, residCNRRM, estKepRR, paramsNRRM, residNRRM)
    end
  else
    if doTime
      return (paramsCNRRM, residCNRRM, estKepRR, runtimeCNRRM)
    else
      return (paramsCNRRM, residCNRRM, estKepRR)
    end
  end
end

# This function will do non-linear fitting for the unconstrained reference region model
function fitNRRM{N}(Ct::Array{Float64,N}, Crr::Vector{Float64}, t::Vector{Float64}; doTime=false)
  DCEMRI.@dprint "fitting DCE data"
  (nT,nV) = size(Ct)
  @assert nT == length(t)

  mask = ones(nV) .> 0
  idxs = find(mask)
  nidxs = length(idxs)

  resid = Inf*ones(nV)
  params = zeros(4, nV)

  DCEMRI.@dprint "attempting Unconstrained Reference Region Tofts-Kety model"
  if doTime
    tic()
    runTime=0.0
  end
  p0 = [0.01, 0.01, 0.01]
  fr1(x,p) = RRM(x, p, Crr)
  p, r, dof = DCEMRI.nlsfit(fr1, Ct, idxs, t, p0)
  r = squeeze(sumabs2(r, 1), 1) / dof
  for k in idxs
    resid[k] = r[k]
    params[1,k] = p[1,k] # Relative KTrans
    params[2,k] = p[2,k]*p[1,k]/p[3,k] # Relative ve
    params[3,k] = p[3,k] # kep or (KTrans/ve)
    params[4,k] = p[2,k] #kep,RR
  end
  runTime = toc()
  if doTime
    return (params, resid, runTime)
  else
    return (params, resid)
  end
end

function fitLRRM(Ct::Matrix{Float64}, Crr::Vector{Float64}, t::Vector{Float64}, getResid=0)
	# Linear reference region model (based on Tofts-Kety)
	# Reference: J Cardenas-Rodriguez et al. MRI 31 (2013) 497-507
	# Input:
	#		t = time points of DCE acquisition. Vector of length nT.
	#		Ct = cocentration-time data. Matrix of dimensions nT-by-N, where N is #voxels.
	#   Crr = concentration-time curve of reference region. Vector of length nT.
	#		getResid = whether to return residuals in results. (1=True, 0=False)
	# Output:
	#		p = Fit parameters. 3-by-N matrix, where N is number of voxels
	#			The 3 elements of p are: [relativeKTrans, relativeVe, kep]
	#   resid = Residuals for each voxel. Vector of length N.
  (nT, nV) = size(Ct)
  resid = zeros(nV)
  trapzCrr = zeros(nT)
  p=zeros(3,nV)
  trapzCrr=cumtrapz(t,Crr)
  for i in 1:nV
    M = zeros(nT, 3);
    curCt = Ct[:,i]
    M[:,1] = Crr
    M[:,2] = trapzCrr
    M[:,3] = -cumtrapz(t,curCt)
    p[:,i] = M\curCt
		resid[i] = norm(curCt-M*p[:,i])
		# p = [kTrans_TOI/kTrans_RR, kTrans_TOI/ve_RR, kTrans_TOI/ve_TOI]
		# desire: p = [kTrans_TOI/ktrans_RR, ve_TOI/ve_RR, kTrans_TOI/ve_TOI]
  end
  p[2,:] = p[2,:]./p[3,:]
  if getResid==0
    return p
  elseif getResid==1
    return p,resid
  end
end

function fitCLRRM(t::Vector{Float64}, Ct::Matrix{Float64}, Crr::Vector{Float64}; getResid=0, estKepRR=1.0)
	# Constrained Linear Reference Region Model
	# Reference: Ahmed & Levesque. (2016), MRM, http://doi.org/10.1002/mrm.26530
	# Input:
	#		t = time points of DCE acquisition. Vector of length nT.
	#		Ct = cocentration-time data. Matrix of dimensions nT-by-N, where N is #voxels.
	#   Crr = concentration-time curve of reference region. Vector of length nT.
	#		getResid = whether to return residuals in results. (1=True, 0=False)
  #   estKepRR = estimate for the kepRR which will be held fixed
	# Output:
	#		p = Fit parameters. 3-by-N matrix, where N is number of voxels
	#			The 3 elements of p are: [relativeKTrans, relativeVe, kep]
	#   resid = Residuals for each voxel. Vector of length N.
  (nT, nV) = size(Ct)
  resid = zeros(nV)
  trapzCrr = zeros(nT)
  p=zeros(2,nV)
  pOut = zeros(3,nV)
  trapzCrr=cumtrapz(t,Crr)
  for i in 1:nV
    M = zeros(nT, 2);
    curCt = Ct[:,i]
    M[:,1] = Crr + estKepRR* trapzCrr
    M[:,2] = -cumtrapz(t,curCt)
    p[1:2,i] = M\curCt
		resid[i] = norm(curCt-M*p[:,i])
  end
  pOut[3,:] = p[2,:]
  pOut[2,:] = estKepRR*p[1,:]
  pOut[1,:] = p[1,:]
  # right now, pOut = [kTrans_TOI/kTrans_RR, kTrans_TOI/ve_RR, kTrans_TOI/ve_TOI]
  # we desire: pOut = [kTrans_TOI/ktrans_RR, ve_TOI/ve_RR, kTrans_TOI/ve_TOI]
  pOut[2,:] = pOut[2,:]./pOut[3,:]
  if getResid==0
    return pOut
  elseif getResid==1
    return pOut,resid
  end
end

function trapz{T<:Real}(x::Vector{T}, y::Vector{T})
	# Trapezoidal integration
	# Reference: https://groups.google.com/d/msg/julia-users/CNYaDUYog8w/lveUAdjkjmIJ
	# Input:
	# 	x = x values
	# 	y = y values of f(x)
	# Output:
	#		r = integral of f(x) from x[1] to x[end] using trapezoidal rule
	if (length(y) != length(x))
	   error("Vectors must be of same length")
	end
	r = 0.0
	for i in 2:length(x)
	    r += (x[i] - x[i-1]) * (y[i] + y[i-1])
	end
	return r/2.0
end

function cumtrapz{T<:Real}(x::Vector{T}, y::Vector{T})
	# Cumulative Trapezoidal integration
	# Input:
	# 	x = x values
	# 	y = y values of f(x)
	# Output:
	#		r = vector, where r(k) is integral of f(x) from x[1] to x[k]
  if (length(y) != length(x))
	   error("Vectors must be of same length")
	end
  r = zeros(length(x))
  for k in 1:length(x)
    r[k]=trapz(x[1:k],y[1:k])
  end
  return r
end

# This is a simple way of downsampling by omitting data points
function downsample(x::Vector{Float64}, dx::Int)
	samplePoints = collect(1:dx:length(x))
	numPoints = length(samplePoints)
	newX = zeros(numPoints)
	for i in 1:numPoints
		newX[i] = x[samplePoints[i]]
	end
	return(newX)
end

function downsample(x::Array{Float64}, dx::Int)
	(nT,nV) = size(x)
	samplePoints = collect(1:dx:nT)
	numPoints = length(samplePoints)
	newX = zeros(numPoints,nV)
	for i in 1:numPoints
		newX[i,:] = x[samplePoints[i], :]
	end
	return(newX)
end

max_stars_repo_name: MPUmri/CRRM

max_issues_repo_name: MPUmri/CRRM

max_forks_repo_name: MPUmri/CRRM


text:
__precompile__(true)

module Stego

using BitStreamLib

function lsb(img::Array{UInt8, 3}, data)
	data_bits::BitStream = BitStream(data)

	length_bits::BitStream = BitStream(Int32(length(data_bits)))

	img .&= 0xfe

	img_ind::Int = 1

	for i in length_bits
		img[img_ind] |= i ? 0x01 : 0x00
		img_ind += 1
	end

	for i in data_bits
		img[img_ind] |= i ? 0x01 : 0x00
		img_ind += 1
	end
end

function lsb(img::Array{UInt8, 3})
	length_bits::BitStream = BitStream()
	img_ind::Int = 1
	for i in 1 : 32
		push!(length_bits, img[i] & 0x01 != 0x00)
		img_ind += 1
	end
	len::Int32 = getvalue(Int32, length_bits)
	data_bits::BitStream = BitStream()

	for i in 1 : len
		push!(data_bits, img[img_ind] & 0x01 != 0x00)
		img_ind += 1
	end
	return data(data_bits)
end

function simplesteg(img::Array{UInt8, 3}, data, bits::Int32)
	data_bits::BitStream = BitStream(data)

	length_bits::BitStream = BitStream(Int32(length(data_bits)))

	img[1:16] .&= 0b11111100

	img_ind::Int = 1
	sub_ind::Int32 = 1

	for i in length_bits
		img[img_ind] |= i ? 0x01 : 0x00
		img_ind += 1
	end
end

end


text:
struct User
    id::Int64
    login::String
end

max_stars_repo_name: tosone/julia-spider

max_issues_repo_name: tosone/julia-spider

max_forks_repo_name: tosone/julia-spider


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Breadth-first search / traversal

#################################################
#
#  Breadth-first visit
#
#################################################
"""
    BreadthFirst

## Conventions in Breadth First Search and Depth First Search
### VertexColorMap
- color == 0    => unseen
- color < 0     => examined but not closed
- color > 0     => examined and closed

### EdgeColorMap
- color == 0    => unseen
- color == 1    => examined
"""
mutable struct BreadthFirst <: AbstractGraphVisitAlgorithm end

function breadth_first_visit_impl!(
    g::AbstractGraph,                       # the graph
    queue::Vector,                       # an (initialized) queue that stores the active vertices
    vertexcolormap::AbstractVertexMap,      # an (initialized) color-map to indicate status of vertices (-1=unseen, otherwise distance from root)
    edgecolormap::AbstractEdgeMap,          # an (initialized) color-map to indicate status of edges
    visitor::AbstractGraphVisitor,          # the visitor
    dir::Symbol                             # direction [:in,:out]
    )

    fneig = dir == :out ? out_neighbors : in_neighbors
    while !isempty(queue)
        u = shift!(queue)
        open_vertex!(visitor, u)
        u_color = vertexcolormap[u]

        for v in fneig(g, u)
            v_color = get(vertexcolormap, v, 0)
            v_edge = Edge(u,v)
            e_color = get(edgecolormap, v_edge, 0)
            examine_neighbor!(visitor, u, v, u_color, v_color, e_color) || return
            edgecolormap[v_edge] = 1
            if v_color == 0
                vertexcolormap[v] = u_color - 1
                discover_vertex!(visitor, v) || return
                push!(queue, v)
            end
        end
        close_vertex!(visitor, u)
        vertexcolormap[u] *= -1
    end
end

function traverse_graph!(
    g::AbstractGraph,
    alg::BreadthFirst,
    source,
    visitor::AbstractGraphVisitor;
    vertexcolormap::AbstractVertexMap = Dict{eltype(g), Int}(),
    edgecolormap::AbstractEdgeMap = DummyEdgeMap(),
    queue = Vector{eltype(g)}(),
    dir = :out)

    for s in source
        vertexcolormap[s] = -1
        discover_vertex!(visitor, s) || return
        push!(queue, s)
    end

    breadth_first_visit_impl!(g, queue, vertexcolormap, edgecolormap
            , visitor, dir)
end

#################################################
#
#  Useful applications
#
#################################################

###########################################
# Constructing BFS trees                  #
###########################################

"""
    TreeBFSVisitorVector{T}

A type for representing a BFS traversal of the graph as a parents array.
"""
mutable struct TreeBFSVisitorVector{T<:Integer} <: AbstractGraphVisitor
    tree::Vector{T}
end

function TreeBFSVisitorVector(n::Integer)
    return TreeBFSVisitorVector(fill(zero(n), n))
end

"""
    tree(parents)

Convert a parents array into a directed graph.
"""
function tree(parents::AbstractVector{T}) where T<:Integer
    n = T(length(parents))
    t = DiGraph(n)
    for i in one(T):n
        parent = parents[i]
        if parent > zero(T)  && parent != i
            add_edge!(t, parent, i)
        end
    end
    return t
end

tree(parents::TreeBFSVisitorVector) = tree(parents.tree)

function examine_neighbor!(visitor::TreeBFSVisitorVector, u::Integer, v::Integer,
                            ucolor::Int, vcolor::Int, ecolor::Int)
    if u != v && vcolor == 0
        visitor.tree[v] = u
    end
    return true
end

function bfs_tree!(visitor::TreeBFSVisitorVector{T},
        g::AbstractGraph,
        s::Integer;
        vertexcolormap = Dict{T,Int}(),
        queue = Vector{T}()) where T<:Integer
    # this version of bfs_tree! allows one to reuse the memory necessary to compute the tree
    # the output is stored in the visitor.tree array whose entries are the vertex id of the
    # parent of the index. This function checks if the scratch space is too small for the graph.
    # and throws an error if it is too small.
    # the source is represented in the output by a fixed point v[root] == root.
    # this function is considered a performant version of bfs_tree for useful when the parent
    # array is more helpful than a DiGraph struct, or when performance is critical.
    nvg = nv(g)
    length(visitor.tree) >= nvg || error("visitor.tree too small for graph")
    visitor.tree[s] = s
    traverse_graph!(g, BreadthFirst(), s, visitor; vertexcolormap=vertexcolormap, queue=queue)
end

"""
    bfs_tree(g, s)

Provide a breadth-first traversal of the graph `g` starting with source vertex `s`,
and return a directed acyclic graph of vertices in the order they were discovered.

### Implementation Notes
This function is a high level wrapper around [`bfs_tree!`](@ref); use that function for more performance.
"""
function bfs_tree(g::AbstractGraph, s::Integer)
    nvg = nv(g)
    visitor = TreeBFSVisitorVector(nvg)
    bfs_tree!(visitor, g, s)
    return tree(visitor)
end

############################################
# Connected Components with BFS            #
############################################
"""
    ComponentVisitorVector{T}

A type of `AbstractGraphVisitor` that represents connected components with
BFS starting from a given seed.
"""
mutable struct ComponentVisitorVector{T<:Integer} <: AbstractGraphVisitor
    labels::Vector{T}
    seed::T
end

function examine_neighbor!(visitor::ComponentVisitorVector, u::Integer, v::Integer,
                            ucolor::Int, vcolor::Int, ecolor::Int)
    if u != v && vcolor == 0
        visitor.labels[v] = visitor.seed
    end
    return true
end

############################################
# Test graph for bipartiteness             #
############################################
mutable struct BipartiteVisitor <: AbstractGraphVisitor
    bipartitemap::Vector{UInt8}
    is_bipartite::Bool
end

BipartiteVisitor(n::Integer) = BipartiteVisitor(zeros(UInt8,n), true)

function examine_neighbor!(visitor::BipartiteVisitor, u::Integer, v::Integer,
        ucolor::Int, vcolor::Int, ecolor::Int)
    if vcolor == 0
        visitor.bipartitemap[v] = (visitor.bipartitemap[u] == 1) ? 2 : 1
    else
        if visitor.bipartitemap[v] == visitor.bipartitemap[u]
            visitor.is_bipartite = false
        end
    end
    return visitor.is_bipartite
end

"""
    is_bipartite(g[, v])

Return `true` if graph `g` is [bipartite](https://en.wikipedia.org/wiki/Bipartite_graph).
If a node `v` is specified, only the connected component to which it belongs is considered.
"""
function is_bipartite(g::AbstractGraph)
    T = eltype(g)
    if !is_directed(g)
        cc = filter(x->length(x)>2, connected_components(g))
    else
        cc = filter(x->length(x)>2, weakly_connected_components(g))
    end
    vmap = Dict{T,Int}()
    for c in cc
        _is_bipartite(g,c[1], vmap=vmap) || return false
    end
    return true
end

is_bipartite(g::AbstractGraph, v::Integer) = _is_bipartite(g, v)

_is_bipartite(g::AbstractGraph, v::Integer; vmap = Dict{eltype(g),Int}()) = _bipartite_visitor(g, v, vmap=vmap).is_bipartite

function _bipartite_visitor(g::AbstractGraph, s::Integer; vmap=Dict{eltype(g),Int}())
    nvg = nv(g)
    visitor = BipartiteVisitor(nvg)
    for v in keys(vmap) #have to reset vmap, otherway problems with digraphs
        vmap[v] = 0
    end
    traverse_graph!(g, BreadthFirst(), s, visitor, vertexcolormap=vmap)
    return visitor
end

@doc_str """
    bipartite_map(g)

For a bipartite graph `g`, return a vector `c` of size ``|V|`` containing
the assignment of each vertex to one of the two sets (``c_i == 1`` or c_i == 2``).
If `g` is not bipartite, return an empty vector.
"""
function bipartite_map(g::AbstractGraph)
    cc = connected_components(g)
    visitors = [_bipartite_visitor(g, x[1]) for x in cc]
    !all([v.is_bipartite for v in visitors]) && return zeros(Int, 0)
    m = zeros(Int, nv(g))
    for i in vertices(g)
        m[i] = any(v->v.bipartitemap[i] == 1, visitors) ? 2 : 1
    end
    m
end

###########################################
# Get the map of the (geodesic) distances from vertices to source #
###########################################

"""
    gdistances!(g, source, dists)

Fill `dists` with the geodesic distances of vertices in `g` from `source`.
`dists` should be a vector of length `nv(g)`. Return `dists`.
For vertices in disconnected components the default distance is -1.
"""
function gdistances!(g::AbstractGraph, source, dists)
    T = eltype(g)
    n = nv(g)
    fill!(dists, -1)
    queue = Vector{T}(n)
    for i in 1:length(source)
        queue[i] = source[i]
        dists[source[i]] = 0
    end
    head = 1
    tail = length(source)
    while head <= tail
        current = queue[head]
        distance = dists[current] + 1
        head += 1
        for j in out_neighbors(g, current)
            if dists[j] == -1
                dists[j] = distance
                tail += 1
                queue[tail] = j
            end
        end
    end
    return dists
end

"""
    gdistances(g, source)

Return a vector filled with the geodesic distances of vertices in  `g` from
`source`. If `source` is a collection of vertices each element should be unique.
For vertices in disconnected components the default distance is -1.
"""
gdistances(g::AbstractGraph, source) = gdistances!(g, source, Vector{Int}(nv(g)))

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Revise
using DynamicalSystems
using Basins
using Plots

# Basin bifurcation in quasiperiodically forced systems Ulrike Feudel, Annette Witt, Ying-Cheng Lai, and Celso Grebogi PRE 28, 1998
function chaotic_map(dz, z, p, n)
    xn = z[1]
    θ = z[2]
    a = p[1]
    ω = (sqrt(5.0) - 1.0) / 2.0
    r = p[2]
    f(x) = r * x * (1.0 - x)
    Mn(n) = reduce(∘, fill(f, n))
    M = Mn(3)
    dz[1] = M(xn) + a * cos(2 * π * θ)
    dz[2] = mod(θ + ω, 1.0)
    return
end

# dummy function to keep the initializator happy
function chaotic_map_J(J, z0, p, n)
    return
end
ds = DiscreteDynamicalSystem(
    chaotic_map,
    [1.0, 0.0],
    [0.0015, 3.833],
    chaotic_map_J,
)
integ = integrator(ds)

θ = range(0.0, 1.0, length = 250)
xg = range(0.0, 1.0, length = 250)
integ.p[2] = 3.833
integ.p[1] = 0.0015

@time bsn = Basins.basins_map2D(xg, θ, integ)
#@time bsn=Basins.basin_general_ds(xg, θ, integ; dt=1, Ncheck=10)

@show bd = box_counting_dim(xg, θ, bsn)
#@show α = 2 - bd

D,De,e,f = uncertainty_exponent(bsn, integ)
println("uncertainty exponent: α=",D, " ± ",De)

plot(θ, xg, bsn.basin, seriestype = :heatmap)

#
# u0=[0.5, 0.];
# v = trajectory(ds,1000,u0)
#
# plot!(v[:,2],v[:,1],seriestype=:scatter)
#
# att1=Array{Float64,1}[]
# att2=Array{Float64,1}[]
# att3=Array{Float64,1}[]
#
# for v in bsn.attractors
# if v[1] == 1.
#     push!(att1,[v[2],v[3]])
# elseif v[1] == 2.
#     push!(att2,[v[2],v[3]])
# #elseif v[1] == 3.
# #    push!(att3,[v[2],v[3]])
# end
#
# end
#
# # att1 = Dataset(att1);
# # att2 = Dataset(att2);
# # att3 = Dataset(att3);
#
#
#
#
# function escape_function(u,u1,u2,u3)
#
#     m1 = minimum([norm(v - u) for v in u1])
#     m2 = minimum([norm(v - u) for v in u2])
#     #m3 = minimum([norm(v - u) for v in u3])
#     m3=1
#     if m1 < 5e-3
#         return 1
#     elseif m2 < 5e-3
#         return 2
#     elseif m3 < 5e-3
#         return 3
#     end
#     return 0
# end
#
#
# basin_t = zeros(length(xg),length(θ))
# for (k,x) in enumerate(xg), (n,y) in enumerate(θ)
#     reinit!(integ,[x,y])
#
#     while escape_function(integ.u,att1,att2,0) == 0
#         step!(integ)
#         integ.t > 200 && break
#     end
#
#     @show basin_t[k,n]=escape_function(integ.u,att1,att2,att3);
# end
#
# plot(θ,xg, basin_t, seriestype=:heatmap)
#

#
# # Estimation using box counting
# ind  = findall(iseven.(bsn.basin) .== true)
# basin_test = deepcopy(bsn.basin)
# [basin_test[k] =basin_test[k]+1 for k in ind ]
# bd = box_counting_dim(xg,θ, basin_test)
# @show α = 2 - bd
# Value in the paper is α=0.07
# But the basins are not reproductible. At least not with this method
#plot(θ,xg,bsn.basin, seriestype=:heatmap)
#
#
# r=3.7
# f(x)=r*x*(1-x)
# Mn(n) = reduce(∘, fill(f, n))
# fnn = Mn(4)
# x=0:0.01:1
# plot(x,fnn.(x))
# plot!([0,1],[0,1])
#


text:
using DataFrames
df = DataFrame()
df[:pb]=rand(1:50, 50);
df[:t]=rand(1:50, 50);
df[:h]=[0.0 for i in 1:50];

using Gadfly
plot(df, x = :pb, y = :t, color = :h, Geom.rectbin)


text:
module Drawers

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export drawer

register_normal_element("q__drawer", context = @__MODULE__)

function drawer(args...; kwargs...)
  q__drawer(args...; kw(kwargs)...)
end

end


text:
using LightXML
using Test

tests = ["parse", "create", "cdata", "pi"]

for t in tests
    fpath = "$t.jl"
    println("running $fpath ...")
    include(fpath)
end


text:
mutable struct ItemUse
    min::Vector{Int64}
    max::Vector{Int64}
    ItemUse(min, max) = new(min, max)
    ItemUse() = new(Int64[], Int64[])
end


text:
# =========================================================================================
for T in [:PotentialType, :SingleLayer, :DoubleLayer]
    @eval @doc """
        abstract type PotentialType end
        struct SingleLayer <: PotentialType end
        struct DoubleLayer <: PotentialType end

    Enum-like representation of single and double layer potentials
    """ $T
end
abstract type PotentialType end
struct SingleLayer <: PotentialType end
struct DoubleLayer <: PotentialType end

# =========================================================================================
for T in [:LocalityType, :LocalES, :NonlocalES]
    @eval @doc """
        abstract type LocalityType end
        struct NonlocalES <: LocalityType end
        struct LocalES    <: LocalityType end

    Enum-like representation of locality assumption:
     * *Local electrostatics*:
       Complete independence of solvent molecules
     * *Nonlocal electrostatics*:
       Allow solvent molecule correlation effects (with area-of-effect radius λ)
    """ $T
end
abstract type LocalityType end
struct NonlocalES <: LocalityType end
struct LocalES    <: LocalityType end

# =========================================================================================
"""
    φmol{T}(
        ξ        ::Vector{T},
        charges  ::Vector{Charge{T}};
        # kwargs
        tolerance::T                 = T(1e-10)
    )

Computes and returns the molecular potential of the given system of point charges in a
structureless medium for the given observation point ξ:

```math
φ_{mol}(ξ) = \\frac{1}{4π ε_0 ε_Ω} \\sum_i \\frac{qᵢ}{|rᵢ-ξ|}
```

If ``|rᵢ-ξ|`` is smaller than the given `tolerance`, the value is replaced by `tolerance`
for the affected charge.

!!! note
    The return value is premultiplied by ``4π ⋅ ε₀ ⋅ ε_Ω``

## Return type
`T`

## Aliases
    φmol{T}(
        Ξ        ::Vector{Vector{T}},
        charges  ::Vector{Charge{T}};
        # kwargs
        tolerance::T                 = T(1e-10)
    )

Computes the molecular potentials for a list of observation points.

    φmol{T}(
        model    ::Model{T, Triangle{T}};
        # kwargs
        tolerance::T                     = T(1e-10)
    )

Computes the molecular potentials for the given surface model, using each triangle center
as observation point.
"""
function φmol(
        ξ        ::Vector{T},
        charges  ::Vector{Charge{T}};
        tolerance::T=T(1e-10)
    ) where T
    # devectorized version of
    # sum([q.val / euclidean(ξ, q.pos) for q in charges])
    ret = zero(T)
    for q in charges
        ret += q.val / max(euclidean(ξ, q.pos), tolerance)
    end
    ret
end

function φmol(
        Ξ        ::Vector{Vector{T}},
        charges  ::Vector{Charge{T}};
        tolerance::T=T(1e-10)
    ) where T
    [φmol(ξ, charges, tolerance=tolerance) for ξ in Ξ]
end

function φmol(
        model    ::Model{T, Triangle{T}};
        tolerance::T=T(1e-10)
    ) where T
    [φmol(ξ.center, model.charges, tolerance=tolerance) for ξ in model.elements]
end

# =========================================================================================
"""
    ∂ₙφmol{T}(ξ::Triangle{T}, charges::Vector{Charge{T}})

Computes and returns the normal derivative of the given system's molecular potential in a
structureless medium, using the given triangle's center as observation point and the
triangle's normal as reference normal.

```math
∂ₙφ_{mol}(ξ) = -\\frac{1}{4π ε_0 ε_Ω} \\sum_i \\frac{qᵢ}{|rᵢ-ξ|³} (rᵢ-ξ) ⋅ n
```

!!! note
    The return value is premultiplied by ``4π ⋅ ε₀ ⋅ ε_Ω``

## Return type
`T`

## Aliases
    ∂ₙφmol{T}(model::Model{T, Triangle{T}})

Computes the normal derivatives of the molecular potentials for the given surface model,
using each triangle center and normal as observation point.
"""
function ∂ₙφmol(
        ξ      ::Triangle{T},
        charges::Vector{Charge{T}}
    ) where T
    ret = zero(T)
    for q in charges
        ret -= q.val * ddot(ξ.center, q.pos, ξ.normal) / euclidean(ξ.center, q.pos)^3
    end
    ret
end

function ∂ₙφmol(model::Model{T, Triangle{T}}) where T
    [∂ₙφmol(ξ, model.charges) for ξ in model.elements]
end

# =========================================================================================
"""
    ∇φmol{T}(ξ::Vector{T}, charges::Vector{Charge{T}})

Computes and returns the gradient of the given system's molecular potential in a
structureless medium for the given observation point ξ.

```math
∇φ_{mol}(ξ) = -\\frac{1}{4π ε_0 ε_Ω} \\sum_i \\frac{qᵢ}{|rᵢ-ξ|³} (rᵢ-ξ)
```

!!! note
    The return value is premultiplied by ``4π ⋅ ε₀ ⋅ ε_Ω``

## Return type
`Vector{T}`

## Aliases
    ∇φmol{T}(Ξ::Vector{Vector{T}})

Computes the molecular potential gradients for a list of observation points.
"""
function ∇φmol(
        ξ      ::Vector{T},
        charges::Vector{Charge{T}}
    ) where T
    # TODO devectorize!
    -sum([q.val * (ξ - q.pos) / euclidean(ξ, q.pos)^3 for q in charges])
end

function ∇φmol(
        Ξ      ::Vector{Vector{T}},
        charges::Vector{Charge{T}}
    ) where T
    [∇φmol(ξ, charges) for ξ in Ξ]
end


text:
# This file includes code originally part of Julia
# Licensed under MIT License, see LICENSE.md

# The plan is to rewrite all of the functionality to not use the utf8proc library,
# and to use tables loaded up on initialization, as with StringLiterals.jl

# Currently, for characters outside the Latin1 subset, this depends on the following C calls:

# utf8proc_errmsg
# utf8proc_charwidth
# utf8proc_category
# utf8proc_category_string
# utf8proc_tolower
# utf8proc_toupper
# utf8proc_totitle

# For grapheme segmentation, this currently depends on the following 2 C calls:

# utf8proc_grapheme_break
# utf8proc_grapheme_break_stateful

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

utf8proc_error(result) =
    error(unsafe_string(ccall(:utf8proc_errmsg, Cstring, (Cssize_t,), result)))

utf8proc_charwidth(ch) = Int(ccall(:utf8proc_charwidth, Cint, (UInt32,), ch))

@inline utf8proc_cat(ch) = ccall(:utf8proc_category, Cint, (UInt32,), ch)
@inline utf8proc_cat_abbr(ch) =
    unsafe_string(ccall(:utf8proc_category_string, Cstring, (UInt32,), ch))

#@inline _lowercase_u(ch) = ccall(:utf8proc_tolower, UInt32, (UInt32,), ch)
#@inline _uppercase_u(ch) = ccall(:utf8proc_toupper, UInt32, (UInt32,), ch)
#@inline _titlecase_u(ch) = ccall(:utf8proc_totitle, UInt32, (UInt32,), ch)

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

# iterators for grapheme segmentation

is_grapheme_break(c1::CodeUnitTypes, c2::CodeUnitTypes) =
    !(c1 <= 0x10ffff && c2 <= 0x10ffff) ||
    ccall(:utf8proc_grapheme_break, Bool, (UInt32, UInt32), c1, c2)

# Stateful grapheme break required by Unicode-9 rules: the string
# must be processed in sequence, with state initialized to Ref{Int32}(0).
# Requires utf8proc v2.0 or later.
is_grapheme_break!(state::Ref{Int32}, c1::CodeUnitTypes, c2::CodeUnitTypes) =
    ((c1 <= 0x10ffff && c2 <= 0x10ffff)
     ? ccall(:utf8proc_grapheme_break_stateful, Bool, (UInt32, UInt32, Ref{Int32}), c1, c2, state)
     : (state[] = 0; true))

is_grapheme_break(c1::Chr, c2::Chr) = is_grapheme_break(codepoint(c1), codepoint(c2))
is_grapheme_break(state::Ref{UInt32}, c1::Chr, c2::Chr) =
    is_grapheme_break(state, codepoint(c1), codepoint(c2))


text:
abstract type StaticOffsetRange{T,P<:StaticRange{T},F} <: AbstractRange{T} end

@pure indexoffset(r::StaticOffsetRange{T,P,F}) where {T,P,F<:SInteger} = F()::F

first(r::StaticOffsetRange) = first(parent(r))
step(r::StaticOffsetRange) = step(parent(r))
last(r::StaticOffsetRange) = last(parent(r))

firstindex(r::StaticOffsetRange) = firstindex(parent(r)) + indexoffset(r)
lastindex(r::StaticOffsetRange) = lastindex(parent(r)) + indexoffset(r)

@inline Base.@propagate_inbounds function getindex(r::StaticOffsetRange, i::Integer)
    Base.@_inline_meta
    parent(r)[i - indexoffset(r)]
end

Base.show(io::IO, r::StaticOffsetRange) = print(io, parent(r))

struct OffsetSRange{T,P,F} <: StaticOffsetRange{T,P,F}
    function OffsetSRange{T,P,F}() where {T,P,F}
        check_params(OffsetSRange{T,P,F})
        new{T,P,F}()
    end
end

@pure parent(r::OffsetSRange{T,P}) where {T,P} = P()::P

mutable struct OffsetMRange{T,P,F} <: StaticOffsetRange{T,P,F}
    parent::P
    offset::F
end

indexoffset(r::OffsetMRange{T,P,F}) where {T,P,F<:BaseInteger} = r.offset
parent(r::OffsetMRange{T,P}) where {T,P} = r.parent::P

function check_params(::Type{OffsetSRange{T,P,F}}) where {T,P,F}
    isstatic(P) || error("OffsetSRange only supports static ranges, got range of type $(nameof(T))")
end