text:
function alphaj(a::Complex128)
    return sqrt(1 - abs(a)^2)
end

function dictionaryElement(z::Complex128,a::Complex128)
    return alphaj(a)/(z-a)
end

function blaschkeElement(z::Complex128,a::Complex128)
    return (1-conj(a)*z)/(z-a)
end

function invBlaschkeElement(z::Complex128,a::Complex128)
    return (z-a)/(1-conj(a)*z)
end

function blaschkeProd(z::Complex128,aj::Vector{Complex128})
    r = one(z)
    for i=1:length(aj)
        r *= blaschkeElement(z,aj[i])
    end
    return r
end

function blaschkeProdExceptQ(z::Complex128,aj::Vector{Complex128},q::Int)
    r = one(z)
    for i=1:length(aj)
        if i!=q
            r *= blaschkeElement(z,aj[i])
        end
    end
    return r
end

function impulseTakenakaMalmquist(n,aj::Vector{Complex128})
    if n<=zero(n)
        return zero(Complex128)
    end
    temp1 = zero(Complex128)
    temp2 = zero(Complex128)
    for i=1:(length(aj)-1)
        temp1 += aj[i]^(n-1)*dictionaryElement(aj[i],aj[end])*blaschkeProdExceptQ(aj[i],aj[1:(length(aj)-1)],i)*(1-abs(aj[i])^2)
    end
    temp2 += aj[end]^(n-1)*alphaj(aj[end])*blaschkeProd(aj[end],aj[1:(length(aj)-1)])
    return temp1 + temp2
end


text:
"""
    hit_and_run(
        N::Int,
        opt_model;
        keepevery = 100,
        samplesize = 1000,
        random_objective = false,
    )

Perform a basic hit and run sampling for `N` iterations on a constrained JuMP
model in `opt_model`. See "Robert L. Smith Efficient Monte Carlo Procedures for
Generating Points Uniformly Distributed over Bounded Regions. Operations
Research 32 (6) 1296-1308 https://doi.org/10.1287/opre.32.6.1296" for more
details.

The process generates `samplesize` samples, and logs the sample state each
`keepevery` iterations.

Warm up points are generated by minimizing and maximizing reactions as in
[`flux_variability_analysis`](@ref), unless the `random_objective` is `true`,
in which case a randomly weighted objective is used for warmup.

Note that `N` needs to be greater than sample size, and should be greater than
the dimensionality of the sampled space (i.e., at least same as the number of
reactions).

# Example
```
using COBREXA
using JuMP
using Tulip

model = load_model(StandardModel, "e_coli_core.json")
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
glucose = findfirst(model.reactions, "EX_glc__D_e")

opt_model = flux_balance_analysis(model, Tulip.Optimizer;
    modifications=[change_objective(biomass),
    modify_constraint(glucose, -12, -12),
    change_optimizer_attribute("IPM_IterationsLimit", 500)])

biomass_index = model[biomass]
λ = JuMP.value(opt_model[:x][biomass_index])
modify_constraint(biomass, 0.99*λ, λ)(model, opt_model)

samples = hit_and_run(100_000, opt_model; keepevery=10, samplesize=5000)
```
"""
function hit_and_run(
    N::Int,
    opt_model;
    keepevery = _constants.sampling_keep_iters,
    samplesize = _constants.sampling_size,
    random_objective = false,
)

    lbs, ubs = get_bound_vectors(opt_model) # get actual ub and lb constraints, can't use model function because the user may have changed them in the function arguments

    wpoints = _get_warmup_points(opt_model; random_objective = random_objective)

    nwpts = size(wpoints, 2) # number of warmup points generated
    samples = zeros(size(wpoints, 1), samplesize) # sample storage
    current_point = zeros(size(wpoints, 1))
    current_point .= wpoints[:, rand(1:nwpts)] # pick random initial point

    sample_num = 0
    samplelength = 0
    updatesamplesizelength = true
    for n = 1:N

        # direction = random point - current point
        if updatesamplesizelength
            direction_point = (@view wpoints[:, rand(1:nwpts)]) - (@view current_point[:]) # use warmup points to find direction in warmup phase
        else
            direction_point =
                (@view samples[:, rand(1:(samplelength))]) - (@view current_point[:]) # after warmup phase, only find directions in sampled space
        end

        λmax = Inf
        λmin = -Inf
        for i in eachindex(lbs)
            δlower = lbs[i] - current_point[i]
            δupper = ubs[i] - current_point[i]
            # only consider the step size bound if the direction of travel is non-negligible
            if direction_point[i] < -_constants.tolerance
                lower = δupper / direction_point[i]
                upper = δlower / direction_point[i]
            elseif direction_point[i] > _constants.tolerance
                lower = δlower / direction_point[i]
                upper = δupper / direction_point[i]
            else
                lower = -Inf
                upper = Inf
            end
            lower > λmin && (λmin = lower) # max min step size that satisfies all bounds
            upper < λmax && (λmax = upper) # min max step size that satisfies all bounds
        end

        if λmax <= λmin || λmin == -Inf || λmax == Inf # this sometimes can happen
            @warn "Infeasible direction at iteration $(n)..."
            continue
        end

        λ = rand() * (λmax - λmin) + λmin # random step size
        current_point .= current_point .+ λ .* direction_point # will be feasible

        if n % keepevery == 0
            sample_num += 1
            samples[:, sample_num] .= current_point
            if sample_num >= samplesize
                updatesamplesizelength = false # once the entire memory vector filled, stop using warm up points
                sample_num = 0 # reset, start replacing the older samples
            end
            updatesamplesizelength && (samplelength += 1) # lags sample_num because the latter is a flag as well
        end

    end

    return samples
end

ext: jl


text:
@generated apply(a,b) = :( apply(a,b,Val( $( hasmethod(*,Tuple{a,b}) ) )) )
apply(a,b,::Val{true}) = a*b
apply(a,b,::Val{false}) = a(b)

"""
    tucker(C,B)

Apply `B[k]` to the `k`th dimension of tensor `C`.

# Examples
```
julia> C = rand(2); B = (rand(3,2),);
       tucker(C,B) == B[1]*C
true

julia> C = rand(2,3); B = (rand(4,2),rand(5,3));
       tucker(C,B) == B[1]*C*B[2]'
true
```
"""
function tucker end

@generated function tucker(
    C::AbstractArray{<:Any,N},
    B::NTuple{N,Any}
) where {N}
    quote
        C_0 = C
        Base.Cartesian.@nexprs $N k->begin
            tmp_k = reshape(C_{k-1},(size(C_{k-1},1),prod(Base.tail(size(C_{k-1})))))
            tmp_k = apply(B[k],tmp_k)
            C_k = Array(reshape(transpose(tmp_k),(Base.tail(size(C_{k-1}))...,size(tmp_k,1))))
            # https://github.com/JuliaLang/julia/issues/30988
        end
        return $(Symbol("C_",N))
    end
end
tucker(C::AbstractArray{<:Any,1}, B::NTuple{1,Any}) = apply(B[1],C)
tucker(C::AbstractArray{<:Any,2}, B::NTuple{2,Any}) = transpose(apply(B[2],transpose(apply(B[1],C))))


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the acoustic perturbation equations

equations = AcousticPerturbationEquations2D(v_mean_global=(0.0, 0.0), c_mean_global=0.0,
                                            rho_mean_global=0.0)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-20.6, 0.0) # minimum coordinates (min(x), min(y))
coordinates_max = (30.6, 51.2) # maximum coordinates (max(x), max(y))

initial_condition = initial_condition_monopole # does not use the global mean values given above
boundary_conditions = (x_neg=boundary_condition_zero,
                       x_pos=boundary_condition_zero,
                       y_neg=boundary_condition_monopole,
                       y_pos=boundary_condition_zero)

# Create a uniformely refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6,
                n_cells_max=100_000,
                periodicity=false)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 24.0
tspan =  (0.0, 24.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=0.8)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks, maxiters=50)

# Print the timer summary
summary_callback()


text:
#Newtonian eqn for impurity
function NewtonImp(phi, V, x, xdot, t, dt, M, lambda, xstep)
    for i = 1:length(x)
        xdot[i] = xdot[i] + real(dt / M * lambda * sum(V[i+1] .* conj(phi) .* phi) * xstep)
        x[i] = x[i] + real(dt * xdot[i])
    end
end

#if abspath(PROGRAM_FILE) == @__FILE__
#    include("helper.jl")
#    using BenchmarkTools
#    phi = ones(100,100)+0im*ones(100,100)
#    V = [ones(100,100)+0im*ones(100,100),ones(100,100)+0im*ones(100,100),
#    ones(100,100)+0im*ones(100,100)]
#    xarray = [i for i in range(-5 / 2.0, 5 / 2.0, length = 100)]
#    V = V_fct(V, [0,0], 0, xarray, 1)
#    x = [0.0 ,0.0]
#    xdot = [0.0, 0.0]
#    @btime NewtonImp(phi, V, x, xdot, 0.0, 0.1, 1.0, 1.0, 1.0)
#end

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
function reverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function goTo(side)
    while !isborder(robot, side)
        putmarker!(robot)
        move!(robot, side)
    end
end

function set(side1, side2, dir)
    while (!isborder(robot, side1) || !isborder(robot, side2))
        while !isborder(robot, side1)
            move!(robot, side1)
            push!(dir, side1)
        end
        while !isborder(robot, side2)
            move!(robot, side2)
            push!(dir, side2)
        end
    end
end

function f()
    dir = []
    set(Nord, West, dir)
    for i in [Sud, Ost, Nord, West]
        goTo(i)
    end
    while length(dir) > 0
        move!(robot, reverse(pop!(dir)))
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
function inflation_var_setup(d)
    f = [ symbols("ϕ$i",real=true) for i in 1:d ]
    fp= [ symbols("π$i",real=true) for i in 1:d ]
    return (f,fp)
end

function sum_separableQ(d,gv,f)
    for i in 1:d
        for j in 1:d
            if i != j && occursin(string(convert(Expr,f[i])),string(convert(Expr,gv[j])))
                return false
            end
        end
    end
    return true
end
function make_funcs(quotes)
    funcs = Dict{String,Function}()
    for (k,v) in quotes
        funcs[k] = @RuntimeGeneratedFunction(v)
    end
    return funcs
end

function inflation_setup(d,V,G,params=[];options::SetupOptions = SetupOptions())
    apply_hca = false
    f,fp = inflation_var_setup(d)
    v = V(f,params)
    g = G(f,params)
    if options.allow_caching
        hashed = hash((d,Quote(v),Quote.(g),Quote.(params),options.background_only,cache_version))
        fname = "funcs$hashed.bin"
        println("Looking for cached symbolic functions in $fname")
        if isfile(fname)
            println("Loading cached symbolic ϵV, Mab, etc...")
            open(fname,"r") do f
                quotes = deserialize(f)
            end
            if options.return_quotes
                return quotes
            else
                funcs = make_funcs(quotes)
                return funcs
            end
        end
    end
    println("Solving for symbolic ϵV, Mab, etc and compiling...")
    if d == 1
        ginv = [g[1,1]^(-1)]
    else
        ginv = options.Ginv(G,f,params)
    end
    if options.simplification > 1
        println("simplifying g, ginv")
        g = SymPy.simplify.(g)
        ginv = SymPy.simplify.(ginv)
    end
    gv = SymPy.Sym[diff(v,f[i]) for i in 1:d ]
    if options.attempt_HCA
        apply_hca = sum_separableQ(d,gv,f)
    end
    if apply_hca
        println("Detected sum-separable potential, applying HCA")
    end
    ggv = [diff(v, f[i], f[j]) for i in 1:d, j in 1:d]
    if options.simplification >= 2
        println("simplifying ggv")
        ggv = simp.(ggv)
    end
    ev = epsv(gv,ginv,v)
    eh = epsh(fp,g)
    hub = hubble(v,eh)
    et = eta(fp,eh,hub,gv)
    om = omega(eh,ev,et)
    om3 = om
    if d > 1
        om3 = omega3(fp,gv,g,ginv,hub,eh)
    end
    println("starting h")
    #h,hl = options.christoffel(f,g,ginv)
    h,hl = options.christoffel(f,g,ginv)
    println("starting pi_eom")
    pi_eom = make_pi_eom(fp,h,ginv,gv,hub,eh;simplification=options.simplification)
    println(sympy.count_ops.(pi_eom))
    cse_pi_eom = sympy.cse(pi_eom)
    println(sympy.count_ops.(cse_pi_eom))
    if !options.background_only
        println("starting Rpp")
        Rpp = Rππ_fast(f,fp,h,hl,g,ginv,simplification=options.simplification)
        println(sympy.count_ops.(Rpp))
        #println(sympy.count_ops.(simplify(pi_eom[1])))
        geo_eom = make_geodesic_eom(fp,h)
        println("starting mass matrix")
        #pi_jac = make_pi_jac(f,fp,hub,eh,pi_eom)
        Covggv = covggv(ggv,h,gv,options.simp,options.simplification)
        println(sympy.count_ops.(Covggv))
        Mab = mab(Covggv,Rpp,fp,hub,gv,eh,g,ginv,options.simp,options.simplification)
        println(sympy.count_ops.(Mab))
        cse_mab = sympy.cse(Mab)
        println("cse Mab")
        println(sympy.count_ops.(cse_mab))

    #Sig = sig(fp,g)
    #S = s(ginv,g,Sig)
    #om2,S2 = omega2(fp,g,hub,h,pi_eom,Sig)
    #Mσσ,Mss = masses(g,Mab,Sig,S)
        println("done")
        println("making pert eom")

        loga = symbols("loga",real=true)
        logk = symbols("logk",real=true)
        t = symbols("t",positive=true)
        Gam = Array{SymPy.Sym}(undef,2*d,2*d)
        for i in 1:2*d
            for j in 1:2*d
                Gam[i,j] = symbols("g$(i)o$j",real=true)
            end
        end
        GT = Array{SymPy.Sym}(undef,2,2)
        for i in 1:2
            for j in 1:i
                GT[i,j] = symbols("gt$(i)o$j",real=true)
                GT[j,i] = GT[i,j]
            end
        end

        #pert_eom = make_pert_eom(t,Gam,fp,Mab,h,eh,loga,logk,hub)
        pert_uab = make_pert_uab(t,fp,Mab,h,eh,loga,logk,hub)
        pert_curv = make_pert_curv(fp,h)
        #pert_jac = make_pert_jac(Gam,fp,pert_eom)
        tensor_eom = make_tensor_eom(t,GT,eh,loga,logk,hub)

        #println(sympy.count_ops.(pert_eom))
        #println("cse pert_eom")
        #cse_pert_eom = sympy.cse(pert_eom)
        #println(sympy.count_ops.(cse_pert_eom[2][end]))
        #println("done")
        println("pert uab")
        println(sympy.count_ops.(pert_uab))
        println("cse pert_uab")
        cse_pert_uab = sympy.cse(pert_uab)
        println(sympy.count_ops.(cse_pert_uab[2][end]))
        println("done")
        println("cse pert_curv")
        cse_pert_curv = sympy.cse(pert_curv)
        println(sympy.count_ops.(cse_pert_curv[2][end]))
        println("done")
        println("tensor eom ops:")
        println(sympy.count_ops.(tensor_eom))

        println("done making pert eom")
        #tensor_jac = make_tensor_jac(GT,tensor_eom)

        #fd = hub*fp
        #if d ==1
        #    sig = sqrt((fd*g*fd)[1,1])
        #else
        #    sig = sqrt((fd'*g*fd)[1,1])
        #end
        #t = fd / sig
        #w = Array{Expr}(undef,d)
        #for a in 1:d
        #    w[a] = 0
        #    for b in 1:d
        #        w[a] += 1/sig*(ginv[a,b]-t[a]*t[b])*gv[b]
        #    end
        #end
        #if d == 1
        #    s = w / sqrt((w*g*w)[1,1])
        #else
        #    s = w / sqrt((w'*g*w)[1,1])
        #end
        #
    end
    if apply_hca
        va = [ integrate(gv[i],f[i]) for i in 1:d]
        va .+= (v - sum(va))/d
        ea = [gvi^2 / v^2 / 2 for gvi in gv ]
        etaa = [ggvi / v for ggvi in diag(ggv)]
        Pz_hca = pz_hca(v,ea,va)
        Ns_hca = ns_hca(v,ea,etaa,eh,va)
        println("Analytic Pζ (HCA):\n $Pz_hca")
        println("Analytic ns (HCA):\n $Ns_hca")
    end
    println("quoting")
    vars = [f;fp;params]
    quotes = Dict{String,Any}()
    quotes["Eh"] = QuoteFn("_Eh",eh,vars)
    println("pi_eom")
    if d > 1
        quotes["Pi_eom"] = QuoteFnArrCSE("_Pieom",cse_pi_eom,vars)
    else
        quotes["Pi_eom"] = QuoteFnArr("_Pieom",pi_eom,vars)
    end
    quotes["V"] = QuoteFn("_Pot",v,vars)
    println("ev")
    quotes["Ev"] = QuoteFn("_Ev",ev,vars)
    quotes["Et"] = QuoteFn("_Et",et,vars)
    println("om")
    quotes["Om"] = QuoteFn("_Om",om,vars)
    #sympy.count_ops(om3)
    cse_om3 = sympy.cse(om3)
    #sympy.count_ops.(cse_om3)
    quotes["OmAlt"] = QuoteFnCSE("_Om3",cse_om3,vars)
    quotes["H"] = QuoteFn("_H",hub,vars)
    println("g")
    quotes["ginv"] = QuoteFnArr("_ginv",SymPy.Sym[ginv...],vars)
    quotes["g"] = QuoteFnArr("_g",SymPy.Sym[g...],vars)
    if !options.background_only
        println("ggv")
        quotes["Vab"] = QuoteFnArr("_ggv",ggv,vars)
        quotes["Geo_eom"] = QuoteFnArr("_Geoeom",geo_eom,vars)
        println("mab")
        quotes["Mab"] = QuoteFnArrCSE("_Mab",cse_mab,vars)
        #mσσquote = QuoteFn(Mσσ,vars)
        #mssquote = QuoteFn(Mss,vars)
        #chquote = :(function Chr($(vars...)) return [$(h...)] end)
        println("pert eom")
        #peomquote = QuoteFnArr("_Perteom",[pert_eom...],[t,loga,logk,Gam...,vars...])
        #peomquote = QuoteFnArrCSE("_Perteom",cse_pert_eom,[t,loga,logk,Gam...,vars...])
        quotes["Pert_uab"] = QuoteFnArrCSE("_Pertuab",cse_pert_uab,[t,loga,logk,vars...])
        quotes["Pert_curv"] = QuoteFnArrCSE("_Pertcurv",cse_pert_curv,[vars...])
        #pjacquote = :(function Pertjac($t,$loga,$logk,$(Gam...),$(vars...)) return [$(pert_jac...)] end )
        quotes["Tensor_eom"] = QuoteFnArr("_Tensoreom",[tensor_eom...],[t,loga,logk,unique(GT)...,vars...])
    end
    if apply_hca
        quotes["Pz_hca"] = QuoteFn("_Pzhca",Pz_hca,vars)
        quotes["ns_hca"] = QuoteFn("_Nshca",Ns_hca,vars)
    end
    #quotes = [ehquote,eomquote,vquote,ggvquote,evquote,etquote,omquote,hquote,mabquote,mσσquote,mssquote,ginvquote,gquote,peomquote,teomquote,pzhcaquote,nshcaquote]
    println("done")
    if options.allow_caching
        open(fname,"w") do f
            serialize(f,quotes)
        end
    end
    if options.return_quotes
        return quotes
    else
        lambdified = make_funcs(quotes)
        return lambdified
    end
end

function make_pi_eom(Pi,h,ginv,va,H,e;simplification=0)
    d = length(Pi)
    curv = Array{Any}(undef,d)
    vupa = Array{Any}(undef,d)
    preeom = Array{Any}(undef,d)
    eom = Array{Any}(undef,d)
    for i in 1:d
        println(i)
        curv[i] = 0
        vupa[i] = 0
        preeom[i] = 0
        for j in 1:d
            vupa[i] += ginv[i,j]*va[j]/(H^2)
            for k in 1:d
                curv[i] += h[i,j,k]*(Pi[j]*Pi[k])
            end
        end
        preeom[i] = vupa[i]+curv[i]
        if simplification >= 3
            preeom[i] = simp(vupa[i]+curv[i])
        end
    end
    for i in 1:d
        eom[i] = -preeom[i]
        eom[i] -= (3-e)*Pi[i]
    end
    #println([-( (6-:p1^2)/:f1 + (3-:p1^2 / 2) * :p1 )])
    #return [-( (6-:p1^2)/:f1 + (3-:p1^2 / 2) * :p1 )]
    return eom
end

function make_geodesic_eom(Pi,h)
    d = length(Pi)
    eom = Array{SymPy.Sym}(undef,d)
    for i in 1:d
        eom[i] = -Pi'*h[i,:,:]*Pi
    end
    eom
end

function make_pi_jac(Phi,Pi,H,eh,eom)
    vars = [Phi;Pi]#;H]
    d = length(Phi)
    l = length(vars)
    J = Array{Any}(undef,l,l)
    J[1:d,1:d] = zeros(d,d)
    J[1:d,d+1:2*d] = Matrix{Float64}(I,d,d)
    #J[end,1:d] = zeros(d)
    #J[end,d+1:2*d] = zeros(d)
    #J[end,end] = -eh
    for i in d+1:2*d
        for j in 1:l
            J[i,j] = df(eom[i-d],vars[j])
        end
    end
    return J
end

function make_pert_eom(t,G,p,mab,h,eh,loga,logk,hub)
    d = length(p)
    u = Array{SymPy.Sym}(undef,2*d,2*d)
    curv = Matrix{SymPy.Sym}(undef,2*d,2*d)
    u[1:d,1:d] = zeros((d,d))
    u[1:d,d+1:end] = Matrix{Float64}(I,d,d)
    #u[d+1:end,1:d] = -(Matrix{Int64}(I,d,d)*exp(2*(logk-t-loga-log(hub)))) - mab
    #u[d+1:end,d+1:end] = Matrix{Int64}(I,d,d)*(eh-3)
    for i in d+1:2*d
        for j in 1:d
            u[i,j] = -mab[i-d,j]
            if i-d == j
                u[i,j] -= exp(2*(logk-t-loga-log(hub)))
            end
        end
        for j in d+1:2*d
            u[i,j] = 0
            if i == j
                u[i,j] = eh-3
            end
        end
    end
    for i in 1:2*d, j in 1:2*d
        curv[i,j] = 0
        for m in 1:d, n in 1:d
            if i <= d
                curv[i,j] += h[i,m,n]*G[m,j]*p[n]
            else
                curv[i,j] += h[i-d,m,n]*G[m+d,j]*p[n]
            end
        #    #if j <= d
        #    #    curv[i,j] += h[j,m,n]*G[i,m]*p[n]
        #    #else
        #    #    curv[i,j] += h[j-d,m,n]*G[i,m+d]*p[n]
        #    #end
        end
        if simplification >= 2
            curv[i,j] = simp(curv[i,j])
            println(i,j)
        end
    end
    #println(u)
    #println(curv)
    return u*G .- curv
end

function make_pert_uab(t,p,mab,h,eh,loga,logk,hub)
    d = length(p)
    u = Array{SymPy.Sym}(undef,2*d,2*d)
    curv = Matrix{SymPy.Sym}(undef,2*d,2*d)
    u[1:d,1:d] = zeros((d,d))
    u[1:d,d+1:end] = Matrix{Float64}(I,d,d)
    u[d+1:end,1:d] = -(Matrix{Int64}(I,d,d)*exp(2*(logk-t-loga-log(hub)))) - mab
    u[d+1:end,d+1:end] = Matrix{Int64}(I,d,d)*(eh-3)

    return u
end
function make_pert_curv(p,h)
    d = length(p)
    hp = Array{SymPy.Sym}(undef,d,d)
    for i in 1:d, j in 1:d
        hp[i,j] = sum([h[i,j,k]*p[k] for k in 1:d])
    end
    return hp
end

function make_pert_jac(G,p,pert_eom)
    d = length(p)
    vars = [G...]
    eom = [pert_eom...]
    l = length(vars)
    J = Array{Any}(undef,l,l)
    for i in 1:l
        for j in 1:l
            J[i,j] = df(eom[i],vars[j])
        end
    end
    return J
end

function make_tensor_eom(t,G,eh,loga,logk,hub)
    w = [ 0 1; -exp(2*(logk-t-loga-log(hub))) (eh-3)]
    dG = [ sum(w[a,c]*G[c,b] + w[b,c]*G[a,c] for c in 1:2) for a in 1:2 for b in 1:2]
    return dG
end

function make_tensor_jac(G,teom)
    vars = [G...]
    eom = [teom...]
    l = length(vars)
    J = Array{Any}(undef,l,l)
    for i in 1:l
        for j in 1:l
            J[i,j] = df(eom[i],vars[j])
        end
    end
    return J
end

function hubble(v,e)
    return sqrt(v/(3-e))
end
#function chrπ(Phi,g,ginv=inv(g))
#    """
#    returns the matrix ``Γ^a_{bc} π^c``, since this can be more efficiently calculated
#    than the christoffel symbols themselves
#    """
#    d = length(Pi)
#    gp = Array{Sym}(undef,d,d)
#    for j in 1:d, k in 1:d
#        gp[j,k] = 0
#        for l in 1:d, m in 1:d
#            gp[j,k] += 1/2*(ginv[j,m]*(diff(g[m,l],Phi[k])+diff(g[k,m],Phi[l])-diff(g[k,l],Phi[m])))*Pi[l]
#        end
#    end
#    return gpp
#
#end
#
#function chrππ(Phi,Pi,g,ginv=inv(g))
#    """
#    returns the vector ``Γ^a_{bc} π^b π^c``, since this can be more efficiently calculated
#    than the christoffel symbols themselves
#    """
#    d = length(Pi)
#    gpp = Array{Sym}(undef,d)
#    for j in 1:d
#        gpp[j] = 0
#        for k in 1:d, l in 1:d, m in 1:d
#            gpp[j] += 1/2*(ginv[j,m]*(diff(g[m,l],Phi[k])+diff(g[k,m],Phi[l])-diff(g[k,l],Phi[m])))*Pi[k]*Pi[l]
#        end
#    end
#    return gpp
#end

function christoffel(Phi,g,ginv=inv(g),simplification=1)
    """
    returns christoffel symbols as ``Γ^a_{bc}``
    """
    d = length(Phi)
    h  = Array{SymPy.Sym}(undef,d,d,d)
    hl = Array{SymPy.Sym}(undef,d,d,d)
    for j in 1:d, k in 1:d, l in k:d
        h[j,k,l]=0
        println(j,k,l)
        for m in 1:d
            h[j,k,l] += 1/2*(ginv[j,m]*(diff(g[m,l],Phi[k])+diff(g[k,m],Phi[l])-diff(g[k,l],Phi[m]) ))
        end
        hl[j,k,l] = 1/2*(diff(g[j,l],Phi[k])+diff(g[k,j],Phi[l])-diff(g[k,l],Phi[j]) )
        hl[j,k,l] = cancel(hl[j,k,l])
        h[j,k,l]  = cancel(h[j,k,l])
        h[j,l,k]  = h[j,k,l]
        hl[j,l,k] = hl[j,k,l]
    end
    return h,hl
end

function Rππ_fast(Phi,Pi,h,hl,g,ginv;simplification=1)
    """
    compute a common contraction of the Riemann curvature tensor ``R^α_{βγδ} π^β π^γ``
    try to use some of the riemann tensor symmetries to help, internally use R_{αβγδ}
    """
    function copy_symm!(R,i,j,k,l)
        # antisymm on first and last pair, and symm on swapping pairs
        R[i,j,l,k] = -R[i,j,k,l]
        R[j,i,k,l] = -R[i,j,k,l]
        R[j,i,l,k] = R[i,j,k,l]
        R[k,l,i,j] = R[i,j,k,l]
        R[l,k,i,j] = -R[i,j,k,l]
        R[k,l,j,i] = -R[i,j,k,l]
        R[l,k,j,i] = R[i,j,k,l]
    end

    d = length(Phi)
    R = Array{SymPy.Sym}(undef,d,d,d,d)
    for i in 1:d
        R[i,i,:,:] = zeros(Int,d,d)
        R[:,:,i,i] = zeros(Int,d,d)
    end
    iterations = 0
    for i in 1:d, j in 1:d, k in 1:d, l in 1:d
        if !isassigned(R,i,j,k,l)
            # use lower index version for nice symmetries
            R[i,j,k,l] = -(diff(hl[l,j,k],Phi[i]) - diff(hl[l,i,k],Phi[j]))
            R[i,j,k,l] -= sum([(h[p,i,k]*hl[p,j,l] - h[p,j,k]*hl[p,i,l]) for p in 1:d])
            if simplification >= 2
                R[i,j,k,l] = simp(R[i,j,k,l])
            end
            # take care of (skew-)symmetries
            copy_symm!(R,i,j,k,l)
            #first Bianchi identity
            if !isassigned(R,k,l,i,j) && isassigned(R,j,k,l,i)
                R[k,l,i,j] = -R[i,j,k,l] - R[j,k,l,i]
                copy_symm!(R,k,l,i,j)
            elseif !isassigned(R,j,k,l,i) && isassigned(R,k,l,i,j)
                R[j,k,l,i] = -R[i,j,k,l] - R[k,l,i,j]
                copy_symm!(R,j,k,l,i)
            end
            iterations += 1
            println(iterations) # few over ideal, need to check correctness
        end
    end
    println("calculating Rππ")
    Rpp = Array{SymPy.Sym}(undef,d,d)
    for i in 1:d, l in 1:d
        Rpp[i,l] = 0
        for j in 1:d, k in 1:d, m in 1:d
            Rpp[i,l] += ginv[i,m]*R[m,j,k,l]*Pi[j]*Pi[k]
        end
        if simplification >= 2
            Rpp[i,l] = simp(Rpp[i,l])
        end
    end
    return Rpp
end

function riemann(Phi,h)
    """
    inefficiently compute Riemann curvature tensor as ``R^α_{βγδ}``
    """
    d = length(Phi)
    R = Array{Any}(undef,d,d,d,d)
    for i in 1:d, j in 1:d, k in 1:d, l in 1:d
        R[i,j,k,l] = diff(h[i,l,j],Phi[k]) - diff(h[i,k,j],Phi[l])
        for m in 1:d
            R[i,j,k,l] += h[i,k,m]*h[m,l,j] - h[i,l,m]*h[m,k,j]
        end
    end
    return R
end

function ricci_tensor(R)
    """
    compute Ricci tensor ``R_{μν} = R^α_{μαν}``
    """
    Ric = Array{Expr}(0,d,d)
    q = CartesianIndices(Ric)
    for (i,j) in q
        for k in 1:d
            Ric[i,j] += R[k,i,k,j]
        end
    end
    return Ric
end
function ricci_scalar(Ric,ginv)
    """
    compute Ricci scalar ``R = g^{μν} R_{μν}``
    """
    return tr(ginv*Ric)
end

function epsh(Pi,g)
    e = 0
    d = length(Pi)
    for i in 1:d
        for j in 1:d
            e += Pi[i] * g[i,j] * Pi[j] / 2
        end
    end
    return e
end

function epsv(gv,ginv,v)
    d = length(gv)
    ev = 0
    for i in 1:d
        for j in 1:d
            ev += ginv[i,j]*(gv[i]*gv[j]) / (2*v^2)
        end
    end
    return ev
end

function eta(Pi,e,H,va)
    d = length(Pi)
    vp = 0
    for i in 1:d
        vp += va[i]*Pi[i]
    end
    return -6 + 2*e - vp/(H^2 * e)
end

function omega(e,ev,eta)
    return -e^2 - 6*ev + 9*ev/e - 1/4*(6 + eta)^2 + e*(6 + ev + eta)
end

function omega2(fp,g,H,h,pi_eom,sig)
    d = length(pi_eom)
    Σ(x) = reduce(Algebra.:+,x)
    P(x) = reduce(Algebra.:*,x)
    fpdot = pi_eom[:] # careful with the copy here, missing this is catastrophic
    signorm = 0
    for i in 1:d
        for j in 1:d
            signorm += sig[i] * g[i,j] * sig[j]
        end
    end
    signorm = sqrt(signorm)
    t = sig ./ 2
    w = fpdot
    w1 = t*(Σ([P([fp[i],fpdot[j],g[i,j]]) + P([fpdot[i],fp[j],g[i,j]]) for i in 1:d for j in 1:d]))
    w2 = t*Σ([P([fp[i],fp[j],fp[k],(h[l,k,i]*g[l,j]+h[l,k,j]*g[i,l])]) for i in 1:d for j in 1:d for k in 1:d for l in 1:d])
    w .-= w2
    w .-= w1
    w *= (H/signorm)
    Om = 0
    wnorm = 0
    for i in 1:d
        for j in 1:d
            wnorm += w[i]*g[i,j]*w[j]
        end
    end
    Om = wnorm / H^2
    wnorm = (wnorm)^(:q)
    s = [wa / wnorm for wa in w]
    s = [sub("q=1/2",sa) for sa in s]
    return Om,s
end

function omega3(fp,gv,g,ginv,hub,eh)
    vu = ginv*gv
    vu = vu ./ sqrt(vu'*g*vu)
    Vv = sqrt(gv'*ginv*gv)
    perpuu = ginv .- vu'.*vu
    piperp = perpuu*g*fp
    pnorm = sqrt(piperp'*g*piperp)
    return (pnorm * Vv / (2*eh*hub^2))^2
end
function covggv(ggv,h,gv,simp,simplification)
    d = length(gv)
    Covggv = Array{SymPy.Sym}(undef,d,d)
    for i in 1:d
        for j in i:d
            Covggv[i,j] = ggv[i,j]
            for k in 1:d
                Covggv[i,j] -= h[k,i,j]*gv[k]
            end
            if simplification >= 1
                Covggv[i,j] = simp(Covggv[i,j])
            end
            Covggv[j,i] = Covggv[i,j]
        end
    end
    return Covggv
end

function mab(Covggv,Rpp,fp,hub,gv,eh,g,ginv,simp,simplification=1)
    d = length(gv)
    Mab = Array{SymPy.Sym}(undef,d,d) # actually Mab/H^2
    sa = Array{SymPy.Sym}(undef,d)
    for w in 1:d
        for b in 1:d
            for a in 1:d
                sa[a] = 0
                if ginv[w,a] != 0
                    s1 = ginv[w,a] * (Covggv[a,b]/hub^2)
                    s2 = -Rpp[w,b]
                    #s2 = -ginv[w,a]*sum([(g[a,m]*R[m,l,i,b])*(fp[l]*fp[i]) for l in 1:d for i in 1:d for m in 1:d])
                    s3 = -ginv[w,a]*(3+eh)*sum([g[a,m]*fp[m] for m in 1:d])*sum([g[b,n]*fp[n] for n in 1:d])
                    sum1 = sum([g[b,n]*fp[n] for n in 1:d])
                    sum2 = sum([g[a,m]*fp[m] for m in 1:d])
                    s4 = prod([ginv[w,a],2,eh,sum1,sum2])
                    s5 = -prod([ginv[w,a],sum2,(-gv[b]/hub^2 + sum1*(eh-3))])
                    s6 = -prod([ginv[w,a],sum1,(-gv[a]/hub^2 + sum2*(eh-3))])
                    if simplification >= 3
                        sa[a] = simp(s6+s5+s4+s3+s2+s1)
                    else
                        sa[a] = s6+s5+s4+s3+s2+s1
                    end
                end
            end
            Mab[w,b] = sum(sa)
        end
    end
    return Mab
end

function pz_hca(v,ea,va)
    # see 1303.3611 for these
    ua = [vai / v for vai in va]
    return v / (24 * pi^2) * sum((ua.^2 ./ ea))
end
function ns_hca(v,ea,etaa,eh,va)
    d = length(va)
    # see 1303.3611 for these
    ua = [vai / v for vai in va]
    num = (1-sum([etaa[i] * ua[i]^2 / (2*ea[i]) for i in 1:d]))
    den = sum([ua[i]^2 / ea[i] for i in 1:d])
    ns = 1 - 2*eh - 4*num/den
    return ns
end
function sig(fp,g)
    d = length(fp)
    norm = 0
    for i in 1:d
        for j in 1:d
            norm += fp[i]*g[i,j]*fp[j]
        end
    end
    norm = sqrt(norm)
    return [p / norm for p in fp]
end
# this is wrong, only works for 2-field
function s(ginv,g,sig)
    d = length(sig)
    sij = Array{Any}(undef,d,d)
    for i in 1:d
        for j in 1:d
            sij[i,j] = ginv[i,j] - sig[i]*sig[j]
        end
    end
    s = sij[:,d]
    norm = 0
    for i in 1:d
        for j in 1:d
            norm += s[i]*g[i,j]*s[j]
        end
    end
    return [si / norm for si in s]
end

function masses(g,Mab,Sig,S)
    d = length(Sig)
    Mσσ = 0
    for i in 1:d
        for j in 1:d
            for k in 1:d
                Mσσ += g[i,k]*Mab[k,j]*Sig[i]*Sig[j]
            end
        end
    end
    Mss = 0
    if d > 1
        for i in 1:d
            for j in 1:d
                for k in d:1
                    Mss += (g[i,k]*Mab[k,j])*S2[i]*S2[j]
                end
            end
        end
    end
    return Mσσ, Mss
end


text:
using Test
using SynthControl, CSV, DataFrames, Dates

@testset "brexit" begin
  df = load_brexit()

  tp = TreatmentPanel("United Kingdom" => Date(2016, 7, 1), df;
                       outcome = :realgdp, id_var = :country, t_var = :quarter)

  s = SynthControlModel(tp)

  @test !(isfitted(s))

  fit!(s)

  @test isfitted(s)

  @test s.tp.no_comps == 22

  @test s.tp.no_pretreat_t == 29

  @test size(s.tp.comps) == (22, 30)

end


text:
#=
License: MIT
Description: Definition of Store for operational semantics
Copyright (c) 2020 thautwarm <twshere@outlook.com> and contributors
=#
export Store

ImD = Base.ImmutableDict
struct Store{K, O}
    unbox :: ImD{K, O}
end

# given σ ∈ Store{K, O} where K, O
# 'σ[k => v]' add a scoped pair k=>v and returns a new store
# where the old one is not affected.
function Store{K, O}() where {K, O}
    Store(ImD{K, O}())
end

Base.getindex(s::Store{K, O}, x::Vector{Pair{K, A}}) where {K, O, A<:O} =
    let im = s.unbox
        for each in x
            im = ImD(im, each)
        end
        Store(im)
    end

Base.getindex(s::Store{K, O}, x::Pair{K, A}) where {K, O, A<:O} = s[[x]]
Base.getindex(s::Store{K, O}, x::K) where {K, O} = s.unbox[x]

Base.get(s::Store, k, default) = get(s.unbox, k, default)
Base.in(k::K, s::Store{K}) where K = haskey(s.unbox, k)


text:
import Base:
	length,
	isempty

"""
We treat TextNode as an empty collection. While it might seem like it could act on the text content, this makes
the behaviour of TextNode and Element unexpectedly different.
"""
Base.length(::TextNode) = 1
Base.length(e::Element) = i + map(Base.length, e.children) |> sum

Base.isempty(::TextNode) = true
Base.isempty(e::Element) = isempty(e.children)

"""
I'm tempted to add AbstractArray functions pointing at Element.children, and Dict structures at Element.props,
but that's quite a bit of boilerplate.
Should be able to do it with a macro!
"""

children(e::Element) = e.children
props(e::Element) = e.props

# Dict ops that take the collection as the first key, and modifies a collection in place
# These are generic on Element, allowing more per-element type overrides to do their own thing
props_ops = [
	:haskey,
	:get,
	:get!
	:getkey,
	:delete!,
	:pop!,
	:keys,
	:values,
	:pairs
	:merge!
	:sizehint!
	:keytype,
	:valtype
]
@eval begin
	$(
		map( op -> :(import Base.$op), props_ops)...
	)
	$(
		map( op -> :(Base.$op(e::Element, x...) = Base.$op(e.props, x...) ), props_ops )...
	)
end

# Ops that take two Dicts should be implemented for two Elements, where it makes sense?


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Markdown

using Compat

## Backport `Pair`. ##

immutable Pair{A,B}
    first::A
    second::B
end

Base.start(p::Pair) = 1
Base.done(p::Pair, i) = i>2
Base.next(p::Pair, i) = (getfield(p,i), i+1)

Base.indexed_next(p::Pair, i::Int, state) = (getfield(p,i), i+1)

Base.hash(p::Pair, h::UInt) = hash(p.second, hash(p.first, h))

==(p::Pair, q::Pair) = (p.first==q.first) & (p.second==q.second)
Base.isequal(p::Pair, q::Pair) = isequal(p.first,q.first) & isequal(p.second,q.second)

Base.isless(p::Pair, q::Pair) = ifelse(!isequal(p.first,q.first), isless(p.first,q.first),
                                                             isless(p.second,q.second))
Base.getindex(p::Pair,i::Int) = getfield(p,i)
Base.getindex(p::Pair,i::Real) = getfield(p, convert(Int, i))
Base.reverse(p::Pair) = Pair(p.second, p.first)

## End backport. ##

import Base: writemime

typealias String AbstractString

include("parse/config.jl")
include("parse/util.jl")
include("parse/parse.jl")

include("Common/Common.jl")
include("GitHub/GitHub.jl")
include("IPython/IPython.jl")
include("Julia/Julia.jl")

include("render/plain.jl")
include("render/html.jl")
include("render/latex.jl")

include("render/terminal/render.jl")

export readme, license, @md_str, @doc_str, @md_mstr, @doc_mstr

parse(markdown::String; flavor = julia) = parse(IOBuffer(markdown), flavor = flavor)
parse_file(file::String; flavor = julia) = parse(readall(file), flavor = flavor)

readme(pkg::String; flavor = github) = parse_file(Pkg.dir(pkg, "README.md"), flavor = flavor)
readme(pkg::Module; flavor = github) = readme(string(pkg), flavor = flavor)

license(pkg::String; flavor = github) = parse_file(Pkg.dir(pkg, "LICENSE.md"), flavor = flavor)
license(pkg::Module; flavor = github) = license(string(pkg), flavor = flavor)

function mdexpr(s, flavor = :julia)
    md = parse(s, flavor = symbol(flavor))
    esc(toexpr(md))
end

function docexpr(s, flavor = :julia)
    quote
        let md = $(mdexpr(s, flavor))
            md.meta[:path] = @__FILE__
            md.meta[:module] = current_module()
            md
        end
    end
end

macro md_str(s, t...)
    mdexpr(s, t...)
end

macro doc_str(s, t...)
    docexpr(s, t...)
end

macro md_mstr(s, t...)
    mdexpr(Base.triplequoted(s), t...)
end

macro doc_mstr(s, t...)
    docexpr(Base.triplequoted(s), t...)
end

function Base.display(d::Base.REPL.REPLDisplay, md::Vector{MD})
    for md in md
        display(d, md)
    end
end

end


text:
using SpatialFields
using Base.Test
import MultiPoly

function test_polynomial()
  x, y = MultiPoly.generators(MultiPoly.MPoly{Float64}, :x, :y)
  poly = x^2 + 3x*y + y + y^2 + 10

  field = PolynomialScalarField(poly)

  @test isapprox(field([3., 7.]), MultiPoly.evaluate(poly, 3., 7.))

  vec_field = SpatialFields.gradient(field)
  x = [100, 200]
  @test isapprox(vec_field(x),
                 Float64[MultiPoly.evaluate(MultiPoly.diff(poly, v),
                                     x[1], x[2]) for v in [:x, :y]])
end

test_polynomial()


text:
#####
##### Status Types
#####

# We use a type based strategy for controlling access to a CachedArray.
# Transitions from one status to another (e.g. NotBusy -> ReadOnly) will potentially
# trigger state updates in the manager/policy and record optional telemetry.
abstract type AbstractStatus end
struct NotBusy <: AbstractStatus end
struct ReadOnly <: AbstractStatus end
struct ReadWrite <: AbstractStatus end

_sym(::Type{NotBusy}) = :NotBusy
_sym(::Type{ReadOnly}) = :ReadOnly
_sym(::Type{ReadWrite}) = :ReadWrite

const Readable = Union{ReadOnly,ReadWrite}
const Writable = ReadWrite

######
###### CachedArray
######

struct CachedArray{T,N,S<:AbstractStatus,M} <: DenseArray{T,N}
    object::Object{M}
    dims::NTuple{N,Int}

    # Inner constructor - do a type check and make sure the finalizer is attached.
    function CachedArray{T,N}(
        object::Object{M}, dims::NTuple{N,Int}, status::S = NotBusy()
    ) where {T,N,S,M}
        if !isbitstype(T)
            throw(ArgumentError("Cannot construct a `CachedArray` from non-isbits types!"))
        end
        return new{T,N,S,M}(object, dims)
    end
end

# Aliases
const ReadableCachedArray{T,N} = CachedArray{T,N,<:Readable,<:Any}
const UnreadableCachedArray{T,N} = CachedArray{T,N,NotBusy,<:Any}
const WritableCachedArray{T,N} = CachedArray{T,N,<:Writable,<:Any}
const UnwritableCachedArray{T,N} = CachedArray{T,N,<:Union{NotBusy,ReadOnly},<:Any}
const BusyCachedArray{T,N} = CachedArray{T,N,<:Union{ReadOnly,ReadWrite},<:Any}

isreadonly(::CachedArray) = false
isreadonly(::UnwritableCachedArray) = true

# Cached API
object(A::CachedArray) = A.object
metadata(A::CachedArray) = metadata(object(A))
datapointer(A::CachedArray) = datapointer(object(A))
manager(A::CachedArray) = manager(object(A))
unsafe_free(A::CachedArray) = unsafe_free(object(A))

# Escape hatch to ALWAYS get a pointer, regardless of the status of the array.
# Should only be called directly by `Base.pointer` or by the `CacheManager`.
@inline unsafe_pointer(A::CachedArray{T}) where {T} = Ptr{T}(pointer(object(A)))

# Don't let functions normally take pointers to un-acquired CachedArrays.
Base.pointer(A::CachedArray) = unsafe_pointer(A)
function Base.pointer(A::UnreadableCachedArray)
    return error("Cannot take a pointer to to a `NotBusy` CachedArray")
end

# Consructors
CachedArray(x::Array{T,N}, manager; kw...) where {T,N} = CachedArray{T,N}(x, manager; kw...)
function CachedArray{T,N}(
    x::Array{T,N}, manager; status = ReadWrite(), priority = PreferLocal
) where {T,N}
    object = alloc(manager, sizeof(x))
    unsafe_copyto!(Ptr{T}(pointer(object)), pointer(x), length(x))
    return CachedArray{T,N}(object, size(x), status)
end

function CachedArray{T}(::UndefInitializer, manager, i::Integer...; kw...) where {T}
    return CachedArray{T}(undef, manager, convert.(Int, i); kw...)
end

function CachedArray{T}(
    ::UndefInitializer,
    manager,
    dims::NTuple{N,Int};
    status = NotBusy(),
    priority = PreferLocal,
) where {T,N}
    object = alloc(manager, prod(dims) * sizeof(T), priority)
    return CachedArray{T,N}(object, dims, status)
end

#####
##### Array Interface
#####

# Go through the "pointer" API to capture "NotBusy" arrays.
Base.unsafe_convert(::Type{Ptr{T}}, A::CachedArray{T}) where {T} = pointer(A)
Base.sizeof(A::CachedArray) = prod(size(A)) * sizeof(eltype(A))
@inline Base.size(A::CachedArray) = A.dims

Base.elsize(::CachedArray{T}) where {T} = sizeof(T)
Base.elsize(::Type{<:CachedArray{T}}) where {T} = sizeof(T)

function Base.getindex(A::ReadableCachedArray, i::Int)
    @boundscheck checkbounds(A, i)
    return LoadStore.unsafe_custom_load(pointer(A), i)
end

function Base.setindex!(A::WritableCachedArray, v, i::Int)
    @boundscheck checkbounds(A, i)
    return LoadStore.unsafe_custom_store!(pointer(A), v, i)
end

Base.IndexStyle(::Type{<:CachedArray}) = Base.IndexLinear()

function Base.similar(
    A::CachedArray,
    ::Type{T} = eltype(A),
    dims::Tuple{Vararg{Int,N}} = size(A);
    status = ReadWrite(),
    priority = ForceLocal,
) where {T,N}
    return CachedArray{T}(undef, manager(A), dims; status = status, priority = priority)
end

function Base.iterate(A::ReadableCachedArray, i::Int = 1)
    i > length(A) && return nothing
    return (@inbounds A[i], i + 1)
end

# For alias detection
Base.dataids(A::CachedArray) = (UInt(pointer(A)),)

function Base.reshape(x::CachedArray{T,M,S}, dims::NTuple{N,Int}) where {T,N,M,S}
    function throw_dmrsa(dims, len)
        return throw(
            DimensionMismatch(
                "new dimensions $(dims) must be consistent with array size $len"
            ),
        )
    end

    if prod(dims) != length(x)
        throw_dmrsa(dims, length(x))
    end

    if N == M && dims == size(x)
        return x
    end

    # Keep the same underlying object.
    return CachedArray{T,N}(object(x), dims, S())
end

#####
##### Speedup functions
#####

function Base.copy(x::CachedArray)
    y = similar(x)
    unsafe_copyto!(pointer(y), pointer(x), length(x))
    return y
end

#####
##### Broadcasting
#####

# Hijack broadcasting so we prioritize CachedArrays.
function Base.BroadcastStyle(::Type{T}) where {T<:CachedArray}
    return Broadcast.ArrayStyle{CachedArray}()
end

function Base.BroadcastStyle(::Type{T}) where {U<:CachedArray,T<:SubArray{<:Any,<:Any,U}}
    return Broadcast.ArrayStyle{CachedArray}()
end

# TODO - right now, we traverse through the BC object to find the first intance of an
# CachedArray.
#
# A more general solution would gather all such arrays and check that all the CachedArrays
# are the same.
function Base.similar(
    bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{CachedArray}}, ::Type{ElType}
) where {ElType}
    cached = findT(CachedArray, bc)
    return similar(cached, ElType, axes(bc))
end

#####
##### Conversion Functions
#####

const __fnmap = [:NotBusy => :release, :ReadOnly => :readable, :ReadWrite => :writable]

const __updates = Dict(
    # Implications of making an array notbusy.
    # 1. TODO:
    :release => [],
    # Implications of making an array readable.
    # 1. TODO: Usage should be updated
    :readable => [],
    # Implications of making an array writable.
    # 1. Array should be marked as dirty.
    # 2. TODO: Usage should be updated
    :writable => [:(setdirty!(x))],
)

for (typ, fn) in __fnmap
    # No-op if already correct type.
    @eval $fn(x::CachedArray{<:Any,<:Any,$typ}) = x
    @eval function $fn(x::CachedArray{T,N}) where {T,N}
        if !isnull(unsafe_pointer(x.object))
            # Optional Telemetry
            @telemetry manager(x) begin
                telemetry_change(gettelemetry(manager(x)), getid(metadata(x)), $(QuoteNode(typ)))
            end
            # unpack any potential policy updates.
            $(__updates...)
        end
        return CachedArray{T,N}(object(x), size(x), $typ())
    end

    # Define `Base.convert` in terms of the "release, readable, writable" methods
    # defined above.
    @eval function Base.convert(
        ::Type{CachedArray{T,N,$typ,M}}, x::CachedArray{T,N,S,M}
    ) where {T,N,S,M}
        return $fn(x)
    end
end


text:
##
using Revise
using GLMakie
import Pkg
Pkg.activate(".")
using Images
using FileIO

function preprocess_image(img_path)
    img = load(img_path)
    return @. Float32(Images.Gray(img))
end
##
desmin, actin = map(preprocess_image, ["./Data/SR2_desmin_image.jpg", "./Data/SR3_actin_image.jpg"])
##
let
    fig = Figure(resolution = (1500, 800))
    ax = GLMakie.Axis(fig[1, 1])
    image!(ax, actin)
    fig
end
desmin = desmin[end-2987+1:end, end-3004+1:end]
##
let
    desmin_vec, actin_vec = map(arr -> reshape(arr, length(arr)), [desmin, actin])
    fig = Figure(resolution = (1500, 800))
    ax = GLMakie.Axis(fig[1, 1])
    scatter!(ax, desmin_vec, actin_vec)
    fig
end


text:
#=
If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.
Find the sum of all the multiples of 3 or 5 below 1000.
=#

function calc()
  filter(x -> x % 3 == 0 || x % 5 == 0, 1:999) |> sum |> println
end
@time calc()

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# -------- Bayes ELO Objective ----------------------------------------------- #

elosigm(d) = 1 / (1 + exp(d/400))

function likelihood(elo1, elo2, res; draw, start)

  p1 = elosigm(elo2 - elo1 - start + draw)
  p2 = elosigm(elo1 - elo2 + start + draw)

  res == 1 ? p1 : (res == -1 ? p2 : 1 - p1 - p2)

end

function loglikelihood(elos, results; kwargs...)

  sum(CartesianIndices(results)) do c
    (i,j,k) = Tuple(c)
    results[i,j,k] * log(likelihood(elos[i], elos[j], k-2; kwargs...))
  end

end

logprior(elos; mean, std) = sum(elo -> -0.5 * ((elo - mean) / std)^2, elos)

function objective(elos, results; mean, std, kwargs...)

  loglikelihood(elos, results; kwargs...) +
  logprior(elos; mean = mean, std = std)

end

# -------- ELO Estimation ---------------------------------------------------- #

function estimate( results  # result array (see function compete below)
                 ; start = nothing
                 , draw = nothing
                 , mean = 0.
                 , std = 1500. )

  # Get the number of players
  k = size(results, 1)

  # Upper and lower bounds and start values for (start, draw)
  lb = [-4std, 1]
  ub = [ 4std, 4std]
  st = [ 0., 1.]

  # If start or draw are given, we do not estimate them but set
  # upper bound = lower bound = start value = given value
  !isnothing(start) && (st[1] = lb[1] = ub[1] = start)
  !isnothing(draw)  && (st[2] = lb[2] = ub[2] = draw)

  # Create the optimizer object
  opt = NLopt.Opt(:LN_BOBYQA, k+2)

  # Define the minimization objective function compatible to NLopt
  f = (x, _) -> -objective( x[1:k]
                          , results
                          ; draw = x[k+2]
                          , start = x[k+1]
                          , mean = mean
                          , std = std )

  # Set the objective and set reasonable lower/upper bounds
  NLopt.min_objective!(opt, f)
  NLopt.lower_bounds!(opt, [fill(mean - 4std, k); lb])
  NLopt.upper_bounds!(opt, [fill(mean + 4std, k); ub])

  # Conduct the optimization
  value, params, ret = NLopt.optimize(opt, [fill(mean, k); st] )

  # Return the found values
  params

end

# -------- Competitions ------------------------------------------------------ #

function match_repetitions(n, pairings)

  n < pairings && @warn "Not every pair competes"

  matches = zeros(Int, pairings)
  matches .+= floor(Int, n / pairings)
  matches[1:(n % pairings)] .+= 1

  shuffle(matches)

end

pairings(l, a) = l * l - (l-a)*(l-a) - a

"""
  compete(players, n [, active]; <keyword arguments>)

Create game results for a number of `n` games between players at positions `i`
and `j` in the collection `players` if at least one of `i` or `j` is in
`active`.

The playings start with `game`, which is infered automatically from
`players` if possible. `callback()` is called after each of the `n` matches.

The return value is a result array of dimensions `(l, l, 3)`, where `l
= length(players)`. The entry at indices `[i,j,k]` stands for the number of
outcomes `k` (`0`: loss, `1`: draw, `2`: win) when player `i` played against
`j`. When both `i` and `j` are inactive, no games are played. The entries at
`[i, j, :]` are `0` in this case.

# Arguments
- `game`: Initial game state. Inferred automatically if possible.
- `callback`: Function called after each of the `n` matches.
- `distributed = false`: Whether to conduct the self-playings on several
   processes in parallel. If `true`, all available workers are used. Alternatively,
   a list of worker ids can be passed.
- `tickets`: Number of chunks in which the workload is distributed if
   `distributed != false`. By default, it is set to the number of workers.
"""
function compete( players
                , n :: Int
                , active = 1:length(players)
                ; game = derive_gametype(players)()
                , callback = () -> nothing
                , distributed = false
                , tickets = nothing )

  # If the competitions takes place in parallel, get the list of workers
  # and cede to the corresponding function in distributed.jl
  if distributed != false

    workers = distributed == true ? Distributed.workers() : distributed
    tickets = isnothing(tickets) ? length(workers) : tickets

    return compete_distributed( players, n, active
                              ; game = game
                              , callback = callback
                              , workers = workers
                              , tickets = tickets )

  end

  # Get the number of pairings that compete
  m = pairings(length(players), length(active))

  matches = match_repetitions(n, m)

  results = zeros(Int, length(players), length(players), 3)
  players = enumerate(players)

  l = 1

  for (i, p1) in players, (j, p2) in players

    if i != j && (i in active || j in active)

      asyncmap(1:matches[l]) do _
        k = pvp(p1, p2, game = game) + 2 # convert -1, 0, 1 to indices 1, 2, 3
        results[i, j, k] += 1
        callback()
      end

      l += 1

    end

  end

  results

end

# -------- Rankings ---------------------------------------------------------- #

"""
Structure that contains the ranking information of a group of competing players.
Based on the results of a competition, elo values and the start advantage as
well as draw bandwidth is determined.
"""
struct Ranking

  players :: Vector{Player}
  results :: Array{Int, 3}

  elos :: Array{Float64}

  draw :: Float64
  start :: Float64

end

"""
    Ranking(players, results; <keyword arguments>)
    Ranking(players, n [, active]; <keyword arguments>)

Get a ranking of `players` based on `results` from a previous call of the
function `compete`. Alternatively, results are generated on the fly via
`compete` for `n` games with all corresponding keyword arguments.

# Arguments
The remaining arguments alter how the Bayesian elo values are calculated.
- `mean = 0`: Mean value of the prior elo distribution.
- `std = 1500`: Standard deviation of the prior elo distribution.
- `start = nothing`: If not `nothing`, the start advantage is fixed to this value.
- `draw = nothing`: If not `nothing`, the draw bandwidth is fixed to this value.
"""
function Ranking(players, results :: Array{Int, 3}; kwargs...)

  est = estimate(results; kwargs...)

  k = size(results, 1)
  elos  = est[1:k]
  start = est[k+1]
  draw  = est[k+2]

  Ranking(players, results, elos, draw, start)

end

function Ranking( players
                , args...
                ; game = derive_gametype(players)()
                , callback = () -> nothing
                , distributed = false
                , tickets = nothing
                , kwargs...)

  results = compete( players
                   , args...
                   ; game = game
                   , callback = callback
                   , distributed = distributed
                   , tickets = tickets )

  Ranking(players, results; kwargs...)

end

# -------- Visualize Rankings ------------------------------------------------ #

function Base.summary(rk :: Ranking, matrix = false)

  # The player with highest elo will come first
  perm = sortperm(rk.elos) |> reverse

  ranks = map(1:length(perm), rk.players[perm], rk.elos[perm]) do i, p, elo
    Printf.@sprintf "%2d. %10.2f  %15s" i elo name(p)
  end

  # TODO: variable widths of the matrix entries, based on the length of the
  # wins/losses
  if matrix

    mat = rk.results[perm, perm, 3] - rk.results[perm, perm, 1]

    for (i,r) in enumerate(ranks)
      vals = map(mat[i,:], 1:length(perm)) do x, j
        i == j ? "    ●" : Printf.@sprintf("%5.0f", x)
      end
      ranks[i] = r * "  " * join(vals)
    end

  end

  start = Printf.@sprintf "start advantage: %.2f" rk.start
  draw  = Printf.@sprintf "draw bandwidth:  %.2f" rk.draw

  join([ranks; ""; start; draw], "\n")

end

# -------- Distributed Contesting -------------------------------------------- #

function compete_distributed(args...; kwargs...)

  sum(with_workers(compete, args...; kwargs...))

end

max_stars_repo_name: roSievers/jtac


text:
# We explore the batch and fidelity here

@everywhere include("../INCLUDEME.jl")

@everywhere using WFLearning

@everywhere function main(n; nbatch=20, epochs=200, learning_rate=0.1)
    pmap(nlayers->task("layers", n, nlayers, nbatch, epochs, learning_rate), 5:15)
end

main(6, nbatch=20, epochs=200, learning_rate=0.1)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
# using CuArrays  # Enable GPU
using Ju
using Flux
using Plots
gr()

env = CartPoleEnv()
input_length = length(sample(observationspace(env)))
output_length = length(actionspace(env))

ϵ = exp_decay(init=1.0, λ=1.0, decay_step=500, clip=0.01)
γ, buffer_size, batch_size = 0.99, 10^3, 32

Qₐ = NeuralNetworkQ(
        Chain(Dense(input_length, 128, relu),
              Dense(128, 128, relu),
              Dense(128, output_length)),
        ADAM())

Qₜ = NeuralNetworkQ(
        Chain(Dense(input_length, 128, relu),
              Dense(128, 128, relu),
              Dense(128, output_length)),
        ADAM())

Flux.loadparams!(Qₜ.model, Qₐ.ps)

learner = DoubleDQN(Qₐ, Qₜ, EpsilonGreedySelector(ϵ); update_period=1, target_update_period=100, γ=γ, batch_size=batch_size)
buffer = CircularSARDBuffer(buffer_size; state_type=Array{Float64, 1}, state_size=(input_length,))
agent = Agent(learner, buffer)

function record_loss()
    loss = []
    function f(env, agent)
        push!(loss, agent.learner.loss)
    end
    function f()
        loss
    end
end

callbacks = (stop_at_step(10000), rewards_of_each_episode(), record_loss())
train!(env, agent;callbacks=callbacks)


text:
int_rules_1_2_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.1.1*(a+b*x+c*x^2)^p =#
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => (2 * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / ((2 * ~p + 1) * (~b + 2 * ~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~p, -1))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~b / 2 + ~c * ~x) / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) * Antiderivative(1 / (~b / 2 + ~c * ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~b + 2 * ~c * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p) / (2 * ~c * (2 * ~p + 1)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ(~p, -1 / 2))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (1 / (~c) ^ ~p) * Antiderivative(Simp((~b / 2 - q / 2) + ~c * ~x, ~x) ^ ~p * Simp(~b / 2 + q / 2 + ~c * ~x, ~x) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~p, 0) && PerfectSquareQ((~b) ^ 2 - 4 * ~a * ~c)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~p, 0) && (EqQ(~a, 0) || Not(PerfectSquareQ((~b) ^ 2 - 4 * ~a * ~c)))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~b + 2 * ~c * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p) / (2 * ~c * (2 * ~p + 1)) - ((~p * ((~b) ^ 2 - 4 * ~a * ~c)) / (2 * ~c * (2 * ~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~p, 0) && IntegerQ(4 * ~p)))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ (3 / 2), ~x) => (-2 * (~b + 2 * ~c * ~x)) / (((~b) ^ 2 - 4 * ~a * ~c) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~b + 2 * ~c * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) - ((2 * ~c * (2 * ~p + 3)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && (NeQ(~p, -3 / 2) && IntegerQ(4 * ~p))))
    @apply_utils Antiderivative(1 / (~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => log(~x) / ~b - log(RemoveContent(~b + ~c * ~x, ~x)) / ~b <-- FreeQ([~b, ~c], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~c / q) * Antiderivative(1 / Simp((~b / 2 - q / 2) + ~c * ~x, ~x), ~x) - (~c / q) * Antiderivative(1 / Simp(~b / 2 + q / 2 + ~c * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (PosQ((~b) ^ 2 - 4 * ~a * ~c) && PerfectSquareQ((~b) ^ 2 - 4 * ~a * ~c)))
    @apply_utils Antiderivative(1 / (~a + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => With([q = 1 - 4 * Simplify((~a * ~c) / (~b) ^ 2)], (-2 / ~b) * Subst(Antiderivative(1 / (q - (~x) ^ 2), ~x), ~x, 1 + (2 * ~c * ~x) / ~b) <-- RationalQ(q) && (EqQ(q ^ 2, 1) || Not(RationalQ((~b) ^ 2 - 4 * ~a * ~c)))) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => -2 * Subst(Antiderivative(1 / Simp(((~b) ^ 2 - 4 * ~a * ~c) - (~x) ^ 2, ~x), ~x), ~x, ~b + 2 * ~c * ~x) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => (1 / (2 * ~c * ((-4 * ~c) / ((~b) ^ 2 - 4 * ~a * ~c)) ^ ~p)) * Subst(Antiderivative(Simp(1 - (~x) ^ 2 / ((~b) ^ 2 - 4 * ~a * ~c), ~x) ^ ~p, ~x), ~x, ~b + 2 * ~c * ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && GtQ(4 * ~a - (~b) ^ 2 / ~c, 0)
    @apply_utils Antiderivative(1 / sqrt(~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => 2 * Subst(Antiderivative(1 / (1 - ~c * (~x) ^ 2), ~x), ~x, ~x / sqrt(~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~b, ~c], ~x)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => 2 * Subst(Antiderivative(1 / (4 * ~c - (~x) ^ 2), ~x), ~x, (~b + 2 * ~c * ~x) / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~b * ~x + ~c * (~x) ^ 2) ^ ~p / ((-(~c) * (~b * ~x + ~c * (~x) ^ 2)) / (~b) ^ 2) ^ ~p) * Antiderivative(((-(~c) * ~x) / ~b - ((~c) ^ 2 * (~x) ^ 2) / (~b) ^ 2) ^ ~p, ~x) <-- FreeQ([~b, ~c], ~x) && (RationalQ(~p) && 3 <= Denominator(~p) <= 4)
    #= Antiderivative(((~a')+(~b')*(~x)+(~c')*(~x)^2)^(~p),~x) :=  (a+b*x+c*x^2)^p/(-c*(a+b*x+c*x^2)/(b^2-4*a*c))^p*Antiderivative((-a*c/(b^2-4*a*c) -b*c*x/(b^2-4*a*c)-c^2*x^2/(b^2-4*a*c))^p,x) <-- FreeQ([a,b,c],x) &&  NeQ(b^2-4*a*c,0) && RationalQ(p) && 3<=Denominator(p)<=4 =#
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([d = Denominator(~p)], ((d * sqrt((~b + 2 * ~c * ~x) ^ 2)) / (~b + 2 * ~c * ~x)) * Subst(Antiderivative((~x) ^ (d * (~p + 1) - 1) / sqrt(((~b) ^ 2 - 4 * ~a * ~c) + 4 * ~c * (~x) ^ d), ~x), ~x, (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (1 / d)) <-- 3 <= d <= 4) <-- FreeQ([~a, ~b, ~c], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && RationalQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (-((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (q * (~p + 1) * (((q - ~b) - 2 * ~c * ~x) / (2q)) ^ (~p + 1))) * HypergeometricFunctions._₂F₁(-(~p), ~p + 1, ~p + 2, (~b + q + 2 * ~c * ~x) / (2q))) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(4 * ~p)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~u + ~(c') * (~u) ^ 2) ^ ~p, ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.1 (a+b x+c x^2)^p.jl


text:
"""
    hydrolysis_cos(temp[, pH=7., seawater::Bool=false])

COS abiotic hydrolysis rate in natural waters. Applicable range: temperature
5--30 C and pH 4--10.

# Parameters

* `temp`: Temperature in Kelvin.
* `pH`: pH value. Default is 7.
* `seawater` (*optional*): If `true`, calculate COS hydrolysis rate in the
  seawater [RC93]. If `false` (default), caculate COS hydrolysis rate in the
  fresh water [E89].

# Returns

The abiotic hydrolysis rate constant of COS [s^-1].

# References

[E89] Elliott, S., Lu, E., and Rowland, F. S. (1989). Rates and
    mechanisms for the hydrolysis of carbonyl sulfide in natural waters.
    *Environ. Sci. Tech.*, 23(4), 458--461.
    <https://doi.org/10.1021/es00181a011>

[RC93] Radford-Knoery, J., and Cutter, G. A. (1993). Determination of
    carbonyl sulfide and hydrogen sulfide species in natural waters using
    specialized collection procedures and gas chromatography with flame
    photometric detection. *Anal. Chem.*, 65(8), 976--982.
    <https://doi.org/10.1021/ac00056a005>

# Examples

```jldoctest
julia> @printf "%g" hydrolysis_cos(288.15)
6.60184e-06

julia> @printf "%g" hydrolysis_cos(288.15, 9.)
3.64899e-05

julia> @printf "%g" hydrolysis_cos(288.15, 8.2, true)
1.40695e-05
```
"""
function hydrolysis_cos(temp, pH=7., seawater::Bool=false)
    T_ref = 298.15
    c_OH = 10. ^ (pH - water_dissoc(temp))
    t_diff = 1. / temp - 1. / T_ref

    # note: the parameter values used here differ from those in the original
    # publications due to a small correction for water pK_w
    if seawater
        # RC93 equation
        return 1.63838819502e-5 * exp(-6444.02904777 * t_diff) +
            11.7115101829 * exp(-2427.27401921 * t_diff) * c_OH
    else
        # E89 equation
        return 2.11834513803e-5 * exp(-10418.3722377 * t_diff) +
            14.16881179 * exp(-6469.11889197 * t_diff) * c_OH
    end
end


text:
print("hello,world!\n")
x = 3
y = 0x233
v = x + y
print(y)
typeof(v)
s = "zhangqishneg"
print(s)
print("fsf$s")
print("fsf\$s")

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
### Note this code is outdated
using Plots
using Convex
using Mosek
using EllZeroTrendFiltering

using Interpolations

g = snp500_data()[1:2000]

N = length(g)

# Transition costs (for recvoering the y-values)
ℓ = compute_discrete_transition_costs(g)
V_N = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)

M = 50

#f = open("data.jld", "w+"); serialize(f, Λ); close(f)
#@time Λ =  pwq_dp_constrained(ℓ, V_N, M)

sols_opt = []
for m=1:M
    I_opt, Y_opt, f_opt = recover_solution(Λ[m, 1], ℓ, V_N)
    push!(sols_opt, (I_opt, Y_opt, f_opt))
end

# Iterative re-weighting according to Kim 2009, et. al
function iterative_l1_reweighting(g, max_cost, ϵ)
    N = length(g)
    H =  spdiagm((ones(N-2), -2*ones(N-2), ones(N-2)), (0,1,2))

    W = spdiagm(ones(N-2))

    x_sol = zeros(size(N))
    f_sol = NaN

    for iter=1:10
        x = Variable(N)
        problem = minimize(norm((W*H)*x, 1))
        problem.constraints += sumsquares(x - g) < max_cost

        solve!(problem, MosekSolver(), verbose=false)
        println(problem.status)

        x_sol = evaluate(x)[:]
        f_sol = sum((x_sol - g).^2)

        W = spdiagm(ones(N-2) ./ (ϵ .+ H*x_sol))
    end

    I_sol = [1; find( abs.(H*x_sol) .> 0.01ϵ); N]

    return (x_sol, I_sol, f_sol)
end

plot(g)
opt_f = []

ϵ_vec = [1e-5, 1e-6]
max_cost_vec = 10 .^ range(log10(0.1), stop=log10(8), length=2)

sols_l1_mat = []
for (i, max_cost)=enumerate(max_cost_vec)
    I = []
    Y = []
    f = []
    for (j, ϵ) = enumerate(ϵ_vec)
        _, I_new, f_new = iterative_l1_reweighting(g, max_cost, ϵ)

        if !isempty(I)
            if I != I_new
                println((i, j))
                println(I)
                println(I_new)
                println(f, ": ", f_new)
                _ = readline()
            end
        end

        I = I_new

        print(I)

        # Find optimal y values and optimal cost
        Y, f = EllZeroTrendFiltering.find_optimal_y_values(ℓ, V_N, I)

    end
    push!(sols_l1_mat, (I, Y, f))
end

I_sols_opt_vec, Y_sols_opt_vec, f_sols_opt_vec = collect(zip(sols_opt...))
f_sols_opt_vec = [f_sols_opt_vec...]

I_sols_l1_vec, Y_sols_l1_vec, f_sols_l1_vec = collect(zip(sols_l1_mat...))
f_sols_l1_vec = [f_sols_l1_vec...]

m_sols_l1_vec = [length.(I_sols_l1_vec)...]
m_sols_opt_vec = [length.(I_sols_opt_vec)...]

plot(m_sols_l1_vec, f_sols_l1_vec);
plot!(m_sols_opt_vec, f_sols_opt_vec, xlims=(0,30), ylims=(0,10));
gui()

plot(g);
#plot!(x_sol, lw=2.5)
plot!(I_vec[m], Y_vec[m], lw=2.5)

plot(H*x_sol);
plot!(I_vec[m][2:end], 0.1*diff(Y_vec[m]) ./ diff(I_vec[m]), m=:o, lw=0)

writedlm("cost_vs_m_opt.csv", [["m"; m_sols_opt_vec] ["cost"; f_sols_opt_vec]], ";")
writedlm("cost_vs_m_l1.csv", [["m"; m_sols_l1_vec] ["cost"; f_sols_l1_vec]] , ";")

writedlm(".csv", [["I"; I_sols_opt_vec[10]] ["Y"; Y_sols_opt_vec[10]]] , ";")

# Computation times for ℓ_0 trend filtering
@time fit_pwl_regularized(g, 0.2)
@time fit_pwl_regularized(g, 0.1)
@time fit_pwl_regularized(g, 0.05)
@time fit_pwl_regularized(g, 0.01)
@time fit_pwl_regularized(g, 0.2)


text:
module MobilePlatforms

# mobile robot definitions and structs; quadrotor, car, bicycle, diff drive, holonomic ...

end # (sub)module MobilePlatforms


text:
using Graphs
# Optional dependencies from Requires.jl:
using .GLMakie, .GraphMakie

function GraphMakie.graphplot(ta::TA, args...; show_guards=true, show_invariants=false, show_resets=true, show_messages=(MSG_IN,MSG_OUT), nlabels_align=(:center,:bottom), nlabels_distance=25, elabels_distance=5, kwargs...)
    g = SimpleDiGraph(length(ta.states))

    nlabels = map(ta.states) do state
        inv = show_invariants ? "\n"*repr("text/plain", ta.invariants[state]) : ""
        repr("text/plain", state)*inv
    end

    elabels=Dict{Graphs.SimpleGraphs.SimpleEdge,String}()
    foreach(ta.arcs) do arc
        id1 = findfirst(==(arc.source), ta.states)
        id2 = findfirst(==(arc.target), ta.states)
        e = Graphs.SimpleGraphs.SimpleEdge(id1,id2)
        add_edge!(g, e)

        msg = arc.message.direction in show_messages ? repr("text/plain", arc.message) : ""
        clks = show_resets ? join(String.(arc.resets) .* ":=0", ", ") : ""
        guard = show_guards ? repr("text/plain", arc.guard) : ""

        elabels[e]=join(filter(!isempty,(msg,clks,guard)), "\n")
    end

    f,ax,p= graphplot(g, args...; nlabels, elabels=getindex.(Ref(elabels),edges(g)), nlabels_align, nlabels_distance, elabels_distance, kwargs...)

    deregister_interaction!(ax, :rectanglezoom)
    register_interaction!(ax, :ndrag, NodeDrag(p))
    return f,ax,p
end

function GraphMakie.graphplot(ts::TTS, args...; kwargs...)
    g = SimpleDiGraph(length(ts.states))

    nlabels = map(ts.states) do (state, zone)
        "State: "*repr("text/plain", state)*"\nZone: "*repr("text/plain", zone)
    end

    elabels=Dict{Graphs.SimpleGraphs.SimpleEdge,String}()
    foreach(ts.transitions) do (from, (msg, to))
        id1 = findfirst(==(from), ts.states)
        id2 = findfirst(==(to), ts.states)
        e = Graphs.SimpleGraphs.SimpleEdge(id1,id2)
        add_edge!(g, e)
        elabels[e]=repr("text/plain", msg)
    end

    f,ax,p= graphplot(g, args...; nlabels, elabels=getindex.(Ref(elabels),edges(g)), kwargs...)

    deregister_interaction!(ax, :rectanglezoom)
    register_interaction!(ax, :ndrag, NodeDrag(p))
    return f,ax,p
end


text:
"""
Author: Michael Palumbo
Created: December 2019
Contact: mlp95@psu.edu
Based on code by Alex Wise (aw@psu.edu)
Refactors & optimized by Eric Ford
"""

#module FitGaussianToCCF

#import ..RVFromCCF: AbstractMeasureRvFromCCF, MeasureRvFromCCFQuadratic
# packages
using LsqFit
using LinearAlgebra

#export MeasureRvFromCCFGaussian

"""  Functor to estimate RV based on fitting a Gaussian quadratic near minimum of the CCF. """
struct MeasureRvFromCCFGaussian <: AbstractMeasureRvFromCCF
    frac_of_width_to_fit::Float64
    measure_width_at_frac_depth::Float64
    init_guess_ccf_σ::Float64
end

"""
Construct functor to estimate RV based on the CCF.
Optional Arguments:
* `frac_of_width_to_fit`: (0.5)
* `measure_width_at_frac_depth`: (0.5)
* `init_guess_ccf_σ``: (2000m/s)
"""
function MeasureRvFromCCFGaussian(; frac_of_width_to_fit::Real = default_frac_of_width_to_fit,
                                    measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth,
                                    init_guess_ccf_σ::Real = default_init_guess_ccf_σ )
    @assert 0.25 <= measure_width_at_frac_depth <= 0.75
    @assert 0.1 <= frac_of_width_to_fit <= 5.0  # TODO: FIgure out appropriate range
    @assert -30e3 <= init_guess_ccf_σ <= 30e3
    MeasureRvFromCCFGaussian(frac_of_width_to_fit,measure_width_at_frac_depth,init_guess_ccf_σ)
end

@. gaussian_line_helper(x, p) = p[4] + p[3] * exp(-0.5*((x-p[1])/p[2])^2)

function (mrv::MeasureRvFromCCFGaussian)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
        # find the min and fit only the part near the minimum of the CCF
        amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)

        # make initial guess parameters
        μ = vels[amin]
        σ = mrv.init_guess_ccf_σ                          # TODO: Make sure this is robust.
        minccf, maxccf = extrema(ccf)
        amp = minccf - maxccf
        y0 = maxccf
        p0 = [μ, σ, amp, y0]

        # fit and return the mean of the distribution
        result = curve_fit(gaussian_line_helper, view(vels,inds), view(ccf,inds), p0)

        if result.converged
           rv = coef(result)[1]
           sigma_rv = stderror(result)[1]
           rvfit = (rv=rv, σ_rv=sigma_rv)
        else
           @warn "Fit of Gaussian to CCF Failed.  Reverting to fit quadratic to CCF."
           quad_fit_to_ccf = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=mrv.frac_of_width_to_fit,measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
           rvfit = quad_fit_to_ccf(vels,ccf)
        end
        return rvfit
end

function (mrv::MeasureRvFromCCFGaussian)(vels::A1, ccf::A2, ccf_var::A3 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} }
        # find the min and fit only the part near the minimum of the CCF
        amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)

        # make initial guess parameters
        μ = vels[amin]
        σ = mrv.init_guess_ccf_σ                          # TODO: Make sure this is robust.
        minccf, maxccf = extrema(ccf)
        amp = minccf - maxccf
        y0 = maxccf
        p0 = [μ, σ, amp, y0]

        # fit and return the mean of the distribution
        result = curve_fit(gaussian_line_helper, view(vels,inds), view(ccf,inds), 1.0 ./ view(ccf_var,inds),  p0)

        if result.converged
           rv = coef(result)[1]
           sigma_rv = stderror(result)[1]
           rvfit = (rv=rv, σ_rv=sigma_rv)
        else
           @warn "Fit of Gaussian to CCF Failed.  Reverting to fit quadratic to CCF."
           # println("# Tried fitting between ", extrema(view(vels,inds)), " over which CCF varied within ", extrema(view(ccf,inds)) )
           quad_fit_to_ccf = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=mrv.frac_of_width_to_fit,measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
           rvfit = quad_fit_to_ccf(vels,ccf,ccf_var)
        end
        return rvfit
end

using PDMats
using Optim
using ForwardDiff

function (mrv::MeasureRvFromCCFGaussian)(vels::A1, ccf::A2, ccf_covar::A3 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,2} }
        # find the min and fit only the part near the minimum of the CCF
        amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)

        # make initial guess parameters
        μ = vels[amin]
        σ = mrv.init_guess_ccf_σ                          # TODO: Make sure this is robust.
        minccf, maxccf = extrema(ccf)
        amp = minccf - maxccf
        y0 = maxccf
        p0 = [μ, σ, amp, y0]

        # fit and return the mean of the distribution
        #result = curve_fit(gaussian_line_helper, view(vels,inds), view(ccf,inds), Matrix(inv(view(ccf_covar,inds,inds))),  p0, show_trace=true)
        #chol_ccf_covar = cholesky(view(ccf_covar,inds,inds))
        #lmfit_helper(p) = chol_ccf_covar * ( gaussian_line_helper(view(vels,inds), p) - view(ccf,inds) )
        #result = LsqFit.lmfit(lmfit_helper,p0,inv(chol_ccf_covar)) # ; kwargs...)

        #=
        if result.converged
           rv = coef(result)[1]
           sigma_rv = stderror(result)[1]
           rvfit = (rv=rv, σ_rv=sigma_rv)
        else
           @warn "Fit of Gaussian to CCF Failed.  Reverting to fit quadratic to CCF."
           # println("# Tried fitting between ", extrema(view(vels,inds)), " over which CCF varied within ", extrema(view(ccf,inds)) )
           quad_fit_to_ccf = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=mrv.frac_of_width_to_fit,measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
           rvfit = quad_fit_to_ccf(vels,ccf,ccf_covar)
        end
        =#

        vels_float64 = convert.(Float64,view(vels,inds))
        pd_ccf_covar = PDMat( Symmetric(view(ccf_covar,inds,inds) ) )
        #optim_helper = TwiceDifferentiable(p-> invquad(pd_ccf_covar ,  ( gaussian_line_helper(vels_float64, p) - view(ccf,inds)  ) ), p0, autodiff=:forward)
        optim_helper(p) = invquad(pd_ccf_covar ,  ( gaussian_line_helper(vels_float64, p) - view(ccf,inds)  ) )
        result = Optim.optimize(optim_helper, p0, ConjugateGradient() )
        if Optim.converged(result)
                fit_param = Optim.minimizer(result)
                rv = fit_param[1]
                numerical_hessian = ForwardDiff.hessian(optim_helper,fit_param)
                var_cov_matrix = inv(numerical_hessian)
                sigma_rv = diag(var_cov_matrix)[1]
                rvfit = (rv=rv, σ_rv=sigma_rv, covar_rv=var_cov_matrix )
        else
                @warn "Fit of Gaussian to CCF Failed.  Reverting to fit quadratic to CCF."
                println("# Tried fitting between ", extrema(view(vels,inds)), " over which CCF varied within ", extrema(view(ccf,inds)) )
                quad_fit_to_ccf = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=mrv.frac_of_width_to_fit,measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
                rvfit = quad_fit_to_ccf(vels,ccf,diag(ccf_covar))
        end

        return rvfit
end

#end # module FitGaussianToCCF


text:
struct RandomQss <: QuerySynthesisStrategy
    epsilon::Union{Float64, Vector{Float64}}
    function RandomQss(; optimizer=nothing, epsilon=0.1)
        check_epsilon(epsilon)
        new(epsilon)
    end
end

function get_query_objects(qs::RandomQss, data::Array{T, 2}, labels::Vector{Symbol}, history::Vector{Array{T, 2}})::Array{T, 2} where T <: Real
    return rand_in_hyper_rect(extrema_arrays(data[:, labels .!= :Lout])..., qs.epsilon)
end


text:
module FinancialFFT
using Requires, FFTW # for conditional dependencies
function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("deps/cuda_dependencies.jl")
    @require DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" include("deps/dual_dependencies.jl")
    @require HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" include("deps/hyper_dependencies.jl")
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" include("deps/forwarddiff_dependencies.jl")
    @require TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" include("deps/taylorseries_dependencies.jl")
end

include("methods/abstracts.jl")
include("methods/utils.jl")
include("methods/charexp.jl")
include("methods/carrmadan.jl")
include("methods/carrmadan_lewis.jl")
include("methods/lewis.jl")

export pricer

end # module


text:
module BOSSArrays

using Requests
using Blosc

export BOSSArray

function __init__()
    Blosc.set_num_threads(cld(Sys.CPU_CORES,2))
end

include("types.jl")
include("base.jl")

end # end of module


text:
function loglik_drilling_history()

    jrng = j1_indexrange(data,i)
    J = length(jrng)

    @views LLmj   = LLJ[:,1:J]
    @views gradmj = gradJ[:,:,1:J]

    fillcols!(LLmj, log.(Prj))
    zero!(gradmj)

    for j in enumerate( jrng )
        for obs in observations(data, j)
            for (m,uv) in enumerate(uvs) # TODO threaded
                LLmj[m,j] += loglik_drilling_in_t(gradmj[:,m,j])
            end
        end
    end

    logsumexp_and_softmax2!(LLm, tmp, LLj)

    gradm[k,m] = sum(gradmj[k,m,:] .* LLj[m,:])

    for obs in observations(data,i)
        for (m,uv) in enumerate(uvs)  # TODO threaded
            LLm[m] += loglik_drilling!(grad,m,uv,drillcomp)
        end # m/ψ
    end # time

end

function loglik_drilling_in_t()

    ubv = view(bigubv, drng.+1, m)

    # payoffs
    @inbounds @simd for d in drng
        ubv[d+1] = flow(u,d,obs) + β * dynamic
    end

    logPrd = ubv[d_obs+1]

    # if no gradient, return
    dograd || return logPrd - logsumexp(ubv)

    logPrd -= logsumexp_and_softmax!(ubv)

    for d in drng
        wt = d==d_obs ? 1-logPrd : -logPrd
        # theta
        for k in idx_drill_cost(model)
            grad[k,m] += wt * (flowdθ(FF, k, θt, σ, wp, s_idx, d, z, ψ, itype...) + prim.β * isev.dEV[z..., ψ, k, sp_idx, itypidx...] )
        end
        # sigma
        if s_idx <= end_ex0(wp)
            dpsi = flowdψ(FF, θt, σ, wp, s_idx, d, z, ψ, itype...) + prim.β * gradient_d(nplus1(Val{NZ}), isev.EV, z..., ψ, sp_idx, itypidx...)::T
            dsig = flowdσ(FF, θt, σ, wp, s_idx, d, z, ψ, itype...) + prim.β * isev.dEVσ[z..., ψ, sp_idx, itypidx...]
            grad[idx_drill_rho(model),m] += wt * (dpsi*_dψ1dθρ(uv..., ρ, σ) + dsig)
        end
    end # d

    return logPrd
end


text:
LinearAlgebra.det(d::Union{WignerDMatrix, WignerdMatrix}) = one(eltype(d))

function LinearAlgebra.tr(d::WignerdMatrix{T}) where {T}
	j = sphericaldegree(d)
	N = twojp1(j)
	iszero(d.beta) && return T(N)
	ωby2 = d.beta/2
	if iszero(sin(ωby2))
		# check if odd or even multiple of π
		if isone(cos(ωby2))
			# even multiple of π
			x = N
		else
			# odd multiple of π
			x = (-1)^(N+1)*N
		end
	else
		x = sin(N*ωby2)/sin(ωby2)
	end
	T(x)
end
function LinearAlgebra.tr(d::WignerdMatrix{T,<:ZeroRadians}) where {T}
	j = sphericaldegree(d)
	T(twojp1(j))
end
function LinearAlgebra.tr(d::WignerdMatrix{T, <:TwoNPiRadians}) where {T}
	j = sphericaldegree(d)
	N = twojp1(j)
	t = wignerdmatrixelement(T, j, -j, -j, d.beta)
	N*t
end
function LinearAlgebra.tr(d::WignerdMatrix{T,<:Union{Irrational{:π}, TwoNPlusOnePiRadians}}) where {T}
	j = sphericaldegree(d)
	if isinteger(j)
		res = wignerdmatrixelement(T, j, 0, 0, d.beta)
	else
		res = zero(T)
	end
	return res
end

function LinearAlgebra.tr(D::WignerDMatrix{T}) where {T}
	j = sphericaldegree(D)
	α,β,γ = eulerangles(D)
	cosωby2 = cos(β/2)cos((α+γ)/2)
	if isone(cosωby2)
		return T(twojp1(j))
	elseif isone(-cosωby2)
		N = twojp1(j)
		return (-1)^(N+1) * N
	end
	ωby2 = acos(cosωby2)
	x = sin((2j+1)ωby2)/sin(ωby2)
	T(x)
end
function LinearAlgebra.tr(D::WignerDMatrix{T, <:Any,
	WignerdMatrix{<:Any, <:Union{Irrational{:π}, TwoNPlusOnePiRadians}}}) where {T}

	j = sphericaldegree(D)
	if isinteger(j)
		res = Complex(tr(parent(D)))
	else
		res = zero(T)
	end
	return res
end

function LinearAlgebra.inv(D::WignerDMatrix)
	j = sphericaldegree(D)
	α, β, γ = eulerangles(D)
	WignerDMatrix(j, -γ, -β, -α)
end

function LinearAlgebra.inv(D::WignerdMatrix)
	j = sphericaldegree(D)
	_, β, _ = eulerangles(D)
	WignerdMatrix(j, -β)
end

function Base.:(*)(d1::WignerdMatrix, d2::WignerdMatrix)
	sphericaldegree(d1) == sphericaldegree(d2) ||
		throw(ArgumentError("d-matrices with j = $(sphericaldegree(d1)) "*
			"and j = $(sphericaldegree(d2)) can not be multiplied"))
	j = sphericaldegree(d1)

	WignerdMatrix(j, d1.beta + d2.beta)
end

function Base.:(*)(D1::WignerDMatrix, D2::WignerDMatrix)
	SpinMatrix(collect(D1) * collect(D2))
end

function Base.:(*)(D1::WignerDMatrix, d2::WignerdMatrix)
	D1 * WignerDMatrix(0, d2, 0)
end
function Base.:(*)(d1::WignerdMatrix, D2::WignerDMatrix)
	WignerDMatrix(0, d1, 0) * D2
end

for f in [:issymmetric, :ishermitian, :isdiag, :isbanded]
	@eval LinearAlgebra.$f(d::AbstractWignerMatrix) = LinearAlgebra.$f(parent(d))
end

function LinearAlgebra.diag(d::AbstractWignerMatrix)
	N = size(d,1)
	v = Vector{eltype(d)}(undef, N)
	j = sphericaldegree(d)
	@inbounds for (i,m) in enumerate(-j:j)
		v[i] = d[m,m]
	end
	return v
end

function LinearAlgebra.diag(d::WignerdMatrix{T,<:ZeroRadians}) where {T}
	ones(T, size(d,1))
end

function LinearAlgebra.diag(d::WignerdMatrix{T, <:TwoNPiRadians}) where {T}
	j = sphericaldegree(d)
	t = wignerdmatrixelement(T, j, -j, -j, d.beta)
	fill(t, size(d,1))
end

function LinearAlgebra.diag(d::WignerdMatrix{T,
	<:Union{Irrational{:π}, TwoNPlusOnePiRadians}}) where {T}

	N = size(d,1)
	v = zeros(T, N)
	j = sphericaldegree(d)
	if isodd(N) # integer j
		v[div(N,2) + 1] = wignerdmatrixelement(T, j, 0, 0, d.beta)
	end
	return v
end

function LinearAlgebra.diag(d::WignerDMatrix{T, <:Any,
	<:WignerdMatrix{<:Any, <:Union{Irrational{:π}, TwoNPlusOnePiRadians}}}) where {T}

	N = size(d,1)
	v = zeros(T, N)
	j = sphericaldegree(d)
	_,β,_ = eulerangles(d)
	if isodd(N) # integer j
		v[div(N,2) + 1] = wignerdmatrixelement(typeofreal(T), j, 0, 0, β)
	end
	return v
end