text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Binomial" begin
    t = BinomialTest(26, 78)
    @test pvalue(t) ≈ 0.004334880883507431
    @test pvalue(t, tail=:left) ≈ 0.002167440441753716
    @test pvalue(t, tail=:right) ≈ 0.9989844298129187
    @test default_tail(t) == :both
    @test_ci_approx confint(t) (0.23058523962930383, 0.4491666887959782)
    @test_ci_approx confint(t, tail=:left) (0.0, 0.4313047758370174)
    @test_ci_approx confint(t, tail=:right) (0.2451709633730693, 1.0)
    @test_ci_approx confint(t, method=:wald) (0.22871819521037956, 0.43794847145628707)
    @test_ci_approx confint(t, tail=:left, method=:wald) (0.0, 0.42112912485444692)
    @test_ci_approx confint(t, tail=:right, method=:wald) (0.24553754181221971, 1.0)
    @test_ci_approx confint(t, method=:wilson) (0.23872670036358601, 0.44358590287381217)
    @test_ci_approx confint(t, tail=:left, method=:wilson) (0.0, 0.42541288951088108)
    @test_ci_approx confint(t, tail=:right, method=:wilson) (0.25242832328277831, 1.0)
    @test_ci_approx confint(t, method=:jeffrey) (0.23626570247518358, 0.44251318323879296)
    @test_ci_approx confint(t, tail=:left, method=:jeffrey) (0.0, 0.42466492683653623)
    @test_ci_approx confint(t, tail=:right, method=:jeffrey) (0.25098836986261724, 1.0)
    @test_ci_approx confint(t, method=:agresti_coull) (0.2384423809121706, 0.44387022232522744)
    @test_ci_approx confint(t, tail=:left, method=:agresti_coull) (0.0, 0.42558712894362222)
    @test_ci_approx confint(t, tail=:right, method=:agresti_coull) (0.25225408385003706, 1.0)
    @test_ci_approx confint(t, method=:arcsine) (0.23366209634204066,0.44117918327686334)
    @test_ci_approx confint(t, tail=:left, method=:arcsine) (0.0,0.4235046425920888)
    @test_ci_approx confint(t, tail=:right, method=:arcsine) (0.2489264087216164,1.0)

    show(IOBuffer(), t)

    t = BinomialTest([trues(6); falses(3)])
    @test pvalue(t) ≈ 0.5078125000000002
    @test pvalue(t, tail=:left) ≈ 0.91015625
    @test pvalue(t, tail=:right) ≈ 0.2539062500000001
    @test_ci_approx confint(t) (0.2992950562085405, 0.9251453685803082)
    @test_ci_approx confint(t, tail=:left) (0.0, 0.9022531865607242)
    @test_ci_approx confint(t, tail=:right) (0.3449413659437032, 1.0)
    @test_ci_approx confint(t, method=:wald) (0.35868803903340479, 0.97464529429992841)
    @test_ci_approx confint(t, tail=:left, method=:wald) (0.0,0.92513047859481645)
    @test_ci_approx confint(t, tail=:right, method=:wald) (0.40820285473851681,1.0)
    @test_ci_approx confint(t, method=:wilson) (0.35420213558039609,0.87941618161308899)
    @test_ci_approx confint(t, tail=:left, method=:wilson) (0.0,0.85802909820500495)
    @test_ci_approx confint(t, tail=:right, method=:wilson) (0.39825972868840931,1.0)
    @test_ci_approx confint(t, method=:jeffrey) (0.34779179347226591,0.89578677833922582)
    @test_ci_approx confint(t, tail=:left, method=:jeffrey) (0.0,0.86830830610561005)
    @test_ci_approx confint(t, tail=:right, method=:jeffrey) (0.39604343455469687,1.0)
    @test_ci_approx confint(t, method=:agresti_coull) (0.350905767251112,0.88271254994237336)
    @test_ci_approx confint(t, tail=:left, method=:agresti_coull) (0.0,0.86049746046629294)
    @test_ci_approx confint(t, tail=:right, method=:agresti_coull) (0.39579136642712159,1.0)
    @test_ci_approx confint(t, method=:arcsine) (0.345812446615087,0.9188773496172281)
    @test_ci_approx confint(t, tail=:left, method=:arcsine) (0.0,0.8879439981269358)
    @test_ci_approx confint(t, tail=:right, method=:arcsine) (0.3965293068864491,1.0)
    show(IOBuffer(), t)

    t = BinomialTest(0, 100, 0.01)
    @test pvalue(t) ≈ 0.7320646825464591
    show(IOBuffer(), t)

    t = BinomialTest(100, 100, 0.99)
    @test pvalue(t) ≈ 0.7320646825464584
    show(IOBuffer(), t)
end

@testset "SignTest" begin
    x = [55, 58, 61, 61, 62, 62, 62, 63, 63, 64, 66, 68, 68, 69, 69, 69, 70, 71, 72, 72]
    @test pvalue(SignTest(x)) ≈ 1.907348632812499e-6
    @test pvalue(SignTest(x), tail=:left) ≈ 1.0
    @test pvalue(SignTest(x), tail=:right) ≈ 9.536743164062495e-7
    @test pvalue(SignTest(x, 70)) ≈ 0.004425048828125003
    @test pvalue(SignTest(x, 70), tail=:left) ≈ 0.0022125244140625013
    @test pvalue(SignTest(x, 70), tail=:right) ≈ 0.9996356964111328
    @test default_tail(SignTest(x)) == :both
    @test_ci_approx confint(SignTest(x, 70)) (62, 69)
    @test_ci_approx confint(SignTest(x, 70), level=0.9998) (61, 71)
    show(IOBuffer(), SignTest(x, 70))

    x = [9, 2, 7, 5]
    y = [7, 2, 6, 4]
    @test pvalue(SignTest(x, y)) ≈ 0.25
    @test_ci_approx confint(SignTest(x, y)) (0, 2)
    show(IOBuffer(), SignTest(x, y))

    # www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
    x = [-4.7, 3.7, 22.4, 13.6, 8.7, 9.1, -7.8, 10.8, 15.6, 23.5, 14.4, 20.2, 6.5, 10.1, -6.9]
    @test pvalue(SignTest(x)) ≈ 0.03515625
    @test pvalue(SignTest(x), tail=:left) ≈ 0.996307373046875
    @test pvalue(SignTest(x), tail=:right) ≈ 0.017578125000000007
    @test_ci_approx confint(SignTest(x), level=0.99995) (-7.8, 23.5)
    @test_ci_approx confint(SignTest(x), level=0.9999) (-6.9, 22.4)
    @test_ci_approx confint(SignTest(x), level=0.993) (-4.7, 20.2)
    @test_ci_approx confint(SignTest(x), level=0.965) (3.7, 15.6)
    @test_ci_approx confint(SignTest(x), level=0.882) (6.5, 14.4)
    @test_ci_approx confint(SignTest(x), level=0.7) (8.7, 13.6)
    @test_ci_approx confint(SignTest(x), level=0.6) (9.1, 10.8)
    show(IOBuffer(), SignTest(x))
end


text:
# get all faces, including
#   bboxes, landmarks, features
using Glob
using PyCall
using JLD2
using FileIO

const FR = pyimport("evolveface")
const PIL = pyimport("PIL")

# get all templates
template_files = glob("data/*.jpg")
template_features = hcat(map(template_files) do f
    img = PIL.Image.open(f)
    bounding_boxes, landmarks = FR.detect_faces(img)
    features = FR.extract_feature_IR50A(img, landmarks)
    @assert size(features)[1] == 1
    return vec(features)
end...)

@save "out/template.jld2" template_features

map(glob("data/vid/*.jpg")) do f
    println(f)
    img = PIL.Image.open(f)
    bounding_boxes, landmarks = FR.detect_faces(img)
    size(landmarks)[1] == 0 && return
    features = FR.extract_feature_IR50A(img, landmarks)
    res = (bounding_boxes, landmarks, features)
    @save "out/$(splitext(basename(f))[1]).jld2" res
end

ext: jl

max_stars_repo_name: BillKerman/FaceNetCustomized

max_issues_repo_name: BillKerman/FaceNetCustomized

max_forks_repo_name: BillKerman/FaceNetCustomized


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %%
# Load the data, don't forget to drop everything missing.
d = DataFrame(CSV.File(datadir("exp_raw/cherry_blossoms.csv"), missingstring = "NA"))
dropmissing!(d, [:temp, :doy])

# scatter(d.temp, d.doy, alpha = 0.3)

# %%
# First, a pure linear model.
flin(temp, α, β) = α .+ β .* temp

@model function cherry_linear(temp, doy)
    # β: With a higher temperature I expect the doy to be earlier. Maybe 1 month per 10°C?
    α ~ Normal(100, 10)
    β ~ Normal(-3, 2)
    σ ~ Exponential(1)
    μ = flin(temp, α, β)
    doy ~ MvNormal(μ, σ)
end

m = cherry_linear(d.temp, d.doy)

# %%
# Sample from the prior. Looks pretty reasonable.
prior = sample(m, Prior(), 100) |> DataFrame

scatter(d.year,  d.doy, alpha = 0.3)
for r in eachrow(prior)
    p = flin(d.temp, r.α, r.β)
    plot!(d.year, p, alpha = 0.2)
end
plot!(legend = false, ylims = (0, 200))

# %%
# Fit and show the fit.
q = quap(m)
post = DataFrame(rand(q.distr, 1000)', ["α", "β", "σ"])

mu = flin(d.temp, post.α', post.β')
mu_m, mu_lower, mu_upper = meanlowerupper(mu)
sim = rand.(Normal.(mu, post.σ'))
sim_m, sim_lower, sim_upper = meanlowerupper(sim)

scatter(d.year, d.doy, alpha = 0.3, legend = false)
plot!(d.year, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m))
plot!(d.year, sim_lower, fillrange = sim_upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %%
# Second, a linear model enhanced with extra wiggle.
using BSplines: BSplineBasis, basismatrix

num_knots = 15
knot_list = quantile(d.year, range(0, 1, length = num_knots))

basis = BSplineBasis(4, knot_list)
B = basismatrix(basis, d.year)

# %%
# Fit the model.
@model function spline(temp, doy, B = B)
    α ~ Normal(100, 10)
    β ~ Normal(-3, 2)
    w ~ filldist(Normal(0, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = flin(temp, α, β) .+ B * w
    doy ~ MvNormal(μ, σ)
end

sample(spline(d.temp, d.doy), NUTS(), 1000)

qspl = quap(spline(d.temp, d.doy))

w_str = ["w[$i]" for i in 1:length(basis)]
post = DataFrame(rand(qspl.distr, 1000)', ["α"; "β"; w_str; "σ"])

smu = fspl(d.temp, post.α', post.β', Array(post[!, w_str])', B)
smu_m, smu_lower, smu_upper = meanlowerupper(smu)

# %%
# Display both solutions together.
scatter(d.year, d.doy, alpha = 0.3, label = "data")
plot!(d.year, smu_m, ribbon = (smu_m .- smu_lower, smu_upper .- smu_m), label = "linear")
plot!(d.year, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m), label = "linear+spline")


text:
"Logger context selector"
abstract ContextSelector

"Return current context"
context(ctxsel::ContextSelector, fqmn::AbstractString, cfgloc::AbstractString) = throw(AssertionError("Function 'context' is not implemented for type $(typeof(ctxsel))"))
context(ctxsel::ContextSelector, cfgloc::AbstractString="") = context(ctxsel, string(current_module()), cfgloc)

"Return all contexts"
contexts(ctxsel::ContextSelector) = throw(AssertionError("Function 'contexts' is not implemented for type $(typeof(ctxsel))"))

"Remove specified logger context"
delete!(ctxsel::ContextSelector, ctx::LoggerContext) = throw(AssertionError("Function 'contexts' is not implemented for type $(typeof(ctxsel))"))

####################################
# Implementations of ContextSelector
####################################

"""Single context selector

Simple logger context selector that keeps only one context inside itself.
"""
type SingleContextSelector <: ContextSelector
    default::LoggerContext
    SingleContextSelector() = new(LoggerContext("Default"))
end
context(ctxsel::SingleContextSelector, fqmn::AbstractString, cfgloc::AbstractString="") = ctxsel.default
contexts(ctxsel::SingleContextSelector) = Dict(:Default=>ctxsel.default)

"""Module context selector choses `LoggerContext` based on module name

All `LoggerContext` instances are stored inside a module where `Log4jl` configuration was performed in a constant with unique name `__Log4jl_LC<ID>__`
"""
type ModuleContextSelector <: ContextSelector
    default::LoggerContext
    ModuleContextSelector() = new(LoggerContext("Default"))
end

function context(ctxsel::ModuleContextSelector,
                 fqmn::AbstractString,
                 cfgloc::AbstractString="")

    # return default context if module doesn't exist
    isempty(fqmn) && return ctxsel.default

    # get module for context selector
    ctx_module = getmodule(fqmn)

    # return context object in the module
    mlc_var = mlcvar(string(fqmn))
    if isconst(ctx_module, mlc_var)
        trace(LOGGER, "Found context: $mlc_var")
        return getfield(ctx_module, mlc_var)
    end

    ctx = ctxsel.default

    # check configuration source
    cfgloc = locateconfig(cfgloc, ctx_module)
    if isempty(cfgloc)
        trace(LOGGER, "Invalid configuration source")

        # look for parent context
        parent_ctx_module = module_parent(ctx_module)
        found = false
        for vs in names(parent_ctx_module, true)
            if startswith(string(vs), LOG4JL_CTX_PREFIX)
                ctx = getfield(parent_ctx_module, vs)
                trace(LOGGER, "Found parent context: $vs")
                found = true
                break
            end
        end
        !found && trace(LOGGER, "No parent context found.")
    else
        trace(LOGGER, "Creating new context in the module: $ctx_module")

        # otherwise create a new context
        ctx = LoggerContext(string(mlc_var), cfgloc)

        # define constant in module with LoggerContext object
        ccall(:jl_set_const, Void, (Any, Any, Any), ctx_module,  mlc_var, ctx)
    end

    # return default context if non found
    return ctx
end

function contexts(ctxsel::ModuleContextSelector)

    function submodules(m::Module, mdls::Vector{Symbol}= Symbol[], ctxs = Dict{Symbol,LoggerContext}())
        for vs in names(m, true)
            !isdefined(m, vs) && continue
            v = getfield(m, vs)
            if isa(v, Module) && vs ∉ [:Main, :Core, :Base]
                fqmn = symbol(v)
                fqmn in mdls && continue
                push!(mdls, fqmn)
                mlc_var = mlcvar(string(fqmn))
                if isdefined(v, mlc_var)
                    ctxs[fqmn] = getfield(v, mlc_var)
                end
                submodules(v, mdls, ctxs)
            end
        end
        return mdls, ctxs
    end

    mdls, ctxs = submodules(Main)
    ctxs[:Default] = ctxsel.default
    return ctxs
end


text:
# Script that plots the time series of the model versus observations
# - observations are automatically extracted from tides and currents
# - only time series where both model and observations are sufficiently complete
#   are plotted
# - Plots daily maximum water levels from the timeseries in addition
#   to full timeseries and lowpass filtered timeseries
using NetCDF
using Plots
using Dates
using HTTP
using CSV
using Statistics
using TimeSeries

using FFTW
include("fourfilt.jl")

# using GR backend, behaves like matlab
gr()

nodatalength = 0.0 # decimal percentage of period (if non-zero then the lowpass filter does not work...)
# low pass filter length
          #days
lpwindow = 28*3600.0*24.0 # seconds
dpiset = 300

# filenames
filename = "ATT/Mv13.61.nc"
sn = ncread(filename, "station_name")
# read parts of the netcdf file, convert to datetimes
time = ncread(filename, "time")
zeta = ncread(filename, "zeta")
zeta[zeta .< -999] .= NaN
lon = ncread(filename, "x")
lat = ncread(filename, "y")
dt = DateTime(2004,6,15) + Dates.Second.(time)
# Make sure only from July 1st
zeta = zeta[:,dt .>= DateTime(2004,7,1)]
deleteat!(dt,dt .< DateTime(2004,7,1))
tmstep = (Dates.value(dt[2]-dt[1]))/1000 #seconds

# This is for reading observations through tidesandcurrents API
base = "https://tidesandcurrents.noaa.gov/api/datagetter?"
suf = "product=hourly_height&datum=MSL&units=metric&time_zone=gmt&application=ports_screen&format=csv"
sd  = join(["begin_date=",Dates.format(dt[1],"yyyymmdd"),"&"])
ed  = join(["end_date=",Dates.format(dt[end],"yyyymmdd"),"&"])

# legend names
linenames = [:Observed, :Model]

# initializing the scatter arrays
xp = Float64[]; yp = Float64[]; mv = Float64[]; ov = Float64[]

# loop over all stations
for ii = 1:size(sn,2)

    # only get stations in our AOI
    #if (lon[ii] > -75 || lon[ii] < -83 ||
    #lat[ii] > 36.5 || lat[ii] < 24); continue; end
    if (lon[ii] > -67 || lon[ii] < -100 ||
        lat[ii] > 46 || lat[ii] < 24); continue; end

    # if modelled doesn't have enough values
    if (sum(isnan.(zeta[ii,:])) > length(dt)*nodatalength); continue; end

    # get and plot the station data
    sta = join(["station=",rstrip(join(sn[:,ii])),"&"])
    res = HTTP.get(join([base,sd,ed,sta,suf]))
    mycsv = CSV.read(res.body)
    # if got error or missing time
    if (ismissing(mycsv[1,1]) || mycsv[1,1][1:5] == "Error"); continue; end
    dtn = DateTime.(mycsv[:,1],"yyyy-mm-dd HH:MM")
    # if observed doesn't have enough non-missing values
    if (sum(ismissing.(mycsv[:,2])) > length(dtn)*nodatalength); continue; end

    # Display the station IDs and indices
    display(ii)
    display(join(sn[:,ii]))

    # Process observed results
    zo = mycsv[:,2]
    zo[ismissing.(mycsv[:,2])] .= NaN
    # remove low-pass filtered timeseries
    tostep = (Dates.value(dtn[2]-dtn[1]))/1000 #seconds
    zofilt = fourfilt(zo,tostep,Inf,lpwindow)
    zoaf = zo .- zofilt
    # convert observed data into the timeseries array
    to = TimeArray(dtn,zoaf)
    # calculate daily maximum
    to = collapse(to,day,first,maximum)
    rename!(to, linenames[1])

    # Process modeled results
    zm = zeta[ii,:]
    # remove low-pass filtered timeseries
    zmfilt = fourfilt(zm,tmstep,Inf,lpwindow)
    zmaf = zm .- zmfilt
    # convert model data into the timeseries array
    tm = TimeArray(dt,zmaf)
    # calculate daily maximum
    tm = collapse(tm,day,first,maximum)
    rename!(tm, linenames[2])

    # plot the observations and model
    plt1 = plot(to,
        title = join(["Daily Maximum Water Levels at Station ID: ",
            rstrip(join(sn[:, ii]))]
        ),
        legend = :topright,
        ylabel = "Water Level (LMSL) [m]",
        xlabel = "DateTime [UTC]"
    )
    plot!(tm)

    plt2 = plot(dtn,zo)
    plot!(dt,zm)
    plot!(dtn,zofilt)
    plot!(dt,zmfilt)
    plt = plot(plt1,plt2,dpi=dpiset,layout = (2,1))

    display(plt) # display in plot pane

    # calculate max WL error
    vtm = values(tm)
    vto = values(to)
    vtm = maximum(vtm[.!isnan.(vtm)])
    vto = maximum(vto[.!isnan.(vto)])

    # add the location plot to show error at stations
    push!(xp,lon[ii])
    push!(yp,lat[ii])
    push!(mv,vtm)
    push!(ov,vto)

    savefig(join(["ATT/DailyWLPlots_Mv13/" rstrip(join(sn[:, ii]))]))

end

# write out the maximum water levels over this time period
using DataFrames
df = DataFrame(Lon = xp, Lat = yp, MaxWL_Obs = ov, MaxWL_Mod = mv)
CSV.write("ATT/DailyWLPlots_Mv13/MaxWLs.csv",df)

# I'm having trouble using scatter here because of clash between Plots and GR package
#plt = scatter(xp,yp,200*ones(length(xp)),error) #,legend=:none)
#display(plt) # display in plot pane
#savefig(join(["ATT/DailyWLPlots_Mv13/MaxWLError"))


text:
module Grids

export AbstractGrid
export RegularRectilinearGrid
export RegularLatLonGrid
export VerticallyStretchedLatLonGrid, LoadVerticallyStretchedLatLonGrid
export Periodic, Bounded, short_show
export replace_grid_storage
export ΔxC, ΔyC, ΔzC, ΔxF, ΔyF, ΔzF, Ax, Ay, Az, volume

using CUDA
using Adapt

using PlanktonIndividuals.Architectures

"""
    AbstractGrid{TX, TY, TZ}
Abstract type for grids with elements of type `Float64` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{TX, TY, YZ} end

"""
    AbstractTopology
Abstract type for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic
Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded
Grid topology for bounded dimensions.
"""
struct Bounded <: AbstractTopology end

import Base: show

include("regular_rectilinear_grid.jl")
include("regular_lat_lon_grid.jl")
include("vertically_stretched_lat_lon_grid.jl")
include("utils.jl")
include("areas_volumes_spacings.jl")

end


text:
export iauTrxp
"""
Multiply a p-vector by the transpose of an r-matrix.

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

Status:  vector/matrix support function.

Given:
   r        double[3][3]   r-matrix
   p        double[3]      p-vector

Returned:
   trp      double[3]      r * p

Note:
   It is permissible for p and trp to be the same array.

Called:
   iauTr        transpose r-matrix
   iauRxp       product of r-matrix and p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauTrxp(double r[3][3], double p[3], double trp[3])

function iauTrxp(r::AbstractMatrix{<:Real}, p::AbstractVector{<:Real})
   trp = zeros(Float64, 3)

   ccall((:iauTrxp, libsofa_c), Cvoid,
       (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 2}, r'),
       convert(Array{Float64, 1}, p),
       trp)

   return SVector{3}(trp)
end


text:
module TestCategories
using Test

using Catlab, Catlab.Theories
using Catlab.CategoricalAlgebra.Sets, Catlab.CategoricalAlgebra.Categories

# Instances
###########

# Categories
#-----------

C = TypeCat(FreeCategory.Ob, FreeCategory.Hom)
@test Ob(C) == TypeSet(FreeCategory.Ob)
@test sprint(show, C) == "TypeCat($(FreeCategory.Ob), $(FreeCategory.Hom))"

x, y = Ob(FreeCategory, :x, :y)
f = Hom(:f, x, y)
@test ob(C, x) == x
@test hom(C, f) == f

# Functors
#---------

F = id(C)
@test (dom(F), codom(F)) == (C, C)
@test startswith(sprint(show, F), "id(TypeCat(")
@test F(x) == x
@test F(f) == f

function InstanceFunctor(T::TypeCat{Ob,Hom}) where {Ob, Hom}
  Functor(x -> functor((Ob,Hom), x), f -> functor((Ob,Hom), f),
          TypeCat{FreeCategory.Ob, FreeCategory.Hom}(), T)
end
F = InstanceFunctor(TypeCat(FreeCategory2.Ob, FreeCategory2.Hom))
@test dom(F) == C

x′, y′ = Ob(FreeCategory2, :x, :y)
f′ = Hom(:f, x′, y′)
@test F(x) == x′
@test F(f) == f′

F_op = op(F)
@test F_op isa Categories.OppositeFunctor
@test (dom(F_op), codom(F_op)) == (op(dom(F)), op(codom(F)))
@test F_op(x) == x′
@test F_op(f) == f′

# Transformations
#----------------

α = id(F)
@test (dom(α), codom(α)) == (F, F)
@test component(α, x) == id(F(x))

end

max_forks_repo_name: epatters/CompCat


text:
include("header.jl")
include("simulator.jl")
include("visual.jl")

function get_conds!(conds)
    for i_exp in 1:size(conds)[1]
        local T0 = conds[i_exp, 1]
        global P = conds[i_exp, 2] * one_atm
        local phi = conds[i_exp, 3]
        conds[i_exp, 4] = f_idt(T0, P, phi, p; dT=dT)
        @printf("(%d) T %.1f [K] P %.1f [atm] phi %.1f IDT0 = %.2e [s] \n",
                i_exp, conds[i_exp, 1], conds[i_exp, 2],
                conds[i_exp, 3], conds[i_exp, 4])
    end
    return conds
end

conds = zeros(20, 4)
conds[:, 1] .= range(700.0, stop=1600, length=size(conds)[1])
conds[:, 2] .= 40.0
conds[:, 3] .= 0.8

# dir_master = conf["master"]
# @load string(ckpt_path, "/../../", dir_master, "/conds.bson") conds

gas = CreateSolution("mechanism/mech_41s168r.yaml");
ns = gas.n_species;
nr = gas.n_reactions;
p = zeros(nr * 3);
conds = get_conds!(conds)
conds_master = deepcopy(conds)

gas = CreateSolution("mechanism/mech_34s121r.yaml");
ns = gas.n_species;
nr = gas.n_reactions;
p = zeros(nr * 3);

conds = get_conds!(conds)
conds_sk = deepcopy(conds)

# gas = CreateSolution("mechanism/gri30_sk23_op.yaml");
# ns = gas.n_species;
# nr = gas.n_reactions;
# p = zeros(nr * 3);

@load string(ckpt_path, "/../p.bson") p
conds = get_conds!(conds)
conds_sk_op = deepcopy(conds)

if std(conds[:, 2]) < 1.e-3
    plt = Plots.plot(1000 ./ conds_master[:, 1], conds_master[:, 4],
                lw=1, ls=:solid, label="Nordin_41s")
else
    plt = Plots.scatter(1000 ./ conds_master[:, 1], conds_master[:, 4],
                label="Nordin_41s")
end
Plots.scatter!(plt, 1000 ./ conds_sk[:, 1], conds_sk[:, 4],
            label="sk34")
Plots.scatter!(plt, 1000 ./ conds_sk_op[:, 1], conds_sk_op[:, 4],
            label="sk34_op")
xlabel!(plt, "1000 / T [K]")
ylabel!(plt, "IDT [s]")
plot!(plt, yscale=:log10, legend=:left)
title!(plt, @sprintf("p=%.1f atm phi=%.1f",
        conds[1, 2], conds[1, 3]))
png(plt, string(fig_path, "/regression_T"))

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
# Given a list of atoms, create a cartesian displamentent in one of them (A)
function create_displacement(BS::BasisSet, A::Int, i::Int, h)

    Aplus = deepcopy(BS.atoms)
    Aplus[A].xyz[i] += h

    Aminus = deepcopy(BS.atoms)
    Aminus[A].xyz[i] -= h

    bs_plus = BasisSet(BS.name, Aplus)
    bs_minus = BasisSet(BS.name, Aminus)

    return bs_plus, bs_minus
end

function ∇FD_overlap(BS::BasisSet, A, i, h = 1e-5)
    return ∇FD_1e(BS, "overlap", A, i, h)
end

function ∇FD_kinetic(BS::BasisSet, A, i, h = 1e-5)
    return ∇FD_1e(BS, "kinetic", A, i, h)
end

function ∇FD_nuclear(BS::BasisSet, A, i, h = 1e-5)
    return ∇FD_1e(BS, "nuclear", A, i, h)
end

function ∇FD_1e(BS::BasisSet, compute::String, A, i, h)

    bs_plus, bs_minus = create_displacement(BS, A, i, h)

    Xplus = ao1e(bs_plus, compute)
    Xminus = ao1e(bs_minus, compute)

    return (Xplus - Xminus) ./ (2*h/Molecules.bohr_to_angstrom)
end

function ∇FD_ERI_2e4c(BS::BasisSet, A, i, h=1e-5)

    bs_plus, bs_minus = create_displacement(BS, A, i, h)

    Xplus = ERI_2e4c(bs_plus)
    Xminus = ERI_2e4c(bs_minus)

    return (Xplus - Xminus) ./ (2*h/Molecules.bohr_to_angstrom)

end

function ∇FD_sparseERI_2e4c(BS::BasisSet, A, i, h=1e-5)

    bs_plus, bs_minus = create_displacement(BS, A, i, h)

    Iplus, Xplus = sparseERI_2e4c(bs_plus)
    Iminus, Xminus = sparseERI_2e4c(bs_minus)

    @assert Iplus == Iminus

    Xout = (Xplus - Xminus) ./ (2*h/Molecules.bohr_to_angstrom)

    return Iplus, Xout
end

function ∇FD_ERI_2e3c(BS::BasisSet, auxBS::BasisSet, A, i, h=1e-5)

    bs_plus, bs_minus = create_displacement(BS, A, i, h)
    Abs_plus, Abs_minus = create_displacement(auxBS, A, i, h)

    Xplus  = ERI_2e3c(bs_plus, Abs_plus)
    Xminus = ERI_2e3c(bs_minus, Abs_minus)

    return (Xplus - Xminus) ./ (2*h/Molecules.bohr_to_angstrom)

end

function ∇FD_ERI_2e2c(BS::BasisSet, A, i, h=1e-5)

    bs_plus, bs_minus = create_displacement(BS, A, i, h)

    Xplus = ERI_2e2c(bs_plus)
    Xminus = ERI_2e2c(bs_minus)

    return (Xplus - Xminus) ./ (2*h/Molecules.bohr_to_angstrom)

end


text:
"""
    gbset(A::GBArray, option, value)
    gbset(option, value)

Set an option either for a specific GBArray, or globally. The commonly used options are:
    - `:format = [:byrow | :bycol]`: The global default or array specific
    column major or row major ordering.
    - `:nthreads = [Integer]`: The global number of OpenMP threads to use.
    - `:burble = [Bool]`: Print diagnostic output.
    - `:sparsity_control = [:full | :bitmap | :sparse | :hypersparse]`: Set the sparsity of a
    single GBArray.
"""
gbset

function gbset(option, value)
    option = option_toconst(option)
    value = option_toconst(value)
    libgb.GxB_Global_Option_set(option, value)
    return nothing
end

function gbget(option)
    option = option_toconst(option)
    return libgb.GxB_Global_Option_get(option)
end

function gbset(A::GBMatrix, option, value)
    option = option_toconst(option)
    value = option_toconst(value)
    libgb.GxB_Matrix_Option_set(A, option, value)
    return nothing
end

function gbget(A::GBMatrix, option)
    option = option_toconst(option)
    return libgb.GxB_Matrix_Option_get(A, option)
end

function gbset(A::GBVector, option, value)
    option = option_toconst(option)
    value = option_toconst(value)
    libgb.GxB_Matrix_Option_set(A, option, value)
    return nothing
end

function gbget(A::GBVector, option)
    option = option_toconst(option)
    return libgb.GxB_Matrix_Option_get(A, option)
end

function format(A::GBVecOrMat)
    t = gbget(A, SPARSITY_STATUS)
    f = gbget(A, FORMAT)
    return (GBSparsity(t), GBFormat(f))
end

const HYPER_SWITCH = libgb.GxB_HYPER_SWITCH
const BITMAP_SWITCH = libgb.GxB_BITMAP_SWITCH
const FORMAT = libgb.GxB_FORMAT
const SPARSITY_STATUS = libgb.GxB_SPARSITY_STATUS
const SPARSITY_CONTROL = libgb.GxB_SPARSITY_CONTROL
const BASE1 = libgb.GxB_PRINT_1BASED
const NTHREADS = libgb.GxB_GLOBAL_NTHREADS
const BURBLE = libgb.GxB_BURBLE

const BYROW = libgb.GxB_BY_ROW
const BYCOL = libgb.GxB_BY_COL

#only translate if it's a symbol
option_toconst(option) = option
function option_toconst(sym::Symbol)
    sym === :format && return FORMAT
    sym === :nthreads && return NTHREADS
    sym === :burble && return BURBLE
    sym === :byrow && return BYROW
    sym === :bycol && return BYCOL
    sym === :sparsity_status && return SPARSITY_STATUS
    sym === :sparsity_control && return SPARSITY_CONTROL
    sym === :full && return GBDENSE
    sym === :bitmap && return GBBITMAP
    sym === :sparse && return GBSPARSE
    sym === :hypersparse && return GBHYPER
end

"""
Sparsity options for GraphBLAS. values can be summed to produce additional options.
"""
@enum GBSparsity::Int32 begin
    GBDENSE = 8 #libgb.GxB_FULL
    GBBITMAP = 4 #libgb.GxB_BITMAP
    GBSPARSE = 2 #libgb.GxB_SPARSE
    GBHYPER = 1 #libgb.GxB_HYPERSPARSE
    GBANYSPARSITY = 15 #libgb.GxB_ANY_SPARSITY
    GBDENSE_OR_BITMAP = 12 #libgb.GxB_FULL + libgb.GXB_BITMAP
    GBSPARSE_OR_HYPER = 3 #libgb.GxB_SPARSE + libgb.GXB_HYPERSPARSE
    #... Probably don't need others
end


text:
using LinearAlgebra
using Pathfinder

function rand_pd_mat(T, n)
    U = qr(randn(T, n, n)).Q
    return Matrix(Symmetric(U * rand_pd_diag_mat(T, n) * U'))
end

rand_pd_diag_mat(T, n) = Diagonal(rand(T, n))

# defined for testing purposes
function Pathfinder.rand_and_logpdf(rng, dist, ndraws)
    x = rand(rng, dist, ndraws)
    if x isa AbstractVector
        logpx = Distributions.logpdf.(dist, x)
        x = collect(x')
    else
        logpx = Distributions.logpdf(dist, x)
    end
    return x, logpx
end

# banana distribution
_ϕb(b, x) = [x[1]; x[2] + b * (x[1]^2 - 100); x[3:end]]
function _fb(b, x)
    y = _ϕb(b, x)
    Σ = Diagonal([100; ones(length(y) - 1)])
    return -dot(y, inv(Σ), y) / 2
end
logp_banana(x) = _fb(0.03, x)


text:
module SchumakerSpline

using Plots
using Dates

include("SchumakerFunctions.jl")

export Schumaker_ExtrapolationSchemes, Curve, Linear, Constant
export Schumaker, evaluate, evaluate_integral
include("roots_optima_intercepts.jl")
export find_derivative_spline, find_roots, find_optima, get_crossover_in_interval, get_intersection_points
include("map_to_shape.jl")
export reshape_values
include("splice_splines.jl")
export splice_splines
include("algebra.jl")
export +,-,*,/
include("plotting.jl")
export plot
include("higher_dimensions.jl")
export Schumaker2d
end


text:
export domainplot, coefficientplot, complexplot

## Fun routines

function _maybe_sort(p::AbstractVector{<:Real}, v)
    σ = sortperm(p)
    p[σ], v[σ]
end

_maybe_sort(p, v) = (p,v)

function plotptsvals(f::Fun)
    if isinf(dimension(space(f)))
        f=pad(f,3ncoefficients(f)+50)
    else
        f=pad(f,dimension(space(f)))
    end
    return _maybe_sort(points(f),values(f))
end

function plotptsvals(f::Fun{<:JacobiWeight})
    f=pad(f,3ncoefficients(f)+50)
    s=space(f)
    pts,vals=points(f),values(f)
    # add endpoints so that singularity is viewable
    if s.β ≥ 0
        pts=push!(pts,leftendpoint(domain(f)))
        vals=push!(vals,first(f))
    end
    if s.α ≥ 0
        pts=insert!(pts,1,rightendpoint(domain(f)))
        vals=insert!(vals,1,last(f))
    end

    pts,vals
end

## Recipes

@recipe function f(g::Fun{S,T}) where {S,T<:Real}
    plotptsvals(g)
end
@recipe function f(g::Fun{S,Complex{T}}) where {S,T<:Real}
    x,v=plotptsvals(g)
    x,Vector{T}[real(v),imag(v)]
end

@recipe function f(x::Fun{V,T},y::Fun{S,T}) where {S,V,T<:Real}
    M=3max(ncoefficients(x),ncoefficients(y))+50
    values(pad(x,M)),values(pad(y,M))
end

@recipe function f(x::AbstractVector{T},g::Fun{S,T}) where {S,T<:Real}
    x,g.(x)
end

@recipe function f(x::AbstractVector{T},g::Fun{S,Complex{T}}) where {S,T<:Real}
    v=g.(x)
    x,Vector{T}[real(v),imag(v)]
end

@recipe function f(G::AbstractVector{<:Fun})
    x=Vector{Float64}[]
    v=Vector{Float64}[]
    for g in G
        xx,vv=plotptsvals(g)
        push!(x,xx)
        push!(v,vv)
    end
    _maybe_sort(x,v)
end

@recipe function f(x::AbstractVector{T},G::AbstractVector{F}) where {T<:Real,F<:Fun}
    v=Vector{Float64}[]
    for g in G
        push!(v,g.(x))
    end
    x,v
end

function complexplotvals(f::Fun)
    vals =plotptsvals(f)[2]
    if isa(domain(f),PeriodicDomain)
        real([vals;vals[1]]),imag([vals;vals[1]])
    else
        real(vals),imag(vals)
    end
end

@recipe function f(::Type{Val{:complexplot}},h)
    complexplotvals(h)
end

@recipe function f(dd::Domain)
    x,y = complexplotvals(Fun(identity,dd))
    # @series (x,y)
    # @series begin
    #     primary := false
    #     ∂(dd)
    # end
end

@recipe function f(dd::UnionDomain)
    @series component(dd,1)
    for k=2:ncomponents(dd)
        @series begin
            primary := false
            component(dd,k)
        end
    end
end

@recipe function f(G::AbstractVector{F}) where F<:Domain
    x=Vector{Float64}[]
    v=Vector{Float64}[]
    for g in G
        xx,vv=complexplotvals(g)
        push!(x,xx)
        push!(v,vv)
    end
    x,v
end

# @userplot DomainPlot
#
# @recipe function f(D::DomainPlot)
#     @assert length(D.args)==1
#     g=D.args[1]
#     @assert isa(g,Fun)
#
#     domain(g)
# end
#
#
#
# @userplot CoefficientPlot
#
# @recipe function f(C::CoefficientPlot)
#     @assert length(C.args)==1
#     g=C.args[1]
#     @assert isa(g,Fun)
#
#     seriestype --> :scatter
#     yscale --> :log10
#     markersize := max(round(Int,5 - log10(ncoefficients(g))),1)
#     xlims --> (0,1.01ncoefficients(g))
#     abs(g.coefficients)
# end

@recipe function f(g::Fun{S,T}) where {S<:ArraySpace,T<:Real}
    components(g)
end

@recipe function f(g::Fun{S,T}) where {S<:PiecewiseSpace,T<:Real}
    p=components(g)
    for k=1:length(p)
        @series begin
            primary := (k==1)
            p[k]
        end
    end
end

@recipe function f(g::Fun{S,T}) where {S<:DiracSpace,T<:Real}
    pts=space(g).points
    n=length(pts)
    ws=pad(g.coefficients,length(pts))

    xlims --> (minimum(pts)-1.,maximum(pts)+1.)

    @series begin
        primary --> true
        fill(1.0,2)*pts[1],[0,1]*ws[1]
    end

    if length(ws) > 1
        @series begin
            primary := false
            fill(1.0,2)*pts[2:end]',[0,1]*ws[2:end]'
        end
    end

    @series begin
        primary := false
        linestyle := :dot
        fill(1.0,2)*pts',[1,2]*ws'
    end
end

@recipe function f(g::Fun{S,T}) where {S<:PointSpace,T<:Real}
    pts=space(g).points
    n=length(pts)
    ws=pad(g.coefficients,length(pts))

    xlims --> (minimum(pts)-1.,maximum(pts)+1.)

    @series begin
        primary --> true
        fill(1.0,2)*pts[1],[0,1]*ws[1]
    end
    if length(ws) > 1
        @series begin
            primary := false
            fill(1.0,2)*pts[2:end]',[0,1]*ws[2:end]'
        end
    end
end

@recipe function f(g::Fun{S,T}) where {S<:HeavisideSpace,T<:Real}
    pts=domain(g).points
    n=length(pts)
    ws=pad(g.coefficients,dimension(space(g)))

    lnsx=Vector{Float64}()
    lnsy=Vector{Float64}()
    dtsx=Vector{Float64}()
    dtsy=Vector{Float64}()
    for k=1:n-1
        push!(lnsx,pts[k])
        push!(lnsy,ws[k])
        push!(lnsx,pts[k+1])
        push!(lnsy,ws[k])

        if k != n-1
            # dotted line, with NaN to separate
            push!(dtsx,pts[k+1])
            push!(dtsx,pts[k+1])
            push!(dtsx,pts[k+1])

            push!(dtsy,ws[k])
            push!(dtsy,ws[k+1])
            push!(dtsy,NaN)

            # extra point for NaN
            push!(lnsx,pts[k+1])
            push!(lnsy,NaN)
        end
    end

    @series begin
        primary --> true
        lnsx,lnsy
    end

    @series begin
       primary := false
       linestyle := :dot
       fillrange := nothing
       dtsx,dtsy
    end
end

###
# Multivariate
###

@recipe function f(g::ProductFun{S,V,SV}) where {S<:UnivariateSpace,
                    V<:UnivariateSpace,
        SV<:TensorSpace}
    g=chop(g,10e-10)
    g=pad(g,max(size(g,1),20),max(size(g,2),20))
    vals=values(g)

    seriestype --> :surface

    if norm(imag(vals),Inf)>10e-9
        @warn "Imaginary part is non-neglible.  Only plotting real part."
    end

    # sort the points
    x = points(factor(space(g),1),size(vals,1))
    y = points(factor(space(g),2),size(vals,2))
    px = sortperm(x)
    py = sortperm(y)

    x[px],y[py],transpose(real(vals))[py,px]
end

@recipe function f(g::LowRankFun{S,V,SV}) where {S<:UnivariateSpace,
                    V<:UnivariateSpace,
        SV<:TensorSpace}
    g=chop(g,10e-10)
    vals=values(g)

    seriestype --> :surface

    if norm(imag(vals),Inf)>10e-9
        @warn "Imaginary part is non-neglible.  Only plotting real part."
    end

    # sort the points
    x = points(factor(space(g),1),size(vals,1))
    y = points(factor(space(g),2),size(vals,2))
    px = sortperm(x)
    py = sortperm(y)

    x[px],y[py],transpose(real(vals))[py,px]
end

@recipe function f(x::AbstractVector,y::AbstractVector,g::MultivariateFun)
    seriestype --> :surface
    x,y,real(g.(x,transpose(y)))
end

@recipe function f(g::Fun{TS,T}) where {TS<:BivariateSpace,T<:Real}
    g = pad(g,ncoefficients(g)+100)
    pts = points(g)
    seriestype --> :surface
    first.(pts),last.(pts),values(g)
end

@recipe function f(g::Fun{TS,T}) where {TS<:AbstractProductSpace,T<:Real}
    ProductFun(g)
end

@recipe function f(g::Fun{TS,T}) where {TS<:AbstractProductSpace,T<:Complex}
    ProductFun(g)
end

@recipe function f(x::AbstractVector,y::AbstractVector,g::Fun{TS,T}) where {TS<:AbstractProductSpace,T<:Real}
    x,y,ProductFun(g)
end

@recipe function f(x::AbstractVector,y::AbstractVector,g::Fun{TS,T}) where {TS<:AbstractProductSpace,T<:Complex}
    x,y,ProductFun(g)
end


text:
function get_all_interactions(n::MangalNetwork)
    page_size = 200
    network_interactions = MangalInteraction[]
    interactions_to_get = count(MangalInteraction, "network_id" => n.id)
    pages_to_do = convert(Int64, ceil(interactions_to_get/page_size))
    for page in 1:pages_to_do
        paging_query = ["page" => "$(page-1)", "count" => "page_size"]
        append!(network_interactions, Mangal.interactions(n, paging_query...))
    end
    return network_interactions
end

import Base.convert

"""
    convert(::Type{UnipartiteNetwork}, n::MangalNetwork; resolution::Type=MangalNode)

Returns a `UnipartiteNetwork` object representation of the `MangalNetwork`
passed as its first argument. The optional keyword `resolution` (can be
`MangalNode` or `MangalReferenceTaxon`) is determined to use which level of
aggregation should be used. The default (`MangalNode`) is raw data, and
`MangalReferenceTaxon` is the cleaned version.
"""
function convert(::Type{UnipartiteNetwork}, n::MangalNetwork; resolution::Type=MangalNode)
    resolution ∈ [MangalNode, MangalReferenceTaxon] || throw(ArgumentError("The resolution argument can only be MangalNode or MangalReferenceTaxon - you used $(resolution)"))

    network_interactions = get_all_interactions(n)
    all_object_nodes = resolution[]

    for i in network_interactions
        push!(all_object_nodes, resolution == MangalNode ? i.from : i.from.taxon)
        push!(all_object_nodes, resolution == MangalNode ? i.to : i.to.taxon)
    end

    object_nodes = unique(all_object_nodes)
    A = zeros(Bool, (length(object_nodes),length(object_nodes)))
    N = UnipartiteNetwork(A, object_nodes)
    for i in network_interactions
        if i.strength != 0
            t_from = resolution == MangalReferenceTaxon ? i.from.taxon : i.from
            t_to   = resolution == MangalReferenceTaxon ? i.to.taxon : i.to
            N[t_from, t_to] = true
        end
    end

    return N
end


text:
"""
    InverseWishart(ν, Ψ)
```julia
ν::Real           degrees of freedom (greater than p - 1)
Ψ::AbstractPDMat  p x p scale matrix
```
The [inverse Wishart distribution](http://en.wikipedia.org/wiki/Inverse-Wishart_distribution)
generalizes the inverse gamma distribution to ``p\\times p`` real, positive definite
matrices ``\\boldsymbol{\\Sigma}``. If ``\\boldsymbol{\\Sigma}\\sim IW_p(\\nu,\\boldsymbol{\\Psi})``,
then its probability density function is

```math
f(\\boldsymbol{\\Sigma}; \\nu,\\boldsymbol{\\Psi}) =
\\frac{\\left|\\boldsymbol{\\Psi}\\right|^{\\nu/2}}{2^{\\nu p/2}\\Gamma_p(\\frac{\\nu}{2})} \\left|\\boldsymbol{\\Sigma}\\right|^{-(\\nu+p+1)/2} e^{-\\frac{1}{2}\\operatorname{tr}(\\boldsymbol{\\Psi}\\boldsymbol{\\Sigma}^{-1})}.
```

``\\mathbf{H}\\sim W_p(\\nu, \\mathbf{S})`` if and only if ``\\mathbf{H}^{-1}\\sim IW_p(\\nu, \\mathbf{S}^{-1})``.
"""
struct InverseWishart{T<:Real, ST<:AbstractPDMat} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    Ψ::ST     # scale matrix
    c0::T     # log of normalizing constant
end

#### Constructors

function InverseWishart(df::T, Ψ::AbstractPDMat{T}) where T<:Real
    p = dim(Ψ)
    df > p - 1 || error("df should be greater than dim - 1.")
    c0 = _invwishart_c0(df, Ψ)
    R = Base.promote_eltype(T, c0)
    prom_Ψ = convert(AbstractArray{R}, Ψ)
    InverseWishart{R, typeof(prom_Ψ)}(R(df), prom_Ψ, R(c0))
end

function InverseWishart(df::Real, Ψ::AbstractPDMat)
    T = Base.promote_eltype(df, Ψ)
    InverseWishart(T(df), convert(AbstractArray{T}, Ψ))
end

InverseWishart(df::Real, Ψ::Matrix) = InverseWishart(df, PDMat(Ψ))

InverseWishart(df::Real, Ψ::Cholesky) = InverseWishart(df, PDMat(Ψ))

function _invwishart_c0(df::Real, Ψ::AbstractPDMat)
    h_df = df / 2
    p = dim(Ψ)
    h_df * (p * typeof(df)(logtwo) - logdet(Ψ)) + logmvgamma(p, h_df)
end

#### Properties

insupport(::Type{InverseWishart}, X::Matrix) = isposdef(X)
insupport(d::InverseWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

dim(d::InverseWishart) = dim(d.Ψ)
size(d::InverseWishart) = (p = dim(d); (p, p))
size(d::InverseWishart, i) = size(d)[i]
rank(d::InverseWishart) = dim(d)
params(d::InverseWishart) = (d.df, d.Ψ, d.c0)
@inline partype(d::InverseWishart{T}) where {T<:Real} = T

### Conversion
function convert(::Type{InverseWishart{T}}, d::InverseWishart) where T<:Real
    P = convert(AbstractArray{T}, d.Ψ)
    InverseWishart{T, typeof(P)}(T(d.df), P, T(d.c0))
end
function convert(::Type{InverseWishart{T}}, df, Ψ::AbstractPDMat, c0) where T<:Real
    P = convert(AbstractArray{T}, Ψ)
    InverseWishart{T, typeof(P)}(T(df), P, T(c0))
end

#### Show

show(io::IO, d::InverseWishart) = show_multline(io, d, [(:df, d.df), (:Ψ, Matrix(d.Ψ))])

#### Statistics

function mean(d::InverseWishart)
    df = d.df
    p = dim(d)
    r = df - (p + 1)
    if r > 0.0
        return Matrix(d.Ψ) * (1.0 / r)
    else
        error("mean only defined for df > p + 1")
    end
end

mode(d::InverseWishart) = d.Ψ * inv(d.df + dim(d) + 1.0)

#  https://en.wikipedia.org/wiki/Inverse-Wishart_distribution#Moments
function cov(d::InverseWishart, i::Integer, j::Integer, k::Integer, l::Integer)
    p, ν, Ψ = (dim(d), d.df, Matrix(d.Ψ))
    ν > p + 3 || throw(ArgumentError("cov only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*(2Ψ[i,j]*Ψ[k,l] + (ν-p-1)*(Ψ[i,k]*Ψ[j,l] + Ψ[i,l]*Ψ[k,j]))
end

function var(d::InverseWishart, i::Integer, j::Integer)
    p, ν, Ψ = (dim(d), d.df, Matrix(d.Ψ))
    ν > p + 3 || throw(ArgumentError("var only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*(ν - p + 1)*Ψ[i,j]^2 + (ν - p - 1)*Ψ[i,i]*Ψ[j,j]
end

#### Evaluation

function _logpdf(d::InverseWishart, X::AbstractMatrix)
    p = dim(d)
    df = d.df
    Xcf = cholesky(X)
    # we use the fact: tr(Ψ * inv(X)) = tr(inv(X) * Ψ) = tr(X \ Ψ)
    Ψ = Matrix(d.Ψ)
    -0.5 * ((df + p + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) - d.c0
end

#### Sampling

_rand!(rng::AbstractRNG, d::InverseWishart, A::AbstractMatrix) =
    (A .= inv(cholesky!(_rand!(rng, Wishart(d.df, inv(d.Ψ)), A))))


text:
# some transitive groups of degree 16

const _16T1 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 4, 1, 6, 3, 8, 5, 10, 7, 12, 9, 14, 11, 16, 13, 15,
                       3, 1, 5, 2, 7, 4, 9, 6, 11, 8, 13, 10, 15, 12, 16, 14,
                       4, 6, 2, 8, 1, 10, 3, 12, 5, 14, 7, 16, 9, 15, 11, 13,
                       5, 3, 7, 1, 9, 2, 11, 4, 13, 6, 15, 8, 16, 10, 14, 12,
                       6, 8, 4, 10, 2, 12, 1, 14, 3, 16, 5, 15, 7, 13, 9, 11,
                       7, 5, 9, 3, 11, 1, 13, 2, 15, 4, 16, 6, 14, 8, 12, 10,
                       8, 10, 6, 12, 4, 14, 2, 16, 1, 15, 3, 13, 5, 11, 7, 9,
                       9, 7, 11, 5, 13, 3, 15, 1, 16, 2, 14, 4, 12, 6, 10, 8,
                       10, 12, 8, 14, 6, 16, 4, 15, 2, 13, 1, 11, 3, 9, 5, 7,
                       11, 9, 13, 7, 15, 5, 16, 3, 14, 1, 12, 2, 10, 4, 8, 6,
                       12, 14, 10, 16, 8, 15, 6, 13, 4, 11, 2, 9, 1, 7, 3, 5,
                       13, 11, 15, 9, 16, 7, 14, 5, 12, 3, 10, 1, 8, 2, 6, 4,
                       14, 16, 12, 15, 10, 13, 8, 11, 6, 9, 4, 7, 2, 5, 1, 3,
                       15, 13, 16, 11, 14, 9, 12, 7, 10, 5, 8, 3, 6, 1, 4, 2,
                       16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 ],
                      (16, 16))

const _16T2 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 1, 6, 7, 8, 3, 4, 5, 12, 13, 14, 9, 10, 11, 16, 15,
                       3, 6, 1, 9, 10, 2, 12, 13, 4, 5, 15, 7, 8, 16, 11, 14,
                       4, 7, 9, 11, 1, 12, 14, 2, 15, 3, 5, 16, 6, 8, 10, 13,
                       5, 8, 10, 1, 11, 13, 2, 14, 3, 15, 4, 6, 16, 7, 9, 12,
                       6, 3, 2, 12, 13, 1, 9, 10, 7, 8, 16, 4, 5, 15, 14, 11,
                       7, 4, 12, 14, 2, 9, 11, 1, 16, 6, 8, 15, 3, 5, 13, 10,
                       8, 5, 13, 2, 14, 10, 1, 11, 6, 16, 7, 3, 15, 4, 12, 9,
                       9, 12, 4, 15, 3, 7, 16, 6, 11, 1, 10, 14, 2, 13, 5, 8,
                       10, 13, 5, 3, 15, 8, 6, 16, 1, 11, 9, 2, 14, 12, 4, 7,
                       11, 14, 15, 5, 4, 16, 8, 7, 10, 9, 1, 13, 12, 2, 3, 6,
                       12, 9, 7, 16, 6, 4, 15, 3, 14, 2, 13, 11, 1, 10, 8, 5,
                       13, 10, 8, 6, 16, 5, 3, 15, 2, 14, 12, 1, 11, 9, 7, 4,
                       14, 11, 16, 8, 7, 15, 5, 4, 13, 12, 2, 10, 9, 1, 6, 3,
                       15, 16, 11, 10, 9, 14, 13, 12, 5, 4, 3, 8, 7, 6, 1, 2,
                       16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 ],
                       (16, 16))

const _16T3 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 1, 6, 7, 8, 3, 4, 5, 12, 13, 14, 9, 10, 11, 16, 15,
                       3, 6, 1, 9, 10, 2, 12, 13, 4, 5, 15, 7, 8, 16, 11, 14,
                       4, 7, 9, 1, 11, 12, 2, 14, 3, 15, 5, 6, 16, 8, 10, 13,
                       5, 8, 10, 11, 1, 13, 14, 2, 15, 3, 4, 16, 6, 7, 9, 12,
                       6, 3, 2, 12, 13, 1, 9, 10, 7, 8, 16, 4, 5, 15, 14, 11,
                       7, 4, 12, 2, 14, 9, 1, 11, 6, 16, 8, 3, 15, 5, 13, 10,
                       8, 5, 13, 14, 2, 10, 11, 1, 16, 6, 7, 15, 3, 4, 12, 9,
                       9, 12, 4, 3, 15, 7, 6, 16, 1, 11, 10, 2, 14, 13, 5, 8,
                       10, 13, 5, 15, 3, 8, 16, 6, 11, 1, 9, 14, 2, 12, 4, 7,
                       11, 14, 15, 5, 4, 16, 8, 7, 10, 9, 1, 13, 12, 2, 3, 6,
                       12, 9, 7, 6, 16, 4, 3, 15, 2, 14, 13, 1, 11, 10, 8, 5,
                       13, 10, 8, 16, 6, 5, 15, 3, 14, 2, 12, 11, 1, 9, 7, 4,
                       14, 11, 16, 8, 7, 15, 5, 4, 13, 12, 2, 10, 9, 1, 6, 3,
                       15, 16, 11, 10, 9, 14, 13, 12, 5, 4, 3, 8, 7, 6, 1, 2,
                       16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 ],
                       (16, 16))

const _16T4 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 6, 7, 1, 8, 4, 12, 13, 14, 3, 5, 10, 11, 16, 9, 15,
                       3, 7, 9, 10, 1, 12, 14, 2, 5, 15, 4, 16, 6, 8, 11, 13,
                       4, 1, 10, 6, 11, 2, 3, 5, 15, 12, 13, 7, 8, 9, 16, 14,
                       5, 8, 1, 11, 9, 13, 2, 14, 3, 4, 15, 6, 16, 7, 10, 12,
                       6, 4, 12, 2, 13, 1, 10, 11, 16, 7, 8, 3, 5, 15, 14, 9,
                       7, 12, 14, 3, 2, 10, 16, 6, 8, 9, 1, 15, 4, 13, 5, 11,
                       8, 13, 2, 5, 14, 11, 6, 16, 7, 1, 9, 4, 15, 12, 3, 10,
                       9, 14, 5, 15, 3, 16, 8, 7, 1, 11, 10, 13, 12, 2, 4, 6,
                       10, 3, 15, 12, 4, 7, 9, 1, 11, 16, 6, 14, 2, 5, 13, 8,
                       11, 5, 4, 13, 15, 8, 1, 9, 10, 6, 16, 2, 14, 3, 12, 7,
                       12, 10, 16, 7, 6, 3, 15, 4, 13, 14, 2, 9, 1, 11, 8, 5,
                       13, 11, 6, 8, 16, 5, 4, 15, 12, 2, 14, 1, 9, 10, 7, 3,
                       14, 16, 8, 9, 7, 15, 13, 12, 2, 5, 3, 11, 10, 6, 1, 4,
                       15, 9, 11, 16, 10, 14, 5, 3, 4, 13, 12, 8, 7, 1, 6, 2,
                       16, 15, 13, 14, 12, 9, 11, 10, 6, 8, 7, 5, 3, 4, 2, 1 ],
                       (16, 16))

const _16T5 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 6, 7, 1, 8, 11, 5, 12, 4, 3, 15, 16, 10, 9, 14, 13,
                       3, 7, 9, 10, 1, 5, 4, 2, 13, 14, 8, 6, 15, 16, 12, 11,
                       4, 1, 10, 9, 7, 2, 3, 5, 14, 13, 6, 8, 16, 15, 11, 12,
                       5, 8, 1, 7, 6, 12, 2, 11, 3, 4, 16, 15, 9, 10, 13, 14,
                       6, 11, 5, 2, 12, 15, 8, 16, 1, 7, 14, 13, 3, 4, 9, 10,
                       7, 5, 4, 3, 2, 8, 1, 6, 10, 9, 12, 11, 14, 13, 16, 15,
                       8, 12, 2, 5, 11, 16, 6, 15, 7, 1, 13, 14, 4, 3, 10, 9,
                       9, 4, 13, 14, 3, 1, 10, 7, 15, 16, 2, 5, 12, 11, 6, 8,
                       10, 3, 14, 13, 4, 7, 9, 1, 16, 15, 5, 2, 11, 12, 8, 6,
                       11, 15, 8, 6, 16, 14, 12, 13, 2, 5, 9, 10, 7, 1, 4, 3,
                       12, 16, 6, 8, 15, 13, 11, 14, 5, 2, 10, 9, 1, 7, 3, 4,
                       13, 10, 15, 16, 9, 3, 14, 4, 12, 11, 7, 1, 6, 8, 5, 2,
                       14, 9, 16, 15, 10, 4, 13, 3, 11, 12, 1, 7, 8, 6, 2, 5,
                       15, 14, 12, 11, 13, 9, 16, 10, 6, 8, 4, 3, 5, 2, 1, 7,
                       16, 13, 11, 12, 14, 10, 15, 9, 8, 6, 3, 4, 2, 5, 7, 1 ],
                       (16, 16))

const _16T8 = reshape([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                       2, 6, 9, 1, 10, 4, 11, 12, 5, 3, 8, 7, 15, 16, 14, 13,
                       3, 7, 6, 8, 1, 5, 4, 2, 13, 14, 15, 16, 10, 9, 12, 11,
                       4, 1, 10, 6, 9, 2, 12, 11, 3, 5, 7, 8, 16, 15, 13, 14,
                       5, 8, 1, 7, 6, 3, 2, 4, 14, 13, 16, 15, 9, 10, 11, 12,
                       6, 4, 5, 2, 3, 1, 8, 7, 10, 9, 12, 11, 14, 13, 16, 15,
                       7, 5, 13, 3, 14, 8, 15, 16, 1, 6, 2, 4, 12, 11, 9, 10,
                       8, 3, 14, 5, 13, 7, 16, 15, 6, 1, 4, 2, 11, 12, 10, 9,
                       9, 11, 4, 12, 2, 10, 1, 6, 15, 16, 14, 13, 3, 5, 7, 8,
                       10, 12, 2, 11, 4, 9, 6, 1, 16, 15, 13, 14, 5, 3, 8, 7,
                       11, 10, 15, 9, 16, 12, 14, 13, 2, 4, 6, 1, 7, 8, 5, 3,
                       12, 9, 16, 10, 15, 11, 13, 14, 4, 2, 1, 6, 8, 7, 3, 5,
                       13, 15, 8, 16, 7, 14, 3, 5, 12, 11, 9, 10, 6, 1, 4, 2,
                       14, 16, 7, 15, 8, 13, 5, 3, 11, 12, 10, 9, 1, 6, 2, 4,
                       15, 14, 12, 13, 11, 16, 9, 10, 7, 8, 5, 3, 4, 2, 1, 6,
                       16, 13, 11, 14, 12, 15, 10, 9, 8, 7, 3, 5, 2, 4, 6, 1 ],
                       (16, 16))


text:
function can_use_greedy(coins::Array{Int, 1})::Bool
	divisible::Bool = true
	for i in 1 : length(coins) - 1
		divisible &= (coins[i] % coins[i + 1]) == 0
	end

	return divisible
end

function min_coins_greedy(coins::Array{Int, 1}, value::Int)::Int
	amt_coins::Int = 0

	for i in coins
		amt_coins += fld(value, i)
		value %= i
	end

	return amt_coins
end

function min_coins_dynamic(coins::Array{Int, 1}, value::Int)::Int
	amt_coins = fld(value, coins[length(coins)])

	for i in length(coins) : -1 : 1
		temp = min_coins_greedy(coins[i : length(coins)], value)
		if temp < amt_coins
			amt_coins = temp
		end
	end

	return amt_coins
end


text:
## Julia 1.1.0

# println(" -> $()")

# dot
println("\n\n---------------- dot ----------------")
println("[1, 2, 3] .^ 3 -> $([1, 2, 3] .^ 3)")

# NaN
println("\n\n---------------- NaN ----------------")
println("NaN == NaN -> $(NaN == NaN)")
println("NaN != NaN -> $(NaN != NaN)")
println("NaN < NaN -> $(NaN < NaN)")
println("NaN > NaN -> $(NaN > NaN)")

# Comparisons
x = 1
println("\n\n---------------- Comparisons ----------------")
println("isequal(NaN, NaN32) -> $(isequal(NaN, NaN32))")
for T in [isfinite, isinf, isnan]
    println("$(T)($(x)) -> $(T(x))")
end

v(x) = (println(x); x)
println("v(1) < v(2) <= v(3) -> $(v(1) < v(2) <= v(3))")

# Rounding
x = 1.1
println("\n\n---------------- Rounding ----------------")
println("round($(x)) -> $(round(x))")
println("round(UInt32, $(x)) -> $(round(UInt32, x))")
println("floor($(x)) -> $(floor(x))")
println("floor(UInt32, $(x)) -> $(floor(UInt32, x))")
println("ceil($(x)) -> $(ceil(x))")
println("ceil(UInt32, $(x)) -> $(ceil(UInt32, x))")
println("trunc($(x)) -> $(trunc(x))")
println("trunc(UInt32, $(x)) -> $(trunc(UInt32, x))")

# Division
x = 17
y = 13
println("\n\n---------------- Division ----------------")
for T in [div, fld, cld, rem, mod, mod1, fldmod, gcd, lcm]
    println("$(T)($(x), $(y)) -> $(T(x, y))")
end
println("mod2pi($(x)) -> $(mod2pi(x))")

# Power, Logs, Roots
x = 17.0
y = 13
println("\n\n---------------- Power, Logs, Roots ----------------")
for T in [sqrt, cbrt, exp, expm1, log, log2, log1p, exponent, significand]
    println("$(T)($(x)) -> $(T(x))")
end
println("hypot($(x), $(y)) -> $(hypot(x, y))")
println("ldexp($(x), 3) -> $(ldexp(x, 3))")
println("log(3.14, $(x)) -> $(log(3.14, x))")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
# Pkg.clone(pwd())
# Pkg.build("Plots")

Pkg.clone("ImageMagick")
Pkg.build("ImageMagick")

Pkg.clone("GR")
Pkg.build("GR")

# Pkg.clone("https://github.com/JuliaPlots/RecipesBase.jl.git")
# Pkg.clone("https://github.com/tbreloff/ExamplePlots.jl.git");
Pkg.clone("https://github.com/JuliaPlots/PlotReferenceImages.jl.git")

# Pkg.clone("https://github.com/JunoLab/Blink.jl.git")
# Pkg.build("Blink")
# import Blink
# Blink.AtomShell.install()
# Pkg.clone("https://github.com/spencerlyon2/PlotlyJS.jl.git")

Pkg.checkout("RecipesBase")
Pkg.clone("VisualRegressionTests")

ENV["PYTHON"] = ""
Pkg.add("PyPlot")
Pkg.build("PyPlot")

Pkg.test("Plots"; coverage=false)


text:
using Random
using IntervalSets

import JuliaRL

module ContMountainCarConst
const vel_limit = (-0.07, 0.07)
const pos_limit = (-1.2, 0.5)
const pos_initial_range = (-0.6, 0.4)

const ActionRange = -1..1

end

"""
ContMountainCar(pos=0.0, vel=0.0, normalized::Bool=false)
ContMountainCar(rng::AbstractRNG, normalized::Bool=false)

    The mountain car environment.

"""
mutable struct ContMountainCar <: RLCore.AbstractEnvironment
    pos::Float64
    vel::Float64
    actions::Interval{:closed, :closed, Float64}
    normalized::Bool
    function ContMountainCar(pos=0.0, vel=0.0, normalized::Bool=false)
        mcc = ContMountainCarConst
        @boundscheck (pos >= mcc.pos_limit[1] && pos <= mcc.pos_limit[2])
        @boundscheck (vel >= mcc.vel_limit[1] && vel <= mcc.vel_limit[2])
        new(pos, vel, ContMountainCarCons.ActionRange, normalized)
    end
end

RLCore.get_actions(env::ContMountainCar) = env.actions
valid_action(env::ContMountainCar, action) = action ∈ env.actions

function RLCore.reset!(env::ContMountainCar, rng::AbstractRNG; kwargs...)
    env.pos = (rand(rng)*(ContMountainCarConst.pos_initial_range[2]
                          - ContMountainCarConst.pos_initial_range[1])
               + ContMountainCarConst.pos_initial_range[1])
    env.vel = 0.0
end

function RLCore.reset!(env::ContMountainCar,
                        start_state::T;
                        kwargs...) where {T<:AbstractArray}
    if env.normalized
        env.pos = start_state[1]
        env.vel = start_state[2]
    else
        pos_limit = ContMountainCarConst.pos_limit
        vel_limit = ContMountainCarConst.vel_limit
        env.pos = (start_state[1]*(pos_limit[2] - pos_limit[1])) + pos_limit[1]
        env.vel = (start_state[2]*(vel_limit[2] - vel_limit[1])) + vel_limit[1]
    end
end

function RLCore.environment_step!(env::ContMountainCar,
                                   action;
                                   rng=Random.GLOBAL_RNG, kwargs...)

    @boundscheck valid_action(env, action)
    env.vel =
        clamp(env.vel + (action)*0.001 - 0.0025*cos(3*env.pos),
              ContMountainCarConst.vel_limit...)
    env.pos =
        clamp(env.pos + env.vel,
              ContMountainCarConst.pos_limit...)
end

function RLCore.get_reward(env::ContMountainCar) # -> determines if the agent_state is terminal
    if env.pos >= ContMountainCarConst.pos_limit[2]
        return 0
    end
    return -1
end

function RLCore.is_terminal(env::ContMountainCar) # -> determines if the agent_state is terminal
    return env.pos >= ContMountainCarConst.pos_limit[2]
end

function RLCore.get_state(env::ContMountainCar)
    if env.normalized
        return get_normalized_state(env)
    else
        return [env.pos, env.vel]
    end
end

function get_normalized_state(env::ContMountainCar)
    pos_limit = ContMountainCarConst.pos_limit
    vel_limit = ContMountainCarConst.vel_limit
    return [(env.pos - pos_limit[1])/(pos_limit[2] - pos_limit[1]),
            (env.vel - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end


text:
# Statistical Models

abstract StatisticalModel

coef(obj::StatisticalModel) = error("coef is not defined for $(typeof(obj)).")
coeftable(obj::StatisticalModel) = error("coeftable is not defined for $(typeof(obj)).")
confint(obj::StatisticalModel) = error("coefint is not defined for $(typeof(obj)).")
deviance(obj::StatisticalModel) = error("deviance is not defined for $(typeof(obj)).")
loglikelihood(obj::StatisticalModel) = error("loglikelihood is not defined for $(typeof(obj)).")
nobs(obj::StatisticalModel) = size(model_response(obj), 1)
stderr(obj::StatisticalModel) = sqrt(diag(vcov(obj)))
vcov(obj::StatisticalModel) = error("vcov is not defined for $(typeof(obj)).")
fit(obj::StatisticalModel, data...) = error("fit is not defined for $(typeof(obj)).")

abstract RegressionModel <: StatisticalModel

fitted(obj::RegressionModel) = error("fitted is not defined for $(typeof(obj)).")
model_response(obj::RegressionModel) = error("model_response is not defined for $(typeof(obj)).")
residuals(obj::RegressionModel) = error("residuals is not defined for $(typeof(obj)).")
predict(obj::RegressionModel) = error("predict is not defined for $(typeof(obj)).")
predict!(obj::RegressionModel) = error("predict! is not defined for $(typeof(obj)).")

## coefficient tables with specialized show method

## Nms are the coefficient names, corresponding to rows in the table
type CoefTable
    mat::Matrix
    colnms::Vector
    rownms::Vector
    pvalcol::Integer
    function CoefTable(mat::Matrix,colnms::Vector,rownms::Vector,pvalcol::Int=0)
        nr,nc = size(mat)
        0 <= pvalcol <= nc || error("pvalcol = $pvalcol should be in 0,...,$nc]")
        length(colnms) in [0,nc] || error("colnms should have length 0 or $nc")
        length(rownms) in [0,nr] || error("rownms should have length 0 or $nr")
        new(mat,colnms,rownms,pvalcol)
    end
end

## format numbers in the p-value column
function format_pvc(pv::Number)
    0. <= pv <= 1. || error("p-values must be in [0.,1.]")
    if pv >= 1e-4
        return @sprintf("%.4f", pv)
    else
        return @sprintf("<1e%2.2d", iceil(max(nextfloat(log10(pv)), -99)))
    end
end

function show(io::IO, ct::CoefTable)
    mat = ct.mat; nr,nc = size(mat); rownms = ct.rownms; colnms = ct.colnms; pvc = ct.pvalcol
    if length(rownms) == 0
        rownms = [lpad("[$i]",ifloor(log10(nr))+3)::String for i in 1:nr]
    end
    rnwidth = max(4,maximum([length(nm) for nm in rownms]) + 1)
    rownms = [rpad(nm,rnwidth) for nm in rownms]
    widths = [length(cn)::Int for cn in colnms]
    str = [sprint(showcompact,mat[i,j]) for i in 1:nr, j in 1:nc]
    if pvc != 0                         # format the p-values column
        for i in 1:nr
            str[i,pvc] = format_pvc(mat[i,pvc])
        end
    end
    for j in 1:nc
        for i in 1:nr
            lij = length(str[i,j])
            if lij > widths[j]
                widths[j] = lij
            end
        end
    end
    widths .+= 1
    println(io," " ^ rnwidth *
            join([lpad(string(colnms[i]), widths[i]) for i = 1:nc], ""))
    for i = 1:nr
        print(io, rownms[i])
        for j in 1:nc
            print(io, lpad(str[i,j],widths[j]))
        end
        println(io)
    end
end


text:
# Mad Libs in Julia

function play_mad_libs()

    print("Enter a verb (action): ")
    verb1 = readline()

    print("Enter an adjective (descriptive word): ")
    adj1 = readline()

    print("Enter a noun (person place or thing): ")
    noun1 = readline()

    print("Enter another noun (person place or thing): ")
    noun2 = readline()

    print("Enter a catchphrase (something like 'hands up!'): ")
    phrase1 = readline()

    base_sentence = "John $verb1 down the street one night, playing with his $adj1 $noun1. When all of a /
    sudden, a $noun2 jumped out at him and said $phrase1"
    print("\n\n", base_sentence)
end

play_mad_libs()


text:
export RandomShuffle
"""
    RandomShuffle() <: Surrogate

A random constrained surrogate, generated by shifting values around.

Random shuffle surrogates preserve the mean, variance and amplitude
distribution of the original signal. Properties not preserved are *any
temporal information*, such as the power spectrum and hence linear
correlations.

The null hypothesis this method can test for is whether the data
are uncorrelated noise, possibly measured via a nonlinear function.
Specifically, random shuffle surrogate can test
the null hypothesis that the original signal is produced by independent and
identically distributed random variables[^Theiler1991, ^Lancaster2018].

*Beware: random shuffle surrogates do not cover the case of correlated noise*[^Lancaster2018].

[^Theiler1991]: J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
"""
struct RandomShuffle <: Surrogate end

function surrogenerator(x::AbstractVector, rf::RandomShuffle)
    return SurrogateGenerator(rf, x, nothing)
end

function (rf::SurrogateGenerator{<:RandomShuffle})()
    n = length(rf.x)
    rf.x[sample(1:n, n, replace = false)]
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
llvm_const{T<:Signed}(::Type{T}, r::Integer) = LLVM.ConstInt(TypeRef(T), r, true)
llvm_const{T<:Unsigned}(::Type{T}, r::Integer) = LLVM.ConstInt(TypeRef(T), r, false)
llvm_const{T<:Real}(::Type{T}, r::Real) = LLVM.ConstReal(TypeRef(T), r)

function intrinsic(mod::ModuleRef, name::AbstractString, rt::Type, params::Type...)
  fn = LLVM.GetNamedFunction(mod, name)
  if fn.ptr != C_NULL return fn end
  fn_type = LLVM.FunctionType(TypeRef(rt), collect(TypeRef, params))
  return LLVM.AddFunction(mod, name, fn_type)
end

function build_grid_stride_loop(mod::ModuleRef, name::AbstractString, eltype::Type, outs::Vector{ANode}, ins::Vector{ANode})
  pt = LLVM.PointerType(LLVM.TypeRef(eltype), 0)
  fun = LLVM.AddFunction(mod, name, LLVM.FunctionType(TypeRef(Void), [TypeRef(Int32); fill(pt, length(outs)+length(ins))]))
  LLVM.AddNamedMetadataOperand(mod, "nvvm.annotations",
    LLVM.MDNode(ValueRef[fun, LLVM.MDString("kernel"), llvm_const(UInt32, 1)]))

  sorted = Seep.toposort(outs, ins)

  limit = LLVM.GetFirstParam(fun)
  LLVM.SetValueName(limit, "limit")

  inptrs = Dict{ANode,ValueRef}()
  outptrs = Dict{ANode,ValueRef}()

  param = LLVM.GetNextParam(limit)
  for i in outs
    outptrs[i] = param
    LLVM.SetValueName(param, string(i.name, "_out"))
    param = LLVM.GetNextParam(param)
  end

  for i in ins
    inptrs[i] = param
    LLVM.SetValueName(param, string(i.name, "_in"))
    param = LLVM.GetNextParam(param)
  end

  builder = LLVM.CreateBuilder()
  entry_block = LLVM.AppendBasicBlock(fun, "entry")
  test_block = LLVM.AppendBasicBlock(fun, "test")
  loop_block = LLVM.AppendBasicBlock(fun, "loop")
  exit_block = LLVM.AppendBasicBlock(fun, "exit")

  # entry:
  LLVM.PositionBuilderAtEnd(builder, entry_block)
  tid = LLVM.BuildCall(builder, intrinsic(mod, "llvm.nvvm.read.ptx.sreg.tid.x", Int32), ValueRef[], "tid")
  ntid = LLVM.BuildCall(builder, intrinsic(mod, "llvm.nvvm.read.ptx.sreg.ntid.x", Int32), ValueRef[], "ntid")
  ctaid = LLVM.BuildCall(builder, intrinsic(mod, "llvm.nvvm.read.ptx.sreg.ctaid.x", Int32), ValueRef[], "ctaid")
  nctaid = LLVM.BuildCall(builder, intrinsic(mod, "llvm.nvvm.read.ptx.sreg.nctaid.x", Int32), ValueRef[], "nctaid")
  first_index = LLVM.BuildAdd(builder, LLVM.BuildMul(builder, ctaid, ntid, ""), tid, "first_index")
  stride = LLVM.BuildMul(builder, ntid, nctaid, "stride")
  LLVM.BuildBr(builder, test_block)

  # test:
  LLVM.PositionBuilderAtEnd(builder, test_block)
  index = LLVM.BuildPhi(builder, TypeRef(Int32), "index")
  LLVM.AddIncoming(index, [first_index], [entry_block])
  c_if = LLVM.BuildICmp(builder, LLVM.IntULE, index, limit, "continue_if")
  LLVM.BuildCondBr(builder, c_if, loop_block, exit_block)

  # loop:
  LLVM.PositionBuilderAtEnd(builder, loop_block)
  dict = Dict{ANode,ValueRef}()
  for (a,ptr) in inptrs
    addr = LLVM.BuildInBoundsGEP(builder, ptr, [index], string(a.name, "_in[index]"))
    dict[a] = LLVM.BuildLoad(builder, addr, a.name)
  end

  for a in sorted
    dict[a] = Seep.llvm_val(a, mod, builder, eltype, map(x->dict[x], a.input)...)
    if haskey(outptrs, a)
      addr = LLVM.BuildInBoundsGEP(builder, outptrs[a], [index], string(a.name, "_out[index]"))
      LLVM.BuildStore(builder, dict[a], addr)
    end
  end
  next_index = LLVM.BuildAdd(builder, index, stride, "next_index")
  LLVM.AddIncoming(index, [next_index], [loop_block])
  LLVM.BuildBr(builder, test_block)

  # exit:
  LLVM.PositionBuilderAtEnd(builder, exit_block)
  LLVM.BuildRetVoid(builder)

  LLVM.dispose(builder)
end

const module_counter = [0]
function create_module()
  module_counter[1] += 1
  module_name = "SeepLlvmModule$(module_counter[1])"
  mod = LLVM.ModuleRef(module_name)

  LLVM.SetTarget(mod, "nvptx$(8*sizeof(Int))-nvidia-cuda")
  LLVM.SetDataLayout(mod, join(["e", "p:64:64:64",
    "i1:8:8", "i8:8:8", "i16:16:16", "i32:32:32", "i64:64:64",
    "f32:32:32", "f64:64:64",
    "v16:16:16", "v32:32:32", "v64:64:64", "v128:128:128",
    "n16:32:64"], "-"))

  return mod
end

function get_bitcode(mod::ModuleRef)
  if false
    return mktempdir() do tempdir
      t = joinpath(tempdir, "bitcode.bc")
      LLVM.WriteBitcodeToFile(mod, t)
      bitcode = open(readbytes, t)
      rm(t)
      bitcode
    end
  else
    t = "bitcode.bc"
    LLVM.WriteBitcodeToFile(mod, t)
    run(`llvm-dis $t`)
    bitcode = open(readbytes, "bitcode.ll")
    return bitcode
  end
end

function get_ptx(bitcode)
  p = NVVM.Program()
  NVVM.errcheck(NVVM.AddModuleToProgram(p, open(readbytes, "/opt/cuda/nvvm/libdevice/libdevice.compute_20.10.bc"), "libdevice"))
  NVVM.errcheck(NVVM.AddModuleToProgram(p, bitcode, "seep"))
  ret = NVVM.CompileProgram(p, ASCIIString[])
  log = NVVM.GetProgramLog(p)
  if log != "\0"
    info(log)
  end

  if ret != 0
    error(bytestring(NVVM.GetErrorString(ret)))
  end

  open(io->write(io, bytestring(NVVM.GetCompiledResult(p))), "bitcode.ptx", "w")

  ptx = bytestring(NVVM.GetCompiledResult(p))
  NVVM.DestroyProgram([p])
  return ptx
end

function check_ready(readyvec::Set{ANode}, readynovec::Set{ANode}, live::Set{ANode}, inputs, nodes)
  for node in nodes
    if node in live continue end
    if all(x->x in live, inputs[node])
      push!(elementwise(node) ? readyvec : readynovec, node)
    end
  end
end

function check_dead(live::Set{ANode}, outputs, nodes)
  for node in nodes
    if any(x->!(x in live), outputs[node]) continue end
    delete!(live, node)
  end
end

function llvm_instance(pool::StoragePool, nodes::Vector{ANode})
  T = atype(pool)
  stored_nodes = Set{ANode}(nodes)

  unsorted = connected_nodes(nodes)
  live = Set{ANode}()
  storage = Dict{ANode,Any}()
  userstorage = Dict{ANode,Any}()
  refcnt = Dict{ANode,RefCounter}()

  outputs,inputs = output_map(unsorted)

  for n in nodes
    if sym(n) == :noop continue end
    incr(refcnt, output_node(n))
  end

  constants = ANode[]
  c_block = :(begin end)

  for n in unsorted
    for b in n.input
      incr(refcnt, output_node(b))
    end

    if !elementwise(n)
      push!(stored_nodes, n)
      union!(stored_nodes, n.input)
    end

    if sym(n) == :input
      push!(live, n)
      userstorage[n] = storage[n] = allocate(pool, n.size...)
    end

    if sym(n) == :load && isa(arg(n), T)
      storage[n] = arg(n)
      push!(live, n)
    end

    if sym(n) == :const
      push!(live, n)
      if isa(arg(n), T)
        storage[n] = arg(n)
      else
        incr(refcnt, n)
        t = storage[n] = allocate(pool, n.size...)
        push!(c_block.args, :(copy!($t, map($(eltype(T)), $(arg(n))))))
        push!(constants, n)
      end
    end
  end

  @gensym initialize_constants
  @eval $(initialize_constants)() = $c_block

  readyvec = Set{ANode}()
  readynovec = Set{ANode}()
  check_ready(readyvec, readynovec, live, inputs, unsorted)

  block = :(begin end)
  if !isempty(constants)
    block = :(begin
      if $(pool).token[1] != $(QuoteNode(initialize_constants))
        $(pool).token[1] = $(QuoteNode(initialize_constants))
        $(initialize_constants)()
      end
    end)
  end

  mod = create_module()
  function_counter = 0

  node_counter = length(live)
  push!(block.args, :(device_synchronize()))
  need_sync = false
  free_nodes = ANode[]

  streams = [CuStream(1) for i in 1:50]
  current_stream = CUDAdrv.default_stream()
  while true
    #println("live: " * join(sort(collect(Int, map(x->x.id, live))), ", "))
    #println("readyvec: " * join(sort(collect(Int, map(x->x.id, readyvec))), ", "))
    #println("readynovec: " * join(sort(collect(Int, map(x->x.id, readynovec))), ", "))

    if !isempty(readynovec)
      ready = collect(readynovec)
      empty!(readynovec)
      i = 1
      #@show node.id

      if need_sync
        push!(block.args, :(device_synchronize()))
        need_sync = false
        for aa in free_nodes
          deallocate(pool, storage[aa])
        end
        empty!(free_nodes)
      end

      sort!(ready, by=length, rev=true)
      for node in ready
        if sym(node) == :noop continue end
        input = map(x->storage[x], node.input)

        if current_stream != streams[(i-1)%50+1] || true
          current_stream = streams[(i-1)%50+1]
          push!(block.args, :(STREAM[1] = $current_stream))
          push!(block.args, :(CUBLAS.cublasSetStream_v2(CUBLAS.cublashandle[1], STREAM[1].handle)))
        end

        if sym(node) == :input || sym(node) == :const
          error("oops");
        elseif sym(node) == :load
          if !haskey(storage, node)
            storage[node] = allocate(pool, node.size...)
          end
          push!(block.args, :(copy!($(storage[node]), arg($node), stream=$current_stream)))
        elseif mutates(node)
          storage[node] = input[1]
          push!(block.args, :(do_forward!($node, $(input...))))
        else
          if !haskey(storage, node)
            storage[node] = allocate(pool, node.size...)
          end
          push!(block.args, :(do_forward!($node, $(storage[node]), $(input...))))
        end

        for a in node.input
          aa = output_node(a)
          if decr(refcnt, aa) && arg(aa) != storage[aa]
            push!(free_nodes, aa)
          end
        end
        i += 1
      end

      if length(ready) > 1
        need_sync = true
      end

      for node in ready
        push!(live, node)
        check_ready(readyvec, readynovec, live, inputs, outputs[node])
        check_dead(live, outputs, inputs[node])
        node_counter += 1
      end
    elseif !isempty(readyvec)
      vnodes = ANode[]
      out_nodes = Set{ANode}()

      maxlen = maximum(length, readyvec)
      while true
        a = collect(filter(x->length(x) == maxlen, readyvec))
        if isempty(a) break end

        setdiff!(readyvec, a)
        append!(vnodes, a)
        union!(live, a)
        for node in a
          check_ready(readyvec, readynovec, live, inputs, outputs[node])
          if node in stored_nodes
            push!(out_nodes, node)
            if !haskey(storage, node)
              storage[node] = allocate(pool, node.size...)
            end
          end
          check_dead(live, outputs, node.input)

          for i in node.input
            aa = output_node(i)
            if decr(refcnt, aa) && aa in stored_nodes && arg(aa) != storage[aa]
              push!(free_nodes, aa)
            end
          end

          node_counter += 1
        end
      end

      for node in vnodes
        if !iszero(refcnt[node])
          push!(out_nodes, node)
          push!(stored_nodes, node)
          if !haskey(storage, node)
            storage[node] = allocate(pool, node.size...)
          end
        end
        @assert length(node) == length(vnodes[1])
      end

      ccs = connected_components(vnodes, inputs, outputs)
      opcount(cc::Vector{ANode}) = sum(length, cc)
      order = Base.Order.By(opcount)
      Collections.heapify!(ccs, order)

      while length(ccs) > 1
        a = Collections.heappop!(ccs, order)
        b = Collections.heappop!(ccs, order)
        if opcount(a) + opcount(b) < 1024*1024
          append!(a, b)
          Collections.heappush!(ccs, a, order)
        else
          push!(ccs, a)
          push!(ccs, b)
          break
        end
      end

      sort!(ccs, by=opcount, rev=true)

      if need_sync
        push!(block.args, :(device_synchronize()))
        need_sync = false
        for aa in free_nodes
          deallocate(pool, storage[aa])
        end
        empty!(free_nodes)
      end

      i = 1
      for cc in ccs
        function_counter += 1
        function_name = "function_$function_counter"
        in_nodes = Set{ANode}()
        for c in cc
          union!(in_nodes, c.input)
        end
        setdiff!(in_nodes, cc)

        in_node_vec = collect(ANode, in_nodes)
        out_node_vec = collect(ANode, filter(c->c in out_nodes, cc))

        args = [length(cc[1]); map(x->storage[output_node(x)], [out_node_vec; in_node_vec])]
        arg_types = (Int32, fill(Ptr{eltype(T)}, length(args)-1)...)
        build_grid_stride_loop(mod, function_name, eltype(T), out_node_vec, in_node_vec)
        push!(block.args, :(cudacall(g._fun[$function_counter], 256, 128, $arg_types, $(args...); stream=$(streams[(i-1)%50+1]))))
        i += 1
      end

      if length(ccs) > 1
        need_sync = true
      end
    else
      break
    end
  end

  current_stream = CUDAdrv.default_stream()
  push!(block.args, :(STREAM[1] = $current_stream))
  push!(block.args, :(CUBLAS.cublasSetStream_v2(CUBLAS.cublashandle[1], STREAM[1].handle)))
  push!(block.args, :(device_synchronize()))

  for aa in free_nodes
    deallocate(pool, storage[aa])
  end

  @assert node_counter == length(unsorted)

  for n in nodes
    if sym(n) == :noop continue end
    aa = output_node(n)
    userstorage[n] = storage[aa]
    if decr(refcnt, aa) && arg(aa) != storage[aa]
      deallocate(pool, storage[aa])
    end
  end

  for i in constants
    if decr(refcnt, i)
      deallocate(pool, storage[i])
    end
  end

  @assert all(iszero, values(refcnt))

  bitcode = get_bitcode(mod)
  LLVM.dispose(mod)

  info("Compiling to PTX")
  ptx = @time get_ptx(bitcode)

  m = CuModuleData(ptx)

  @gensym SeepLlvmInstance
  type_block = quote
    immutable $SeepLlvmInstance <: Instance
      _mod::CuModule
      _fun::Vector{CuFunction}
      _stream::Vector{CuStream}
      _pool::StoragePool
      _storage::Dict{ANode,$T}
    end
    Base.show(io::IO, ::$SeepLlvmInstance) = print(io, "Seep LLVM Instance with $($(length(unsorted))) nodes.")
    @compat (g::$SeepLlvmInstance)() = $block
    $SeepLlvmInstance($m, $(CuFunction[CuFunction(m, "function_$i") for i in 1:function_counter]), $streams, $pool, $userstorage,)
  end

  for (n,a) in userstorage
    if haskey(userstorage, n) && n.name != ""
      if length(filter(x->x.name == n.name, unsorted)) > 1
        warn("There are multiple nodes named $(n.name).")
        continue
      end
      push!(type_block.args[2].args[3].args, :($(symbol(n.name))::$(typeof(a))))
      push!(type_block.args[end].args, a)
    end
  end

  info("Evaluating type")
  return @time eval(type_block)
end


text:
@testset "Simplify" begin
  Qx, x = QQ["x"]
  K, a = number_field(x^2 - x + 1289)
  Kt, t = K["t"]
  L, = number_field(t^2 - 5, cached = false)
  k, = @inferred Hecke.simplified_absolute_field(L, cached = false)
  @test is_isomorphic(k, number_field(x^4 + 513*x^2 + 67081)[1])

  K, a = number_field(x^3 + 42*x - 192, "a", cached = false);
  Kt, t = K["t"]
  L, = number_field(t^2 + 3, cached = false)
  k, = @inferred  Hecke.simplified_absolute_field(L)

  K, a = number_field(x^6 - 2*x^5 + x^4 - 398*x^3 - 15961*x^2 + 41496*x + 829621)
  Kx, x = K["x"]
  g = x^2 + 715563140//10487238697*a^5 - 6095690353//10487238697*a^4 +
  48931518460//10487238697*a^3 - 538842733308//10487238697*a^2 -
  6596745411420//10487238697*a + 4782497108432//806710669
  L, b = number_field(g, "b");
  k, = @inferred Hecke.simplified_absolute_field(L)

  Qx, x = QQ["x"]
  K, a = number_field(x^3 + 42*x - 192, "a")
  Kt, t = K["t"]
  L, = number_field([t^2 + 3], cached = false);
  k, = @inferred absolute_simple_field(L, cached = false, simplify = true)

  Qx, x = QQ["x"]
  f = x^2 - x + 743
  K, a = number_field(f, "a")
  Kt, t = K["t"]
  L, = number_field([t^2 - a + 27])
  k, = @inferred absolute_simple_field(L, cached = false, simplify = true)
end


text:
using Weave

cd("D:/animated-adventure-in-mathematics/")

documentsListRaw = ["ElementaryFnTrans",
    "ConicalCurvesNContours",
    "ParametricSurfaces",
    "RandomPts",
    "TrigFnNTransformations"]
documentInList = documentsListRaw.* "/" .* documentsListRaw.*".jmd"
formats = ["github",
            "md2pdf"][1]#,
            #"pandoc2html"]

for i in 1:length(documentInList)#, j in formats

    j = "github"

    err = ""

    try
    print("Generating document:"*documentInList[i]*",$j \n")

    cd("D:/animated-adventure-in-mathematics/"*documentsListRaw[i])

    weave(documentsListRaw[i]*".jmd",doctype = j)#,fig_path = documentsList[i])
    catch
        err
    end
    print(err)
    #cd("D:/visjl/")
    #weave(i,doctype = "md2pdf")
end


text:
# ==================================================================================================================
function solveDiffusion(EL2NOD, θStokes, rStokes, Δt,κ, dQdT, ρ, Cp, shear_heating,  PhaseID,
    T,∂Ωt,tfix,tfree,CMidx,to)

    @timeit to "Thermal diffusion" begin
        @timeit to "Assembly" CC, KK, Rhs =
            assembly_thermal_cylindric(EL2NOD, θStokes, rStokes,
                    Δt, κ, dQdT, ρ, Cp, shear_heating, PhaseID, CMidx)

        @timeit to "BCs" begin
            # KK, Rhs = _prepare_matrices(KK, Rhs, CC, T, Δt)
            _prepare_matrices!(KK, Rhs, CC, T, Δt)
            _apply_bcs!(T,KK,Rhs,∂Ωt,tfix)
        end

        # _CG!(T,KK,Rhs,tfree) # faster depending on MKL/BLAS setup and nthreads?
        @timeit to "Solve" T[tfree] .= KK[tfree,tfree]\Rhs[tfree]

    end
    return T,to
end
# ==================================================================================================================

@inline function _apply_bcs!(T,KK,Rhs,∂Ω,vfix)
    @views T[∂Ω] .= vfix # write prescribed temperature values
    Rhs .-=  KK[:,∂Ω] * vfix # move BC to rhs
end

@inline function _prepare_matrices(KK ,Rhs, CC, T, Δt)
    α  = 1.0
    if α != 1.0
        Rhs .+= (KK - (Δt*(1-α))*CC) * T
    else
        Rhs .+= KK*T
    end
    KK .+= @. (Δt*α)*CC
    return KK, Rhs
end

@inline function _prepare_matrices!(KK, Rhs, CC, T, Δt)
    α  = 1.0
    if α != 1.0
        Rhs .+= (KK - (Δt*(1-α))*CC) * T
    else
        Rhs .+= KK*T
    end
    KK .+= @. (Δt*α)*CC
    return KK, Rhs
end

function assembly_thermal_cylindric(EL2NOD, θStokes, rStokes,
    Δt, κ, dQdt, ρ, Cp, shear_heating, PhaseID,CMidx)

    # ============================================ MODEL AND BLOCKING PARAMETERS
    ndim    = 2
    nip     = 6
    nnodel  = size(EL2NOD,1)
    nel     = size(EL2NOD,2)
    nelblk  = min(nel, 1200)
    # nelblk  = min(nel, nblk)
    nblk    = ceil(Int,nel/nelblk)
    il      = one(nelblk); iu = nelblk

    # ======================== STORAGE FOR DATA OF ALL ELEMENT MATRICES/VECTORS
    n1 = Int(nnodel*(nnodel+1)/2)
    C_all = fill(0.0,n1,nel)
    M_all = fill(0.0,n1,nel);
    Rhs_all = fill(0.0,nnodel,nel);

    # =========== PREPARE INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    ni, nn,nn3 = Val(nip), Val(nnodel), Val(3)
    N, ∇N, dN3ds, w_ip = _get_SF(ni,nn)
    N3, = _get_SF(ni,nn3)
    detJ_PL = Vector{Float64}(undef,nelblk)
    R_31 = similar(detJ_PL) # =  detJa_PL*dxi_dth
    R_21 = similar(detJ_PL) # = -detJa_PL*deta_dth
    Th_31 = similar(detJ_PL) # = -detJa_PL*dxi_dr
    Th_21 = similar(detJ_PL) # =  detJa_PL*deta_dr
    NxN = Array{Float64,2}(undef,nnodel,nnodel)
    dNdx = Array{Float64,2}(undef,nelblk,size(∇N[1],2))
    dNdy = similar(dNdx)

    # =====================================================================
    # STORAGE FOR DATA OF ELEMENTS IN BLOCK
    # =====================================================================
    C_blk = fill(0.0,nelblk, n1)
    M_blk = fill(0.0,nelblk, n1)
    Rhs_blk = fill(0.0,nelblk, nnodel)
    ω  = Vector{Float64}(undef,Int(nelblk))
    invJx_double = Array{Float64,2}(undef,Int(nelblk),ndim) # storage for x-components of Jacobi matrix
    invJz_double = similar(invJx_double) # storage for z-components of Jacobi matrix

    # =========================================================================
    # BLOCK LOOP - MATRIX COMPUTATION
    # =========================================================================
    for ib ∈ 1:nblk

        #===========================================================================
        CALCULATE JACOBIAN, ITS DETERMINANT AND INVERSE
        ===========================================================================#
        #=
        NOTE: For triangular elements with non-curved edges the Jacobian is
              the same for all integration points (i.e. calculated once
              before the integration loop). Further, linear 3-node shape are
              sufficient to calculate the Jacobian.
        =#
        VCOORD_th = view(θStokes,:,il:iu)
        VCOORD_r = view(rStokes,:,il:iu)
        J_th = VCOORD_th'*dN3ds'
        J_r = VCOORD_r'*dN3ds'

        _fill_R_J_thermal!(J_th, J_r, VCOORD_th, VCOORD_r,
                   R_31, R_21, Th_31, Th_21, detJ_PL)

        # ------------------------ NUMERICAL INTEGRATION LOOP (GAUSS QUADRATURE)
        _ip_loop!(C_blk, M_blk, Rhs_blk,
                  Δt, κ, dQdt, ρ, Cp, shear_heating, EL2NOD, PhaseID, il:iu,
                  N3, N, ∇N, w_ip, nip, nnodel,
                  dNdx, dNdy, NxN,
                  ω, invJx_double, invJz_double, detJ_PL,
                  R_21, R_31, Th_21, Th_31,
                  VCOORD_th, VCOORD_r)

        # --------------------------------------- WRITE DATA INTO GLOBAL STORAGE
        @views begin
            C_all[:,il:iu] = C_blk'
            M_all[:,il:iu] = M_blk'
            Rhs_all[:,il:iu] = Rhs_blk'
        end

        emptyblocks!(C_blk,
                     M_blk,
                     Rhs_blk)

        # -------------------------------- READJUST START, END AND SIZE OF BLOCK
        il += nelblk;
        if ib == nblk-1
           # ------------ Account for different number of elements in last block
           nelblk = nel-iu  # new block size
           # ------------------------- Reallocate at the edge the block sequence
           detJ_PL          = Vector{Float64}(undef,nelblk)
           R_31             = similar(detJ_PL) # =  detJa_PL*dxi_dth
           R_21             = similar(detJ_PL) # = -detJa_PL*deta_dth
           Th_31            = similar(detJ_PL) # = -detJa_PL*dxi_dr
           Th_21            = similar(detJ_PL) # =  detJa_PL*deta_dr
           dNdx             = Array{Float64,2}(undef,nelblk,size(∇N[1],2))
           dNdy             = similar(dNdx)
           C_blk            = fill(0.0, nelblk, n1)
           M_blk            = fill(0.0, nelblk, n1)
           Rhs_blk          = fill(0.0, nelblk, nnodel)
           ω                = Vector{Float64}(undef, Int(nelblk))
           invJx_double     = Array{Float64,2}(undef, Int(nelblk),ndim) # storage for x-components of Jacobi matrix
           invJz_double     = similar(invJx_double) # storage for z-components of Jacobi matrix

        end
        iu += nelblk

    end # end block loop

    #===========================================================================
    ASSEMBLY OF GLOBAL SPARSE MATRICES AND RHS-VECTOR
    ===========================================================================#
    C = vec(C_all)
    M = vec(M_all)
    R = vec(Rhs_all)
    E = vec(view(EL2NOD,1:nnodel,:))
    CC, MM, Rhs = _thermal_matrices(C,M,R,E,CMidx)

    return CC,MM,Rhs

end # END OF ASSEMBLY FUNCTION

#===============================================================================
 ------------------------------------------------------------------------------
                        BEGINING OF ASSEMBLY SUBFUNCTIONS
 ------------------------------------------------------------------------------
===============================================================================#
@inline function _thermal_matrices(C,M,R,E,CM_i,CM_j)
    tmp   = sparse(CM_i, CM_j, C)
    CC    = tmp + tril(tmp,-1)'; # other half is given by symmetry
    tmp   = sparse(CM_i, CM_j, M);
    MM    = tmp + tril(tmp,-1)'; # other half is given by symmetry
    Rhs   = accumarray(E, R)

    return CC, MM, Rhs
end

@inline function _thermal_matrices(C,M,R,E,CMidx)
    tmp   = sparse(CMidx._i, CMidx._j, C)
    CC    = tmp + tril(tmp,-1)'; # other half is given by symmetry
    tmp   = sparse(CMidx._i, CMidx._j, M);
    MM    = tmp + tril(tmp,-1)'; # other half is given by symmetry
    Rhs   = accumarray(E, R)

    return CC, MM, Rhs
end

function _create_triplets(EL2NOD,nnodel)

    indx_j  = Array{Int64}(undef,nnodel,nnodel)
    indx_i  = Array{Int64}(undef,nnodel,nnodel)
    dummy   = collect(1:nnodel)
    @inbounds for i ∈ 1:3
        indx_j[i,:] = dummy
        indx_i[:,i] = dummy'
    end

    # indx_i      = indx_j';
    indx_i      = tril(indx_i); indxx_i = vec(indx_i); filter!(x->x>0,indxx_i)
    indx_j      = tril(indx_j); indxx_j = vec(indx_j); filter!(x->x>0,indxx_j)

    CM_i         = copy(EL2NOD[indxx_i,:]); CMM_i = vec(CM_i)
    CM_j         = copy(EL2NOD[indxx_j,:]); CMM_j = vec(CM_j)
    indx         = CMM_i .< CMM_j
    tmp          = copy(CMM_j[indx])
    CMM_j[indx]  = CM_i[indx]
    CMM_i[indx]  = tmp

    return CMM_i,CMM_j

end

function _ip_loop!(C_blk, M_blk, Rhs_blk,
    Δt,κ, dQdt, ρ, Cp, shear_heating, EL2NOD, PhaseID,els,
    N3,N, ∇N, w_ip, nip,  nnodel,
    dNdx,dNdy, NxN,
    ω,invJx_double, invJz_double, detJ_PL,
    R_21,R_31,Th_21,Th_31,
    VCOORD_th,VCOORD_r)

    th_ip   = similar(detJ_PL)
    r_ip    = similar(detJ_PL)

    @inbounds for ip=1:nip

        N_ip  = N[ip]

        #=======================================================================
        PROPERTIES OF ELEMENTS AT ip-TH EVALUATION POINT
        =======================================================================#
        Cond_blk = _element_conductivity(κ,PhaseID,els)
        dQdt_blk = _element_dQdt(dQdt,PhaseID,els)
        RhoCp_blk = _element_RhoCp(ρ,Cp,EL2NOD,PhaseID, els, N_ip)
        sh_blk = view(shear_heating,els,ip)

        #===========================================================================================================
        CALCULATE 2nd JACOBIAN (FROM CARTESIAN TO POLAR COORDINATES --> curved edges), ITS DETERMINANT AND INVERSE
        ===========================================================================================================
        NOTE: For triangular elements with curved edges the Jacobian needs to be computed at each integration
               point (inside the integration loop).
        TODO: fix type stability error =#
        th_ip   = VCOORD_th'*N3[ip]'
        r_ip    = VCOORD_r'*N3[ip]'

        _derivative_weights!(∇N[ip],ω,dNdx,dNdy,w_ip[ip],th_ip,r_ip,
            R_21,R_31,Th_21,Th_31,detJ_PL,invJx_double,invJz_double)

        #=======================================================================
        3-node triangle allows lumping to improve stability
        =======================================================================#
        if nnodel == 3 # lumping
            mul!(NxN,N_ip',N_ip)
            @inbounds for ii ∈ 1:nnodel
                NxN[ii,ii] = NxN[ii,1] + NxN[ii,2]+ NxN[ii,3]
            end
            @inbounds for ii ∈ 1:nnodel,jj ∈ 1:nnodel
                if ii != jj
                    NxN[ii,jj] = 0.0
                end
            end
        else
            mul!(NxN,N_ip',N_ip)
        end

       _fill_blk_matrices!(C_blk,M_blk,Rhs_blk, sh_blk,
                ω,dNdx,dNdy,NxN,N_ip, nnodel,
                Cond_blk,dQdt_blk,RhoCp_blk,Δt)

    end # end integration point loop
end # ENF OF IP_LOOP FUNCTION

function sumrows(a::Array{Float64,2})
    nrow,ncol = size(a,1), size(a,2)
    b         = fill(0.0,nrow)
    for i ∈ 1:nrow,j ∈ 1:ncol
            b[i] += a[i,j]
    end
    return b
end

function _fill_R_J_thermal!(J_th,J_r,VCOORD_th,VCOORD_r,
    R_31,R_21,Th_31,Th_21,detJ_PL)

    VCOORD_r1    = view(VCOORD_r,1,:)
    VCOORD_r2    = view(VCOORD_r,2,:)
    VCOORD_r3    = view(VCOORD_r,3,:)
    VCOORD_th1   = view(VCOORD_th,1,:)
    VCOORD_th2   = view(VCOORD_th,2,:)
    VCOORD_th3   = view(VCOORD_th,3,:)

    @inbounds for i ∈ 1:length(R_31)
        detJ_PL[i] = J_th[i,1]*J_r[i,2] - J_th[i,2]*J_r[i,1]
        R_31[i]    = VCOORD_r3[i]  - VCOORD_r1[i]  # =  detJa_PL*dxi_dth
        R_21[i]    = VCOORD_r2[i]  - VCOORD_r1[i]  # = -detJa_PL*deta_dth
        Th_31[i]   = VCOORD_th3[i] - VCOORD_th1[i] # = -detJa_PL*dxi_dr
        Th_21[i]   = VCOORD_th2[i] - VCOORD_th1[i] # =  detJa_PL*deta_dr
    end

    return R_31,R_21,Th_31,Th_21,detJ_PL
end

function _derivative_weights!(∇N,ω,dNdx,dNdy,w_ip,th_ip,r_ip,
    R_21,R_31,Th_21,Th_31,detJ_PL,invJx_double,invJz_double)

    # --- Start loop (faster and more memory friendly than single - line declarations)
    for i ∈ 1:length(th_ip)
        # -------------------------------- i-th extractions for some extra speed
        detJ    = detJ_PL[i]
        invdetJ = 1/detJ
        sin_ip  = sin(th_ip[i])
        cos_ip  = cos(th_ip[i])
        R_31_i  = R_31[i]
        R_21_i  = R_21[i]
        Th_31_i = Th_31[i]
        Th_21_i = Th_21[i]
        r_i     = r_ip[i]
        invr_i = 1/r_i
        # --------------------------------------- Inverse of the double Jacobian
        invJx_double[i,1] = ( R_31_i*cos_ip*invr_i - Th_31_i*sin_ip ) * invdetJ
        invJx_double[i,2] = (-R_21_i*cos_ip*invr_i + Th_21_i*sin_ip ) * invdetJ
        invJz_double[i,1] = (-R_31_i*sin_ip*invr_i - Th_31_i*cos_ip ) * invdetJ
        invJz_double[i,2] = ( R_21_i*sin_ip*invr_i + Th_21_i*cos_ip ) * invdetJ
        # ---------------------------- Numerical integration of element matrices
        ω[i]  = r_i * detJ * w_ip
        # ω[i]  = detJ * w_ip
    end

    # --------------------------------------- Derivatives wrt global coordinates
    mul!(dNdx,invJx_double, ∇N)
    mul!(dNdy,invJz_double, ∇N)

end

@inline function _fill_blk_matrices!(C_blk,M_blk,Rhs_blk, sh_blk,
    ω, dNdx, dNdy, NxN, N_ip, nnodel,
    Cond_blk::Float64, dQdt_blk::Float64, RhoCp_blk, Δt)

    indx = 1
    nn = size(M_blk,1)
    @inbounds for i ∈ 1:nnodel
        for j ∈ i:nnodel
            # Conductivity matrices of all elements in block
            # C = B'*kappa*B, with B = [ dN1dx dN2dx dN3dx dN4dx dN5dx dN6dx;
            #                            dN1dy dN2dy dN3dy dN4dy dN5dy dN6dy ]
            NxN_ij = NxN[i,j]
            @fastmath  @simd for k ∈ 1:nn
                C_blk[k,indx] +=  ω[k] *
                    Cond_blk*(dNdx[k,i]*dNdx[k,j] + dNdy[k,i]*dNdy[k,j])

                M_blk[k,indx] +=  ω[k]* RhoCp_blk[k] * NxN_ij #N_ip[i]*N_ip[j] # NxN_ij
            end
            indx += 1;
        end
    end

    # RIGHT HAND SIDE
    if dQdt_blk !=0.0
        Rhs_blk .+= (N_ip'*(Δt* dQdt_blk*ω)')';
    end

    # Rhs_blk .+= (Δt* sh_blk.*ω)*N_ip  # shear heating

end

@inline function _fill_blk_matrices!(C_blk,M_blk,Rhs_blk,
    ω,dNdx,dNdy,NxN,N_ip, nnodel,
    Cond_blk::Array{Float64,1},dQdt_blk::Array{Float64,1},RhoCp_blk,Δt)

    indx   = 1
    nn     = size(M_blk,1)
    @inbounds for i ∈ 1:nnodel
        for j ∈ i:nnodel
            # Conductivity matrices of all elements in block
            # C = B'*kappa*B, with B = [ dN1dx dN2dx dN3dx dN4dx dN5dx dN6dx;
            #                            dN1dy dN2dy dN3dy dN4dy dN5dy dN6dy ]
            NxN_ij = NxN[i,j]
            @fastmath  @simd for k ∈ 1:nn
                C_blk[k,indx] +=  ω[k] *
                    Cond_blk* (dNdx[k,i]*dNdx[k,j] + dNdy[k,i]*dNdy[k,j])

                M_blk[k,indx] +=  ω[k]* RhoCp_blk[k] * NxN_ij
            end
            indx += 1;
        end
    end

    # RIGHT HAND SIDE
    if sum(dQdt_blk) !=0.0
        Rhs_blk .+= @. (N_ip'*(Δt* dQdt_blk*ω)')';
    end

end

# @inline function  _element_conductivity(κ,PhaseID,els)
#     # Use the constant conductivity defined in PHYSICS.K
#     length(PhaseID)==1 ? Cond_blk = κ[PhaseID]' : Cond_blk = vec(κ[PhaseID[els]]')

#     return  Cond_blk
# end

_element_conductivity(κ,PhaseID,els) = length(PhaseID)== 1 ? κ[PhaseID]' : vec(κ[PhaseID[els]]')

@inline function _element_RhoCp(ρ,Cp,EL2NOD,PhaseID,els,Nip)

    Dens_blk = vec(ρ[view(EL2NOD,:,els)]'*Nip')

    # Use the constant conductivity defined in PHYSICS.K
    length(PhaseID)==1 ? Cp_blk = Cp[PhaseID]' : Cp_blk = vec(Cp[PhaseID[els]]')
    # Cp_blk = ifelse(length(PhaseID)==1,
    #                 Cp[PhaseID]',
    #                 vec(Cp[PhaseID[els]]'))
    # if length(PhaseID)==1
    #     # Dens_blk    = ρ[PhaseID]'
    #     Cp_blk      = Cp[PhaseID]'
    # else
    #     # Dens_blk    = vec(ρ[PhaseID[els]]')
    #     Cp_blk      = vec(Cp[PhaseID[els]]')
    # end

    # Make sure a column-vector is returned
    return  Dens_blk.*Cp_blk
end

@inline function _element_dQdt(dQdt,PhaseID,els)
    # Use the constant conductivity defined in PHYSICS.K
    length(PhaseID)==1 ? dQdt_blk = dQdt[PhaseID]' : dQdt_blk = dQdt[PhaseID[els]]'
    return dQdt_blk
end

@inline function emptyblocks!(C_blk,M_blk,Rhs_blk)
    @inbounds for i in eachindex(C_blk)
        C_blk[i] = 0.0
        M_blk[i] = 0.0
    end
    @inbounds for i in eachindex(Rhs_blk)
        Rhs_blk[i] = 0.0
    end
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
export noop

import Base: show

"""
    noop()

Creates a noop operator, which does nothing, but breaks operator composition type inference checking procedure for Julia's compiler.
It might be useful for very long chain of operators, because Julia tries to statically infer data types at compile-time for the whole chain and
can run into StackOverflow issues.

```jldoctest
using Rocket

# foo() block to enforce local scope for scope variable
function foo()
    source = from(1:5)

    for i in 1:1000
        source = source |> map(Int, d -> d + 1) |> noop()
    end

    subscribe!(source, logger())
end

foo()
;

# output

[LogActor] Data: 1001
[LogActor] Data: 1002
[LogActor] Data: 1003
[LogActor] Data: 1004
[LogActor] Data: 1005
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`logger`](@ref), [`map`](@ref)
"""
noop() = NoopOperator()

struct NoopOperator <: InferableOperator end

operator_right(operator::NoopOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::NoopOperator, source) where L
    return proxy(L, source, NoopProxy())
end

struct NoopProxy <: ActorSourceProxy end

actor_proxy!(::Type{L}, proxy::NoopProxy, actor)   where L = NoopActor{L}(actor)
source_proxy!(::Type{L}, proxy::NoopProxy, source) where L = NoopSource{L}(source)

struct NoopActor{L} <: Actor{L}
    actor
end

on_next!(actor::NoopActor{L}, data::L) where L = next!(actor.actor, data)
on_error!(actor::NoopActor, err)               = error!(actor.actor, err)
on_complete!(actor::NoopActor)                 = complete!(actor.actor)

struct NoopSource{L} <: Subscribable{L}
    source
end

function on_subscribe!(source::NoopSource, actor)
    return subscribe!(source.source, actor)
end

Base.show(io::IO, ::NoopOperator)          = print(io, "NoopOperator()")
Base.show(io::IO, ::NoopProxy)             = print(io, "NoopProxy()")
Base.show(io::IO, ::NoopActor{L})  where L = print(io, "NoopActor($L)")
Base.show(io::IO, ::NoopSource{L}) where L = print(io, "NoopSource($L)")