text:
abstract type MaskedArray{T, N} <: AbstractArray{T, N} end

@inline isactive(x::MaskedArray, i::Int) = (@_propagate_inbounds_meta; getmask(x)[i])
@inline checkactive(::Type{Bool}, x::MaskedArray, i::Int) = (@_propagate_inbounds_meta; isactive(x, i)) # checkbounds as well
@inline checkactive(x::MaskedArray, i::Int) = (@_propagate_inbounds_meta; checkactive(Bool, x, i) || throw(UndefRefError()))

@inline function Base.getindex(x::MaskedArray, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        checkactive(x, i)
        x.data[i]
    end
end

@inline function Base.setindex!(x::MaskedArray, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        x.data[i] = v
        x.mask[i] = true
    end
    x
end

@inline function Base.setindex!(x::MaskedArray, ::Nothing, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.mask[i] = false
    x
end

@inline function unsafe_getindex(x::MaskedArray, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[i]
end

@inline function unsafe_setindex!(x::MaskedArray, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[i] = v
    v
end

getmask(x::MaskedArray) = x.mask

# `f` should be `f(mask)`
function findentry(f, x::MaskedArray)
    i = f(vec(x.mask)) # `vec` is not needed for BitArray
    i === nothing ? nothing : @inbounds x[i]
end

struct MaskedDenseArray{T, N} <: MaskedArray{T, N}
    data::Array{T, N}
    mask::Array{Bool, N} # faster than BitArray?
end

@inline function MaskedDenseArray{T}(::UndefInitializer, dims::NTuple{N, Int}) where {T, N}
    data = Array{T}(undef, dims)
    mask = fill!(similar(data, Bool), false)
    MaskedDenseArray(data, mask)
end
@inline MaskedDenseArray{T}(u::UndefInitializer, dims::Vararg{Int, N}) where {T, N} =
    MaskedDenseArray{T}(u, dims)
@inline function MaskedDenseArray(data::Array)
    mask = fill!(similar(data, Bool), false)
    MaskedDenseArray(data, mask)
end

Base.IndexStyle(::Type{<: MaskedDenseArray}) = IndexLinear()
Base.size(x::MaskedDenseArray) = size(x.data)

@inline Base.isassigned(x::MaskedDenseArray, i::Int) = isassigned(x.data, i)

struct MaskedStructArray{T, N, Ttuple} <: MaskedArray{T, N}
    data::StructArray{T, N, Ttuple, Int}
    mask::Array{Bool, N} # faster than BitArray?
end

@inline function MaskedStructArray{T}(::UndefInitializer, dims::NTuple{N, Int}) where {T, N}
    data = StructArray{T}(undef, dims)
    mask = fill!(similar(data, Bool), false)
    MaskedStructArray(data, mask)
end
@inline MaskedStructArray{T}(u::UndefInitializer, dims::Vararg{Int, N}) where {T, N} =
    MaskedStructArray{T}(u, dims)
@inline function MaskedStructArray(data::StructArray)
    mask = fill!(similar(data, Bool), false)
    MaskedStructArray(data, mask)
end

Base.propertynames(x::MaskedStructArray) = (:data, :mask, propertynames(x.data)...)
@inline function Base.getproperty(x::MaskedStructArray, name::Symbol)
    name == :data && return getfield(x, :data)
    name == :mask && return getfield(x, :mask)
    MaskedDenseArray(getproperty(getfield(x, :data), name), getfield(x, :mask))
end

Base.IndexStyle(::Type{<: MaskedStructArray}) = IndexLinear()
Base.size(x::MaskedStructArray) = size(x.data)

@inline Base.isassigned(x::MaskedStructArray, i::Int) = isassigned(x.data, i)

# https://discourse.julialang.org/t/poor-time-performance-on-dict/9656/14
struct FastHashInt; i::Int; end
Base.:(==)(x::FastHashInt, y::FastHashInt) = x.i == y.i
Base.hash(x::FastHashInt, h::UInt) = xor(UInt(x.i), h)

struct MaskedHashArray{T, N} <: MaskedArray{T, N}
    data::Dict{FastHashInt, T}
    mask::Array{Bool, N}
    dims::NTuple{N, Int}
end

@inline function MaskedHashArray{T}(::UndefInitializer, dims::NTuple{N, Int}) where {T, N}
    data = Dict{FastHashInt, T}()
    mask = fill!(Array{Bool}(undef, dims), false)
    MaskedHashArray(data, mask, dims)
end
@inline MaskedHashArray{T}(u::UndefInitializer, dims::Vararg{Int, N}) where {T, N} =
    MaskedHashArray{T}(u, dims)

Base.IndexStyle(::Type{<: MaskedHashArray}) = IndexLinear()
Base.size(x::MaskedHashArray) = x.dims

# Base.keys(x::MaskedHashArray) = keys(x.data) # disable because of using FastHashInt
Base.values(x::MaskedHashArray) = values(x.data)
Base.haskey(x::MaskedHashArray, i) = haskey(x.data, FastHashInt(i))
Base.delete!(x::MaskedHashArray, i) = delete!(x.data, FastHashInt(i))

@inline function Base.get(x::MaskedHashArray, i::Int, default)
    get(x.data, FastHashInt(i), default)
end

@inline function Base.getindex(x::MaskedHashArray, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        checkactive(x, i)
        x.data[FastHashInt(i)]
    end
end

@inline function Base.setindex!(x::MaskedHashArray, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        x.data[FastHashInt(i)] = v
        x.mask[i] = true
    end
    x
end
@inline function Base.setindex!(x::MaskedHashArray, ::Nothing, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.mask[i] = false
    x
end

@inline function unsafe_getindex(x::MaskedHashArray, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[FastHashInt(i)]
end

@inline function unsafe_setindex!(x::MaskedHashArray, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[FastHashInt(i)] = v
    v
end


text:
using Test
using LinearAlgebra
using Nystrom
using StaticArrays
using Random
Random.seed!(1)

@testset "DIM" begin
    # construct interior solution
    rtol=5e-2
    @testset "Greens identity (interior) 2d" begin
        Geometry.clear_entities!()
        Ω   = ParametricSurfaces.Circle() |> Geometry.Domain
        Γ   = boundary(Ω)
        M   = ParametricSurfaces.meshgen(Γ,(5,))
        mesh = NystromMesh(view(M,Γ),order=3)
        xout = SVector(3,3)
        ops = (
            Laplace(;dim=2),
            Helmholtz(;dim=2,k=1.2),
            Elastostatic(;dim=2,μ=2,λ=3)
        )
        for pde in ops
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(qnode,xout)*c
            dudn = (qnode) -> AdjointDoubleLayerKernel(pde)(qnode,xout)*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)
            γ₁u_norm = norm(norm.(γ₁u,Inf),Inf)
            # single and double layer
            S     = SingleLayerOperator(pde,mesh)
            Smat  = S |> Matrix
            D     = DoubleLayerOperator(pde,mesh)
            Dmat = D |> Matrix
            e0   = Nystrom.error_interior_green_identity(Smat,Dmat,γ₀u,γ₁u)/γ₀u_norm
            Sdim  = Nystrom.assemble_dim(S)
            Ddim  = Nystrom.assemble_dim(D)
            e1    = Nystrom.error_interior_green_identity(Sdim,Ddim,γ₀u,γ₁u)/γ₀u_norm
            @testset "Single/double layer $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
            # adjoint double-layer and hypersingular
            K     = AdjointDoubleLayerOperator(pde,mesh)
            Kmat     = K |> Matrix
            H     = HyperSingularOperator(pde,mesh)
            Hmat = H |> Matrix
            e0    = Nystrom.error_interior_derivative_green_identity(Kmat,Hmat,γ₀u,γ₁u)/γ₁u_norm
            Kdim  = Nystrom.assemble_dim(K)
            Hdim  = Nystrom.assemble_dim(H)
            e1   = Nystrom.error_interior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)/γ₁u_norm
            @testset "Adjoint double-layer/hypersingular $(string(pde))" begin
                @test norm(e0,Inf) > rtol
                @test norm(e1,Inf) < rtol
            end
        end
    end

    @testset "Greens identity (interior) 3d" begin
        Geometry.clear_entities!()
        Ω   = ParametricSurfaces.Sphere(;radius=1) |> Geometry.Domain
        Γ   = boundary(Ω)
        M   = ParametricSurfaces.meshgen(Γ,(2,2))
        mesh = NystromMesh(view(M,Γ),order=3)
        xout = SVector(3,3,3)
        ops = (
            Laplace(;dim=3),
            Helmholtz(;dim=3,k=1.2),
            Elastostatic(;dim=3,μ=2,λ=3),
            Maxwell(;k=2.)
        )
        for pde in ops
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xout,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xout,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)
            γ₁u_norm = norm(norm.(γ₁u,Inf),Inf)
            # single and double layer
            S     = SingleLayerOperator(pde,mesh)
            Smat  = S |> Matrix
            D     = DoubleLayerOperator(pde,mesh)
            Dmat = D |> Matrix
            e0   = Nystrom.error_interior_green_identity(Smat,Dmat,γ₀u,γ₁u)/γ₀u_norm
            Sdim  = Nystrom.assemble_dim(S)
            Ddim  = Nystrom.assemble_dim(D)
            e1    = Nystrom.error_interior_green_identity(Sdim,Ddim,γ₀u,γ₁u)/γ₀u_norm
            @testset "Single/double layer $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
            # adjoint double-layer and hypersingular
            pde isa Maxwell && continue
            K     = AdjointDoubleLayerOperator(pde,mesh)
            Kmat     = K |> Matrix
            H     = HyperSingularOperator(pde,mesh)
            Hmat = H |> Matrix
            e0   = Nystrom.error_interior_derivative_green_identity(Kmat,Hmat,γ₀u,γ₁u)/γ₁u_norm
            Kdim  = Nystrom.assemble_dim(K)
            Hdim  = Nystrom.assemble_dim(H)
            e1   = Nystrom.error_interior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)/γ₁u_norm
            @testset "Adjoint double-layer/hypersingular $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
        end
    end

    @testset "Greens identity (exterior) 2d" begin
        Geometry.clear_entities!()
        Ω   = ParametricSurfaces.Circle() |> Geometry.Domain
        Γ   = boundary(Ω)
        M   = ParametricSurfaces.meshgen(Γ,(7,))
        mesh = NystromMesh(view(M,Γ),order=3)
        xin = SVector(0.1,0.2)
        ops = (
            Laplace(;dim=2),
            Helmholtz(;dim=2,k=1.2),
            Elastostatic(;dim=2,μ=2,λ=3)
        )
        for pde in ops
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xin,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xin,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)
            γ₁u_norm = norm(norm.(γ₁u,Inf),Inf)
            # single and double layer
            S     = SingleLayerOperator(pde,mesh)
            Smat  = S |> Matrix
            D     = DoubleLayerOperator(pde,mesh)
            Dmat = D |> Matrix
            e0   = Nystrom.error_exterior_green_identity(Smat,Dmat,γ₀u,γ₁u)/γ₀u_norm
            Sdim  = Nystrom.assemble_dim(S)
            Ddim  = Nystrom.assemble_dim(D)
            e1    = Nystrom.error_exterior_green_identity(Sdim,Ddim,γ₀u,γ₁u)/γ₀u_norm
            @testset "Single/double layer $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
            # adjoint double-layer and hypersingular
            K     = AdjointDoubleLayerOperator(pde,mesh)
            Kmat     = K |> Matrix
            H     = HyperSingularOperator(pde,mesh)
            Hmat = H |> Matrix
            e0   = Nystrom.error_exterior_derivative_green_identity(Kmat,Hmat,γ₀u,γ₁u)/γ₁u_norm
            Kdim  = Nystrom.assemble_dim(K)
            Hdim  = Nystrom.assemble_dim(H)
            e1   = Nystrom.error_exterior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)/γ₁u_norm
            @testset "Adjoint double-layer/hypersingular $(string(pde))" begin
                @test norm(e0,Inf) > rtol
                @test norm(e1,Inf) < rtol
            end
        end
    end

    @testset "Greens identity (exterior) 3d" begin
        Geometry.clear_entities!()
        Ω   = ParametricSurfaces.Sphere(;radius=3) |> Geometry.Domain
        Γ   = boundary(Ω)
        M   = ParametricSurfaces.meshgen(Γ,(4,4))
        mesh = NystromMesh(view(M,Γ),order=3)
        xs = SVector(0.1,-0.1,0.2)
        ops = (
            Laplace(;dim=3),
            Helmholtz(;dim=3,k=1.2),
            Elastostatic(;dim=3,μ=2,λ=3),
            Maxwell(;k=1)
        )
        for pde in ops
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xs,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xs,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)
            γ₁u_norm = norm(norm.(γ₁u,Inf),Inf)
            # single and double layer
            S     = SingleLayerOperator(pde,mesh)
            Smat  = S |> Matrix
            D     = DoubleLayerOperator(pde,mesh)
            Dmat = D |> Matrix
            e0   = Nystrom.error_exterior_green_identity(Smat,Dmat,γ₀u,γ₁u)/γ₀u_norm
            Sdim  = Nystrom.assemble_dim(S)
            Ddim  = Nystrom.assemble_dim(D)
            e1    = Nystrom.error_exterior_green_identity(Sdim,Ddim,γ₀u,γ₁u)/γ₀u_norm
            @testset "Single/double layer $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
            # adjoint double-layer and hypersingular
            pde isa Maxwell && continue
            K     = AdjointDoubleLayerOperator(pde,mesh)
            Kmat     = K |> Matrix
            H     = HyperSingularOperator(pde,mesh)
            Hmat = H |> Matrix
            e0   = Nystrom.error_exterior_derivative_green_identity(Kmat,Hmat,γ₀u,γ₁u)/γ₁u_norm
            Kdim  = Nystrom.assemble_dim(K)
            Hdim  = Nystrom.assemble_dim(H)
            e1   = Nystrom.error_exterior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)/γ₁u_norm
            @testset "Adjoint double-layer/hypersingular $(string(pde))" begin
                @test norm(e0,Inf) > norm(e1,Inf)
                @test norm(e1,Inf) < rtol
            end
        end
    end
end

ext: jl

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text: ispangram(input::AbstractString) = length(Set(collect((m.match for m = eachmatch(r"[a-z]", lowercase(input)))))) == 26


text:
module gRPC

using Sockets
using ProtoBuf
using HTTP2
using HTTP2.Session

import Base: close, run
import Sockets: TCPServer
import ProtoBuf: call_method
import HTTP2: Headers
import HTTP2.Session: HTTPConnection, EvtGoaway, EvtRecvHeaders, EvtRecvData

include("error.jl")
include("common.jl")
include("server.jl")
include("client.jl")

export gRPCController
export gRPCClient, close, stub
export gRPCServer, close, run
export gRPCException, wrap

end # module


text:
###############################################################################
#
# Optimize leaf invstment Vcmax+Jmax and Leaf area
#
###############################################################################
"""
    optimize_leaf!(
                node::SPACSimple{FT},
                photo_set::AbstractPhotoModelParaSet{FT},
                weather::Array{FT,2},
                printing::Bool
    ) where {FT<:AbstractFloat}

Optimize leaf area (LAI within 0-20) and photosynthetic capacity (within
    5-200), given
- `node` [`SPACSimple`] type struct
- `photo_set` [`AbstractPhotoModelParaSet`] type struct
- `weather` Weather profile in a growing season
- `printing` Optional. If true, printing progress
"""
function optimize_leaf!(
            node::SPACSimple{FT},
            photo_set::AbstractPhotoModelParaSet{FT},
            weather::Array{FT,2},
            printing::Bool = false
) where {FT<:AbstractFloat}
    # 1. use the opt_laba and opt_vmax to initialize
    @inline f(x) = (tmp_node = deepcopy(node);
                    leaf_allocation!(tmp_node, photo_set, x[1], x[2]);
                    tmp_prof = annual_profit(tmp_node, photo_set, weather);
                    if printing
                        @info "\tLABA=$(x[1])\tVMAX=$(x[2])\tPROF=$(tmp_prof)";
                    end;
                    return tmp_prof);

    ms = ReduceStepMethodND{FT}(
                x_mins=FT[1,5],
                x_maxs=FT[node.gaba*20, 200],
                x_inis=FT[node.opt_laba, node.opt_vmax],
                Δ_inis=FT[100,10]);
    st = SolutionToleranceND{FT}(FT[0.9, 0.09], 50);
    lv = find_peak(f, ms, st);

    # 2. update the opt_laba and opt_vmax to node
    leaf_allocation!(node, photo_set, lv[1], lv[2]);
    node.opt_laba = lv[1];
    node.opt_vmax = lv[2];

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module SimpleGF2
using LinearAlgebra

import Base: show, +, -, *, /, xor, &, |, abs
import Base: convert, promote_rule, isless, real, rand

export GF2

struct GF2 <: Number
  val::UInt8
  function GF2(x::T) where T<:Integer
    new(x&1)
  end
end

convert(::Type{GF2}, x::T) where T<:Integer = GF2(x)
convert(::Type{T}, x::GF2) where T<:Integer = T(x.val)

+(x::GF2,y::GF2) = GF2(x.val+y.val)    # addition
-(x::GF2,y::GF2) = GF2(x.val+y.val)    # subtraction
-(x::GF2) = x                          # unary minus
*(x::GF2,y::GF2) = GF2(x.val & y.val)  # multiplication
function /(x::GF2,y::GF2)              # division
  if y.val == 0x00
    error("Division by zero in GF(2)")
  end
  return x
end

# Bitwise operators
xor(x::GF2,y::GF2) = GF2(xor(x.val,y.val))   # XOR
(&)(x::GF2,y::GF2) = GF2(x.val&y.val)   # AND
(|)(x::GF2,y::GF2) = GF2(x.val|y.val)   # OR

# isless defined giving all relations
isless(x::GF2,y::GF2) = isless(x.val,y.val)
isless(x::GF2,y::T) where T<:Real = isless(x.val,y)
isless(x::T,y::GF2) where T<:Real = isless(x,y.val)

abs(x::GF2) = x
real(x::GF2) = x.val

# random values and matrices
rand(::Type{GF2}) = GF2(rand(Int))
#rand(::Type{GF2},dims::Integer...) = map(GF2,rand(Int,dims...))
rand(::Type{GF2},dims::Integer...) = GF2.(rand(Int,dims...))

promote_rule(::Type{GF2}, ::Type{Int} ) = GF2

function show(io::IO, x::GF2)
  if x.val == 0x00
    print(io,"GF2(0)")
  else
    print(io,"GF2(1)")
  end
end

include("solving.jl")
# include("GF2inv.jl") -- now moved to "solving"

end # end of module


text:
include("sd0.jl")

function var0(a)
    # create variable to store ouput
    out = 0.0

    # sd
    s = sd0(a)

    # calculation
    out = s^2

    # return
    return out

end


text:
# This file is a part of DataDeps.jl. License is MIT.
module DataDeps

if VERSION < v"1.3"
    # Load in `deps.jl`, complaining if it does not exist
    const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
    if !isfile(depsjl_path)
        error("DataDeps not installed properly, run `] build DataDeps`, restart Julia and try again")
    end
    include(depsjl_path)

    # Emulate `p7zip_jll` behaviour
    p7zip(f) = f(p7zip_exe)
else
    using p7zip_jll
end

using HTTP
using Reexport
@reexport using SHA

export DataDep, ManualDataDep
export register, resolve, @datadep_str, preupload_check
export unpack

include("errors.jl")
include("types.jl")

include("util.jl")
include("registration.jl")

include("locations.jl")
include("verification.jl")

include("resolution.jl")
include("resolution_automatic.jl")
include("resolution_manual.jl")

include("preupload.jl")

include("fetch_helpers.jl")
include("post_fetch_helpers.jl")

end # module


text:
module OpenBoundaryCondition
    export DVector, Distance, Shift!

    function DVector(obj1::Base.Main.DefectObj, obj2::Base.Main.DefectObj, universe::Base.Main.Universe)
        dr = obj2.position - obj1.position
        dr
    end

    function Distance(obj1::Base.Main.DefectObj, obj2::Base.Main.DefectObj, universe::Base.Main.Universe)
        dVector = DVector(obj1, obj2, universe)
        distance = 0.0
        for d in 1:universe.box.nDimension
            distance += dVector[d]^2
        end
        return distance = distance^0.5
    end

    function Shift!(obj::Base.Main.DefectObj, dr::Base.Main.Vector{Int64}, universe::Base.Main.Universe)
        obj.position += dr
        for d in 1:universe.box.nDimension
            if obj.position[d] < universe.box.boundary[d,1]
                Base.Main.Introduce!(universe, [Base.Main.Disappear()], obj)
            elseif obj.position[d] >= universe.box.boundary[d,2]
                Base.Main.Introduce!(universe, [Base.Main.Disappear()], obj)
            end
        end
    end
end

module PPPBoundaryCondition
    export DVector, Distance, Shift!

    function DVector(obj1::Base.Main.DefectObj, obj2::Base.Main.DefectObj, universe::Base.Main.Universe)
        dr = obj2.position - obj1.position
        for d in 1:universe.box.nDimension
            if dr[d] >= universe.box.halfLength[d]
                dr[d] -= universe.box.length[d]
            elseif dr[d] < -universe.box.halfLength[d]
                dr[d] += universe.box.length[d]
            end
        end
        dr
    end

    function Distance(obj1::Base.Main.DefectObj, obj2::Base.Main.DefectObj, universe::Base.Main.Universe)
        dVector = DVector(obj1, obj2, universe)
        distance = 0.0
        for d in 1:universe.box.nDimension
            distance += dVector[d]^2
        end
        return distance = distance^0.5
    end

    function Shift!(obj::Base.Main.DefectObj, dr::Base.Main.Vector{Int64}, universe::Base.Main.Universe)
        obj.position += dr
        for d in 1:universe.box.nDimension
            if obj.position[d] < universe.box.boundary[d,1]
                obj.position[d] += universe.box.length[d]
            elseif obj.position[d] >= universe.box.boundary[d,2]
                obj.position[d] -= universe.box.length[d]
            end
        end
    end
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mc-ke

max_forks_repo_name: permissionx/mc-ke


text:
using Base.Test

include("cgop.jl")

tcls1 = [
 parse("[x,y].[+P(x),+Q(x,y)]"),
 parse("[y,x].[+P(a),-Q(x,y)]"),
 parse("[x,y].[+Q(x,x),-P(x)]"),
 parse("[].[-P(b)]"),
 parse("[x].[-P(x)]")
]

tcls2 = [
 parse("[x,y].[+P(x),+Q(x,y)]"),
 parse("[y,x].[+P(a),-Q(x,y)]"),
 parse("[x,y].[+Q(x,x),-P(x)]"),
 parse("[].[-P(b),+Q(b,b)]"),
 parse("[x].[-P(x),+Q(x,x)]")
]

tcls3 = [
 parse("[x,y].[+Q(x,x),+R(x),-P(x),+R(x)]"),
 parse("[y,x].[+P(a),-R(a),-Q(x,y),-R(a)]"),
 parse("[x].[-P(b),+R(b),+R(x),+Q(b,b)]"),
 parse("[x,y].[+P(x),+Q(x,y),-R(x),-R(y)]"),
 parse("[y,x].[+P(a),-Q(x,y),-R(a)]"),
 parse("[x,y].[-P(x),+Q(x,x),+R(x,f(a)),+R(x,y)]")
]

@testset "makecdb" begin
 cdb1 = []

 cdb2 = makecdb(cdb1,tcls1[1])
 @test cdb2[1] == ([:x,:y], parse("[+P(x),+Q(x,y)]").args)
 cdb3 = makecdb(cdb2,tcls1[4])
 @test cdb3[2] == ([], parse("[-P(b)]").args)

end

@testset "makeldb" begin
 ldb1 = Dict()
@show tcls1[2]
 ldb2 = makeldb(ldb1, 2, tcls1[2])
@show ldb1
@show ldb2
 @test ldb2[(2,1)] == parse("+P(a)")
 @test ldb2[(2,2)] == parse("-Q(x,y)")
end

@testset "makepgr" begin
 ldb = Dict()
 ldb = makeldb(ldb, 1, tcls1[1])
 gr1 = makepg(ldb)

end

@testset "makedb" begin
 resdb= makedb(tcls1)

 @test length(resdb.cdb) == 5
 @test length(values(resdb.ldb)) == 8
 @test length(keys(resdb.pgr)) == 2
 @test length(values(resdb.pgr)) == 2

 @test getvars(resdb.cdb,1) == [:x,:y]
 @test getvars(resdb.cdb,4) == []

 @test litcounts(resdb.cdb) == [2,2,2,1,1]

 @test getlit(resdb.cdb,(2,2)) == (:-, parse("Q(x,y)"))

end

@testset "putcls" begin
 adb1= makedb(tcls1)
 cdb1=adb1.cdb
 @test putcls(([:x,:y],parse("[+(P(x)),+Q(x,y)]")),cdb1) == 1
 @test putcls(([:w,:y],parse("[+(P(w)),+Q(w,y)]")),cdb1) == 1
 @test putcls(([:w,:y],parse("[+(P(w)),+Q(x,y)]")),cdb1) == size(cdb1,1)+1
 @test putcls(([:x,:y],parse("[+(R(x)),+Q(x,y)]")),cdb1) == size(cdb1,1)+1

end

@testset "findunit" begin
 rdb1= makedb(tcls1)
 rdb2= makedb(tcls2)

 @test findunit(rdb1.cdb) == [4,5]
 @test findunit(rdb2.cdb) == []
end

@testset "putcdb" begin

end

@testset "putldb" begin

end

@testset "putpgr" begin

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text: println("Gaussian Elimination!")


text:
module ExprTools

using MacroTools
using MacroTools: @q

export combinedef, combinearg, splitdef, splitarg

@static if VERSION >= v"0.6.0"
    include("utilsv6.jl")

    """
    combinedef(dict::Dict)

    `combinedef` is the inverse of `splitdef`. It takes a splitdef-like dict
    and returns a function definition. """
    function combinedef(dict::Dict)
        rtype = get(dict, :rtype, :Any)
        # We have to combine params and whereparams because f{}() where {} = 0 is
        # a syntax error unless as a constructor.
        all_params = [get(dict, :params, [])..., get(dict, :whereparams, [])...]
        :(function $(dict[:name]){$(all_params...)}($(dict[:args]...);
                                                    $(dict[:kwargs]...))::$rtype
          $(dict[:body])
          end)
    end
else
    longdef1_where(ex) = ex
    function splitwhere(fdef)
        @assert(@capture(longdef1(fdef),
                         function (fcall_ | fcall_)
                         body_ end),
                "Not a function definition: $fdef")
        return fcall, body, nothing
    end

    """
        combinedef(dict::Dict)

    `combinedef` is the inverse of `splitdef`. It takes a splitdef-like dict
    and returns a function definition. """
    function combinedef(dict::Dict)
        rtype = get(dict, :rtype, :Any)
        params = get(dict, :params, [])
        :(function $(dict[:name]){$(params...)}($(dict[:args]...);
                                                $(dict[:kwargs]...))::$rtype
              $(dict[:body])
          end)
    end
end

function longdef1(ex)
  @match ex begin
    (f_(args__) = body_) => @q function $f($(args...)) $body end
    (f_(args__)::rtype_ = body_) => @q function $f($(args...))::$rtype $body end
    ((args__,) -> body_) => @q function ($(args...),) $body end
    (arg_ -> body_) => @q function ($arg,) $body end
    _ => longdef1_where(ex)
  end
end

doc"""    splitdef(fdef)

Match any function definition

```julia
function name{params}(args; kwargs)::rtype where {whereparams}
   body
end
```

and return `Dict(:name=>..., :args=>..., etc.)`. The definition can be rebuilt by
calling `MacroTools.combinedef(dict)`, or explicitly with the appropriate quoted
expression (good luck).
"""
function splitdef(fdef)
    error_msg = "Not a function definition: $fdef"
    fcall, body, whereparams = splitwhere(fdef)
    @assert(@capture(fcall, ((func_(args__; kwargs__)) |
                             (func_(args__; kwargs__)::rtype_) |
                             (func_(args__)) |
                             (func_(args__)::rtype_))),
            error_msg)
    @assert(@capture(func, (fname_{params__} | fname_)), error_msg)
    di = Dict(:name=>fname, :args=>args,
              :kwargs=>(kwargs===nothing ? [] : kwargs), :body=>body)
    if rtype !== nothing; di[:rtype] = rtype end
    if whereparams !== nothing; di[:whereparams] = whereparams end
    if params !== nothing; di[:params] = params end
    di
end

"""
    combinearg(arg_name, arg_type, is_splat, default)

`combinearg` is the inverse of `splitarg`. """
function combinearg(arg_name, arg_type, is_splat, default)
    a = arg_name===nothing ? :(::$arg_type) : :($arg_name::$arg_type)
    a2 = is_splat ? Expr(:..., a) : a
    return default === nothing ? a2 : Expr(:kw, a2, default)
end

macro splitcombine(fundef)
    dict = splitdef(fundef)
    esc(rebuilddef(striplines(dict)))
end

"""
    splitarg(arg)

Match function arguments (whether from a definition or a function call) such as
`x::Int=2` and return `(arg_name, arg_type, is_splat, default)`. `arg_name` and
`default` are `nothing` when they are absent. For example:

```julia
> map(splitarg, (:(f(a=2, x::Int=nothing, y, args...))).args[2:end])
4-element Array{Tuple{Symbol,Symbol,Bool,Any},1}:
 (:a, :Any, false, 2)
 (:x, :Int, false, :nothing)
 (:y, :Any, false, nothing)
 (:args, :Any, true, nothing)
```
"""
function splitarg(arg_expr)
    splitvar(arg) =
        @match arg begin
            ::T_ => (nothing, T)
            name_::T_ => (name::Symbol, T)
            x_ => (x::Symbol, :Any)
        end
    (is_splat = @capture(arg_expr, arg_expr2_...)) || (arg_expr2 = arg_expr)
    if @capture(arg_expr2, arg_ = default_)
        @assert default !== nothing "splitarg cannot handle `nothing` as a default. Use a quoted `nothing` if possible. (MacroTools#35)"
        return (splitvar(arg)..., is_splat, default)
    else
        return (splitvar(arg_expr2)..., is_splat, nothing)
    end
end

end # module


text:
test_in = b"""
    Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor
    incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud
    exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat
    """

@testset "simple_compression" begin
    comp = lz4_compress(test_in, 5)
    @test comp != test_in

    decomp = lz4_decompress(comp)
    @test decomp == test_in

    comp = lz4_hc_compress(test_in, 12)
    @test comp != test_in
    decomp = lz4_decompress(comp)
    @test decomp == test_in

    huge_vector = rand(UInt8, CodecLz4.LZ4_MAX_INPUT_SIZE + 1)
    @test_throws CodecLz4.LZ4Exception lz4_hc_compress(huge_vector)
    @test_throws CodecLz4.LZ4Exception lz4_compress(huge_vector)
    @test_throws CodecLz4.LZ4Exception lz4_decompress(comp, floor(Int, length(test_in)/2))
    @test_throws CodecLz4.LZ4Exception lz4_decompress(test_in)
end


text:
function reShoot!(sys::System, ens::Ensemble, i::Int64)
    for k in eachindex(ens)
        sample = genstate(sys.data, sys.times[i])
        for i in eachindex(sample)
            ens[k].x[i] = round(Int,sample[i]/ens[k].θ[4+i])
        end
        ens[k].t = sys.times[i]
    end
end

function addFirstEns(sys::System, tauVar::TauVar)
    ens = addEnsemble(sys, method = :Sobol)
    npoints = length(sys.times)
    if sys.routine.ssa == :Direct
        for i in 2:npoints
            if sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            if i > 2 && sys.routine.splitting
                split!(sys, ens)
            end
            simIntDM!(sys, ens, i)
            getCost!(sys, ens, i)
        end
    elseif sys.routine.ssa == :Tau
        for i in 2:npoints
            if sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            if i > 2 && sys.routine.splitting
                split!(sys, ens)
            end
            simIntTau!(sys, ens, i, tauVar)
            getCost!(sys, ens, i)
        end
    end
    ens
end

function addEns(sys::System, tauVar::TauVar, n::Int; gn::Int=0)
    ens = addEnsemble(sys, nSamples = n, method = :Dist, gn=gn)
    npoints = length(sys.times)
    if sys.routine.ssa == :Direct
        for i in 2:npoints
            if sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            if i > 2 && sys.routine.splitting
                split!(sys, ens)
            end
            simIntDM!(sys, ens, i)
            getCost!(sys, ens, i)
        end
    elseif sys.routine.ssa == :Tau
        for i in 2:npoints
            if sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            if i > 2 && sys.routine.splitting
                split!(sys, ens)
            end
            simIntTau!(sys, ens, i, tauVar)
            getCost!(sys, ens, i)
        end
    end
    ens
end

function split!(sys::System, ens::Ensemble)
    inds = Vector{Int}(undef, 0)
    costs = Vector{Float64}(undef, length(ens))
    for k in eachindex(ens)
        costs[k] = getCostVal(ens[k])
    end
    δ = quantile(costs, sys.routine.nElite / length(ens))
    append!(inds, findall(costs .<= δ))
    s = sample(inds, length(ens))
    new_ens = Ensemble()
    ens = setPath(sys,ens,s)
end

function setPath(sys, ens, s)
    new_ens = Ensemble()
    if sys.state.i == 1
        sob = SobolSeq(length(sys.model.bounds[1]), sys.model.bounds[1],  sys.model.bounds[2])
        for i in eachindex(s)
            par = next!(sob)
            push!(new_ens, copyPath2(ens[s[i]], par))
        end
    else
        if sys.routine.sampling == :Normal
            dist = MvNormal(sys.state.θ, sys.state.σ2)
            for i in eachindex(s)
                par = rand(dist)
                push!(new_ens, copyPath2(ens[s[i]], par))
            end
        elseif sys.routine.sampling == :log
            dist = MvLogNormal(sys.state.θ, sys.state.σ2)
            for i in eachindex(s)
                par = rand(dist)
                push!(new_ens, copyPath2(ens[s[i]], par))
            end
        end
    end
    new_ens
end

function copyPath2(p::Path, par::Vector{Float64})
    p2 = copyPath(p)
    p2.θ = copy(par)
    p2
end

function getEliteData(sys::System, ens::Ensemble)
    eliteSet = Ensemble()
    mincost = Inf
    costs = Vector{Float64}(undef, length(ens))
    for k in eachindex(ens)
        costs[k] = getCostVal(ens[k])
        costs[k] < mincost && (mincost = costs[k])
    end
    δ = quantile(costs, sys.routine.nElite / length(ens))
    for xx in ens[costs .<= δ]
        push!(eliteSet, xx)
    end
    return eliteSet, δ, mincost
end

function getEliteGroups(sys::System, ens::Ensemble)
    ng = ens[end].group
    groupcosts = zeros(ng)
    for g in 1:ng
        mincost = Inf
        for k in eachindex(ens)
            (g == ens[k].group) && (groupcosts[g] += getCostVal(ens[k]))
        end
    end
    δ = quantile(groupcosts, sys.routine.nElite / ng)
    gelite = collect(1:ng)[groupcosts .<= δ]

    groupSet = Ensemble()

    for xx in ens
        if in(xx.group, gelite)
            push!(groupSet, xx)
        end
    end

    return groupSet, δ, minimum(groupcosts)
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
# Frequency
using Base.Dates

abstract Frequency

type NoFrequency <: Frequency end
type Once <: Frequency end
type Annual <: Frequency end
type Semiannual <: Frequency end
type EveryFourthMonth <: Frequency end
type Quarterly <: Frequency end
type Bimonthly <: Frequency end
type Monthly <: Frequency end
type EveryFourthWeek <: Frequency end
type Biweekly <: Frequency end
type Weekly <: Frequency end
type Daily <: Frequency end
type OtherFrequency <: Frequency end

# default value if no freq specified
value(::Frequency) = -1

value(::NoFrequency) = -1
value(::Once)			 	= 0
value(::Annual)			= 1
value(::Semiannual)		= 2
value(::EveryFourthMonth)  = 3
value(::Quarterly)		 	= 4
value(::Bimonthly)		 	= 6
value(::Monthly)			= 12
value(::EveryFourthWeek)  	= 13
value(::Biweekly)		 	= 26
value(::Weekly)			= 52
value(::Daily)			 	= 365
value(::OtherFrequency)   	= 999

period(::Annual) = Year(1)
period(::Semiannual) = Month(6)


text:
include("sn.jl")
include("sn_bigr.jl")

#r=10.
r=100.
#r=110.
#r=1.1
#r=0.01
l_max = 20
n_max = l_max^2+2*l_max
sn_big = zeros(BigFloat,n_max+1);
snew = zeros(typeof(r),n_max+1);
snew_big = zeros(BigFloat,n_max+1);
#b=sqrt((r+.5)^2+0.1^2)
#b=1e-3
#b=0.99
#b=0.9
#b=0.55
b = 1+r-1e-8
#b = 1-r+1e-8
#b = 110.99
#b = 1.0
#b=100.
#b = r+1-1e-18
diff_frac = zeros(n_max+1)
diff_rel  = zeros(n_max+1)
for i=1:length(b)
  s_n_bigr!(l_max,r,b[i],snew)
  s_n_bigr!(l_max,big(r),big(b[i]),snew_big)
  s_n!(l_max,big(r),big(b[i]),sn_big)
  l = 0; n = 0; m = 0
  for n=0:n_max
    mu = l-m; nu = l+m
    if snew[n+1] != 0.0
      diff = abs(snew[n+1]/convert(Float64,snew_big[n+1])-1.)
      diff_frac[n+1]=diff
      diff_rel[n+1]=abs(snew[n+1]-convert(Float64,snew_big[n+1]))
#      if diff > 0.01
      println("n: ",n," l: ",l," m: ",m," mu: ",mu," nu: ",nu," s: ",snew[n+1],
#       " s_old: ",convert(Float64,sn_big[n+1])," d_old: ",snew[n+1]/convert(Float64,sn_big[n+1])-1.,
       " s_old: ",convert(Float64,sn_big[n+1]),
       " s_big: ",convert(Float64,snew_big[n+1])," d_old: ",snew[n+1]/convert(Float64,sn_big[n+1])-1.,
       " d_big: ",snew[n+1]/convert(Float64,snew_big[n+1])-1.)
#      end
    end
    m +=1
    if m > l
      l += 1
      m = -l
    end
  end
end

using PyPlot
clf()
semilogy(diff_frac,"o",label="fractional",c="blue")
semilogy(diff_rel,"o",label="relative",c="orange")
legend(loc="upper left")
xlabel("n")
ylabel("Error")
axis([-5,n_max,1e-25,1e-2])

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
const JSONRPCID = Union{AbstractString, Integer}

@enum(
    ErrorCode,
    ParseError = -32700,
    InvalidRequest = -32600,
    MethodNotFound = -32601,
    InvalidParams = -32602,
    InternalError = -32603,
    # -32000 to -32099: Server error
    MethodException = -32001,
    SerializationError = -32002,
    UserError = -40000,
)

abstract type Response end

struct ResultResponse <: Response
    id::JSONRPCID
    result
    jsonrpc::Symbol
    ResultResponse(id, result) = new(id, result, Symbol("2.0"))
end

struct Error
    code::Int
    message::String
    data
end

Error(code::ErrorCode, message, data) = Error(Int(code), message, data)
Error(code, message) = Error(code, message, nothing)

usererror(message::AbstractString, data=nothing) =
    Error(UserError, message, data)

struct ErrorResponse <: Response
    id::Union{JSONRPCID, Nothing}
    error::Error
    jsonrpc::Symbol
    ErrorResponse(id, error) = new(id, error, Symbol("2.0"))
end

ErrorResponse(code::ErrorCode) =
    ErrorResponse(nothing, Error(code, string(code)))

ErrorResponse(id, code::ErrorCode) =
    ErrorResponse(id, Error(code, string(code)))

ErrorResponse(id, code::ErrorCode, message::String, data=nothing) =
    ErrorResponse(id, Error(code, message, data))

function jsonrpc(dispatcher, rf, acc, input::AbstractDict{String})
    @debug "jsonrpc" dispatcher rf acc input

    id = get(input, "id", nothing)
    if !(id isa String || id isa Int)
        id = nothing
    end
    @debug "`id` field is parsed" id
    if get(input, "jsonrpc", nothing) != "2.0"
        return rf(acc, ErrorResponse(id, InvalidRequest))
    end
    @debug "`jsonrpc` field is correct"

    fname = get(input, "method", nothing)
    if !(fname isa String)
        return rf(acc, ErrorResponse(id, InvalidRequest))
    end
    @debug "`method` field is correct"

    params = get(input, "params", nothing)
    if !(params isa Vector || params isa Dict || params isa Nothing)
        return rf(acc, ErrorResponse(id, InvalidRequest))
    end
    @debug "`params` field is correct"

    @debug "Request parsed" id fname params
    id === nothing && return acc  # notification
    @debug "Request is valid"

    f = get(dispatcher, fname, nothing)
    if f === nothing
        return rf(acc, ErrorResponse(id, MethodNotFound))
    end
    @debug "Method is found"

    result, success = try
        f(params), true
    catch err
        @debug "Failed" f params exception=(err, catch_backtrace())
        err, false
    end
    @debug "got" result success
    if success
        if result isa Exception
            return rf(acc, ErrorResponse(id, UserError, sprint(showerror, result)))
        elseif result isa Error
            return rf(acc, ErrorResponse(id, result))
        else
            return rf(acc, ResultResponse(id, result))
        end
    else
        return rf(acc, ErrorResponse(id, MethodException, sprint(showerror, result)))
    end
end

function jsonrpc(dispatcher, rf, acc, input::AbstractVector)
    isempty(input) && return rf(acc, ErrorResponse(InvalidRequest))
    results = foldl([], init=input) do acc, input
        jsonrpc(dispatcher, rf, acc, input)
    end
    isempty(results) && return acc  # notification
    return rf(acc, results)
end

jsonrpc(dispatcher, rf, acc, response::Response) = rf(acc, response)

jsonrpc(dispatcher, rf, acc, _) =
    rf(acc, ErrorResponse(InvalidRequest))


text:
using Oceananigans.Grids: Flat, Bounded, y_domain
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: _top_tke_flux, CATKEVDArray

import Oceananigans.Grids: validate_size, validate_halo, HRegRectilinearGrid
import Oceananigans.TurbulenceClosures: time_discretization, calculate_diffusivities!, with_tracers
import Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∂ⱼ_τ₃ⱼ, ∇_dot_qᶜ
import Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: top_tke_flux
import Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U

#####
##### Implements a "single column model mode" for HydrostaticFreeSurfaceModel
#####

const YZSliceGrid = Union{AbstractGrid{<:AbstractFloat, <:Flat, <:Bounded, <:Bounded},
                          AbstractGrid{<:AbstractFloat, <:Flat, <:Periodic, <:Bounded}}

@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid::YZSliceGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i]
    return ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid::YZSliceGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i]
    return ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∇_dot_qᶜ(i, j, k, grid::YZSliceGrid, closure_array::ClosureArray, c, tracer_index, args...)
    @inbounds closure = closure_array[i]
    return ∇_dot_qᶜ(i, j, k, grid, closure, c, tracer_index, args...)
end

struct SliceEnsembleSize
    ensemble :: Int
    Ny :: Int
    Nz :: Int
    Hy :: Int
    Hz :: Int
end

SliceEnsembleSize(; size, ensemble=0, halo=(1, 1)) = SliceEnsembleSize(ensemble, size[1], size[2], halo[1], halo[2])

validate_size(TX, TY, TZ, e::SliceEnsembleSize) = tuple(e.ensemble, e.Ny, e.Nz)
validate_halo(TX, TY, TZ, e::SliceEnsembleSize) = tuple(0, e.Hy, e.Hz)

#####
##### CATKEVerticalDiffusivity helpers
#####

""" Compute the flux of TKE through the surface / top boundary. """
@inline function top_tke_flux(i, j, grid::YZSliceGrid, clock, fields, parameters, closure_array::CATKEVDArray, buoyancy)
    top_tracer_bcs = parameters.top_tracer_boundary_conditions
    top_velocity_bcs = parameters.top_velocity_boundary_conditions
    @inbounds closure = closure_array[i]

    return _top_tke_flux(i, j, grid, closure.surface_TKE_flux, closure,
                         buoyancy, fields, top_tracer_bcs, top_velocity_bcs, clock)
end

@inline function hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid::YZSliceGrid,
                                                               val_tracer_index::Val{tracer_index},
                                                               advection,
                                                               closure_array::CATKEVDArray, args...) where tracer_index

    @inbounds closure = closure_array[i]
    return hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid, val_tracer_index, advection, closure, args...)
end

#####
##### Arrays of Coriolises
#####

const CoriolisVector = AbstractVector{<:AbstractRotation}

@inline x_f_cross_U(i, j, k, grid::YZSliceGrid, coriolis::CoriolisVector, U) = @inbounds x_f_cross_U(i, j, k, grid, coriolis[i], U)
@inline y_f_cross_U(i, j, k, grid::YZSliceGrid, coriolis::CoriolisVector, U) = @inbounds y_f_cross_U(i, j, k, grid, coriolis[i], U)
@inline z_f_cross_U(i, j, k, grid::YZSliceGrid, coriolis::CoriolisVector, U) = @inbounds z_f_cross_U(i, j, k, grid, coriolis[i], U)

function FFTImplicitFreeSurfaceSolver(arch, grid::YZSliceGrid, settings)

    grid isa HRegRectilinearGrid ||
        throw(ArgumentError("FFTImplicitFreeSurfaceSolver requires horizontally-regular rectilinear grids."))

    # Construct a "horizontal grid". We support either x or y being Flat, but not both.
    TY = topology(grid, 2)

    sz = SliceEnsembleSize(size=(grid.Ny, 1), ensemble=grid.Nx, halo=(grid.Hy, 0))

    horizontal_grid = RectilinearGrid(; topology = (Flat, TY, Flat),
                                        size = sz,
                                        halo = grid.Hy,
                                        y = y_domain(grid))

    solver = FFTBasedPoissonSolver(arch, horizontal_grid)
    right_hand_side = solver.storage

    return FFTImplicitFreeSurfaceSolver(solver, grid, horizontal_grid, right_hand_side)
end


text:
function parseInput(file)
    allNumbers::Vector{Int} = []
    open(file) do file
        for line in eachline(file)
            push!(allNumbers, parse(Int,line))
        end
    end
    return allNumbers
end

function testNumberAtIndex(allNumbers, numberIndex, preambleLength)
    result = 0
    numberToTest = allNumbers[numberIndex]
    lastIx = numberIndex-1
    startIx = numberIndex-preambleLength
    for i in startIx:1:lastIx
        for j in startIx:1:lastIx
            a = allNumbers[i]
            b = allNumbers[j]
            # println("a+b=", a+b)
            if a != b && a + b == numberToTest
                result += 1
            end
        end
    end
    return result > 0
end

function findWeakness(file, preambleLength)
    allNumbers = parseInput(file)
    for n in preambleLength+1:length(allNumbers)
        # println(allNumbers[n])
        if !testNumberAtIndex(allNumbers,n,preambleLength)
            println("Found weaknes in ", allNumbers[n], " at index ", n)
        end
    end
end

findWeakness("data.txt", 25)


text:
# This file contains the definition of a function
# to generate an example of the input/output file

"""
    example(file::AbstractString)

Generate an example of the input/output file.

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Input.example(file)

# output

```

See: [`InputStruct`](@ref).
"""
function example(file::AbstractString)

    # Strip the string
    file = strip(file)

    # Check if it's a directory
    if isdir(file)
        throw(ScatsInputIsADir(file))
    end

    open(file, "w") do io
        # Print
        println(io, "Sample size")
        println(io, "230")
        println(io)
        println(io, "Sample step")
        println(io, RT(1))
        println(io)
        println(io, "Significance level")
        println(io, RT(0.01))
        println(io)
        println(io, "Time array")
        println(io, Array{RT}([i for i in 0:229]))
        println(io)
        println(io, "Values array")
        println(io, Array{RT}([i for i in 0:229]))
    end

end


text:
plt = Observable{Any}()
ui = dom"div"(plt)
# Ploting recon
map!(t-> begin
        @manipulate for slice = 1:size(image,3)
            aux = log.(abs.(kspace[:,:,slice]).+1)
            plot_image(aux,zmin=0,zmax=.1*maximum(aux[:]);darkmode,title="Reconstruction ($slice/$(size(image,3)))")
        end
    end
    , plt, img_obs)
content!(w, "div#content", ui)


text:
function mean(x::Vector{T}) where {T}
    m = sum(x) / length(x)
    return m
end

function std(x::Vector{T}) where {T}
    m = mean(x)
    z = x .- m
    m = sqrt( sum(z.^2) / (length(x) - 1) )
    return m
end

function center(x::Vector)
    m = mean(x)
    y = copy(x)
    y .-= m
    return y
end

function cov(x::Vector{T}) where {T}
    return std(x)
end

function cov(x::Vector{T}, y::Vector{T}) where {T}
    mx = mean(x)
    my = mean(y)
    L = length(x)

    return 1 / (L - 1) * dot(x .- mx, y .- my)
end

function cov(x::Matrix{T}) where {T}
    # covariance in column way
    L = size(x,1)
    return 1/(L-1) * (x' * x)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
mutable struct Data
	geneIDs
	PPDE
	logFC

	function Data(geneIDs, PPDE, logFC)
		new(geneIDs, PPDE, logFC)
	end
end

"""
Read DE data from a delimited file
"""
function readData(pathToData; delimiter='\t', geneIDhead="GeneID", PPDEhead="PPDE", logFChead="logFC")
	(data, cNames) = readdlm(pathToData, delimiter, header=true)
	geneIDs = data[:, findfirst(isequal(geneIDhead), cNames[:])][:]
	PPDE = data[:, findfirst(isequal(PPDEhead), cNames[:])][:]
	FC = data[:, findfirst(isequal(logFChead), cNames[:])][:]

	return Data(geneIDs, PPDE, FC)
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
RpcApi.SimpleOrganism[RpcApi.SimpleOrganism(0x0000000000090ddf,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_33::Int64 = 3761244718335218341
        function func_31(var_29::String="kxAJlYap",var_30::Int16=6327)
            o.mem[var_30] = var_30
            Creature.stepLeft(o)
            return var_29
        end
        local var_28::String = "BI2aHdRl"
        local var_20::String = "ozWRZlps"
        function func_7(var_6::String="ZKRF9L8v")
            local var_37::String = "AaTltb79"
            Creature.stepDown(o)
            return var_6
        end
        var_33 = (Int64)(var_33) | (Int64)(var_33)
        var_33 = Creature.energyDown(o)
        var_33 = Creature.energyUp(o)
        var_2 = if 127 >= 0.385419
                Int8(round(0.385419))
            else
                127
            end
        for i::Int8 = 1:div(var_2,16)
        end
        var_2 = if 127 >= 0.388149
                Int8(round(0.388149))
            else
                127
            end
        var_2 = if 127 >= abs(30028)
                Int8(30028)
            else
                127
            end
        Creature.eatDown(c,o,Int(var_2))
        var_33 = Int(Creature.idRight(o))
        var_33 = if isempty(var_20)
                0
            else
                1
            end
        Creature.eatDown(c,o,Int(var_2))
        Creature.eatDown(c,o,Int(var_2))
        Creature.eatRight(c,o,Int(var_2))
        var_33 = Creature.energyUp(o)
        Creature.eatUp(c,o,Int(var_2))
        Creature.eatLeft(c,o,Int(var_2))
        var_33 = Creature.energyUp(o)
        Creature.eatRight(c,o,Int(var_2))
        Creature.eatUp(c,o,Int(var_2))
        var_33 = Creature.energyUp(o)
        var_20 = string(24637)
        Creature.eatUp(c,o,Int(var_2))
        var_20 = string(var_20)
        var_33 = Int(Creature.idLeft(o))
        Creature.eatUp(c,o,Int(var_2))
        Creature.eatRight(c,o,Int(var_2))
        Creature.eatLeft(c,o,Int(var_2))
        var_33 = Creature.energyUp(o)
        var_33 = Creature.energyLeft(o)
        var_20 = func_7(var_28)
        var_33 = Creature.energyLeft(o)
        var_33 = Creature.energyDown(o)
        var_20 = var_20 * var_20
        var_2 = if 127 >= abs(var_2)
                Int8(var_2)
            else
                127
            end
        var_2 = if 127 >= abs(var_33)
                Int8(var_33)
            else
                127
            end
        for i::Int8 = 1:div(var_2,16)
            local var_36::Int8 = 29
            var_36 = if 127 >= abs(1571)
                    Int8(1571)
                else
                    127
                end
            Creature.eatDown(c,o,Int(var_36))
        end
        var_13 = Int16(Creature.eatDown(c,o,Int(var_2)))
        var_20 = string(-3373)
        Creature.eatRight(c,o,Int(var_2))
        var_33 = Int(Creature.idLeft(o))
        Creature.eatDown(c,o,Int(var_2))
        return true
    end),[98,38,30,73,83,1,72,0,51],47.4821,0,31.5481,298580,0x0000006b,Dict(542=>542,540=>540,546=>546,545=>545,544=>544,6327=>6327,541=>541,543=>543),[557,490],13144,51,53)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
predict_links(S::Int64, n::Int64, chain)
    S: number of species
    n: number of links that will be predicted
    chain: Monte-Carlo chain of the flexible links model fitted to empirical data
Returns n values of predicted numbers of links
"""
function predict_links(S::Int64, n::Int64, chain)
    predicted_links = zeros(Int64, n) # create vector for the predicted links
    p = get_params(chain[200:end,:,:]) # get fitted parameters of the flexible links model
    # simulate n number of links
    for i in 1:n
        # select parameter values randomly
        p_i = rand(1:length(p.μ))
        μ, ϕ = p.μ[p_i], p.ϕ[p_i]
        # predict number of links
        predicted_links[i] = rand(BetaBinomial(S^2-(S-1), μ*exp(ϕ), (1-μ)*exp(ϕ))) + (S-1)
    end
    return predicted_links
  end


text:
"""
    fit_tensors(directions, adc)

Fit effective diffusion tensors to directionalized ADCs.
The six components of the symmetric diffusion tensors are fitted by least
squares to the gradient directions and resulting ADCs.
"""
function fit_tensors(directions, adc)

    # Sizes
    ncompartment, nsequence, ndirection = size(adc)

    # Get vectors of gradient interactions
    g = directions
    G = zeros(ndirection, 6)
    G[:, 1] = g[1, :] .^ 2
    G[:, 2] = g[2, :] .^ 2
    G[:, 3] = g[3, :] .^ 2
    G[:, 4] = 2 .* g[1, :] .* g[2, :]
    G[:, 5] = 2 .* g[1, :] .* g[3, :]
    G[:, 6] = 2 .* g[2, :] .* g[3, :]

    tensors = Array{SMatrix{3,3},2}(undef, ncompartment, nsequence)
    for icmpt = 1:ncompartment, iseq = 1:nsequence

        # Deduce vector of effective diffusion tensor components usinag a least squares fit
        # of the directions to the computed ADCs
        Dvec = G \ adc[icmpt, iseq, :]

        # Tensorize component vectors (symmetric)
        tensors[icmpt, iseq] = SMatrix{3,3}(
            [
                Dvec[1, :] Dvec[4, :] Dvec[5, :]
                Dvec[4, :] Dvec[2, :] Dvec[6, :]
                Dvec[5, :] Dvec[6, :] Dvec[3, :]
            ],
        )

    end

    tensors
end


text:
export blendMode

function blendMode(mode)
	if mode == "REPLACE"
		glBlendEquation(GL_FUNC_ADD)
		glBlendFunc(GL_ONE, GL_ZERO)
	elseif mode == "BLEND"
		# glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
		glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE)
	elseif mode == "ADD"
		# glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
		glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE, GL_ONE, GL_ONE)
	elseif mode == "SUBTRACT"
		# glBlendEquationSeparate(GL_FUNC_REVERSE_SUBTRACT, GL_FUNC_ADD);
		glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE, GL_ONE, GL_ONE)
	elseif mode == "LIGHTEST"
		# glBlendEquationSeparate(GL_FUNC_MAX, GL_FUNC_ADD);
		glBlendFuncSeparate(GL_ONE, GL_ONE, GL_ONE, GL_ONE)
	elseif mode == "DARKEST"
		# glBlendEquationSeparate(GL_FUNC_MIN, GL_FUNC_ADD);
		glBlendFuncSeparate(GL_ONE, GL_ONE, GL_ONE, GL_ONE)
	elseif mode == "EXCLUSION"
		# glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
		glBlendFuncSeparate(GL_ONE_MINUS_DST_COLOR, GL_ONE_MINUS_SRC_COLOR, GL_ONE, GL_ONE)
	elseif mode == "MULTIPLY"
		# glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
		glBlendFuncSeparate(GL_ZERO, GL_SRC_COLOR, GL_ONE, GL_ONE)
	elseif mode == "SCREEN"
		# glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
		glBlendFuncSeparate(GL_ONE_MINUS_DST_COLOR, GL_ONE, GL_ONE, GL_ONE)
	end
end

#createGraphics


text:
using Test, YaoSym

@testset "apply" begin
    @vars θ γ η
    reg2 = rand_state(4)
    reg1 = ArrayReg(Basic.(state(reg2)))
    for G in [X, ConstGate.T, H]
        g = control(4, 3, 2=>G)
        @test copy(reg1) |> g ≈ copy(reg2) |> g
        # print reg2 errors
    end

    reg2 = ArrayReg(bit"0011")
    reg1 = ArrayReg(Basic.(state(reg2)))
    for GC in [Rx, Ry, Rz, shift, phase, θ->rot(SWAP,θ)]
        G = GC(θ)
        locs = nqubits(G) == 1 ? 2 : (3,1)
        g1 = put(4, locs=>GC(θ))
        g2 = put(4, locs=>GC(π/2))
        regs = state(copy(reg1) |> g1)
        res = subs.(regs, θ, Basic(π)/2)
        @test  res ≈ state(copy(reg2) |> g2)
    end

    CRk(::Type{T}, i::Int, j::Int, k::Int) where T = control([i, ], j=>shift(2*T(π)/(1<<k)))
    CRot(::Type{T}, n::Int, i::Int) where T = chain(n, i==j ? put(i=>H) : CRk(T, j, i, j-i+1) for j = i:n)
    qft(::Type{T}, n::Int) where T = chain(n, CRot(T, n, i) for i = 1:n)
    res = ket"11" |> qft(Basic, 2)
    l1 = length(string(res.state.nzval[2]))
    res.state.nzval .= simplify_expi.(res.state.nzval)
    @test length(string(res.state.nzval[2])) < l1
    @test res ≈ ArrayReg(bit"11") |> qft(Float64, 2)
end


text:
ENV["NEMO_THREADED"] = 1

using Arblib, Test, LinearAlgebra, SpecialFunctions

@testset "Arblib" begin
    include("arb_types-test.jl")
    include("types-test.jl")
    include("arbcall-test.jl")
    include("precision-test.jl")
    include("setters-test.jl")
    include("constructors-test.jl")
    include("predicates-test.jl")
    include("show-test.jl")
    include("promotion.jl")
    include("examples.jl")
    include("arithmetic.jl")
    include("rand.jl")
    include("float.jl")
    include("interval.jl")
    include("vector.jl")
    include("matrix.jl")
    include("eigen.jl")
    include("calc_integrate.jl")
    include("ref-test.jl")
    include("poly.jl")
    include("series.jl")
    include("special-functions.jl")
    include("threading.jl")
end


text:
"""
    BilinearUpsample(factor::Tuple{Integer,Integer})

Create an upsampling layer that uses bilinear interpolation.

The width and height dimensions grow by the `factor` tuple.

# Examples
```jldoctest; setup = :(using Flux: BilinearUpsample; using Random; Random.seed!(0))
julia> b = Flux.BilinearUpsample((2, 2))
BilinearUpsample(2, 2)

julia> b(rand(2, 2, 1, 1))
4×4×1×1 Array{Float64,4}:
[:, :, 1, 1] =
 0.823648  0.658877  0.329336  0.164566
 0.845325  0.675933  0.337149  0.167757
 0.888679  0.710044  0.352773  0.174138
 0.910357  0.7271    0.360586  0.177329```
"""
struct BilinearUpsample{T<:Integer}
    factor::Tuple{T,T}
end

function (b::BilinearUpsample)(x::AbstractArray)
    W, H, C, N = size(x)

    newW = W * b.factor[1]
    newH = H * b.factor[2]

    out = similar(x, (newW, newH, C, N))

    @inbounds for w = 1:newW, h = 1:newH
        w₀ = (w - 0.5) / b.factor[1] + 0.5
        h₀ = (h - 0.5) / b.factor[2] + 0.5

        w1 = floor(Int, w₀)
        w2 = w1 + 1
        h1 = floor(Int, h₀)
        h2 = h1 + 1

        i1 = clamp(w1, 1, W)
        i2 = clamp(w2, 1, W)
        j1 = clamp(h1, 1, H)
        j2 = clamp(h2, 1, H)

        @views out[w, h, :, :] =
            (
                x[i1, j1, :, :] * (w2 - w₀) * (h2 - h₀) +
                x[i1, j2, :, :] * (w2 - w₀) * (h₀ - h1) +
                x[i2, j1, :, :] * (w₀ - w1) * (h2 - h₀) +
                x[i2, j2, :, :] * (w₀ - w1) * (h₀ - h1)
            ) / ((w2 - w1) * (h2 - h1))
    end

    out
end

function Base.show(io::IO, b::BilinearUpsample)
    print(io, "BilinearUpsample(", b.factor[1], ", ", b.factor[2], ")")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 273c7c85-8178-44a7-99f0-581754aeb8c8
begin
	"""
	Mark a vector as being append-only: let Firebasey know that it can diff this array simply by comparing lengths, without looking at its contents.

	It was made specifically for logs: Logs are always appended, OR the whole log stream is reset. AppendonlyMarker is like SubArray (a view into another array) except we agree to only ever append to the source array. This way, firebase can just look at the index and diff based on that.
	"""
	struct AppendonlyMarker{T} <: AbstractVector{T}
		mutable_source::Vector{T}
		length_at_time_of_creation::Int
	end
	AppendonlyMarker(arr) = AppendonlyMarker(arr, length(arr))

	# Poor mans vector-proxy
	# I think this is enough for Pluto to show, and for msgpack to pack
	function Base.size(arr::AppendonlyMarker)
		return (arr.length_at_time_of_creation,)
	end
	Base.getindex(arr::AppendonlyMarker, index::Int) = arr.mutable_source[index]
	Base.iterate(arr::AppendonlyMarker, args...) = Base.iterate(arr.mutable_source, args...)
end

# ╔═╡ 183cef1f-bfe9-42cd-8239-49e9ed00a7b6
md"""
## AppendonlyMarker(s)

Example of how to solve performance problems with Firebasey:
We make a new type with a specific diff function.
It might be very specific per problem, but that's fine for performance problems (I think).
It also keeps the performance solutions as separate modules/packages to whatever it is you're actually modeling.
"""

# ╔═╡ 971709de-074e-49cf-8bd4-9c675b037dfd
md"## `@skip_as_script`"

# ╔═╡ 9fce9aa9-d3c6-4134-9692-8a8756fa3cff
function is_inside_pluto(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ db9167c4-7ba7-42e1-949b-0ad18e2d7b25
"""
	@skip_as_script expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex)
	if is_inside_pluto(__module__)
		esc(ex)
	else
		nothing
	end
end

# ╔═╡ 35d3bcd7-af51-466a-b4c4-cc055e74d01d
@skip_as_script appendonly_1, appendonly_2 = let
	array_1 = [1,2,3,4]
	appendonly_1 = AppendonlyMarker(array_1)
	push!(array_1, 5)
	appendonly_2 = AppendonlyMarker(array_1)

	appendonly_1, appendonly_2
end;

# ╔═╡ 1017f6cc-58ac-4c7b-a6d0-a03f5e387f1b
@skip_as_script appendonly_1_large, appendonly_2_large = let
	large_array_1 = [
		Dict{String,Any}(
			"x" => 1,
			"y" => [1,2,3,4],
			"z" => "hi",
		)
		for i in 1:10000
	];
	appendonly_1 = AppendonlyMarker(large_array_1)
	push!(large_array_1, Dict("x" => 5))
	appendonly_2 = AppendonlyMarker(large_array_1)

	appendonly_1, appendonly_2
end;

# ╔═╡ 33facb62-af28-4e94-946d-545637f320e9
"The opposite of `@skip_as_script`"
macro only_as_script(ex) is_inside_pluto(__module__) ? nothing : esc(ex) end

# ╔═╡ 092c4b11-8b75-446f-b3ad-01fa858daebb
# Only define this in Pluto - assume we are `using Test` otherwise
begin
	@skip_as_script begin
		import Pkg
		Pkg.activate(mktempdir())
		Pkg.add(Pkg.PackageSpec(name="PlutoTest"))
		using PlutoTest
	end
	# Do nothing inside pluto (so we don't need to have Test as dependency)
	# test/Firebasey is `using Test` before including this file
	@only_as_script begin
		if !isdefined(@__MODULE__, Symbol("@test"))
			macro test(e...) nothing; end
			macro test_throws(e...) nothing; end
			macro test_broken(e...) nothing; end
			macro testset(e...) nothing; end
		end
	end
end

# ╔═╡ 06492e8d-4500-4efe-80ee-55bf1ee2348c
@skip_as_script @test length([AppendonlyMarker([1,2,3])...]) == 3

# ╔═╡ 2284ae12-5b8c-4542-81fa-c4d34f2483e7
@skip_as_script @test length([AppendonlyMarker([1,2,3], 1)...]) == 1

# ╔═╡ dc5cd268-9cfb-49bf-87fb-5b7db4fa6e3c
md"## Import Firebasey when running inside notebook"

# ╔═╡ 0c2f23d8-8e98-47b7-9c4f-5daa70a6c7fb
# OH how I wish I would put in the time to refactor with fromFile or SOEMTGHINLAS LDKJ JULIA WHY ARE YOU LIKE THIS GROW UP
@skip_as_script Firebasey = let
	wrapper_module = Module()
	Core.eval(wrapper_module, :(module Firebasey
			include("Firebasey.jl")
		end
	))
	wrapper_module.Firebasey
end

# ╔═╡ 2903d17e-c6fd-4cea-8585-4db26a00b0e7
function Firebasey.diff(a::AppendonlyMarker, b::AppendonlyMarker)
	if a.mutable_source !== b.mutable_source
		[Firebasey.ReplacePatch([], b)]
	else
		if a.length_at_time_of_creation > b.length_at_time_of_creation
			throw(ErrorException("Not really supposed to diff AppendonlyMarker with the original being longer than the next version (you know, 'append only' and al)"))
		end

		map(a.length_at_time_of_creation+1:b.length_at_time_of_creation) do index
			Firebasey.AddPatch([index], b.mutable_source[index])
		end
	end
end

# ╔═╡ 8537488d-2ff9-42b7-8bfc-72d43fca713f
@skip_as_script @test array_diff(appendonly_1, appendonly_2) == [Firebasey.AddPatch([5], 5)]

# ╔═╡ 70179239-357a-424d-bac3-3a1431aff536
var"@track" = Firebasey.var"@track"

# ╔═╡ 37fe8c10-09f0-4f72-8cfd-9ce044c78c13
@skip_as_script @track for _ in 1:1000 array_diff(appendonly_1_large, appendonly_2_large) end

# ╔═╡ Cell order:
# ╟─183cef1f-bfe9-42cd-8239-49e9ed00a7b6
# ╠═273c7c85-8178-44a7-99f0-581754aeb8c8
# ╠═2903d17e-c6fd-4cea-8585-4db26a00b0e7
# ╠═35d3bcd7-af51-466a-b4c4-cc055e74d01d
# ╠═1017f6cc-58ac-4c7b-a6d0-a03f5e387f1b
# ╠═06492e8d-4500-4efe-80ee-55bf1ee2348c
# ╠═2284ae12-5b8c-4542-81fa-c4d34f2483e7
# ╠═8537488d-2ff9-42b7-8bfc-72d43fca713f
# ╠═37fe8c10-09f0-4f72-8cfd-9ce044c78c13
# ╟─971709de-074e-49cf-8bd4-9c675b037dfd
# ╟─9fce9aa9-d3c6-4134-9692-8a8756fa3cff
# ╟─db9167c4-7ba7-42e1-949b-0ad18e2d7b25
# ╟─33facb62-af28-4e94-946d-545637f320e9
# ╟─dc5cd268-9cfb-49bf-87fb-5b7db4fa6e3c
# ╠═0c2f23d8-8e98-47b7-9c4f-5daa70a6c7fb
# ╠═70179239-357a-424d-bac3-3a1431aff536
# ╠═092c4b11-8b75-446f-b3ad-01fa858daebb


text:
using ApproxFun, Plots

# In this example, we numerically compute with Galerkin orthogonal polynomials.
# If H is a separable Hilbert space and B contains ν linear functional constraints,
# then the Galerkin orthogonal polynomials are those orthogonal polynomials with
# respect to H_B := \{ u \in H : Bu = 0 \}.
# For example, if H = L^2([-1,1]) and we wish to enforce that u(-1) = u'(-1) =
# u'(1) = 0, then the 2-normalized orthonormal polynomials may be created by
# taking the QR factorization of the connection coefficients between the
# normalized Legendre polynomials \sqrt{n+1/2}P_n(x) and a degree-graded
# polynomial basis that satisfies the constraints.

N = 501

NS = NormalizedLegendre()
B = [Evaluation(NS, -1, 0); Evaluation(NS, -1, 1); Evaluation(NS, 1, 1)]
ν = size(B, 1)
QS = QuotientSpace(B)
A = Conversion(QS, NS)
Q, R = qr(A[1:N+ν,1:N])

p = Vector{Fun}(undef, N)
for n = 1:N
    en = [zeros(n-1);1.0;zeros(N+ν-n)]
    p[n] = Conversion(NS, Legendre())*Fun(NS, (Q*en)[1:n+ν])
    pad!(p[n], 3n+50)
end

# The mass matrix with entries M_{i,j} = \int_{-1}^1 p_i(x) p_j(x) dx,
# is numerically the identity, as we expect.

M = [innerproduct(p[i], p[j]) for i = 1:N, j = 1:N]
@show norm(M-I) ≤ 4*norm(M)*eps()
@show opnorm(M-I) ≤ 4*sqrt(N)*opnorm(M)*eps()

# Although the polynomials now satisfy deg(p_n) = n+ν instead of deg(p_n) = n,
# the index appears to determine the number of roots of the Galerkin orthogonal
# polynomials. In particular, p_n has exactly n roots in (-1,1).

pl = plot(p[1]; legend = false)
for n in 2:10
    plot!(p[n])
end
plot(pl)

# This need not be the case! If the second constraint is modified to impose
# u'(0) = 0, then the odd-index polynomials have one more root than their index
# and the even-index polynomials have one fewer.

B = [Evaluation(NS, -1, 0); Evaluation(NS, 0, 1); Evaluation(NS, 1, 1)]
QS = QuotientSpace(B)
A = Conversion(QS, NS)
Q, R = qr(A[1:N+ν,1:N])

q = Vector{Fun}(undef, N)
for n = 1:N
    en = [zeros(n-1);1.0;zeros(N+ν-n)]
    q[n] = Conversion(NS, Legendre())*Fun(NS, (Q*en)[1:n+ν])
    pad!(q[n], 3n+50)
end

pl = plot(q[1]; legend = false)
for n in 2:10
    plot!(q[n])
end
plot(pl)

# More information on Galerkin orthogonal polynomials is available in
#
#   P. W. Livermore. Galerkin orthogonal polynomials, J. Comp. Phys., 229:2046–2060, 2010.
#
# The banded QR factorization of the connection coefficients is a linear
# complexity algorithm to represent the Galerkin orthogonal polynomials in H_B
# in terms of the polynomial basis for H. This is described in
#
#   J. L. Aurentz and R. M. Slevinsky. On symmetrizing the ultraspherical spectral method for self-adjoint problems, arXiv:1903.08538, 2019.
#


text:
using MendelBase
using SnpArrays
using IHT
using DataFrames
using Distances
using StatsBase

export MendelIHT

include("MendelIHT_utilities.jl")

"""
This is the wrapper function for the Iterative Hard Thresholding analysis option in Open Mendel.
"""
function MendelIHT(control_file = ""; args...)
    const MENDEL_IHT_VERSION :: VersionNumber = v"0.1.0"
    #
    # Print the logo. Store the initial directory.
    #
    print(" \n \n")
    println("     Welcome to OpenMendel's")
    println("      IHT analysis option")
    println("        version ", MENDEL_IHT_VERSION)
    print(" \n \n")
    println("Reading the data.\n")
    initial_directory = pwd()
    #
    # The user specifies the analysis to perform via a set of keywords.
    # Start the keywords at their default values.
    #
    keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
    #
    # Define some keywords unique to this analysis option.
    #
    keyword["data_type"] = ""
    keyword["predictors"] = ""
    #
    # Process the run-time user-specified keywords that will control the analysis.
    # This will also initialize the random number generator.
    #
    process_keywords!(keyword, control_file, args)
    #
    # Check that the correct analysis option was specified.
    #
    lc_analysis_option = lowercase(keyword["analysis_option"])
    if (lc_analysis_option != "" &&
      lc_analysis_option != "iht")
     throw(ArgumentError(
       "An incorrect analysis option was specified.\n \n"))
    end
    keyword["analysis_option"] = "Iterative Hard Thresholding"
    #
    # Read the genetic data from the external files named in the keywords.
    #
    (pedigree, person, nuclear_family, locus, snpdata,
    locus_frame, phenotype_frame, pedigree_frame, snp_definition_frame) =
    read_external_data_files(keyword)
    #
    # Execute the specified analysis.
    #
    println(" \nAnalyzing the data.\n")
##
    phenotype = convert(Array{Float64,1}, pedigree_frame[:Trait])
    k = keyword["predictors"]
    result = L0_reg(snpdata, phenotype, k)
    return result
##
    # execution_error = iht_gwas(person, snpdata, pedigree_frame, keyword)
    # if execution_error
    #   println(" \n \nERROR: Mendel terminated prematurely!\n")
    # else
    #   println(" \n \nMendel's analysis is finished.\n")
    # end
    # #
    # # Finish up by closing, and thus flushing, any output files.
    # # Return to the initial directory.
    # #
    # close(keyword["output_unit"])
    # cd(initial_directory)
    # return nothing
end #function IterHardThreshold

"""
This function performs IHT on GWAS data.
"""
function L0_reg(
    x        :: SnpData,
    y        :: Vector{Float64},
    k        :: Int;
    v        :: IHTVariable = IHTVariables(x, y, k),
    mask_n   :: Vector{Int} = ones(Int, size(y)),
    tol      :: Float64 = 1e-4,
    max_iter :: Int     = 100,
    max_step :: Int     = 50,
)

    # start timer
    tic()

    # first handle errors
    k        >= 0            || throw(ArgumentError("Value of k must be nonnegative!\n"))
    max_iter >= 0            || throw(ArgumentError("Value of max_iter must be nonnegative!\n"))
    max_step >= 0            || throw(ArgumentError("Value of max_step must be nonnegative!\n"))
    tol      >  eps(Float64) || throw(ArgumentError("Value of global tol must exceed machine precision!\n"))

    # initialize return values
    mm_iter   = 0                 # number of iterations of L0_reg
    mm_time   = 0.0               # compute time *within* L0_reg
    next_loss = oftype(tol,Inf)   # loss function value

    # initialize floats
    current_obj = oftype(tol,Inf) # tracks previous objective function value
    the_norm    = 0.0             # norm(b - b0)
    scaled_norm = 0.0             # the_norm / (norm(b0) + 1)
    μ           = 0.0             # Landweber step size, 0 < tau < 2/rho_max^2

    # initialize integers
    i       = 0                   # used for iterations in loops
    mu_step = 0                   # counts number of backtracking steps for mu

    # initialize booleans
    converged = false             # scaled_norm < tol?

    #convert bitarrays to Float64 genotype matrix, standardize each SNP, and add intercept
    snpmatrix = convert(Array{Float64,2}, x.snpmatrix)
    # snpmatrix = use_A2_as_minor_allele(x.snpmatrix) #to compare results with using PLINK
    snpmatrix = StatsBase.zscore(snpmatrix, 1)
    snpmatrix = [ones(size(snpmatrix, 1)) snpmatrix]

    #
    # Begin IHT calculations
    #
    fill!(v.xb, 0.0) #initialize β = 0 vector, so Xβ = 0
    copy!(v.r, y)    #redisual = y-Xβ = y  CONSIDER BLASCOPY!
    v.r[mask_n .== 0] .= 0 #bit masking? idk why we need this yet

    # calculate the gradient v.df = -X'(y - Xβ) = X'(-1*(Y-Xb)). Future gradient
    # calculations are done in iht!. Note the negative sign will be cancelled afterwards
    # when we do b+ = P_k( b - μ∇f(b)) = P_k( b + μ(-∇f(b))) = P_k( b + μ*v.df)

    # Can we use v.xk instead of snpmatrix?
    At_mul_B!(v.df, snpmatrix, v.r)

    for mm_iter = 1:max_iter
        # save values from previous iterate
        copy!(v.b0, v.b)   # b0 = b   CONSIDER BLASCOPY!
        copy!(v.xb0, v.xb) # Xb0 = Xb  CONSIDER BLASCOPY!
        loss = next_loss

        #calculate the step size μ. Can we use v.xk instead of snpmatrix?
        (μ, μ_step) = iht!(v, snpmatrix, y, k, nstep=max_step, iter=mm_iter)

        # iht! gives us an updated x*b. Use it to recompute residuals and gradient
        v.r .= y .- v.xb
        v.r[mask_n .== 0] .= 0 #bit masking, idk why we need this yet

        At_mul_B!(v.df, snpmatrix, v.r) # v.df = X'(y - Xβ) Can we use v.xk instead of snpmatrix?

        # update loss, objective, gradient, and check objective is not NaN or Inf
        next_loss = sum(abs2, v.r) / 2
        !isnan(next_loss) || throw(error("Objective function is NaN, aborting..."))
        !isinf(next_loss) || throw(error("Objective function is Inf, aborting..."))

        # track convergence
        the_norm    = chebyshev(v.b, v.b0) #max(abs(x - y))
        scaled_norm = the_norm / (norm(v.b0, Inf) + 1)
        converged   = scaled_norm < tol

        if converged
            mm_time = toq()   # stop time

            println("IHT converged in " * string(mm_iter) * " iterations")
            println("It took " * string(mm_time) * " seconds to converge")
            println("The estimated model is stored in 'beta'")
            println("There are " * string(countnz(v.b)) * " non-zero entries of β")
            return IHTResults(mm_time, next_loss, mm_iter, copy(v.b))
        end

        if mm_iter == max_iter
            mm_time = toq() # stop time
            throw(error("Did not converge!!!!! The run time for IHT was " *
                string(mm_time) * "seconds"))
        end
    end
end #function L0_reg

"""
Calculates the IHT step β+ = P_k(β - μ ∇f(β)).
Returns step size (μ), and number of times line search was done (μ_step).

This function updates: b, xb, xk, gk, xgk, idx
"""
function iht!(
    v         :: IHTVariable,
    snpmatrix :: Matrix{Float64},
    y         :: Vector{Float64},
    k         :: Int;
    iter      :: Int = 1,
    nstep     :: Int = 50,
)
    # compute indices of nonzeroes in beta and store them in v.idx (also sets size of v.gk)
    _iht_indices(v, k)

    # fill v.xk, which stores columns of snpmatrix corresponding to non-0's of b
    v.xk[:, :] .= snpmatrix[:, v.idx]

    # fill v.gk, which store only k largest components of gradient (v.df)
    # fill_perm!(v.gk, v.df, v.idx)  # gk = g[v.idx]
    v.gk .= v.df[v.idx]

    # now compute X_k β_k and store result in v.xgk
    A_mul_B!(v.xgk, v.xk, v.gk)

    # warn if xgk only contains zeros
    all(v.xgk .≈ 0.0) && warn("Entire active set has values equal to 0")

    #compute step size and notify if step size too small
    μ = norm(v.gk, 2)^2 / norm(v.xgk, 2)^2
    isfinite(μ) || throw(error("Step size is not finite, is active set all zero?"))
    μ <= eps(typeof(μ)) && warn("Step size $(μ) is below machine precision, algorithm may not converge correctly")

    #Take the gradient step and compute ω. Note in order to compute ω, need β^{m+1} and xβ^{m+1} (eq5)
    _iht_gradstep(v, μ, k)
    ω = compute_ω!(v, snpmatrix) #is snpmatrix required? Or can I just use v.x

    #compute ω and check if μ < ω. If not, do line search by halving μ and checking again.
    μ_step = 0
    for i = 1:nstep
        #exit loop if μ < ω where c = 0.01 for now
        if _iht_backtrack(v, ω, μ); break; end

        #if μ >= ω, step half and warn if μ falls below machine epsilon
        μ /= 2
        μ <= eps(typeof(μ)) && warn("Step size equals zero, algorithm may not converge correctly")

        # recompute gradient step
        copy!(v.b, v.b0)
        _iht_gradstep(v, μ, k)

        # re-compute ω based on xβ^{m+1}
        A_mul_B!(v.xb, snpmatrix, v.b)
        ω = sqeuclidean(v.b, v.b0) / sqeuclidean(v.xb, v.xb0)

        μ_step += 1
    end

    return (μ, μ_step)
end


text:
function readInput(filename::String)::Vector{Int}
	open(filename) do file
		return parse.(Int, split(readline(file), ","))
	end
end

function centroid(vec::Vector{<:Real}; weight_fun::Function)::Int
	min_ = typemax(Int)
	for i in minimum(vec):maximum(vec)
		min_ = min(min_, sum(weight_fun.(abs.(vec .- i))))
	end
	return min_
end

part1(vec) = sum(abs.(vec .- round(Int, mean(vec))))
part2(vec) = centroid(vec, weight_fun=n -> sum(1:n))

vec = readInput("input.txt")
println(part1(vec))
println(part2(vec))


text:
# Model Configuration
## obs likelihood model configuration
Random.seed!(313)
Nhs = [256,256,256,256] # hidden layer dims
dilations = [1,3,5,7]
Nks = 3 .*ones(Int,length(Nhs))         # conv kernel dims
strides = ones(Int,length(Nhs))
strides[end] = subsample                # output subsampling
dropout = zeros(length(Nhs))
fout = identity                         # function in last layer

# training opts
λ1 = 1f-5      # l1 output regularization
lr = 1f-3      # learning rate
Nb = 16        # batch size
epochs_cur = 5
epochs     = 15
opt = ADAM(lr)
curriculum_training = true


text:
# Work around differences in julia 0.2 and 0.3 set constructors.
function set(T::Type, itr)
    S = Set{T}()
    union!(S, itr)
end

function set(itr)
    set(Any, itr)
end

# Is this usable data?
function isconcrete{T<:Number}(x::T)
    !isna(x) && isfinite(x)
end

function isconcrete(x)
    !isna(x)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    n = 0
    for x in xs
        if isconcrete(x)
            n += 1
        end
    end
    n
end

function concrete_minimum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_min = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_max = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function nonzero_length(xs)
    n = 0
    for x in xs
        if x != 0
            n += 1
        end
    end
    n
end

function default_formatter(xs)
    [string(x) for x in xs]
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit{T}(a::T, b::T)
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!{T}(a::T, b::T)
    for field in T.names
        aval = getfield(a, field)
        bval = getfield(b, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if aval === nothing || aval === string || aval == default_formatter
            setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = is(u, nothing)
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has{T,N}(xs::AbstractArray{T,N}, y::T)
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

Maybe(T) = Union(T, Nothing)

function lerp(x::Float64, a::Float64, b::Float64)
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::String)
    s = replace(s, r"<[^>]*>", "")
    s = replace(s, r"\s", "_")
    s
end

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable{T <: FloatingPoint}(x::T) = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    @assert n > 1

    step = (b - a) / (n - 1)
    xs = Array(typeof(a + step), n)
    for i in 1:n
        xs[i] = a + (i-1) * step
    end

    df = DataFrame(xs, map(f, xs))
    # NOTE: 'colnames!' is the older deprecated name. 'names!' was also defined
    # but threw an error.
    try
        names!(df, [:x, :f_x])
    catch
        colnames!(df, ["x", "f_x"])
    end
    df
end


text:
###
# Tweedie with Covariates
###
module tw
using JuMP, Gurobi, MLBase

###
#Solves QP for optimal weights
#assuming bound on h : norm(h)^2 <= hmax dualized to lam1
#assuming bound on (a,b) : a^2 + b^2 <= bmax^2 dualized to lam 2
#has been dualized into parameter lambda.
#k is a two parameter kernel: k(y1, y2)
###
#quadratic portion of obj.
function form_omega(ys, ws, k, k1)
	const n = length(ys)
	out = zeros(2n + 2, 2n + 2)
	for ix = 1:n
		KKei = vcat(1, ys[ix],
					[ k(ys[jx], ys[ix]) for jx = 1:n],
					[k1(ys[jx], ys[ix]) for jx = 1:n]
					)
		out += ws[ix] * KKei * KKei'
	end
	return out
end

#Hnorm portion of the regularizer
function form_graham(ys, k, k1, k12, offset)
	const n = length(ys)
	K   = [  k(ys[ix], ys[jx]) for ix = 1:n, jx=1:n]
	K1  = [ k1(ys[ix], ys[jx]) for ix = 1:n, jx=1:n]
	K12 = [k12(ys[ix], ys[jx]) for ix = 1:n, jx=1:n]
	out = 	[K  K1';
	 		 K1  K12]
	out += offset * eye(2n)
	return out
end

function form_lincoeff(ys, ws, k1, k12)
	const n = length(ys)
	K1_T =  [ k1(ys[jx], ys[ix]) for ix = 1:n, jx=1:n]
	K12_T = [k12(ys[jx], ys[ix]) for ix = 1:n, jx=1:n]
	return 2 * vcat(0, 1, K1_T * ws, K12_T * ws)
end

#col_vals is the vector (b, a, alphas, betas)
function eval_score(col_vals, ys, new_y, k, k1)
	const n = length(ys)
	out = col_vals[1] + col_vals[2] * new_y
	out += sum(col_vals[2 + ix] * k(ys[ix], new_y) for ix = 1:n)
	out += sum(col_vals[2 + n + ix] * k1(ys[ix], new_y) for ix = 1:n)
	return out
end

#Call to Gurobi to do the lifting in the lazy way.
#A smarter solution would combine proj gradient descent with a sweeping alg.
#lam 1 is hnorm regularizer norm_H(h)^2
#lam 2 is two norm regularize a^2 + b^2
function fit_cond_score(ys, ws, k, k1, k12, lam1::Float64, lam2::Float64;
					force_inc = true, output=false, TOL=1e-5)
	const n = length(ys)
	@assert n >=2 "Require at least 2 data points"

	omega = form_omega(ys, ws, k, k1)
	b = form_lincoeff(ys, ws, k1, k12)
	graham = form_graham(ys, k, k1, k12, TOL)

	m = Model(solver=GurobiSolver(OutputFlag=Int(output)))
	@variable(m, alphas[1:n])
	@variable(m, betas[1:n])
	@variable(m, null_coef[1:2])  #null[1] + null[2] * Y
	col_var = vcat(null_coef, alphas, betas)

	#Monotonicity constraints
	if force_inc
		indx = sortperm(ys)
		for ix = 2:n
			vs_ix = null_coef[2] * ys[indx[ix]] +
					sum(alphas[jx] * k(ys[jx], ys[indx[ix]]) for jx = 1:n) +
					sum(betas[jx] * k1(ys[jx], ys[indx[ix]]) for jx = 1:n)
			vs_p = null_coef[2] * ys[indx[ix - 1]] +
					sum(alphas[jx] * k(ys[jx], ys[indx[ix - 1]]) for jx = 1:n) +
					sum(betas[jx] * k1(ys[jx], ys[indx[ix - 1]]) for jx = 1:n)
			@constraint(m, vs_p <= vs_ix)
		end
	end

	#objective
	obj = col_var' * omega * col_var + b' * col_var
	@objective(m, Min, obj +
						lam1 * col_var[3:end]' * graham * col_var[3:end] +
						lam2 * null_coef' * null_coef)
	status = solve(m)
	#@assert status== :Optimal "Error solving: $status"
	#println(status)

	col_vals = getvalue(col_var)
	vals = [col_vals[1] + col_vals[2] * ys[ix] +
			sum(col_vals[2 + jx] * k(ys[jx], ys[ix]) for jx = 1:n) +
			sum(col_vals[2 + n + jx] * k1(ys[jx], ys[ix]) for jx = 1:n)
			for ix = 1:n]

	return vals, col_vals
end

#Neglect monotonicity and solve a linear system directly.
#lam1 is hnorm regularizer
#lam2 is null_space regularizer
function fit_cond_score_ls(ys, ws, k, k1, k12, lam1::Float64, lam2::Float64; TOL=1e-5)
	const n = length(ys)

	Q = form_omega(ys, ws, k, k1)
	b = form_lincoeff(ys, ws, k1, k12)
	graham_h = form_graham(ys, k, k1, k12, TOL)

	Q += [lam2 * ones(2, 2)  	zeros(2, 2n);
		  zeros(2n, 2) 			lam1 * graham_h]

	#min x' Q x + b'x ->  Qx  = - .5 * b
	col_vals = \(Q, -.5 * b)

	vals = [col_vals[1] + col_vals[2] * ys[ix] +
			sum(col_vals[2 + jx] * k(ys[jx], ys[ix]) for jx = 1:n) +
			sum(col_vals[2 + n + jx] * k1(ys[jx], ys[ix]) for jx = 1:n)
			for ix = 1:n]

	return vals, col_vals
end

#kw is now the kernel for the weights
#kw: R -> R
function fit_NW_ws(xs, ys, new_x, kw)
	const n = length(xs)
	#use the standard scaling for conenience.
	const h = n^(-1/6)
	ws = [kw(norm(xs[ix] - new_x)/h) for ix = 1:n] #extra h scaling drops out
	ws /= sum(ws)
	return ws
end

#####
#A tuned model consists of
#	training data (xs, ys)
#	a procedure that given new_x -> weights
#	Kernel parameters for the RKHS
#	lambdas

# Evaluation entails
#	computing the new weights
#	solving the estimation problem for conditional score at new pt
#	outputting estimate
######

####
#Assumes
#   weight procedure is fixed: fit_ws(xs, ys, new_x) -> ws
#	Kernel parameters are fixed
#	lam2 fixed, small
#Only Tunes lam1
########
function tune_lam1(xs, ys, fit_ws, k, k1, k12, lam1_grid, lam2::Float64, TOL=1e-5)
	const num_lams = length(lam1_grid)
	const n = length(ys)
	cv_vals = zeros(2, num_lams) #mean and std
	for (ilam, lam) in enumerate(lam1_grid)
		#Nothing to learn in the estimation phase!
		function est_fun(indices)
			return xs[indices], ys[indices]
		end

		function eval_fun(fit, test_indices)
			(xs_train, ys_train) = fit
			#compute MSE on the trained values
			MSE = 0.
			for indx in test_indices
				#VG Check how test set is passed
				new_x, new_y = xs[indx], ys[indx]
				ws = fit_ws(xs_train, ys_train, new_x)
				vals, col_vals = fit_cond_score(ys_train, ws, k, k1, k12, lam, lam2, TOL=TOL)
				MSE += (new_y - (new_y + new_x[1]^2 * eval_score(col_vals, ys_train, new_y, k, k1)))^2
			end
			println(MSE)
			return MSE
		end

		#now evaluate the cross_val.
		folds_results = cross_validate(est_fun, eval_fun, n, Kfold(n, 5))
		cv_vals[1, ilam] = mean(folds_results)
		cv_vals[2, ilam] = std(folds_results)
	end

	#For now, return everything
	cv_vals
end

#Some helper functions
k_gauss(x, y, c)   =  exp(-c * (x - y)^2 )
k1_gauss(x, y, c)  = -2c * k_gauss(x, y, c) * (x-y)
k12_gauss(x, y, c) =  k_gauss(x, y, c) * (2c - 4c^2 *(x-y)^2)

end #ends module


text:
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1,T2,W} = promote_type(T1,T2)
Base.@pure function _ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1, T2, W}
    T_canon = promote_type(T1,T2)
    if pick_vector_width(T_canon) < W
        T1
    else
        T_canon
    end
end
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: Integer, T2 <: Integer,W} = _ff_promote_rule(T1,T2,Val{W}())
ff_promote_rule(::Type{T1}, ::Type{T2}, ::Val{W}) where {T1 <: FloatingTypes, T2<:FloatingTypes,W} = _ff_promote_rule(T1,T2,Val{W}())

Base.promote_rule(::Type{V}, ::Type{T2}) where {W,T1,T2<:NativeTypes,V<:AbstractSIMDVector{W,T1}} = Vec{W,ff_promote_rule(T1,T2,Val{W}())}
Base.promote_rule(::Type{V1}, ::Type{V2}) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2}} = Vec{W,promote_type(T1,T2)}

_assemble_vec_unroll(::Val{N}, ::Type{V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = VecUnroll{N,W,T,V}
Base.promote_rule(::Type{VecUnroll{N,W,T1,V}}, ::Type{T2}) where {N,W,T1,V,T2<:NativeTypes} = _assemble_vec_unroll(Val{N}(), promote_type(V,T2))
Base.promote_rule(::Type{<:VecUnroll{N,W,T,V1}}, ::Type{V2}) where {N,W,T,V1,V2<:AbstractSIMDVector{W}} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))
Base.promote_rule(::Type{VecUnroll{N,W,T1,V1}}, ::Type{VecUnroll{N,W,T2,V2}}) where {N,W,T1,T2,V1,V2} = _assemble_vec_unroll(Val{N}(), promote_type(V1,V2))

Base.promote_rule(::Type{Bit}, ::Type{T}) where {T <: Number} = T

# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{T2}) where {M,N,W,T1,T2<:NativeTypes} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{Vec{W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecUnroll{M,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecTile{M,N,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}


text:
using BPM

using PyPlot
rc("figure", figsize=(3.5, 3.5))
rc("font", size=10.0)
rc("lines", linewidth=1.5)
rc("lines", markersize=3.0)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
rc("figure.subplot", left=0.15, bottom=0.05, top=0.97, right=.93)
rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
import Atmosphere
# function runme()
meshgrid(x,y) = (repmat(x',length(y),1),repmat(y,1,length(x)))

# INITIALIZE
nturbines = 8
x_test = ones(nturbines)*0.0 # x-locations of turbines (m)
y_test = ones(nturbines)*0.0 # y-locations of turbines (m)
obs_test = [0., 20.0, 1.5] # x-, y-, and z-location of the observer (m)
winddir_test = 0. # wind direction (deg)
rpm_test = ones(nturbines)*5000.0 # rotation rate of the tubrines (rpm)
windvel_test = ones(nturbines)*1.0 # wind velocity (m/s)
h_test = 15.24 # height of the turbine hub (m)
noise_corr = 0.8697933840957954 # correction factor for noise

alpha =  [5.10514, 5.10514, 6.8703, 6.96041, 6.07495, 4.58821, 3.20019, 3.20019]#[10.0,5.0,1.0]
R_station = [0.0778218, 0.177878, 0.277935, 0.377991, 0.478048, 0.578104, 0.678161, 0.778218]#[0.0,0.125,0.25]*2
p_chord = [0.155659, 0.19989, 0.210314, 0.172972, 0.125581, 0.0843676, 0.0460056, 0.0100891]#[.01,.05,.01]
AR = R_station[end]/mean(p_chord) #TODO: actual AR #10.0 #TODO: actual AR
psi = 14.0 # solid angle (deg)

area = pi*R_station[end]^2 * nturbines

Rtip_orig = R_station[end]

rho,mu,a = Atmosphere.atmospherefit(0.0) #Use average altitude for each evaluation

blades = 2

tipspeed = collect(linspace(0.1,0.8,10))
omega_correct_speed = tipspeed*a/R_station[end]
RPM_correct_speed = omega_correct_speed * 60 / (2*pi)
numprops = [2,4,8,16,32]#collect(linspace(2,32,5))
nBlades = collect(linspace(2,5,length(numprops)))
numprops2 = [2,4,8,16,32]
db_test1 = zeros(length(numprops),length(tipspeed))
checkrpm = zeros(length(numprops),length(tipspeed))

levels = [60,65,70,72.5]

X,Y=meshgrid(tipspeed,numprops)

for j = 1:length(tipspeed)
    for i = 1:length(numprops)
        nturbines = 8#round(Int,numprops[i])
        Rtip_scale = sqrt(area/nturbines/pi)/Rtip_orig
        blades = 4#round(Int,nBlades[i])
        R_station = [0.0778218, 0.177878, 0.277935, 0.377991, 0.478048, 0.578104, 0.678161, 0.778218]*Rtip_scale
        p_chord = [0.155659, 0.19989, 0.210314, 0.172972, 0.125581, 0.0843676, 0.0460056, 0.0100891]*Rtip_scale

        omega_correct_speed = tipspeed[j]*a/R_station[end]
        RPM_correct_speed = omega_correct_speed * 60 / (2*pi)

        x_test = ones(nturbines)*0.0 # x-locations of turbines (m)
        y_test = ones(nturbines)*0.0 # y-locations of turbines (m)
        windvel_test = ones(nturbines)*1.0 # wind velocity (m/s)
        AR = R_station[end]/mean(p_chord) #TODO: actual AR #10.0 #TODO: actual AR

        rpm_test = ones(nturbines)*RPM_correct_speed # rotation rate of the tubrines (rpm)
        checkrpm[i,j] = RPM_correct_speed
        db_test1[i,j] = BPM.turbinepos(x_test, y_test, obs_test, winddir_test, windvel_test, rpm_test, blades, h_test, R_station, p_chord, p_chord*0.25, alpha, mu/rho, a, psi, AR, noise_corr)
        println(mean(p_chord))
    end
    # println(tipspeed[j])

end

println(maximum(db_test1))

close("all")
figname = "Noise_Tipspeed_Contour_scaledchord"
fig = figure(figname)
CS = PyPlot.contourf(Y,X,db_test1,interpolation = "spline16",origin = "lower",cmap = PyPlot.cm[:Blues],levels = linspace(28,maximum(db_test1),1000))#linspace(28,77,1000))
PyPlot.colorbar(orientation = "horizontal",extend = "both",label = "SPL (dBa)",pad = 0.2,aspect = 20,ticks = [30,40,50,60,70])
CS2 = PyPlot.contour(Y,X,db_test1,colors = "white",levels = levels)#linspace(28,77,1000))
clabel( CS2,fmt = "%1.1f")
xticks(numprops2)
PyPlot.xlabel("Propellers (#)")
PyPlot.ylabel("Tipspeed (mach)")
# xlim([Inlet_pos,Outlet_pos])
# ylim([Bot_pos,Top_pos])
PyPlot.savefig("../MDAO_paper/figures/$figname.png",transparent = true,dpi = 600)

db_test2 = zeros(length(numprops),length(tipspeed))

for j = 1:length(tipspeed)
    for i = 1:length(numprops)
        nturbines = round(Int,numprops[i])
        Rtip_scale = sqrt(area/nturbines/pi)/Rtip_orig
        blades = round(Int,nBlades[i])
        R_station = [0.0778218, 0.177878, 0.277935, 0.377991, 0.478048, 0.578104, 0.678161, 0.778218]*Rtip_scale
        p_chord = [0.155659, 0.19989, 0.210314, 0.172972, 0.125581, 0.0843676, 0.0460056, 0.0100891]#*Rtip_scale

        omega_correct_speed = tipspeed[j]*a/R_station[end]
        RPM_correct_speed = omega_correct_speed * 60 / (2*pi)

        x_test = ones(nturbines)*0.0 # x-locations of turbines (m)
        y_test = ones(nturbines)*0.0 # y-locations of turbines (m)
        windvel_test = ones(nturbines)*1.0 # wind velocity (m/s)
        AR = R_station[end]/mean(p_chord) #TODO: actual AR #10.0 #TODO: actual AR

        rpm_test = ones(nturbines)*RPM_correct_speed # rotation rate of the tubrines (rpm)
        checkrpm[i,j] = RPM_correct_speed
        db_test2[i,j] = BPM.turbinepos(x_test, y_test, obs_test, winddir_test, windvel_test, rpm_test, blades, h_test, R_station, p_chord, p_chord*0.25, alpha, mu/rho, a, psi, AR, noise_corr)

    end
    println(tipspeed[j])
end

figname = "Noise_Tipspeed_Contour_constchord"
fig = figure(figname)
CS = PyPlot.contourf(Y,X,db_test2,interpolation = "spline16",origin = "lower",cmap = PyPlot.cm[:Blues],levels = linspace(28,maximum(db_test1),1000))#linspace(28,77,1000))
PyPlot.colorbar(orientation = "horizontal",extend = "both",label = "SPL (dBa)",pad = 0.2,aspect = 20,ticks = [30,40,50,60,70])
CS2 = PyPlot.contour(Y,X,db_test2,colors = "white",levels = levels)#linspace(28,77,1000))
clabel( CS2,fmt = "%1.1f")
xticks(numprops2)
PyPlot.xlabel("Propellers (#)")
PyPlot.ylabel("Tipspeed (mach)")
# xlim([Inlet_pos,Outlet_pos])
# ylim([Bot_pos,Top_pos])
PyPlot.savefig("../MDAO_paper/figures/$figname.png",transparent = true,dpi = 600)

# figname = "CheckRPM"
# figure(figname)
# CS = PyPlot.contourf(Y,X,checkrpm,interpolation = "spline16",origin = "lower",cmap = PyPlot.cm[:Blues],levels = linspace(minimum(checkrpm),maximum(checkrpm),1000))
# PyPlot.colorbar(orientation = "horizontal",extend = "both",label = "SPL (dba)",pad = 0.2,aspect = 20,ticks = round.(linspace(minimum(checkrpm),maximum(checkrpm),8),1))
# xticks(numprops2)
# PyPlot.xlabel("Propellers (#)")
# PyPlot.ylabel("Tipspeed (mach)")
# # xlim([Inlet_pos,Outlet_pos])
# # ylim([Bot_pos,Top_pos])
# # PyPlot.savefig("../MDAO_paper/figures/$figname.pdf",transparent = true)


text:
using Base.Test
# using Lint

# lintpkg("AutomotiveDrivingModels")

using AutomotiveDrivingModels

include(Pkg.dir("AutomotiveDrivingModels", "test", "core", "test_core.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "test", "test_actions.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "test_driver_models.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "test", "test_roadway_generation.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "test_minkowski.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "test", "test_features.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "test", "test_feature_extractors.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "test", "evaluation", "test_evaluation.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "test", "behaviors", "test_behaviors.jl"))