text:
#=
Source:
  Ntaimo, L. and S. Sen, "The 'million-variable' march for stochastic combinatorial optimization," Journal of Optimization, 2005.

Input:
  nJ: number of servers
  nI: number of clients
  nS: number of scenarios

Sets:
  sI: clients
  sJ: servers
  sZ: zones

Variables (1st Stage):
  x[j]: 1 if a server is located at site j, 0 otherwise

Variables (2nd Stage):
  y[i,j]: 1 if client i is served by a server at location j, 0 otherwise
  y0[j]: any overflows that are not served due to limitations in server capacity

Parameters (general):
  c[j]: cost of locating a server at location j
  q[i,j]: revenue from client i being served by server at location j
  q0[j]: overflow penalty
  d[i,j]: client i resource demand from server at location j
  u: server capacity
  v: an upper bound on the total number of servers that can be located
  w[z]: minimum number of servers to be located in zone z
  Jz[z]: the subset of server locations that belong to zone z
  p[s]: probability of occurrence for scenario s

Parameters (scenario):
  h[i,s]: 1 if client i is present in scenario s, 0 otherwise
=#
using DualDecomposition
using JuMP, Ipopt, GLPK
using Random

const DD = DualDecomposition

function main_sslp(nJ::Int, nI::Int, nS::Int, seed::Int=1)
    Random.seed!(seed)

    sJ = 1:nJ
    sI = 1:nI
    sS = 1:nS

    c = rand(40:80,nJ)
    q = rand(0:25,nI,nJ,nS)
    q0 = ones(nJ)*1000
    d = q
    u = 1.5*sum(d)/nJ
    v = nJ
    h = rand(0:1,nI,nS)
    Pr = ones(nS)/nS

    # This creates a Lagrange dual problem for each scenario s.
    function create_scenario_model(s::Int64)
        model = Model(GLPK.Optimizer)

        @variable(model, x[j=sJ], Bin)
        @variable(model, y[i=sI,j=sJ], Bin)
        @variable(model, y0[j=sJ] >= 0)

        @objective(model, Min,
              sum(c[j]*x[j] for j in sJ)
            - sum(q[i,j,s]*y[i,j] for i in sI for j in sJ)
            + sum(q0[j]*y0[j] for j in sJ))

        @constraint(model, sum(x[j] for j in sJ) <= v)
        @constraint(model, [j=sJ], sum(d[i,j,s]*y[i,j] for i in sI) - y0[j] <= u*x[j])
        @constraint(model, [i=sI], sum(y[i,j] for j in sJ) == h[i,s])

        return model
    end

    # Create DualDecomposition instance.
    algo = DD.LagrangeDual()

    # Add Lagrange dual problem for each scenario s.
    models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in sS)
    for s in sS
        DD.add_block_model!(algo, s, models[s])
    end

    coupling_variables = Vector{DD.CouplingVariableRef}()
    for s in sS
        model = models[s]
        xref = model[:x]
        for i in sJ
            push!(coupling_variables, DD.CouplingVariableRef(s, i, xref[i]))
        end
    end

    # Set nonanticipativity variables as an array of symbols.
    DD.set_coupling_variables!(algo, coupling_variables)

    # Solve the problem with the solver; this solver is for the underlying bundle method.
    DD.run!(algo, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
end

main_sslp(10,50,50)
# main_sslp(10,50,100)
# main_sslp(10,50,500)
# main_sslp(10,50,1000)
# main_sslp(10,50,2000)


text:
using Phylogenetics
using Base.Test

include("test/operations.jl")
include("test/io.jl")

ext: jl


text:
# ~/zjul/corona/wrld/jmh/src/jcor.jl  - Grand Loop Version from 006
# 200424 007
#
#
using CSV, DataFrames, Dates, Printf, Plots

# Specify which countries to extract and give population in millions:
# (from: https://www.worldometers.info/world-population/population-by-country/)
clist  = ["Argentina" "Belgium" "Brazil" "Norway" "Panama" "Spain" "Sweden" "Thailand"]
plist  = [45.195774,  11.589623, 212.559417, 5.421241,4.314767, 46.754778, 10.099265, 67.799978]

# Select which data source to read: cfm, dth or rcv
data = "dth"    #"dth" # "rcv" #  "dth" "cfm"

# Grand Loop:
for data in ["cfm" "dth" "rcv"]
    println("running jcor: ", data, "...")
    # end  # at EOF

    if length(ARGS) >= 1 data=ARGS[1] end  # Command line selection
    if data == "cfm"
        df_dta_src = CSV.read("../data/confirmed.csv")
        th = 100; ths="100";
        plotstring="Confirmed Cases";
    elseif data == "rcv"
        df_dta_src = CSV.read("../data/recovered.csv")
        th = 30; ths="30";  # Error in data for Norway end on 32
        plotstring="Recovered";
    elseif  data == "dth"
        df_dta_src = CSV.read("../data/deaths.csv")
        th = 10; ths="10";
        plotstring="Deaths";
    else
        println("Incorrect data selection: should be one of: cfm dth rcv")
        #@goto EOF  # readline()
    end

    countries = df_dta_src[:,2]
    lastdate = String(names(df_dta_src)[end])

    # Extracting the epidemic data:
    cindx = Array{Int64}(undef, length(clist))
    dta   = zeros(size(df_dta_src)[2]-4, length(clist)); # Data starts in column 4
    for i in 1:length(clist)
        indx = findall(x->x==clist[i], df_dta_src[:,2])
        dta[:,i] = convert(Array, df_dta_src[indx,5:end])
    end
    df_dta = DataFrame(dta)
    for i in 1:length(clist)
         rename!(df_dta, names(df_dta)[i] => Symbol(clist[i]));
    end

    # Converting the Dates:
    colstring = Array{Union{Nothing, String}}(nothing, length(names(df_dta_src))-4,2)
    coldates  = Array{Date,length(names(df_dta_src))-4}
    for i in 1:length(names(df_dta_src))-4 colstring[i,1] = String(names(df_dta_src)[i+4])*"20" end #  adding 20 to get 2020
    colstring2=colstring
    for i=1:length(colstring[:,1]) colstring2[i,1]=replace(colstring[i,1], "/" => "-") end
    coldates=Date.(colstring2[:,1], Dates.DateFormat("mm-dd-yyyy"))
    df_00 = DataFrame([coldates])
    rename!(df_00, :x1 => :Date)
    CSV.write("dfs/df_00.csv", df_00)

    # Joining Dates and Data:
    df_dta = hcat(df_00, df_dta)
    CSV.write("dfs/df_"*data*".csv", df_dta)

    # Finding first occurrence of more than 'th' confirmed cases
    firsts = zeros(Int, size(df_dta)[2]-1)
    for j in 1:size(df_dta)[2]-1 firsts[j] = findfirst(x->x>th, dta[:,j]) end

    # Extracting above threshold data:
    mtx = Array{Union{Missing, Float64}}(missing,size(df_dta)[1]-minimum(firsts)+1, size(df_dta)[2]-1)
    for j in 1:size(mtx)[2]
        for i in 1:size(dta)[1]-firsts[j]+1   mtx[i,j] = dta[firsts[j]-1+i,j] end
    end
    df_dta_th = DataFrame(mtx)
    for i in 1:length(clist)
         rename!(df_dta_th, names(df_dta_th)[i] => Symbol(clist[i]));
    end
    CSV.write("dfs/df_"*data*"_"*ths*".csv", df_dta_th)

    # Per Million inhabitants:
    mtx2 = copy(mtx)
    for j in 1:length(clist) mtx2[:,j] = mtx2[:,j]./plist[j] end
    df_dta_thm = DataFrame(mtx2)
    for i in 1:length(clist)
        rename!(df_dta_thm, names(df_dta_thm)[i] => Symbol(clist[i]));
    end
    CSV.write("dfs/df_"*data*"_"*ths*"m.csv", df_dta_thm)

    # PLOTTING:
    updated = @sprintf "%10s" df_dta.Date[end]

    plot(df_dta[:,1], dta,
         label = clist,
         legend=:outertopright, xrotation=45,
         xlabel = "End of data: "*updated, xguidefontsize=8,
         title = "Corona 2020\n(1)"*plotstring*" - Selected Countries")
    savefig("figs/"*data*"_01.png")

    labels1 = String[]
    for i in 2:length(names(df_dta)) push!(labels1, String(names(df_dta)[i])) end
    labels = reshape(labels1,1,:)
    nmlast = [findfirst( x-> ismissing(x), mtx[:,i]) for i in 1:size(mtx)[2]]
    for i in 1:length(nmlast)
        if nmlast[i] == nothing
            nmlast[i] = size(mtx,1)
        else
            nmlast[i]= nmlast[i] - 1
        end
    end

    plot(mtx, yaxis=:log,
         label = clist,
         legend=:outertopright, xrotation=45,
         ylabel = "Logaritmic Scale", xguidefontsize=8,
         xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated,
         annotations=([(nmlast[i], mtx[nmlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(nmlast)]),
         title  = "Corona 2020\n(2)"*plotstring*" - Selected Countries")
    savefig("figs/"*data*"_02.png")

    plot(mtx2,
         label = clist,
         legend=:outertopright, xrotation=45,
         xlabel="Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
         annotations=([(nmlast[i], mtx2[nmlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(nmlast)]),
         title = "Corona 2020\n(3)"*plotstring*" per Million Population\nSelected Countries")
    savefig("figs/"*data*"_03.png")

end # Grand Loop

# @label EOF
# EOF

max_stars_repo_name: dyrehaugen/jdt

max_issues_repo_name: dyrehaugen/jdt

max_forks_repo_name: dyrehaugen/jdt


text:
##
## Join / merge
##

function join_idx(left, right, max_groups)
    ## adapted from Wes McKinney's full_outer_join in pandas (file: src/join.pyx).

    # NA group in location 0

    left_sorter, where, left_count = groupsort_indexer(left, max_groups)
    right_sorter, where, right_count = groupsort_indexer(right, max_groups)

    # First pass, determine size of result set
    tcount = 0
    rcount = 0
    lcount = 0
    for i in 1:(max_groups + 1)
        lc = left_count[i]
        rc = right_count[i]

        if rc > 0 && lc > 0
            tcount += lc * rc
        elseif rc > 0
            rcount += rc
        else
            lcount += lc
        end
    end

    # group 0 is the NA group
    tposition = 0
    lposition = 0
    rposition = 0

    left_pos = 0
    right_pos = 0

    left_indexer = Array(Int, tcount)
    right_indexer = Array(Int, tcount)
    leftonly_indexer = Array(Int, lcount)
    rightonly_indexer = Array(Int, rcount)
    for i in 1:(max_groups + 1)
        lc = left_count[i]
        rc = right_count[i]
        if rc == 0
            for j in 1:lc
                leftonly_indexer[lposition + j] = left_pos + j
            end
            lposition += lc
        elseif lc == 0
            for j in 1:rc
                rightonly_indexer[rposition + j] = right_pos + j
            end
            rposition += rc
        else
            for j in 1:lc
                offset = tposition + (j-1) * rc
                for k in 1:rc
                    left_indexer[offset + k] = left_pos + j
                    right_indexer[offset + k] = right_pos + k
                end
            end
            tposition += lc * rc
        end
        left_pos += lc
        right_pos += rc
    end

    ## (left_sorter, left_indexer, leftonly_indexer,
    ##  right_sorter, right_indexer, rightonly_indexer)
    (left_sorter[left_indexer], left_sorter[leftonly_indexer],
     right_sorter[right_indexer], right_sorter[rightonly_indexer])
end

function PooledDataVecs(df1::AbstractDataFrame,
                        df2::AbstractDataFrame)
    # This method exists to allow merge to work with multiple columns.
    # It takes the columns of each DataFrame and returns a DataArray
    # with a merged pool that "keys" the combination of column values.
    # The pools of the result don't really mean anything.
    dv1, dv2 = PooledDataVecs(df1[1], df2[1])
    refs1 = dv1.refs + 1   # the + 1 handles NA's
    refs2 = dv2.refs + 1
    ngroups = length(dv1.pool) + 1
    for j = 2:ncol(df1)
        dv1, dv2 = PooledDataVecs(df1[j], df2[j])
        for i = 1:length(refs1)
            refs1[i] += (dv1.refs[i]) * ngroups
        end
        for i = 1:length(refs2)
            refs2[i] += (dv2.refs[i]) * ngroups
        end
        ngroups = ngroups * (length(dv1.pool) + 1)
    end
    pool = [1:ngroups]
    (PooledDataArray(refs1, pool), PooledDataArray(refs2, pool))
end

function PooledDataArray(df::AbstractDataFrame)
    # This method exists to allow another way for merge to work with
    # multiple columns. It takes the columns of the DataFrame and
    # returns a DataArray with a merged pool that "keys" the
    # combination of column values.
    # Notes:
    #   - I skipped the sort to make it faster.
    #   - Converting each individual one-row DataFrame to a Tuple
    #     might be faster.
    refs = zeros(POOLED_DATA_VEC_REF_TYPE, nrow(df))
    poolref = Dict{AbstractDataFrame, Int}()
    pool = Array(Uint64, 0)
    j = 1
    for i = 1:nrow(df)
        val = df[i,:]
        if has(poolref, val)
            refs[i] = poolref[val]
        else
            push!(pool, hash(val))
            refs[i] = j
            poolref[val] = j
            j += 1
        end
    end
    return PooledDataArray(refs, pool)
end

function merge(df1::AbstractDataFrame, df2::AbstractDataFrame, bycol, jointype)

    dv1, dv2 = PooledDataVecs(df1[bycol], df2[bycol])
    left_indexer, leftonly_indexer, right_indexer, rightonly_indexer =
        join_idx(dv1.refs, dv2.refs, length(dv1.pool))

    if jointype == "inner"
        return cbind(df1[left_indexer,:], without(df2, bycol)[right_indexer,:])
    elseif jointype == "left"
        left = df1[[left_indexer,leftonly_indexer],:]
        right = rbind(without(df2, bycol)[right_indexer,:],
                      nas(without(df2, bycol), length(leftonly_indexer)))
        return cbind(left, right)
    elseif jointype == "right"
        left = rbind(without(df1, bycol)[left_indexer, :],
                     nas(without(df1, bycol), length(rightonly_indexer)))
        right = df2[[right_indexer,rightonly_indexer],:]
        return cbind(left, right)
    elseif jointype == "outer"
        mixed = cbind(df1[left_indexer, :], without(df2, bycol)[right_indexer, :])
        leftonly = cbind(df1[leftonly_indexer, :],
                         nas(without(df2, bycol), length(leftonly_indexer)))
        leftonly = leftonly[:, colnames(mixed)]
        rightonly = cbind(nas(without(df1, bycol), length(rightonly_indexer)),
                          df2[rightonly_indexer, :])
        rightonly = rightonly[:, colnames(mixed)]
        return rbind(mixed, leftonly, rightonly)
    end
end

merge(df1::AbstractDataFrame, df2::AbstractDataFrame, bycol) = merge(df1, df2, bycol, "inner")

function merge(df1::AbstractDataFrame, df2::AbstractDataFrame)
    s1 = Set{ByteString}()
    for coln in colnames(df1)
        add!(s1, coln)
    end
    s2 = Set{ByteString}()
    for coln in colnames(df2)
        add!(s2, coln)
    end
    bycol = first(collect(intersect(s1, s2)))
    merge(df1, df2, bycol, "inner")
end


text:
@testset "push-pull" begin
    Signals.async_mode(false)
    A = Signal(1)
    B = Signal(x -> x+1, A)

    C = Signal(x -> x+1, A)
    D = Signal(x -> x+1, C)
    E = Signal((x,y) -> x+y, B, D)
    @test E[] == 5
    A(2)
    @test E[] == 7
    A[] = 3
    @test E[] == 7      # still unaffected
    @test E() == 9      # properly pulled
    A(1.0)
    @test E[] === 5.0   # dynamic change of type
end

@testset "strict_push" begin
    Signals.async_mode(true)
    A_strict = Signal(1; strict_push = true)
    B_strict = Signal(A_strict; state = 0) do a, state
        state.x += a
        a + 1
    end

    A_soft = Signal(1)
    B_soft = Signal(A_soft; state = 0) do a, state
        state.x += a
        a + 1
    end

    A_soft(0); sleep(1) # restart the eventlopp

    @test B_soft.state.x == 1
    A_soft(10); A_soft(100);
    yield()
    @test B_soft.state.x == 101

    @test B_strict.state.x == 1
    A_strict(10); A_strict(100);
    yield()
    @test B_strict.state.x == 111
end


text:
function read_bus_ids(line)
    buses = Int64[]
    for id in split(line, ",")
        if id != "x"
            push!(buses, parse(Int64, id))
        end
    end
    buses
end

function waiting_time(departure, buses)
    min_wt = typemax(Int64)
    min_bus_id = nothing
    for bus_id in buses
        remainder = departure % bus_id
        if remainder == 0
            return 0
        end
        wt = bus_id - remainder
        if wt < min_wt
            min_wt = wt
            min_bus_id = bus_id
        end
    end
    min_wt * min_bus_id
end

function main()
    lines = readlines("../resources/input_13.txt")
    departure = parse(Int64, lines[1])
    buses = read_bus_ids(lines[2])
    result = waiting_time(departure, buses)
    println(result)
end

main()


text:
module HyperParameters

import Base: convert, eltype, show, ==, *, /, +, -, ^, besselk, exp, gamma, tanh

export Bound, Interval, leftbounded, rightbounded, unbounded, checkbounds, Variable, fixed,
    Argument, HyperParameter

#== Bound Type ==#

immutable Bound{T<:Real}
    value::T
    isopen::Bool
end
Bound{T<:Real}(value::T, isopen::Bool) = Bound{T}(value, isopen)

function Bound{T<:Real}(value::T, boundtype::Symbol)
    if boundtype == :open
        return Bound(value, true)
    elseif boundtype == :closed
        return Bound(value, false)
    else
        error("Bound type $boundtype not recognized")
    end
end

eltype{T<:Real}(::Bound{T}) = T
convert{T<:Real}(::Type{Bound{T}}, B::Bound) = Bound(convert(T, B.value), B.isopen)

#== Interval Type ==#

immutable Interval{T<:Real}
    left::Nullable{Bound{T}}
    right::Nullable{Bound{T}}
    function Interval(left::Nullable{Bound{T}}, right::Nullable{Bound{T}})
        if !isnull(left) && !isnull(right)
            rbound = get(right)
            lbound = get(left)
            if lbound.isopen || rbound.isopen
                lbound.value <  rbound.value || error("Invalid bounds")
            else
                lbound.value <= rbound.value || error("Invalid bounds")
            end
        end
        new(left, right)
    end
end
Interval{T<:Real}(left::Nullable{Bound{T}}, right::Nullable{Bound{T}}) = Interval{T}(left, right)
Interval{T<:Real}(left::Bound{T}, right::Bound{T}) = Interval(Nullable(left), Nullable(right))

eltype{T}(::Interval{T}) = T

function leftbounded{T<:Real}(value::T, boundtype::Symbol)
    if boundtype == :open
        Interval(Nullable(Bound(value, true)),  Nullable{Bound{T}}())
    elseif boundtype == :closed
        Interval(Nullable(Bound(value, false)), Nullable{Bound{T}}())
    else
        error("Bound type $boundtype not recognized")
    end
end
leftbounded{T<:Real}(left::Bound{T}) = Interval(Nullable(left), Nullable{Bound{T}}())

function rightbounded{T<:Real}(value::T, boundtype::Symbol)
    if boundtype == :open
        Interval(Nullable{Bound{T}}(), Nullable(Bound(value, true)))
    elseif boundtype == :closed
        Interval(Nullable{Bound{T}}(), Nullable(Bound(value, false)))
    else
        error("Bound type $boundtype not recognized")
    end
end
rightbounded{T<:Real}(right::Bound{T}) = Interval(Nullable{Bound{T}}(), Nullable(right))

unbounded{T<:Real}(::Type{T}) = Interval(Nullable{Bound{T}}(), Nullable{Bound{T}}())

function convert{T<:Real}(::Type{Interval{T}}, I::Interval)
    if isnull(I.left)
        isnull(I.right) ? unbounded(T) : rightbounded(convert(Bound{T}, get(I.right)))
    else
        if isnull(I.right)
            leftbounded(convert(Bound{T}, get(I.left)))
        else
            Interval(convert(Bound{T}, get(I.left)), convert(Bound{T}, get(I.right)))
        end
    end
end

function description_string{T}(I::Interval{T})
    interval =  string("Interval{", T, "}")
    if isnull(I.left)
        if isnull(I.right)
            string(interval, "(-∞,∞)")
        else
            string(interval, "(-∞,", get(I.right).value, get(I.right).isopen ? ")" : "]")
        end
    else
        left = string(get(I.left).isopen ? "(" : "[",  get(I.left).value, ",")
        if isnull(I.right)
            string(interval, left, "∞)")
        else
            string(interval, left, get(I.right).value, get(I.right).isopen ? ")" : "]")
        end
    end
end
function show{T}(io::IO, I::Interval{T})
    print(io, description_string(I))
end

function checkbounds{T<:Real}(I::Interval{T}, x::T)
    if isnull(I.left)
        if isnull(I.right)
            true
        else
            ub = get(I.right)
            ub.isopen ? (x < ub.value) : (x <= ub.value)
        end
    else
        lb = get(I.left)
        if isnull(I.right)
            lb.isopen ? (lb.value < x) : (lb.value <= x)
        else
            ub = get(I.right)
            if ub.isopen
                lb.isopen ? (lb.value < x < ub.value) : (lb.value <= x < ub.value)
            else
                lb.isopen ? (lb.value < x <= ub.value) : (lb.value <= x <= ub.value)
            end
        end
    end
end

#== Variable Type ==#

immutable Variable{T<:Real}
    value::T
    isfixed::Bool
end
Variable{T<:Real}(value::T, isfixed::Bool=false) = Variable{T}(value, false)
Variable{T<:Real}(value::Variable{T}) = value
fixed{T<:Real}(v::T) = Variable{T}(v, true)

eltype{T<:Real}(::Variable{T}) = T

convert{T<:Real}(::Type{Variable{T}}, var::Variable) = Variable(convert(T, var.value), var.isfixed)

typealias Argument{T<:Real} Union{T,Variable{T}}

#== HyperParameter Type ==#

type HyperParameter{T<:Real}
    value::T
    bounds::Interval{T}
    isfixed::Bool
    function HyperParameter(value::T, bounds::Interval{T}, isfixed::Bool)
        checkbounds(bounds, value) || error("Value $(value) must be in range " * string(bounds))
        new(value, bounds, isfixed)
    end
end
function HyperParameter{T<:Real}(x::T, bounds::Interval{T} = NullBound(T), isfixed::Bool=false)
    HyperParameter{T}(x, bounds, isfixed)
end
function HyperParameter{T<:Real}(x::Variable{T}, bounds::Interval{T} = NullBound(T))
    HyperParameter(x.value, bounds, x.isfixed)
end

function convert{T<:Real}(::Type{HyperParameter{T}}, θ::HyperParameter)
    HyperParameter{T}(convert(T, θ.value), convert(Interval{T}, θ.bounds), θ.isfixed)
end

function show{T}(io::IO, θ::HyperParameter{T})
    print(io, "HyperParameter{" * string(T) * "}(", θ.value, ") ∈ ", θ.bounds)
end

@inline *(a::Real, v::HyperParameter) = *(a, v.value)
@inline *(v::HyperParameter, a::Real) = *(v.value, a)

@inline /(a::Real, v::HyperParameter) = /(a, v.value)
@inline /(v::HyperParameter, a::Real) = /(v.value, a)

@inline +(a::Real, v::HyperParameter) = +(a, v.value)
@inline +(v::HyperParameter, a::Real) = +(v.value, a)

@inline -(v::HyperParameter) = -(v.value)
@inline -(a::Real, v::HyperParameter) = -(a, v.value)
@inline -(v::HyperParameter, a::Real) = -(v.value, a)

@inline ^(a::Real, v::HyperParameter)          = ^(a, v.value)
@inline ^(v::HyperParameter, a::Integer)       = ^(v.value, a)
@inline ^(v::HyperParameter, a::AbstractFloat) = ^(v.value, a)

@inline besselk(v::HyperParameter, x::Real) = besselk(v.value, x)
@inline exp(v::HyperParameter)  = exp(v.value)
@inline gamma(v::HyperParameter) = gamma(v.value)
@inline tanh(v::HyperParameter) = tanh(v.value)

@inline ==(a::Real, v::HyperParameter) = ==(a, v.value)
@inline ==(v::HyperParameter, a::Real) = ==(v.value, a)
@inline ==(v1::HyperParameter, v2::HyperParameter) = ==(v1.value, v2.value)

end # End HyperParameter


text:
"""
    mutable struct UnitCell

    Cell(dim, crystal, params)

A unit cell representing a real-space and inverse-space basis in 1, 2 or 3 dimensions.

A specified crystal system controls the geometry of the cell. These include:

*1D: (:lamellar)
*2D: (:square, :rectangular, :hexagonal, :oblique)
*3D: (:cubic, :tetragonal, :orthorhombic, :hexagonal, :monoclinic, :triclinic)
"""
mutable struct UnitCell
    dim     :: Int
    crystal :: Symbol
    volume  :: Float64
    params  :: Vector{Float64}

    R       :: SMat3D{Float64}  # R[:,i] = Bravais lattice basis vector a_i
    G       :: SMat3D{Float64}  # G[:,i] = Reciprocal lattice basis vector b_i
    dR      :: Array{Float64,3} # Derivatives of a_i w.r.t. cell parameters
    dG      :: Array{Float64,3} # Derivatives of b_i
    dRR     :: Array{Float64,3} # Derivatives of a_i ⋅ a_j
    dGG     :: Array{Float64,3} # Derivatives of b_i ⋅ b_j

    # Grids for k-vectors
    ksq     :: FieldGrid{Float64}
    dksq    :: Vector{FieldGrid{Float64}}

    function UnitCell(dim::Integer, crystal::Symbol, params::AbstractVector{<:Real})
        cell = new()
        cell.dim = dim
        cell.crystal = crystal
        cell.params = params

        npr = nparams(cell)
        cell.dR = zeros(3, 3, npr)
        cell.dG = zeros(3, 3, npr)
        cell.dRR = zeros(3, 3, npr)
        cell.dGG = zeros(3, 3, npr)

        # Add the k-grids
        cell.ksq = FieldGrid{Float64}(undef, 0, 0, 0)
        cell.dksq = []

        # Update the basis information to make current before return
        update!(cell)

        return cell
    end
end

UnitCell(dim::Integer, crystal::Symbol, params...) = UnitCell(dim, crystal, [params...])

#==============================================================================#

Base.show(io::IO, cell::UnitCell) = @printf(io, "Cell(dim = %d, V = %.2f)", cell.dim, cell.volume)

nparams(cell::UnitCell) = length(cell.params)

function paramstring(cell::UnitCell)
    str_params = map(x -> @sprintf("%.3e",x), cell.params)
    return "[" * join(str_params, ", ") * "]"
end

"""
    setup!(cell, sys)

Allocate k-grids of an appropriate size for the system.
"""
function setup!(cell::UnitCell, sys::AbstractSystem)
    dim = ndims(sys.dims)
    if dim != cell.dim
        error("Number of non-singleton dimensions in System and Cell do not match!")
    end

    # Setup the cell
    cell.ksq = zeros(Float64, ksize(sys.dims))
    cell.dksq = [zeros(Float64, ksize(sys.dims)) for k = 1:nparams(cell)]

    ksq!(cell)
    dksq!(cell)

    return nothing
end

"""
    update!(cell)

Update the cell basis matrices and wavevector grids
given the current set of parameters (lengths/angles).

Must be called after every update to cell parameters during a variable cell calculation.
"""
function update!(cell::UnitCell)
    # Construct a new basis matrix from cell params
    basis!(cell)

    # Calculate cell basis derivatives: dR, dG, dRR, dGG
    dbasis!(cell)
    ddbasis!(cell)

    # Update the ksq and dksq grids
    ksq!(cell)
    dksq!(cell)

    return nothing
end

#==============================================================================#

function basis!(cell::UnitCell)
    cell.R = make_basis(cell.dim, cell.crystal, cell.params)
    cell.G = 2*pi*transpose(inv(cell.R))
    cell.volume = det(cell.R)
    return nothing
end

function dbasis!(cell::UnitCell)
    # Calculates the variation of the real and reciprocal basis matrices
    #   w.r.t. each of the (6) cell parameters, a/b/c/alpha/beta/gamma
    # Analytically:
    #   dR_ijk = d(R_ij) / d(theta_k)
    #   dG_ijk = d(G_ij) / d(theta_k)
    fill!(cell.dR, 0.0)
    fill!(cell.dG, 0.0)

    # Calculate the derivatives of the real basis matrix w.r.t each parameter
    #   R == 3x3 matrix
    #   grad(R) = 3x3x6 array, gradient of cell dimensions adds a rank
    if cell.dim == 1
        if cell.crystal == :lamellar
            cell.dR[1,1,1] = 1.0
        end
    elseif cell.dim == 2
        if cell.crystal == :square
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,1] = 1.0
        elseif cell.crystal == :rectangular
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,2] = 1.0
        elseif cell.crystal == :hexagonal
            cell.dR[1,1,1] = 1.0
            cell.dR[1,2,1] = 0.5
            cell.dR[2,2,1] = sqrt(0.75)
        elseif cell.crystal == :oblique
            a, b, gamma = cell.params
            cell.dR[1,1,1] = 1.0
            cell.dR[1,2,2] = cos(gamma)
            cell.dR[2,2,2] = sin(gamma)
            cell.dR[1,2,3] = -b * sin(gamma)
            cell.dR[2,2,3] = b * cos(gamma)
        end
    elseif cell.dim == 3
        if cell.crystal == :cubic
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,1] = 1.0
            cell.dR[3,3,1] = 1.0
        elseif cell.crystal == :tetragonal
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,1] = 1.0
            cell.dR[3,3,2] = 1.0
        elseif cell.crystal == :orthorhombic
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,2] = 1.0
            cell.dR[3,3,3] = 1.0
        elseif cell.crystal == :hexagonal
            cell.dR[1,1,1] = 1.0
            cell.dR[1,2,1] = 0.5
            cell.dR[2,2,1] = sqrt(0.75)
            cell.dR[3,3,2] = 1.0
        elseif cell.crystal == :monoclinic
            a, b, c, beta = cell.params
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,2] = 1.0
            cell.dR[1,3,3] = cos(beta)
            cell.dR[3,3,3] = sin(beta)
            cell.dR[1,3,4] = -c*sin(beta)
            cell.dR[3,3,4] = c*cos(beta)
        elseif cell.crystal == :triclinic
            a, b, c, alpha, beta, gamma = cell.params
            cx, cy, cz = cell.R[1,3], cell.R[2,3], cell.R[3,3]
            cell.dR[1,1,1] = 1.0
            cell.dR[1,2,2] = cos(gamma)
            cell.dR[2,2,2] = sin(gamma)
            cell.dR[1,3,3] = cos(beta)
            cell.dR[2,3,3] = cy / c
            cell.dR[3,3,3] = cz / c
            # The dirty ones w/ angle derivatives, 3rd column gets messy
            cell.dR[2,3,4] = -c * sin(alpha) / sin(gamma)
            cell.dR[3,3,4] = c^2 * (cos(alpha) - cos(beta)*cos(gamma)) * sin(alpha) / sin(gamma)^2 / cz
            cell.dR[1,3,5] = -c * sin(beta)
            cell.dR[2,3,5] = c * sin(beta) / tan(gamma)
            cell.dR[3,3,5] = c^2 * (-cos(alpha) / tan(gamma) + cos(beta) / sin(gamma)) * sin(beta) / sin(gamma) / cz
            cell.dR[1,2,6] = -b * sin(gamma)
            cell.dR[2,2,6] = b * cos(gamma)
            cell.dR[2,3,6] = c * (-cos(alpha) / tan(gamma) + cos(beta) / sin(gamma)) / sin(gamma)
            cell.dR[3,3,6] = c^2 * ((2 + cos(2*alpha) + cos(2*beta)) / tan(gamma) - cos(alpha) * cos(beta)
                 * (3 + cos(2*gamma)) / sin(gamma) )  / sin(gamma)^2 / cz / 2
        elseif cell.crystal == :pnna
            cell.dR[1,1,1] = 2.0
            cell.dR[2,2,1] = sqrt(3.0)
            cell.dR[3,3,1] = 1.0
        elseif cell.crystal == :fddd1
            cell.dR[1,1,1] = 1.0
            cell.dR[2,2,1] = sqrt(3.0)
            cell.dR[3,3,1] = 2.0*sqrt(3.0)
        end
    end

    # Update the dG matrix using the dR matrix calculated above
    #   d(b_u) = -(1/2pi) * sum_v [ b_u * d(a_v) b_v ] (Jian's thesis & 2003 stress paper)
    #   So we can express reciprocal variation in terms of real variation dR matrix
    G = cell.G
    dR = cell.dR
    for k = 1:nparams(cell)
        for i = 1:cell.dim
            for j = 1:cell.dim
                for l = 1:cell.dim
                    for m = 1:cell.dim
                        cell.dG[i,j,k] -= G[i,l] * dR[m,l,k] * G[m,j]
                    end
                end
            end
        end
    end
    cell.dG ./= 2*pi

    return nothing
end

function ddbasis!(cell::UnitCell)
    # Calculate derivatives of basis vector dot products w.r.t. cell params
    # Analytically:
    #   dRR_ijk = d(a_i ⋅ a_j) / d(theta_k)
    #   dGG_ijk = d(b_i ⋅ b_j) / d(theta_k)
    fill!(cell.dRR, 0.0)
    fill!(cell.dGG, 0.0)

    # Updates dRR and dGG based on the current basis matrix
    R = transpose(cell.R) # Transpose works to correct indexing from PSCF
    G = cell.G
    dR = cell.dR
    dG = cell.dG
    for k = 1:nparams(cell)
        for i = 1:cell.dim
            for j = 1:cell.dim
                for l = 1:cell.dim
                    cell.dRR[i,j,k] += R[i,l] * dR[l,j,k] + R[j,l] * dR[l,i,k]
                    cell.dGG[i,j,k] += G[l,i] * dG[l,j,k] + G[l,j] * dG[l,i,k]
                end
            end
        end
    end

    return nothing
end

#==============================================================================#

"""
    ksq!(cell)

Update the |k^2| grid given the current reciprocal cell basis `Gbasis`.
The wavevectors are shifted in accordance with the output of a real FFT.
"""
function ksq!(cell::UnitCell)
    Nx, Ny, Nz = size(cell.ksq)
    G = cell.G

    for idx in CartesianIndices(cell.ksq)
        ix, iy, iz = idx[1]-1, idx[2]-1, idx[3]-1

        jx = ix  # Assumes first dim is half-space
        jy = iy <= Ny/2 ? iy : Ny-iy
        jz = iz <= Nz/2 ? iz : Nz-iz
        jvec = @SVector [jx, jy, jz]
        kvec = G * jvec

        @inbounds cell.ksq[idx] = dot(kvec, kvec)
    end

    return nothing
end

"""
    dksq!(cell)

Update the derivatives of |k^2| for each cell parameter.
"""
function dksq!(cell::UnitCell)
    Nx, Ny, Nz = size(cell.ksq)
    G = cell.G
    dGG = cell.dGG

    #Threads.@threads for k = 1:length(cell.dksq)
    for k = 1:length(cell.dksq)
        dksq = cell.dksq[k]
        fill!(dksq, 0.0)

        for idx in CartesianIndices(dksq)
            ix, iy, iz = idx[1]-1, idx[2]-1, idx[3]-1

            jx = ix  # Assumes first dim is half-space
            jy = iy <= Ny/2 ? iy : Ny-iy
            jz = iz <= Nz/2 ? iz : Nz-iz
            jvec = @SVector [jx, jy, jz]

            for l = 1:cell.dim
                for m = 1:cell.dim
                    @inbounds dksq[idx] += jvec[l] * jvec[m] * dGG[l,m,k]
                end
            end

        end
    end

    return nothing
end


text:
@testset "Misc Integrals" begin
    @reset
    aoints = Fermi.Integrals.IntegralHelper()
    @test_throws Fermi.Options.FermiException Fermi.Integrals.compute!(aoints, "invalid")

    @set {
        basis = cc-pvtz
        jkfit auto
        rifit auto
    }

    jkfit = Fermi.Integrals.JKFIT()
    @test Fermi.string_repr(jkfit) == "JKFIT: cc-pvtz-jkfit"

    rifit = Fermi.Integrals.RIFIT()
    @test Fermi.string_repr(rifit) == "RIFIT: cc-pvtz-rifit"

    @reset
    @set precision "xablau"
    @test_throws Fermi.Options.FermiException Fermi.Integrals.IntegralHelper()

    @set precision "single"
    ints = Fermi.Integrals.IntegralHelper()
    S = ints["S"]
    @test S[1,1] |> typeof === Float32

    ints["T"]
    ints["V"] = rand(5,5)
    @test begin
        "S" in keys(ints.cache) &&
        "T" in keys(ints.cache) &&
        "V" in keys(ints.cache)
    end

    Fermi.Integrals.delete!(ints, "V")
    @test begin
        "S" in keys(ints.cache) &&
        "T" in keys(ints.cache) &&
        !("V" in keys(ints.cache))
    end
    Fermi.Integrals.delete!(ints)

    @test begin
        !("S" in keys(ints.cache)) &&
        !("T" in keys(ints.cache)) &&
        !("V" in keys(ints.cache))
    end

    @reset
    mol = Molecule(molstring="He 0.0 0.0 0.0")
    bset = Fermi.Integrals.BasisSet("sto-3g", mol)
    ints = Fermi.Integrals.IntegralHelper(bset)
    ints["S"]
    @test Fermi.string_repr(ints) == " ⇒ Fermi IntegralHelper\n ⋅ Data Type:                 Float64\n ⋅ Basis:                     sto-3g\n ⋅ ERI:                       SparseERI\n ⋅ Orbitals:                  AtomicOrbitals\n ⋅ Stored Integrals:          S "

end

@testset "1-e AO -> MO" begin
    @reset
    @set {
        printstyle none
        df true
        basis cc-pvdz
        jkfit cc-pvtz-jkfit
    }
    aoints = Fermi.Integrals.IntegralHelper()
    wfn = @energy aoints => rhf
    Cμi = wfn.orbitals.C
    moints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)

    Fermi.Integrals.compute!(moints, aoints, "S")
    @test moints.cache["S"] ≈ Cμi' * aoints["S"] * Cμi

    Fermi.Integrals.compute!(moints, aoints, "T")
    @test moints.cache["T"] ≈ Cμi' * aoints["T"] * Cμi

    Fermi.Integrals.compute!(moints, aoints, "V")
    @test moints.cache["V"] ≈ Cμi' * aoints["V"] * Cμi
end

@testset "DF AO -> MO" begin

    @reset
    @set {
        printstyle none
        df true
        basis cc-pvdz
        jkfit cc-pvtz-jkfit
        drop_occ 2
    }

    @molecule {
        C       -1.1170269915      4.4925791088      0.0000000582
        C        0.3969669792      4.3675789789     -0.0000001928
        H       -1.5729510532      3.5466015111      0.3608481999
        H       -1.4776736551      4.7005983883     -1.0293753400
        H       -1.4262457918      5.3234861671      0.6685275007
        H        0.7061857730      3.5366719134     -0.6685276289
        H        0.7576136432      4.1595597057      1.0293752066
        H        0.8528910462      5.3135565710     -0.3608483431

    }

    wfn = @energy rhf
    Cμi = wfn.orbitals.C
    Co = Cμi[:, 3:9]
    Cv = Cμi[:,10:end]

    moints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals, eri_type=Fermi.Integrals.RIFIT())
    aoints = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.RIFIT())

    # contract B(prs) = C(μr) * B(pμν) * C(νs) from scratch
    function test_contract(C1, C2, B)
        p,μ,ν = size(B)
        r = size(C1, 2)
        s = size(C2, 2)

        Bμpν = permutedims(B, (2,1,3))
        Bμpν = reshape(Bμpν, (μ, p*ν))
        Brpν = C1' * Bμpν
        Brpν = reshape(Brpν, (r*p, ν))
        Brps = Brpν * C2
        Brps = reshape(Brps, (r,p,s))
        return permutedims(Brps, (2,1,3))
    end

    # Contract (ijkl) = B(pij) * B(pkl)
    function test_contract2(B1, B2)
        p,i,j = size(B1)
        _,k,l = size(B2)
        B1 = reshape(B1, p,i*j)
        B2 = reshape(B2, p,k*l)
        out = B1' * B2
        return reshape(out, (i,j,k,l))
    end

    # Full Bpq
    @test moints["ERI"] ≈ test_contract(Cμi, Cμi, aoints["ERI"])

    # Boo
    @test moints["BOO"] ≈ test_contract(Co, Co, aoints["ERI"])

    # Bov
    @test moints["BOV"] ≈ test_contract(Co, Cv, aoints["ERI"])

    # Bvv
    @test moints["BVV"] ≈ test_contract(Cv, Cv, aoints["ERI"])

    Boo = moints["BOO"]
    Bov = moints["BOV"]
    Bvv = moints["BVV"]

    oooo = test_contract2(Boo, Boo)
    ooov = test_contract2(Boo, Bov)
    oovv = test_contract2(Boo, Bvv)
    ovov = test_contract2(Bov, Bov)
    ovvv = test_contract2(Bov, Bvv)
    vvvv = test_contract2(Bvv, Bvv)

    @test moints["OOOO"] ≈ oooo
    @test moints["OOOV"] ≈ ooov
    @test moints["OOVV"] ≈ oovv
    @test moints["OVOV"] ≈ ovov
    @test moints["OVVV"] ≈ ovvv
    @test moints["VVVV"] ≈ vvvv

    # Test again, cleaning the cache such that Bpq arrays need to be compute on the fly
    delete!(moints)

    @test moints["OOOO"] ≈ oooo
    @test moints["OOOV"] ≈ ooov
    @test moints["OOVV"] ≈ oovv
    @test moints["OVOV"] ≈ ovov
    @test moints["OVVV"] ≈ ovvv
    @test moints["VVVV"] ≈ vvvv
end

@testset "Chonky AO -> MO" begin

    @reset
    @set {
        printstyle none
        basis cc-pvdz
        drop_occ 1
    }

    @molecule {
        C        0.3340594014      4.3347983152      0.0216449677
        H        1.4398420378      4.2538788474     -0.0166555421
        H       -0.0088087137      5.0751637968     -0.7300667682
        H        0.0222217713      4.6642047575      1.0340770936
        H       -0.1170174852      3.3459458686     -0.2007749095
    }

    wfn = @energy rhf
    C  = wfn.orbitals.C
    Co = C[:, 2:5]
    Cv = C[:,6:end]

    moints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals, eri_type=Fermi.Integrals.Chonky())
    aoints = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.Chonky())

    # contract pqrs = C1(μp) * C2(νq) * C3(ρr) * C4(σs) * μνρσ
    function test_contract(C1, C2, C3, C4, μνρσ)

        μ,ν,ρ,σ = size(μνρσ)
        _,s = size(C4)
        _,r = size(C3)
        _,q = size(C2)
        _,p = size(C1)

        # μνρ,σ * σs
        μνρs = reshape(reshape(μνρσ, (μ*ν*ρ, σ)) * C4, (μ,ν,ρ,s))

        # μνs,ρ * ρr
        μνsr = reshape(reshape(permutedims(μνρs, (1,2,4,3)), (μ*ν*s, ρ)) * C3, (μ,ν,s,r))

        # μrs,ν * νq
        μrsq = reshape(reshape(permutedims(μνsr, (1,4,3,2)), (μ*r*s, ν)) * C2, (μ,r,s,q))

        # rsq,μ * μp
        pqrs = reshape(C1' * reshape(permutedims(μrsq, (1,4,2,3)), (μ, q*r*s)), (p,q,r,s))

        return pqrs
    end

    # Full ERI
    @test moints["ERI"] ≈ test_contract(C, C, C, C, aoints["ERI"])
    delete!(moints)

    # oooo
    @test moints["OOOO"] ≈ test_contract(Co, Co, Co, Co, aoints["ERI"])

    # ooov
    @test moints["OOOV"] ≈ test_contract(Co, Co, Co, Cv, aoints["ERI"])

    # oovv
    @test moints["OOVV"] ≈ test_contract(Co, Co, Cv, Cv, aoints["ERI"])

    # ovov
    @test moints["OVOV"] ≈ test_contract(Co, Cv, Co, Cv, aoints["ERI"])

    # ovvv
    @test moints["OVVV"] ≈ test_contract(Co, Cv, Cv, Cv, aoints["ERI"])

    # vvvv
    @test moints["VVVV"] ≈ test_contract(Cv, Cv, Cv, Cv, aoints["ERI"])
end

@testset "Sparse AO -> MO" begin

    @reset
    @set {
        printstyle none
        basis cc-pvdz
        drop_occ 1
        drop_vir 3
    }

    @molecule {
        N       -1.0811450850      4.5346506907     -0.6162384680
        H       -0.0380366676      4.5141429777     -0.6171266138
        H       -1.4191680136      3.9897206949      0.2067368203
        H       -1.4191724649      4.0628676478     -1.4832192333
    }

    wfn = @energy rhf

    mo1 = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals, eri_type=Fermi.Integrals.Chonky())
    mo2 = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals, eri_type=Fermi.Integrals.Chonky())
    aosp = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.SparseERI())
    aock = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.Chonky())

    # oooo
    Fermi.Integrals.compute_OOOO!(mo1, aosp)
    Fermi.Integrals.compute_OOOO!(mo2, aock)
    @test mo1["OOOO"] ≈ mo2["OOOO"]

    # ooov
    Fermi.Integrals.compute_OOOV!(mo1, aosp)
    Fermi.Integrals.compute_OOOV!(mo2, aock)
    @test mo1["OOOV"] ≈ mo2["OOOV"]

    # oovv
    Fermi.Integrals.compute_OOVV!(mo1, aosp)
    Fermi.Integrals.compute_OOVV!(mo2, aock)
    @test mo1["OOVV"] ≈ mo2["OOVV"]

    # ovov
    Fermi.Integrals.compute_OVOV!(mo1, aosp)
    Fermi.Integrals.compute_OVOV!(mo2, aock)
    @test mo1["OVOV"] ≈ mo2["OVOV"]

    # ovvv
    Fermi.Integrals.compute_OVVV!(mo1, aosp)
    Fermi.Integrals.compute_OVVV!(mo2, aock)
    @test mo1["OVVV"] ≈ mo2["OVVV"]

    # vvvv
    Fermi.Integrals.compute_VVVV!(mo1, aosp)
    Fermi.Integrals.compute_VVVV!(mo2, aock)
    @test mo1["VVVV"] ≈ mo2["VVVV"]
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using JointMoments
using Base.Test

include("data.jl")
include("precomputed.jl")
include("tolerance.jl")

ε = 1e-4

tests = ["center",
         "replicate",
         "tensors",
         "collapse",
         "statistics"]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
function plot_fw(wf::WaterFlooding)
    fw(sw) = fractional_flow(sw, wf.kr, wf.μw, wf.μo)
    swi = wf.Si
    swj = wf.Sj
    s̃w = wf.S̃
    swr = wf.kr.swr
    sor = wf.kr.sor
    sw = collect(swr:0.01:1 - sor)

    plot(sw, fw(sw), lw=3, label=false)
    plot!(lims=(0, 1), title="Water Fractional Flow")
    plot!([swi s̃w swj], [fw(swi) fw(s̃w) fw(swj)], labels=["Swi" "S̃w" "Swj"], seriestype=:scatter)
    plot!([swi, s̃w], [fw(swi), fw(s̃w)], ls=:dash, label="Shock")
    plot!(xlabel="Water Saturation, Sw", ylabel="Water Fractional Flow, fw")
end

function plot_fw(μw::Float64, μo::Float64, krs::RelPerms...)
    plot()

    for kr in krs
        fw(sw) = fractional_flow(sw, kr, μw, μo)
        sw = collect(kr.swr:0.01:1 - kr.sor)
        display(plot!(sw, fw(sw), lw=3))
    end
    plot!(lims=(0, 1))
    plot!(xlabel="Water Saturation, Sw", ylabel="Water Fractional Flow, fw")
end

function plot_sw_profile(wf::WaterFlooding, time::Float64...)
    plot()
    for t in time
        t_label = round(t, digits=2)
        x, sw = saturation_profile(wf, t)
        display(plot!(x, sw, lw=3.0, fill=(0, 0.1), label="Buckley Leverett", title="Injected PV = $t_label", lims=(0,1)))
    end
    plot!(lims=(0, 1))
    plot!(xlabel="Position, x", ylabel="Water Saturation, Sw")
end

function animate_sw_profile(wf::WaterFlooding)
    dfw(sw) = fw_derivative(sw, wf.kr, wf.μw, wf.μo)
    tmax = 1 / dfw(wf.S̃)
    times = collect(range(0, 3 * tmax, length=300))

    anim = @animate for t in times
        plot_sw_profile(wf::WaterFlooding, t)
    end

    gif(anim, fps=10)
end

# Tracer

function plot_tracer_fw(wf, tracer)
    fw(s) = fractional_flow(s, wf.kr, wf.μw, wf.μo)
    plot_fw(wf)
    plot!([0, tracer.sc], [0, fw(tracer.sc)], marker=:circle, ls=:dot, label="Tracer")
end

function plot_sw_profile(wf::WaterFlooding, c::Tracer, time::Float64)
    plot_sw_profile(wf, time)
    x, s = tracer_profile(wf, c, time)
    xc = c.v * time

    plot!(x, s, lw=0, fill=(0, 0.1, :yellowgreen), label=false)
    plot!([xc, xc], [0, c.sc], lw=2, ls=:dash,  label="Tracer", color=:yellowgreen)

end

function animate_sw_profile(wf::WaterFlooding, c::Tracer)
    dfw(sw) = fw_derivative(sw, wf.kr, wf.μw, wf.μo)
    tmax = 1 / dfw(wf.S̃)
    times = collect(range(0, 3 * tmax, length=300))

    anim = @animate for t in times
        plot_sw_profile(wf, c, t)
    end

    gif(anim, fps=10)
end

# * Polymer Flood

function plot_sw_profile(pf::PolymerFlooding, time::Float64...)
    plot()
    for t in time
        x, s = saturation_profile(pf, t)
        display(plot!(x, s, lw=3.0, fill=(0, 0.1), label="Time $t"))
    end
    plot!(lims=(0, 1))
    plot!(xlabel="Position, x", ylabel="Water Saturation, Sw")
end

function plot_sw_profile(pf::PolymerFlooding, c::Tracer, time::Float64)
    x, s = tracer_profile(pf, c, time)
    xp, sp = polymer_profile(pf, time)
    xc = c.v * time

    plot_sw_profile(pf, time)

    plot!(x, s, lw=0, fill=(0, 0.1, :yellowgreen), label=false)
    plot!([xc, xc], [0, c.sc], lw=3, ls=:dot,  label="Tracer", color=:yellowgreen)

    plot!(xp, sp, lw=0, fill=(0, 0.2, :deeppink2), label=false)
    plot!([xp[1], xp[1]], [0, pf.Sb], lw=3, ls=:dot,  label="Polymer", color=:deeppink2)

end

function animate_sw_profile(pf::PolymerFlooding, c::Tracer)
    dfw(sw) = fw_derivative(sw, pf.kr, pf.μw, pf.μo)
    tmax = 1 / pf.Vb
    times = collect(range(0, 2 * tmax, length=200))

    anim = @animate for t in times
        plot_sw_profile(pf, c, t)
    end

    gif(anim, fps=10)
end


text:
# Autogenerated wrapper script for CUDA_jll for x86_64-w64-mingw32
export compute_sanitizer, libcublas, libcublasLt, libcudadevrt, libcudart, libcufft, libcufftw, libcupti, libcurand, libcusolver, libcusolverMg, libcusparse, libdevice, libnppc, libnppial, libnppicc, libnppidei, libnppif, libnppig, libnppim, libnppist, libnppisu, libnppitc, libnpps, libnvblas, libnvtoolsext, libnvvm, nvdisasm

JLLWrappers.@generate_wrapper_header("CUDA")
JLLWrappers.@declare_executable_product(compute_sanitizer)
JLLWrappers.@declare_library_product(libcublas, "cublas64_11.dll")
JLLWrappers.@declare_library_product(libcublasLt, "cublasLt64_11.dll")
JLLWrappers.@declare_file_product(libcudadevrt)
JLLWrappers.@declare_library_product(libcudart, "cudart64_110.dll")
JLLWrappers.@declare_library_product(libcufft, "cufft64_10.dll")
JLLWrappers.@declare_library_product(libcufftw, "cufftw64_10.dll")
JLLWrappers.@declare_library_product(libcupti, "cupti64_2021.3.0.dll")
JLLWrappers.@declare_library_product(libcurand, "curand64_10.dll")
JLLWrappers.@declare_library_product(libcusolver, "cusolver64_11.dll")
JLLWrappers.@declare_library_product(libcusolverMg, "cusolverMg64_11.dll")
JLLWrappers.@declare_library_product(libcusparse, "cusparse64_11.dll")
JLLWrappers.@declare_file_product(libdevice)
JLLWrappers.@declare_library_product(libnppc, "nppc64_11.dll")
JLLWrappers.@declare_library_product(libnppial, "nppial64_11.dll")
JLLWrappers.@declare_library_product(libnppicc, "nppicc64_11.dll")
JLLWrappers.@declare_library_product(libnppidei, "nppidei64_11.dll")
JLLWrappers.@declare_library_product(libnppif, "nppif64_11.dll")
JLLWrappers.@declare_library_product(libnppig, "nppig64_11.dll")
JLLWrappers.@declare_library_product(libnppim, "nppim64_11.dll")
JLLWrappers.@declare_library_product(libnppist, "nppist64_11.dll")
JLLWrappers.@declare_library_product(libnppisu, "nppisu64_11.dll")
JLLWrappers.@declare_library_product(libnppitc, "nppitc64_11.dll")
JLLWrappers.@declare_library_product(libnpps, "npps64_11.dll")
JLLWrappers.@declare_library_product(libnvblas, "nvblas64_11.dll")
JLLWrappers.@declare_library_product(libnvtoolsext, "nvToolsExt64_1.dll")
JLLWrappers.@declare_library_product(libnvvm, "nvvm64_40_0.dll")
JLLWrappers.@declare_executable_product(nvdisasm)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        compute_sanitizer,
        "bin\\compute-sanitizer.exe",
    )

    JLLWrappers.@init_library_product(
        libcublas,
        "bin\\cublas64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcublasLt,
        "bin\\cublasLt64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libcudadevrt,
        "lib\\cudadevrt.lib",
    )

    JLLWrappers.@init_library_product(
        libcudart,
        "bin\\cudart64_110.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcufft,
        "bin\\cufft64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcufftw,
        "bin\\cufftw64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcupti,
        "bin\\cupti64_2021.3.0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcurand,
        "bin\\curand64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolver,
        "bin\\cusolver64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolverMg,
        "bin\\cusolverMg64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusparse,
        "bin\\cusparse64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libdevice,
        "share\\libdevice\\libdevice.10.bc",
    )

    JLLWrappers.@init_library_product(
        libnppc,
        "bin\\nppc64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppial,
        "bin\\nppial64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppicc,
        "bin\\nppicc64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppidei,
        "bin\\nppidei64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppif,
        "bin\\nppif64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppig,
        "bin\\nppig64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppim,
        "bin\\nppim64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppist,
        "bin\\nppist64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppisu,
        "bin\\nppisu64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnppitc,
        "bin\\nppitc64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnpps,
        "bin\\npps64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvblas,
        "bin\\nvblas64_11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvtoolsext,
        "bin\\nvToolsExt64_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvvm,
        "bin\\nvvm64_40_0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nvdisasm,
        "bin\\nvdisasm.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct GPUOffloadFactorization <: AbstractFactorization end

function SciMLBase.solve(cache::LinearCache, alg::GPUOffloadFactorization; kwargs...)
    if cache.isfresh
        fact = do_factorization(alg, CUDA.CuArray(cache.A), cache.b, cache.u)
        cache = set_cacheval(cache, fact)
    end

    copyto!(cache.u,cache.b)
    y = Array(ldiv!(cache.cacheval, CUDA.CuArray(cache.u)))
    SciMLBase.build_linear_solution(alg,y,nothing,cache)
end

function do_factorization(alg::GPUOffloadFactorization, A, b, u)
    A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("LU is not defined for $(typeof(A))")

    if A isa AbstractDiffEqOperator
        A = A.A
    end
    fact = qr(CUDA.CuArray(A))
    return fact
end

function LinearAlgebra.ldiv!(x::CUDA.CuArray,_qr::CUDA.CUSOLVER.CuQR,b::CUDA.CuArray)
  _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
  x .= vec(_x)
  CUDA.unsafe_free!(_x)
  return x
end
# make `\` work
LinearAlgebra.ldiv!(F::CUDA.CUSOLVER.CuQR, b::CUDA.CuArray) = (x = similar(b); ldiv!(x, F, b); x)

export GPUOffloadFactorization


text:
using DataFrames, Gadfly
set_default_plot_size(6.6inch, 3.3inch)

expand_grid(xv, yv, n::Int) = vcat([[x y z] for x in xv, y in yv, z in 1:n]...)

a = expand_grid(6.0:10, 1.0:4, 3)

D= DataFrame(x= a[:,1], y=a[:,2], z=a[:,3].*a[:,1].*a[:,2], g = string.(floor.(Int, a[:,3])) )
coord = Coord.cartesian(xmin=6, xmax=10, ymin=1, ymax=4)

plot(D, xgroup=:g, x=:x, y=:y, color=:z,
    Geom.subplot_grid(coord, layer(z=:z, Geom.contour(levels=7))),
    Scale.color_continuous(minvalue=0, maxvalue=120)
)


text:
using ..CryptoGroups.Fields: F2PB, F2GNB, tobits
using CryptoUtils: sqrt_mod_prime

function int2octet(x::Integer)

    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end

    return hex2bytes(hex)
end

function int2octet(x::Integer, N::Int)
    k = div(N, 8, RoundUp)

    bytes = int2octet(x)

    pad = UInt8[0 for i in 1:(k - length(bytes))]

    return UInt8[pad..., bytes...]
end

octet2int(x::Vector{UInt8}) = parse(BigInt, bytes2hex(x), base=16)

octet2int(x::String) =  octet2int(hex2bytes(x))

function octet2bits(x::Vector{UInt8})
    bv = BitVector(u << -i % Bool for u in x for i in 7:-1:0)
    return bv
end

octet2bits(x::Vector{UInt8}, N::Int) = octet2bits(x)[end - N + 1:end]

octet2bits(x::String, N::Int) = octet2bits(hex2bytes(x), N)

function bits2uint8(x::BitVector)

    s = UInt8(0)

    for (i, j) in enumerate(reverse(x))

        if j == 1
            s += UInt8(2)^(i - 1)
        end
    end

    return s
end

function bits2octet(_x::BitVector)

    x = copy(_x)

    if mod(length(x), 8) != 0

        padding = BitVector(0 for i in 1:(8 - mod(length(x), 8)))
        prepend!(x, padding)

    end

    # For now assuming that x is in the length of octets

    N = div(length(x), 8, RoundUp)

    b = reshape(x, 8, N)

    bytes = UInt8[bits2uint8(b[:, i]) for i in 1:N]

    return bytes
end

function decompress(x::BigInt, ỹ::Bool, spec::ECP)

    p = modulus(spec)

    α = mod(x^3 + a(spec) * x + b(spec), p)

    β = sqrt_mod_prime(α, p)

    # Rightmost bit of β. Perhaps this works:
    # Alternativelly I could do circshift
    if mod(β, 2) == ỹ
        y = β
    else
        y = p - β
    end

    return y
end

function decompress(x::BitVector, ỹ::Bool, spec::EC2N)
    @warn "Not implemented, returning nothing..."
    return nothing
end

function point(x_octet::Vector{UInt8}, ỹ::Bool, spec::ECP)
    x = octet2int(x_octet)
    y = decompress(x, ỹ, spec)
    return (x, y)
end

function point(x::Vector{UInt8}, ỹ::Bool, spec::EC2N)
    error("Not implmented")
end

# I could also add validation that the point is on curve.
point(x::Vector{UInt8}, y::Vector{UInt8}, spec::ECP) = (octet2int(x), octet2int(y))

point(x::Vector{UInt8}, y::Vector{UInt8}, spec::EC2N) = (octet2bits(x, bitlength(spec)), octet2bits(y, bitlength(spec)))

function point(x_octet::Vector{UInt8}, y_octet::Vector{UInt8}, ỹ::Bool, spec::ECP)

    x, y = point(x_octet, y_octet, spec)

    y′ = decompress(x, ỹ, spec)

    if !(y == y′)
        @warn "Contradictory compression."
    end

    return (x, y)
end

function point(x_octet::Vector{UInt8}, y_octet::Vector{UInt8}, ỹ::Bool, spec::EC2N)

    x, y = point(x_octet, y_octet, spec)

    @warn "Sign bit is ignored."

    return (x, y)
end

function point(po::Vector{UInt8}, spec::Spec)

    pc = po[1]

    if pc == 0
        return zero(spec)
    elseif pc == 2
        ỹ = false
        x = po[2:end]
        return point(x, ỹ, spec)
    elseif pc == 3
        ỹ = true
        x = po[2:end]
        return point(x, ỹ, spec)
    elseif pc == 4 || pc == 6 || pc == 7

        l = div(length(po) - 1, 2)
        x = po[2:(l+1)]
        y = po[(l+2):(2l+1)]

        if pc == 4
            return point(x, y, spec)
        elseif pc == 6
            return point(x, y, false, spec)
        elseif pc == 7
            return point(x, y, true, spec)
        end
    else
        error("Wrong PC: $pc")
    end
end

point(po::String, spec::Spec) = point(hex2bytes(po), spec)

function _compressed_octet(x::Vector{UInt8}, y::Vector{UInt8}, ỹ::Bool)

    if ỹ == false
        return UInt8[2, x...]
    elseif ỹ == true
        return UInt8[3, x...]
    end

end

function _hybrid_octet(x::Vector{UInt8}, y::Vector{UInt8}, ỹ::Bool)

    if ỹ == false
        return UInt8[6, x..., y...]
    elseif ỹ == true
        return UInt8[7, x..., y...]
    end

end

function _uncompressed_octet(x::Vector{UInt8}, y::Vector{UInt8})
     return UInt8[4, x..., y...]
end

function octet(x::BigInt, y::BigInt, N::Int; mode::Symbol = :uncompressed)

    _x = int2octet(x, N)
    _y = int2octet(y, N)

    if mode == :uncompressed

        return _uncompressed_octet(_x, _y)

    elseif mode in [:compressed, :hybrid]

        ỹ = mod(y, 2) % Bool

        if mode == :compressed

            return _compressed_octet(_x, _y, ỹ)

        elseif mode == :hybrid

            return _hybrid_octet(_x, _y, ỹ)

        end

    else
        error("Unrecognized mode $mode")
    end

end

octet(x::BigInt, y::BigInt, spec::ECP; mode::Symbol = :uncompressed) = octet(x, y, bitlength(modulus(spec)); mode)

_field(basis::PB) = F2PB(basis.f)
_field(basis::GNB) = F2GNB{basis.m, basis.T}

function octet(x::BitVector, y::BitVector, basis::BinaryBasis; mode::Symbol = :uncompressed)

    _x = bits2octet(x)
    _y = bits2octet(y)

    if mode == :uncompressed

        return _uncompressed_octet(_x, _y)

    elseif mode in [:compressed, :hybrid]

        F = _field(basis)

        # Reverse order is just a guess

        #x_ = F(x)
        #y_ = F(y)

        x_ = F(reverse(x))
        y_ = F(reverse(y))

        @warn "Calculation of ỹ could be wrong due to insufficient tests."

        z = y_ * inv(x_)

        ỹ = tobits(z)[end]

        if mode == :compressed

            return _compressed_octet(_x, _y, ỹ)

        elseif mode == :hybrid

            return _hybrid_octet(_x, _y, ỹ)

        end

    else
        error("Unrecognized mode $mode")
    end
end

octet(x::BitVector, y::BitVector, spec::EC2N; mode::Symbol = :uncompressed) = octet(x, y, spec.basis; mode)


text:
# MIT License
#
# Copyright © 2020 Michael J Simms. All rights reserved.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

module Peaks

import Statistics
export find_peaks

# Defines a point. X values are integers. Y values are floating point.
mutable struct GraphPoint
    x::UInt
    y::Float64
end

# List of points.
#struct GraphLine <: Array{GraphPoint} end

# Defines a peak. A peak is described by three points: a left trough, a peak, and a right trough.
mutable struct GraphPeak
    left_trough::GraphPoint
    peak::GraphPoint
    right_trough::GraphPoint
    area::Float64
end

# List of points.
#struct GraphPeakList <: Array{GraphPeak} end

function compute_area(data, current_peak::GraphPeak)
    current_peak.area = 0.0

    if current_peak.left_trough.x < current_peak.right_trough.x
        for i in current_peak.left_trough.x+1:current_peak.right_trough.x
            b = data[i] + data[i - 1]
            current_peak.area += (0.5 * b)
        end
    end
end

# Returns a list of the peaks that rise above the threshold.

function find_peaks_over_threshold(data::Array{Float64}, threshold = 0.0)
    peaks = []
    current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)
    x = 1

    for y in data

        if y < threshold

            # Have we found a peak? If so, add it and start looking for the next one.
            if current_peak.right_trough.x > 0

                # Still descending
                if y <= current_peak.right_trough.y
                    current_peak.right_trough.x = x
                    current_peak.right_trough.y = y

                # Rising
                else
                    compute_area(data, current_peak)
                    push!(peaks, current_peak)
                    current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)
                end

            # Are we looking for a left trough?
            elseif current_peak.left_trough.x == 0
                current_peak.left_trough.x = x
                current_peak.left_trough.y = y

            # If we have a left trough and an existing peak, assume this is the right trough - for now.
            elseif (current_peak.peak.x > current_peak.left_trough.x) && (current_peak.left_trough.x > 0)
                current_peak.right_trough.x = x
                current_peak.right_trough.y = y

            else
                current_peak.left_trough.x = x
                current_peak.left_trough.y = y
            end

        elseif current_peak.left_trough.x > 0 # Left trough is set.

            # Are we looking for a peak or is this bigger than the current peak?
            if current_peak.peak.x == 0 || y >= current_peak.peak.y
                current_peak.peak.x = x
                current_peak.peak.y = y
                current_peak.right_trough.x = 0
                current_peak.right_trough.y = 0.0
            end

        elseif current_peak.right_trough.x > 0 # Right trough is set.
            compute_area(data, current_peak)
            push!(peaks, current_peak)
            current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)

        else # Nothing is set, but the value is above the threshold.
            current_peak.left_trough.x = x
            current_peak.left_trough.y = y
        end

        x = x + 1
    end

    peaks
end

# Returns a list of all statistically significant peaks in the given waveform.
# These are defined as peaks that rise more than one standard deviation above the mean for at least three points on the x axis.

function find_peaks_over_stddev(data::Array{Float64}, sigmas = 1.0)
    mean = Statistics.mean(data)
    stddev = sigmas * Statistics.std(data)
    threshold = mean + stddev
    peaks = find_peaks_over_threshold(data, threshold)
    peaks
end

# Utility functions

function average(data::Array{GraphPoint})
    sum = 0.0

    for item in data
        sum = sum + item.y
    end

    avg = sum / length(data)
    avg
end

function variance(data::Array{GraphPoint}, mean)
    numerator = 0.0

    for item in data
        numerator = numerator + ((item.y - mean) * (item.y - mean))
    end

    variance = numerator / length(data)
    variance
end

function standard_deviation(data::Array{GraphPoint}, mean)
    var = variance(data, mean)
    dev = sqrt(var)
    dev
end

# Returns a list of the peaks that rise above the threshold.

function find_peaks_over_threshold(data::Array{GraphPoint}, threshold = 0.0)
    peaks = []
    current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)
    x = 1

    for point in data

        if point.y < threshold

            # Have we found a peak? If so, add it and start looking for the next one.
            if current_peak.right_trough.x > 0

                # Still descending
                if y <= current_peak.right_trough.y
                    current_peak.right_trough.x = x
                    current_peak.right_trough.y = y

                # Rising
                else
                    compute_area(data, current_peak)
                    push!(peaks, current_peak)
                    current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)
                end

            # Are we looking for a left trough?
            elseif current_peak.left_trough.x == 0
                current_peak.left_trough = point

            # If we have a left trough and an existing peak, assume this is the right trough - for now.
            elseif (current_peak.peak.x > current_peak.left_trough.x) && (current_peak.left_trough.x > 0)
                current_peak.right_trough = point

            else
                current_peak.left_trough = point
            end

        elseif current_peak.left_trough.x > 0 # Left trough is set.

            # Are we looking for a peak or is this bigger than the current peak?
            if current_peak.peak.x == 0 || y >= current_peak.peak.y
                current_peak.peak = point
                current_peak.right_trough.x = 0
                current_peak.right_trough.y = 0.0
            end

        elseif current_peak.right_trough.x > 0 # Right trough is set.
            compute_area(data, current_peak)
            push!(peaks, current_peak)
            current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)

        else # Nothing is set, but the value is above the threshold.
            current_peak.left_trough = point
        end

        x = x + 1
    end

    peaks
end

# Returns a list of all statistically significant peaks in the given waveform.
# These are defined as peaks that rise more than one standard deviation above the mean for at least three points on the x axis.

function find_peaks_over_stddev(data::Array{GraphPoint}, sigmas = 1.0)
    mean = Statistics.mean(data)
    stddev = sigmas * standard_deviation(data, mean)
    threshold = mean + stddev
    peaks = find_peaks_over_threshold(data, threshold)
    peaks
end

# Returns a list of all statistically significant peaks in the given waveform, with the specified minimum area.

function find_peaks_of_size(data::Array{GraphPoint}, minpeakarea, sigmas = 1.0)
    peaks = []

    mean = average(data)
    stddev = sigmas * standard_deviation(data, mean)
    threshold = mean + stddev
    current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)
    x = 1

    for point in data

        if point.y < threshold

            # Have we found a peak? If so, add it and start looking for the next one.
            if current_peak.right_trough.x > 0
                compute_area(data, current_peak)
                if current_peak.area >= minpeakarea
                    push!(peaks, current_peak)
                end
                current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)

            # Are we looking for a left trough?
            elseif current_peak.left_trough.x == 0
                current_peak.left_trough = point

            # If we have a left trough and an existing peak, assume this is the right trough - for now.
            elseif (current_peak.peak.x > current_peak.left_trough.x) && (current_peak.left_trough.x > 0)
                current_peak.right_trough = point

            else
                current_peak.left_trough = point
            end

        elseif current_peak.left_trough.x > 0 # Left trough is set.

            # Are we looking for a peak or is this bigger than the current peak?
            if current_peak.peak.x == 0 || y >= current_peak.peak.y
                current_peak.peak = point
                current_peak.right_trough.x = 0
                current_peak.right_trough.y = 0.0
            end

        elseif current_peak.right_trough.x > 0 # Right trough is set.
            compute_area(data, current_peak)
            if current_peak.area >= minpeakarea
                push!(peaks, current_peak)
            end
            current_peak = GraphPeak(GraphPoint(0,0.0), GraphPoint(0,0.0), GraphPoint(0,0.0), 0.0)

        else # Nothing is set, but the value is above the threshold.
            current_peak.left_trough = point
        end

        x = x + 1
    end

    peaks
end

end

max_stars_repo_name: msimms/PeakFinder

max_issues_repo_name: msimms/PeakFinder

max_forks_repo_name: msimms/PeakFinder


text:
################################################################################
#
#    NfAbs/Hilbert.jl: Hilbert symbols
#
# This file is part of Hecke.
#
# Copyright (c) 2015-2019: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2019 Markus Kirschmer
#
################################################################################

export quadratic_defect, hilbert_symbol

function quadratic_defect(a::Union{Rational{<:Integer},IntegerUnion,fmpq}, p::IntegerUnion)
  return quadratic_defect(fmpq(a), fmpz(p))
end

@doc doc"""
    quadratic_defect(a::Union{NumFieldElem,Rational,fmpq}, p) -> Union{Inf, PosInf}

Returns the valuation of the quadratic defect of the element $a$ at $p$, which
can either be prime object or an infinite place of the parent of $a$.
"""
function quadratic_defect(a, p) end

function quadratic_defect(a::fmpq, p::fmpz)
  if iszero(a)
    return inf
  end

  v = valuation(a, p)

  if isodd(v)
    return v
  end

  a = a//QQ(p)^v

  if isodd(p)
    return jacobi_symbol(mod(a, p), p) == 1 ? inf : v
  end

  a = mod(a, 8)

  if a == 1
    return inf
  elseif a == 5
    return v + 2
  end

  return v + 1
end

function _quadratic_defect_unit(a::NumFieldElem, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  @assert valuation(a, p) == 0 && isdyadic(p)
  o = order(p)
  f = nf(o)
  parent(a) != f && error("incompatible arguments")
  k, h = ResidueField(o, p)
  hex = extend(h, f)
  ok, s = issquare_with_sqrt(hex(a))

  a = a * (hex\(s))^-2;
  w = isone(a) ? inf : valuation(a-1, p)
  ee = 2*valuation(o(2), p);
  pi = f(uniformizer(p))

  while w < ee && iseven(w)
    ok, s = issquare_with_sqrt(hex((a-1) * pi^-w))
    a =  a * (1+(hex\(s))*pi^(div(w,2)))^-2;
    w = isone(a) ? inf : valuation(a-1, p)
  end

  if w > ee return inf, one(f)
  elseif w == ee
    kx, x = PolynomialRing(k, cached = false)
    d = x^2 + x + hex((a-1)//4)
    if !isirreducible(d)
      return inf, one(f)
    end
  end
  return w, a
end

function quadratic_defect(a::NumFieldElem, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  if iszero(a)
    return inf
  end

  v = valuation(a, p)

  if isodd(v)
    return v
  end

  o = order(p)
  f = nf(o)
  pi = f(uniformizer(p))
  a = a//pi^v

  if !isdyadic(p)
    k, h = ResidueField( o, p )
    hex = extend(h, f)
    ok, s = issquare_with_sqrt(hex(a))
    return ok ? inf : v
  end

  w, aa = _quadratic_defect_unit(a, p)
  return w + v
end

@doc doc"""
    quadratic_defect(a::NfOrdElem, p) -> Union{Inf, PosInf}

Returns the valuation of the quadratic defect of the element $a$ at $p$, which
can either be prime object or an infinite place of the parent of $a$.
"""
function quadratic_defect(a::Union{NfRelOrdElem, NfAbsOrdElem}, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  return quadratic_defect(elem_in_nf(a), p)
end

################################################################################
#
#  Hilbert symbol
#
################################################################################

function hilbert_symbol(a::NumFieldElem, b::IntegerUnion, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  return hilbert_symbol(a, parent(a)(b), p)
end

function hilbert_symbol(a::IntegerUnion, b::NumFieldElem, p::Union{NfAbsOrdIdl, NfRelOrdIdl})
  return hilbert_symbol(b, a, p)
end

@doc Markdown.doc"""
    hilbert_symbol(a::NumFieldElem, b::NumFieldElem, p::NfOrdIdl) -> Int

Returns the local Hilbert symbol $(a,b)_p$.
"""
function hilbert_symbol(a::T, b::T, p::Union{NfAbsOrdIdl, NfRelOrdIdl}) where {T <: NumFieldElem}
  (iszero(a) || iszero(b)) && error("arguments must be non-zero")
  o = order(p)
  f = nf(o)
  f !== parent(a) && error("Incompatible parents")
  pi = f(uniformizer(p))
  v = valuation(a,p)
  w = valuation(b,p)

  if isodd(v)
    if iseven(w)
      a,b,v,w = b,a,w,v
    else
      a = -a*b
      v = v+w
    end
  end

  # now v is even, make a a local unit
  a = a // pi^v

  if !isdyadic(p)
    return isodd(w) && iszero(quadratic_defect(a, p)) ? -1 : 1
  end

  v, a = _quadratic_defect_unit(a, p)

  b = b // pi^(2*fld(w, 2)) # fld(x, y) = div(x, y, RoundDown), i.e., round towards -oo
  w = mod(w, 2)
  # turn b into a uniformizer
  if isfinite(v) && isodd(v)
    if w == 0
      b = (1-a)*b // pi^(v-1)
      w = 1
      @assert valuation(b, p) == 1
    end
    k, h = ResidueField(o, p)
    h = extend(h, f)
  end

  # lift v until it becomes infinite or even
  while isfinite(v) && isodd(v)
    t = pi^(v-1)*b
    ok, s = issquare_with_sqrt( h( (a-1) // t) )
    a = a * (1 - (h\(s))^2 * t)
    v, a = _quadratic_defect_unit(a, p)
  end

  return isfinite(v) && isodd(w) ? -1 : 1
end

function hilbert_symbol(a::NumFieldElem, b::IntegerUnion, p::Plc)
  return hilbert_symbol(a, parent(a)(b), p)
end

function hilbert_symbol(a::IntegerUnion, b::NumFieldElem, p::Plc)
  return hilbert_symbol(b, a, p)
end

@doc Markdown.doc"""
    hilbert_symbol(a::nf_elem, b::nf_elem, p::InfPlc) -> Int

Returns the local Hilbert symbol $(a,b)_p$.
"""
function hilbert_symbol(a::T, b::T, p::Plc) where {T <: NumFieldElem}
  return iscomplex(p) || ispositive(a, p) || ispositive(b, p) ? 1 : -1
end

function hilbert_symbol(a::IntegerUnion, b::IntegerUnion, p::IntegerUnion)
  return hilbert_symbol(fmpz(a), fmpz(b), fmpz(p))
end

function hilbert_symbol(a::IntegerUnion, b::IntegerUnion, p::PosInf)
  return hilbert_symbol(a, b, 0)
end

@doc Markdown.doc"""
    hilbert_symbol(a::fmpz, b::fmpz, p::fmpz) -> Int

Returns the local Hilbert symbol $(a,b)_p$.
"""
function hilbert_symbol(a::fmpz, b::fmpz, p::fmpz)
  if p <= 0
    return (a < 0 && b < 0) ? -1 : 1
  end
  v = valuation(a,p)
  w = valuation(b,p)

  if isodd(v)
    if iseven(w)
      a,b,v,w = b,a,w,v
    else
      a,v = -a*b,v+w
    end
  end

  a = divexact(a, p^v)

  if isodd(p)
    return (isodd(w) && jacobi_symbol( mod(a, p), p) == -1) ? -1 : 1
  end

  a = mod(a, 8)
  if a == 1
    return 1
  elseif a == 5
    return isodd(w) ? -1 : 1
  end

  if isodd(w)
    b = divexact(b*(1-a), p^(w+1))
  else
    b = divexact(b, p^w)
  end
  @assert valuation(b, p) == 0
  b = mod(b, 8)

  return (b == 3) || (b == 7) ? -1 : 1
end

function hilbert_symbol(a::Union{fmpq,fmpz,Integer,Rational{<:Integer}},
                        b::Union{fmpq,fmpz,Integer,Rational{<:Integer}},
                        p::IntegerUnion)
  return hilbert_symbol(fmpq(a), fmpq(b), fmpz(p))
end

@doc Markdown.doc"""
    hilbert_symbol(a::Union{fmpq,fmpz,Int,Rational{Int}}, b::Union{fmpq,fmpz,Int,Rational{Int}}, p::Union{fmpz,Int}) -> {-1, +1}

Returns the local Hilbert symbol $(a,b)_p$.
"""
function hilbert_symbol(a::fmpq, b::fmpq, p::fmpz)
  a = FlintQQ(a)
  b = FlintQQ(b)
  hilbert_symbol(numerator(a) * denominator(a), numerator(b) * denominator(b), p)
end


text:
MODULES  = (:BGK,
            :Landau1P)

for mod_name in MODULES
    include(joinpath(dirname(@__FILE__),"..","..","src","BoltzmannSolver","$mod_name.jl"))
end

function test_1v_damping_step(mod_name, Δt, f₀, f₁, v, exp_norm_dev, args...; plot_needed=false, kwargs...)
    if plot_needed
        plt = plot(v, f₀, label = "initial")
        plot!(v, f₁, label = "expected")
    end

    f = copy(f₀)

    if mod_name == :BGK
        damp! = eval(:($mod_name)).generate_solver(f, v, Δt, args...; kwargs...)
    else
        damp! = eval(:($mod_name)).generate_solver(f₀, f, v, Δt, args...; kwargs...)
    end
    for k = 1:100
        damp!()
        # plot!(v, f)
        f₀ = copy(f)
    end

    if plot_needed
        plot!(v, f, label = "$mod_name $args")
    end

    println("$mod_name $args $(values(values(kwargs))): $(norm(f - f₁))")
    @test norm(f - f₁) ≈ 0 atol=exp_norm_dev

    if plot_needed
        display(plt)
    end
end

@testset "Test 1V Boltzmann solvers" begin
    Δt = 0.1
    v = collect(-4:0.1:4)

    f₀ = @. exp(-v^2)
    f₁ = @. exp(-v^2)

    test_1v_damping_step(:BGK, Δt, f₀, f₁, v, 1e-4, 1e-2)
    test_1v_damping_step(:Landau1P, Δt, f₀, f₁, v, 3e-3, 1e-2)

    a = 3e-1
    n = 0.5*sqrt(π)*(a+2.)
    T = 0.25*sqrt(π)*(3a+2.)/n
    f₀ = @. exp(-v^2)*(1.0 + a*v^2)
    f₁ = @. n/sqrt(2π*T)*exp(-v^2/(2T))

    test_1v_damping_step(:BGK, Δt, f₀, f₁, v, 2e-3, 1e-1)
end


text:
function dpll(fla::Formula)::SatSolverResult
  pa = PartialAssignment(fla.variable_count)
  dpll_rec(fla, pa)
end

function dpll_rec(fla::Formula, pa::PartialAssignment)::SatSolverResult
  up = unit_propagation!(fla, pa)
  if false == up
    return UnsatisfiableFormula()
  end

  # Find an undecided variable.
  undecided_var = findfirst((v) -> 0 == v, pa)
  if nothing == undecided_var
    return SatisfiableFormula(pa)
  end

  new_pa = copy(pa)
  new_pa[undecided_var] = 1
  sol_true = dpll_rec(fla, new_pa)
  if is_formula_satisfiable(sol_true)
    return sol_true
  end

  new_pa = copy(pa)
  new_pa[undecided_var] = -1
  sol_false = dpll_rec(fla, new_pa)
  if is_formula_satisfiable(sol_false)
    return sol_false
  end

  UnsatisfiableFormula()
end


text:
using DataFrames: DataFrame

@testset "drive" begin
    @testset "basic" begin
        @system SDrive(Controller) begin
            a => [2, 4, 6] ~ drive
        end
        s = instance(SDrive)
        @test s.a' == 2
        update!(s)
        @test s.a' == 4
        update!(s)
        @test s.a' == 6
    end

    @testset "unit" begin
        @system SDriveUnit(Controller) begin
            a => [2, 4, 6] ~ drive(u"m")
        end
        s = instance(SDriveUnit)
        @test s.a' == 2u"m"
        update!(s)
        @test s.a' == 4u"m"
        update!(s)
        @test s.a' == 6u"m"
    end

    @testset "type" begin
        @system SDriveType(Controller) begin
            a => [2, 4, 6] ~ drive::int(u"m")
        end
        s = instance(SDriveType)
        @test s.a' === 2u"m"
        update!(s)
        @test s.a' === 4u"m"
        update!(s)
        @test s.a' === 6u"m"
    end

    @testset "parameter" begin
        @system SDriveParameter(Controller) begin
            a ~ drive(parameter)
        end
        a = [2, 4, 6]
        c = :0 => :a => a
        s = instance(SDriveParameter; config=c)
        @test s.a' == 2
        update!(s)
        @test s.a' == 4
        update!(s)
        @test s.a' == 6
    end

    @testset "provide" begin
        @system SDriveProvide(Controller) begin
            p => DataFrame(index=(0:2)u"hr", a=[2,4,6], x=1:3, c=(4:6)u"m") ~ provide
            a ~ drive(from=p)
            b ~ drive(from=p, by=:x)
            c ~ drive(from=p, u"m")
        end
        s = instance(SDriveProvide)
        @test s.a' == 2
        @test s.b' == 1
        @test s.c' == 4u"m"
        update!(s)
        @test s.a' == 4
        @test s.b' == 2
        @test s.c' == 5u"m"
        update!(s)
        @test s.a' == 6
        @test s.b' == 3
        @test s.c' == 6u"m"
    end

    @testset "produce" begin
        @system SDriveProduce begin
            i(context.clock.tick) ~ preserve::int
            a => 0:3 ~ drive::int
        end
        @system SDriveProduceController(Controller) begin
            p => produce(SDriveProduce) ~ produce
        end
        s = instance(SDriveProduceController)
        p = s.p'
        @test isempty(p)
        update!(s)
        @test p[1].i' == 0 && p[1].a' == 1
        update!(s)
        @test p[1].i' == 0 && p[1].a' == 2
        @test p[2].i' == 1 && p[2].a' == 1
        update!(s)
        @test p[1].i' == 0 && p[1].a' == 3
        @test p[2].i' == 1 && p[2].a' == 2
        @test p[3].i' == 2 && p[3].a' == 1
    end

    @testset "error" begin
        @test_throws LoadError @eval @system SDriveErrorMissingFrom(Controller) begin
            a ~ drive(by=:a)
        end

        @test_throws LoadError @eval @system SDriveErrorProvideParameter(Controller) begin
            p => DataFrame(index=(0:2)u"hr", a=[2,4,6]) ~ provide
            a ~ drive(from=p, parameter)
        end

        @test_throws LoadError @eval @system SDriveErrorProvideBody(Controller) begin
            p => DataFrame(index=(0:2)u"hr", a=[2,4,6]) ~ provide
            a => [1, 2, 3] ~ drive(from=p)
        end
    end
end


text:
using RayTracer, Images, BenchmarkTools, ArgParse

function render(scene, screen_size, gillum)
    cam = Camera(
        Vec3(0.5f0, 0.2f0, -0.5f0),
        Vec3(0.0f0, 0.1f0,  0.0f0),
        Vec3(0.0f0, 1.0f0,  0.0f0),
        45.0f0,
        1.0f0,
        screen_size...
    )
    origin, direction = get_primary_rays(cam)
    light = PointLight(
        Vec3(1.0f0, 1.0f0, 1.0f0),
        2.0f8,
        Vec3(10.0f0, 10.0f0, 5.0f0)
    )
    image_packed = raytrace(
        origin,
        direction,
        scene,
        light,
        origin,
        gillum
    )
    return zeroonenorm(
               reshape(hcat(image_packed.x, image_packed.y, image_packed.z),
                      (3, screen_size...)))
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--mesh_path", "-m"
            help = "path to the mesh being rendered"
            arg_type = String
            required = true
        "--global_illumination", "-g"
            help = "Activate global illumination (>= 2 for no)"
            arg_type = Int
            default = 2
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Running Benchmarks for Mesh:\n
             1. Path to the Mesh - $(parsed_args["mesh_path"])\n
             2. Global Illumination - $(parsed_args["global_illumination"])")
    println("The benchmarking time involves the following:\n
             1. Instantiating the Camera\n
             2. Rendering\n
             3. Reshaping and Normalizing the Image")
    for screen_size in [(32, 32), (64, 64), (128, 128), (512, 512), (1024, 1024)]
        @info "Screen Size --> $screen_size"
        scene = load_obj(parsed_args["mesh_path"])
        print("Time without BVH --> ")
        @btime render($scene, $screen_size, $(parsed_args["global_illumination"]))
        bvh = BoundingVolumeHierarchy(scene)
        print("Time with BVH --> ")
        @btime render($bvh, $screen_size, $(parsed_args["global_illumination"]))
    end
    println("Benchmarking Completed")
end

main()


text:
function TNmarginal(model::MetNet, state::AbstractMetState, ider)
    ider = rxnindex(model, ider)
    μ_ = μ(model, state, ider)
    σ_ = σ(model, state, ider)
    return Truncated(Normal(μ_, sqrt(σ_)), lb(model, ider), ub(model, ider))
end


text:
#########################################################################
# Project #4: The Second-Order Moller-Plesset
#             Perturbation Theory (MP₂) Energy
#########################################################################

using PyCall: pyimport

using TensorOperations: @tensor

########################################################
#3: TRANSFORM THE TWO-ELECTRON INTEGRALS TO THE MO BASIS
########################################################

function ao2mo_noddy(nao,C,eri,mol)
   eri_mo = zeros(length(eri))
   @inbounds (
   for q in 1:nao, p in q:nao
      for r in 1:p
         lim = ifelse(p == r, q, r)
         for s in 1:lim
            pqrs = get_4idx(p,q,r,s)
            for i in 1:nao, j in 1:nao, k in 1:nao, l in 1:nao
               ijkl = get_4idx(i,j,k,l)
               eri_mo[pqrs] += C[i,p]*C[j,q]*eri[ijkl]*C[k,r]*C[l,s]
            end
         end
      end
   end )
   eri_mo
end

function ao2mo_smart(nao,C,eri,mol)
   """
   (ij|kl) => (ij|ks) => (ij|rs) => (iq|rs) => (pq|rs)
   """
   M      = nao*(nao+1)÷2
   X = zeros(nao, nao)
   tmp = zeros(M,M)
   @inbounds (
   for j in 1:nao, i in j:nao
      ij = index(i,j)
      for l in 1:nao, k in l:nao
         ijkl = get_4idx(i,j,k,l)
         X[k,l] = X[l,k] = eri[ijkl]
      end

      Y = zeros(nao, nao)
      Y = C' * X
      X = zeros(nao, nao)
      X = Y * C
      for l in 1:nao, k in l:nao
         kl = index(k,l)
         tmp[kl,ij] = X[k,l]
      end
   end )

   eri_mo = zeros(M*(M+1)÷2)

   @inbounds (
   for l in 1:nao, k in l:nao
      kl = index(k,l)
      X = zeros(nao,nao)
      Y = zeros(nao,nao)
      for j in 1:nao, i in j:nao
         ij = index(i,j)
         X[i,j] = X[j,i] = tmp[kl,ij]
      end
      Y = zeros(nao, nao)
      Y = C' *  X
      X = zeros(nao, nao)
      X = Y * C
      for j in 1:nao, i in j:nao
         klij = get_4idx(k,l,i,j)
         eri_mo[klij] = X[i,j]
      end
   end )
   eri_mo
end
#eri_mo   = pyscf.ao2mo(nao,C,eri,mol)

########################################################
#4: COMPUTE THE MP₂ ENERGY
########################################################
function short_mp2(e,C,eri,mol)
   nocc  = mol.nelectron ÷ 2
   nao   = mol.nao_nr() |> n -> convert(Int,n) # why not converted ??!
   # Convert AO -> MO basis
   eri_mo = zeros(nao,nocc,nao,nocc)
   emp2 = 0.0
   occupied =      1:nocc
   virtual  = nocc+1:nao
   pyscf = pyimport("pyscf")
   eri = pyscf.ao2mo.restore(1,eri,nao) # reshape
   Cocc = @view C[:,1:nocc]
   @tensor eri_mo[a,i,b,j] := Cocc[p,a] * C[q,i] * Cocc[r,b] * C[s,j] * eri[p,q,r,s]
   @inbounds (
   for i in occupied, a in virtual,
                 j in occupied, b in virtual
       eiajb, eibja = eri_mo[i,a,j,b],  eri_mo[i,b,j,a]
       emp2 += eiajb * (2.0eiajb-eibja) / (e[i]+e[j]-e[a]-e[b])
   end )
   emp2
end

function mymp2(e,C,eri,mol; alg::Symbol=:ao2mo_smart)
   # alg ∈ [ao2mo_noddy, ao2mo_smart]
   ao2mo    = getfield(@__MODULE__, alg) # get function
   nocc     = mol.nelectron ÷ 2
   nao      = mol.nao_nr()       |> n -> convert(Int,n) # why not converted ??!
   emp2     = 0.0
   occupied =      1:nocc
   virtual  = nocc+1:nao
   eri_mo   = ao2mo(nao,C,eri,mol)
   @inbounds (
   for i in occupied, a in virtual,
                 j in occupied, b in virtual
      iajb = get_4idx(i,a,j,b)
      ibja = get_4idx(i,b,j,a)
      emp2 += eri_mo[iajb] * (2.0eri_mo[iajb]-eri_mo[ibja]) / (e[i]+e[j]-e[a]-e[b])
   end )

   # In the basis of spin-molecular orbitals:
   #emp2 += 0.25eri_mo[ijab]*eri_mo[ijab]/(e[i]+e[j]-e[a]-e[b])
   emp2
end


text:
module WGLMakie

using Hyperscript
using JSServe, Observables, AbstractPlotting
using Colors, GeometryBasics
using ShaderAbstractions, LinearAlgebra
using GeometryBasics: GeometryBasics

using JSServe: Application, Session, evaljs, linkjs
using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session
using JSServe: JSObject, onload, uuidstr
using JSServe.DOM
using ShaderAbstractions: VertexArray, Buffer, Sampler, AbstractSampler
using ShaderAbstractions: InstancedProgram
import AbstractPlotting.FileIO
using StaticArrays
using GeometryBasics: decompose_uv
using ImageMagick: ImageMagick

using FreeTypeAbstraction
using AbstractPlotting: get_texture_atlas, glyph_uv_width!, SceneSpace, Pixel
using AbstractPlotting: attribute_per_char, glyph_uv_width!, layout_text

using ImageTransformations

struct WebGL <: ShaderAbstractions.AbstractContext end
struct WGLBackend <: AbstractPlotting.AbstractBackend end

const THREE = JSServe.Dependency(:THREE,
                                 ["https://cdn.jsdelivr.net/gh/mrdoob/three.js/build/three.js"])
const WGL = JSServe.Dependency(:WGLMakie, [joinpath(@__DIR__, "wglmakie.js")])

struct ThreeDisplay <: AbstractPlotting.AbstractScreen
    context::JSObject
end
JSServe.session(td::ThreeDisplay) = JSServe.session(td.context)

function Base.insert!(td::ThreeDisplay, scene::Scene, plot::AbstractPlot)
    js_scene = serialize_three(scene, plot)
    td.context.add_plot(js_scene)
    return
end

"""
    get_plot(td::ThreeDisplay, plot::AbstractPlot)

Gets the ThreeJS object representing the plot object.
"""
function get_plot(td::ThreeDisplay, plot::AbstractPlot)
    return td.context.get_plot(string(objectid(plot)))
end

include("serialization.jl")
include("events.jl")

include("particles.jl")
include("lines.jl")
include("meshes.jl")
include("imagelike.jl")
include("display.jl")

function activate!()
    b = WGLBackend()
    AbstractPlotting.register_backend!(b)
    AbstractPlotting.set_glyph_resolution!(AbstractPlotting.Low)
    AbstractPlotting.current_backend[] = b
    return AbstractPlotting.inline!(true) # can't display any different atm
end

function __init__()
    # Activate WGLMakie as backend!
    return activate!()
end

end # module


text:
include("getdigit_n.jl")
include("myprintf.jl")
include("simpson_noparallel.jl")
using Printf
using .GetDigit_n
using .MyPrintf
using .Simpson_noparallel

const DIGIT = 12
const MAXN = 280000000

function circle_4(n)
    return BigFloat(4) * Simpson_noparallel.simpson(x::BigFloat -> sqrt(BigFloat(1) - x * x), BigFloat(0), BigFloat(1), n)
end

function main()
    n = GetDigit_n.getn_even(circle_4, MAXN, DIGIT, false)
    precision = GetDigit_n.getdigit(circle_4, n, DIGIT, false)
    precision = GetDigit_n.getdigit_2(circle_4, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, circle_4(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
OpenXML(filename) =
    Lift(f -> parsexml(readstring(f)), convert(Combinator, filename))

XMLRoot() =
    Lift(getroot, It)

XMLRoot(name; kw...) =
    XMLRoot() >> ThenFilter(_xmlchildfilter(name; kw...))

XMLName() =
    Lift(getname, It)

LiftQuery{V}(fn, argtypes::Tuple{Vararg{Type}}, restype::Type{XMLElementVector{V}}, qs::Vector{Query}) =
    let q = invoke(LiftQuery, (Any, Tuple{Vararg{Type}}, Type, Vector{Query}), fn, argtypes, restype, qs)
        q >> Query(DecodeVectorSig(), Input(XMLElementVector{V}), Output(XMLElement) |> setoptional() |> setplural())
    end

XMLParent() =
    Lift(getparent, It)

XMLChild() =
    Lift(getchildren, It)

XMLChild(name; kw...) =
    XMLChild() >> ThenFilter(_xmlchildfilter(name; kw...))

XMLDescendant() =
    Lift(getdescendants, It)

XMLDescendant(name; kw...) =
    XMLDescendant() >> ThenFilter(_xmlchildfilter(name; kw...))

function _xmlchildfilter(name; kw...)
    C = XMLName() .== name
    for (k, v) in kw
        if v == true
            C = C & Exists(XMLAttribute(string(k)))
        elseif v == false
            C = C & ~Exists(XMLAttribute(string(k)))
        else
            C = C & (XMLAttribute(string(k)) .== v)
        end
    end
    C
end

XMLText() =
    Lift(gettext, It)

XMLAttribute(attr) =
    Lift(getattribute, It, convert(Combinator, attr))

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto