text:
module Bio3DView

export
    Style,
    Surface,
    IsoSurface,
    Box,
    Axes,
    CameraAngle,
    viewfile,
    viewstring,
    viewstruc,
    viewpdb

using Requires

# Counter for viewers so they can be named individually
element_count = 0

function __init__()
    @require BioStructures="de9282ab-8554-53be-b2d6-f6c222edabfc" include("biostructures.jl")
    @require Blink="ad839575-38b3-5650-b840-f874b8c74a25" include("blink.jl")
end

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited

path_lib = normpath(@__DIR__, "..", "js")
path_3dmol = joinpath(path_lib, "3Dmol-nojquery-min.js")
path_jquery = joinpath(path_lib, "jquery-3.3.1.min.js")
js_3dmol = read(path_3dmol, String)
js_jquery = read(path_jquery, String)

"""
    Style(style_type)
    Style(style_type, options)

A style for a molecular visualisation, with an optional `Dict` of options.
Examples are `Style("cartoon")` and
`Style("cartoon", Dict("color"=> "spectrum", "ribbon"=> true, "thickness"=> 1.0))`.
"""
struct Style
    name::String
    options::Dict{String, Any}
end

Style(s::AbstractString) = Style(s, Dict())

# Default style for molecular visualisation, depends on file format
function defaultstyle(format::AbstractString)
    if format == "pdb"
        return Style("cartoon", Dict("color"=> "spectrum"))
    elseif format in ("sdf", "xyz", "mol2")
        return Style("stick")
    else
        throw(ArgumentError("Not a valid file format: \"$fmt\""))
    end
end

"""
    Surface()
    Surface(options)

A style for a molecular VDW surface visualisation, with an optional `Dict` of
options.
An example is `Surface(Dict("opacity"=> 0.8, "colorscheme"=> "whiteCarbon"))`.
"""
struct Surface
    options::Dict{String, Any}
end

Surface() = Surface(Dict())

"""
    IsoSurface(voldata, isoval)

Data and styling for an isosurface visualisation.
Arguments are the filepath with the volume data in "cube" format and the value
to view the isosurface at.
Optional keyword arguments are `color`, `opacity`, `wireframe` and `smoothness`.
"""
struct IsoSurface
    voldata::String
    isoval::Float64
    color::String
    opacity::Float64
    wireframe::Bool
    smoothness::Float64
end

function IsoSurface(voldata::AbstractString,
                isoval::Real;
                color::AbstractString="blue",
                opacity::Real=0.8,
                wireframe::Bool=false,
                smoothness::Real=5)
    return IsoSurface(voldata, isoval, color, opacity, wireframe, smoothness)
end

"""
    Box(center, dimensions)

Data and styling for a box visualisation.
Arguments are a `Vector{Float64}` of the center coordinates and a
`Vector{Float64}` of the box size in each dimension.
Optional keyword arguments are `color` and `wireframe`.
"""
struct Box
    center::Vector{Float64}
    dims::Vector{Float64}
    color::String
    wireframe::Bool
end

Box(center, dims; color::AbstractString="black", wireframe::Bool=true) = Box(
        center, dims, color, wireframe)

"""
    Axes(length, radius)
    Axes()

Styling for visualisation of the coordinate axes.
Arguments are numbers for the length and the radius of the arrows.
The optional keyword argument `colors` is a list of 3 colors as `String`s for
the 3 axes.
"""
struct Axes
    len::Float64
    radius::Float64
    colors::Vector{String}
end

function Axes(len::Real=1.0,
                radius::Real=0.1;
                colors::Vector{String}=["red", "green", "blue"])
    return Axes(len, radius, colors)
end

"""
    CameraAngle(posx, posy, posz, zoom, qx, qy, qz, qw)

A custom perspective to view the molecule from.
Arguments are x/y/z translation, zoom, and x/y/z/w rotation quaternion.
Default is `CameraAngle(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0)`.
"""
struct CameraAngle
    posx::Float64
    posy::Float64
    posz::Float64
    zoom::Float64
    qx::Float64
    qy::Float64
    qz::Float64
    qw::Float64
end

"""
    viewfile(file)
    viewfile(file, format)

View a molecular structure from a file.
Displays in a popup window, or in the output cell for an IJulia notebook.
Arguments are the filepath and the format ("pdb", "sdf", "xyz" or "mol2").
If not provided, the format is guessed from the file extension, e.g.
"myfile.xyz" is treated as being in the xyz format.
Optional keyword arguments are `style`, `surface`, `isosurface`, `box`,
`vtkcell`, `axes`, `cameraangle`, `height`, `width`, `html` and `debug`.
"""
function viewfile(f::AbstractString,
                format::AbstractString=lowercase(split(f, ".")[end]);
                style::Style=defaultstyle(format),
                kwargs...)
    if !isfile(f)
        throw(ArgumentError("Cannot find file \"$f\""))
    end
    return view("data-type='$format'", read(f, String); style=style, kwargs...)
end

"""
    viewstring(str, format)

View a molecular structure contained in a string.
Displays in a popup window, or in the output cell for an IJulia notebook.
Arguments are the molecule string and the format ("pdb", "sdf", "xyz" or
"mol2").
Optional keyword arguments are `style`, `surface`, `isosurface`, `box`,
`vtkcell`, `axes`, `cameraangle`, `height`, `width`, `html` and `debug`.
"""
function viewstring(s::AbstractString,
                format::AbstractString;
                style::Style=defaultstyle(format),
                kwargs...)
    return view("data-type='$format'", s; style=style, kwargs...)
end

"""
    viewpdb(pdbid)

View a structure from the Protein Data Bank (PDB).
Displays in a popup window, or in the output cell for an IJulia notebook.
Argument is the four letter PDB ID, e.g. "1AKE".
Requires an internet connection to work.
Optional keyword arguments are `style`, `surface`, `isosurface`, `box`,
`vtkcell`, `axes`, `cameraangle`, `height`, `width`, `html` and `debug`.
"""
function viewpdb(p::AbstractString;
                style::Style=defaultstyle("pdb"),
                kwargs...)
    if !occursin(r"^[a-zA-Z0-9]{4}$", p)
        throw(ArgumentError("Not a valid PDB ID: \"$p\""))
    end
    return view("data-pdb='$p'"; style=style, kwargs...)
end

# Get the script to add an isosurface from an IsoSurface object
function isosurfacestring(iso::IsoSurface)
    return "data = `$(read(iso.voldata, String))`\n" *
        "var voldata = new \$3Dmol.VolumeData(data, \"cube\");\n" *
        "viewer.addIsosurface(voldata, {isoval: $(iso.isoval), " *
        "color: \"$(iso.color)\", alpha: $(iso.opacity), " *
        "wireframe: $(iso.wireframe), smoothness: $(iso.smoothness)});\n"
end

# Get the script to add a box from a Box object
function boxstring(box::Box)
    return "viewer.addBox({" *
        "center:{x:$(box.center[1]),y:$(box.center[2]),z:$(box.center[3])}, " *
        "dimensions:{w:$(box.dims[1]),h:$(box.dims[2]),d:$(box.dims[3])}, " *
        "color:'$(box.color)', " *
        "wireframe:$(box.wireframe)});\n"
end

# Get the script to add a unit cell from a vtk format file
function vtkcellstring(f::AbstractString)
    coords = []
    lines = []
    mode = "none"
    open(f) do file
        for l in eachline(file)
            if startswith(l, "POINTS")
                mode = "points"
            elseif startswith(l, "LINES")
                mode = "lines"
            elseif mode == "points"
                push!(coords, parse.(Float64, split(l)))
            elseif mode == "lines"
                # Account for zero-based indexing
                push!(lines, parse.(Int, split(l))[2:end] .+ 1)
            end
        end
    end
    o = ""
    for (i, j) in lines
        o *= "viewer.addLine({\n" *
        "start:{x:$(coords[i][1]), y:$(coords[i][2]), z:$(coords[i][3])},\n" *
        "end:{x:$(coords[j][1]), y:$(coords[j][2]), z:$(coords[j][3])}\n});\n"
    end
    return o
end

# Get the script to add arrows showing the coordinate axes
function axesstring(a::Axes)
    o = ""
    for (x, y, z, col) in ((a.len, 0, 0, a.colors[1]), (0, a.len, 0, a.colors[2]),
                (0, 0, a.len, a.colors[3]))
        o *= "viewer.addArrow({\nstart:{x:0, y:0, z:0},\n" *
                "end:{x:$x, y:$y, z:$z},\ncolor:'$col',\nradius:$(a.radius)\n});\n"
    end
    return o
end

# Get the script to change the view to a custom angle
function cameraanglestring(a::CameraAngle)
    return "viewer.setView([$(a.posx), $(a.posy), $(a.posz), $(a.zoom), " *
            "$(a.qx), $(a.qy), $(a.qz), $(a.qw)]);\n"
end

# Generate HTML to view a molecule
function view(tag_str::AbstractString,
                data_str::AbstractString="";
                style::Style,
                surface::Union{Surface, Nothing}=nothing,
                isosurface::Union{IsoSurface, Nothing}=nothing,
                box::Union{Box, Nothing}=nothing,
                vtkcell::Union{AbstractString, Nothing}=nothing,
                axes::Union{Axes, Nothing}=nothing,
                cameraangle::Union{CameraAngle, Nothing}=nothing,
                height::Integer=540,
                width::Integer=540,
                html::Bool=false,
                debug::Bool=false)
    global element_count
    element_count += 1
    if length(data_str) > 0
        data_id = "3dmol_data_$element_count"
        data_div = "<textarea style='display:none;' id='$data_id'>$data_str</textarea>"
        tag_str *= " data-element='$data_id'"
    else
        data_div = ""
    end
    if surface != nothing
        surface_tag = "data-surface='$(tagstring(surface))'"
    else
        surface_tag = ""
    end
    viewer_id = "3dmol_viewer_$element_count"
    div_str = "

"
    script_str = ""
    ijulia_html = "\n" *
            "\n$data_div\n$div_str\n$script_str\n"
    # Return stand-alone HTML only
    if html
        return "\n<meta charset=\"UTF-8\">\n\n" *
                "$ijulia_html\n"
    end
    if isijulia()
        return HTML(ijulia_html)
    else
        if !isdefined(Main, :Blink)
            throw("Cannot run the command as you do not appear to be in an IJulia " *
                    "notebook, Blink is not loaded and html is not set to true")
        end
        if Sys.iswindows()
            req_path = replace(path_jquery, "\\" => "\\\\")
        else
            req_path = path_jquery
        end
        # The first part gets jQuery to work with Electron
        blink_html = "\n" *
            "\n" *
            "\n$data_div\n$div_str\n$script_str\n"
        viewblink(blink_html, height, width, debug)
    end
end

# Convert Style instance to a tag string
# Style format is "cartoon:color=red,ribbon=true"
function tagstring(s::Style)
    o = s.name
    if length(s.options) > 0
        o *= ":"
        for k in keys(s.options)
            o *= "$k=$(string(s.options[k])),"
        end
        # Strip trailing comma
        o = o[1:end - 1]
    end
    return o
end

# Convert Surface instance to a tag string
# Surface format is "opacity:0.7;color:white"
function tagstring(s::Surface)
    o = ""
    if length(s.options) > 0
        for k in keys(s.options)
            o *= "$k:$(string(s.options[k]));"
        end
        # Strip trailing comma
        o = o[1:end - 1]
    end
    return o
end

end


text:
@enum PairedReadOrientation::UInt64 FwRv=1 RvFw=2

struct PairedReads <: ReadDatastore{LongSequence{DNAAlphabet{4}}}
    filename::String         # Filename datastore was opened from.
    name::String             # Name of the datastore. Useful for other applications.
    defaultname::String      # Default name, useful for other applications.
    readsize::UInt64         # Maximum size of any read in this datastore.
    chunksize::UInt64        # Number of chunks of sequence data per read.
    fragsize::UInt64         # Fragment size of library.
    readpos_offset::UInt64
    size::UInt64
    orientation::PairedReadOrientation
    stream::IOStream
end

@inline orientation(prds::PairedReads) = prds.orientation
@inline maxseqlen(prds::PairedReads) = prds.readsize
@inline name(prds::PairedReads) = prds.name
@inline stream(prds::PairedReads) = prds.stream

###
### PairedReads Header format
###

# | Field             | Value  | Type        |
# |:-----------------:|:------:|:-----------:|
# | Magic number      | 0x05D5 | UInt16      | 2
# | Datastore type    | 0x0001 | UInt16      | 2
# | Version number    | 0x0001 | UInt16      | 2
# | Default name      | N/A    | String      | 7
# | Maximum read size | N/A    | UInt64      | 8
# | Chunk size        | N/A    | UInt64      | 8
# | Fragment size     | N/A    | UInt64      | 8
# | Orientation       | N/A    | Orientation | 8
# | Number of Reads   | N/A    | UInt64      | 8

const PairedDS_Version = 0x0001

function PairedReads(rdrx::FASTQ.Reader, rdry::FASTQ.Reader,
                     outfile::String, name::String,
                     minsize::UInt64, maxsize::UInt64,
                     fragsize::UInt64, orientation::PairedReadOrientation)

    # Create and allocate the sequence and record objects.
    lread = FASTQ.Record()
    rread = FASTQ.Record()
    lseq = LongSequence{DNAAlphabet{4}}(maxsize)
    rseq = LongSequence{DNAAlphabet{4}}(maxsize)

    #chunksize::UInt64 = BioSequences.seq_data_len(DNAAlphabet{4}, maxsize)
    chunksize::UInt64 = length(BioSequences.encoded_data(lseq))

    fd = open(outfile, "w")

    # Write magic no, datastore type, version number.
    sizepos = write(fd, ReadDatastoreMAGIC, PairedDS, PairedDS_Version) +
    # Write the default name of the datastore.
    writestring(fd, name) +
    # Write the read size, and chunk size.
    write(fd, maxsize, chunksize, fragsize, orientation)
    # Write space for size variable (or number of read pairs).
    readpos = write(fd, UInt64(0)) + sizepos

    pairs = discarded = truncated = 0

    @info "Building paired read datastore from FASTQ files"
    @info "Writing paired reads to datastore"

    while !eof(rdrx) && !eof(rdry)
        # Read in the two records.
        read!(rdrx, lread)
        read!(rdry, rread)

        llen = UInt64(FASTQ.seqlen(lread))
        rlen = UInt64(FASTQ.seqlen(rread))
        # If either read is too short, discard them both.
        if llen < minsize || rlen < minsize
            discarded += 1
            continue
        end

        if llen > maxsize
            truncated += 1
            ln = maxsize
        else
            ln = llen
        end

        if rlen > maxsize
            truncated += 1
            rn = maxsize
        else
            rn = rlen
        end

        pairs += 1

        # Copy FASTQ records to sequence variables, thus encoding
        # them in 2 bit format.
        copyto!(lseq, 1, lread, 1, ln)
        copyto!(rseq, 1, rread, 1, rn)

        # Write the two sizes reads one after the other:
        # For each sequence, write the size, and then the datachunk.

        write(fd, ln)
        write(fd, lseq.data)

        write(fd, rn)
        write(fd, rseq.data)
    end

    nreads::UInt64 = pairs * 2

    seek(fd, sizepos)
    write(fd, nreads)

    close(fd)

    @info "Done writing paired read sequences to datastore"
    @info string(discarded, " read pairs were discarded due to a too short sequence")
    @info string(truncated, " reads were truncated to ", maxsize, " base pairs")
    @info string("Created paired sequence datastore with ", pairs, " sequence pairs")

    stream = open(outfile, "r+")
    return PairedReads(outfile, name, name, maxsize, chunksize, fragsize,
                       readpos, nreads, orientation, stream)
end

function Base.open(::Type{PairedReads}, filename::String)
    fd = open(filename, "r")
    magic = read(fd, UInt16)
    dstype = reinterpret(Filetype, read(fd, UInt16))
    version = read(fd, UInt16)

    @assert magic == ReadDatastoreMAGIC
    @assert dstype == PairedDS
    @assert version == PairedDS_Version

    default_name = readuntil(fd, '\0')

    readsize = read(fd, UInt64)
    chunksize = read(fd, UInt64)
    fragsize = read(fd, UInt64)
    orientation = reinterpret(PairedReadOrientation, read(fd, UInt64))
    nreads = read(fd, UInt64)
    readpos_offset = position(fd)

    return PairedReads(filename, default_name, default_name,
                       readsize, chunksize, fragsize,
                       readpos_offset, nreads, orientation, fd)
end

@inline Base.length(prds::PairedReads) = prds.size

Base.summary(io::IO, prds::PairedReads) = print(io, "Paired Read Datastore '", prds.name, "': ", length(prds), " reads")

function Base.show(io::IO, prds::PairedReads)
    summary(io, prds)
end

bytes_per_read(prds::PairedReads) = (prds.chunksize + 1) * sizeof(UInt64)
@inline unsafe_read_offset_in_file(prds::PairedReads, idx::Integer) = prds.readpos_offset + (bytes_per_read(prds) * (idx - 1))

@inline _inbounds_index_of_sequence(prds::PairedReads, idx::Integer) = prds.readpos_offset + (bytes_per_read(prds) * (idx - 1))

@inline function Base.getindex(prds::PairedReads, idx::Integer)
    @boundscheck checkbounds(prds, idx)
    seq = eltype(prds)(prds.readsize)
    return inbounds_load_sequence!(prds, idx, seq)
end

ext: jl


text:
# bool
@test ismissing(rcopy(R"NA"))
v = [true ,true ,missing, false]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]

# integer
@test ismissing(rcopy(R"as.integer(NA)"))
v = [1 ,2 ,missing, 4]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]

v = Union{Int32, Missing}[1 ,2 ,missing, 4]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]

# float
@test ismissing(rcopy(R"as.double(NA)"))
v = [1.0 ,2.0 ,missing, 4.0]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]

v = Union{Float64, Missing}[1.0 ,2.0 ,missing, 4.0]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]

#string
@test ismissing(rcopy(R"as.character(NA)"))
v = ["a" ,"b" ,missing, "d"]
r = RObject(v)
@test isna(r, 3)
@test !isna(r, 1)
@test ismissing(rcopy(r)[3])
@test rcopy(r)[1:2] == v[1:2]


text:
using LoadQML

load_qml_app(joinpath(dirname(@__FILE__), "colors.qml"))


text:
using CompositeGrids, Test, StaticArrays, LinearAlgebra, Printf, Random, Statistics, JLD2, FileIO
# import Test: @test, @testset

if isempty(ARGS)
    include("chebyshev.jl")
    include("grid.jl")
    include("interpolate.jl")
    include("io.jl")
else
    include(ARGS[1])
end


text:
import Dolo

model = Dolo.Model(joinpath(Dolo.pkg_path, "examples","models","rbc_mc.yaml"))

Dolo.time_iteration(model)
dr0 = Dolo.time_iteration(model, tol_η=1e-6; grid=Dict(:n=>[5])).dr
sol = Dolo.time_iteration(model, maxit=1000; grid=Dict(:n=>[100]))

# faster version: limit the number of steps in the newton solver
@time sol = Dolo.time_iteration(model, maxit=1000; solver=Dict(:maxit=>2), grid=Dict(:n=>[100]))

# or use explicit formulas given in the model
@time sol = Dolo.time_iteration(model, maxit=1000; solver=Dict(:type=>:direct), grid=Dict(:n=>[100]))


text:
# Functions for performing shrinkage

function getprobabilitiesshrinkage(frequencies::Array{Int}, lambda::Void)
	function getlambda(n::Int, normalizedvalues, target)
		# Unbiased estimator of variance of u
		varu = normalizedvalues .* (1 - normalizedvalues) / (n - 1)
		msp = sum((normalizedvalues - target).^2) # misspecification ???

		# Estimate shrinkage intensity
		lambda = msp == 0 ? 1 : sum(varu) / msp

		# Make lambda be between 0 and 1 inclusive
		return lambda > 1 ? 1 : (lambda < 0 ? 0 : lambda)
	end

	target = 1 / length(frequencies) # Target is uniform distribution
	n = sum(frequencies)
	normalizedvalues = frequencies / n
	lambda = n == 1 || n == 0 ? 1 : getlambda(n, normalizedvalues, frequencies)

	return lambda * target + (1 - lambda) * normalizedvalues
end
function getprobabilitiesshrinkage(frequencies::Array{Int}, lambda::Real)
	target = 1 / length(frequencies) # Target is uniform distribution
	n = sum(frequencies)
	normalizedcounts = frequencies / n

	return lambda * target + (1 - lambda) * normalizedcounts
end


text:
function create_prior(observations, K::Int, dt::Float64)
    a = rand(Float64, K)
    a /= sum(a)
    A = rand(Float64, K, K)
    @inbounds for i = 1:K
        A[i, :] /= sum(A[i, :])
    end
    R = kmeans(abs2.((filter(!isnothing, observations))'), 3, tol = 1e-6; maxiter = 10000)
    D = R.centers # get the cluster centers
    D ./= 4dt
    # D .-= (abs2(er) / dt)
    D = reverse(sort(Array{Float64,1}(D[:])))
    a, A, D
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using PyPlot
if false
w = [3, 4]
fig, axs = PyPlot.subplots(1, 1)
annotate(
L"$\theta: $" * "$w" * ".\n" *
L"$\lambda: $" * "$range" * ".\n"*
"quad: " * "$quadType",
 xy=[.1;.2],
    fontsize=80.0
)
end
rangeθ = [2, 3]
plt, axs = PyPlot.subplots(2, 2)
axs[2, 2].annotate(L"$\theta: $" * "$rangeθ" * ".\n",  xy=[.1;.2],
fontsize=80.0)

    #xycoords="axes fraction",
	# xytext=[-10,10],
	# textcoords="offset points",
    #
    #
    #ha="center",
    #va="bottom"

    #	xy=[.4;.5],

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# London Tube Lines

#import Pkg; Pkg.add("VegaLite")
#Pkg.add("VegaDatasets")
using VegaLite, VegaDatasets

londonBoroughs = dataset("londonBoroughs")
londonCentroids = dataset("londonCentroids")
londonTubeLines = dataset("londonTubeLines")

@vlplot(
    width=700, height=500,
    config={
        view={
            stroke=:transparent
        }
    }
) +
@vlplot(
    data={
        values=londonBoroughs,
        format={
            type=:topojson,
            feature=:boroughs
        }
    },
    mark={
        :geoshape,
        stroke=:white,
        strokeWidth=2
    },
    color={value="#d1ffc2"}
) +
@vlplot(
    data=londonCentroids,
    transform=[{
        calculate="indexof (datum.name,' ') > 0  ? substring(datum.name,0,indexof(datum.name, ' ')) : datum.name",
        as=:bLabel
    }],
    mark=:text,
    longitude="cx:q",
    latitude="cy:q",
    text="bLabel:n",
    size={value=8},
    opacity={value=0.6}
) +
@vlplot(
    data={
        values=londonTubeLines,
        format={
            type=:topojson,
            feature=:line
        }
    },
    mark={
        :geoshape,
        filled=false,
        strokeWidth=2
    },
    color={
        "id:n",
        legend={
            title=nothing,
            orient="bottom-right",
            offset=0
        },
        scale={
            domain=[
                "Bakerloo",
                "Central",
                "Circle",
                "District",
                "DLR",
                "Hammersmith & City",
                "Jubilee",
                "Metropolitan",
                "Northern",
                "Piccadilly",
                "Victoria",
                "Waterloo & City"
            ],
            range=[
                "rgb(137,78,36)",
                "rgb(220,36,30)",
                "rgb(255,206,0)",
                "rgb(1,114,41)",
                "rgb(0,175,173)",
                "rgb(215,153,175)",
                "rgb(106,114,120)",
                "rgb(114,17,84)",
                "rgb(0,0,0)",
                "rgb(0,24,168)",
                "rgb(0,160,226)",
                "rgb(106,187,170)"
            ]
        }
    }
)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
type MaxoutLayer <: Layer
    inputActivation::Vol
    outputActivation::Vol
    groupSize::Int
    argmax::Array{Int64, 3}

    function MaxoutLayer(groupSize::Int, dimensions)
        (inputX, inputY, inputDepth) = dimensions
        outputDepth = floor(inputDepth / groupSize)
        outputDimensions = (inputX, inputY, outputDepth)
        new(zeros(dimensions), zeros(outputDimensions), groupSize, zeros(Int64, outputDimensions))
    end
end

function forward(self::MaxoutLayer, inputActivation::Vol)
    self.inputActivation = inputActivation
    (outputX, outputY, outputDepth) = size(self.outputActivation)
    for d in 1:outputDepth, y in 1:outputY, x in 1:outputX
        maxValue, maxIndex = findmax([V[x, y, ds] for ds in (d - 1) * groupSize + 1:d * groupSize])
        self.outputActivation.dw[x, y, d] = maxValue
        self.argmax[x, y, d] = maxIndex
    end
    return self.outputActivation
end

function backward(self::MaxoutLayer)
    fill!(self.inputActivation.dw, 0.0)
    (outputX, outputY, outputDepth) = size(self.outputActivation)
    for d in 1:outputDepth, y in 1:outputY, x in 1:outputX
        chainGrad = self.outputActivation.dw[x, y, d]
        maxIndex = self.argmax[x, y, d]
        maxInputIndex = (d-1) * groupSize + maxIndex
        self.inputActivation.dw[x, y, maxinputIndex] = chainGrad
    end
end

gradientStates(self::MaxoutLayer) = []


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libsuitesparseconfig"], :suitesparseconfig),
    LibraryProduct(prefix, String["libamd"], :amd),
    LibraryProduct(prefix, String["libbtf"], :btf),
    LibraryProduct(prefix, String["libcamd"], :camd),
    LibraryProduct(prefix, String["libccolamd"], :ccolamd),
    LibraryProduct(prefix, String["libcolamd"], :colamd),
    LibraryProduct(prefix, String["libcholmod"], :cholmod),
    LibraryProduct(prefix, String["libldl"], :ldl),
    LibraryProduct(prefix, String["libklu"], :klu),
    LibraryProduct(prefix, String["libumfpack"], :umfpack),
    LibraryProduct(prefix, String["librbio"], :rbio),
    LibraryProduct(prefix, String["libspqr"], :spqr),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaLinearAlgebra/SuiteSparseBuilder/releases/download/v5.2.0-0.2.20"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/SuiteSparse.aarch64-linux-gnu.tar.gz", "f4b08b4b1dee6e6e7929ea2059309eefdc395554ee7d771a62e98779929aad1d"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/SuiteSparse.arm-linux-gnueabihf.tar.gz", "d80c1718d8ce2835877907a2366a2b412f7920c8b79d0eaa9b4f5f4f3c464e2f"),
    Linux(:i686, :glibc) => ("$bin_prefix/SuiteSparse.i686-linux-gnu.tar.gz", "0eaa7a6eff01ce3c9b4a2d08350fba875f5559d349e018d3a839fcfe10964533"),
    Windows(:i686) => ("$bin_prefix/SuiteSparse.i686-w64-mingw32.tar.gz", "7284c271c9ba7731bb88ab12d912173479c568093f09ffc1e8b352b68b54b8c4"),
    MacOS(:x86_64) => ("$bin_prefix/SuiteSparse.x86_64-apple-darwin14.tar.gz", "6e7e4a5590ac7458026b2dafb49db6c5adaae3a30e98e8b627b8b97185d3abe1"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/SuiteSparse.x86_64-linux-gnu.tar.gz", "e1a0ce686e696be79d62ad10160137911510c7618612a5f3019d002ef6903945"),
    Windows(:x86_64) => ("$bin_prefix/SuiteSparse.x86_64-w64-mingw32.tar.gz", "4e6e104ba32dfff14eb7a68bc7557c9fc9d4434c935c7ca6cfbc7c243dd5d25a"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
import Base: -, +, *, /, ∘, zero, conj, real, imag, abs

### On scalar grid data ####

# Set it to negative of itself
function (-)(p_in::ScalarGridData)
  p = deepcopy(p_in)
  p.data .= -p.data
  return p
end

function (-)(p1::T,p2::T) where {T <: ScalarGridData}
   return T(p1.data .- p2.data)
 end

function (+)(p1::T,p2::T) where {T <: ScalarGridData}
  return T(p1.data .+ p2.data)
end

# Multiply and divide by a constant
function (*)(p::T,c::Number) where {T<:ScalarGridData}
  return T(c*p.data)
end

function (/)(p::T,c::Number) where {T<:ScalarGridData}
  return T(p.data ./ c)
end

function product!(out::Nodes{T, NX, NY},
                  p::Nodes{T, NX, NY},
                  q::Nodes{T, NX, NY}) where {T, NX, NY}

    inds = node_inds(T, (NX, NY))
    @inbounds for y in 1:inds[2], x in 1:inds[1]
        out[x,y] = p[x,y] * q[x,y]
    end
    out
end

function product(p::Nodes{T, NX, NY}, q::Nodes{T, NX, NY}) where {T, NX, NY}
    product!(Nodes(T, p), p, q)
end

function (∘)(p::Nodes{T, NX, NY}, q::Nodes) where {T, NX, NY}
    product!(Nodes(T, p), p, q)
end

(*)(c::Number,p::T) where {T<:GridData} = *(p,c)

### On vector grid data ####

function (-)(p_in::VectorGridData)
  p = deepcopy(p_in)
  p.u .= -p.u
  p.v .= -p.v
  return p
end

function (-)(p1::T,p2::T) where {T<:VectorGridData}
  return T(p1.u - p2.u, p1.v - p2.v)
end

function (+)(p1::T,p2::T) where {T<:VectorGridData}
  return T(p1.u + p2.u, p1.v + p2.v)
end

function (*)(p::T,c::Number) where {T<:VectorGridData}
  return T(c*p.u,c*p.v)
end

function (/)(p::T,c::Number) where {T<:VectorGridData}
  return T(p.u / c, p.v / c)
end

"""
    product!(out::Edges/Nodes,p::Edges/Nodes,q::Edges/Nodes)

Compute the Hadamard (i.e. element by element) product of edge or nodal
(primal or dual) data `p` and `q` and return the result in `out`.

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> out = p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> product!(out,p,q)
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function product!(out::Edges{T, NX, NY},
                  p::Edges{T, NX, NY},
                  q::Edges{T, NX, NY}) where {T, NX, NY}

    uinds, vinds = edge_inds(T, (NX, NY))
    @inbounds for y in 1:uinds[2], x in 1:uinds[1]
        out.u[x,y] = p.u[x,y] * q.u[x,y]
    end

    @inbounds for y in 1:vinds[2], x in 1:vinds[1]
        out.v[x,y] = p.v[x,y] * q.v[x,y]
    end
    out
end

"""
    product(p::Edges/Nodes,q::Edges/Nodes) --> Edges/Nodes

Compute the Hadamard product of edge or nodal (primal or dual) data `p` and `q` and return
the result. This operation can also be carried out with the `∘` operator:

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> p∘q
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function product(p::Edges{T, NX, NY}, q::Edges{T, NX, NY}) where {T, NX, NY}
    product!(Edges(T, p), p, q)
end

function (∘)(p::Edges{T, NX, NY}, q::Edges) where {T, NX, NY}
    product!(Edges(T, p), p, q)
end

### ON TENSORS ####

function (-)(p_in::EdgeGradient)
  p = deepcopy(p_in)
  p.dudx .= -p.dudx
  p.dvdy .= -p.dvdy
  p.dudy .= -p.dudy
  p.dvdx .= -p.dvdx
  return p
end

function (-)(p1::T,p2::T) where {T <: EdgeGradient}
  return T(p1.dudx - p2.dudx, p1.dvdy - p2.dvdy, p1.dudy - p2.dudy, p1.dvdx - p2.dvdx)
end

function (+)(p1::T,p2::T) where {T <: EdgeGradient}
  return T(p1.dudx + p2.dudx, p1.dvdy + p2.dvdy, p1.dudy + p2.dudy, p1.dvdx + p2.dvdx)
end

function (*)(p::T,c::Number) where {T <: EdgeGradient}
  return T(c*p.dudx,c*p.dvdy,c*p.dudy,c*p.dvdx)
end

function (/)(p::T,c::Number) where {T <: EdgeGradient}
  return T(p.dudx / c, p.dvdy / c, p.dudy / c, p.dvdx / c)
end

#### ON ALL TYPES ####

zero(::Type{T}) where {T <: GridData} = T()

#### ON COMPLEX GRID DATA

for f in (:conj,)
    @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
        Acopy = deepcopy(A)
        Acopy .= broadcast($f,Acopy)
        return Acopy
    end
end

for f in (:real, :imag, :abs)
  @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
      Acopy = similar(A,element_type=Float64)
      Acopy .= broadcast($f,A)
      return Acopy
  end
end


text:
n = 1000

xs = convert(Array{Any}, readlines())

# Parse input
for (i, x) in enumerate(xs)
    dash = findfirst('-', x)
    colon = findnext(':', x, dash)

    mn = parse(Int, x[begin:dash - 1])
    mx = parse(Int, x[dash + 1:colon - 2])
    c = x[colon - 1]
    pass = x[colon + 2:end]

    xs[i] = (mn, mx, c, pass)
end

# Part 1
cnt = 0
for (mn, mx, c, pass) in xs
    occ = count("$c", pass)

    if mn <= occ <= mx
        global cnt

        cnt += 1
    end
end

println("Part 1 : $cnt")

# Part 2
cnt = 0
for (mn, mx, c, pass) in xs
    # ⊻ is the bitwise xor
    if (pass[mn] == c) ⊻ (pass[mx] == c)
        global cnt

        cnt += 1
    end
end

println("Part 2 : $cnt")


text:
#= Visualize the hierarchical structure given the position of nodes
returned by function getpos =#

    function hierviz(N, plateList, Pos)

    Val = sum(convert(Array, value.(Pos)), dims = 1)
    #Val = sum(Pos, dims = 1)

    nNodes = 0
    for i in 1:length(N)
        nNodes += length(N[i])
    end

    opt = round.((Val .+ 1 .+ nNodes) ./2)

    nodeList = collect(1:nNodes)
    newNodeList = collect(1:nNodes)
    nPlate = length(plateList)

    for i in 1:nNodes
        newNodeList[Int(opt[i])] = nodeList[i]
    end

    Drawing(1800, 1000)
    origin()
    background("white")
    sethue("black")
    fontsize(20)

        radius = 30
        xloca = (opt.-1).*3 .*radius

    for i in 1:nNodes
        u = Point(xloca[i], 0)
        text("u$i", u, halign=:center, valign=:middle)
        circle.(u, radius, :stroke)
    end

    for i in 1:nPlate
        for j in 1:(nNodes - length(plateList[i]) + 1)
            if sort(newNodeList[j:(j + length(plateList[i]) - 1)]) == sort(plateList[i])
                randomhue()
                box(Point(xloca[newNodeList[j]] - radius - 20*rand(), - radius - 20*rand()),
                    Point(xloca[newNodeList[j + length(plateList[i]) - 1]] + radius + 20*rand(), radius + 20* rand()), :stroke)
            end
        end
    end

    finish()
    preview()

end


text:
using QuadGK
using Statistics
using GNSSTools
using ProgressMeter
using JLD
using Random
using Optim
using PyPlot
pygui(true)
include("kalman_filter_gain_error_functions.jl")

# k₁ = 0.017299
# k₂ = 0.024087
# B2_truth = (k₁ + k₂/k₁)/4
# H2 = init_H(k₁, k₂)
# B2 = get_B(H2)

# H3 = init_H(0.5, 2.0, 4)
# B3 = get_B(H3)

# h₀ = 1e-19
# h₋₂ = 2e-20
# T = 0.001  # s
# qₐ = 0.1
# CN0 = 30 # dB⋅Hz
# sig_freq = L1_freq
# state_num = 2
# # R = [1e-4]
# R = [1/(2*10^(CN0/10)*T)*(1 + 1/(2*10^(CN0/10)*T))]
# A = calcA(T, state_num)
# C = calcC(T, state_num)
# Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num)
# K = dkalman(A, C, Q, GNSSTools.Diagonal(R))
# if length(K) == 2
#     k1 = K[1]
#     k2 = K[2]
#     B_truth = (k1 + k2/k1)/4
#     B = get_B(init_H(k1, k2))
#     println([B, B_truth, B_truth-B])
# else
#     k1 = K[1]
#     k2 = K[2]
#     k3 = K[3]
#     B = get_B(init_H(k1, k2, k3))
#     println(B)
# end

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

# t_length = 1
# f_s = 5e6
# # f_s = 2.048e6
# # f_s = 25e6
# # h_parms = h_parms_tcxo[1]
# h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # TCXO
# # h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO

# sigtype = define_l1ca_code_type(; sig_freq=1*L1_freq)
# # sigtype = define_l5_code_type(t_length; sig_freq=1*L5_freq, channel="both")
# # sigtype = definesignaltype(definecodetype(l5i_codes, L5_chipping_rate),
#                         #    L5_freq, "I")
# # signal = definesignal(sigtype, 5e6, t_length;
# #                       prn=1, code_start_idx=1 000, f_d=163000, fd_rate=449,
# #                       CN0=45, include_phase_noise=false,
# #                       include_thermal_noise=true, receiver_h_parms=h_parms);
# signal = definesignal(sigtype, f_s, t_length;
#                       prn=1, code_start_idx=1000, f_d=1000, fd_rate=0, f_if=1.25e6,
#                       CN0=25, include_phase_noise=true,
#                       include_thermal_noise=true, include_adc=true,
#                       include_carrier=true, receiver_h_parms=h_parms)
# generatesignal!(signal);
# # signal.noexp = false
# # generatereplica!(signal);

# acqresults, trackresults, corr_result, SNR_est = process(signal,
#                                                          sigtype,
#                                                          1, "I"; fd_range=5000,
#                                                          fine_acq_method=:fft,
#                                                          return_corrresult=true,
#                                                          h₀=h_parms[3],
#                                                          h₋₂=h_parms[1],
#                                                          q_mult=1,
#                                                          cov_mult=1,
#                                                          R_mult=1,
#                                                          q_a=0.1,
#                                                          dynamickf=true,
#                                                          state_num=3,
#                                                          dll_b=0.1,
#                                                          acquisition_T=1e-3,
#                                                          fine_acq_T=20e-3,
#                                                          tracking_T=4e-3,
#                                                          M=20,
#                                                          σᵩ²=missing)

# figure()
# if size(trackresults.K)[1] == 3
#     Bs = get_B.(init_H.(trackresults.K[1,:], trackresults.K[2,:],
#                         trackresults.K[3,:]))
# else
#     Bs = get_B.(init_H.(trackresults.K[1,:], trackresults.K[2,:]))
# end
# plot(trackresults.t, Bs)
# ylim([0, 10])
# xlabel("Time (s)")
# ylabel("Equivalent Bandwidth (Hz)")

# figure()
# plot(trackresults.t, trackresults.K[1,:], label="K₁")
# plot(trackresults.t, trackresults.K[2,:], label="K₂")
# if size(trackresults.K)[1] == 3
#     plot(trackresults.t, trackresults.K[3,:], label="K₃")
# end
# ylim([0,10])
# xlabel("Time (s)")
# ylabel("Kₙ (unitless)")
# legend()

# h₀ = 1e-19
# h₋₂ = 2e-20
# T = 0.001  # s
# qₐ = 0.1
# CN0 = 36 # dB⋅Hz
# sig_freq = L1_freq
# state_num = 2
# # R = [1e-4]
# R = [1/(2*10^(CN0/10)*T)]
# A = calcA(T, state_num)
# C = calcC(T, state_num)
# Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num)
# K = dkalman(A, C, Q, GNSSTools.Diagonal(R))
# if length(K) == 2
#     k1 = K[1]
#     k2 = K[2]
#     B_truth = (k1 + k2/k1)/4
#     B = get_B(init_H(k1, k2))
#     println([B, B_truth, B_truth-B])
# else
#     k1 = K[1]
#     k2 = K[2]
#     k3 = K[3]
#     B = get_B(init_H(k1, k2, k3))
#     println(B)
# end

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

# """
#     sim_KF_B(state_num)

# Estimates the steady state KF gain for various C/N₀ and qₐ.
# """
# function sim_KF_B(CN0s, qas, state_num; iterations=25, T=1e-3)
#     t_length = 1
#     f_s = 5e6
#     h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # TCXO
#     # h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO
#     sigtype = define_l1ca_code_type(; sig_freq=1*L1_freq)
#     signal = definesignal(sigtype, f_s, t_length;
#                         prn=1, code_start_idx=1, f_d=1000, fd_rate=0,
#                         f_if=1.25e6, CN0=36, include_phase_noise=true,
#                         include_thermal_noise=true, include_adc=true,
#                         include_carrier=true, receiver_h_parms=h_parms)
#     Bs = zeros(length(qas), length(CN0s))
#     p = Progress(iterations*length(CN0s)*length(qas), 1, "Processing...")
#     for j in 1:length(CN0s)
#         CN0 = CN0s[j]
#         for i in 1:length(qas)
#             for iteration in 1:iterations
#                 ϕ₀ = rand(0:0.0001:2π)
#                 f_d = rand(-5e3:0.1:5e5)
#                 # code_start_idx = rand(1:1:floor(Int, T*f_s))
#                 code_start_idx = 1
#                 definesignal!(signal; CN0=CN0, new_phase_noise=true,
#                              new_thermal_noise=true, phi=ϕ₀, f_d=f_d,
#                              code_start_idx=code_start_idx)
#                 generatesignal!(signal)
#                 # σᵩ² = 1/(2*10^(CN0s[j]/10)*T)  # rad²
#                 σᵩ² = GNSSTools.phase_noise_variance(CN0, T)
#                 (acqresults, trackresults, corr_result,
#                  SNR_est) = process(signal,
#                                     sigtype,
#                                     1, "I"; fd_range=5000,
#                                     fine_acq_method=:fft,
#                                     return_corrresult=true,
#                                     h₀=h_parms[3],
#                                     h₋₂=h_parms[1],
#                                     q_mult=1,
#                                     σω=1000,
#                                     cov_mult=1,
#                                     R_mult=1,
#                                     q_a=qas[i],
#                                     dynamickf=true,
#                                     state_num=state_num,
#                                     dll_b=0.1,
#                                     acquisition_T=1e-3,
#                                     fine_acq_T=20e-3,
#                                     tracking_T=T,
#                                     M=20,
#                                     show_plot=false,
#                                     σᵩ²=σᵩ²)
#                 if size(trackresults.K)[1] == 3
#                     Bs[i,j] += get_B(init_H(trackresults.K[1,end],
#                                             trackresults.K[2,end],
#                                             trackresults.K[3,end]))
#                 else
#                     Bs[i,j] += get_B(init_H(trackresults.K[1,end],
#                                             trackresults.K[2,end]))
#                 end
#                 # Bs[i,j] += get_B(init_H(trackresults.K[:,end]...))
#                 next!(p)
#             end
#             Bs[i,j] = Bs[i,j] / iterations
#         end
#     end
#     return (Bs, CN0s, qas)
# end

# CN0s = Array(range(25, 50, step=1))
# # qas = [0.0, 0.1, 1, 10, 100]
# qas = [0, 1.0]
# Bs_2state, CN0s, qas_ = sim_KF_B(CN0s, 0.0, 2; iterations=10, T=1e-3)
# Bs_3state, CN0s, qas = sim_KF_B(CN0s, qas, 3; iterations=10, T=1e-3)

# save("ch3_kf_pll_bandwidth_3-state.jld", "Bs_2state", Bs_2state)
# save("ch3_kf_pll_bandwidth_3-state.jld", "Bs_3state", Bs_3state)
# # results = load("ch4_acquisition_performance.jld", "results")

# # figure(figsize=(7.5,2.5))
# # ax1 = subplot(1,2,1)
# # for i in 1:length(qas)
# #     ax1.plot(CN0s, Bs_2state[i,:], label="qₐ = $(qas[i])")
# # end
# # xlabel("C/N₀ (dB⋅Hz)")
# # ylabel("Bandwidth (Hz)")
# # title("(a)")
# # ax1.legend()
# figure(figsize=(7.5,3))
# ax2 = subplot(1,1,1)
# for i in reverse(1:length(qas))
#     ax2.plot(CN0s, Bs_3state[i,:], label="3-state, qₐ = $(qas[i])")
# end
# ax2.plot(CN0s, Bs_2state[1,:], "k-", label="2-state")
# xlabel("C/N₀ (dB⋅Hz)")
# ylabel("Bandwidth (Hz)")
# ax2.legend()
# subplots_adjust(wspace=0.4, bottom=0.15, left=0.08, right=0.92, top=0.95)
# savefig(string(directory, "figures/ch3_kf_pll_bandwidth.svg"), dpi=300)

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

function estimate_BKσ(T, CN0, state_num, qₐ, h_parms, sig_freq)
    h₀ = h_parms[3]
    h₋₁ = h_parms[2]
    h₋₂ = h_parms[1]
    A = calcA(T, state_num)
    C = calcC(T, state_num)
    Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num)
    σᵩ²_unfiltered_expected = phase_noise_variance(CN0, T)
    R = [σᵩ²_unfiltered_expected]
    K = dkalman(A, C, Q, GNSSTools.Diagonal(R))
    # B = get_B(init_H(K...))
    B = calc_B(K...)
    σᵩ²_filtered_expected = phase_noise_variance(B, CN0, T, h₀, h₋₁, h₋₂)
    return (K, B, σᵩ²_unfiltered_expected, σᵩ²_filtered_expected)
end

function get_kf_steady_state_info(h_parms, Ts, CN0s, qₐ, state_num)
    Ks = zeros(length(h_parms), length(Ts), length(CN0s), state_num)
    Bs = zeros(length(h_parms), length(Ts), length(CN0s))
    σᵩ²s_unfiltered = zeros(length(h_parms), length(Ts), length(CN0s))
    σᵩ²s_filtered = zeros(length(h_parms), length(Ts), length(CN0s))
    p = Progress(length(h_parms)*length(CN0s)*length(Ts), 1, "Processing...")
    for i in 1:length(h_parms)
        h_parm = h_parms[i]
        for j in 1:length(Ts)
            T = Ts[j]
            for k in 1:length(CN0s)
                CN0 = CN0s[k]
                K, B, σᵩ²_unfiltered, σᵩ²_filtered = estimate_BKσ(T, CN0,
                                                                 state_num, qₐ,
                                                                 h_parm,
                                                                 sig_freq)
                Ks[i,j,k,:] = K
                Bs[i,j,k] = B
                σᵩ²s_unfiltered[i,j,k] = σᵩ²_unfiltered
                σᵩ²s_filtered[i,j,k] = σᵩ²_filtered
                next!(p)
            end
        end
    end
    return (Ks, Bs, σᵩ²s_unfiltered, σᵩ²s_filtered)
end

sig_freq = L1_freq
Ts = [1e-3, 2e-3, 5e-3, 10e-3, 20e-3]  # seconds
CN0s = Array(range(20, 50, step=1))

tcxo_h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # TCXO
ocxo_h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO
h_parms = [tcxo_h_parms, ocxo_h_parms]

fig = figure(figsize=(7.5, 3.5))
cmap = get_cmap("viridis")

# 2-state KF
qₐ = 0
state_num = 2
Ks, Bs, σᵩ²s_unfiltered, σᵩ²s_filtered = get_kf_steady_state_info(h_parms, Ts,
                                                                 CN0s, qₐ,
                                                                 state_num);
# ax1 = fig.add_subplot(2,2,1)
ax1 = fig.add_subplot(1,1,1)
for i in 1:length(Ts)
    color = cmap(float(i)/length(Ts))
    ax1.plot(CN0s, Bs[1,i,:], "-", color=color,
             label="T = $(floor(Int, Ts[i]*1000))ms; TCXO")
    # ax1.plot(CN0s, Bs[2,i,:], ":", color=color,
            #  label="T = $(floor(Int, Ts[i]*1000))ms; OCXO")
end
for i in 1:length(Ts)
    color = cmap(float(i)/length(Ts))
    ax1.plot(CN0s, Bs[2,i,:], ":", color=color,
             label="T = $(floor(Int, Ts[i]*1000))ms; OCXO")
end

xlabel("C/N₀ (dB⋅Hz)")
ylabel("B (Hz)")
# title("(a)")
legend(ncol=2)

# # 3-state KF
# qₐ = 0.1
# state_num = 3
# Ks, Bs, σᵩ²s_unfiltered, σᵩ²s_filtered = get_kf_steady_state_info(h_parms, Ts,
#                                                                  CN0s, qₐ,
#                                                                  state_num);
# ax2 = fig.add_subplot(2,2,2)
# for i in 1:length(Ts)
#     color = cmap(float(i)/length(Ts))
#     ax2.plot(CN0s, Bs[1,i,:], "-", color=color)
#     ax2.plot(CN0s, Bs[2,i,:], ":", color=color,
#              label="T = $(floor(Int, Ts[i]*1000))ms; OCXO")
# end
# xlabel("C/N₀ (dB⋅Hz)")
# ylabel("B (Hz)")
# title("(b)")
# legend()

# # 3-state KF
# qₐ = 1
# state_num = 3
# Ks, Bs, σᵩ²s_unfiltered, σᵩ²s_filtered = get_kf_steady_state_info(h_parms, Ts,
#                                                                  CN0s, qₐ,
#                                                                  state_num);
# ax3 = fig.add_subplot(2,2,3)
# for i in 1:length(Ts)
#     color = cmap(float(i)/length(Ts))
#     ax3.plot(CN0s, Bs[1,i,:], "-", color=color)
#     ax3.plot(CN0s, Bs[2,i,:], ":", color=color)
# end
# xlabel("C/N₀ (dB⋅Hz)")
# ylabel("B (Hz)")
# title("(c)")

# # 3-state KF
# qₐ = 100
# state_num = 3
# Ks, Bs, σᵩ²s_unfiltered, σᵩ²s_filtered = get_kf_steady_state_info(h_parms, Ts,
#                                                                  CN0s, qₐ,
#                                                                  state_num);
# ax4 = fig.add_subplot(2,2,4)
# for i in 1:length(Ts)
#     color = cmap(float(i)/length(Ts))
#     ax4.plot(CN0s, Bs[1,i,:], "-", color=color)
#     ax4.plot(CN0s, Bs[2,i,:], ":", color=color)
# end
# xlabel("C/N₀ (dB⋅Hz)")
# ylabel("B (Hz)")
# title("(d)")

subplots_adjust(wspace=0.25, hspace=0.45, bottom=0.15,
                left=0.08, right=0.92, top=0.94)

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
savefig(string(directory, "figures/ch3_kf_pll_bandwidth_numerical_2-state.svg"), dpi=300)


text:
# Helper structs to dispatch known decision variable construction
struct NoSpecifiedConstraint <: MOI.AbstractSet end
Base.copy(set::NoSpecifiedConstraint) = set

struct DecisionSet{S <: MOI.AbstractSet} <: MOI.AbstractScalarSet
    constraint::S

    function DecisionSet(; constraint::MOI.AbstractSet = NoSpecifiedConstraint())
        S = typeof(constraint)
        return new{S}(constraint)
    end
end
struct KnownSet <: MOI.AbstractScalarSet end

struct SingleDecisionSet{T, S} <: MOI.AbstractScalarSet
    decision::Decision{T}
    constraint::S
end
Base.copy(set::SingleDecisionSet) = reuse(set, set.decision)

struct SingleKnownSet{T} <: MOI.AbstractScalarSet
    known::Decision{T}
end
Base.copy(set::SingleKnownSet) = reuse(set, set.known)

set_constraint(set::SingleDecisionSet) = set.constraint
set_constraint(::SingleKnownSet) = NoSpecifiedConstraint()

struct FreeDecision <: MOI.AbstractScalarSet end

struct MultipleDecisionSet{T, S} <: MOI.AbstractVectorSet
    decisions::Vector{Decision{T}}
    constraint::S
end
MOI.dimension(set::MultipleDecisionSet) = length(set.decisions)
Base.copy(set::MultipleDecisionSet) = reuse(set, set.decisions)

struct MultipleKnownSet{T} <: MOI.AbstractVectorSet
    knowns::Vector{Decision{T}}
end
MOI.dimension(set::MultipleKnownSet) = length(set.knowns)
Base.copy(set::MultipleKnownSet) = reuse(set, set.knowns)

MOIU.variable_function_type(::Type{<:SingleDecisionSet}) = SingleDecision
MOIU.variable_function_type(::Type{<:MultipleDecisionSet}) = VectorOfDecisions
MOIU.variable_function_type(::Type{<:SingleKnownSet}) = SingleKnown
MOIU.variable_function_type(::Type{<:MultipleKnownSet}) = VectorOfKnowns

function JuMP.in_set_string(print_mode, set::SingleDecisionSet)
    if set.constraint == NoSpecifiedConstraint()
        return string(JuMP._math_symbol(print_mode, :in), " Decisions")
    else
        return string(JuMP._math_symbol(print_mode, :in), " Decisions($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function JuMP.in_set_string(print_mode, set::SingleKnownSet)
    return string(JuMP._math_symbol(print_mode, :in), " Known(value = $(set.known.value))")
end

function JuMP.in_set_string(print_mode, ::MultipleDecisionSet)
    if set.constraint == NoSpecifiedConstraint()
        return string(JuMP._math_symbol(print_mode, :in), " Decisions")
    else
        return string(JuMP._math_symbol(print_mode, :in), " Decisions($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function JuMP.in_set_string(print_mode, set::MultipleKnownSet)
    return string(JuMP._math_symbol(print_mode, :in), " Known(values = $([k.value for k in set.knowns]))")
end

function reuse(set::SingleDecisionSet, decision::Decision)
    return SingleDecisionSet(decision, copy(set.constraint))
end

function reuse(set::SingleKnownSet, decision::Decision)
    return SingleKnownSet(decision)
end

function reuse(set::MultipleDecisionSet, decisions::Vector{<:Decision})
    return MultipleDecisionSet(decisions, copy(set.constraint))
end

function reuse(set::MultipleKnownSet, decisions::Vector{<:Decision})
    return MultipleKnownSet(decisions)
end

function VariableRef(model::Model, index::MOI.VariableIndex, ::Union{SingleDecisionSet, MultipleDecisionSet})
    return DecisionRef(model, index)
end

function VariableRef(model::Model, index::MOI.VariableIndex, ::Union{SingleKnownSet, MultipleKnownSet})
    return KnownRef(model, index)
end

function set(variable::JuMP.ScalarVariable, ::Type{DecisionRef}, constraint::MOI.AbstractSet)
    return SingleDecisionSet(Decision(variable.info, Float64), constraint)
end

function set(variable::JuMP.ScalarVariable, ::Type{KnownRef}, ::NoSpecifiedConstraint)
    return SingleKnownSet(KnownDecision(variable.info, Float64))
end

function set_decision!(decisions::Decisions, index::MOI.VariableIndex, set::SingleDecisionSet)
    set_decision!(decisions, index, set.decision)
    return nothing
end

function set_decision!(decisions::Decisions, index::MOI.VariableIndex, set::SingleKnownSet)
    set_known_decision!(decisions, index, set.known)
    return nothing
end

function set_decision!(decisions::Decisions, index::MOI.VariableIndex, set_index::Int, set::MultipleDecisionSet)
    set_decision!(decisions, index, set.decisions[set_index])
    return nothing
end

function set_decision!(decisions::Decisions, index::MOI.VariableIndex, set_index::Int, set::MultipleKnownSet)
    set_known_decision!(decisions, index, set.knowns[set_index])
    return nothing
end


text:
using Test
using FileIO, Images, Colors
using CSV

# !!! commented out because of local repository of test images. !!!

# base_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/"
# dataset_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/dataset/"
# landmark_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/landmarks/"
# loc_table = CSV.read(base_path * "location_table.csv")
#
# train_rows = loc_table[loc_table[:status] .== "training", :]
#
# @testset "LAP_julia.jl" begin
#     TEST_SIZE = 2
#     @testset "Tests on random images $k" for k = 1:TEST_SIZE
#         random_rows = train_rows[rand(1:size(train_rows, 1), TEST_SIZE), :]
#         source_loc = (random_rows[k, :])[Symbol("Source image")]
#         target_loc = (random_rows[k, :])[Symbol("Target image")]
#         source_in = load(dataset_path * source_loc)
#         target_in = load(dataset_path * target_loc)
#         gray_source = Gray.(source_in)
#         gray_target = Gray.(target_in)
#
#         @testset "rescale!" begin
#             both((x,y)) = [x[:]; y[:]]
#             @test maximum(both(LAP_julia.helpers.rescale!(gray_source, gray_target))) == 1
#             @test minimum(both(LAP_julia.helpers.rescale!(gray_source, gray_target))) == 0
#         end
#     end
# end

# @testset "Registration Algorithms" begin
#     TimerOutputs.enable_debug_timings(LAP_julia.registration)
#     img, imgw, flow = gen_init(:lena, flow_args=[15, 140])
#
#     @testset "PFLAP" begin
#         timer = TimerOutput("Registration");
#         @timeit timer "polyfilter lap" begin
#             flow_est, source_reg = pflap(img, imgw, display=false, timer=timer)
#         end
#         print_timer(timer)
#
#         flow_est, source_reg = pflap(img, imgw, display=false)
#         showflow(flow)
#         showflow(flow_est)
#         @test angle_mae(flow, flow_est) < 10

@testset "helpers" begin
    @testset "pad_images $k, $l, $m, $n" for k = 1:3, l = 1:3, m = 1:3, n = 1:3
        im1 = ones(k, l)
        im2 = ones(m, n)
        @test size(LAP_julia.pad_images(im1, im2)[1]) <= (k+m-1, l+n-1)
        @test size(LAP_julia.pad_images(im1, im2)[1]) == ((k<m) ? m : k, (l<n) ? n : l)
    end
end

@testset "inpaint" begin
    u = ones(10,10) .+ 2im .* ones(10,10);
    nan_u = u
    nan_u[3:8, 3:8] .= NaN .+ NaN .* 1im;
    LAP_julia.inpaint.inpaint_nans!(nan_u)
    @test all(nan_u .== u)
    @test any(isnan.(nan_u)) == false
end


text:
mutable struct Weight <: FieldVector{2,Float64}
    d::Float64
    e::Float64
    Weight() = new(0.0, 0.0)
    Weight(d, e) = new(d, e)
end

const Base.zero(::Type{Weight}) = Weight(0.0, 0.0)
const Base.abs(w::Weight) = abs(w.d) + abs(w.e) # define abs(Weight)

function evalAllG!(G, K, T0idx, varT, evalG; kwargs...)
    for g in G
        tin, tout = g.Tpair
        g.weight = evalG(K, varT[T0idx+tin], varT[T0idx+tout], kwargs...)
    end
end

# function phase(varT, Tpair, isF)
#     tInL, tOutL, tInR, tOutR = varT[Tpair[INL]], varT[Tpair[OUTL]], varT[Tpair[INR]],
#     varT[Tpair[OUTR]]
#     if (isF)
#         return cos(π / β * ((tInL + tOutL) - (tInR + tOutR)))
#     else
#         return cos(π / β * ((tInL - tOutL) + (tInR - tOutR)))
#     end
# end

function eval(para, ver4::Ver4, varK, varT, legK, evalG::Function, evalV::Function, fast = false; kwargs...)
    KinL, KoutL, KinR, KoutR = legK
    spin = para.spin
    T0idx = para.firstTauIdx
    Kidx = para.firstLoopIdx + ver4.loopidxOffset

    if ver4.loopNum == 0
        qd = KinL - KoutL
        qe = KinL - KoutR
        if para.interactionTauNum == 1
            sign = para.isFermi ? -1 : 1
            ver4.weight[1].d = -evalV(qd)
            ver4.weight[1].e = (-evalV(qe)) * sign
        else
            Tidx = para.firstTauIdx + ver4.TidxOffset
            τIn, τOut = varT[Tidx], varT[Tidx+1]
            error("not implemented!")
            # elseif ver4.interactionTauNum == 2
            # vd, wd, ve, we = vertexDynamic(para, qd, qe, τIn, τOut)

            # ver4.weight[1].d = vd
            # ver4.weight[1].e = ve
            # ver4.weight[2].d = wd
            # ver4.weight[2].e = 0.0
            # ver4.weight[3].d = 0.0
            # ver4.weight[3].e = we
        end
        return
    end

    # LoopNum>=1
    for w in ver4.weight
        w.d, w.e = 0.0, 0.0 # initialize all weights
    end
    G = ver4.G
    K = varK[:, Kidx]

    evalAllG!(G[1], K, T0idx, varT, evalG, kwargs...)

    PhaseFactor = 1.0 / (2π)^para.loopDim
    Kt, Ku, Ks = similar(K), similar(K), similar(K) #Kt, Ku and Ks will be re-created later, slow in performance

    for c in ver4.chan
        if c == T
            @. Kt = KoutL + K - KinL
            evalAllG!(G[Int(c)], Kt, T0idx, varT, evalG, kwargs...)
            # println("initializating", G[Int(c)])
        elseif c == U
            # can not be in box!
            @. Ku = KoutR + K - KinL
            evalAllG!(G[Int(c)], Ku, T0idx, varT, evalG, kwargs...)
        elseif c == S
            # S channel, and cann't be in box!
            @. Ks = KinL + KinR - K
            evalAllG!(G[Int(c)], Ks, T0idx, varT, evalG, kwargs...)
        else
            error("not impossible!")
        end
    end
    for b in ver4.bubble
        c = b.chan
        Factor = SymFactor[Int(c)] * PhaseFactor
        if para.isFermi == false
            Factor = abs(Factor)
        end

        if c == T
            eval(para, b.Lver, varK, varT, [KinL, KoutL, Kt, K], evalG, evalV; kwargs...)
            eval(para, b.Rver, varK, varT, [K, Kt, KinR, KoutR], evalG, evalV; kwargs...)
        elseif c == U
            eval(para, b.Lver, varK, varT, [KinL, KoutR, Ku, K], evalG, evalV; kwargs...)
            eval(para, b.Rver, varK, varT, [K, Ku, KinR, KoutL], evalG, evalV; kwargs...)
        elseif c == S
            # S channel
            eval(para, b.Lver, varK, varT, [KinL, Ks, KinR, K], evalG, evalV; kwargs...)
            eval(para, b.Rver, varK, varT, [K, KoutL, Ks, KoutR], evalG, evalV; kwargs...)
        else
            error("not implemented")
        end

        rN = length(b.Rver.weight)
        gWeight = 0.0
        for (l, Lw) in enumerate(b.Lver.weight)
            for (r, Rw) in enumerate(b.Rver.weight)
                map = b.map[(l-1)*rN+r]

                gWeight = map.G0.weight * map.Gx.weight * Factor

                if fast && ver4.level == 1
                    # gWeight *= phase(varT, ver4.Tpair[map.ver])
                    w = ver4.weight[1]
                else
                    w = ver4.weight[map.vidx]
                end

                if c == T
                    w.d += gWeight * (Lw.d * Rw.d * spin + Lw.d * Rw.e + Lw.e * Rw.d)
                    w.e += gWeight * Lw.e * Rw.e
                elseif c == U
                    w.d += gWeight * Lw.e * Rw.e
                    w.e += gWeight * (Lw.d * Rw.d * spin + Lw.d * Rw.e + Lw.e * Rw.d)
                elseif c == S
                    # S channel,  see the note "code convention"
                    w.d += gWeight * (Lw.d * Rw.e + Lw.e * Rw.d)
                    w.e += gWeight * (Lw.d * Rw.d + Lw.e * Rw.e)
                else
                    error("not implemented")
                end

            end
        end

    end
end


text:
function evaluate_subset(df, ex::Expr)
    if ex.head == :call
        return Expr(ex.head, ex.args[1], (evaluate_subset(df, ex.args[i]) for i in 2:length(ex.args))...)
    else
        return Expr(ex.head, (evaluate_subset(df, ex.args[i]) for i in 1:length(ex.args))...)
    end
end
evaluate_subset(df, ex::Symbol) = df[!, ex]
evaluate_subset(df, ex)  = ex


text:
import Base: broadcast, broadcast!

function broadcasted(::typeof(+), dest::CuArray{T}, srcs::AbstractCuArray{T}...) where T
    for src in srcs
        dest === src && continue
        @assert ndims(dest) >= ndims(src)
        _broadcast!(+, dest, src)
    end
    dest
end

function _broadcast!(::typeof(+), dest::CuArray{T}, src::CuArray{T}) where T
    if length(dest) == length(src)
        BLAS.axpy!(T(1), src, dest)
    else
        CUDNN.addto!(1, src, 1, dest)
    end
    dest
end

@generated function _broadcast!(::typeof(+), dest::CuArray{T,N}, src::AbstractCuArray{T}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void broadcast_add(Array<$Ct,$N> dest, Array<$Ct,$N> src) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= dest.length()) return;

        int sub[$N];
        dest.ind2sub(sub, idx);
        dest(sub) += src(sub);
    }""")
    quote
        gdims, bdims = cudims(length(dest))
        $k(gdims, bdims, dest, src)
        dest
    end
end

function broadcast!(::typeof(*), y::CuArray{T}, x1::CuArray{T}, x2::CuArray{T}) where T
    throw("Not implemented yet.")
end