text:
using Flux
using Base.Test
using DataFlow
using Keras

Keras.new_type(a, b) = Keras.new_type(a, nothing, b)

include("pool_tests.jl")
include("activation_test.jl")
include("normalize_tests.jl")

ext: jl


text:
module Day24

using AdventOfCode2019

function day24(input::String = readInput(joinpath(@__DIR__, "..", "data", "day24.txt")))
    layout = parse_input(input)
    p1 = part1(layout)
    p2 = part2(layout)
    return [p1, p2]
end

function part1(layout::Array{Bool,2})
    state = copy(layout)
    tmp = similar(state)
    ratings = [biodiversity_rating(state)]
    while true
        next_state!(tmp, state)
        copy!(state, tmp)
        br = biodiversity_rating(state)
        if br in ratings
            return br
        end
        push!(ratings, br)
    end
end

function part2(layout::Array{Bool,2}; niter::Int = 200)
    layers = initialize_layers(layout)
    for i = 1:niter
        next_state!(layers)
    end
    return sum.(values(layers)) |> sum
end

function parse_input(input::String)
    layout = zeros(Bool, 5, 5)
    sp = split(input, "\n")
    for (i, row) in enumerate(sp)
        for (j, c) in enumerate(row)
            if sp[i][j] == '#'
                layout[i, j] = true
            end
        end
    end
    return layout
end

function in_bounds(i::Int, j::Int, layout::Array{Bool,2})
    return i > 0 && j > 0 && i <= size(layout, 1) && j <= size(layout, 2)
end

function number_of_infested_neighbours(i::Int, j::Int, layout::Array{Bool,2}; exclude_middle = false)
    n = 0
    for (k, l) in ((-1, 0), (1, 0), (0, -1), (0, 1))
        if in_bounds(i + k, j + l, layout) && layout[i + k, j + l]
            if exclude_middle && ceil.(Int, size(layout) ./ 2) == (i + k, j + l)
                continue
            end
            n += 1
        end
    end
    return n
end

function number_of_infested_neighbours(i::Int, j::Int, l::Int, layers::Dict{Int,Array{Bool,2}})
    if i == 3 && j == 3
        return 0
    end
    n = number_of_infested_neighbours(i, j, layers[l]; exclude_middle = true)
    if i == 1 && haskey(layers, l - 1) && layers[l - 1][2, 3]
        n += 1
    end
    if j == 1 && haskey(layers, l - 1) && layers[l - 1][3, 2]
        n += 1
    end
    if i == 5 && haskey(layers, l - 1) && layers[l - 1][4, 3]
        n += 1
    end
    if j == 5 && haskey(layers, l - 1) && layers[l - 1][3, 4]
        n += 1
    end
    if i == 2 && j == 3 && haskey(layers, l + 1)
        for j = 1:5
            if layers[l + 1][1, j]
                n += 1
            end
        end
    end
    if i == 4 && j == 3 && haskey(layers, l + 1)
        for j = 1:5
            if layers[l + 1][5, j]
                n += 1
            end
        end
    end
    if i == 3 && j == 2 && haskey(layers, l + 1)
        for i = 1:5
            if layers[l + 1][i, 1]
                n += 1
            end
        end
    end
    if i == 3 && j == 4 && haskey(layers, l + 1)
        for i = 1:5
            if layers[l + 1][i, 5]
                n += 1
            end
        end
    end
    return n
end

function next_state!(dest::Array{Bool,2}, source::Array{Bool,2})
    copy!(dest, source)
    for i = 1:size(source, 1)
        for j = 1:size(source, 2)
            n = number_of_infested_neighbours(i, j, source)
            if source[i, j]
                if n != 1
                    dest[i, j] = false
                end
            else
                if n == 1 || n == 2
                    dest[i, j] = true
                end
            end
        end
    end
end

function next_state!(layers::Dict{Int,Array{Bool,2}})
    add_two_layers!(layers)
    source = deepcopy(layers)
    for l in keys(source)
        for i = 1:5
            for j = 1:5
                n = number_of_infested_neighbours(i, j, l, source)
                if source[l][i, j]
                    if n != 1
                        layers[l][i, j] = false
                    end
                else
                    if n == 1 || n == 2
                        layers[l][i, j] = true
                    end
                end
            end
        end
    end
end

function biodiversity_rating(layout::Array{Bool,2})
    rating = 0
    for i in eachindex(layout)
        if layout'[i]
            rating += 2^(i - 1)
        end
    end
    return rating
end

function initialize_layers(layout::Array{Bool,2})
    layers = Dict{Int,Array{Bool,2}}()
    layers[0] = copy(layout)
    return layers
end

function add_two_layers!(layers::Dict{Int,Array{Bool,2}})
    mi = minimum(keys(layers))
    ma = maximum(keys(layers))
    layers[mi - 1] = zeros(Bool, size(layers[mi]))
    layers[ma + 1] = zeros(Bool, size(layers[ma]))
end

end # module


text:
module ELMNIST

using  ExtremeLearningMachine
import MLDatasets
import Random

include("dataset.jl")
include("model.jl")

export ELM,
       Linear, ReLU, Sigmoid, Square, Tanh,
       add_data!, solve, predict,
       prepare_model, calculate_error,
       download, testdata, traindata

end # module


text:
@testset verbose = true "Models" begin
    @testset "Poisson" begin
        mark_dist = Dists.Categorical([0.1, 0.3, 0.6])
        pp = PoissonProcess(5.0, mark_dist)
        h = rand(pp, 0.0, 100.0)
        pp_est = fit(PoissonProcess{Dists.Categorical}, h)
        error = Dists.probs(mark_distribution(pp_est)) - Dists.probs(mark_distribution(pp))
        @test maximum(error) < 0.1
    end
end


text:
################
# Division
################

"""
    Base.:/(a::AbstractUncertainValue, b::AbstractUncertainValue) -> UncertainValue

Division operator for pairs of uncertain values.

Computes the element-wise quotients between for a default of `n = 10000` realizations of `a` and
`b`, then returns an uncertain value based on a kernel density estimate to the distribution
of the element-wise quotients.

Use the `/(a, b, n)` syntax to tune the number (`n`) of draws.
"""
function Base.:/(a::AbstractUncertainValue, b::AbstractUncertainValue; n::Int = 30000)
    UncertainValue(resample(a, n) ./ resample(b, n))
end

"""
    Base.:/(a::Real, b::AbstractUncertainValue) -> UncertainValue

Division operator for between scalars and uncertain values.

Computes the element-wise quotients between `a` and `b` for a default of `n = 10000` realizations
of `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise quotients.

Use the `/(a, b, n)` syntax to tune the number (`n`) of draws.
"""
Base.:/(a::Real, b::AbstractUncertainValue; n::Int = 30000) =
    UncertainValue(a ./ resample(b, n))

"""
    Base.:/(a::AbstractUncertainValue, b::Real) -> UncertainValue

Division operator for between uncertain values and scalars.

Computes the element-wise quotients between `a` and `b` for a default of `n = 10000` realizations
of `a`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise quotients.

Use the `/(a, b, n)` syntax to tune the number (`n`) of draws.
"""
Base.:/(a::AbstractUncertainValue, b::Real; n::Int = 30000) =
    UncertainValue(resample(a, n) ./ b)

"""
    Base.:/(a::AbstractUncertainValue, b::AbstractUncertainValue, n::Int) -> UncertainValue

Division operator for pairs of uncertain values.

Computes the element-wise quotients between `a` and `b` for `n` realizations
of `a` and `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise quotients.

Call this function using the `/(a, b, n)` syntax.
"""
function Base.:/(a::AbstractUncertainValue, b::AbstractUncertainValue, n::Int)
    UncertainValue(resample(a, n) ./ resample(b, n))
end
"""
    Base.:/(a::Real, b::AbstractUncertainValue, n::Int) -> UncertainValue

Division operator for scalar-uncertain value pairs.

Computes the element-wise quotients between `a` and `b` for `n` realizations
of `b`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise quotients.

Call this function using the `/(a, b, n)` syntax.
"""
Base.:/(a::Real, b::AbstractUncertainValue, n::Int) =
    UncertainValue(a ./ resample(b, n))

"""
    Base.:/(a::AbstractUncertainValue, b::Real, n::Int) -> UncertainValue

Division operator for scalar-uncertain value pairs.

Computes the element-wise quotients between `a` and `b` for `n` realizations
of `a`, then returns an uncertain value based on a kernel density estimate to the
distribution of the element-wise quotients.

Call this function using the `/(a, b, n)` syntax.
"""
Base.:/(a::AbstractUncertainValue, b::Real, n::Int) =
    UncertainValue(resample(a, n) ./ b)

#####################################################################################
# Special cases
#####################################################################################

import ..UncertainValues: CertainValue

##################
# `CertainValue`s
#################
"""
    Base.:/(a::Union{CertainValue, Real}, b::Union{CertainValue, Real})

Division of certain values with themselves or scalars acts as regular division, but
returns the result wrapped in a `CertainValue` instance.
"""
Base.:/(a::Union{CertainValue, Real}, b::Union{CertainValue, Real})

Base.:/(a::CertainValue, b::CertainValue) = CertainValue(a.value / b.value)
Base.:/(a::CertainValue, b::Real) = CertainValue(a.value / b)
Base.:/(a::Real, b::CertainValue) = CertainValue(a / b.value)


text:
using Test

@testset "alltests" begin
  include("typeinfer.jl")
  include("ciid.jl")
  include("condition.jl")
  include("distributions.jl")
  include("core.jl")
  include("intervene.jl")
  include("logpdf.jl")
  include("multivariate.jl")
  include("namedtuple.jl")
  include("solution.jl")
  include("tagging.jl")
  include("typeinfer.jl")
  include("var.jl")
end


text:
using FLoops
using StaticArrays, StatsBase, Combinatorics

function birthday_problem_floop(t, ncores )

v = 0
  @floop ThreadedEx(basesize=t÷ncores) for _ in 1:t
       months =  @SVector [rand(1:12) for i in 1:20]
       counts =  @SVector [sum(months.==i) for i=1:12]
       success = sum(counts.==2) == 4  &&  sum(counts.==3)==4
       @reduce(v += success)
  end

return v/t

end

 birthday_problem_floop(10_000_000,4)

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
##

using MAT

"""
    read_SPM(pth::String,fil::String="interp_precomputed.mat")

Reads pre-computed interpolation (sparse matrix) from
`pth*"interp_precomputed.mat"`.
"""
function read_SPM(pth::String,fil::String="interp_precomputed.mat")
    #vars = matread(pth*"interp_precomputed.mat")
    file = matopen(pth*fil)
    interp=read(file, "interp")
    lon=read(file, "lon")
    lat=read(file, "lat")
    SPM=interp["SPM"]
    #println(keys(interp))
    close(file)
    return SPM,lon,lat
end

##

using JLD

"""
    prep_MTRX()

Repackage interpolation matrix, mask, etc to `.jld` file.
"""
function prep_MTRX()
    GCMGridSpec()
    GCMGridLoad()
    msk2d=mask(view(MeshArrays.hFacC,:,:,1),NaN,0)
    msk3d=mask(MeshArrays.hFacC,NaN,0)
    msk2d=convert2gcmfaces(msk2d)
    msk3d=convert2gcmfaces(msk3d)

    dirIn=""
    MTRX,lon,lat=read_SPM(dirIn)
    lon=vec(lon[:,1])
    lat=vec(lat[1,:])

    fid = open("GRID_LLC90/RC.data")
    dep=Array{Float64,1}(undef,50)
    read!(fid,dep)
    close(fid)
    dep = -hton.(dep)

    siz2d=(length(lon),length(lat))
    siz3d=(length(lon),length(lat),50)

    save(dirIn*"MTRX.jld", "MTRX", MTRX, "lon", lon, "lat", lat, "dep", dep,
    "msk2d", msk2d, "msk3d", msk3d, "siz2d", siz2d, "siz3d", siz3d)
end


text:
# BED Writer
# ==========

"""
    BED.Writer(output::IO)

Create a data writer of the BED file format.

# Arguments:
* `output`: data sink
"""
struct Writer <: BioCore.IO.AbstractWriter
    output::IO
end

function BioCore.IO.stream(writer::Writer)
    return writer.output
end

function Base.write(writer::Writer, record::Record)
    return write(writer.output, record, '\n')
end


text:
export mysum1!,
       musum2!,
       dotproduct!,
       dotproduct1!,
       OneStepForward!,
       OneStepAdjoint!,
       MultiStepForward,
       MultiStepAdjoint,
       MultiStepAdjoint_spt

include("OneStepForwardAdjoint.jl")
include("MultiStepForward.jl")
include("MultiStepAdjoint.jl")

max_stars_repo_name: WneleiGao/ElasticWave


text:
include("CauchyWeight.jl")
include("Geometry.jl")
include("evaluation.jl")
include("skewProductFun.jl")
include("convolutionProductFun.jl")

# GreensFun

export GreensFun

# In GreensFun, K must be <: BivariateFun, since a kernel
# could consist of a Vector of ProductFun's and LowRankFun's.
# If there are GreensFun's in the Vector, then we recurse.

struct GreensFun{K<:BivariateFun,T} <: BivariateFun{T}
    kernels::Vector{K}
    function GreensFun{K,T}(Kernels::Vector{K}) where {K,T}
        if greensfun_checkdomains(Kernels)
            if greensfun_checkgreensfun(Kernels)
                return GreensFun(vcat(map(kernels,Kernels)...))
            end
            new{K,T}(Kernels)
        else
            error("Cannot create GreensFun: all kernel domains must equal $(domain(Kernels[1]))")
        end
    end
end

function greensfun_checkdomains(Kernels::Vector{K}) where K
    d = domain(Kernels[1])
    ret = true
    for i in 2:length(Kernels)
        ret *= (domain(Kernels[i]) == d)
    end
    ret
end

function greensfun_checkgreensfun(Kernels::Vector{K}) where K
    ret = true
    for i in 1:length(Kernels)
        ret *= !(typeof(Kernels[i]) <: GreensFun)
    end
    !ret
end

GreensFun(kernels::Vector{K}) where {K<:MultivariateFun} =
    GreensFun{eltype(kernels),mapreduce(cfstype,promote_type,kernels)}(kernels)

GreensFun(F::K) where {K<:MultivariateFun} = GreensFun(K[F])

convert(::Type{GreensFun{K,T}},A::GreensFun{K,T}) where {K<:BivariateFun,T} = A
convert(::Type{GreensFun{K,T}},A::GreensFun) where {K<:BivariateFun,T} = error("Cannot convert GreensFun")

Base.length(G::GreensFun) = length(G.kernels)
Base.transpose(G::GreensFun) = GreensFun(mapreduce(transpose,+,G.kernels))
convert(::Type{GreensFun},F::Union{ProductFun,LowRankFun}) = GreensFun(F)
rank(G::GreensFun) = error("Not all kernels are low rank approximations.")

domain(G::GreensFun) = domain(first(G.kernels))
(G::GreensFun)(x,y)=evaluate(G,x,y)
evaluate(G::GreensFun,x,y) = mapreduce(f->evaluate(f,x,y),+,G.kernels)
kernels(B::BivariateFun) = B
kernels(G::GreensFun) = G.kernels

rank(G::GreensFun{L}) where {L<:LowRankFun} = mapreduce(rank,+,G.kernels)
slices(G::GreensFun{L}) where {L<:LowRankFun} = mapreduce(x->x.A,vcat,G.kernels),mapreduce(x->x.B,vcat,G.kernels)
slices(G::GreensFun{L},k::Int) where {L<:LowRankFun} = slices(G)[k]
LowRankIntegralOperator(G::GreensFun{L}) where {L<:LowRankFun} = LowRankIntegralOperator(slices(G)...)

Base.promote_rule(::Type{GreensFun{K,T}},::Type{GreensFun{K1,T1}}) where {K,T,K1,T1} = GreensFun{promote_type(K,K1),promote_type(T,T1)}

getindex(⨍::Operator,G::GreensFun) = mapreduce(f->⨍[f],+,G.kernels)

# avoid ambiguity
for TYP in (:(ApproxFun.DefiniteLineIntegralWrapper),:DefiniteLineIntegral)
    @eval function getindex(⨍::$TYP,B::Matrix{F}) where F<:BivariateFun
        m,n = size(B)
        wsp = domainspace(⨍)
        @assert m == length(wsp.spaces)
        ⨍j = DefiniteLineIntegral(component(wsp,1))
        ret = Array{Operator{promote_type(eltype(⨍j),map(cfstype,B)...)}}(undef,m,n)
        for j=1:n
            ⨍j = DefiniteLineIntegral(component(wsp,j))
            for i=1:m
                ret[i,j] = ⨍j[B[i,j]]
            end
        end
        ops=promotespaces(ret)
        InterlaceOperator(ops,
                          PiecewiseSpace(map(domainspace,ops[1,:])),
                          PiecewiseSpace(map(rangespace,ops[:,1])))
    end
end

# Algebra with BivariateFun's

+(F::GreensFun,G::GreensFun) = GreensFun([F.kernels;G.kernels])
-(F::GreensFun,G::GreensFun) = GreensFun([F.kernels;-G.kernels])
+(G::GreensFun,B::BivariateFun) = GreensFun([G.kernels;kernels(B)])
-(G::GreensFun,B::BivariateFun) = GreensFun([G.kernels;-kernels(B)])
+(B::BivariateFun,G::GreensFun) = GreensFun([kernels(B);G.kernels])
-(B::BivariateFun,G::GreensFun) = GreensFun([kernels(B);-G.kernels])

# work around 0.4 bug

# Custom operations on Arrays required to infer type of resulting Array{GreensFun}

for op in (:+,:-)
    @eval begin
        function $op(A::Array{GreensFun{F,T1}},B::Array{GreensFun{G,T2}}) where {F<:BivariateFun,T1,G<:BivariateFun,T2}
            C = similar(A, GreensFun{promote_type(F,G),promote_type(T1,T2)}, promote_shape(size(A),size(B)))
            for i=1:length(A)
                @inbounds C[i] = $op(A[i],B[i])
            end
            return C
        end
    end
end

## TODO: Get ProductFun & LowRankFun in different CauchyWeight spaces to promote to GreensFun.
#=
+{S<:UnivariateSpace,V<:UnivariateSpace,O1,O2}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,CauchyWeight{O2}}) = GreensFun([F,G])
+{S<:UnivariateSpace,V<:UnivariateSpace,O1,SS}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,SS}) = GreensFun([F,G])
+{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,O1}(F::ProductFun{S,V,SS},G::ProductFun{S,V,CauchyWeight{O1}}) = GreensFun([F,G])

-{S<:UnivariateSpace,V<:UnivariateSpace,O1,O2}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,CauchyWeight{O2}}) = GreensFun([F,-G])
-{S<:UnivariateSpace,V<:UnivariateSpace,O1,SS}(F::ProductFun{S,V,CauchyWeight{O1}},G::ProductFun{S,V,SS}) = GreensFun([F,-G])
-{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,O1}(F::ProductFun{S,V,SS},G::ProductFun{S,V,CauchyWeight{O1}}) = GreensFun([F,-G])
=#

## Constructors

GreensFun(f::Function,args...;kwds...) = GreensFun(dynamic(f),args...;kwds...)
GreensFun(f::Function,g::Function,args...;kwds...) = GreensFun(dynamic(f),dynamic(g),args...;kwds...)

function GreensFun(f::DFunction,ss::SS;method::Symbol=:lowrank,kwds...) where SS<:AbstractProductSpace
    if method == :standard
        F = ProductFun(f,ss,kwds...)
    elseif method == :convolution
        F = convolutionProductFun(f,ss,kwds...)
    elseif method == :unsplit
        # Approximate imaginary & smooth part.
        F1 = skewProductFun((x,y)->imag(f(x,y)),ss.space;kwds...)
        # Extract diagonal value.
        d = domain(ss)
        xm,ym = mean([first(d[1]),last(d[1])]),mean([first(d[2]),last(d[2])])
        F1m = F1(xm,ym)
        # Set this normalized part to be the singular part.
        F1 = ProductFun(-coefficients(F1)/F1m/2,ss)
        # Approximate real & smooth part after singular extraction.
        m,n = size(F1)
        if typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Chebyshev,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n)+1)
        elseif typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Laurent,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n))
        end
        F = [F1;F2]
    elseif method == :lowrank
        F = LowRankFun(f,ss;method=:standard,kwds...)
    elseif method == :Cholesky
        F = LowRankFun(f,ss;method=method,kwds...)
    end
    GreensFun(F)
end

function GreensFun(f::DFunction,g::DFunction,ss::SS;method::Symbol=:unsplit,kwds...) where SS<:AbstractProductSpace
    if method == :unsplit
        # Approximate Riemann function of operator.
        G = skewProductFun(g,ss.space;kwds...)
        # Normalize and set to the singular part.
        F1 = ProductFun(-coefficients(G)/2,ss)
        # Approximate real & smooth part after singular extraction.
        m,n = size(F1)
        if typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Chebyshev,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n)+1)
        elseif typeof(ss.space) <: TensorSpace && all(k->typeof(ss.space.spaces[k]) <: Laurent,1:2)
            F2 = skewProductFun((x,y)->f(x,y) - F1(x,y),ss.space,nextpow2(m),nextpow2(n))
        end
        F = [F1;F2]
    end
    GreensFun(F)
end

# Array of GreensFun on TensorSpace of PiecewiseSpaces

function GreensFun(f::DFunction,ss::AbstractProductSpace{Tuple{PWS1,PWS2}};method::Symbol=:lowrank,tolerance::Symbol=:absolute,kwds...) where {PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}
    M,N = ncomponents(factor(ss,1)),ncomponents(factor(ss,2))
    @assert M == N
    G = Array{GreensFun}(undef,N,N)
    if method == :standard
        for i=1:N,j=1:N
            G[i,j] = GreensFun(f,component(ss,i,j);method=method,kwds...)
        end
    elseif method == :convolution
        for i=1:N,j=i:N
            G[i,j] = GreensFun(f,component(ss,i,j);method=method,kwds...)
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    elseif method == :unsplit
        maxF = Array{Number}(undef,N)
        for i=1:N
          G[i,i] = GreensFun(f,component(ss,i,i);method=method,kwds...)
          maxF[i] = one(real(mapreduce(cfstype,promote_type,G[i,i].kernels)))/2π
        end
        for i=1:N,j=i+1:N
            G[i,j] = GreensFun(f,component(ss,i,j).space;method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    elseif method == :lowrank
        for i=1:N,j=1:N
            G[i,j] = GreensFun(f,component(ss,i,j);method=method,kwds...)
        end
    elseif method == :Cholesky
        if tolerance == :relative
            for i=1:N
                G[i,i] = GreensFun(f,component(ss,i,i);method=method,kwds...)
                for j=i+1:N
                    G[i,j] = GreensFun(f,component(ss,i,j);method=:lowrank,kwds...)
                end
            end
        elseif tolerance == :absolute
            maxF = Array{Number}(undef,N)
            for i=1:N
                F,maxF[i] = LowRankFun(f,component(ss,i,i);method=method,retmax=true,kwds...)
                G[i,i] = GreensFun(F)
            end
            for i=1:N,j=i+1:N
                G[i,j] = GreensFun(f,component(ss,i,j);method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
            end
        end
        for i=2:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    end
    G
end

function GreensFun(f::DFunction,g::DFunction,ss::AbstractProductSpace{Tuple{PWS1,PWS2}};method::Symbol=:unsplit,tolerance::Symbol=:absolute,kwds...) where {PWS1<:PiecewiseSpace,PWS2<:PiecewiseSpace}
    M,N = ncomponents(factor(ss.space,1)),ncomponents(factor(ss.space,2))
    @assert M == N
    G = Array{GreensFun}(undef,N,N)
    if method == :unsplit
        maxF = Array{Number}(undef,N)
        for i=1:N
            G[i,i] = GreensFun(f,g,ss[i,i];method=method,kwds...)
            maxF[i] = one(real(mapreduce(cfstype,promote_type,G[i,i].kernels)))/2π
        end
        for i=1:N
            for j=i+1:N
                G[i,j] = GreensFun(f,ss[i,j].space;method=:lowrank,tolerance=(tolerance,max(maxF[i],maxF[j])),kwds...)
            end
        end
        for i=1:N,j=1:i-1
            G[i,j] = transpose(G[j,i])
        end
    end
    G
end


text:
using Test
using Distributions
using DataStructures

import GlobalSensitivityAnalysis: ishigami, sample

##
## 1. Utils
##

# DeltaData
N = 100
parameters = OrderedDict(
    :param1 => Normal(0,1),
    :param2 => LogNormal(5, 20),
    :param3 => TriangularDist(0, 4, 1)
)

data1 = DeltaData(params = parameters, N = N)
@test data1.params[:param1] == parameters[:param1]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param3] == parameters[:param3]
@test data1.N == N

data2 = DeltaData()
@test data2.params == nothing
@test data2.N == 1000

data3 = DeltaData(params = parameters)
data4 = DeltaData(N = 100)

##
## 2. Sample Delta
##

samples = sample(data1)
@test size(samples, 2) == length(data1.params)
@test size(samples, 1) == data1.N

samples3 = sample(data3)
@test size(samples3, 2) == length(data3.params)
@test size(samples3, 1) == data3.N

@test_throws ErrorException sample(data2) # params are nothing
@test_throws ErrorException sample(data4) # params are nothing

# ##
# ## 3a. Analyze Delta - Indices (TODO is there anything else we should test indices for?)
# ##

Y1 = ishigami(samples)
results = analyze(data1, samples, Y1)
for Si in results[:firstorder]
    @test Si <= 1
end
for Si in results[:delta]
    @test Si <= 1
end
for CI in results[:firstorder_conf]
    @test CI > 0
end
for CI in results[:delta_conf]
    @test CI > 0
end

Y3 = ishigami(samples3)
results = analyze(data3, samples3, Y3)
for Si in results[:firstorder]
    @test Si <= 1
end
for Si in results[:delta]
    @test Si <= 1
end
for CI in results[:firstorder_conf]
    @test CI > 0
end
for CI in results[:delta_conf]
    @test CI > 0
end

# ##
# ## 3b. Analyze Delta - Optional Keyword Args
# ##

data = DeltaData(
    params = OrderedDict(:x1 => Normal(1, 0.2),
        :x2 => Uniform(0.75, 1.25),
        :x3 => LogNormal(0, 0.5)),
    N = 1000
)
samples = sample(data)
Y = ishigami(samples)
results = analyze(data, samples, Y)

results = analyze(data, samples, Y; progress_meter = false) # no progress bar should show
results = analyze(data, samples, Y; num_resamples = 100, conf_level = 0.90, progress_meter = false) # keywords

@test length(analyze(data, samples, Y; N_override = 10)) == 5
results_override = analyze(data, samples, Y; N_override = data.N)
results_original = analyze(data, samples, Y)
@test results_override[:firstorder] == results_original[:firstorder]
@test_throws ErrorException analyze(data, samples, Y; N_override = data.N + 1) # N_override > N

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
"""
Convenience code to deal with finding files, reading parameters, etc.

Author: Eric Ford
Created: Sept 2020
"""

default_paths_to_search = [".",joinpath(".","examples"),pwd(),joinpath(pwd(),"examples") ]

""" `read_data_paths( ; path_to_search, filename, verbose )`
Looks for data_paths.jl and includes it to set data_paths

Warning:  Malicious users could insert arbitrary code into data_paths.jl.  Don't be a malicous user.
"""
function read_data_paths(; paths_to_search::Union{String,AbstractVector{String}} = default_paths_to_search, filename::String = "data_paths.jl", verbose::Bool = true)
   if verbose   println("# Looking for data paths and config files in default_paths_to_search.")    end
   idx_path = findfirst(isfile,map(d->joinpath(d,filename),paths_to_search))
   if isnothing(idx_path)
      @error("Can't find ", filename, "  in any of: ", paths_to_search)
   end
   data_paths_jl = paths_to_search[idx_path]
   if verbose   println("# Found ", data_paths_jl, " directory.")   end
   if !isdir(data_paths_jl)
         @error("Can't access instrument's base data directory ", data_paths_jl, ".")
   end

   path_to_data_paths = joinpath(pwd(),data_paths_jl,filename)
   if isfile(path_to_data_paths)
      code_to_include_param = quote
         include($path_to_data_paths)
      end
   else
      println("# Did not locate ", path_to_data_paths )
   end
   return code_to_include_param
end

#=
""" make_manifest(data_path::String, target_subdir::String, Inst::Module; [opts] )
Returns a dataframe containing a list of files to be read and some metadata (e.g., observation times)

# Optional arguements
- verbose = true

"""
function make_manifest(data_path::String, target_subdir::String, Inst::Module; verbose::Bool = true)
   #=
   if Inst == EXPRES
      #@assert isdefined(Main,:expres_data_path)
      data_path = joinpath(expres_data_path,target_subdir)
   elseif Inst == NEID
      #  @assert isdefined(Main,:solar_data_path)
      data_path = joinpath(solar_data_path,target_subdir)
   else
      @error "Specified instrument didn't match a module name, ", string(Inst), "."
   end
   =#
   if !isdir(data_path)
      @error "Can't access target data directory ", data_path, "."
   end

   if verbose  println("# Creating manifest of files to process.")    end
   df_files = Inst.make_manifest(data_path)
   if size(df_files,1) < 1
      @error("Did not find any files in ", data_path ,".")
   end
   #=
   df_files = df_files |> @filter( _.target == target_fits ) |> DataFrame
   if size(df_files,1) < 1
      @error("Did not find any files with target field matching ", target_fits,".")
   end
   if size(df_files,1) > max_spectra_to_use
      df_files = df_files |> @take(max_spectra_to_use) |> DataFrame
   end
   =#
   return df_files
end
=#

"""   `code_to_include_param_jl( ; path_to_search, filename, verbose )`

Returns a Code object.  After `res = code_toread_param_jl( path_to_search )`,
execute `eval(res)` to actually include the param.jl file.
This is useful since it allows variables to be placed into caller's namespace.

Warning:  Malicious users could insert arbitrary code into param.jl.  Don't be a malicous user.
"""
function code_to_include_param_jl(;paths_to_search::Union{String,AbstractVector{String}} = default_paths_to_search, filename::String = "param.jl", verbose::Bool = true)
   if verbose   println("# Looking for param.jl file to set configuration parameters.")    end
   idx_path = findfirst(isfile,map(d->joinpath(d,filename),paths_to_search))
   if isnothing(idx_path)
      @error(" Couldn't find $filename in $default_paths_to_search")
      return Expr()
   end
   data_path = paths_to_search[idx_path]

   path_to_param = joinpath(pwd(),data_path,filename)
   if isfile(path_to_param )
      if verbose    println("# Reading parameter values from ", path_to_param)   end
      code_to_include_param = quote
         include($path_to_param)
      end
   else
      println("# Did not locate ", filename, " in ", data_path, ".") #  Returning code to set default values for lots of things.")
      # Set defaults in case not in param.jl
      code_to_include_param = quote
         #=
         df_files_use = df_files
         espresso_filename = joinpath(pkgdir(RvSpectMLBase),"data","masks","G2.espresso.mas")
         ccf_mid_velocity = 0
         tophap_ccf_mask_scale_factor=1.6
         =#
      end
   end
   return code_to_include_param
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
subject = 600851475143
upperlimit = Int(floor(sqrt(subject)))
candidates = []
primenumbers = []

for i = 3:upperlimit
  if subject % i == 0
    push!(candidates, i)
  end
end

for candidate in candidates
  isprime = true
  for i = 3:Int(floor(sqrt(candidate)))
    if candidate % i == 0
      isprime = false
      break
    end
  end
  if isprime
    push!(primenumbers, candidate)
  end
end
println(maximum(primenumbers))

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
#!/usr/bin/env julia

try
	using LibSerialPort
catch
	println(stderr, "I require LibSerialPort.jl to work.")
	println(stderr, "Open 'julia' and enter into Pkg mode by hitting ]")
	println(stderr, "Then enter 'add LibSerialPort' to install the dependencies.")
	exit(1)
end

function main(args:: AbstractVector{String} = ARGS)
	if isempty(args)
		println("Usage: $(PROGRAM_FILE |> basename) <port> [clearscreen]")
		println("Available ports:")
		list_ports()
		return
	end

	clearscreen = parse(Bool, get(args, 2, "false"))

	rssi_data = Int[]
	framesize_data = Int[]
	frameid_data = Int[]
	bmp_pres_data = Float32[]
	bmp_temp_data = Float32[]
	lm35_temp_data = Float32[]
	dht11_temp_data = Float32[]
	dht11_humi_data = Float32[]

	open(args[1], 115200) do sp
		# in case we start listening in the middle of a frame, discard it
		line = readline(sp)
		while line != "-- END"
			line = readline(sp)
		end

		# start reading frames indefinitely
		while true
			readline(sp) # consume "-- START"

			rssi = parse(Int, readline(sp))
			framesize = parse(Int, readline(sp))
			frameid = parse(Int, readline(sp))
			bmp_pres = parse(Float32, readline(sp))
			bmp_temp = parse(Float32, readline(sp))
			lm35_temp = parse(Float32, readline(sp))
			dht11_temp = parse(Float32, readline(sp))
			dht11_humi = parse(Float32, readline(sp))

			push!(rssi_data, rssi)
			push!(framesize_data, framesize)
			push!(frameid_data, frameid)
			push!(bmp_pres_data, bmp_pres)
			push!(bmp_temp_data, bmp_temp)
			push!(lm35_temp_data, lm35_temp)
			push!(dht11_temp_data, dht11_temp)
			push!(dht11_humi_data, dht11_humi)

			clearscreen && print("\e[2J\e[;H")
			lines = [
				"Frame $frameid | $framesize bytes | RSSI $rssi",
				"BMP: $bmp_pres hPa, $bmp_temp °C",
				"LM35: $lm35_temp °C",
				"DHT11: $dht11_humi%, $dht11_temp °C"
			]
			maxlen = maximum(length, lines)

			println('|', '='^(maxlen+2), '|')
			for l in lines
				println("| ", rpad(l, maxlen), " |")
			end
			println('|', '='^(maxlen+2), '|')
			println()

			readline(sp) # consume "-- END"
		end
	end
end

try
	Base.exit_on_sigint(false)
	main()
catch e
	!isa(e, InterruptException) && rethrow()
	println()
end

max_stars_repo_name: CanOpeners/code-scraps

max_issues_repo_name: CanOpeners/code-scraps

max_forks_repo_name: CanOpeners/code-scraps


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    spheredir(θ, φ)

Returns the 3D direction given polar angle `θ` and
azimuthal angle `φ` in degrees according to the ISO
convention.
"""
function spheredir(theta, phi)
  θ, φ = deg2rad(theta), deg2rad(phi)
  Vec(sin(θ)*cos(φ), sin(θ)*sin(φ), cos(θ))
end

"""
    aniso2distance(semiaxes, angles; convention=TaitBryanExtr)

Return the distance associated with the ellipsoid with given `semiaxes`,
`angles` and `convention`. For more information, see `Meshes.Ellipsoid`.
"""
aniso2distance(semiaxes, angles; convention=TaitBryanExtr) =
  metric(Ellipsoid(semiaxes, angles, convention=convention))


text:
module Jackknife

using Statistics

"""
    Reducer(f)

Type that holds a function and enforces at runtime that the function returns either a
`Real` value or `missing`. This is used internally to sanity check the input for functions
which accept a function argument.
"""
struct Reducer
    f::Function
end
(r::Reducer)(x) = r.f(x)::Union{Real,Missing}

const FuncOrReducer = Union{Function,Reducer}

"""
    Jackknife.leaveoneout(estimator::Function, x::AbstractVector{<:Real})

Estimate the parameter ``\\hat{\\theta}`` using the estimating function `estimator`
for each subsample of `x`, systematically omitting each index one at a time. The
result is a vector of length `length(x)-1` of parameter estimates.
"""
leaveoneout(estimator::Function, x::AbstractVector{T}) where {T<:Union{Real,Missing}} =
    leaveoneout(Reducer(estimator), x)

function leaveoneout(estimator::Reducer, x::AbstractVector{T}) where T<:Union{Real,Missing}
    length(x) > 1 || throw(ArgumentError("The sample must have size > 1"))
    inds = eachindex(x)
    return map(i->estimator(view(x, filter(!isequal(i), inds))), inds)
end

"""
    Jackknife.variance(estimator::Function, x)

Compute the jackknife estimate of the variance of `estimator`, which is given as a
function that computes a point estimate when passed a real-valued vector `x`.
"""
function variance(estimator::FuncOrReducer, x::AbstractVector{T}) where T<:Union{Real,Missing}
    θ = leaveoneout(estimator, x)
    n = length(x)
    return var(θ) * (n - 1)^2 / n
end

"""
    Jackknife.bias(estimator::Function, x)

Compute the jackknife estimate of the bias of `estimator`, which is given as a
function that computes a point estimate when passed a real-valued vector `x`.
"""
function bias(estimator::FuncOrReducer, x::AbstractVector{T}) where {T<:Union{Real,Missing}}
    θ = leaveoneout(estimator, x)
    return (length(x) - 1) * (mean(θ) - estimator(x))
end

"""
    Jackknife.estimate(estimator::Function, x)

Compute the bias-corrected jackknife estimate of the parameter `estimator(x)`.
"""
function estimate(estimator::FuncOrReducer, x::AbstractVector{T}) where T<:Union{Real,Missing}
    θ = leaveoneout(estimator, x)
    n = length(x)
    return n * estimator(x) - (n - 1) * mean(θ)
end

end # module


text:
struct Mean <: MLJBase.AggregationMode end

function (::Mean)(dicts::Array{Dict{Any, Any}})
    valuesDict = Dict()
    grps = keys(dicts[1])
    n = length(dicts)
    for grp in grps
        valSum = 0
        for dict in dicts
            valSum += dict[grp]
        end
        if typeof(valSum)==Int64
            valuesDict[grp] = valSum # When the individual metric value is an integer(eg. for true_positive), then aggregation is via sum
        else
            valuesDict[grp] = valSum/n # When the metric is a float, then it is a rate(eg. for true_positive_rate). So, aggregation is via mean
        end
    end
    return valuesDict
end

function (::Mean)(dict::Dict{Any,Any})
    return dict
end

struct MetricWrapper <: MLJBase.Measure
    measure::MLJBase.Measure
    grp::Symbol
end

function MetricWrapper(measure::MLJBase.Measure; grp=:class)
    MetricWrapper(measure, grp)
end

function MetricWrappers(measures::Array{MLJBase.Measure, 1}; grp=:class)
    wrappedMetrics = []
    for measure in measures
        push!(wrappedMetrics, MetricWrapper(measure; grp=grp))
    end
    return wrappedMetrics
end

function (FM::MetricWrapper)(ŷ, X, y)
    if typeof(ŷ) <: UnivariateFiniteArray
        ŷ = StatsBase.mode.(ŷ)
    end
    if typeof(y) <: UnivariateFiniteArray
        y = StatsBase.mode.(y)
    end
    grps = X[:, FM.grp]
    ft = fair_tensor(categorical(ŷ), categorical(y), categorical(grps))

    valuesDict = Dict()
    valuesDict["overall"] = FM.measure(ft)
    for grp in levels(grps)
        valuesDict[grp] = FM.measure(ft; grp=grp)
    end
    return valuesDict
end

MLJBase.name(::Type{<:MetricWrapper}) = "MetricWrapper"
# MLJBase.target_scitype(::Type{<:MetricWrapper}) = AbstractArray{Multiclass{2},1}
MLJBase.supports_weights(::Type{<:MetricWrapper}) = false # for now
# MLJBase.prediction_type(::Type{<:MetricWrapper}) = :deterministic # Not specifying it to have check_measures false
MLJBase.orientation(::Type{<:MetricWrapper}) = :other # other options are :score, :loss
MLJBase.reports_each_observation(::Type{<:MetricWrapper}) = false
MLJBase.aggregation(::Type{<:MetricWrapper}) = Mean()
MLJBase.is_feature_dependent(::Type{<:MetricWrapper}) = true

# To display the name of actual measure while printing instead of the name MetricWrapper
Base.show(stream::IO, m::MetricWrapper) = print(stream, m.measure)


text:
const defcolors = ["#1F77B4", "#FF7F0E", "#2CA02C", "#D62728", "#9467BD",
                   "#8C564B", "#E377C2", "#7F7F7F", "#BCBD22", "#17BECF"]

sv = IScatterSpectrum.ScatterVolume(230e6, 1e-6, 50000e-9, 0.0)
p  = IScatterSpectrum.Plasma(1.5e11, 2000., 1000., sv)
f  = 5.0:5:5000
freq = [0; f]
tr   = GenericTrace[]
push!(tr, scatter(x=freq,
                  y=[IScatterSpectrum.pwrspec(p, sv);
                     [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                  line_color=defcolors[1],
                  name="Te=2000 K"))
p  = IScatterSpectrum.Plasma(1.5e11, 1000., 1000., sv)
push!(tr, scatter(x=freq,
                  y=[IScatterSpectrum.pwrspec(p, sv);
                     [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                  line_color=defcolors[2],
                  name="Te=1000 K"))
p  = IScatterSpectrum.Plasma(1.5e11, 1500., 1000., sv)
push!(tr, scatter(x=freq,
                  y=[IScatterSpectrum.pwrspec(p, sv);
                     [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                  line_color=defcolors[3],
                  name="Te=1500 K"))
p  = IScatterSpectrum.Plasma(1.5e11, 3000., 1000., sv)
push!(tr, scatter(x=freq,
                  y=[IScatterSpectrum.pwrspec(p, sv);
                     [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                  line_color=defcolors[4],
                  name="Te=3000 K"))
layout = Layout(xaxis_title="Frequency in Hz",
                yaxis_title="Power",
                font_size=18,
                legend_x=0.85,
                title="Collisionless Incoherent Scatter Spectrum")
plt = plot(tr, layout)
PlotlyJS.savefig(plt, "/home/scb/code/julia/pkg/v0.6/IScatterSpectrum/test/CollisionlessPlasma.pdf")
PlotlyJS.savefig(plt, "/home/scb/code/julia/pkg/v0.6/IScatterSpectrum/test/CollisionlessPlasma.png")
tr2 = [tr[k] for k in [1, 2, 4]]
sv = IScatterSpectrum.ScatterVolume(230e6, 1e-6, 50000e-9, 70*π/180)
p  = IScatterSpectrum.Plasma(1.5e11, 2000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[1],
                   line_dash="dash",
                   name="α=70 deg"))
p  = IScatterSpectrum.Plasma(1.5e11, 1000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[2],
                   line_dash="dash",
                   showlegend=false))
p  = IScatterSpectrum.Plasma(1.5e11, 3000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[4],
                   line_dash="dash",
                   showlegend=false))
sv = IScatterSpectrum.ScatterVolume(230e6, 1e-6, 50000e-9, 30*π/180)
p  = IScatterSpectrum.Plasma(1.5e11, 2000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[1],
                   line_dash="dot",
                   name="α=30 deg"))
p  = IScatterSpectrum.Plasma(1.5e11, 1000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[2],
                   line_dash="dot",
                   showlegend=false))
p  = IScatterSpectrum.Plasma(1.5e11, 3000., 1000., sv)
push!(tr2, scatter(x=freq,
                   y=[IScatterSpectrum.pwrspec(p, sv);
                      [IScatterSpectrum.pwrspec(freq, p, sv) for freq in f]],
                   line_color=defcolors[4],
                   line_dash="dot",
                   showlegend=false))
plt = plot(tr2, layout)
PlotlyJS.savefig(plt, "/home/scb/code/julia/pkg/v0.6/IScatterSpectrum/test/CollisionlessPlasmaMagField.pdf")
PlotlyJS.savefig(plt, "/home/scb/code/julia/pkg/v0.6/IScatterSpectrum/test/CollisionlessPlasmaMagField.png")


text:
# Read the image format written by RADMC3D.

# Read the ascii text in `image.out` and parse things into a 3 dimensional matrix (x, y, lambda)

# The first four lines are format information
# iformat # = 1 (2 is local observer)
# im_nx   im_ny #number of pixels in x and y directions
# nlam # number of images at different wavelengths
# pixsize_x  pixsize_y # size of the pixels in cm
# lambda[1] ... lambda[nlam + 1] # wavelengths (um) correspending to images
# pixels, ordered from left to right (increasing x) in the inner loop, and from bottom to top (increasing y) in the outer loop. And wavelength is the outermost loop.

module image

export imread, imToSky, imToSpec, SkyImage, blur, -
export tauread

using ..constants
import Images # The Images.jl package, not affiliated w/ DiskJockey
using Dierckx

import Base.- # extend this for Image

# Define an image type, which can store the data as well as pixel spacing

abstract Image

# RawImage reflects the RADMC convention that both x and y are increasing
# with array index. This means that to display the image as RADMC intends it,
# you must set the first array element to the lower left corner.
type RawImage <: Image
    data::Array{Float64, 3} # [ergs/s/cm^2/Hz/ster]
    pixsize_x::Float64 # [cm]
    pixsize_y::Float64 # [cm]
    lams::Vector{Float64} # [μm]
end

# SkyImage is a holder that has both RA and DEC increasing with array index
# This convention is necessary for the FFT step
# However, to display this image in the traditional sky convention (North up,
# East to the left), you must set the first array element to the lower left
# corner *and* flip the array along the RA axis: `fliplr(data)` or flipdim(data, 2)
type SkyImage <: Image
    data::Array{Float64, 3} # [Jy/pixel]
    ra::Vector{Float64} # [arcsec]
    dec::Vector{Float64} # [arcsec]
    lams::Vector{Float64} # [μm]
end

# TausurfImage is designed to hold the projected output
type TausurfImage <: Image
    data::Array{Float64, 3} # cm above/behind central projected plane of disk
    pixsize_x::Float64 # [cm]
    pixsize_y::Float64 # [cm]
    lams::Vector{Float64} # [μm]
end

# Encapsulates the 3D positions of the tau=1 surface, for all pixels
type Tausurf3D
    data::Array{Float64, 4} # [cm], x, y, z positions for multiple wavelengths
end

function -(img1::SkyImage, img2::SkyImage)
    # @assert img1.lams == img2.lams "Images must have the same wavelengths."
    @assert isapprox(img1.ra, img2.ra) "Images must have same RA coordinates."
    @assert isapprox(img1.dec, img2.dec) "Images must have same DEC coordinates."
    data = img1.data - img2.data
    return SkyImage(data, img1.ra, img1.dec, img1.lams)
end

# SkyImage constructor for just a single frame
SkyImage(data::Matrix{Float64}, ra::Vector{Float64}, dec::Vector{Float64}, lam::Float64) =
SkyImage(reshape(data, tuple(size(data)..., 1)), ra, dec, [lam])

# Read the image file (default=image.out) and return it as an Image object, which contains the fluxes in Jy/pixel,
# the sizes and locations of the pixels in arcseconds, and the wavelengths corresponding to the images
function imread(file="image.out")

    fim = open(file, "r")
    iformat = parse(Int, readline(fim))
    im_nx, im_ny = split(readline(fim))
    im_nx = parse(Int, im_nx)
    im_ny = parse(Int, im_ny)
    nlam = parse(Int, readline(fim))
    pixsize_x, pixsize_y = split(readline(fim))
    pixsize_x = parse(Float64, pixsize_x)
    pixsize_y = parse(Float64, pixsize_y)

    # Read the wavelength array
    lams = Array(Float64, nlam)
    for i=1:nlam
        lams[i] = parse(Float64, readline(fim))
    end

    # Create an array with the proper size, and then read the file into it
    data = Array(Float64, (im_ny, im_nx, nlam))

    # According to the RADMC manual, section A.15, the pixels are ordered
    # left to right (increasing x) in the inner loop, and from bottom to top
    # (increasing y) in the outer loop.
    # Basically, pack the array in order but display with origin=lower.

    # Because of the way an image is stored as a matrix, we actually pack the
    # array indices as data[y, x, lam]
    # radmc3dPy achieves something similar by keeping indices the x,y but
    # swaping loop order (radmcPy/image.py:line 675)
    for k=1:nlam
        readline(fim) # Junk space
        for j=1:im_ny
            for i=1:im_nx
                data[j,i,k] = parse(Float64, readline(fim))
            end
        end
    end

    close(fim)

    # According to the RADMC3D manual, the units are *intensity* [erg cm−2 s−1 Hz−1 ster−1]

    return RawImage(data, pixsize_x, pixsize_y, lams)
end

"""Like imread, but for tausurf. Values where there is no surface set to NaN."""
function tauread(file="image.out")
    fim = open(file, "r")
    iformat = parse(Int, readline(fim))
    im_nx, im_ny = split(readline(fim))
    im_nx = parse(Int, im_nx)
    im_ny = parse(Int, im_ny)
    nlam = parse(Int, readline(fim))
    pixsize_x, pixsize_y = split(readline(fim))
    pixsize_x = parse(Float64, pixsize_x)
    pixsize_y = parse(Float64, pixsize_y)

    # Read the wavelength array
    lams = Array(Float64, nlam)
    for i=1:nlam
        lams[i] = parse(Float64, readline(fim))
    end

    # Create an array with the proper size, and then read the file into it
    data = Array(Float64, (im_ny, im_nx, nlam))

    for k=1:nlam
        readline(fim) # Junk space
        for j=1:im_ny
            for i=1:im_nx
                val = parse(Float64, readline(fim))

                # If RADMC3D signaled that there is no tau=1 surface here, set height to NaN
                if isapprox(val, -1e91)
                    data[j,i,k] = NaN #
                else
                    data[j,i,k] = val
                end
            end
        end
    end

    close(fim)

    # According to the RADMC3D manual, the units are *intensity* [erg cm−2 s−1 Hz−1 ster−1]

    return TausurfImage(data, pixsize_x, pixsize_y, lams)

end

# Assumes dpc is parsecs
function imToSky(img::RawImage, dpc::Float64)

    # The RawImage is oriented with North up and East increasing to the left.
    # this means that for the RawImage, the delta RA array goes from + to -

    # However, the SkyImage actually requires RA (ll) in increasing form.
    # Therefore we flip along the RA axis, fliplr(data) or flipdim(data, 2)

    #println("Min and max intensity ", minimum(img.data), " ", maximum(img.data))
    #println("Pixel size ", img.pixsize_x)
    #println("Steradians subtended by each pixel ",  img.pixsize_x * img.pixsize_y / (dpc * pc)^2)

    #convert from ergs/s/cm^2/Hz/ster to to Jy/ster
    conv = 1e23 # [Jy/ster]

    # Conversion from erg/s/cm^2/Hz/ster to Jy/pixel at 1 pc distance.
    # conv = 1e23 * img.pixsize_x * img.pixsize_y / (dpc * pc)^2

    # Flip across RA dimension
    dataJy = flipdim(img.data, 2) .* conv

    (im_ny, im_nx) = size(dataJy)[1:2] #y and x dimensions of the image

    # The locations of pixel centers in cm
    # if n_x = 16, goes [-7.5, -6.5, ..., -0.5, 0.5, ..., 6.5, 7.5] * pixsize
    xx = ((Float64[i for i=0:im_nx-1] + 0.5) - im_nx/2.) * img.pixsize_x
    yy = ((Float64[i for i=0:im_ny-1] + 0.5) - im_ny/2.) * img.pixsize_y

    # The locations of the pixel centers in relative arcseconds
    # Note both RA and DEC increase with array index.
    ra = xx./(AU * dpc)
    dec = yy./(AU * dpc)

    return SkyImage(dataJy, ra, dec, img.lams)

end

# Following Images.jl, give the number of arcseconds in each dimension on how to Gaussian
# blur the channel maps. Unfortunately only aligned Gaussians are allowed so far, no rotation.
function blur(img::SkyImage, sigma)
    # convert sigma in arcseconds into pixels
    # sigma is a length 2 array with the [sigma_y, sigma_x] blurring scales

    # measure image size in arcsecs
    width = img.ra[end] - img.ra[1] # [arcsec]
    npix = size(img.data)[1]
    nchan = size(img.data)[3]

    println("Image width: $width [arcsec], npix: $npix, nchan: $nchan, sigma: $sigma [arcsec]")

    pixel_arcsec = npix / width #

    println("Pixel_arcsec: $pixel_arcsec")

    sigma *= pixel_arcsec # [pixels]

    println("Pixel sigma: $sigma [pixels]")

    data_blur = Array(Float64, size(img.data)...)
    # go through each channel
    for i=1:nchan
        # Now, load this into an Images.jl frame
        img_Images = Images.Image(img.data[:,:,i])
        data_blur[:,:,i] = Images.imfilter_gaussian(img_Images, sigma)
    end

    return SkyImage(data_blur, img.ra, img.dec, img.lams)

end

# Take an image and integrate all the frames to create a spatially-integrated spectrum
function imToSpec(img::SkyImage)

    # pixels in SkyImage are Jy/ster

    # convert from Jy/str to Jy/pixel using str/pixel
    dRA = abs(img.ra[2] - img.ra[1]) * arcsec
    dDEC = abs(img.dec[2] - img.dec[1]) * arcsec

    # Add up all the flux in the pixels to create the spectrum
    flux = squeeze(sum(img.data .* dRA .* dDEC, (1, 2)), (1,2))
    spec = hcat(img.lams, flux) #First column is wl, second is flux

    return spec
end

# Calculate the integrated line flux
function integrateSpec(spec::Matrix{Float64}, lam0::Float64)
    #First column is wl, second is flux

    wl = spec[:,1]
    fl = spec[:,2]

    # Convert wl to kms
    vs = c_kms * (spec[:,1] .- lam0)/lam0

    if vs[2] - vs[1] < 0
        reverse!(vs)
        reverse!(fl)
    end

    # println(vs)
    # println(fl)

    spl = Spline1D(vs, fl)
    tot = integrate(spl, vs[1], vs[end])

    return tot
end

# Storage for zeroth moment map
type ZerothMoment <: Image
   data::Array{Float64, 2} # [Jy · Hz / pixel]
   ra::Vector{Float64} # [arcsec]
   dec::Vector{Float64} # [arcsec]
end

function convert(::Type{ZerothMoment}, img::SkyImage)
    # Sum along the frequency axis
    data = squeeze(sum(img.data, (3)), 3)
    return ZerothMoment(data, img.ra, img.dec)
end

end #model

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
function initialize(obj::FieldModeHandler, arg0::jboolean, arg1::FieldAbsoluteDate)
    return jcall(obj, "initialize", void, (jboolean, FieldAbsoluteDate), arg0, arg1)
end


text:
abstract type Optim end

mutable struct SGD <: Optim
    """
    確率的勾配降下法
    Stocastic Gradient Descent
    """
    lr
    function SGD(; lr = 0.01)
        new(lr)
    end
end

function update!(self::SGD, params, grads)
    for i = 1:lastindex(params)
        for j = 1:lastindex(params[i])
            params[i][j] = params[i][j] .- self.lr .* grads[i][j] # .-=は破壊的変更
        end
    end
end

mutable struct Adam <: Optim
    """
    Adam
    (http://arxiv.org/abs/1412.6980v8)
    """
    lr
    β_1
    β_2
    iter
    m
    v
    function Adam(; lr = 0.001, β_1 = 0.9, β_2 = 0.999)
        self = new()
        self.lr = lr
        self.β_1 = β_1
        self.β_2 = β_2
        self.iter = 0
        self.m = nothing
        self.v = nothing
        return self
    end
end

function update!(self::Adam, params, grads)
    if self.m === nothing
        self.m, self.v = [], []
        for i = 1:lastindex(params)
            for j = 1:lastindex(params[i])
                push!(self.m, zero(params[i][j]))
                push!(self.v, zero(params[i][j]))
            end
        end
    end
    self.iter += 1
    lr_t = self.lr * sqrt(1.0 - self.β_2^self.iter) / (1.0 - self.β_1^self.iter)

    for i = 1:lastindex(params)
        for j = 1:lastindex(params[i])
            self.m[i] += (1 - self.β_1) .* (grads[i][j] .- self.m[i])
            self.v[i] += (1 - self.β_2) .* (grads[i][j] .^ 2 .- self.v[i])
            params[i][j] -= lr_t .* self.m[i] ./ (sqrt.(self.v[i]) .+ 1e-7)
        end
    end
end


text:
## Definition of new vehicle state, vehicle, scene, and action

# Struct that contains noise for position and velocity (assumed to be in s, t coords)
@with_kw struct Noise
    pos::VecE2 = VecE2(0,0)
    vel::Float64 = 0
end

abstract type NoisyState end

# State for a pedestrian that has noise
@with_kw struct NoisyPedState <: NoisyState
    veh_state::VehicleState # Noise on position and velocity
    noise::Noise # Noise on position and velocity
end

# State for vehicles with blinkers, lane intention and noise
@with_kw struct BlinkerState <: NoisyState
    veh_state::VehicleState # position and velocity
    blinker::Bool # Whether or not the blinker is on
    goals::Array{Int} # The list of possible goals that this vehicle can have
    noise::Noise # Noise on position and velocity
end

# Creates a copy of the state  with the new specified vehicle state
update_veh_state(s::NoisyPedState, veh_state::VehicleState) = NoisyPedState(veh_state, s.noise)
update_veh_state(s::BlinkerState, veh_state::VehicleState) = BlinkerState(veh_state, s.blinker, s.goals, s.noise)

# Position and velocity functions
AutomotiveSimulator.posf(s::NoisyState) = posf(s.veh_state)
AutomotiveSimulator.posg(s::NoisyState) = posg(s.veh_state)
AutomotiveSimulator.vel(s::NoisyState)  = vel(s.veh_state)
AutomotiveSimulator.velf(s::NoisyState) = velf(s.veh_state)
AutomotiveSimulator.velg(s::NoisyState) = velg(s.veh_state)

# Get the blinker of a state or entity
blinker(s) = false # default behavior
blinker(s::BlinkerState) = s.blinker
blinker(veh::Entity) = blinker(veh.state)

# Get the goals of a state or entity
goals(s) = [laneid(s)] # default behavior
goals(s::BlinkerState) = s.goals
goals(veh::Entity) = goals(veh.state)

# Get the noise of a state or entity
noise(veh) = Noise() # default behavior
noise(s::NoisyState) = s.noise
noise(veh::Entity) = noise(veh.state)

# Makes a copy of the scene with the noise added to the vehicles in the state
function noisy_scene(scene::Scene, roadway::Roadway)
    noisy_scene = Scene(Entity)
    for (i,veh) in enumerate(scene)
        push!(noisy_scene, noisy_entity(veh, roadway))
    end
    noisy_scene
end

# Makes a copy of the entity but with the noise applied to the VehicleState
function noisy_entity(ent, roadway::Roadway)
    f = posf(ent)
    Δs, Δt = noise(ent).pos
    noisy_f = Frenet(roadway[f.roadind.tag], f.s + Δs, f.t + Δt, f.ϕ)
    noisy_g = posg(noisy_f, roadway)
    noisy_v = vel(ent) + noise(ent).vel
    noisy_vs = VehicleState(noisy_g, noisy_f, noisy_v)
    Entity(update_veh_state(ent.state, noisy_vs), ent.def, ent.id)
end

const PEDESTRIAN_DEF = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)

# Constructor for making a blinker vehicle
function BlinkerVehicle(;roadway::Roadway, lane::Int, s::Float64, v::Float64, id::Int, goals::Array{Int}, blinker::Bool, t::Float64 = 0., ϕ::Float64 = 0., noise::Noise = Noise())
    f = Frenet(roadway, lane, s, t, ϕ)
    bs = BlinkerState(VehicleState(f, roadway, v), blinker, goals, noise)
    Entity(bs, VehicleDef(), id)
end

# Constructor for making a noisy pedestrian
function NoisyPedestrian(;roadway::Roadway, lane::Int, s::Float64, v::Float64, id::Int, t::Float64=0., ϕ::Float64 = 0., noise::Noise = Noise())
    f = Frenet(roadway_lane(roadway, lane), s, t, ϕ)
    ps = NoisyPedState(VehicleState(f, roadway, v), noise)
    Entity(ps, PEDESTRIAN_DEF, id)
end

# Instructions for rendering the Blinker Vehicle
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, veh::Entity{BlinkerState, VehicleDef, Int64})
    reg_veh = Entity(veh.state.veh_state, veh.def, veh.id)
    add_renderable!(rendermodel, FancyCar(car=reg_veh))
    #TODO: Add noisy ghost
    li = laneid(veh)
    bo = BlinkerOverlay(on = blinker(veh), veh = reg_veh, right=Tint_signal_right[li])
    add_renderable!(rendermodel, bo)
    return rendermodel
end

# Instructions for rendering the noisy pedestrian
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, ped::Entity{NoisyPedState, VehicleDef, Int64})
    reg_ped = Entity(ped.state.veh_state, ped.def, ped.id)
    add_renderable!(rendermodel, FancyPedestrian(ped=reg_ped))
    #TODO: Add noisy ghost
    return rendermodel
end

## Definition of BlinkerVehicleControl action with helpers and propagate
abstract type Disturbance end

# Define a new action that conrols a blinker vehicle
@with_kw struct BlinkerVehicleControl <: Disturbance
    a::Float64 = 0 # Acceleration along lane (set by the model)
    da::Float64 = 0 # Acceleration disturbance along lane
    toggle_goal::Bool = false # Switches vehicle turn intention
    toggle_blinker::Bool = false # Toggles the turn signal
    noise::Noise = Noise() #
end

# Define a new action that controls a pedestrian
@with_kw struct PedestrianControl <: Disturbance
    a::VecE2 = VecE2(0., 0.)
    da::VecE2 = VecE2(0., 0.)
    noise::Noise = Noise()
end

# The function that propagates the PedestrianControl action
function AutomotiveSimulator.propagate(ped::Entity{NoisyPedState, D, I}, action::PedestrianControl, roadway::Roadway, Δt::Float64) where {D, I}
    starting_lane = laneid(ped)
    vs_entity = Entity(ped.state.veh_state, ped.def, ped.id)
    a_lat_lon = reverse(action.a + action.da)
    vs = propagate(vs_entity, LatLonAccel(a_lat_lon...), roadway, Δt)
    vs = VehicleState(vs.posG, vs.posF, clamp(vs.v, -3., 3.)) # Max pedestrian speed
    nps = NoisyPedState(set_lane(vs, laneid(ped), roadway), action.noise)
    @assert starting_lane == laneid(nps)
    nps
end

# The function that propogates the BlinkerVehicleControl action
function AutomotiveSimulator.propagate(veh::Entity{BlinkerState, D, I}, action::BlinkerVehicleControl, roadway::Roadway, Δt::Float64) where {D,I}
    # set the new goal
    vs = veh.state.veh_state
    if action.toggle_goal
        gs = goals(veh)
        curr_index = findfirst(gs .== laneid(veh))
        @assert !isnothing(curr_index)
        new_goal = gs[curr_index % length(gs) + 1]
        if can_have_goal(veh, new_goal, roadway)
            vs = set_lane(vs, new_goal, roadway)
        end
    end

    starting_lane = laneid(vs)

    # Update the kinematics of the vehicle (don't allow v < 0)
    vs_entity = Entity(vs, veh.def, veh.id)
    vs = propagate(vs_entity, LaneFollowingAccel(action.a + action.da), roadway, Δt)

    # Set the blinker state and return
    new_blink = action.toggle_blinker ? !blinker(veh) : blinker(veh)
    bs = BlinkerState(vs, new_blink, goals(veh), action.noise)
    @assert starting_lane == laneid(bs)
    bs
end

## Define a pedestrian control module

# Define the wrapper for the adversarial pedestrian
@with_kw mutable struct AdversarialPedestrian <: DriverModel{PedestrianControl}
    idm::IntelligentDriverModel = IntelligentDriverModel(v_des= 1.0)
    next_action::PedestrianControl = PedestrianControl()
    ignore_idm = false
end

# Sample an action from AdversarialPedestrian model
function Base.rand(rng::AbstractRNG, model::AdversarialPedestrian)
    na = model.next_action
    if !model.ignore_idm
        return PedestrianControl((model.idm.a, 0), na.da, na.noise)
    else
        return na
    end
end

# Observe function for AdversarialPedestrian model
function AutomotiveSimulator.observe!(model::AdversarialPedestrian, scene::Scene, roadway::Roadway, egoid::Int64)
    ego = get_by_id(scene, egoid)
    ego_v = vel(ego)
    track_longitudinal!(model.idm, ego_v, NaN, NaN)
    model
end

## Definition of Blindspot
struct Blindspot
    θ::Float64 # Offset of the blindspot from vehicle direction
    ϕ::Float64 # Width of the blindspot
end

function in_blindspot(egopos::VecSE2{Float64}, blindspot::Blindspot, otherpos::VecSE2{Float64})
    Δy = otherpos.y - egopos.y
    Δx = otherpos.x - egopos.x
    ψ = atan(Δy, Δx)
    θ = blindspot.θ + egopos.θ
    (ψ > θ - blindspot.ϕ / 2.) && (ψ < θ  + blindspot.ϕ / 2.)
end

struct RenderableBlindspot
    pos::VecSE2
    blindspot::Blindspot
    length::Float64
    color::Colorant
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, rb::RenderableBlindspot)
    θ = rb.blindspot.θ + rb.pos.θ
    low = θ - rb.blindspot.ϕ / 2.
    high = θ + rb.blindspot.ϕ / 2.
    pt1 = VecE2{Float64}(rb.length*cos(low), rb.length*sin(low))
    pt2 = VecE2{Float64}(rb.length*cos(high), rb.length*sin(high))
    p = VecE2(rb.pos)
    pts = [p, p+pt1, p+pt2, p]

    add_instruction!(
        rendermodel, AutomotiveVisualization.render_closed_line,
        (pts,  rb.color, 0.1, rb.color),
        coordinate_system=:scene
    )
    return rendermodel
end

## Definition of the T-Intersection DriverModel with observe! and rand()

# Define a driving model for a T-intersection IDM model
@with_kw mutable struct TIDM <: DriverModel{BlinkerVehicleControl}
    idm::IntelligentDriverModel = IntelligentDriverModel(v_des = 15.) # underlying idm
    noisy_observations::Bool = false # Whether or not this model gets noisy observations
    ttc_threshold = 7 # threshold through intersection
    next_action::BlinkerVehicleControl = BlinkerVehicleControl() # The next action that the model will do (for controllable vehicles)

    # Describes the intersection and rules of the road
    yields_way::Dict{Int64, Array{Int64}} = Dict() # Lane priorities
    intersection_enter_loc::Dict{Int64, VecSE2} = Dict() # Entry location of intersection
    intersection_exit_loc::Dict{Int64, VecSE2} = Dict()  # Exit location of intersection
    goals::Dict{Int64, Array{Int64}} = Dict() # Possible goals of each lane
    should_blink::Dict{Int64, Bool} = Dict()  # Wether or not the blinker should be on

    blindspot::Union{Blindspot, Nothing} = nothing # Blindspot of the vehicle
end

# Sample an action from TIDM model
function Base.rand(rng::AbstractRNG, model::TIDM)
    na = model.next_action
    BlinkerVehicleControl(model.idm.a, na.da, na.toggle_goal, na.toggle_blinker, na.noise)
 end

# Observe function for TIDM
function AutomotiveSimulator.observe!(model::TIDM, input_scene::Scene, roadway::Roadway, egoid::Int64)
    # If this model is susceptible to noisy observations, adjust all the agents by noise
    scene = model.noisy_observations ? noisy_scene(input_scene, roadway) : input_scene

    # Get the ego and the ego lane
    ego = get_by_id(scene, egoid)
    ego_v = vel(ego)
    li = laneid(ego)

    # Compute ego headway to intersection and time to cross
    intrsxn_Δs = distance_to_point(ego, roadway, model.intersection_enter_loc[li])
    intrsxn_exit_Δs = distance_to_point(ego, roadway, model.intersection_exit_loc[li])
    time_to_cross = time_to_cross_distance_const_acc(ego, model.idm, distance_to_point(ego, roadway, model.intersection_exit_loc[li]))

    # Get headway to the forward car
    fore = find_neighbor(scene, roadway, ego, targetpoint_ego = VehicleTargetPointFront(), targetpoint_neighbor = VehicleTargetPointRear())
    fore_v, fore_Δs = isnothing(fore.ind) ? (NaN, Inf) : (vel(scene[fore.ind]), fore.Δs)

    # Check to see if ego car has right of way
    has_right_of_way = true
    lanes_to_yield_to = model.yields_way[li]
    vehicles_to_yield_to = []
    for (i,veh) in enumerate(scene)
        # if the vehicle is the ego vehicle then move on
        veh.id == egoid && continue

        # Check if the other entity is in the blind spot. If so, move one
        !isnothing(model.blindspot) && in_blindspot(posg(ego), model.blindspot, posg(veh)) && continue

        # Check to see if the vehicle is in the ego's lane without the same laneid
        if laneid(ego) != laneid(veh)
            Δs_inlane = compute_inlane_headway(ego, veh, roadway)
            if Δs_inlane < fore_Δs
                fore_Δs = Δs_inlane
                fore_v = vel(veh)
            end
        end

        # If the vehicle is in a lane that the ego should yield to, store it
        if lane_belief(veh, model, roadway) in lanes_to_yield_to
            has_right_of_way = false
            push!(vehicles_to_yield_to, veh)
        end
    end

    # TODO: Reasoning to take place if you enter an intersection
    # if i_headway <-r
        # has_right_of_way = true
    # end

    next_idm = track_longitudinal!(model.idm, ego_v, fore_v, fore_Δs)
    # If the vehicle does not have right of way then stop before the intersection
    if !has_right_of_way && intrsxn_exit_Δs > 0
        Nv = length(vehicles_to_yield_to)
        in_intersection = Vector{Bool}(undef, Nv)
        exit_time = Vector{Float64}(undef, Nv)
        enter_time = Vector{Float64}(undef, Nv)
        Δs_in_lane = Vector{Float64}(undef, Nv)
        for (i, veh) in zip(1:Nv, vehicles_to_yield_to)
            _exit = distance_to_point(veh, roadway, model.intersection_exit_loc[laneid(veh)])
            _enter = distance_to_point(veh, roadway, model.intersection_enter_loc[laneid(veh)])
            in_intersection[i] = (_enter < 2 && _exit > 0)
            exit_time[i] = time_to_cross_distance_const_vel(veh, _exit)
            enter_time[i] = time_to_cross_distance_const_vel(veh, _enter)
            Δs_in_lane[i] = compute_inlane_headway(ego, veh, roadway)
        end

        # The intersection is clear of car i if, it exited the intersection in the past, or
        # it will enter the intersection after you have crossed it, or
        # it will have exited a while before you crossed
        intersection_clear = (exit_time .<= 0) .| (enter_time .> time_to_cross) .| (exit_time .+ model.ttc_threshold .< time_to_cross)
        intersection_clear = intersection_clear .& (.!in_intersection)
        intersection_clear = intersection_clear .& (Δs_in_lane .> intrsxn_exit_Δs)
        if !all(intersection_clear)
            # yield to oncoming traffic
            minΔs_to_yield = minimum(Δs_in_lane[.!intersection_clear]) # headways of cars you care about
            # If there isn't a leading car, or if it is past the intersection, use intersection point, otherwise use car
            if isinf(fore_Δs) || minΔs_to_yield < fore_Δs || (intrsxn_Δs > 0 && intrsxn_Δs < fore_Δs)
                next_idm = track_longitudinal!(model.idm, ego_v, 0., min(minΔs_to_yield, intrsxn_Δs))
            end
        end

    end
    model.idm = next_idm
    model
end

# function find_nearest_car_forward(ego, scene, roadway)
#     for (i,veh) in enumerate(scene)
#         if veh.id != egoid &&
# end

function compute_inlane_headway(ego, veh, roadway)
    elane = laneid(ego)
    v = Entity(set_lane(veh.state.veh_state, elane, roadway), veh.def, veh.id)
    f = posf(v)
    Δs = f.s - v.def.length/2. - posf(ego).s - targetpoint_delta(VehicleTargetPointFront(), ego)
    tbound = get_lane(roadway, ego.state.veh_state).width / 2.
    (abs(f.t) <= tbound && Δs > 0) ? Δs : Inf
end

# Get the belief of the lane of the specificed vehicle
function lane_belief(veh::Entity, model::TIDM, roadway::Roadway)
    possible_lanes = model.goals[laneid(veh)]
    length(possible_lanes) == 1 && return possible_lanes[1]
    possible_lanes = possible_lanes[[can_have_goal(veh, l, roadway) for l in possible_lanes]]
    length(possible_lanes) == 1 && return possible_lanes[1]

    blinker_match = findfirst([blinker(veh) == model.should_blink[l] for l in possible_lanes])
    if !isnothing(blinker_match)
        return possible_lanes[blinker_match]
    else
        return model.goals[laneid(veh)][1]
    end
end

## Definition of helper functions for managing vehicles on the road
AutomotiveSimulator.get_lane(roadway::Roadway, state::BlinkerState) = roadway[posf(state).roadind.tag]
AutomotiveSimulator.get_lane(roadway::Roadway, state::NoisyPedState) = roadway[posf(state).roadind.tag]

# Get the laneid of a vehicle
laneid(veh) = posf(veh).roadind.tag.segment

# Returns a vehicle with the same state, but projected onto the desired lane
function set_lane(state::VehicleState, laneid::Int, roadway::Roadway)
    desired_lane = roadway_lane(roadway, laneid)
    posF = Frenet(posg(state), desired_lane, roadway)
    VehicleState(posF, roadway, vel(state))
end

# Returns a new vehicle that is a copy of the provided vehicle except with the specified id
set_veh_id(veh::Entity, id::Int) = Entity(veh.state, veh.def, id)

# Check whether a car is allowed to switch goals (requires posF.t to be small)
function can_have_goal(veh::Entity, goal::Int, roadway::Roadway; ϵ = 1e-6)
    desired_lane = roadway[goal].lanes[1]
    posF = Frenet(posg(veh), desired_lane, roadway)
    return abs(posF.t) <= ϵ
end

# Computes the time it takes to cover a given distance, assuming the maximum acceleration of the provided idm
function time_to_cross_distance_const_acc(veh::Entity, idm::IntelligentDriverModel, ds::Float64)
    v = vel(veh)
    d = v^2 + 2*idm.a_max*ds
    d < 0 && return 0 # We have already passed the point we are trying to get to
    vf = min(idm.v_des, sqrt(d))
    2*ds/(vf + v)
end

# Compute the time it takes to travel a given distance at the current vehicle velocity
time_to_cross_distance_const_vel(veh::Entity, ds::Float64) = ds / vel(veh)

# Computes the distance between the vehicle and the provided point along the roadway
# If the pt is not on the roadway then it is the distance to the projection onto the road
function distance_to_point(veh::Entity, roadway::Roadway, pt::VecSE2)
    s0 = posf(veh).s
    lane = get_lane(roadway, veh)
    s_end = Frenet(pt, lane, roadway).s
    s_end - s0
end

# computes the closest distance between the ego vehicle and another vehicle on the road
function min_dist(scene, egoid)
    minval = 10000
    (length(scene) == 0 || !(egoid in scene)) && return minval
    ego_pos = posg(get_by_id(scene, egoid).state)
    for (i,veh) in enumerate(scene)
        if veh.id != egoid
            pos = posg(veh)
            dist = sqrt((pos.x - ego_pos.x)^2 + (pos.y - ego_pos.y)^2)
            if dist < minval
                minval = dist
            end
        end
    end
    minval
end

# Define a callback to remove cars that have reached the end of the lane
struct CleanSceneCallback
    egoid
end

# Check if the ego vehicle is in collision with another car
function ego_collides(egoid, scene)
    if !(egoid in scene)
        return false
    end

    ego = get_by_id(scene, egoid)
    for (i,veh) in enumerate(scene)
        if egoid != veh.id && collision_checker(ego, veh) && vel(ego) > 0.1
            return true
        end
    end
    false
end

# Check if any two vehicles collide
function any_collides(scene)
    for (i,veh) in enumerate(scene)
        for (j,veh2) in enumerate(scene)
            j <= i && continue
            collision_checker(veh, veh2) && return true
        end
    end
    false
end

# Check if a given car has reached the end of the roadway
function end_of_road(veh, roadway, set_end)
    s = posf(veh).s

    lane = get_lane(roadway, veh)
    s_end = lane.curve[end].s
    return s >= min(s_end, set_end)
end

# Removes cars that have reached the end of the lane
function AutomotiveSimulator.run_callback(c::CleanSceneCallback, scenes::Array{Scene}, actions::Nothing, roadway::R, models::Dict{I,M}, tick::Int) where {F,I,R,M<:DriverModel}
    for (i,veh) in enumerate(scenes[tick])
        if end_of_road(veh, roadway)
           deleteat!(scenes[tick], findfirst(veh.id, scenes[tick]))
       end
    end
    return ego_collides(scenes[tick], c.egoid)
end

## This stuff is used when scene has different types
Base.in(id, scene::Scene) = !isnothing(findfirst(id, scene))

function Base.findfirst(id, scene::Scene)
    for entity_index in 1 : scene.n
        scene.entities[entity_index].id == id && return entity_index
    end
    nothing
end

function AutomotiveSimulator.get_by_id(scene::Scene, id)
    entity_index = findfirst(id, scene)
    isnothing(entity_index) && throw(BoundsError(scene, [id]))
    scene[entity_index]
end

## Neural Network driving model
@with_kw mutable struct PolicyModel <: DriverModel{LaneFollowingAccel}
    policy
    state = nothing
end

# Sample an action from TIDM model
function Base.rand(rng::AbstractRNG, model::PolicyModel)
    action(model.policy, model.state)
end

# Observe function for TIDM
function AutomotiveSimulator.observe!(model::PolicyModel, input_scene::Scene, roadway::Roadway, egoid::Int64)
    model.state = input_scene
end


text:
# p22.jl - 5th eigenvector of Airy equation u_xx = lambda*x*u

clf();
for N = 12:12:48
    (D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
    (lam,V) = eig(D2,diagm(x[2:N]));      # generalized ev problem
    ii = find(lam.>0);
    V = V[:,ii]; lam = lam[ii];
    ii = sortperm(lam)[5]; lambda = lam[ii];
    v = [0;V[:,ii];0]; v = v/v[Int(N/2+1)]*airyai(0);
    xx = -1:.01:1; vv = polyval(polyfit(x,v),xx);
    subplot(2,2,N/12); plot(xx,vv); grid(true)
    title("N = $N     eig = $(signif(lambda,13))");
end


text:
using Test
using RCall
using EmpiricalBayes
using StatsBase
using Distributions

R"load('./notebooks/datasets/prostz.Rdata')"

@rget prostz
R"library(deconvolveR)"

R"tau <- seq(-3.6,3.6, by = 0.2)"
@rget tau
R"result <- deconv(tau = tau, X = prostz, deltaAt = 0, family = 'Normal', pDegree = 5)"
R"result_mat <- result$stats"
@rget result_mat

#result_mat[18:20, 2:3]
#plot(tau, result_mat[:,2])
#tau_nz = copy(tau)
#deleteat!(tau_nz, 19)
#gs_nz = result_mat[:,2]
#deleteat!(gs_nz, 19)
#plot(tau_nz, gs_nz)

R"result_cov <- result$cov.g"
@rget result_cov

g_prob = result_mat[:,2]
g_bias = result_mat[:,6]

prob_greater_1 = sum(g*(abs(t) >1) for (t,g) in zip(tau,g_prob) )
prob_greater_1 = sum(g*(abs(t) >1) for (t,g) in zip(tau,g_prob) )

greater_1_linearization = one(Float64).*(abs.(tau) .>1)
greater_1_linearization'*(g_prob .+ g_bias)
greater_1_linearization'*g_prob

prob_greater_1_ci = 1.96*sqrt(greater_1_linearization'*result_cov*greater_1_linearization)

greater_eq_2_linearization = one(Float64).*(abs.(tau) .>=2)
greater_eq_2_linearization'*g_prob
prob_greater_2_ci = 1.96*sqrt(greater_eq_2_linearization'*result_cov*greater_eq_2_linearization)

# Why is this factor 2 off efron's result?

R"result_P <- result$P"
@rget result_P

se_1 = diag(result_cov).^0.5
se_2 = result_mat[:, 3]

@test se_1 == se_2

R"class(result) <- 'deconvolveR'"
R"class(result)"

a = R"result"

brad_jl = fit(BradDeconvolveR, prostz; c0=1.0)

@test greater_eq_2_linearization'*g_prob ≈ estimate(brad_jl, PriorTailProbability(2.0);debias=false)

(l1,r1) = EmpiricalBayes.confint(brad_jl, PriorTailProbability(2.0), 0.95; debias=false)
r1-l1

estimate(brad_jl, PosteriorTarget(PosteriorMeanNumerator(2.0));debias=false)

confint(brad_jl, PosteriorTarget(PosteriorMeanNumerator(2.0)), 0.95, ;debias=false)

## With LFSR Numerator

## First a case in which result is really what we do not want
## which happens when there is a large delta at 0.

estimate(brad_jl, PosteriorTarget(LFSRNumerator(0.0));debias=true)
confint(brad_jl, PosteriorTarget(LFSRNumerator(0.0)), 0.95, ;debias=false)

lfsr_riesz = riesz_representer.(LFSRNumerator(2.0), tau)
lfsr_riesz_at_0 = copy(lfsr_riesz)
lfsr_riesz_at_0[19] /=2
lfsr_riesz'*g_prob
lfsr_riesz_at_0'*g_prob

#
brad_jl_no0 = fit(BradDeconvolveR, prostz; c0=1.0,
                     deltaAt = :nothing)

lfsr_riesz_at_0'*brad_jl_no0.g_prior
lfsr_riesz'*brad_jl_no0.g_prior

# Difference still not negligible, let us use a much larger grid.

tau2 = collect(-3.6:0.002:3.6)
brad_jl_tau2 = fit(BradDeconvolveR, prostz; c0=0.0,
                     deltaAt = :nothing,  prior_grid=tau2)

lfsr_riesz = riesz_representer.(LFSRNumerator(2.0), tau2)
lfsr_riesz_at_0 = copy(lfsr_riesz)
lfsr_riesz_at_0[181] /=2

lfsr_riesz_at_0'*brad_jl_tau2.g_prior
lfsr_riesz'*brad_jl_tau2.g_prior

confint(brad_jl_tau2, LFSRNumerator(2.0), 0.95, ;debias=false)
# Discretization error here is larger than confindence interval length. Ugh.

# Let us..

tau2
# OK need to rework this from scratch more or less..
# Idea: Make grid finer and properly integrate from a bit left to right.

# TODO: 1) Implement the below
# 2) Check if method gives correct coverage in Normal/Normal example
# 3) Check how method behaves in our other two examples which we want to simulate from..
function riesz_representer(target::LinearInferenceTarget, b::BradDeconvolveR)

end

# Afou ftasw kai isws koimh8w ligo ->
# prwta na dw ola ta experiment settings pou 8elw na tre3w
# Meta na dw pws leitourgei tou Brad, toulaxiston kapoies fores 8a 8elw swsto coverage
# Isws kai ena grhgoro peirama na dw ti coverage petyxainoume...

# Maybe even demonstrate everything on 1 simulation example....


text:
function dimension(mesh::UniformMesh{dim}) where dim
    return dim
end

function number_of_elements(mesh)
    return mesh.total_number_of_elements
end

function linear_to_cartesian_index(mesh::UniformMesh{2},idx::Z) where {Z<:Integer}

    1 <= idx <= mesh.total_number_of_elements || throw(BoundsError(mesh,idx))
    elx = ceil(Int,idx/mesh.nelements[2])
    ely = idx - (elx-1)*mesh.nelements[2]
    return elx,ely

end

function element(mesh::UniformMesh{2},idx)

    elx,ely = linear_to_cartesian_index(mesh,idx)
    xL = mesh.x0 + [(elx-1)*mesh.element_size[1],(ely-1)*mesh.element_size[2]]
    xR = mesh.x0 + [(elx)*mesh.element_size[1],(ely)*mesh.element_size[2]]
    return xL,xR
end

function element(mesh::UniformMesh{1},idx)
    xL = mesh.x0 + (idx-1)*mesh.element_size
    xR = xL + mesh.element_size
    return xL,xR
end

function neighbors(mesh::UniformMesh{2},idx)

    elx,ely = linear_to_cartesian_index(mesh,idx)
    bdryflag = mesh.bdryflag

    bottom = ely == 1 ? bdryflag : idx-1
    right = elx == mesh.nelements[1] ? bdryflag : idx+mesh.nelements[2]
    top = ely == mesh.nelements[2] ? bdryflag : idx+1
    left = elx == 1 ? bdryflag : idx-mesh.nelements[2]

    return [bottom,right,top,left]

end

function neighbors(mesh::UniformMesh{1},idx)
    nelements = number_of_elements(mesh)
    @assert 1 <= idx <= nelements
    left = idx == 1 ? 0 : idx-1
    right = idx == nelements ? 0 : idx+1
    return [left,right]
end

function faces_per_cell(mesh::UniformMesh{2})
    return 4
end

function faces_per_cell(mesh::UniformMesh{1})
    return 2
end


text:
module SpaceMethods

export AbstractSpaceMethod, AbstractDGMethod, odefun!

"""
    AbstractSpaceMethod

Supertype for spatial discretizations
"""
abstract type AbstractSpaceMethod end

"""
    AbstractDGMethod <: AbstractSpaceMethod

Supertype for discontinuous Galerkin spatial discretizations
"""
abstract type AbstractDGMethod <: AbstractSpaceMethod end

"""
    odefun!(disc::AbstractSpaceMethod, dQ, Q, t; increment)

Evaluates the right-hand side of the spatial discretization defined by `disc` at
time `t` with state `Q`.
The result is either added into `dQ` if `increment` is true or stored in `dQ` if it is false.
Namely, the semi-discretization is of the form
``
  \\dot{Q} = F(Q, t)
``
and after the call `dQ += F(Q, t)` if `increment == true`
or `dQ = F(Q, t)` if `increment == false`

!!! note

    There is no generic implementation of this function. This must be
    implemented for each subtype of `AbstractSpaceMethod`
"""
function odefun! end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
nz=43; nx=35
ny=randn(nz, nx);

my1=GeoPhyInv.get_rhovxI(ny)
my2=GeoPhyInv.get_rhovzI(ny)

myp1 = randn(size(my1));
myp2 = randn(size(my2));

# spray
nyp=zeros(nz, nx);

GeoPhyInv.grad_modrr_sprayrr!(nyp, myp1, myp2)
GeoPhyInv.grad_modrr_sprayrrvx!(nyp, myp1)
GeoPhyInv.grad_modrr_sprayrrvz!(nyp, myp2)

# dot product test
@test LinearAlgebra.dot(vcat(my1,my2), vcat(myp1,myp2)) ≈ LinearAlgebra.dot(ny, nyp)


text:
struct Closure{T}
    params::T
end
(s::Closure)(x::Vector) = s.apply(x, s.params)

mutable struct StepClosure{S,F,G,H} <: AbstractStep
    selections::S
    params::Union{Nothing, Dict{Symbol,Closure}}
    get_params::F
    apply::G
    broadcast::Bool
    trained::Bool
    skip::Bool
    prehook::H
end

function transform!(s::StepClosure, df)
    if s.broadcast
        for (k,f) in s.params
            df[k] = s.apply.(df[k], Ref(f.params))
        end
    else
        for (k,f) in s.params
            df[k] = s.apply(df[k], f.params)
        end
    end
    df
end
function transform(s::StepClosure, df)
    transform!(s, copy(df))
end
function fit!(s::StepClosure, df)
    s.params = Dict(col => Closure(s.get_params(s.prehook(df[col]))) for col in getselectionkeys(df, s.selections))
    s.trained = true
end
function step_closure!(r::Recipe, s...; get_params=nothing, apply=nothing, broadcast::Bool=true, skip::Bool=false, prehook=identity)
    get_params === nothing || apply === nothing && throw(error("`get_params` or `apply` kwargs not set for `step_closure!`"))
    push!(r.steps,
          StepClosure([s...],
                       nothing,
                       get_params,
                       apply,
                       broadcast,
                       false,
                       skip,
                       prehook))
end


text:
const ID_LENGTH = 32
const ENDIANNESS = 16911887
const SDF_VERSION = 1
const SDF_REVISION = 4
const SDF_MAGIC = "SDF1"

const BLOCKTYPE_SCRUBBED = Int32(-1)
const BLOCKTYPE_NULL = Int32(0)
const BLOCKTYPE_PLAIN_MESH = Int32(1)
const BLOCKTYPE_POINT_MESH = Int32(2)
const BLOCKTYPE_PLAIN_VARIABLE = Int32(3)
const BLOCKTYPE_POINT_VARIABLE = Int32(4)
const BLOCKTYPE_CONSTANT = Int32(5)
const BLOCKTYPE_ARRAY = Int32(6)
const BLOCKTYPE_RUN_INFO = Int32(7)
const BLOCKTYPE_SOURCE = Int32(8)
const BLOCKTYPE_STITCHED_TENSOR = Int32(9)
const BLOCKTYPE_STITCHED_MATERIAL = Int32(10)
const BLOCKTYPE_STITCHED_MATVAR = Int32(11)
const BLOCKTYPE_STITCHED_SPECIES = Int32(12)
const BLOCKTYPE_SPECIES = Int32(13)
const BLOCKTYPE_PLAIN_DERIVED = Int32(14)
const BLOCKTYPE_POINT_DERIVED = Int32(15)
const BLOCKTYPE_CONTIGUOUS_TENSOR = Int32(16)
const BLOCKTYPE_CONTIGUOUS_MATERIAL = Int32(17)
const BLOCKTYPE_CONTIGUOUS_MATVAR = Int32(18)
const BLOCKTYPE_CONTIGUOUS_SPECIES = Int32(19)
const BLOCKTYPE_CPU_SPLIT = Int32(20)
const BLOCKTYPE_STITCHED_OBSTACLE_GROUP = Int32(21)
const BLOCKTYPE_UNSTRUCTURED_MESH = Int32(22)
const BLOCKTYPE_STITCHED = Int32(23)
const BLOCKTYPE_CONTIGUOUS = Int32(24)
const BLOCKTYPE_LAGRANGIAN_MESH = Int32(25)
const BLOCKTYPE_STATION = Int32(26)
const BLOCKTYPE_STATION_DERIVED = Int32(27)
const BLOCKTYPE_DATABLOCK = Int32(28)
const BLOCKTYPE_NAMEVALUE = Int32(29)

const DATATYPE_NULL = Int32(0)
const DATATYPE_INTEGER4 = Int32(1)
const DATATYPE_INTEGER8 = Int32(2)
const DATATYPE_REAL4 = Int32(3)
const DATATYPE_REAL8 = Int32(4)
const DATATYPE_REAL16 = Int32(5)
const DATATYPE_CHARACTER = Int32(6)
const DATATYPE_LOGICAL = Int32(7)
const DATATYPE_OTHER = Int32(8)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis: vandermonde_matrix
using Test

polynomial = :(1 + u + v + u*v)
coordinates = ((-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0))
V = vandermonde_matrix(polynomial, coordinates)
@debug "create Vandermonde matrix" polynomial coordinates V
V_expected = [
    1.0 -1.0 -1.0  1.0
    1.0  1.0 -1.0 -1.0
    1.0  1.0  1.0  1.0
    1.0 -1.0  1.0 -1.0]
@test isapprox(V, V_expected)

polynomial = :(1 + u + v)
coordinates = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
V = vandermonde_matrix(polynomial, coordinates)
V_expected = [
    1.0 0.0 0.0
    1.0 1.0 0.0
    1.0 0.0 1.0]
@test isapprox(V, V_expected)