text:
using MPI
using StaticArrays
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging

using ClimateMachine.BalanceLaws:
    BalanceLaw,
    Prognostic,
    Auxiliary,
    GradientFlux,
    UpwardIntegrals,
    DownwardIntegrals

import ClimateMachine.BalanceLaws:
    vars_state,
    integral_load_auxiliary_state!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    update_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    boundary_conditions,
    boundary_state!,
    compute_gradient_argument!,
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ClimateMachine.DGMethods: init_ode_state
using ClimateMachine.Mesh.Geometry: LocalGeometry

struct IntegralTestSphereModel{T} <: BalanceLaw
    Rinner::T
    Router::T
end

function update_auxiliary_state!(
    dg::DGModel,
    m::IntegralTestSphereModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)
    reverse_indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)

    return true
end

vars_state(::IntegralTestSphereModel, ::UpwardIntegrals, T) = @vars(v::T, r::T)
vars_state(::IntegralTestSphereModel, ::DownwardIntegrals, T) =
    @vars(v::T, r::T)
vars_state(m::IntegralTestSphereModel, ::Auxiliary, T) = @vars(
    int::vars_state(m, UpwardIntegrals(), T),
    rev_int::vars_state(m, UpwardIntegrals(), T),
    r::T,
    a::T
)

vars_state(::IntegralTestSphereModel, ::Prognostic, T) = @vars()
vars_state(::IntegralTestSphereModel, ::GradientFlux, T) = @vars()

flux_first_order!(::IntegralTestSphereModel, _...) = nothing
flux_second_order!(::IntegralTestSphereModel, _...) = nothing
source!(::IntegralTestSphereModel, _...) = nothing
boundary_conditions(::IntegralTestSphereModel) = (nothing,)
boundary_state!(_, ::Nothing, ::IntegralTestSphereModel, _...) = nothing
init_state_prognostic!(::IntegralTestSphereModel, _...) = nothing
wavespeed(::IntegralTestSphereModel, _...) = 1

function nodal_init_state_auxiliary!(
    m::IntegralTestSphereModel,
    aux::Vars,
    tmp::Vars,
    g::LocalGeometry,
)

    x, y, z = g.coord
    aux.r = hypot(x, y, z)
    θ = atan(y, x)
    ϕ = asin(z / aux.r)
    # Exact integral
    aux.a = 1 + cos(ϕ)^2 * sin(θ)^2 + sin(ϕ)^2
    aux.int.v = exp(-aux.a * aux.r^2) - exp(-aux.a * m.Rinner^2)
    aux.int.r = aux.r - m.Rinner
    aux.rev_int.v = exp(-aux.a * m.Router^2) - exp(-aux.a * aux.r^2)
    aux.rev_int.r = m.Router - aux.r
end

@inline function integral_load_auxiliary_state!(
    m::IntegralTestSphereModel,
    integrand::Vars,
    state::Vars,
    aux::Vars,
)
    integrand.v = -2 * aux.r * aux.a * exp(-aux.a * aux.r^2)
    integrand.r = 1
end

@inline function integral_set_auxiliary_state!(
    m::IntegralTestSphereModel,
    aux::Vars,
    integral::Vars,
)
    aux.int.v = integral.v
    aux.int.r = integral.r
end

@inline function reverse_integral_load_auxiliary_state!(
    m::IntegralTestSphereModel,
    integral::Vars,
    state::Vars,
    aux::Vars,
)
    integral.v = aux.int.v
    integral.r = aux.int.r
end

@inline function reverse_integral_set_auxiliary_state!(
    m::IntegralTestSphereModel,
    aux::Vars,
    integral::Vars,
)
    aux.rev_int.v = integral.v
    aux.rev_int.r = integral.r
end

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

using Test
function test_run(mpicomm, topl, ArrayType, N, FT, Rinner, Router)
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )
    dg = DGModel(
        IntegralTestSphereModel(Rinner, Router),
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    dQdt = similar(Q)

    exact_aux = copy(dg.state_auxiliary)
    dg(dQdt, Q, nothing, 0.0)
    (int_r_ind, rev_int_r_ind) = varsindices(
        vars_state(dg.balance_law, Auxiliary(), FT),
        ("int.r", "rev_int.r"),
    )

    # We should be exact for the integral of ∫_{R_{inner}}^{r} 1
    if N != 0
        @test exact_aux[:, int_r_ind, :] ≈ dg.state_auxiliary[:, int_r_ind, :]
        @test exact_aux[:, rev_int_r_ind, :] ≈
              dg.state_auxiliary[:, rev_int_r_ind, :]
    else
        # FIXME: With the N = 0 case, the value of the integral actually will be
        # the value of the computed integral is at the cell faces, NOT the cell
        # centers, so we will be off in the direct comparison. So instead we
        # only compare the last value for now of the FORWARD integral. The
        # reverse integral and other cases will be fixed up later.
        nvertelem = topl.stacksize
        nhorzelem = div(length(topl.elems), nvertelem)
        naux = size(exact_aux, 2)
        aux = reshape(
            dg.state_auxiliary.data,
            (N + 1, N + 1, N + 1, naux, nvertelem, nhorzelem),
        )
        # Only the last value of the forward integral is currently correct
        @test all(Router - Rinner .≈ aux[:, :, end, int_r_ind, end, :])
        # FIXME: Reverse integral is currently off by one cell width
        Δ = (Router - Rinner) / nvertelem
        @test all(Router - Rinner - Δ .≈ aux[:, :, 1, rev_int_r_ind, 1, :])
        # All the `JcV` (line integral metrics) values should be `Δ / 2`
        @test all(Δ .≈ 2grid.vgeo[:, Grids._JcV, :])
    end

    euclidean_distance(exact_aux, dg.state_auxiliary)
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    base_Nhorz = 4
    base_Nvert = 2
    Rinner = 1 // 2
    Router = 1

    expected_result = Dict()
    expected_result[1] = [
        1.5934735012225074e-02
        4.0030667455285352e-03
        1.0020652111566574e-03
        2.5059856392475887e-04
    ]
    expected_result[4] = [
        4.662884229467401e-7,
        7.218989778540723e-9,
        1.1258613174916711e-10,
        1.7587739986848968e-12,
    ]
    lvls = integration_testing ? length(expected_result[4]) : 1

    for N in (0, 1, 4)
        for FT in (Float64,) # Float32)
            err = zeros(FT, lvls)
            for l in 1:lvls
                @info (ArrayType, FT, "sphere", N, l)
                Nhorz = 2^(l - 1) * base_Nhorz
                Nvert = 2^(l - 1) * base_Nvert
                Rrange = grid1d(FT(Rinner), FT(Router); nelem = Nvert)
                topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
                err[l] = test_run(
                    mpicomm,
                    topl,
                    ArrayType,
                    N,
                    FT,
                    FT(Rinner),
                    FT(Router),
                )
                if N != 0
                    @test expected_result[N][l] ≈ err[l] rtol = 1e-3 atol =
                        eps(FT)
                end
            end
            if integration_testing && N != 0
                @info begin
                    msg = "polynomialorder order $N"
                    for l in 1:(lvls - 1)
                        rate = log2(err[l]) - log2(err[l + 1])
                        msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                    end
                    msg
                end
            end
        end
    end
end

nothing


text:
resolutions(::Type{EarthEnv{HabitatHeterogeneity}}) = ("1km", "5km", "25km")
defres(::Type{EarthEnv{HabitatHeterogeneity}}) = "25km"
layers(::Type{EarthEnv{HabitatHeterogeneity}}) = (
    :cv, :evenness, :range, :shannon, :simpson, :std, :Contrast, :Correlation,
    :Dissimilarity, :Entropy, :Homogeneity, :Maximum, :Uniformity, :Variance
)

"""
    getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, [layer::Union{Tuple,Integer}]; res::Int=25) => String
    getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layer::Integer, res::Int=25) => String

Download EarthEnv habitat heterogeneity data, choosing `layer` from:
`$(layers(EarthEnv{HabitatHeterogeneity}))` and `res` from
`$(resolutions(EarthEnv{HabitatHeterogeneity}))`.

Without a layer argument, all layers will be downloaded and a tuple of paths returned.
If the data is already downloaded the path will be returned without the getraster.
"""
function getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res::String=defres(T))
    getraster(T, layer, res)
end
function getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol, res::String)
    _check_layer(T, layer)
    _check_res(T, res)
    path = rasterpath(T, layer; res)
    url = rasterurl(T, layer; res)
    return _maybe_download(url, path)
end

rastername(::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res::String=defres(T)) =
    "$(layer)_$(res).tif"
rasterpath(T::Type{<:EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res::String=defres(T)) =
    joinpath(rasterpath(T), string(res), rastername(T, layer; res))
rasterurl(T::Type{EarthEnv{HabitatHeterogeneity}}, layer; res::String=defres(T)) =
    joinpath(rasterurl(T), "$(res)/$(layer)_01_05_$(res)_$(_getprecision(layer, res)).tif")

_pathsegment(::Type{HabitatHeterogeneity}) = "habitat_heterogeneity"

function _getprecision(layer, res)
    precision = "uint16"
    # WELLCOME TO HELL
    precision = ((res == 10) && (layer == :cv)) ? "uint32" : precision
    precision = (layer == :Contrast) ? "uint32" : precision
    precision = (layer == :Dissimilarity) ? "uint32" : precision
    precision = (layer == :Variance) ? "uint32" : precision
    precision = ((res == 25) && (layer == :Entropy)) ? "uint32" : precision
    # /HELL
end

ext: jl


text:
# Autogenerated wrapper script for Imath_jll for powerpc64le-linux-gnu
export libImath

JLLWrappers.@generate_wrapper_header("Imath")
JLLWrappers.@declare_library_product(libImath, "libImath-3_1.so.29")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libImath,
        "lib/libImath-3_1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## This must come before the "module" (until we find out why)
using DeepConvert

module PrimeSieve

using Compat
using Primes

import DeepConvert: @mkdeepconvert, @bigint

include("../deps/deps.jl")
import Base: convert

if isdefined(Base, :ccall)
    import Base: ccall
end

export genprimes, nprimes
export snthprimea
export ntcountprimes
export countprimes, scountprimes

export printprimes
export printprimes2
export printprimes3
export printprimes4
export printprimes5
export printprimes6

export primesievesize, primetest, primesieve_num_threads, isprime

##

# Convert numbers to Int128 or UInt64, hopefully the subexpressions
# have not overflowed.  Eg. 10^19.
# Unquoted expressions pass through
@mkdeepconvert(conv128,Int128)
@mkdeepconvert(convu64,UInt64)
@mkdeepconvert(convint,Int64)
macro i128_str(s) conv128(s) end

#include("deepconvert.jl")
include("nextprime.jl")
include("primesieve_c.jl")
include("primecount_c.jl")
include("tables.jl")
include("wrappers.jl")
include("tuples.jl")
include("primeit.jl")
include("randprime.jl")
include("msieve.jl")
include("primesieve_c_extra.jl")
include("gaussian.jl")
include("jacobi.jl")

end # module PrimeSieve


text:
using Compose, DataFrames, Gadfly, RDatasets, Statistics
set_default_plot_size(21cm, 8cm)

salaries = dataset("car","Salaries")
salaries.Salary /= 1000.0
salaries.Discipline = ["Discipline $(x)" for x in salaries.Discipline]
df = by(salaries, [:Rank,:Discipline], :Salary=>mean, :Salary=>std)
df.ymin, df.ymax = df.Salary_mean.-df.Salary_std, df.Salary_mean.+df.Salary_std
df.label = string.(round.(Int, df.Salary_mean))

p1 = plot(df, x=:Discipline, y=:Salary_mean, color=:Rank,
    Scale.x_discrete(levels=["Discipline A", "Discipline B"]),
    ymin=:ymin, ymax=:ymax, Geom.errorbar, Stat.dodge,
    Geom.bar(position=:dodge),
    Scale.color_discrete(levels=["Prof", "AssocProf", "AsstProf"]),
    Guide.colorkey(title="", pos=[0.76w, -0.38h]),
    Theme(bar_spacing=0mm, stroke_color=c->"black")
)
p2 = plot(df, x=:Salary_mean, y=:Discipline, color=:Rank,
    Scale.y_discrete(levels=["Discipline A", "Discipline B"]),
    label=:label, Geom.label(position=:centered), Stat.dodge(position=:stack, axis=:y),
    Geom.bar(position=:stack, orientation=:horizontal),
    Guide.yticks(orientation=:vertical), Guide.ylabel(nothing)
)
hstack(p1, p2)


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Create a decision tree to recursively split encounters in the ACASX Problem. SA algorithm.
Example usage: config=configure(ACASX_SA_Tree,"normal","nvn_dasc"); acasx_sa_tree(;config...)
"""
module ACASX_SA_Tree

export configure, acasx_sa_tree

import Compat.ASCIIString
using DecisionTrees
using ExprSearch.SA
using Datasets
using RLESUtils, Obj2Dict, Configure
using Reexport

using GrammarExpts
using ACASXProblem, SA_Tree_Logs
using DerivTreeVis, DecisionTreeVis
import Configure.configure

include("dtree_callbacks.jl")

const CONFIGDIR = joinpath(dirname(@__FILE__), "..", "config")

configure(::Type{Val{:ACASX_SA_Tree}}, configs::AbstractString...) = configure_path(CONFIGDIR, configs...)

function train_dtree{T}(psa_params::PSAESParams, problem::ACASXClustering, Dl::DFSetLabeled{T},
                        maxdepth::Int64, loginterval::Int64)

  logs = default_logs()
  add_folder!(logs, "members", [ASCIIString, ASCIIString, Int64], ["members_true", "members_false", "decision_id"])

  num_data = length(Dl)
  T1 = Bool #predict_type
  T2 = Int64 #label_type

  p = DTParams(num_data, maxdepth, T1, T2)

  dtree = build_tree(p, Dl, problem, psa_params, logs, loginterval) #userargs...

  return dtree, logs
end

function acasx_sa_tree(;outdir::AbstractString="./ACASX_SA_Tree",
                       seed=1,
                       logfileroot::AbstractString="acasx_sa_tree_log",

                       runtype::Symbol=:nmacs_vs_nonnmacs,
                       data::AbstractString="dasc",
                       manuals::AbstractString="dasc_manual",
                       clusterdataname::AbstractString="josh1",

                       maxsteps::Int64=20,
                       T1::Float64=12.184,
                       alpha::Float64=0.99976,
                       n_epochs::Int64=50,
                       n_starts::Int64=1,
                       n_threads::Int64=1,
                       maxdepth::Int64=1,

                       loginterval::Int64=100,
                       vis::Bool=true,
                       plotpdf::Bool=true,
                       limit_members::Int64=10)
  mkpath(outdir)

  problem = ACASXClustering(runtype, data, manuals, clusterdataname)

  observer = Observer()
  par_observer = Observer()

  sa_params = SAESParams(maxsteps, T1, alpha, n_epochs, n_starts, observer)
  psa_params = PSAESParams(n_threads, sa_params, par_observer)

  Dl = problem.Dl
  dtree, logs = train_dtree(psa_params, problem, Dl, maxdepth, loginterval)

  #add to log
  push!(logs, "parameters", ["seed", seed, 0])
  push!(logs, "parameters", ["runtype", runtype, 0])
  push!(logs, "parameters", ["clusterdataname", clusterdataname, 0])

  outfile = joinpath(outdir, "$(logfileroot).txt")
  save_log(outfile, logs)

  #visualize
  if vis
    decisiontreevis(dtree, Dl, joinpath(outdir, "$(logfileroot)_vis"), limit_members,
                    FMT_PRETTY, FMT_NATURAL; plotpdf=plotpdf)
  end

  return dtree, logs
end

end #module


text:
function NiellMappingFunctionModel(arg0::jdouble)
    return NiellMappingFunctionModel((jdouble,), arg0)
end

function NiellMappingFunctionModel(arg0::jdouble, arg1::TimeScale)
    return NiellMappingFunctionModel((jdouble, TimeScale), arg0, arg1)
end

function get_parameters_drivers(obj::NiellMappingFunctionModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function mapping_factors(obj::NiellMappingFunctionModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::NiellMappingFunctionModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(MPI.Get_library_version)})
    precompile(Tuple{typeof(MPI.__init__)})
end


text:
using Test, ARFIMA, Statistics, Random, StatsBase

N = 1000
σ = 2.3

@testset "Size Validation" begin
    x = arfima(N, σ, nothing) # white noise
    @test abs(std(x) - σ) < 1
    @test length(x) == N

    x = arfima(N, σ, nothing, SVector(0.2, 0.2))
    @test length(x) == N
    @test all(!isinf, x)
    x = arfima(N, σ, nothing, SVector(0.2,0.2), SVector(0.1))
    @test length(x) == N
    @test all(!isinf, x)
    x = arfima(N, σ, 1, SVector(0.02,0.02), SVector(0.01))
    @test length(x) == N
    @test all(!isinf, x)
    x = arfima(N, σ, 0.25, SVector(0.02,0.02), SVector(0.01))
    @test length(x) == N
    @test all(!isinf, x)
end

@testset "ARMA" begin
    # white noise
    rng = Random.MersenneTwister(1)
    noise = arma(rng, N, σ, nothing, nothing);
    @test std(noise) ≈ σ atol=0.01

    # AR
    @testset "AR Noise" for φ in [SVector(zeros(1)...), SVector(zeros(3)...), SVector(zeros(5)...)]
        rng = Random.MersenneTwister(1)
        x = arma(rng, N, σ, φ, nothing); # produce white noise
        @test autocor(noise, 1:2) ≈ autocor(x, 1:2) atol=0.05
        @test mean(noise - x) ≈ 0.0 atol=0.05
    end

    rng = Random.MersenneTwister(1)
    φ = 0.8
    ar = arma(rng, 100*N, σ, SVector(φ), nothing);
    @test [φ, φ^2] ≈ autocor(ar, 1:2) atol=0.01

    # MA
    @testset "MA Noise" for θ in [SVector(zeros(1)...), SVector(zeros(3)...), SVector(zeros(5)...)]
        rng = Random.MersenneTwister(1)
        x = arma(rng, N, σ, nothing, θ); # produce white noise
        @test autocor(noise, 1:2) ≈ autocor(x, 1:2) atol=0.05
        @test mean(noise - x) ≈ 0.0 atol=0.05
    end

    rng = Random.MersenneTwister(1)
    θ = 0.96
    ma = arma(rng, 100*N, σ, nothing, SVector(θ));
    @test [-θ/(1+θ^2),0] ≈ autocor(ma, 1:2) atol=0.01

    # ARMA
    @testset "ARMA Noise" for θ in [SVector(zeros(1)...), SVector(zeros(2)...)],
                              φ in [SVector(zeros(1)...), SVector(zeros(3)...), SVector(zeros(5)...)]
        rng = Random.MersenneTwister(1)
        x = arma(rng, N, σ, φ, θ); # produce white noise
        @test autocor(noise, 1:2) ≈ autocor(x, 1:2) atol=0.05
        @test mean(noise - x) ≈ 0.0 atol=0.05
    end

    rng = Random.MersenneTwister(1)
    φ = 0.7
    θ = 0.2
    x = arma(rng, 100*N, σ, SVector(φ), SVector(θ));
    ρ(k) = k == 0 ? (1+φ*θ)*(φ+θ)/(1+2φ*θ+θ^2) : φ*ρ(k-1)
    @test ρ.(1:3) ≈ autocor(x, 1:3) atol=0.05
end


text:
# register.jl
# Register NIfTI volumes based on header information

using NIfTI, Grid

# Affine matrices that add or subtract one
# Needed because Julia and Grid.jl assume zero-based indexing
const ADD1 = [1 0 0 1
              0 1 0 1
              0 0 1 1
              0 0 0 1]
const SUB1 = [1 0 0 -1
              0 1 0 -1
              0 0 1 -1
              0 0 0 1]

# Round to an integer if the second argument is an integer type
mayberound{T<:Integer}(x, ::Type{T}) = round(T, x)
mayberound(x, ::DataType) = x

function apply_rotation!(out, A, data)
    A = ADD1*A*SUB1

    # Get pointers to each pixel dim of data
    s = stride(data, 4)
    ptrs = [pointer_to_array(pointer(data)+(i-1)*s, (size(data, 1), size(data, 2), size(data, 3)))
            for i = 1:size(data, 4)]

    # Allocate B, which represents voxels in targ, and C, which
    # represents voxels in mov
    B = Array(Float32, 4)
    B[4] = 1
    C = Array(Float32, 4)
    Csub = pointer_to_array(pointer(C), 3)

    ic = InterpGridCoefs(ptrs[1], InterpLinear)
    for i = 1:size(out, 3)
        B[3] = i
        for j = 1:size(out, 2)
            B[2] = j
            for k = 1:size(out, 1)
                B[1] = k
                # Convert voxels in targ coords to voxels in mov coords
                A_mul_B!(C, A, B)

                # Interpolate voxels in targ coords from voxels in mov
                # coords
                set_position(ic, BCnan, false, Csub)
                v = interp(ic, ptrs[1])
                if !isnan(v)
                    out[k, j, i, 1] = mayberound(v, eltype(out))
                    for l = 2:size(data, 4)
                        out[k, j, i, l] = mayberound(interp(ic, ptrs[l]), eltype(out))
                    end
                end
            end
        end
    end
end

function register{T<:FloatingPoint}(targ::NIVolume, mov::NIVolume{T}, outtype::DataType=T)
    # Construct the affine matrix that maps voxels of targ to voxels of
    # mov
    targaffine = getaffine(targ.header)
    movaffine = getaffine(mov.header)
    A = inv(movaffine)*targaffine

    # Allocate output
    out = zeros(outtype, size(targ, 1), size(targ, 2), size(targ, 3), size(mov, 4))

    # Perform rotation
    apply_rotation!(out, A, mov.raw)

    # Set header of mov
    newheader = deepcopy(mov.header)
    setaffine(newheader, targaffine)
    newheader.pixdim = (newheader.pixdim[1], targ.header.pixdim[2:4]..., newheader.pixdim[5:end]...)
    NIVolume(newheader, mov.extensions, out)
end
register(targ::NIVolume, mov::NIVolume) =
    register(targ, NIVolume(mov.header, mov.extensions, map(Float32, mov.raw)), eltype(mov))
register(targ::AbstractString, mov::AbstractString, out::AbstractString) =
    niwrite(out, register(niread(targ, mmap=true), niread(mov, mmap=true)))


text:
function assert_option(x)
    return is_option(x) || error("$(typeof(x)) is not an option type")
end

@option struct ToDictOption
    include_defaults::Bool = true
    exclude_nothing::Bool = false
end

"""
    TOMLStyle::ToDictOption

Predefined option for TOML compatible [`to_dict`](@ref) option.
"""
const TOMLStyle = ToDictOption(; include_defaults=true, exclude_nothing=true)

"""
    YAMLStyle::ToDictOption

Predefined option for YAML compatible [`to_dict`](@ref) option.
"""
const YAMLStyle = ToDictOption(; include_defaults=true, exclude_nothing=false)

"""
    JSONStyle::ToDictOption

Predefined option for JSON compatible [`to_dict`](@ref) option.
"""
const JSONStyle = ToDictOption(; include_defaults=true, exclude_nothing=false)

"""
    to_dict(x; include_defaults=true, exclude_nothing=false) -> OrderedDict

Convert an object `x` to an `OrderedDict`.

# Kwargs

- `include_defaults`: include the default value, default is `true`.
- `exclude_nothing`: exclude fields that have value `nothing`,
    this supersedes `include_defaults` when they are both `true`.

# Format Compatibilty

When mapping an option struct from Julia to TOML/YAML/JSON/etc. format,
there are some subtle semantic compatibilty one need to deal with, we provide
some convenient predefined conversion option constants as [`TOMLStyle`](@ref),
[`YAMLStyle`](@ref), [`JSONStyle`](@ref).

!!! tips
    `to_dict` does not export fields that are of the same values as the defaults.
    In most cases, this should be the default behaviour, and users should not use
    `include_defaults`, however,  this can be overridden by changing `include_defaults`
    to `true`.
"""
function to_dict(x; kw...)
    return to_dict(x, ToDictOption(; kw...))
end

"""
    to_dict(x, option::ToDictOption) -> OrderedDict

Convert an object `x` to an `OrderedDict` with `ToDictOption` specified.

# Example

```julia
to_dict(x, TOMLStyle) # TOML compatible
to_dict(x, YAMLStyle) # YAML compatible
to_dict(x, JSONStyle) # JSON compatible
```
"""
function to_dict(x, option::ToDictOption)
    assert_option(x)
    return to_dict(typeof(x), x, option)
end

# disambiguity
to_dict(x::Type, ::ToDictOption) = error("$x is not an option type")

"""
    to_dict(::Type{T}, x, option::ToDictOption) where T

Convert `x` when `x` is inside an option type `T`. `option`
is a set of options to determine the conversion behaviour. this can
be overloaded to change the behaviour of `to_dict(x; kw...)`.

    to_dict(::Type{T}, x) where T

One can also use the 2-arg version when `x` is not or does not
contain an option type for convenience.

# Example

The following is a builtin overload to handle list of options.

```julia
function Configurations.to_dict(::Type{T}, x::Vector, option::ToDictOption) where T
    if is_option(eltype(x))
        return map(p->to_dict(T, p, include_defaults), x)
    else
        return x
    end
end
```

The following overloads the 2-arg `to_dict` to convert all `VersionNumber` to
a `String` for all kinds of option types.

```julia
Configurations.to_dict(::Type, x::VersionNumber) = string(x)
```
"""
function to_dict(::Type{T}, x, option::ToDictOption) where {T}
    if is_option(x)
        return _option_to_dict(x, option)
    else
        return to_dict(T, x) # fall through 2-arg version
    end
end

to_dict(::Type, x) = x

# handle list of options as builtin
function to_dict(::Type{T}, x::Vector, option::ToDictOption) where {T}
    return map(x) do each
        d = to_dict(T, each, option)
        if eltype(x) isa Union && is_option(each)
            FieldType = typeof(each)
            alias = type_alias(FieldType)
            idx = find_reflect_field(FieldType)
            if alias !== nothing && idx === nothing
                return OrderedDict{String,Any}(alias => d)
            end
        end
        return d
    end
end

function _option_to_dict(x, option::ToDictOption)
    assert_option(x)

    d = OrderedDict{String,Any}()
    T = typeof(x)
    for name in fieldnames(T)
        type = fieldtype(T, name)
        value = getfield(x, name)
        name_str = string(name)
        if option.exclude_nothing && value === nothing
            continue
        end

        if option.include_defaults || value != field_default(T, name)
            field_dict = to_dict(T, value, option)

            # 1. option type contains field of type Reflect
            # 2. option type contains field of Union{options...}
            # 3. other types
            #
            # NOTE:
            # we always add an alias if it's a Union
            # of multiple option types
            if type === Reflect
                # we don't implement this via `to_dict`
                # because we want it error when `Reflect`
                # is used as a normal type
                # `Reflect` should be only used to denote
                # a field contains the type info as a `String`.
                if type_alias(T) === nothing
                    d[name_str] = full_typename(T)
                else
                    d[name_str] = type_alias(T)
                end
            elseif is_option(value) && is_union_of_multiple_options(type)
                if contains_reflect_type(typeof(value))
                    d[name_str] = field_dict
                    continue
                end

                alias = type_alias(typeof(value))
                if alias === nothing
                    error("please define an alias for option type $(typeof(value))")
                end

                d[name_str] = OrderedDict{String,Any}(alias => field_dict)
            else
                d[name_str] = field_dict
            end
        end
    end
    return d
end

"""
    is_option_maybe(::Type{T}) where T

`T` is an option struct or if `T` is an union, one of the types
is an option struct.
"""
function is_option_maybe(::Type{T}) where {T}
    is_option(T) && return true
    T isa Union || return false
    return is_option_maybe(T.a) || is_option_maybe(T.b)
end

"""
    to_toml([f::Function], io::IO, option; sorted=false, by=identity, kw...)

Convert an instance `option` of option type to TOML and write it to `IO`. See [`to_dict`](@ref)
for other valid keyword options. See also `TOML.print` in the stdlib for the explaination of
`sorted`, `by` and `f`.

# Exclude `nothing`

In TOML specification, there is [no null type](https://github.com/toml-lang/toml/issues/802). One
should exclude the field if it is not specified (of value `nothing` in Julia). In `to_toml` the option
`exclude_nothing` is always `true`.

In most cases, `nothing` is used with another type to denote optional or not specified field,
thus one should always put a default
value `nothing` to the option struct, e.g

One should define

```julia
@option struct OptionX
    a::Union{Nothing, Int} = nothing
    b::Maybe{Int} = nothing
end
```

Here `Maybe{T}` is a convenient alias of `Union{Nothing, T}`.
"""
function to_toml(f, io::IO, x; sorted::Bool=false, by=identity, kw...)
    is_option(x) || error("argument is not an option type")
    d = to_dict(x; exclude_nothing=true, kw...)
    return TOML.print(f, io, d; sorted=sorted, by=by)
end

"""
    to_toml([f::Function], filename::String, option; sorted=false, by=identity, kw...)

Convert an instance `option` of option type to TOML and write it to `filename`. See also `TOML.print`.
"""
function to_toml(f, filename::String, x; sorted::Bool=false, by=identity, kw...)
    open(filename, "w+") do io
        to_toml(f, io, x; sorted=sorted, by=by, kw...)
    end
end

function to_toml(io::IO, x; sorted::Bool=false, by=identity, kw...)
    return to_toml(identity, io, x; sorted=sorted, by=by, kw...)
end

function to_toml(filename::String, x; sorted::Bool=false, by=identity, kw...)
    return to_toml(identity, filename, x; sorted=sorted, by=by, kw...)
end

"""
    to_toml(x; sorted=false, by=identity, kw...)

Convert an instance `x` of option type to TOML and write it to `String`. See also `TOML.print`.

`to_toml` does not export fields that are of the same values as the defaults. This can be
overridden by changing `include_defaults` to `true`.
"""
function to_toml(x; sorted::Bool=false, by=identity, kw...)
    return sprint(x) do io, x
        to_toml(io, x; sorted=sorted, by=by, kw...)
    end
end


text:
#!/usr/bin/env julia

# Usage example:
#      gen_namelist.jl global args  > fort.4
# will merge global.toml and args.toml and convert to namelist format

using Harmonie_namelists

io = stdout   # io = open("fort.4")

dicts = read_namelist.(ARGS)
totdict = merge_namelists(dicts)
replace_env!(totdict)
dict2namelist(io,totdict)


text:
push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, SimpleSDMLayers
using GBIF
using Statistics
using Literate

# Literate files
for ENDING in ["examples", "sdm"]
    expl = joinpath("docs", "src", ENDING)
    corefiles = [joinpath(expl, f) for f in readdir(expl)]
    filter!(f -> endswith(f, "jl"), corefiles)
    vignetteconfig = Dict(
        "repo_root_url" => "https://github.com/EcoJulia/SimpleSDMLayers.jl",
        "flavor" => Literate.DocumenterFlavor(),
        "credit" => false
    )
    for corefile in corefiles
        Literate.markdown(corefile, expl; config=vignetteconfig)
    end
end

makedocs(
    sitename = "Simple SDM Layers",
    modules = [SimpleSDMLayers],
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Types" => "man/types.md",
            "Indexing" => "man/indexing.md",
            "Clipping" => "man/clip.md",
            "Operations on values" => "man/changevalues.md",
            "Overloads" => "man/overloads.md",
            "Other operations" => "man/operations.md",
            "Data access" => "man/data.md",
            "IO" => "man/io.md"
        ],
        "General examples" => [
            "Introduction: elevation data" => "examples/elevation.md",
            "Geometry for clipping" => "examples/geometry.md",
            "Sliding window analysis" => "examples/slidingwindow.md",
            "Landcover data" => "examples/landcover.md",
            "Multivariate mapping" => "examples/multivariate.md"
        ],
        "SDM case studies" => [
            "GBIF integration" => "sdm/gbif.md",
            "Variable selection (VIF)" => "sdm/vif.md",
            "Building the BIOCLIM model" => "sdm/bioclim.md",
            "Pseudo-absences" => "sdm/pseudoabsences.md",
            "Dealing with future data" => "sdm/future.md",
            "BRTs and climate change" => "sdm/brt.md"
        ]
    ]
)

run(`find . -type f -size +40M -delete`)

deploydocs(
    repo = "github.com/EcoJulia/SimpleSDMLayers.jl.git",
    push_preview = true,
    devbranch = "main"
)


text:
module LIBRA

using UnitSystems, Parameters, LinearAlgebra, FFTW
using Dierckx, Arpack, Infiltrator, Statistics
export C_e, Negative, Constants, Positive, Seperator, Flux, C_se, Phi_s, Phi_e, Phi_se, DRA
export RealisationAlgorthim, TransferFun, flatten, R, F, CellDef, Sim_Model, D_Linear, _bisection, cell

include("Functions/Transfer/C_e.jl")
include("Functions/Transfer/C_se.jl")
include("Functions/Transfer/Flux.jl")
include("Functions/Transfer/Phi_s.jl")
include("Functions/Transfer/Phi_e.jl")
include("Functions/Transfer/Phi_se.jl")
include("Methods/DRA.jl")
include("Functions/Sim_Model.jl")

const F,R = faraday(Metric), universal(SI2019) #Faraday Constant / Universal Gas Constant
const Debug = 0 #Print Variables for Debugging

function D_Linear(CellData, ν_neg, ν_pos, σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, κ_eff_Sep)
    D = Array{Float64}(undef,0,1)
    Dt = Array{Float64}(undef,0,1)
    D_ = Array{Float64}(undef,5,1)
    q = Int64(1)
    tfs = CellData.Transfer.tfs
    for i in 1:size(tfs,1)
        if tfs[i,1] == C_e
            Dt =  zeros(length(tfs[i,3]))
        elseif tfs[i,1] == Phi_e
            Dt = Array{Float64}(undef,0,1)
            for pt in tfs[i,3]
                if pt <= CellData.Neg.L+eps()
                D_[q]  =  @. (CellData.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg*pt/CellData.Neg.L)) - pt*ν_neg*sinh(ν_neg) + CellData.Neg.L*(cosh(ν_neg)-cosh(ν_neg*(CellData.Neg.L-pt)/CellData.Neg.L)))/(CellData.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg)*ν_neg) #Lee. Eqn. 4.22 @ ∞
                elseif pt <= CellData.Neg.L + CellData.Sep.L + eps()
                D_[q] = @. (CellData.Neg.L - pt)/(CellData.Const.CC_A*κ_eff_Sep) + (CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/(CellData.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) #Lee. Eqn. 4.23 @ ∞
                else
                D_[q] = @. -CellData.Sep.L/(CellData.Const.CC_A*κ_eff_Sep) + CellData.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg)/(CellData.Const.CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg) + (CellData.Pos.L*(-σ_eff_Pos*cosh(ν_pos) + σ_eff_Pos*cosh((CellData.Const.Ltot-pt)*ν_pos/CellData.Pos.L) +  κ_eff_Pos*(cosh((pt-CellData.Neg.L-CellData.Sep.L)*ν_pos/CellData.Pos.L)-1)) - (pt-CellData.Neg.L-CellData.Sep.L)*κ_eff_Pos*sinh(ν_pos)*ν_pos)/(CellData.Const.CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos*sinh(ν_pos))
                end
                Dt = [Dt; D_[q]]
                q = q+1
            end

        elseif tfs[i,1] == C_se
            Dt = zeros(length(tfs[i,3]))
        elseif tfs[i,2] == "Pos"
            σ_eff = σ_eff_Pos
            κ_eff = κ_eff_Pos
            if tfs[i,1] == Phi_se
                Dt = @. -1*CellData.Pos.L/(CellData.Const.CC_A*ν_pos*sinh(ν_pos))*((1/κ_eff)*cosh(ν_pos*tfs[i,3])+(1/σ_eff)*cosh(ν_pos*(tfs[i,3]-1))) # Contribution to D as G->∞
            elseif tfs[i,1] == Flux
                Dt = @. -1*ν_pos*(σ_eff*cosh(ν_pos*tfs[i,3])+κ_eff*cosh(ν_pos*(tfs[i,3]-1)))/(CellData.Pos.as*F*CellData.Pos.L*CellData.Const.CC_A*(κ_eff+σ_eff)*sinh(ν_pos))
            elseif tfs[i,1] == Phi_s
                Dt = @. -1*(-CellData.Pos.L*(κ_eff*(cosh(ν_pos)-cosh(tfs[i,3]-1)*ν_pos))-CellData.Pos.L*(σ_eff*(1-cosh(tfs[i,3]*ν_pos)+tfs[i,3]*ν_pos*sinh(ν_pos))))/(CellData.Const.CC_A*σ_eff*(κ_eff+σ_eff)*ν_pos*sinh(ν_pos)) # Contribution to D as G->∞
            end
        elseif tfs[i,2] == "Neg"
            σ_eff = σ_eff_Neg
            κ_eff = κ_eff_Neg
            if tfs[i,1] == Phi_se
                Dt = @. CellData.Neg.L/(CellData.Const.CC_A*ν_neg*sinh(ν_neg))*((1/κ_eff)*cosh(ν_neg*tfs[i,3])+(1/σ_eff)*cosh(ν_neg*(tfs[i,3]-1))) # Contribution to D as G->∞
            elseif tfs[i,1] == Flux
                Dt = @. ν_neg*(σ_eff*cosh(ν_neg*tfs[i,3])+κ_eff*cosh(ν_neg*(tfs[i,3]-1)))/(CellData.Neg.as*F*CellData.Neg.L*CellData.Const.CC_A*(κ_eff+σ_eff)*sinh(ν_neg))
            elseif tfs[i,1] == Phi_s
                Dt = @. (-CellData.Neg.L*(κ_eff*(cosh(ν_neg)-cosh(tfs[i,3]-1)*ν_neg))-CellData.Neg.L*(σ_eff*(1-cosh(tfs[i,3]*ν_neg)+tfs[i,3]*ν_neg*sinh(ν_neg))))/(CellData.Const.CC_A*σ_eff*(κ_eff+σ_eff)*ν_neg*sinh(ν_neg)) # Contribution to D as G->∞
            end
        end
        D = [D; Dt]
    end
    return D
end

function flatten end
flatten() = ()
flatten(a::Tuple) = Tuple(a)
flatten(a) = (a,)
flatten(a::Tuple, b...) = tuple(a..., flatten(b...)...)
flatten(a, b...) = tuple(a, flatten(b...)...)
flatten_tuple(x::Tuple) = flatten(x...)

"""
    bisection(f, a, b; fa = f(a), fb = f(b), ftol, wtol)

Bisection algorithm for finding the root ``f(x) ≈ 0`` within the initial bracket
`[a,b]`.

Returns a named tuple

`(x = x, fx = f(x), isroot = ::Bool, iter = ::Int, ismaxiter = ::Bool)`.

Terminates when either

1. `abs(f(x)) < ftol` (`isroot = true`),
2. the width of the bracket is `≤wtol` (`isroot = false`),
3. `maxiter` number of iterations is reached. (`isroot = false, maxiter = true`).

which are tested for in the above order. Therefore, care should be taken not to make `wtol` too large.

"""
function bisection(f,CellData, a::Real, b::Real; fa::Real = f(a), fb::Real = f(b),
                   ftol = √eps(), wtol = 0, maxiter = 100)
    @assert fa * fb ≤ 0 "initial values don't bracket zero"
    @assert isfinite(a) && isfinite(b)
    _bisection(f, CellData, float.(promote(a, b, fa, fb, ftol, wtol))..., maxiter)
end

function _bisection(f,CellData, a, b, fa, fb, ftol, wtol, maxiter)
    iter = 0
    abs(fa) < ftol && return (x = a, fx = fa, isroot = true, iter = iter, ismaxiter = false)
    abs(fb) < ftol && return (x = b, fx = fb, isroot = true, iter = iter, ismaxiter = false)
    while true
        iter += 1
        m = middle(a, b)
        fm = f(CellData,m)
        abs(fm) < ftol && return (x = m, fx = fm, isroot = true, iter = iter, ismaxiter = false)
        abs(b-a) ≤ wtol && return (x = m, fx = fm, isroot = false, iter = iter, ismaxiter = false)
        if fa * fm > 0
            a, fa = m, fm
        else
            b, fb = m, fm
        end
        iter == maxiter && return (x = m, fx = fm, isroot = false, iter = iter, ismaxiter = true)
    end
end

function cell(CellType)
    if CellType == "Doyle_94"
        CellType = string(CellType,".jl")
        include(joinpath(dirname(pathof(LIBRA)), "Data/Doyle_94", CellType))
    elseif CellType == "LG_M50"
        CellType = string(CellType,".jl")
        include(joinpath(dirname(pathof(LIBRA)), "Data/Chen_2020", CellType))
    end
    return CellData
end

end # module

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
using Gallium
using Base.Test
include("Hooking.jl")
# Breakpointing needs to be run in a separate process since it overrides methods
# for testing purposes
run(`$(Base.julia_cmd()) -f $(joinpath(dirname(@__FILE__),"breakpointing.jl"))`)

# #self# needs to be captured properly (for clousures)
fclosure,gclosure = let y = 2
    (()->(@Base._noinline_meta; y),()->(@Base._noinline_meta; y=3))
end

global hit_breakpoint = false
@conditional breakpoint(fclosure, Tuple{}) (global hit_breakpoint = true; false)

@test fclosure()==2
@test hit_breakpoint

# Breakpointing functions with keyword arguments
hit_breakpoint = false
@noinline fkeywords(x; a = 2) = x + a
@conditional breakpoint(fkeywords, Tuple{Int64}) (global hit_breakpoint = true; false)

@test fkeywords(1) == 3
@test hit_breakpoint

# Floating point value recovery
hit_breakpoint = false
function f131(x)
    y=x+1
    z=2y
    return z
end
@conditional breakpoint(f131) (global hit_breakpoint = true; (@test x==3.0); false)
f131(3.0)
@test hit_breakpoint

# Floating point value recovery in breakpoint()
@noinline function do_xmmbptest()
    # Manually do what's done in breakpoint() to look at the stack
    RC = Hooking.getcontext()
    # -1 to skip breakpoint (getcontext is inlined)
    stack, RCs = Gallium.stackwalk(RC; fromhook = false)
    frame = collect(filter(x->isa(x,Gallium.JuliaStackFrame),stack))[end-1]
    id = findfirst(sym->sym==:y,frame.linfo.def.lambda_template.slotnames)
    # For now only test this on OS X since it doesn't work elsewhere
    # TODO: Figure out why
    if is_apple()
        @test 2.0 == get(frame.env.locals[id])
    end
end

function xmmbptest(x)
    y = x+1.0
    do_xmmbptest()
    # Make sure to keep this live
    y
end
@test xmmbptest(1.0) == 2.0

# Floating point value recovery in breakpoint_on_error()
using Gallium
function xmmerrtest(x)
    y = x+1.0
    error(y)
    y
end
hit_breakpoint = false
Gallium.breakpoint_on_error()
push!(Gallium.bp_on_error_conditions,:(is_apple() && @test y == 2.0; global hit_breakpoint = true; false))
try; xmmerrtest(1.0); end
@test hit_breakpoint
empty!(Gallium.bp_on_error_conditions)
Gallium.breakpoint_on_error(false)

include("fakememory.jl")


text:
# These tests don't do much except ensure that the package loads,
# and does something sensible if it does.
using Loess
using Base.Test

srand(100)
xs = 10 .* rand(100)
ys = sin(xs) .+ 0.5 * rand(100)

model = loess(xs, ys)

us = collect(minimum(xs):0.1:maximum(xs))
vs = predict(model, us)

@test minimum(vs) >= -1.1
@test maximum(vs) <= +1.1


text:
using PowerSystems

######################################
############ Generators ##############
######################################

######## Machine Data #########

machine_classic() = BaseMachine(
    0.0, #R
    0.2995, #Xd_p
    0.7087,
)  #MVABase

machine_oneDoneQ() = OneDOneQMachine(
    0.0, #R
    1.3125, #Xd
    1.2578, #Xq
    0.1813, #Xd_p
    0.25, #Xq_p
    5.89, #Td0_p
    0.6, #Tq0_p
)

machine_simple_marconato() = SimpleMarconatoMachine(
    0.0,
    1.3125, #Xd
    1.2578, #Xq
    0.1813, #Xd_p
    0.25, #Xq_p
    0.14, #Xd_pp
    0.18, #Xq_pp
    5.89, #Td0_p
    0.6, #Tq0_p
    0.5, #Td0_pp
    0.023, #Tq0_pp
    0.0,
) #MVABase

machine_marconato() = MarconatoMachine(
    0.0,
    1.3125, #Xd
    1.2578, #Xq
    0.1813, #Xd_p
    0.25, #Xq_p
    0.14, #Xd_pp
    0.18, #Xq_pp
    5.89, #Td0_p
    0.6, #Tq0_p
    0.5, #Td0_pp
    0.023, #Tq0_pp
    0.0,
) #MVABase

machine_anderson() = AndersonFouadMachine(
    0.0, #R
    0.8979, #Xd
    0.646, #Xq
    0.2995, #Xd_p
    0.646, #Xq_p
    0.23, #Xd_pp
    0.4, #Xq_pp
    3.0, #Td0_p
    0.1, #Tq0_p
    0.01, #Td0_pp
    0.033, #Tq0_pp
)

machine_simple_anderson() = SimpleAFMachine(
    0.0, #R
    0.8979, #Xd
    0.646, #Xq
    0.2995, #Xd_p
    0.646, #Xq_p
    0.23, #Xd_pp
    0.4, #Xq_pp
    3.0, #Td0_p
    0.1, #Tq0_p
    0.01, #Td0_pp
    0.033, #Tq0_pp
)

machine_genrou() = RoundRotorExponential(
    R = 0.0,
    Td0_p = 8.0,
    Td0_pp = 0.03,
    Tq0_p = 0.4,
    Tq0_pp = 0.05,
    Xd = 1.8,
    Xq = 1.7,
    Xd_p = 0.3,
    Xq_p = 0.55,
    Xd_pp = 0.25,
    Xl = 0.2,
    Se = (0.0, 0.0),
)

#Not available yet
#=
machine_kundur() = SimpleFullMachine(
    0.003, #R on Example 3.1 and 4.1 of Kundur
    0.0006, #R_f
    0.0284, #R_1d or RD in Machowski
    0.0062, #R_1q or RQ on Machowski
    1.81, #L_d
    1.76, #L_q
    1.66, #L_ad or k*M_f or k*M_D in Machowski
    1.61, #L_aq or k*M_Q in Machowski
    1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
    1.825, #L_ff
    0.1713, #L_1d or L_D in Machowski
    0.7525, #L_1q or L_Q in Machowski
    555.0,
) #MVABase

machine_full_kundur() = FullMachine(
    0.003, #R on Example 3.1 and 4.1 of Kundur
    #0.0006, #R_f
    0.003, #R_f
    0.0284, #R_1d or RD in Machowski
    0.0062, #R_1q or RQ on Machowski
    1.81, #L_d
    1.76, #L_q
    1.66, #L_ad or k*M_f or k*M_D in Machowski
    1.61, #L_aq or k*M_Q in Machowski
    1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
    1.825, #L_ff
    0.1713, #L_1d or L_D in Machowski
    0.7525, #L_1q or L_Q in Machowski
)
=#

machine_multi_ref() = BaseMachine(
    0.0, #R
    0.2995, #Xd_p
    1.0901,
)  #MVABase

machine_multi() = BaseMachine(
    0.0, #R
    0.2995, #Xd_p
    0.9516,
)  #MVABase

######## Shaft Data #########

shaft_damping() = SingleMass(
    3.148, #H
    2.0,
) #D

shaft_no_damping() = SingleMass(
    3.01, #H (M = 6.02 -> H = M/2)
    0.0,
) #D

shaft_genrou() = SingleMass(H = 6.175, D = 0.05)

shaft_fivemass() = FiveMassShaft(
    3.01, #5.148, #H
    0.3348, #H_hp
    0.7306, #H_ip
    0.8154, #H_lp
    0.0452, #H_ex,
    0.0, #2.0, #D
    0.5180, #D_hp
    0.2240, #D_ip
    0.2240, #D_lp
    0.1450, #D_ex
    0.0518, #D_12
    0.0224, #D_23
    0.0224, #D_34
    0.0145, #D_45
    33.07, #K_hp
    28.59, #K_ip
    44.68, #K_lp
    21.984,
) #K_ex

######## PSS Data #########

pss_none() = PSSFixed(0.0)

######## TG Data #########

tg_none() = TGFixed(1.0) #eff

tg_type1() = TGTypeI(
    0.02, #R
    0.1, #Ts
    0.45, #Tc
    0.0, #T3
    12.0, #T4
    50.0, #T5
    (min = 0.3, max = 1.2), #P_lims
)

tg_type2() = TGTypeII(
    0.05, #R
    2.0, #T1
    1.0, #T2
    (min = 0.1, max = 1.5), #τ_lims
)

########  AVR Data #########

avr_none() = AVRFixed(0.0)

avr_propr() = AVRSimple(500.0) #Kv

avr_fixed() = AVRFixed(1.05) #Emf

avr_type1() = AVRTypeI(
    20.0, #Ka - Gain
    0.01, #Ke
    0.063, #Kf
    0.2, #Ta
    0.314, #Te
    0.35, #Tf
    0.001, #Tr
    (min = -5.0, max = 5.0),
    0.0039, #Ae - 1st ceiling coefficient
    1.555,
) #Be - 2nd ceiling coefficient

avr_type2() = AVRTypeII(
    200.0, #K0 - Gain
    4.0, #T1 - 1st pole
    1.0, #T2 - 1st zero
    0.006, #T3 - 2nd pole
    0.06, #T4 - 2nd zero
    0.0001, #Te - Field current time constant
    0.0001, #Tr - Measurement time constant
    (min = -50.0, max = 50.0),
    0.0, #Ae - 1st ceiling coefficient
    0.0,
) #Be - 2nd ceiling coefficient

######################################
############# Inverters ##############
######################################

###### Converter Data ######
converter_low_power() = AverageConverter(rated_voltage = 690.0, rated_current = 2.75)

converter_high_power() = AverageConverter(rated_voltage = 138.0, rated_current = 100.0)

###### DC Source Data ######
dc_source_lv() = FixedDCSource(voltage = 600.0) #Not in the original data, guessed.
dc_source_hv() = FixedDCSource(voltage = 1500.0) #Not in the original data, guessed.

###### Filter Data ######
filt() = LCLFilter(lf = 0.08, rf = 0.003, cf = 0.074, lg = 0.2, rg = 0.01)
filt_gfoll() = LCLFilter(lf = 0.009, rf = 0.016, cf = 2.5, lg = 0.002, rg = 0.003)

###### PLL Data ######
pll() = KauraPLL(
    ω_lp = 500.0, #Cut-off frequency for LowPass filter of PLL filter.
    kp_pll = 0.084,  #PLL proportional gain
    ki_pll = 4.69,   #PLL integral gain
)

reduced_pll() = ReducedOrderPLL(
    ω_lp = 1.32 * 2 * pi * 50, #Cut-off frequency for LowPass filter of PLL filter.
    kp_pll = 2.0,  #PLL proportional gain
    ki_pll = 20.0,   #PLL integral gain
)

no_pll() = PSY.FixedFrequency()

###### Outer Control ######
function outer_control()
    function virtual_inertia()
        return VirtualInertia(Ta = 2.0, kd = 400.0, kω = 20.0)
    end
    function reactive_droop()
        return ReactivePowerDroop(kq = 0.2, ωf = 1000.0)
    end
    return OuterControl(virtual_inertia(), reactive_droop())
end

function outer_control_nopll()
    function virtual_inertia()
        return VirtualInertia(Ta = 2.0, kd = 0.0, kω = 20.0)
    end
    function reactive_droop()
        return ReactivePowerDroop(kq = 0.2, ωf = 1000.0)
    end
    return OuterControl(virtual_inertia(), reactive_droop())
end

function outer_control_droop()
    function active_droop()
        return PSY.ActivePowerDroop(Rp = 0.05, ωz = 2 * pi * 5)
    end
    function reactive_droop()
        return ReactivePowerDroop(kq = 0.2, ωf = 1000.0)
    end
    return OuterControl(active_droop(), reactive_droop())
end

function outer_control_gfoll()
    function active_pi()
        return ActivePowerPI(Kp_p = 2.0, Ki_p = 30.0, ωz = 0.132 * 2 * pi * 50)
    end
    function reactive_pi()
        return ReactivePowerPI(Kp_q = 2.0, Ki_q = 30.0, ωf = 0.132 * 2 * pi * 50.0)
    end
    return OuterControl(active_pi(), reactive_pi())
end

######## Inner Control ######
inner_control() = VoltageModeControl(
    kpv = 0.59,     #Voltage controller proportional gain
    kiv = 736.0,    #Voltage controller integral gain
    kffv = 0.0,     #Binary variable enabling the voltage feed-forward in output of current controllers
    rv = 0.0,       #Virtual resistance in pu
    lv = 0.2,       #Virtual inductance in pu
    kpc = 1.27,     #Current controller proportional gain
    kic = 14.3,     #Current controller integral gain
    kffi = 0.0,     #Binary variable enabling the current feed-forward in output of current controllers
    ωad = 50.0,     #Active damping low pass filter cut-off frequency
    kad = 0.2,
)

current_mode_inner() = CurrentModeControl(
    kpc = 0.37,     #Current controller proportional gain
    kic = 0.7,     #Current controller integral gain
    kffv = 1.0,     #Binary variable enabling the voltage feed-forward in output of current controllers
)

####### Devices #######

function dyn_gen_second_order(generator)
    return PSY.DynamicGenerator(
        name = get_name(generator),
        ω_ref = 1.0, # ω_ref,
        machine = machine_oneDoneQ(), #machine
        shaft = shaft_no_damping(), #shaft
        avr = avr_type1(), #avr
        prime_mover = tg_none(), #tg
        pss = pss_none(), #pss
    )
end

function inv_case78(static_device)
    return DynamicInverter(
        name = get_name(static_device),
        ω_ref = 1.0, # ω_ref,
        converter = converter_high_power(), #converter
        outer_control = outer_control(), #outer control
        inner_control = inner_control(), #inner control voltage source
        dc_source = dc_source_lv(), #dc source
        freq_estimator = pll(), #pll
        filter = filt(), #filter
    )
end


text:
function grad{T <: Number}(
    f::Function,
    x::T;
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> T
    if method == :finite
        return finitegrad(f, x, direction)
    elseif method == :ad
        return autograd(f, x)
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
end

function grad!{T <: Number}(
    gr::Vector{T},
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Vector{T}
    if method == :finite
        finitegrad!(gr, f, x, direction)
    elseif method == :ad
        # TODO: Raise a warning here that memory is being allocated?
        autograd!(gr, Array(Dual{T}, length(x)), f, x)
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
    return gr
end

function grad!{T <: Number}(
    gr::Vector{T},
    x_d::Vector{Dual{T}},
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Vector{T}
    if method == :finite
        finitegrad!(gr, f, x, direction)
    elseif method == :ad
        autograd!(gr, x_d, f, x)
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
    return gr
end

function grad{T <: Number}(
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Vector{T}
    n = length(x)
    gr = Array(T, n)
    if method == :ad
        x_d = Array(Dual{T}, n)
        grad!(gr, x_d, f, x, method = method, direction = direction)
    else
        grad!(gr, f, x, method = method, direction = direction)
    end
    return gr
end


text:
function rowcol(s)
    # Find row and col of the seatid

    rowpart = s[1:7]
    rowpart = replace(rowpart, "B" => "1")
    rowpart = replace(rowpart, "F" => "0")
    colpart = s[8:end]
    colpart = replace(colpart, "R" => "1")
    colpart = replace(colpart, "L" => "0")
    rowbin, colbin = parse(UInt8, rowpart, base=2), parse(UInt8, colpart, base=2)
    Int.((rowbin & UInt(127), colbin & UInt(7)))
end

seatid((row, col)) = row * 8 + col
# Every seat also has a unique seat ID: multiply the row by 8, then add the column

# First Half

puzzeinput = readlines("inputs/day05.txt")
maximum(map(seatid ∘ rowcol, puzzeinput))

# Second Half

function missingid(strs)

    rowcols = map(rowcol, strs)
    filter!(rc -> rc[1] > 0 && rc[1] < 127, rowcols)
    # Removes very front and very back rows, because your seat wasn't
    # at the very front or back

    sorted_ids = map(seatid, rowcols) |> sort
    for (i, id) in enumerate(sorted_ids)
        if sorted_ids[i + 1] ≠ id + 1
            # Missing ID found

            return id + 1
        end
    end
end

puzzeinput = readlines("inputs/day05.txt")
missingid(puzzeinput)


text:
# EventHistory.jl
# Event History Analysis in Julia
#
# Copyright (C) 2013   Klaus K. Holst
#
# Authors: Klaus Kähler Holst <kkho@biostat.ku.dk>
# Keywords: event history analysis, survival analysis,
#           cox regression, competing risks models
##################################################

module EventHistory

using Distributions, DataFrames
import Base.show
import DataFrames: ModelFrame, ModelMatrix, CoefTable, coefnames
import StatsBase: coef, coeftable, confint, vcov, predict
import Calculus: deparse
using StatsBase: StatisticalModel, RegressionModel

type EventHistoryModel <: RegressionModel
    model::AbstractString
    formula::Formula
    eventtype::DataType
    coef::Vector{Float64}
    coefmat::CoefTable
    IC::Matrix{Float64}
    invhess::Matrix{Float64}
    vcov::Matrix{Float64}
    opt::Vector
    grad::Vector{Float64}
    X::Matrix{Float64}
    eventtime::Matrix
    chaz::Matrix{Float64}
end

function show(io::IO, obj::EventHistoryModel)
    ## outcome = deparse(obj.call.args[2])
    ## covars = deparse(obj.call.args[3])
    ## print(io,"\nModel: ", obj.model,",", obj.eventtype, "; ",
    ##       string("$outcome ~ $covars"))
    print(io,"\nModel: ", obj.model, "; ", obj.formula,"\n")
    n = size(obj.eventtime,1)
    events::Int = sum(obj.eventtime[:,2])
    print(io,"\nn=",n,", events=",events,"\n\n")
    print(io,obj.coefmat)
end

include("event.jl")
include("phreg.jl") ## Cox Proportional Hazards Model
#include("lifetable.jl") ## Life-table construction
#include("aalen.jl") ## Aalens Additive Model
include("optim.jl") ## Optimization routines

export phreg, Event, predict ## lifetable, aalen

end # module


text:
const zeroid = "0"^40

function readpktline(io)
    len = parse(Int, String(read(io,4)), base=16)
    len == 0 && return nothing
    chomp(String(read(io, len-4)))
end

function fetch_refs(io::IO)
    reflist = Pair{String, SHA1}[]
    capabilities = String[]

    refpkt1 = readpktline(io)
    if refpkt1 == "version=1"
        refpkt1 = readpktline(io)
    end
    if refpkt1 === nothing
        return reflist, capabilities
    end

    refpkt, capstr = split(refpkt1, '\0', limit=2)
    capabilities = split(capstr, ' ')

    if refpkt == "$zeroid capabilities^{}"
        refpkt = readpktline(io)
	@assert refpkt === nothing
        return reflist, capabilities
    end

    while refpkt !== nothing
	hashstr, ref = split(refpkt,' ',limit=2)
	push!(reflist, ref => SHA1(hashstr))
	refpkt = readpktline(io)
    end
    return reflist, capabilities
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
"""
    parse_difflist(data, offset, bytes_per_sample_id, has_genotypes)

Parses a single `DiffList` from `data`, starting with `offset + 1`-th byte.
"""
function parse_difflist(data::AbstractVector{UInt8},
    offset::Integer,
    bytes_per_sample_id::Integer,
    has_genotype::Bool)
    # length
    len, offset = decode_single(pointer(data); offset=offset)
    # length-zero list
    sample_id_dtype = bytes_to_UInt[bytes_per_sample_id]
    if len == 0
        return DiffList{sample_id_dtype, Nothing}(
            0, C_NULL, C_NULL, has_genotype, nothing, C_NULL), offset
    end
    # sample id bases
    n_groups = ceil_int(len, 0x000040) # 64 in decimal
    sample_id_bases = reinterpret(sample_id_dtype, view(data,
        Int(offset + 1) : Int(offset + n_groups * bytes_per_sample_id)))

    offset += n_groups * bytes_per_sample_id

    # sizes of the final components
    final_component_sizes = view(data, Int(offset + 1): Int(offset + n_groups - 1))
    offset += n_groups - 1

    # genotypes
    if has_genotype
        genotype_bytes = ceil_int(len, 4)
        # genotype_data = view(data, offset + 1 : offset + genotype_bytes)
        genotypes = BitsVector(pointer(data, offset + 1), 2, len)
        # BitsVector{typeof(genotype_data)}(Ref(genotype_data), 2, len)
        offset += genotype_bytes
    else
        genotypes = nothing
    end

    # final component: differences of indices
    final_component_size = sum(final_component_sizes) + 63 * length(final_component_sizes)
    last_group_size = len % 64
    last_group_size = last_group_size == 0 ? 64 : last_group_size
    # increment has one less value
    last_group_size -= 1
    last_incr_size = size_n(data, last_group_size, offset + final_component_size)
    final_component_size += last_incr_size
    final_component = pointer(data, offset + 1)#view(data, Int(offset + 1) : Int(offset + final_component_size))

    # TODO: defer computation of offset of the last block to parse_difflist_sampleids!
    offset += final_component_size

    DiffList{typeof(sample_id_bases),typeof(genotypes)}(
        len, sample_id_bases, pointer(final_component_sizes),
        has_genotype,
        genotypes, final_component), offset
end

"""
    parse_difflist_sampleids!(idx, idx_incr, dl, gid, sid_incr_offset=nothing)

Parses a `gid`-th group of 64 sample ids in a `DiffList` `dl` into `idx`.
"""
function parse_difflist_sampleids!(idx::AbstractArray, idx_incr::AbstractArray,
        dl::DiffList, gid::Integer, sid_incr_offset::Union{Nothing, UInt} = nothing)
    # TODO: compute offset of the last block in this function.
    n_groups = ceil_int(dl.len, 64)
    @assert gid <= n_groups
    if sid_incr_offset === nothing
        sid_incr_offset = zero(UInt)
        @inbounds for i in 1:(gid - 1)
            sid_incr_offset += unsafe_load(dl.last_component_sizes, i) + 0x3f
            #dl.last_component_sizes[][i] + 0x3f # offset by 63 (0x3f)
        end
    end
    idx_incr[1] = 0
    # offset by one to make it 1-based.
    baseidx = dl.sample_id_bases[gid] + 1
    if gid != n_groups
        decode_multiple!(pointer(idx_incr, 2), dl.sample_id_increments;
            count = 63, offset = UInt(sid_incr_offset))
        cumsum!(idx_incr, idx_incr)
        idx .= baseidx .+ idx_incr
    else
        count = dl.len % 64
        count = count == 0 ? 0x00000040 : count
        decode_multiple!(pointer(idx_incr, 2), dl.sample_id_increments;
            count = count - 1, offset = UInt(sid_incr_offset))
        idx[1] = baseidx
        @inbounds for i in 2:count
            idx[i] = idx[i-1] + idx_incr[i]
        end
        @inbounds for i in (count + 1):length(idx)
            idx[i] = 0
        end
        # cumsum!(view(idx_incr, 1:count), view(idx_incr, 1:count))
        # idx[1:count] .= baseidx .+ view(idx_incr, 1:count)
        # fill!(@view(idx[(count + 1):end]), 0)
    end
    idx
end


text:
export quadform

function quadform(x::Value, A::AbstractExpr)
    return x' * A * x
end

function quadform(x::AbstractExpr, A::Value)
    if length(size(A)) != 2 || size(A, 1) != size(A, 2)
        error("Quadratic form only takes square matrices")
    end
    if !issymmetric(A)
        error("Quadratic form only defined for symmetric matrices")
    end
    V = eigvals(Symmetric(Matrix(A)))

    if all(V .>= 0)
        factor = 1
    elseif all(V .<= 0)
        factor = -1
    else
        error("Quadratic forms supported only for semidefinite matrices")
    end

    P = real(sqrt(Matrix(factor * A)))
    return factor * square(norm2(P * x))
end


text:
module FlexibleLicenses

using Dates
import PkgTemplates: FilePlugin, Template, customizable, default_file, defaultkw, destination, prompt, source, view

"""
    default_format_year(2021)

Omitting year is gaining popularity (e.g. Facebook),
but I believe it breaks some license identification tools.
"""
default_format_year(y) = "$y-present"

"""
    OwnerFlexibleLicense(; owner="Yoyodynene, Inc.", name="MIT", path=nothing, destination="LICENSE")
Creates a license file. Provides the ability to specify an owner other than the package authors.
## Keyword Arguments
- `name::AbstractString`: Name of a license supported by PkgTemplates.
  Available licenses can be seen
  [here](https://github.com/invenia/PkgTemplates.jl/tree/master/templates/licenses).
- `path::Union{AbstractString, Nothing}`: Path to a custom license file.
  This keyword takes priority over `name`.
- `destination::AbstractString`: File destination, relative to the repository root.
  For example, `"LICENSE.md"` might be desired.
"""
struct FlexibleLicense <: FilePlugin
    path::String
    destination::String
    owner::Union{Nothing,String}
    format_year::Function
end

function FlexibleLicense(;
    name::AbstractString="MIT",
    path::Union{AbstractString, Nothing}=nothing,
    destination::AbstractString="LICENSE",
    owner::Union{Nothing,String}=nothing,
    format_year::Function=default_format_year,
)
    if path === nothing
        path = if name == "UNLICENSED"
            joinpath(@__DIR__, "..", "templates", "UNLICENSED")
        else
            default_file("licenses", name)
        end
        isfile(path) || throw(ArgumentError("FlexibleLicense '$(basename(path))' is not available"))
    end
    return FlexibleLicense(path, destination, owner, format_year)
end

defaultkw(::Type{FlexibleLicense}, ::Val{:path}) = nothing
defaultkw(::Type{FlexibleLicense}, ::Val{:name}) = "MIT"
defaultkw(::Type{FlexibleLicense}, ::Val{:destination}) = "LICENSE"
defaultkw(::Type{FlexibleLicense}, ::Val{:owner}) = nothing
defaultkw(::Type{FlexibleLicense}, ::Val{:format_year}) = default_format_year

source(p::FlexibleLicense) = p.path
destination(p::FlexibleLicense) = p.destination
view(p::FlexibleLicense, t::Template, ::AbstractString) = Dict(
    "AUTHORS" => p.owner !== nothing ? p.owner : join(t.authors, ", "),
    "YEAR" => p.format_year(year(today())),
)

function prompt(::Type{FlexibleLicense}, ::Type, ::Val{:name})
    options = readdir(default_file("licenses"))
    pushfirst!("UNLICENSED")
    # Move MIT to the top.
    deleteat!(options, findfirst(==("MIT"), options))
    pushfirst!(options, "MIT")
    menu = RadioMenu(options; pagesize=length(options))
    println("Select a license:")
    idx = request(menu)
    return options[idx]
end

customizable(::Type{FlexibleLicense}) = (:name => String,)

end # module


text:
module Pixell

using WCS
using WCS_jll
import WCS: AbstractWCSTransform
using FITSIO
using FFTW
using Printf
import Unitful, UnitfulAngles
import Unitful: uconvert, ustrip
using StaticArrays
using DSP: unwrap, unwrap!
import FastTransforms: chebyshevjacobimoments1, clenshawcurtisweights
using Libsharp
import Libsharp: sharp_execute!
import Healpix: Alm, map2alm, alm2map, alm2cl

include("enmap.jl")
include("enmap_geom.jl")
include("enmap_ops.jl")
include("transforms.jl")
include("plot.jl")
include("transform_distance.jl")

export Enmap, CarClenshawCurtis, getwcs
export geometry, fullsky_geometry, slice_geometry
export pix2sky, pix2sky!, sky2pix, sky2pix!, skyarea, pixareamap, pixareamap!
export read_map, write_map
export Alm, map2alm, alm2cl, alm2map, alm2map!
export distance_transform, ApproxSeqSDT, ExactSeqSDT, BruteForceSDT

# set up some shortcuts for common angles
const radian = Unitful.rad
const degree = Unitful.°
const arcminute = UnitfulAngles.arcminute

function __init__()
    Enplot.register_colors!()
end

end


text:
# This file was generated, do not modify it. # hide
coefs, intercept = fitted_params(mtm.fitresult.fitresult.machine)
round.(coefs, sigdigits=2)


text:
using Random
using Plots
using Revise
using EvoTrees

# @load "blog/model_linear.bson" model
# @load "data/model_linear_8.bson" model
# @load "data/model_gaussian_5.bson" model

model = EvoTrees.load("data/model_linear_4.bson");
var_names = ["var_$i" for i in 1:100]
plot(model)
plot(model, 2)
plot(model, 3, var_names)
plot(model.trees[2])
plot(model.trees[2], var_names)

typeof(tree_layout[1])
BezierCurve(tree_layout[1])
mutable struct BCurve{T <: GeometryBasics.Point}
    control_points::Vector{T}
end

function (bc::BCurve)(t::Real)
    p = zero(P2)
    n = length(bc.control_points) - 1
    for i in 0:n
        p += bc.control_points[i + 1] * binomial(n, i) * (1 - t)^(n - i) * t^i
    end
    p
end


text:
# Memoization for calls to TauPy
"Raypath cache holding previously-computed ray paths"
const RAY_CACHE = Dict{UInt64,Any}()
"Counter for number of calls to each item in raypath cache"
const RAY_NCALLS = Dict{UInt64,Int64}()
"Maximum size in MB of the raypath cache"
const MAX_SIZE_MB_CACHE = Ref(1024.0)

"Return the key for the ray cache for a particular set of arguments"
cache_args_key(args...) = hash(args)

"""
    update_cache!(key, val)

Update the raypath cache, dropping least-used values when the cache is full, and
return the value in the cache.
"""
function update_cache!(key, val)
    if haskey(RAY_CACHE, key)
        RAY_NCALLS[key] += 1
        return val
    end
    RAY_CACHE[key] = val
    RAY_NCALLS[key] = 1
    while Base.summarysize(RAY_CACHE)/1024^2 > MAX_SIZE_MB_CACHE[]
        minkey = argmin(RAY_NCALLS)
        delete!(RAY_CACHE, minkey)
        delete!(RAY_NCALLS, minkey)
    end
    val
end

"""
    get_cache(key) -> ::Union{Vector{Phase},Vector{PhaseGeog}}

Get a value from the ray path cache and update the call count.
"""
get_cache(key) = update_cache!(key, RAY_CACHE[key])

"""
    clear_cache!()

Remove all values from the ray path cache.
"""
clear_cache!() = (empty!(RAY_CACHE); empty!(RAY_NCALLS))

"""
    set_cache_size_mb!(max_size)

Set the maximum size of the ray path cache.
"""
set_cache_size_mb!(max_size) = MAX_SIZE_MB_CACHE[] = max_size

max_stars_repo_name: anowacki/TauPy


text:
using Statistics

@testset "Auxiliary Functions Test" begin

    @testset "check constant columns" begin

		try
			Perceptrons.check_constant_cols([1.0 1; 1 2;1 3])
		catch
			@test true
		end

		try
			Perceptrons.check_constant_cols([1.0 1 1][:,:])
		catch
			@test true
		end

		try
			Perceptrons.check_constant_cols([1.0 2 3][:,:])
		catch
			@test true
		end

		try
			Perceptrons.check_constant_cols([1.0; 1; 1][:,:])
		catch
			@test true
		end

		@test Perceptrons.check_constant_cols([1.0 1;2 2;3 3])
		@test Perceptrons.check_constant_cols([1.0;2;3][:,:])

	end

	@testset "centralize" begin

		X        = [1; 2 ;3.0][:,:]
		X        = Perceptrons.centralize_data(X,mean(X,dims=1),std(X,dims=1))
		@test all(X .== [-1,0,1.0])

	end

	@testset "decentralize" begin

		Xo        = [1; 2; 3.0][:,:]
		Xn        = [-1,0,1.0][:,:]
		Xn        = Perceptrons.decentralize_data(Xn,mean(Xo,dims=1),std(Xo,dims=1))
		@test all(Xn .== [1; 2; 3.0])

	end

	@testset "checkparams" begin

         try
			 Perceptrons.check_params("linear")
		 catch
			 @test true
		 end
		 try
			 Perceptrons.check_params("x")
		 catch
			 @test true
		 end

	end

	@testset "checkdata" begin

		 try
			 Perceptrons.check_data(zeros(0,0), 0)
		 catch
			 @test true
		 end
		 try
			 Perceptrons.check_data(zeros(1,1), 10)
		 catch
			 @test true
		 end
		 @test Perceptrons.check_data(ones(1,1), 1)

	end

	@testset "check binary labels" begin

		try
			Perceptrons.check_linear_binary_labels([1,0,2])
		catch
			@test true
		end
		try
			Perceptrons.check_linear_binary_labels([1,-1])
		catch
			@test true
		end

		@test Perceptrons.check_linear_binary_labels([1,0])

	end

end;


text:
using Cxx
using Graphs
using GraphViz

import Base: start, next, done
import Graphs: is_directed, num_vertices, vertices,
    implements_vertex_list, implements_adjacency_list, out_neighbors,
    vertex_index

include(Pkg.dir("Cxx","test","llvmincludes.jl"))

immutable LLVMGraph{GraphType,NodeType} <: AbstractGraph{NodeType,Edge{NodeType}}
    g::GraphType
end

function LLVMGraph{GraphType}(g::GraphType)
    LLVMGraph{GraphType,pcpp(@cxx llvm::GraphTraits{$(GraphType)}::NodeType)}(g)
end

is_directed(g::LLVMGraph) = true
println(:(@cxx GraphTraits{$(Expr(:$,:GraphType))}::size(g.g)))
println(macroexpand(:(@cxx GraphTraits{$(Expr(:$,:GraphType))}::size(g.g))))
num_vertices{GraphType}(g::LLVMGraph{GraphType}) = @cxx llvm::GraphTraits{$GraphType}::size(g.g)
implements_vertex_list(::LLVMGraph) = true
implements_vertex_map(::LLVMGraph) = true
implements_adjacency_list(::LLVMGraph) = true

immutable NodeIterator{GT}
    G::GT
end

iterator_type{GraphType}(::Type{NodeIterator{GraphType}}) = @cxx llvm::GraphTraits{$(GraphType)}::nodes_iterator

start{GraphType}(I::NodeIterator{GraphType}) =
    @cxx llvm::GraphTraits{$GraphType}::nodes_begin(I.G)

function next{GraphType}(I::NodeIterator{GraphType},s)
    @assert isa(s, iterator_type(typeof(I)))
    v = icxx"&*$s;"
    icxx"(void)++$s;"
    (v,s)
end

done{GraphType}(I::NodeIterator{GraphType},s) =
    icxx"$s == llvm::GraphTraits<$GraphType>::nodes_end($(I.G));"

vertices{GraphType,NodeType}(g::LLVMGraph{GraphType,NodeType}) =
    NodeIterator{GraphType}(g.g)

immutable EdgeIterator{GT,NT}
    N::NT
end

iterator_type{GT,NT}(::Type{EdgeIterator{GT,NT}}) = @cxx llvm::GraphTraits{$(GT)}::ChildIteratorType

start{GT,NT}(I::EdgeIterator{GT,NT}) =
    @cxx llvm::GraphTraits{$GT}::child_begin(I.N)

function next{GT,NT}(I::EdgeIterator{GT,NT},s)
    @assert isa(s, iterator_type(typeof(I)))
    v = icxx"*$s;"
    icxx"(void)++$s;"
    (v,s)
end

done{GT,NT}(I::EdgeIterator{GT,NT},s) =
    icxx"$s == llvm::GraphTraits<$GT>::child_end($(I.N));"

out_neighbors{GT,NT}(v::NT,g::LLVMGraph{GT,NT}) =
    EdgeIterator{GT,NT}(v)

vertex_index(v,g) = reinterpret(Uint64, v.ptr)

buf = IOBuffer()

code_llvmf(f,t) = pcpp"llvm::Function"(ccall(:jl_get_llvmf, Ptr{Void}, (Any,Any,Bool), f, t, false))
function code_graph(f,args)
    v = @cxx std::string()
    os = @cxx llvm::raw_string_ostream(v)
    graphf = code_llvmf(f,args)
    @cxx llvm::WriteGraph(os,graphf)
    @cxx os->flush()
    bytestring((@cxx v->data()), (@cxx v->length()))
end

gt = code_graph(factorize,(typeof(rand(4,4)),))
lf = code_llvmf(factorize,(typeof(rand(4,4)),))

to_dot(LLVMGraph(lf),buf)
display(Graph(takebuf_string(buf)))


text:
"""
histogram(array; bins = 100)
# Description
return arrays for plotting histogram
"""
function histogram(
    array;
    bins = minimum([100, length(array)]),
    normalize = true,
)
    tmp = zeros(bins)
    down, up = extrema(array)
    down, up = down == up ? (down - 1, up + 1) : (down, up) # edge case
    bucket = collect(range(down, up, length = bins + 1))
    normalization = normalize ? length(array) : 1
    for i in eachindex(array)
        # normalize then multiply by bins
        val = (array[i] - down) / (up - down) * bins
        ind = ceil(Int, val)
        # handle edge cases
        ind = maximum([ind, 1])
        ind = minimum([ind, bins])
        tmp[ind] += 1 / normalization
    end
    return (bucket[2:end] + bucket[1:(end - 1)]) .* 0.5, tmp
end


text:
function check(result::Some)::(Union{Nothing,T} where {T<:Any})

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(something(result), Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object, catch_backtrace())
        @error(error_message)

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end

function process_csv_api_data(api_call_raw_data::String)::DataFrame

    # create a data table from the CSV data -
    tmp_data_table = CSV.read(IOBuffer(api_call_raw_data), DataFrame)

    # sort the table according to the timestamps -
    idx_sort = sortperm(tmp_data_table[:, 1])

    # create a sorted data table -
    sorted_data_table = tmp_data_table[idx_sort, :]

    # return -
    return sorted_data_table
end

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity