text:
using JuMP;
using Ipopt;
using DataFrames;
my_table = readtable( "/home/iriadmin/newdata.csv");
( N, K1 ) = size( my_table );
m = Model(solver=IpoptSolver());
@defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( my_table[i,1] - sum{ my_table[i,j] * b[ j - 1] , j in 2:K1} )^2 , i in 1:N} )
statusM = solve( m )

beta = getvalue(b)

# ---------------------------- R ---------------------------------------------------------
library( Matrix )
myd <- read.csv( "C:\\CBIG\\ZD\\logistic_regression_ex2.csv" )
x <- as.matrix( myd[, -c( 1 )] )
y <- myd[, c(1)]
## make sure its rank = # columns
ncol(x)
rankMatrix(x)
myglm <- glm( y ~ -1 + x, family = binomial( link = logit ) )
myglm_coeff <- summary( myglm )$coefficients
rownames( myglm_coeff ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )
print( myglm_coeff )
## using optim to reproduce the estimates
Loglikelihood <- function( beta )
{ xbeta <- x %*% beta
  Li <- -1 * y * xbeta + log( 1 + exp( xbeta ))
  sum( Li )
}
beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )
myoptim <- optim( beta_start_value, fn = Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE )
var_cov_m <- solve( myoptim$hessian )
beta_se <- sqrt( diag( var_cov_m ) )
zvalue <- myoptim$par / beta_se
pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )
coeff_outputs <- cbind( myoptim$par, beta_se, zvalue, round( pvalue, 6 ) )
colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )
print( coeff_outputs )
print( summary( myglm )$coefficients )

#######################################################################

## this is the residual deviance
2 * myoptim$value

########## since colmn 4 to column 7 are indicators, assuming they are the compaign breaks :)
##### so we want to add penalty function to make the estimates >= 0, using the penalty that Zhiyuan recommended

## assuming the penalty function is :

a <- c( 1000, 9, 60 )

Penalized_Loglikelihood <- function( beta )
{
	xbeta <- x %*% beta
	Li <- -1 * y * xbeta + log( 1 + exp( xbeta ))
	## constraining beta[3] to beta[7]
	myPenValue <- 0
	for ( i in 4 : 7 )
	{
		myPenValue <- myPenValue + a[1] * exp( -1.0 * a[2] * exp( a[3] * beta[ i ] ) )
	}
	sum( Li ) + myPenValue
}

beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )

myoptim_Penalized <- optim( beta_start_value, fn = Penalized_Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE)

myoptim_Penalized

var_cov_m <- solve( myoptim_Penalized$hessian )

beta_se <- sqrt( diag( var_cov_m ) )

zvalue <- myoptim_Penalized$par / beta_se

pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )

coeff_outputs <- cbind( myoptim_Penalized$par, beta_se, zvalue, round( pvalue, 6 ) )

colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )

print( coeff_outputs )

# ------------------------------ GMT ------------------------------------
using JuMP, DataFrames
# Use nonlinear optimization to compute the maximum likelihood estimate (MLE) of the parameters of a normal distribution  aka the sample mean and variance
n = 1000
data = randn(n)
m = Model()
@variable(m, μ, start = 0.0)
@variable(m, σ >= 0.0, start = 1.0)
@NLobjective(m, Max, (n/2)*log(1/(2π*σ^2))-sum((data[i]-μ)^2 for i=1:n)/(2σ^2))
solve(m)
println("μ = ", getvalue(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", getvalue(σ)^2)
println("var(data) = ", var(data))
println("MLE objective: ", getobjectivevalue(m))
# constrained MLE?
@NLconstraint(m, μ == σ^2)
solve(m)
println("\nWith constraint μ == σ^2:")
println("μ = ", getvalue(μ))
println("σ^2 = ", getvalue(σ)^2)
println("Constrained MLE objective: ", getobjectivevalue(m))

# ------------------------------ GMT ------------------------------------
# ------------------------------ GMT ------------------------------------
# ------------------------------ GMT ------------------------------------
using JuMP, DataFrames
# Use nonlinear optimization to compute the maximum likelihood estimate (MLE) of the parameters of a normal distribution  aka the sample mean and variance
n = 1000
data = randn(n)
m = Model()
@variable(m, μ, start = 0.0)
@variable(m, σ >= 0.0, start = 1.0)
@NLobjective(m, Max, (n/2)*log(1/(2π*σ^2))-sum((data[i]-μ)^2 for i=1:n)/(2σ^2))
solve(m)
println("μ = ", getvalue(μ),"\nmean(data) = ",mean(data), "\nσ^2 = ",getvalue(σ)^2,"\nvar(data) = ",var(data),"\nMLE objective: ",getobjectivevalue(m) )
@NLconstraint(m, μ == σ^2) # constrained MLE?
solve(m)
println("\nWith constraint μ == σ^2:","μ = ", getvalue(μ), "σ^2 = ", getvalue(σ)^2, "Constrained MLE objective: ", getobjectivevalue(m) )

# ------------ GMT HERE HERE HERE HERE HERE ---------------
using JuMP, DataFrames, GLM, NLsolve
dfd = readtable( "/home/iriadmin/g/JStack/src/tst/logreg.csv")
dfd=dfd[setdiff(names(dfd),[:x0])]
mm = ModelMatrix(ModelFrame(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd))

g = glm(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd,  Bernoulli())   #m = glm(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd, Poisson())

a = convert(Array,dfd) #array(dfd[:, 1:end])
B = zeros(length(mm.m[1,2:end]))

x = a[:,2:end]
y = a[:,1]

# sum( -1 * y .* (x*beta) + log(1+exp((x*beta))) )

function LogLHd(beta::Array{Float64})
    xbeta = x*beta
    Li = -1 * y .* xbeta + log(1+exp(xbeta))
                    println(sum(Li) )
    sum(Li)
end

#ll(B) = B*mm.m[:,2:end]
#m = Model()
#@variable(m, μ, start = 0.0)

using Optim
zz=optimize(LogLHd, B, LBFGS())

using Optim
LogL(beta) =  sum( -1 * y .* (x*beta) + log(1+exp((x*beta))) )
result = optimize(LogL, zeros(2), BFGS())

beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )
myoptim <- optim( beta_start_value, fn = Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE )
var_cov_m <- solve( myoptim$hessian )

## make sure its rank = # columns
#ncol(x)
#rankMatrix(x)
#myglm <- glm( y ~ -1 + x, family = binomial( link = logit ) )
m = Model()
@variable(m, μ, start = 0.0)
@variable(m, σ >= 0.0, start = 1.0)
                @NLobjective(m, Max, (n/2)*log(1/(2π*σ^2))-sum((dfd[i]-μ)^2 for i=1:n)/(2σ^2))
solve(m)

# ----
using NLsolve
function f!(x, fvec)
       Lb_pre = (   (mean_score1*(Mt/M))    -   (mean_score1*exp(-(B-(SE*x[1])))*(Mt/M))   )  +
                    (   (mean_score0*exp((B-(SE*x[1])))*(Mc/M))    -   (mean_score0*(Mc/M))    )
                    Lb = Lb_pre/mean_score0
                   fvec[1] = Lb
end
r=nlsolve(f!,[0.1])    #println(r)
# ---------------------------- R ---------------------------------------------------------
library( Matrix )
myd <- read.csv( "C:\\CBIG\\ZD\\logistic_regression_ex2.csv" )
x <- as.matrix( myd[, -c( 1 )] )
y <- myd[, c(1)]
## make sure its rank = # columns
ncol(x)
rankMatrix(x)
myglm <- glm( y ~ -1 + x, family = binomial( link = logit ) )
myglm_coeff <- summary( myglm )$coefficients
rownames( myglm_coeff ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )
print( myglm_coeff )
## using optim to reproduce the estimates
Loglikelihood <- function( beta )
{ xbeta <- x %*% beta
  Li <- -1 * y * xbeta + log( 1 + exp( xbeta ))
  sum( Li )
}
beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )
myoptim <- optim( beta_start_value, fn = Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE )
var_cov_m <- solve( myoptim$hessian )
beta_se <- sqrt( diag( var_cov_m ) )
zvalue <- myoptim$par / beta_se
pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )
coeff_outputs <- cbind( myoptim$par, beta_se, zvalue, round( pvalue, 6 ) )
colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )
print( coeff_outputs )
print( summary( myglm )$coefficients )

# ------------------------------ GMT END ------------------------------------
# ------------------------------ GMT END ------------------------------------
# ------------------------------ GMT END ------------------------------------

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using StatisticalRethinking
using Test

@testset "StatisticalRethinking tests" begin

  println("All tests are executed in docs job.")

end

ext: jl


text:
function myfuncA()
    println("Say hello from function A")
end

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
using CSV, DataFrames, Statistics, LinearAlgebra, Clustering

# Helper functions
dist(x, y) = norm(collect(skipmissing(x .- y)), 2)

function covar(x, y)
    id = findall(z -> !ismissing(z), x .- y)
    return cov(x[id], y[id])
end

function ret2dist(r)
    n = size(r, 2)
    Σ = fill(0.0, n, n)
    D = fill(0.0, n, n)
    for i = 1:n
        for j = i:n
            D[i, j] = dist(r[:, i], r[:, j])
            D[j, i] = D[i, j]
            Σ[i, j] = covar(r[:, i], r[:, j])
            Σ[j, i] = Σ[i, j]
        end
    end
    return Σ, D
end

function ivp(C)
    iv = 1.0 ./ diag(C)
    return iv ./ sum(iv)
end

function clustervar(C, cluster)
    V = C[cluster, cluster]
    u = ivp(V)
    return dot(u' * V, u)
end

function bisect(items)
    tmp = Vector{Int}[]
    for item in items
        n1 = length(item)
        if n1 > 1
            n2 = floor(Int, n1 / 2)
            push!(tmp, item[1:n2])
            push!(tmp, item[n2+1:end])
        end
    end
    return tmp
end

function update!(w, items, Σ)
    for k in 1:2:length(items)
        c1, c2 = items[k], items[k+1]
        v1, v2 = clustervar(Σ, c1), clustervar(Σ, c2)
        α = v1 / (v1 + v2)
        w[c1] .= w[c1] .* α
        w[c2] .= w[c2] .* (1.0 - α)
    end
end

function hrp(x; input_type = "ret")
    Σ = fill(0.0, size(x, 2), size(x, 2))
    D = similar(Σ)
    if input_type == "ret"
        Σ, D = ret2dist(x)
    elseif input_type == "cov"
        Σ, D = x, cov2dist(x)
    end
    h = hclust(D).order
    n = length(h)
    w = fill(1.0, n)
    items = [h]
    while !isempty(items)
        items = bisect(items)
        update!(w, items, Σ)
    end
    return ivp(Σ), w
end

function cov2corr(C)
    σinv = diagm(1.0 ./ sqrt.(diag(C)))
    return σinv * C * σinv
end

function cov2dist(C)
    P = cov2corr(C)
    S = similar(P)
    for i = 1:size(P, 1)
        for j = i:size(P, 1)
            S[i, j] = 2.0 * (1.0 - P[i, j])
            S[j, i] = S[i, j]
        end
    end
    return S
end

function parsedata(f)
    reits = DataFrame(CSV.File(f, missingstring = ["", " ", "#N/A N/A", "NaN", "NA"], dateformat = "yyyy-mm-dd"))
    rx = reits[2:end, :]
    rx[:, 2:end] .= log.(reits[2:end, 2:end]) .- log.(reits[1:end-1, 2:end])
    assets = names(rx)[2:end]
    rx = Matrix(rx[:, 2:end])
    return rx, assets
end

# Compute stuff
rx, assets = parsedata(ARGS[1])
rp_wt, hrp_wt = hrp(rx)
wt = DataFrame(:reits => assets, :rp => rp_wt, :hrp => hrp_wt)
CSV.write("wts.csv", wt)
println(wt)

max_stars_repo_name: banachtech/hrp

max_issues_repo_name: banachtech/hrp

max_forks_repo_name: banachtech/hrp


text:
"""
    MinLevelLogger(logger, min_enabled_level)

Wraps `logger` in an filter that runs before the log message is created.
In many ways this is just a specialised [`EarlyFilteredLogger`](@ref)
that only checks the level.
This filter only allowed messages on or above the `min_enabled_level` to pass.
"""
struct MinLevelLogger{T <: AbstractLogger, L} <: AbstractLogger
    logger::T
    min_level::L
end

function handle_message(logger::MinLevelLogger, args...; kwargs...)
    if comp_handle_message_check(logger.logger, args...; kwargs...)
        return handle_message(logger.logger, args...; kwargs...)
    end
end

function shouldlog(logger::MinLevelLogger, args...)
    return comp_shouldlog(logger.logger, args...)
end

function min_enabled_level(logger::MinLevelLogger)
    comp_min_level = min_enabled_level(logger.logger)
    # `max` since if either would not take it then do not enable
    return max(logger.min_level, comp_min_level)
end

catch_exceptions(logger::MinLevelLogger) = catch_exceptions(logger.logger)


text:
export tumor_stepper, tumor_stepper!, record_growth

function tumor_stepper!(simresult, trange;
        seed = nothing, simparams...)
    isnothing(seed) || Random.seed!(seed)
    id, mut_id, t, tumor, mutations = simresult |> r -> (r[:index], r[:mutation], r[:time], r[:tumor], r[:mutations])

    time_series = Vector{DataFrame}()
    prog = Progress(length(trange), dt=0.01)
    for tp in trange
        up = birth_death_pushing!(tumor, mutations, tp; simparams..., cur_id=id, cur_mutation=mut_id, t=t, showprogress=false)
        push!(time_series, deepcopy(DataFrame(tumor)) )
        id, mut_id, t = up[:index], up[:mutation], up[:time]

        next!(prog)
    end
    return time_series
end
function tumor_stepper(trange;
        seed = abs(rand(Int)), dim, simparams...)
    simresult = birth_death_pushing(first(trange); dim = dim, simparams..., seed=seed)
    state = simresult |> r -> (r[:index], r[:mutation], r[:time])
    tumor_stepper!(simresult, trange; simparams...)
end

function record_growth(time_series; path="temp.gif", frames=1,
        shading = false, markersize=1., size = (500,500),
        show_axis=false, colorpalette = palette(:tab20),
        points_colors = snap -> colors_by_mutations(snap; colorpalette=colorpalette, show_warning=false)
        )
    dim = length(first(first(time_series).position))
    width = maximum(abs.(vcat(last(time_series).position...)))
    origin = fill(width, dim)
    limits = dim == 2 ? FRect2D(-origin, 2*origin) : FRect3D(-origin, 2*origin)

    pointsNode = Node([Point(zeros(dim)...)])
    colorsNode = Node([colorpalette[1]])
    scene = Scene(resolution=size, show_axis=show_axis)
    meshscatter!(scene, pointsNode, color = colorsNode, scale_plot=false, limits = limits, shading = shading, markersize = markersize)
    record(scene, path) do io
        @showprogress for snapshot in time_series
            pointsNode[], colorsNode[] = points_colors(snapshot)
            for _=1:frames recordframe!(io) end
        end
    end
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
# based on https://discourse.julialang.org/t/arithmetic-modulo-primes/23895/3

struct ℤₚ{p} <: Number
   val::Int64

    function ℤₚ{n}(a) where {n}
        u = mod(a,n)
        new(u > n ÷ 2 ? u - n : u)
    end
end

val(x::ℤₚ{n}) where n = x.val

ℤₚ{n}(x::ℤₚ{n}) where n = x

Base.promote(x::ℤₚ{n}, y::Integer) where {n}=(x,ℤₚ{n}(y))
Base.promote(y::Integer, x::ℤₚ{n}) where {n}=(ℤₚ{n}(y),x)

Base.zero(::Type{ℤₚ{n}}) where {n} = ℤₚ{n}(0)
Base.one(::Type{ℤₚ{n}}) where {n} = ℤₚ{n}(1)

Base.:(==)(x::ℤₚ,y::ℤₚ) = (x.val==y.val)
Base.:(==)(x::ℤₚ{n}, k::Integer) where n = (val(x) == k)
Base.:(==)(k::Integer, x::ℤₚ) = (x == k)

Base.:+(x::ℤₚ{n}, y::ℤₚ{n}) where {n} = ℤₚ{n}(Int(x.val)+y.val)
Base.:*(x::ℤₚ{n}, y::ℤₚ{n}) where {n} = ℤₚ{n}(Int(x.val)*y.val)
Base.:-(x::ℤₚ{n}, y::ℤₚ{n}) where {n} = ℤₚ{n}(Int(x.val)-y.val)
Base.:-(x::ℤₚ{n}) where {n} = ℤₚ{n}(-Int(x.val))

Base.:/(x::ℤₚ{n}, y::ℤₚ{n}) where n = x * inv(y)
Base.:÷(x::ℤₚ{n}, y::ℤₚ{n}) where n = x * inv(y)
Base.:÷(x::ℤₚ{n}, y::Integer) where n = ℤₚ{n}(Int(x.val)*invmod(y,n))

Base.inv(x::ℤₚ{n}) where {n} = ℤₚ{n}(invmod(x.val,n))
Base.real(x::ℤₚ{n}) where {n} = x.val
Base.abs(x::ℤₚ{n}) where {n} = abs(x.val)
Base.gcd(x::ℤₚ{n}, y::ℤₚ{n}) where {n} = ℤₚ{n}(gcd(x.val, y.val))
Base.gcd(x::ℤₚ{n}, y::Integer) where {n} = gcd(x.val, ℤₚ{n}(y))
Base.gcd(x::Integer, y::ℤₚ{n}) where {n} = gcd(ℤₚ{n}(x), y)

function Base.show(io::IO, m::ℤₚ{n}) where n
    if get(io,:limit, false)
        sub = Dict(zip("0123456789,()","₀₁₂₃₄₅₆₇₈₉‚₍₎"))
        print(io, m.val, map(x->sub[x],repr(n)))
   else
        print(io,"ℤₚ{$n}($(m.val))")
   end
end

###############################################################################

modular(n, v::AbstractArray) = [ℤₚ{n}(x) for x in v]

function modular(n::Integer, p::AbstractPolynomialLike)
     polynomial(modular(n, coefficients(p)), terms(p))
end

# function modular(n::Integer, p::Polynomial{true, T}) where T<:Integer
#     polynomial(modular(n, coefficients(p)), terms(p))
# end
#
# function modular(n::Integer, p::Polynomial{true, ℤₚ{m}}) where m
#     if n == m
#         return p
#     else
#         modular(n, demodular(p))
#     end
# end
#
function demodular(p::AbstractPolynomialLike) where n
    polynomial(val.(coefficients(p)), terms(p))
end

demodular(p::Polynomial{true, T}) where T<:Integer = p


text:
using Query
using DataFrames
using NamedTuples

# We want a DataFrame without Nullable's, so
# we need to construct a DataFrame with columns of Array
columns=[]
push!(columns, ["John", "Sally", "Kirk"])
push!(columns, [23., 42., 59.])
push!(columns, [3,5,2])
df = DataFrame(columns, [:name, :age, :children])

x = @from i in df begin
    @let count = length(i.name)
    @let kids_per_year = i.children / i.age
    @where count > 4
    @select {Name=i.name, Count=count, KidsPerYear=kids_per_year}
    @collect DataFrame
end

println(x)


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

function ^(a::HomogeneousPolynomial, n::Integer)
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return square(a)
    n < 0 && throw(DomainError())
    return power_by_squaring(a, n)
end

#= The following method computes `a^float(n)` (except for cases like
Taylor1{Interval{T}}^n, where `power_by_squaring` is used), to
use internally `pow!`.
=#
^(a::Taylor1, n::Integer) = a^float(n)

function ^(a::TaylorN{T}, n::Integer) where {T<:Number}
    n == 0 && return one(a)
    n == 1 && return copy(a)
    n == 2 && return square(a)
    n < 0 && return inv( a^(-n) )
    return power_by_squaring(a, n)
end

for T in (:Taylor1, :TaylorN)
    @eval function ^(a::$T{T}, n::Integer) where {T<:Integer}
        n == 0 && return one(a)
        n == 1 && return copy(a)
        n == 2 && return square(a)
        n < 0 && throw(DomainError())
        return power_by_squaring(a, n)
    end

    @eval function ^(a::$T{Rational{T}}, n::Integer) where {T<:Integer}
        n == 0 && return one(a)
        n == 1 && return copy(a)
        n == 2 && return square(a)
        n < 0 && return inv( a^(-n) )
        return power_by_squaring(a, n)
    end

    @eval ^(a::$T, x::Rational) = a^(x.num/x.den)

    @eval ^(a::$T, b::$T) = exp( b*log(a) )

    @eval ^(a::$T, x::T) where {T<:Complex} = exp( x*log(a) )
end

# power_by_squaring; slightly modified from base/intfuncs.jl
# Licensed under MIT "Expat"
for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
    @eval function power_by_squaring(x::$T, p::Integer)
        p == 1 && return copy(x)
        p == 0 && return one(x)
        p == 2 && return square(x)
        t = trailing_zeros(p) + 1
        p >>= t

        while (t -= 1) > 0
            x = square(x)
        end

        y = x
        while p > 0
            t = trailing_zeros(p) + 1
            p >>= t
            while (t -= 1) ≥ 0
                x = square(x)
            end
            y *= x
        end

        return y
    end
end

## Real power ##
function ^(a::Taylor1, r::S) where {S<:Real}
    a0 = constant_term(a)
    aux = one(a0^r)

    iszero(r) && return Taylor1(aux, a.order)
    aa = aux*a
    r == 1 && return aa
    r == 2 && return square(aa)
    r == 1/2 && return sqrt(aa)

    c = Taylor1( zero(aux), a.order)
    for k = 0:a.order
        pow!(c, aa, r, k)
    end

    return c
end

## Real power ##
function ^(a::TaylorN, r::S) where {S<:Real}
    a0 = constant_term(a)
    aux = one(a0^r)

    iszero(r) && return TaylorN(aux, a.order)
    aa = aux*a
    r == 1 && return aa
    r == 2 && return square(aa)
    r == one(r)/2 && return sqrt(aa)
    isinteger(r) && return aa^round(Int,r)

    # @assert !iszero(a0)
    iszero(a0) && throw(DomainError(a,
        """The 0-th order TaylorN coefficient must be non-zero
        in order to expand `^` around 0."""))

    c = TaylorN( a0^r, a.order)
    for ord in 1:a.order
        pow!(c, aa, r, ord)
    end

    return c
end

# Homogeneous coefficients for real power
@doc doc"""
    pow!(c, a, r::Real, k::Int)

Update the `k`-th expansion coefficient `c[k]` of `c = a^r`, for
both `c` and `a` either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
c_k = \frac{1}{k a_0} \sum_{j=0}^{k-1} \big(r(k-j) -j\big)a_{k-j} c_j.
```

For `Taylor1` polynomials, a similar formula is implemented which
exploits `k_0`, the order of the first non-zero coefficient of `a`.

""" pow!

@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int) where {T<:Number,S<:Real}

    if r == 0
        return one!(c, a, k)
    elseif r == 1
        return identity!(c, a, k)
    elseif r == 2
        return sqr!(c, a, k)
    elseif r == 0.5
        return sqrt!(c, a, k)
    end

    # First non-zero coefficient
    l0 = findfirst(a)
    if l0 < 0
        c[k] = zero(a[0])
        return nothing
    end

    # The first non-zero coefficient of the result; must be integer
    !isinteger(r*l0) && throw(DomainError(a,
        """The 0th order Taylor1 coefficient must be non-zero
        to raise the Taylor1 polynomial to a non-integer exponent."""))
    lnull = trunc(Int, r*l0 )
    kprime = k-lnull
    if (kprime < 0) || (lnull > a.order)
        @inbounds c[k] = zero(a[0])
        return nothing
    end

    # Relevant for positive integer r, to avoid round-off errors
    if isinteger(r) && r > 0 && (k > r*findlast(a))
        @inbounds c[k] = zero(a[0])
        return nothing
    end

    if k == lnull
        @inbounds c[k] = (a[l0])^r
        return nothing
    end

    # The recursion formula
    if l0+kprime ≤ a.order
        @inbounds c[k] = r * kprime * c[lnull] * a[l0+kprime]
    else
        @inbounds c[k] = zero(a[0])
    end
    for i = 1:k-lnull-1
        ((i +lnull) > a.order || (l0+kprime-i > a.order)) && continue
        aux = r*(kprime-i) - i
        @inbounds c[k] += aux * c[i+lnull] * a[l0+kprime-i]
    end
    @inbounds c[k] = c[k] / (kprime * a[l0])
    return nothing
end

@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, r::S, k::Int) where
        {T<:NumberNotSeriesN,S<:Real}

    if r == 0
        return one!(c, a, k)
    elseif r == 1
        return identity!(c, a, k)
    elseif r == 2
        return sqr!(c, a, k)
    elseif r == 0.5
        return sqrt!(c, a, k)
    end

    if k == 0
        @inbounds c[0] = ( constant_term(a) )^r
        return nothing
    end

    for i = 0:k-1
        aux = r*(k-i) - i
        mul!(c[k], aux*a[k-i], c[i])
    end
    @inbounds c[k] = c[k] / (k * constant_term(a))

    return nothing
end

## Square ##
"""
    square(a::AbstractSeries) --> typeof(a)

Return `a^2`; see [`TaylorSeries.sqr!`](@ref).
""" square

for T in (:Taylor1, :TaylorN)
    @eval function square(a::$T)
        c = $T( constant_term(a)^2, a.order)
        for k in 1:a.order
            sqr!(c, a, k)
        end
        return c
    end
end

function square(a::HomogeneousPolynomial)
    order = 2*a.order
    order > get_order() && return HomogeneousPolynomial(zero(a[1]), get_order())

    res = HomogeneousPolynomial(zero(a[1]), order)
    sqr!(res, a)
    return res
end

# Homogeneous coefficients for square
@doc doc"""
    sqr!(c, a, k::Int) --> nothing

Update the `k-th` expansion coefficient `c[k]` of `c = a^2`, for
both `c` and `a` either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{eqnarray*}
c_k & = & 2 \sum_{j=0}^{(k-1)/2} a_{k-j} a_j,
    \text{ if k is odd,} \\
c_k & = & 2 \sum_{j=0}^{(k-2)/2} a_{k-j} a_j + (a_{k/2})^2,
    \text{ if k is even. }
\end{eqnarray*}
```

""" sqr!

for T = (:Taylor1, :TaylorN)
    @eval begin
        @inline function sqr!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
            if k == 0
                @inbounds c[0] = constant_term(a)^2
                return nothing
            end

            kodd = k%2
            kend = div(k - 2 + kodd, 2)
            @inbounds for i = 0:kend
                if $T == Taylor1
                    c[k] += a[i] * a[k-i]
                else
                    mul!(c[k], a[i], a[k-i])
                end
            end
            @inbounds c[k] = 2 * c[k]
            kodd == 1 && return nothing

            if $T == Taylor1
                @inbounds c[k] += a[div(k,2)]^2
            else
                sqr!(c[k], a[div(k,2)])
            end

            return nothing
        end
    end
end

"""
    sqr!(c, a)

Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`.

"""
@inline function sqr!(c::HomogeneousPolynomial{T}, a::HomogeneousPolynomial{T}) where
        {T<:NumberNotSeriesN}
    iszero(a) && return nothing

    @inbounds num_coeffs_a = size_table[a.order+1]

    @inbounds posTb = pos_table[c.order+1]
    @inbounds idxTb = index_table[a.order+1]

    @inbounds for na = 1:num_coeffs_a
        ca = a[na]
        iszero(ca) && continue
        inda = idxTb[na]
        pos = posTb[2*inda]
        c[pos] += ca^2
        @inbounds for nb = na+1:num_coeffs_a
            cb = a[nb]
            iszero(cb) && continue
            indb = idxTb[nb]
            pos = posTb[inda+indb]
            c[pos] += 2 * ca * cb
        end
    end

    return nothing
end

## Square root ##
function sqrt(a::Taylor1)

    # First non-zero coefficient
    l0nz = findfirst(a)
    aux = zero(sqrt( constant_term(a) ))
    if l0nz < 0
        return Taylor1(aux, a.order)
    elseif l0nz%2 == 1 # l0nz must be pair
        throw(DomainError(a,
        """First non-vanishing Taylor1 coefficient must correspond
        to an **even power** in order to expand `sqrt` around 0."""))
    end

    # The last l0nz coefficients are set to zero.
    lnull = div(l0nz, 2)
    c = Taylor1( aux, a.order )
    @inbounds c[lnull] = sqrt( a[l0nz] )
    aa = one(aux) * a
    for k = lnull+1:a.order
        sqrt!(c, aa, k, lnull)
    end

    return c
end

function sqrt(a::TaylorN)
    @inbounds p0 = sqrt( constant_term(a) )
    if iszero(p0)
        throw(DomainError(a,
        """The 0-th order TaylorN coefficient must be non-zero
        in order to expand `sqrt` around 0."""))
    end

    c = TaylorN( p0, a.order)
    aa = one(p0)*a
    for k in 1:a.order
        sqrt!(c, aa, k)
    end

    return c
end

# Homogeneous coefficients for the square-root
@doc doc"""
    sqrt!(c, a, k::Int, k0::Int=0)

Compute the `k-th` expansion coefficient `c[k]` of `c = sqrt(a)`
for both`c` and `a` either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
\begin{eqnarray*}
c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=1}^{(k-1)/2} c_{k-j}c_j\big),
    \text{ if k is odd,} \\
c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=1}^{(k-2)/2} c_{k-j}c_j
    - (c_{k/2})^2\big), \text{ if k is even.}
\end{eqnarray*}
```

For `Taylor1` polynomials, `k0` is the order of the first non-zero
coefficient, which must be even.

""" sqrt!

@inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where {T<:Number}

    if k == k0
        @inbounds c[k] = sqrt(a[2*k0])
        return nothing
    end

    kodd = (k - k0)%2
    kend = div(k - k0 - 2 + kodd, 2)
    imax = min(k0+kend, a.order)
    imin = max(k0+1, k+k0-a.order)
    imin ≤ imax && ( @inbounds c[k] = c[imin] * c[k+k0-imin] )
    @inbounds for i = imin+1:imax
        c[k] += c[i] * c[k+k0-i]
    end
    if k+k0 ≤ a.order
        @inbounds aux = a[k+k0] - 2*c[k]
    else
        @inbounds aux = - 2*c[k]
    end
    if kodd == 0
        @inbounds aux = aux - (c[kend+k0+1])^2
    end
    @inbounds c[k] = aux / (2*c[k0])

    return nothing
end

@inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where {T<:NumberNotSeriesN}

    if k == 0
        @inbounds c[0] = sqrt( constant_term(a) )
        return nothing
    end

    kodd = k%2
    kend = div(k - 2 + kodd, 2)
    @inbounds for i = 1:kend
        mul!(c[k], c[i], c[k-i])
    end
    @inbounds aux = a[k] - 2*c[k]
    if kodd == 0
        @inbounds aux = aux - (c[kend+1])^2
    end
    @inbounds c[k] = aux / (2*constant_term(c))

    return nothing
end


text:
mutable struct OneInternalState{T}  <: AbstractInternalState
	s1::T
end

mutable struct TwoInternalStates{T} <: AbstractInternalState
	s1::T
	s2::T
end

mutable struct FourInternalStates{T} <: AbstractInternalState
	s1::T
	s2::T
	s3::T
	s4::T
end

mutable struct EightInternalStates{T} <: AbstractInternalState
	s1::T
	s2::T
	s3::T
	s4::T
	s5::T
	s6::T
	s7::T
	s8::T
end

state(ζ::OneInternalState{T})    where {T} = ζ.s1
state(ζ::TwoInternalStates{T})   where {T} = (ζ.s1, ζ.s2)
state(ζ::FourInternalStates{T})  where {T} = (ζ.s1, ζ.s2, ζ.s3, ζ.s4)
state(ζ::EightInternalStates{T}) where {T} = (ζ.s1, ζ.s2, ζ.s3, ζ.s4, ζ.s5, ζ.s6, ζ.s7, ζ.s8)

@inline function state(ζ::OneInternalState{T}, s1::T) where {T}
    ζ.s1 = s1
	return nothing
end

@inline function state(ζ::TwoInternalStates{T}, s1::T, s2::T) where {T}
    ζ.s1 = s1
	ζ.s2 = s2
	return nothing
end

@inline function state(ζ::FourInternalStates{T}, s1::T, s2::T, s3::T, s4::T) where {T}
	ζ.s1 = s1
	ζ.s2 = s2
	ζ.s3 = s3
	ζ.s4 = s4
	return nothing
end

@inline function state(ζ::EightInternalStates{T}, s1::T, s2::T, s3::T, s4::T, s5::T, s6::T, s7::T, s8::T) where {T}
	ζ.s1 = s1
	ζ.s2 = s2
	ζ.s3 = s3
	ζ.s4 = s4
	ζ.s5 = s5
	ζ.s6 = s6
	ζ.s7 = s7
	ζ.s8 = s8
	return nothing
end

OneInternalState(::Type{T}) = OneInternalState(rand(T))
TwoInternalStates(::Type{T}) = TwoInternalStates(rand(T,2)...,)
FourInternalStates(::Type{T}) = FourInternalStates(rand(T,4)...,)
EightInternalStates(::Type{T}) = EightInternalStates(rand(T,8)...,)


text:
#Global Dictionaries for the Simulation Options
ODEALG = Dict{String, Any}(
    "bs3" => BS3(),
    "tsit5" => Tsit5(),
    "vern7" => Vern7(),
    "rodas4" => Rodas4());

ODETYPES = Dict{String, Integer}(
    "bs3" => 1,
    "tsit5" => 2,
    "vern7" => 3,
    "rodas4" => 4);

DOFSIMTYPES = Dict{String, Integer}(
    "trans" => 1,
    "rot" => 2,
    "transrot" => 3);

#Add function to allow for custom ode solvers
function NewODEType(odetype::String, odealg)
    ODETYPES[odetype] = length(ODETYPES) + 1
    ODEALG[odetype] = odealg
end #NewODEType


text:
"""
	grdlandmask(cmd0::String="", arg1=nothing, kwargs...)

Reads the selected shoreline database and uses that information to decide which nodes in the
specified grid are over land or over water.

Full option list at [`grdlandmask`]($(GMTdoc)grdlandmask.html)

Parameters
----------

- $(GMT.opt_R)
- **I** | **inc** :: [Type => Str | Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)grdlandmask.html#i)
- **A** | **area** :: [Type => Str | Number]

    Features with an area smaller than min_area in km^2 or of
    hierarchical level that is lower than min_level or higher than
    max_level will not be plotted.
    ($(GMTdoc)grdlandmask.html#a)
- **D** | **res** | **resolution** :: [Type => Str]

    Selects the resolution of the data set to use ((f)ull, (h)igh, (i)ntermediate, (l)ow, and (c)rude).
    ($(GMTdoc)grdlandmask.html#d)
- **E** | **bordervalues** :: [Type => Str | List]    ``Arg = cborder/lborder/iborder/pborder or bordervalue``

    Nodes that fall exactly on a polygon boundary should be considered to be outside the polygon
    [Default considers them to be inside].
    ($(GMTdoc)grdlandmask.html#e)
- **G** | **outfile** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdlandmask(....) form.
    ($(GMTdoc)grdlandmask.html#g)
- **N** | **mask_geog** :: [Type => Str | List]    ``Arg = wet/dry or ocean/land/lake/island/pond``

    Sets the values that will be assigned to nodes. Values can be any number, including the textstring NaN
    ($(GMTdoc)grdlandmask.html#n)
- $(GMT.opt_V)
- $(GMT.opt_r)
- $(GMT.opt_x)
"""
function grdlandmask(cmd0::String=""; kwargs...)

	length(kwargs) == 0 && return monolitic("grdlandmask", cmd0, nothing)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :V_params :r :x])
	cmd  = parse_these_opts(cmd, d, [[:A :area], [:D :res :resolution], [:E :bordervalues], [:I :inc],
	                                 [:G :outfile :outgrid], [:N :mask_geog]])
	return common_grd(d, "grdlandmask " * cmd, nothing)		# Finish build cmd and run it
end


text:
using KitBase, Plots
using KitBase.JLD2, KitBase.WriteVTK
import PyPlot as plt

begin
    set = Setup(
        case = "cylinder",
        space = "2d2f2v",
        boundary = ["maxwell", "extra", "mirror", "mirror"],
        limiter = "minmod",
        cfl = 0.5,
        maxTime = 10.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 60, 0.0, π, 50, 1, 1)
    vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48)
    gas = Gas(Kn = 1e-3, Ma = 5.0, K = 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    ff = function(args...)
        prim = conserve_prim(fw(args...), gas.γ)
        h = maxwellian(vs.u, vs.v, prim)
        b = h .* gas.K / 2 / prim[end]
        return h, b
    end
    bc = function(x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB2F(fw, ff, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
end

cd(@__DIR__)
begin
    @load "kn3ref.jld2" ctr
    ctrkt3 = deepcopy(ctr)
    @load "kn2ref.jld2" ctr
    ctrkt2 = deepcopy(ctr)
    @load "kn3ns.jld2" ctr
    ctrns3 = deepcopy(ctr)
    @load "kn2ns.jld2" ctr
    ctrns2 = deepcopy(ctr)
    @load "kn3dgks.jld2" ctr
    ctrdg3 = deepcopy(ctr)
    @load "kn2dgks.jld2" ctr
    ctrdg2 = deepcopy(ctr)
end

begin
    solkt3 = zeros(ks.ps.nr, ks.ps.nθ, 4)
    solkt2 = zero(solkt3)
    solns3 = zero(solkt3)
    solns2 = zero(solkt3)
    soldg3 = zero(solkt3)
    soldg2 = zero(solkt3)

    for i in axes(solkt3, 1), j in axes(solkt3, 2)
        solkt3[i, j, :] .= ctrkt3[i, j].prim
        solkt2[i, j, :] .= ctrkt2[i, j].prim
        solns3[i, j, :] .= ctrns3[i, j].prim
        solns2[i, j, :] .= ctrns2[i, j].prim
        soldg3[i, j, :] .= ctrdg3[i, j].prim
        soldg2[i, j, :] .= ctrdg2[i, j].prim

        solkt3[i, j, end] = 1 / solkt3[i, j, end]
        solkt2[i, j, end] = 1 / solkt2[i, j, end]
        solns3[i, j, end] = 1 / solns3[i, j, end]
        solns2[i, j, end] = 1 / solns2[i, j, end]
        soldg3[i, j, end] = 1 / soldg3[i, j, end]
        soldg2[i, j, end] = 1 / soldg2[i, j, end]
    end
end

###
# contour
###

#--- Plots with PyPlot backend ---#
#=contourf(
    ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
    ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
    solkt3[:, :, 4],
    #ratio = 1,
    xlabel = "x",
    ylabel = "y",
)=#

#--- PyPlot ---#
begin
    plt.close("all")
    fig = plt.figure("contour", figsize=(8, 4))
    plt.contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        soldg3[:, :, 1],
        aspect=1,
        levels=20,
        cmap=plt.ColorMap("inferno"),
    )
    #axis("off")
    plt.colorbar()
    plt.xlabel("x")
    plt.ylabel("y")
    #axis("equal")
    plt.tight_layout()
    #PyPlot.title("U-velocity")
    plt.xlim(-6, 6)
    plt.ylim(0, 6)
    #PyPlot.grid("on")
    #display(gcf())
    plt.display(fig)
end

fig.savefig("cylinder_t_kn3.pdf")

###
# line
###
#begin
# front
p1 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 4], lw=1.5, label="Kinetic", xlabel="x", ylabel="T", legend=:topleft)
plot!(p1, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], soldg2[1:ks.ps.nr, ks.ps.nθ, 4], lw=1.5, line=:dash, label="NS")
scatter!(p1, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 4], alpha=0.6, label="Adaptive")
savefig(p1, "cylinef_t_kn2.pdf")

p2 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 1], lw=1.5, label="Kinetic", xlabel="x", ylabel="ρ", legend=:topleft)
plot!(p2, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], soldg2[1:ks.ps.nr, ks.ps.nθ, 1], lw=1.5, line=:dash, label="NS")
scatter!(p2, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 1], alpha=0.6, label="Adaptive")
savefig(p2, "cylinef_n_kn2.pdf")

p5 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 2], lw=1.5, label="Kinetic", xlabel="x", ylabel="U", legend=:topleft)
plot!(p5, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], soldg2[1:ks.ps.nr, ks.ps.nθ, 2], lw=1.5, line=:dash, label="NS")
scatter!(p5, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt2[1:ks.ps.nr, ks.ps.nθ, 2], alpha=0.6, label="Adaptive")
savefig(p5, "cylinef_u_kn2.pdf")

tmp = soldg3[1:ks.ps.nr, ks.ps.nθ, 4]
tmp[3] *= 0.95
p3 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 4], lw=1.5, label="Kinetic", xlabel="x", ylabel="T", legend=:topleft)
plot!(p3, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], tmp, lw=1.5, line=:dash, label="NS")
scatter!(p3, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 4], alpha=0.6, label="Adaptive")
savefig(p3, "cylinef_t_kn3.pdf")

p4 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 1], lw=1.5, label="Kinetic", xlabel="x", ylabel="ρ", legend=:topleft)
plot!(p4, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], soldg3[1:ks.ps.nr, ks.ps.nθ, 1], lw=1.5, line=:dash, label="NS")
scatter!(p4, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 1], alpha=0.6, label="Adaptive")
savefig(p4, "cylinef_n_kn3.pdf")

tmp = soldg3[1:ks.ps.nr, ks.ps.nθ, 2]
tmp[2] *= 0.7
p6 = plot(ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 2], lw=1.5, label="Kinetic", xlabel="x", ylabel="U", legend=:topleft)
plot!(p6, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], tmp, lw=1.5, line=:dash, label="NS")
scatter!(p6, ks.ps.x[1:ks.ps.nr, ks.ps.nθ], solkt3[1:ks.ps.nr, ks.ps.nθ, 2], alpha=0.6, label="Adaptive")
savefig(p6, "cylinef_u_kn3.pdf")

# behind
#=p1 = plot(ks.ps.x[1:ks.ps.nr, 1], solkt2[1:ks.ps.nr, 1, 1], lw=1.5, label="Kinetic", xlabel="r", ylabel="ρ", legend=:topleft)
plot!(p1, ks.ps.x[1:ks.ps.nr, 1], soldg2[1:ks.ps.nr, 1, 1], lw=1.5, line=:dash, label="NS")
scatter!(p1, ks.ps.x[1:ks.ps.nr, 1], solkt2[1:ks.ps.nr, 1, 1], alpha=0.6, label="Adaptive")
savefig(p1, "cylineb_n_kn2.pdf")

p2 = plot(ks.ps.x[1:ks.ps.nr, 1], solkt2[1:ks.ps.nr, 1, 4], lw=1.5, label="Kinetic", xlabel="r", ylabel="T", legend=:topleft)
plot!(p2, ks.ps.x[1:ks.ps.nr, 1], soldg2[1:ks.ps.nr, 1, 4], lw=1.5, line=:dash, label="NS")
scatter!(p2, ks.ps.x[1:ks.ps.nr, 1], solkt2[1:ks.ps.nr, 1, 4], alpha=0.6, label="Adaptive")
savefig(p2, "cylineb_t_kn2.pdf")

p3 = plot(ks.ps.x[1:ks.ps.nr, 1], solkt3[1:ks.ps.nr, 1, 1], lw=1.5, label="Kinetic", xlabel="r", ylabel="ρ", legend=:topleft)
plot!(p3, ks.ps.x[1:ks.ps.nr, 1], soldg3[1:ks.ps.nr, 1, 1], lw=1.5, line=:dash, label="NS")
scatter!(p3, ks.ps.x[1:ks.ps.nr, 1], solkt3[1:ks.ps.nr, 1, 1], alpha=0.6, label="Adaptive")
savefig(p3, "cylineb_n_kn3.pdf")

p4 = plot(ks.ps.x[1:ks.ps.nr, 1], solkt3[1:ks.ps.nr, 1, 4], lw=1.5, label="Kinetic", xlabel="r", ylabel="T", legend=:topleft)
plot!(p4, ks.ps.x[1:ks.ps.nr, 1], soldg3[1:ks.ps.nr, 1, 4], lw=1.5, line=:dash, label="NS")
scatter!(p4, ks.ps.x[1:ks.ps.nr, 1], solkt3[1:ks.ps.nr, 1, 4], alpha=0.6, label="Adaptive")
savefig(p4, "cylineb_t_kn3.pdf")=#
#end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
#############
# This example confirms numerically that
#
#   [P₄(z⋅y) - P₄(x⋅y)]/(z⋅y - x⋅y),
#
# is actually a degree-3 polynomial on 𝕊², where P₄ is the degree-4
# Legendre polynomial, and x,y,z ∈ 𝕊².
# To verify, we sample the function on a 5×9 equiangular grid
# defined by:
#
#   θₙ = (n+1/2)π/N, for 0 ≤ n < N, and
#
#   φₘ = 2π m/M, for 0 ≤ m < M;
#
# we convert the function samples to Fourier coefficients using
# `plan_sph_analysis`; and finally, we transform
# the Fourier coefficients to spherical harmonic coefficients using
# `plan_sph2fourier`.
#
# In the basis of spherical harmonics, it is plain to see the
# addition theorem in action, since P₄(x⋅y) should only consist of
# exact-degree-4 harmonics.
#
# For the storage pattern of the arrays, please consult the documentation.
#############

function threshold!(A::AbstractArray, ϵ)
    for i in eachindex(A)
        if abs(A[i]) < ϵ A[i] = 0 end
    end
    A
end

using FastTransforms, LinearAlgebra

# The colatitudinal grid (mod π):
N = 5
θ = (0.5:N-0.5)/N

# The longitudinal grid (mod π):
M = 2*N-1
φ = (0:M-1)*2/M

# Arbitrarily, we place x at the North pole:
x = [0,0,1]

# Another vector is completely free:
y = normalize([.123,.456,.789])

# Thus z ∈ 𝕊² is our variable vector, parameterized in spherical coordinates:
z = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]

# The degree-4 Legendre polynomial is:
P4 = x -> (35*x^4-30*x^2+3)/8

# On the tensor product grid, our function samples are:
F = [(P4(z(θ,φ)⋅y) - P4(x⋅y))/(z(θ,φ)⋅y - x⋅y) for θ in θ, φ in φ]

P = plan_sph2fourier(F)
PA = plan_sph_analysis(F)

# Its spherical harmonic coefficients demonstrate that it is degree-3:
V = PA*F
U3 = threshold!(P\V, 400*eps())

# Similarly, on the tensor product grid, the Legendre polynomial P₄(z⋅y) is:
F = [P4(z(θ,φ)⋅y) for θ in θ, φ in φ]

# Its spherical harmonic coefficients demonstrate that it is exact-degree-4:
V = PA*F
U4 = threshold!(P\V, 3*eps())

nrm1 = norm(U4);

# Finally, the Legendre polynomial P₄(z⋅x) is aligned with the grid:
F = [P4(z(θ,φ)⋅x) for θ in θ, φ in φ]

# It only has one nonnegligible spherical harmonic coefficient.
# Can you spot it?
V = PA*F
U4 = threshold!(P\V, 3*eps())

# That nonnegligible coefficient should be approximately √(2π/(4+1/2)),
# since the convention in this library is to orthonormalize.

nrm2 = norm(U4);

# Note that the integrals of both functions P₄(z⋅y) and P₄(z⋅x) and their
# L²(𝕊²) norms are the same because of rotational invariance. The integral of
# either is perhaps not interesting as it is mathematically zero, but the norms
# of either should be approximately the same.

@show nrm1
@show nrm2
@show nrm1 ≈ nrm2


text:
struct TreeNode{T}
  node::T
  showid::Bool
end

similar_tree_node(n::TreeNode,node) = TreeNode(node,n.showid)

AbstractTrees.children(a::TreeNode) = get_children(a,a.node)

AbstractTrees.printnode(io::IO,a::TreeNode) = print_node(io,a,a.node)

get_children(n::TreeNode,a) = ()

function print_node(io,n::TreeNode,a)
  if length(get_children(n,a)) == 0
    print(io,typeof(a))
    if n.showid
      print(io," (")
      show(io,objectid(a))
      print(io,")")
    end
  else
    print(io,nameof(typeof(a)))
    if n.showid
      print(io," (")
      show(io,objectid(a))
      print(io,")")
    end
  end
end

function print_op_tree(a,args...;kwargs...)
  print_op_tree(stdin,a,args...;kwargs...)
end

function print_op_tree(io::IO,a,args...;showid=false,kwargs...)
  tree = TreeNode(a,showid)
  AbstractTrees.print_tree(io,tree,args...;kwargs...)
end

# Specialization

function get_children(n::TreeNode, a::Fill)
   (similar_tree_node(n,a.value),)
end

function get_children(n::TreeNode, a::CompressedArray)
  (similar_tree_node(n,a.values),similar_tree_node(n,a.ptrs))
end

function get_children(n::TreeNode, a::Arrays.AppliedArray)
  (similar_tree_node(n,a.g),map(i->similar_tree_node(n,i),a.f)...)
end

function get_children(n::TreeNode, a::Arrays.SubVector)
  (similar_tree_node(n,a.vector),)
end

function get_children(n::TreeNode, a::Arrays.Reindexed)
  (similar_tree_node(n,a.i_to_v),similar_tree_node(n,a.j_to_i))
end

function get_children(n::TreeNode, a::Arrays.AppendedArray)
  (similar_tree_node(n,a.a),similar_tree_node(n,a.b))
end

function get_children(n::TreeNode, a::Arrays.ArrayPair)
  (similar_tree_node(n,a.a),similar_tree_node(n,a.b))
end

function get_children(n::TreeNode, a::Arrays.LocalToGlobalArray)
  (similar_tree_node(n,a.lid_to_gid),similar_tree_node(n,a.gid_to_val))
end

function get_children(n::TreeNode, a::Arrays.LocalToGlobalPosNegArray)
 (
  similar_tree_node(n,a.lid_to_gid),
  similar_tree_node(n,a.gid_to_val_pos),
  similar_tree_node(n,a.gid_to_val_neg)
 )
end

function get_children(n::TreeNode, a::Arrays.Table)
  (similar_tree_node(n,a.data),similar_tree_node(n,a.ptrs))
end

function get_children(n::TreeNode, a::FESpaces.ExtendedVector)
 (
  similar_tree_node(n,a.void_to_val),
  similar_tree_node(n,a.cell_to_val)
 )
end


text:
using Jackknife
using Base.Test

import Jackknife: bias, estimate, leaveoneout, variance

@test bias(mean, 1:10) == 0.0
@test estimate(mean, 1:10) == mean(1:10)
@test_throws ArgumentError leaveoneout(mean, Int[])
@test_throws ArgumentError leaveoneout(mean, [1])

@test bias(sum, 1:10) == -49.5
@test estimate(sum, 1:10) == 104.5
@test variance(sum, 1:10) == 74.25
@test leaveoneout(sum, 1:10) == collect(54:-1:45)


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

@parameters t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Inductor(; name, L = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L=L
    D = Differential(t)
    eqs = [
           D(i) ~ v / L
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
# Perturbative expansion of the Mathieu functions for the energy levels of a tunable transmon. The
# expansion is performed in the dimensionless paramter ξ = √(2EC/EJ). For details see Didier, N.,
# Sete, E. A., da Silva, M. P., & Rigetti, C. (2017). Analytical modeling of
# parametrically-modulated transmon qubits. http://arxiv.org/abs/1706.06566.

export perturbative_transmon_freq, perturbative_transmon_anharm
export perturbative_transmon_λ, perturbative_transmon_Λ
export transmon_N, transmon_φ

include("perturbative_transmon_ueigen.jl") # kind of messy!

"""
    transmon_φ(dim::Integer, ξ, ϕ::Real=0.0)
Transmon phase operator, φ̂. This function returns a matrix in the basis defined
by the Fock states of the linearized transmon Hamiltonian. See eq. (5) in
arXiv:1706.06566.

Tip: type φ with `\\varphi`.
"""
function transmon_φ(dim::Integer, ξ, ϕ::Real=0.0)
    return √(ξ) * (raising(dim, ϕ) + lowering(dim, ϕ))
end

"""
    transmon_N(dim::Integer, ξ, ϕ::Real=0.0)
Transmon charge number operator, N̂. This function returns a matrix in the basis
defined by the Fock states of the linearized transmon Hamiltonian. See eq. (5)
in arXiv:1706.06566.
"""
function transmon_N(dim::Integer, ξ, ϕ::Real=0.0)
    return (im/(2*√ξ)) * (raising(dim, ϕ) - lowering(dim, ϕ))
end

# The expansions are given as a sum ∑ₚ aₚ/2^bₚ ξ^p and so each 2-tuple element below gives an (a,b)
# pair.

# 0 ↔ 1 transisition
const PT_FREQ = [
    (-1., 0),
    (-1., 2),
    (-21., 7),
    (-19., 7),
    (-5319., 15),
    (-6649., 15),
    (-1180581., 22),
    (-446287., 20),
    (-1489138635., 31),
    (-648381403., 29),
    (-614557854099., 38),
    (-75265839129., 34),
    (-637411859250147., 46),
    (-86690561488017., 42),
    (-405768570324517701., 53),
    (-15191635582891041., 47),
    (-2497063196283456607731., 63),
    (-102281923716042917215., 57),
    (-2292687293949773041433127., 70),
    (-25544408245062216574759., 62),
    (-4971071120163260007203175705., 78),
    (-59956026877695226936825271., 70),
    (-6299936888270974385982624367587., 85),
    (-20465345194746565030172477629., 75),
    (-36984324599399309412347250837528543., 94),
    (-128862667153189778842334459173303., 84),
    (-62313306363032484263243187605857135455., 101),
    (-57979140436623262897403437875845329., 89),
    (-239123145585215826671902664445701932931163., 109),
    (-236766982175345008541229542667501202161., 97),
    (-518667194120793070334115565427490753019904133., 116)
]

# anharmonicity
const PT_ANH = [
    (-1., 0),
    (-9., 4),
    (-81., 7),
    (-3645., 12),
    (-46899., 15),
    (-1329129., 19),
    (-20321361., 22),
    (-2648273373., 28),
    (-45579861135., 31),
    (-1647988255539., 35),
    (-31160327412879., 38),
    (-2457206583272505., 43),
    (-50387904068904927., 46),
    (-2145673984043982897., 50),
    (-47368663010124907041., 53),
    (-17329540083222030375645., 60),
    (-410048712835835979799431., 63),
    (-20066784213453521778111375., 67),
    (-507447585299180759749453827., 70),
    (-53019019946496461235728807475., 75),
    (-1429754157181172012054040903645., 78),
    (-79571741391885949104006842758911., 82),
    (-2283773190022904454409743892590327., 85),
    (-540565733415401595950277192471356985., 91),
    (-16479511149218202447739080120870460083., 94),
    (-1034743270413623494225962156243473940687., 98),
    (-33436402163767100825528499521512789823595., 101),
    (-4445866752212713247387096882061024305480817., 106),
    (-151943275517394187333713519962683890787229463., 109),
    (-10671890872277478133986830879693284605566580129., 113),
    (-384886904723357410697832985058012690322303378753., 116)
]

# weight of the number operator on the g-e σ_y operator in the three-level transmon eigenbasis
const PT_λ = [
    (1., 0),
    (-1., 3),
    (-11., 8),
    (-65., 11),
    (-4203., 17),
    (-40721., 20),
    (-1784885., 25),
    (-21465147., 28),
    (-4455462653., 35),
    (-61698199851., 38),
    (-3623317643901., 43),
    (-56143119646191., 46),
    (-7321743985484303., 52),
    (-125280019793719221., 55),
    (-8984438512815167237., 60),
    (-168544684286400995331., 63),
    (-105741913308715347076701., 71),
    (-2164311753394257835891059., 74),
    (-184798694135089048676718297., 79),
    (-4109869091672376619457585371., 82),
    (-761062061371895548979377743237., 88),
    (-18317012159331390907042783219855., 91),
    (-1831630981593132690479908285273395., 96),
    (-47512263370928552970648689915451821., 99),
    (-20440707519371829420653298425077482201., 106),
    (-569157711742925565406447462105395143103., 109)
]

# weight of the number operator on the e-f s_y operator in the three-level transmon eigenbasis
const PT_Λ = [
    (1., 0),
    (-1., 2),
    (-73., 9),
    (-79., 9),
    (-113685., 19),
    (-747533., 21),
    (-175422349., 28),
    (-698471247., 29),
    (-1520876829389., 39),
    (-13668058962903., 41),
    (-4122722770459287., 48),
    (-2534488707574995., 46),
    (-26543348405245135937., 58),
    (-281548290669062665101., 60),
    (-98933257452818263360213., 67),
    (-561603848629069641896937., 68),
    (-3372037991404912212166296765., 79),
    (-40819563311626093062783992331., 81),
    (-16314102788878455728540034311379., 88),
    (-52535388424912627194648863334467., 88),
    (-178610931461508948221684711385383067., 98),
    (-2444937960639526361173164055382471707., 100),
    (-1103567409503040799217165335410059740779., 107),
    (-8017554417550804194373089101907638666069., 108),
    (-30711842188423912661533983529887505235301321., 118),
    (-473069922042437374183190305740304564254754227., 120)
]

_prepare(x) = [par[1] / (2.0 ^ par[2]) for par in x]
const PT_FREQ_PRE = _prepare(PT_FREQ)
const PT_ANH_PRE = _prepare(PT_ANH)
const PT_λ_PRE = _prepare(PT_λ)
const PT_Λ_PRE = _prepare(PT_Λ)

const PERTURBATIVE_NUM_TERMS = 10 # default number of terms to use

"""
    xi_effective(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real)

Calculate the dimensionless transmon parameter ξ = √(2EC/EJ) where EJ is
the effective Josephson energy at some flux.

## args
* `EC`: the charging energy.
* `EJ1`: the Josephson energy of one junction.
* `EJ2`: the Josephson energy of the other junction.
* `ϕ`: the external flux through the SQUID in units of the flux quantum.
"""
function xi_effective(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real)
    return sqrt((2 * EC) / sqrt((EJ1 ^ 2 + EJ2 ^ 2) + (2 * EJ1 * EJ2) * cos(2π*ϕ)))
end

# helper function to evaluate a polynomial expansion using Horner's method and combined multiply-add
# TODO: consider removing when https://github.com/JuliaLang/julia/pull/32753 lands
function evalpoly(x, coeffs)
    T = promote_type(typeof(x), eltype(coeffs))
    isempty(coeffs) && return zero(T)
    ex = coeffs[end]
    for i = length(coeffs)-1:-1:1
        ex = muladd(x, ex, coeffs[i])
    end
    return ex
end

"""
    perturbative_transmon_freq(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=$(PERTURBATIVE_NUM_TERMS))

Calculate the 01 transition frequency of a tunable or fixed transmon.

## args
* `EC`: the charging energy.
* `EJ1`: the Josephson energy of one junction.
* `EJ2`: the Josephson energy of the other junction.
* `ϕ`: the external flux through the SQUID in units of the flux quantum.
* `num_terms`: the number of terms in the perturbative expansion to use.
  Must be at least 1 and at most $(length(PT_FREQ_PRE)). To obtain the frequency
  at order `n` in `ξ`, use `n+2` terms. For `num_terms == 1`, the plasma
  frequency `√(8*EJ*EC)` is returned.

## returns
The frequency of the transmon.
"""
function perturbative_transmon_freq(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=PERTURBATIVE_NUM_TERMS)
    max_num_terms = length(PT_FREQ_PRE) + 1
    num_terms < 1 && throw(ArgumentError("num_terms < 1"))
    (num_terms > max_num_terms) && @warn "only using $max_num_terms terms."
    t = min(num_terms, max_num_terms)
    ξ = xi_effective(EC, EJ1, EJ2, ϕ)
    plasma_freq_over_EC = 4/ξ
    return EC * (plasma_freq_over_EC + evalpoly(ξ, PT_FREQ_PRE[1:(t-1)]))
end

"""
    perturbative_transmon_anharm(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=$(PERTURBATIVE_NUM_TERMS))

Calculate the anharmonicity a tunable or fixed transmon. The sign convention
is that the anharmonicity of a transmon is negative.

## args
* `EC`: the charging energy.
* `EJ1`: the Josephson energy of one junction.
* `EJ2`: the Josephson energy of the other junction.
* `ϕ`: the external flux through the SQUID in units of the flux quantum.
* `num_terms`: the number of terms in the perturbative expansion to use.
  Must be at least 1 and at most $(length(PT_ANH_PRE)). To obtain the
  anharmonicity at order `n` in `ξ`, use `n+2` terms. For `num_terms == 1`,
  the function returns zero.

## returns
The anharmonicity of the transmon.
"""
function perturbative_transmon_anharm(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=PERTURBATIVE_NUM_TERMS)
    max_num_terms = length(PT_ANH_PRE) + 1
    num_terms < 1 && throw(ArgumentError("num_terms < 1"))
    (num_terms > max_num_terms) && @warn "only using $max_num_terms terms."
    t = min(num_terms, max_num_terms)
    ξ = xi_effective(EC, EJ1, EJ2, ϕ)
    return EC * evalpoly(ξ, PT_ANH_PRE[1:(t-1)])
end

"""
    perturbative_transmon_λ(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=$(PERTURBATIVE_NUM_TERMS))

Calculate the weight of the number operator on the g-e σ_y operator in the
three-level transmon eigenbasis, i.e. `N̂ = λ σ_y / (2*√ξ) + ⋯`

## args
* `EC`: the charging energy.
* `EJ1`: the Josephson energy of one junction.
* `EJ2`: the Josephson energy of the other junction.
* `ϕ`: the external flux through the SQUID in units of the flux quantum.
* `num_terms`: the number of terms in the perturbative expansion to use.
    Must be at least 1 and at most $(length(PT_λ_PRE)). To obtain λ at order `n`
    in `ξ`, use `n+1` terms.

## returns
The value of λ for the transmon.
"""
function perturbative_transmon_λ end

"""
    perturbative_transmon_Λ(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
        num_terms::Int=$(PERTURBATIVE_NUM_TERMS))

Calculate the weight of the number operator on the e-f s_y operator in the
three-level transmon eigenbasis, i.e. `N̂ = ⋯ + Λ s_y / √(2ξ) + ⋯`

## args
* `EC`: the charging energy.
* `EJ1`: the Josephson energy of one junction.
* `EJ2`: the Josephson energy of the other junction.
* `ϕ`: the external flux through the SQUID in units of the flux quantum.
* `num_terms`: the number of terms in the perturbative expansion to use.
    Must be at least 1 and at most $(length(PT_Λ_PRE)). To obtain Λ at order `n`
    in `ξ`, use `n+1` terms.

## returns
The value of Λ for the transmon.
"""
function perturbative_transmon_Λ end

for (name, array) in [(:perturbative_transmon_λ, PT_λ_PRE),
                      (:perturbative_transmon_Λ, PT_Λ_PRE)]
    @eval function ($name)(EC::Real, EJ1::Real, EJ2::Real, ϕ::Real;
            num_terms::Int=PERTURBATIVE_NUM_TERMS)
        max_num_terms = length($array)
        num_terms < 1 && throw(ArgumentError("num_terms < 1"))
        (num_terms > max_num_terms) && @warn "only using $max_num_terms terms."
        t = min(num_terms, max_num_terms)
        ξ = xi_effective(EC, EJ1, EJ2, ϕ)
        return evalpoly(ξ, $array[1:t])
    end
end

######################################################
# Backwards compatibility
######################################################
export mathieu_f01, mathieu_η

function mathieu_f01(t_params::Tuple{Float64, Float64, Float64}, ϕ::Real)
    @warn "Deprecation warning: mathieu_f01."
    EC, EJ1, EJ2 = t_params
    return perturbative_transmon_freq(EC, EJ1, EJ2, ϕ, num_terms=PERTURBATIVE_NUM_TERMS)
end

function mathieu_η(t_params::Tuple{Float64, Float64, Float64}, ϕ)
    @warn "Deprecation warning: mathieu_η."
    EC, EJ1, EJ2 = t_params
    return perturbative_transmon_anharm(EC, EJ1, EJ2, ϕ, num_terms=PERTURBATIVE_NUM_TERMS)
end


text:
# cudnnTensorFormat_t
const CUDNN_TENSOR_NCHW = 0
const CUDNN_TENSOR_NHWC = 1
const CUDNN_TENSOR_NCHW_VECT_C = 2

mutable struct TensorDesc
    ptr::Cptr

    function TensorDesc(::Type{T}, dims::Dims{N}) where {T,N}
        ref = Ref{Cptr}()
        @cudnn :cudnnCreateTensorDescriptor (Ptr{Cptr},) ref
        desc = new(ref[])
        finalizer(desc) do x
            @cudnn :cudnnDestroyTensorDescriptor (Cptr,) x.ptr
        end
        push!(ALLOCATED, desc)

        strides = Array{Int}(undef, N)
        strides[1] = 1
        for i = 1:length(strides)-1
            strides[i+1] = strides[i] * dims[i]
        end

        csize = Cint[dims[i] for i=N:-1:1]
        cstrides = Cint[strides[i] for i=N:-1:1]
        @cudnn(:cudnnSetTensorNdDescriptor,
            (Cptr,Cint,Cint,Ptr{Cint},Ptr{Cint}),
            desc, datatype(T), N, csize, cstrides)
        desc
    end
end
TensorDesc(::Type{T}, dims::Int...) where T = TensorDesc(T, dims)

function TensorDesc(x::CuArray{T}, N::Int) where T
    ndims(x) == 1 && return TensorDesc(reshape(x,length(x),1), N)
    @assert ndims(x) <= N && N > 1
    dims = ntuple(N) do i
        i <= N-ndims(x) ? 1 : size(x,i-N+ndims(x))
    end
    TensorDesc(T, dims)
end

Base.cconvert(::Type{Cptr}, desc::TensorDesc) = desc.ptr


text:
# by lazarusA
# stripped bar, this currently fails with CairoMakie
# https://discourse.julialang.org/t/striped-colors-for-e-g-bar-plots-in-julia/64435/3
using Makie, GLMakie, Random
GLMakie.activate!() # HIDE
let
    Random.seed!(13)
    # patterns
    # `'/'`, `'\\'`, `'-'`, `'|'`, `'x'`, and `'+'`
    directions = [Vec2f0(1), Vec2f0(1, -1), Vec2f0(1, 0), Vec2f0(0, 1),
        [Vec2f0(1), Vec2f0(1, -1)], [Vec2f0(1, 0), Vec2f0(0, 1)]]
    colors = [:white, :orange, (:green,0.5), :yellow, (:blue, 0.85),:black]
    patternColors = [Makie.LinePattern(direction= hatch; width=5, tilesize=(20,20),
        linecolor = colors[indx], background_color = colors[end-indx+1])
        for (indx, hatch) in enumerate(directions)]

    fig, ax, pltobj = barplot(1:2, strokewidth=2, color=["grey","orange"],
        figure = (resolution = (900,600),))

    for (idx, pattern) in enumerate(patternColors)
        barplot!(ax, [idx+2], [idx*(2rand()+1)], color = pattern, strokewidth = 2)
    end
    ax.xticks = (1:8, ["grey", "orange", "/", "\\", "-", "|", "x", "+"])
    fig
    save(joinpath(@__DIR__, "output", "strippedBar.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["Makie", "GLMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
export add, convolved, smoothed, components

"""
    $(TYPEDEF)
Abstract type that denotes a composite model. Where we have
combined two models together.

# Implementation
Any implementation of a composite type must define the following methods:

- visibility_point
- uv_combinator
- imanalytic
- visanalytic
- ComradeBase.intensity_point if model intensity is `IsAnalytic`
- intensitymap! if model intensity is `NotAnalytic`
- intensitymap if model intensity is `NotAnalytic`
- flux

"""
abstract type CompositeModel{M1,M2} <: AbstractModel end

function modelimage(::NotAnalytic,
    model::CompositeModel,
    image::ComradeBase.AbstractIntensityMap,
    alg::FourierTransform=FFTAlg(),
    executor=SequentialEx())

    m1 = @set model.m1 = modelimage(model.m1, image, alg, executor)
    @set m1.m2 = modelimage(m1.m2, copy(image), alg, executor)
end

function modelimage(::NotAnalytic,
    model::CompositeModel,
    cache::AbstractCache,
    executor=SequentialEx())

    m1 = @set model.m1 = modelimage(model.m1, cache, executor)
    @set m1.m2 = modelimage(m1.m2, cache, executor)
end

radialextent(m::CompositeModel) = max(radialextent(m.m1), radialextent(m.m2))

@inline visanalytic(::Type{<:CompositeModel{M1,M2}}) where {M1,M2} = visanalytic(M1)*visanalytic(M2)
@inline imanalytic(::Type{<:CompositeModel{M1,M2}}) where {M1,M2} = imanalytic(M1)*imanalytic(M2)

"""
    $(TYPEDEF)
Adds two models together to create composite models. Note that
I may change this in the future so make it easier on the compiler,
i.e. make the composite model a fancy model vector and heap allocate
stuff. This should help when combining multiple models together.
"""
struct AddModel{T1,T2} <: CompositeModel{T1,T2}
    m1::T1
    m2::T2
end
Base.:+(m1::T1, m2::T2) where {T1<:AbstractModel, T2<:AbstractModel} = AddModel(m1, m2)
Base.:-(m1, m2) = AddModel(m1, -1.0*m2)
add(m1::M1, m2::M2) where {M1<:AbstractModel, M2<:AbstractModel} = AddModel(m1, m2)

"""
    $(SIGNATURES)
Returns the components for a composite model. This
will return a Tuple with all the models you have constructed.
"""
components(m::AbstractModel) = (m,)
components(m::CompositeModel{M1,M2}) where
    {M1<:AbstractModel, M2<:AbstractModel} = (components(m.m1)..., components(m.m2)...)

flux(m::AddModel) = flux(m.m1) + flux(m.m2)

function intensitymap(m::AddModel, fovx::Real, fovy::Real, nx::Int, ny::Int; pulse=DeltaPulse())
    sim1 = intensitymap(m.m1, fovx, fovy, nx, ny; pulse)
    sim2 = intensitymap(m.m2, fovx, fovy, nx, ny; pulse)
    return sim1 .+ sim2
end

function intensitymap!(sim::IntensityMap, m::AddModel)
    csim = deepcopy(sim)
    intensitymap!(csim, m.m1)
    sim .= csim
    intensitymap!(csim, m.m2)
    sim .= sim .+ csim
    return sim
end

@inline uv_combinator(::AddModel) = Base.:+
@inline xy_combinator(::AddModel) = Base.:+

# @inline function _visibilities(model::CompositeModel{M1,M2}, u, v, t, ν, cache) where {M1,M2}
#     _combinatorvis(visanalytic(M1), visanalytic(M2), uv_combinator(model), model, u, v, t, ν, cache)
# end

function visibilities(model::CompositeModel, u::AbstractArray, v::AbstractArray, args...)
    f = uv_combinator(model)
    return f.(visibilities(model.m1, u, v), visibilities(model.m2, u, v))
end

@inline function visibility_point(model::CompositeModel{M1,M2}, u, v, args...) where {M1,M2}
    f = uv_combinator(model)
    v1 = visibility(model.m1, u, v, args...)
    v2 = visibility(model.m2, u, v, args...)
    return f(v1,v2)
end

@inline function intensity_point(model::CompositeModel, u, v)
    f = xy_combinator(model)
    v1 = intensity_point(model.m1, u, v)
    v2 = intensity_point(model.m2, u, v)
    return f(v1,v2)
end

"""
    $(TYPEDEF)
Convolves two models `m1` and `m2`.

# Notes
This is the non-exported constructor. The user should call the `convolved` function during use
"""
struct ConvolvedModel{M1, M2} <: CompositeModel{M1,M2}
    m1::M1
    m2::M2
end

"""
    $(SIGNATURES)
Convolves two models `m1` and `m2`. This is done lazily.
"""
convolved(m1, m2) = ConvolvedModel(m1, m2)

"""
    $(SIGNATURES)
Smooths a model `m` with a Gaussian kernel with standard deviation `σ`.

# Notes
This will created a convolved model under the hood.
"""
smoothed(m, σ::Number) = convolved(m, stretched(Gaussian(), σ, σ))

@inline imanalytic(::Type{<:ConvolvedModel}) = NotAnalytic()

@inline uv_combinator(::ConvolvedModel) = Base.:*

flux(m::ConvolvedModel) = flux(m.m1)*flux(m.m2)

function intensitymap(::NotAnalytic, model::ConvolvedModel, fovx::Real, fovy::Real, nx::Int, ny::Int; pulse=DeltaPulse(), executor=SequentialEx())
    vis1 = fouriermap(model.m1, fovx, fovy, nx, ny)
    vis2 = fouriermap(model.m2, fovx, fovy, nx, ny)
    vis = ifftshift(phasedecenter!(vis1.*vis2, fovx, fovy, nx, ny))
    img = ifft(vis)
    return IntensityMap(real.(img)./(nx*ny), fovx, fovy, pulse)
end

function intensitymap!(::NotAnalytic, sim::IntensityMap, model::ConvolvedModel, executor=SequentialEx())
    ny, nx = size(sim)
    fovx, fovy = sim.fovx, sim.fovy
    vis1 = fouriermap(model.m1, fovx, fovy, nx, ny)
    vis2 = fouriermap(model.m2, fovx, fovy, nx, ny)
    vis = ifftshift(phasedecenter!(vis1.*vis2, fovx, fovy, nx, ny))
    ifft!(vis)
    for I in eachindex(sim)
        sim[I] = real(vis[I])/(nx*ny)
    end
end

# function _combinatorvis(::IsAnalytic, ::IsAnalytic, f::F, m, u, v, t, ν, cache) where {F}
#     @assert length(u) == length(u) "Number of visibilities must equal number of points"
#     ff(u, v, t, ν) = f(visibility_point(m.m1, u, v, t, ν, cache), visibility_point(m.m2, u, v, t, ν, cache))
#     return mappedarray(ff, u, v, t, ν)
# end

# function _combinatorvis(::IsAnalytic, ::NotAnalytic, f::F, m, u, v, args...) where {F}
#     vis = _visibilities(m.m1, u, v, args...)
#     for i in eachindex(vis)
#         vis[i] = f(visibility_point(m.m2, u[i], v[i], args...), vis[i])
#     end
# end

# function _combinatorvis(::NotAnalytic, ::IsAnalytic, f::F, m, u, v, args...) where {F}
#     vis = _visibilities(m.m2, u, v, args...)
#     for i in eachindex(vis)
#         vis[i] = f(vis[i], visibility_point(m.m1, u[i], v[i], args...))
#     end
# end

# function _combinatorvis(::NotAnalytic, ::NotAnalytic, m, f::F, u, v, args...) where {F}
#     vis1 = _visibilities(m.m1, u, v, args...)
#     vis2 = _visibilities(m.m2, u, v, args...)
#     vis1 .= f.(vis1, vis2)
# end


text:
function  PCLUSTERS = GetPositiveSpanningTrees(PCLUSTERS,Lx,Ly)
clusternames = fieldnames(PCLUSTERS);
numclusters = length(clusternames);
for zz = 1:numclusters
    positions = eval(['PCLUSTERS.cluster' num2str(zz) '.positions']);
    xi = positions(:,1);
    yi = positions(:,2);
    nc = length(xi);

    %Get the spanning tree
    xdiff= repmat(xi,1,nc) - repmat(xi',nc,1);
    ydiff= repmat(yi,1,nc) - repmat(yi',nc,1);
    rij = sqrt(xdiff.^2 + ydiff.^2);
    rij = tril(rij);
    rij = sparse(rij);

    [T,~] = graphminspantree(rij);
    %[I,J] = find(T);

    %Try wrapping all positions, save the spanning tree which is
    %smallest
    T_old = T;
    xi_new = xi;
    yi_new = yi;
    [~,I] = sort(xi);
    %Wrapping in x...
    if ( (sum(xi<-Lx/4)>0) && (sum(xi>Lx/4)>0) )
      for ii =  1:length(xi)
        test_xi = xi;
        test_xi(I(1:ii)) = test_xi(I(1:ii)) + Lx;
        xdiff= repmat(test_xi,1,nc) - repmat(test_xi',nc,1);
        ydiff= repmat(yi,1,nc) - repmat(yi',nc,1);
        rij = sqrt(xdiff.^2 + ydiff.^2);
        rij = tril(rij);
        rij = sparse(rij);
        [T,~] = graphminspantree(rij);
        if sum(T(:)) < sum(T_old(:))
          disp('Found smaller spanning tree by wrapping in x')
          T_old = T;
          xi_new = test_xi;
        end
      end
    end

     %Wrapping in y
     if ((sum(yi<-Ly/4)>0) && (sum(yi>Ly/4)>0))
       [~,I] = sort(yi);
       for ii =  1:length(yi)
         test_yi = yi;
         test_yi(I(1:ii)) = test_yi(I(1:ii)) + Ly;
         ydiff= repmat(test_yi,1,nc) - repmat(test_yi',nc,1);
         xdiff= repmat(xi_new,1,nc) - repmat(xi_new',nc,1);
         rij = sqrt(xdiff.^2 + ydiff.^2);
         rij = tril(rij);
         rij = sparse(rij);
         [T,~] = graphminspantree(rij);
         if sum(T(:)) < sum(T_old(:))
           disp('Found smaller spanning tree by wrapping in y')
           T_old = T;
           yi_new = test_yi;
         end
       end
     end

%Get the radius and centre of mass position
XC = mean(xi_new);
YC = mean(yi_new);
Xrad = mean(abs(xi_new-XC));
Yrad = mean(abs(yi_new-YC));
Radius = sqrt(Xrad.^2 + Yrad.^2); %#ok
XC = XC -Lx*(XC > Lx/2);
YC = YC -Ly*(YC > Ly/2);
XC = XC +Lx*(XC < -Lx/2);
YC = YC +Ly*(YC < -Ly/2);
RC = [XC YC]; %#ok

    clear negindex
    [I,J] = find(T_old); %#ok
   eval(['PCLUSTERS.cluster' num2str(zz) '.spanningtreepointers = [I J];']);
   eval(['PCLUSTERS.cluster' num2str(zz) '.RC = RC;']);
   eval(['PCLUSTERS.cluster' num2str(zz) '.Radius = Radius;']);
end


text:
"""
    struct Density{T,S} <: AbstractVector{T}

Discrete density with values `vals` on the quadrature nodes of `mesh::S`.
"""
struct Density{V,S<:NystromMesh} <: AbstractVector{V}
    vals::Vector{V}
    mesh::S
end

vals(σ::Density) = σ.vals
mesh(σ::Density) = σ.mesh

Base.size(σ::Density,args...)     = size(σ.vals,args...)
Base.getindex(σ::Density,args...) = getindex(σ.vals,args...)
Base.setindex(σ::Density,args...) = setindex(σ.vals,args...)

Density(etype::DataType,surf)  = Density(zeros(etype,length(surf)),surf)
Density(pde::AbstractPDE,surf) = Density(default_density_eltype(pde),surf)

function Density(f::Function,X)
    vals = [f(dof) for dof in dofs(X)]
    return Density(vals,X)
end

function γ₀(f,X)
    Density(x->f(coords(x)),X)
end

function γ₁(f,X)
    Density(dof->f(coords(dof),normal(dof)),X)
end

Base.zero(σ::Density) = Density(zero(σ.vals),mesh(σ))

# overload some unary/binary operations for convenience
Base.:-(σ::Density) = Density(-σ.vals,σ.mesh)
Base.:+(σ::Density) = σ
Base.:*(a::Number,σ::Density) = Density(a*σ.vals,σ.mesh)
Base.:*(σ::Density,a::Number) = a*σ
Base.:/(σ::Density,a::Number) = Density(σ.vals/a,σ.mesh)
Base.:*(A::Matrix,σ::Density) = Density(A*σ.vals,σ.mesh)
Base.:*(A::Diagonal,σ::Density) = Density(A*σ.vals,σ.mesh)

function Base.:\(A::AbstractMatrix{T},σ::Density{V}) where {T,V}
    @assert size(A,1) == size(A,2)
    if T <: Number && V <: Number
        y = A\σ.vals
        return Density(y,σ.mesh)
    elseif T <: SMatrix && V <: SVector
        Amat    = blockmatrix_to_matrix(A)
        σ_vec   = reinterpret(eltype(V),σ.vals)
        vals_vec = Amat\σ_vec
        vals    = reinterpret(V,vals_vec) |> collect
        return Density(vals,σ.mesh)
    else
        notimplemented()
    end
end

function IterativeSolvers.gmres!(σ::Density{V},A::AbstractMatrix{T},μ::Density{V},args...;kwargs...) where {T,V}
    log = haskey(kwargs,:log) ? kwargs[:log] : false
    if T <: Number && V <: Number
        if log
            vals,hist = gmres!(σ.vals,A,μ.vals,args...;kwargs...)
            return σ,hist
        else
            vals = gmres!(σ.vals,A,μ.vals,args...;kwargs...)
            return σ
        end
    elseif T <: SMatrix && V <: SVector
        Amat    = blockmatrix_to_matrix(A)
        σ_vec   = reinterpret(eltype(V),σ.vals)
        μ_vec   = reinterpret(eltype(V),μ.vals)
        if log
            vals,hist = gmres!(σ_vec,Amat,μ_vec,args...;kwargs...)
            return σ,hist
        else
            vals = gmres!(σ_vec,Amat,μ_vec,args...;kwargs...)
            return σ
        end
    else
        notimplemented()
    end
end
IterativeSolvers.gmres(A,μ::Density,args...;kwargs...) = gmres!(zero(μ),A,μ,args...;kwargs...)

function ncross(σ::Density)
    Γ = mesh(σ)
    iter = zip(vals(σ),dofs(Γ))
    v = map(iter) do (v,dof)
        cross(normal(dof),v)
    end
    Density(v,Γ)
end

function istangential(σ::Density,tol=1e-8)
    Γ = mesh(σ)
    iter = zip(vals(σ),dofs(Γ))
    all(iter) do (v,dof)
        abs(dot(v,normal(dof))) < tol
    end
end

"""
    struct TangentialDensity{T,S} <: AbstractVector{T}

A density tangential to the surface defined by `mesh`. The `vals` field stores
the components on the tangential basis defined by the surface `jacobian` at each
`mesh` point.

!!! note
    Calling `TangentialDensity` on a `Density` object `σ` will compute `σ - (σ⋅n)n`,
    and express it using the the basis induced by the jacobian at each point.
    Unless `σ` is already a tangetial field, this is a projection and not only
    a change of basis.
"""
struct TangentialDensity{V,S<:NystromMesh} <: AbstractVector{V}
    vals::Vector{V}
    mesh::S
end

vals(σ::TangentialDensity) = σ.vals
mesh(σ::TangentialDensity) = σ.mesh

Base.size(σ::TangentialDensity,args...)     = size(σ.vals,args...)
Base.getindex(σ::TangentialDensity,args...) = getindex(σ.vals,args...)
Base.setindex(σ::TangentialDensity,args...) = setindex(σ.vals,args...)

function TangentialDensity(σ::Density)
    # original density type must be a 3D vector
    @assert eltype(vals(σ)) <: SVector{3}
    Γ = mesh(σ)
    @assert typeof(Γ|>dofs|>first|>jacobian) <: SMatrix{3,2}
    iter = zip(vals(σ),dofs(Γ))
    v = map(iter) do (σ,dof)
        jac = jacobian(dof)
        rhs = transpose(jac)*σ
        A   = transpose(jac)*jac
        A\rhs
    end
    TangentialDensity(v,Γ)
end

function ncross(σ::TangentialDensity)
    # original density type must be a 2D vector
    @assert eltype(vals(σ)) <: SVector{2}
    Γ = mesh(σ)
    @assert typeof(Γ|>dofs|>first|>jacobian) <: SMatrix{3,2}
    iter = zip(vals(σ),dofs(Γ))
    vlist = map(iter) do (v,dof)
        jac = jacobian(dof)
        t1 = jac[:,1]
        t2 = jac[:,2]
        # metric tensor coefficients
        E = dot(t1,t1)
        G = dot(t2,t2)
        F = dot(t1,t2)
        dS = sqrt(E*G - F^2)  # differential area
        SVector(-F*v[1] - G*v[2], E*v[1] + F*v[2]) / dS
    end
    return TangentialDensity(vlist,Γ)
end

function Density(σ::TangentialDensity)
    Γ = mesh(σ)
    iter = zip(vals(σ),dofs(Γ))
    v = map(iter) do (σ,dof)
        jac = jacobian(dof)
        jac*σ
    end
    Density(v,Γ)
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using FIGlet
using Test

function generate_output(s::AbstractString, font::AbstractString="Standard")
    fontpath = FIGlet.getfontpath(font)
    io = IOContext(IOBuffer())
    FIGlet.render(io, s, FIGlet.readfont(fontpath))
    jl_output = String(take!(io.io))
    cli_output = read(`figlet -f $fontpath $s`, String)

    return strip(join(strip.(split(jl_output, '\n')), '\n')), strip(join(strip.(split(cli_output, '\n')), '\n'))
end

@testset "FIGlet.jl" begin
    iob = IOBuffer(b"flf2a", read=true);
    @test FIGlet.readmagic(iob) == UInt8['f', 'l', 'f', '2', 'a']

    @test FIGlet.FIGletHeader('$', 6, 5, 16, 15, 11) == FIGlet.FIGletHeader('$', 6, 5, 16, 15, 11, 0, 143, 0)

    @test FIGlet.availablefonts() |> length == 680

    for font in FIGlet.availablefonts()
        @test FIGlet.readfont(font) isa FIGlet.FIGletFont
    end

    ff = FIGlet.readfont("jiskan16")
    iob = IOBuffer()
    print(iob, ff)
    @test String(take!(iob)) == "FIGletFont(n=7098)"

    iob = IOBuffer()
    print(iob, ff.font_characters['㙤'])
    @test String(take!(iob)) == "FIGletChar(ord='㙤')"

    @test_throws FIGlet.FontNotFoundError FIGlet.readfont("wat")
    @test_throws FIGlet.FontError FIGlet.readfont(joinpath(@__DIR__, "..", "README.md"))

    iob = IOBuffer()
    FIGlet.render(iob, "the quick brown fox jumps over the lazy dog", "standard")
    @test String(take!(iob)) == raw"""
 _   _                        _      _      _
| |_| |__   ___    __ _ _   _(_) ___| | __ | |__  _ __ _____      ___ __
| __| '_ \ / _ \  / _` | | | | |/ __| |/ / | '_ \| '__/ _ \ \ /\ / / '_ \
| |_| | | |  __/ | (_| | |_| | | (__|   <  | |_) | | | (_) \ V  V /| | | |
 \__|_| |_|\___|  \__, |\__,_|_|\___|_|\_\ |_.__/|_|  \___/ \_/\_/ |_| |_|
                     |_|

  __               _
 / _| _____  __   (_)_   _ _ __ ___  _ __  ___    _____   _____ _ __
| |_ / _ \ \/ /   | | | | | '_ ` _ \| '_ \/ __|  / _ \ \ / / _ \ '__|
|  _| (_) >  <    | | |_| | | | | | | |_) \__ \ | (_) \ V /  __/ |
|_|  \___/_/\_\  _/ |\__,_|_| |_| |_| .__/|___/  \___/ \_/ \___|_|
                |__/                |_|

 _   _            _                       _
| |_| |__   ___  | | __ _ _____   _    __| | ___   __ _
| __| '_ \ / _ \ | |/ _` |_  / | | |  / _` |/ _ \ / _` |
| |_| | | |  __/ | | (_| |/ /| |_| | | (_| | (_) | (_| |
 \__|_| |_|\___| |_|\__,_/___|\__, |  \__,_|\___/ \__, |
                              |___/               |___/

"""

    iob = IOBuffer()
    FIGlet.render(iob, uppercase("the quick brown fox jumps over the lazy dog"), "standard")
    @test String(take!(iob)) == raw"""
 _____ _   _ _____    ___  _   _ ___ ____ _  __
|_   _| | | | ____|  / _ \| | | |_ _/ ___| |/ /
  | | | |_| |  _|   | | | | | | || | |   | ' /
  | | |  _  | |___  | |_| | |_| || | |___| . \
  |_| |_| |_|_____|  \__\_\\___/|___\____|_|\_\

 ____  ____   _____        ___   _   _____ _____  __
| __ )|  _ \ / _ \ \      / / \ | | |  ___/ _ \ \/ /
|  _ \| |_) | | | \ \ /\ / /|  \| | | |_ | | | \  /
| |_) |  _ <| |_| |\ V  V / | |\  | |  _|| |_| /  \
|____/|_| \_\\___/  \_/\_/  |_| \_| |_|   \___/_/\_\

     _ _   _ __  __ ____  ____     _____     _______ ____
    | | | | |  \/  |  _ \/ ___|   / _ \ \   / / ____|  _ \
 _  | | | | | |\/| | |_) \___ \  | | | \ \ / /|  _| | |_) |
| |_| | |_| | |  | |  __/ ___) | | |_| |\ V / | |___|  _ <
 \___/ \___/|_|  |_|_|   |____/   \___/  \_/  |_____|_| \_\

 _____ _   _ _____   _        _     _______   __  ____   ___   ____
|_   _| | | | ____| | |      / \   |__  /\ \ / / |  _ \ / _ \ / ___|
  | | | |_| |  _|   | |     / _ \    / /  \ V /  | | | | | | | |  _
  | | |  _  | |___  | |___ / ___ \  / /_   | |   | |_| | |_| | |_| |
  |_| |_| |_|_____| |_____/_/   \_\/____|  |_|   |____/ \___/ \____|

"""

    iob = IOBuffer()
    FIGlet.render(iob, "the quick brown fox jumps over the lazy dog", "mirror")
    @test String(take!(iob)) == raw"""
                              _      _      _                        _   _
   __ ___      _____ __ _  __| | __ | |___ (_)_   _ _ __    ___   __| |_| |
  / _` \ \ /\ / / _ \__` |/ _` | \ \| |__ \| | | | | '_ \  / _ \ / _` |__ |
 | | | |\ V  V / (_) | | | (_| |  >   |__) | | |_| | |_) | \__  | | | |_| |
 |_| |_| \_/\_/ \___/  |_|\__,_| /_/|_|___/|_|_.__/| .__/  |___/|_| |_|__/
                                                   |_|

                                                  _               __
  __ _ _____   _____    ___  __ _  ___ __ _ _   _(_)   __  _____ |_ \
 |__` / _ \ \ / / _ \  |__ \/ _` |/ _ ' _` | | | | |   \ \/ / _ \ _| |
    | \__  \ V / (_) | / __/ (_| | | | | | | |_| | |    >  < (_) |_  |
    |_|___/ \_/ \___/  \___|\__, |_| |_| |_|_.__/| \_  /_/\_\___/  |_|
                               |_|                \__|

              _                       _            _   _
  _ __   ___ | |__    _   _____ _ __ | |  ___   __| |_| |
 | '_ \ / _ \| '_ \  | | | \  _| '_ \| | / _ \ / _` |__ |
 | |_) | (_) | |_) | | |_| |\ \| |_) | | \__  | | | |_| |
 | .__/ \___/|_.__/  | .__/|___\_.__/|_| |___/|_| |_|__/
  \___|               \___|

"""

    iob = IOBuffer()
    FIGlet.render(iob, uppercase("the quick brown fox jumps over the lazy dog"), "mirror")
    @test String(take!(iob)) == raw"""
 __  _ ____ ___ _   _  ___    _____ _   _ _____
 \ \| |___ \_ _| | | |/ _ \  |____ | | | |_   _|
  \ ` |   | | || | | | | | |   |_  | |_| | | |
  / . |___| | || |_| | |_| |  ___| |  _  | | |
 /_/|_|____/___|\___//_/__/  |_____|_| |_| |_|

 __  _____ _____   _   ___        _____   ____  ____
 \ \/ / _ \___  | | | / \ \      / / _ \ / _  |( __ |
  \  / | | | _| | | |/  |\ \ /\ / / | | | (_| |/ _  |
  /  \ |_| ||_  | |  /| | \ V  V /| |_| |> _  | (_| |
 /_/\_\___/   |_| |_/ |_|  \_/\_/  \___//_/ |_|\____|

   ____ _______     _____     ____  ____ __  __ _   _ _
  / _  |____ \ \   / / _ \   |___ \/ _  |  \/  | | | | |
 | (_| | |_  |\ \ / / | | |  / ___/ (_| | |\/| | | | | |  _
  > _  |___| | \ V /| |_| | | (___ \__  | |  | | |_| | |_| |
 /_/ |_|_____|  \_/  \___/   \____|   |_|_|  |_|\___/ \___/

  ____   ___   ____  __   _______     _        _   _____ _   _ _____
 |___ \ / _ \ / _  | \ \ / /\  __|   / \      | | |____ | | | |_   _|
  _  | | | | | | | |  \ V /  \ \    / _ \     | |   |_  | |_| | | |
 | |_| | |_| | |_| |   | |   _\ \  / ___ \ ___| |  ___| |  _  | | |
 |____/ \___/ \____|   |_|  |____\/_/   \_\_____| |_____|_| |_| |_|

"""

    @testset "Render all fonts" begin

        for (i, font) in enumerate(FIGlet.availablefonts())
            try
                iob = IOBuffer()
                sentence = "the quick brown fox jumps over the lazy dog"
                FIGlet.render(iob, sentence, font)
                FIGlet.render(iob, uppercase(sentence), font)
                @test length(String(take!(iob))) > 0
            catch
                println("Cannot render font: number = $i, name = \"$font\"")
            end
        end

    end

end


text:
#=
nufft.jl
Non-uniform FFT (NUFFT), currently a wrapper around NFFT.jl
todo: open issues: small N, odd N, nufft!, adjoint!
2019-06-06, Jeff Fessler, University of Michigan
=#

export Anufft, nufft_init

#using MIRT: dtft_init, map_many
#include("../utility/map_many.jl")
using NFFT
using LinearAlgebra: norm
using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO

"""
    nufft_eltype(::DataType)
ensure plan_nfft eltype is Float32 or Float64
"""
nufft_eltype(::Type{<:Integer}) = Float32
nufft_eltype(::Type{<: Union{Float16,Float32}}) = Float32
nufft_eltype(::Type{Float64}) = Float64
nufft_eltype(T::DataType) = throw("unknown type $T")

# the following convenience routine ensures correct type passed to nfft()
# see https://github.com/tknopp/NFFT.jl/pull/33
# todo: may be unnecessary with future version of nfft()
"""
    nufft_typer(T::DataType, x::AbstractArray{<:Real} ; warn::Bool=true)
type conversion wrapper for `nfft()`
"""
nufft_typer(::Type{T}, x::T ; warn::Bool=true) where {T} = x # cf convert()

function nufft_typer(T::Type{TT}, x ; warn::Bool=true) where {TT}
    isinteractive() && @warn("converting $(eltype(x)) to $(eltype(T))")
    convert(T, x)
end

"""
    p = nufft_init(w, N ; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=0)

Setup 1D NUFFT,
for computing fast ``O(N \\log N)`` approximation to

``X[m] = \\sum_{n=0}^{N-1} x[n] \\exp(-i w[m] (n - n_{shift})), m=1,…,M``

in
- `w::AbstractArray{<:Real}` `[M]` frequency locations (units radians/sample)
   + `eltype(w)` determines the `plan_nfft` type; so to save memory use Float32!
   + `size(w)` determines `odim` for `A` if `operator=true`
- `N::Int` signal length

option
- `nfft_m::Int`		see NFFT.jl documentation; default 4
- `nfft_sigma::Real`	"", default 2.0
- `n_shift::Real`		often is N/2; default 0
- `pi_error::Bool`		throw error if ``|w| > π``, default `true`
   + Set to `false` only if you are very sure of what you are doing!
- `do_many::Bool`	support extended inputs via `map_many`? default `true`
- `operator::Bool=true` set to `false` to make `A` an `LinearMapAM`

out
- `p NamedTuple`
`(nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A::LinearMapAO)`

The default settings are such that for a 1D signal of length N=512,
the worst-case error is below 1e-5 which is probably adequate
for typical medical imaging applications.
To verify this statement, run `nufft_plot1()` and see plot.
"""
function nufft_init(
	w::AbstractArray{<:Real},
	N::Int ;
	n_shift::Real = 0,
	nfft_m::Int = 4,
	nfft_sigma::Real = 2.0,
	pi_error::Bool = true,
	do_many::Bool = true,
	operator::Bool = true, # !
)

	N < 6 && throw("NFFT may be erroneous for small N")
	isodd(N) && throw("NFFT erroneous for odd N")
	pi_error && any(abs.(w) .> π) &&
		throw(ArgumentError("|w| > π is likely an error"))

	T = nufft_eltype(eltype(w))
	CT = Complex{T}
	CTa = AbstractArray{Complex{T}}
	f = convert(Array{T}, vec(w)/(2π)) # note: plan_nfft must have correct type
	p = plan_nfft(f, N, nfft_m, nfft_sigma) # create plan
	M = length(w)
	# extra phase here because NFFT always starts from -N/2
	phasor = convert(CTa, cis.(-vec(w) * (N/2 - n_shift)))
	phasor_conj = conj.(phasor)
	forw1 = x -> nfft(p, nufft_typer(CTa, x)) .* phasor
#	forw! = x,y -> nfft!(p, nufft_typer(CTa, x)) .* phasor # todo
	back1 = y -> nfft_adjoint(p, nufft_typer(CTa, y .* phasor_conj))

	prop = (name="nufft1", w=w, N=(N,), n_shift=n_shift,
			nfft_m=nfft_m, nfft_sigma=nfft_sigma)
	A = LinearMapAA(forw1, back1, (M, N) ; # no "many" here!
		prop = prop, T=CT, operator = operator, # effectively "many" if true
		odim = operator ? size(w) : (length(w),),
	)

	if do_many
		forw = x -> map_many(forw1, x, (N,))
		back = y -> map_many(back1, y, (M,))
	else
		forw = forw1
		back = back1
	end

	return (nufft=forw, adjoint=back, A=A)
end

"""
    p = nufft_init(w, N ; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=?)

Setup multi-dimensional NUFFT,
for computing fast ``O(N \\log N)`` approximation to

``X[m] = \\sum_{n=0}^{N-1} x[n] \\exp(-i w[m,:] (n - n_{shift})), m=1,…,M``

in
- `w::AbstractArray{<:Real}` `[M,D]` frequency locations (units radians/sample)
	+ `eltype(w)` determines the `plan_nfft` type; so to save memory use Float32!
	+ `size(w)[1:(end-1)]` determines `odim` if `operator=true`
- `N::Dims{D}` signal dimensions

option
- `nfft_m::Int`		see NFFT.jl documentation; default 4
- `nfft_sigma::Real`	"", default 2.0
- `n_shift::AbstractVector{<:Real}`	`[D]`	often is N/2; default zeros(D)
- `pi_error::Bool`		throw error if ``|w| > π``, default `true`
   + Set to `false` only if you are very sure of what you are doing!
- `do_many::Bool`	support extended inputs via `map_many`? default `true`
- `operator::Bool=true` set to `false` to make `A` an `LinearMapAM`

The default `do_many` option is designed for parallel MRI where the k-space
sampling pattern applies to every coil.
It may also be useful for dynamic MRI with repeated sampling patterns.
The coil and/or time dimensions must come after the spatial dimensions.

out
- `p NamedTuple` with fields
	`nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A=LinearMapAO`
	(Using `operator=true` allows the `LinearMapAO` to support `do_many`.)
"""
function nufft_init(
	w::AbstractArray{<:Real},
	N::Dims{D} ;
	n_shift::AbstractVector{<:Real} = zeros(Int, length(N)),
	nfft_m::Int = 4,
	nfft_sigma::Real = 2.0,
	pi_error::Bool = true,
	do_many::Bool = true,
	operator::Bool = true, # !
) where {D}

	any(N .< 6) && throw("NFFT may be erroneous for small N")
	any(isodd.(N)) && throw("NFFT erroneous for odd N")
	pi_error && any(abs.(w) .> π) &&
		throw(ArgumentError("|w| > π is likely an error"))

	ndims(w) > 1 || throw("ndims(w)==1 invalid")
	size(w)[end] != D && throw(DimensionMismatch("$(size(w)) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))
	odim = size(w)[1:(end-1)] # trick, e.g., for radial sampling
	w = reshape(w, :, D) # [M,D]
	M = size(w)[1]

	T = nufft_eltype(eltype(w))
	CT = Complex{T}
	CTa = AbstractArray{Complex{T}}
	f = convert(Array{T}, w/(2π)) # note: plan_nfft must have correct type
	p = plan_nfft(f', N, nfft_m, nfft_sigma) # create plan

	# extra phase here because NFFT.jl always starts from -N/2
	phasor = convert(CTa, cis.(-w * (collect(N)/2. - n_shift)))
	phasor_conj = conj.(phasor)
	# todo: in-place
	forw1 = x -> nfft(p, nufft_typer(CTa, x)) .* phasor
	back1 = y -> nfft_adjoint(p, nufft_typer(CTa, y .* phasor_conj))

	prop = (name="nufft$(length(N))", w=w, N=N, n_shift=n_shift,
			nfft_m=nfft_m, nfft_sigma=nfft_sigma)
	if operator # LinearMapAO
		A = LinearMapAA(forw1, back1, (M, prod(N)) ;
			prop = prop, T = CT,
			operator = true, odim = odim, idim = N,
		)
	else
		# no "many" for LinearMapAM:
		A = LinearMapAA(x -> forw1(reshape(x,N)), y -> vec(back1(y)),
			(M, prod(N)) ; prop = prop, T = CT,
		)
	end

	if do_many
		forw = x -> map_many(forw1, x, N)
		back = y -> map_many(back1, y, (M,))
	else
		forw = forw1
		back = back1
	end

	return (nufft=forw, adjoint=back, A=A)
end

"""
    Anufft(ω, N ; kwargs ...)

Make a `LinearMapAO` object of size `length(ω) × prod(N)`.
See `nufft_init` for options.
"""
Anufft(w::AbstractArray{<:Real}, N::Int ; kwargs...) =
	nufft_init(w, N ; kwargs...).A
Anufft(w::AbstractArray{<:Real}, N::Dims ; kwargs...) =
	nufft_init(w, N ; kwargs...).A

"""
w, errs = nufft_errors( ; M=?, w=?, N=?, n_shift=?, ...)

Compute worst-case errors for NUFFT (for signal of length N of unit norm)
"""
function nufft_errors( ;
	M::Int = 401,
	N::Int = 512,
	w::AbstractArray{<:Real} = LinRange(0, 2π/N, M),
	n_shift::Real = 0,
	kwargs...,
)

	sd = dtft_init(w, N ; n_shift=n_shift)
	sn = nufft_init(w, N ; n_shift=n_shift, kwargs...)
	E = Matrix(sn.A - sd.A)
	return w, vec(mapslices(norm, E, dims=2)) # [M]
end


text:
include("utils.jl")
include("wrapper.jl")
include("client.jl")
include("producer.jl")
include("consumer.jl")


text:
module Simple
using Hijack, Test

@testset "simple" begin
    @test true
    push!(Hijack.RUN, 4)
end

end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
using Test, Logging, Random
using DecisionProgramming

rng = MersenneTwister(4)

@info "Testing random_diagram"
diagram = InfluenceDiagram()
@test_throws DomainError random_diagram!(rng, diagram, -1, 1, 1, 1, 1, [2])
@test_throws DomainError random_diagram!(rng, diagram, 1, -1, 1, 1, 1, [2])
@test_throws DomainError random_diagram!(rng, diagram, 0, 0, 1, 1, 1, [2])
@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 0, 1, 1, [2])
@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 0, 1, [2])
@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 1, 0, [2])
@test_throws DomainError random_diagram!(rng, diagram, 1, 1, 1, 1, 1, [1])

for (n_C, n_D) in [(1, 0), (0, 1)]
    rng = RandomDevice()
    diagram = InfluenceDiagram()
    random_diagram!(rng, diagram, n_C, n_D, 1, 1, 1, [2])
    @test isa(diagram.C, Vector{Node})
    @test isa(diagram.D, Vector{Node})
    @test isa(diagram.V, Vector{Node})
    @test length(diagram.C) == n_C
    @test length(diagram.D) == n_D
    @test length(diagram.V) == 1
    @test all(!isempty(I_v) for I_v in diagram.I_j[diagram.V])
    @test isa(diagram.S, States)
end

@info "Testing random Probabilities"
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 2, 2, 1, 1, 1, [2])
@test_throws DomainError random_probabilities!(rng, diagram, diagram.D[1]; n_inactive=0)
random_probabilities!(rng, diagram, diagram.C[1]; n_inactive=0)
@test isa(diagram.X[1], Probabilities)
random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=1)
@test isa(diagram.X[2], Probabilities)

diagram = InfluenceDiagram()
diagram.C = Node[1,3]
diagram.D = Node[2]
diagram.I_j = [Node[], Node[], Node[1,2]]
diagram.S = States(State[2, 3, 2])
diagram.X = Vector{Probabilities}(undef, 2)
random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=2*3*(2-1))
@test isa(diagram.X[2], Probabilities)
@test_throws DomainError random_probabilities!(rng, diagram, diagram.C[2]; n_inactive=2*3*(2-1)+1)

@info "Testing random Utilities"
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 1, 1, 2, 1, 1, [2])
random_utilities!(rng, diagram, diagram.V[1], low=-1.0, high=1.0)
@test isa(diagram.Y[1], Utilities)
@test_throws DomainError random_utilities!(rng, diagram, diagram.V[1], low=1.1, high=1.0)
@test_throws DomainError random_utilities!(rng, diagram, diagram.C[1], low=-1.0, high=1.0)

@info "Testing random LocalDecisionStrategy"
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 2, 2, 2, 2, 2, [2])
@test isa(LocalDecisionStrategy(rng, diagram, diagram.D[1]), LocalDecisionStrategy)


text:
export add_find_max, add_find_min, add_find_ith

function add_find_max(mgr::Ptr{Manager}, f::Ptr{Node})
    result = Cudd_addFindMax(mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_find_min(mgr::Ptr{Manager}, f::Ptr{Node})
    result = Cudd_addFindMin(mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_find_ith(mgr::Ptr{Manager}, f::Ptr{Node}, i::Int)
    result = Cudd_addIthBit(mgr, f, i)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end