text:
#######################
## Logistic Circuits
#######################

abstract type LogisticΔNode{O} <: DecoratorΔNode{O} end
abstract type LogisticLeafNode{O} <: LogisticΔNode{O} end
abstract type LogisticInnerNode{O} <: LogisticΔNode{O} end

struct LogisticLiteral{O} <: LogisticLeafNode{O}
    origin::O
end

struct Logistic⋀{O} <: LogisticInnerNode{O}
    origin::O
    children::Vector{<:LogisticΔNode{<:O}}
end

mutable struct Logistic⋁{O} <: LogisticInnerNode{O}
    origin::O
    children::Vector{<:LogisticΔNode{<:O}}
    thetas::Array{Float64, 2}
end

const LogisticΔ{O} = AbstractVector{<:LogisticΔNode{O}}

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension

@inline GateType(::Type{<:LogisticLiteral}) = LiteralGate()
@inline GateType(::Type{<:Logistic⋀}) = ⋀Gate()
@inline GateType(::Type{<:Logistic⋁}) = ⋁Gate()

#####################
# constructors and conversions
#####################

function Logistic⋁(::Type{O}, origin, children, classes::Int) where {O}
    Logistic⋁{O}(origin, children, Array{Float64, 2}(undef, (length(children), classes)))
end

const LogisticCache = Dict{ΔNode, LogisticΔNode}

function LogisticΔ(circuit::Δ, classes::Int, cache::LogisticCache = LogisticCache())

    sizehint!(cache, length(circuit)*4÷3)

    O = grapheltype(circuit) # type of node in the origin

    pc_node(::LiteralGate, n::ΔNode) = LogisticLiteral{O}(n)
    pc_node(::ConstantGate, n::ΔNode) = error("Cannot construct a logistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")

    pc_node(::⋀Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Logistic⋀{O}(n, children)
    end

    pc_node(::⋁Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Logistic⋁(O, n, children, classes)
    end

    map(circuit) do node
        pcn = pc_node(GateType(node), node)
        cache[node] = pcn
        pcn
    end
end

#####################
# methods
#####################

import LogicCircuits: literal, children # make available for extension

@inline literal(n::LogisticLiteral)::Lit  = literal(n.origin)
@inline children(n::LogisticInnerNode) = n.children
@inline classes(n::Logistic⋁) = size(n.thetas)[2]

num_parameters(n::Logistic⋁) = num_children(n) * classes(n)
num_parameters(c::LogisticΔ) = sum(n -> num_parameters(n), ⋁_nodes(c))

num_parameters_perclass(n::Logistic⋁) = num_children(n)
num_parameters_perclass(c::LogisticΔ) = sum(n -> num_parameters_perclass(n), ⋁_nodes(c))

"Return the first origin that is a Logistic circuit node"
logistic_origin(n::DecoratorΔNode)::LogisticΔNode = origin(n,LogisticΔNode)

"Return the first origin that is a Logistic circuit"
logistic_origin(c::DecoratorΔ)::LogisticΔ = origin(c, LogisticΔNode)

# TODO Learning

# Class Conditional Probability
function class_conditional_likelihood_per_instance(fc::FlowΔ,
                                                    classes::Int,
                                                    batch::PlainXData{Bool})
    lc = origin(origin(fc))
    @assert(lc isa LogisticΔ)
    pass_up_down(fc, batch)
    likelihoods = zeros(num_examples(batch), classes)
    for n in fc
        orig = logistic_origin(n)
        if orig isa Logistic⋁
            # For each class. orig.thetas is 2D so used eachcol
            for (idx, thetaC) in enumerate(eachcol(orig.thetas))
                foreach(n.children, thetaC) do c, theta
                    likelihoods[:, idx] .+= prod_fast(downflow(n), pr_factors(origin(c))) .* theta
                end
            end
        end
    end
    likelihoods
end

"""
Calculate conditional log likelihood for a batch of samples with evidence P(c | x).
(Also returns the generated FlowΔ)
"""
function class_conditional_likelihood_per_instance(lc::LogisticΔ,
                                                    classes::Int,
                                                    batch::PlainXData{Bool})
    opts = (max_factors = 2, compact⋀=false, compact⋁=false)
    fc = FlowΔ(lc, num_examples(batch), Float64, opts)
    (fc, class_conditional_likelihood_per_instance(fc, classes, batch))
end


text:
using OrbitalTrajectories

using Test

@testset "SpiceUtils" begin
    @test OrbitalTrajectories.SpiceUtils.mass_fraction("jupiter", "europa") ≈ 2.528e-5 atol=1e-5
    @test OrbitalTrajectories.SpiceUtils.mass_fraction("mars", "phobos") ≈ 1.0659e16 / (1.0659e16 + 6.4171e23) atol=1e-5
end

ext: jl


text:
module Wut

include("BitPat.jl")

abstract type AbstractEncoder
end

# implemented by all concrete encoders:
#   getWidth(a::ConcreteEncoder)
#   getDescription(a::ConcreteEncoder)
#   encodeIntoArray(a::ConcreteEncoder, inputData, output)

#  function setLearning(a::AbstractEncoder, learningEnabled):
#  function setFieldStats(a::AbstractEncoder, fieldName, fieldStatistics):
#  function getScalarNames(a::AbstractEncoder, parentFieldName=''):
#  function getDecoderOutputFieldTypes(a::AbstractEncoder):
#  function setStateLock(a::AbstractEncoder,lock):
#  function _getInputValue(a::AbstractEncoder, obj, fieldName):
#  function getEncoderList(a::AbstractEncoder):
#  function getScalars(a::AbstractEncoder, inputData):
#  function getEncodedValues(a::AbstractEncoder, inputData):
#  function getBucketIndices(a::AbstractEncoder, inputData):
#  function scalarsToStr(a::AbstractEncoder, scalarValues, scalarNames=None):
#  function getFieldDescription(a::AbstractEncoder, fieldName):
#  function encodedBitDescription(a::AbstractEncoder, bitOffset, formatted=False):
#  function pprintHeader(a::AbstractEncoder, prefix=""):
#  function pprint(a::AbstractEncoder, output, prefix=""):
#  function decode(a::AbstractEncoder, encoded, parentFieldName=''):
#  function decodedToStr(a::AbstractEncoder, decodeResults):
#  function getBucketValues(a::AbstractEncoder):
#  function getBucketInfo(a::AbstractEncoder, buckets):
#  function topDownCompute(a::AbstractEncoder, encoded):
#  function closenessScores(a::AbstractEncoder, expValues, actValues, fractional=True):

function encode(a::AbstractEncoder, v)
    encodeIntoArray(a, v, BitPat(getWidth(a)))
end
export encode

function getDisplayWidth(a::AbstractEncoder)
     getWidth(a) + length(getDescription(a)) - 1
end
export getDisplayWidth

function decode(e::AbstractEncoder, b::BitPat; parentFieldName="")

    parentName = (parentFieldName == "") ? getName(e) : @sprintf("%s.%s", parentFieldName, getName(e))
    x=Dict(parentName => [])
    ("abstract",x,[parentName])
end
export decode

include("ScalarEncoder.jl")

#include("AdaptiveScalarEncoder.jl")
include("CategoryEncoder.jl")
#include("CoordinateEncoder.jl")
include("DateEncoder.jl")
#include("DeltaEncoder.jl")
#include("GeospatialCoordinateEncoder.jl")
#include("LogEncoder.jl")
#include("MultiEncoder.jl")
#include("PassThroughEncoder.jl")
include("RandomDistributedScalarEncoder.jl")
#include("SDRCategoryEncoder.jl")
#include("ScalarSpaceEncoder.jl")
#include("SparsePassThroughEncoder.jl")

end # module


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct SingleCageInductionMachine <: DynamicInjection
        name::String
        R_s::Float64
        R_r::Float64
        X_ls::Float64
        X_lr::Float64
        X_m::Float64
        H::Float64
        A::Float64
        B::Float64
        base_power::Float64
        ext::Dict{String, Any}
        C::Float64
        τ_ref::Float64
        B_shunt::Float64
        X_ad::Float64
        X_aq::Float64
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of 5-states three-phase single cage induction machine with quadratic torque-speed relationship.

# Arguments
- `name::String`
- `R_s::Float64`: Armature stator resistance, validation range: `(0, nothing)`
- `R_r::Float64`: Rotor resistance, validation range: `(0, nothing)`
- `X_ls::Float64`: Stator Leakage Reactance, validation range: `(0, nothing)`
- `X_lr::Float64`: Rotor Leakage Reactance, validation range: `(0, nothing)`
- `X_m::Float64`: Stator-Rotor Mutual Reactance, validation range: `(0, nothing)`
- `H::Float64`: Motor Inertia Constant [s], validation range: `(0, nothing)`
- `A::Float64`: Torque-Speed Quadratic Term, validation range: `(0, 1)`
- `B::Float64`: Torque-Speed Linear Term, validation range: `(0, 1)`
- `base_power::Float64`: Base power, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `C::Float64`: Torque-Speed Constant Term
- `τ_ref::Float64`: Reference torque parameter
- `B_shunt::Float64`: Susceptance Initialization Corrector Term
- `X_ad::Float64`: Equivalent d-axis reactance
- `X_aq::Float64`: Equivalent q-axis reactance
- `states::Vector{Symbol}`: The states are:
	ψ_qs: stator flux in the q-axis,
	ψ_ds: stator flux in the d-axis,
	ψ_qr: rotor flux in the q-axis,
	ψ_dr: rotor flux in the d-axis,
	ωr: Rotor speed [pu],
- `n_states::Int`: SingleCageInductionMachine has 5 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SingleCageInductionMachine <: DynamicInjection
    name::String
    "Armature stator resistance"
    R_s::Float64
    "Rotor resistance"
    R_r::Float64
    "Stator Leakage Reactance"
    X_ls::Float64
    "Rotor Leakage Reactance"
    X_lr::Float64
    "Stator-Rotor Mutual Reactance"
    X_m::Float64
    "Motor Inertia Constant [s]"
    H::Float64
    "Torque-Speed Quadratic Term"
    A::Float64
    "Torque-Speed Linear Term"
    B::Float64
    "Base power"
    base_power::Float64
    ext::Dict{String, Any}
    "Torque-Speed Constant Term"
    C::Float64
    "Reference torque parameter"
    τ_ref::Float64
    "Susceptance Initialization Corrector Term"
    B_shunt::Float64
    "Equivalent d-axis reactance"
    X_ad::Float64
    "Equivalent q-axis reactance"
    X_aq::Float64
    "The states are:
	ψ_qs: stator flux in the q-axis,
	ψ_ds: stator flux in the d-axis,
	ψ_qr: rotor flux in the q-axis,
	ψ_dr: rotor flux in the d-axis,
	ωr: Rotor speed [pu],"
    states::Vector{Symbol}
    "SingleCageInductionMachine has 5 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SingleCageInductionMachine(name, R_s, R_r, X_ls, X_lr, X_m, H, A, B, base_power, ext=Dict{String, Any}(), )
    SingleCageInductionMachine(name, R_s, R_r, X_ls, X_lr, X_m, H, A, B, base_power, ext, PowerSystems.calculate_IM_torque_params(A, B), 1.0, 0.0, (1.0 / X_m + 1.0 / X_ls + 1.0 / X_lr)^(-1), X_ad, [:ψ_qs, :ψ_ds, :ψ_qr, :ψ_dr, :ωr], 5, InfrastructureSystemsInternal(), )
end

function SingleCageInductionMachine(; name, R_s, R_r, X_ls, X_lr, X_m, H, A, B, base_power, ext=Dict{String, Any}(), C=PowerSystems.calculate_IM_torque_params(A, B), τ_ref=1.0, B_shunt=0.0, X_ad=(1.0 / X_m + 1.0 / X_ls + 1.0 / X_lr)^(-1), X_aq=X_ad, states=[:ψ_qs, :ψ_ds, :ψ_qr, :ψ_dr, :ωr], n_states=5, internal=InfrastructureSystemsInternal(), )
    SingleCageInductionMachine(name, R_s, R_r, X_ls, X_lr, X_m, H, A, B, base_power, ext, C, τ_ref, B_shunt, X_ad, X_aq, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function SingleCageInductionMachine(::Nothing)
    SingleCageInductionMachine(;
        name="init",
        R_s=0,
        R_r=0,
        X_ls=0,
        X_lr=0,
        X_m=0,
        H=0,
        A=0.0,
        B=0.0,
        base_power=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`SingleCageInductionMachine`](@ref) `name`."""
get_name(value::SingleCageInductionMachine) = value.name
"""Get [`SingleCageInductionMachine`](@ref) `R_s`."""
get_R_s(value::SingleCageInductionMachine) = value.R_s
"""Get [`SingleCageInductionMachine`](@ref) `R_r`."""
get_R_r(value::SingleCageInductionMachine) = value.R_r
"""Get [`SingleCageInductionMachine`](@ref) `X_ls`."""
get_X_ls(value::SingleCageInductionMachine) = value.X_ls
"""Get [`SingleCageInductionMachine`](@ref) `X_lr`."""
get_X_lr(value::SingleCageInductionMachine) = value.X_lr
"""Get [`SingleCageInductionMachine`](@ref) `X_m`."""
get_X_m(value::SingleCageInductionMachine) = value.X_m
"""Get [`SingleCageInductionMachine`](@ref) `H`."""
get_H(value::SingleCageInductionMachine) = value.H
"""Get [`SingleCageInductionMachine`](@ref) `A`."""
get_A(value::SingleCageInductionMachine) = value.A
"""Get [`SingleCageInductionMachine`](@ref) `B`."""
get_B(value::SingleCageInductionMachine) = value.B
"""Get [`SingleCageInductionMachine`](@ref) `base_power`."""
get_base_power(value::SingleCageInductionMachine) = value.base_power
"""Get [`SingleCageInductionMachine`](@ref) `ext`."""
get_ext(value::SingleCageInductionMachine) = value.ext
"""Get [`SingleCageInductionMachine`](@ref) `C`."""
get_C(value::SingleCageInductionMachine) = value.C
"""Get [`SingleCageInductionMachine`](@ref) `τ_ref`."""
get_τ_ref(value::SingleCageInductionMachine) = value.τ_ref
"""Get [`SingleCageInductionMachine`](@ref) `B_shunt`."""
get_B_shunt(value::SingleCageInductionMachine) = value.B_shunt
"""Get [`SingleCageInductionMachine`](@ref) `X_ad`."""
get_X_ad(value::SingleCageInductionMachine) = value.X_ad
"""Get [`SingleCageInductionMachine`](@ref) `X_aq`."""
get_X_aq(value::SingleCageInductionMachine) = value.X_aq
"""Get [`SingleCageInductionMachine`](@ref) `states`."""
get_states(value::SingleCageInductionMachine) = value.states
"""Get [`SingleCageInductionMachine`](@ref) `n_states`."""
get_n_states(value::SingleCageInductionMachine) = value.n_states
"""Get [`SingleCageInductionMachine`](@ref) `internal`."""
get_internal(value::SingleCageInductionMachine) = value.internal

"""Set [`SingleCageInductionMachine`](@ref) `R_s`."""
set_R_s!(value::SingleCageInductionMachine, val) = value.R_s = val
"""Set [`SingleCageInductionMachine`](@ref) `R_r`."""
set_R_r!(value::SingleCageInductionMachine, val) = value.R_r = val
"""Set [`SingleCageInductionMachine`](@ref) `X_ls`."""
set_X_ls!(value::SingleCageInductionMachine, val) = value.X_ls = val
"""Set [`SingleCageInductionMachine`](@ref) `X_lr`."""
set_X_lr!(value::SingleCageInductionMachine, val) = value.X_lr = val
"""Set [`SingleCageInductionMachine`](@ref) `X_m`."""
set_X_m!(value::SingleCageInductionMachine, val) = value.X_m = val
"""Set [`SingleCageInductionMachine`](@ref) `H`."""
set_H!(value::SingleCageInductionMachine, val) = value.H = val
"""Set [`SingleCageInductionMachine`](@ref) `A`."""
set_A!(value::SingleCageInductionMachine, val) = value.A = val
"""Set [`SingleCageInductionMachine`](@ref) `B`."""
set_B!(value::SingleCageInductionMachine, val) = value.B = val
"""Set [`SingleCageInductionMachine`](@ref) `base_power`."""
set_base_power!(value::SingleCageInductionMachine, val) = value.base_power = val
"""Set [`SingleCageInductionMachine`](@ref) `ext`."""
set_ext!(value::SingleCageInductionMachine, val) = value.ext = val
"""Set [`SingleCageInductionMachine`](@ref) `C`."""
set_C!(value::SingleCageInductionMachine, val) = value.C = val
"""Set [`SingleCageInductionMachine`](@ref) `τ_ref`."""
set_τ_ref!(value::SingleCageInductionMachine, val) = value.τ_ref = val
"""Set [`SingleCageInductionMachine`](@ref) `B_shunt`."""
set_B_shunt!(value::SingleCageInductionMachine, val) = value.B_shunt = val
"""Set [`SingleCageInductionMachine`](@ref) `X_ad`."""
set_X_ad!(value::SingleCageInductionMachine, val) = value.X_ad = val
"""Set [`SingleCageInductionMachine`](@ref) `X_aq`."""
set_X_aq!(value::SingleCageInductionMachine, val) = value.X_aq = val


text:
include("A2C_utils.jl")
using Flux, Statistics, Random, Distributions, StatsFuns, Zygote

T = 100
Nhidden = 50
Nstates = 1
Nstate_rep = 1
Naction = 2
Nin = Naction+1+1 #no state information for this agent, just actions, reward, time
Nout = Naction+1
βv = 0.05f0
βe = 0.05f0
γ = 0.9f0

function forward(m, x)
    #this just adds a softmax to the policy output
    y = m(x)
    logπ = y[1:2]
    logπ = logπ .- StatsFuns.logsumexp(logπ) #softmax
    return [logπ; y[3]]
end

function env(y, t, state, ps)
    ### input the agent output; output action, reward, and new input
    πt = exp.(y[1:2]) #probability of pulling arms
    d = Binomial(1, Float64(πt[2])) #probability of second arm
    a = (rand(d)+1) # 0 -> 1, 1 -> 2
    rew = rand(Binomial(1, Float64(ps[a]))) #draw reward from corresponding arm
    ahot = zeros(2)
    ahot[a] = 1.
    x = [ahot; rew; t] #input is action, reward, timestep
    ps = ps #reward probabilities do not change
    s_new = zeros(1) #state is empty
    return Int(a), Float32(rew), Float32.(x), Float32.(ps), Float32.(s_new)
end
Zygote.@nograd env #don't take gradients of this

function initialize(ps, state)
    if maximum(ps) == 0 #if not provided
        p = rand()*0.95+0.025
        ps = [p; 1-p]
    end
    state = Float32.(zeros(1)) #no state
    return Float32.(ps), Float32.(state)
end
Zygote.@nograd initialize #don't take gradients of this

function model_eval()
    #evaluation function; compute mean reward (random is 0.5; oracle is 0.745)
    means = []
    all_as = []
    all_ps = [[i/50.; 1-i/50.] for i = 1:49]
    for ps = all_ps
        L, ys, rews, as, ps = run_episode(m, ps = ps, hidden = false)
        means = [means; mean(rews)]
        all_as = [all_as; mean(as)]
    end
    return mean(means), mean(all_as)
end

max_stars_repo_name: KrisJensen/metaRL

max_issues_repo_name: KrisJensen/metaRL

max_forks_repo_name: KrisJensen/metaRL


text:
# a simple LP example
#
#   maximize 3x + 2y
#
#   s.t. 1 x + 1 y == 20
#        x \in [0, 10], y \in [0, 5]
#

using Xpress, Compat.Test

# @testset "Basics 1" begin
    model = Xpress.Model("lp_01", :minimize)

    # add variables
    add_cvar!(model, 3.0, 0., 10.)  # x
    add_cvar!(model, 2.0,  0., 5.)  # y

    # add constraints
    add_constr!(model, [1., 1.], '=', 20.)

    @test Xpress.get_obj(model) == [3, 2]
    @test Xpress.get_constrmatrix(model) == sparse([1 1])
    @test Xpress.get_rhs(model) == [20]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [10, 5]
    @test Xpress.get_sense(model) == [:(==)]

    # perform optimization
    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :infeasible

    lrp = Float64[0.0]
    grp = Float64[0.1]
    lbp = Float64[0.0,0.0]
    ubp = Float64[0.0,0.0]

    @test 0 == Xpress.repairweightedinfeasibility(model, lrp, grp, lbp, ubp, phase2 = Cchar('d'))
    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

@show Xpress.get_objval(model)
@show Xpress.get_solution(model)
@show Xpress.get_dual(model)

    @test 0 == Xpress.repairweightedinfeasibility(model, lrp, grp, lbp, ubp, phase2 = Cchar('x'))
    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

@show Xpress.get_objval(model)
@show Xpress.get_solution(model)
@show Xpress.get_dual(model)

    @test 0 == Xpress.repairweightedinfeasibility(model, lrp, grp, lbp, ubp, phase2 = Cchar('f'))
    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

@show Xpress.get_objval(model)
@show Xpress.get_solution(model)
@show Xpress.get_dual(model)

    @test 0 == Xpress.repairweightedinfeasibility(model, lrp, grp, lbp, ubp, phase2 = Cchar('n'))
    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

@show Xpress.get_objval(model)
@show Xpress.get_solution(model)
@show Xpress.get_dual(model)
;


text:
"""
Nasdaq

This provides instance types to represent the stocks listed on the Nasdaq stock exchange

See README.md for the full documentation

Copyright 2020, Eric Forgy, Scott P. Jones and other contributors

Licensed under MIT License, see LICENSE.md
"""
module Nasdaq

using Assets, HTTP, JSON3

const _nasdaq_data = Dict{Symbol,Tuple{ListedEquity,String}}()

allsymbols() = keys(_nasdaq_data)
allpairs() = pairs(_nasdaq_data)

#__init__() = load()

function load(; host="localhost", port=8000, force=false, name="Nasdaq", symbols=[], currency=:USD)
    exch = lowercase(name)
    jsondir = joinpath(@__DIR__, "..", "data")
    jsonfile = joinpath(jsondir, string(exch,"-",join(symbols,"-"),".json"))
    if !force && isfile(jsonfile)
        body = read(jsonfile)
    else
        body = HTTP.get("http://$host:$port/api/$exch", query=Dict(:symbols => join(symbols,","))).body
        ispath(jsondir) || mkpath(jsondir)
        open(io -> write(io, body), jsonfile, "w")
    end
    response = JSON3.read(body, Dict{Symbol,String})
    exchsym = Symbol(uppercase(name))
    stocks = ListedEquity[]
    for key in keys(response)
        symbol = Symbol(key)
        push!(stocks, ListedEquity(exchsym, symbol, currency))
        _nasdaq_data[symbol] = (ListedEquity(exchsym, symbol, currency), response[symbol])
    end
    stocks
end

end # module Nasdaq


text:
@inline _linear_endpoints(z::Real) = (-z, z)
@inline _linear_endpoints(z::AbstractInterval) = endpoints(z)

@inline _radial_endpoints(r::T) where {T <: Real} = (zero(T), r)
@inline _radial_endpoints(r::AbstractInterval) = endpoints(r)

@inline _in_angular_interval_closed(α::Real, α_int::Nothing, tol::Real = 0) = true
@inline _in_angular_interval_closed(α::Real, α_int::AbstractInterval{T}, tol::Real = 0) where {T} = mod(α - (α_int.left-tol), T(2π)) ≤ (α_int.right+tol) - (α_int.left-tol)

@inline function _in_angular_interval_closed(α::T, α_int::Tuple{T,T}; csgtol::T = csg_default_tol(T)) where {T}
    m = mod(α - α_int[1], T(2π))
    d = (α_int[2] - α_int[1])
    m ≤ d + csgtol
end
@inline function _in_angular_interval_open(α::T, α_int::Tuple{T,T}; csgtol::T = csg_default_tol(T)) where {T}
    csgtol < mod(α - α_int[1], T(2π)) < (α_int[2] - α_int[1]) - csgtol
end

@inline _in_angular_interval_open(α::T, α_int::Nothing) where {T<:Real} = 0 < mod(α, T(2π)) < T(2π)
@inline _in_angular_interval_open(α::Real, α_int::AbstractInterval{T}) where {T} = 0 < mod(α - α_int.left, T(2π)) < width(α_int)


text:
module TreeNodes

Base.Experimental.@optlevel 1

using Redux
using CImGui

# helpers
"""
    TreeNode(f::Function, label::AbstractString) -> Bool
Return `true` when the node is open.
"""
function TreeNode(f::Function, label::AbstractString)
    is_open = CImGui.TreeNode(label)
    if is_open
        f()
        CImGui.TreePop()
    end
    return is_open
end

"""
    TreeNodeEx(f::Function, label::AbstractString, flags=0) -> Bool
Return `true` when the node is open.
"""
function TreeNodeEx(f::Function, label::AbstractString, flags=0)
    is_open = CImGui.TreeNodeEx(label, flags)
    if is_open
        f()
        CImGui.TreePop()
    end
    return is_open
end

end # module


text:
macro test_binary_operator(op, vals)
    quote
        @test $(op)(false, false) == $(vals)[1]
        @test $(op)(true, false) == $(vals)[2]
        @test $(op)(false, true) == $(vals)[3]
        @test $(op)(true, true) == $(vals)[4]
    end
end

@testset "Expressions" begin
    @testset "Built-in Operators" begin
        @test isdefined(Deductive, :¬)
        @test isdefined(Deductive, :∧)
        @test isdefined(Deductive, :∨)
        @test isdefined(Deductive, :→)
        @test isdefined(Deductive, :⟷)

        @test ¬(true) == false
        @test ¬(false) == true

        @test_binary_operator ∨ (false, true, true, true)
        @test_binary_operator ∧ (false, false, false, true)

        @test_binary_operator → (true, false, true, true)
        @test_binary_operator ⟷ (true, false, false, true)
    end

    @testset "LogicalSymbol Core" begin
        @test (LogicalSymbol <: Deductive.AbstractExpression)

        a = LogicalSymbol(:a)

        my_metadata = "AbCdEfG"
        with_metadata = LogicalSymbol(:b, my_metadata)

        # conversion to julia symbol
        @test Symbol(a) == :a
        @test Symbol(with_metadata) == :b

        # tree / node tools
        @test !istree(a)
        @test isnode(a)
        @test !istree(with_metadata)
        @test isnode(with_metadata)

        # metadata getters
        @test isnothing(metadata(a))
        @test !isnothing(metadata(with_metadata))
        @test metadata(with_metadata) == my_metadata

        # constituent variables
        @test length(variables(a)) == 1
        @test first(variables(a)) === a
        @test length(variables(with_metadata)) == 1
        @test first(variables(with_metadata)) == with_metadata

        # constituent symbols (always none)
        @test isempty(operations(a))
        @test isempty(operations(with_metadata))

        # symbol comparison
        @test isequal(a, a)
        @test isequal(with_metadata, with_metadata)
        @test isequal(a, LogicalSymbol(:a))
        @test isequal(with_metadata, LogicalSymbol(:b, my_metadata))
        @test !isequal(a, LogicalSymbol(:b))
        @test !isequal(with_metadata, LogicalSymbol(:b))
        @test !isequal(a, LogicalSymbol(:a, 0))
        @test !isequal(a, LogicalSymbol(:a, (nothing,)))
    end

    @testset "LogicalOperation Core" begin
        neg_op = Deductive.LogicalOperation(x -> !x, :¬, 1, false, false)
        and_op = Deductive.LogicalOperation(:∧, 2, true, true) do x, y
            x && y
        end
        ternary_op = Deductive.LogicalOperation(:Δ, 3, true, true) do x, y, z
            x && y && z
        end

        @testset "Operation Names" begin
            @test Symbol(neg_op) == :¬
            @test Symbol(and_op) == :∧
            @test Symbol(ternary_op) == :Δ
        end

        @testset "isunary and isbinary" begin
            @test isunary(neg_op)
            @test !isbinary(neg_op)
            @test !isunary(and_op)
            @test isbinary(and_op)
            @test !isunary(ternary_op)
            @test !isbinary(ternary_op)
        end

        @testset "Equality Checking" begin
            # two operators are equal if their names, arguments, and associative / commutative properties are the same
            @test isequal(neg_op, Deductive.LogicalOperation((x) -> !x, :¬, 1, false, false))
            @test !isequal(neg_op, Deductive.LogicalOperation((x) -> !x, :¬, 2, false, false))
            @test !isequal(neg_op, Deductive.LogicalOperation((x) -> !x, :-, 1, false, false))
            @test isequal(neg_op, Deductive.LogicalOperation((x) -> x, :¬, 1, false, false))
            @test !isequal(neg_op, Deductive.LogicalOperation((x) -> x, :¬, 1, true, false))
            @test !isequal(neg_op, Deductive.LogicalOperation((x) -> x, :¬, 1, false, true))
        end

        @testset "Material Operation" begin
            @test neg_op(false)
            @test !neg_op(true)
            @test and_op(true, true)
            @test !and_op(true, false)
            @test ternary_op(true, true, true)
            @test !ternary_op(false, true, true)
            @test !ternary_op(true, false, true)
            @test !ternary_op(true, true, false)
            @test ternary_op(BitVector((1, 1, 1)))
            @test !ternary_op(BitVector((1, 0, 1)))
        end

        @testset "isassociative and iscommutative" begin
            # isassociative for builtins
            @test isassociative(∧)
            @test isassociative(∨)
            @test isassociative(⟷)
            @test !isassociative(¬)
            @test !isassociative(→)

            # iscommutative for builtins
            @test iscommutative(∧)
            @test iscommutative(∨)
            @test iscommutative(⟷)
            @test !iscommutative(¬)
            @test !iscommutative(→)
        end
    end

    @testset "LogicalExpression Core" begin
        @test (Deductive.LogicalExpression <: Deductive.AbstractExpression)

        a, b, c, d = LogicalSymbol.([:a, :b, :c, :d])
        st = a ∧ b

        @testset "Basics" begin
            @test st isa Deductive.LogicalExpression
            @test istree(st)
            @test !isnode(st)
            @test operation(st) === ∧
            @test all(arguments(st) .== (a, b))
            @test isnothing(metadata(st))
        end

        @testset "Equality Checking" begin
            @test isequal(st, st)
            @test isequal(st, a ∧ b)
            @test isequal(a ∧ b, a ∧ b)
            @test !isequal(a ∧ b, a ∨ b)
            @test !isequal(¬a, ¬¬a)
            @test !isequal((a ∧ b) ∧ c, a ∧ (b ∧ c))
            @test !isequal(a ∧ b, b ∧ a)
            @test !isequal(a ∧ b, a ∧ c)
        end

        @testset "Variables Tree" begin
            @test variables(¬a) == Set([a])
            @test variables(¬¬a) == Set([a])
            @test variables(a ∧ b) == Set([a, b])
            @test variables(¬(a ∨ b)) == Set([a, b])
            @test variables(¬(a ∨ b) → a) == Set([a, b])
            @test variables(¬(a ∨ b) → c) == Set([a, b, c])
        end

        @testset "Operations Tree" begin
            @test operations(¬a) == Set([¬])
            @test operations(¬¬a) == Set([¬])
            @test operations(a ∧ b) == Set([∧])
            @test operations(a ∧ ¬b) == Set([∧, ¬])
            @test operations(a ∧ ¬(b ∨ c)) == Set([∧, ∨, ¬])
            @test operations(¬(a ∨ b) → a) == Set([∨, ¬, →])
        end

        @testset "Binary Operator Methods" begin
            @test left(a ∧ b) == a
            @test right(a ∧ b) == b
            @test left((a ∧ b) ∨ (b ∧ c)) == (a ∧ b)
            @test right((a ∧ b) ∨ (b ∧ c)) == (b ∧ c)
        end

        @testset "Associative Expressions" begin
            @test isassociative(a ∧ b)
            @test isassociative(a ∨ b)
            @test isassociative(a ∨ b ∨ c)
            @test isassociative(a ⟷ b ⟷ c)
            @test !isassociative(¬a)
            @test !isassociative(a → b)
            @test !isassociative(a → b → c)
        end

        @testset "Commutative Expressions" begin
            @test iscommutative(a ∧ b)
            @test iscommutative(a ∨ b)
            @test iscommutative(a ∨ b ∨ c)
            @test iscommutative(a ⟷ b ⟷ c)
            @test !iscommutative(¬a)
            @test !iscommutative(a → b)
            @test !iscommutative(a → b → c)
        end

        @testset "Associative Orderings" begin
            @test associative_ordering(a) == [a]
            @test associative_ordering(b) == [b]
            @test associative_ordering(a ∧ b) == [a, b]
            @test associative_ordering(b ∧ a) == [b, a]
            @test associative_ordering(¬b ∧ a) == [b, a]
            @test associative_ordering((¬b ∨ c) ∧ a) == [b, c, a]
            @test associative_ordering((¬a ∨ (b ∨ c)) ∧ d) == [a, b, c, d]
        end

        @testset "Associative Equality" begin
            @test isequal_associative(a, a)
            @test isequal_associative(b, b)
            @test !isequal_associative(a, b)
            @test isequal_associative(a ∧ b, a ∧ b)
            @test !isequal_associative(a ∧ b, b ∧ a)
            @test !isequal_associative(a ∧ b, a ∨ b)
            @test !isequal_associative(a ∧ ¬b, a ∧ ¬b)
        end

        @testset "Associative Tree Count" begin
            @test Deductive.associative_tree_count(0) == 1
            @test Deductive.associative_tree_count(1) == 1
            @test Deductive.associative_tree_count(2) == 2
            @test Deductive.associative_tree_count(3) == 5
            @test Deductive.associative_tree_count(4) == 14
            @test Deductive.associative_tree_count(5) == 42
            @test Deductive.associative_tree_count(6) == 132
        end

        # mutability tests
        @testset "Copy Expression" begin
            # Base.copy tests
            copied_st = copy(st)
            @test copied_st == st   # still the same expression
            @test copied_st !== st  # but not the same reference
            @test first(arguments(copied_st)) == first(arguments(copied_st))
            @test first(arguments(copied_st)) === first(arguments(copied_st))

            # Base.deepcopy tests
            st2 = (a ∧ b) ∧ c
            dc_st2 = deepcopy(st2)
            @test dc_st2 == st2
            @test dc_st2 !== st2
            @test first(arguments(dc_st2)) == first(arguments(st2))
            @test first(arguments(dc_st2)) !== first(arguments(st2))
            @test arguments(dc_st2)[2] == arguments(st2)[2]
            @test arguments(dc_st2)[2] === arguments(st2)[2]  # the second arg is a symbol so it shouldn't be copied
        end

        @testset "Parent Expressions" begin
            st2 = a ∧ b
            @test isempty(parents(st2))
            st3 = a ∧ st2
            @test st3 ∈ parents(st2)
            st4 = a ∧ st3
            @test st4 ∈ parents(st3)

            # parent removal upon argument reassignment
            st4.arguments[2] = b
            @test st4 ∉ parents(st3)
            st3.arguments[2] = st4
            @test st3 ∈ parents(st4)
        end

        @testset "Mutate Operation" begin
            # simple expression
            st2 = a ∧ b
            @test operation(st2) == ∧
            @test operations(st2) == Set{Deductive.LogicalOperation}([∧])
            st2.operation = ∨
            @test operation(st2) == ∨
            @test operations(st2) == Set{Deductive.LogicalOperation}([∨])
            st2.operation = →
            @test operation(st2) == →
            @test operations(st2) == Set{Deductive.LogicalOperation}([→])
            try
                st2.operation = ¬  # cannot set operation of different argument count
                @test false
            catch e
                @test true
            end

            # nested expression
            st3 = (a ∧ b) → (c ∨ d)
            @test operation(st3) == →
            @test operations(st3) == Set{Deductive.LogicalOperation}([→, ∧, ∨])
            st3.operation = ⟷
            @test operation(st3) == ⟷
            @test operations(st3) == Set{Deductive.LogicalOperation}([⟷, ∧, ∨])
            st3.operation = ∧
            @test operation(st3) == ∧
            @test operations(st3) == Set{Deductive.LogicalOperation}([∧, ∨])
            st3.operation = →
            arguments(st3)[1].operation = →
            @test operation(first(arguments(st3))) == →
            @test operations(st3) == Set{Deductive.LogicalOperation}([→, ∨])
            arguments(st3)[2].operation = →
            @test operation(last(arguments(st3))) == →
            @test operations(st3) == Set{Deductive.LogicalOperation}([→])
        end

        @testset "Mutate Arguments" begin
            # simple expression
            st2 = a ∧ b
            @test arguments(st2) == Vector{Deductive.AbstractExpression}([a, b])
            @test variables(st2) == Set{Deductive.LogicalSymbol}([a, b])
            st2.arguments[1] = b
            @test arguments(st2) == Vector{Deductive.AbstractExpression}([b, b])
            @test variables(st2) == Set{Deductive.LogicalSymbol}([b])
            st2.arguments[2] = a
            @test arguments(st2) == Vector{Deductive.AbstractExpression}([b, a])
            @test variables(st2) == Set{Deductive.LogicalSymbol}([a, b])
            st2.arguments[1:2] = [a, a]
            @test arguments(st2) == Vector{Deductive.AbstractExpression}([a, a])
            @test variables(st2) == Set{Deductive.LogicalSymbol}([a])

            # nested expression
            st3 = (a ∧ b) → (c ∨ d)
            arguments(st3)[1].arguments[1] = b
            @test arguments(st3) == Vector{Deductive.AbstractExpression}([b ∧ b, c ∨ d])
            @test variables(st3) == Set{Deductive.LogicalSymbol}([b, c, d])
            arguments(st3)[2].arguments[2] = a
            @test arguments(st3) == Vector{Deductive.AbstractExpression}([b ∧ b, c ∨ a])
            @test variables(st3) == Set{Deductive.LogicalSymbol}([b, c, a])
        end
    end
end


text:
module ThreadGantt

import IOCapture: capture
using PlotlyJS

export @workunit, plotgantt, capture

include("main.jl")

end


text:
using LightGraphs
using Plots
using StatsBase
using DataFrames

function computeEnergy()
    global Network
    global nodes
    global Data

    #Compute sytem energy
    E = 0
    for i in 1:length(nodes)
        Data.E[end] += nodes[i].energy
    end
end

#Computes node energy by assessing the change in preference of an inneighbor
function dE(ID,oldVal,newVal)
    global Network
    global nodes
    global T
    global J

    if newVal == nodes[ID].vote && oldVal == nodes[ID].vote
        ε = nodes[ID].energy
    elseif newVal == nodes[ID].vote && oldVal != nodes[ID].vote
        ε = nodes[ID].energy - J
    elseif newVal != nodes[ID].vote && oldVal == nodes[ID].vote
        ε = nodes[ID].energy + J
    else
        ε = nodes[ID].energy
    end

    return ε
end

#Slower: recomputes node energy thoroughly
function dE(ID)
    global Network
    global nodes
    global T
    global J

    #Goes through inneighbors and computes Potts node energy (ε)
    options = inneighbors(Network,ID)
    ε = 0
    for i in 1:length(options)
        if nodes[ID].vote == nodes[options[i]].vote
            ε += -J
        end
    end

    #Assign node energy
    return ε
end

#=Like dE(), but returns energy difference instead of assigning it automatically:
function ΔE(ID,newVal)
    global Network
    global nodes
    global T
    global J

    #Goes through inneighbors and computes Potts node energy (ε)
    options = inneighbors(Network,ID)
    ε = 0
    for i in 1:length(options)
        if nodes[ID].vote == nodes[options[i]].vote
            ε += -J
        end
    end

    #Assign node energy
    nodes[ID].energy = ε
end=#

#=(Not ready yet)
function generateEnergyDistribution()
    #Determine maximum energy
    E_max = findmax(Data.E)[1]
    global distrib = zeros(E_max+1)

    #Count energy levels
    for i in 1:length(Data.E)
        distrib()
    end

    #Normalize?

    return 0
end
=#

function Entropy()

end

function stateVariable()

end

function trackPreference()
    global nodes
    global Data

    #Count popularity of each candidate
    for i in 1:length(nodes)
        #Unrolled loop for each candidate
        #(unfortunately it makes it difficult to change candidate number)
        if nodes[i].vote == 1
            Data.c1[end] += 1
        elseif nodes[i].vote == 2
            Data.c2[end] += 1
        elseif nodes[i].vote == 3
            Data.c3[end] += 1
        elseif nodes[i].vote == 4
            Data.c4[end] += 1
        elseif nodes[i].vote == 5
            Data.c5[end] += 1
        elseif nodes[i].vote == 6
            Data.c6[end] += 1
        elseif nodes[i].vote == 7
            Data.c7[end] += 1
        elseif nodes[i].vote == 8
            Data.c8[end] += 1
        elseif nodes[i].vote == 9
            Data.c9[end] += 1
        elseif nodes[i].vote == 10
            Data.c10[end] += 1
        end
    end
end

function trackPreference(old, new)
    global nodes
    global Data
    #Copy previous distribution
    Data.c1[end] = Data.c1[end-1]
    Data.c2[end] = Data.c2[end-1]
    Data.c3[end] = Data.c3[end-1]
    Data.c4[end] = Data.c4[end-1]
    Data.c5[end] = Data.c5[end-1]
    Data.c6[end] = Data.c6[end-1]
    Data.c7[end] = Data.c7[end-1]
    Data.c8[end] = Data.c8[end-1]
    Data.c9[end] = Data.c9[end-1]
    Data.c10[end] = Data.c10[end-1]

    #Remove old value
    if old == 1
        Data.c1[end] -= 1
    elseif old == 2
        Data.c2[end] -= 1
    elseif old == 3
        Data.c3[end] -= 1
    elseif old == 4
        Data.c4[end] -= 1
    elseif old == 5
        Data.c5[end] -= 1
    elseif old == 6
        Data.c6[end] -= 1
    elseif old == 7
        Data.c7[end] -= 1
    elseif old == 8
        Data.c8[end] -= 1
    elseif old == 9
        Data.c9[end] -= 1
    elseif old == 10
        Data.c10[end] -= 1
    end

    #Add new value
    if new == 1
        Data.c1[end] +=1
    elseif new == 2
        Data.c2[end] +=1
    elseif new == 3
        Data.c3[end] +=1
    elseif new == 4
        Data.c4[end] +=1
    elseif new == 5
        Data.c5[end] +=1
    elseif new == 6
        Data.c6[end] +=1
    elseif new == 7
        Data.c7[end] +=1
    elseif new == 8
        Data.c8[end] +=1
    elseif new == 9
        Data.c9[end] +=1
    elseif new == 10
        Data.c10[end] +=1
    end
end

#Plotting
#TODO: get it from storage
function plotAnalysis(t,dir)
    global Data
    global nodes
    global Network

    #Energy evolution
    plot1=plot(1:t+1, Data.E[1:t+1]#=/Data.E[1]=#, legend=false)
    xlabel!("Time")
    ylabel!("E/E_0")
    #xlims!(51000,52000)
    #ylims!(2,2.5)
    title!("Energy")
    png(string("Data/",dir,"/","Energy"))

    #Energy distribution
    energyDistribution=counts(-Data.E)
    plot2=plot(1:length(energyDistribution), energyDistribution, legend=false)
    title!("Energy distribution")
    png(string("Data/",dir,"/","EnergyDistribution"))

    #Inneighbor histogram
    noInneighbors = []
    for i in 1:length(nodes)
        push!(noInneighbors, length(inneighbors(Network, i)))
    end
    plot3=histogram(noInneighbors, bins = 15)
    title!("Inneighbor histogram")
    png(string("Data/",dir,"/","InneighborHistogram"))

    #Preference evolution
    #=UNDER CONSTRUCTION
    plot3=plot(1:t+1,Preference,legend=true)
    xlabel!("Time")
    ylabel!("E/E_0")
    #xlims!(51000,52000)
    #ylims!(2,2.5)
    title!("Energy")
    png("Data\\Energy")
    =#
end

max_stars_repo_path: Scripts/Prototype 2/analysis.jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_path: Scripts/Prototype 2/analysis.jl

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_path: Scripts/Prototype 2/analysis.jl

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
include("header.jl")

@testset "math" begin
    o = (:delta=>0.0001,:atol=>0.01,:rtol=>0.01)
    ϵ = 0.1
    abs_gt_0(x)=sign(x) * (abs(x) + ϵ)
    abs_gt_1(x)=sign(x) * (abs(x) + 1 + ϵ)
    abs_lt_1(x)=rand() * (2-2ϵ) - (1-ϵ)
    val_lt_1(x)=rand() * (1-2ϵ) + ϵ
    val_gt_1(x)=abs(x) + 1 + ϵ
    val_gt_0(x)=abs(x) + ϵ
    val_gt_m1(x)=abs(x) - 1 + ϵ
    val_cot(x)=(abs(cot(x)) > 10 ? x+pi/2 : x) # avoid multiples of pi
    val_cotd(x)=(abs(cotd(x)) > 10 ? x+90 : x) # avoid multiples of 180
    val_tan(x)=(abs(tan(x)) > 10 ? x+pi/2 : x) # avoid pi/2, 3pi/2 etc.
    val_tand(x)=(abs(tand(x)) > 10 ? x+90 : x) # avoid 90, 270, etc.
    val_logb(x)=(x = abs(x) + ϵ; abs(x-1) < ϵ ? x+0.5 : x)     # avoid log base close to 1
    id = identity
    val(v) = (x->v)

    @test randcheck(acos,abs_lt_1;o...)
    @test randcheck(acosd,abs_lt_1;o...)
    @test randcheck(acosh,val_gt_1;o...)
    @test randcheck(acot,abs_gt_0;o...)
    @test randcheck(acotd,abs_gt_0;o...)
    @test randcheck(acoth,abs_gt_1;o...)
    @test randcheck(acsc,abs_gt_1;o...)
    @test randcheck(acscd,abs_gt_1;o...)
    @test randcheck(acsch,abs_gt_0;o...)
    @test randcheck(asec,abs_gt_1;o...)
    @test randcheck(asecd,abs_gt_1;o...)
    @test randcheck(asech,val_lt_1;o...)
    @test randcheck(asin,abs_lt_1;o...)
    @test randcheck(asind,abs_lt_1;o...)
    @test randcheck(asinh;o...)
    @test randcheck(atan;o...)
    @test randcheck(atand;o...)
    @test randcheck(atanh,abs_lt_1;o...)
    @test randcheck(cbrt;o...)
    @test randcheck(clamp,identity,val(-0.5),val(0.5); args=1,o...)
    #@test randcheck(clamp!) # overwrites
    @test randcheck(cos;o...)
    @test randcheck(cosc;o...)
    @test randcheck(cosd;o...)
    @test randcheck(cosh;o...)
    @test randcheck(cospi;o...)
    @test randcheck(cot,val_cot;o...)
    @test randcheck(cotd,val_cotd;o...)
    @test randcheck(coth,abs_gt_0;o...)
    @test randcheck(csc,val_cot;o...)
    @test randcheck(cscd,val_cotd;o...)
    @test randcheck(csch,abs_gt_0;o...)
    @test randcheck(deg2rad;o...)
    @test randcheck(exp;o...)
    @test randcheck(exp10,abs_lt_1;o...)
    @test randcheck(exp2;o...)
    @test randcheck(expm1;o...)
    @test randcheck(exponent;o...)
    #@test randcheck(frexp) # returns 2 values
    @test randcheck(hypot,id,id;o...)
    @test randcheck(ldexp,id,val(rand(-5:5)); args=1,o...)
    @test randcheck(log,val_gt_0;o...)
    @test randcheck(log,val_logb,val_gt_0;o...)
    @test randcheck(log10,val_gt_0;o...)
    @test randcheck(log1p,val_gt_m1;o...)
    @test randcheck(log2,val_gt_0;o...)
    @test randcheck(max,id,id;o...)
    @test randcheck(min,id,id;o...)
    #@test randcheck(minmax) # returns 2 values
    @test randcheck(mod2pi;o...)
    #@test randcheck(modf)   # returns 2 values
    @test randcheck(rad2deg;o...)
    @test randcheck(rem2pi,id,val(RoundNearest); args=1,o...)
    @test randcheck(sec,val_tan;o...)
    @test randcheck(secd,val_tand;o...)
    @test randcheck(sech;o...)
    @test randcheck(significand;o...)
    @test randcheck(sin;o...)
    @test randcheck(sinc;o...)
    #@test randcheck(sincos) # returns 2 values
    @test randcheck(sind;o...)
    @test randcheck(sinh;o...)
    @test randcheck(sinpi;o...)
    @test randcheck(sqrt,val_gt_0;o...)
    @test randcheck(tan, val_tan;o...)
    @test randcheck(tand, val_tand;o...)
    @test randcheck(tanh;o...)
end


text:
@adjoint function evaluate(s::SqEuclidean, x::AbstractVector, y::AbstractVector)
  δ = x .- y
  sum(abs2, δ), Δ -> begin
    x̄ = (2 * Δ) .* δ
    (nothing, x̄, -x̄)
  end
end

@adjoint function evaluate(s::Euclidean, x::AbstractVector, y::AbstractVector)
  D = x.-y
  δ = sqrt(sum(abs2,D))
  δ, Δ -> begin
    x̄ = Δ .* D / (δ + eps(δ))
    (nothing, x̄, -x̄)
  end
end

@adjoint function evaluate(s::DotProduct, x::AbstractVector, y::AbstractVector)
  dot(x,y), Δ -> begin
    (nothing, Δ.*y, Δ.*x)
  end
end


text:
using MAT

@testset "find2dPIS" begin
    # MATLAB output
    test = matread("resources/find2dPIS.mat")

    A = [1 0; 1 1]
    B = [1 -1; 1 -2]
    C = [0 1]
    D = [0 0]
    ABCD = [A B; C D]
    u = 1 / pi

    s = find2dPIS(u, ABCD, false)
    @test test["s"] == s
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
# runs all scripts in the facet_verifications directory
script_dir = "./script/facet_verifications/"
println("running all facet verification scripts")

expected_num_scripts = 6

pass_list = []
for file in readdir(script_dir)
    if occursin(r"run_all\.jl$", file)
        continue
    elseif occursin(r"\.jl$", file)
        println("running $file")
        @time test_pass = include(file)
        push!(pass_list, test_pass)
        println("test pass? : ", test_pass)
    end
end

if expected_num_scripts == length(pass_list)
    println("Ran all expected scripts")
else
    throw(ErrorException("Not all scripts ran"))
end

if all(pass_list)
    println("All facets verified.")
else
    throw(ErrorException("An failure occurred while running facet verifications"))
end


text:
using MatterPower
using DoubleExponentialFormulas
using Plots, LaTeXStrings
using Dierckx
using DelimitedFiles, Printf

# %% Specify a redshift
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
ired = 1

# %% Read in linear and non-linear power spectra computed by CLASS and spline interpolate in k
filename = @sprintf("data/matterpower05_z%1d_pk.dat", ired)
d = readdlm(filename, comments = true)
pklin_class = Spline1D(d[:, 1], d[:, 2])
filename = @sprintf("data/matterpower05_z%1d_pk_nl.dat", ired)
d = readdlm(filename, comments = true)
pknl_class = Spline1D(d[:, 1], d[:, 2])

# %% Define a function to return a linear matter power spectrum (in units of Mpc^3/h^3)
# as a function of the comoving wavenumber, kovh, in units of h/Mpc.
# Here is an example using Einstein & Hu's analytical transfer function

# Cosmological parameters
As, ns, kpivot = 2.097e-9, 0.9652, 0.05
Ωm, ΩΛ, Ωb = 0.315, 0.685, 0.049
Ωk = 1 - Ωm - ΩΛ
h0 = 0.674
ωm, fb = Ωm * h0^2, Ωb / Ωm

# Tabulate linear growth factor as a function of scale factor, a
sol = setup_growth(Ωm, ΩΛ)
a = 1 / (1 + redshift[ired])
D1 = sol(a)[1]

pk(kovh) =
   D1^2 *
   As *
   (kovh * h0 / kpivot)^(ns - 1) *
   (2 * kovh^2 * 2998^2 / 5 / Ωm)^2 *
   t_nowiggle(kovh * h0, ωm, fb)^2 *
   2 *
   π^2 / kovh^3

# %% Alternatively you may read in pre-computed data and define a spline function
# using Dierckx
# pk = Spline1D(tabulated_kovh, tabulated_power_spectrum)

# %% Compute the non-linear power spectrum using halofit
# Get three parameters [kσ, neff, C] needed for the halofit transform
p = kσ, neff, C = setup_halofit(pk)

# Calculate Ωm(z) = Ωm(z=0)(1+z)^3 / E^2(z), where E^2(z) = H^2(z) / H0^2
z = redshift[ired]
E2 = Ωm * (1 + z)^3 + Ωk * (1 + z)^2 + ΩΛ
Ωmz = Ωm * (1 + z)^3 / E2

# Define a function to return a halofit non-linear power spectrum
pknl(kovh) = halofit(pk, p, Ωmz, kovh) # Mpc^3/h^3

# %% Compute correlation functions
kmin, kmax = 1e-4, Inf
R = Float64.(10:2:200)
ξ1 = zeros(length(R))
ξ2 = zeros(length(R))
ξ3 = zeros(length(R))
ξ4 = zeros(length(R))
for i = 1:length(R)
   func1(k) = k^2 * pk(k) * sin(k * R[i]) / (k * R[i]) / 2 / π^2
   ξ1[i], err = quaddeo(func1, R[i], 0, kmin, kmax)
   func2(k) = k^2 * pknl(k) * sin(k * R[i]) / (k * R[i]) / 2 / π^2
   ξ2[i], err = quaddeo(func2, R[i], 0, kmin, kmax)
   func3(k) = k^2 * pklin_class(k) * sin(k * R[i]) / (k * R[i]) / 2 / π^2
   ξ3[i], err = quaddeo(func3, R[i], 0, kmin, kmax)
   func4(k) = k^2 * pknl_class(k) * sin(k * R[i]) / (k * R[i]) / 2 / π^2
   ξ4[i], err = quaddeo(func4, R[i], 0, kmin, kmax)
end

# %% Plot results and save to "xi.pdf"
p = plot(
   R,
   R .^ 2 .* ξ4,
   m = 2,
   c = :red,
   lab = "Halofit ξ(R) with BAO",
   ylab = L"R^2\xi(R)",
   xlab = L"R~[Mpc/h]",
   title = "Redshift = $z",
)
p =
   plot!(R, R .^ 2 .* ξ2, ls = :dash, c = :red, lab = "Halofit ξ(R) w/o BAO")
p = plot!(R, R .^ 2 .* ξ3, m = 2, c = 1, lab = "Linear ξ(R) with BAO")
p = plot!(R, R .^ 2 .* ξ1, ls = :dash, c = 1, lab = "Linear ξ(R) w/o BAO")
savefig("xi.pdf")
display(p)


text:
export char_poly

"""
`char_poly(A)` returns the characteristic polynomial of
the exact matrix `A`.
"""
function char_poly(A::AbstractMatrix{T}) where {T<:TypeX}
    r, c = size(A)
    @assert r == c "Matrix must be square"

    xI = zeros(SimplePolynomial, r, r)
    x = getx(T)
    for i = 1:r
        xI[i, i] = x
    end

    f = detx(xI - A)
    p = integerize(numerator(f))
end


text:
# Create mutable structs to cache function pointers
macro create_func_ptr_struct(name, syms)
    e = :(mutable struct $name end)
    for s in eval(syms)
        push!(e.args[3].args, :($s::Ptr{Nothing}))
    end
    push!(e.args[3].args, :(
        $name() = new( fill(C_NULL,length($syms))... )
    ) )
    e
end

include("libgr_syms.jl")
include("libgrm_syms.jl")
include("libgr3_syms.jl")

@create_func_ptr_struct LibGR_Ptrs libGR_syms
@create_func_ptr_struct LibGRM_Ptrs libGRM_syms
@create_func_ptr_struct LibGR3_Ptrs libGR3_syms

const libGR_ptrs = LibGR_Ptrs()
const libGRM_ptrs = LibGRM_Ptrs()
const libGR3_ptrs = LibGR3_Ptrs()

const libs_loaded = Ref(false)

"""
    load_libs(always = false)

    Load shared GR libraries from either GR_jll or from GR tarball.
    always is a boolean flag that is passed through to
"""
function load_libs(always::Bool = false)
    if gr_provider[] == "BinaryBuilder"
        try
            @eval GR import GR_jll
            libGR_handle[] = GR_jll.libGR_handle
            libGR3_handle[] = GR_jll.libGR3_handle
            libGRM_handle[] = GR_jll.libGRM_handle
            grdir[] = joinpath(dirname(GR_jll.libGR_path), "..")
        catch err
            @error "Error importing GR_jll:" err
            ENV["GRDIR"] = ""
            depsfile_succeeded[] = false
            __init__()
            load_libs()
            return
        end
    else
        # Global grdir should be set in __init__

        flag = occursin("site-packages", grdir[])
        loadpath = grdir[]
        if flag
            ENV["GKS_FONTPATH"] = grdir[]
        elseif os != :Windows
            loadpath = joinpath(loadpath, "lib")
        else
            loadpath = joinpath(loadpath, "bin")
        end
        push!(Base.DL_LOAD_PATH, loadpath)

        libGR_handle[] = Libdl.dlopen(libGR)
        libGR3_handle[] = Libdl.dlopen(libGR3)
        libGRM_handle[] = Libdl.dlopen(libGRM)
    end
    @debug "Library handles" libGR_handle[] libGR3_handle[] libGRM_handle[]

    libs_loaded[] = true

    check_env[] = true
    init(always)
end

function get_func_ptr(handle::Ref{Ptr{Nothing}}, ptrs::Union{LibGR_Ptrs, LibGRM_Ptrs, LibGR3_Ptrs}, func::Symbol, loaded=libs_loaded[])
    if !loaded
        load_libs(true)
    end
    s = getfield(ptrs, func)
    if s == C_NULL
        s = Libdl.dlsym(handle[], func)
        setfield!(ptrs, func, s)
    end
    return getfield(ptrs,func)
end

libGR_ptr(func) = get_func_ptr(libGR_handle, libGR_ptrs, func)
libGRM_ptr(func) = get_func_ptr(libGRM_handle, libGRM_ptrs, func)
libGR3_ptr(func) = get_func_ptr(libGR3_handle, libGR3_ptrs, func)

precompile(get_func_ptr, (Base.RefValue{Ptr{Nothing}},LibGR_Ptrs, Symbol) )
precompile(get_func_ptr, (Base.RefValue{Ptr{Nothing}},LibGRM_Ptrs, Symbol) )
precompile(get_func_ptr, (Base.RefValue{Ptr{Nothing}},LibGR3_Ptrs, Symbol) )
precompile(libGR_ptr, (Symbol,))
precompile(libGRM_ptr, (Symbol,))
precompile(libGR3_ptr, (Symbol,))
precompile(load_libs, (Bool,))


text:
function analyze(c::Constructor)
    arg_list = get_args(c)

    args = Any[]
    for arg in arg_list
        proxy = get_proxy(spelling(cu_type(arg)))
        (proxy == Union{}) ? push!(args, arg) : push!(args, proxy)
    end
    return WrappedConstructor(spelling(c),
                         c,
                         cindex.getCursorLexicalParent(c),
                         args)
end
analyze(d::Destructor) = WrappedDestructor(d, cindex.getCursorLexicalParent(d))
analyze(m::CXXMethod) = analyze_method(m)

function analyze(c::ClassDecl)

    constructors = WrappedConstructor[]
    methods = WrappedMethod[]
    for node in children(c)
        if isa(node, CXXMethod)
            push!(methods, analyze(node))
        elseif isa(node, Constructor)
            push!(constructors, analyze(node))
        end
    end

    WrappedClass(
    spelling(c),
    constructors,
    methods,
    c
    )
end

function analyze(x::TranslationUnit)
    classnodes = WrappedClass[]
    for node in children(x)
        if isa(node, ClassDecl)
            push!(classnodes, analyze(node))
        end
    end
    WrappedHeader(x, classnodes)
end


text:
"""
Orbit visualization tests. Do plots run without error?
"""
module VisualizationTests

using LinearAlgebra, Unitful, UnitfulAstro, GeneralAstrodynamics, Test
using DifferentialEquations
using Plots

@testset verbose=false "Trajectory Plots" begin

    # CR3BP trajectory
    L2Halo = halo(SunEarth; Az = 30_000u"km", L = 2)
    traj   = propagate(L2Halo...)

    try
        plot(traj)
        @test true
    catch e
        @test false
    end

    # Discrete trajectory
    traj = map(t -> Orbit(traj, t), solution(traj).t)
    try
        plot(traj)
        @test true
    catch e
        @test false
    end

    # Close all plots
    Plots.closeall()

end

@testset verbose=false "Manifold Plots" begin

    # CR3BP manifold
    L2Halo = halo(SunEarth; Az = 30_000u"km", L = 2)

    try
        plot(manifold(L2Halo...))
        @test true
    catch e
        @test false
    end

    # Close all plots
    Plots.closeall()

end

@testset verbose=false "Energy Plots" begin

    # CR3BP trajectory
    L2Halo = halo(EarthMoon; Az = 15_000u"km", L = 2)

    # CR3BP zero velocity curves
    try
        zerovelocityplot(L2Halo.orbit)
        @test true
    catch e
        @test false
    end

    # Close all plots
    Plots.closeall()

end

end # module


text:
const DATA_DOC = """
    Module Data

The module Data is created in the module where `@init_parallel_kernel` is called from. It provides the following types:

--------------------------------------------------------------------------------
    Data.Number

The type of numbers used by @zeros, @ones and @rand and in all array types of module `Data` (selected with argument `numbertype` of [`@init_parallel_kernel`](@ref)).

--------------------------------------------------------------------------------
    Data.Array{ndims}

Expands to `Data.Array{numbertype, ndims}`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the datatype `Data.Array` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (Array for Threads and CUDA.CuArray or CUDA.CuDeviceArray for CUDA; [`@parallel`](@ref) and [`@parallel_indices`](@ref) convert CUDA.CuArray automatically to CUDA.CuDeviceArray in kernels when required).

--------------------------------------------------------------------------------
!!! note "Advanced"
        Data.DeviceArray{ndims}

    Expands to `Data.DeviceArray{numbertype, ndims}`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the datatype `Data.DeviceArray` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (Array for Threads and CUDA.CuDeviceArray for CUDA).

    !!! warning
        This datatype is not intended for explicit manual usage. [`@parallel`](@ref) and [`@parallel_indices`](@ref) convert CUDA.CuArray automatically to CUDA.CuDeviceArray in kernels when required.
"""

function Data_cuda(numbertype::DataType)
    :(baremodule Data # NOTE: there cannot be any newline before 'module Data' or it will create a begin end block and the module creation will fail.
        import CUDA
        Number         = $numbertype
        Array{N}       = CUDA.CuArray{$numbertype, N}
        DeviceArray{N} = CUDA.CuDeviceArray{$numbertype, N}
    end)
end

function Data_threads(numbertype::DataType)
    :(baremodule Data # NOTE: there cannot be any newline before 'module Data' or it will create a begin end block and the module creation will fail.
        import Base
        Number         = $numbertype
        Array{N}       = Base.Array{$numbertype, N}
        DeviceArray{N} = Base.Array{$numbertype, N}
    end)
end

function Data_none()
    :(baremodule Data # NOTE: there cannot be any newline before 'module Data' or it will create a begin end block and the module creation will fail.
    end)
end


text:
# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function record!(it::Int64, issp::Int64, iss::Int64, pac, pap)

    for ipw in pac.activepw
        rinterpolatew = pap[ipw].ss[issp].rinterpolatew
        for rfield in pac.rfields
            pw = view(pap[ipw].w1[:t][rfield],:)
            recs = pap[ipw].ss[issp].records[rfield][it]
            w=rinterpolatew[rfield]
            mul!(recs, w, pw)
        end
    end
end


text:
using Distributed

@everywhere function host_id()
    return [gethostname() myid()]
end

vs = Any[]
for worker in workers()
    s = @spawnat worker host_id()
    push!(vs, s)
end

vhost = Array{String, 1}()
vid = Array{Int, 1}()
for (index, s) in enumerate(vs)
    a, b = fetch(s)
    push!(vhost, a)
    push!(vid, b)
end

idx_id = sortperm(vid)
vhost = vhost[idx_id]
vid = vid[idx_id]
onehost = unique(vhost)
oneid = vid[indexin(onehost, vhost)]
for index in 1:length(oneid)
    println("$(onehost[index]) $(oneid[index])")
end

max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples


text:
function GetExprMetric(file)
    df = DataFrame(CSV.File(file,types=Dict(1=>String)))
    exprs = []
    for (d,f) ∈ eachrow(df[df.Coefficient .!= "0",[:Index,:Coefficient]])
        #println((d,f))
        d = d == "-" ? "" : d
        d = replace(d,"th"=>"θ")
        ps = "ψ"*d
        push!(exprs,"($f)*$(ps)(r,x)")
    end
    func_def = "$(join(exprs," + "))"
    println(func_def)
    Meta.parse(func_def)
end

function GetExprDerivative(file)
    df = DataFrame(CSV.File(file,types=Dict(1=>String)))
    exprs = []
    for (d,f) ∈ eachrow(df[df.Coefficient .!= "0",:])
        #println((d,f))
        d = d == "-" ? "" : d
        d = replace(d,"th"=>"θ")
        ps = "ψ"*d
        push!(exprs,"($f)*$(ps)(r,x)")
    end
    func_def = "$(join(exprs," + "))"
    println(func_def)
    Meta.parse(func_def)
end

function MakeδT(file; additional_params = ())
    thisexpr = GetExprMetric(file)
    @eval δT(r,x,a,m,ω,s,ψ11,ψ12,ψ22,ψ1,ψ2,ψ,$(additional_params...); M=1) = $thisexpr
    δT
end

#δT = MakeδT("TuekolskyShifts.csv")

function Make∂ωT(file; additional_params = ())
    thisexpr2 = GetExprDerivative(file)
    @eval ∂ωT(r,x,a,m,ω,s,ψ11,ψ12,ψ22,ψ1,ψ2,ψ,$(additional_params...); M=1) = $thisexpr2
    ∂ωT
end

#∂ωT = Make∂ωT("TuekolskyFrequencyDerivative.csv")

max_stars_repo_name: Potatoasad/KerrQNMShifts

max_issues_repo_name: Potatoasad/KerrQNMShifts

max_forks_repo_name: Potatoasad/KerrQNMShifts


text:
"""

# d_separation

$(SIGNATURES)

Computes the d_separation between 2 sets of nodes conditioned on a third set.

### Required arguments
```julia
d_separation(
* `d::DAG`                             : DAG
* `f::SymbolList`                      : First set
*  s::SymbolList`                      : Second set
)
```

### Keyword arguments
```julia
* `c::SymbolListOrNothing=nothing`     : Conditioning set
* `debug=false`                        : Trace execution
```

### Returns
```julia
* `res::Bool`                          : Boolean result of test
```

# Extended help

### Example

### d_separation between mechanics and statistics, conditioning on algebra
```julia
using StructuralCausalModels, CSV

df = DataFrame!(CSV.File(scm_path("..", "data", "marks.csv"));

d = OrderedDict(
  :mechanics => [:vectors, :algebra],
  :vectors => [:algebra],
  :analysis => [:algebra],
  :statistics => [:algebra, :analysis]
);

dag = DAG("marks", d, df);
d_separation(marks, [:statistics], [:mechanics]; c=[:algebra]))
```
### Acknowledgements

Original author:                       Giovanni M. Marchetti

Translated to Julia:                   Rob J Goedman

### License

The R package ggm is licensed under License: GPL-2.

The Julia translation is licenced under: MIT.

Part of the API, exported.
"""
function d_separation(d::DAG, f::SymbolList, s::SymbolList;
  c::SymbolListOrNothing=nothing, debug=false)

  e = induced_covariance_graph(d, vcat(f, s), c; debug=debug)
  sum(e[f, s]) == 0

end

export
  d_separation


text:
include("registry.jl")
include("types.jl")
include("print.jl")

include("lib.jl")

include("substitude.jl")
include("match.jl")
include("simplify.jl")


text:
export NodeNil, is_nil

import .Nodes.is_nil

struct NodeNil <: Ranger.AbstractNode
    id::UInt32
    name::String

    function NodeNil()
        new(0, "_Nil_")
    end
end

function is_nil(node::NodeNil)
    node.id == 0
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using QuantumOpticsControlExtensions
using Test

@testset "QuantumOpticsControlExtensions.jl" begin
    # Write your tests here.
end


text:
@testset "others" begin
    @testset "TuringWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dW1 = Wishart(dim + 4, A)
        dW2 = TuringWishart(dim + 4, A)
        dW3 = TuringWishart(dW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, meanlogdet, entropy, cov, var)
            @test F(dW1) == F(dW2) == F(dW3)
        end
        @test Matrix(mode(dW1)) == mode(dW2) == mode(dW3)
        xw = rand(dW2)
        @test insupport(dW1, xw)
        @test insupport(dW2, xw)
        @test insupport(dW3, xw)
        @test logpdf(dW1, xw) == logpdf(dW2, xw) == logpdf(dW3, xw)
    end

    @testset "TuringInverseWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dIW1 = InverseWishart(dim + 4, A)
        dIW2 = TuringInverseWishart(dim + 4, A)
        dIW3 = TuringInverseWishart(dIW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, cov, var)
            @test F(dIW1) == F(dIW2) == F(dIW3)
        end
        @test Matrix(mode(dIW1)) == mode(dIW2) == mode(dIW3)
        xiw = rand(dIW2)
        @test insupport(dIW1, xiw)
        @test insupport(dIW2, xiw)
        @test insupport(dIW3, xiw)
        @test logpdf(dIW1, xiw) == logpdf(dIW2, xiw) == logpdf(dIW3, xiw)
    end

    @testset "TuringMvNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                C = Matrix{Float64}(I, 3, 3)
                d1 = TuringMvNormal(m, C)
            elseif TD <: TuringDiagMvNormal
                C = ones(3)
                d1 = TuringMvNormal(m, C)
            else
                C = 1.0
                d1 = TuringMvNormal(m, C)
            end
            d2 = MvNormal(m, C)

            @testset "$F" for F in (length, size)
                @test F(d1) == F(d2)
            end

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test isapprox(logpdf(d1, x1), logpdf(d2, x1), rtol = 1e-6)
            @test isapprox(logpdf(d1, x2), logpdf(d2, x2), rtol = 1e-6)
        end
    end

    @testset "TuringMvLogNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                C = Matrix{Float64}(I, 3, 3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            elseif TD <: TuringDiagMvNormal
                C = ones(3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            else
                C = 1.0
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            end
            d2 = MvLogNormal(MvNormal(m, C))

            @test length(d1) == length(d2)

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test isapprox(logpdf(d1, x1), logpdf(d2, x1), rtol = 1e-6)
            @test isapprox(logpdf(d1, x2), logpdf(d2, x2), rtol = 1e-6)

            x2[:, 1] .= -1
            @test isinf(logpdf(d1, x2)[1])
            @test isinf(logpdf(d2, x2)[1])
        end
    end

    @testset "TuringUniform" begin
        @test logpdf(TuringUniform(), 0.5) == 0
        if AD == "All" || AD == "Tracker"
            @test logpdf(TuringUniform(), param(0.5)) == 0
        end
    end

    if AD == "All" || AD == "Tracker"
        @testset "Semicircle" begin
            @test Tracker.data(logpdf(Semicircle(1.0), param(0.5))) == logpdf(Semicircle(1.0), 0.5)
        end
    end

    @testset "TuringPoissonBinomial" begin
        d1 = TuringPoissonBinomial([0.5, 0.5])
        d2 = PoissonBinomial([0.5, 0.5])
        @test quantile(d1, 0.5) == quantile(d2, 0.5)
        @test minimum(d1) == minimum(d2)
    end

    @testset "Inverse of pi" begin
        @test 1/pi == inv(pi)
    end

    @testset "Cholesky" begin
        A = rand(3, 3)'; A = A + A' + 3I;
        C = cholesky(A; check = true)
        factors, info = DistributionsAD.turing_chol(A, true)
        @test factors == C.factors
        @test info == C.info
    end

    function test_reverse_mode_ad( f, ȳ, x...; rtol=1e-6, atol=1e-6)
        # Perform a regular forwards-pass.
        y = f(x...)

        # Use finite differencing to compute reverse-mode sensitivities.
        x̄s_fdm = FDM.j′vp(central_fdm(5, 1), f, ȳ, x...)

        if AD == "All" || AD == "Zygote"
            # Use Zygote to compute reverse-mode sensitivities.
            y_zygote, back_zygote = Zygote.pullback(f, x...)
            x̄s_zygote = back_zygote(ȳ)

            # Check that Zygpte forwards-pass produces the correct answer.
            @test isapprox(y, y_zygote, atol=atol, rtol=rtol)

            # Check that Zygote reverse-mode sensitivities are correct.
            @test all(zip(x̄s_zygote, x̄s_fdm)) do (x̄_zygote, x̄_fdm)
                isapprox(x̄_zygote, x̄_fdm; atol=atol, rtol=rtol)
            end
        end

        if AD == "All" || AD == "ReverseDiff"
            test_rd = length(x) == 1 && y isa Number
            if test_rd
                # Use ReverseDiff to compute reverse-mode sensitivities.
                if x[1] isa Array
                    x̄s_rd = similar(x[1])
                    tp = ReverseDiff.GradientTape(x -> f(x), x[1])
                    ReverseDiff.gradient!(x̄s_rd, tp, x[1])
                    x̄s_rd .*= ȳ
                    y_rd = ReverseDiff.value(tp.output)
                    @assert y_rd isa Number
                else
                    x̄s_rd = [x[1]]
                    tp = ReverseDiff.GradientTape(x -> f(x[1]), [x[1]])
                    ReverseDiff.gradient!(x̄s_rd, tp, [x[1]])
                    y_rd = ReverseDiff.value(tp.output)[1]
                    x̄s_rd = x̄s_rd[1] * ȳ
                    @assert y_rd isa Number
                end

                # Check that ReverseDiff forwards-pass produces the correct answer.
                @test isapprox(y, y_rd, atol=atol, rtol=rtol)

                # Check that ReverseDiff reverse-mode sensitivities are correct.
                @test isapprox(x̄s_rd, x̄s_fdm[1]; atol=atol, rtol=rtol)
            end
        end

        if AD == "All" || AD == "Tracker"
            # Use Tracker to compute reverse-mode sensitivities.
            y_tracker, back_tracker = Tracker.forward(f, x...)
            x̄s_tracker = back_tracker(ȳ)

            # Check that Tracker forwards-pass produces the correct answer.
            @test isapprox(y, Tracker.data(y_tracker), atol=atol, rtol=rtol)

            # Check that Tracker reverse-mode sensitivities are correct.
            @test all(zip(x̄s_tracker, x̄s_fdm)) do (x̄_tracker, x̄_fdm)
                isapprox(Tracker.data(x̄_tracker), x̄_fdm; atol=atol, rtol=rtol)
            end
        end
    end
    _to_cov(B) = B + B' + 2 * size(B, 1) * Matrix(I, size(B)...)

    @testset "logsumexp" begin
        x, y = rand(3), rand()
        test_reverse_mode_ad(logsumexp, y, x; rtol=1e-8, atol=1e-6)
    end

    @testset "zygote_ldiv" begin
        A = rand(3, 3)'; A = A + A' + 3I;
        B = copy(A)
        Ȳ = rand(3, 3)
        @test DistributionsAD.zygote_ldiv(A, B) == A \ B
        test_reverse_mode_ad((A,B)->DistributionsAD.zygote_ldiv(A,B), Ȳ, A, B)
    end

    @testset "logdet" begin
        rng, N = MersenneTwister(123456), 7
        y, B = randn(rng), randn(rng, N, N)
        test_reverse_mode_ad(B->logdet(cholesky(_to_cov(B))), y, B; rtol=1e-8, atol=1e-6)
        test_reverse_mode_ad(B->logdet(cholesky(Symmetric(_to_cov(B)))), y, B; rtol=1e-8, atol=1e-6)
    end

    @testset "fill" begin
        if AD == "All" || AD == "Tracker"
            @test fill(param(1.0), 3) isa TrackedArray
        end
        rng = MersenneTwister(123456)
        test_reverse_mode_ad(x->fill(x, 7), randn(rng, 7), randn(rng))
        test_reverse_mode_ad(x->fill(x, 7, 11), randn(rng, 7, 11), randn(rng))
        test_reverse_mode_ad(x->fill(x, 7, 11, 13), rand(rng, 7, 11, 13), randn(rng))
    end
    @testset "Tracker, Zygote and ReverseDiff + MvNormal" begin
        rng, N = MersenneTwister(123456), 11
        B = randn(rng, N, N)
        m, A = randn(rng, N), B' * B + I

        # Generate from the TuringDenseMvNormal
        d = TuringDenseMvNormal(m, A)
        x = rand(d)

        # Check that the logpdf agrees with MvNormal.
        d_ref = MvNormal(m, PDMat(A))
        @test logpdf(d, x) ≈ logpdf(d_ref, x)

        test_reverse_mode_ad((m, B, x)->logpdf(MvNormal(m, _to_cov(B)), x), randn(rng), m, B, x)
        test_reverse_mode_ad((m, B, x)->logpdf(TuringMvNormal(m, _to_cov(B)), x), randn(rng), m, B, x)
        test_reverse_mode_ad((m, B, x)->logpdf(TuringMvNormal(m, Symmetric(_to_cov(B))), x), randn(rng), m, B, x)
    end

    @testset "Entropy" begin
        sigmas = exp.(randn(10))
        d1 = TuringDiagMvNormal(zeros(10), sigmas)
        d2 = MvNormal(zeros(10), sigmas)

        @test isapprox(entropy(d1), entropy(d2), rtol = 1e-6)
    end

    @testset "Params" begin
        m = rand(10)
        sigmas = randexp(10)

        d = TuringDiagMvNormal(m, sigmas)
        @test params(d) == (m, sigmas)

        d = TuringScalMvNormal(m, sigmas[1])
        @test params(d) == (m, sigmas[1])
    end
end


text:
# Test raw wrapper with direct library calls

@testset "create small problem and solve" begin
    scip__ = Ref{Ptr{SCIP.SCIP_}}() # SCIP**
    rc = SCIP.SCIPcreate(scip__)
    @test rc == SCIP.SCIP_OKAY

    scip_ = scip__[] # dereference to SCIP*
    @test scip_ != C_NULL

    rc = SCIP.SCIPincludeDefaultPlugins(scip_)
    @test rc == SCIP.SCIP_OKAY

    # disable output
    rc = SCIP.SCIPsetIntParam(scip_, "display/verblevel", 0)
    @test rc == SCIP.SCIP_OKAY

    # create problem
    rc = SCIP.SCIPcreateProbBasic(scip_, "")
    @test rc == SCIP.SCIP_OKAY

    # add variable: x >= 0, objcoef: 1
    var__ = Ref{Ptr{SCIP.SCIP_VAR}}()
    rc = SCIP.SCIPcreateVarBasic(scip_, var__, "x",
                                 0.0, SCIP.SCIPinfinity(scip_),
                                 1.0, SCIP.SCIP_VARTYPE_CONTINUOUS)
    @test rc == SCIP.SCIP_OKAY
    var_ = var__[]
    @test var_ != C_NULL
    rc = SCIP.SCIPaddVar(scip_, var_)
    @test rc == SCIP.SCIP_OKAY

    # add constraint: 2x >= 3   ( really: 3 <= 2 * x <= inf )
    cons__ = Ref{Ptr{SCIP.SCIP_CONS}}()
    rc = SCIP.SCIPcreateConsBasicLinear(scip_, cons__, "c", 0, C_NULL, C_NULL,
                                        3.0, SCIP.SCIPinfinity(scip_))
    @test rc == SCIP.SCIP_OKAY
    cons_ = cons__[]
    @test cons_ != C_NULL
    rc = SCIP.SCIPaddCoefLinear(scip_, cons_, var_, 2.0)
    @test rc == SCIP.SCIP_OKAY
    rc = SCIP.SCIPaddCons(scip_, cons_)
    @test rc == SCIP.SCIP_OKAY

    # solve problem
    rc = SCIP.SCIPsolve(scip_)
    @test rc == SCIP.SCIP_OKAY

    # check solution value
    sol_ = SCIP.SCIPgetBestSol(scip_)
    @test sol_ != C_NULL
    val = SCIP.SCIPgetSolVal(scip_, sol_, var_)
    @test val ≈ 3.0 / 2.0

    # release variables and solver
    rc = SCIP.SCIPreleaseCons(scip_, cons__)
    @test rc == SCIP.SCIP_OKAY
    rc = SCIP.SCIPreleaseVar(scip_, var__)
    @test rc == SCIP.SCIP_OKAY
    rc = SCIP.SCIPfree(scip__)
    @test rc == SCIP.SCIP_OKAY
end

@testset "SCIP_CALL macro (@SC)" begin
    # should do nothing
    @SCIP.SC SCIP.SCIP_OKAY

    @test_throws AssertionError @SCIP.SC SCIP.SCIP_ERROR
end


text:
#=
`show` methods for phylogenetic trees
=#

function Base.show(io::IO, br::AbstractBranch)
    msg = br.label == "" ? "" : " \"" * br.label * "\" "
    if br.length === nothing
        msg *= "without length"
    else
        msg *= "of length " * string(br.length)
    end
    print(io::IO, "Branch ", msg)

    return nothing
end # function Base.show

function Base.show(io::IO, t::AbstractTree)
    if length(t.label) == 0
        print(io, "Tree")
    else
        print(io, "Tree \"", t.label, "\"")
    end

    return nothing
end # function Base.show

function Base.show(io::IO, ::MIME"text/plain", t::AbstractTree)
    if length(t.label) == 0
        print(io, "Tree")
    else
        print(io, "Tree \"", t.label, "\"")
    end
    print(":\n   ", t.rooted ? "Rooted" : "Unrooted")

    return nothing
end # function Base.show

function Base.show(io::IO, p::AbstractNode)
    if length(p.label) == 0
        print(io, "Node")
    else
        print(io, "Node \"", p.label, "\"")
    end

    return nothing
end # function Base.show

function Base.show(io::IO, ::MIME"text/plain", p::AbstractNode)
    if length(p.label) == 0
        print(io, "Node")
    else
        print(io, "Node \"", p.label, "\"")
    end
    print(":\n   ", n_neighbour(p), " neighbour(s)")

    return nothing
end # function Base.show


text:
module AnalyticBSplines

using Polynomials
using BSplinesQuasi
import BSplinesQuasi: AbstractKnotSet

include("intervals.jl")
include("ppoly.jl")

function gen_basis(t::AbstractKnotSet{k,ml,mr,T}, bl=one(T), br=one(T)) where {k,ml,mr,T}
    d(a,b) = b == 0 ? 0 : a/b
    B = [[PPoly([((t[j],t[j+1]) => [one(T)])])
          for j = 1:length(t)-1]]
    for k = 2:order(t)
        Bₖ = [d(Poly([-t[i],1]),(t[i+k-1]-t[i]))*B[end][i] +
              d(Poly([t[i+k],-1]),(t[i+k]-t[i+1]))*B[end][i+1]
              for i = 1:length(t)-k]
        push!(B,Bₖ)
    end
    B[end][1] *= bl
    B[end][end] *= br
    B
end

function scalar_op(Bₖ, t::AbstractKnotSet{k,ml,mr,T}, f::Poly=Poly([one(T)])) where {k,ml,mr,T}
    B = Bₖ[end]
    n = length(B)
    O = Matrix{T}(undef, n,n)
    l(v::Rational) = denominator(v) == 1 ? numerator(v) : 1.0*v
    l(v) = v
    for i = 1:n
        for j = 1:n
            O[i,j] = polyint(B[i]*f*B[j], l(first(t)),l(last(t)))
        end
    end
    O
end

function derop(L, R, t::AbstractKnotSet{k,ml,mr,T}, o::Integer) where {k,ml,mr,T}
    m = length(L[end])
    n = length(R[end])
    O = Matrix{T}(undef, m,n)
    l(v::Rational) = denominator(v) == 1 ? numerator(v) : 1.0*v
    l(v) = v
    d,r = divrem(o,2)
    a,b = d,d+r
    for i = 1:m
        for j = 1:n
            O[i,j] = (iseven(a) ? 1 : -1)*polyint(polyder(L[end][i],a)*polyder(R[end][j],b), l(first(t)),l(last(t)))
        end
    end
    O
end

derop(B, t, o) = derop(B, B, t, o)

export gen_basis, scalar_op, derop

end # module


text:
# ------------------------------------------------------------------------------------------
# # 认识 Julia
#
# 这个 notebook 用于速成 Julia 语法，让你了解 Julia 是多么轻量级、多么好用
# —— 就像你最喜欢的其他高级语言一样！
#
#
# ## 目录
# - [字符串（Strings）](#字符串)
# - [数据结构（Data structures）](#数据结构)
#     - [元组（Tuples）](#元组（Tuples）)
#     - [字典（Dictionaries）](#字典（Dictionaries）)
#     - [数组（Arrays）](#数组（Arrays）)
# - [循环（Loops）](#循环)
#     - [`for` 循环](#for-循环)
#     - [`while` 循环](#while-循环)
# - [条件判断（Conditionals）](#条件判断)
#     - [`if` 分支](#if-分支)
#     - [三元运算符](#三元运算符)
# - [函数（Functions）](#函数)
#     - [如何声明一个函数](#如何声明一个函数)
#     - [Julia 中的鸭子类型（Duck-typing）](#Julia-中的鸭子类型（Duck-typing）)
#     - [调用函数](#调用函数)
#     - [原地修改与非原地修改函数（Mutating vs. non-mutating functions）](#原地修改与非原地修改函数（Mutating-vs.-non-
# mutating-functions）)
#     - [一些高阶函数](#一些高阶函数)
#     - [广播（broadcast）](#广播（broadcast）)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 字符串
# ------------------------------------------------------------------------------------------

string1 = "How many cats "

string2 = "is too many cats?"

string(string1, string2)

😺 = 10
println("I don't know but $😺 are too few!")

# ------------------------------------------------------------------------------------------
# 注意：Julia 允许用 Unicode 字符，例如😺
#
# 这就让我们可以把代码写成这个样子
# ------------------------------------------------------------------------------------------

😺 = 1
😀 = 0
😞 = -1

😺 + 😞 == 😀

# ------------------------------------------------------------------------------------------
# ## 数据结构
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 元组（Tuples）
#
# 要创建创建一个元组，就用小括号依次 `( )` 括起一系列元素
#
#
# 语法：
# ```julia
# (item1, item2, ...)
# ```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = ("penguins", "cats", "sugargliders")

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# ###  字典（Dictionaries）
#
# 想存键值对，用字典。通过 `Dict()` 构造器来创建字典。
#
# 语法：
# ```julia
# Dict(key1 => value1, key2 => value2, ...)
# ```
#
# 比如电话簿需要保存人名和电话号码的成对关系
# ------------------------------------------------------------------------------------------

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

myphonebook["Jenny"]

# ------------------------------------------------------------------------------------------
# ### 数组（Arrays）
#
# 不同于元组，数组可变；不同于字典，数组有序。
# 要创建创建一个数组，用方括号依次 `[ ]` 括起一系列元素。
#
# 语法：
# ```julia
# [item1, item2, ...]
# ```
#
# 例如：
# ------------------------------------------------------------------------------------------

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mixture = [1, 1, 2, 3, "Ted", "Robyn"]

# ------------------------------------------------------------------------------------------
# 数组可以包含其他数据结构，也可以包含其他数组变成高维数组。
# ------------------------------------------------------------------------------------------

numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

rand(4, 3)

# ------------------------------------------------------------------------------------------
# ## 循环
#
# ### `for` 循环
#
# `for` 循环语法：
# ```julia
# for <var> in <loop iterable>
#     <loop body>
# end
# ```
# ------------------------------------------------------------------------------------------

for n in 1:10
    println(n)
end

# ------------------------------------------------------------------------------------------
# ### `while` 循环
#
# `while` 循环语法：
# ```julia
# while <condition>
#     <loop body>
# end
# ```
# ------------------------------------------------------------------------------------------

n = 0
while n < 10
    n += 1
    println(n)
end

# ------------------------------------------------------------------------------------------
# ## 条件判断
#
# ### `if` 分支
#
# Julia 中 `if` 语法为：
#
# ```julia
# if <condition 1>
#     <option 1>
# elseif <condition 2>
#     <option 2>
# else
#     <option 3>
# end
# ```
#
# 这让我们可以按条件执行不同的语句
# ------------------------------------------------------------------------------------------

x, y = 1, 2
if x > y
    x
else
    y
end

# ------------------------------------------------------------------------------------------
# ### 三元运算符
#
# 我们可以用三元运算符语法重写上一个代码块
#
# ```julia
# a ? b : c
# ```
#
# 它等价于
#
# ```julia
# if a
#     b
# else
#     c
# end
# ```
# ------------------------------------------------------------------------------------------

(x > y) ? x : y

# ------------------------------------------------------------------------------------------
# ## 函数
#
# 将要谈论的主题：
# 1. 如何声明一个函数
# 2. Julia 中的鸭子类型（Duck-typing）
# 3. 原地修改与非原地修改函数（Mutating vs. non-mutating functions）
# 4. 高阶函数
#
# ### 如何声明一个函数
#
# #### 方法一：用 `function` 和 `end` 关键字
# ------------------------------------------------------------------------------------------

function f(x)
    x^2
end

# ------------------------------------------------------------------------------------------
# #### 方法二：用 `=`
# ------------------------------------------------------------------------------------------

f2(x) = x^2

# ------------------------------------------------------------------------------------------
# #### 方法三：匿名函数
# ------------------------------------------------------------------------------------------

f3 = x -> x^2

# ------------------------------------------------------------------------------------------
# ### 调用函数
# ------------------------------------------------------------------------------------------

f(42)

f2(42)

f3(42)

# ------------------------------------------------------------------------------------------
# ### Julia 中的鸭子类型（Duck-typing）
# > *“假如它嘎嘎像个鸭鸭，那它就是个鸭鸭”*
#
# 在未指定函数参数的类型时，不管输入参数是什么类型，只要函数体可以正常执行，
# 那么 Julia 函数就可以应用在这种类型的参数上
# 比如上面的 `f` 函数可以接受一个矩阵（matrix）作为参数
# ------------------------------------------------------------------------------------------

A = rand(3, 3)

f(A)

# ------------------------------------------------------------------------------------------
# 然而，`f` 函数无法接受一个向量（vector）。
# 不同于 `A^2` 是已经定义好的矩阵乘方，当 `v` 是向量时的 `v^2` 是未定义的，会导致错误。
# ------------------------------------------------------------------------------------------

v = rand(3)

f(v)

# ------------------------------------------------------------------------------------------
# ### 原地修改与非原地修改函数（Mutating vs. non-mutating functions）
#
# 约定：函数名以 `!` 结束的原地修改传入的变量，函数名非以 `!` 结束的则不会改变传入的变量。
#
# 例如下面的 `sort` 和 `sort!`
# ------------------------------------------------------------------------------------------

v = [3, 5, 2]

sort(v)

v

# ------------------------------------------------------------------------------------------
# `sort(v)`会返回一个新的、将 `v` 排序后的数组，`v` 本身不会改变
#
# 而当执行 `sort!(v)` 后，`v` 将会被原地排序修改
# ------------------------------------------------------------------------------------------

sort!(v)

v

# ------------------------------------------------------------------------------------------
# ### 一些高阶函数
#
# #### `map`
#
# `map` 是 Julia 中的高阶函数。“高阶”指的是它会*接受一个函数*作为它的传入参数。
# `map` 将传入的函数作用于传入的数据结构的每一个元素。例如执行
#
# ```julia
# map(f, [1, 2, 3])
# ```
# 将得到一个输出数组，数组的每一个元素为 `f` 作用在数组 `[1, 2, 3]` 的对应元素
# ```julia
# [f(1), f(2), f(3)]
# ```
# ------------------------------------------------------------------------------------------

map(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 这里我们把向量 `[1, 2, 3]` 的每一个元素都平方了，而不是对向量 `[1, 2, 3]` 本身进行平方。
#
# `map` 可以接受匿名函数作为参数，像这样
# ------------------------------------------------------------------------------------------

map(x -> x^3, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 这就完成了对数组 `[1, 2, 3]` 所有元素立方的计算！
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 广播（broadcast）
#
# 广播（`broadcast`）和 `map` 一样也是高阶函数。
# 广播是一种广义化的 `map`，所以 `map` 能做的它都能做。并不止如此，广播的调用方法和 `map` 一样。
# ------------------------------------------------------------------------------------------

broadcast(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 我们又一次通过函数 `f` 把向量 `[1, 2, 3]` 的每一个元素都平方了
# ——这次是用“广播”函数 `f` 的方式！
#
# 广播（`broadcast`）有个语法糖：调用要广播的函数时，在函数名和输入变量之间加一个 `.`。
#
# 比如：
# ```julia
# broadcast(f, [1, 2, 3])
# ```
# 就等价于
# ```julia
# f.([1, 2, 3])
# ```
# ------------------------------------------------------------------------------------------

f.([1, 2, 3])

# ------------------------------------------------------------------------------------------
# 再一次强调 `map` 和广播方式与如下直接调用函数的区别。
#
# ```julia
# f([1, 2, 3])
# ```
# 我们可以求向量中每个元素的平方，但是不能求向量的平方！
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 为了再讲清楚这一点，我们再看一下以下调用的区别
#
# ```julia
# f(A)
# ```
# 与
# ```julia
# f.(A)
# ```
# 其中 `A` 是个矩阵：
# ------------------------------------------------------------------------------------------

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

# ------------------------------------------------------------------------------------------
# 和之前一样，对于矩阵 `A` 来说
# ```
# f(A) = A^2 = A * A
# ```
#
# 另一方面，
# ------------------------------------------------------------------------------------------

B = f.(A)

# ------------------------------------------------------------------------------------------
# 结果包含了矩阵 `A` 中每个元素的平方
#
# 广播的点语法让比较复杂的元素间运算的表示更自然、更贴近数学表示。
#
# 例如：
# ------------------------------------------------------------------------------------------

C = A .+ 2 .* f.(A) ./ A

# ------------------------------------------------------------------------------------------
# 而不是
# ------------------------------------------------------------------------------------------

broadcast(x -> x + 2 * f(x) / x, A)

# ------------------------------------------------------------------------------------------
# 并且这么写仍然能获得和 C 语言一样高效的代码！
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
const NCube{N} = Tuple{Vararg{AbstractHypercubeTransform, N}}

struct TupleHC{T} <: VectorHC
    transformations::T
    dimension::Int
    function TupleHC(transformations::T) where {T<:NCube}
        new{T}(transformations, _sum_dimensions(transformations))
    end
    function TupleHC(transformations::T
                    ) where {N, S <: NCube, T <: NamedTuple{N,S}}
        new{T}(transformations, _sum_dimensions(transformations))
    end
end
ascube(transformations::NCube) = TupleHC(transformations)
ascube(dists::Tuple{Vararg{Union{Dists.Distribution, MT.AbstractMeasure},N}}) where {N}= ascube(ascube.(dists))

ascube(transformations::NamedTuple{N, <:NCube}) where {N} = TupleHC(transformations)

function ascube(dists::NamedTuple{N, <:Tuple{Vararg{Union{Dists.Distribution, MT.AbstractMeasure},M}}}) where {N,M}
    ascube(NamedTuple{N}(ascube.(values(dists))))
end
dimension(tt::TupleHC) = tt.dimension

"""
    $(SIGNATURES)
Helper function that steps through the transformation tuple
"""
function transform_tuple(tt::NCube, x, index)
    return _transform_tuple(x, index, tt)
end

# ends the tuple stepper
_transform_tuple(x::AbstractVector, index, ::Tuple{}) = (), index

# steps throught
function _transform_tuple(x::AbstractVector, index, ts)
    tfirst = first(ts)
    yfirst, index1 = _step_transform(tfirst, x, index)
    yrest, index2  = _transform_tuple(x, index1, Base.tail(ts))
    return (yfirst, yrest...), index2
end

function _step_transform(c::TupleHC{<:Tuple}, x, index)
    transform_tuple(c.transformations, x, index)
end

function _step_transform(c::TupleHC{<:NamedTuple{N}}, x, index) where {N}
    transformations = c.transformations
    y, index′ = transform_tuple(values(transformations), x, index)
    NamedTuple{keys(transformations)}(y), index′
end


text:
function call(::Types.Login, args::Vector, storage, conn::IO)
    ip_addr = string(first(getpeername(conn)))
    name = args[1]

    if (name in [value.name for (_, value) in storage.active_clients])
        return (write(conn, "This account is already used\n"))
    end

    arr = Requests.get_account(storage.sql_conn, name)

    if (length(arr) == 0)
        return (write(conn, "Invalid username or password\n"))
    end
    password = String(args[2])
    arr = map(x -> string(x), arr)
    if (Bcrypt.CompareHashAndPassword(arr[2], password) == false)
        return (write(conn, "Invalid username or password\n"))
    end

    storage.active_clients[ip_addr].id = arr[1]
    storage.active_clients[ip_addr].name = name
    storage.active_clients[ip_addr].is_auth = true
    storage.active_clients[ip_addr].current_channel_id = "1"
end


text:
using KnuthBendix
using Test

include("abstract_words.jl")

@testset "KnuthBendix.jl" begin
   include("words.jl")
   include("bufferwords.jl")
   include("alphabets.jl")
   include("orderings.jl")
   include("rewriting.jl")
   include("automata.jl")
   include("kbs1.jl")
   include("kbs.jl")

   include("kbmag_parsing.jl")

   include("rws_examples.jl")
   include("test_examples.jl")
end


text:
using Test
using Nystrom
using StaticArrays
using LinearAlgebra

@testset "Kernels" begin
    pde  = Helmholtz(;dim=3,k=1)
    G    = SingleLayerKernel(pde)
    dG    = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == ComplexF64
    @test Nystrom.return_type(dG) == ComplexF64
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Laplace(;dim=3)
    G    = SingleLayerKernel(pde)
    dG    = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == Float64
    @test Nystrom.return_type(dG) == Float64
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Elastostatic(;μ=1,λ=2.0,dim=3)
    G    = SingleLayerKernel(pde)
    dG   = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == SMatrix{3,3,Float64,9}
    @test Nystrom.return_type(dG) == SMatrix{3,3,Float64,9}
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)

    pde  = Maxwell(;k=1.0)
    G    = SingleLayerKernel(pde)
    dG   = DoubleLayerKernel(pde)
    @test Nystrom.kernel_type(G) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(dG) == Nystrom.DoubleLayer()
    @test Nystrom.return_type(G) == SMatrix{3,3,ComplexF64,9}
    @test Nystrom.return_type(dG) == SMatrix{3,3,ComplexF64,9}
    @test Nystrom.combined_field_coefficients(G) == (0,-1)
    @test Nystrom.combined_field_coefficients(dG) == (1,0)
end

@testset "Maxwell identities" begin
    # essentially vector calculus identities verified by Maxwell's diadic greens
    # function
    op = Maxwell(2)
    T  = Nystrom.default_density_eltype(op)
    G   = SingleLayerKernel(op)
    dG  = DoubleLayerKernel(op)
    E  = rand(T)
    x  = rand(SVector{3})
    y  = rand(SVector{3})
    jac = rand(SMatrix{3,2})
    xdof = NystromDOF(x,1.,jac,0,0)
    ydof = NystromDOF(y,1.,jac,0,0)
    # [ n(y) × (∇ʸ × G(x,y) ) ]ᵗ E = [ ∇ʸ × G ]ᵗ (E × n(y)) = -[ ∇ˣ G ] (n(y) × E)
    lhs = dG(xdof,ydof)*E
    rhs = Nystrom._curl_y_green_tensor_maxwell(x,y,op.k) * cross(E,normal(ydof))
    @test lhs ≈ rhs
    # [ ∇ʸ × G ]ᵗ = ∇ˣ × G
    lhs =  Nystrom._curl_y_green_tensor_maxwell(x,y,op.k) |> transpose
    rhs =  Nystrom._curl_x_green_tensor_maxwell(x,y,op.k)
    @test lhs == rhs
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

module FEMBasis

import Calculus
using LinearAlgebra
using Tensors
# Reexport Vec
export Vec

const Vecish{N, T} = Union{NTuple{N, T}, Vec{N, T}}

abstract type AbstractBasis{dim} end
# Forward methods on instances to types
Base.length(B::T) where T<:AbstractBasis = length(T)
Base.size(B::T) where T<:AbstractBasis = size(T)
eval_basis!(B::T, N, xi) where T<:AbstractBasis = eval_basis!(T, N, xi)
eval_dbasis!(B::T, dN, xi) where T<:AbstractBasis = eval_dbasis!(T, dN, xi)
eval_basis(B::AbstractBasis{dim}, xi) where {dim} = eval_basis!(B, zeros(length(B)), xi)
eval_dbasis(B::AbstractBasis{dim}, xi) where {dim} = eval_dbasis!(B, zeros(Vec{dim}, length(B)), xi)

include("subs.jl")
include("vandermonde.jl")

include("create_basis.jl")
export get_reference_element_coordinates, eval_basis!, eval_dbasis!

include("lagrange_segments.jl")
export Seg2, Seg3
include("lagrange_quadrangles.jl")
export Quad4, Quad8, Quad9
include("lagrange_triangles.jl")
export Tri3, Tri6, Tri7
include("lagrange_tetrahedrons.jl")
export Tet4, Tet10
include("lagrange_hexahedrons.jl")
export Hex8, Hex20, Hex27
include("lagrange_wedges.jl")
export Wedge6, Wedge15
include("lagrange_pyramids.jl")
export Pyr5

include("nurbs.jl")
include("nurbs_segment.jl")
export NSeg
include("nurbs_surface.jl")
export NSurf
include("nurbs_solid.jl")
export NSolid
include("math.jl")
export interpolate, interpolate!, jacobian
export grad, grad!
export BasisInfo

end


text:
## Exercise 4-10
## Write an appropriately general set of functions that can draw shapes as in Turtle pies.
## https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#fig04-3
println("Ans: ")

println("End.")