text:
# auto-generated tests from julia-repl docstrings
using Test, ModuleElts
function mytest(a::String,b::String)
  omit=a[end]==';'
  a=replace(a,"\\\\"=>"\\")
  a=repr(MIME("text/plain"),eval(Meta.parse(a)),context=:limit=>true)
  if omit a="nothing" end
  a=replace(a,r" *(\n|$)"s=>s"\1")
  a=replace(a,r"\n$"s=>"")
  b=replace(b,r" *(\n|$)"s=>s"\1")
  b=replace(b,r"\n$"s=>"")
  i=1
  while i<=lastindex(a) && i<=lastindex(b) && a[i]==b[i]
    i=nextind(a,i)
  end
  if a!=b print("exec=$(repr(a[i:end]))\nmanl=$(repr(b[i:end]))\n") end
  a==b
end
@testset "ModuleElts" begin
@test mytest("a=ModuleElt(:xy=>1,:yx=>-1)",":xy-:yx")
@test mytest("repr(a)","\"ModuleElt([:xy => 1, :yx => -1])\"")
@test mytest("ModuleElt([:xy=>1//2,:yx=>-1])","(1//2):xy+(-1//1):yx")
@test mytest("a-a","0")
@test mytest("a*99","99:xy-99:yx")
@test mytest("a//2","(1//2):xy+(-1//2):yx")
@test mytest("a/2","0.5:xy-0.5:yx")
@test mytest("a+ModuleElt(:yx=>1)",":xy")
@test mytest("a[:xy]","1")
@test mytest("a[:xx]","0")
@test mytest("haskey(a,:xx)","false")
@test mytest("first(a)",":xy => 1")
@test mytest("collect(a)","2-element Vector{Pair{Symbol, Int64}}:\n :xy => 1\n :yx => -1")
@test mytest("collect(keys(a))","2-element Vector{Symbol}:\n :xy\n :yx")
@test mytest("collect(values(a))","2-element Vector{Int64}:\n  1\n -1")
@test mytest("length(a)","2")
@test mytest("eltype(a)","Pair{Symbol, Int64}")
@test mytest("a=ModuleElt(:yy=>1, :yx=>2, :xy=>3, :yy=>-1;check=false)",":yy+2:yx+3:xy-:yy")
@test mytest("a=ModuleElt(:yy=>1, :yx=>2, :xy=>3, :yy=>-1)","3:xy+2:yx")
@test mytest("a+ModuleElt([:z=>1.0])","3.0:xy+2.0:yx+1.0:z")
end


text:
module RandomFunctions

include("collatz.jl")
export collatz_steps, collatz_steps_02, max_stop_time, max_stop_time_02, max_stop_time_03

include("edit_distance.jl")
export edit_distance

include("sieve.jl")
export sieve_01,sieve_02,sieve_03

end

ext: jl


text:
module ExportNURBS

using BasicBSpline
using Luxor
using IntervalSets
using Colors
using ColorVectorSpace

export save_png, save_svg, save_pov

include("_common.jl")
include("_luxor.jl")
include("_povray.jl")

end # module


text:
# A mini balance law for computing vorticity using the DG kernels.
#

using StaticArrays

using ..BalanceLaws
using ..VariableTemplates
using ..MPIStateArrays

import ..BalanceLaws:
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    eq_tends,
    flux,
    source,
    wavespeed,
    boundary_conditions,
    boundary_state!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    transform_post_gradient_laplacian!,
    init_state_auxiliary!,
    init_state_prognostic!,
    update_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

# A mini balance law that is used to take the gradient of u and v
# to obtain vorticity.
struct VorticityModel <: BalanceLaw end

# output vorticity vector
vars_state(::VorticityModel, ::Prognostic, FT) = @vars(Ω_bl::SVector{3, FT})
# input velocity vector
vars_state(::VorticityModel, ::Auxiliary, FT) = @vars(u::SVector{3, FT})
vars_state(::VorticityModel, ::Gradient, FT) = @vars()
vars_state(::VorticityModel, ::GradientFlux, FT) = @vars()

# required for each balance law
function init_state_auxiliary!(
    m::VorticityModel,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
) end
function init_state_prognostic!(
    ::VorticityModel,
    state::Vars,
    aux::Vars,
    localgeo,
    t,
) end
function nodal_update_auxiliary_state!(
    m::VorticityModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end;
function flux_first_order!(
    ::VorticityModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    u = aux.u
    @inbounds begin
        flux.Ω_bl = @SMatrix [
            0 u[3] -u[2]
            -u[3] 0 u[1]
            u[2] -u[1] 0
        ]
    end
end
function compute_gradient_argument!(
    ::VorticityModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
flux_second_order!(::VorticityModel, _...) = nothing
source!(::VorticityModel, _...) = nothing
boundary_conditions(::VorticityModel) = ntuple(i -> nothing, 6)
boundary_state!(nf, ::Nothing, ::VorticityModel, _...) = nothing


text:
using Printf

using Oceananigans
using Oceananigans.Advection
using Oceananigans.Fields
using Oceananigans.Utils
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.AbstractOperations
using Oceananigans.Utils

# Some physical constants
const f₀ = 1e-4  # Coriolis parameter [s⁻¹]
const ρ₀ = 1027  # Density of seawater [kg/m³]
const cₚ = 4000  # Specific heat capacity of seawater at constant pressure [J/(kg·K)]

# Set up grid
N, L = 256, 100
topology = (Periodic, Periodic, Bounded)
grid = RegularCartesianGrid(topology=topology, size=(N, N, N), x=(0, L), y=(0, L), z=(-L, 0))

# Set up wind stress
τ = 0.02  # Wind stress [N/m^2]
u_top_bc = FluxBoundaryCondition(τ / ρ₀)

Q = 0  # Upward surface heat flux [W/m²]
T_top_bc = FluxBoundaryCondition(Q / (ρ₀ * cₚ))

∂T∂z = 0.01  # Initial temperature stratification [°C/m]
T_bot_bc = GradientBoundaryCondition(∂T∂z)

u_bcs = UVelocityBoundaryConditions(grid, top=u_top_bc)
T_bcs = TracerBoundaryConditions(grid, top=T_top_bc, bottom=T_bot_bc)

# Set up model
model = IncompressibleModel(
            architecture = GPU(),
                    grid = grid,
             timestepper = :RungeKutta3,
               advection = WENO5(),
                coriolis = FPlane(f=f₀),
                 tracers = :T,
                buoyancy = SeawaterBuoyancy(constant_salinity=true),
                 closure = AnisotropicMinimumDissipation(),
     boundary_conditions = (u=u_bcs, T=T_bcs)
)

# Set initial conditions
ε(σ) = σ * randn() # Gaussian noise with mean 0 and variance μ
T₀(x, y, z) = 20 + ∂T∂z * z + ε(1e-10) * exp(4z/L)
set!(model, T=T₀)

# Set up adaptive time stepping
wizard = TimeStepWizard(cfl=0.3, Δt=3.0, max_change=1.2, max_Δt=30.0)
cfl = AdvectiveCFL(wizard)

function print_simulation_stats(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time
    cfl = simulation.parameters.cfl

    progress = 100 * (t / simulation.stop_time)

    u_max = maximum(abs, model.velocities.u.data.parent)
    v_max = maximum(abs, model.velocities.v.data.parent)
    w_max = maximum(abs, model.velocities.w.data.parent)
    T_min = minimum(abs, model.tracers.T.data.parent)
    T_max = maximum(abs, model.tracers.T.data.parent)

    @info @sprintf("[%05.2f%%] i: %d, t: %s, umax: (%6.3g, %6.3g, %6.3g) m/s, T: (min=%6.4f, max=%6.4f), CFL: %6.4g, next Δt: %s",
                   progress, i, prettytime(t), u_max, v_max, w_max, T_min, T_max, cfl(model), prettytime(simulation.Δt.Δt))
end

# Set up simulation
simulation = Simulation(model, Δt=wizard, stop_time=4days, progress=print_simulation_stats,
                        iteration_interval=50, parameters=(cfl=cfl,))

#####
##### Set up output writers
#####

# Code credit: https://discourse.julialang.org/t/collecting-all-output-from-shell-commands/15592
function execute(cmd::Cmd)
    out, err = Pipe(), Pipe()

    process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
    close(out.in)
    close(err.in)

    return (stdout = out |> read |> String, stderr = err |> read |> String, code = process.exitcode)
end

global_attributes = Dict(
    "ClimateSurrogates.jl git commit SHA1" => execute(`git rev-parse HEAD`).stdout |> strip,
    "Reference density" => "$ρ₀ kg/m³",
    "Specific_heat_capacity" => "$cₚ J/kg/K",
    "Coriolis parameter" => "$f₀ s⁻¹",
    "Gravitational acceleration" => "$(model.buoyancy.gravitational_acceleration) m/s²",
    "Heat flux" => "$Q W/m²",
    "Initial stratification" => "$∂T∂z °C/m",
    "closure" => "$(model.closure)",
    "coriolis" => "$(model.coriolis)"
)

(u, v, w), T = model.velocities, model.tracers.T

L⁻²∫u_dxdy = AveragedField(u, dims=(1, 2))
L⁻²∫v_dxdy = AveragedField(v, dims=(1, 2))
L⁻²∫T_dxdy = AveragedField(T, dims=(1, 2))

L⁻²∫uu_dxdy = AveragedField(u * u, dims=(1, 2))
L⁻²∫vv_dxdy = AveragedField(v * v, dims=(1, 2))
L⁻²∫ww_dxdy = AveragedField(w * w, dims=(1, 2))
L⁻²∫uv_dxdy = AveragedField(u * v, dims=(1, 2))
L⁻²∫uw_dxdy = AveragedField(u * w, dims=(1, 2))
L⁻²∫vw_dxdy = AveragedField(v * w, dims=(1, 2))

L⁻²∫uT_dxdy = AveragedField(u * T, dims=(1, 2))
L⁻²∫vT_dxdy = AveragedField(v * T, dims=(1, 2))
L⁻²∫wT_dxdy = AveragedField(w * T, dims=(1, 2))

profiles = Dict(
    "u"  => L⁻²∫u_dxdy,
    "v"  => L⁻²∫v_dxdy,
    "T"  => L⁻²∫T_dxdy,
    "uu" => L⁻²∫uu_dxdy,
    "vv" => L⁻²∫vv_dxdy,
    "ww" => L⁻²∫ww_dxdy,
    "uv" => L⁻²∫uv_dxdy,
    "uw" => L⁻²∫uw_dxdy,
    "vw" => L⁻²∫vw_dxdy,
    "uT" => L⁻²∫uT_dxdy,
    "vT" => L⁻²∫vT_dxdy,
    "wT" => L⁻²∫wT_dxdy
)

nc_filepath = "wind_mixing_horizontal_averages_$(τ)Nm2.nc"
simulation.output_writers[:profiles] =
    NetCDFOutputWriter(model, profiles, filepath=nc_filepath,
                       global_attributes=global_attributes,
                       schedule=TimeInterval(10minutes), verbose=true)

run!(simulation)


text:
module ExtraTests

module UtilityTests
using Base.Test
using ExpressionPatterns
using ExpressionPatterns.Helper

cd("../src")
@test isfile(ExpressionPatterns.path("Analyzer"))
cd("../test")

@test clean_code(quote x end).args == [:x]

end

module PrintingTests
using Base.Test
using ExpressionPatterns
using ExpressionPatterns.Dispatch
using ExpressionPatterns.Dispatch.Reflection

@metafunction f(x,y) = ()
@metafunction g()[g] = ()
@metafunction h()    = ()
@metafunction h(x)   = ()

totree(ex) = ExpressionPatterns.Analyzer.analyze(ex, PrintingTests)
showcpt(io, ex) = show(IOContext(io, :compact=>true), ex)

@test sprint(show, totree(:(10)))     == "pattern`10`"
@test sprint(show, totree(:(x+y)))    == "pattern`x + y`"
@test sprint(show, totree(:(*{x})))   == "pattern`(:*){x}`"
@test sprint(show, totree(:(:?{x})))  == "pattern`(:?){x}`"

@test sprint(show, totree(:(:T{x,Int64})))   == "pattern`(:type){x, Int64}`"
@test sprint(show, totree(:(:L{x})))         == "pattern`(:equals){x}`"
@test sprint(show, totree(:(:B{+})))         == "pattern`(:binding){+}`"
@test sprint(show, totree(:(:P{x, iseven}))) == "pattern`(:predicate){x, iseven}`"

@test sprint(show, (@metamethods f)[1])   == "[---]<x y>"
@test sprint(show, (@metamethods g)[1])   == "[g]<>"

import ExpressionPatterns.Dispatch.Reflection.metafuntable
@test sprint(showcpt, metafuntable(:f, PrintingTests)) == "metafunction f"
@test sprint(show,    metafuntable(:f, PrintingTests)) == "metafunction f with 1 method:\n [---]<x y>\n"
@test sprint(show,    metafuntable(:h, PrintingTests)) == "metafunction h with 2 methods:\n [---]<>\n [---]<x>\n"

end

module MakeDocsTest

cd("../src/Docs")
include("../src/Docs/makedocs.jl")
cd("../../test")

end

module TreeComparison
using Base.Test
using ExpressionPatterns.Analyzer
using ExpressionPatterns.Matching.Comparison

totree(x) = analyze(x, TreeComparison)

@test totree(:(a,*{x}))   ⊆ totree(:(*{x},))
@test totree(:(*{x},))    ⊇ totree(:(a,*{x}))
@test totree(:(*{x},1))   ⊇ totree(:(*{x},2,1))
@test totree(:(*{x},2,1)) ⊆ totree(:(*{x},1))

end

end


text:
# creating a module
module Prism_Mod

# exporting a function for aliasing
# export only works with `using` keyword
export volume

struct Prism
    length::Real
    width::Real
    height::Real

    function Prism()
        new(1, 1, 1)
    end

    function Prism(l::Real, w::Real, h::Real)

        if l < 0 || w < 0 || h < 0
            error("Can't have negative values for lengths.")
        elseif w < l
            error("Can't have shorter width than length.")
        else
            new(l, w, h)
        end
    end
end

function volume(p::Prism)
    v = p.length * p.width * p.height
end

function surface_area(p::Prism)
    s1 = p.length * p.width
    s2 = p.length * p.height
    s3 = p.width * p.height

    sa = 2 * (s1 + s2 + s3)
end

end # end of module

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
import SolverTraces: base_exp

function test_base_exp(v, b, e)
    bv,be = base_exp(v)
    @test bv ≈ b
    @test be == e
end

@testset "Tolerance" begin
    test_base_exp(0.0, 0.0, 0)
    test_base_exp(3.5e2, 3.5, 2)
    test_base_exp(-3.5e2, -3.5, 2)
    test_base_exp(0.1, 1.0, -1)
    test_base_exp(-0.01π, -π, -2)

    t = Tolerance(1e-3)
    t.current=1.0
    test_esc_str_equal(format(t, 1), "\e[38;2;255;0;0m1.00×10⁰  \e[0m(10⁻³)")

    t.current=1e-1
    test_esc_str_equal(format(t, 1), "\e[38;2;255;170;0m1.00×10⁻¹ \e[0m(10⁻³)")

    t.current=1e-2
    test_esc_str_equal(format(t, 1), "\e[38;2;170;255;0m1.00×10⁻² \e[0m(10⁻³)")

    t.current=1e-3
    test_esc_str_equal(format(t, 1), "\e[38;2;0;255;0m1.00×10⁻³ \e[0m(10⁻³)")

    t = Tolerance(1e-3, print_target=false)
    t.current=1.0
    test_esc_str_equal(format(t, 1), "\e[38;2;255;0;0m1.00×10⁰  \e[0m")

    t.current=1e-1
    test_esc_str_equal(format(t, 1), "\e[38;2;255;170;0m1.00×10⁻¹ \e[0m")

    t.current=1e-2
    test_esc_str_equal(format(t, 1), "\e[38;2;170;255;0m1.00×10⁻² \e[0m")

    t.current=1e-3
    test_esc_str_equal(format(t, 1), "\e[38;2;0;255;0m1.00×10⁻³ \e[0m")

    t.current=Inf
    test_esc_str_equal(format(t, 1), "\e[38;2;255;0;0mInf       \e[0m")

    t.current=-Inf
    test_esc_str_equal(format(t, 1), "\e[38;2;255;0;0m-Inf      \e[0m")

    t.current=NaN
    test_esc_str_equal(format(t, 1), "\e[38;2;255;0;0mNaN       \e[0m")
end


text:
function [ varm,var,varp ] = asselin( varm,var,varp,asscoef,nx,nz )
%UNTITLED10 Summary of this function goes here
%   Detailed explanation goes here

% Apply Asselin filtering for 2dt filtering
for k = 2 : nz-1
  for i = 2 : nx-1

var[i,k] = var[i,k] + asscoef*(varp[i,k]-2*var[i,k]+varm[i,k]);

  end
end

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
@inline S_inner_to_outer(S_in, u, xi, phi0) =
    1/u + xi - u * phi0*phi0/3 * (1 - u * xi) + u*u*u * S_in

@inline S_u_inner_to_outer(S_u_in, S_in, u, xi, phi0) =
    -1/(u*u) - phi0*phi0/3 * (1 - u * xi) + u * xi * phi0*phi0/3 +
    3 * u*u * S_in + u*u*u * S_u_in

@inline F_inner_to_outer(F_in, u) = u*u * F_in

@inline F_u_inner_to_outer(F_u_in, F_in, u) = 2 * u * F_in + u*u * F_u_in

@inline Sd_inner_to_outer(Sd_in, u, xi, phi0) =
    1/(2*u*u) + xi/u + xi*xi/2 - phi0*phi0/6 + u*u * Sd_in

@inline Bd_inner_to_outer(Bd_in, u) = u*u*u * Bd_in

@inline phid_inner_to_outer(phid_in, u, phi0) =
    -phi0/2 + u*u * phi0*phi0*phi0 * phid_in

@inline A_inner_to_outer(A_in, u, xi, phi0) =
    1/(u*u) + 2*xi/u + xi*xi - 2/3 * phi0*phi0 + u*u * A_in

@inline A_u_inner_to_outer(A_u_in, A_in, u, xi, phi0) =
    -2/(u*u*u) - 2*xi/(u*u) + 2*u * A_in + u*u * A_u_in


text:
export iauS2p
"""
Convert spherical polar coordinates to p-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   theta   double       longitude angle (radians)
   phi     double       latitude angle (radians)
   r       double       radial distance

Returned:
   p       double[3]    Cartesian coordinates

Called:
   iauS2c       spherical coordinates to unit vector
   iauSxp       multiply p-vector by scalar

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauS2p(double theta, double phi, double r, double p[3])
function iauS2p(theta::Real, phi::Real, r::Real)
   p = zeros(Float64, 3)

   ccall((:iauS2p, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, theta),
       convert(Float64, phi),
       convert(Float64, r),
       p)

   return SVector{3}(p)
end


text:
@inline function user_is_authorized_to_register(username, repo; auth)
    gh_user = GitHub.owner(username; auth)
    gh_repo = GitHub.repo(repo; auth)
    return user_is_authorized_to_register(gh_user, gh_repo; auth)
end

@inline function user_is_authorized_to_register(gh_user::GitHub.Owner, gh_repo::GitHub.Repo; auth)
    repo_owner = gh_repo.owner
    if repo_owner.login == gh_user.login
        return true
    end
    if repo_owner.typ == "Organization"
        if GitHub.check_membership(repo_owner, gh_user; auth)
            return true
        end
    end
    return false
end

@inline function _get_pr_title_and_body(commit_message::AbstractString)
    lines = split(strip(commit_message), '\n')
    first_line = lines[1]
    remaining_lines = lines[2:end]
    pr_title = strip(first_line)
    pr_body = strip(join(remaining_lines, "\n"))
    return pr_title, pr_body
end

"""
    register(; user, issue_title, kwargs...)

Register a new package or a new version of an existing package.
"""
@inline function register(; issue_title = nothing, kwargs...)
    kwargs_dict = Dict(kwargs)
    user = User!(kwargs_dict)
    config = Config!(kwargs_dict)
    parse_trigger_issue_title!(kwargs_dict, issue_title)
    new_kwargs = NamedTuple(kwargs_dict)
    registration_request = RegistrationRequest(; new_kwargs...)
    return register(config, user, registration_request)
end

@inline function register(config::Config, user::User, registration_request::RegistrationRequest)
    auth = GitHub.authenticate(config)
    if user_is_authorized_to_register(user.username, registration_request.repo; auth)
        whoami = GitHub.whoami(; auth)
        @info("Authenticated to GitHub as $(whoami.login)")
        @info("User is authorized to register", user, registration_request.repo)
        package_url_unauthenticated = "https://$(config.gh_hostname_for_clone)/$(registration_request.repo).git"
        registry_url_authenticated = "https://$(whoami.login):$(config.github_token[])@$(config.gh_hostname_for_clone)/$(config.registry).git"
        with_cloned_repo(package_url_unauthenticated) do cloned_package_path
            with_cloned_repo(registry_url_authenticated) do cloned_registry_path
                cd(cloned_package_path)
                if !isempty(strip(registration_request.branch))
                    run(`git checkout $(strip(registration_request.branch))`)
                end
                run(`git config user.name "JuliaRegistrator"`)
                run(`git config user.email "contact@julialang.org"`)
                cd(cloned_registry_path)
                run(`git checkout master`)
                pr_branch_name = get_pr_branch_name(cloned_package_path)
                run(`git checkout -B $(pr_branch_name)`)
                run(`git config user.name "Registrator"`)
                run(`git config user.email "ci@juliacomputing.com"`)

                cloned_registry = ClonedRegistry(cloned_registry_path)
                cloned_package = ClonedPackage(cloned_package_path)
                @info "" cloned_registry.path cloned_package.path
                cd(cloned_registry_path)
                foo = LocalRegistry.register(
                    cloned_package,
                    cloned_registry;
                    commit = true,
                    push = false,
                    repo = package_url_unauthenticated,
                )
                cd(cloned_registry_path)
                commit_message = read(`git log -1 --pretty=%B`, String)
                pr_title, pr_body = _get_pr_title_and_body(commit_message)
                pr_body *= "\nCreated by: @$(user.username)"
                registry_gh_repo = GitHub.repo(config.registry)
                pr_params = Dict()
                pr_params["head"] = pr_branch_name
                pr_params["base"] = "master"
                pr_params["body"] = pr_body
                pr_params["maintainer_can_modify"] = true
                pr_params["title"] = pr_title
                run(`git push -f -u origin $(pr_branch_name)`)
                pr = create_or_update_pull_request(registry_gh_repo, pr_params; auth)
            end
        end
    else
        msg = """
        The user is not authorized to make new releases of this package.

        In order to use the serverless Registrator, you need to either be the
        owner of the repository, or you need to be a member of the organization
        that owns the repository. (Please note: if you are a member of the organization,
        your membership must be public.)

        If you do not meet either of those criteria, you will need to use either
        the Registrator comment bot or the Registrator web interface.
        """
        throw(UserNotAuthorizedError(msg, user, registration_request.repo))
    end
    return nothing
end


text:
export GLMCopulaCSObs, GLMCopulaCSModel

# abstract type QCobs end

struct GLMCopulaCSObs{T <: BlasReal, D, Link} # <: QCobs
    # data
    n::Int
    p::Int
    y::Vector{T}
    X::Matrix{T}
    V::SymmetricToeplitz{T}
    vec::Vector{T}
    # working arrays
    ∇CSV::SymmetricToeplitz{T}
    ∇β::Vector{T}   # gradient wrt β
    ∇resβ::Matrix{T}# residual gradient matrix d/dβ_p res_ij (each observation has a gradient of residual is px1)
    ∇ρ::Vector{T}
    ∇σ2::Vector{T}
    Hβ::Matrix{T}   # Hessian wrt β
    Hρ::Matrix{T}   # Hessian wrt ρ
    Hσ2::Matrix{T}   # Hessian wrt ρ
    Hρσ2::Matrix{T}  # Hessian wrt ρ, σ2
    Hβσ2::Vector{T}  # Hessian wrt β and σ2
    Hβρ::Vector{T}
    res::Vector{T}  # residual vector res_i
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    xtx::Matrix{T}  # Xi'Xi
    storage_n::Vector{T}
    storage_p1::Vector{T}
    storage_np::Matrix{T}
    storage_pp::Matrix{T}
    added_term_numerator::Matrix{T}
    added_term2::Matrix{T}
    η::Vector{T}    # η = Xβ systematic component
    μ::Vector{T}    # μ(β) = ginv(Xβ) # inverse link of the systematic component
    varμ::Vector{T} # v(μ_i) # variance as a function of the mean
    dμ::Vector{T}   # derivative of μ
    d::D            # distribution()
    link::Link      # link function ()
    wt::Vector{T}   # weights wt for GLM.jl
    w1::Vector{T}   # working weights in the gradient = dμ/v(μ)
    w2::Vector{T}   # working weights in the information matrix = dμ^2/v(μ)
end

function GLMCopulaCSObs(
    y::Vector{T},
    X::Matrix{T},
    d::D,
    link::Link) where {T <: BlasReal, D, Link}
    n, p = size(X, 1), size(X, 2)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    V = SymmetricToeplitz(ones(T, n))
    vec = Vector{T}(undef, n)
    ∇CSV = SymmetricToeplitz(ones(T, n))
    ∇β  = Vector{T}(undef, p)
    ∇resβ  = Matrix{T}(undef, n, p)
    ∇ρ  = Vector{T}(undef, 1)
    ∇σ2  = Vector{T}(undef, 1)
    Hβ  = Matrix{T}(undef, p, p)
    Hρ  = Matrix{T}(undef, 1, 1)
    Hσ2  = Matrix{T}(undef, 1, 1)
    Hρσ2 = Matrix{T}(undef, 1, 1)
    Hβσ2 = zeros(T, p)
    Hβρ = zeros(T, p)
    res = Vector{T}(undef, n)
    t   = [tr(V)/2]
    q   = Vector{T}(undef, 1)
    xtx = transpose(X) * X
    storage_n = Vector{T}(undef, n)
    storage_p1 = Vector{T}(undef, p)
    storage_np = Matrix{T}(undef, n, p)
    storage_pp = Matrix{T}(undef, p, p)
    added_term_numerator = Matrix{T}(undef, n, p)
    added_term2 = Matrix{T}(undef, p, p)
    η = Vector{T}(undef, n)
    μ = Vector{T}(undef, n)
    varμ = Vector{T}(undef, n)
    dμ = Vector{T}(undef, n)
    wt = Vector{T}(undef, n)
    fill!(wt, one(T))
    w1 = Vector{T}(undef, n)
    w2 = Vector{T}(undef, n)
    # constructor
    GLMCopulaCSObs{T, D, Link}(n, p, y, X, V, vec, ∇CSV, ∇β, ∇resβ, ∇ρ, ∇σ2, Hβ, Hρ, Hσ2, Hρσ2, Hβσ2, Hβρ,
       res, t, q, xtx, storage_n, storage_p1, storage_np, storage_pp, added_term_numerator, added_term2,
        η, μ, varμ, dμ, d, link, wt, w1, w2)
end

"""
    GLMCopulaCSModel
GLMCopulaCSModel(gcs)
Gaussian copula compound symmetric (CS) model, which contains a vector of
`GLMCopulaCSObs` as data, model parameters, and working arrays.
"""
struct GLMCopulaCSModel{T <: BlasReal, D, Link} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GLMCopulaCSObs{T, D, Link}}
    Ytotal::T
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    ρ::Vector{T}            # autocorrelation parameter
    σ2::Vector{T}           # autoregressive noise parameter
    θ::Vector{T}
    # working arrays
    ∇β::Vector{T}   # gradient of beta from all observations
    ∇ρ::Vector{T}           # gradient of rho from all observations
    ∇σ2::Vector{T}          # gradient of sigmasquared from all observations
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    Hβ::Matrix{T}    # Hessian from all observations
    Hρ::Matrix{T}    # Hessian from all observations
    Hσ2::Matrix{T}    # Hessian from all observations
    Hρσ2::Matrix{T}
    Hβσ2::Vector{T}
    Hβρ::Vector{T}
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    TR::Matrix{T}
    QF::Matrix{T}         # n-by-1 matrix with qik = res_i' Vi res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    d::Vector{D}
    link::Vector{Link}
    penalized::Bool
end

function GLMCopulaCSModel(gcs::Vector{GLMCopulaCSObs{T, D, Link}}; penalized::Bool = false) where {T <: BlasReal, D, Link}
    n, p = length(gcs), size(gcs[1].X, 2)
    β   = Vector{T}(undef, p)
    τ   = [1.0]
    ρ = [1.0]
    σ2 = [1.0]
    θ   = Vector{T}(undef, 1)
    ∇β  = Vector{T}(undef, p)
    ∇ρ  = Vector{T}(undef, 1)
    ∇σ2  = Vector{T}(undef, 1)
    XtX = zeros(T, p, p) # sum_i xi'xi
    Hβ  = Matrix{T}(undef, p, p)
    Hρ  = Matrix{T}(undef, 1, 1)
    Hσ2  = Matrix{T}(undef, 1, 1)
    Hρσ2 = Matrix{T}(undef, 1, 1)
    Hβσ2 = zeros(T, p)
    Hβρ = zeros(T, p)
    Ainv    = zeros(T, p + 2, p + 2)
    Aevec   = zeros(T, p + 2, p + 2)
    M       = zeros(T, p + 2, p + 2)
    vcov    = zeros(T, p + 2, p + 2)
    ψ       = Vector{T}(undef, p + 2)
    TR  = Matrix{T}(undef, n, 1) # collect trace terms
    QF  = Matrix{T}(undef, n, 1)
    Ytotal = 0.0
    ntotal = 0.0
    d = Vector{D}(undef, n)
    link = Vector{Link}(undef, n)
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        Ytotal  += sum(gcs[i].y)
        BLAS.axpy!(one(T), gcs[i].xtx, XtX)
        TR[i, :] = gcs[i].t
        d[i] = gcs[i].d
        link[i] = gcs[i].link
    end
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, 1)
    storage_θ = Vector{T}(undef, 1)
    GLMCopulaCSModel{T, D, Link}(gcs, Ytotal, ntotal, p, β, τ, ρ, σ2, θ,
        ∇β, ∇ρ, ∇σ2, XtX, Hβ, Hρ, Hσ2, Hρσ2, Hβσ2, Hβρ, Ainv, Aevec,  M, vcov, ψ,
        TR, QF, storage_n, storage_m, storage_θ, d, link, penalized)
end

# use ToeplitzMatrices
"""
    get_V!(ρ, gc)
Forms the CS covariance structure given ρ (correlation parameter), gc (single cluster observation) object.
"""
function get_V!(ρ, gc::Union{GLMCopulaCSObs{T, D, Link}, NBCopulaCSObs{T, D, Link}, GaussianCopulaCSObs{T}}) where {T, D, Link}
    gc.vec[1] = 1.0
    @inbounds for i in 2:gc.n
        gc.vec[i] = ρ
    end
    gc.V.vc .= gc.vec
end

"""
    get_∇V!(gc)
Forms the first derivative of CS covariance structure wrt to ρ, given ρ (correlation parameter)
"""
function get_∇V!(gc::Union{GLMCopulaCSObs{T, D, Link}, NBCopulaCSObs{T, D, Link}, GaussianCopulaCSObs{T}}) where {T, D, Link}
    gc.vec[1] = 0.0
    @inbounds for i in 2:gc.n
        gc.vec[i] = 1.0
    end
    gc.∇CSV.vc .= gc.vec
end

function loglikelihood!(
    gc::GLMCopulaCSObs,
    β::Vector{T},
    ρ::T,
    σ2::T,
    needgrad::Bool = false,
    needhess::Bool = false;
    penalized::Bool = false
    ) where {T <: BlasReal}
    # n, p = size(gc.X, 1), size(gc.X, 2)
    needgrad = needgrad || needhess
    if needgrad
        fill!(gc.∇β, 0)
        get_∇V!(gc)
    end
    needhess && fill!(gc.Hβ, 0)
    fill!(gc.∇β, 0.0)
    update_res!(gc, β)
    standardize_res!(gc)
    fill!(gc.∇resβ, 0.0) # fill gradient of residual vector with 0
    std_res_differential!(gc) # this will compute ∇resβ

    # form V
    get_V!(ρ, gc)

    #evaluate copula loglikelihood
    mul!(gc.storage_n, gc.V, gc.res, one(T), zero(T)) # storage_n = V[k] * res
    if needgrad
        BLAS.gemv!('T', σ2, gc.∇resβ, gc.storage_n, 1.0, gc.∇β) # stores ∇resβ*Γ*res (standardized residual)
    end

    q = dot(gc.res, gc.storage_n)
    c1 = 1 + (0.5 * gc.n * σ2)
    # if q < 0
    #     q = 0
    # end
    c2 = 1 + (0.5 * σ2 * q)
    # loglikelihood
    logl = QuasiCopula.component_loglikelihood(gc)
    logl += -log(c1)
    # @show logl
    logl += log(c2)
    # add L2 ridge penalty
    if penalized
        logl -= 0.5 * (σ2)^2
    end
    # @show logl
    if needgrad
        inv1pq = inv(c2)
        # gradient with respect to rho
        mul!(gc.storage_n, gc.∇CSV, gc.res, one(T), zero(T)) # storage_n = ∇CSV * res
        q2 = dot(gc.res, gc.storage_n) #
        # gc.∇ρ .= inv(c2) * 0.5 * σ2 * transpose(gc.res) * gc.∇CSV * gc.res
        gc.∇ρ .= inv1pq * 0.5 * σ2 * q2
        # gradient with respect to sigma2
        gc.∇σ2 .= 0.5 * ((q * inv1pq) - gc.n * inv(c1))
        if penalized
            gc.∇σ2 .-= σ2
        end
      if needhess
            gc.Hρ .= -(inv(c2) * 0.5 * σ2 * q2)^2

            # hessian for sigma2
            gc.Hσ2 .= 0.25 * gc.n^2 * inv(c1)^2 - inv1pq^2 * (0.25 * q^2)
            if penalized
                gc.Hσ2 .-= 1.0
            end
            # hessian cross term for rho and sigma2
            # gc.Hρσ2 .= 0.5 * q2 * inv1pq - 0.25 * σ2 * q * q2 * inv1pq^2
            gc.Hρσ2 .= (0.5 * q2 * c2 - 0.25 * σ2 * q * q2) * inv1pq^2

            # hessian cross term for beta and sigma2
            # gc.Hβσ2 .= inv1pq * gc.∇β - 0.5 * q * inv1pq^2 * σ2 * gc.∇β
            gc.Hβσ2 .= ((c2 * gc.∇β * inv(σ2)) - (gc.∇β * 0.5 * q)) * inv1pq^2

            #  hessian cross term for beta and rho
            # gc.Hβρ .= inv1pq * σ2 * transpose(gc.∇resβ) * gc.∇ARV * gc.res - 0.5 * σ2^2 * inv1pq^2 * q2 * transpose(gc.∇resβ) * gc.V * gc.res
            BLAS.gemv!('T', σ2, gc.∇resβ, gc.storage_n, 1.0, gc.storage_p1) # stores ∇resβ*Γ*res (standardized residual)

            # gc.Hβρ .= ((c2 * gc.storage_p1) .- (0.5 * gc.∇β * σ2 * q2)) * inv1pq^2

            BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, 0.0, gc.Hβ) # only lower triangular
            fill!(gc.added_term_numerator, 0.0) # fill gradient with 0
            fill!(gc.added_term2, 0.0) # fill hessian with 0
            # gc.V .= get_AR_cov(n, ρ, σ2, gc.V)
            mul!(gc.added_term_numerator, gc.V, gc.∇resβ, one(T), zero(T)) # storage_n = V[k] * res
            BLAS.gemm!('T', 'N', σ2, gc.∇resβ, gc.added_term_numerator, one(T), gc.added_term2)
            gc.added_term2 .*= inv1pq
            gc.Hβ .+= gc.added_term2
            gc.Hβ .+= QuasiCopula.glm_hessian(gc)
      end
      gc.∇β .= gc.∇β .* inv1pq
      gc.res .= gc.y .- gc.μ
      gc.∇β .+= QuasiCopula.glm_gradient(gc)
    end
    logl
end

function loglikelihood!(
    gcm::GLMCopulaCSModel,
    needgrad::Bool = false,
    needhess::Bool = false
    )
    if needgrad
        fill!(gcm.∇β, 0.0)
        fill!(gcm.∇ρ, 0.0)
        fill!(gcm.∇σ2, 0.0)
    end
    if needhess
        fill!(gcm.Hβ, 0.0)
        fill!(gcm.Hρ, 0.0)
        fill!(gcm.Hσ2, 0.0)
        fill!(gcm.Hρσ2, 0.0)
        fill!(gcm.Hβσ2, 0.0)
        # @show gcm.Hβσ2
        fill!(gcm.Hβρ, 0)
    end
    logl = zeros(Threads.nthreads())
    Threads.@threads for i in eachindex(gcm.data)
        @inbounds logl[Threads.threadid()] += loglikelihood!(gcm.data[i], gcm.β,
         gcm.ρ[1], gcm.σ2[1], needgrad, needhess; penalized = gcm.penalized)
     end
     @inbounds for i in eachindex(gcm.data)
         if needgrad
             gcm.∇β .+= gcm.data[i].∇β
             gcm.∇ρ .+= gcm.data[i].∇ρ
             gcm.∇σ2 .+= gcm.data[i].∇σ2
         end
         if needhess
             gcm.Hβ .+= gcm.data[i].Hβ
             gcm.Hρ .+= gcm.data[i].Hρ
             gcm.Hσ2 .+= gcm.data[i].Hσ2
             gcm.Hρσ2 .+= gcm.data[i].Hρσ2
             gcm.Hβσ2 .+= gcm.data[i].Hβσ2
         end
     end
    return sum(logl)
end


text:
# This file was generated, do not modify it. # hide
lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])


text:
include("../btg.jl")
include("constrained_guess.jl")
include("sample.jl")
include("../computation/finitedifference.jl")
include("visualize.jl")

using Plots
using Random
using Ipopt, JuMP
import JuMP.MathOptInterface
const MOI = MathOptInterface

include("optUCB.jl")
include("updateBTG.jl")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
immutable FDist <: ContinuousUnivariateDistribution
    ndf::Float64
    ddf::Float64
    function FDist(d1::Real, d2::Real)
    	if d1 > 0 && d2 > 0
    		new(float64(d1), float64(d2))
    	else
    		error("Both numerator and denominator degrees of freedom must be positive")
    	end
    end
end

@_jl_dist_2p FDist f

insupport(d::FDist, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x

mean(d::FDist) = 2.0 < d.ddf ? d.ddf / (d.ddf - 2.0) : NaN

function var(d::FDist)
	if 4.0 < d.ddf
		return 2.0 * d.ddf^2 *
		       (d.ndf + d.ddf - 2.0) /
		       (d.ndf * (d.ddf - 2.0)^2 * (d.ddf - 4.0))
	else
		return NaN
	end
end


text:
@testset "Single Criterion Decision Making Functions" begin
    @testset "Single Criterion Decision Making tools" begin
        @testset "Laplace" begin

            tol = 0.00001

            mat = [
                3000 2750 2500 2250
                1500 4750 8000 7750
                2000 5250 8500 11750
            ]

            dm = makeDecisionMatrix(mat)

            result = laplace(dm)

            @test isa(result, LaplaceResult)

            @test isapprox(result.expected_values, [2625.0, 5500, 6875], atol = tol)

            @test result.bestIndex == 3
        end

        @testset "Maximin" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = maximin(dm)

            @test isa(result, MaximinResult)

            @test isapprox(result.rowmins, [18, 18, 24, 20], atol = tol)

            @test result.bestIndex == 3
        end

        @testset "Maximax" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = maximax(dm)

            @test isa(result, MaximaxResult)

            @test isapprox(result.rowmaxs, [26, 34, 34, 30], atol = tol)

            @test result.bestIndex in [2, 3]
        end

        @testset "Minimax" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = minimax(dm)

            @test isa(result, MinimaxResult)

            @test isapprox(result.rowmaxs, [26, 34, 34, 30], atol = tol)

            @test result.bestIndex == 1
        end

        @testset "Minimin" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = minimin(dm)

            @test isa(result, MiniminResult)

            @test isapprox(result.rowmins, [18, 18, 24, 20], atol = tol)

            @test result.bestIndex in [1, 2]
        end

        @testset "Savage" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = savage(dm)

            @test isa(result, SavageResult)

            @test result.bestIndex == 4

        end

        @testset "Hurwicz" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            dm = makeDecisionMatrix(mat)

            result = hurwicz(dm)

            @test isa(result, HurwiczResult)

            @test result.bestIndex == 3

        end

        @testset "Maximum likelihood for single criterion" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            weights = [0.2, 0.5, 0.2, 0.1]

            dm = makeDecisionMatrix(mat)

            result = mle(dm, weights)

            @test isa(result, MLEResult)

            @test result.bestIndex == 2

            @test isapprox(result.scores, [24, 29.2, 27, 27], atol = tol)
        end

        @testset "Expected Regret" begin

            tol = 0.00001

            mat = [
                26 26 18 22
                22 34 30 18
                28 24 34 26
                22 30 28 20
            ]

            weights = [0.2, 0.5, 0.2, 0.1]

            dm = makeDecisionMatrix(mat)

            result = expectedregret(dm, weights)

            @test isa(result, ExpectedRegretResult)

            @test result.bestIndex == 2

            @test isapprox(result.scores, [8, 2.8, 5, 5], atol = tol)
        end
    end
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
###################################################################################################
# MULTIPLE HYPOTHESIS CORRECTION
###################################################################################################
"""
    `add_qvalues_to_path(PATH)`

    Function that takes in the differential output and adds BH q-value.

    # Examples
    ```julia-repl
    julia> CpelNano.add_qvalues_to_path(path)
    ```
"""
function add_qvalues_to_path(path::String)::Nothing

    # Leave if no data
    filesize(path) > 0 || return nothing

    # Get data
    all_data = readdlm(path, '\t', Any)
    qvals = fill(NaN, size(all_data)[1])

    # Multiple hypothesis testing correction
    ind = .!isnan.(all_data[:,5])
    if sum(ind) > 0
        qvals[ind] = MultipleTesting.adjust(convert(Vector{Float64}, all_data[ind,5]), BenjaminiHochberg())
    end

    # Append to output matrix
    all_data = hcat(all_data, qvals)

    # Write to temp output
    temp_path = path * ".tmp"
    open(temp_path, "w") do io
        writedlm(io, all_data, '\t')
    end

    # Move to original file
    mv(temp_path, path, force=true)

    # Return
    return nothing

end
"""
    `mult_hyp_corr(CONFIG)`

    Function that takes in all the differential output and adds BH q-value in each one.

    # Examples
    ```julia-repl
    julia> CpelNano.mult_hyp_corr(config)
    ```
"""
function mult_hyp_corr(config::CpelNanoConfig)::Nothing

    # Add q-values
    add_qvalues_to_path(config.out_diff_files.tmml_file)
    add_qvalues_to_path(config.out_diff_files.tnme_file)
    config.matched || add_qvalues_to_path(config.out_diff_files.tcmd_file)

    # Return
    return nothing

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(command_main),Array{String,1}})
end


text:
x = create_var(1, 1)
y = create_var(3, 1)
s = sum_squares(x) + sum_squares(y)
test_expr(s)

s = 4 * s
test_expr(s)

s = s * 2
test_expr(s)

s = s / 4
test_expr(s)

s = 0
s += sum_squares(x)
s += sum_squares(y)
s += 5
s += sum_squares(9)
s += sum_squares()
test_expr(s)

s = sum_squares(x)
@test_throws ErrorException  s + (-3)

info("All sum squares tests passed")


text:
using DataFrames, DataFramesMeta, StatsPlots, StatsBase, CSV

function main()
    df_c = CSV.read(joinpath(homedir(), "Desktop", "c.csv"), DataFrame)
    df_bf = CSV.read(joinpath(homedir(), "Desktop", "bf.csv"), DataFrame)

    df_c_int8 = @where(df_c, :int_type .== "Int8")
    df_c_int128 = @where(df_c, :int_type .== "Int128")
    df_bf_int8 = @where(df_bf, :int_type .== "Int8")

    df_c_int_types = [first(i.int_type) for i in groupby(df_c, :int_type)]
    df_c_runtime_avgs = [mean(i.run_time_ms) for i in groupby(df_c, :int_type)]
    df_bf_int_types = [first(i.int_type) for i in groupby(df_bf, :int_type)]
    df_bf_runtime_avgs = [mean(i.run_time_ms) for i in groupby(df_bf, :int_type)]

    println(df_bf_int_types)
    println(df_bf_runtime_avgs)

    plt = bar(
        # df_c.int_type,
        # df_c.run_time_ms,
        # categories,
        # unique(df_c.int_type),
        df_c_int_types,
        # Int[get(frequencies, c, 0) for c in categories],
        # df_c.run_time_ms,
        df_c_runtime_avgs,
        yaxis=:log,
        # legend = false,
        label = "Clever",
        xticks = (0.5:length(df_c_int_types), df_c_int_types),
        xrotation = 45,
        fontfamily = font("Times"),
        xlabel = "Int type",
        ylabel = "Run time in ms",
        title = "Run time analysis on varying integer types"
    )
    bar!(
        # df_bf.int_type,
        df_bf_int_types,
        # [first(i.int_type) for i in groupby(df_bf, :int_type)],
        df_bf_runtime_avgs,
        # df_c.run_time_ms,
        # [mean(i.run_time_ms) for i in groupby(df_bf, :int_type)],
        label = "Brute force"
    )

    savefig(plt, joinpath(homedir(), "Desktop", "test.pdf"))

    return nothing
end

main()

max_stars_repo_name: jakewilliami/netlify

max_issues_repo_name: jakewilliami/netlify

max_forks_repo_name: jakewilliami/netlify


text:
mutable struct AdaptiveCholeskyFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LayoutMatrix{T}
    data::CachedMatrix{T,DM,M}
    ncols::Int
end

size(U::AdaptiveCholeskyFactors) = size(U.data.array)
bandwidths(A::AdaptiveCholeskyFactors) = (0,bandwidth(A.data,2))

function AdaptiveCholeskyFactors(::SymmetricLayout{<:AbstractBandedLayout}, S::AbstractMatrix{T}) where T
    A = parent(S)
    l,u = bandwidths(A)
    data = BandedMatrix{T}(undef,(0,0),(l,u)) # pad super
    AdaptiveCholeskyFactors(CachedArray(data,A), 0)
end
AdaptiveCholeskyFactors(A::AbstractMatrix{T}) where T = AdaptiveCholeskyFactors(MemoryLayout(A), A)
MemoryLayout(::Type{AdaptiveCholeskyFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()

function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int) where T
    if n > F.ncols
        _,u = bandwidths(F.data.array)
        resizedata!(F.data,n+u,n+u);
        ncols = F.ncols
        kr = ncols+1:n
        factors = view(F.data.data,kr,kr)
        banded_chol!(factors, UpperTriangular)
        # multiply remaining columns
        kr2 = max(n-u+1,kr[1]):n
        U1 = UpperTriangular(view(F.data.data,kr2,kr2))
        B = view(F.data.data,kr2,n+1:n+u)
        ldiv!(U1',B)
        muladd!(-one(T),B',B,one(T),view(F.data.data,n+1:n+u,n+1:n+u))
        F.ncols = n
    end
    F
end

function getindex(F::AdaptiveCholeskyFactors, k::Int, j::Int)
    partialcholesky!(F, max(k,j))
    F.data.data[k,j]
end

adaptivecholesky(A) = Cholesky(AdaptiveCholeskyFactors(A), :U, 0)

ArrayLayouts._cholesky(::SymmetricLayout{<:AbstractBandedLayout}, ::NTuple{2,OneToInf{Int}}, A) = adaptivecholesky(A)

function colsupport(F::AdaptiveCholeskyFactors, j)
    partialcholesky!(F, maximum(j)+bandwidth(F,2))
    colsupport(F.data.data, j)
end

function rowsupport(F::AdaptiveCholeskyFactors, j)
    partialcholesky!(F, maximum(j)+bandwidth(F,2))
    rowsupport(F.data.data, j)
end

colsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = rowsupport(parent(F), j)
rowsupport(F::AdjOrTrans{<:Any,<:AdaptiveCholeskyFactors}, j) = colsupport(parent(F), j)

function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',<:AdaptiveLayout},<:PaddedLayout})
    A,B = M.A,M.B
    T = eltype(M)
    COLGROWTH = 1000 # rate to grow columns
    tol = floatmin(real(T))

    require_one_based_indexing(B)
    mA, nA = size(A)
    mB = length(B)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    sB = B.datasize[1]
    l,_ = bandwidths(A)

    jr = 1:min(COLGROWTH,nA)

    P = parent(parent(A))

    @inbounds begin
        while first(jr) < nA
            j = first(jr)
            cs = colsupport(A, last(jr))
            cs_max = maximum(cs)
            kr = j:cs_max
            resizedata!(B, min(cs_max,mB))
            if (j > sB && maximum(abs,view(B.data,j:last(colsupport(P,j)))) ≤ tol)
                break
            end
            partialcholesky!(P, min(cs_max,nA))
            U_N = UpperTriangular(view(P.data.data, jr, jr))
            ldiv!(U_N', view(B.data, jr))
            jr1 = last(jr)-l+1:last(jr)
            jr2 = last(jr)+1:last(jr)+l
            muladd!(-one(T), view(P.data.data, jr1,jr2)', view(B.data,jr1), one(T), view(B.data,jr2))
            jr = last(jr)+1:min(last(jr)+COLGROWTH,nA)
        end
    end
    B
end

function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveCholeskyFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}})
    n = B.datasize[1]
    partialcholesky!(parent(R), n)
    ldiv!(UpperTriangular(view(parent(R).data.data,oneto(n),oneto(n))), view(B.data,oneto(n)))
    B
end


text:
export AbstractShape
export level, ndir, project, center, bounds

abstract type AbstractShape{K} end

Base.in(x::AbsVecReal, s::AbstractShape{K}) where {K} = level(x,s) ≤ 0

# Return the value of the level set function of the shape.
level(x::AbsVecReal, s::AbstractShape{K}) where {K} = level(SVec{K}(x), s)

# Return the center position of the shape, where the level set function is minimized.
function center end

# Return the minimum axis-aligned bounding box of the shape.
function bounds end

# Return the outward normal direction, even for x inside the shape.
ndir(x::AbsVecReal, s::AbstractShape{K}) where {K} = ndir(SVec{K}(x), s)

# Project x onto s along the direction normal to the shape.  The result is only
# approximate.
project(x::AbsVecReal, s::AbstractShape{K}) where {K} = project(SVec{K}(x), s)

include("translated.jl")
include("simple/simple.jl")
include("composite/composite.jl")

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
import StatsBase: mean

"""
    NormalisedPredictiveAsymmetryTest(predictive_test::CausalityTest; f::Number = 1.0)

The parameters for a normalised predictive asymmetry causality test [1]. For the
non-normalised version, see [`PredictiveAsymmetryTest`](@ref).

## Mandatory keywords

- **`predictive_test`**: An instance of a predictive causality test that explicitly
    uses prediction lags (e.g. [`VisitationFrequencyTest`](@ref) or
    [`TransferOperatorGridTest`](@ref)).
- **`f::Number`**: The normalisation factor.

## About the prediction lags

The prediction lags in the predictive causality test must consist of `n` negative
integers and `n` positive integers that are symmetric around zero.

In other words, negative lags  must exactly match the positive lags but with opposite
sign. The zero lag can be included, but will be ignored, so it is possible to give
ranges too.

## Examples

Let's define a a transfer entropy test with which to compute transfer entropies
symmetrically around lag zero. We'll use that as input to the predictive asymmetry
test.

```julia
bin = RectangularBinning(4) # divide each axis into 4 equal-length intervals
ηs = [-5:-1; 1:5] # exclude the zero lag (it is irrelevant for the asymmetry)
test_visitfreq = VisitationFrequencyTest(binning = bin, ηs = ηs)

# Normalising predictive asymmetry values to f = 1.0 times the mean transfer entropy.
NormalisedPredictiveAsymmetryTest(predictive_test = test_visitfreq, f = 1.0)
```

## References

1. Diego, David, Kristian Agasøster Haaga, Jo Brendryen, and Bjarte Hannisdal.
    A simple test for causal asymmetry in complex systems. In prep.
"""
Base.@kwdef mutable struct NormalisedPredictiveAsymmetryTest{TEST, N} <: AbstractPredictiveAsymmetryTest where TEST
    predictive_test::TEST
    f::Number

    # If a threshold is given, use it.
    function NormalisedPredictiveAsymmetryTest(test::T; f::Number) where {T <: TransferEntropyCausalityTest}
        # Check that prediction lags are okay
        verified_prediction_lags(test.ηs)
        N = length(test.ηs[test.ηs .> 0])
        return new{T, N}(test, f)
    end

    function NormalisedPredictiveAsymmetryTest(test::T; f::Number) where {T <: TransferEntropyCausalityTest}
        # Check that prediction lags are okay
        verified_prediction_lags(test.ηs)
        N = length(test.ηs[test.ηs .> 0])
        return new{T, N}(test, f)
    end
end

"""
    lagnormalised_statistic(values, ηs, f)

In a cumulative manner, normalise the `values` for the statistic
to some fraction `f` of their mean over prediction lags `ηs`.
"""
function lagnormalised_statistic(values, ηs, f)
    length(values) == length(ηs) || error("length(values) must equal length(ηs)")
    avg_vals = zeros(Float64, length(ηs[ηs .> 0]))

    for (i, η) in enumerate(ηs[ηs .> 0])
        startidx = findfirst(ηs .== -η)
        stopidx  = findlast(ηs .== η)
        avg_vals[i] = mean(values[startidx:stopidx])*f
    end

    return avg_vals
end

"""
    normalised_predictive_asymmetry(source, target, test::PredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

Compute the predictive asymmetry from `source` to `target` using the provided predictive test `p`.
The test can be a [`VisitationFrequencyTest`](@ref), [`TransferOperatorGridTest`](@ref) or
[`NearestNeighbourMITest`](@ref) and has to be defined for prediction lags symmetrically
around zero.

## Example

We'll use the [`logistic4`](@ref) system, which consists of three variables
and is unidirectionally coupled x -> y -> z.

```julia
# Example data
sys = logistic()
npts = 300
orbit = trajectory(sys, npts, Ttr = 300);
x, y, z = columns(orbit);

# Test setup
η_max = 10
ηs = [-η_max:-1; 1:η_max]
bin = RectangularBinning(floor(Int, npts^(1/4)))
test = VisitationFrequencyTest(ηs = ηs, binning = bin)
pa_test = PredictiveAsymmetryTest(test)

# Analysis
pas_xy = causality(x, y, pa_test)
pas_yx = causality(y, x, pa_test)
pas_yz = causality(y, z, pa_test)
pas_zy = causality(z, y, pa_test)

[pas_xy pas_yx pas_yz pas_zy]
````

If the test correctly picked up the correct directionality, then the values of
of `asymmetries` corresponding to the causal interactions (1st and 3rd column)
should be mostly positive, and those corresponding to non-causal interactions
should be mostly negative.
"""
function normalised_predictive_asymmetry(source, target,
        p::NormalisedPredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

    # Update the test parameters so that we have symmetric prediction lags
    test = update_ηs(p.predictive_test)

    # Predictions from source to target
    preds = causality(source, target, test)

    # Normalised predictions
    normalised_preds = lagnormalised_statistic(preds, test.ηs, p.f)

    # The number of predictive asymmetries is half the number of prediction lags,
    # which is encoded in the type parameter `N`
    ηs = test.ηs
    As = zeros(Float64, N)

    for (i, η) in enumerate(ηs[ηs .> 0])
        lag_idxs = findfirst(ηs .== -η):findfirst(ηs .== η)
        As[i] = causalbalance(ηs[lag_idxs], preds[lag_idxs]) ./ normalised_preds[i]
    end

    return return_predictive_asymmetry(p.predictive_test.ηs, As, N)
end

function normalised_predictive_asymmetry(source, target, cond,
    p::NormalisedPredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

    # Update the test parameters so that we have symmetric prediction lags
    test = update_ηs(p.predictive_test)

    # Predictions from source to target
    preds = causality(source, target, cond, test)

    # Normalised predictions
    normalised_preds = lagnormalised_statistic(preds, test.ηs, p.f)

    # The number of predictive asymmetries is half the number of prediction lags,
    # which is encoded in the type parameter `N`
    ηs = test.ηs
    As = zeros(Float64, N)

    for (i, η) in enumerate(ηs[ηs .> 0])
        lag_idxs = findfirst(ηs .== -η):findfirst(ηs .== η)
        As[i] = causalbalance(ηs[lag_idxs], preds[lag_idxs]) ./ normalised_preds[i]
    end

    return return_predictive_asymmetry(p.predictive_test.ηs, As, N)
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, p::NormalisedPredictiveAsymmetryTest{CT}) where {T<:Real, CT}
    normalised_predictive_asymmetry(source, target, p)
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, cond::AbstractVector{T}, p::NormalisedPredictiveAsymmetryTest{CT}) where {T<:Real, CT}
    normalised_predictive_asymmetry(source, target, cond, p)
end

export NormalisedPredictiveAsymmetryTest


text:
using HDF5

@defstruct HDF5OutputLayer Layer (
  name :: AbstractString = "hdf5-output",
  (bottoms :: Vector{Symbol} = [], length(bottoms) > 0),
  (datasets :: Vector{Symbol} = Symbol[], length(datasets) == 0 || length(datasets) == length(bottoms)),
  (filename :: AbstractString = "", !isempty(filename)),
  force_overwrite :: Bool = false,
)
@characterize_layer(HDF5OutputLayer,
  is_sink => true
)

type HDF5OutputLayerState <: LayerState
  layer  :: HDF5OutputLayer
  file   :: HDF5File
  buffer :: Vector{Array}
  dsets  :: Vector{Any}
  index  :: Int
end

function setup(backend::Backend, layer::HDF5OutputLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  if isfile(layer.filename)
    if !layer.force_overwrite
      error("HDF5OutputLayer: output file '$(layer.filename)' already exists")
    else
      @warn("HDF5OutputLayer: output file '$(layer.filename)' already exists, overwriting")
    end
  end
  if length(layer.datasets) == 0
    datasets = layer.bottoms
  else
    datasets = layer.datasets
  end

  file = h5open(layer.filename, "w")
  buffer = Array{Array}(length(inputs))
  dsets  = Array{Any}(length(inputs))
  for i = 1:length(inputs)
    data_type = eltype(inputs[i])
    dims = size(inputs[i])
    dsets[i] = d_create(file, string(datasets[i]), datatype(data_type),
        dataspace(dims, max_dims=tuple(dims[1:end-1]..., -1)), "chunk", dims)
    buffer[i] = Array{data_type}(dims)
  end

  return HDF5OutputLayerState(layer, file, buffer, dsets, 1)
end

function forward(backend::Backend, state::HDF5OutputLayerState, inputs::Vector{Blob})
  for i = 1:length(inputs)
    copy!(state.buffer[i], inputs[i])
    dims = size(state.buffer[i])
    batch_size = dims[end]

    # extend the HDF5 dataset
    set_dims!(state.dsets[i], tuple(dims[1:end-1]..., state.index*batch_size))

    # write data
    idx = map(x -> 1:x, dims[1:end-1])
    state.dsets[i][idx...,(state.index-1)*batch_size+1:state.index*batch_size] = state.buffer[i]
  end
  state.index += 1
end

function backward(backend::Backend, state::HDF5OutputLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end

function shutdown(backend::Backend, state::HDF5OutputLayerState)
  close(state.file)
end


text:
struct Graph500Results
    scale::Integer
    nbfs::Integer
    kernel_1_time::Float64
    kernel_2_time::Vector{Float64}
    kernel_2_nedge::Vector{Float64}
end

### required for initialising Dict{Graph{T},Graph500Results}() constructor
Graph500Results() = Graph500Results(0, 0, 0.0, zeros(0), zeros(0))

doc"""
    driver(scale, edgefactor; replicate, seed)

Run Graph500 benchmark(http://graph500.org/?page_id=12) on a graph with `2^scale` vertices and
`2^scale*edgefactor` edges. Returns a `ResultStatistic` structure

### Optional Arguments
- `replicate=true`: If true, returns a graph generated with intitial seed = `seed`
"""
function driver(
  scale::Integer,
  edgefactor::Integer;
  replicate::Bool=false,
  seed::AbstractArray=[])

    println("Generating graph.")
    ij = kronecker_generator(scale, edgefactor; replicate=replicate, seed=seed)   # generate edge list

    println("Starting Kernel 1.")

    g, kernel_1_time = @timed kernel_1(ij)    # kernel1

    T = eltype(g)
    search_key = key_sampling(g)    # sample min(nv(g),64) vertices for dfs
    nbfs  = length(search_key)      # number of bfs
    kernel_2_time =  zeros(Float64, nbfs)   # stores per bfs time
    kernel_2_nedge = zeros(Float64, nbfs)   # stores number of edges traversed per bfs

    ### kernel_2
    @showprogress 1 "Kernel 2: " for (k, key) in enumerate(search_key)
        parent_, kernel_2_time[k] = @timed kernel_2(g, key)
        err = validate(g, parent_, ij, key)
        err <= 0 && error("BFS ", k, " from search key ", key, " failed to validate: ", err)

        ### Volume/2
        total_edges_traversed = 0
        for i in 1:length(parent_)
            if(parent_[i] > 0)
                total_edges_traversed += indegree(g, parent_[i])
            end
        end
        kernel_2_nedge[k] = Float64(total_edges_traversed)/2
    end

    ### Result statistics
    result_state = Graph500Results(scale, nbfs, kernel_1_time, kernel_2_time, kernel_2_nedge)
    result = Dict{Graph{T}, Graph500Results}()
    result[g] = result_state
    output(result_state, nv(g), ne(g))
    return result
end


text:
include("params.jl")
include("mb.jl")

struct EfficientNet{S, B, H, P, F}
  stem::S
  blocks::B

  head::H
  pooling::P
  top::F
end
Flux.@functor EfficientNet

function EfficientNet(
  model_name, block_params, global_params; in_channels, n_classes, pretrain,
)
  pad, bias = SamePad(), false
  out_channels = round_filter(32, global_params)
  stem = Chain(
    Conv((3, 3), in_channels=>out_channels; bias, stride=2, pad),
    BatchNorm(out_channels, swish))

  blocks = MBConv[]
  for bp in block_params
    in_channels = round_filter(bp.in_channels, global_params)
    out_channels = round_filter(bp.out_channels, global_params)
    repeat = global_params.depth_coef ≈ 1 ?
      bp.repeat : ceil(Int64, global_params.depth_coef * bp.repeat)

    push!(blocks, MBConv(
      in_channels, out_channels, bp.kernel, bp.stride;
      expansion_ratio=bp.expansion_ratio))
    for _ in 1:(repeat - 1)
      push!(blocks, MBConv(
        out_channels, out_channels, bp.kernel, 1;
        expansion_ratio=bp.expansion_ratio))
    end
  end
  blocks = Chain(blocks...)

  head_out_channels = round_filter(1280, global_params)
  head = Chain(
    Conv((1, 1), out_channels=>head_out_channels; bias, pad),
    BatchNorm(head_out_channels, swish))

  top = Dense(head_out_channels, n_classes)
  model = EfficientNet(stem, blocks, head, AdaptiveMeanPool((1, 1)), top)
  pretrain && loadpretrain!(model, "EfficientNet" * model_name)
  model
end

"""
    EfficientNet(block_params, global_params; in_channels = 3)

Construct an EfficientNet model
([reference](https://arxiv.org/abs/1905.11946)).

# Arguments
- `model_name::String`: Name of the model. Accepts `b0`-`b8` names.
- `in_channels::Int`: Number of input channels. Default is `3`.
- `n_classes::Int`: Number of output classes. Default is `1000`.
- `pretrain::Bool`: Whether to load ImageNet pretrained weights.
  Default is `false`.
"""
EfficientNet(
  model_name::String; in_channels::Int = 3,
  n_classes::Int = 1000, pretrain::Bool = false,
) = EfficientNet(
  model_name, get_efficientnet_params(model_name)...;
  in_channels, n_classes, pretrain)

(m::EfficientNet)(x) = m.top(Flux.flatten(m.pooling(m.head(m.blocks(m.stem(x))))))


text:
module Spinach
using Reexport
import LLVM
# Write your package code here.
    @reexport import LLVM
    include("create_ast.jl")
    include("codegen.jl")
    export generate_IR, global_machine, run_module
end


text:
"""
    functional_cartography(N::T, L::Dict{E,Int64}) where {T<:BinaryNetwork, E<:AllowedSpeciesTypes}

This function will take the output of a modularity analysis (*i.e.* a network
and a partition), and return a dictionary where every species is associated to
its functional role, as defined in Olesen et al (2005). The first element is the
within-module degree z-score, and the second is the participation coefficient.
"""
function functional_cartography(N::T, L::Dict{E,Int64}) where {T<:BinaryNetwork, E<:AllowedSpeciesTypes}
    κ = Dict{last(eltype(N)), Int64}()
    S = Dict{Int64,Tuple{Float64, Float64}}()
    carto = Dict{last(eltype(N)),Tuple{Float64, Float64}}()

    for s in species(N)
        sp_module = L[s]
        k_i = 0
        if s ∈ species(N; dims=2)
            for s2 in N[:,s]
                L[s2] == sp_module && (k_i += 1)
            end
        end
        if s ∈ species(N; dims=1)
            for s2 in N[s,:]
                L[s2] == sp_module && (k_i += 1)
            end
        end
        κ[s] = k_i
    end

    for s in values(L)
        sp = filter(k -> k.second == s, L)
        k_i = [get(κ, x, 0) for x in keys(sp)]
        S[s] = (mean(k_i), std(k_i))
    end

    # Participation coefficient
    for s in species(N)
        ks = zeros(Int64, length(unique(values(L))))
        if s ∈ species(N; dims=2)
            for s2 in N[:,s]
                ks[L[s2]] += 1
            end
        end
        if s ∈ species(N; dims=1)
            for s2 in N[s,:]
                ks[L[s2]] += 1
            end
        end
        carto[s] = ((κ[s] - S[L[s]][1])/S[L[s]][2], 1-sum((ks./degree(N)[s]).^2))
    end

    return carto
end


text:
# using Dates
# using DelimitedFiles
# using CSV
# using Plots
using Dates

#calculates the potential Evaporation based on the equation by Hargreaves-Samani

function epot_hargreaves(Temp_min::Float64, Temp::Float64, Temp_max::Float64, KT::Float64, Radiation::Float64)
    #Radiation in mm/day
    # KT 0.162 for interior regions where land mass dominates, and
    # 0.190 for coastal regions, where air masses are influenced by a nearby water body
    # 0.17 for Salt Lake City
    Rs = KT * ((Temp_max - Temp_min)^0.5) * Radiation #global solar radiation in mm/day
    Epot = 0.0135 * Rs * (Temp + 17.8)  #potential Evaporation in mm/day
    #Epot = 0.0135 * Rs * ((Temp_max + Temp_min)/2 + 17.8)
    return Epot, Rs

end

function radiation(Latitude::Float64, Day::Int)
    # formula derived from Crop evapotranspiration-Guidelines for computing crop water requirements- FAO Irrigation and drainage paper 56
    # Latitude should be in the form 13°44' = 13+44/66
    #Latitude = floor(Latitude) + mod(Latitude,1) * 10 /60
    Solar_Constant = 0.0820 #MJ/m2/min
    Distance_Earth_Sun = 1 + 0.033 * cos(2 * pi / 365 * Day)
    Latitude = Latitude * pi / 180  #to convert degree to rad
    Solar_Declination = 0.409 * sin((2 * pi/ 365) * Day -1.39) # in rad
    Sunsethour_angle = acos(-tan(Latitude) * tan(Solar_Declination)) # in rad
    Radiation = 24 * 60 / pi * Solar_Constant * Distance_Earth_Sun * (Sunsethour_angle * sin(Latitude) * sin(Solar_Declination) + cos(Latitude) * cos(Solar_Declination) * sin(Sunsethour_angle))
    Radiation = 0.408 * Radiation #to convert it from MJ/m2/day to mm/day
    return Radiation::Float64
end

function getEpot(Temp_min::Array{Float64, 1}, Temp::Array{Float64, 1}, Temp_max::Array{Float64, 1}, KT::Float64, Timeseries::Array{Date, 1}, Latitude::Float64)
    Evaporation = Float64[]
    Radiation = Float64[]
    for (t, current_date) in enumerate(Timeseries)
        #current_date = Date(current_date, dateformat"y,m,d")
        Day = Dates.dayofyear(current_date)
        Current_Radiation = radiation(Latitude, Day)
        Current_Evaporation, Current_Rs = epot_hargreaves(Temp_min[t], Temp[t], Temp_max[t], KT, Current_Radiation)
        push!(Evaporation, Current_Evaporation)
        push!(Radiation, Current_Radiation)
    end
    return Evaporation::Array{Float64,1}, Radiation::Array{Float64,1}
end

# data from https://www.worlddata.info/europe/austria/sunset.php

# Thornthwaite Method
# the monthly mean temperature needs to be calculated
function monthlytemp(Timeseries, Temp::Array{Float64, 1})
    #calculates the mean monthly temperatures of months in the timeseries
    @assert length(Timeseries[:,1]) == length(Temp)
    sum_Temp = 0
    monthly_Temp = Float64[]
    for (i, current_date) in enumerate(Timeseries)
        sum_Temp += Temp[i]
        #current_date = Date(current_date, dateformat"y,m,d")
        day = Dates.day(current_date)
        if day == Dates.daysinmonth(current_date)
            monthly_Temp_current = sum_Temp / Dates.daysinmonth(current_date)
            push!(monthly_Temp, monthly_Temp_current)
            sum_Temp = 0
        end
    end
    #print(length(monthly_Temp))
    return monthly_Temp::Array{Float64, 1}
end

function heatindex(monthly_Temp::Array{Float64, 1})
    #calculates the annual heat index
    @assert mod(length(monthly_Temp), 12) == 0
    Heatindex = Float64[]
    sum = 0
    for (index, current_Temp) in enumerate(monthly_Temp)
        I = Complex(current_Temp/ 5).^1.514
        I = real.(I)
        @assert imag.(I) == 0im
        sum += I
        #at the end of each year the annual heat index is given as an output
        if index > 2 && mod(index,12)  == 0
            push!(Heatindex, sum)
            sum = 0
        end
    end
    return Heatindex::Array{Float64, 1}
end

function epot_thornthwaite(Temp_daily::Float64, Heatindex::Float64, daysinmonth::Int64, sunhours::Float64)
    #calculates the potential evaporation by Thornthwaite method
    # Thornthwaite method assumes no evaporation below zero degrees (COMPUTER PROGRAM FOR ESTIMATING EVAPOTRANSPIRATION USING THE THORNTHWAITE METHOD)
    a = 675 * 10^(-9) * Heatindex^3 - 771 * 10^(-7) * Heatindex^2 + 1792 * 10^(-5) * Heatindex + 0.49239
    if Temp_daily < 0
        Temp_daily = 0
    end
    E_gross = 16 * (10 * Temp_daily / Heatindex)^a
    Epot = E_gross * (sunhours * daysinmonth / 360)
    Epot_daily = Epot/daysinmonth
    # if 16 * (10 * Temp_daily / Heatindex) <= 0
    #     Epot_daily = 0
    # else
    #     E_gross = 16 * (10 * Temp_daily / Heatindex)^a
    #     Epot = E_gross * (sunhours * daysinmonth / 360)
    #     Epot_daily = Epot/daysinmonth
    # end
    return Epot_daily
end

function getEpot_thornthwaite(Temp::Array{Float64, 1}, Timeseries::Array{Date, 1}, sunhours::Array{Float64, 1})
    # assertion for that the timeseries contain whole years
    @assert length(Timeseries) == length(Temp)
    Evaporation = Float64[]
    # calculate the mean monthly temperatures of the timeseries
    Temp_month = monthlytemp(Timeseries, Temp)
    #calculate the annual heat indices of the timeseries
    annual_heatindex = heatindex(Temp_month)
    first_year = Dates.year(Timeseries[1])
    for (i, current_date) in enumerate(Timeseries)
        #current_date = Date(current_date, dateformat"y,m,d")
        #get the number of days of the current month
        daysinmonth = Dates.daysinmonth(current_date)
        #print(daysinmonth)
        datetuple = Dates.yearmonthday(current_date)
        #print(datetuple[2])
        #get the sunhours of the current month
        sunhours_current = sunhours[datetuple[2]]
        # get the monthly evaporation for current month in the timeseries
        j = 12 * abs(datetuple[1]- first_year) + datetuple[2]
        Temp_monthly_current = Temp_month[j]
        # get the annual heat index of the current year
        current_annual_heatindex = annual_heatindex[abs(datetuple[1]-first_year) + 1]
        Epot = epot_thornthwaite(Temp_monthly_current, current_annual_heatindex, daysinmonth, sunhours_current)
        #Epot = Epot/daysinmonth
        push!(Evaporation, Epot)
    end
    return Evaporation::Array{Float64,1}
end

function getEpot_Daily_thornthwaite(Temp::Array{Float64, 1}, Timeseries::Array{Date, 1}, sunhours::Array{Float64, 1})
    # assertion for that the timeseries contain whole years
    @assert length(Timeseries) == length(Temp)
    Evaporation = Float64[]
    #calculate the annual heat indices of the timeseries
    Temp_month = monthlytemp(Timeseries, Temp)
    annual_heatindex = heatindex(Temp_month)
    first_year = Dates.year(Timeseries[1])
    for (i, current_date) in enumerate(Timeseries)
        #current_date = Date(current_date, dateformat"y,m,d")
        #get the number of days of the current month
        daysinmonth = Dates.daysinmonth(current_date)
        #print(daysinmonth)
        datetuple = Dates.yearmonthday(current_date)
        #print(datetuple[2])
        #get the sunhours of the current month
        sunhours_current = sunhours[datetuple[2]]
        # get the annual heat index of the current year
        current_annual_heatindex = annual_heatindex[abs(datetuple[1]-first_year) + 1]
        Epot = epot_thornthwaite(Temp[i], current_annual_heatindex, daysinmonth, sunhours_current)
        #Epot = Epot/daysinmonth
        push!(Evaporation, Epot)
    end
    return Evaporation::Array{Float64,1}
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
module EverestRcon

loaded = true
try
    using HTTP, YAML

catch e
    # HTTP not installed
    # Causes problems on Windows 7 due to github TLS versions

    global loaded = false
end

if loaded
    queryurl(url, query) = url * HTTP.escapeuri(query)

    function request(url::String, query::Dict{String, Any}=Dict{String, Any}())
        try
            res = HTTP.get(queryurl(url, query), readtimeout=1)

            return true, String(res.body)

        catch e
            if isa(e, HTTP.ExceptionRequest.StatusError)
                return false, String(e.response.body)

            else
                println(Base.stderr, e)
                return false, "Connection timed out"
            end
        end
    end

    respawn(baseUrl::String) = request("$baseUrl/respawn?")

    focus(baseUrl::String) = request("$baseUrl/focus")

    reload(baseUrl::String, area::String) = request("$baseUrl/reloadmap?", Dict{String, Any}("area" => area))
    reload(baseUrl::String, area::String, side::String) = request("$baseUrl/reloadmap?", Dict{String, Any}("area" => area, "side" => side))

    # Ahorn doesn't support sides yet, ignore for now
    teleport(baseUrl::String, area::String, room::String; force::Bool=false) = request("$baseUrl/tp?", Dict{String, Any}("area" => area, "room" => room, "forcenew" => force))
    teleport(baseUrl::String, area::String, room::String, x::Integer, y::Integer; force::Bool=false) = request("$baseUrl/tp?", Dict{String, Any}("area" => area, "room" => room, "forcenew" => force, "x" => x, "y" => y))

    teleportToRoom(baseUrl::String, level::String; force::Bool=false) = request("$baseUrl/tp?", Dict{String, Any}("level" => level, "forcenew" => force))
    teleportToRoom(baseUrl::String, level::String, x::Integer, y::Integer; force::Bool=false) = request("$baseUrl/tp?", Dict{String, Any}("level" => level, "forcenew" => force, "x" => x, "y" => y))

    function session(baseUrl::String)
        success, data = request("$baseUrl/session")

        if success
            return true, YAML.load(data)

        else
            return false, data
        end
    end
end

end


text:
"""
    testgeometry(geom)

Test whether the required interface for your `geom` has been implemented correctly.
"""
function testgeometry(geom)
    @assert isgeometry(geom) "$geom doesn't implement `isgeometry`."
    type = geomtrait(geom)
    @assert !isnothing(type) "$geom doesn't implement `geomtrait`."

    if type == PointTrait()
        n = ncoord(geom)
        if n >= 1  # point could be empty
            getcoord(geom, 1)  # point always needs at least 2
        end
    else
        n = ngeom(geom)
        if n >= 1  # geometry could be empty
            g2 = getgeom(geom, 1)
            subtype = subtrait(type)
            if !isnothing(subtype)
                issub = geomtrait(g2) isa subtype
                !issub && error("Implemented hierarchy for $geom type is incorrect. Subgeometry should be a $subtype")
            end
            @assert testgeometry(g2) "Subgeometry implementation is not valid."
        end
    end
    return true
end

"""
    testfeature(feature)

Test whether the required interface for your `feature` has been implemented correctly.
"""
function testfeature(feature)
    @assert isfeature(typeof(feature)) "$feature doesn't implement `isfeature`."
    @assert trait(feature) isa AbstractFeatureTrait "$feature does not return an `AbstractFeatureTrait` for `geomtrait`."
    @assert geomtrait(feature) == nothing
    geom = geometry(feature)
    if !isnothing(geom)
        @assert isgeometry(geom) "geom $geom from $feature doesn't implement `isgeometry`."
        @assert coordinates(feature) == coordinates(geometry(feature))
    end

    props = properties(feature)
    if !isnothing(props)
        @assert first(propertynames(props)) isa Symbol "`propertynames` of $props does not return an iterable of `Symbol`"
        map(n -> getproperty(props, n), propertynames(props))
    end
    ext = extent(feature)
    @assert ext isa Union{Nothing,Extent}
    return true
end

"""
    testfeaturecollection(featurecollection)

Test whether the required interface for your `featurecollection` has been implemented correctly.
"""
function testfeaturecollection(fc)
    @assert isfeaturecollection(typeof(fc)) "$fc doesn't implement `isfeaturecollection`."
    @assert trait(fc) isa AbstractFeatureCollectionTrait "$fc does not return an `AbstractFeatureCollectionTrait` for `geomtrait`."
    @assert geomtrait(fc) == nothing
    @assert isa(nfeature(fc), Integer) "feature collection $fc doesn't return an `Integer` from `nfeatures`."
    if nfeature(fc) > 0
        @assert isfeature(getfeature(fc, 1)) "For $fc `getfeature(featurecollection, 1)` does not return an object where `isfeature(obj) == true`."
        @assert isfeature(getfeature(fc, nfeature(fc))) "For $fc `getfeature(featurecollection, nfeatures(featurecollection))` does not return an object where `isfeature(obj) == true`."
    else
        @warn "`nfeatures == 0` for feature collection, cannot test some properties"
    end
    @assert coordinates(fc) == coordinates.(getfeature(fc))
    return true
end


text:
type Module
	speed::Float64
	energy::Float64
	production::Float64
end

function stats(pse)
	if pse[1] > 0
		return [[0.04 -0.15 0.05]; [0.06 -0.15 0.6]; [0.1 -0.15 0.8]][pse[1],:]
	end
	if pse[2] > 0
		return [[0 0.2 0.5]; [0 0.3 0.6]; [0 0.5 0.7]][pse[2],:]
	end
	if pse[3] > 0
		return [[0 0 -0.3]; [0 0 -0.4]; [0 0 -0.5]][pse[3],:]
	end
end

type Assembler
	name::AbstractString
	kw::Int64
	speed::Float64
	slots::Int64
	#slots::Vector{Module}
end

function stats(a::Assembler, pse)
	([1 a.speed a.kw]' * [1+pse[1] 1+pse[2] 1+pse[3]])[1,:]
end

function string(a::Assembler, prods, speeds, effs)

	#@sprintf "%s : %0.2fkW %0.2f speed %0.2f J per Work" a.name
end

macro taken()
	:(sum([sum(prod) sum(speed) sum(eff)]))
end

macro rem()
	:(0:max(0, a.slots-@taken))
end

prod = [0 0 0]
speed = [0 0 0]
eff = [0 0 0]
for a in [Assembler("A1", 90, 0.5, 0) Assembler("A2", 150, 0.75, 2) Assembler("A3", 210, 1.25, 4)]
	for prod[1] in @rem
		for prod[2] in @rem
			for prod[3] in @rem
				for speed[1] in @rem
					for speed[2] in @rem
						for speed[3] in @rem
							for eff[1] in @rem
								for eff[2] in @rem
									for eff[3] in @rem
										@printf "%s - " a.name
										stat = [0 0 0]
										for p in 1:3
											if prod[p] > 0
												stat = prod[p] * stats([p 0 0])
												@printf "\t%d Prod %d" prod[p] p
											end
										end
										for s in 1:3
											if speed[s] > 0
												stat = speed[s] * stats([0 s 0])
												@printf "\t%d Speed %d" speed[s] s
											end
										end
										for e in 1:3
											if eff[e] > 0
												stat = eff[e] * stats([0 0 e])
												@printf "\t%d Eff %d" eff[e] e
											end
										end
										@printf "\t%s\n" stats(a, stat)
									end
									eff[3] = 0
								end
								eff[2] = 0
							end
							eff[1] = 0
						end
						speed[3] = 0
					end
					speed[2] = 0
				end
				speed[1] = 0
			end
			prod[3] = 0
		end
		prod[2] = 0
	end
	prod[1] = 0
end

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
# julia ncurses interface
# TODO:
# - examine the option to use Clang.jl's C wrapper generator
#   * so far I can spot a few reasons why this option wouldn't work

module NCurses

using Base
using BinDeps
using Compat

# The following loads "../deps/build.jl" and forces re-install of library
# this explains why other modules prefer to load deps.jl instead.
# @BinDeps.load_dependencies "../deps/deps.jl" [:ncurses] # not working as expected

if isfile(joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl"))
  include("../deps/deps.jl")
else
  error("NCurses not properly installed. Please run Pkg.build(\"NCurses\")")
end

# think if it is needed to do this, or should we leave it to each application
# follow-up : no we certainly do not want to do this init here, as this
# returns the screns' main pointer which should be available to the application.
# Now, thinking about moving this to a module singleton object.
function __init__()
    err = ccall((:initscr, ncurses), Ptr{Void}, ())
    err != C_NULL || error("error initializing NCurses module")
    atexit() do
        ccall((:endwin, ncurses), Ptr{Void}, ())
    end
end

function version()
    nc_ver = ccall((:curses_version, ncurses) , Ptr{UInt8}, ())
    if nc_ver != C_NULL
        return bytestring(nc_ver)
      else
        return ""
    end
end

include("common.jl")
include("curses.jl")
include("text.jl")
include("window.jl")
include("screen.jl")
include("form.jl")
include("menu.jl")
include("panel.jl")
include("event.jl")
include("other.jl")

end #module


text:
module NumericalFluxes

export Rusanov, CentralGradPenalty, CentralNumericalFluxDiffusive,
       CentralNumericalFluxNonDiffusive

using StaticArrays, LinearAlgebra
using GPUifyLoops: @unroll
using CLIMA.VariableTemplates
import ..DGmethods: BalanceLaw, Grad, Vars, vars_state, vars_diffusive,
                    vars_aux, vars_gradient, boundary_state!, wavespeed,
                    flux_nondiffusive!, flux_diffusive!, diffusive!, num_state,
                    num_gradient, gradvariables!

"""
    GradNumericalPenalty

Any `P <: GradNumericalPenalty` should define methods for:

   diffusive_penalty!(gnf::P, bl::BalanceLaw, diffF, nM, QM, QdiffM, QauxM, QP,
                      QdiffP, QauxP, t)
   diffusive_boundary_penalty!(gnf::P, bl::BalanceLaw, l_Qvisc, nM, l_GM, l_QM,
                               l_auxM, l_GP, l_QP, l_auxP, bctype, t)

"""
abstract type GradNumericalPenalty end

function diffusive_penalty! end
function diffusive_boundary_penalty! end

"""
    CentralGradPenalty <: GradNumericalPenalty

"""
struct CentralGradPenalty <: GradNumericalPenalty end

function diffusive_penalty!(::CentralGradPenalty, bl::BalanceLaw,
    diff_penalty::Vars{D}, n::SVector,
    transform⁻::Vars{T}, state⁻::Vars{S}, aux⁻::Vars{A},
    transform⁺::Vars{T}, state⁺::Vars{S}, aux⁺::Vars{A},
    t) where {D,T,S,A}

  G = n .* (parent(transform⁺) .- parent(transform⁻))' ./ 2
  diffusive!(bl, diff_penalty, Grad{T}(G), state⁻, aux⁻, t)
end

function diffusive_boundary_penalty!(nf::CentralGradPenalty, bl::BalanceLaw,
    diff_penalty::Vars{D}, n::SVector,
    transform⁻::Vars{T}, state⁻::Vars{S}, aux⁻::Vars{A},
    transform⁺::Vars{T}, state⁺::Vars{S}, aux⁺::Vars{A},
    bctype, t, state1⁻::Vars{S}, aux1⁻::Vars{A}) where {D,T,S,A}

  boundary_state!(nf, bl, state⁺, aux⁺, n, state⁻, aux⁻,
    bctype, t, state1⁻, aux1⁻)

  gradvariables!(bl, transform⁺, state⁺, aux⁺, t)

  diffusive_penalty!(nf, bl, diff_penalty, n,
    transform⁻, state⁻, aux⁻,
    transform⁺, state⁺, aux⁺,
    t)
end

"""
    NumericalFluxNonDiffusive

Any `N <: NumericalFluxNonDiffusive` should define the a method for

    numerical_flux_nondiffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QauxM, QP,
                                 QauxP, t)

where
- `F` is the numerical flux array
- `nM` is the unit normal
- `QM`/`QP` are the minus/positive state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_nondiffusive!(nf::N, bl::BalanceLaw, F, nM, QM,
                                          QauxM, QP, QauxP, bctype, t)

"""
abstract type NumericalFluxNonDiffusive end

function numerical_flux_nondiffusive! end

function numerical_boundary_flux_nondiffusive!(nf::NumericalFluxNonDiffusive,
    bl::BalanceLaw, fluxᵀn::Vars{S}, n::SVector,
    state⁻::Vars{S}, aux⁻::Vars{A}, state⁺::Vars{S}, aux⁺::Vars{A},
    bctype, t, state1⁻::Vars{S}, aux1⁻::Vars{A}) where {S,A}

  boundary_state!(nf, bl, state⁺, aux⁺, n,
    state⁻, aux⁻, bctype, t, state1⁻, aux1⁻)

  numerical_flux_nondiffusive!(nf, bl, fluxᵀn, n, state⁻, aux⁻, state⁺, aux⁺, t)
end

"""
    Rusanov <: NumericalFluxNonDiffusive

The Rusanov (aka local Lax-Friedrichs) numerical flux.

# Usage

    Rusanov()

Requires a `flux_nondiffusive!` and `wavespeed` method for the balance law.
"""
struct Rusanov <: NumericalFluxNonDiffusive end

update_penalty!(::Rusanov, ::BalanceLaw, _...) = nothing

function numerical_flux_nondiffusive!(nf::Rusanov,
    bl::BalanceLaw, fluxᵀn::Vars{S}, n::SVector,
    state⁻::Vars{S}, aux⁻::Vars{A}, state⁺::Vars{S}, aux⁺::Vars{A}, t) where {S,A}

  numerical_flux_nondiffusive!(CentralNumericalFluxNonDiffusive(),
    bl, fluxᵀn, n, state⁻, aux⁻, state⁺, aux⁺, t)

  Fᵀn = parent(fluxᵀn)
  λ⁻ = wavespeed(bl, n, state⁻, aux⁻, t)
  λ⁺ = wavespeed(bl, n, state⁺, aux⁺, t)
  λ = max(λ⁻, λ⁺)
  λΔQ = λ * (parent(state⁻) - parent(state⁺))

  # TODO: should this operate on ΔQ or λΔQ?
  update_penalty!(nf, bl, n, λ,
    Vars{S}(λΔQ), state⁻, aux⁻, state⁺, aux⁺, t)

  Fᵀn .+= λΔQ/2
end

"""
    CentralNumericalFluxNonDiffusive() <: NumericalFluxNonDiffusive

The central numerical flux for nondiffusive terms

# Usage

    CentralNumericalFluxNonDiffusive()

Requires a `flux_nondiffusive!` method for the balance law.
"""
struct CentralNumericalFluxNonDiffusive <: NumericalFluxNonDiffusive end

function numerical_flux_nondiffusive!(::CentralNumericalFluxNonDiffusive,
    bl::BalanceLaw, fluxᵀn::Vars{S}, n::SVector,
    state⁻::Vars{S}, aux⁻::Vars{A}, state⁺::Vars{S}, aux⁺::Vars{A}, t) where {S,A}

  FT = eltype(fluxᵀn)
  nstate = num_state(bl,FT)
  Fᵀn = parent(fluxᵀn)

  F⁻ = similar(Fᵀn, Size(3, nstate))
  fill!(F⁻, -zero(FT))
  flux_nondiffusive!(bl, Grad{S}(F⁻), state⁻, aux⁻, t)

  F⁺ = similar(Fᵀn, Size(3, nstate))
  fill!(F⁺, -zero(FT))
  flux_nondiffusive!(bl, Grad{S}(F⁺), state⁺, aux⁺, t)

  Fᵀn .+= (F⁻ + F⁺)' * (n/2)
end

"""
    NumericalFluxDiffusive

Any `N <: NumericalFluxDiffusive` should define the a method for

    numerical_flux_diffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QdiffM, QauxM, QP,
                              QdiffP, QauxP, t)

where
- `F` is the numerical flux array
- `nM` is the unit normal
- `QM`/`QP` are the minus/positive state arrays
- `QdiffM`/`QdiffP` are the minus/positive diffusive state arrays
- `QdiffM`/`QdiffP` are the minus/positive auxiliary state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_diffusive!(nf::N, bl::BalanceLaw, F, nM, QM, QdiffM,
                                       QauxM, QP, QdiffP, QauxP, bctype, t)

"""
abstract type NumericalFluxDiffusive end

function numerical_flux_diffusive! end

function numerical_boundary_flux_diffusive!(nf::NumericalFluxDiffusive,
    bl::BalanceLaw, fluxᵀn::Vars{S}, n::SVector,
    state⁻::Vars{S}, diff⁻::Vars{D}, aux⁻::Vars{A},
    state⁺::Vars{S}, diff⁺::Vars{D}, aux⁺::Vars{A},
    bctype, t,
    state1⁻::Vars{S}, diff1⁻::Vars{D}, aux1⁻::Vars{A}) where {S,D,A}

  boundary_state!(nf, bl, state⁺, diff⁺, aux⁺,
    n, state⁻, diff⁻, aux⁻, bctype, t,
    state1⁻, diff1⁻, aux1⁻)

  numerical_flux_diffusive!(nf, bl, fluxᵀn, n,
    state⁻, diff⁻, aux⁻, state⁺, diff⁺, aux⁺, t)
end

"""
    CentralNumericalFluxDiffusive <: NumericalFluxDiffusive

The central numerical flux for diffusive terms

# Usage

    CentralNumericalFluxDiffusive()

Requires a `flux_diffusive!` for the balance law.
"""
struct CentralNumericalFluxDiffusive <: NumericalFluxDiffusive end

function numerical_flux_diffusive!(::CentralNumericalFluxDiffusive,
  bl::BalanceLaw, fluxᵀn::Vars{S}, n::SVector,
  state⁻::Vars{S}, diff⁻::Vars{D}, aux⁻::Vars{A},
  state⁺::Vars{S}, diff⁺::Vars{D}, aux⁺::Vars{A}, t) where {S,D,A}

  FT = eltype(fluxᵀn)
  nstate = num_state(bl,FT)
  Fᵀn = parent(fluxᵀn)

  F⁻ = similar(Fᵀn, Size(3, nstate))
  fill!(F⁻, -zero(FT))
  flux_diffusive!(bl, Grad{S}(F⁻), state⁻, diff⁻, aux⁻, t)

  F⁺ = similar(Fᵀn, Size(3, nstate))
  fill!(F⁺, -zero(FT))
  flux_diffusive!(bl, Grad{S}(F⁺), state⁺, diff⁺, aux⁺, t)

  Fᵀn .+= (F⁻ + F⁺)' * (n/2)
end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Result (ModiaMath/Result/_module.jl)
#

"""
    mutable struct RawResult

Hold result data in a Float64 matrix (first column time, i-th column variable i).
Typical usage:

```
    raw = RawResult(100,3) # Initial storage: 100 time points, 3 variables
    storeRawResult!(raw, [1.0, 2.0, 1.0])
    storeRawResult!(raw, [2.0, 3.0, 5.0])
      ...
    res = getDictResult(res, ["time", "r[1]", "r[2]"]

    plot(res["time"], res["r[1]"])
```

If initial storage is not sufficient, it is automatically doubled.
"""
mutable struct RawResult
    nt::Int
    data::Matrix{Float64}
    names::Vector{String}

    function RawResult(nt::Int, names::Vector{String})
        @assert(nt > 0)
        @assert(length(names) > 0)
        nv = length(names)
        new(0, zeros(nt, nv), names)
    end
end

nResults(res::RawResult) = res.nt

function storeRawResult!(res::RawResult, v::Vector{Float64})
    res.nt += 1
    i = res.nt
    ntMax = size(res.data, 1)

    if i > ntMax
        # Result storage is full, double the current storage
        res.data = [res.data; zeros(ntMax, size(res.data, 2))]
    end

    for j = 1:length(v)
        res.data[i,j] = v[j]
    end

    return nothing
end

function storeRawResult!(res::RawResult, model, t::Float64, x::Vector{Float64}, derx::Vector{Float64}, w::Vector{Float64})
    @assert(length(x) == length(derx))
    res.nt += 1
    i = res.nt
    ntMax = size(res.data, 1)

    if i > ntMax
        # Result storage is full, double the current storage
        res.data = [res.data; zeros(ntMax, size(res.data, 2))]
    end

    # Copy data in raw result storage
    res.data[i,1] = t
    jstart = 1
    nx = length(x)

    for j = 1:nx
        res.data[i,jstart + j] = x[j]
    end

    jstart = jstart + nx
    for j = 1:nx
        res.data[i,jstart + j] = derx[j]
    end

    jstart = jstart + nx
    for j = 1:length(w)
        res.data[i,jstart + j] = w[j]
    end

    return nothing
end


text:
"""
Module Help

Provides functions for displaying help about the commands during speech recognition.

# Functions
- [`help`](@ref)

To see a description of a function type `?<functionname>`.
"""
module Help

import ..JustSayIt: command, command_names, next_token, pretty_cmd_string, PyKey

const COMMANDS_KEYWORD = "commands"

"Show help for your commands or a spoken command or module."
function help()
    valid_input = [COMMANDS_KEYWORD, command_names()...]
    keyword = next_token(valid_input)
    if keyword == COMMANDS_KEYWORD
        cmd_length_max = maximum(length.(command_names()))
        @info join(["", "Your commands:",
                    map(sort([command_names()...])) do x
                       join((x, pretty_cmd_string(command(x))), " "^(cmd_length_max+1-length(x)) * "=> ")
                    end...
                    ], "\n")
    elseif keyword in command_names()
        cmd = command(keyword)
        if isa(cmd, Function)
            @info "Command $keyword" ""=Base.Docs.doc(cmd)
        elseif isa(cmd, PyKey)
            @info "Command $keyword\n   =   Keyboard key $(pretty_cmd_string(cmd))"
        else
            @info "Command $keyword\n   =   Keyboard shortcut $(pretty_cmd_string(cmd))"
        end
    else
        @info "Keyword not recognized."
    end
    return
end

end # module Help


text:
# `clip-01`

using StatisticalRethinking
#gr(size=(600,600));

ProjDir = @__DIR__

# ### snippet 9.1

# Metropolis algorithm

function generate_walk(N::Int64)
  num_weeks = N
  positions = zeros(Int64, num_weeks);
  current = 10
  d = Uniform(0, 1)

  for i in 1:N
    positions[i] = current  # Record current position
    proposal = current + sample([-1, 1], 1)[1] # Generate proposal
    proposal = proposal < 1  ? 10 : proposal
    proposal = proposal > 10  ? 1 : proposal
    prob_move = proposal/current  # Move?
    current = rand(d, 1)[1] <  prob_move ? proposal : current
  end

  positions

end

# Generate the visits

N = 100000
walk = generate_walk(N)

# Plot the first 100 weeks and a histogram of weeks per iosland

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = plot(walk[1:100], leg=false, xlabel="Week", ylabel="Island", title="First 100 steps")
p[2] = histogram(walk, leg=false, xlabel="Island", ylabel="Number of weeks",
  title="$N steps")
plot(p..., layout=(1, 2))

# End of `clip-01.jl`