text:
get_s3_bucket_arn(cluster_name) = get_cluster(cluster_name).s3_bucket_arn
get_s3_bucket_path(cluster_name) =
    replace(get_cluster_s3_bucket_arn(), "arn:aws:s3:::" => "s3://")
function get_s3fs_bucket_path(cluster_name)
    arn = get_cluster_s3_bucket_arn()
    joinpath(
        homedir(),
        ".banyan",
        "mnt",
        "s3",
        arn[findfirst("arn:aws:s3:::", arn).stop+1:end],
    )
end

#########################
# MOUNTED S3 FILESYSTEM #
#########################

tmp_paths = []

# You should first call `download_remote_path` to get a remote path, then
# call `with_downloaded_path_for_reading` for the local path to read from.

function download_remote_path(remotepath)
    global tmp_paths

    # This returns an AbstractPath that can be manipulated using any of the
    # `Base.Filesystem` methods. To actually use it for reading, you should
    # first pass the returned `AbstractPath` through
    # `with_downloaded_path_for_reading`. That will ensure that if the path
    # is still referring to an S3 object, it will be first downloaded locally.
    # All files will be deleted afterwards if they are in a temporary
    # directory.

    if startswith(remotepath, "s3://")
        download_remote_s3_path(remotepath)
        # This will either return an `S3Path` still referring to a remote
        # location or it will return a local S3FS path for reading from.
    elseif startswith(remotepath, "http://") || startswith(remotepath, "https://")
        tmp_path = tempname() * splitext(remotepath)[2]
        push!(tmp_paths, tmp_path)
        Downloads.download(remotepath, tmp_path)
    else
        throw(
            ArgumentError(
                "$remotepath does not start with either http:// or https:// or s3://",
            ),
        )
    end
end

failed_to_use_s3fs = false

function download_remote_s3_path(path)
    global failed_to_use_s3fs

    # Get information about requested object
    s3path = S3Path(path, config = get_aws_config())
    bucket = s3path.bucket
    key = s3path.key
    # bucket = "banyan-cluster-data-myfirstcluster"
    mount = joinpath(homedir(), ".banyan", "mnt", "s3", bucket)

    # If specified, do not allow even attempting to use S3FS.
    if haskey(ENV, "BANYAN_USE_S3FS") && ENV["BANYAN_USE_S3FS"] == "0"
        failed_to_use_s3fs = true
    end

    # Don't attempt to use S3FS if we have ever failed to use S3FS.
    if !failed_to_use_s3fs
        # Ensure path to mount exists
        no_mount = false
        try
            if !isdir(mount)
                mkpath(mount)
                # TODO: Ensure that no directory really means there is no mount
                no_mount = true
                @warn "Attempting to remount S3FS because no directory found at $mount"
            end
            if !ismount(mount)
                no_mount = true
                @warn "Attempting to remount S3FS because no mount found at $mount"
            end
        catch
            no_mount = true
            @warn "Attempting to remount S3FS because attempting to stat the directory at $mount failed"
        end

        # Ensure something is mounted
        if no_mount
            try
                run(`umount -fq $mount`)
            catch e
                @warn "Failed to re-mount S3FS with error: $e. You may try to force unmounting with \`umount -fq $mount\` and then re-run."
            end

            # TODO: Store buckets from different accounts/IAMs/etc. seperately
            try
                ACCESS_KEY_ID = get_aws_config()[:creds].access_key_id
                SECRET_ACCESS_KEY = get_aws_config()[:creds].secret_key
                passwd_s3fs_contents = ACCESS_KEY_ID * ":" * SECRET_ACCESS_KEY
                HOME = homedir()
                region = get_aws_config_region()
                run(pipeline(`echo $passwd_s3fs_contents`, "$HOME/.passwd-s3fs"))
                run(`chmod 600 $HOME/.passwd-s3fs`)
                run(`s3fs $bucket $mount -o url=https://s3.$region.amazonaws.com -o endpoint=$region -o passwd_file=$HOME/.passwd-s3fs`)
            catch e
                @warn "Failed to mount S3 bucket \"$bucket\" at $mount using s3fs with error: $e. You may ensure s3fs is in PATH or mount manually."
                failed_to_use_s3fs = true
            end
        end
    end

    if failed_to_use_s3fs
        # If there is still no mount, return an S3Path
        s3path
    else
        # Return local path to object
        joinpath(mount, key)
    end
end

function with_downloaded_path_for_reading(func::Function, downloaded_path; for_writing=false)
    global tmp_paths

    # There are 3 cases here: `downloaded_path` is an S3Path or a local S3FS
    # path or an http:// file that has been downloaded to tempdir() (i.e., /tmp).

    temp_downloaded_path = if downloaded_path isa S3Path
        temp_downloaded_path = Path(tempname() * splitext(downloaded_path)[2])
        push!(tmp_paths, temp_downloaded_path)
        if !for_writing
            cp(downloaded_path, temp_downloaded_path, force=true)
        end
        func(string(temp_downloaded_path))
        if for_writing
            cp(temp_downloaded_path, downloaded_path, force=true)
        end
        string(temp_downloaded_path)
    else
        func(string(downloaded_path))
        string(downloaded_path)
    end
end

function cleanup_tmp()
    global tmp_paths
    for tmp_path in tmp_paths
        try
            rm(tmp_path, recursive=true, force=true)
        catch e
            @warn "Unable to delete $tmp_path: $e"
        end
    end
    empty!(tmp_paths)
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
module ForcedFlowFreeSlip

using Printf

using Oceananigans
using Oceananigans.BoundaryConditions
using Oceananigans.OutputWriters
using Oceananigans.Fields

# Functions that define the forced flow problem

 ξ(t) = 1 + sin(t^2)
ξ′(t) = 2t * cos(t^2)

 f(x, t) =   cos(x - ξ(t))
fₓ(x, t) = - sin(x - ξ(t))

Fᵘ(x, y, t) = (4 * f(x, t) - 2 * ξ′(t) * fₓ(x, t)) * cos(y)
Fᵛ(x, y, t) = sin(2y) / 2

u(x, y, t) =   f(x, t) * cos(y)
v(x, y, t) = -fₓ(x, t) * sin(y)

#####
##### x, z
#####

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time
    progress = 100 * (i / simulation.stop_iteration)
    @info @sprintf("[%05.2f%%] iteration: %d, time: %.5e", progress, i, t)
    return nothing
end

const DATA_DIR = joinpath(@__DIR__, "..", "data")

function setup_xz_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir=DATA_DIR)

    grid = RegularRectilinearGrid(size=(Nx, 1, Nx), x=(0, 2π), y=(0, 1), z=(0, π),
                                topology=(Periodic, Periodic, Bounded))

    model = IncompressibleModel(architecture = architecture,
                                        grid = grid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = nothing,
                                     closure = IsotropicDiffusivity(ν=1),
                                     forcing = (u = (x, y, z, t) -> Fᵘ(x, z, t),
                                                w = (x, y, z, t) -> Fᵛ(x, z, t))
                                )

    set!(model, u = (x, y, z) -> u(x, z, 0),
                w = (x, y, z) -> v(x, z, 0))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress=print_progress, iteration_interval=40)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, model.velocities;
                                                          dir = dir, force = true,
                                                          field_slicer = nothing,
                                                          prefix = @sprintf("forced_free_slip_xz_Nx%d_Δt%.1e", Nx, Δt),
                                                          schedule = TimeInterval(stop_iteration * Δt / 2))

    return simulation
end

function setup_and_run_xz(; setup...)
    simulation = setup_xz_simulation(; setup...)
    @info "Running forced flow free slip simulation with Nx = Nz = $(setup[:Nx]), Δt = $(setup[:Δt])"
    @time run!(simulation)
    return nothing
end

#####
##### x, y
#####

function setup_xy_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir=DATA_DIR)

    grid = RegularRectilinearGrid(size=(Nx, Nx, 1), x=(0, 2π), y=(0, π), z=(0, 1),
                                topology=(Periodic, Bounded, Bounded))

    model = IncompressibleModel(architecture = architecture,
                                        grid = grid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = nothing,
                                     closure = IsotropicDiffusivity(ν=1),
                                     forcing = (u = (x, y, z, t) -> Fᵘ(x, y, t),
                                                v = (x, y, z, t) -> Fᵛ(x, y, t))
                                )

    set!(model, u = (x, y, z) -> u(x, y, 0),
                v = (x, y, z) -> v(x, y, 0))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress=print_progress, iteration_interval=40)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, model.velocities;
                                                          dir = dir, force = true,
                                                          field_slicer = nothing,
                                                          prefix = @sprintf("forced_free_slip_xy_Nx%d_Δt%.1e", Nx, Δt),
                                                          schedule = TimeInterval(stop_iteration * Δt / 2))

    return simulation
end

function setup_and_run_xy(; setup...)
    simulation = setup_xy_simulation(; setup...)
    @info "Running forced flow free slip simulation with Nx = Ny = $(setup[:Nx]), Δt = $(setup[:Δt])"
    @time run!(simulation)
    return nothing
end

end # module

ext: jl


text:
using CompScienceMeshes
using BEAST

dir = dirname(pathof(BEAST))
Γ = readmesh(joinpath(dir,"../examples/sphere2.in"))

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

bels, bad = BEAST.assemblydata(X)
tels, tad = BEAST.assemblydata(Y)

Y = BEAST.RTBasis(Y.geo.mesh, Y.fns, Y.pos)

Δt, Nt = 0.6, 200
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
T3 = timebasisshiftedlagrange(Δt, Nt, 3)
δ = timebasisdelta(Δt, Nt)

DL = TDMaxwell3D.doublelayer(speedoflight=1.0)
SL = TDMaxwell3D.singlelayer(speedoflight=1.0, numdiffs=1)

# BEAST.allocatestorage(K,W,V,Val{:bandedstorage},BEAST.LongDelays{:ignore})

function as(op, W, V)

    X = spatialbasis(V)
    Y = spatialbasis(W)

    aux = BEAST.EmptyRP(1.0)

    M = numfunctions(Y)
    N = numfunctions(X)

    K0 = fill(typemax(Int), M, N)
    K1 = zeros(Int, M, N)

    function store(v,m,n,k)
        K0[m,n] = min(K0[m,n],k)
        K1[m,n] = max(K1[m,n],k)
    end

    BEAST.assemble_chunk!(aux, W, V, store)

end


text:
data = split(read(stdin, String), "\n\n", keepempty=false)

mutable struct Tile
    name
    l
    u
    r
    d
    rotation
    borders
end

new_tile(name, l, u, r, d; borders=[]) = Tile(name, l, u, r, d, 0, borders)

# Reverse rotation (ltr)
tiles = []
for tile in data
    tile = split(tile, '\n', keepempty=false)
    name = parse(Int, tile[1][6:end - 1])

    tile = [tile[i][j] for i = 1:length(tile), j = 1:length(tile[1])]
    l = tile[end:-1:2, 1]
    u = tile[2, :]
    r = tile[2:end, end]
    d = tile[end, end:-1:begin]

    push!(tiles, new_tile(name, l, u, r, d, borders=[l, u, r, d]))
end

# Find borders
count = 0
# ids[tileid] = #edge
ids = Dict()
for tile in tiles
    global count

    borders = [tile.l, tile.u, tile.r, tile.d]
    for (borderi, border) in enumerate(borders)
        match = false
        for tile2 in tiles
            if tile2 == tile
                continue
            end

            borders2 = [tile2.l, tile2.u, tile2.r, tile2.d]

            for border2 in borders2
                # Found matching border
                if border2 == border || border2 == border[end:-1:begin]
                    match = true
                    break
                end
            end

            if match
                break
            end
        end

        if !match
            count += 1

            if haskey(ids, tile.name)
                push!(ids[tile.name], borderi)
            else
                ids[tile.name] = [borderi]
            end
        end
    end
end

corners = [tile for (tile, v) in ids if length(v) == 2]

@assert length(corners) == 4

prod = reduce(*, corners)

println("Part 1 : $prod")

# TODO : Remove
# # Top left corner
# first = [t for t in tiles if t.name == corners[end]][begin]
# println("First : $(first.name)")

# current = first
# for tile in tiles
#     global current

#     if tile == current
#         continue
#     end

#     match = false
#     for (borderi, border) in enumerate(tile.borders)
#         if border == current.borders[3]
#             # Not flipped
#             println(tile.name)
#             match = true
#             break
#         elseif border == current.borders[3][end:-1:begin]
#             # Flipped
#             println("$(tile.name) flipped $borderi")
#             match = true
#             break
#         end
#     end

#     if match
#         break
#     end
# end

# TODO : Remove
# Old part 1
# visited = repeat([false], length(tiles))

# # map[i, j] = tile at i, j
# map = Dict()
# map[(0, 0)] = tiles[begin]
# visited[begin] = true

# println(tiles[begin])
# println(tiles[end])

# mini = -0xffff
# minj = -0xffff
# maxi = 0xffff
# maxj = 0xffff

# # Left
# # rot_id : left = 1, down = 2, ...
# function solve(i, j, ij_new, rot_id)
#     global map, visited, tiles
#     global mini, minj, maxi, maxj

#     if !(mini <= i <= maxi) || !(minj <= j <= maxj)
#         return false
#     end

#     current = map[(i, j)]
#     # if current == nothing
#     #     return false
#     # end

#     cedges = [current.l, current.u, current.r, current.d]
#     if !haskey(map, ij_new)
#         for (tilei, (visit, tile)) in enumerate(zip(visited, tiles))
#             if !visit
#                 matches = false
#                 newedges = [tile.l, tile.u, tile.r, tile.d]

#                 # Test each rotation possibility
#                 for rot = 0:3
#                     if newedges[1 + rot] == cedges[1 + (rot_id + current.rotation - 1) % 4][end:-1:begin] || newedges[1 + rot] == cedges[1 + (rot_id + current.rotation - 1) % 4]
#                         matches = true
#                         tile.rotation = rot
#                         break
#                     end
#                 end

#                 if matches
#                     visited[tilei] = true
#                     map[ij_new] = tile
#                     return true
#                 end
#             end
#         end
#     end

#     # map[ij_new] = nothing

#     if rot_id == 1
#         minj = max(minj, j)
#     elseif rot_id == 2
#         mini = max(mini, i)
#     elseif rot_id == 3
#         maxj = min(maxj, j)
#     elseif rot_id == 4
#         maxi = min(maxi, i)
#     end

#     return false
# end

# function solveall(i, j)
#     solve(i, j, (i, j - 1), 1) && solveall(i, j - 1)
#     solve(i, j, (i - 1, j), 2) && solveall(i - 1, j)
#     solve(i, j, (i, j + 1), 3) && solveall(i, j + 1)
#     solve(i, j, (i + 1, j), 4) && solveall(i + 1, j)
# end

# # function solvefrom(i, j)
# #     global map, visited, tiles

# #     solveall(i, j)
# # end

# # solve(0, 0, (-1, 0), 2)
# # solve(0, 0, (-1, 0), 2)
# # solvefrom(0, 0, (-1, 0), 2)
# solveall(0, 0)

# # last = (0, 0)
# # for i = 1:dim(tiles, 1)
# #     for j = 1:dim(tiles, 2)
# #         solveall(last)
# #     end
# # end

# @assert all(visited) "Error : All tiles not visited"

# println()
# println(map)

# # for (k, v) in map
# #     if v != nothing
# #         (i, j) = k

# #         global mini = min(mini, i)
# #         global minj = min(minj, j)
# #         global maxi = max(maxi, i)
# #         global maxj = max(maxj, j)
# #     end
# # end

# println(mini)
# println(minj)
# println(maxi)
# println(maxj)

# # println("> $(map[(1, 0)])")

# prod = map[(mini, minj)] * map[(maxi, minj)] * map[(maxi, maxj)] * map[(mini, maxj)]
# println("Part 1 : $prod")


text:
module Photopolarimetry

using Light, Images, ImageMagick, Colors, FixedPointNumbers

export path2stokes, getfnames, build_opt, getsz, stokes2mats, mats2images, stokes2profile, stokes2ellipse

const N = 1000
const N2 = round(Int,N/2)
const prop = Dict(:spatialorder => ["x","y"], :colorspace => "RGB", :pixelspacing => [1,1])

function getfnames(path::AbstractString)
	files = readdir(path)
	filter!(r"^[^\.].*\.(?:NEF|CR2)$",files)
	@assert length(files) == 4 || length(files) == 6
	return map(x -> joinpath(path,x), files)
end
getsz(fname) = map(x -> round(Int,parse(Int,x)/2),split(split(readall(`dcraw -i -v $fname`),'\n')[14])[[3,5]])
function build_opt(SZ,flip, flop, rotate, cropright, cropbottom, cropleft, croptop, scale)
	toflip = flip ? `-flip` : ``
	toflop = flop ? `-flop` : ``
	sz = SZ*scale/100
	w0 = sz[1]*cosd(rotate) +  sz[2]*cosd(90 - abs(rotate))
	w = round(Int, w0*cropright/100)
	h0 = sz[1]*sind(abs(rotate)) +  sz[2]*sind(90 - abs(rotate))
	h = round(Int, h0*cropbottom/100)
	x = round(Int, w0*cropleft/100)
	y = round(Int, h0*croptop/100)
	return ([w,h], `-scale $scale% -rotate $rotate +repage -crop $(w)x$h+$x+$y $toflip $toflop `)
end
#@everywhere
function myimread(fname::AbstractString,sz::Vector{Int},opt::Cmd)
	stream, _ = open(pipeline(`dcraw -w -H 0 -o 0 -h -4 -c $fname`,`convert - $opt -channel G -separate  Gray:-`))
	read(stream, UInt16, 1, sz...)
end
function loaddata{T <: AbstractString}(fnames::Vector{T},sz::Vector{Int},opt::Cmd)
	n = length(fnames)
	#x = Shared
	x = Array(UInt16,n,sz...)
	#@sync @parallel
	for i = 1:n
		x[i,:,:] = myimread(fnames[i],sz,opt)
	end
	return x
end
function data2stokes(x::Array{UInt16,3}, angleoffset::Real)
	n,w,h = size(x)
	S = Array(Stokes,w,h)
	m = zeros(n)
	for i = 1:w*h
		for j = 1:n
			m[j] = (x[j,i] + 1)/65536
		end
		tmp = Stokes(measurements2stokes(m...)...)
		S[i] = rotate(tmp,angleoffset)
	end
	return S
end
function path2stokes(path::AbstractString,opts...; angleoffset::Real = 0)
	fnames = getfnames(path)
	SZ = getsz(fnames[1])
	sz, opt = build_opt(SZ, opts...)
	x = loaddata(fnames,sz,opt)
	return data2stokes(x, angleoffset)
end
function path2stokes(path::AbstractString; angleoffset::Real = 0)
	fnames = getfnames(path)
	SZ = getsz(fnames[1])
	x = loaddata(fnames,SZ,` `)
	return data2stokes(x, angleoffset)
end

roundit(x::Float64) = round(Int, x*(N - 1) + 1)
function normalize(p::Polar)
	if isnan(p.I)
		return (1,1,rand(1:N),N2)
	else
		I = p.I > 1 ? N : roundit(p.I)
		dolp = p.dolp > 1 ? N : roundit(p.dolp)
		aop = I == 1 || dolp == 1 ? rand(1:N) : roundit(p.aop/pi + 0.5)
		tmp = (p.docp + 1)/2
		docp = tmp < 0 ? 1 : tmp > 1 ? N : roundit(tmp)
	end
	return (I,dolp,aop,docp)
end
function colorwheel(r::Int)
	row = Int[]
	col = Int[]
	ind = Int[]
	r2 = r*r
	for i in -r:r, j in -r:r
		if i*i + j*j <= r2
			push!(row,i + r + 1)
			push!(col,j + r + 1)
			tmp = atan2(j,i)/pi
			a = tmp > 0 ? roundit(tmp) : roundit(1 + tmp)
			push!(ind,a)
		end
	end
	return (row,col,ind)
end
function stokes2mats(s::Array{Stokes,2})
	sz = size(s)
	I = zeros(Int,sz)
	dolp = zeros(Int,sz)
	aop = zeros(Int,sz)
	docp = zeros(Int,sz)
	for i in eachindex(s)
		I[i], dolp[i], aop[i], docp[i] = normalize(Polar(s[i]))
	end
	return (I, dolp, aop, docp)
end
function mats2images(I::Array{Int,2},dolp::Array{Int,2},aop::Array{Int,2},docp::Array{Int,2})
	sz = size(I)
	row,col,ind = colorwheel(round(Int,0.1*mean(sz)))
	buff = round(Int,0.025*mean(sz))
	for (r,c,i) in zip(row,col,ind)
		aop[sz[1] - r - buff, c + buff] = i
	end
	Iimg = ImageCmap(I, linspace(RGB(0,0,0),RGB(1,1,1),N); prop...)
	colorbar = repeat(reshape(collect(round(Int,linspace(N,1,sz[2]))),1,sz[2]),outer = [round(Int,0.1sz[1]),1])
	dolpimg = ImageCmap(cat(1,dolp,colorbar),linspace(RGB(0,0,0),RGB(1,1,1),N); prop...)
	aopimg = convert(Image, ImageCmap(aop,convert(Array{RGB{FixedPointNumbers.UFixed{UInt8,8}},1}, linspace(HSV(0,1,1), HSV(360,1,1),N)); prop...))
	docpimg = ImageCmap(cat(1,docp,colorbar), [linspace(RGB(1,0,0),RGB(0,0,0),Int(N/2)); linspace(RGB(0,0,0), RGB(0,1,0),Int(N/2))]; prop...)
	return (Iimg, dolpimg, aopimg, docpimg)
end

profile(a::Stokes,ψ) = (a.s0 + ψ*a.s1)/(a.s0 - ψ*a.s1)
normalizeprofile(x,Sp) = (x - 1/Sp)/(Sp - 1/Sp)

function stokes2profile(s::Array{Stokes,2}, Sp::Number)
	sz = size(s)
	p = zeros(Int,sz)
	ψ = (Sp - 1)/(Sp + 1)
	for i in eachindex(s)
		p[i] = isnan(s[i].s0) ? rand(1:N) : roundit(normalizeprofile(profile(s[i],ψ),Sp))
	end
	colorbar = repeat(reshape(collect(round(Int,linspace(N,1,sz[2]))),1,sz[2]),outer = [round(Int,0.1sz[1]),1])
	ImageCmap(cat(1,p,colorbar),linspace(RGB(0,0,0),RGB(1,1,1),N); prop...)
end

import Base: /, *
immutable Centimeter
	x::Real
end
*(x::Centimeter,y::Real) = x.x*y
*(y::Real,x::Centimeter) = x*y
/(x::Centimeter,y::Real) = x.x/y
/(y::Real,x::Centimeter) = x/y
immutable Pixel
	x::Integer
end
*(x::Pixel,y::Real) = x.x*y
*(y::Real,x::Pixel) = x*y
/(x::Pixel,y::Real) = x.x/y
/(y::Real,x::Pixel) = x/y
immutable Length
	cm::Centimeter
	px::Pixel
end
Length(x::Centimeter,RATIO::Real) = Length(x,Pixel(round(Int,x/RATIO)))
Length(x::Pixel,RATIO::Real) = Length(Centimeter(x*RATIO),x)
findwindow(x::Real,RATIO::Real) = round(x/2/RATIO)*RATIO
function sense(el::PolEllipse,window::Length,maxpol::Real)
	p = round(el.polarization/maxpol,4)
	a = round(window.cm.x*el.a/maxpol,4)
	b = round(window.cm.x*el.b/maxpol,4)
	theta = round(rad2deg(el.angle),4)
	color = el.lefthand ? "dolcpcolor" : "dorcpcolor"
	return (color,p,theta,a,b)
end
function point2measurement(p1::Int,p2::Int,s::Array{Stokes,2},w::Int)
	m = Stokes(0.,0.,0.,0.)
	for j = -w:w, k = -w:w
		tmp = s[p1 + j, p2 + k]
		m += isnan(tmp.s0) ? Stokes(0.,0.,0.,0.) : tmp
	end
	return m/4*w*w
end
function grid(s::Array{Stokes,2},sz::Tuple{Length,Length},window::Length,maxpol::Real,RATIO::Real, minpol::Real, eccentricity::Real, name::AbstractString, linearcolor::RGB)
	open("ellipse.tikz","w") do o
		println(o,"""\\begin{tikzpicture}
			\\begin{axis}[xmin=0,xmax=1,ymin=0,ymax=1, ticks=none, axis lines=none, width = $(sz[1].cm.x) cm, height = $(sz[2].cm.x) cm, scale only axis]
				\\addplot graphics [xmin=0,xmax=1,ymin=0,ymax=1]{$name.jpg};
			\\end{axis}""")

		for i = (1 + window.px.x):2window.px.x:(sz[1].px.x - window.px.x), j = (1 + window.px.x):2window.px.x:(sz[2].px.x - window.px.x)
			ss = point2measurement(i,sz[2].px.x - j + 1,s,window.px.x)
			isnan(ss.s0) && continue
			el = PolEllipse(ss)
			el.polarization < minpol && continue
			#isnan(el.a) < minpol && continue
			color,p,theta,a,b = sense(el,window,maxpol)
			x = RATIO*(i - .5)
			y = RATIO*(j - .5)
			#y = RATIO*(sz[2].px.x - (j - 1.5))
			if a/b < eccentricity
				x1 = x + b*cosd(theta + 180 + 90)
				y1 = y + b*sind(theta + 180 + 90)
				x2 = x1 + 2b*cosd(theta + 90)
				y2 = y1 + 2b*sind(theta + 90)
				#println(o,"\\draw [$color,draw opacity = $p] ($x1 cm,$y1 cm) -- ($x2 cm, $y2 cm);")
				println(o,"\\draw [dolpcolor] ($x1 cm,$y1 cm) -- ($x2 cm, $y2 cm);")
			else
				#println(o,"\\draw [$color,draw opacity = $p] ($x cm,$y cm) circle [x radius=$a cm, y radius=$b cm, rotate=$theta];")
				println(o,"\\draw [$color] ($x cm,$y cm) circle [x radius=$a cm, y radius=$b cm, rotate=$theta];")
			end
		end
		println(o,"""\\draw[step = $(2window.cm.x) cm, gray, very thin, draw opacity=0.25] (0,0) grid ($(sz[1].cm.x) cm, $(sz[2].cm.x) cm);
	\\end{tikzpicture}""")
	end
end
function stokes2ellipse(s::Array{Stokes,2}; dpi::Int = 96, window_in_cm::Real = 0.15, maxpol::Real = 1, minpol::Real = 0, eccentricity::Real = 0.01, name = "RGB", linearcolor::RGB = RGB(0.0,1.0,0.0))
	SZ = size(s)
	const RATIO = 2.54/dpi
	sz = (Length(Pixel(SZ[1]),RATIO),Length(Pixel(SZ[2]),RATIO))
	window = Length(Centimeter(findwindow(window_in_cm,RATIO)),RATIO)
	grid(s,sz,window,maxpol,RATIO, minpol, eccentricity,name,linearcolor)
end

end


text:
using Documenter, ConstrainedLasso

# makedocs(
# #      doctest   = false,
#       format    = :html,
#       clean     = true,
#       sitename  = "ConstrainedLasso.jl",
#       modules   = [ConstrainedLasso],
#       pages     = [
#           "Home"     => "index.md",
# 	  "Interface" => "interface.md",
#           "References" => "references.md"
#       ]
#
# )
#
# deploydocs(
#       repo    = "github.com/Hua-Zhou/ConstrainedLasso.jl.git",
#       target  = "build",
#       julia   = "0.6",
#       deps    = nothing,
#       make    = nothing
# )

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs()

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo   = "github.com:Hua-Zhou/ConstrainedLasso.jl.git",
    julia  = "0.6.2",
    osname = "osx"
)


text:
using SimpleWeightedGraphs
using LinearAlgebra

include("tsptwtrajectorystate.jl")

struct TsptwFeaturization <: AbstractFeaturization end

"""
    TsptwStateRepresentation{F, TS}

This is the Tsptw representation used by Quentin Cappart in Combining Reinforcement Learning and Constraint Programming
for Combinatorial Optimization (https://arxiv.org/pdf/2006.01610.pdf).
"""
mutable struct TsptwStateRepresentation{F, TS} <: FeaturizedStateRepresentation{F, TS}
    dist::Matrix
    timeWindows::Matrix
    pos::Matrix
    nodeFeatures::Union{Nothing, AbstractMatrix{Float32}}
    variableIdx::Union{Nothing, Int64}
    possibleValuesIdx::Union{Nothing, AbstractVector{Int64}}
end

function TsptwStateRepresentation{F, TS}(model::CPModel; action_space=nothing) where {F, TS}
    dist, timeWindows, pos = get_tsptw_info(model)
    sr = TsptwStateRepresentation{F, TS}(dist, timeWindows, pos, nothing, nothing, nothing)
    sr.nodeFeatures = featurize(sr)
    return sr
end

TsptwStateRepresentation(model::CPModel) = TsptwStateRepresentation{TsptwFeaturization, TsptwTrajectoryState}(model)

function TsptwTrajectoryState(sr::TsptwStateRepresentation{F, TsptwTrajectoryState}) where F
    # TODO change this once the InitializingPhase is fixed
    if isnothing(sr.variableIdx)
        sr.variableIdx = 1
    end
    if isnothing(sr.possibleValuesIdx)
        throw(ErrorException("Unable to build a DefaultTrajectoryState, when the possible values vector is nothing."))
    end

    n = size(sr.dist, 1)
    adj = ones(Int, n, n) - I
    edgeFeatures = build_edge_feature(adj, sr.dist)
    fg = GeometricFlux.FeaturedGraph(adj; nf=sr.nodeFeatures, ef=edgeFeatures)
    return TsptwTrajectoryState(fg, sr.variableIdx, sr.possibleValuesIdx)
end

function get_tsptw_info(model::CPModel)
    dist, timeWindows, pos, grid_size = model.adhocInfo

    max_d = Base.maximum(dist)
    max_tw = Base.maximum(timeWindows)

    dist = dist ./ max_d
    timeWindows = timeWindows ./ max_tw
    pos = pos ./ grid_size

    return dist, timeWindows, pos
end

function update_representation!(sr::TsptwStateRepresentation, model::CPModel, x::AbstractIntVar)
    sr.possibleValuesIdx = collect(x.domain)

    i = 1
    while x.id != "a_"*string(i)
        i += 1
    end
    if SeaPearl.isbound(model.variables["v_"*string(i)])
        sr.variableIdx = SeaPearl.assignedValue(model.variables["v_"*string(i)])
    end
    update_features!(sr, model)
    return sr
end

function build_edge_feature(adj::AbstractMatrix, weighted_adj::AbstractMatrix)
    adj_list = GeometricFlux.adjacency_list(adj)
    n = length(adj_list)
    return hcat([build_edge_feature_aux(i, adj_list[i], weighted_adj) for i in 1:n]...)
end

function build_edge_feature_aux(i, js, weighted_adj::AbstractMatrix)
    hcat([weighted_adj[i, j] for j = js]...)
end

"""
    function featurize(sr::TsptwStateRepresentation{TsptwFeaturization})

Create nodeFeatures for every node of the graph. Supposed to be overwritten.
Tsptw behavior is to call `Tsptw_featurize`.
"""
function featurize(sr::FeaturizedStateRepresentation{TsptwFeaturization, TS}) where TS
    n = size(sr.dist, 1)
    nodeFeatures = zeros(Float32, 6, n)
    for i in 1:n
        if !isnothing(sr.possibleValuesIdx) && !(i in sr.possibleValuesIdx)
            nodeFeatures[5, i] = 1.
        end
        if i == sr.variableIdx
            nodeFeatures[6, i] = 1.
        end
    end

    nodeFeatures[1:2, :] = transpose(sr.pos)
    nodeFeatures[3:4, :] = transpose(sr.timeWindows)
    return nodeFeatures
end

function update_features!(sr::FeaturizedStateRepresentation{TsptwFeaturization, TS}, model::CPModel) where TS
    sr.nodeFeatures[5, :] .= 0
    if !isnothing(sr.possibleValuesIdx)
        sr.nodeFeatures[5, sr.possibleValuesIdx] .= 1
    end
    sr.nodeFeatures[6, :] .= 0
    if !isnothing(sr.variableIdx)
        sr.nodeFeatures[6, sr.variableIdx] = 1
    end
    return
end

"""
    feature_length(gen::TsptwGenerator, ::Type{TsptwStateRepresentation{TsptwFeaturization}})

Returns the length of the feature vector, useful for SeaPearl to choose the size of the container
"""
feature_length(::Type{TsptwStateRepresentation{TsptwFeaturization, TS}}) where TS = 6


text:
function coeffAero(f::Fragment, m::Float64)
    if m<0.1
        cd = f.cd_sub
    elseif m<0.75
        cd = 0.2/0.65*(m-0.75)+(f.cd_sub+0.2)
    elseif m<1.5
        cd = 0.45/(1.5-0.75)*(m-1.5)+(f.cd_sub+0.65)
    elseif m<3.0
        cd = -0.15/1.5*(m-3.0)+(f.cd_sub+0.5)
    else
        cd = f.cd_sub+0.5
    end
    return cd

end


text:
using Test, LinearAlgebra
using DiffEqSensitivity, StochasticDiffEq
using ForwardDiff
using Random

@info "SDE Non-Diagonal Noise Adjoints"

seed = 100
Random.seed!(seed)

tstart = 0.0
tend = 0.1
dt = 0.005
trange = (tstart, tend)
t = tstart:dt:tend
tarray = collect(t)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

# non-diagonal noise
@testset "Non-diagonal noise tests" begin
    Random.seed!(seed)

    u₀ = [0.75,0.5]
    p = [-1.5,0.05,0.2, 0.01]

    dtnd = tend/1e3

    # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
    function f_nondiag!(du,u,p,t)
      du[1] = p[1]*u[1] + p[2]*u[2]
      du[2] = p[2]*u[1] + p[1]*u[2]
      nothing
    end

    function g_nondiag!(du,u,p,t)
      du[1,1] = p[3]*u[1] + p[4]*u[2]
      du[1,2] = p[3]*u[1] + p[4]*u[2]
      du[2,1] = p[4]*u[1] + p[3]*u[2]
      du[2,2] = p[4]*u[1] + p[3]*u[2]
      nothing
    end

    function f_nondiag(u,p,t)
      dx = p[1]*u[1] + p[2]*u[2]
      dy = p[2]*u[1] + p[1]*u[2]
      [dx,dy]
    end

    function g_nondiag(u,p,t)
      du11 = p[3]*u[1] + p[4]*u[2]
      du12 = p[3]*u[1] + p[4]*u[2]
      du21 = p[4]*u[1] + p[3]*u[2]
      du22 = p[4]*u[1] + p[3]*u[2]

      [du11  du12
        du21  du22]
    end

    function f_nondiag_analytic(u0,p,t,W)
      A = [[p[1], p[2]] [p[2], p[1]]]
      B = [[p[3], p[4]] [p[4], p[3]]]
      tmp = A*t + B*W[1] + B*W[2]
      exp(tmp)*u0
    end

    noise_matrix = similar(p,2,2)
    noise_matrix .= false

    Random.seed!(seed)
    prob = SDEProblem(f_nondiag!,g_nondiag!,u₀,trange,p,noise_rate_prototype=noise_matrix)
    sol = solve(prob, EulerHeun(), dt=dtnd, save_noise=true)

    noise_matrix = similar(p,2,2)
    noise_matrix .= false
    Random.seed!(seed)
    proboop = SDEProblem(f_nondiag,g_nondiag,u₀,trange,p,noise_rate_prototype=noise_matrix)
    soloop = solve(proboop,EulerHeun(), dt=dtnd, save_noise=true)

    res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint())

    @info res_sde_p

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint())

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test_broken isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test_broken isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint())

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test_broken isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test_broken isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    function compute_grads_nd(sol)
      xdis = sol(tarray)

      mat1 = Matrix{Int}(I, 2, 2)
      mat2 = ones(2,2)-mat1

      tmp1 = similar(p)
      tmp1 *= false

      tmp2 = similar(xdis.u[1])
      tmp2 *= false

      for (i, u) in enumerate(xdis)
        tmp1[1]+=xdis.t[i]*u'*mat1*u
        tmp1[2]+=xdis.t[i]*u'*mat2*u
        tmp1[3]+=sum(sol.W(xdis.t[i])[1])*u'*mat1*u
        tmp1[4]+=sum(sol.W(xdis.t[i])[1])*u'*mat2*u

        tmp2 += u.^2
      end

      return tmp2 ./ xdis.u[1], tmp1
    end

    res1, res2 = compute_grads_nd(soloop)

    @test isapprox(res1, res_sde_u0, rtol=1e-5)
    @test isapprox(res2, res_sde_p', rtol=1e-4)

end

@testset "diagonal but mixing noise tests" begin

  Random.seed!(seed)
  u₀ = [0.75,0.5]
  p = [-1.5,0.05,0.2, 0.01]
  dtmix = tend/1e3

  # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
  function f_mixing!(du,u,p,t)
    du[1] = p[1]*u[1] + p[2]*u[2]
    du[2] = p[2]*u[1] + p[1]*u[2]
    nothing
  end

  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end

  function f_mixing(u,p,t)
    dx = p[1]*u[1] + p[2]*u[2]
    dy = p[2]*u[1] + p[1]*u[2]
    [dx,dy]
  end

  function g_mixing(u,p,t)
    dx = p[3]*u[1] + p[4]*u[2]
    dy = p[3]*u[1] + p[4]*u[2]
    [dx,dy]
  end

  Random.seed!(seed)
  prob = SDEProblem(f_mixing!,g_mixing!,u₀,trange,p)

  soltsave = collect(trange[1]:dtmix:trange[2])
  sol = solve(prob, EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave )

  Random.seed!(seed)
  proboop = SDEProblem(f_mixing,g_mixing,u₀,trange,p)
  soloop = solve(proboop,EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  #oop

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @info res_sde_p

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=false,noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true, noise=DiffEqSensitivity.ZygoteNoise()))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=false, noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(),noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_p

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise(), noisemixing=true))

  @test_broken isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test_broken isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=false,noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5) # would pass with 1e-4 but last noise value is off

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true, noise=DiffEqSensitivity.ZygoteNoise()))

  @test_broken isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test_broken isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=false, noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(),noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  function GSDE(p)
    Random.seed!(seed)
    tmp_prob = remake(prob,u0=eltype(p).(prob.u0),p=p,
                  tspan=eltype(p).(prob.tspan))
    _sol = solve(tmp_prob,EulerHeun(),dt=dtmix,adaptive=false,saveat=Array(t))
    A = convert(Array,_sol)
    res = g(A,p,nothing)
  end

  res_sde_forward = ForwardDiff.gradient(GSDE,p)

  @test isapprox(res_sde_p', res_sde_forward, rtol=1e-5)

  function GSDE2(u0)
    Random.seed!(seed)
    tmp_prob = remake(prob,u0=u0,p=eltype(p).(prob.p),
                      tspan=eltype(p).(prob.tspan))
    _sol = solve(tmp_prob,EulerHeun(),dt=dtmix,adaptive=false,saveat=Array(t))
    A = convert(Array,_sol)
    res = g(A,p,nothing)
  end

  res_sde_forward = ForwardDiff.gradient(GSDE2,u₀)

  @test isapprox(res_sde_forward, res_sde_u0, rtol=1e-5)

end

@testset "mixing noise inplace/oop tests" begin
  Random.seed!(seed)
  u₀ = [0.75,0.5]
  p = [-1.5,0.05,0.2, 0.01]
  dtmix = tend/1e3

  # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
  function f_mixing!(du,u,p,t)
    du[1] = p[1]*u[1] + p[2]*u[2]
    du[2] = p[2]*u[1] + p[1]*u[2]
    nothing
  end

  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end

  function f_mixing(u,p,t)
    dx = p[1]*u[1] + p[2]*u[2]
    dy = p[2]*u[1] + p[1]*u[2]
    [dx,dy]
  end

  function g_mixing(u,p,t)
    dx = p[3]*u[1] + p[4]*u[2]
    dy = p[3]*u[1] + p[4]*u[2]
    [dx,dy]
  end

  Random.seed!(seed)
  prob = SDEProblem(f_mixing!,g_mixing!,u₀,trange,p)
  soltsave = collect(trange[1]:dtmix:trange[2])
  sol = solve(prob, EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  Random.seed!(seed)
  proboop = SDEProblem(f_mixing,g_mixing,u₀,trange,p)
  soloop = solve(proboop,EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  @test sol.u ≈ soloop.u atol = 1e-14

  # BacksolveAdjoint

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @test res_sde_u0 ≈ res_sde_u02 atol = 1e-14
  @test res_sde_p ≈ res_sde_p2 atol = 1e-14

  @show res_sde_u0

  adjproboop = SDEAdjointProblem(soloop,BacksolveAdjoint(noisemixing=true),dg!,tarray, nothing)
  adj_soloop = solve(adjproboop,EulerHeun(); dt=dtmix, tstops=soloop.t, adaptive=false)

  @test adj_soloop[end][length(p)+length(u₀)+1:end] == soloop.u[1]
  @test - adj_soloop[end][1:length(u₀)] == res_sde_u0
  @test adj_soloop[end][length(u₀)+1:end-length(u₀)] == res_sde_p'

  adjprob = SDEAdjointProblem(sol,BacksolveAdjoint(noisemixing=true, checkpointing=true),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-15

  adjprob = SDEAdjointProblem(sol,BacksolveAdjoint(noisemixing=true, checkpointing=false),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-8

  # InterpolatingAdjoint

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test res_sde_u0 ≈ res_sde_u02 atol = 1e-8
  @test res_sde_p ≈ res_sde_p2 atol = 5e-8

  @show res_sde_u0

  adjproboop = SDEAdjointProblem(soloop,InterpolatingAdjoint(noisemixing=true),dg!,tarray, nothing)
  adj_soloop = solve(adjproboop,EulerHeun(); dt=dtmix, tstops=soloop.t, adaptive=false)

  @test - adj_soloop[end][1:length(u₀)] ≈ res_sde_u0  atol = 1e-14
  @test adj_soloop[end][length(u₀)+1:end] ≈ res_sde_p' atol = 1e-14

  adjprob = SDEAdjointProblem(sol,InterpolatingAdjoint(noisemixing=true, checkpointing=true),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈ adj_sol[end]  rtol=1e-8

  adjprob = SDEAdjointProblem(sol,InterpolatingAdjoint(noisemixing=true, checkpointing=false),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-8
end


text:
mutable struct KalmanFilter
    estimate_position::Float64
    estimate_velocity::Float64
    p_xx::Float64
    p_xv::Float64
    p_vv::Float64
    position_variance::Float64
    velocity_variance::Float64
    r::Float64
    has_velocity::Bool

    function KalmanFilter()
        return KalmanFilter(0.0, false)
    end

    function KalmanFilter(position::Float64, video::Bool)
        local estimate_position, estimate_velocity, has_velocity
        local p_xx, p_xv, p_vv
        local position_variance, velocity_variance, r

        estimate_position = position
        estimate_velocity = 0.0
        has_velocity = false

        if !video
            # normal settings used for tracking X-Y coordinates
            p_xx = 0.1
            p_xv = 0.1
            p_vv = 0.1
            position_variance = 0.01
            velocity_variance = 0.01
            r = 100000.0
        else
            # special setting used by video streaming
            p_xx = 0.1 * position
            p_xv = 1.0
            p_vv = 1.0
            position_variance = 0.1 * position
            velocity_variance = 0.01 * position / 1000.0
            r = 0.01 * position
        end

        return new(
            estimate_position,
            estimate_velocity,
            p_xx,
            p_xv,
            p_vv,
            position_variance,
            velocity_variance,
            r,
            has_velocity,
        )
    end
end

function reset(filter::KalmanFilter, position::Float64)
    filter.estimate_position = position
    filter.estimate_velocity = 0.0
    filter.has_velocity = false
end

function update(filter::KalmanFilter, position::Float64, deltat::Float64)
    if deltat <= 0.0
        return
    end

    if !filter.has_velocity
        filter.estimate_velocity = (position - filter.estimate_position) / deltat
        filter.estimate_position = position
        filter.has_velocity = true
    else
        # Temporal update (predictive)
        filter.estimate_position += filter.estimate_velocity * deltat

        # Update covariance
        filter.p_xx += deltat * (2.0 * filter.p_xv + deltat * filter.p_vv)
        filter.p_xv += deltat * filter.p_vv

        filter.p_xx += deltat * filter.position_variance
        filter.p_vv += deltat * filter.velocity_variance

        # Observational update (reactive)
        vi = 1.0 / (filter.p_xx + filter.r)
        kx = filter.p_xx * vi
        kv = filter.p_xv * vi

        filter.estimate_position += (position - filter.estimate_position) * kx
        filter.estimate_velocity += (position - filter.estimate_position) * kv

        filter.p_xx *= 1.0 - kx
        filter.p_xv *= 1.0 - kx
        filter.p_vv -= kv * filter.p_xv
    end
end

function predict(filter::KalmanFilter, position::Float64, deltat::Float64)::Float64
    return position + filter.estimate_velocity * deltat
end


text:
@generated function staticschema(::Type{T}) where {T}
    name_tuple = Expr(:tuple, [QuoteNode(f) for f in fieldnames(T)]...)
    type_tuple = Expr(:curly, :Tuple, [Expr(:call, :fieldtype, :T, i) for i in 1:fieldcount(T)]...)
    Expr(:curly, :NamedTuple, name_tuple, type_tuple)
end

@generated function staticschema(::Type{T}) where {T<:Tuple}
    name_tuple = Expr(:tuple, [QuoteNode(Symbol("x$f")) for f in fieldnames(T)]...)
    type_tuple = Expr(:curly, :Tuple, [Expr(:call, :fieldtype, :T, i) for i in 1:fieldcount(T)]...)
    Expr(:curly, :NamedTuple, name_tuple, type_tuple)
end

staticschema(::Type{T}) where {T<:NamedTuple} = T

getnames(::Type{NamedTuple{names, types}}) where {names, types} = names
gettypes(::Type{NamedTuple{names, types}}) where {names, types} = types

function fields(::Type{T}) where {T}
    getnames(staticschema(T))
end


text:
# This file consists of code that is simplified/modified/inspired by/from
# the DSP.jl package. The purpose of these functions is to speed up repeated
# power spectrum computations for real-valued 1D vectors, by utilizing fft-plans
# and in-place computations.
#
# The DSP module is distributed under the MIT license.
# Copyright 2012-2021 DSP.jl contributors
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software withspectrum restriction, including withspectrum limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
using DSP
using LinearAlgebra: mul!

"""
    prepare_spectrum(s, plan) → y::Vector

Pre-allocate a vector that will hold the one-sided power spectrum for a signal `s`,
taking into account the number of points that will be used for the Fourier transform.
"""
function prepare_spectrum(s, plan)
    return plan * s
end

"""
    _powerspectrum_from_fft!(
        spectrum::AbstractArray{T},
        ℱ::AbstractVector{Complex{T}},
        nfft::Int,
        r::Real,
        offset::Int = 0
        ) where T

Compute one-sided power spectrum from a pre-computed Fourier transform `ℱ` of a signal,
into a pre-allocated `spectrum vector``. `nfft` is the number of points used for
the Fourier transform. `n` is the length of the signal,
"""
function _powerspectrum_from_fft!(
        spectrum::AbstractArray{T},
        ℱ::AbstractVector{Complex{T}},
        nfft::Int,
        n::Real,
        offset::Int = 0
        ) where T

    l = length(ℱ)
    m1, m2 = convert(T, 1 / n), convert(T, 2 / n)

    for i = 2:l-1
       @inbounds spectrum[offset + i] += abs2(ℱ[i]) * m2
    end
    @inbounds spectrum .+= abs2.(ℱ) .* m2
    @inbounds spectrum[offset + l] += abs2(ℱ[end]) * ifelse(iseven(nfft), m1, m2)

    return spectrum
end

"""
    powerspectrum_onesided!(spectrum, signal, forward)

Let `n = DSP.nextfastfft(length(s))`. Modifies `spectrum` in-place, so `spectrum`
is reset to all zeros every time this function is called.

- `forward` is a forward fft-plan for the signal.
- `spectrum` is a pre-allocated `AbstractVector{<:Real}` that will hold the spectrum.
    Has the length of the Fourier transform resulting from `forward * signal`.
- `signal` is the signal, a `AbstractVector{<:Real}` of length `n`.
"""
function powerspectrum!(ℱp, spectrum, signal, forward)
    # Fourier transform of the signal, based on pre-computed plan `forward`.
    mul!(ℱp, forward, signal) #ℱp = forward * signal

    l = length(signal)

    # Reset spectrum, since we're doing multiple in-place additions to it.
    spectrum .= 0.0

    # In-place computation of spectrum based on the transform `ℱ`
    n = nextfastfft(l)
    _powerspectrum_from_fft!(spectrum, ℱp, n, l)

    return spectrum
end


text:
abstract type SrMat{T} <: LinearMap{T} end

Base.eltype(S::SrMat{T}) where T = T
LinearAlgebra.ishermitian(S::SrMat) = true

# Lazy Algebra: addition/subtraction
add!(S::SrMat, α::Number) = (S.shift += α; return S)
subtract!(S::SrMat, α) = add!(S, -α)
shift!(S::SrMat, α) = (S.shift = α; return S)

#Base.:+(S::SrMat, α::Number) = add!(copy(S), α)
Base.:+(S::SrMat, α::LinearAlgebra.UniformScaling) = add!(copy(S), α.λ)

# Constructor
SrMatrix(T, O::AbstractMatrix, α::Number, args...) = begin
    if T <: Real
        S = SrMatrixR(O, α, args...)
    else
        S = SrMatrixC(O, α, args...)
    end
    return S
end

mutable struct SrMatrixC{T,Tv<:AbstractMatrix{T},Ct,Cr,Pc} <: SrMat{T}
    O::Tv
    shift::T

    # cache
    ṽ::Ct
    res::Cr
    parallel_cache::Pc
end

SrMatrixC(O::AbstractMatrix, α::Number=0, par_cache=NotParallel()) =
    SrMatrixC(O, convert(eltype(O), α),
              similar(O, size(O,2)), similar(O, size(O,1)),
              par_cache)

Base.copy(S::SrMatrixC) = SrMatrixC(S.O, S.shift,
                                    S.ṽ, S.res,
                                    S.parallel_cache)

Base.size(S::SrMatrixC) = (size(S.O,1), size(S.O,1))
Base.isreal(S::SrMatrixC) = isreal(S.O)
LinearAlgebra.issymmetric(S::SrMatrixC) = isreal(S)

function init!(S::SrMatrixC, O)
    # Check size
    if size(S.O) != size(O)
        S.O = conj(O)
        S.ṽ = similar(O, size(O,2))
    else
        S.O .= conj.(O)
    end

    return S
end

function LinearAlgebra.mul!(C::AbstractVector, S::SrMatrixC{T}, F::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
    O     = S.O
    shift = S.shift

    #cache
    ṽ   = S.ṽ
    res = S.res

    # calculation
    mul!(ṽ, O', F)
    mul!(res, O, ṽ)
    workers_sum!(res, S.parallel_cache)

    N = size(O,2) * num_workers(S.parallel_cache)

    if β == 0
        C .= α .* (F .* shift .+ res ./ N)
    else
        C .= α .* (F .* shift .+ res ./ N) .+ β .* C
    end
end

##
mutable struct SrMatrixR{T,Tv<:AbstractMatrix{T},Ct,Cr,Pc} <: SrMat{T}
    Oᵣ::Tv
    Oᵢ::Tv
    shift::T

    # cache
    ṽ::Ct
    res::Cr
    parallel_cache::Pc
end

SrMatrixR(O::AbstractMatrix, α::Number=0, par_cache=NotParallel()) =
    SrMatrixR(real(O),
              imag(O),
              convert(eltype(O), α),
              similar(O, real(eltype(O)), size(O,2)),
              similar(O, size(O,1)),
              par_cache)

Base.copy(S::SrMatrixR) = SrMatrixR(S.Oᵣ, S.Oᵢ, S.shift,
                                    S.ṽ, S.res, S.parallel_cache)

Base.size(S::SrMatrixR) = (size(S.Oᵣ,1), size(S.Oᵣ,1))

init!(S::SrMatrixR, O) = begin
    # Check size
    # Check size
    if size(S.Oᵣ) != size(O)
        S.Oᵣ = real(O)
        S.Oᵢ = imag(O)
        S.ṽ  = similar(O, real(eltype(O)), size(O,2))
    else
        S.Oᵣ .= real.(O)
        S.Oᵢ .= imag.(O)
    end

    return S
end

function LinearAlgebra.mul!(C::AbstractVector, S::SrMatrixR{T}, F::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
    Oᵣ     = S.Oᵣ
    Oᵢ     = S.Oᵢ
    shift = S.shift

    #cache
    ṽ   = S.ṽ
    res = S.res

    # calculation real part
    mul!(ṽ, Oᵣ', F)
    mul!(res, Oᵣ, ṽ)

    # calculation imaginary part
    mul!(ṽ, Oᵢ', F)
    mul!(res, Oᵢ, ṽ, one(eltype(ṽ)), one(eltype(ṽ)))
    workers_sum!(res, S.parallel_cache)

    N = size(Oᵣ,2) * num_workers(S.parallel_cache)

    if β == 0
        C .= α .* (F .* shift .+ res ./ N)
    else
        C .= α .* (F .* shift .+ res ./ N) .+ β .* C
    end
end


text:
# %% imports

using JLD
using EvolutionaryStrategies

include("cartpole_fitness.jl")
#include("swingup_fitness.jl")
#include("mountcar_fitness.jl")
include("cambrianrun.jl")

# %% create snestopgen
init = 0
save("Julia/Actors/topgenxNES.jld","init", init)
# %% xNES optimisation

for opti in 1:50
    cfg = get_config("Julia/Actors/cartpole_esconfig.yml")
    es = xNES(cfg,fitness)
    a = run!(es, opti)
    println("optim $(opti)")
end

# %% plot the distribution of gen convergence

f = load("Julia/Actors/topgenxNES.jld")
gens = Int64[]
for i in 1:length(f)-1
    push!(gens, f["$(i)_run"])
end

histogram(f, bins = 0:10:300,
            label = "First R>195 individual (mean = 82.68)",
            color = :green,
            xlabel = "Generation",
            ylabel = "Number of optimisations")

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the string alignment.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Alignment to the right" begin
    str = "Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test"

    expected = "                    Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test"
    aligned_str = align_string(str, 40, :r)

    @test aligned_str === expected
    @test printable_textwidth(aligned_str) == 40

    aligned_str = align_string(str, 40, :r; fill = true)

    @test aligned_str === expected
    @test printable_textwidth(aligned_str) == 40

    aligned_str = align_string(str, 10, :r)
    @test str == aligned_str
end

@testset "Alignment to the center" begin
    str = "Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test"
    expected = "          Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test"
    aligned_str = align_string(str, 40, :c)

    @test aligned_str === expected
    @test printable_textwidth(aligned_str) == 30

    expected = "          Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test          "
    aligned_str = align_string(str, 40, :c; fill = true)

    @test aligned_str === expected
    @test printable_textwidth(aligned_str) == 40

    aligned_str = align_string(str, 10, :c)
    @test str == aligned_str
end

@testset "Alignment to the left" begin
    str = "Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test"
    aligned_str = align_string(str, 40, :l)

    @test aligned_str === str
    @test printable_textwidth(aligned_str) == 20

    expected = "Test 😃 \e[38;5;231;48;5;243mTest \e[0m😅 Test                    "
    aligned_str = align_string(str, 40, :l; fill = true)

    @test aligned_str == expected
    @test printable_textwidth(aligned_str) == 40

    aligned_str = align_string(str, 10, :l)
    @test str == aligned_str
end

@testset "Multiple lines alignment" begin
    str = """
    We have \e[38;5;231;48;5;243mhere\e[0m 😅😃 the first line
    We now have the 😊 \e[38;5;231;48;5;243msecond\e[0m 😃 line"""

    expected_str = "                  We have \e[38;5;231;48;5;243mhere\e[0m 😅😃 the first line\n                 We now have the 😊 \e[38;5;231;48;5;243msecond\e[0m 😃 line"
    aligned_str = align_string_per_line(str, 50, :r)

    @test aligned_str == expected_str
end

@testset "Corner cases" begin
    str = """
    We have \e[38;5;231;48;5;243mhere\e[0m 😅😃 the first line
    We now have the 😊 \e[38;5;231;48;5;243msecond\e[0m 😃 line"""

    aligned_str = align_string_per_line(str, -1, :l)

    @test aligned_str == str
end


text:
struct Console
end

function print(msg::Console)
    print(msg)
end

max_stars_repo_name: aole/RogueJ

max_issues_repo_name: aole/RogueJ

max_forks_repo_name: aole/RogueJ


text:
using Multigraphs, LightGraphs, SparseArrays
using Test

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

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

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

@testset "di_multigraph.jl" begin
    include("di_multigraph_adjlist.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 06240e71-df45-4f63-b797-4042432e6716
using FFTW, Plots, SampledSignals

# ╔═╡ f8ba933c-a1e3-11eb-2263-17f14870c1fa
md"# Basic Waveforms
First let's see what basic waveforms can be made directly."

# ╔═╡ bfacded1-6007-4e18-aba1-1b7325878dd2
begin
    fs = 8000
	freqs = SampleBuf(Array(LinRange(0, 2*pi*220, fs)), fs)
	t = Array(0:1/100:1) * 2 * pi
end

# ╔═╡ 65dbcdfe-2dbd-4b28-ab59-fb12425913ca
function plotfreqs(tone, title)
	freqs = fft(tone) |> fftshift
	domain = FFTW.fftfreq(length(tone), fs) |> fftshift
    plot(domain, 10*log10.(abs.(freqs)/fs), title=title,
		xlimit=(0, fs/2), ylimit=(-50, 0))
end

# ╔═╡ 4bd75a55-4457-4625-a184-a2f7490cbf0e
md"## Sinusoids
The builtin sin() and cos() functions produce sinusoidal waves, they are roughly
equivalent for this, being phase shifted versions of each other."

# ╔═╡ 664b9502-13de-4fdf-971c-66cf0ea47c3b
plot(t, sin.(t))

# ╔═╡ 5e805acc-c6b1-4de5-8bc1-5bdcbd84ca04
sine_tone = sin.(freqs)

# ╔═╡ 2a51d397-322b-414b-9d9f-a1e8ba6b2364
plotfreqs(sine_tone, "Sine Spectrum")

# ╔═╡ dd813478-66a1-4d3e-9872-322e7b521e11
md"## Modulo waveforms
### Sawtooth
Using the periodic behavior of the remainder operator we get a sawtooth wave.
However, making a sawtooth directly like this gives aliasing
since this equivalent to sampling an analog sawtooth wave with all harmonics.
Harmonics higher than half the sampling rate cannot accurately be represented
and appear as lower frequencies when reproducing the sound.
This aliasing shows up as noise in the signal that is visible on the frequency graph.
"

# ╔═╡ a2357e48-116a-4df1-a52d-100533bad25d
function saw(t)
	t/pi % 2 - 1
end

# ╔═╡ 50824ffc-622d-4a00-acce-558cecf79e55
plot(t, saw.(t))

# ╔═╡ 3ae066cb-b5f2-4913-8a73-20348a401026
saw_tone = saw.(freqs)

# ╔═╡ 7bd89cfd-9894-4644-b0e4-cae8a63cda9a
plotfreqs(saw_tone, "Aliased Saw Spectrum")

# ╔═╡ dfa8eb03-2955-4bb6-b9ca-eb7cd6583e57
md"
### Triangle
The jagged shape of the sawtooth wave, with it's sharp discontinuity,
is why it's higher harmonics are so loud.
If we use two ramps to remove the discontinuity we get a triangle wave
with harmonics that fall off more quickly.
The aliasing is reduced since the harmonics are softer by the time they reach the frequency limit. It also no longer has even harmonics, since the second half is the first half upside down."

# ╔═╡ a9f4a512-c8b4-474a-a93c-2cef99db534c
function tri(t)
	abs(t/pi%2-1)*2-1
end

# ╔═╡ f6f28913-cbfa-48f8-98ab-0d5ed267b350
plot(t, tri.(t))

# ╔═╡ 231185ea-8913-490d-8fe8-261983552457
tri_tone = tri.(freqs)

# ╔═╡ b627d7af-55d1-4d97-8f15-7f0d172f9b5a
plotfreqs(tri_tone, "Triangle Spectrum")

# ╔═╡ 99708158-d74e-4311-a10f-d0ac4f79d4f2
md" ### Polynomial waveforms
We can use polynomials of a sawtooth wave to get various periodic waveforms.
If the two endpoints are equal the wave will be continuous.
If the derivitives at the endpoints are equal the waveform will be smooth.
Smooth waveforms have harmonics that fall off even faster than continuous waveforms.

Using a parabola we get harmonic roll-off similar to a triangle wave
but including the even harmonics:
"

# ╔═╡ 1f127661-47b0-4498-9fc2-857f314b6869
function poly2(t)
	"""
	Parabolic waveform
	since t^2 gives output in [0, 1] we use 2t^2-1 to get output in [-1,1]

	Note that since f' = 2x + b is linear we can't make a smooth function
	since no two points can be equal
	"""
	t = t/pi %2 -1
	2*t*t - 1
end

# ╔═╡ 99f21cb1-b50f-4e8b-9d97-8bf5297afbec
poly2_tone = poly2.(freqs)

# ╔═╡ 8802e5fe-8416-4c67-a1f0-b99bb6af3e3f
plot(t, poly2.(t))

# ╔═╡ b561d76e-e2a7-4c3d-9677-42f8b5661b7c
# Kind of like a triangle wave but with even harmonics
plotfreqs(poly2_tone, "Poly2 Spectrum")

# ╔═╡ 09d7e4f3-fed8-4fc8-b598-03ecc8e75d6f
md"With a cubic we can get a smooth function.
This has harmonics that fall off quickly enough that there is no noticable aliasing,
even at this low of a samplerate.
"

# ╔═╡ aff77549-c6da-4329-8e92-5f1fb561ece6
function poly3(t)
	"""
	f(x) = x^3 + bx^2 + cx
	function is continuous when
	1 + b + c = f(1) = f(-1) = -1 + b - c
	2c = -2, c = -1

	f'(x) = 3x^2 + 2bx - 1
	function is smooth when
    3 + 2b - 1 = f'(1) = f'(-1) = 3 - 2b - 1
	4b = 0, b = 0

	f(x) = x^3 - x
	f'(x) = 3x^2-1 = 0  when x = +-1/sqrt(3)
	f(-1/sqrt(3)) = -1/sqrt(3)/3+1/sqrt(3) = 2/(3sqrt(3))
	We scale by the inverse 3sqrt(3)/2
	"""
	t = t/pi%2-1
	x = 3*sqrt(3)/2*t*(t*t - 1)
end

# ╔═╡ e91b624c-cef9-4f74-8bf0-76a484c9c895
poly3_tone = poly3.(freqs)

# ╔═╡ 6880cf94-23ee-4a96-8b9f-0417fc239374
plot(t, poly3.(t))

# ╔═╡ 4caa26e8-c5ce-4d30-b80b-d3ed7659efa0
# Here there is no aliasing >-40db
plotfreqs(poly3_tone, "Poly3 Spectrum")

# ╔═╡ 1701a11b-b17b-499a-bf5f-0d721c026590
md"As we smooth the function more by setting higher derivitives equal
the polynomials approach a pure sine tone"

# ╔═╡ e060b886-4af2-4692-8882-79ceabe64ff1
function normalize(tone)
	""" Normalizing polynomial waveforms directly is annoying
	so we will use this function instead of calculating minimum and maximums"""
	m = minimum(tone)
	M = maximum(tone)
	2*(tone.-m)/(M-m) .-1
end

# ╔═╡ 85236c40-9c9c-4e00-b8e8-54787a5760af
function poly4(t, b)
	"""
	f = x^4 + bx^3 + cx^2 + dx
	1+b+c+d = f(1) = f(-1) = 1-b+c-d
	2(b+d) = 0, b+d = 0, d = -b
	-4+3b+2c-b = f'(1) = f'(-1) = -4+3b-2c-b
	2(4+2c) = 0, 4+2c = 0, c = -2
	f = x^4 + bx^3 - 2x^2 - bx
	"""
	t = t/pi%2-1
	t = (((t + b)*t -2)*t - b)*t
end

# ╔═╡ f1a60f5c-53c6-4533-b3ed-e715fb67d7ba
poly4_tone = normalize(poly4.(freqs, 0))

# ╔═╡ 788ec89e-c1d0-498b-8e8c-82de0bc4eda6
plot(t, normalize(poly4.(t, 0)))

# ╔═╡ e24388fa-d474-4770-95a0-d34627ae7a04
# No overtones reach the niquist rate, no aliasing at all
plotfreqs(poly4_tone, "Poly4 Spectrum")

# ╔═╡ bf89825a-944b-4771-a619-e4c5d1220f61
function poly5(t)
	"""
	f = ax^5 + bx^4 + cx^3 + dx^2 + ex + f
	f(-1) = 0 = -a + b -c + d + -e + f
	a+c+e = b+d+f
	f(1) = 0 = a+b+c+d+e+f = 2*(a+c+e) => a+c+e = 0
	5a-4b+3c-2d+e = f'(-1) = f'(1) = 5a+4b+3c+2d+e
	<=> 4b + 2d = 0 <=> d = -2b
	b - 2b + f = 0
	f = b
	20a+12b+6c+2d = f''(1) = f''(-1) = -20a +12b -6c +2d
	40a + 12 c = 0
	c = -10/3a
	60a+24b+6c = f'''(1) = f'''(-1) = 60a -24b + 6c
	b = 0
	"""
    t = t/pi%2-1
	a = 1
	c = -10/3*a
	b = 0
	((((a*t + b)*t + c)*t - 2*b)*t-(c + a))*t + b
end

# ╔═╡ 076a215e-1335-4974-b61e-06fa4c9e9295
begin
	plot(t, -normalize(poly5.(t)))
	plot!(t, sin.(t))
end

# ╔═╡ e5ac5d3e-1240-4a3c-bc9a-d8e8e6fec60f
poly5_tone = normalize(poly5.(freqs))

# ╔═╡ 82c3a459-0944-4637-a481-cafd15402afe
# We can see the progression of higher degree waveforms approaching the sine wave spectrum
plotfreqs(poly5_tone, "Poly5 Spectrum")

# ╔═╡ 36e284d1-a982-45c2-a0cc-dae1cc355362
md"## Odd harmonics
To get only odd harmonics
we need the second half of the waveform
to be the negative of the first half.
This sounds like we need an odd function where $f(-x) = -f(x)$,
but this would have the second half in reverse order.
We actually need a function where for some a, $f(x+a) = -f(x)$.
The only functions like that everywhere are already periodic,
since $f(x+2a) = -f(x+a) = f(x)$.
So we can't solve for polynomials like before.

However a waveform with this property can be constructed very easily piecewise
as $g(x) = \{f(x)$ if $x > 0$ else $-f(x+1)\}$ from any function defined on the interval $[0,1]$.
Using the sign of $x$ we can remove the if statement by taking
$g(x) = \text{sign}(x)f(2x-\text{sign}(x))$ on the interval $[-1, 1]$,
for a function $f(x)$ on $[-1, 1]$.

For the resulting function to be connected we need $f(-1) = f(1) = 0$.
For this function to be smooth we need $f'(-1) = -f'(1)$ as well.
Any even function ($f(x) = f(-x)$ everywhere) has both these properties,
provided that it is smooth on the interval and it is translated so that it's
endpoints are $0$.
"

# ╔═╡ 3e73b73b-79a3-4ba6-bfc5-94f39816680f
function para(t)
	""" f(x) = x^2-1
	g(x) = sign(x)f(2x-sign(x))
	= sign(x)((2x-sign(x))^2-1) = sign(x)((4x^2 -4xsign(x) + 1)-1)
	= 4*x(abs(x)-1)
	"""
	t = t/pi%2-1
	4*t*(abs(t)-1)
end

# ╔═╡ 97ec7198-36e8-410b-9602-d6b61f31db6b
begin
	plot(t, para.(t))
	plot!(t, sin.(t))
end

# ╔═╡ 2765c9be-3636-45de-9804-7fecb69f4a7b
para_tone = para.(freqs)

# ╔═╡ 52f4cba9-6eb4-42fe-bef5-8ea1a0f4e18a
plotfreqs(para_tone, "Parabola Spectrum")

# ╔═╡ 0ade700b-6193-4518-8470-bb787f87ba47
function ecubic(t)
	"""
	f(x) = abs(x)^3-1
	g(x) = sign(x)f(2x-sign(x))
         = sign(x)(abs(2x-sign(x))^3-1)
	"""
	t = t/pi%2-1
	s = sign(t)
	t = s*(abs(2*t-s)^3-1)
end

# ╔═╡ f116f443-0643-4cbb-8da8-4f9a312fd82f
plot(t, ecubic)

# ╔═╡ 35acdceb-34f5-4a39-aa06-569ab6f16205
ecubic_tone = ecubic.(freqs)

# ╔═╡ d7f37c37-a912-4324-8ca2-063f7c1f3ebc
plotfreqs(ecubic_tone, "Even Cubic Spectrum")

# ╔═╡ f86cda73-edd8-4905-81e0-83a6529daacc
function quart(t)
	"""
	f(x) = x^4-1
	g(x) = sign(x)((2x-sign(x))^4-1)
	"""
	t = t/pi%2-1
	s = sign(t)
	t = s*((2*t-s)^4-1)
end

# ╔═╡ 89fee0e4-7e7f-400b-9cc2-8b3186ab0384
plot(t, quart.(t))

# ╔═╡ b649ead1-acf8-4390-bb9e-ae876a1c7f48
quart_tone = quart.(freqs)

# ╔═╡ c492d40f-fc0d-4b06-a925-7c75bff32503
plotfreqs(quart_tone, "Quartic Spectrum")

# ╔═╡ d5e54701-36a3-45f9-ac2c-8917b2e3a09d
md"Polynomial methods are cheap to compute and have more tone than a pure sine wave but aren't inherently bandlimited and so can't make brighter sounds similar to analog square or saw waves without aliasing."

# ╔═╡ Cell order:
# ╟─f8ba933c-a1e3-11eb-2263-17f14870c1fa
# ╠═06240e71-df45-4f63-b797-4042432e6716
# ╠═bfacded1-6007-4e18-aba1-1b7325878dd2
# ╠═65dbcdfe-2dbd-4b28-ab59-fb12425913ca
# ╟─4bd75a55-4457-4625-a184-a2f7490cbf0e
# ╠═664b9502-13de-4fdf-971c-66cf0ea47c3b
# ╠═5e805acc-c6b1-4de5-8bc1-5bdcbd84ca04
# ╠═2a51d397-322b-414b-9d9f-a1e8ba6b2364
# ╠═dd813478-66a1-4d3e-9872-322e7b521e11
# ╠═a2357e48-116a-4df1-a52d-100533bad25d
# ╠═50824ffc-622d-4a00-acce-558cecf79e55
# ╠═3ae066cb-b5f2-4913-8a73-20348a401026
# ╠═7bd89cfd-9894-4644-b0e4-cae8a63cda9a
# ╠═dfa8eb03-2955-4bb6-b9ca-eb7cd6583e57
# ╠═a9f4a512-c8b4-474a-a93c-2cef99db534c
# ╠═f6f28913-cbfa-48f8-98ab-0d5ed267b350
# ╠═231185ea-8913-490d-8fe8-261983552457
# ╠═b627d7af-55d1-4d97-8f15-7f0d172f9b5a
# ╟─99708158-d74e-4311-a10f-d0ac4f79d4f2
# ╠═1f127661-47b0-4498-9fc2-857f314b6869
# ╠═99f21cb1-b50f-4e8b-9d97-8bf5297afbec
# ╠═8802e5fe-8416-4c67-a1f0-b99bb6af3e3f
# ╠═b561d76e-e2a7-4c3d-9677-42f8b5661b7c
# ╠═09d7e4f3-fed8-4fc8-b598-03ecc8e75d6f
# ╠═aff77549-c6da-4329-8e92-5f1fb561ece6
# ╠═e91b624c-cef9-4f74-8bf0-76a484c9c895
# ╠═6880cf94-23ee-4a96-8b9f-0417fc239374
# ╠═4caa26e8-c5ce-4d30-b80b-d3ed7659efa0
# ╠═1701a11b-b17b-499a-bf5f-0d721c026590
# ╠═e060b886-4af2-4692-8882-79ceabe64ff1
# ╠═85236c40-9c9c-4e00-b8e8-54787a5760af
# ╠═f1a60f5c-53c6-4533-b3ed-e715fb67d7ba
# ╠═788ec89e-c1d0-498b-8e8c-82de0bc4eda6
# ╠═e24388fa-d474-4770-95a0-d34627ae7a04
# ╠═bf89825a-944b-4771-a619-e4c5d1220f61
# ╠═076a215e-1335-4974-b61e-06fa4c9e9295
# ╠═e5ac5d3e-1240-4a3c-bc9a-d8e8e6fec60f
# ╠═82c3a459-0944-4637-a481-cafd15402afe
# ╠═36e284d1-a982-45c2-a0cc-dae1cc355362
# ╠═3e73b73b-79a3-4ba6-bfc5-94f39816680f
# ╠═97ec7198-36e8-410b-9602-d6b61f31db6b
# ╠═2765c9be-3636-45de-9804-7fecb69f4a7b
# ╠═52f4cba9-6eb4-42fe-bef5-8ea1a0f4e18a
# ╠═0ade700b-6193-4518-8470-bb787f87ba47
# ╠═f116f443-0643-4cbb-8da8-4f9a312fd82f
# ╠═35acdceb-34f5-4a39-aa06-569ab6f16205
# ╠═d7f37c37-a912-4324-8ca2-063f7c1f3ebc
# ╠═f86cda73-edd8-4905-81e0-83a6529daacc
# ╠═89fee0e4-7e7f-400b-9cc2-8b3186ab0384
# ╠═b649ead1-acf8-4390-bb9e-ae876a1c7f48
# ╠═c492d40f-fc0d-4b06-a925-7c75bff32503
# ╟─d5e54701-36a3-45f9-ac2c-8917b2e3a09d

max_stars_repo_name: jnha/music

max_issues_repo_name: jnha/music

max_forks_repo_name: jnha/music


text:
using IndirectShooting
using Test

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


text:
"""
A `DerivedDict` is a dictionary that inherits most of its behaviour from
an underlying dictionary.

The abstract type for derived dictionaries implements the interface of a
dictionary using composition and delegation. Concrete derived dictionaries
may override functions to specialize behaviour. For example, a mapped dictionary
may override the evaluation routine to apply the map first.
"""
abstract type DerivedDict{S,T} <: Dictionary{S,T} end

abstract type HighlySimilarDerivedDict{S,T} <: DerivedDict{S,T} end

###########################################################################
# Warning: derived sets implements all functionality by delegating to the
# underlying set, as if the derived set does not want to change any of that
# behaviour. This may result in incorrect defaults. The concrete set should
# override any functionality that is changed by it.
# Update: more functionality is copied for HighlySimilarDerivedDict, less
# for the derived dictionaries.
###########################################################################

# Assume the concrete set has a field called superdict -- override if it doesn't
superdict(d::DerivedDict) = d.superdict

# The concrete subset should implement similardictionary, as follows:
#
# similardictionary(d::ConcreteDerivedDict, d2::Dictionary) = ConcreteDerivedDict(d2)
#
# This function calls the constructor of the concrete set. We can then
# generically implement other methods that would otherwise call a constructor,
# such as resize and promote_eltype.

similar(d::DerivedDict, ::Type{T}, dims::Int...) where {T} =
    similardictionary(d, similar(superdict(d), T, dims))

resize(d::DerivedDict, dims) = similardictionary(d, resize(superdict(d), dims))

for op in (:coefficienttype,)
    @eval $op(d::DerivedDict) = $op(superdict(d))
end

# Delegation of properties
for op in (:isreal, :isbasis, :isframe, :isdiscrete, :measure)
    @eval $op(d::DerivedDict) = $op(superdict(d))
end

for op in (:isorthogonal, :isbiorthogonal, :isorthonormal)
    @eval $op(d::DerivedDict, μ::Measure) = $op(superdict(d), μ)
end

# Delegation of feature methods
for op in (:hasderivative, :hasantiderivative, :hasinterpolationgrid, :hasextension)
    @eval $op(Φ::DerivedDict) = $op(superdict(Φ))
end
for op in (:hasderivative, :hasantiderivative)
    @eval $op(Φ::DerivedDict, order) = $op(superdict(Φ), order)
end

# hastransform has extra arguments
hasgrid_transform(d::DerivedDict, gb, grid) = hasgrid_transform(superdict(d), gb, grid)

# When getting started with a discrete dictionary, you may want to write:
# hasderivative(d::ConcreteDict) = false
# hasantiderivative(d::ConcreteDict) = false
# hasinterpolationgrid(d::ConcreteDict) = false
# hastransform(d::ConcreteDict) = false
# hastransform(d::ConcreteDict, gb::GridBasis) = false
# hasextension(d::ConcreteDict) = false
# ... and then implement those operations one by one and remove the definitions.

zeros(::Type{T}, d::DerivedDict) where {T} = zeros(T, superdict(d))
tocoefficientformat(a, d::DerivedDict) = tocoefficientformat(a, superdict(d))

# Delegation of methods
for op in (:length, :extensionsize, :size, :interpolation_grid,
        :components, :tail, :ordering, :support, :dimensions)
    @eval $op(d::DerivedDict) = $op(superdict(d))
end

# Delegation of methods with an index parameter
for op in (:size, :element, :support)
    @eval $op(d::DerivedDict, i) = $op(superdict(d), i)
end

approx_length(d::DerivedDict, n::Int) = approx_length(superdict(d), n)

mapped_dict(d::HighlySimilarDerivedDict, map) = similardictionary(d, mapped_dict(superdict(d), map))

dict_in_support(d::DerivedDict, i, x) = in_support(superdict(d), i, x)

#########################
# Indexing and iteration
#########################

IndexStyle(dict::DerivedDict) = IndexStyle(superdict(dict))

native_index(dict::DerivedDict, idx) = native_index(superdict(dict), idx)

linear_index(dict::DerivedDict, idxn) = linear_index(superdict(dict), idxn)

eachindex(dict::DerivedDict) = eachindex(superdict(dict))

linearize_coefficients!(d::DerivedDict, coef_linear::Vector, coef_native) =
    linearize_coefficients!(superdict(d), coef_linear, coef_native)

delinearize_coefficients!(d::DerivedDict, coef_native, coef_linear::Vector) =
    delinearize_coefficients!(superdict(d), coef_native, coef_linear)

approximate_native_size(d::DerivedDict, size_l) = approximate_native_size(superdict(d), size_l)

linear_size(d::DerivedDict, size_n) = linear_size(superdict(d), size_n)

unsafe_eval_element(d::DerivedDict, idx, x) = unsafe_eval_element(superdict(d), idx, x)

unsafe_eval_element_derivative(d::DerivedDict, idx, x, order) =
    unsafe_eval_element_derivative(superdict(d), idx, x, order)

#########################
# Wrapping of operators
#########################

for op in (:transform_dict,)
    @eval $op(d::DerivedDict; options...) = $op(superdict(d); options...)
end

for op in (:derivative_dict, :antiderivative_dict)
    @eval $op(Φ::DerivedDict, order; options...) =
        similardictionary(Φ, $op(superdict(Φ), order; options...))
end

for op in (:extension, :restriction)
    @eval $op(::Type{T}, src::DerivedDict, dest::DerivedDict; options...) where {T} =
        wrap_operator(src, dest, $op(T, superdict(src), superdict(dest); options...))
end

function transform_from_grid(T, d1::GridBasis, d2::DerivedDict, grid; options...)
    op = transform_from_grid(T, d1, superdict(d2), grid; options...)
    wrap_operator(d1, d2, op)
end

function transform_to_grid(T, d1::DerivedDict, d2::GridBasis, grid; options...)
    op = transform_to_grid(T, superdict(d1), d2, grid; options...)
    wrap_operator(d1, d2, op)
end

for op in (:differentiation, :antidifferentiation)
    @eval $op(::Type{T}, d1::DerivedDict, d2::DerivedDict, order; options...) where {T} =
        wrap_operator(d1, d2, $op(T, superdict(d1), superdict(d2), order; options...))
end

function evaluation(::Type{T}, dict::DerivedDict, gb::GridBasis, grid; options...) where {T}
    A = evaluation(T, superdict(dict), gb, grid; options...)
    wrap_operator(dict, gb, A)
end

## Printing

show(io::IO, mime::MIME"text/plain", d::DerivedDict) = composite_show(io, mime, d)

Display.object_parentheses(d::DerivedDict) = false
Display.stencil_parentheses(d::DerivedDict) = false
Display.displaystencil(d::DerivedDict) = [modifiersymbol(d), "(", superdict(d), ")"]

#########################
# Concrete dict
#########################

"""
For testing purposes we define a concrete subset of DerivedDict. This set should
pass all interface tests and be functionally equivalent to the underlying set.
"""
struct ConcreteDerivedDict{S,T} <: HighlySimilarDerivedDict{S,T}
    superdict   ::  Dictionary{S,T}
end

# Implementing similardictionary is all it takes.

similardictionary(dict::ConcreteDerivedDict, dict2::Dictionary) = ConcreteDerivedDict(dict2)


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %%
# Load the data, don't forget to drop everything missing.
d = DataFrame(CSV.File(datadir("exp_raw/cherry_blossoms.csv"), missingstring = "NA"))
d2 = dropmissing(d, :doy)

# %%
using BSplines: BSplineBasis, basismatrix

num_knots = 15
knot_list = quantile(d2.year, range(0, 1, length = num_knots))

basis = BSplineBasis(4, knot_list)
B = basismatrix(basis, d2.year)

@model function spline(doy, B = B)
    # α ~ Normal(100, 10)
    # w ~ filldist(Normal(0, 10), size(B, 2))
    w ~ filldist(Normal(100, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = B * w
    doy ~ MvNormal(μ, σ)
end

qspl = quap(spline(d2.doy))

# %%
w_str = ["w[$i]" for i in 1:length(basis)]
post = DataFrame(rand(qspl.distr, 1000)', [w_str; "σ"])

smu = B * Array(post[!, w_str])'
smu_m, smu_lower, smu_upper = meanlowerupper(smu)

scatter(d2.year, d.doy, alpha = 0.3)
plot!(d2.year, smu_m, ribbon = (smu_m .- smu_lower, smu_upper .- smu_m))

# %%
"""
Not sure what he means by "The first basis functions could substitute for the intercept".
If I'm not mistaken I need to add the intercept to all of the basis functions since the
first one doesn't reach to the end?
"""


text:
export
    DisplayObj,
    Label,
    Panel

"""
DisplayObj

the inheritance of DisplayObj must have the following fields:
* `content::DisplayRow` or `content::DisplayArray`
* `style::Crayon`
* `pos::Vector{Int}` or `pos::Tuple{Int, Int}`
"""
abstract type DisplayObj end

function render(stream::IO, display_obj::DisplayObj)
    T.print(stream, display_obj.style);
    render(stream, display_obj.content, pos=display_obj.pos);
    T.print(stream, string(Crayon(reset=true)))
end

render(display_obj::DisplayObj) = render(T.out_stream, display_obj)

struct Label <: DisplayObj
    content::DisplayRow
    style::Crayon
    pos::Vector{Int}
end

function Label(
    str::String,
    pos::Vector{Int};
    style::Crayon=Crayon(foreground=:green)
)
    content = DisplayRow(textwidth(str))
    content[1:end] = str

    return Label(content, style, pos)
end

struct Panel <: DisplayObj
    content::DisplayArray
    style::Crayon
    pos::Vector{Int}
end

function Panel(
    str::String,
    size::Vector{Int},
    pos::Vector{Int};
    style::Crayon=Crayon(foreground=:blue, bold=true), kwargs...
)

    h, w = size
    content = DisplayArray(Rectangle(h, w;kwargs...))
    str_width = textwidth(str)
    (w-str_width > 2) && (content[1, 2:(str_width+1)] = str)
    return Panel(content, style, pos)
end


text:
import BasisFunctions: approximate
export Fun

###################
# Helper functions
###################

determine_return_type(fun, S) = Base.promote_op(fun, S)
determine_return_type(fun, S::Type{Tuple{A}}) where {A} = Base.promote_op(fun, A)
determine_return_type(fun, S::Type{Tuple{A,B}}) where {A,B} = Base.promote_op(fun, A, B)
determine_return_type(fun, S::Type{Tuple{A,B,C}}) where {A,B,C} = Base.promote_op(fun, A, B, C)
determine_return_type(fun, S::Type{SVector{1,A}}) where {A} = Base.promote_op(fun, A)
determine_return_type(fun, S::Type{SVector{2,A}}) where {A} = Base.promote_op(fun, A, A)
determine_return_type(fun, S::Type{SVector{3,A}}) where {A} = Base.promote_op(fun, A, A, A)
determine_return_type(fun, S::Type{SVector{4,A}}) where {A} = Base.promote_op(fun, A, A, A, A)

promote_dictionary(dict, fun) = _promote_dictionary(dict, fun, determine_return_type(fun, domaintype(dict)))
_promote_dictionary(dict, fun, ::Union{}) = dict
_promote_dictionary(dict, fun, ::Type{T}) where {T} = BasisFunctions.ensure_coefficienttype(promote_type(codomaintype(dict),T), dict)

function directsolver(A::DiagonalOperator; verbose = false, options...)
    verbose && println("Using direct solver: inverse of diagonal matrix")
    inv(A)
end

function directsolver(A; directsolver = :svd, verbose = false, options...)
    if directsolver == :svd
        verbose && println("Using direct solver: SVD")
        SVD_solver(A)
    elseif directsolver == :qr
        verbose && println("Using direct solver: QR")
        QR_solver(A)
	elseif directsolver == :regsvd
        verbose && println("Using direct solver: regularized SVD")
        regularized_SVD_solver(A; verbose=verbose, options...)
    else
        DS = directsolver(A; verbose=verbose, options...)
        verbose && println("Using direct solver: $DS")
        DS
    end
end

function iterativesolver(A; iterativesolver = :lsqr, verbose = false, options...)
    if iterativesolver == :lsqr
        verbose && println("Using iterative solver: LSQR")
        LSQR_solver(A; verbose=verbose, options...)
    elseif iterativesolver == :lsmr
        verbose && println("Using iterative solver: LSMR")
        LSMR_solver(A; verbose=verbose, options...)
    else
        DS = iterativesolver(A)
        verbose && println("Using direct solver: $DS")
        DS
    end
end

####################
# The Fun interface
####################

guess_coefficienttype(dict, fun) = promote_type(codomaintype(dict), determine_return_type(fun, domaintype(dict)))

"""
`Fun` is used to approximate functions in the same way as `approximate` with
the same interface. It discards all outputs of `approximate` and only
returns the function approximation.
"""
function Fun(fun, dict::Dictionary, args...;
        coefficienttype = guess_coefficienttype(dict, fun),
        options...)
	dict2 = BasisFunctions.ensure_coefficienttype(coefficienttype, dict)
    ap = approximationproblem(fun, dict2, args...; options...)
    Fun(ap; options...)
end

function Fun(fun, platform::Platform, args...; options...)
    ap = approximationproblem(fun, platform, args...)
    Fun(ap; options...)
end

Fun(fun, ap::ApproximationProblem; options...) = Fun(_ap(fun, ap); options...)

function Fun(ap::ApproximationProblem; verbose = false, options...)
    if verbose
        println("Fun: using the following dictionary:")
        show(dictionary(ap))
		println()
    end
    F, A, B, C, S = approximate(ap; verbose=verbose, options...)
    F
end

function Fun(ap::AdaptiveApproximation; verbose = false, options...)
    verbose && println("Fun: adaptive approximation using platform $(platform(ap))")
    approximate(ap; verbose=verbose, options...)
end

function Fun(dict::Dictionary, coefficients::AbstractArray; verbose = false, options...)
	verbose && println("Fun: coefficients are given, no approximation problem is solved")
	Expansion(dict, coefficients)
end

############################
# The approximate function
############################

approximate(fun, dict::Dictionary, args...; options...) =
    approximate(_ap(fun, dict, args...); options...)
approximate(fun, platform::Platform, args...; options...) =
    approximate(_ap(fun, platform, args...); options...)

approximate(ap::ApproximationProblem;
            samplingstyle = SamplingStyle(ap),
            solverstyle = SolverStyle(ap), options...) =
    approximate(promote(samplingstyle, solverstyle)..., ap; options...)

# If the sampling has product structure and a single solverstyle is specified,
# we turn the solverstyle into a product style as well.
promote(samplingstyle::SamplingStyle, solverstyle::SolverStyle) = (samplingstyle,solverstyle)
promote(samplingstyle::ProductSamplingStyle, solverstyle::ProductSolverStyle) = (samplingstyle,solverstyle)
promote(samplingstyle::ProductSamplingStyle, solverstyle::SolverStyle) =
	(samplingstyle,ProductSolverStyle(map(x->solverstyle, samplingstyle.styles)))

# Construct the approximation problem and solve it
function approximate(samplingstyle::SamplingStyle, solverstyle::SolverStyle, ap::ApproximationProblem;
            verbose = false, options...)
    dict = dictionary(ap)

    verbose && println("Approximate: using sampling style $samplingstyle")
	verbose && println("Approximate: using solver style $solverstyle")

    S = sampling_operator(samplingstyle, ap; verbose=verbose, options...)
	if verbose
		println()
		println("Approximate: sampling operator with size $(size(S)) is")
		println(S)
		println()
	end
	A = discretization(ap; options...)
	b = sample_data(ap; options...)
    coef = solve(solverstyle, ap, A, b; S=S, verbose=verbose, samplingstyle=samplingstyle, options...)
    F = Expansion(dictionary(ap), coef)
    res = norm(A*coef-b)
    verbose && println("Approximate: ended with residual $res\n")
    F, A, b, coef, S, samplingparameter(ap)
end

export solve
solve(style::SolverStyle, ap, A::DictionaryOperator, B; options...) =
    solver(style, ap, A; B=B, options...) * B
solve(style::SolverStyle, ap, A::AbstractOperator, fun; options...) =
    apply(solver(style, ap, A; options...), fun; options...)


text:
export SymbolicProbabilityEstimator

"""
A probability estimator based on symbolization.
"""
abstract type SymbolicProbabilityEstimator <: ProbabilitiesEstimator end

include("utils.jl")
include("SymbolicPermutation.jl")
include("SymbolicWeightedPermutation.jl")
include("SymbolicAmplitudeAware.jl")

"""
    permentropy(x; τ = 1, m = 3, base = MathConstants.e)
Shorthand for `genentropy(x, SymbolicPermutation(m = m, τ = τ); base)` which
calculates the permutation entropy of order `m` with delay time `τ` used for embedding.
"""
function permentropy(x; τ = 1, m = 3, base = MathConstants.e)
    return genentropy(x, SymbolicPermutation(m = m, τ = τ); base)
end

export permentropy


text:
export UserKNN

"""
    UserKNN(
        data::DataAccessor,
        k::Int,
        normalize::Bool=false
    )

[User-based CF using the Pearson correlation](https://dl.acm.org/citation.cfm?id=312682). `k` represents number of neighbors, and `normalize` specifies if weighted sum of neighbors' rating is normalized.

The technique gives a weight to a user-user pair by the following equation:

```math
s_{u, v} = \\frac{ \\sum_{i \\in \\mathcal{I}^+_{u \\cap v}}  (r_{u, i} - \\overline{r}_u) (r_{v, i} - \\overline{r}_v)}
{ \\sqrt{\\sum_{i \\in \\mathcal{I}^+_{u \\cap v}} (r_{u,i} - \\overline{r}_u)^2} \\sqrt{\\sum_{i \\in \\mathcal{I}^+_{u \\cap v}} (r_{v, i} - \\overline{r}_v)^2} },
```

where ``\\mathcal{I}^+_{u \\cap v}`` is a set of items which were observed by both user ``u`` and ``v``, and ``r_{u, i}`` indicates a ``(u, i)`` element in ``R``. Plus, ``\\overline{r}_u`` and ``\\overline{r}_v`` are respectively mean values of ``r_{u, i}`` and ``r_{v, i}`` over ``i \\in \\mathcal{I}^+_{u \\cap v}``. By using the weights, user-based CF selects the top-``k`` highest-weighted users (i.e., nearest neighbors) of a target user ``u``, and predicts missing ``r_{u, i}`` for ``i \\in \\mathcal{I}^-_u`` as:

```math
r_{u, i} = \\overline{r}_u + \\frac{\\sum^k_{t=1} \\left(r_{\\sigma(t), i} - \\overline{r}_{\\sigma(t)}\\right) \\cdot s_{u,\\sigma(t)} }{ \\sum^k_{t=1} s_{u,\\sigma(t)} },
```

where ``\\sigma(t)`` denotes the ``t``-th nearest-neighborhood user. Ultimately, sorting items ``\\mathcal{I}^-_u`` by the predicted values enables us to make recommendation to a user ``u``.

It should be noted that user-based CF is highly inefficient because gradually increasing massive users and their dynamic tastes require us to frequently recompute the similarities for every pairs of users.
"""
struct UserKNN <: Recommender
    data::DataAccessor
    k::Int
    sim::AbstractMatrix
    normalize::Bool

    function UserKNN(data::DataAccessor, k::Int, normalize::Bool)
        n_user = size(data.R, 1)
        new(data, k, matrix(n_user, n_user), normalize)
    end
end

UserKNN(data::DataAccessor, k::Int) = UserKNN(data, k, false)
UserKNN(data::DataAccessor) = UserKNN(data, 20, false)

isbuilt(recommender::UserKNN) = isfilled(recommender.sim)

function build!(recommender::UserKNN)
    # Pearson correlation

    R = copy(recommender.data.R)

    n_row = size(R, 1)

    for ri in 1:n_row
        for rj in ri:n_row
            # pairwise correlation (i.e., ignore NaNs)
            ij = broadcast(!isalmostzero, R[ri, :]) .& broadcast(!isalmostzero, R[rj, :])

            vi = R[ri, :] .- mean(R[ri, ij])
            vj = R[rj, :] .- mean(R[rj, ij])

            numer = dot(vi[ij], vj[ij])
            denom = sqrt(dot(vi[ij], vi[ij]) * dot(vj[ij], vj[ij]))

            c = (denom == 0) ? 0 : (numer / denom)

            recommender.sim[ri, rj] = c
            if (ri != rj); recommender.sim[rj, ri] = c; end # symmetric
        end
    end
end

function predict(recommender::UserKNN, u::Int, i::Int)
    check_build_status(recommender)

    numer = denom = 0

    pairs = collect(zip(1:size(recommender.data.R)[1], recommender.sim[u, :]))
    # closest neighbor is always target user him/herself, so omit him/her
    ordered_pairs = sort(pairs, by=tuple->last(tuple), rev=true)[2:(recommender.k + 1)]

    for (u_near, w) in ordered_pairs
        v_near = recommender.data.R[u_near, :]

        r = v_near[i]
        if isalmostzero(r); continue; end

        r_ = 0
        if recommender.normalize
            jj = broadcast(!isalmostzero, v_near)
            r_ = mean(v_near[jj])
        end

        numer += (r - r_) * w
        denom += w
    end

    pred = (denom == 0) ? 0 : numer / denom
    if recommender.normalize
        ii = broadcast(!isalmostzero, recommender.data.R[u, :])
        m = mean(recommender.data.R[u, ii])
        pred += isnan(m) ? 0 : m
    end
    pred
end


text:
module Cbc

if isfile(joinpath(Pkg.dir("Cbc"),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Cbc not properly installed. Please run Pkg.build(\"Cbc\")")
end

include("CoinMPInterface.jl")
include("CbcSolverInterface.jl")

using Cbc.CbcMathProgSolverInterface
export CbcSolver

end # module


text:
@enum BracketState begin
    ISOTOPE = Int8(0)
    SYMBOL = Int8(1)
    MULTIPLICITY = Int8(2)
    CHARGE = Int8(3)
    FINISH = Int8(4)
end

struct BracketParseException <: Exception
    msg::String
    bracketcontents::String
    charloc::Int
end
Base.showerror(io::IO, e::BracketParseException) = print(io, "Failed to parse bracket(\"[$(e.bracketcontents)]\") on character ", e.charloc, "!\n", e.msg)

function ParseBracket(s::AbstractString)
    len = length(s)

    state = ISOTOPE
    symbol, isotope, aromatic, ringID, hydrogens, charge = nothing, nothing, false, [], 0, 0

    idx = 1
    while idx <= len
        cursor = s[idx]
        #Handle isotope, explicit hydrogen multiplicity
        if isdigit(cursor)
            if state == ISOTOPE
                isotope, idx = tryparseint16(s, idx)
                state = SYMBOL
            elseif state == MULTIPLICITY
                hydrogens, idx = tryparseint16(s, idx)
                state = CHARGE
            else
               throw(BracketParseException("Invalid SMILES. Charge must be denoted with +/-.", s, idx))
            end
        #Handle Element symbol parsing
        elseif isletter(cursor)
            if state == SYMBOL || state == ISOTOPE
                if isa(symbol, Nothing)
                    symbollist = isuppercase(cursor) ? bracket : bracket_aromatic
                    aromatic = islowercase(cursor)
                    symbol, idx = tryparseelement(s, idx, symbollist)
                    if isa(symbol, Nothing)
                        throw(BracketParseException("Invalid SMILES. Unrecognized chemical symbol.", s, idx))
                    end
                elseif cursor == 'H'
                    state = MULTIPLICITY
                    hydrogens = 1
                    idx = nextind(s, idx)
                else
                    throw(BracketParseException("Invalid SMILES. Bracket should contain one element (with or without Hydrogen).", s, idx))
                end
            else
                throw(BracketParseException("Invalid SMILES. Bracket should contain one element.", s, idx))
            end
        #Handle charge
        elseif ispm(cursor)
            if state == CHARGE || state == MULTIPLICITY || state == SYMBOL || state == ISOTOPE
                charge, idx = tryparsecharge(s, idx)
                state = SYMBOL
            else
                throw(BracketParseException("Invalid SMILES. Charge should follow element symbol.", s, idx))
            end
        #Handle operators
        elseif isspecialoperator(cursor)
            throw(BracketParseException("Invalid SMILES. Special operators (" * join(specialoperators, ",") * ") should not be contained in a bracket.", s, idx))
        elseif isbondoperator(cursor)
            if cursor == '@'
                @warn("Chirality is not yet supported, ignoring @ operator", maxlog=1)
                _, idx = tryparsechirality(s, idx)
            else
                throw(BracketParseException("Invalid SMILES. Bond operators(" * join(setdiff(bondoperators, '@'), ",") * ") should not be contained in a bracket.", s, idx))
            end
        end
    end
    return Element( symbol, isotope, aromatic, ringID, hydrogens, 0, charge )
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

export LAPACKException,
       SingularException,
       PosDefException,
       RankDeficientException,
       ZeroPivotException

struct LAPACKException <: Exception
    info::BlasInt
end

"""
    SingularException

Exception thrown when the input matrix has one or more zero-valued eigenvalues, and is not invertible.
A linear solve involving such a matrix cannot be computed.
The `info` field indicates the location of (one of) the singular value(s).
"""
struct SingularException <: Exception
    info::BlasInt
end

"""
    PosDefException

Exception thrown when the input matrix was not [positive definite](https://en.wikipedia.org/wiki/Definiteness_of_a_matrix).
Some linear algebra functions and factorizations are only applicable to positive definite matrices.
The `info` field indicates the location of (one of) the eigenvalue(s) which is (are) less than/equal to 0.
"""
struct PosDefException <: Exception
    info::BlasInt
end
function Base.showerror(io::IO, ex::PosDefException)
    print(io, "PosDefException: matrix is not ")
    if ex.info == -1
        print(io, "Hermitian")
    else
        print(io, "positive definite")
    end
    print(io, "; Cholesky factorization failed.")
end

struct RankDeficientException <: Exception
    info::BlasInt
end

"""
    ZeroPivotException <: Exception

Exception thrown when a matrix factorization/solve encounters a zero in a pivot (diagonal)
position and cannot proceed.  This may *not* mean that the matrix is singular:
it may be fruitful to switch to a diffent factorization such as pivoted LU
that can re-order variables to eliminate spurious zero pivots.
The `info` field indicates the location of (one of) the zero pivot(s).
"""
struct ZeroPivotException <: Exception
    info::BlasInt
end
function Base.showerror(io::IO, ex::ZeroPivotException)
    print(io, "ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.")
end

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
module MyApp

using Genie, Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = MyApp))

  Genie.genie(; context = @__MODULE__)

  Base.eval(Main, :(const Genie = MyApp.Genie))
  Base.eval(Main, :(using Genie))
end

end

max_stars_repo_name: dvasya/herocker

max_issues_repo_name: dvasya/herocker

max_forks_repo_name: dvasya/herocker


text:
p = 2
function pow_array(x::Vector{Float64})
s = 0.0
for y in x
s = s + y^p
end
return s
end

max_stars_repo_path: Module 2/Chapter04/pow_Array.jl

max_issues_repo_path: Module 2/Chapter04/pow_Array.jl

max_forks_repo_path: Module 2/Chapter04/pow_Array.jl


text:
module Foundation

include("extensions/array.jl")
include("extensions/string.jl")
include("extensions/math.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 9513dfd4-cb93-11ec-03d9-2b4109ca3b97
begin
	using Pkg
	Pkg.activate("../../")
	using PyPlotly
end

# ╔═╡ 170bc566-3b0a-4436-a249-2cb26d891e90
begin
	df = px.data.iris()
	fig = px.scatter(
		df,
		x="sepal_width",
		y="sepal_length",
		color="species",
		size="petal_length",
		hover_data=["petal_width"],
	)
	fig
end

# ╔═╡ Cell order:
# ╠═9513dfd4-cb93-11ec-03d9-2b4109ca3b97
# ╠═170bc566-3b0a-4436-a249-2cb26d891e90