text:
# # BSON
# BSON is a efficient compressed json-like data encoding.
# This BSON specification is copied verbatim from [http://bsonspec.org/spec.html](http://bsonspec.org/spec.html).
# The original document uses pseudo-BNF syntax,
# this document uses rearranged but otherwise quite similar syntax
using CombinedParsers
# The [Notes](#notes) contain information for reassembling the parsed part of the byte array into a JSON-like structure.
# This document adds the julia functions to the parser.
# The result `document` is a pure julia BSON parser.
#
# ## Specification Version 1.1
# BSON is a binary format in which zero or more ordered key/value pairs are stored as a single entity. We call this entity a document.

# The following grammar specifies version 1.1 of the BSON standard. We've written the grammar using a `CombinedParsers` syntax.
# Valid BSON data is represented by the `document` non-terminal.

# ### Basic Types
# The following basic types are used as terminals in the rest of the grammar. Each type must be serialized in little-endian format.
byte = Bytes(1,UInt8)         #  (8-bits)
int32 = Bytes(4,Int32)        # (32-bit signed integer, two's complement)
int64 = Bytes(8,Int64)        # (64-bit signed integer, two's complement)
uint64 = Bytes(8,UInt64)      # (64-bit unsigned integer)
double = Bytes(8,Float64)     # (64-bit IEEE 754-2008 binary floating point)
using DecFP
decimal128 = Bytes(16,Dec128) # (128-bit IEEE 754-2008 decimal floating point)

# ### Non-terminals
# The following specifies the rest of the BSON grammar.
# Note that quoted strings represent terminals, and should be interpreted with C semantics
# (e.g. "\x01" represents the byte 0000 0001).
@syntax t = Sequence("\\x", integer_base(16,2,2)) do v
    CharIn(convert(UInt8,v[2]))
end

element = Either{Any}(Any[])
e_list = Repeat(element)
nothing # hide
# Also note that we use the * operator as shorthand for repetition (e.g. ("\x01"*2) is "\x01\x01").
# When used as a unary operator, * means that the repetition can occur 0 or more times.
# !!! note
#     `*` here is replaced with `Repeat`.

# BSON Document. int32 is the total number of bytes comprising the document.
@syntax document = Sequence(int32, e_list, t"\x00")[2]

# Zero or more modified UTF-8 encoded characters followed by '\x00'. The (byte*) MUST NOT contain '\x00', hence it is not full UTF-8.
@syntax cstring = map(String,Repeat_until(byte, t"\x00"))
@syntax e_name = map(Symbol,cstring)            # Key name

# String - The int32 is the number bytes in the (byte\*) + 1 (for the trailing '\x00'). The (byte\*) is zero or more UTF-8 encoded characters.
@inline string_until_before(v) = Bytes{String}(v-1)
@syntax lstring = (after(string_until_before, String, int32) * t"\x00")[1]

# Binary - The int32 is the number of bytes in the (byte*).
@syntax subtype = Either((
    t"\x00",                                  # Generic binary subtype
    t"\x01",                                  # Function
    t"\x02",                                  # Binary (Old)
    t"\x03",                                  # UUID (Old)
    t"\x04",                                  # UUID
    t"\x05",                                  # MD5
    t"\x06",                                  # Encrypted BSON value
    t"\x80",                                  # User defined
))
@syntax binary = after(Vector{UInt8}, int32) do v
    subtype * Bytes(v,Vector{UInt8})
end

#
@syntax code_w_s = int32 * lstring * document # Code w/ scope — Deprecated
for e in [
    :float64 => t"\x01" * e_name * double,               # 64-bit binary floating point
    :string => t"\x02" * e_name * lstring,                # UTF-8 string
    :embedded => t"\x03" * e_name * document,            # Embedded document
    :array => t"\x04" * e_name * document,               # Array
    t"\x05" * e_name * binary,                           # Binary data
    t"\x06" * e_name * (Always() => missing),            # Undefined (value) — Deprecated
    t"\x07" * e_name * Bytes(12,Vector{UInt8}),          # ObjectId
    :boolean => t"\x08" * e_name * Either(
        t"\x00" => false,                                # Boolean "false"
        t"\x01" => true),                                # Boolean "true"
    :UTC => t"\x09" * e_name * map(int64) do v
    v |> Millisecond |> Dates.UTInstant |> DateTime
    end,                                                 # UTC datetime
    t"\x0A" * e_name * (Always() => nothing),            # Null value
    :regex => t"\x0B" * e_name * map(v->Regex(v...),cstring * cstring),
    ## Regular expression - The first cstring is the regex pattern,
    ## the second is the regex options string.
    ## Options are identified by characters, which must be stored in alphabetical order.
    t"\x0C" * e_name * lstring * Bytes(12,Vector{UInt8}), # DBPointer — Deprecated
    t"\x0D" * e_name * lstring,                           # JavaScript code
    t"\x0E" * e_name * lstring,                           # Symbol. — Deprecated
    t"\x0F" * e_name * code_w_s,                         # JavaScript code w/ scope — Deprecated
    t"\x10" * e_name * int32,                            # 32-bit integer
    t"\x11" * e_name * uint64,                           # Timestamp
    :int => t"\x12" * e_name * int64,                    # 64-bit integer
    t"\x13" * e_name * decimal128,                       # 128-bit decimal floating point
    t"\xFF" * e_name,                                    # Min key
    t"\x7F" * e_name,                                    # Max key
]
    pairp = parser(e)
    push!(element,
          map(if result_type(pairp) <: Tuple && fieldcount(result_type(pairp))==3
              v->v[2]=>v[3]
              else
              v->v[2]=>v[3:end]
              end,
              pairp))
end
nothing # hide

# ### Notes
# - Array - The document for an array is a normal BSON document with integer values for the keys, starting with 0 and continuing sequentially. For example, the array ['red', 'blue'] would be encoded as the document {'0': 'red', '1': 'blue'}. The keys must be in ascending numerical order.
# - UTC datetime - The int64 is UTC milliseconds since the Unix epoch.
# - Timestamp - Special internal type used by MongoDB replication and sharding. First 4 bytes are an increment, second 4 are a timestamp.
# - Min key - Special type which compares lower than all other possible BSON element values.
# - Max key - Special type which compares higher than all other possible BSON element values.
# - Generic binary subtype - This is the most commonly used binary subtype and should be the 'default' for drivers and tools.
# - The BSON "binary" or "BinData" datatype is used to represent arrays of bytes. It is somewhat analogous to the Java notion of a ByteArray. BSON binary values have a subtype. This is used to indicate what kind of data is in the byte array. Subtypes from zero to 127 are predefined or reserved. Subtypes from 128-255 are user-defined.
#      - \x02 Binary (Old) - This used to be the default subtype, but was deprecated in favor of \x00. Drivers and tools should be sure to handle \x02 appropriately. The structure of the binary data (the byte* array in the binary non-terminal) must be an int32 followed by a (byte*). The int32 is the number of bytes in the repetition.
#      - \x03 UUID (Old) - This used to be the UUID subtype, but was deprecated in favor of \x04. Drivers and tools for languages with a native UUID type should handle \x03 appropriately.
#      - \x80-\xFF "User defined" subtypes. The binary data can be anything.
# - Code w/ scope - Deprecated. The int32 is the length in bytes of the entire code_w_s value. The string is JavaScript code. The document is a mapping from identifiers to values, representing the scope in which the string should be evaluated.

# ## Testing examples
using Dates
testdata = Dict(:a => 1, :b => 3.0, :c => "yeay", :d => true)

# The [BSON.jl](https://github.com/JuliaIO/BSON.jl) package provides an optimized and feature-rich BSON serializer and parser.
using BSON
bson("test.bson", testdata)

s = read("test.bson")
parse(document,s, log=true)

# Aim of this page is demonstrating the straightforward implementation of a BSON parser.
# The parser is untested and not feature-complete regarding subtypes etc.
# Parsing logic can be plugged into the `element`s with `map`.

# ### Benchmarking
using BenchmarkTools
# BSON has an optimized parser, used for comparison
@benchmark BSON.parse(IOBuffer(s))
# The `CombinedParser` matching is slower.
@benchmark match(document,s)
# Transformation is doing too much work currently - lazy access optimizations are planned.
@benchmark parse(document,s)

# ### Profiling
# Profiling reveals that most time is lost due to the `Either` state type.
## using Profile
## Profile.clear()
## f(n,a...) = for _ in 1:n; match(a...); end
## @profile f(100000,document,s);
## using ProfileView
## ProfileView.view()

# ![](bson-profile.png)
# This is looked into, I have a strong hunch that the allocation can be removed so `CombinedParser` can come close to `BSON.jl`.

# ### Dates
# Seems, reading a document written with BSON works.
# BSON.jl has a nonstandard way with `Dates`...
testdata = Dict(:a => 1, :b => 3.0, :c => "yeay", :d => true, :e => now())
bson("test.bson", testdata)
s = read("test.bson")
parse(document,s, log=true)

# This parser does not do any of the magic in `BSON.jl` to reassemble into julia types.
# BSON dates work though:
let d = now()
    Any[:e => d] ==
        parse(document,UInt8[0x10,0x0,0x0,0x0, # we actually ignore the 4 byte length of the vector
                             0x09,0x65,0x0,reinterpret(UInt8,[Dates.value(d)])...,0x0], log=true)
end


text:
using Printf
using LinearAlgebra
using Random

include("PWGrid_v04.jl")
include("read_kpts.jl")
include("../common/gen_lattice.jl")

function test_main()
    LatConst = 10.
    LatVecs = gen_lattice_hexagonal( LatConst )
    ecutwfc_Ry = 20.0
    kpts_red = read_kpts("KPATH_HCP_60")

    pw = PWGrid( ecutwfc_Ry*0.5, LatVecs, kpts_red )

    @printf("Number of G-vectors: %d\n", pw.gvec.Ng)
    @printf("Number of Gk-vectors\n")
    Ngw = pw.gkvec.Ngw
    for ik = 1:length(Ngw)
        @printf("%3d %8d\n", ik, Ngw[ik])
    end

    Ns = pw.Ns
    Npoints = prod(Ns) # note that for the current implementation Ns = Ng
    Ngwx = pw.gkvec.Ngwx

    @printf("Ngwx = %d\n", pw.gkvec.Ngwx)
    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])

end

test_main()

ext: jl


text:
module PlainBox
# lines
const horizontal = '\u2500'
const vertical = '\u2502'

# elbows
const up_right_corner = '\u250C'
const up_left_corner = '\u2510'
const down_right_corner = '\u2514'
const down_left_corner = '\u2518'

# t intersections
const vertical_right_t = '\u251C'
const vertical_left_t  = '\u2524'
const horizontal_down_t = '\u252C'
const horizontal_up_t = '\u2534'

const cross = '\u253C'

@doc """
Thin line unicode box drawing characters arranged for easy access.

Elements are arranged so that converting booleans of whether or not the
left, right, top or bottom neighbouring box drawing characters are present
into a linear index selects the correct box drawing character.

Example
=======
```julia
left =   true
right=   false
top  =   false
bottom = true

simple[bottom + 2*top + 4*right + 8*left] # == up_left_corner
```
""" ->
const simple =
[
    ' ',
    ' ',
    ' ',
    vertical,
    ' ',
    up_right_corner,
    down_right_corner,
    vertical_right_t,
    ' ',
    up_left_corner,
    down_left_corner,
    vertical_left_t,
    horizontal,
    horizontal_down_t,
    horizontal_up_t,
    cross
]
end


text:
"""
Ur fileformat
See EUROCONTROL NEST Manual for Ur fileformat description

Example Ur fileformat
AZ	2020/02/01	2020/02/29	791			Portugal S M
EB	2020/02/01	2020/02/29	9114			Belg.-Luxembourg
ED	2020/02/01	2020/02/29	6374			Germany
EE	2020/02/01	2020/02/29	3151			Estonia
EF	2020/02/01	2020/02/29	4366			Finland
EG	2020/02/01	2020/02/29	6512	0.848504	   GBP	United Kingdom
EH	2020/02/01	2020/02/29	6739			Netherlands
EI	2020/02/01	2020/02/29	2461			Ireland
EK	2020/02/01	2020/02/29	5759	7.47133	   DKK	Denmark

Example Ur fileformat
AZ 2017/12/01 2017/12/31 1089
EB 2017/12/01 2017/12/31 6550
ED 2017/12/01 2017/12/31 8268
EF 2017/12/01 2017/12/31 5632
EG 2017/12/01 2017/12/31 9381
EH 2017/12/01 2017/12/31 6709
EI 2017/12/01 2017/12/31 2976
"""
module Ur

export read

# include("utility.jl")
using ..util
using Dates

const yyyymmdd = DateFormat("yyyy/mm/dd")

struct UnitRate
    start_date::Date
    end_date::Date
    unitrate::Float64
    valutaconversion::Float64
    valuta::String
    country::String
    function UnitRate(line::String)
        line_elements = split(line, '\t')
        start_date = format_date(line_elements[2], yyyymmdd)
        end_date = format_date(line_elements[3], yyyymmdd)
        unitrate = parse(Float64, line_elements[4]) / 100.0
        if length(line_elements[5]) > 0
            valutaconversion = parse(Float64, line_elements[5])
            valuta = strip(line_elements[6])
            country = strip(line_elements[7])
        else
            valutaconversion = 1.0
            valuta = "EUR"
            country = strip(line_elements[7])
        end
        new(start_date, end_date, unitrate, valutaconversion, valuta, country)
    end
end

function read(file)
    ur = Dict{String, UnitRate}()
    for line in eachline(file)
        ur[line[1:2]] = UnitRate(line)
    end
    return ur
end

end # module


text:
using DataStructures

"""
Hybrid Actor-critic PF
----------------------

Hybrid model : Actor-critic & Value iteration with Particle filter

This algorithm is a weighted sum of Actor-critic and a model-based algorithm
that implements value iteration and estimates the
transition matrix via a particle filter.
In other words, the model-based structure is the same as in
Forward Learner, but now the world-model is estimated optimally and not via a
delta rule.

See learner/forwardlearnerBF.jl, learner/actorcriticlearner.jl
and learner/hybridlearner.jl for the core implementation.

Note 1
------
With particle filtering we need to run multiple times with different seeds.
As final LL we report the MEAN across seeds.
For signals/results we keep track of the ones of all seeds
and in the end report the one with the MAXIMUM LL.

Note 2
------
get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl

Note 3
------
Temperature of softmax is removed, because this learner in NOT convex and weights
can be in (0, +Inf). The temperature is absorbed in the weights. Otherwise we
it is overparametrized.
"""
function getLL_Hybrid_AC_FWD_BF(data::Array{Int,2},
                            alpha::Float64, eta_critic::Float64,
                            gamma::Float64, lambda_actor::Float64,
                            lambda_critic::Float64,
                            surpriseprobability::Float64,
                            stochasticity::Float64,
                            nparticles::Union{Int, Float64},
                            w_MB::Float64,
                            w_MF::Float64;
                            particlefilterseeds = [1,2,3],
                            init_R = 0.,
                            init_policy=nothing,
                            init_Vs = nothing,
                            use_replacing_trace=true,
                            do_debug_print = false)

    nparticles = converttoInt(nparticles)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_Hybrid_AC_FWD supports only one subject")
    end
    all_puzzles = get_all_puzzles(data)

    action_selection_probs_seeds = Array{Array{Float64, 1}, 1}(undef, length(particlefilterseeds))
    [action_selection_probs_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    LL_seeds = zeros(Float64, length(particlefilterseeds))

    Tsas_seeds = Array{Any, 1}(undef, length(particlefilterseeds))

    # signals:
    TD_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [TD_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    Surprise_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [Surprise_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    V_state_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [V_state_seeds[i_seed] = zeros(Float64, nr_samples + data_stats.nr_rewards) for i_seed in 1:length(particlefilterseeds)]
    γ_surprise_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [γ_surprise_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]
    SPE_seeds =Array{Array{Float64,1}, 1}(undef, length(particlefilterseeds))
    [SPE_seeds[i_seed] = zeros(Float64, nr_samples) for i_seed in 1:length(particlefilterseeds)]

    final_policyparameters_record_seed = Array{Dict{Int,Any}, 1}(undef, length(particlefilterseeds))
    [final_policyparameters_record_seed[i_seed] = Dict{Int,Any}() for i_seed in 1:length(particlefilterseeds)]
    finalQSA_Hybrid_record_seed = Array{Dict{Int,Any}, 1}(undef, length(particlefilterseeds))
    [finalQSA_Hybrid_record_seed[i_seed] = Dict{Int,Any}() for i_seed in 1:length(particlefilterseeds)]
    finalQSA_MB_record_seed = Array{Dict{Int,Any}, 1}(undef, length(particlefilterseeds))
    [finalQSA_MB_record_seed[i_seed] = Dict{Int,Any}() for i_seed in 1:length(particlefilterseeds)]

    for i_seed in 1:length(particlefilterseeds)
        learnerMB = []
        learnerMF = []
        learnerHybrid = []
        modulatingfactor = 0.
        action_counter = 0
        for i_puzzle in all_puzzles
            # Puzzle = independent block of the experiment.
            # The code was made for general purposes, in this experiment we
            # have only ONE puzzle.
            learnerMB = ForwardLearnerBF(nr_states = nr_states,
                            nr_actions = nr_actions,
                            γ = gamma,
                            init_R = 0.,
                            nparticles = nparticles,
                            surpriseprobability = surpriseprobability,
                            stochasticity = stochasticity,
                            seedlearner = particlefilterseeds[i_seed],
                            )
            if init_policy == nothing # providing init_policy w/o init_Vs is considered a bug
                learnerMF = ActorCritic(nr_states = nr_states,
                                nr_actions = nr_actions,
                                α_critic = eta_critic,
                                α_actor = alpha,
                                γ = gamma,
                                λ_critic = lambda_critic,
                                λ_actor = lambda_actor,
                                use_replacing_trace = use_replacing_trace,
                                do_debug_print = do_debug_print)
            else
                learnerMF = ActorCritic(nr_states = nr_states,
                                nr_actions = nr_actions,
                                α_critic = eta_critic,
                                α_actor = alpha,
                                γ = gamma,
                                λ_critic = lambda_critic,
                                λ_actor = lambda_actor,
                                policyparameters = copy(init_policy[current_puzzle_id]),
                                V = copy(init_Vs[current_puzzle_id]),
                                use_replacing_trace = use_replacing_trace,
                                do_debug_print = do_debug_print)
            end
            learnerHybrid = Hybrid(nr_states = nr_states,
                            nr_actions = nr_actions,
                            w_MB = w_MB,
                            w_MF = w_MF
                            )
            all_epi_in_puzzle = get_all_episodes_in_puzzle(data, i_puzzle)
            for i_episode in all_epi_in_puzzle
                episode_data = get_data_by_puzzle_episode(data, i_puzzle, i_episode)
                path_length = size(episode_data)[1]  # = T
                if do_debug_print
                    println("i_puzzle:$i_puzzle | i_episode:$i_episode | path_length:$path_length")
                end
                isendofexperiment = false
                for t in 1:path_length
                    action_counter += 1
                    state = episode_data[t, Int(c_state)]
                    action = episode_data[t, Int(c_action)]
                    reward = episode_data[t, Int(c_reward)]
                    next_state = episode_data[t, Int(c_nextState)]
                    done = reward>0 ? true : false
                    if do_debug_print
                        # println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | p(s,a):$(round(p_action,2)) | p(s,1;+-bias): $(round(p1_pos,2))/$(round(p1_neg,2))")
                        println("t:$t | s:$state | a:$action | R:$reward | s':$next_state | done: $done")
                    end

                    # Hybrid learner
                    update!(learnerHybrid, state, action, reward, learnerMB.QSA, learnerMF.policyparameters) #Q_hybrid =  w_MB*learnerMB.QSA + w_MF*learnerMF.QSA

                    p_actions = getSoftmaxProbs(learnerHybrid.QSA[state, :], 1., [0., 0.])

                    action_selection_probs_seeds[i_seed][action_counter] = copy(p_actions[action])
                    LL_seeds[i_seed] += log(p_actions[action])

                    # ForwardLearner-like update (ie Value iteration),
                    # but with model estimated via particle filter
                    update!(learnerMB, state, action, next_state, reward)
                    γ_surprise_seeds[i_seed][action_counter] = copy(learnerMB.Testimate.γ_surprise[1])
                    Surprise_seeds[i_seed][action_counter] = copy(learnerMB.Testimate.Sgm[1])
                    SPE_seeds[i_seed][action_counter] = copy(learnerMB.Testimate.SPE[1])

                    # Actor-Critic
                    update!(learnerMF, state, action, next_state, reward, p_actions[action], modulatingfactor)
                    TD_seeds[i_seed][action_counter] = copy(learnerMF.td_error[1])

                end # end t
            end # end i_episode
            final_policyparameters_record_seed[i_seed][i_puzzle] = copy(learnerMF.policyparameters)
            finalQSA_Hybrid_record_seed[i_seed][i_puzzle] = copy(learnerHybrid.QSA)
            finalQSA_MB_record_seed[i_seed][i_puzzle] = copy(learnerMB.QSA)
        end # end i_puzzle
        Tsas_seeds[i_seed]= deepcopy(learnerMB.Testimate.Ps1a0s0)
    end # seeds

    LL = get_meanLL_acrossseeds(LL_seeds)

    action_selection_probs_cat = hcat(action_selection_probs_seeds...)
    action_selection_probs = mean(action_selection_probs_cat, dims=2)[:]
    trial_by_trial_p_sa = (action_selection_probs)

    # For signals: return the ones corresponding to the Max LL across seeds
    LL_max, seed_idx = findmax(LL_seeds)
    TD_record = copy(TD_seeds[seed_idx])
    SPE_record = copy(SPE_seeds[seed_idx])
    trial_by_trial_p_sa_max = (action_selection_probs_seeds[seed_idx])

    algo = "hybrid_actorcritic_FWD_BF"

    param_v = [alpha, eta_critic, gamma, lambda_actor, lambda_critic, surpriseprobability, stochasticity, nparticles, w_MB, w_MF]
    param_str = "alpha=$(round(alpha,digits=4)),
                eta_critic=$(round(eta_critic,digits=4)),
                gamma=$(round(gamma,digits=4)),
                lambda_actor=$(round(lambda_actor,digits=4)),
                lambda_critic=$(round(lambda_critic,digits=4)),
                surpriseprobability=$(round(surpriseprobability,digits=4)),
                stochasticity=$(round(stochasticity,digits=4)),
                nparticles=$nparticles,
                w_MB=$(round(w_MB,digits=4)),
                w_MF=$(round(w_MF,digits=4))"

    results_dict = Dict{String,Any}("finalQSA_Hybrid_record" => finalQSA_Hybrid_record_seed[seed_idx],
                                    "final_PolicyParameters_MF" => final_policyparameters_record_seed[seed_idx],
                                    "finalQSA_MB_record" => finalQSA_MB_record_seed[seed_idx],
                                    "final_Tsas" => Tsas_seeds[seed_idx],
                                    "Surprise" => Surprise_seeds[seed_idx],
                                    "gamma_surprise" => γ_surprise_seeds[seed_idx],
                                    "nparticles" => nparticles,
                                    "particlefilterseeds" => particlefilterseeds,
                                    "trial_by_trial_p_sa_max" => trial_by_trial_p_sa_max)

    signals = Step_by_step_signals(data_stats.all_subject_ids[1],
    TD_record, #RPE
    SPE_record, #SPE
    trial_by_trial_p_sa,
    LL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
    return LL, signals
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/11

# %%
versioninfo()

# %%
# resized original version

height, width = 187, 746;
#org_sized = rand(Float32, (2001, 2001)) * 60;
org_sized = rand(Float32, (401, 401)) * 60;
#org_sized = rand(Float32, (101, 101)) * 60;
shadow_time_hrs = zeros(Float32, size(org_sized));
height_mat = rand(Float32, (height, width)) * 100; # originally values getting larger from (0, width//2) to the outside with the distance squared

angle_mat = round.(Int32, 2 .* atand.(0:height-1, (0:width-1)' .-(width/2-1))) .+ 1;

enlarged_size = zeros(eltype(org_sized), size(org_sized, 1) + height, size(org_sized, 2) + width);
enlarged_size[1:size(org_sized, 1), range(Int(width/2), length=size(org_sized, 2))] = org_sized;

weights = rand(Float32, 361)/ 10 .+ 0.01;  # weights originally larger in the middle

function computeSumHours(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    shadowed_segments = zeros(eltype(weights), 361)

    for x2 in 1:width
        for y2 in 1:height
            overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            if overshadowed
                angle = angle_mat[y2, x2]
                if shadowed_segments[angle] == 0.0
                    shadowed_segments[angle] = weights[angle]
                end
            end
        end
    end

    return sum(shadowed_segments)
end

function computeAllLines(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
    for x in 1:size(org_sized, 2) - 1
        for y  in 1:size(org_sized, 1) - 1
            shadow_time_hrs[x, y] = computeSumHours(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
        end
    end
    return shadow_time_hrs
end

@time result = computeAllLines(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
@time result = computeAllLines(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
@time result = computeAllLines(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
r0 = result;

# %%
# revised version

weights = Float32.(weights)
shadowed_segments = zeros(eltype(weights), 361)

function computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    shadowed_segments .= 0

    @inbounds for x2 in 1:width
        for y2 in 1:height
            overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            if overshadowed
                angle = angle_mat[y2, x2]
                if shadowed_segments[angle] == 0.0
                    shadowed_segments[angle] = weights[angle]
                end
            end
        end
    end

    return sum(shadowed_segments)
end

function computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
    for x in 1:size(org_sized, 2) - 1
        for y  in 1:size(org_sized, 1) - 1
            shadow_time_hrs[x, y] = computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
        end
    end
    return shadow_time_hrs
end

@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
r1 = result;

# %%
# revised version 2

weights = Float32.(weights)
shadowed_segments = zeros(eltype(weights), 361)

function computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    shadowed_segments .= 0

    Threads.@threads for x2 in 1:width
        for y2 in 1:height
            @inbounds overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            @inbounds if overshadowed
                angle = angle_mat[y2, x2]
                if shadowed_segments[angle] == 0.0
                    shadowed_segments[angle] = weights[angle]
                end
            end
        end
    end

    return sum(shadowed_segments)
end

function computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
    for x in 1:size(org_sized, 2) - 1
        for y  in 1:size(org_sized, 1) - 1
            shadow_time_hrs[x, y] = computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
        end
    end
    return shadow_time_hrs
end

@show Threads.nthreads()
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
r2 = result;

# %%
# revised version 3

weights = Float32.(weights)
shadowed_segments = zeros(eltype(weights), 361)

function computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    shadowed_segments .= 0

    Threads.@threads for x2 in 1:width
        @inbounds @simd for y2 in 1:height
            overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            if overshadowed
                angle = angle_mat[y2, x2]
                if shadowed_segments[angle] == 0.0
                    shadowed_segments[angle] = weights[angle]
                end
            end
        end
    end

    return sum(shadowed_segments)
end

function computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
    shadowed_segments = zeros(eltype(weights), 361)
    for x in 1:size(org_sized, 2) - 1
        for y  in 1:size(org_sized, 1) - 1
            shadow_time_hrs[x, y] = computeSumHours!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, shadowed_segments)
        end
    end
    return shadow_time_hrs
end

@show Threads.nthreads()
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
@time result = computeAllLines!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights, shadowed_segments)
r3 = result;

# %%
r0 ≈ r1 ≈ r2 ≈ r3

# %%

max_stars_repo_path: 0009/Comparing Numba and Julia for a complex matrix computation.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
#===============================================================================
# Logistic Regression for the Numerai dataset
#
# Author: Carl Cortright
# Date: 1/4/2017
#
===============================================================================#
using ScikitLearn
using DataFrames, DataArrays

@sk_import linear_model: LogisticRegression
@sk_import preprocessing: (LabelBinarizer, StandardScaler)

# Define he training set
train = readtable("/home/carl/Documents/Numerai/Datasets/Jan92017/numerai_training_data.csv")

# Pull out features and targets
features = convert(Array, train[1:50])
targets = convert(Array, train[51])

# Fit a logistic regression
logreg = LogisticRegression(fit_intercept=true)
ScikitLearn.fit!(logreg, features, targets)

accuracy = sum(ScikitLearn.predict(logreg, features) .== targets) / length(targets)
println("accuracy: $accuracy")

max_stars_repo_name: carlcortright/Numerai

max_issues_repo_name: carlcortright/Numerai

max_forks_repo_name: carlcortright/Numerai


text:
# Methods for manipulating MolecularDynamicsTrial objects

include("MolecularDynamicsTrial.jl")
import JSON
import HDF5

# Check the filename to make sure it has a descriptive extension
function mdcheckextension( filename::String, r::MDTrial )
    r::Regex
    if typeof(r) == MolecularDynamicsTrial
        ex = r".mdv"                                    # a verlet trial
    elseif typeof(r) == MDMetropolisTrial
        ex = r".mdm"                                    # a metropolis trial
    end
    ismatch(ex, filename) || error("Needs descriptive extension: ", filename)
end

# Write the trial to JSON
function writejson( filename::String, r::MDTrial )
    mdcheckextension(filename, r)
    isfile( filename ) && rm( filename )
    file = open( filename, "w" )
    write( file, JSON.json( r, 4 ) )
    close( file )
end

# Save the trial using HDF5
function save( filename::String, r::MDTrial )
    mdcheckextension(filename, r)
    print("Saving trial as ",filename,"...")
    fileid = HDF5.h5open( filename, "w" )               # open, overwrite
    for field in names( MolecularDynamicsTrial )
        fileid[ string( field ) ] = r.( field )
    end
    close( fileid )
    println(" done.")
end

# Load the trial using HDF5
function load!( filename::String, r::MolecularDynamicsTrial )
    print("Loading trial from ",filename,"...")
    fileid = HDF5.h5open( filename, "r" )               # open, read-only
    for field in names( fileid )
        r.( symbol( field ) ) = HDF5.read( fileid, field )
    end
    close( fileid )
    println(" done.")
end

# Initialize and load a trial using HDF5
function loadmdtrial( filename::String )
    r = MolecularDynamicsTrial()
    load!( filename, r )
    return r
end


text:
"""
    Base.write(filename::AbstractString, A::AbstractRaster; kw...)

Write an [`AbstractRaster`](@ref) to file, guessing the backend from the file extension.

Keyword arguments are passed to the `write` method for the backend.
"""
function Base.write(
    filename::AbstractString, A::AbstractRaster; kw...
)
    write(filename, _sourcetype(filename), A; kw...)
end
Base.write(A::AbstractRaster) = write(filename(A), A)

"""
    Base.write(filename::AbstractString, s::AbstractRasterStack; suffix, kw...)

Write any [`AbstractRasterStack`](@ref) to file, guessing the backend from the file extension.

## Keywords

- `suffix`: suffix to append to file names. By default the layer key is used.

Other keyword arguments are passed to the `write` method for the backend.

If the source can't be saved as a stack-like object, individual array layers will be saved.
"""
function Base.write(filename::AbstractString, s::AbstractRasterStack; suffix=nothing, ext=nothing, kw...)
    if isnothing(ext)
        base, ext = splitext(filename)
    else
        base = filename
    end
    T = _sourcetype(filename)
    if haslayers(T)
        write(filename, _sourcetype(filename), s; kw...)
    else
        # Otherwise write separate files for each layer
        suffix1 = if suffix === nothing
            divider = Sys.iswindows() ? '\\' : '/'
            # Add an underscore to the key if there is a file name already
            spacer = last(filename) == divider ? "" : "_"
            map(k -> string(spacer, k), keys(s))
        else
            suffix
        end
        @warn string("Cannot write stacks to \"", ext, "\", writing layers as individual files")
        map(keys(s), suffix1) do key, sfx
            fn = string(base, sfx, ext)
            write(fn, _sourcetype(filename), s[key])
        end
    end
end

# Trait for source data that has stack layers
haslayers(T) = false


text:
using VariationalImaging.SumRegsDenoise
using AlgTools.LinOps

###################
# Learning function
###################
function sumregs_learning_function(x::AbstractVector{Float64},data,Δ;Δt=1e-3)
	op₁ = FwdGradientOp()
	op₂ = BwdGradientOp()
	op₃ = CenteredGradientOp()
    u = sumregs_denoise(data[2],x,op₁,op₂,op₃)
    cost = 0.5*norm₂²(u-data[1])
	if Δ > Δt
    	grad = sumregs_gradient(x,op₁,op₂,op₃,u,data[1])
	else
		grad = sumregs_gradient_reg(x,op₁,op₂,op₃,u,data[1])
	end
    return u,cost,grad
end

function sumregs_learning_function(x::AbstractArray{T,3},data,Δ;Δt=1e-3) where T
	op₁ = FwdGradientOp()
	op₂ = BwdGradientOp()
	op₃ = CenteredGradientOp()
	ū = data[1]
	pOp = PatchOp(x[:,:,1],ū[:,:,1])
    u = sumregs_denoise(data[2],x,op₁,op₂,op₃,pOp)
    cost = 0.5*norm₂²(u-ū)
	if Δ > Δt
    	grad = sumregs_gradient(x,op₁,op₂,op₃,pOp,u,ū)
	else
		grad = sumregs_gradient_reg(x,op₁,op₂,op₃,pOp,u,ū)
	end
    return u,cost,grad
end

function sumregs_denoise(data,x::AbstractVector{Float64},op₁::LinOp,op₂::LinOp,op₃::LinOp)
	denoise_params = (
        ρ = 0,
        α₁ = x[1],
        α₂ = x[2],
        α₃ = x[3],
		op₁ = op₁,
		op₂ = op₂,
		op₃ = op₃,
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 5000,
        verbose_iter = 5001,
        save_iterations = false
    )
    st_opt, iterate_opt = initialise_visualisation(false)
    opt_img = sumregs_denoise_pdps(data; iterate=iterate_opt, params=denoise_params)
    finalise_visualisation(st_opt)
    return opt_img
end

function sumregs_denoise(data,x::AbstractArray{T,3},op₁::LinOp,op₂::LinOp,op₃::LinOp,pOp::PatchOp) where T
    x̄ = pOp(x)
	denoise_params = (
        ρ = 0,
        α₁ = x̄[:,:,1],
        α₂ = x̄[:,:,2],
        α₃ = x̄[:,:,3],
		op₁ = op₁,
		op₂ = op₂,
		op₃ = op₃,
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 5000,
        verbose_iter = 5001,
        save_iterations = false
    )
    st_opt, iterate_opt = initialise_visualisation(false)
    opt_img = sumregs_denoise_pdps(data; iterate=iterate_opt, params=denoise_params)
    finalise_visualisation(st_opt)
    return opt_img
end

function sumregs_gradient(x::AbstractVector{Float64},op₁::LinOp,op₂::LinOp,op₃::LinOp,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T
	M,N,O = size(u)
	grad = zeros(size(x))
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = sumregs_gradient(x,op₁,op₂,op₃,u1,u2)
		grad += g
	end
	return grad
end

function sumregs_gradient_reg(x::AbstractVector{Float64},op₁::LinOp,op₂::LinOp,op₃::LinOp,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T
	@info "Using reg gradient"
	M,N,O = size(u)
	grad = zeros(size(x))
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = sumregs_gradient_reg(x,op₁,op₂,op₃,u1,u2)
		grad += g
	end
	return grad
end

function sumregs_gradient_reg(x::AbstractVector{Float64},op₁::LinOp,op₂::LinOp,op₃::LinOp,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	n = size(u,1)
	γ = 1e3

	G₁ = matrix(op₁,n)
	Gu₁ = G₁*u[:]
	nGu₁ = xi(Gu₁)
	act1₁ = nGu₁ .- 1/γ
	act₁ = max.(0,act1₁) .!= 0
	inact₁ = 1 .- act₁
	Act₁ = spdiagm(0=>act₁)
	Inact₁ = spdiagm(0=>inact₁)
	den₁ = Act₁*nGu₁ + inact₁
	Den₁ = spdiagm(0=>1 ./den₁)
	prodGuGu₁ = prodesc(Gu₁./(den₁.^3),Gu₁)
	I = spdiagm(0=>ones(n^2))
	B₁ = γ*Inact₁
	C₁ = (Act₁*(prodGuGu₁-Den₁))

	G₂ = matrix(op₂,n)
	Gu₂ = G₂*u[:]
	nGu₂ = xi(Gu₂)
	act1₂ = nGu₂ .- 1/γ
	act₂ = max.(0,act1₂) .!= 0
	inact₂ = 1 .- act₂
	Act₂ = spdiagm(0=>act₂)
	Inact₂ = spdiagm(0=>inact₂)
	den₂ = Act₂*nGu₂ + inact₂
	Den₂ = spdiagm(0=>1 ./den₂)
	prodGuGu₂ = prodesc(Gu₂./(den₂.^3),Gu₂)
	I = spdiagm(0=>ones(n^2))
	B₂ = γ*Inact₂
	C₂ = (Act₂*(prodGuGu₂-Den₂))

	G₃ = matrix(op₃,n)
	Gu₃ = G₃*u[:]
	nGu₃ = xi(Gu₃)
	act1₃ = nGu₃ .- 1/γ
	act₃ = max.(0,act1₃) .!= 0
	inact₃ = 1 .- act₃
	Act₃ = spdiagm(0=>act₃)
	Inact₃ = spdiagm(0=>inact₃)
	den₃ = Act₃*nGu₃ + inact₃
	Den₃ = spdiagm(0=>1 ./den₃)
	prodGuGu₃ = prodesc(Gu₃./(den₃.^3),Gu₃)
	I = spdiagm(0=>ones(n^2))
	B₃ = γ*Inact₃
	C₃ = (Act₃*(prodGuGu₃-Den₃))

	p = (I+x[1]*G₁'*(B₁-C₁)*G₁+x[2]*G₂'*(B₂-C₂)*G₂+x[3]*G₃'*(B₃-C₃)*G₃)\(ū[:]-u[:])
	grad = [p'*(G₁'*(Act₁*Den₁*Gu₁+γ*Inact₁*Gu₁));p'*(G₂'*(Act₂*Den₂*Gu₂+γ*Inact₂*Gu₂));p'*(G₃'*(Act₃*Den₃*Gu₃+γ*Inact₃*Gu₃))]
	return grad
end

function sumregs_gradient(x::AbstractArray{T,3},op₁::LinOp,op₂::LinOp,op₃::LinOp,pOp::PatchOp,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T

	M,N,O = size(u)
	grad = zeros(size(x))
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = sumregs_gradient(x,op₁,op₂,op₃,pOp,u1,u2)
		grad += g
	end
	return grad
end

function sumregs_gradient_reg(x::AbstractArray{T,3},op₁::LinOp,op₂::LinOp,op₃::LinOp,pOp::PatchOp,u::AbstractArray{T,3},ū::AbstractArray{T,3}) where T

	M,N,O = size(u)
	grad = zeros(size(x))
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = sumregs_gradient_reg(x,op₁,op₂,op₃,pOp,u1,u2)
		grad += g
	end
	return grad
end

function sumregs_gradient_reg(x::AbstractArray{T,3},op₁::LinOp,op₂::LinOp,op₃::LinOp,pOp::PatchOp,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	m,n = size(u)
	γ = 1e8

	I = spdiagm(0=>ones(m*n))

	G₁ = matrix(op₁,m)
	Gu₁ = G₁*u[:]
	nGu₁ = xi(Gu₁)
	act1₁ = nGu₁ .- 1/γ
	act₁ = max.(0,act1₁) .!= 0
	inact₁ = 1 .- act₁
	Act₁ = spdiagm(0=>act₁)
	Inact₁ = spdiagm(0=>inact₁)
	den₁ = Act₁*nGu₁ + inact₁
	Den₁ = spdiagm(0=>1 ./den₁)
	prodGuGu₁ = prodesc(Gu₁./(den₁.^3),Gu₁)
	B₁ = γ*Inact₁
	C₁ = (Act₁*(prodGuGu₁-Den₁))

	G₂ = matrix(op₂,m)
	Gu₂ = G₂*u[:]
	nGu₂ = xi(Gu₂)
	act1₂ = nGu₂ .- 1/γ
	act₂ = max.(0,act1₂) .!= 0
	inact₂ = 1 .- act₂
	Act₂ = spdiagm(0=>act₂)
	Inact₂ = spdiagm(0=>inact₂)
	den₂ = Act₂*nGu₂ + inact₂
	Den₂ = spdiagm(0=>1 ./den₂)
	prodGuGu₂ = prodesc(Gu₂./(den₂.^3),Gu₂)
	B₂ = γ*Inact₂
	C₂ = (Act₂*(prodGuGu₂-Den₂))

	G₃ = matrix(op₃,m)
	Gu₃ = G₃*u[:]
	nGu₃ = xi(Gu₃)
	act1₃ = nGu₃ .- 1/γ
	act₃ = max.(0,act1₃) .!= 0
	inact₃ = 1 .- act₃
	Act₃ = spdiagm(0=>act₃)
	Inact₃ = spdiagm(0=>inact₃)
	den₃ = Act₃*nGu₃ + inact₃
	Den₃ = spdiagm(0=>1 ./den₃)
	prodGuGu₃ = prodesc(Gu₃./(den₃.^3),Gu₃)
	B₃ = γ*Inact₃
	C₃ = (Act₃*(prodGuGu₃-Den₃))

	x₁ = pOp(x[:,:,1])
	x₂ = pOp(x[:,:,2])
	x₃ = pOp(x[:,:,3])

	p = (I + x₁[:] .*G₁'*(B₁-C₁)*G₁ + x₂[:] .*G₂'*(B₂-C₂)*G₂ + x₃[:] .*G₃'*(B₃-C₃)*G₃)\(ū[:]-u[:])
	g₁ = reshape(spdiagm(0=>p)*(G₁'*(Act₁*Den₁*Gu₁+γ*Inact₁*Gu₁)),m,n)
	g₂ = reshape(spdiagm(0=>p)*(G₂'*(Act₂*Den₂*Gu₂+γ*Inact₂*Gu₂)),m,n)
	g₃ = reshape(spdiagm(0=>p)*(G₃'*(Act₃*Den₃*Gu₃+γ*Inact₃*Gu₃)),m,n)

	gx = zeros(pOp.size_in...,3)

	gx[:,:,1] = calc_adjoint(pOp,g₁)
	gx[:,:,2] = calc_adjoint(pOp,g₂)
	gx[:,:,3] = calc_adjoint(pOp,g₃)

	return gx
end

function sumregs_gradient(x::AbstractVector{Float64},op₁::LinOp,op₂::LinOp,op₃::LinOp,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	n = size(u,1)

	# Generate centered gradient matrix fwd
	G₁ = matrix(op₁,n)
	Gu₁ = G₁*u[:]
	nGu₁ = xi(Gu₁)
	act₁ = nGu₁ .< 1e-12
	inact₁ = 1 .- act₁
	Act₁ = spdiagm(0=>act₁)
	Inact₁ = spdiagm(0=>inact₁)

	## Vector with grad norm in inactive components and one in the active
	den₁ = Inact₁*nGu₁+act₁
	Den₁ = spdiagm(0=>1 ./den₁)

	## prod KuKuᵗ/norm³
	prodKuKu₁ = prodesc(Gu₁ ./den₁.^3,Gu₁)

	# Generate centered gradient matrix fwd
	G₂ = matrix(op₂,n)
	Gu₂ = G₂*u[:]
	nGu₂ = xi(Gu₂)
	act₂ = nGu₂ .< 1e-12
	inact₂ = 1 .- act₂
	Act₂ = spdiagm(0=>act₂)
	Inact₂ = spdiagm(0=>inact₂)

	## Vector with grad norm in inactive components and one in the active
	den₂ = Inact₂*nGu₂+act₂
	Den₂ = spdiagm(0=>1 ./den₂)

	## prod KuKuᵗ/norm³
	prodKuKu₂ = prodesc(Gu₂ ./den₂.^3,Gu₂)

	# Generate centered gradient matrix fwd
	G₃ = matrix(op₃,n)
	Gu₃ = G₃*u[:]
	nGu₃ = xi(Gu₃)
	act₃ = nGu₃ .< 1e-12
	inact₃ = 1 .- act₃
	Act₃ = spdiagm(0=>act₃)
	Inact₃ = spdiagm(0=>inact₃)

	## Vector with grad norm in inactive components and one in the active
	den₃ = Inact₃*nGu₃+act₃
	Den₃ = spdiagm(0=>1 ./den₃)

	## prod KuKuᵗ/norm³
	prodKuKu₃ = prodesc(Gu₃ ./den₃.^3,Gu₃)

	Adj = [spdiagm(0=>ones(n^2)) -G₁' -G₂' -G₃';
			Act₁*G₁+Inact₁*x[1]*(Den₁-prodKuKu₁)*G₁ Inact₁+eps()*Act₁ spzeros(2*n^2,2*n^2) spzeros(2*n^2,2*n^2);
			Act₂*G₂+Inact₂*x[2]*(Den₂-prodKuKu₂)*G₂ spzeros(2*n^2,2*n^2) Inact₂+eps()*Act₂ spzeros(2*n^2,2*n^2);
			Act₃*G₃+Inact₃*x[3]*(Den₃-prodKuKu₃)*G₃ spzeros(2*n^2,2*n^2) spzeros(2*n^2,2*n^2) Inact₃+eps()*Act₃]

	Track=[(u[:]-ū[:]);zeros(6*n^2)]
	mult = Adj\Track
	p = @view mult[1:n^2]
	return -[p'*(G₁'*Inact₁*Den₁*Gu₁);p'*(G₂'*Inact₂*Den₂*Gu₂);p'*(G₃'*Inact₃*Den₃*Gu₃)]
end

function sumregs_gradient(x::AbstractArray{T,3},op₁::LinOp,op₂::LinOp,op₃::LinOp,pOp::PatchOp,u::AbstractArray{T,2},ū::AbstractArray{T,2}) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	m,n = size(u)

	# Generate centered gradient matrix fwd
	G₁ = matrix(op₁,n)
	Gu₁ = G₁*u[:]
	nGu₁ = xi(Gu₁)
	act₁ = nGu₁ .< 1e-12
	inact₁ = 1 .- act₁
	Act₁ = spdiagm(0=>act₁)
	Inact₁ = spdiagm(0=>inact₁)

	## Vector with grad norm in inactive components and one in the active
	den₁ = Inact₁*nGu₁+act₁
	Den₁ = spdiagm(0=>1 ./den₁)

	## prod KuKuᵗ/norm³
	prodKuKu₁ = prodesc(Gu₁ ./den₁.^3,Gu₁)

	# Generate centered gradient matrix fwd
	G₂ = matrix(op₂,n)
	Gu₂ = G₂*u[:]
	nGu₂ = xi(Gu₂)
	act₂ = nGu₂ .< 1e-12
	inact₂ = 1 .- act₂
	Act₂ = spdiagm(0=>act₂)
	Inact₂ = spdiagm(0=>inact₂)

	## Vector with grad norm in inactive components and one in the active
	den₂ = Inact₂*nGu₂+act₂
	Den₂ = spdiagm(0=>1 ./den₂)

	## prod KuKuᵗ/norm³
	prodKuKu₂ = prodesc(Gu₂ ./den₂.^3,Gu₂)

	# Generate centered gradient matrix fwd
	G₃ = matrix(op₃,n)
	Gu₃ = G₃*u[:]
	nGu₃ = xi(Gu₃)
	act₃ = nGu₃ .< 1e-12
	inact₃ = 1 .- act₃
	Act₃ = spdiagm(0=>act₃)
	Inact₃ = spdiagm(0=>inact₃)

	## Vector with grad norm in inactive components and one in the active
	den₃ = Inact₃*nGu₃+act₃
	Den₃ = spdiagm(0=>1 ./den₃)

	## prod KuKuᵗ/norm³
	prodKuKu₃ = prodesc(Gu₃ ./den₃.^3,Gu₃)

	x₁ = pOp(x[:,:,1])
	x₂ = pOp(x[:,:,2])
	x₃ = pOp(x[:,:,3])

	Adj = [spdiagm(0=>ones(n^2)) -G₁' -G₂' -G₃';
			Act₁*G₁+Inact₁*spdiagm(0=>[x₁[:];x₁[:]])*(Den₁-prodKuKu₁)*G₁ Inact₁+eps()*Act₁ spzeros(2*n^2,2*n^2) spzeros(2*n^2,2*n^2);
			Act₂*G₂+Inact₂*spdiagm(0=>[x₂[:];x₂[:]])*(Den₂-prodKuKu₂)*G₂ spzeros(2*n^2,2*n^2) Inact₂+eps()*Act₂ spzeros(2*n^2,2*n^2);
			Act₃*G₃+Inact₃*spdiagm(0=>[x₃[:];x₃[:]])*(Den₃-prodKuKu₃)*G₃ spzeros(2*n^2,2*n^2) spzeros(2*n^2,2*n^2) Inact₃+eps()*Act₃]

	Track=[(u[:]-ū[:]);zeros(6*n^2)]
	mult = Adj\Track
	p = @view mult[1:n^2]

	g₁ = reshape(-spdiagm(0=>p[:])*(G₁'*Inact₁*Den₁*Gu₁),m,n)
	g₂ = reshape(-spdiagm(0=>p[:])*(G₂'*Inact₂*Den₂*Gu₂),m,n)
	g₃ = reshape(-spdiagm(0=>p[:])*(G₃'*Inact₃*Den₃*Gu₃),m,n)

	gx = zeros(pOp.size_in...,3)

	gx[:,:,1] = calc_adjoint(pOp,g₁)
	gx[:,:,2] = calc_adjoint(pOp,g₂)
	gx[:,:,3] = calc_adjoint(pOp,g₃)
	return gx
end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round 2009
# Problem A: Alien Language

# Determines whether the pattern pat matches the word w.
function ismatch(pat, w)
	for i = 1:length(w)
		if !in(w[i], pat[i])
			return false
		end
	end
	true
end

# Returns the number of words that the pattern pat matches.
nmatches(pat, words) = count(w -> ismatch(pat, w), words)

# Parses a string a pattern, returning a list of lists of characters. Ordinary
# letters in the pattern result in singleton lists, while letters enclosed in
# parentheses are grouped in the same list.
function parse_pattern(str)
	pat = Vector{Char}[]
	nested = false
	for c in str
		if !nested
			if c == '('
				push!(pat, [])
				nested = true
			else
				push!(pat, [c])
			end
		else
			if c == ')'
				nested = false
			else
				push!(last(pat), c)
			end
		end
	end
	pat
end

function main()
	_, nw, n = int(split(readline()))
	words = map(i -> chomp(readline()), 1:nw)
	for i = 1:n
		pat = parse_pattern(readline())
		nm = nmatches(pat, words)
		println("Case #$i: ", nm)
	end
end

main()


text:
# Generate self-signed certificates for testing

using Dates

using OpenSSL_jll

days_validity = 365*100

openssl = OpenSSL_jll.openssl()

# Required ??
# $opensslcfg

root_ca_cert_file = "root_ca.crt"
root_ca_key_file = "root_ca.key"
root_ca_serial_file = "root_ca.ser"

server_csr_file = "server.csr"
server_key_file = "server.key"
server_cert_file = "server.pem"

# OpenSSL_jll uses nonexistant default config location
opensslcfg = `-config /usr/lib/ssl/openssl.cnf`

# Generate self signed root CA cert + key
run(`$openssl req $opensslcfg
    -nodes
    -days $days_validity
    -x509
    -newkey rsa:2048
    -keyout $root_ca_key_file
    -out $root_ca_cert_file
    -subj "/OU=HttpDownloadsTest/CN=HttpDownloadsCA"
`)

# Generate server private key + cert to be signed
run(`$openssl req $opensslcfg
    -nodes
    -newkey rsa:2048
    -keyout $server_key_file
    -out $server_csr_file
    -subj "/OU=HttpDownloadsTest/CN=localhost"
`)

# Sign the server cert
run(`$openssl x509
    -req -in $server_csr_file
    -days $days_validity
    -CA $root_ca_cert_file
    -CAkey $root_ca_key_file
    -CAcreateserial
    -CAserial $root_ca_serial_file
    -out $server_cert_file
`)

# For the purposes of testing, we only need:
# - server: server cert + private key
# - client: root CA cert
#
# We can remove the other files as we won't be signing any more certificates
# with this CA.
rm(root_ca_key_file)
rm(root_ca_serial_file)
rm(server_csr_file)


text:
push!(LOAD_PATH, pwd())
using Mea.Periodize
using Mea.Transport
using ProfileView

function code_to_profile()
modelvec = Periodize.buildmodelvec("self_ctow0.dat", "statsparams0.json")
Periodize.calcdos2(modelvec)
end

#taken from: https://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/
function benchmark()
    # Any setup code goes here.

    # Run once, to force compilation.
    println("======================= First run:")
    srand(69)
    @time code_to_profile()

    # Run a second time, with profiling.
    println("\n\n======================= Second run:")
    srand(69)
    Profile.init(delay=0.01)
    Profile.clear()
    Profile.clear_malloc_data()
    @profile @time code_to_profile()

    # Write profile results to profile.bin.
    r = Profile.retrieve()
    f = open("profile.bin", "w")
    serialize(f, r)
    close(f)
end

function view_results()
f=open("profile.bin")
r=deserialize(f);
ProfileView.view(r[1], lidict=r[2])

end


text:
push!(LOAD_PATH,"../src/")
using Documenter
using StatsModels, DataFrames, VegaLite
using LinearRegression

makedocs(sitename="LinearRegression.jl", modules = [LinearRegression] ,
    pages = Any[
        "Home" => "index.md",
        "Tutorials" => Any[
            "Basic" => "basic_tutorial.md",
            "Multiple regression" => "multi_tutorial.md",
            "Weighted regression" => "weighted_regression_tutorial.md" ]
    ])

deploydocs(
    repo = "github.com/ericqu/LinearRegression.jl.git",
    push_preview = false,
    devbranch = "main",
    devurl = "dev",
)


text:
function returnSonum8List()
    returnSonumList(8)
end

function returnSonum16List()
    returnSonumList(16)
end

function returnSonumList(nbit::Int)
    if nbit == 8
        return sonum8
    elseif nbit == 16
        return sonum16
    else
        throw(error("Only 8/16bits supported."))
    end
end

function returnSonumBounds(nbit::Int)
    if nbit == 8
        return bounds8
    elseif nbit == 16
        return bounds16
    else
        throw(error("Only 8/16bits supported."))
    end
end


text:
using EmpiricalRisks
using Base.Test
import DualNumbers

### auxiliary functions

function verify_uniloss(loss::UnivariateLoss, f, u::Float64, y::Real)
    # verify inferred types
    for VT in [Float64, Float32]
        @test Base.return_types(EmpiricalRisks.value, Tuple{typeof(loss), VT, VT}) == [VT]
        @test Base.return_types(EmpiricalRisks.deriv, Tuple{typeof(loss), VT, VT}) == [VT]
        @test Base.return_types(EmpiricalRisks.value_and_deriv, Tuple{typeof(loss), VT, VT}) == [Tuple{VT, VT}]
    end

    # verify computation correctness
    r = f(DualNumbers.dual(u, 1.0), y)
    v_r = DualNumbers.realpart(r)
    dv_r = DualNumbers.epsilon(r)

    @test_approx_eq v_r value(loss, u, y)
    @test_approx_eq dv_r deriv(loss, u, y)

    (v, dv) = value_and_deriv(loss, u, y)
    @test_approx_eq v_r v
    @test_approx_eq dv_r dv
end

function verify_uniloss(loss::UnivariateLoss, f, us::AbstractVector, ys::AbstractVector)
    for u in us, y in ys
        verify_uniloss(loss, f, u, y)
    end
end

### test cases

# AbsLoss

_abs(u::DualNumbers.Dual) = DualNumbers.realpart(u) == 0.0 ? DualNumbers.dual(0.0, 0.0) : DualNumbers.abs(u)

verify_uniloss(AbsLoss(),
    (p, y) -> _abs(p - y), -3.0:3.0, -1.0:0.5:1.0)

# SqrLoss

verify_uniloss(SqrLoss(),
    (p, y) -> abs2(p - y) / 2, -3.0:3.0, -1.0:0.5:1.0)

# QuantileLoss

function _quanlossf(t::Float64, u::DualNumbers.Dual, y)
    DualNumbers.realpart(u) > y ? t * (u - y) :
    DualNumbers.realpart(u) < y ? (1.0 - t) * (y - u) :
    DualNumbers.dual(0.0, 0.0)
end

verify_uniloss(QuantileLoss(0.3), (p, y) -> _quanlossf(0.3, p, y), -2.0:0.5:2.0, -1.0:0.5:1.0)
verify_uniloss(QuantileLoss(0.5), (p, y) -> _quanlossf(0.5, p, y), -2.0:0.5:2.0, -1.0:0.5:1.0)

# EpsilonInsLoss

function _epsinsensf(eps::Float64, u::DualNumbers.Dual, y)
    a = abs(DualNumbers.realpart(u) - y)
    a > eps ? _abs(u - y) - eps : DualNumbers.dual(0.0, 0.0)
end

verify_uniloss(EpsilonInsLoss(0.3), (p, y) -> _epsinsensf(0.3, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)
verify_uniloss(EpsilonInsLoss(0.5), (p, y) -> _epsinsensf(0.5, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)

# HuberLoss

function _huberf(h::Float64, u::DualNumbers.Dual, y)
    a = abs(DualNumbers.realpart(u) - y)
    a > h ? h * _abs(u - y) - 0.5 * h^2 : 0.5 * abs2(u - y)
end

verify_uniloss(HuberLoss(0.3), (p, y) -> _huberf(0.3, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)
verify_uniloss(HuberLoss(0.5), (p, y) -> _huberf(0.5, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)

# HingeLoss

_hingef(u::DualNumbers.Dual, y) = y * DualNumbers.realpart(u) < 1.0 ? 1.0 - y * u : DualNumbers.dual(0.0, 0.0)
verify_uniloss(HingeLoss(), _hingef, -2.0:0.5:2.0, [-1.0, 1.0])

# SquaredHingeLoss

_sqrhingef(u::DualNumbers.Dual, y) = y * DualNumbers.realpart(u) < 1.0 ? (1.0 - y * u).^2 : DualNumbers.dual(0.0, 0.0)
verify_uniloss(SqrHingeLoss(), _sqrhingef, -2.0:0.5:2.0, [-1.0, 1.0])

# SmoothedHingeLoss

function _sm_hingef(h::Float64, u::DualNumbers.Dual, y)
    yu = y * DualNumbers.realpart(u)
    yu >= 1.0 + h ? DualNumbers.dual(0.0, 0.0) :
    yu <= 1.0 - h ? 1.0 - y * u :
    abs2(1.0 + h - y * u) / (4 * h)
end

verify_uniloss(SmoothedHingeLoss(0.2), (p, y) -> _sm_hingef(0.2, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)
verify_uniloss(SmoothedHingeLoss(0.5), (p, y) -> _sm_hingef(0.5, p, y), -2.0:0.25:2.0, -1.0:0.5:1.0)

# LogisticLoss

verify_uniloss(LogisticLoss(),
    (p, y) -> log(1 + exp(-y * p)), -2.0:0.5:2.0, [-1.0, -0.5, 0.5, 1.0])


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using DifferentialEquations

# +
function macroParasiteModelFunction(dY,Y,p,t)
    #Host Birth Rate
    a = p[1]
    #Parasite Influence on host birth rate
    b = p[2]
    #Parasite induced host mortality
    α = p[3]
    #Parasite induced decrease in host reproduction
    β = p[4]
    #Intrinsic death rate of parasites
    μ = p[5]
    #dispersion aggregation parameter
    k = p[6]
    #Rate of production of new free-living stages
    λ = p[7]
    #Death Rate of Free-Living Stages
    γ = p[8]

    #Host Population
    H = Y[1]
    #Parasite Population
    P = Y[2]
    #Infective Stages
    W = Y[3]

dY[1] = (a-b)*H - α*P
dY[2] = β*H*W - (μ + α + b) * P - (α*((P^2)/H)*((k+1)/k))
dY[3] = λ*P - (γ*W) - (β*H*W)
end

# -

par=[1.4,1.05,0.0003,0.01,0.5,0.1,10.0,10.0]
init=[100.0,10.0,10.0]
tspan=(0.0,100.0)

macro_odeProblem = ODEProblem(macroParasiteModelFunction,init,tspan,par)

sol=solve(macro_odeProblem);

using Plots

plot(sol,xlabel="Time",yscale=:log10)

max_stars_repo_name: epirecipes/epiRecipes

max_issues_repo_name: epirecipes/epiRecipes

max_forks_repo_name: epirecipes/epiRecipes


text:
#
# separable.jl --
#
# Separable Multidimensional interpolation.
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2021, Éric Thiébaut.
#

# FIXME: use a table for the innermost dimensions!

function apply!(dst::AbstractArray{T,2},
                ker::Kernel{T,S},
                x1::AbstractVector{T},
                x2::AbstractVector{T},
                src::AbstractArray{T,2}) where {T,S}
    apply!(1, Direct, ker, x1, ker, x2, src, 0, dst)
end

function apply!(dst::AbstractArray{T,2},
                ::Type{P},
                ker::Kernel{T,S},
                x1::AbstractVector{T},
                x2::AbstractVector{T},
                src::AbstractArray{T,2}) where {P<:Union{Direct,Adjoint},T,S}
    apply!(1, P, ker, x1, ker, x2, src, 0, dst)
end

function apply!(dst::AbstractArray{T,2},
                ker1::Kernel{T,S1},
                x1::AbstractVector{T},
                ker2::Kernel{T,S2},
                x2::AbstractVector{T},
                src::AbstractArray{T,2}) where {T,S1,S2}
    return apply!(1, Direct, ker1, x1, ker2, x2, src, 0, dst)
end

function apply!(dst::AbstractArray{T,2},
                ::Type{P},
                ker1::Kernel{T,S1},
                x1::AbstractVector{T},
                ker2::Kernel{T,S2},
                x2::AbstractVector{T},
                src::AbstractArray{T,2}) where {P<:Union{Direct,Adjoint},
                                                T,S1,S2}
    return apply!(1, P, ker1, x1, ker2, x2, src, 0, dst)
end

#------------------------------------------------------------------------------
# Direct operations.

@generated function apply!(α::Real,
                           ::Type{Direct},
                           ker1::Kernel{T,S1},
                           x1::AbstractVector{T},
                           ker2::Kernel{T,S2},
                           x2::AbstractVector{T},
                           src::AbstractArray{T,2},
                           β::Real,
                           dst::AbstractArray{T,2}) where {T,S1,S2}

    # Generate pieces of code.
    J1 = [Symbol(:j1_,s) for s in 1:S1]
    W1 = [Symbol(:w1_,s) for s in 1:S1]
    J2 = [Symbol(:j2_,s) for s in 1:S2]
    W2 = [Symbol(:w2_,s) for s in 1:S2]
    code1 = (Meta.generate_getcoefs(J1, W1, :ker1, :lim1, :(x1[i1])),)
    code2 = (Meta.generate_getcoefs(J2, W2, :ker2, :lim2, :(x2[i2])),
             [:( $(W2[i]) *= alpha ) for i in 1:S2]...)
    expr = Meta.generate_interp_expr(:src, J1, W1, J2, W2)

    quote
        # Get dimensions and limits.
        @assert size(dst) == (length(x1), length(x2))
        n1, n2 = size(dst)
        lim1 = limits(ker1, size(src, 1))
        lim2 = limits(ker2, size(src, 2))

        # Apply direct operator.
        if α == 0
            vscale!(dst, β)
        elseif β == 0
            alpha = convert(T, α)
            @inbounds for i2 in 1:n2
                $(code2...)
                for i1 in 1:n1
                    $(code1...)
                    dst[i1,i2] = $expr
                end
            end
        else
            alpha = convert(T, α)
            beta = convert(T, β)
            @inbounds for i2 in 1:n2
                $(code2...)
                for i1 in 1:n1
                    $(code1...)
                    dst[i1,i2] = $expr + beta*dst[i1,i2]
                end
            end
        end
        return dst
    end
end

#------------------------------------------------------------------------------
# In-place adjoint operation.

@generated function apply!(α::Real,
                           ::Type{Adjoint},
                           ker1::Kernel{T,S1},
                           x1::AbstractVector{T},
                           ker2::Kernel{T,S2},
                           x2::AbstractVector{T},
                           src::AbstractArray{T,2},
                           β::Real,
                           dst::AbstractArray{T,2}) where {T<:AbstractFloat,
                                                           S1,S2}

    # Generate pieces of code.
    J1 = [Symbol(:j1_,s) for s in 1:S1]
    W1 = [Symbol(:w1_,s) for s in 1:S1]
    J2 = [Symbol(:j2_,s) for s in 1:S2]
    W2 = [Symbol(:w2_,s) for s in 1:S2]
    temp = [Symbol(:tmp_,s) for s in 1:S2]
    code1 = [Meta.generate_getcoefs(J1, W1, :ker1, :lim1, :(x1[i1]))]
    code2 = [Meta.generate_getcoefs(J2, W2, :ker2, :lim2, :(x2[i2]))]
    for i2 in 1:S2
        push!(code1, :( $(temp[i2]) = $(W2[i2])*val ))
        for i1 in 1:S1
            push!(code1, :(
                dst[$(J1[i1]),$(J2[i2])] += $(W1[i1])*$(temp[i2])
            ))
        end
    end

    quote
        # Get dimensions and limits.
        @assert size(src) == (length(x1), length(x2))
        n1, n2 = size(src)
        lim1 = limits(ker1, size(dst, 1))
        lim2 = limits(ker2, size(dst, 2))

        # Apply adjoint operator.
        vscale!(dst, β)
        if α != 0
            alpha = convert(T, α)
            @inbounds for i2 in 1:n2
                $(code2...)
                for i1 in 1:n1
                    val = src[i1,i2]*alpha
                    $(code1...)
                end
            end
        end
        return dst
    end
end


text:
# Parameters controlling how a plot appears

const title_font_desc = "'PT Sans','Helvetica Neue','Helvetica',sans-serif"
const label_font_desc = "'PT Sans Caption','Helvetica Neue','Helvetica',sans-serif"

# Choose highlight color by darkening the fill color
default_discrete_highlight_color(fill_color::Color) = RGB(1, 1, 1)

default_discrete_highlight_color(fill_color::TransparentColor) = RGBA{Float32}(
        default_discrete_highlight_color(color(fill_color)), fill_color.alpha)

function default_continuous_highlight_color(fill_color::Color)
    c = convert(LCHab, fill_color)
    return LCHab(max(0, c.l - 40), c.c, c.h)
end

default_continuous_highlight_color(fill_color::TransparentColor) = RGBA{Float32}(
        default_continuous_highlight_color(color(fill_color)), fill_color.alpha)

function default_stroke_color(fill_color::Color)
    fill_color = convert(LCHab, fill_color)
    c = LCHab(fill_color.l, fill_color.c, fill_color.h)
    LCHab(c.l - 15, c.c, c.h)
end

default_stroke_color(fill_color::TransparentColor) = RGBA{Float32}(
        default_stroke_color(color(fill_color)), fill_color.alpha)

function default_lowlight_color(fill_color::Color)
    fill_color = convert(LCHab, fill_color)
    c = LCHab(fill_color.l, fill_color.c, fill_color.h)
    LCHab(90, 20, c.h)
end

function default_lowlight_color(fill_color::TransparentColor)
    @warn "For opacity, use `alpha=[a]`, or use `Theme(alphas=[a])` and/or `Scale.alpha_discrete()`, or use `Scale.alpha_continuous()`"
   RGBA{Float32}(Gadfly.default_lowlight_color(color(fill_color)), fill_color.alpha)
end

# Choose a middle color by darkening the fill color
function default_middle_color(fill_color::Color)
    fill_color = convert(LCHab, fill_color)
    LCHab(fill_color.l + 40, fill_color.c, fill_color.h)
end

default_middle_color(fill_color::TransparentColor) = RGBA{Float32}(
        default_middle_color(color(fill_color)), fill_color.alpha)

get_stroke_vector(::Nothing) = []
get_stroke_vector(vec::AbstractVector) = vec
function get_stroke_vector(linestyle::Symbol)
    ldash = 6 * Compose.mm
    dash = 4 * Compose.mm
    dot = 2 * Compose.mm
    gap = 1 * Compose.mm
    linestyle == :solid && return []
    linestyle == :dash && return [dash, gap]
    linestyle == :dot && return [dot, gap]
    linestyle == :dashdot && return [dash, gap, dot, gap]
    linestyle == :dashdotdot && return [dash, gap, dot, gap, dot, gap]
    linestyle == :ldash && return [ldash, gap]
    linestyle == :ldashdash && return [ldash, gap, dash, gap]
    linestyle == :ldashdot && return [ldash, gap, dot, gap]
    linestyle == :ldashdashdot && return [ldash, gap, dash, gap, dot, gap]
    error("unsupported linestyle: ", linestyle)
end

using DocStringExtensions

"""
$(FIELDS)
"""
@varset Theme begin
    "If the color aesthetic is not mapped to anything, this is the color that is used.  (Color)",
    default_color,         ColorOrNothing,  LCHab(70, 60, 240)

    "Size of points in the point, boxplot, and beeswarm geometries.  (Measure)",
    point_size,            Measure,         0.9mm

    "Minimum size of points in the point geometry.  (Measure)",
    point_size_min,        Measure,         0.45mm,
    "point_size_{min,max} to be deprecated. Instead try the new `Scale.size_discrete2` and/or with `style(discrete_sizemap)`"

    "Maximum size of points in the point geometry.  (Measure)",
    point_size_max,        Measure,         1.8mm,
    "point_size_{min,max} to be deprecated. Instead try the new `Scale.size_discrete2` and/or with `style(discrete_sizemap)`"

    "The function `f` in  [`Scale.size_discrete2`](@ref).",
    discrete_sizemap,      Function,        Scale.default_discrete_sizes

    "The function `f` in  [`Scale.size_radius`](@ref) and [`Scale.size_area`](@ref).",
    continuous_sizemap,      Function,        Scale.default_continuous_sizes

    "A new `Theme` field, in development. Currently only works with `Guide.manual_discrete_key` and `Guide.manual_color_key`. (Function)",
    discrete_colormap,      Function,         Scale.default_discrete_colors

    "Shapes of points in the point geometry.  (Function in circle, square, diamond, cross, xcross, utriangle, dtriangle, star1, star2, hexagon, octagon, hline, vline, ltriangle, rtriangle)",
    point_shapes,          Vector{Function},  [Shape.circle, Shape.square, Shape.diamond, Shape.cross, Shape.xcross,
                                               Shape.utriangle, Shape.dtriangle, Shape.star1, Shape.star2,
                                               Shape.hexagon, Shape.octagon, Shape.hline, Shape.vline, Shape.ltriangle, Shape.rtriangle]

    "Width of lines in the line geometry. (Measure)",
    line_width,            Measure,         0.3mm

    "Style of lines in the line geometry. The default palette is `[:solid, :dash, :dot, :dashdot, :dashdotdot, :ldash, :ldashdash, :ldashdot, :ldashdashdot]` which is a Vector{Symbol}, or customize using Vector{Vector{<:Measure}}",
    line_style,            (Vector{<:Union{Symbol, Vector{<:Measure}}}),   [:solid, :dash, :dot, :dashdot, :dashdotdot, :ldash, :ldashdash, :ldashdot, :ldashdashdot]

    "Alpha palette. The default palette is [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]. Customize using a Vector of length one or greater, with 0.0≤values≤1.0",
    alphas,         Vector{Float64}, [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]

    "Background color used in the main plot panel. (Color or Nothing)",
    panel_fill,            ColorOrNothing,  nothing

    "Border color of the main plot panel. (Color or Nothing)",
    panel_stroke,          ColorOrNothing,  nothing

    "Border line width for main plot panel. (Measure)",
    panel_line_width,      Measure,      0.3mm

    "Opacity of the plot background panel. (Float in [0.0, 1.0])",
    panel_opacity,         Float64,         0.0,
    "The keyword argument `panel_opacity` has been deprecated. Instead, provide a e.g. RGBA() color to panel_fill."

    "Background color for the entire plot. If nothing, no background. (Color or Nothing)",
    background_color,      ColorOrNothing,  nothing

    "Padding around the plot. The order of padding is: `plot_padding=[left, right, top, bottom]`. If a vector of length one is provided e.g.  `[5mm]` then that value is applied to all sides. Absolute or relative units can be used. (Vector{<:Measure})",
    plot_padding,          (Vector{<:Measure}),         [5mm]

    "Color of grid lines. (Color or Nothing)",
    grid_color,            ColorOrNothing,  colorant"#D0D0E0"
    "Style of grid lines. (Symbol in :solid, :dash, :dot, :dashdot, :dashdotdot, or Vector of Measures)",
    grid_line_style,       Union{Symbol,Vector},   [0.5mm, 0.5mm]

    "In the D3 backend, mousing over the plot makes the grid lines emphasised by transitioning to this color. (Color or Nothing)",
    grid_color_focused,    ColorOrNothing,  colorant"#A0A0A0"

    "Width of grid lines. (Measure)",
    grid_line_width,       Measure,         0.2mm

    "Context order of the grid lines (default=0). The plot panel background has a context order of -1. (Int)",
    grid_line_order,        Int,            0

    "Font used for minor labels such as tick labels and entries in keys. (String)",
    minor_label_font,      AbstractString,          label_font_desc
    "Font size used for minor labels. (Measure)",
    minor_label_font_size, Measure,         8pt
    "Color used for minor labels. (Color)",
    minor_label_color,     ColorOrNothing,  colorant"#6c606b"

    "Font used for major labels such as guide titles and axis labels. (String)",
    major_label_font,      AbstractString,          title_font_desc
    "Font size used for major labels. (Measure)",
    major_label_font_size, Measure,         11pt
    "Color used for major labels. (Color)",
    major_label_color,     ColorOrNothing,  colorant"#564a55"

    "Font used for labels in Geom.label. (String)",
    point_label_font,      AbstractString,          label_font_desc
    "Font size used for labels. (Measure)",
    point_label_font_size, Measure,         8pt
    "Color used for labels. (Color)",
    point_label_color,     ColorOrNothing,  colorant"#4c404b"

    "Font used for titles of keys. (String)",
    key_title_font,      AbstractString,          title_font_desc
    "Font size used for key titles. (Measure)",
    key_title_font_size, Measure,         11pt
    "Color used for key titles. (Color)",
    key_title_color,     ColorOrNothing,  colorant"#362a35"

    "Font used for key entry labels. (String)",
    key_label_font,      AbstractString,          title_font_desc
    "Font size used for key entry labels. (Measure)",
    key_label_font_size, Measure,         8pt
    "Color used for key entry labels. (Color)",
    key_label_color,     ColorOrNothing,  colorant"#4c404b"

    "How many gradations to show in a continuous color key. (Int)",
    key_color_gradations, Int,            40

    "Spacing between bars in [`Geom.bar`](@ref). (Measure)",
    bar_spacing,           Measure,         -0.05mm

    "Spacing between boxplots in [`Geom.boxplot`](@ref). (Measure)",
    boxplot_spacing,       Measure,         1mm

    "Length of caps on error bars. (Measure)",
    errorbar_cap_length,   Measure,         3mm

    # SEEMS TO BE ONLY USED BY GEOM.ERRORBAR ???
    # Lines are drawn in a slightly different color than fills, e.g. to
    # differentiate histogram bars from error bars.
    stroke_color,          Function,       default_stroke_color

    "Width of lines drawn around plot geometry like points, and boxplot rectangles. (Measure)",
    highlight_width,       Measure,         0.3mm

    "Color used to outline plot geometry. This is a function that alters (e.g. darkens) the fill color of the geometry. (Function)",
    discrete_highlight_color,       Function,        default_discrete_highlight_color
    "Color used to outline plot geometry. This is a function that alters (e.g. darkens) the fill color of the geometry. (Function)",
    continuous_highlight_color,     Function,        default_continuous_highlight_color

    "Color used to draw background geometry, such as `Geom.ribbon` and `Geom.polygon`. This is a function that alters the fill color of the geometry.  (Function)",
    lowlight_color,        Function,        default_lowlight_color

    "Color altering function used to draw the midline in boxplots. (Function)",
    middle_color,          Function,        default_middle_color

    "Width of the middle line in boxplots. (Measure)",
    middle_width,          Measure,         0.6mm

    "One of `:left`, `:center`, `:right` indicating the placement of the title of color key guides. (Symbol)",
    guide_title_position,  Symbol,          :left

    "The shape used in color swatches in the color key guide. Either `:circle` or `:square`  (Symbol)",
    colorkey_swatch_shape, Symbol,          :square

    # will eventually replace `colorkey_swatch_shape`
    "Shape used in keys for swatches (Function as in `point_shapes`)",
    key_swatch_shape,       Union{Function,Nothing},        nothing

    "Default color used in keys for swatches.  Currently works for `Guide.shapekey` and `Guide.sizekey` (Color)",
    key_swatch_color,       ColorOrNothing,        nothing

    "Size of key swatches, will override `Theme(point_size=)`.  Currently works for `Guide.shapekey`  (Measure)",
    key_swatch_size,       Union{Measure,Nothing},        nothing

    "Where key should be placed relative to the plot panel. One of `:left`, `:right`, `:top`, `:bottom`, `:inside` or `:none`. Setting to `:none` disables the key. Setting to `:inside` places the key in the lower right quadrant of the plot. (Symbol)",
    key_position,          Symbol,          :right

    "Color used to stroke bars in bar plots. If a function is given, it's used to transform the fill color of the bars to obtain a stroke color. (Function, Color, or Nothing)",
    bar_highlight,         Union{Nothing, Function, Color},   nothing

    rug_size,             Measure,          2.0mm

    "Number of annealing iterations.  Used by `Geom.label(position=:dynamic)`",
    label_placement_iterations,  Int,       1000

    "Penalty for a label not being contained within the plot frame.  Used by `Geom.label(position=:dynamic)`",
    label_out_of_bounds_penalty, Float64,   10.0

    "Penalty for making a label hidden to avoid overlaps.  Used by `Geom.label(position=:dynamic)`",
    label_hidden_penalty,        Float64,   0.5

    "Probability of proposing a visibility flip during label layout.  Used by `Geom.label(position=:dynamic)`",
    label_visibility_flip_pr,    Float64,   0.2

    "Padding between marker and label.  Used by `Geom.label(position=:dynamic)`",
    label_padding,    Measure,   1mm

    "Maximum number of columns for key entry labels. (Int)",
    key_max_columns,             Int,       4

    "A `DiscreteColorScale` see [`Scale.color_discrete_hue`](@ref)",
    discrete_color_scale,        Scale.DiscreteColorScale, Scale.color_discrete()

    "A `ContinuousColorScale` see [`Scale.color_continuous`](@ref)",
    continuous_color_scale,      Scale.ContinuousColorScale, Scale.color_continuous()
end

### should we export {current,with,pop,push)_theme?
const theme_stack = Theme[Theme()]

"""
    current_theme()

Get the `Theme` on top of the theme stack.
"""
current_theme() = theme_stack[end]

"""
    style(; kwargs...) -> Theme

Return a new `Theme` that is a copy of the current theme as modifed by the
attributes in `kwargs`.  See [Themes](@ref) for available fields.

# Examples

```
style(background_color="gray")
```
"""
style(; kwargs...) = Theme(current_theme(); kwargs...)

"""
    push_theme(t::Theme)

Set the current theme by placing `t` onto the top of the theme stack.
See also [`pop_theme`](@ref) and [`with_theme`](@ref).
"""
function push_theme(t::Theme)
    push!(theme_stack, t)
    nothing
end

"""
    pop_theme() -> Theme

Return to using the previous theme by removing the top item on the theme stack.
See also [`pop_theme`](@ref) and [`with_theme`](@ref).
"""
function pop_theme()
    length(theme_stack) == 1 && error("There default theme cannot be removed")
    pop!(theme_stack)
end

"""
    push_theme(t::Symbol)

Push a `Theme` by its name.  Available options are `:default` and `:dark`.
See also [`get_theme`](@ref).
"""
push_theme(t::Symbol) = push_theme(get_theme(Val{t}()))

"""
    get_theme()

Register a theme by name by adding methods to `get_theme`.

# Examples

```
get_theme(::Val{:mytheme}) = Theme(...)
push_theme(:mytheme)
```
"""
get_theme(::Val{name}) where {name} = error("No theme $name found")

"""
    with_theme(f, theme)

Call function `f` with `theme` as the current `Theme`.
`theme` can be a `Theme` object or a symbol.

# Examples

```
with_theme(style(background_color=colorant"#888888"))) do
    plot(x=rand(10), y=rand(10))
end
```
"""
function with_theme(f, theme)
    push_theme(theme)
    p = f()
    pop_theme()
    p
end

"""
    get_theme(::Val{:default})

A dark foreground on a light background.
"""
get_theme(::Val{:default}) = Theme()

### Dark theme

const dark_theme = let label_color=colorant"#a1a1a1",
    bgcolor=colorant"#222831",
    grid_color=colorant"#575757",
    fgcol1=colorant"#FE4365",
    fgcol2=colorant"#eca25c",
    fgcol3=colorant"#3f9778",
    fgcol4=colorant"#005D7F"

    function border_color(fill_color)
        fill_color = convert(LCHab, fill_color)
        c = LCHab(fill_color.l, fill_color.c, fill_color.h)
        LCHab(60, 20, c.h)
    end

    function gen_dark_colors(n)
      cs = distinguishable_colors(n, [fgcol1, fgcol2,fgcol3],
          lchoices=Float64[58, 45, 72.5, 90],
          transform=c -> deuteranopic(c, 0.1),
          cchoices=Float64[20,40],
          hchoices=[75,51,35,120,180,210,270,310]
      )

      convert(Vector{Color}, cs)
    end

    function lowlight_color(fill_color)
        fill_color = convert(LCHab, fill_color)
        c = LCHab(fill_color.l, fill_color.c, fill_color.h)
        c2 = convert(RGBA, LCHab(c.l, 50, c.h))
        RGBA(c2.r, c2.g, c2.b, .53)
    end

    function dark_theme_discrete_colors(;
        levels=nothing,
        order=nothing,
        preserve_order=true)

        Gadfly.Scale.DiscreteColorScale(
            gen_dark_colors,
            levels=levels,
            order=order,
            preserve_order=preserve_order
        )
    end

    function dark_theme_continuous_colors()
        Scale.color_continuous(
          colormap=Scale.lab_gradient(fgcol4, fgcol2, fgcol1)
        )
    end

    Theme(
          default_color=fgcol1,
          stroke_color=default_stroke_color,
          panel_fill=bgcolor,
          major_label_color=label_color,
          minor_label_color=label_color,
          grid_color=grid_color,
          key_title_color=label_color,
          key_label_color=label_color,
          lowlight_color=lowlight_color,
          background_color=bgcolor,
          discrete_highlight_color=border_color,
          discrete_color_scale=dark_theme_discrete_colors(),
          continuous_color_scale=dark_theme_continuous_colors(),
    )
end

"""
    get_theme(::Val{:dark})

A light foreground on a dark background.
"""
get_theme(::Val{:dark}) = dark_theme


text:
import GPUArrays: is_gpu, is_cpu, name, threads, blocks, global_memory
import GPUArrays: supports_double, local_memory

function devices(filter_funcs...)
    sort(filter(cl.devices()) do dev
        # These drivers have some problems -  TODO figure out if we are the problem
        !contains(cl.info(dev, :version), "(Build 10)") && # intels experimental 2.1 driver
        !any(f-> !f(dev), filter_funcs)
    end, by = x-> !is_gpu(x)) # now gpu devices come first
end

is_gpu(dev::cl.Device) = cl.info(dev, :device_type) == :gpu
is_cpu(dev::cl.Device) = cl.info(dev, :device_type) == :cpu

name(dev::cl.Device) = string("CL ", cl.info(dev, :name))
function cl_version(dev::cl.Device)
    ver_str = cl.info(dev, :version)
    vmatch = match(r"\d.\d", ver_str)
    major, minor = parse.(Int, split(vmatch.match, '.'))
    VersionNumber(major, minor)
end

threads(dev::cl.Device) = cl.info(dev, :max_work_group_size) |> Int
blocks(dev::cl.Device) = cl.info(dev, :max_work_item_size)

global_memory(dev::cl.Device) = cl.info(dev, :global_mem_size) |> Int
local_memory(dev::cl.Device) = cl.info(dev, :local_mem_size) |> Int

function supports_double(dev::cl.Device)
    result = Ref{cl.CL_uint}()
    cl.@check cl.api.clGetDeviceInfo(
        dev.id, cl.CL_DEVICE_PREFERRED_VECTOR_WIDTH_DOUBLE,
        sizeof(cl.CL_uint), result,
        C_NULL
    )
    result[] != 0
end


text:
module PlotAdaptiveGrid1D
import Random: MersenneTwister
export adapted_grid
include("adapted_grid.jl")
end


text:
using Dierckx

"""
    a₀ : Bohr radius (in cm)
"""
const a₀ = ustrip(BohrRadius |> u"cm") # 0.529 Å

"""
    Rₐ(elm::Element)

Classic formula for the atomic screening radius in cm
"""
Rₐ(elm::Element) = a₀ * z(elm)^-0.333333333333

"""
Algorithms implementing the elastic scattering cross-section

    σₜ(::Type{ScreenedRutherford}, elm::Element, E::Float64)
"""
abstract type ElasticScatteringCrossSection end

"""
Basic screened Rutherford algorithm where V(r) = (-Ze²/r)exp(-r/R) where R=a₀Z⁻¹/³ is solved using
the first Born approximation.
"""
abstract type ScreenedRutherfordType <: ElasticScatteringCrossSection end

struct ScreenedRutherford <: ScreenedRutherfordType end

"""
Liljequist's simple refinement of the basic ScreenedRutherford algorithm.

Journal of Applied Physics, 65, 24-31 (1989) as corrected in J. Appl. Phys. 68 (7) 3061-3065 (1990)
"""
struct Liljequist1989 <: ScreenedRutherfordType end

"""
Browning's scattering cross section according to a draft 1994 article
"""
struct Browning1994 <: ScreenedRutherfordType end

"""
Browning's scattering cross section

Appl. Phys. Lett. 58, 2845 (1991); https://doi.org/10.1063/1.104754
"""
struct Browning1991 <: ScreenedRutherfordType end

"""
    ξ(::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64)::Float64

  * E in eV
"""
function ξ(::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64)::Float64
    Rₑ, mc² = ustrip((PlanckConstant * SpeedOfLightInVacuum * RydbergConstant) |> u"eV"), ustrip(ElectronMass * SpeedOfLightInVacuum^2 |> u"eV")
    return 0.5 * π * a₀^2 * (4.0 * z(elm) * ((E + mc²) / (E + 2.0 * mc²)) * (Rₑ / E))^2 # As corrected in Liljequist1989

end

"""
    ϵ(elm::Element, E::Float64)

Screening factor.
"""
ϵ(::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64) = 2.0 * (kₑ(E) * Rₐ(elm))^2

# A spline interpolation based on Liljequist Table III
const LiljequistCorrection = begin
    zs = Float64[4, 6, 13, 20, 26, 29, 32, 38, 47, 50, 56, 64, 74, 79, 82]
    es =
        1000.0 * [
            0.1,
            0.15,
            0.2,
            0.3,
            0.4,
            0.5,
            0.7,
            1,
            1.5,
            2,
            3,
            4,
            5,
            7,
            10,
            15,
            20,
            30,
            40,
            50,
            70,
            100,
        ] # eV
    tbl3 = [
        1.257 1.211 1.188 1.165 1.153 1.147 1.139 1.134 1.129 1.127 1.123 1.121 1.119 1.116 1.113 1.111 1.11 1.11 1.11 1.112 1.115 1.119
        1.506 1.394 1.330 1.257 1.217 1.192 1.163 1.140 1.123 1.115 1.107 1.103 1.100 1.097 1.095 1.093 1.093 1.093 1.094 1.096 1.099 1.104
        3.14 2.589 2.301 1.993 1.824 1.714 1.576 1.458 1.352 1.294 1.23 1.196 1.175 1.15 1.132 1.118 1.111 1.105 1.103 1.103 1.104 1.107
        3.905 3.192 2.823 2.429 2.211 2.065 1.878 1.713 1.558 1.466 1.358 1.295 1.254 1.202 1.161 1.128 1.111 1.094 1.087 1.084 1.081 1.081
        7.061 5.076 4.102 3.207 2.786 2.536 2.24 1.997 1.781 1.655 1.508 1.421 1.363 1.288 1.225 1.17 1.141 1.111 1.097 1.088 1.08 1.076
        10.154 7.184 5.521 3.987 3.31 2.933 2.517 2.199 1.929 1.777 1.603 1.501 1.433 1.344 1.269 1.202 1.166 1.129 1.110 1.098 1.087 1.079
        13.213 8.783 6.638 4.674 3.796 3.309 2.781 2.39 2.068 1.891 1.691 1.575 1.497 1.397 1.31 1.232 1.189 1.144 1.12 1.106 1.09 1.08
        16.187 9.955 7.475 5.351 4.378 3.813 3.176 2.694 2.297 2.081 1.838 1.698 1.605 1.484 1.379 1.282 1.227 1.167 1.134 1.114 1.091 1.075
        16.265 13.22 9.985 6.781 5.391 4.624 3.785 3.156 2.643 2.365 2.057 1.882 1.765 1.615 1.483 1.36 1.289 1.207 1.162 1.133 1.098 1.071
        14.33 13.534 10.829 7.385 5.802 4.936 4.005 3.319 2.763 2.464 2.132 1.944 1.82 1.659 1.518 1.387 1.31 1.222 1.172 1.133 1.101 1.071
        15.97 13.946 11.072 7.779 6.189 5.285 4.293 3.554 2.951 2.624 2.259 2.052 1.913 1.735 1.579 1.433 1.347 1.247 1.189 1.14 1.105 1.068
        25.349 25.137 19.078 11.778 8.611 6.957 5.298 4.19 3.367 2.946 2.493 2.242 2.077 1.866 1.682 1.512 1.411 1.293 1.224 1.151 1.12 1.073
        43.544 47.121 37.079 21.021 13.996 10.545 7.322 5.369 4.072 3.464 2.849 2.524 2.315 2.053 1.828 1.622 1.501 1.358 1.274 1.178 1.145 1.085
        52.567 58.958 46.173 24.508 15.831 11.792 8.101 5.87 4.388 3.7 3.015 2.655 2.426 2.14 1.896 1.674 1.543 1.389 1.298 1.218 1.157 1.09
        43.737 52.101 45.521 25.950 16.794 12.472 8.537 6.160 4.576 3.843 3.115 2.735 2.493 2.193 1.938 1.705 1.568 1.407 1.312 1.247 1.164 1.093
    ]
    # Spline it first at each energy over Z
    zspl = [Spline1D(zs, tbl3[:, ie], k = 3) for ie in eachindex(es)]
    # Then use this to estimate the interpolation over the range Z=1 to 92
    [Spline1D(es, [zspl[ie](Float64(z)) for ie in eachindex(es)], k = 3) for z = 1:92]
end

"""
    σₜᵣ(::Type{ScreenedRutherford}, elm::Element, E::Float64)

The transport cross-section in cm².  The transport cross-section gives the correct transport
mean free path - the mean free path in the direction of initial propagation after an infinite
number of collisions.
"""
function σₜᵣ(::Type{ScreenedRutherford}, elm::Element, E::Float64)
    ϵv = ϵ(ScreenedRutherford, elm, E)
    return ξ(ScreenedRutherford, elm, E) * (log(2.0 * ϵv + 1) - 2.0 * ϵv / (2.0 * ϵv + 1.0))
end
function σₜᵣ(::Type{Liljequist1989}, elm::Element, E::Float64)
    return σₜᵣ(ScreenedRutherford, elm, E) / LiljequistCorrection[z(elm)](E)
end

"""
    σₜ(::Type{ScreenedRutherford}, elm::Element, E::Float64)
    σₜ(::Type{Liljequest1989}, elm::Element, E::Float64)
    σₜ(::Type{Browning1991}, elm::Element, E::Float64)
    σₜ(::Type{Browning1994}, elm::Element, E::Float64)

Total cross section per atom in cm².
"""
function σₜ(::Type{ScreenedRutherford}, elm::Element, E::Float64)
    ϵv = ϵ(ScreenedRutherford, elm, E)
    return ξ(ScreenedRutherford, elm, E) * (2.0 * ϵv^2 / (2.0 * ϵv + 1.0))
end
function σₜ(::Type{Liljequist1989}, elm::Element, E::Float64)
    return σₜ(ScreenedRutherford, elm, E) / LiljequistCorrection[z(elm)](E)
end
function σₜ(::Type{Browning1991}, elm::Element, E::Float64)
    e = 0.001 * E
    u = log10(8.0 * e * z(elm)^-1.33)
    return 4.7e-18 * (z(elm)^1.33 + 0.032 * z(elm)^2) / (
        (e + 0.0155 * (z(elm)^1.33) * sqrt(e)) * (1.0 - 0.02 * sqrt(z(elm)) * exp(-u^2))
    )
end
function σₜ(::Type{Browning1994}, elm::Element, E::Float64)
    e = 0.001 * E
    return 3.0e-18 * z(elm)^1.7 /
           (e + 0.005 * z(elm)^1.7 * sqrt(e) + 0.0007 * z(elm)^2 / sqrt(e))
end

"""
    δσδΩ(::Type{ScreenedRutherford}, θ::Float64, elm::Element, E::Float64)::Float64

The *differential* screened Rutherford cross-section per atom.
"""
function δσδΩ(::Type{ScreenedRutherford}, θ::Float64, elm::Element, E::Float64)::Float64
    return ξ(ScreenedRutherford, elm, E) *
           (1.0 - cos(θ) + ϵ(ScreenedRutherford, elm, E)^-1)^-2
end
function δσδΩ(::Type{Liljequist1989}, θ::Float64, elm::Element, E::Float64)::Float64
    return σ(ScreenedRutherford, θ, elm, E) / LiljequistCorrection[z(elm)](E)
end

"""
    λ(ty::Type{<:ElasticScatteringCrossSection}, θ::Float64, elm::Element, E::Float64)::Float64

The mean free path.  The mean distance between elastic scattering events.
"""
function λ(ty::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64, N::Float64)
    return (σₜ(ty, elm, E) * N)^-1
end
function λ(ty::Type{<:ScreenedRutherfordType}, mat::Material, elm::Element, E::Float64)
    return λ(ty, elm, E, atoms_per_cm³(mat, elm))
end

"""
    Base.rand(ty::Type{<:ScreenedRutherfordType}, mat::Material, E::Float64)::NTuple{3, Float64}

 Returns a randomly selected elastic scattering event description.  The result is ( λ, θ, ϕ ) where
 λ is a randomized mean free path for the first scattering event.  θ is a randomized scattering
 angle on (0.0, π) and ϕ is a randomized azimuthal angle on [0, 2π).

 The algorithm considers scattering by any element in the material and picks the shortest randomized
 path.  This implementation depends on two facts: 1) We are looking for the first scattering event
 so we consider all the elements and pick the one with the shortest path. 2) The process is memoryless.
"""
function Base.rand(
    ty::Type{<:ScreenedRutherfordType},
    mat::Material,
    E::Float64,
)::NTuple{3,Float64}
    elm′, λ′ = elements[119], 1.0e308
    for (i, z) in enumerate(keys(mat))
        l = -λ(ty, mat, z, E) * log(rand())
        (elm′, λ′) = l < λ′ ? (z, l) : (elm′, λ′)
    end
    @assert elm′ != elements[119] "Are there any elements in $mat?  Is the density ($(mat[:Density])) too low?"
    return (λ′, rand(ty, elm′, E), 2.0 * π * rand())
end

"""
    λₜᵣ(ty::Type{<:ElasticScatteringCrossSection}, θ::Float64, elm::Element, E::Float64)::Float64
    λₜᵣ(ty::Type{<:ScreenedRutherfordType}, mat::Material, elm::Element, E::Float64)

The transport mean free path. The mean distance in the initial direction of propagation between
elastic scattering events.

  * N is the number of atoms per cm³
  * E is the electron kinetic energy in eV
"""
function λₜᵣ(
    ty::Type{<:ElasticScatteringCrossSection},
    elm::Element,
    E::Float64,
    N::Float64,
)
    return (σₜᵣ(ty, elm, E) * N)^-1
end
function λₜᵣ(ty::Type{<:ScreenedRutherfordType}, mat::Material, elm::Element, E::Float64)
    return λₜᵣ(ty, elm, E, atoms_per_cm³(mat, elm))
end

"""
    Base.rand(::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64)::Float64

Draw an angle distributed according to the angular dependence of the differential screened Rutherford cross-section.
"""
function Base.rand(ty::Type{<:ScreenedRutherfordType}, elm::Element, E::Float64)::Float64
    Y = rand()
    return acos(1.0 + (Y - 1.0) / (ϵ(ty, elm, E) * Y + 0.5))
end
function Base.rand(ty::Type{Browning1994}, elm::Element, E::Float64)::Float64
    α, R = 7.0e-3 / (0.001 * E), rand()
    return acos(1.0 - 2.0 * α * R / (1.0 + α - R))
end