text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     comment_magics: false
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Julia 1.5.3
#     language: julia
#     name: julia-1.5
# ---

# # Sampling from a probability distribution
#
# This is a test notebook to study the best way to sample a density function.
# It is not directly linked with the GEMPIC package.

import Pkg; Pkg.add(["Plots", "BenchmarkTools", "Sobol"])

using Plots, BenchmarkTools, Sobol

s = SobolSeq(2)

n = 1000
points = zeros(Float64,(2,n))
for i = 1:n
    points[:,i] .= next!(s)
end
scatter(points[1,:], points[2,:])

"""
Solve ``P(x) = \\frac{2\\pi R}{k} ``

```math
P(x) = \\int_0^x (1 + \\alpha cos(k y)) dy
```
```math
P(x) = x + \\frac{\\alpha}{k} sin (k x)
```

"""
function dichotomie(R)

    kx, alpha = 0.5, 0.1
    eps = 1e-12
    a, b = 0, 2π/kx
    while true
        x = 0.5 * (a+b)
        D = ( kx * x + alpha * sin(kx*x)) / 2π
        if ( D<R-eps )
           a = x
        elseif ( D>R+eps )
           b = x
        else
           return x
        end
    end
end

"""
Input r is a random number ``\\in [0,1]``

```math
    f(x) = 1 + \\alpha cos(k x)
```
on some domain ``[0, 2\\pi/k]``

Solve the equation ``P(x)-r=0`` with Newton’s method

```math
    x^{n+1} = x^n – (P(x)-(2\\pi r / k)/f(x)
```

with
```math
P(x) = \\int_0^x (1 + \\alpha cos(k y)) dy
```
```math
P(x) = x + \\frac{\\alpha}{k} sin (k x)
```
"""
function newton(r)
    x0, x1 = 0.0, 1.0
    alpha, k = 0.1, 0.5
    r *= 2 * pi / k
    while (abs(x1-x0) > 1e-12)
        p = x0 + alpha * sin( k * x0) / k
        f = 1 + alpha * cos( k * x0)
        x0, x1 = x1, x0 - (p - r) / f
    end
    x1
end

?newton

@btime newton(0.9)

?dichotomie

@btime dichotomie(0.9)

# +
function landau( nbpart :: Int64)

   xp = Float64[]
   vp = Float64[]

   s = SobolSeq(2)

   for k=0:nbpart-1

      v = sqrt(-2 * log( (k+0.5)/nbpart))
      r1, r2 = next!(s)
      θ = r1 * 2π
      push!(xp,  newton(r2))
      push!(vp,  v * sin(θ))

   end

   xp, vp

end

xp, vp = landau(100000);
# -

p = histogram([xp,vp], normalize=true, bins = 100,  layout=(2,1))
plot!(p[1,1], x-> (1+0.1*cos(0.5*x))/4π, 0., 4π)
plot!(p[2,1], x-> (exp(-x^2/2))/sqrt(2π), -6, 6)

ρ(x) = exp(-x^2/2)
xmin, xmax, nx = -5, 5, 100
x = range( xmin, stop=xmax, length=nx)
f = map(ρ, x)
plot(ρ, xmin, xmax)
scatter!(x, f)

dx = (xmax - xmin)/(nx-1)
@show nx * dx == xmax - xmin
v = (cumsum(f) .* dx) / sqrt(2π)
plot( x, v)

samples  = rand(10000)
for i in eachindex(samples)
    j = findmin(abs.(v .- samples[i] ))[2]
    samples[i] =  xmin + j * dx
end
histogram(samples, normalize=true, bins = 100)
plot!(x, f ./ sqrt(2π))

s = SobolSeq(1)
for i in eachindex(samples)
    rnd = next!(s)[1]
    j = findmin(abs.(v .- rnd ))[2]
    samples[i] =  xmin + j * dx
end
histogram(samples, normalize=true, bins = 100)
plot!(x, f ./ sqrt(2π))

alpha, kx = 0.5, 0.5
xmin, xmax, nx = 0.0, 2π/kx, 10000
dx = (xmax - xmin) / (nx-1)
x = range(xmin, stop=xmax, length=nx) |> collect
f = 1 .+ alpha*cos.(kx .* x)
v = cumsum(f)*dx
plot(x, v)

# Construct the CDF numerically and find the closest value

nbpart = 10000
s  = SobolSeq(1)
xp = Float64[]
for k=0:nbpart-1
   r = next!(s)[1] * 4π
   j = findmin(abs.(v .- r) )[2]
   push!(xp,  x[j])
end

histogram(xp, normalize=true, bins = 100)
plot!(x-> (1+alpha*cos(kx*x))/4π, 0., 4π)


text:
module TestCurrentController

println("\nCurrentController: Demonstrating the ability to simulate mixed domain models")

using Modia
using Modia.Electric
using Modia.Rotational
using Modia.Blocks

# Desired:
#   using ModiaMath: plot
#
# In order that ModiaMath need not to be defined in the user environment, it is included via Modia:
using Modia.ModiaMath: plot

@model CurrentController begin
  k=30   # Gain of PI current controller
  T=0.005   # Time constant of PI current controller
  Jmotor=Inertia(J=0.0025)
  emf=EMF(k=1.016) #, phi=Float(state=false))
  inductor=Inductor(L=0.061)
  resistor=Resistor(R=13.8)
  ground=Ground()
  signalVoltage=SignalVoltage()
  currentSensor=CurrentSensor()
  firstOrder=FirstOrder(k=1,T=0.001)
  feedback=Feedback()
  PI1=PI(T=T, k=k)
  gear=IdealGear(ratio=105)
  load=Inertia(J=100)
  step=Step(height=1,startTime=0)
  spring=SpringDamper(c=5e5, d=500)
@equations begin
  connect(emf.flange, Jmotor.flange_a)
  connect(resistor.n, inductor.p)
  connect(ground.p, emf.n)
  connect(inductor.n, emf.p)
  connect(signalVoltage.p, resistor.p)
  connect(currentSensor.p, ground.p)
  connect(currentSensor.n, signalVoltage.n)
  connect(currentSensor.i, feedback.u2)
  connect(PI1.y, firstOrder.u)
  connect(feedback.y, PI1.u)
  connect(signalVoltage.v, firstOrder.y)
  connect(feedback.u1, step.y)
  connect(Jmotor.flange_b, gear.flange_a)
  connect(load.flange_a, spring.flange_b)
  connect(gear.flange_b, spring.flange_a)
  end
end

#result = simulate(CurrentController, 0.1, storeEliminated=false)

# result = checkSimulation(CurrentController, 0.1, "load.w", 0.07929151315487117, removeSingularities=false)
# plot(result, ("currentSensor.i", "step.y", "load.w"))

@static if VERSION < v"0.7.0-DEV.2005"
  result = checkSimulation(CurrentController, 0.1, "load.w", 0.07927285979038304, removeSingularities=true)
  plot(result, [("currentSensor.i", "step.y"), "load.w"], heading="CurrentController", figure=11)
end

end

ext: jl


text:
using Distributed
using DifferentialEquations
using ParameterizedFunctions
using Plots
using Random
using DataFrames
using CSV

pyplot()

################################################################################
# Model:
# A <-> R <-> I
# ka: R -> A
# kr: A -> R
# ki: R -> I
# kp: I -> R
################################################################################

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]
end

function isoutofdomain(c, k, t)
    if (c[1]<=1) && (c[2]<=1) && (c[3]<=1) && (c[1]>=0) && (c[2]>=0) && (c[3]>=0)
        return false
    else
        return true
    end
end

function run_recruitment(k_time_pairs)
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    c0 = [1.0; 0.0; 0.0]

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, isoutofdomain=isoutofdomain, saveat=0.01)
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    data = [as rs is]

    return [t, data]
end

function get_sim_data(ka, ks, ki)
    kas = Float64[]
    kss = Float64[]
    kis = Float64[]
    kps = Float64[]
    ts = Float64[]
    as = Float64[]
    rs = Float64[]
    is = Float64[]
    k_recruit = Dict(
        'a' => ka,
        's' => ks,
        'i' => ki,
        'p' => 0)
    k_release = Dict(
        'a' => ka,
        's' => 0,
        'i' => 0,
        'p' => 0)
    recruit_pair = [k_recruit, (0.0, 5.0)]
    release_pair = [k_release, (5.0, 35.0)]
    sol = run_recruitment([recruit_pair, release_pair])
    ts = [ts ; sol[1]]
    as = [as ; sol[2][:, 1]]
    rs = [rs ; sol[2][:, 2]]
    is = [is ; sol[2][:, 3]]

    append!(kas, transpose([ka for i in range(1, stop=length(sol[1]))]))
    append!(kss, transpose([ks for i in range(1, stop=length(sol[1]))]))
    append!(kis, transpose([ki for i in range(1, stop=length(sol[1]))]))
    append!(kps, transpose([0  for i in range(1, stop=length(sol[1]))]))
    return [kas kss kis kps ts as rs is]
end

function run_sims_orig()
    k_scale_ki = collect(0:0.025:0.5)
    k_scale_ka = collect(0:0.05:1)
    k_scale_ks = collect(0:0.5:10)
    for ka in k_scale_ka
        for ks in k_scale_ks
            for ki in k_scale_ki
                data = get_sim_data(ka, ks, ki)
                df_dict = Dict(
                    "ka" => data[:, 1],
                    "ks" => data[:, 2],
                    "ki" => data[:, 3],
                    "kp" => data[:, 4],
                    "t" => data[:, 5],
                    "a" => data[:, 6],
                    "r" => data[:, 7],
                    "i" => data[:, 8])
                df = DataFrame(df_dict)
                println("./sims_orig_3sm_5d/ka-" * string(ka) * "_ks-" * string(ks) * "_ki-" * string(ki) *".csv")
                CSV.write("./sims_orig_3sm_5d/ka-" * string(ka) * "_ks-" * string(ks) * "_ki-" * string(ki) *".csv", df)
            end
        end
    end
end

function run_sims_rand(num_samples = 8000)
    function get_df(kvals)
        ka, ks, ki = kvals .* [1, 10, 0.5]
        data = get_sim_data(ka, ks, ki)
        df_dict = Dict(
                       "ka" => data[:, 1],
                       "ks" => data[:, 2],
                       "ki" => data[:, 3],
                       "kp" => data[:, 4],
                       "t" => data[:, 5],
                       "a" => data[:, 6],
                       "r" => data[:, 7],
                       "i" => data[:, 8])
        df = DataFrame(df_dict)
        return df
    end

    kvals = [rand(3) for i in 1:num_samples]
    df = vcat(pmap(get_df, kvals)...)
    CSV.write("./randomized_cr_df.csv", df)
    println("Wrote dataframe for " * string(num_samples) * " samples.")
end


text:
function print_description(res_dir)
    tree = Core.ResDirTree(res_dir);
    println(Core.load_desc(tree));
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
# ---
# title: 906. Super Palindromes
# id: problem906
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/super-palindromes/description/>
# hidden: true
# ---
#
# Let's say a positive integer is a  _superpalindrome_  if it is a palindrome,
# and it is also the square of a palindrome.
#
# Now, given two positive integers `L` and `R` (represented as strings), return
# the number of superpalindromes in the inclusive range `[L, R]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: L = "4", R = "1000"
#     Output: 4
#     **Explanation** : 4, 9, 121, and 484 are superpalindromes.
#     Note that 676 is not a superpalindrome: 26 * 26 = 676, but 26 is not a palindrome.
#
#
#
# **Note:**
#
#   1. `1 <= len(L) <= 18`
#   2. `1 <= len(R) <= 18`
#   3. `L` and `R` are strings representing integers in the range `[1, 10^18)`.
#   4. `int(L) <= int(R)`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module BerghAlgorithms
include("StackyFan.jl")
include("AlgA.jl")
include("AlgC.jl")
# include("AlgDp.jl")

function __init__()

    println("         ")
    println("   ^     ")
    println("  / \\    | WARNING! The BerghAlgorithms module is new and key features are still being implemented and tested.")
    println(" / ! \\   | As such, please be patient as bugs and mathematical inaccuracies are resolved.")
    println("-------  | If you find an issue, please consider reporting it to htalbott@umich.edu.")

end

end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
module JaynesCummings
using Compat

if VERSION < v"0.5.0"
    product(x,y) = reshape(collect(Compat.Iterators.product(x,y)),length(x),length(y))
else
    import Base.product
end

export
gen_initialstate, gen_hamiltonian, gen_displacementop, gen_timeevoarray,
gen_drive_hamiltonian,
calc_qubittimeevo, calc_qubittimeevo_rel, calc_photonnumbers,
calc_densitymatrix_resonator, calc_wignerfunction_resonator,
plot_qubittimeevo, plot_photonnumbers, plot_densitymatrix, plot_wignerfunction,
plot_wignerfunction3d,
ptrace,
ħ, a, a_dagger, σ_minus, σ_plus, σ_z, qeye, ⊗

include("constants.jl")
include("generation.jl")
include("calculation.jl")

load_plotfuns() = include(joinpath(dirname(@__FILE__),"plotting.jl"))
runtests() = include(joinpath(dirname(dirname(@__FILE__)),"test","runtests.jl"))
runperf() = include(joinpath(dirname(dirname(@__FILE__)),"test","perf","runperf.jl"))

end


text:
module CUDATutorials

using CUDA, Literate

function make_dir(path)
    try
        mkdir(path)
    catch e
        e.code == -17 && return true
        throw(e)
    end
end

"""
    generate_notebooks(outpath = pwd())
Output all tutorials as Jupyter Notebooks to `outpath`
"""
function generate_notebooks(outpath = pwd())
    outpath = joinpath(outpath, "notebook")
    make_dir(outpath)
    tutorial_path = joinpath(@__DIR__, "tutorials")
    for (root, dirs, files) in walkdir(tutorial_path)
        for dir in dirs
            make_dir(joinpath(outpath, dir))
        end

        for file in files
            # TODO: non-regex solution ?
            m = match(Regex("^($tutorial_path)/(.+)"), root)

            if file[(end - 2):end] == ".jl"
                Literate.notebook(joinpath(root, file), joinpath(outpath, String(m[2]));
                                    name = file, execute=false, credit = false)
            else
                cp(joinpath(root, file), joinpath(outpath, String(m[2]), file))
            end
        end
    end
end

end #module


text:
# ============================================================
# Hydrodynamic Fluxes
# Reference: http://www.cfdbooks.com/
# ============================================================

"""
    flux_upwind(uL, uR, Ω::T, n::T, dt) where {T<:AbstractVector}

Upwind flux

"""
function flux_upwind(uL, uR, Ω::T, n::T, dt = 1.0) where {T<:AbstractVector}
    ip = dot(Ω, n) # inner product

    if ip > 0
        return dt * ip * uL
    else
        return dt * ip * uR
    end
end

"""
    flux_lax!(fw::AbstractArray{<:Real,1}, wL::AbstractArray{<:Real,1}, wR::AbstractArray{<:Real,1}, γ::Real, dt::Real, dx::Real)

Lax-Friedrichs flux

_P. D. Lax, Weak Solutions of Nonlinear Hyperbolic Equations and Their Numerical Computation,
Commun. Pure and Applied Mathematics, 7, 159-193, 1954._

"""
function flux_lax!(
    fw::X,
    wL::Y,
    wR::Y,
    γ,
    dt,
    dx,
) where {X<:AbstractArray{<:Real,1},Y<:AbstractArray{<:Real,1}}
    fw .= 0.5 * dt .* (euler_flux(wL, γ)[1] + euler_flux(wR, γ)[1] - dx / dt .* (wR - wL))
    return nothing
end

"""
    flux_hll!(fw::AbstractArray{<:Real,1}, wL::AbstractArray{<:Real,1}, wR::AbstractArray{<:Real,1}, γ::Real, dt::Real)

HLL flux for the Euler equations

- @args: variables at left & right sides of interface
- @args: specific heat ratio

"""
function flux_hll!(
    fw::X,
    wL::Y,
    wR::Y,
    γ,
    dt,
) where {X<:AbstractArray{<:Real,1},Y<:AbstractArray{<:Real,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    aL = sound_speed(primL, γ)
    aR = sound_speed(primR, γ)

    λmin = primL[2] - aL
    λmax = primR[2] + aR

    if λmin >= 0.0
        fw .= euler_flux(wL, γ)[1]
    elseif λmax <= 0.0
        fw .= euler_flux(wR, γ)[1]
    else
        factor = 1.0 / (λmax - λmin)

        flux1 = euler_flux(wL, γ)[1]
        flux2 = euler_flux(wR, γ)[1]

        @. fw = factor * (λmax * flux1 - λmin * flux2 + λmax * λmin * (wR - wL))
    end

    fw .*= dt

    return nothing

end

"""
    flux_roe!(fw::AbstractArray{<:Real,1}, wL::AbstractArray{<:Real,1}, wR::AbstractArray{<:Real,1},
    γ::Real, dt::Real, n = [1.0, 0.0]::AbstractArray{<:Real,1})

Roe's flux with entropy fix

_P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference Schemes, Journal of Computational Physics, 43, pp. 357-372._
(_cf. http://cfdbooks.com/cfdcodes.html_)

* @args primL[1:4] = left state (rhoL, uL, vL, pL)
* @args primR[1:4] = right state (rhoR, uR, vR, pR)
* @args γ: specific heat ratio
* @args n[2]: unit face normal (L -> R)

"""
function flux_roe!(
    fw::X,
    wL::Y,
    wR::Y,
    γ,
    dt,
    n = [1.0, 0.0],
) where {X<:AbstractArray{<:Real,1},Y<:AbstractArray{<:Real,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    if length(fw) == 3 # 1D

        # left state
        rhoL = primL[1]
        vL = primL[2]
        pL = 0.5 * primL[1] / primL[end]
        aL = sound_speed(primL, γ)
        HL = (wL[end] + pL) / rhoL

        # right state
        rhoR = primR[1]
        vR = primR[2]
        pR = 0.5 * primR[1] / primR[end]
        aR = sound_speed(primR, γ)
        HR = (wR[end] + pR) / rhoR

        # compute Roe averages
        RT = sqrt(rhoR / rhoL)
        rho = RT * rhoL
        v = (vL + RT * vR) / (1.0 + RT)
        H = (HL + RT * HR) / (1.0 + RT)
        a = sqrt((γ - 1.0) * (H - 0.5 * v^2))

        # differences in primitive variables
        drho = rhoR - rhoL
        du = vR - vL
        dP = pR - pL

        # wave strength (characteristic variables)
        dV = [
            0.5 * (dP - rho * a * du) / a^2,
            -(dP / a^2 - drho),
            0.5 * (dP + rho * a * du) / a^2,
        ]

        # absolute values of the wave speeds (eigenvalues)
        ws = [abs(v - a), abs(v), abs(v + a)]

        # modified wave speeds for nonlinear fields (to remove expansion shocks).
        Da = max(0.0, 4.0 * ((vR - aR) - (vL - aL)))
        if ws[1] < 0.5 * Da
            ws[1] = ws[1]^2 / Da + 0.25 * Da
        end
        Da = max(0.0, 4.0 * ((vR + aR) - (vL + aL)))
        if ws[3] < 0.5 * Da
            ws[3] = ws[3]^2 / Da + 0.25 * Da
        end

        # right eigenvectors
        R = zeros(3, 3)
        R[1, 1] = 1.0
        R[2, 1] = v - a
        R[3, 1] = H - v * a

        R[1, 2] = 1.0
        R[2, 2] = v
        R[3, 2] = 0.5 * v * v

        R[1, 3] = 1.0
        R[2, 3] = v + a
        R[3, 3] = H + v * a

        # compute average flux
        fw .= 0.5 .* (euler_flux(wL, γ)[1] + euler_flux(wR, γ)[1])

        # add matrix dissipation term to complete Roe flux
        for j = 1:3, k = 1:3
            fw[j] -= 0.5 * ws[k] * dV[k] * R[j, k]
        end

    elseif length(fw) == 4 # 2D

        # normal vector
        nx = n[1]
        ny = n[2]

        # tangent vector
        mx = -ny
        my = nx

        #--- primitive and other variables ---#
        # left state
        rhoL = primL[1]
        uL = primL[2]
        vL = primL[3]
        unL = uL * nx + vL * ny
        umL = uL * mx + vL * my
        pL = 0.5 * primL[1] / primL[4]
        aL = sound_speed(primL[4], γ)
        HL = aL^2 / (γ - 1.0) + 0.5 * (uL^2 + vL^2)

        # right state
        rhoR = primR[1]
        uR = primR[2]
        vR = primR[3]
        unR = uR * nx + vR * ny
        umR = uR * mx + vR * my
        pR = 0.5 * primR[1] / primR[4]
        aR = sound_speed(primR[4], γ)
        HR = aR^2 / (γ - 1.0) + 0.5 * (uR^2 + vR^2)

        # Roe averages
        RT = sqrt(rhoR / rhoL)
        rho = RT * rhoL
        u = (uL + RT * uR) / (1.0 + RT)
        v = (vL + RT * vR) / (1.0 + RT)
        H = (HL + RT * HR) / (1.0 + RT)
        a = sqrt((γ - 1.0) * (H - 0.5 * (u^2 + v^2)))
        un = u * nx + v * ny
        um = u * mx + v * my

        # wave strengths
        drho = rhoR - rhoL
        dp = pR - pL
        dun = unR - unL
        dum = umR - umL

        LdU = [
            (dp - rho * a * dun) / (2.0 * a^2),
            rho * dum,
            drho - dp / (a^2),
            (dp + rho * a * dun) / (2.0 * a^2),
        ]

        # wave speed
        ws = abs.([un - a, un, un, un + a])

        # Harten's entropy fix JCP(1983), 49, pp357-393
        # only for the nonlinear fields.
        if ws[1] < 0.2
            ws[1] = 0.5 * (ws[1]^2 / 0.2 + 0.2)
        end
        if ws[4] < 0.2
            ws[4] = 0.5 * (ws[4]^2 / 0.2 + 0.2)
        end

        # right eigenvectors
        Rv = zeros(4, 4)

        Rv[1, 1] = 1.0
        Rv[2, 1] = u - a * nx
        Rv[3, 1] = v - a * ny
        Rv[4, 1] = H - un * a

        Rv[1, 2] = 0.0
        Rv[2, 2] = mx
        Rv[3, 2] = my
        Rv[4, 2] = um

        Rv[1, 3] = 1.0
        Rv[2, 3] = u
        Rv[3, 3] = v
        Rv[4, 3] = 0.5 * (u * u + v * v)

        Rv[1, 4] = 1.0
        Rv[2, 4] = u + a * nx
        Rv[3, 4] = v + a * ny
        Rv[4, 4] = H + un * a

        # dissipation term
        diss = zeros(4)
        for i = 1:4, j = 1:4
            diss[i] += ws[j] * LdU[j] * Rv[i, j]
        end

        # compute fluxes
        fL = [
            rhoL * unL,
            rhoL * unL * uL + pL * nx,
            rhoL * unL * vL + pL * ny,
            rhoL * unL * HL,
        ]
        fR = [
            rhoR * unR,
            rhoR * unR * uR + pR * nx,
            rhoR * unR * vR + pR * ny,
            rhoR * unR * HR,
        ]

        @. fw = 0.5 * dt * (fL + fR - diss)

    end

    return nothing

end


text:
##############
using Revise
using FDMIC_Geodynamics
using LoopVectorization, Printf, Base.Threads, Plots, Revise, LinearAlgebra, Statistics, SparseArrays, LazyGrids
##############
function SetMarkers!( p::Markers, params::ModelParameters, dom::ModelDomain, α )
    R     = params.user[1]
    a_ell = 1.75
    b_ell = 0.25
    # Use this function to set up the model geometry
    for k=1:p.nmark
        x_ell =  (p.x[k])*cos(α) + (p.y[k])*sin(α)
        y_ell = -(p.x[k])*sin(α) + (p.y[k])*cos(α)
        in         = ((x_ell/a_ell)^2 + (y_ell/b_ell)^2) < R^2
        p.phase[k] = (in==0)* 1.0 + (in==1)*2.0
    end
end
##############
function Rheology!( f, Eps, Tau, params, materials, BC, ncx, ncy )
    # Compute effective viscosity for each phase
    f.etac .= 0.0
    for m=1:materials.nphase
        eta      = materials.eta0[m] * Eps.II.^(1.0/materials.n[m] - 1.0)
        f.etac .+= f.phase_perc[m,:,:] .* eta
    end
    @time CentroidsToVertices!( f.etav, f.etac, ncx, ncy, BC )
    println("min Eii: ", minimum(Eps.II), " --- max Eii: ", maximum(Eps.II))
    println("min eta: ", minimum(f.etac), " --- max eta: ", maximum(f.etac))
end
export Rheology!
##############
@views function main( N, α, eta_fact )
    # Main routine
    println("#####################")
    println("###### M2Di.jl ######")
    println("#####################")
    # Name of the experiement
    experiment       = "PeriodicSimpleShear"
    # Domain
    domain           = ModelDomain()
    domain.xmin      = -1.0
    domain.xmax      =  1.0
    domain.ymin      = -1.0
    domain.ymax      =  1.0
    # parameters
    params           = ModelParameters()
    params.Ebg       = -1.0      # reference strain rate (>0 horizontal compression)
    # Scales
    scale            = CharScales()
    scale.L          = 1.0
    scale.T          = 1.0
    scale.t          = 1.0/abs(params.Ebg)
    # Spatial discretisation
    ncx              = N; domain.ncx = ncx
    ncy              = N; domain.ncy = ncy
    # Time discretisation
    params.advection = 0
    params.nt        = 1
    params.dt        = 1e4
    params.Courant   = 0.40  # Courant number
    params.dt_var    = 1     # variable dt
    params.t         = 0.0   # time
    # Boundary conditions
    params.BC_type   = 1     # 1: Pure shear / 2: Simple shear / 3: Simple shear periodic
    params.ALE       = 0     # Deform box
    params.comp      = 0     # Compressible
    # Solver
    #------------------------ PINNED PRESSURE
    params.solver    = 2     # -1: coupled with pinned P, 0: coupled with slight compressibility --- 1: decoupled Powell-Hestenes --- 2: decoupled KSP GCR
    #------------------------ PINNED PRESSURE
    params.gamma     = 1e4   # Penalty factor
    params.Dir_scale = 1.0   # Dirichlet scaling factor
    # Non-linear iterations
    params.niter_nl  = 5     # max. number of non-linear iterations
    params.tol_nl    = 1e-7  # non-linear tolerance
    params.JFNK      = 0     # Jacobian-Free Newton_Krylov
    # Visualisation
    params.show_figs = 1     # activates visualisation...
    params.nout      = 10    # ... every nout
    # Initialize coordinates
    domain.L0 = domain.xmax-domain.xmin
    GenerateMesh!( domain )
    dx, dy = domain.dx, domain.dy
    # Material setup
    params.user[1]        = 1.0/3.0  # layer half thickness
    # Material parameters
    materials        = MaterialParameters()
    materials.nphase = 2
    materials.Tref   = zeros(Float64, materials.nphase)
    materials.eta0   = zeros(Float64, materials.nphase)
    materials.n      = zeros(Float64, materials.nphase)
    # Material 1
    materials.Tref[1]   = 2.0  # reference flow stress
    materials.n[1]      = 1.0
    materials.eta0[1]   = materials.Tref[1] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[1])  # matrix viscosity
    # Material 2
    materials.Tref[2]   = 2.0/eta_fact # reference flow stress
    materials.n[2]      = 1.0
    materials.eta0[2]   = materials.Tref[2] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[2]) # inclusion viscosity
    # BC definition
    BC = BoundaryConditions()
    # Pure shear - full Dirichlet
    BC.Vx.type_W = 1;     BC.Vx.Dir_W = 1.0*ones(ncy+2)
    BC.Vx.type_E = 1;     BC.Vx.Dir_E =-1.0*ones(ncy+2)
    BC.Vx.type_S =13;     BC.Vx.Dir_S = 0.0*ones(ncx+1)
    BC.Vx.type_N =13;     BC.Vx.Dir_N = 0.0*ones(ncx+1)
    BC.Vy.type_W =13;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
    BC.Vy.type_E =13;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
    BC.Vy.type_S = 1;     BC.Vy.Dir_S =-1.0*ones(ncx+2)
    BC.Vy.type_N = 1;     BC.Vy.Dir_N = 1.0*ones(ncx+2)
    # Evaluate BC's: West/East Vx
    for j=1:ncy+2
        BC.Vx.Dir_W[j] = -domain.xmin*params.Ebg
        BC.Vx.Dir_E[j] = -domain.xmax*params.Ebg
    end
    # Evaluate BC's: South/North Vx
    for i=1:ncx+1
        BC.Vx.Dir_S[i] = 0.0
        BC.Vx.Dir_N[i] = 0.0
    end
    # Evaluate BC's: West/East Vy
    for j=1:ncy+1
        BC.Vy.Dir_W[j] = 0.0
        BC.Vy.Dir_E[j] = 0.0
    end
    # Evaluate BC's: South/North Vy
    for i=1:ncx+2
        BC.Vy.Dir_S[i] = domain.ymin*params.Ebg
        BC.Vy.Dir_N[i] = domain.ymax*params.Ebg
    end
    # Allocate tables
    fields           = Fields2D()
    fields.etac      = zeros(Float64, ncx  , ncy  )
    fields.etav      = zeros(Float64, ncx+1, ncy+1)
    fields.Kc        = zeros(Float64, ncx  , ncy  )
    fields.rhoc      = zeros(Float64, ncx  , ncy  )
    fields.Pc0       = zeros(Float64, ncx+0, ncy+0)
    fields.Pc        = zeros(Float64, ncx+0, ncy+0)
    fields.Vx        = zeros(Float64, ncx+1, ncy+2) # !!! GHOST ROWS
    fields.Vy        = zeros(Float64, ncx+2, ncy+1) # !!! GHOST COLUMNS
    fields.Damc      = ones(Float64, ncx  , ncy  )
    fields.Damv      = ones(Float64, ncx+1, ncy+1)
    Tau              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Eps              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    if BC.periodix==0
        fields.Fx    = zeros(Float64, ncx+1, ncy+0)
    else
        fields.Fx    = zeros(Float64, ncx+0, ncy+0)
    end
    fields.Fy        = zeros(Float64, ncx+0, ncy+1)
    fields.Fp        = zeros(Float64, ncx  , ncy  )
    # Initialise markers
    p                = Markers()
    p.nmx            = 4            # 2 marker per cell in x
    p.nmy            = 4            # 2 marker per cell in y
    InitialiseMarkerPositions!( p, domain )
    # Over allocate markers
    p.phase          = zeros(Float64, p.nmark_max)
    p.T              = zeros(Float64, p.nmark_max)
    p.cellx          = zeros(Int64,   p.nmark_max)
    p.celly          = zeros(Int64,   p.nmark_max)
    # Marker 2 cell mapping
    fields.mpc           = zeros(Float64,(ncx  ,ncy  )) # markers per cell
    fields.mpc_th        = zeros(Float64,(nthreads(), ncx  ,ncy   )) # markers per cell
    fields.phase_perc    = zeros(Float64,(materials.nphase, ncx  ,ncy  )) # markers per cell
    fields.phase_perc_th = zeros(Float64,(nthreads(), materials.nphase, ncx  ,ncy   )) # markers per cell
    # Initial configuration
    SetMarkers!( p, params, domain, α )  # Define phase on markers
    SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
    fields.Vxc    = 0.5.*(fields.Vx[1:end-1,2:end-1] .+ fields.Vx[2:end-0,2:end-1])
    fields.Vyc    = 0.5.*(fields.Vy[2:end-1,1:end-1] .+ fields.Vy[2:end-1,2:end-0])
    # Visualisation
    viz_directory = string( "Figures_", experiment )
    if isdir( viz_directory ) == false
        mkdir( viz_directory )
    end
    path = string( "./", viz_directory, "/" )
    anim = Plots.Animation( path, String[] )
    #------------------------ PINNED PRESSURE
    # Set pinned pressure
    fields.Pc[1,1] = 0.0
    #------------------------ PINNED PRESSURE
    # TIME LOOP
    for it=1:params.nt
        @printf("-------------------------------------\n")
        @printf("---------- Time step #%04d ----------\n", it)
        @printf("-------------------------------------\n")
        # History
        fields.Pc0 .= fields.Pc
        # Count markers and reseed
        @time CountMarkers3!( p, fields, materials, domain )
        @time ReseedMarkers2!( p, fields.mpc, domain)
        nmark = p.nmark
        @time CountMarkers3!( p, fields, materials, domain )
        # Define dt
        if params.dt_var == 1
            maxVx     = maximum(abs.(fields.Vx))
            maxVy     = maximum(abs.(fields.Vy))
            params.dt = params.Courant*(maximum([domain.dx,domain.dy])/maximum([maxVx,maxVy]))
            @printf("dt = %2.2e\n", params.dt)
        end
        # Update cell info on markers
        @time LocateMarkers2!(p, domain)
        # Number equations
        println("Numbering")
        @time NumberingStokes!( fields, BC, domain )
        # Non-linear iterations
        for iter=1:params.niter_nl

            @printf("-------------------------------------\n")
            @printf("---------- Iteration #%04d ----------\n", iter)
            @printf("-------------------------------------\n")

            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            @time Rheology!( fields, Eps, Tau, params, materials, BC, ncx, ncy )
            @time Stress!( Tau, Eps, fields )
            @time ResidualsComp!( fields, Tau, Eps, BC, domain, params )
            println("|Fx| = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy| = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp| = ", norm(fields.Fp)/length(fields.Fp))
            if ( (norm(fields.Fx)/length(fields.Fx) < params.tol_nl) & (norm(fields.Fy)/length(fields.Fy) < params.tol_nl) & (norm(fields.Fp)/length(fields.Fp) < params.tol_nl) )
                println("Mechanical solver converged")
                break
            end

            # Assemble Stokes matrices
            println("Assembly")
            @time Kuu, Kup, Kpu = StokesAssembly( BC, fields, params.Dir_scale, domain )
            # Call solver
            println("Solver")
            @time StokesSolver!( fields, params, domain, materials, BC, Kuu, Kup, Kpu )

            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            @time Stress!( Tau, Eps, fields )
            @time ResidualsComp!( fields, Tau, Eps, BC, domain, params )
            println("|Fx| = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy| = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp| = ", norm(fields.Fp)/length(fields.Fp))
        end
        params.t += params.dt
        # Advect particles
        if params.advection == 1
            @time RungeKutta2!( p, fields, params, BC, domain )
            # Deform box
            if params.ALE==1
                dx, dy = PureShearBoxUpdate!( domain, BC, params )
                SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
            end
        end
    end
    # Compute dissipation
    psi1  = Eps.II .* Tau.II # by invariant product
    psi2  = Eps.xx.*Tau.xx .+ Eps.yy.*Tau.yy .+ 2.0.*Eps.xy_c.*Tau.xy_c # by tensor contraction
    psi1_inc =  mean( fields.phase_perc[2,:,:] .* psi1 )
    psi2_inc =  mean( fields.phase_perc[2,:,:] .* psi2 )
    return fields, domain, psi1_inc, psi2_inc
end

####################################################
####################################################
####################################################
α        = deg2rad(65)
ncx, ncy = 200, 200
eta_val  = [1, 1.1, 1.2, 1.5, 2, 3, 4, 5, 10, 50, 100, 500, 1000, 1250, 5000, 7500, 10000, 11250, 11500]
diss1v   = zeros(length(eta_val))
diss2v   = zeros(length(eta_val))
# Call solver
for ieta=1:length(eta_val)
    @time fields, domain, diss1v[ieta], diss2v[ieta] = main( ncx, α, eta_val[ieta] )
end
# Visualize
step  = 10
ieta  = length(eta_val)
@time fields, domain, diss1v[ieta], diss2v[ieta] = main( ncx, α, eta_val[end] ) # one extra solve for visualisation
Vxc   = 0.5*(fields.Vx[1:end-1,:2:end-1] .+ fields.Vx[2:end-0,:2:end-1])
Vyc   = 0.5*(fields.Vy[2:end-1,:1:end-1] .+ fields.Vy[2:end-1,:2:end-0])
Vxq   = Vxc[1:step:end,1:step:end]
Vyq   = Vyc[1:step:end,1:step:end]
# Make figure
p1 = plot( log10.(eta_val), diss1v, label="ψ₁")
p1 = plot!(log10.(eta_val), diss2v, label="ψ₂", xlabel="ηₘ/ηᵢ", ylabel="ψ",legend=:topright)
p2 = heatmap(domain.xc, domain.yc, Array(fields.Pc)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="P")
xc2q, yc2q = ndgrid( domain.xc[1:step:end], domain.yc[1:step:end] )
p2 = quiver!( xc2q[:], yc2q[:], quiver=(Vxq[:], Vyq[:]), lw=0.1, c=:blue)
display( plot(p1, p2) )


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "SlideShow"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
   showdocs([slideshow, slide])
)

$(vskip(4em))

"""

function main(window)
    docpage(intro)
end


text:
"""
 Calculate contribution of activations towards loss according to https://arxiv.org/abs/1611.06440.

High value of `contribution` means that the neuron of that index has a high contribution to the loss score.

Short summary is that the first order taylor approximation of the optimization problem: "which neurons shall I remove to minimize impact on the loss function" boils down to: "the ones which minimize abs(gradient * activation)" (assuming parameter independence).
"""
mutable struct ActivationContribution <: AbstractMutableComp
    layer
    contribution
    aggmethod
end
ActivationContribution(l::AbstractMutableComp) = ActivationContribution(l, zeros(Float32, nout(l)), Ewma(0.05))
ActivationContribution(l) = ActivationContribution(l, missing, Ewma(0.05))

layer(m::ActivationContribution) = layer(m.layer)
layertype(m::ActivationContribution) = layertype(m.layer)
wrapped(m::ActivationContribution) = m.layer
NaiveNASlib.minΔninfactor(m::ActivationContribution) = minΔninfactor(wrapped(m))
NaiveNASlib.minΔnoutfactor(m::ActivationContribution) = minΔnoutfactor(wrapped(m))

function(m::ActivationContribution)(x...)
    act = wrapped(m)(x...)

    return hook(act) do grad
        m.contribution = agg(m.aggmethod, m.contribution, mean_squeeze(abs.(act .* grad).data, actdim(ndims(act))))
        return grad
    end
end

actdim(nd::Integer) = nd - 1

function NaiveNASlib.mutate_outputs(m::ActivationContribution, outputs::AbstractVector{<:Integer})
    m.contribution = select(m.contribution, 1 => outputs)
    mutate_outputs(wrapped(m), outputs)
end

"""
    mean_squeeze(x, dimkeep)

Return mean value of `x` along all dimensions except `dimkeep` as a 1D array (singleton dimensions are removed).
"""
function mean_squeeze(x, dimkeep)
    dims = filter(i -> i != dimkeep, 1:ndims(x))
    return dropdims(mean(x, dims=dims), dims=Tuple(dims))
end

# To peel the onion...
neuron_value(v::AbstractVertex) = neuron_value(base(v))
neuron_value(v::InputSizeVertex) = ones(nout(v))
neuron_value(v::CompVertex) = neuron_value(v.computation)
neuron_value(m::AbstractMutableComp) = neuron_value(wrapped(m))
neuron_value(m::ActivationContribution) = m.contribution
neuron_value(l) = neuron_value(layertype(l), l)

# Default: highest mean of abs of weights + bias. Not a very good metric, but should be better than random
# Maybe do something about state in recurrent layers as well, but CBA to do it right now
neuron_value(t::FluxParLayer, l) = mean_squeeze(abs.(weights(l)), outdim(l)) + abs.(bias(l))

"""
    Ewma(α::Real)

Exponential moving average.

Parameter `α` acts as a forgetting factor, i.e larger values means faster convergence but more noisy estimate.
"""
struct Ewma
    α
    function Ewma(α::Real)
        0 <= α <= 1 || error("α must be between 0 and 1, was $α")
        new(α)
    end
end
agg(m::Ewma, x ,y) = m.α .* x .+ (1 - m.α) .* y
agg(m::Ewma, ::Missing, y) = y


text:
int_rules_5_3_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.3.1*(a+b*arctan(c*x^n))^p =#
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~(n'))) ^ ~(p'), ~x) => ~x * (~a + ~b * atan(~c * (~x) ^ ~n)) ^ ~p - ~b * ~c * ~n * ~p * Antiderivative(((~x) ^ ~n * (~a + ~b * atan(~c * (~x) ^ ~n)) ^ (~p - 1)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (IGtQ(~p, 0) && (EqQ(~n, 1) || EqQ(~p, 1)))
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~(n'))) ^ ~(p'), ~x) => ~x * (~a + ~b * acot(~c * (~x) ^ ~n)) ^ ~p + ~b * ~c * ~n * ~p * Antiderivative(((~x) ^ ~n * (~a + ~b * acot(~c * (~x) ^ ~n)) ^ (~p - 1)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (IGtQ(~p, 0) && (EqQ(~n, 1) || EqQ(~p, 1)))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandIntegrand(((~a + (im * ~b * log(1 - im * ~c * (~x) ^ ~n)) / 2) - (im * ~b * log(1 + im * ~c * (~x) ^ ~n)) / 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandIntegrand(((~a + (im * ~b * log(1 - (im * (~x) ^ -(~n)) / ~c)) / 2) - (im * ~b * log(1 + (im * (~x) ^ -(~n)) / ~c)) / 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative((~a + ~b * acot((~x) ^ -(~n) / ~c)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative((~a + ~b * atan((~x) ^ -(~n) / ~c)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * atan(~c * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && FractionQ(~n))
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * acot(~c * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~p, 1) && FractionQ(~n))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~(n'))) ^ ~p, ~x) => Unintegrable((~a + ~b * atan(~c * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~(n'))) ^ ~p, ~x) => Unintegrable((~a + ~b * acot(~c * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.3 Inverse tangent/5.3.1 (a+b arctan(c x^n))^p.jl


text:
using DataFrames;
using Distributions;

#source_dir = typeof(Base.source_dir()) == Void ? joinpath(Pkg.dir("LocalAlignment"), "src") : Base.source_dir()
source_dir = pwd()

include(joinpath(source_dir, "src", "algorithms", "smithwaterman.jl"));
include(joinpath(source_dir, "src", "utils", "parse.jl"));
include(joinpath(source_dir, "src", "utils", "benchmarking.jl"));

sfp = ARGS[1];
T = ARGS[2];
N = ARGS[3];
pos_a_fp = ARGS[4];
neg_a_fp = ARGS[5];

"""
optimize(S::DataFrame, T::Float64, N::Int, g::Float64, e::Float64, pos_a_fp::String, neg_a_fp::String)
Optimize the scoring matrix S with a Monte Carlo-Markov Chain algorithm with simulated
annealing paramter T over N iterations. Score postive alignments in POS_A_FP and
negative alignments in NEG_A_FP with gap penalty G and extension penalty E. Objective function
is the sum of TPRs for FPRs = [0.0, 0.1, 0.2, 0.3, 0.4].
"""
function optimize(S, T, N, g, e, pos_a_fp, neg_a_fp)

    pos_aligns = parse_alignments(pos_a_fp);
    neg_aligns = parse_alignments(neg_a_fp);

    Scurr = S

    prev_pos_scores = []
    prev_neg_scores = []
    for pos in pos_aligns
        push!(prev_pos_scores, score_alignment(pos[1], pos[2], Scurr, g, e)[1])
    end

    for neg in neg_aligns
        push!(prev_neg_scores, score_alignment(neg[1], neg[2], Scurr, g, e)[1])
    end

    Jprev = eval_objective_function(prev_pos_scores, prev_neg_scores)

    for i in 1:N
        println(string("Iteration: ", i));
        curr_pos_scores = []
        curr_neg_scores = []

        Sp = gen_new_matrix(Scurr)

        for pos in pos_aligns
            push!(curr_pos_scores, score_alignment(pos[1], pos[2], Sp, g, e)[1])
        end

        for neg in neg_aligns
            push!(curr_neg_scores, score_alignment(neg[1], neg[2], Sp, g, e)[1])
        end

        Jcurr = eval_objective_function(curr_pos_scores, curr_neg_scores)

        if Jcurr >= Jprev
            Scurr = Sp
            prev_neg_scores = curr_neg_scores
            prev_pos_scores = curr_pos_scores
            Jprev = Jcurr
        else
            p = e^( (Jcurr - Jprev) * i ^T)
            if rand(Uniform()) < p
                Scurr = Sp
                prev_neg_scores = curr_neg_scores
                prev_pos_scores = curr_pos_scores
                Jprev = Jcurr
            else
                continue;
            end
        end

    end

    Scurr

end

"""
eval_objective_function(pos_scores::Array, neg_scores::Array)
Evaluate the objective function for sequence optimization on positive scores POS_SCORES and
negative scores NEG_SCORES.
"""
function eval_objective_function(pos_scores, neg_scores)

    val = 0.0
    fpr = [0.0, 0.1, 0.2, 0.3]

    for fp in fpr
        thresh = find_tp_threshold(neg_scores, fp)
        val += length(pos_scores[pos_scores .> thresh]) / length(pos_scores)
    end

    return val
end

"""
score_alignment(align1::String, align2::String, S::DataFrame, g::Float64, e::Float64)
Score the alignment stored in ALIGN1 and ALIGN2 according to the substitution score
matrix S, gap penalty G and extension penalty E.
"""
function score_alignment(align1, align2, S, g, e)

    score = 0.0
    penalty = 0.0

    for i in 1:length(align1)

        if align1[i] != '-' && align2[i] != '-'
            score += S[S[:, :AA] .== string(uppercase(align1[i])),
                            Symbol(uppercase(align2[i]))];
        end

    end

    curr_e = 0
    for i in 1:length(align1)

        if align1[i] == '-' || align2[i] == '-'
            if curr_e == 0
                penalty += g
            else
                penalty += e
            end
            curr_e += 1
        else
            curr_e = 0;
        end

    end

    score = score - penalty
    return score

end

"""
parse_alignment(fp::String)
Parse alignment file in FP and return list of alignment pairs.
"""
function parse_alignments(fp)

    alignments = []

    f = open(fp, "r");
    lines = collect(readlines(f))

    for i in 1:3:length(lines)

        align1 = strip(lines[i]);
        align2 = strip(lines[i+1]);
        push!(alignments, [align1, align2]);

    end

    alignments

end

"""
gen_new_matrix(S::DataFrame)
Generate a new symmetric substitution matrix from S by adding random normal noise
to each subsitution value.
"""
function gen_new_matrix(S)

    Sp = deepcopy(S)
    Sp = Array(Sp[:,2:end])
    aa = S[:,:AA]

    for i in 1:size(Sp, 1)
        for j in size(Sp, 1):-1:i
            Sp[i, j] += randn()
            Sp[j, i] = Sp[i,j]
        end
    end

    Sp = DataFrame(Sp);
    Sp = hcat(aa, Sp);
    names!(Sp, names(S));

    Sp

end

S = parse_score_matrix(sfp);

S2 = optimize(S, T, N, pos_a_fp, neg_a_fp);

writetable("opt_scor_matrix.txt", S2, quotemark=' ')


text:
using Statistics: mean
using Oceananigans.Fields: Field, location
using Oceananigans.OutputReaders: FieldTimeSeries, FieldDataset

"""
    coarse_grain(field::Field, new_grid; dims=3)

Coarse grain a `field` onto a `new_grid` along `dims`. Returns a new `Field`.
"""
function coarse_grain(field::Field{X, Y, Center}, new_grid; dims=3) where {X, Y}

    # TODO: Generalize to x and y.
    @assert dims == 3

    # TODO: Generalize `coarse_grain` to non-integer ratios.
    r = field.grid.Nz / new_grid.Nz
    @assert isinteger(r)
    r = Int(r)

    coarse_field = Field(location(field)..., field.architecture, new_grid, field.boundary_conditions)

    coarse_data = zeros(size(coarse_field))
    for K in 1:new_grid.Nz
        k1 = r * (K-1) + 1
        k2 = r * K
        coarse_data[:, :, K] .= mean(field[:, :, k1:k2])
    end

    set!(coarse_field, coarse_data)

    return coarse_field
end

function coarse_grain(field::Field{X, Y, Face}, new_grid; dims=3) where {X, Y}

    # TODO: Generalize to x and y.
    @assert dims == 3

    r = field.grid.Nz / new_grid.Nz

    coarse_field = Field(location(field)..., field.architecture, new_grid, field.boundary_conditions)

    coarse_data = zeros(size(location(field), new_grid))

    # Left-most and right-most faces must match between the two grids.
    coarse_data[:, :, 1] .= field.data[:, :, 1]
    coarse_data[:, :, new_grid.Nz+1] .= field.data[:, :, field.grid.Nz+1]

    for K in 2:new_grid.Nz
        c = 1 + (K-1) * r
        coarse_data[K] = (floor(c)+1 - c) * field.data[floor(Int, c)] + (c - floor(c)) * field.data[floor(Int, c)+1]
    end

    set!(coarse_field, coarse_data)

    return coarse_field
end

function coarse_grain(fts::FieldTimeSeries, new_grid)
    fts_new = FieldTimeSeries(new_grid, location(fts), fts.times; ArrayType=Array{Float32})

    Nt = size(fts, 4)
    for n in 1:Nt
        fts_new.data[:, :, :, n] .= coarse_grain(fts[n], new_grid).data
    end

    return fts_new
end

function coarse_grain(fds::FieldDataset, new_grid)
    coarse_fields = Dict{String, FieldTimeSeries}(
        name => coarse_grain(fts, new_grid)
        for (name, fts) in fds.fields
    )
    return FieldDataset(coarse_fields, fds.metadata, fds.filepath)
end


text:
using PictureSegmentation, Images, TestImages
using Test

function generate_3d_img()
    img=Array{RGB{Float64},3}(undef,100,100,100)
    for i in CartesianIndices(img)
        if (i[1]-50)^2 + (i[2]-50)^2 < (i[3]/4)^2
            img[i]=RGB(i[3]/100,0,1-(i[3]/100))
        else
            img[i]=RGB(0,0,0)
        end
    end
    for i = 10:40
        for j = 10:40
            for k = 10:40
                img[i,j,k] = RGB(0.25,0.66,0.02)
            end
        end
    end
    return img
end

@testset "PictureSegmentation.jl" begin
    include("growcut.jl")
    include("growcut3d.jl")
    include("seedpixel.jl")
end


text:
"""
Routines for single-node inference that aren't specific
to any particular method of inference (e.g MCMC, DeterministicVI),
yet are also not about the statistical model (i.e., the Model
module).
Currently what's in here are routines that effectively
truncate the Gaussians in the model.
TODO: rename this module to something more meaningful
"""
module Infer

import WCS
using StaticArrays

import ..Config
using ..Model
import ..Log

"""
Computes the nearby light sources in the catalog for each of the target
sources.

Arguments:
    target_sources: indexes of astronomical objects in the catalog to infer
    catalog: astronomical objects appearing the images
    images: astronomical images
"""
function find_neighbors(target_sources::Vector{Int64},
                        catalog::Vector{CatalogEntry},
                        images::Vector{Image})
    psf_width_ub = zeros(NUM_BANDS)
    for img in images
        psf_width = Model.get_psf_width(img.psf)
        psf_width_ub[img.b] = max(psf_width_ub[img.b], psf_width)
    end

    epsilon_lb = fill(Inf, NUM_BANDS)
    for img in images
        epsilon = img.sky[div(img.H, 2), div(img.W, 2)]
        epsilon_lb[img.b] = min(epsilon_lb[img.b], epsilon)
    end

    radii_map = zeros(length(catalog))
    for s in 1:length(catalog)
        ce = catalog[s]
        for img in images
            radius_pix = Model.choose_patch_radius(ce, img, width_scale=1.2)
            radii_map[s] = max(radii_map[s], radius_pix)
        end
        @assert radii_map[s] <= 25
    end

    dist(ra1, dec1, ra2, dec2) = (3600 / 0.396) * max(abs(dec2 - dec1), abs(ra2 - ra1))

    neighbor_map = Vector{Int64}[Int64[] for s in target_sources]

    # If this loop isn't super fast in pratice, we can tile (the sky, not the
    # images) or build a spatial index with a library before distributing
    for ts in 1:length(target_sources)
        s = target_sources[ts]
        ce = catalog[s]

        for s2 in 1:length(catalog)
            ce2 = catalog[s2]
            ctrs_dist = dist(ce.pos[1], ce.pos[2], ce2.pos[1], ce2.pos[2])

            if s2 != s && ctrs_dist < radii_map[s] + radii_map[s2]
                push!(neighbor_map[ts], s2)
            end
        end
    end

    neighbor_map
end

function get_sky_patches(images::Vector{Image},
                         catalog::Vector{CatalogEntry};
                         radius_override_pix=NaN)
    N = length(images)
    S = length(catalog)
    patches = Matrix{SkyPatch}(S, N)

    for n = 1:N
        img = images[n]

        for s=1:S
            world_center = catalog[s].pos
            pixel_center = WCS.world_to_pix(img.wcs, world_center)
            wcs_jacobian = Model.pixel_world_jacobian(img.wcs, pixel_center)
            radius_pix = Model.choose_patch_radius(catalog[s], img, width_scale=1.2)
            @assert radius_pix <= 25
            if !isnan(radius_override_pix)
                radius_pix = radius_override_pix
            end

            hmin = max(0, floor(Int, pixel_center[1] - radius_pix - 1))
            hmax = min(img.H - 1, ceil(Int, pixel_center[1] + radius_pix - 1))
            wmin = max(0, floor(Int, pixel_center[2] - radius_pix - 1))
            wmax = min(img.W - 1, ceil(Int, pixel_center[2] + radius_pix - 1))

            # some light sources are so far from some images that they don't
            # overlap at all
            H2 = max(0, hmax - hmin + 1)
            W2 = max(0, wmax - wmin + 1)

            # all pixels are active by default
            active_pixel_bitmap = trues(H2, W2)

            patches[s, n] = SkyPatch(world_center,
                                     radius_pix,
                                     img.psf,
                                     wcs_jacobian,
                                     pixel_center,
                                     SVector(hmin, wmin),
                                     active_pixel_bitmap)
        end
    end

    patches
end

"""
Record which pixels in each patch will be considered when computing the
objective function.

Non-standard arguments:
  noise_fraction: The proportion of the noise below which we will remove pixels.
"""
function load_active_pixels!(config::Config,
                             images::Vector{Image},
                             patches::Matrix{SkyPatch};
                             exclude_nan=true,
                             noise_fraction=0.5)
    S, N = size(patches)

    for n = 1:N, s=1:S
        img = images[n]
        p = patches[s,n]

        # (h2, w2) index the local patch, while (h, w) index the image
        H2, W2 = size(p.active_pixel_bitmap)
        for w2 in 1:W2, h2 in 1:H2
            h = p.bitmap_offset[1] + h2
            w = p.bitmap_offset[2] + w2

            # skip masked pixels
            if isnan(img.pixels[h, w]) && exclude_nan
                p.active_pixel_bitmap[h2, w2] = false
                continue
            end

            # include pixels that are close, even if they aren't bright
            sq_dist = (h - p.pixel_center[1])^2 + (w - p.pixel_center[2])^2
            if sq_dist < config.min_radius_pix^2
                p.active_pixel_bitmap[h2, w2] = true
                continue
            end

            # if this pixel is bright, let's include it
            # (in the future we may want to do something fancier, like
            # fitting an elipse, so we don't include nearby sources' pixels,
            # or adjusting active pixels during the optimization)
            # Note: This is risky because bright pixels are disproportionately likely
            # to get included, even if it's because of noise. Therefore it's important
            # to keep the noise fraction pretty low.
            threshold = img.nelec_per_nmgy[h] * img.sky[h, w] * (1. + noise_fraction)
            p.active_pixel_bitmap[h2, w2] = img.pixels[h, w] > threshold
        end
    end
end

# legacy wrapper
function load_active_pixels!(images::Vector{Image},
                             patches::Matrix{SkyPatch};
                             exclude_nan=true,
                             noise_fraction=0.5)
    load_active_pixels!(
        Config(),
        images,
        patches,
        exclude_nan=exclude_nan,
        noise_fraction=noise_fraction,
    )
end

# The range of image pixels in a vector of patches
function get_active_pixel_range(
    patches::Matrix{SkyPatch}, sources::Vector{Int}, n::Int)
    H_min = minimum([ p.bitmap_offset[1] + 1 for p in patches[sources, n] ])
    W_min = minimum([ p.bitmap_offset[2] + 1 for p in patches[sources, n] ])
    H_max = maximum([ p.bitmap_offset[1] + size(p.active_pixel_bitmap, 1)
                      for p in patches[sources, n] ])
    W_max = maximum([ p.bitmap_offset[2] + size(p.active_pixel_bitmap, 2)
                      for p in patches[sources, n] ])
    H_min, W_min, H_max, W_max
end

# Is a pixel (h, w) in whole-image coordinates an active pixel in the patch p?
function is_pixel_in_patch(h::Int, w::Int, p::SkyPatch)
    hp = h - p.bitmap_offset[1]
    wp = w - p.bitmap_offset[2]
    in_patch =
        (hp > 0) &
        (wp > 0) &
        (hp <= size(p.active_pixel_bitmap, 1)) &
        (wp <= size(p.active_pixel_bitmap, 2))
    if !in_patch
        return false
    else
        return p.active_pixel_bitmap[hp, wp]
    end
end

end  # module


text:
module ReinforcementLearningZoo

const RLZoo = ReinforcementLearningZoo
export RLZoo

using CircularArrayBuffers
using ReinforcementLearningBase
using ReinforcementLearningCore
using Setfield: @set
using StableRNGs
using Logging
using Flux.Losses
using Dates
using Random
using Random: shuffle
using CUDA
using Zygote
using Zygote: ignore
using Flux
using Flux: onehot, normalise
using StatsBase
using StatsBase: sample, Weights, mean
using LinearAlgebra: dot
using MacroTools
using Distributions: Categorical, Normal, logpdf
using StructArrays

include("patch.jl")
include("algorithms/algorithms.jl")

using Requires

# dynamic loading environments
function __init__()
    @require ReinforcementLearningEnvironments = "25e41dd2-4622-11e9-1641-f1adca772921" begin
        include("experiments/rl_envs/rl_envs.jl")
        @require ArcadeLearningEnvironment = "b7f77d8d-088d-5e02-8ac0-89aab2acc977" include("experiments/atari/atari.jl")
        # @require OpenSpiel = "ceb70bd2-fe3f-44f0-b81f-41608acaf2f2" include("experiments/open_spiel/open_spiel.jl")
    end
end

end # module


text:
using CxxWrap, LibGit2

# Path to the CxxWrap dependencies
cxxhome = dirname(dirname(CxxWrap.jlcxx_path))

url = "https://github.com/hildebrandmw/pcm"
branch = "master"
localdir = joinpath(@__DIR__, "pcm")
ispath(localdir) || LibGit2.clone(url, localdir; branch = branch)

juliahome = dirname(Base.Sys.BINDIR)
run(`make JULIA_HOME=$juliahome CXXWRAP_HOME=$cxxhome -j all `)


text:
@testset "pso" begin
    obj = pso.Objective(pso.TestFunctions.sphere, [-5.0, 5.0], 3)
    @testset "Objective" begin
        # obj1 = pso.Objective(pso.TestFunctions.sphere, [-5.0, 5.0], 3)
        @test obj.search_space == [[-5.0, 5.0], [-5.0, 5.0], [-5.0, 5.0]]
    end
    @testset "update_position" begin
        position = [1.0, 2.0]
        pso.update_position!(position, [0.75, -0.25], [0.0, 0.0], [0.0, 0.0])
        @test position == [1.75, 1.75]
    end
    # @testset "update_velocity" begin
    #     velocity = [1.0, 2.0]
    #     Optimizers.update_velocity!([1.0, 2.0], velocity, )
    #     @test velocity == [1.75, 1.75]
    # end
    @testset "confinements" begin
        position = [6.0, 2.0]
        velocity = [1.0, -2.0]
        pso.confinement!(position, velocity, [[-5.0, 5.0], [-5.0, 5.0]])
        @test position == [5.0, 2.0]
        @test velocity == [-0.5, -2.0]
    end
    @testset "evaluate" begin
        position = [1.0, 1.0]
        pos_best = [2.0, -2.0]
        res_best = [4.0, 3.0, pso.TestFunctions.sphere(pos_best), 5.0]
        pso.evaluate!(position, res_best, pos_best, 3, <, pso.TestFunctions.sphere)
        @test position == [1.0, 1.0]
        @test pos_best == [1.0, 1.0]
        @test res_best == [4.0, 3.0, pso.TestFunctions.sphere(position), 5.0]
    end
    @testset "get_localbest" begin
        @test pso.get_localbest([0.2, 0.5, 0.6, 0.1], [1, 3], <) == 1
        @test pso.get_localbest([0.2, 0.5, 0.6, 0.1], [1, 3], >) == 3
        @test pso.get_localbest([0.5, 0.5, 0.5, 0.4, 0.5], [2, 4, 5], <) == 4
    end
    neig = pso.GlobalNeighbourhood(20)
    optimizer = pso.PSO(obj, neig)
    @testset "PSO object" begin
        @test size(optimizer.position) == (20,)
        @test size(optimizer.velocity) == (20,)
        # for el in optimizer.position
        #     @test -5 ≤ el ≤ 5
        # end
        # for el in optimizer.position_matrix.+optimizer.velocity_matrix
        #     @test -5 ≤ el ≤ 5
        # end
    end
    @testset "optimize! function" begin
        pso.optimize!(optimizer, 500)
        result = pso.getoptimum(optimizer)
        @test all(result[1] .< [1e-10, 1e-10, 1e-10])
        @test result[2] < 1e-10
    end
end


text:
##########################################################
#           Implementación de Backpropagation
# Basado en el algoritmo 4 y 5, sección 5.3 de la memoria
##########################################################
export backpropagation!

using Random

VectorOrMatrix = Union{Matrix,Vector}
function descent_gradient(
    neural_network::AbstractOneLayerNeuralNetwork,
    X_train :: VectorOrMatrix,
    Y_train :: VectorOrMatrix,
    activation_function,
    derivative_activation_funcion,
    batch_size :: Int = 32,
    )
    len, _ = size(X_train)
    if len < batch_size
        error("batch_size should be equal or smaller than the size of train dataset,
        but $batch_size  > $len")
    end

    # Derivadas parciales a calcular
    n_1,d = size(neural_network.W1)
    s,n = size(neural_network.W2)

    partial_W1 = zeros(Float64, (n_1,d))
    partial_W2 = zeros(Float64, (s,n))
    derivative_B = zeros(Float64, (s,n))
    # Variables auxiliares para reducir número de operaciones
    sensibilities_img = zeros(Float64, n)

    index = first(randperm(len),batch_size)
    # Para cada muestra del entrenamiento
    for train_index in index
        # 1. Forward propagation almacenando resultado relevantes
        𝛅 = neural_network.W1 * push!(copy(X_train[train_index,:]),1)
        sensibilities_img = map(activation_function, 𝛅)
        derivative_sensibilities_img= map(derivative_activation_funcion, 𝛅)

        forward_propagation_x = neural_network.W2 * sensibilities_img
        error = forward_propagation_x .- Y_train[train_index,:]

        # 2. Parcial de B (W_2)
        for v in 1:n
            for w in 1:s
                partial_W2[w,v] = partial_W2[w,v] + error[w]*sensibilities_img[v]
            end
        end

        # 3. Parcial W_1
        for i in 1:n
            for k in 1:s
            derivative_B[k,i] =
                neural_network.W2[k,i]*derivative_sensibilities_img[i]
            end
        end

        for v in 1:n
            for k in 1:s
                difference_times_derivative = error[k]*derivative_B[k,s]
                partial_W2[k, v] += error[k]*derivative_sensibilities_img[v]
                partial_W1[v,d] += difference_times_derivative

                for u in 1:(d-1)
                    partial_W1[v,u] += difference_times_derivative * X_train[train_index, u]
                end
            end
        end

     end
     return partial_W1, partial_W2
end

"""
    function backpropagation!(
        neural_network::AbstractOneLayerNeuralNetwork,
        X_train :: Matrix,
        Y_train :: Vector,
        activation_function,
        derivative_activation_funcion,
        batch_size :: Int = 32,
        η :: Float64 =  0.1
    )

Compute backpropagation
"""
function backpropagation!(
    neural_network::AbstractOneLayerNeuralNetwork,
    X_train :: Matrix,
    Y_train :: Vector,
    activation_function,
    derivative_activation_funcion,
    batch_size :: Int = 32,
    η :: Float64 =  0.1
)
    ∂w1, ∂w2 = descent_gradient(neural_network,
         X_train, Y_train,
         activation_function, derivative_activation_funcion,
        batch_size)
        neural_network.W1 -= η*∂w1
        neural_network.W2 -= η*∂w2
    return neural_network
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

polyorder = 1
numqp = required_quadrature_order(polyorder)
basis = TensorProductBasis(2, polyorder)
levelset = InterpolatingPolynomial(1, basis)
levelsetcoeffs = [-1.0, -1.0, 1.0, 1.0]
update!(levelset, levelsetcoeffs)

xL, xR = [-1.0, -1.0], [1.0, 1.0]
pquad = area_quadrature(levelset, +1, xL, xR, numqp)
nquad = area_quadrature(levelset, -1, xL, xR, numqp)
squad = surface_quadrature(levelset, xL, xR, numqp)

lambda, mu = 1.0, 2.0
penalty = 1.0
dx = 0.1
e11 = dx / 2.0
e22 = -lambda / (lambda + 2mu) * e11
dy = e22
penalty = 1.0
stiffness = CutCell.plane_strain_voigt_hooke_matrix(lambda, mu)
cellmap = CutCell.CellMap([0.0, 0.0], [2.0, 1.0])

pbf = CutCell.bilinear_form(basis, pquad, stiffness, cellmap)
nbf = CutCell.bilinear_form(basis, nquad, stiffness, cellmap)

normals = repeat([1.0, 0.0], 1, numqp)
top =
    CutCell.incoherent_traction_operator(basis, squad, squad, normals, stiffness, cellmap)

facescale = CutCell.scale_area(cellmap, normals)
components = CutCell.tangents(normals)
mm =
    penalty *
    CutCell.interface_component_mass_matrix(basis, squad, squad, components, facescale)

eta = -1.0
BF = [
    pbf zeros(8, 8)
    zeros(8, 8) nbf
]
T1 = -0.5 * [
    -top -top
    +top +top
]
T2 = eta * T1'

MM = [
    +mm -mm
    -mm +mm
]

K = BF + T1 + T2 + MM
@test norm(K - K') > 1.0

rhs = zeros(16)

CutCell.apply_dirichlet_bc!(K, rhs, [5, 6], 1, 0.0, 2)
CutCell.apply_dirichlet_bc!(K, rhs, [5], 2, 0.0, 2)
CutCell.apply_dirichlet_bc!(K, rhs, [3, 4], 1, dx, 2)

sol = K \ rhs
disp = reshape(sol, 2, :)

testdisp = [
    0.0 0.0 dx dx 0.0 0.0 dx dx
    0.0 dy 0.0 dy 0.0 dy 0.0 dy
]

@test allapprox(disp, testdisp, 1e2 * eps())

###############################################################################
# Use different material properties

L = 1.0
W = 1.0
lambda1, mu1 = 1.0, 2.0
lambda2, mu2 = 3.0, 4.0
penaltyfactor = 1e2
nelmts = 1
dx = L / nelmts
dy = W / nelmts
penalty = penaltyfactor / dx * (lambda1 + mu1 + lambda2 + mu2)*0.5

applydisplacement = 0.01

stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)
polyorder = 1
numqp = 2
basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh([0.0, 0.0], [L, W], [nelmts, nelmts], basis)
levelset = InterpolatingPolynomial(1, basis)
x0 = [0.5, 0.0]
normal = [1.0, 0.0]
levelsetcoeffs =
    CutCell.levelset_coefficients(x -> plane_distance_function(x, normal, x0), mesh)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
interfacecondition =
    CutCell.incoherent_interface_condition(basis,interfacequads,stiffness,cutmesh,penalty)

onbottomboundary(x) = x[2] ≈ 0.0 ? true : false
onrightboundary(x) = x[1] ≈ L ? true : false
ontopboundary(x) = x[2] ≈ W ? true : false
onleftboundary(x) = x[1] ≈ 0.0 ? true : false

leftbc = CutCell.DisplacementComponentCondition(
    x -> 0.0,
    basis,
    facequads,
    stiffness,
    cutmesh,
    onleftboundary,
    [1.0, 0.0],
    penalty,
)
bottombc = CutCell.DisplacementComponentCondition(
    x -> 0.0,
    basis,
    facequads,
    stiffness,
    cutmesh,
    onbottomboundary,
    [0.0, 1.0],
    penalty,
)
rightbc = CutCell.DisplacementComponentCondition(
    x -> applydisplacement,
    basis,
    facequads,
    stiffness,
    cutmesh,
    onrightboundary,
    [1.0, 0.0],
    penalty,
)

sysmatrix = CutCell.SystemMatrix()
sysrhs = CutCell.SystemRHS()

CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
CutCell.assemble_interface_condition!(sysmatrix,interfacecondition,cutmesh)
CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,leftbc,cutmesh)
CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,bottombc,cutmesh)
CutCell.assemble_penalty_displacement_bc!(sysmatrix,sysrhs,rightbc,cutmesh)

matrix = CutCell.make_sparse(sysmatrix,cutmesh)
rhs = CutCell.rhs(sysrhs,cutmesh)

sol = matrix\rhs

displacement = reshape(sol,2,:)

E1 = mu1*(lambda1+mu1)/(lambda1+2mu1)
E2 = mu2*(lambda2+mu2)/(lambda2+2mu2)

uI = E1/(E1+E2)*applydisplacement

u10 = uI + (applydisplacement - uI)/0.5*(-0.5)
Δx = applydisplacement

u12 = -lambda1/(lambda1+2mu1)*2*(Δx - uI)
u22 = -lambda2/(lambda2+2mu2)*2*uI

testdisplacement = [u10   u10   Δx   Δx   0.0   0.0   2uI   2uI
                    0.0   u12   0.0  u12  0.0   u22   0.0   u22]

@test allapprox(displacement,testdisplacement,1e2eps())


text:
# We define a new broadcast style for domains, because they may
# represent continuous sets

"The broadcast style associated with domains"
struct DomainSetStyle <: Base.Broadcast.BroadcastStyle end

Base.BroadcastStyle(::Type{<:Domain}) = DomainSetStyle()

Base.broadcastable(d::Domain) = d

# DomainSetStyle doesn't mix with ArrayStyle
Base.BroadcastStyle(::DomainSetStyle, ::Base.Broadcast.AbstractArrayStyle) = DomainSetStyle()
Base.BroadcastStyle(::Base.Broadcast.AbstractArrayStyle, ::DomainSetStyle) = DomainSetStyle()

import Base.Broadcast: broadcasted

broadcasted(::DomainSetStyle, ::typeof(+), a::Union{Number,AbstractArray}, d::Domain) =
    map_domain(Translation(a), d)
broadcasted(::DomainSetStyle, ::typeof(+), d::Domain, a::Union{Number,AbstractArray}) =
    map_domain(Translation(a), d)

broadcasted(::DomainSetStyle, ::typeof(-), a::Union{Number,AbstractArray}, d::Domain) =
    map_domain(AffineMap(-1, a), d)
broadcasted(::DomainSetStyle, ::typeof(-), d::Domain, a::Union{Number,AbstractArray}) =
    map_domain(Translation(-a), d)
broadcasted(::DomainSetStyle, ::typeof(-), d::Domain{T}) where {T} =
    map_domain(LinearMap{T}(-1), d)

broadcasted(::DomainSetStyle, ::typeof(*), a::Number, d::Domain{T}) where {T} =
    map_domain(LinearMap{T}(a), d)
broadcasted(::DomainSetStyle, ::typeof(*), d::Domain{T}, a::Number) where {T} =
    map_domain(LinearMap{T}(a), d)

broadcasted(::DomainSetStyle, ::typeof(/), d::Domain, a::Number) =
    mapped_domain(LinearMap(a), d)

broadcasted(::DomainSetStyle, ::typeof(\), a::Number, d::Domain) =
    mapped_domain(LinearMap(a), d)

broadcasted(::DomainSetStyle, m::AbstractMap, d::Domain) = map_domain(m, d)

broadcasted(::DomainSetStyle, fun::Function, d::Domain{T}) where {T} =
    convert(Map{T}, fun).(d)

# Intercept broadcast applied to `in`, e.g. in.(A, d).
# This gives domains an opportunity to provide a more efficient implementation
# when invoked with a set of points `A` at once, especially if `A`
# has particular structure. A common case would be a raster of points,
# for plotting purposes.
# This call can be avoided by typing in.(A, Ref(d)) instead.
broadcasted(::DomainSetStyle, ::typeof(in), A, d::Domain) = broadcast_in(A, d)
broadcasted(::DomainSetStyle, ::typeof(approx_in), A, d::Domain, tol) = broadcast_approx_in(A, d, tol)

"Vectorized version of `in`: apply `x ∈ d` to all elements of `A`."
broadcast_in(A, d::Domain) = in.(A, Ref(d))

"Vectorized version of `approx_in`: apply `x ∈ d` to all elements of `A`."
broadcast_approx_in(A, d::Domain, tol) = approx_in.(A, Ref(d), tol)

## Some arithmetics

# Allow unary minus, but use broadcast for the implementation
-(d::Domain) = (-).(d)

# Allow multiplication and division by numbers, like for vectors
*(a::Number, domain::Domain) = a .* domain
*(domain::Domain, a::Number) = domain .* a
/(domain::Domain, a::Number) = domain ./ a
\(a::Number, domain::Domain) = a .\ domain


text:
################################################################################
### Test constructors

cs     = ColumnSchema(:customer_id, "Customer ID", Int, !CATEGORICAL, REQUIRED, UNIQUE, 1:1_000_000)
ts     = TableSchema(:mytable, "My table", [cs], [:customer_id])
schema = Schema(:myschema, "My data set", Dict(:mytable => ts))

@test_throws ErrorException ColumnSchema(:customer_id, "Customer ID", Int, !CATEGORICAL, REQUIRED, UNIQUE, UInt)     # Int != UInt
@test_throws ErrorException ColumnSchema(:new, "Customer is new", Char, CATEGORICAL, REQUIRED, !UNIQUE, ["y", "n"])  # Char != String
@test_throws ErrorException TableSchema(:mytable, "My table", [cs], [:XXXcustomer_id])  # primarykey non-existent

################################################################################
### Compare DataFrame to Schema

# Schema
patientid = ColumnSchema(:patientid, "Patient ID",  UInt,   !CATEGORICAL, REQUIRED,  UNIQUE, UInt)
age       = ColumnSchema(:age,       "Age (years)", Int,    !CATEGORICAL, REQUIRED, !UNIQUE, Int)
dose      = ColumnSchema(:dose,      "Dose size",   String,  CATEGORICAL, REQUIRED, !UNIQUE, ["small", "medium", "large"])
fever     = ColumnSchema(:fever,     "Had fever",   Bool,    CATEGORICAL, REQUIRED, !UNIQUE, Bool)
ts        = TableSchema(:mytable, "My table", [patientid, age, dose, fever], [:patientid])
schema    = Schema(:fever, "Fever schema", Dict(:mytable => ts))

pid2 = ColumnSchema(:pid2, "Patient ID - version 2", UInt, !CATEGORICAL, !REQUIRED, UNIQUE, UInt)
@test_throws ErrorException TableSchema(:mytable, "My table", [pid2, age, dose, fever], [:pid2])  # Primary key not unique

# DataFrame
tbl = DataFrame(
    patientid = [1, 2, 3, 4],
    age       = [11, 22, 33, 444],
    dose      = ["small", "medium", "large", "medium"],
    fever     = [false, true, true, false]
)

# Compare data to schema
outdata, issues_in, issues_out = compare(ts, tbl)
@test size(issues_in, 1) == 4

# Modify data to comply with the schema
categorical!(tbl, [:dose, :fever])  # Ensure :dose and :fever contain categorical data
outdata, issues_in, issues_out = compare(ts, tbl)
@test size(issues_in,  1) == 2
@test size(issues_out, 1) == 0

tbl[!, :patientid] = convert(Vector{UInt}, tbl[!, :patientid])
outdata, issues_in, issues_out = compare(ts, tbl)
@test size(issues_in, 1) == 0

# Modify schema: Forbid tbl[:age] having values of 120 or above
schema.tables[:mytable].colname2colschema[:age].validvalues = 0:120

# Compare again
outdata, issues_in, issues_out = compare(schema.tables[:mytable], tbl)
@test size(issues_in, 1)  == 1
@test size(issues_out, 1) == 1
@test ismissing(outdata[4, :age]) == true

# Fix data
tbl[4, :age] = 44
outdata, issues_in, issues_out = compare(schema.tables[:mytable], tbl)
@test size(issues_in, 1) == 0

# Add a new column to the schema
zipcode = ColumnSchema(:zipcode, "Zip code", Int, !CATEGORICAL, !REQUIRED, !UNIQUE, 10000:99999)
schema.tables[:mytable].colname2colschema[:zipcode] = zipcode
push!(schema.tables[:mytable].columnorder, :zipcode)
@test schema.tables[:mytable].columnorder[end] == :zipcode
@test haskey(schema.tables[:mytable].colname2colschema, :zipcode)
@test schema.tables[:mytable].colname2colschema[:zipcode] == zipcode

# Write the updated schema to disk
#schemafile = joinpath(dirname(pathof(Schemata)), "..", "test/schemata/fever_updated.yaml")
#writeschema(schemafile, schema)
#schema_from_disk = readschema(schemafile)
#@test schema == schema_from_disk

# Add a corresponding (non-compliant) column to the data
tbl[!, :zipcode] = ["11111", "22222", "33333", "NULL"];  # CSV file was supplied with "NULL" values, forcing eltype to be String.
outdata, issues_in, issues_out = compare(schema.tables[:mytable], tbl)
@test size(issues_in, 1)  == 2
@test size(issues_out, 1) == 0

# Add a new column to the schema
dosedate = ColumnSchema(:date, "Dose date", Date, !CATEGORICAL, !REQUIRED, !UNIQUE, Date)
schema.tables[:mytable].colname2colschema[:date] = dosedate
push!(schema.tables[:mytable].columnorder, :date)

# Add a corresponding (compliant) column to the data
tbl[!, :date] = ["2017-12-01", "2017-12-01", "2017-12-11", "2017-12-09"];
outdata, issues_in, issues_out = compare(schema.tables[:mytable], tbl)
@test size(issues_in, 1)  == 4
@test size(issues_out, 1) == 0

################################################################################
# Test CustomParser

# Define custom parser
function my_zdt_custom_parser(s::T, tz::String) where {T <: AbstractString}
    occursin(':', s) && return ZonedDateTime(DateTime(s[1:16]), TimeZone(tz))  # Example: s="2020-12-31T09:30:59+10:00"
    dt = Date(eval(Meta.parse(s)))  # Examples: s="today()", s="2020-11-01"
    ZonedDateTime(DateTime(dt), TimeZone(tz))
end

my_zdt_custom_parser(dttm::DateTime, tz::String) = ZonedDateTime(dttm, TimeZone(tz))

# Dict for ColumnSchema constructor, obtained after reading yaml
d = Dict("name"          => "zdt", "description" => "Test custom parser for TimeZones.ZonedDateTime",
         "datatype"      => "ZonedDateTime",
         "iscategorical" => false, "isrequired" => true, "isunique" => true,
         "validvalues"   => "(today()-Year(2), Hour(1), today()-Day(1))",  # Ensure that the range has sufficient resolution
         "parser"        => Dict("function" => "my_zdt_custom_parser", "args"=>["Australia/Melbourne"]))

# Need to eval datatype and parser.function in the same scope that they were defined (and before constructing the ColumnSchema).
# Schemata.jl can't see the datatype and parser.function until it receives them from the current scope.
d["datatype"] = eval(Meta.parse(d["datatype"]))
d["parser"]["function"] = eval(Meta.parse(d["parser"]["function"]))

# Now the schema constructors can be used
cs = ColumnSchema(d)
ts = TableSchema(:mytable, "My table", [cs], [:zdt])

tbl = DataFrame(zdt=[DateTime(today() - Day(7)) + Hour(i) for i = 1:3])
target = [ZonedDateTime(tbl[i, :zdt], TimeZone("Australia/Melbourne")) for i = 1:3]
outdata, issues_in, issues_out = compare(ts, tbl);
@test outdata[!, :zdt] == target

tbl = DataFrame(zdt=[string(DateTime(today() - Day(7)) + Hour(i)) for i = 1:3])  # String type
outdata, issues_in, issues_out = compare(ts, tbl);
@test outdata[!, :zdt] == target

tbl = DataFrame(zdt=[string(ZonedDateTime(DateTime(today() - Day(7)) + Hour(i), TimeZone("Australia/Melbourne"))) for i = 1:3])  # String type
outdata, issues_in, issues_out = compare(ts, tbl);
@test outdata[!, :zdt] == target

################################################################################
# Test conversion of String arg to DateFormat when the datatype is Date.

d = Dict("name"        => "dt",
         "description" => "Test Date parser with a prebuilt DateFormat stored in args, rather than construct a DateFormat for each parsing",
         "datatype"    => "Date", "iscategorical" => false, "isrequired" => true, "isunique" => true,
         "validvalues" => "(01JAN2019, Day(1), 31DEC2019)",
         "parser"      => Dict("function" => "Date", "args"=>["dduuuyyyy"]))

cs  = ColumnSchema(d)
ts  = TableSchema(:mytable, "My table", [cs], [:dt])
tbl = DataFrame(dt=["02DEC2019", "03DEC2019", "04DEC2019"])
outdata, issues_in, issues_out = compare(ts, tbl);
@test outdata[!, :dt] == [Date(2019, 12, 1) + Day(i) for i = 1:3]

################################################################################
# Test intra-row constraints
function test_row_constraints()
    filename = joinpath(dirname(pathof(Schemata)), "..", "test/schemata/row_constraints.yaml")
    schema   = readschema(filename)
    d = DataFrame(
                  patientid = UInt.([1,2,3]),
                  dob=Date.(["1992-10-01", "1988-03-23", "1983-11-18"]),
                  date_of_marriage=[Date("2015-09-13"), missing, Date("1981-11-01")]
                 )
    compare(schema.tables[:dates], d)
end
outdata, issues_in, issues_out = test_row_constraints()
@test size(issues_in, 1) == 1


text:
using IntervalArithmetic
using Base.Test

@testset "Operations on boxes" begin
    A = IntervalBox(1..2, 3..4)
    B = IntervalBox(0..2, 3..6)

    @test 2*A == IntervalBox(2..4, 6..8)
    @test A + B == IntervalBox(1..4, 6..10)
    @test dot(A, B) == @interval(9, 28)

    @test A ⊆ B
    @test A ∩ B == A
    @test A ∪ B == B

    X = IntervalBox(1..2, 3..4)
    Y = IntervalBox(3..4, 3..4)

    @test isempty(X ∩ Y)
    @test X ∪ Y == IntervalBox(1..4, 3..4)

    X = IntervalBox(2..4, 3..5)
    Y = IntervalBox(3..5, 4..17)
    @test X ∩ Y == IntervalBox(3..4, 4..5)

    v = [@interval(i, i+1) for i in 1:10]
    V = IntervalBox(v...)
    @test length(V) == 10

    Y = IntervalBox(1..2)  # single interval
    @test isa(Y, IntervalBox)
    @test length(Y) == 1
    @test Y == IntervalBox( (Interval(1., 2.),) )

end

# @testset "@intervalbox tests" begin
#     @intervalbox f(x, y) = (x + y, x - y)
#
#     X = IntervalBox(1..1, 2..2)
#     @test f(X) == IntervalBox(3..3, -1 .. -1)
#
#     @intervalbox g(x, y) = x - y
#     @test isa(g(X), IntervalBox)
#
#     @test emptyinterval(X) == IntervalBox(∅, ∅)
#
# end

@testset "setdiff for IntervalBox" begin
    X = IntervalBox(2..4, 3..5)
    Y = IntervalBox(3..5, 4..6)
    @test setdiff(X, Y) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(2..3, 3..5) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(2..5, 3..4),
                              IntervalBox(2..5, 5..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(4..6, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(4..5, 3..4),
                              IntervalBox(4..5, 5..6),
                              IntervalBox(2..4, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(3..4, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(3..4, 5..6),
                              IntervalBox(2..3, 3..6),
                              IntervalBox(4..5, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(2..4, 10..20)
    @test setdiff(X, Y) == typeof(X)[X]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, -10..10)
    @test setdiff(X, Y) == typeof(X)[]

    X = IntervalBox(1..4, 3..6, 7..10)
    Y = IntervalBox(2..3, 4..5, 8..9)
    @test setdiff(X, Y) == [ IntervalBox(2..3, 4..5, 7..8),
                              IntervalBox(2..3, 4..5, 9..10),
                              IntervalBox(2..3, 3..4, 7..10),
                              IntervalBox(2..3, 5..6, 7..10),
                              IntervalBox(1..2, 3..6, 7..10),
                              IntervalBox(3..4, 3..6, 7..10) ]

end

@testset "mid, diam, × for IntervalBox" begin
    X = (0..2) × (3..5)
    @test length(X) == 2
    @test X == IntervalBox(Interval(0, 2), Interval(3, 5))

    @test diam(X) == 2
    @test mid(X) == [1, 4]

    Y = X × (4..8)
    @test isa(Y, IntervalBox)
    @test length(Y) == 3
    @test Y == IntervalBox(Interval(0, 2), Interval(3, 5), Interval(4, 8))
    @test diam(Y) == 4
end

@testset "Constructing multidimensional IntervalBoxes" begin
    @test IntervalBox(1..2, Val{1}) == IntervalBox(1..2)
    @test IntervalBox(1..2, Val{2}) == (1..2) × (1..2)
    @test IntervalBox(1..2, Val{5}) == (1..2) × (1..2) × (1..2) × (1..2) × (1..2)

    @test IntervalBox(1..2, 3) == IntervalBox(1..2, Val{3})

end


text:
using Grid
using PyPlot

f(x) = 2 .* cos(6x) .+ sin(14x) .+ 2.5
c_grid = 0:.2:1

Af = CoordInterpGrid(c_grid, f(c_grid), BCnil, InterpLinear)

f_grid = linspace(0, 1, 150)

fig = figure()
ax = fig[:add_subplot](111)

ax[:plot](f_grid, f(f_grid), "b-", lw=2, alpha=0.8, label="true function")
ax[:plot](f_grid, Af[f_grid], "g-", lw=2, alpha=0.8,
          label="linear approximation")

ax[:vlines]([c_grid], [c_grid] * 0, f(c_grid), linestyle="dashed", alpha=0.5)
ax[:legend](loc="upper center")


text:
#= Periodic Table
Chemical constants defining properties of atoms on the periodic table.

Only extends till element 108 (Hassium) due to experimental data limitations.

Data taken from PubChem, US NIH.

National Center for Biotechnology Information. PubChem Database. Periodic Table of Elements, https://pubchem.ncbi.nlm.nih.gov/periodic-table/ (accessed on Mar. 4, 2020)

=#

const periodic_table = Dict{String, Int}("H" => 1, "He" => 2, "Li" => 3, "Be" => 4, "B" => 5,
                                         "C" => 6, "N" => 7, "O" => 8, "F" => 9, "Ne" => 10,
                                         "Na" => 11, "Mg" => 12, "Al" => 13, "Si" => 14, "P" => 15,
                                         "S" => 16, "Cl" => 17, "Ar" => 18, "K" => 19, "Ca" => 20,
                                         "Sc" => 21, "Ti" => 22, "V" => 23, "Cr" => 24, "Mn" => 25,
                                         "Fe" => 26, "Co" => 27, "Ni" => 28, "Cu" => 29, "Zn" => 30,
                                         "Ga" => 31, "Ge" => 32, "As" => 33, "Se" => 34, "Br" => 35,
                                         "Kr" => 36, "Rb" => 37, "Sr" => 38, "Y" => 39, "Zr" => 40,
                                         "Nb" => 41, "Mo" => 42, "Tc" => 43, "Ru" => 44, "Rh" => 45,
                                         "Pd" => 46, "Ag" => 47, "Cd" => 48, "In" => 49, "Sn" => 50,
                                         "Sb" => 51, "Te" => 52, "I" => 53, "Xe" => 54, "Cs" => 55,
                                         "Ba" => 56, "La" => 57, "Ce" => 58, "Pr" => 59, "Nd" => 60,
                                         "Pm" => 61, "Sm" => 62, "Eu" => 63, "Gd" => 64, "Tb" => 65,
                                         "Dy" => 66, "Ho" => 67, "Er" => 68, "Tm" => 69, "Yb" => 70,
                                         "Lu" => 71, "Hf" => 72, "Ta" => 73, "W" => 74, "Re" => 75,
                                         "Os" => 76, "Ir" => 77, "Pt" => 78, "Au" => 79, "Hg" => 80,
                                         "Tl" => 81, "Pb" => 82, "Bi" => 83, "Po" => 84, "At" => 85,
                                         "Rn" => 86, "Fr" => 87, "Ra" => 88, "Ac" => 89, "Th" => 90,
                                         "Pa" => 91, "U" => 92, "Np" => 93, "Pu" => 94, "Am" => 95,
                                         "Cm" => 96, "Bk" => 97, "Cf" => 98, "Es" => 99, "Fm" => 100,
                                         "Md" => 101, "No" => 102, "Lr" => 103, "Rf" => 104, "Db" => 105,
                                         "Sg" => 106, "Bh" => 107, "Hs" => 108)

const periodic_table_reverse = Dict{Int, String}(values(periodic_table) .=> keys(periodic_table))

const periodic_table_mass = [1.0080, 4.00260, 7.0, 9.012183, 10.81, 12.011, 14.007, 15.999, 18.99840316,
                             20.180, 22.9897693, 24.305, 26.981538, 28.085, 30.97376200, 32.07, 35.45,
                             39.9, 39.098, 40.08, 44.95591, 47.87, 50.941, 51.996, 54.93804, 55.84,
                             58.93319, 58.693, 63.55, 65.4, 69.72, 72.63, 74.92159, 78.97, 79.90, 83.80,
                             85.468, 87.6, 88.9058, 91.22, 92.9064, 96.0, 97.90721, 101.1, 102.9055, 106.4,
                             107.868, 112.41, 114.82, 118.71, 121.76, 127.6, 126.9045, 131.29, 132.9054520,
                             137.33, 138.9055, 140.12, 140.9077, 144.24, 144.91276, 150.4, 151.96, 157.2,
                             158.92535, 162.50, 164.93033, 167.26, 168.93422, 173.04, 174.967, 178.5,
                             180.9479, 183.8, 186.21, 190.2, 192.22, 195.08, 196.96657, 200.59, 204.383,
                             207, 208.9804, 208.98243, 209.98715, 222.01758, 223.01973, 226.02541,
                             227.02775, 232.038, 231.0359, 238.0289, 237.04817, 244.06420, 243.06138,
                             247.07035, 247.07031, 251.07959, 252.0830, 257.09511, 258.09843, 259.10100,
                             262.110, 267.122, 268.126, 271.134, 274.144, 277.152]

const periodic_table_electrons = ["1s1", "1s2", "[He]2s1", "[He]2s2", "[He]2s2 2p1", "[He]2s2 2p2",
                                  "[He]   2s2 2p3", "[He]2s2 2p4", "[He]2s2 2p5", "[He]2s2 2p6", "[Ne]3s1",
                                  "[Ne]3s2", "[Ne]3s2 3p1", "[Ne]3s2 3p2", "[Ne]3s2 3p3", "[Ne]3s2 3p4",
                                  "[Ne]3s2 3p5", "[Ne]3s2 3p6", "[Ar]4s1", "[Ar]4s2", "[Ar]4s2 3d1",
                                  "[Ar]4s2 3d2", "[Ar]4s2 3d3", "[Ar]3d5 4s1", "[Ar]4s2 3d5",
                                  "[Ar]4s2 3d6", "[Ar]4s2 3d7", "[Ar]4s2 3d8", "[Ar]4s1 3d10",
                                  "[Ar]4s2 3d10", "[Ar]4s2 3d10 4p1", "[Ar]4s2 3d10 4p2",
                                  "[Ar]4s2 3d10 4p3", "[Ar]4s2 3d10 4p4", "[Ar]4s2 3d10 4p5",
                                  "[Ar]4s2 3d10 4p6", "[Kr]5s1", "[Kr]5s2"," [Kr]5s2 4d1", "[Kr]5s2 4d2",
                                  "[Kr]5s1 4d4", "[Kr]5s1 4d5", "[Kr]5s2 4d5", "[Kr]5s1 4d7",
                                  "[Kr]5s1 4d8", "[Kr]4d10", "[Kr]5s1 4d10"," [Kr]5s2 4d10",
                                  "[Kr]5s2 4d10 5p1"," [Kr]5s2 4d10 5p2", "[Kr]5s2 4d10 5p3",
                                  "[Kr]5s2 4d10 5p4", "[Kr]5s2 4d10 5p5", "[Kr]5s2 4d10 5p6", "[Xe]6s1",
                                  "[Xe]6s2", "[Xe]6s2 5d1", "[Xe]6s2 4f1 5d1"," [Xe]6s2 4f3",
                                  "[Xe]6s2 4f4", "[Xe]6s2 4f5", "[Xe]6s2 4f6", "[Xe]6s2 4f7",
                                  "[Xe]6s2 4f7 5d1", "[Xe]6s2 4f9", "[Xe]6s2 4f10", "[Xe]6s2 4f11",
                                  "[Xe]6s2 4f12", "[Xe]6s2 4f13", "[Xe]6s2 4f14", "[Xe]6s2 4f14 5d1",
                                  "[Xe]6s2 4f14 5d2", "[Xe]6s2 4f14 5d3", "[Xe]6s2 4f14 5d4",
                                  "[Xe]6s2 4f14 5d5", "[Xe]6s2 4f14 5d6", "[Xe]6s2 4f14 5d7",
                                  "[Xe]6s1 4f14 5d9", "[Xe]6s1 4f14 5d10", "[Xe]6s2 4f14 5d10",
                                  "[Xe]6s2 4f14 5d10 6p1", "[Xe]6s2 4f14 5d10 6p2",
                                  "[Xe]6s2 4f14 5d10 6p3", "[Xe]6s2 4f14 5d10 6p4",
                                  "[Xe]6s2 4f14 5d10 6p5", "[Xe]6s2 4f14 5d10 6p6", "[Rn]7s1",
                                  "[Rn]7s2", "[Rn]7s2 6d1", "[Rn]7s2 6d2", "[Rn]7s2 5f2 6d1",
                                  "[Rn]7s2 5f3 6d1", "[Rn]7s2 5f4 6d1", "[Rn]7s2 5f6", "[Rn]7s2 5f7",
                                  "[Rn]7s2 5f7 6d1", "[Rn]7s2 5f9", "[Rn]7s2 5f10", "[Rn]7s2 5f11",
                                  "[Rn] 5f12 7s2", "[Rn]7s2 5f13", "[Rn]7s2 5f14", "[Rn]7s2 5f14 6d1",
                                  "[Rn]7s2 5f14 6d2", "[Rn]7s2 5f14 6d3", "[Rn]7s2 5f14 6d4",
                                  "[Rn]7s2 5f14 6d5", "[Rn]7s2 5f14 6d6"]

const periodic_table_electronegativity = [2.2, missing, 0.98, 1.57, 2.04, 2.55, 3.04, 3.44, 3.98, missing,
                                          0.93, 1.31, 1.61, 1.9, 2.19, 2.58, 3.16, missing, 0.82, 1.0,
                                          1.36, 1.54, 1.63, 1.66, 1.55, 1.83, 1.88, 1.91, 1.9, 1.65, 1.81,
                                          2.01, 2.18, 2.55, 2.96, 3.0, 0.82, 0.95, 1.22, 1.33, 1.6, 2.16,
                                          1.9, 2.2, 2.28, 2.2, 1.93, 1.69, 1.78, 1.96, 2.05, 2.1, 2.66,
                                          2.6, 0.79, 0.89, 1.1, 1.12, 1.13, 1.14, missing, 1.17, missing,
                                          1.2, missing, 1.22, 1.23, 1.24, 1.25, missing, 1.27, 1.3, 1.5,
                                          2.36, 1.9, 2.2, 2.2, 2.28, 2.54, 2.0, 1.62, 2.33, 2.02, 2.0, 2.2,
                                          missing, 0.7, 0.9, 1.1, 1.3, 1.5, 1.38, 1.36, 1.28, 1.3, 1.3,
                                          1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, missing, missing, missing,
                                          missing, missing]