text:
__precompile__(true)

"""
    NKFtool
[Julia](https://julialang.org) package to guess and convert
encodings of Japanese characters.
This is a wrapper
to Network kanji filter, aka [nkf](https://osdn.net/projects/nkf/),
and provides the interface to use nkf command installed in your system from Julia.

NKFtool requires nkf to be installed in your system.
Major Unix-like operating systems offer its precompiled version
in their package management systems.
"""
module NKFtool

export nkf_version, nkf_help, nkf_guess, nkf_convert

try
    success(`nkf -v`)
    global foundNKF = true
catch
    global foundNKF = false
end

if (!foundNKF)
    println("NKFtool cannot be loaded: nkf is not available on this system.")
end

function __init__()

    try
        global NKF_VERSION = begin
            output = IOBuffer()
            run(pipeline(`nkf --version`, stdout = output))
            String(take!(output))
            # close(output)
        end
    catch
		global NKF_VERSION = ""
    end

    try
        global NKF_HELP = begin
            output = IOBuffer()
            run(pipeline(`nkf --help`, stdout = output))
            String(take!(output))
            # close(output)
        end
    catch
		global NKF_HELP = ""
    end
end

"""
    nkf_version()

Print the version message of nkf command.
The same as the output of `nkf --version` command.

# Examples
```julia-repl
julia> nkf_version() |> print
Network Kanji Filter Version 2.1.5 (2018-12-15)
Copyright (C) 1987, FUJITSU LTD. (I.Ichikawa).
Copyright (C) 1996-2018, The nkf Project.
```
"""
function nkf_version()
    NKF_VERSION
end

"""
    nkf_help()

Print the help message of nkf command.
The same as the output of `nkf --help` command.

# Examples
```julia-repl
julia> nkf_help() |> print
```
"""
function nkf_help()
    NKF_HELP
end

"""
    nkf_guess(from::String)

Try to guess the encoding of the input text `from`, and
return a string representing its encoding,
which is just the result of the command line `echo <from> | nkf -g`.

# Examples
```julia-repl
julia> nkf_guess(raw"こんにちわ")
"UTF-8"

julia> nkf_convert( raw"こんにちわ", "-j") |> nkf_guess
"ISO-2022-JP"

julia> nkf_convert( raw"こんにちわ", "-e") |> nkf_guess
"EUC-JP"

julia> nkf_convert( raw"こんにちわ", "-s") |> nkf_guess
"Shift_JIS"
```
"""
function nkf_guess(from::String)
    nkf_guess(IOBuffer(from))
end

"""
    nkf_guess(from::IO)

Try to guess the encoding of the input stream `from`, and
return a string representing its encoding.

# Examples
```julia-repl
julia> nkf_guess(IOBuffer(raw"こんにちわ"))
"UTF-8"

julia> open("hello_sjis.txt","w") do f
           print(f, nkf_convert(raw"こんにちわ", "-s"))
       end
       #
       encoding=open("hello_sjis.txt") do f
           nkf_guess(f)
       end
"Shift_JIS"
```
"""
function nkf_guess(from::IO)
    output = IOBuffer()
    run(pipeline(`nkf -g`, stdin = from, stdout = output))
    chomp(String(take!(output)))
end

"""
    nkf_convert(from::String, options="-w -m0")

Convert the input string `from` to the encoding
specified by the option directive `options`,
and return the output text string,
which is just the result of the command line `echo <from> | nkf <options>`.

# Arguments
- `from::String`: the input string
- `options::String`: the option directive to be passed to nkf command.

    * Output encoding
        * `-j` : ISO-2022-JP
        * `-s` : Shift_JIS
        * `-e` : EUC-JP
        * `-w[8[0],{16,32}[{B,L}[0]]]` : UTF with options
    * Input encoding
        * `-J` : ISO-2022-JP
        * `-S` : Shift_JIS
        * `-E` : EUC-JP
        * `-W[8,[16,32][B,L]]` : UTF with option
    * MIME decode : `-m[BQSN0]`
        * B:base64
        * Q:quoted
        * S:strict
        * N:nonstrict
        * 0:no decode
    * MIME encode : `-M[BQ]`
        * B:base64
        * Q:quoted

# Examples
```julia-repl
julia> nkf_convert(raw"こんにちわ", "-w -m0")
"こんにちわ"

julia> using Base64

julia> nkf_convert( raw"こんにちわ", "-j") |> base64encode
"GyRCJDMkcyRLJEEkbxsoQg=="

julia> String(base64decode(ans)) |> nkf_convert
"こんにちわ"
```
"""
function nkf_convert(from::String, options="-w m0")
    nkf_convert(IOBuffer(from), options)
end

"""
    nkf_convert(from::IO, options="-w -m0")

Convert the input stream `from` to the encoding
specified by the option directive `options`,
and return the output text string,
which is just the result of the command line `cat <from> | nkf <options>`

# Arguments
- `text::String`: the input string
- `options::String`: the directibr to be passed to nkf command. See [`nkf_convert(from::String, options="-w -m0")`](@ref)

# Examples
```julia-repl
julia> open("hello_sjis.txt","w") do f
           print(f, nkf_convert(raw"こんにちわ", "-s"))
       end
       #
       hello_utf=open("hello_sjis.txt") do f
           nkf_convert(f, "-w -m0")
       end
"こんにちわ"
```
"""
function nkf_convert(from::IO, options="-w -m0")
    output = IOBuffer()
    run(pipeline(`nkf $options`, stdin = from, stdout = output))
    String(take!(output))
end

end # module


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*csc(x))^-1*cos(x)^4, x) == :(-1//3*a^-1*cos(x)^3+-1//8*x*a^-1+-1//8*a^-1*cos(x)*sin(x)+1//4*a^-1*cos(x)^3*sin(x))
@test integrate((a+a*csc(x))^-1*cos(x)^3, x) == :((1/2)*a^-1*sin(x)^2+-1//3*a^-1*sin(x)^3)
@test integrate((a+a*csc(x))^-1*cos(x)^2, x) == :(-1*a^-1*cos(x)+-1//2*x*a^-1+(1/2)*a^-1*cos(x)*sin(x))
@test integrate((a+a*csc(x))^-1*cos(x), x) == :(a^-1*sin(x)+-1*a^-1*log(1+sin(x)))
@test integrate((a+a*csc(x))^-1*sec(x), x) == :((1/2)*a^-1*sec(x)^2+(1/2)*a^-1*arctanh(sin(x))+-1//2*a^-1*sec(x)*tan(x))
@test integrate((a+a*csc(x))^-1*sec(x)^2, x) == :(-1//3*a^-1*tan(x)^3+1//3*a^-1*sec(x)^3)
@test integrate((a+a*csc(x))^-1*sec(x)^3, x) == :(1//4*a^-1*sec(x)^4+1//8*a^-1*arctanh(sin(x))+-1//4*a^-1*sec(x)^3*tan(x)+1//8*a^-1*sec(x)*tan(x))
@test integrate((a+a*csc(x))^-1*sec(x)^4, x) == :(-1//3*a^-1*tan(x)^3+-1//5*a^-1*tan(x)^5+1//5*a^-1*sec(x)^5)
@test integrate((a+b*csc(x))^-1*cos(x)^4, x) == :(-1//8*a^-4*(8*b*(a^2+-1*b^2)+-1*a*(-4*b^2+3*a^2)*sin(x))*cos(x)+-1//12*a^-2*cos(x)^3*(4b+-3*a*sin(x))+1//8*x*a^-5*(3*a^4+8*b^4+-12*a^2*b^2)+2*b*a^-5*(a^2+-1*b^2)^3//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*cos(x)^3, x) == :(-1//3*a^-1*sin(x)^3+a^-3*(a^2+-1*b^2)*sin(x)+(1/2)*b*a^-2*sin(x)^2+-1*b*a^-4*(a^2+-1*b^2)*log(b+a*sin(x)))
@test integrate((a+b*csc(x))^-1*cos(x)^2, x) == :((1/2)*x*a^-3*(a^2+-2*b^2)+-1//2*a^-2*(2b+-1*a*sin(x))*cos(x)+2*b*a^-3*(a^2+-1*b^2)^(1/2)*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*cos(x), x) == :(a^-1*sin(x)+-1*b*a^-2*log(b+a*sin(x)))
@test integrate((a+b*csc(x))^-1*sec(x), x) == :((-2b+2a)^-1*log(1+sin(x))+-1*(2a+2b)^-1*log(1+-1*sin(x))+-1*b*(a^2+-1*b^2)^-1*log(b+a*sin(x)))
@test integrate((a+b*csc(x))^-1*sec(x)^2, x) == :(-1*(a^2+-1*b^2)^-1*(b+-1*a*sin(x))*sec(x)+2*a*b*(a^2+-1*b^2)^-3//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*sec(x)^3, x) == :(-1*(-2*b^2+2*a^2)^-1*sec(x)^2*(b+-1*a*sin(x))+-1//4*a*(a+b)^-2*log(1+-1*sin(x))+1//4*a*(a+-1b)^-2*log(1+sin(x))+-1*b*a^2*(a^2+-1*b^2)^-2*log(b+a*sin(x)))
@test integrate((a+b*csc(x))^-1*sec(x)^4, x) == :(-1*(-3*b^2+3*a^2)^-1*sec(x)^3*(b+-1*a*sin(x))+-1//3*(a^2+-1*b^2)^-2*(3*b*a^2+-1*a*(b^2+2*a^2)*sin(x))*sec(x)+2*b*a^3*(a^2+-1*b^2)^-5//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))

ext: jl

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.1.3 (d cos)^n (a+b csc)^m.jl


text:
"""Three variables come into play in a wavelet filterbank:
* the spatial variable `V`, which may be multi-dimensional ;
* the scale variable `γ::V`, which is one-dimensional and non-negative ;
* the isometric variable `θ::V`, which belongs to a ""point group"".
A point group is an algebraic group of transformations which have a single
fixed point. Reflections of the real line and rotations of the plane are
examples of such point groups.
Conversely, translations and scalings are not point groups.
"""
abstract AbstractPointGroup

"""The `TrivialGroup` is a finite group of a single neutral element."""
immutable TrivialGroup <: AbstractPointGroup end

"""The `ReflectionGroup` is a finite group of two elements, which represents
the reflections of the real line."""
immutable ReflectionGroup <: AbstractPointGroup end

"""The `RotationGroup` is a continuous group, which represents the rotations of
the plane. Its orientations are discretized uniformly."""
immutable RotationGroup <: AbstractPointGroup
    n_orientations::Int
end

typealias LineGroups Union{TrivialGroup,ReflectionGroup}
typealias PlaneGroups Union{TrivialGroup,RotationGroup}

"""Returns the number of orientations associated to a point group.
Equal to `1` for `TrivialGroup`, `2` for `ReflectionGroup`,
and `group.n_orientations` for `RotationGroup`."""
get_n_orientations(group::TrivialGroup) = 1
get_n_orientations(group::ReflectionGroup) = 2
get_n_orientations(group::RotationGroup) = group.n_orientations


text:
chargers = sort(parse.(Int, split(strip(read("input.txt", String)), "\n")))

diffs(chargers) = map(t -> t[2] - t[1], zip(chargers[1:end - 1], chargers[2:end]))

joltages = [0; chargers; chargers[end] + 3]
differences = diffs(joltages)

count(==(1), differences) * count(==(3), differences)

examplechargers = [0, 1, 4, 5, 6, 7, 10, 11, 12, 15, 16, 19, 22]
examplediffs = [1, 3, 1, 1, 1, 3, 1, 1, 3, 1, 3, 3]
exampleanswer = 8

otherexample = [0, 1, 2, 3, 4, 7, 8, 9, 10, 11, 14, 17, 18, 19, 20, 23, 24, 25, 28, 31,
32, 33, 34, 35, 38, 39, 42, 45, 46, 47, 48, 49, 52]
otherdiffs = [1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 3, 1, 1, 1, 3, 1, 1, 3, 3, 1, 1, 1, 1, 3, 1, 3, 3, 1, 1, 1, 1, 3]
otheranswer = 19208

"""Splits a full list of differences into sublists for easier calculation"""
function calculatebranches(diffs)
    iter = Iterators.Stateful(diffs)
    arglist = []
    while !isempty(iter)
        args = [Iterators.takewhile(<(3), iter)...]
        push!(arglist, args)
        Iterators.takewhile(==(3), iter)
    end
    mapreduce(calculatesubbranch, *, filter(!isempty, arglist))
end

"""Takes a sublist of diffs that are less than 3 and calculates the possible permutations of the subbranch"""
function calculatesubbranch(diffs)
    if length(diffs) == 1
        return 1
    elseif length(diffs) == 2
        return 2
    elseif length(diffs) == 3
        return 4
    elseif length(diffs) == 4
        return 7
    end
end

@assert calculatebranches(examplediffs) == exampleanswer
@assert calculatebranches(otherdiffs) == otheranswer
calculatebranches(differences) # holy shit it actually works

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
#!/opt/julia0.6/julia
using DFWannier
using CUDAnative
using CuArrays
CUDAnative.@profile begin
  T=Float64
  test_wfc1 =  DFWannier.host2gpu(read_xsf_file("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/wan_00003.xsf",Atom{T}(0.0,0.0,0.1,0.1),T))
  test_wfc2=  DFWannier.host2gpu(read_xsf_file("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/wan_00004.xsf",Atom{T}(0.0,0.0,0.1,0.1),T))
  dims = size(test_wfc1.values)
  cu_dims = CuArray(UInt32[dims...])
  Lx = CuArray(zeros(Complex{T},dims))
  Ly = CuArray(zeros(Complex{T},dims))
  Lz = CuArray(zeros(Complex{T},dims))
  n1 = CuArray(zeros(Complex{T},dims))
  n2 = CuArray(zeros(Complex{T},dims))
  grid = Array(test_wfc1.grid)
  origin = grid[1,1,1]
  a = grid[2,1,1] .- origin
  b = grid[1,2,1] .- origin
  c = grid[1,1,2] .- origin
  V = CuArray(inv([[a...] [b...] [c...]]))
  indices = Array{Tuple{Tuple{Int32,Int32,Int32},Tuple{Int32,Int32,Int32}},1}()
  t_coeff = Array{Complex{T}}(27)
  t=1
  k=T[0.2,0.2,0.3]
  for R1=-1:1,R2=-1:1,R3=-1:1
    R= R1*test_wfc1.cell[1]+R2*test_wfc1.cell[2]+R3*test_wfc1.cell[3]
    ind1,ind2 = DFWannier.find_start(test_wfc1,R,27)
    push!(indices,((Int32(ind1[1]),Int32(ind1[2]),Int32(ind1[3])),(Int32(ind1[1]-ind2[1]),Int32(ind1[2]-ind2[2]),Int32(ind1[3]-ind2[3]))))
    t_coeff[t] = Complex{T}(exp(dot(-2*pi*k,[R1,R2,R3])*1im))
    t+=1
  end
  coefficients = CuArray(t_coeff)
  cu_indices = CuArray(indices)
  k_wfcs = Array{Wfc3D_gpu{T},1}(2)
  for (n,wfc) in enumerate([test_wfc1,test_wfc2])
    k_wfcs[n] = Wfc3D_gpu(wfc.grid,CuArray(zeros(Complex{T},size(wfc.values))),wfc.cell,wfc.atom)
  end
  # DFWannier.construct_bloch_sums([test_wfc1,test_wfc2],k_wfcs,k,cu_indices,coefficients)
  begin
  correct = calculate_angmom(test_wfc1,test_wfc2,V,CuArray([test_wfc1.atom.center.x,test_wfc1.atom.center.y,test_wfc1.atom.center.z]),dims,Lx,Ly,Lz,n2,n2)
  test = calculate_angmom(test_wfc1,test_wfc2,V,CuArray([test_wfc1.atom.center.x,test_wfc1.atom.center.y,test_wfc1.atom.center.z]),dims,Lx,Ly,Lz,n2,n2)
  @time for i =1:2000
    test = calculate_angmom(test_wfc1,test_wfc2,V,CuArray([test_wfc1.atom.center.x,test_wfc1.atom.center.y,test_wfc1.atom.center.z]),dims,Lx,Ly,Lz,n2,n2)

  end
  assert(test==correct)
  end
end


text:
module Elvur

using StructTypes
using JSON3
using SQLite
using DataFrames
using Logging

io = open("log.txt", "w+")
logger = SimpleLogger(io)
global_logger(logger)

include("types/types.jl")

allitems() = AllItems(JSON3.read(read("D:\\Code\\.julia\\dev\\Elvur\\data\\items.json", String), Vector{Item}))

include("data/data.jl")
include("main.jl")

end


text:
include("lpc.jl")

infile = "in.wav"
if length(ARGS) >= 1
    infile = ARGS[1]
end

pcm, fs = wavread(infile)
xs = pcm[:, 1]

out, n, frameTypes = lpc_encode_decode(xs, fs)

outfile = splitext(infile)[1] * "_lp.wav"
dump_voicing_labels(out, fs, n, frameTypes, outfile)

outfile = splitext(infile)[1] * "_labelled.wav"
dump_voicing_labels(xs, fs, n, frameTypes, outfile)


text:
###################################################
## networkposition.jl
##     Location information with respect to the network (projection onto the network...)
## Authors: Sébastien Martin, 2017
###################################################

"""
    Types that represent positions in the network, always contain.
"""
abstract type NetworkPosition end

" Returns the exact latitude of the position"
lat(pos::NetworkPosition) = pos.lat
" Returns the exact longitude of the position"
lon(pos::NetworkPosition) = pos.lon
" Returns the x position in the ENU coordinate system of the network"
x(pos::NetworkPosition) = pos.x
" Returns the y position in the ENU coordinate system of the network"
y(pos::NetworkPosition) = pos.y
" Return the id of the best node approximation of the position"
node(pos::NetworkPosition) = error("The `node` function has not been written for this NetworkPosition")
" Returns the distance from the node given by the `node` function"
distToProj(n::Network, pos::NetworkPosition) = error("The `distToProj` function has not been written for this NetworkPosition")

distanceGeo(p1::NetworkPosition, p2::NetworkPosition) =
    distanceGeo(lon(p1), lat(p1), lon(p2), lat(p2))

"""
    A simple nearest-node projection to the network. Fast to compute.
"""
struct NodePosition <: NetworkPosition
    "Latitude"
    lat::Float64
    "Longitude"
    lon::Float64
    "x position (planar representation)"
    x::Float64
    "y position (planar representation)"
    y::Float64
    "node id in the routing graph (0 if not computed)"
    node::Int
end

node(pos::NodePosition) = pos.node > 0 ? pos.node : error("Node not computed")

function distToProj(n::Network, pos::NodePosition)
    node = n.nodes[pos.node]
    distanceGeo(node.lon, node.lat, pos.lon, pos.lat)
end

"""
    KD-tree wrapper to project real locations onto network, (on nearest nodes only)
"""
struct NodeProjector
    "The routing network"
    network::Network
    "The KD-tree of nodes"
    tree::KDTree
end
function NodeProjector(n::Network)
    # Constructing tree
    nodePos = Array(Float64,(2,length(n.nodes)))
    for (i,node) in enumerate(n.nodes)
       nodePos[1,i] = node.x
       nodePos[2,i] = node.y
    end
    return NodeProjector(n,KDTree(nodePos))
end

"""
    Returns the nearest node id with respect to a x/y position
"""
nearestNode(proj::NodeProjector, x, y) = knn(proj.tree,[x,y],1)[1][1]

"""
    Returns the Network Position, given latitude and longitude
"""
function NetworkPosition(proj::NodeProjector, lat, lon)
    x, y = toENU(lon, lat, proj.network)
    return NodePosition(lat, lon, x, y, nearestNode(proj, x, y))
end

"""
    A more advanced projection to the network. project on a road and give position on this road.
"""
struct RoadPosition <: NetworkPosition
    "Latitude"
    lat::Float64
    "Longitude"
    lon::Float64
    "x position (planar representation)"
    x::Float64
    "y position (planar representation)"
    y::Float64
    "node id of the origin of the road"
    nodeO::Int
    "node id of the destination of the road"
    nodeD::Int
    "fraction representing the position in between the node. 0=nodeO, 1.=nodeD"
    fraction::Float64
end

function node(pos::RoadPosition)
    if pos.fraction <=0.5
        return pos.nodeO
    else
        return pos.nodeD
    end
end
"""
    Each road is discretized into regularly spaced sub-nodes, used by the RoadProjector
"""
struct RoadSubNode
    "Origin node id"
    nodeO::Int
    "Destination node id"
    nodeD::Int
    "x position"
    x::Float64
    "y position"
    y::Float64
end

"""
    More fancy projector, that project onto the exact nearest road.
"""
struct RoadProjector
    "The routing network"
    network::Network
    "The KD-tree for all RoadSubNodes"
    tree::KDTree
    "The list of all subnodes"
    subnodes::Vector{RoadSubNode}
    "The discretization parameter, in meters: the max distance between 2 subnodes of same road"
    maxDistance::Float64
end

function RoadProjector(n::Network; maxDistance::Float64=50., noHighway::Bool=false)
    if noHighway # remove highways from network
        n = roadTypeSubset(n, 3:6)
    end
    # Creating the subnodes.
    subnodes = RoadSubNode[]
    uniqueArcs = Set{Tuple{Int, Int}}()
    for r in values(n.roads)
        od = minmax(r.orig, r.dest)
        if !(od in uniqueArcs)
            push!(uniqueArcs, od)
            nSubNodes = max(1,ceil(Int, r.distance/maxDistance))+1
            nodeO = n.nodes[r.orig]
            nodeD = n.nodes[r.dest]
            for i = 1:nSubNodes
                frac = (i-1)/(nSubNodes-1)
                x = nodeO.x + frac*(nodeD.x-nodeO.x)
                y = nodeO.y + frac*(nodeD.y-nodeO.y)
                push!(subnodes, RoadSubNode(r.orig, r.dest, x, y))
            end
        end
    end
    # Constructing tree
    nodePos = Array{Float64}(2,length(subnodes))
    for (i,subnode) in enumerate(subnodes)
       nodePos[1,i] = subnode.x
       nodePos[2,i] = subnode.y
    end
    return RoadProjector(n,KDTree(nodePos), subnodes, maxDistance)
end

"""
    Returns the Network Position, given latitude and longitude
    Use the subnodes to project on the road quickly.
"""
function NetworkPosition(proj::RoadProjector, lat, lon)
    x, y = toENU(lon, lat, proj.network)
    id, dist = knn(proj.tree,[x,y],1)
    closestSubNode = proj.subnodes[id[1]]

    candidateSubNodes = inrange(proj.tree, [x, y], dist[1] + proj.maxDistance)

    candidateRoads = Tuple{Int, Int}[] # using a vector as very small sets are inefficient
    for subNodeId in candidateSubNodes
        subnode = proj.subnodes[subNodeId]
        od = (subnode.nodeO, subnode.nodeD)
        if !(od in candidateRoads)
            push!(candidateRoads, od)
        end
    end
    length(candidateRoads) > 0 || error("Something is wrong here.")
    bestOD = (0, 0)
    bestFrac = -1.
    bestDist = Inf
    otherSide = false
    for od in candidateRoads
        dist, frac, side = roadProjection(proj.network.nodes[od[1]], proj.network.nodes[od[2]], x, y)
        if dist < bestDist
            bestOD = od
            bestFrac = frac
            bestDist = dist
            otherSide = side
        end
    end
    (o,d) = bestOD
    if otherSide && haskey(proj.network.roads, (d,o)) # Select the road direction
        bestFrac = 1.-bestFrac
        (o,d) = (d,o)
    end

    return RoadPosition(lat, lon, x, y, o, d, bestFrac)
end

"""
    Project a position onto a road segment.
"""
function roadProjection(nO::Node, nD::Node, x::Float64, y::Float64)
    O = Float64[nO.x, nO.y]
    D = Float64[nD.x, nD.y]
    p = Float64[x,y]

    OD = D - O
    ODsquared = dot(OD,OD)
    ODsquared == 0. && error("origin and destination are the same")

    Op = p - O
    # from http://stackoverflow.com/questions/849211/
    # Consider the line extending the segment, parameterized as O + frac (D - O)
    # We find projection of point p onto the line.
    # It falls where frac = [(p-O) . (D-O)] / |D-O|^2\
    frac = dot(Op,OD)/ODsquared

    if (frac < 0.0)
        frac = 0.
    elseif (frac > 1.0)
        frac = 1.
    end

    proj = O + frac * OD

    rotOD = [-OD[2], OD[1]] # rotate 90 degrees counter-clockwise
    otherSide = dot(rotOD, p-proj) > 0 # test the side of the street
    dist = norm(p-proj)

    return dist, frac, otherSide
end

function distToProj(n::Network, pos::RoadPosition)
    nO = n.nodes[pos.nodeO]
    nD = n.nodes[pos.nodeD]
    projlat = nO.lat + pos.fraction*(nD.lat-nO.lat)
    projlon = nO.lon + pos.fraction*(nD.lon-nO.lon)
    return distanceGeo(projlon, projlat, pos.lon, pos.lat)
end

function projXY(n::Network, pos::RoadPosition)
    nO = n.nodes[pos.nodeO]
    nD = n.nodes[pos.nodeD]
    x = nO.x + pos.fraction*(nD.x-nO.x)
    y = nO.y + pos.fraction*(nD.y-nO.y)
    return x,y
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
using OrdinaryDiffEq, DiffEqDevTools, Sundials, ModelingToolkit, ODEInterfaceDiffEq,
      Plots, DASSL, DASKR
using LinearAlgebra

ModelingToolkit.@parameters begin
      k₁=18.7
      k₂=0.58
      k₃=0.09
      k₄=0.42
      kbig=34.4
      kla=3.3
      ks=115.83
      po2=0.9
      hen=737
end

@variables begin
      t
      y₁(t)  = 0.444
      y₂(t)  = 0.00123
      y₃(t)  = 0.0
      y₄(t)  = 0.007
      y₅(t)  = 1.0
      y₆(t)  = 115.83*0.444*0.007 # ks*y₁*y₄
end

D = Differential(t)

r₁  = k₁ * (y₁^4.)*sqrt(abs(y₂))
r₂  = k₂ * y₃ * y₄
r₃  = k₂/kbig * y₁ * y₅
r₄  = k₃*y₁*(y₄^2)
r₅  = k₄*(y₆^2)*sqrt(abs(y₂))
fin = kla*(po2/hen-y₂)

eqs = [
      D(y₁) ~ -2. * r₁ + r₂ - r₃ - r₄
      D(y₂) ~ -0.5 * r₁ - r₄ - 0.5*r₅ + fin
      D(y₃) ~ r₁ - r₂ + r₃
      D(y₄) ~ -r₂ + r₃ - 2. * r₄
      D(y₅) ~ r₂ - r₃ + r₅
      0. ~ ks * y₁ * y₄ - y₆
]

ModelingToolkit.@named sys = ModelingToolkit.ODESystem(eqs)

tspan = (0.0, 180.0)
mmprob = ODEProblem(sys, [], tspan)
sol = solve(mmprob, Rodas4(),abstol=1/10^14,reltol=1/10^14)

odaeprob = ODAEProblem(structural_simplify(sys),[],tspan)
ode_ref_sol = solve(odaeprob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14);

du = mmprob.f(mmprob.u0,mmprob.p,0.0)
du0 = D.(states(sys)) .=> du
daeprob = DAEProblem(sys,du0,[],tspan)
ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);

probs = [mmprob,daeprob,odaeprob]
refs = [ref_sol,ref_sol,ode_ref_sol];

plot(ode_ref_sol, vars = [y₁,y₂,y₃,y₄,y₅,y₆])

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
]

wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (2:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (3:5);
setups = [Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 2, :alg=>DASSL.dassl()),
          Dict(:prob_choice => 2, :alg=>DASKR.daskr()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()), # too slow
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (6:8)
reltols = 1.0 ./ 10.0 .^ (2:4);
setups = [Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 3, :alg=>Rosenbrock23()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>CVODE_BDF()),
          Dict(:prob_choice => 3, :alg=>TRBDF2()),
          Dict(:prob_choice => 3, :alg=>KenCarp4()),
          ]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:12)
reltols = 1.0 ./ 10.0 .^ (4:9)

setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
          Dict(:prob_choice => 3, :alg=>Rodas5()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          #Dict(:prob_choice => 1, :alg=>FBDF()),
          Dict(:prob_choice => 1, :alg=>QNDF()),
          Dict(:prob_choice => 1, :alg=>rodas()),
          Dict(:prob_choice => 1, :alg=>radau()),
          Dict(:prob_choice => 1, :alg=>RadauIIA5()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          ]

wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(probs,abstols,reltols,setups;error_estimate = :l2,
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:9)
reltols = 1.0 ./ 10.0 .^ (4:6)

setups = [Dict(:prob_choice => 1, :alg=>Rodas5()),
          Dict(:prob_choice => 3, :alg=>Rodas5()),
          Dict(:prob_choice => 1, :alg=>Rodas4()),
          Dict(:prob_choice => 3, :alg=>Rodas4()),
          Dict(:prob_choice => 2, :alg=>DFBDF()),
          Dict(:prob_choice => 2, :alg=>IDA()),
          Dict(:prob_choice => 2, :alg=>DASKR.daskr()),
          ]
gr()
wp = WorkPrecisionSet(probs,abstols,reltols,setups;
                      save_everystep=false,appxsol=refs,maxiters=Int(1e5),numruns=10)
plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
module Assignment05

export isreverse,
       complement,
       isreversecomplement,
       isreversecomplement,
       reverse_complement,
       find_dna

include("assignment.jl")

end


text:
module MaxHelpingHandAllSideTentacles

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/AllSideTentacles" AllSideTentacles(only::String="*", exclude::String="", side::String="Right", color::String="", offset::Number=0)

placements = AllSideTentacles

function Ahorn.canFgBg(effect::AllSideTentacles)
    return true, true
end

function Ahorn.editingOptions(effect::AllSideTentacles)
    return Dict{String, Any}(
        "side" => String["Top", "Bottom", "Left", "Right"]
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using ..SymEngine: Basic
export simplify_expi

function Base.iszero(x::Basic)
    isempty(free_symbols(x)) && iszero(N(x))
end

SymEngine.free_symbols(syms::Union{Real,Complex}) = Basic[]
SymEngine.free_symbols(syms::AbstractArray{T}) where {T<:Union{Real,Complex}} = Basic[]
SymEngine.free_symbols(syms::AbstractArray{T}) where {T<:Union{Basic,SymEngine.BasicType}} =
    SymEngine.free_symbols(Matrix(syms))
function rand_assign(syms...)
    fs = union(free_symbols.(syms)...)
    Dict(zip(fs, randn(length(fs))))
end

function _basic_approx(x, y; atol = 1e-8)
    diff = x - y
    assign = rand_assign(x, y)
    length(assign) > 0 && (diff = subs.(diff, Ref.((assign...,))...))
    nres = ComplexF64.(diff)
    all(isapprox.(nres, 0; atol = atol))
end

Base.:≈(x::AbstractArray{<:Basic}, y::AbstractArray; atol = 1e-8) = _basic_approx(x, y, atol = atol)
Base.:≈(x::AbstractArray, y::AbstractArray{<:Basic}; atol = 1e-8) = _basic_approx(x, y, atol = atol)
Base.:≈(x::AbstractArray{<:Basic}, y::AbstractArray{<:Basic}; atol = 1e-8) =
    _basic_approx(x, y, atol = atol)
function Base.Complex{T}(a::Basic) where {T}
    a = SymEngine.evalf(a)
    T(real(a)) + im * T(imag(a))
end

# simplification rules
function simplify_expi(expr::Basic)
    subs(
        expr,
        exp(im * Basic(π) / 2) => Basic(im),
        exp(im * Basic(π)) => -1,
        exp(im * Basic(π) * 3 / 2) => -Basic(im),
    )
end

Broadcast.broadcasted(
    ::Broadcast.AbstractArrayStyle{2},
    ::typeof(_pretty_basic),
    a::PermMatrix,
) = PermMatrix(a.perm, _pretty_basic.(a.vals))

Broadcast.broadcasted(
    ::LinearAlgebra.StructuredMatrixStyle{Diagonal},
    ::typeof(_pretty_basic),
    a::IMatrix{N,T},
) where {N, T} = IMatrix{N,Basic}()


text:
using Test

@testset "Summary" begin

include("cma_es.jl")
include("collect_points.jl")
include("continuous_time_rnn.jl")
include("dummy_app.jl")

end


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing SLiDE objects.
"""
module SLiDE

#################################################################################
# IMPORTS
import CSV
# import InvertedIndices
import DataFrames;          using DataFrames
import Dates
import DelimitedFiles;      using DelimitedFiles
import Ipopt;               using Ipopt
import JSON;                using JSON
import JuMP;                using JuMP
import Logging;             using Logging
# import PowerSimulations;
# import Printf
# import Query;               using Query
# import Revise
import Statistics;          using Statistics
# import Test;
import XLSX
import YAML

import InteractiveUtils
const IU = InteractiveUtils

import InfrastructureSystems
const IS = InfrastructureSystems

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

# First, generate structs to ensure all exports are possible.
const SLIDE_DIR = abspath(joinpath(dirname(Base.find_package("SLiDE")), ".."))
export SLIDE_DIR

# EXPORTS
export Add
export Combine
export Describe
export Drop
export Group
export OrderedGroup
export Map
export Melt
export Operate
export Order
export Match
export Rename
export Replace
export Deselect
export Stack

export SetInput
export CSVInput
export GAMSInput
export XLSXInput
export DataInput

export FileInput

export Parameter

export EconomicSystemsType
export CGE
export DataStream
export Check
export Edit
export File

# UTILITIES
export append
export convert_type
export datatype
export dropnan!
export dropnan
export dropzero!
export dropzero
export dropvalue!
export dropvalue
export ensurearray
export ensuretuple
export ensurefinite
export find_oftype
# export isarray
# export istype
export permute
# export add_permutation!
export findindex
export findvalue
export findunits
export indexjoin
export convertjoin
export propertynames_with

# EDIT
export edit_with
export fill_zero
export fill_with
export extrapolate_region
export extrapolate_year
export filter_with

# CALCULATE
export combine_over
export transform_over

# READ
export read_file
export read_from
export load_from
export write_yaml
export run_yaml

# CHECK
export compare_summary
export compare_keys
export compare_values
export verify_over
export benchmark_against

# BUILD
export build
export partition
export calibrate
export share
export disagg

export module_energy!
export module_elegen!
export module_co2emis!

# MODEL
export model_input

#################################################################################
# INCLUDES
"""
* See [PowerSystems.jl](https://github.com/NREL/PowerSystems.jl/blob/master/src/PowerSystems.jl).*
Supertype for all SLiDE types.
All subtypes must include a InfrastructureSystemsInternal member.
Subtypes should call InfrastructureSystemsInternal() by default, but also must
provide a constructor that allows existing values to be deserialized.
"""
# abstract type EconomicSystemsType <: IS.InfrastructureSystemsType end
abstract type EconomicSystemsType end

abstract type DataStream <: EconomicSystemsType end
abstract type Edit <: DataStream end
abstract type File <: DataStream end
abstract type Check <: DataStream end

abstract type CGE <: EconomicSystemsType end

# CONSTANTS
include("definitions.jl")
export SUB_ELAST
export TRANS_ELAST
export MODEL_LOWER_BOUND

# TYPES
include(joinpath("parse", "generated_check", "includes.jl"))
include(joinpath("parse", "generated_edit", "includes.jl"))
include(joinpath("parse", "generated_file", "includes.jl"))

include(joinpath("model", "generated_cge", "includes.jl"))

# UTILITIES
include(joinpath("utils", "calc.jl"))
include(joinpath("utils", "fill_zero.jl"))
include(joinpath("utils", "indexjoin.jl"))
include(joinpath("utils", "label.jl"))
include(joinpath("utils", "utils.jl"))

include(joinpath("parse", "edit_with.jl"))
include(joinpath("parse", "filter_with.jl"))
include(joinpath("parse", "load_from.jl"))
include(joinpath("parse", "read_file.jl"))
include(joinpath("parse", "run_yaml.jl"))
include(joinpath("parse", "check_data.jl"))

include(joinpath("build", "build.jl"))
include(joinpath("build", "partition.jl"))
include(joinpath("build", "calibrate.jl"))
include(joinpath("build", "share", "share.jl"))
include(joinpath("build", "share", "share_cfs.jl"))
include(joinpath("build", "share", "share_gsp.jl"))
include(joinpath("build", "share", "share_pce.jl"))
include(joinpath("build", "share", "share_sgf.jl"))
include(joinpath("build", "share", "share_utd.jl"))
include(joinpath("build", "disagg", "disagg_region.jl"))

include(joinpath("model", "model_input.jl"))

function __init__()

    last_updated = Dates.DateTime("2021-01-20T00:00:00.0")
    data = joinpath(SLiDE.SLIDE_DIR, "data")
    if isdir(data) && last_updated > Dates.unix2datetime(ctime(data))
        @warn("SLiDE input data has been updated.
            Remove or rename $data and
            rebuild SLiDE (] build) to avoid compatibility issues.")
    end

    # See: http://pages.cs.wisc.edu/~ferris/path/LICENSE
    # https://docs.julialang.org/en/v1/base/base/#Base.ENV
    # https://docs.julialang.org/en/v1/manual/modules/index.html -> init
    Base.ENV["PATH_LICENSE_STRING"] =
        "2830898829&Courtesy&&&USR&45321&5_1_2021&1000&PATH&GEN&31_12_2025&0_0_0&6000&0_0"
end

end # module

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using Revise
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

dir="./test/data/NAMD"
atoms = readPDB("$dir/structure.pdb")

@time options = Options(lastframe=1,seed=321,StableRNG=true,nthreads=1,silent=true)

@time tmao = Selection(select(atoms,"resname TMAO"),natomspermol=14)

# cross-correlation

@time protein = Selection(select(atoms,"protein"),nmols=1)

@time traj = Trajectory("$dir/trajectory.dcd",protein,tmao)

@time samples = CM.Samples(md=(traj.solvent.nmols-1)/2,random=options.n_random_samples)

@time R = Result(traj,options)

@time framedata = CM.FrameData(traj,R)

@time CM.nextframe!(traj)

@time RNG = CM.init_random(options)

@time CM.mddf_frame!(1,framedata,options,RNG,R)

# self

@time traj = Trajectory("$dir/trajectory.dcd",tmao)

@time samples = CM.Samples(md=(traj.solvent.nmols-1)/2,random=options.n_random_samples)

@time R = Result(traj,options)

@time framedata = CM.FrameData(traj,R)

@time CM.nextframe!(traj)

@time RNG = CM.init_random(options)

@time CM.mddf_frame_self!(1,framedata,options,RNG,R)


text:
plot(((Dates.value.(sa.times)).*(10^(-9))).-Dates.value(sa.times[1]).*(10^(-9)), sa.losses, label = "twostage loss over time")
plot(sa.times, sa.losses, label = "l2 loss over time")

savefig("l2_conv.pdf")
plot(((Dates.value.(sa_l2.times)).*(10^(-9))).-Dates.value(sa_l2.times[1]).*(10^(-9)), sa_l2.losses, label = "l2 loss over time")

using BSON,Dates,Plots
#bson("DiffEqFlux.jl/test/results/verification_5000_epochs.bson", Dict(:times => sa.times,:l2s => sa.l2s, :losses => sa.losses))

bson("DiffEqFlux.jl/test/results/mix/sa_mix_"*name*"_epochs.bson", Dict(:times => sa.times, :losses => sa.losses))

l2 = BSON.load("DiffEqFlux.jl/test/results/l2/sa_l2_1000_epochs.bson")
mix = BSON.load("DiffEqFlux.jl/test/results/mix/sa_mix_3000_8_1_40_20_epochs.bson")
two = BSON.load("DiffEqFlux.jl/test/results/twostage/sa_two_169555_epochs.bson")

l2_a =(Dates.value.(l2[:times]).*(10^(-9))).-Dates.value(l2[:times][1]).*(10^(-9))
mix_a =(Dates.value.(mix[:times]).*(10^(-9))).-Dates.value(mix[:times][1]).*(10^(-9))
two_a =(Dates.value.(two[:times]).*(10^(-9))).-Dates.value(two[:times][1]).*(10^(-9))
precompiler_time = mix_a[5001]-mix_a[5000]
for i in 1:length(mix_a)
    if i>5000
        mix_a[i]= mix_a[i]-precompiler_time
    end
end
savefig(string(name,"long.pdf"))
l2_b = l2[:losses]
mix_b = mix[:losses]
two_b = two[:losses]

plot(two_a[1:207:end],log.(two_b[1:207:end].+1),linecolor ="#ff8533",label = "collocation", grid=false, xlabel = "time [sec]", ylabel = " error" )
plot!(l2_a[1:820],log.(l2_b[1:820].+1),label = "L2",linecolor = "#668cff")
plot(mix_a[1:5:5000],log.(mix_b[1:5:5000].+1),label = "two stage",linecolor = "#ffd39b")

length(two_a[1:207:end])
length(l2_a[1:820])
length(mix_a[1:10:8494])
u=((Dates.value.(sa.times)).*(10^(-9))).-Dates.value(sa.times[1]).*(10^(-9))
v=sa.losses

for i in 1:length(mix_a)
    if mix_a[i]>283
        print(i)
        break
    end
end

plot(u[1:207:end],log.(v[1:207:end].+1),linecolor ="#ff8533",label = "collocation", grid=false, xlabel = "time [sec]", ylabel = " error" )
plot!(a[1:820],log.(b[1:820].+1),label = "L2",linecolor = "#668cff")

u[1:207:end]

savefig("280sec_loss_drop_mix.pdf")

plot(u[1:2:end],log.(v[1:2:end].+1),linecolor ="#ff8533",label = "mix", grid=false, xlabel = "time [sec]", ylabel = " error" )

plot(sa.times[1:8000],sa.losses[1:8000])
plot(sa.times[8000:end],sa.losses[8000:end])

log(sa.losses[820]+1)

mix_a =(Dates.value.(sa.times).*(10^(-9))).-Dates.value(sa.times[1]).*(10^(-9))
mix_b = sa.losses

mix_a[8000]
460-mix_a[8001]

## collocation check with L2 number

time =(Dates.value.(sa.times).*(10^(-9))).-Dates.value(sa.times[1]).*(10^(-9))
losses = sa.losses
l2s = sa.l2s
end_epoch = 3800
scatter(l2_a[1:17],log.(l2_b[1:17].+1),markersize=3,label = "L2 model: loss")
plot!(two_a[1:2:end_epoch],log.(losses[1:2:end_epoch].+1),linecolor ="#ff8533",xlabel="time",ylabel="error",label = "Collocation model: loss")
plot!(two_a[1:2:end_epoch],log.(l2s[1:2:end_epoch].+1),linecolor ="green",label = "Collocation model: L2 to observation")
l2 = BSON.load("DiffEqFlux.jl/test/results/l2/sa_l2_1000_epochs.bson")
l2_a =(Dates.value.(l2[:times]).*(10^(-9))).-Dates.value(l2[:times][1]).*(10^(-9))
l2_b = l2[:losses]

two = BSON.load("DiffEqFlux.jl/test/results/twostage/sa_two_169555_epochs.bson")
two_a =(Dates.value.(two[:times]).*(10^(-9))).-Dates.value(two[:times][1]).*(10^(-9))
two_b = two[:losses]
plot!(two_a[1:207:end],log.(two_b[1:207:end].+1),linecolor ="#ff8533",label = "collocation", grid=false, xlabel = "time [sec]", ylabel = " error" )

savefig("col_los_vs_real_obs_3800_with_l2model.pdf")

two_a[3800]

l2_a[17]


text:
module LinearWrapper

using JavaCall

const BaseMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultivariateOptimizer
const BaseOptimizer = @jimport org.hipparchus.optim.BaseOptimizer
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const Enum = @jimport java.lang.Enum
const GoalType = @jimport org.hipparchus.optim.nonlinear.scalar.GoalType
const JString = @jimport java.lang.String
const LinearConstraint = @jimport org.hipparchus.optim.linear.LinearConstraint
const LinearConstraintSet = @jimport org.hipparchus.optim.linear.LinearConstraintSet
const LinearObjectiveFunction = @jimport org.hipparchus.optim.linear.LinearObjectiveFunction
const LinearOptimizer = @jimport org.hipparchus.optim.linear.LinearOptimizer
const List = @jimport java.util.List
const MultivariateOptimizer = @jimport org.hipparchus.optim.nonlinear.scalar.MultivariateOptimizer
const NonNegativeConstraint = @jimport org.hipparchus.optim.linear.NonNegativeConstraint
const Object = @jimport java.lang.Object
const OptimizationData = @jimport org.hipparchus.optim.OptimizationData
const Optional = @jimport java.util.Optional
const PivotSelectionRule = @jimport org.hipparchus.optim.linear.PivotSelectionRule
const PointValuePair = @jimport org.hipparchus.optim.PointValuePair
const RealVector = @jimport org.hipparchus.linear.RealVector
const Relationship = @jimport org.hipparchus.optim.linear.Relationship
const SimplexSolver = @jimport org.hipparchus.optim.linear.SimplexSolver
const SimplexTableau = @jimport org.hipparchus.optim.linear.SimplexTableau
const SolutionCallback = @jimport org.hipparchus.optim.linear.SolutionCallback

export LinearConstraint
export LinearConstraintSet
export LinearObjectiveFunction
export LinearOptimizer
export NonNegativeConstraint
export PivotSelectionRule
export Relationship
export SimplexSolver
export SimplexTableau
export SolutionCallback
export compare_to
export compute_objective_value
export describe_constable
export do_optimize
export equals
export get_class
export get_coefficients
export get_constant_term
export get_constraints
export get_convergence_checker
export get_declaring_class
export get_evaluations
export get_goal_type
export get_iterations
export get_lower_bound
export get_max_evaluations
export get_max_iterations
export get_relationship
export get_solution
export get_start_point
export get_upper_bound
export get_value
export hash_code
export is_restricted_to_non_negative
export is_solution_optimal
export name
export normalize_constraints
export notify
export notify_all
export opposite_relationship
export optimize
export ordinal
export to_string
export value
export value_of
export values
export wait

include("linear_constraint.jl")
include("linear_constraint_set.jl")
include("linear_objective_function.jl")
include("linear_optimizer.jl")
include("non_negative_constraint.jl")
include("pivot_selection_rule.jl")
include("relationship.jl")
include("simplex_solver.jl")
include("simplex_tableau.jl")
include("solution_callback.jl")

end


text:
using Test

@setset "cut_mirror" begin
    , kG  = gen_kGrid(4,2)
    , kG2 = gen_kGrid(6,2)
    @test all(expand_mirror(cut_mirror(kG)) .== collect(kG))
    @test all(expand_mirror(cut_mirror(kG2)) .== collect(kG2))
end


text:
module LForaging

using Reexport
#add CSV@0.8.5
@reexport using DataFrames, CategoricalArrays, StatsPlots, BrowseTables, StatsBase, Bootstrap, MixedModels, Optim, Distributions
#CSV,
import Statistics: median, std

include("Analysis_fun.jl")
include("preprocess.jl")
include("RewRateUpdates.jl")

export survivalrate_algorythm, cumulative_algorythm, hazardrate_algorythm, function_analysis
export median, std
export preprocess_pokes, trial_info!, bout_info!, travel_info, leave_info!
export preprocess_bouts
export get_rew_rate

end # module


text:
include("view.jl")
using WaveSim

# Run the simulator and display results.
function main()
  # Define simulation parameters (use many default values, see WaveSimParameters).
  sim_params = WaveSimParameters(
    focus_depth = 0.03,  # [m]
    steer_angle = 10.0,  # [deg]
    aperture_size = 0.02,  # [m]
    temporal_res = 0.1e-6,  # [s]
    spatial_res = [256, 512]  # [pixels]
  )

  # Compute focusing delays for the elements of the phased array.
  trans_delays = WaveSim.delays_from_focus_and_steer(sim_params)

  # Run the simulation.
  images = WaveSim.wavesim(trans_delays, sim_params)
  beam_energy_map, transmit_time_map = WaveSim.beam_energy_map_and_transmit_time_map(images, sim_params)

  # Display results.
  imshowall(images, beam_energy_map, transmit_time_map, sim_params)

  return beam_energy_map, transmit_time_map, images
end

beam_energy_map, transmit_time_map, images = main()


text:
# TODO : these are not actually CSF but proximal to ventricular structures
# - Angular bundle
# - Subventricular zone
#
# TODO - foramens (monro etc.)

const CerebrospinalFluid = AnatomicalStructure("cerebrospinal fluid")
const CSF = CerebrospinalFluid

const Ventricles = AnatomicalStructure("ventricles")

const LateralVentricle = BilateralStructure(" lateral ventricle")
const InferiorLateralVentricle = BilateralStructure("inferior lateral ventricle")
const InferiorHorn = InferiorLateralVentricle

const PosteriorLateralVentricle = BilateralStructure("posterior lateral ventricle")
const PosteriorHorn = PosteriorLateralVentricle

const AnteriorLateralVentricle = BilateralStructure("anterior lateral ventricle")
const AnteriorHorn = AnteriorLateralVentricle

const BodyOfLateralVentricle = AnatomicalStructure("body of lateral ventricle")

const ThirdVentricle = AnatomicalStructure("third ventricle")

const FourthVentricle = AnatomicalStructure("fourth ventricle")

# note: this is commonly called the fifth ventricle (it shouldn't be though)
const CavumSeptumPellucidum = AnatomicalStructure("cavum septum pellucidum")


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ a80cb638-1dde-11eb-08dd-593b9630a62e
using Pkg, DrWatson

# ╔═╡ a80cedd8-1dde-11eb-2ac5-ff913973a868
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots, Images
end

# ╔═╡ 3bc81896-1dde-11eb-22ff-3372c06964a9
md"## Listing 1.13"

# ╔═╡ 924b85f8-1ddf-11eb-025f-6f6f1ba174c0
Pkg.add("ImageIO")

# ╔═╡ a80d6ee8-1dde-11eb-3a36-c1d068d268f0
begin
	img = load(datadir("stars.png"))
	gImg = red.(img)*0.299 + green.(img)*0.587 + blue.(img)*0.114
	rows, cols = size(img)
end

# ╔═╡ a8187068-1dde-11eb-290c-919a177c2f04
"Highest intensity pixel: ", findmax(gImg)

# ╔═╡ a818e0ca-1dde-11eb-2d10-6be0b50f7358
function boxBlur(image,x,y,d)
    if x<=d || y<=d || x>=cols-d || y>=rows-d
        return image[x,y]
    else
        total = 0.0
        for xi = x-d:x+d
            for yi = y-d:y+d
                total += image[xi,yi]
            end
        end
        return total/((2d+1)^2)
    end
end

# ╔═╡ a821fab6-1dde-11eb-36fd-db97e9de0689
begin
	blurImg = [boxBlur(gImg,x,y,5) for x in 1:cols, y in 1:rows]

	yOriginal, xOriginal = argmax(gImg).I
	yBoxBlur, xBoxBlur   = argmax(blurImg).I

	heatmap(gImg, c=:Greys, yflip=true)
end

# ╔═╡ a8223dee-1dde-11eb-3345-630187ae9ff6
p1 = scatter!((xOriginal, yOriginal), ms=60, ma=0, msw=4, msc=:red)

# ╔═╡ a82b3c52-1dde-11eb-1e09-f3ad27bbd9ee
heatmap(blurImg, c=:Greys, yflip=true)

# ╔═╡ a82bc988-1dde-11eb-17a5-67a907d95a5e
p2 = scatter!((xBoxBlur, yBoxBlur), ms=60, ma=0, msw=4, msc=:red)

# ╔═╡ a831c28e-1dde-11eb-2f74-53e1521437e3
plot(p1, p2, size=(800, 400), ratio=:equal, xlims=(0,cols), ylims=(0,rows),
	colorbar_entry=false, border=:none, legend=:none)

# ╔═╡ 0349ba5c-1de2-11eb-359e-cfe15670884f
md"## End of listing 1.13"

# ╔═╡ Cell order:
# ╟─3bc81896-1dde-11eb-22ff-3372c06964a9
# ╠═a80cb638-1dde-11eb-08dd-593b9630a62e
# ╠═a80cedd8-1dde-11eb-2ac5-ff913973a868
# ╠═924b85f8-1ddf-11eb-025f-6f6f1ba174c0
# ╠═a80d6ee8-1dde-11eb-3a36-c1d068d268f0
# ╠═a8187068-1dde-11eb-290c-919a177c2f04
# ╠═a818e0ca-1dde-11eb-2d10-6be0b50f7358
# ╠═a821fab6-1dde-11eb-36fd-db97e9de0689
# ╠═a8223dee-1dde-11eb-3345-630187ae9ff6
# ╠═a82b3c52-1dde-11eb-1e09-f3ad27bbd9ee
# ╠═a82bc988-1dde-11eb-17a5-67a907d95a5e
# ╠═a831c28e-1dde-11eb-2f74-53e1521437e3
# ╟─0349ba5c-1de2-11eb-359e-cfe15670884f


text:
module ContactImplicitMPC

using BenchmarkTools
using Colors
using CoordinateTransformations
using FFMPEG
using FileIO
using ForwardDiff
using GeometryBasics
using IfElse
using InteractiveUtils
using JLD2
using LinearAlgebra
using Logging
using QDLDL
using MeshCat
using MeshCatMechanisms
using Meshing
using MeshIO
using Parameters
using Plots
using Random
using Rotations
using SparseArrays
using StaticArrays
using Symbolics
using Test

# Utilities
include("utils.jl") #

# Solver
include("solver/gn.jl")
include("solver/lu.jl")
include("solver/ldl.jl")
include("solver/qr.jl")
include("solver/schur.jl")
include("solver/cones.jl") #

# Environment
include("simulator/environment.jl")

# Dynamics
include("dynamics/model.jl")

# Simulator
include("simulation/index.jl")

# Solver
include("solver/interior_point.jl")

# Simulator
include("simulation/contact_methods.jl")
include("simulation/simulation.jl")
include("simulator/trajectory.jl")

include("dynamics/code_gen_dynamics.jl")
include("dynamics/fast_methods_dynamics.jl")

# Models
include("dynamics/quaternions.jl")
include("dynamics/mrp.jl")
include("dynamics/euler.jl")

include("dynamics/particle_2D/model.jl")
include("dynamics/particle/model.jl")
include("dynamics/hopper_2D/model.jl")
include("dynamics/hopper_3D/model.jl")
include("dynamics/quadruped/model.jl")
include("dynamics/flamingo/model.jl")
include("dynamics/pushbot/model.jl")
include("dynamics/rigidbody/model.jl")

# Simulator
include("simulator/policy.jl")
include("simulator/disturbances.jl")
include("simulator/simulator.jl")

# Simulation
include("simulation/environments/flat.jl")
include("simulation/environments/piecewise.jl")
include("simulation/environments/quadratic.jl")
include("simulation/environments/slope.jl")
include("simulation/environments/sinusoidal.jl")
include("simulation/environments/stairs.jl")

include("simulation/residual_approx.jl")
include("simulation/code_gen_simulation.jl")

# Controller
include("controller/linearized_step.jl")
include("controller/implicit_dynamics.jl")
include("controller/objective.jl")
include("controller/linearized_solver.jl")
include("controller/newton.jl")
include("controller/newton_indices.jl")
include("controller/newton_residual.jl")
include("controller/newton_jacobian.jl")
include("controller/mpc_utils.jl")
include("controller/policy.jl")
include("controller/newton_structure_solver/methods.jl")

# Visuals
include("dynamics/visuals.jl")
include("dynamics/visual_utils.jl")

export
    World,
    LinearizedCone,
    NonlinearCone,
    ContactModel,
    Dimensions,
    BaseMethods,
    DynamicsMethods,
    ResidualMethods,
    Environment,
    environment_2D,
    environment_3D,
    environment_2D_flat,
    environment_3D_flat,
    get_model,
    SparseStructure,
    LinearizedStep,
    bil_addition!,
    r_linearized!,
    rz_linearized!,
    ImplicitTraj,
    linearization!,
    implicit_dynamics!,
    TrackingObjective,
    TrackingVelocityObjective,
    second_order_cone_product,
    generate_base_expressions,
    RLin,
    RZLin,
    RθLin,
    ContactTraj,
    Simulation,
    num_var,
    num_data,
    get_simulation,
    get_trajectory,
    interior_point,
    InteriorPointOptions,
    interior_point_solve!,
    r!,
    rz!,
    rθ,
    generate_base_expressions,
    save_expressions,
    instantiate_base!,
    generate_dynamics_expressions,
    save_expressions,
    instantiate_dynamics!,
    environment_3D_flat,
    friction_dim,
    dim,
    sqrt_quat,
    cayley_map,
    L_multiply,
    R_multiply,
    R2,
    R3,
    FrictionCone,
    rotation,
    module_dir,
    open_loop_disturbances,
    disturbances,
    open_loop_policy,
    policy,
    linearized_mpc_policy,
    NewtonOptions,
    LinearizedMPCOptions,
    SimulatorOptions,
    simulator,
    simulate!,
    generate_residual_expressions,
    instantiate_residual!,
    ϕ_func,
    tracking_error,
    repeat_ref_traj,
    Schur,
    schur_factorize!,
    schur_solve!,
    LinearSolver,
    LUSolver,
    lu_solver,
    factorize!,
    linear_solve!,
    OptimizationSpace,
    OptimizationSpace13,
    index_q2,
    index_γ1,
    index_b1,
    index_ψ1,
    index_s1,
    index_η1,
    index_s2,
    index_q0,
    index_q1,
    index_u1,
    index_w1,
    index_μ,
    index_h,
    index_dyn,
    index_imp,
    index_mdp,
    index_fri,
    index_bimp,
    index_bmdp,
    index_bfri,
    linearization_var_index,
    linearization_term_index,
    index_ort,
    index_soc,
    num_var,
    num_data,
    num_bilinear,
    index_equr,
    index_ortr,
    index_socr

end # module


text:
module JML

# package code goes here
include("functions.jl")

end # module


text:
# https://github.com/JuliaPlots/UnicodePlots.jl

# don't warn on unsupported... there's just too many warnings!!
warn_on_unsupported_args(::UnicodePlotsBackend, plotattributes) = nothing

# ------------------------------------------------------------------------------------------
const _canvas_map = (
    ascii = UnicodePlots.AsciiCanvas,
    block = UnicodePlots.BlockCanvas,
    braille = UnicodePlots.BrailleCanvas,
    density = UnicodePlots.DensityCanvas,
    dot = UnicodePlots.DotCanvas,
    heatmap = UnicodePlots.HeatmapCanvas,
    lookup = UnicodePlots.LookupCanvas,
)

# do all the magic here... build it all at once,
# since we need to know about all the series at the very beginning
function unicodeplots_rebuild(plt::Plot{UnicodePlotsBackend})
    plt.o = UnicodePlots.Plot[]

    for sp in plt.subplots
        xaxis = sp[:xaxis]
        yaxis = sp[:yaxis]
        xlim = collect(axis_limits(sp, :x))
        ylim = collect(axis_limits(sp, :y))

        # we set x/y to have a single point,
        # since we need to create the plot with some data.
        # since this point is at the bottom left corner of the plot,
        # it shouldn't actually be shown
        x = Float64[xlim[1]]
        y = Float64[ylim[1]]

        # create a plot window with xlim/ylim set,
        # but the X/Y vectors are outside the bounds
        canvas_type = if (ct = _canvas_type[]) == :auto
            isijulia() ? :ascii : :braille
        else
            ct
        end

        kw = (
            compact = true,
            title = texmath2unicode(sp[:title]),
            xlabel = texmath2unicode(xaxis[:guide]),
            ylabel = texmath2unicode(yaxis[:guide]),
            xscale = xaxis[:scale],
            yscale = yaxis[:scale],
            border = isijulia() ? :ascii : :solid,
            xlim = xlim,
            ylim = ylim,
        )

        o = UnicodePlots.Plot(x, y, _canvas_map[canvas_type]; kw...)
        for series in series_list(sp)
            o = addUnicodeSeries!(sp, o, kw, series, sp[:legend_position] != :none)
        end

        for ann in sp[:annotations]
            x, y, val = locate_annotation(sp, ann...)
            o = UnicodePlots.annotate!(
                o,
                x,
                y,
                texmath2unicode(val.str);
                color = up_color(val.font.color),
                halign = val.font.halign,
                valign = val.font.valign,
            )
        end

        push!(plt.o, o)  # save the object
    end
end

up_color(col::UnicodePlots.UserColorType) = col
up_color(col::RGBA) =
    (c = convert(ARGB32, col); map(Int, (red(c).i, green(c).i, blue(c).i)))
up_color(col) = :auto

# add a single series
function addUnicodeSeries!(
    sp::Subplot{UnicodePlotsBackend},
    up::UnicodePlots.Plot,
    kw,
    series,
    addlegend::Bool,
)
    st = series[:seriestype]

    # get the series data and label
    x, y = if st == :straightline
        straightline_data(series)
    elseif st == :shape
        shape_data(series)
    else
        float(series[:x]), float(series[:y])
    end

    # special handling (src/interface)
    if st == :histogram2d
        kw[:xlim][:] .= kw[:ylim][:] .= 0
        return UnicodePlots.densityplot(x, y; kw...)
    elseif st == :heatmap
        rng = range(0, 1, length = length(UnicodePlots.COLOR_MAP_DATA[:viridis]))
        cmap = [(red(c), green(c), blue(c)) for c in get(get_colorgradient(series), rng)]
        return UnicodePlots.heatmap(
            series[:z].surf;
            zlabel = sp[:colorbar_title],
            colormap = cmap,
            kw...,
        )
    elseif st == :spy
        return UnicodePlots.spy(series[:z].surf; kw...)
    end

    series_kw = (;)

    # now use the ! functions to add to the plot
    if st in (:path, :straightline, :shape)
        func = UnicodePlots.lineplot!
    elseif st == :scatter || series[:markershape] != :none
        func = UnicodePlots.scatterplot!
        series_kw = (; marker = series[:markershape])
    else
        error("Series type $st not supported by UnicodePlots")
    end

    label = addlegend ? series[:label] : ""

    for (n, segment) in enumerate(series_segments(series, st; check = true))
        i, rng = segment.attr_index, segment.range
        lc = get_linecolor(series, i)
        up = func(
            up,
            x[rng],
            y[rng];
            color = up_color(lc),
            name = n == 1 ? label : "",
            series_kw...,
        )
    end

    for (xi, yi, str, fnt) in EachAnn(series[:series_annotations], x, y)
        up = UnicodePlots.annotate!(
            up,
            xi,
            yi,
            str;
            color = up_color(fnt.color),
            halign = fnt.halign,
            valign = fnt.valign,
        )
    end

    return up
end

# ------------------------------------------------------------------------------------------

# since this is such a hack, it's only callable using `png`...
# should error during normal `show`
function png(plt::Plot{UnicodePlotsBackend}, fn::AbstractString)
    fn = addExtension(fn, "png")

    @static if Sys.isapple()
        # make some whitespace and show the plot
        println("\n\n\n\n\n\n")
        gui(plt)
        sleep(0.5)
        # use osx screen capture when my terminal is maximized
        # and cursor starts at the bottom (I know, right?)
        run(`screencapture -R50,600,700,420 $fn`)
        return
    elseif Sys.islinux()
        run(`clear`)
        gui(plt)
        run(`import -window $(ENV["WINDOWID"]) $fn`)
        return
    end

    error(
        "Can only savepng on MacOS or Linux with UnicodePlots " *
        "(though even then I wouldn't do it)",
    )
end

# ------------------------------------------------------------------------------------------
Base.show(plt::Plot{UnicodePlotsBackend}) = show(stdout, plt)
Base.show(io::IO, plt::Plot{UnicodePlotsBackend}) = _show(io, MIME("text/plain"), plt)

function _show(io::IO, ::MIME"text/plain", plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    nr, nc = size(plt.layout)
    lines_colored = Array{Union{Nothing,Vector{String}}}(undef, nr, nc)
    lines_uncolored = copy(lines_colored)
    l_max = zeros(Int, nr)
    w_max = zeros(Int, nc)
    buf = IOBuffer()
    cbuf = IOContext(buf, :color => true)
    re_col = r"\x1B\[[0-9;]*[a-zA-Z]"
    sps = 0
    for r in 1:nr
        lmax = 0
        for c in 1:nc
            l = plt.layout[r, c]
            if l isa GridLayout && size(l) != (1, 1)
                @error "UnicodePlots: complex nested layout is currently unsupported !"
            else
                if get(l.attr, :blank, false)
                    lines_colored[r, c] = lines_uncolored[r, c] = nothing
                else
                    sp = plt.o[sps += 1]
                    show(cbuf, sp)
                    colored = String(take!(buf))
                    uncolored = replace(colored, re_col => "")
                    lines_colored[r, c] = lc = split(colored, "\n")
                    lines_uncolored[r, c] = lu = split(uncolored, "\n")
                    lmax = max(length(lc), lmax)
                    w_max[c] = max(maximum(length.(lu)), w_max[c])
                end
            end
        end
        l_max[r] = lmax
    end
    empty = String[' '^w for w in w_max]
    for r in 1:nr
        for n in 1:l_max[r]
            for c in 1:nc
                pre = c == 1 ? '\0' : ' '
                lc = lines_colored[r, c]
                if lc === nothing || length(lc) < n
                    print(io, pre, empty[c])
                else
                    lu = lines_uncolored[r, c]
                    print(io, pre, lc[n], ' '^(w_max[c] - length(lu[n])))
                end
            end
            println(io)
        end
        r < nr && println(io)
    end
    nothing
end

function _display(plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    map(display, plt.o)
    nothing
end


text:
include("euler/euler.jl")
using .Calculus: is_pentagonal

function compute()::Int
    pentagonals = [1]
    j = 1
    for i ∈ pentagonals
        for k ∈ pentagonals
            if is_pentagonal(i - k) && is_pentagonal(abs(i - 2k))
                return abs(i - 2k)
            end
        end
        j += 3
        append!(pentagonals, i + j)
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Aqua,
      DataFrames,
      ModelParameters,
      Setfield,
      StaticArrays,
      Test,
      Unitful

import BenchmarkTools

@testset "param setproperties" begin
    param = Param(1; a=2.0, b="3", c='4')
    @set! param.val = 2
    @test param.val == 2
    @set! param.a = "99"
    @test param.a == "99"
end

@testset "param math" begin
    # We don't have to test everything, that is for AbstractNumbers.jl
    @test 2 * Param(5.0; bounds=(5.0, 15.0)) == 10.0
    @test Param(5.0; bounds=(5.0, 15.0)) + 3 == 8.0
    @test Param(5.0; bounds=(5.0, 15.0))^2 === 25.0
    @test Param(5; bounds=(5.0, 15.0))^2 === 25
end

struct S1{A,B,C,D,E,F}
    a::A
    b::B
    c::C
    d::D
    e::E
    f::F
end

struct S2{H,I,J}
    h::H
    i::I
    j::J
end

s2 = S2(
    Param(99),
    7,
    Param(100.0; bounds=(50.0, 150.0))
)

s1 = S1(
   Param(1.0; bounds=(5.0, 15.0)),
   Param(2.0; bounds=(5.0, 15.0)),
   Param(3.0; bounds=(5.0, 15.0)),
   Param(4.0),
   (Param(5.0; bounds=(5.0, 15.0)), Param(6.0; bounds=(5.0, 15.0))),
   s2,
)

pars = ModelParameters.params(s1)

@testset "params are correctly flattened from an object" begin
    @test length(pars) == 8
    @test all(map(p -> isa(p, Param), pars))
end

@testset "missing fields are added to Model Params" begin
    m = Model(s1);
    @test all(map(p -> propertynames(p) == (:val, :bounds), params(m)))
end

@testset "getproperties returns column tuples of param fields" begin
    m = Model(s1)
    @test m[:component] === (S1, S1, S1, S1, Tuple, Tuple, S2, S2)
    @test m[:fieldname] === (:a, :b, :c, :d, 1, 2, :h, :j)
    @test m[:val] === (1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 99, 100.0)
    @test m[:bounds] === ((5.0, 15.0), (5.0, 15.0), (5.0, 15.0), nothing,
                       (5.0, 15.0), (5.0, 15.0), nothing, (50.0, 150.0))
end

@testset "setindex updates and adds param fields" begin
    m = Model(s1)
    m[:val] = m[:val] .* 2
    @test m[:val] == (2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 198, 200.0)
    m[:newfield] = ntuple(x -> x, 8)
    @test m[:newfield] == ntuple(x -> x, 8)
end

@testset "show" begin
    m = Model(s1)
    sh = sprint(show, MIME"text/plain"(), m)
    @test occursin("Model with parent", sh)
    @test occursin("S1", sh)
    @test occursin("Param", sh)
    @test occursin("┌──────", sh)
    @test occursin("component", sh)
    @test occursin("100.0", sh)
end

@testset "strip params from model" begin
    m = Model(s1);
    stripped = stripparams(m)
    @test params(stripped) == ()
    @test stripped.c == 3.0
    @test stripped.f.j == 100.0
end

@testset "Tables interface" begin
    m = Model(s1);
    s = Tables.schema(m)
    @test keys(m) == s.names == (:component, :fieldname, :val, :bounds)
    @test s.types == (
        Union{DataType,UnionAll},
        Union{Int64,Symbol},
        Union{Float64,Int64},
        Union{Nothing,Tuple{Float64,Float64}},
    )
    @test Tables.rows(m) isa Tables.RowIterator
    @test Tables.columns(m) isa Model

    df = DataFrame(m)
    @test all(df.component .== m[:component])
    @test all(df.fieldname .== m[:fieldname])
    @test all(df.val .== m[:val])
    @test all(df.bounds .== m[:bounds])

    df.val .*= 3
    ModelParameters.update!(m, df)
    @test m[:val] == (3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 297, 300.0)
    df.val ./= 3
    newm = @inferred ModelParameters.update(m, df)
    @test newm[:val] == (1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 99.0, 100.0)
end

@testset "use Unitful units, with StaticModel" begin
    s1 = S1(
       Param(1.0; bounds=(5.0, 15.0)),
       Param(2.0; bounds=(5.0, 15.0), units=u"s"),
       Param(3.0; bounds=(5.0, 15.0), units=u"K"),
       Param(4.0; units=u"m"),
       Param(5.0; bounds=(5.0, 15.0)),
       Param(6.0),
    )
    s2 = S2( s1,
        Param(7.0; bounds=(50.0, 150.0), units=u"m*s^2"),
        Param(8.0),
    )
    m = StaticModel(s2)
    sh = sprint(show, m)
    @test occursin("bounds", sh)
    @test occursin("units", sh)
    @test m[:units] == (nothing, u"s", u"K", u"m", nothing, nothing, u"m*s^2", nothing)
    # Values have units now
    @test withunits(m) == (1.0, 2.0u"s", 3.0u"K", 4.0u"m", 5.0, 6.0, 7.0u"m*s^2", 8.0)
    @test withunits(m, :bounds) ==
        ((5.0, 15.0), (5.0, 15.0) .* u"s", (5.0, 15.0) .* u"K", nothing,
         (5.0, 15.0), nothing, (50.0, 150.0) .* u"m*s^2", nothing)
    @test stripunits(m, (1.0, 2.0u"s", 3.0u"K", 4.0u"m", 5.0, 6.0, 7.0u"m*s^2", 8.0)) ==
        (1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
    @test stripunits(m, ((5.0, 15.0), (5.0, 15.0) .* u"s", (5.0, 15.0) .* u"K", nothing,
                         (5.0, 15.0), nothing, (50.0, 150.0) .* u"m*s^2", nothing)) ==
                        ((5.0, 15.0), (5.0, 15.0), (5.0, 15.0), nothing,
                         (5.0, 15.0), nothing, (50.0, 150.0), nothing)
end

@testset "parameters in StaticArrays" begin
    s2 = S2(
        SA[Param(99)],
        7,
        SA[Param(100.0) Param(200.0)]
    ) |> Model
    @test params(s2) === (Param(99), Param(100.0), Param(200.0))
    s2[:val] = s2[:val] .+ 1.0
    @test params(s2) === (Param(100.0), Param(101.0), Param(201.0))
end

@testset "type stable update" begin
    s1 = S1(
       Param(1.0; bounds=(5.0, 15.0)),
       Param(2.0; bounds=(5.0, 15.0), units=u"s"),
       Param(3.0; bounds=(5.0, 15.0), units=u"K"),
       Param(4.0; units=u"m"),
       Param(5.0; bounds=(5.0, 15.0)),
       Param(6.0),
    )
    s2 = S2(
        s1,
        Param(7.0; bounds=(50.0, 150.0), units=u"m*s^2"),
        Param(8.0),
    )
    m = Model(s2)
    ps = collect(m[:val]).*2.0
    new_s2 = @inferred update(s2, ps)
    @test all(Model(new_s2)[:val] .== ps)
    b = BenchmarkTools.@benchmark update($s2, $ps)
    @test b.allocs == 0
end

@testset "parameter grouping" begin
    s1 = S1(
       Param(1.0; bounds=(5.0, 15.0), group=:A),
       Param(2.0; bounds=(5.0, 15.0), units=u"s", group=:A),
       Param(3.0; bounds=(5.0, 15.0), units=u"K", group=:A),
       Param(4.0; units=u"m", group=:B),
       Param(5.0; bounds=(5.0, 15.0), group=:B),
       Param(6.0, group=:B),
    )
    s2 = S2(
        s1,
        Param(7.0; bounds=(50.0, 150.0), units=u"m*s^2", group=:A),
        Param(8.0, group=:B),
    )
    m = Model(s2)
    # groupparams
    groupedparams = groupparams(m, :group)
    @test haskey(groupedparams, :A)
    @test length(groupedparams.A) == 4
    @test haskey(groupedparams, :B)
    @test length(groupedparams.B) == 4
    groupedparams = groupparams(m, :group, :fieldname)
    @test haskey(groupedparams, :A)
    @test haskey(groupedparams, :B)
    @test groupedparams.A.a == [s1.a]
    @test groupedparams.A.b == [s1.b]
    @test groupedparams.A.c == [s1.c]
    @test groupedparams.B.d == [s1.d]
    @test groupedparams.B.e == [s1.e]
    @test groupedparams.B.f == [s1.f]
    @test groupedparams.A.i == [s2.i]
    @test groupedparams.B.j == [s2.j]
    # flat
    groupedvals = mapflat(p -> p.val, groupedparams)
    @test groupedvals.A.a == [s1.a.val]
    @test groupedvals.A.b == [s1.b.val]
    @test groupedvals.A.c == [s1.c.val]
    @test groupedvals.B.d == [s1.d.val]
    @test groupedvals.B.e == [s1.e.val]
    @test groupedvals.B.f == [s1.f.val]
    @test groupedvals.A.i == [s2.i.val]
    @test groupedvals.B.j == [s2.j.val]
    # convert to tuples; uses maptype kwarg to recurse exclusively on NamedTuples, not arrays
    tuplegroups = mapflat(Tuple, groupedparams; maptype=NamedTuple)
    @test isa(tuplegroups.A.a, Tuple)
    @test isa(tuplegroups.A.b, Tuple)
    @test isa(tuplegroups.A.c, Tuple)
    @test isa(tuplegroups.B.d, Tuple)
    @test isa(tuplegroups.B.e, Tuple)
    @test isa(tuplegroups.B.f, Tuple)
    @test isa(tuplegroups.A.i, Tuple)
    @test isa(tuplegroups.B.j, Tuple)
end


text:
module DataFramesMeta

using Reexport

using MacroTools

@reexport using DataFrames

@reexport using Chain

# Basics:
export @with,
       @subset, @subset!,
       @orderby,
       @by, @combine,
       @transform, @select, @transform!, @select!,
       @eachrow, @eachrow!,
       @byrow,
       @based_on, @where # deprecated

include("parsing.jl")
include("macros.jl")
include("linqmacro.jl")
include("eachrow.jl")

end # module


text:
import Distributions
const Dists = Distributions
import MeasureTheory
const MT = MeasureTheory
using HypercubeTransform
using Test

@testset "TupleHC" begin
    a = Dists.Normal()
    b = Dists.Uniform()
    c = Dists.Product([Dists.Normal(), Dists.Normal()])
    d = Dists.Gamma()
    pos = [0.5, 0.5, 0.5, 0.5, 0.0]

    tc = ascube((a, b, c, d))
    nc = ascube((;a,b, c, d))
    @test dimension(tc) == length(pos)
    println(transform(tc, pos))
    @test transform(tc, pos) == values(transform(nc, pos))

    @inferred NTuple{3, Float64} transform(tc, pos)
end


text:
__precompile__()

module LibExpat
using Pkg.Artifacts
using Expat_jll

import Base: getindex, show, parse

include("lX_common_h.jl")
include("lX_defines_h.jl")
include("lX_expat_h.jl")
#include("lX_exports_h.jl")

export ETree, xp_parse, xpath, @xpath_str

# streaming
export XPCallbacks, XPStreamHandler,
       parse, stop, pause, resume, free, parsefile

const SymbolAny = Tuple{Symbol,Any}

mutable struct ETree
    # XML Tag
    name::AbstractString
    # Dict of tag attributes as name-value pairs
    attr::Dict{AbstractString,AbstractString}
    # List of child elements.
    elements::Vector{Union{ETree,AbstractString}}
    parent::ETree

    ETree() = ETree("")
    function ETree(name)
        pd=new(
            name,
            Dict{AbstractString, AbstractString}(),
            Union{ETree,AbstractString}[])
        pd.parent=pd
        pd
    end
end

Base.@deprecate_binding ParsedData ETree

function show(io::IO, pd::ETree)
    print(io,'<',pd.name)
    for (name,value) in pd.attr
        print(io,' ',name,'=','"',replace(value, '"' => """),'"')
    end
    if length(pd.elements) == 0
        print(io,'/','>')
    else
        print(io,'>')
        for ele in pd.elements
            if isa(ele, ETree)
                show(io, ele)
            else
                print(io, replace(ele, '<' => "<"))
            end
        end
        print(io,'<','/',pd.name,'>')
    end
end

string_value(pd::ETree) = String(take!(string_value(pd,IOBuffer())))
function string_value(pd::ETree, str::IOBuffer)
    for node in pd.elements
        if isa(node, AbstractString)
            write(str, node::AbstractString)
        elseif isa(node,ETree)
            string_value(node::ETree, str)
        end
    end
    str
end

mutable struct XPHandle
  parser::Union{XML_Parser,Nothing}
  pdata::ETree
  in_cdata::Bool

  XPHandle(p) = new(p, ETree(""), false)
end

function xp_geterror(p::Union{XML_Parser,Nothing})
    ec = XML_GetErrorCode(p)

    if ec != 0
        return ( unsafe_string(XML_ErrorString(XML_GetErrorCode(p))),
                XML_GetCurrentLineNumber(p),
                XML_GetCurrentColumnNumber(p) + 1,
                XML_GetCurrentByteIndex(p) + 1
            )
     else
        return ("", 0, 0, 0)
     end

end

function xp_geterror(xph::XPHandle)
    return xp_geterror(xph.parser)
end

function xp_close(xph::XPHandle)
  if (xph.parser != nothing)    XML_ParserFree(xph.parser) end
  xph.parser = nothing
end

function start_cdata(p_xph::Ptr{Nothing})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    xph.in_cdata = true
    return
end

function end_cdata(p_xph::Ptr{Nothing})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    xph.in_cdata = false
    return
end

function cdata(p_xph::Ptr{Nothing}, s::Ptr{UInt8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle

    txt = unsafe_string(s, Int(len))
    push!(xph.pdata.elements, txt)

    return
end

function comment(p_xph::Ptr{Nothing}, data::Ptr{UInt8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = unsafe_string(data)
    return
end

function default(p_xph::Ptr{Nothing}, data::Ptr{UInt8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = unsafe_string(data)
    return
end

function default_expand(p_xph::Ptr{Nothing}, data::Ptr{UInt8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = unsafe_string(data)
    return
end

function attrs_in_to_dict(attrs_in::Ptr{Ptr{UInt8}})
    attrs = Dict{AbstractString,AbstractString}()

    if (attrs_in != C_NULL)
        i = 1
        attr = unsafe_load(attrs_in, i)
        while (attr != C_NULL)
            k = unsafe_string(attr)

            i += 1
            attr = unsafe_load(attrs_in, i)

            attr == C_NULL && error("Attribute does not have a name!")
            v = unsafe_string(attr)

            attrs[k] = v

            i += 1
            attr = unsafe_load(attrs_in, i)
        end
    end

    return attrs
end

function start_element(p_xph::Ptr{Nothing}, name::Ptr{UInt8}, attrs_in::Ptr{Ptr{UInt8}})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    name = unsafe_string(name)

    new_elem = ETree(name)
    new_elem.parent = xph.pdata

    push!(xph.pdata.elements, new_elem)

    merge!(new_elem.attr, attrs_in_to_dict(attrs_in))
    xph.pdata = new_elem

    return
end

function end_element(p_xph::Ptr{Nothing}, name::Ptr{UInt8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    txt = unsafe_string(name)

    xph.pdata = xph.pdata.parent

    return
end

function start_namespace(p_xph::Ptr{Nothing}, prefix::Ptr{UInt8}, uri::Ptr{UInt8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    prefix = unsafe_string(prefix)
    uri = unsafe_string(uri)
    return
end

function end_namespace(p_xph::Ptr{Nothing}, prefix::Ptr{UInt8})
    xph = unsafe_pointer_to_objref(p_xph)::XPHandle
    prefix = unsafe_string(prefix)
    return
end

# Unsupported callbacks: External Entity, NotationDecl, Not Stand Alone, Processing, UnparsedEntityDecl, StartDocType
# SetBase and GetBase

function xp_parse(txt::AbstractString)
    xph = nothing
    xph = xp_make_parser()

    try
        rc = XML_Parse(xph.parser, txt, sizeof(txt), 1)
        (rc != XML_STATUS_OK) && error("Error parsing document : $rc")

        # The root element will only have a single child element in a well formed XML
        return xph.pdata.elements[1]
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(xph)
        rethrow("$e, $err, $line, $column, $pos")

    finally
        (xph != nothing) && xp_close(xph)
    end
end

function find(pd::ETree, path::T) where T<:AbstractString
    # What are we looking for?
    what = :node
    attr = ""

    pathext = split(path, "#")
    if (length(pathext)) > 2 error("Invalid path syntax")
    elseif (length(pathext) == 2)
        if (pathext[2] == "string")
            what = :string
        else
            error("Unknown extension : [$(pathext[2])]")
        end
    end

    xp = SymbolAny[]
    if path[1] == '/'
        # This will treat the incoming pd as the root of the tree
        push!(xp, (:root,:element))
        pathext[1] = pathext[1][2:end]
        #else - it will start searching the children....
    end

    nodes = split(pathext[1], "/")
    idx = false
    descendant = :child
    for n in nodes
        idx = false
        if length(n) == 0
            if descendant == :descendant
                error("too many / in a row")
            end
            descendant = :descendant
            continue
        end
        # Check to see if it is an index into an array has been requested, else default to 1
        m =  match(r"([\:\w]+)\s*(\[\s*(\d+)\s*\])?\s*(\{\s*(\w+)\s*\})?", n)

        if ((m == nothing) || (length(m.captures) != 5))
            error("Invalid name $n")
        else
            node = m.captures[1]
            push!(xp, (descendant,:element))
            descendant = :child
            push!(xp, (:name, convert(T,node)))

            if m.captures[5] != nothing
                if (n == nodes[end])
                    what = :attr
                end
                attr = m.captures[5]
                push!(xp, (:filter, (:attribute, convert(T,attr))))
            end

            if m.captures[3] != nothing
                push!(xp, (:filter, (:number, Base.parse(Int, m.captures[3]))))
                idx = true
            end
        end
    end

    pd = xpath(pd, XPath{T,Vector{ETree}}((:xpath, xp)))
    if what == :node
        if idx
            if length(pd) == 1
                return pd[1]
            else
                return nothing
            end
        else
            # If caller did not specify an index, return a list of leaf nodes.
            return pd
        end
    elseif length(pd) == 0
        return nothing

    elseif length(pd) != 1
        error("More than one instance of $pd, please specify an index")

    else
        pd = pd[1]
        if what == :string
            return string_value(pd)

        elseif what == :attr
            return pd.attr[attr]

        else
            error("Unknown request type")
        end
    end

    return nothing
end

include("xpath.jl")

include("streaming.jl")

function xp_make_parser(sep='\0')
    cb_start_cdata = @cfunction(start_cdata, Nothing, (Ptr{Nothing},))
    cb_end_cdata = @cfunction(end_cdata, Nothing, (Ptr{Nothing},))
    cb_cdata = @cfunction(cdata, Nothing, (Ptr{Nothing}, Ptr{UInt8}, Cint))
    cb_comment = @cfunction(comment, Nothing, (Ptr{Nothing}, Ptr{UInt8}))
    cb_default = @cfunction(default, Nothing,  (Ptr{Nothing}, Ptr{UInt8}, Cint))
    cb_default_expand = @cfunction(default_expand, Nothing, (Ptr{Nothing}, Ptr{UInt8}, Cint))
    cb_start_element = @cfunction(start_element, Nothing, (Ptr{Nothing}, Ptr{UInt8}, Ptr{Ptr{UInt8}}))
    cb_end_element = @cfunction(end_element, Nothing, (Ptr{Nothing}, Ptr{UInt8}))
    cb_start_namespace = @cfunction(start_namespace, Nothing, (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}))
    cb_end_namespace = @cfunction(end_namespace, Nothing, (Ptr{Nothing}, Ptr{UInt8}))

    p::XML_Parser = (sep == '\0') ? XML_ParserCreate(C_NULL) : XML_ParserCreateNS(C_NULL, sep);
    if (p == C_NULL) error("XML_ParserCreate failed") end

    xph = XPHandle(p)
    p_xph = pointer_from_objref(xph)
    XML_SetUserData(p, p_xph);

    XML_SetCdataSectionHandler(p, cb_start_cdata, cb_end_cdata)
    XML_SetCharacterDataHandler(p, cb_cdata)
    XML_SetCommentHandler(p, cb_comment)
    XML_SetDefaultHandler(p, cb_default)
    XML_SetDefaultHandlerExpand(p, cb_default_expand)
    XML_SetElementHandler(p, cb_start_element, cb_end_element)
#    XML_SetExternalEntityRefHandler(p, f_ExternaEntity)
    XML_SetNamespaceDeclHandler(p, cb_start_namespace, cb_end_namespace)
#    XML_SetNotationDeclHandler(p, f_NotationDecl)
#    XML_SetNotStandaloneHandler(p, f_NotStandalone)
#    XML_SetProcessingInstructionHandler(p, f_ProcessingInstruction)
#    XML_SetUnparsedEntityDeclHandler(p, f_UnparsedEntityDecl)
#    XML_SetStartDoctypeDeclHandler(p, f_StartDoctypeDecl)

    return xph
end

end


text:
module SciMLTutorials

using Weave, Pkg, IJulia, InteractiveUtils, Markdown

repo_directory = joinpath(@__DIR__,"..")
cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")
default_builds = (:script,:html,:github)

function weave_file(folder,file,build_list=default_builds)
  target = joinpath(repo_directory, "tutorials", folder, file)
  @info("Weaving $(target)")

  if isfile(joinpath(repo_directory, "tutorials", folder, "Project.toml"))
    @info("Instantiating", folder)
    Pkg.activate(joinpath(repo_directory,"tutorials", folder))
    Pkg.instantiate()
    Pkg.build()

    @info("Printing out `Pkg.status()`")
    Pkg.status()
  end

  args = Dict{Symbol,String}(:folder=>folder,:file=>file)
  if :script ∈ build_list
    println("Building Script")
    dir = joinpath(repo_directory,"script",basename(folder))
    mkpath(dir)
    tangle(target; out_path=dir)
  end
  if :html ∈ build_list
    println("Building HTML")
    dir = joinpath(repo_directory,"html",basename(folder))
    mkpath(dir)
    weave(target,doctype = "md2html",out_path=dir,args=args,css=cssfile,fig_ext=".svg")
  end
  if :pdf ∈ build_list
    println("Building PDF")
    dir = joinpath(repo_directory,"pdf",basename(folder))
    mkpath(dir)
    try
      weave(target,doctype="md2pdf",out_path=dir,template=latexfile,args=args)
    catch ex
      @warn "PDF generation failed" exception=(ex, catch_backtrace())
    end
  end
  if :github ∈ build_list
    println("Building Github Markdown")
    dir = joinpath(repo_directory,"markdown",basename(folder))
    mkpath(dir)
    weave(target,doctype = "github",out_path=dir,args=args)
  end
  if :notebook ∈ build_list
    println("Building Notebook")
    dir = joinpath(repo_directory,"notebook",basename(folder))
    mkpath(dir)
    Weave.convert_doc(target,joinpath(dir,file[1:end-4]*".ipynb"))
  end
end

function weave_all(build_list=default_builds)
  for folder in readdir(joinpath(repo_directory,"tutorials"))
    folder == "test.jmd" && continue
    weave_folder(folder,build_list)
  end
end

function weave_folder(folder,build_list=default_builds)
  for file in readdir(joinpath(repo_directory,"tutorials",folder))
    # Skip non-`.jmd` files
    if !endswith(file, ".jmd")
      continue
    end

    try
      weave_file(folder,file,build_list)
    catch e
      @error(e)
    end
  end
end

function tutorial_footer(folder=nothing, file=nothing)
    display(md"""
    ## Appendix

    These tutorials are a part of the SciMLTutorials.jl repository, found at: <https://github.com/SciML/SciMLTutorials.jl>.
    For more information on high-performance scientific machine learning, check out the SciML Open Source Software Organization <https://sciml.ai>.

    """)
    if folder !== nothing && file !== nothing
        display(Markdown.parse("""
        To locally run this tutorial, do the following commands:
        ```
        using SciMLTutorials
        SciMLTutorials.weave_file("$folder","$file")
        ```
        """))
    end
    display(md"Computer Information:")
    vinfo = sprint(InteractiveUtils.versioninfo)
    display(Markdown.parse("""
    ```
    $(vinfo)
    ```
    """))

    display(md"""
    Package Information:
    """)

    proj = sprint(io -> Pkg.status(io=io))
    mani = sprint(io -> Pkg.status(io=io, mode = Pkg.PKGMODE_MANIFEST))

    md = """
    ```
    $(chomp(proj))
    ```

    And the full manifest:

    ```
    $(chomp(mani))
    ```
    """
    display(Markdown.parse(md))
end

function open_notebooks()
  Base.eval(Main, Meta.parse("import IJulia"))
  weave_all((:notebook,))
  path = joinpath(repo_directory,"notebook")
  IJulia.notebook(;dir=path)
end

end


text:
module Plugins

using ..Types: uuid

abstract type AbstractPlugin end

using Base: UUID

struct Notary
    Signer::Function # allows to create signatures when necessery. id is obtained with verify on random data.
    Signature::Function
    verify::Function # verify returns id of the signature.
    hash::Function
end

struct Cypher
    G # a group for DiffieHellman key exchange
    rng::Function
    secureio::Function
end

struct CypherSuite{T} end ### This type belongs to plugins
CypherSuite(uuid::UUID) = CypherSuite{uuid.value}
CypherSuite(m::Module) = CypherSuite(uuid(m))

Notary(cyphersuite::UUID,crypto::Symbol) = Notary(CypherSuite(cyphersuite),crypto)::Notary
Cypher(cyphersuite::UUID,crypto::Symbol) = Cypher(CypherSuite(cyphersuite),crypto)::Cypher

struct Plugin{T} end

import Base.getindex

getindex(::Type{Plugin},uuid::UUID) = Plugin{uuid.value}
getindex(::Type{Plugin},x::Module) = Plugin[uuid(x)]

export Notary, Cypher, CypherSuite, Plugin

end


text:
"""
Data structure representing the estimated state of a state-space model
after filtering/smoothing.

#### Fields
- observations : Array of observations. Each column is  a single observation vector.
- state : Array of (possibly multivariate) Distributions. Each distribution
represents the estimate and uncertainty of the hidden state at that time.
- loglik : The log-likelihood of the model, i.e. the probabilty of the observations
given the model and state estimates.
- smoothed : Bool. Whether the state has been smoothed, or just filtered.
"""
type FilteredState{T, D<:ContinuousMultivariateDistribution}
	observations::Array{T, 2}
	state::Array{D}
	input::Array{T, 2}
	times::Vector{T}
	loglik::T
	smoothed::Bool
end

function show{T}(io::IO, fs::FilteredState{T})
	n = length(fs.state)
	dobs = size(fs.observations, 1)
	dstate = length(fs.state[1])
	s = ifelse(fs.smoothed, "smoothed", "filtered")
	println("FilteredState{$T}")
	println("$n $s estimates of $dstate-D process from $dobs-D observations")
	println("Log-likelihood: $(fs.loglik)")
end

"""
Returns the log-likelihood of the FilteredState object.
"""
loglikelihood(fs::FilteredState) = fs.loglik

for op in (:mean, :var, :cov, :cor, :rand)
	@eval begin
		function ($op){T}(s::FilteredState{T})
			result = Array(T, length(s.state[1]), length(s.state))
			for i in 1:length(s.state)
				result[:, i] = ($op)(s.state[i])
			end
			return result
		end
	end
end


text:
function findtracks(A,t,p,l);
    #FINDTRACKS creates track associative array.

        # Find docs that have person
        DocIDwPer = getrow(A[:,p])

        # Find docs that have person and location.
        DocIDwPerLoc = getrow(A[DocIDwPer,l])

        # Find docs that have person, location and time.
        DocIDwPerLocTime = getrow(A[DocIDwPerLoc,t])

        # Limit to these documents.
        AA = A[DocIDwPerLocTime,:]

        # Get person sub array.
        Aper = AA[DocIDwPerLocTime,p]
        TrackPer,DocAper = find(Aper')

        # Get location sub array.
        Aloc = AA[DocIDwPerLocTime,l]
        EntAloc,DocAloc = find(Aloc')

        # Get Single location per document- order by actors
        uLocDocs = unique(DocAloc)
        uLocDocIdx = indexin(uLocDocs,DocAloc) # get first index of unique docs (correct for getting highest)
        uDocLocs = EntAloc[uLocDocIdx] # single locations per document
        uLocDocIdxinDocAper = indexin(DocAper,uLocDocs) # locations of unique loc docs in per doc list
        TrackLoc = uDocLocs[uLocDocIdxinDocAper]

        # Get time sub array.
        Atime = AA[DocIDwPerLocTime,t]
        EntAtime,DocAtime = find(Atime')

        # Get Single time per document- order by actors
        uTimeDocs = unique(DocAtime)
        uTimeDocIdx = indexin(uTimeDocs,DocAtime) # get first index of unique docs (correct for getting highest)
        uDocTimes = EntAtime[uTimeDocIdx] # single times per document
        uTimeDocIdxinDocAper = indexin(DocAper,uTimeDocs) # locations of unique time docs in per doc list
        TrackTime = uDocTimes[uTimeDocIdxinDocAper]

        Tr = Assoc(TrackTime,TrackPer,TrackLoc);

    end


text:
h5p_set_alignment(plist_id, threshold, alignment) =
    ccall((:H5Pset_alignment, libhdf5), Herr, (Hid, Hsize, Hsize), plist_id, threshold, alignment)

h5p_get_alignment(plist_id, threshold, alignment) =
    ccall((:H5Pget_alignment, libhdf5), Herr, (Hid, Ptr{Hsize}, Ptr{Hsize}), plist_id, threshold, alignment)

function get_alignment(plist_id)
    threshold = Ref{Hsize}()
    alignment = Ref{Hsize}()
    h5p_get_alignment(plist_id, threshold, alignment)
    return threshold[], alignment[]
end

@init HDF5.hdf5_prop_get_set["alignment"] = (nothing, h5p_set_alignment, H5P_FILE_ACCESS)


text:
const TESTDIR = dirname(@__FILE__)

include(joinpath(TESTDIR, "ExampleTests.jl"))


text:
using Revise
using Soss

m = @model x begin
    α ~ Cauchy()
    β ~ Normal()
    σ ~ HalfNormal()
    y ~ For(x) do xj
        Normal(α + β * xj, σ)
    end
end;

truth = rand(m(x=randn(10)));
using Plots
xx = range(extrema(truth.x)...,length=100)
plot(xx, truth.α .+ truth.β .* xx, legend=false)
scatter!(truth.x,truth.y)

post = dynamicHMC(m(X=truth.X), (y=truth.y,));

particles(post)

pred = predictive(m,:β)

function fit(m::Model, verbosity::Integer, X, y)
    fitresult = dynamicHMC(m(X=X), (y=y,))
    cache = nothing
    report = nothing
    (fitresult, cache, report)
end

struct MLJpredictor <: Distributions.Sampleable{Univariate, Continuous}
    pred :: Soss.Model
    Xrow :: Matrix
    βs   :: Vector{NamedTuple}
end

function Base.rand(p::MLJpredictor)
    args = merge(rand(p.βs), (X=p.Xrow,))
    rand(p.pred(args)).y
end

function predict(m, fitresult, Xnew)
    pred = predictive(m, setdiff(variables(m),[:X,:y])...)
    map(eachrow(Xnew)) do x
        X = reshape(x, 1, :)
        MLJpredictor(pred, X, fitresult)
    end
end

fitresult = fit(m, 0, truth.X, truth.y)[1];
p = predict(m,fitresult, rand(4,3));

map(rand, p)

# mutable struct MLJmodel <: MLJ.Probabilistic
#     sampler #function that can do MCMC
#     m :: Soss.Model
# end

# MLJBase.fit(model::SomeSupervisedModel, verbosity::Integer, X, y) -> fitresult, cache, report
# MLJBase.predict(model::SomeSupervisedModel, fitresult, Xnew) -> yhat

# function fit(mljModel, verbosity, X, y)
#     # This line is wrong, we'll need to grab the right NamedTuple keys
#     (samples, report) = mljModel.sampler(mljModel.m(X=X), (y=y,))
#     cache = ...
#     (samples, cache, report)
# end

# function Base.rand(p::MLJpredictor)
#     β = rand(p.βs)

#     # This line is wrong, we'll need to grab the right NamedTuple keys
#     rand(p.pred(X=p.Xrow, β=β))
# end

# function predict(mljModel,βs,Xnew)
#     pred = predictive(mljModel.m, mljModel.X)
#     [MLJpredictor(pred, Xrow, βs) for Xrow in eachrow(Xnew)]
# end


text:
export MisfitParam,getMisfitParam

const ArrayUnionLocRemote = Union{Array{RemoteChannel},Array{Future},Array{Float64}}
const UnionLocRemote      = Union{RemoteChannel,Future,AbstractFloat}

"""
mutable struct jInv.InverseSolve.MisfitParam

Type storing information about one term in the misfit

F(m) = sum_i^n phi_i(pFor(model(m)),dobs,Wd)

Fields:

    pFor::ForwardProbType  - forward problem
    Wd                     - inverse standard deviation
    dobs                   - observed data
    misfit::Function       - misfit function
    modelfun::Function     - model function (evaluated locally)
    gloc                   - mapping from inverse to forward mesh (including active cells projection).

Constructors:

getMisfitParam(pFor,Wd,dobs,misfit,model,gloc=getGlobalToLocal(1.0))

getMisfitParam(pForRFs::Array{RemoteChannel}, Wd::Array, dobs::Array, misfit::Function,
                            Iact,sigmaBack::Vector,
                            Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}}=ones(length(pForRFs)),
                            modelfun::Function=identityMod,fname="")

"""
mutable struct  MisfitParam
    pFor::ForwardProbType
    Wd
    dobs
    misfit::Function
    modelfun::Function
    gloc
end

"""
function jInv.InverseSolve.getMisfitParam

Required Input:

Input for a call:
    getMisfitParam(pFor,Wd,dobs,misfit,model,gloc=getGlobalToLocal(1.0))

    pFor::ForwardProbType  - forward problem
    Wd                     - inverse standard deviation of data
    dobs                   - observed data
    misfit::Function       - Misfit function
    modelfun::Function     - model function. If all misfits have same model function,
                             put modelfun in inverseParam and identityMod here for efficiency.
    gloc                   - mapping from inverse to forward mesh (default: identity). See globalToLocal.jl

Input for a call:
    getMisfitParam(pForRFs::Array{RemoteChannel}, Wd::Array, dobs::Array, misfit::Function,
                            Iact,sigmaBack::Vector, Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}}=ones(length(pForRFs)),
                            modelfun::Function=identityMod,fname="")

    pForRFs::Array{RemoteChannel}                      - Array of remote references for forward problem parameters on each worker.
    Wd::Array                                - Array of inverse standard deviation for the data on each worker.
    dobs::Array                                - Array of observed data for each worker.
    misfit::Function                            - Misfit function
    Iact                                        - Projector to active cells.
    sigmaBack::Vector                            - Background model ("frozen" cells).
    Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}}    - mapping from inverse to forward mesh (default: identity)
    modelfun::Function=identityMod                         - model function. If all misfits have same model function,
                                          put modelfun in inverseParam and identityMod here for efficiency.
    fname=""                                 - (optional)

"""
function getMisfitParam(pFor::ForwardProbType, Wd, dobs, misfit::Function, modelfun::Function=identityMod, gloc=getGlobalToLocal(1.0))
    return MisfitParam(pFor,Wd,dobs,misfit,modelfun,gloc);
end

function getMisfitParam(pForRFs::Array{RemoteChannel}, Wd::Array, dobs::Array, misfit::Function,
                            Iact,sigmaBack::Union{Vector,AbstractFloat,AbstractModel},
                            Mesh2MeshRFs::ArrayUnionLocRemote=ones(length(pForRFs)),
                            modelfun::Function=identityMod,fname="")
    pMis                = Array{RemoteChannel}(undef, length(pForRFs));
    # the next references are for not to broadcast Iact and sigmaBackground
    IactRFs            = Array{Future}(undef,maximum(workers()));
    sigmaBackgroundRFs = Array{Future}(undef,maximum(workers()));

    @sync begin
        for p=workers()
            @async begin
                # communicate Iact and sigmaBack
                IactRFs[p]            = remotecall_wait(identity,p,Iact)            # send active cells projection matrix
                sigmaBackgroundRFs[p] = remotecall_wait(identity,p,sigmaBack) # send background conductivity
                for idx=1:length(pForRFs)
                    if p==pForRFs[idx].where
                        pMis[idx] = initRemoteChannel(getMisfitParam,p,pForRFs[idx],Wd[idx],dobs[idx],misfit,IactRFs[p],sigmaBackgroundRFs[p],Mesh2MeshRFs[idx],modelfun,fname);
                    end
                end
            end
        end
    end
    return pMis;
end

function getMisfitParam(pForRF::RemoteChannel, Wd, dobs, misfit::Function,
                            Iact::Future,sigmaBack::Future,
                            Mesh2MeshRF::UnionLocRemote = 1.0,
                            modelfun::Function=identityMod,fname="")
    # Single worker version of getMisfitParam.
    worker = pForRF.where;
    if !isa(Mesh2MeshRF,AbstractFloat)
        Mesh2Mesh     = fetch(Mesh2MeshRF);
        finalize(Mesh2MeshRF)  # to prevent memory leak
        if worker!=Mesh2MeshRF.where
            error("getMisfitParam::Mesh2Mesh and pFor not on the same worker");
        end
    else
        Mesh2Mesh = Mesh2MeshRF;
    end
    pFor         = fetch(pForRF);
    Iact         = fetch(Iact);
    sigmaBack     = fetch(sigmaBack);
    Wd              = fetch(Wd);
    dobs            = fetch(dobs);
    pMis         = getMisfitParam(pFor,Wd,dobs,misfit,modelfun, prepareGlobalToLocal(Mesh2Mesh,Iact,sigmaBack,fname));
    return pMis;
end

import jInv.Utils.clear!
function clear!(P::MisfitParam;clearPFor::Bool=true, clearData::Bool=false,clearMesh2Mesh::Bool=false)
    if clearPFor
        clear!(P.pFor);
    end
    if clearData
        P.Wd   = 0;
        P.dobs = 0;
    end
    if clearMesh2Mesh
        P.gloc = getGlobalToLocal(1.0);
    end
    return P;
end


text:
module Space

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

export space

Genie.Renderer.Html.register_normal_element("q__separator", context = @__MODULE__)

function space(args...; wrap::Function = StippleUI.DEFAULT_WRAPPER, kwargs...)
  wrap() do
    q__separator(args...; kwargs...)
  end
end

end