text:
using ImportMacros

@import ChargedParticleDynamics.GuidingCenter4d.TokamakFastDeeplyTrapped as problem

# run configuration
const Δt    = 5.
const ntime = 1250000

# plot configuration
const nplot = 1
const splot = true
const XLIM  = ( 2.10,  3.00)
const YLIM  = (-0.80, +0.80)

# output configuration
const coll_title = "Guiding Center 4D Fast Deeply Trapped"


text:
# Function for calculating the sample autocovariance of a vector for lags from
# 0 to maxlag, returning a vector of length maxlag+1
function autocov(sample::Vector{Float64}, maxLag::Int)
  n = size(sample, 1)

  if n <= maxLag
    throw("The length of the input vector must be at least maxlag+1")
  end

  sample = sample-mean(sample)

  acv = Array(Float64, maxLag+1)

  for i = 0:maxLag
     acv[i+1] = dot(sample[1:n-i], sample[1+i:n])
  end

  acv = acv/n

  return acv
end

# Function for calculating the standard Normal cumulative density function more
# accurately for small values of its argument
function logOfNormalCdf(x::Union(Float64, Vector{Float64}))

  c::Float64 = -6.5
  a = [1, -1, 3, -15, 105, -945, 10395, -135135, 2027025, -34459425, 654729075]

  n = size(x, 1)
  y = Array(Float64, n)

  for i in 1:n
    if x[i]>c
      y[i] = log(cdf(Normal(0., 1.), x[i]))
    else
      y[i] = log(-dot(a, x[i].^(-2*(1:11)+1)))-x[i]^2/2-log(2*pi)/2
    end
  end

  return y
end

# Function for estimating the variance of a single MCMC chain using the initial
# monotone sequence estimator of Geyer, see
# Practical Markov Chain Monte Carlo, Charles J. Geyer, Statistical Science,
# Vol. 7, No. 4. (1992), pp. 473-483, doi:10.2307/2246094
function imse(sample::Union(Array{Float64, 1}, Array{Float64, 2}), maxLag::Int)
  nData = size(sample, 1)
  nPars = size(sample, 2)

  k = convert(Int, floor((maxLag-1)/2))

  # Preallocate memory
  acv = Array(Float64, maxLag+1, nPars)
  g = Array(Float64, k+1, nPars)
  m = (k+1)*ones(nPars)
  variance = Array(Float64, nPars)

  # Calculate empirical autocovariance
  for i = 1:nPars
    acv[:, i] = autocov(sample[:, i], maxLag)
  end

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for i = 1:nPars
    for j = 0:k
      g[j+1, i] = acv[2*j+1, i]+acv[2*j+2, i]
      if g[j+1, i] <= 0
        m[i] = j
        break
      end
    end
  end

  # Create the monotone sequence of g
  for i = 1:nPars
    if m[i] > 1
      for j = 2:m[i]
        if g[j, i] > g[j-1, i]
          g[j, i] = g[j-1, i]
        end
      end
    end
  end

  # Calculate the initial monotone sequence estimator
  for i = 1:nPars
    variance[i] = (-acv[1, i]+2*sum(g[1:m[i], i]))/nData
  end

  return variance
end

ext: jl


text:
struct ForgeAPIError <: GitForge.ForgeError
    message::String
end

struct ForgeAPINotImplemented <: GitForge.ForgeError
    func::Function
    args::Tuple
    noted::Bool
    api
end

Base.showerror(io::IO, e::ForgeAPINotImplemented) =
    print(io, typeof(e.api), " has not implemented the function $(e.func) for arguments $(typeof.((e.api, e.func, e.args...)))\n  use: @not_implemented(::$(typeof(e.api)), ::typeof($(e.func)), $(join(map(a-> "::$(typeof(a))", e.args), ", ")))")

not_implemented(api, func, args...; noted = true) = throw(ForgeAPINotImplemented(func, args, noted, api))

"""
A forge is an online platform for Git repositories.
The most common example is [GitHub](https://github.com).

`Forge` subtypes can access their respective web APIs.
"""
abstract type Forge end

"""
    Endpoint(
        method::Symbol,
        url::$AStr;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::Dict=Dict(),
        allow_404::Bool=false,
    ) -> Endpoint

Contains information on how to call an endpoint.

## Arguments
- `method::Symbol`: HTTP request method to use.
- `url::$AStr`: Endpoint URL, relative to the base URL.

## Keywords
- `headers::Vector{<:Pair}=HTTP.Header[]`: Request headers to add.
- `query::Dict=Dict()`: Query string parameters to add.
- `allow_404::Bool=false`: Sends responses  with 404 statuses to the postprocessor.
"""
struct Endpoint
    method::Symbol
    url::String
    headers::Vector{<:Pair}
    query::Dict
    allow_404::Bool

    function Endpoint(
        method::Symbol, url::AStr;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::Dict=Dict(),
        allow_404::Bool=false,
    )
        return new(method, url, headers, query, allow_404)
    end
end

"""
    base_url(::Forge) -> String

Returns the base URL of all API endpoints.
"""
base_url(::Forge) = ""

"""
    request_headers(::Forge, ::Function) -> Vector{Pair}

Returns the headers that should be added to each request.
"""
request_headers(::Forge, ::Function) = HEADERS

"""
    request_query(::Forge, ::Function) -> Dict

Returns the query string parameters that should be added to each request.
"""
request_query(::Forge, ::Function) = Dict()

"""
    request_kwargs(::Forge, ::Function) -> Dict{Symbol}

Returns the extra keyword arguments that should be passed to `HTTP.request`.
"""
request_kwargs(::Forge, ::Function) = Dict()

"""
    postprocessor(::Forge, ::Function) -> PostProcessor

Returns the [`PostProcessor`](@ref) to be used.
"""
postprocessor(::Forge, ::Function) = DoNothing()

"""
    into(::Forge, ::Function) -> Type

Returns the type that the [`PostProcessor`](@ref) should create from the response.

"""
into(::Forge, ::Function) = Any

"""
    endpoint(::Forge, ::Function, args...) -> Endpoint

Returns an [`Endpoint`](@ref) for a given function.
Trailing arguments are usually important for routing.
For example, [`get_user`](@ref) can take some ID parameter which becomes part of the URL.
"""
endpoint(forge::T, func::Function, args...) where T <: Forge =
    throw(ForgeAPINotImplemented(func, args, false, forge))

macro not_implemented(api::Expr, func::Expr, rest...)
    api_name = length(api.args) == 2 ? api.args[1] : :api
    func_name = func.args[1].args[end]
    rest = [rest[n].head !== :(::) ? rest[n] : :($(Symbol("a$n"))::$(rest[n].args[end])) for n in 1:length(rest)]
    rest_names = map(a-> a.head !== :(::) ? a : a.args[1], rest)
    :(endpoint(api::$(esc((api.args[end]))), $(esc(func)), $(rest...)) = throw(ForgeAPINotImplemented($func_name, ($(rest_names...),), true, $api_name)))
end

"""
    has_rate_limits(::Forge, ::Function) -> Bool

Returns whether or not the forge server uses rate limiting.
"""
has_rate_limits(::Forge, ::Function) = false

"""
    rate_limit_check(::Forge, ::Function) -> Bool

Returns whether or not there is an active rate limit.
If one is found, [`on_rate_limit`](@ref) is called to determine how to react.
"""
rate_limit_check(::Forge, ::Function) = false

"""
    on_rate_limit(::Forge, ::Function) -> OnRateLimit

Returns an [`OnRateLimit`](@ref) that determines how to react to an exceeded rate limit.
"""
on_rate_limit(::Forge, ::Function) = ORL_THROW

"""
    rate_limit_wait(::Forge, ::Function)

Wait for a rate limit to expire.
"""
rate_limit_wait(::Forge, ::Function) = nothing

"""
    rate_limit_period(::Forge, ::Function) -> Period

Compute the amount of time until a rate limit expires.
"""
rate_limit_period(::Forge, ::Function) = nothing

"""
    rate_limit_update!(::Forge, ::Function, ::HTTP.Response)

Update the rate limiter with a new response.
"""
rate_limit_update!(::Forge, ::Function, ::HTTP.Response) = nothing


text:
module TcrCore

##############################################################################
##
## Dependencies
##
##############################################################################

using Distributions, HypothesisTests, StatsBase

##############################################################################
##
## Exported methods and types (in addition to everything reexported above)
##
##############################################################################

export
    #Agents:
    Agent,
    setupAgentsWithFixedAccuracy,
    setupRandomAgents,
    evaluateCandidateByAgent,
    evaluateCandidateByAgentPrecise,

    #Items:
    getCandidate,
    registryMean,

    #Actions:
    simpleVote,
    binaryTokenVote,
    proRataTokenVote,
    selectRandomAgent,
    selectRandomAgentWithMinBalance,
    noneRedistribution,
    onlyChallengerRewardRedistribution,
    fullRedistribution,
    challenge,
    application,

    #Simulations:
    simSimple,
    simChallenge,
    simToken

    #Research
    # tokensAndDiversifications,
    # tokensAndDiversificationsEffects

##############################################################################
##
## Load files
##
##############################################################################

include("environment/items.jl")
include("environment/agents.jl")
include("actions/vote.jl")
include("actions/select-agent.jl")
include("actions/redistribute.jl")
include("actions/composite-actions.jl")
include("simulations/simulations.jl")
include("experiments/benchmarks.jl")
include("experiments/sample-scenarios.jl")
include("experiments/research.jl")

end # module TcrSimulations

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
using SQLite

"""
    select_columns(db, table)
Return an array with names of columns in the given table
"""
function select_columns(db::SQLite.DB, table)
    cols_query = SQLite.query(db, "SHOW COLUMNS FROM $table;")
    cols_query[2]
end

"""
    select_all_tables_mysql(con)

Return an array of all tables in a given MySQL database
"""
function select_all_tables(db::SQLite.DB)
    tables_query = SQLite.query(db, "SELECT name FROM sqlite_master WHERE type='table'")
    tables_query[1].values
end

"""
    query(db, query_code)

Execute a SQLite command
"""
function db_query(db::SQLite.DB, query_code)
    sel = SQLite.query(db, query_code)
    return sel
end

"""
    insert_row!(db, tablename, values)
Insert a row of values into the specified table for a given a SQLite database handle

### Arguments:

* `db::MySQLDB`: Database object (connection and map)
* `data_values::Dict{ASCIIString, Any}`: Array of (string) values
* `verbose`: Print debugginh info
"""
function insert_row!{T}(db::SQLite.DB, tablename, data_values::Dict{Symbol, T},
    verbose = false)

    cols_string, vals_string = assemble_cols_and_vals(data_values)
    lastid = -1
    try
        q = db_query(db, "INSERT INTO `$tablename` ($cols_string) values $vals_string")
    catch e
        if verbose
            Base.showerror(STDOUT, e)
            println("\n")
            println("Warning! Could not insert values $vals_string into table $tablename")
        end
        return -1
    end

    lastid_query = db_query(db, "SELECT last_insert_rowid()")
    lastid = get(lastid_query[1][1], -1)
    if verbose
        println("Row successfully inserted into table: $tablename")
    end

    return lastid
end


text:
using UUIDs,Plots,JLD2

function sin_generator(freq, t)
    sin(2*pi*t*freq)
end

function draw_graph(time_dat_path::String,
                    w_c_dat_path::String,
                    w_s_dat_path::String,
                    w_fm_dat_path::String)

    timeline::Array{Float64} = []
    w_fm_vals::Array{Float64} = []
    w_c_vals::Array{Float64} = []
    w_s_vals::Array{Float64} = []

    open(time_dat_path,"r") do f_t
        open(w_c_dat_path,"r") do f_c
            open(w_s_dat_path,"r") do f_s
                open(w_fm_dat_path,"r") do f_fm

                    #plot_graph=plot(timeline, w_s_vals, xlims=(0.0,1/1000), size=(2000,480), fmt =:png)
                    #plot!(plot_graph, timeline, w_fm_vals)
                    #plot!(plot_graph, timeline, w_c_vals)

                    #plot_graph=plot(timeline, w_fm_vals, xlims=(0.0003,0.0003+0.000005), size=(2000,700), fmt =:png)
                    #plot!(plot_graph, timeline, w_s_vals)
                    #plot!(plot_graph, timeline, w_c_vals)
                    #savefig(plot_graph,"fm_signal.png")

                    #plot_graph=plot(timeline, w_fm_vals, xlims=(0.0007,0.0007+0.000005), size=(2000,700), fmt =:png)
                    #plot!(plot_graph, timeline, w_s_vals)
                    #plot!(plot_graph, timeline, w_c_vals)
                    #savefig(plot_graph,"fm_signal_2.png")

                    #plot_graph=plot(timeline, w_fm_vals, xlims=(0.0005,0.0005+0.000005), size=(2000,700), fmt =:png)
                    #plot!(plot_graph, timeline, w_s_vals)
                    #plot!(plot_graph, timeline, w_c_vals)
                    #savefig(plot_graph,"fm_signal_3.png")
                end
            end
        end
    end
end

function main()

    uuid = uuid4()

    println(string(uuid))

    mkdir(string(uuid))

    time_dat_path = string(uuid) * "/time.dat"
    w_c_dat_path  = string(uuid) * "/w_c.dat"
    w_s_dat_path  = string(uuid) * "/w_s.dat"
    w_fm_dat_path = string(uuid) * "/w_fm.dat"

    # f_s_max = 22KHz
    # channel spacing 200kHz
    # frequnency deviation delta_f=75kHz
    # should be carrier_freq > delta_f
    carrier_freq  =      7780*1000 #7.78MHz
    sampling_rate =  100*1000*1000
    Δt = 1/sampling_rate
    Δf = 75*1000
    accumulate = 0

    f_t  = open(time_dat_path, "w")
    f_c  = open(w_c_dat_path,  "w")
    f_s  = open(w_s_dat_path,  "w")
    f_fm = open(w_fm_dat_path, "w")

    for step in 1:100*1000*1000
        t = step/sampling_rate
        w_s = sin_generator(1000, t)
        accumulate += w_s*Δt
        #println(2*pi*Δf*accumulate)
        w_c = sin(2*pi*carrier_freq*t)
        w_fm = sin(2*pi*carrier_freq*t+2*pi*Δf*accumulate)
        write(f_t, t)
        write(f_c, w_c)
        write(f_s, w_s)
        write(f_fm, w_fm)
    end

    close(f_fm)
    close(f_s)
    close(f_c)
    close(f_t)

    jldopen(ARGS[1],"w") do file
        file["samplingrate"] = sampling_rate

        vals = []
        open(time_dat_path,"r") do f
            while !eof(f_t)
                push!(vals, read(f, Float64))
            end
        end
        file["times"] = vals

        vals = []
        open(w_s_dat_path,"r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_s"] = vals

        vals = []
        open(w_c_dat_path,"r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_c"] = vals

        vals = []
        open(w_fm_dat_path,"r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_fm"] = vals
    end

    rm(string(uuid), force=true, recursive=true)
end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
# This file was generated, do not modify it. # hide
leaf_range = range(lgb, :min_data_in_leaf, lower=1, upper=50)

curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=leaf_range, resolution=50,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Min data in leaf", fontsize=14)
ylabel("RMSE", fontsize=14)

plt.savefig(joinpath(@OUTPUT, "lgbm_hp3.svg")) # hide


text:
using SpectralDistances, DynamicAxisWarping, DSP, Distances

nfft = 32
T = 200
q = sin.(1:T)
y = [zeros(10nfft);q;zeros(10nfft)] .+ 0.0001 .* randn.()

Q = spectrogram(q,nfft)
Y = spectrogram(y,nfft)

n, m = size(Q.power)
rad = 5

## Discrete grid transport
res = dtwnn(
    log1p.(Q.power),
    log1p.(Y.power),
    DiscreteGridTransportDistance(Cityblock(), eltype(Q.power), n, n),
    rad,
    saveall = true,
)

@test res.cost < 0.01
@test res.loc == 21
@test all(isfinite, res.dists)
isinteractive() && plot(res.dists)

fm  = TimeWindow(inner = LS(na = 2, λ=1e-3), n = nfft, noverlap = nfft÷2)

Qm  = fm(q) |> change_precision(Float64)
Ym  = fm(y) |> change_precision(Float64);

res = dtwnn(
    Qm,
    Ym,
    EuclideanRootDistance(p = 1, weight = simplex_residueweight),
    rad,
    saveall = true,
)
@test res.cost < 0.01
@test  20 <= res.loc <= 22
@test all(isfinite, res.dists)
isinteractive() && plot(res.dists)
##
res_dtwotrd = dtwnn(
    Qm,
    Ym,
    OptimalTransportRootDistance(p = 1, weight = simplex_residueweight),
    rad,
    saveall = true,
    tol = 1e-3,
)
@test res.cost < 0.01
@test  20 <= res.loc <= 22
@test all(isfinite, res.dists)
isinteractive() && plot(res.dists)

dist = TimeDistance(
    inner = OptimalTransportRootDistance(p = 1, β = 0.5, weight = simplex_residueweight),
    tp    = 1,
    c     = 0.1,
)
res = distance_profile(dist, Qm, Ym, tol=1e-3, check_interval=10)
@test minimum(res) < 1
@test  20 <= argmin(res) <= 22
@test all(isfinite, res)
isinteractive() && plot(res)


text:
##
## Using RadiativeTransfer.Scattering to perform Mie computations
##

using Revise
using RadiativeTransfer
using RadiativeTransfer.Scattering

# ? Should we make a wrapper around LogNormal?
using Distributions
using FastGaussQuadrature
# using BenchmarkTools

##
## STEP 1: Create the Aerosol
##

# Aerosol particle distribution and properties
μ  = (π/3) * 0.55 /(2π)                # Log mean radius
σ  = 1.1               # Log stddev of radius
nᵣ = 1.5                # Real part of refractive index
nᵢ = 0.0                # Imag part of refractive index
r_max = 30.0            # Maximum radius
nquad_radius = 2500     # Number of quadrature points for integrating of size dist.

size_distribution = make_log_normal_size_dist(μ, σ)

# Create the aerosol
aero = Aerosol(size_distribution, nᵣ, nᵢ)

##
## STEP 2: Create the Mie Calculations model
##

λ = 0.55   # Incident wavelength
polarization_type = Stokes_I()
truncation_type = δBGE(10, 10)
wigner_file_path = "/home/rjeyaram/RadiativeTransfer/src/Scattering/Mie/wigner_values.jld"

model_NAI2 = make_mie_model(NAI2(), aero, λ, polarization_type, truncation_type, r_max, nquad_radius)
model_PCW = make_mie_model(PCW(), aero, λ, polarization_type, truncation_type, r_max, nquad_radius, wigner_file_path)

##
## STEP 3: Perform the Mie Calculations
##

# @btime aerosol_optics_NAI2 = compute_aerosol_optical_properties(model_NAI2);
# @btime aerosol_optics_PCW = compute_aerosol_optical_properties(model_PCW);

# Compute the derivatives!
aerosol_optics_NAI2 = compute_aerosol_optical_properties(model_NAI2, autodiff=false);
aerosol_optics_PCW = compute_aerosol_optical_properties(model_PCW);

# Compute scattering matrix from optical properties
μ_quad, w_μ = gausslegendre(1000)
scattering_matrix = Scattering.reconstruct_phase(aerosol_optics_NAI2.greek_coefs, μ_quad);

# aerosol_optics_NAI2.greek_coefs.α ≈ aerosol_optics_PCW.greek_coefs.α

# aerosol_optics_NAI2.greek_coefs ≈ aerosol_optics_PCW.greek_coefs
aerosol_optics_NAI2 ≈ aerosol_optics_PCW


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Guile"
version = v"2.2.4"

# Collection of sources required to build Guile
sources = [
    "https://ftp.gnu.org/gnu/guile/guile-2.2.4.tar.gz" =>
    "33b904c0bf4e48e156f3fb1d0e6b0392033bd610c6c9d9a0410c6e0ea96a3e5c",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd guile-2.2.4/
./configure LDFLAGS=-L$prefix/lib CFLAGS=-I$prefix/include --prefix=$prefix --host=$targt
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libguile", :libguile),
    ExecutableProduct(prefix, "guile", :guile)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libiconv-v1.15-0/build_Libiconv.v1.15.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Gettext-v0.19.8-0/build_Gettext.v0.19.8.jl",
    "https://github.com/yukota/LibtoolBuilder/releases/download/v2.4.6-1/build_Libtool.v2.4.6.jl",
    "https://github.com/yukota/libunistringBuilder/releases/download/v0.9.10-2/build_libunistring.v0.9.10.jl",
    "https://github.com/yukota/libgcBuilder/releases/download/v8.0.2-0/build_libgc.v8.0.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libffi-v3.2.1-0/build_Libffi.v3.2.1.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: yukota/GuileBuilder

max_issues_repo_name: yukota/GuileBuilder

max_forks_repo_name: yukota/GuileBuilder


text:
using Reexport

@reexport module IntegrationDynamicalSystems
    import ..CausalityTests
    import CausalityToolsBase: causality
    import CausalityToolsBase: CausalityTest
    import CausalityToolsBase
    import DynamicalSystems: DynamicalSystem, DiscreteDynamicalSystem, trajectory, ContinuousDynamicalSystem

    ############################
    # Dynamical system resampling
    ############################
    include("DynamicalSystemResampling.jl")

    ############################
    # Dynamical system setup
    ############################
    include("DynamicalSystemSetup.jl")
    include("DiscreteSystemSetup.jl")

    ##########################
    # Continuous system setup
    ##########################
    include("ContinuousSystemSetup.jl")

    ############################################
    # causality tests for discrete systems
    ############################################

    """
        causality(system::DiscreteDynamicalSystem, setup::DiscreteSystemSetup,
            test::CausalityTest)

    Apply the causality `test` to the given discrete `system` using the provided `setup` parameters.

    ## Example

    ```julia
    # Use one of the built-in dynamical systems that ship with CausalityTools,
    # `logistic2_unidir`, which is a system of two chaotic logistic maps `x` and `y` with
    # coupling `x` to `y`, with `c_xy` controlling the coupling strength.
    sys = logistic2_unidir(c_xy = 1.0)

    # We treat the first variable as the source and the second variable as the target.
    # Iterate the system until we have time series with 200 observations.
    setup = DiscreteSystemSetup(source = 1, target = 2, n_pts = 200)

    # Define some causality tests (with default values, you'd want to tweak the
    # parameters to your needs)
    test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
    test_cm = CrossMappingTest(n_reps = 10)
    test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
    test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
    test_jdd = JointDistanceDistributionTest()
    test_jddt = JointDistanceDistributionTTest()
    predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
    test_pa = PredictiveAsymmetryTest(predictive_test = predtest)

    # Compute the causality statistics using the different tests
    causality(sys, setup, test_ccm)
    causality(sys, setup, test_cm)
    causality(sys, setup, test_vf)
    causality(sys, setup, test_tog)
    causality(sys, setup, test_jdd)
    causality(sys, setup, test_jddt)
    ```
    """
    function causality(system::DiscreteDynamicalSystem, setup::DiscreteSystemSetup{NoResampling},
            test::CausalityTest)

        # Iterate the system and create orbit
        dt, Ttr, n_pts = setup.dt, setup.Ttr, setup.n_pts
        orbit = trajectory(system, n_pts, dt = dt, Ttr = Ttr)

        # Causality analysis
        source = orbit[:, setup.source]
        target = orbit[:, setup.target]

        # Resampling scheme is `NoResampling`, so call causality without any constraints.
        causality(source, target, test)
    end

    """
        causality(system::ContinuousDynamicalSystem, setup::ContinuousSystemSetup,
            test::CausalityTest)

    Apply the causality `test` to the given continuous `system` using the provided `setup` parameters.

    ## Example

    Let's re-create the `lorenz_lorenz` system that ships with `CausalityTools` and analyse the
    dynamical influence between some of its components.

    ```julia
    function eom_lorenz_lorenz(u, p, t)
        c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
        x1, x2, x3, y1, y2, y3 = (u...,)

        dx1 = -a₁*(x1 - x2) + c_yx*(y1 - x1)
        dx2 = -x1*x3 + a₂*x1 - x2
        dx3 = x1*x2 - a₃*x3
        dy1 = -b₁*(y1 - y2) + c_xy*(x1 - y1)
        dy2 = -y1*y3 + b₂*y1 - y2
        dy3 = y1*y2 - b₃*y3

        return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
    end

    function lorenz_lorenz(; u0 = rand(6),
            c_xy = 0.2, c_yx = 0.0,
            a₁ = 10, a₂ = 28, a₃ = 8/3,
            b₁ = 10, b₂ = 28, b₃ = 9/3)
        ContinuousDynamicalSystem(eom_lorenz_lorenz, u0, [c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃])
    end

    # Create an instance of the system with default parameters
    sys = lorenz_lorenz()

    # Analysis setup. We'll check for an influence from variable x1 to y1, which are the
    # 1st and 4st components of the orbit. We'll generate orbits consisting of 500 points.
    setup = ContinuousSystemSetup(source = 1, target = 4, n_pts = 500)

    # Predictive asymmetry causality test
    predtest = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = -5:5)
    test_pa = PredictiveAsymmetryTest(predictive_test = predtest)

    causality(sys, setup, test_pa)
    ```
    """
    function causality(system::ContinuousDynamicalSystem, setup::ContinuousSystemSetup{NoResampling},
            test::CausalityTest)

        dt, Ttr, n_pts, sample_step = setup.dt, setup.Ttr, setup.n_pts, setup.sample_step

        # the system is recorded at times t0:dt:T
        T = n_pts*dt*sample_step

        orbit = trajectory(system, T, dt = dt, Ttr = Ttr*dt)[1:sample_step:end-1, :]
        #(percent_noise > 0) ? add_noise(o, percent_noise) : o

        # Causality analysis
        source = orbit[:, setup.source]
        target = orbit[:, setup.target]

        # Resampling scheme is `NoResampling`, so call causality without any constraints.
        causality(source, target, test)
    end

    export causality, ContinuousSystemSetup, DiscreteSystemSetup
end

"""
    IntegrationDynamicalSystems

Module providing an interface for causality detection for dynamical systems as defined
in DynamicalSystems.jl.
"""
IntegrationDynamicalSystems


text:
using LinearAlgebra
using SparseArrays

using FastGaussQuadrature

"Construct the Element-To-Element (EToE) and
Element-To-Face (EtoF) maps
"
function Connect1D(EToV)
    K = size(EToV)[1] # Number of cells
    Nfaces = 2

    TotalFaces = Nfaces*K  # Total number of faces
    Nv = K+1  # Number of vertices

    # Local face to local vertex connections
    vn = [1,2]

    FToVrows = 1:TotalFaces
    FToVcols = zeros(TotalFaces)
    FToVvals = ones(Int, TotalFaces)

    row = 1
    for k=1:K
        for face=1:Nfaces
            FToVcols[row] = EToV[k,vn[face]]
            row += 1
        end
    end
    FToV = sparse(FToVrows, FToVcols, FToVvals)

    # Global face to global face sparse array
    FToF = FToV * FToV' - I

    faces1, faces2, _ = findnz(FToF) # Find non-zero elements

    # Convert face global number to element and face numbers
    element1 = floor.(Int, (faces1 .- 1)/Nfaces) .+ 1
    face1 = mod.( faces1 .- 1, Nfaces ) .+ 1
    element2 = floor.(Int, (faces2 .- 1)/Nfaces) .+ 1
    face2 = mod.( faces2 .- 1, Nfaces ) .+ 1

    # Rearrange into Nelements x Nfaces sized arrays

    # This initialises all faces to connect to themselves
    # Faces which connect to other cells are then overwritten
    EToE = (1:K) * ones(Int, 1, Nfaces)
    EToF = ones(Int, K,1)*(1:Nfaces)'

    # this index code doesn't work as expected
    #ind = sub2ind([K,Nfaces], element1, face1)
    #EToE[ind] = element2
    #EToF[ind] = face2

    for i = 1:length(element1)
        EToE[element1[i], face1[i]] = element2[i]
        EToF[element1[i], face1[i]] = face2[i]
    end

    return EToE, EToF
end

"
K    Number of cells
Nfp  Number of nodes per face
Np   Number of nodes per cell
"
function BuildMaps1D(K::Int, Nfaces::Int, Nfp::Int, Np::Int, Fmask, EToE, EToF, x)

    nodeids = reshape(1:K*Np, Np, K)

    vmapM   = zeros(Int, Nfp, Nfaces, K)
    vmapP   = zeros(Int, Nfp, Nfaces, K)

    for k1=1:K
        for f1=1:Nfaces
            # find index of face nodes with respect to volume node ordering
            vmapM[:,f1,k1] .= nodeids[Fmask[f1], k1]
        end
    end

    for k1=1:K
        for f1=1:Nfaces
            # find neighbor
            k2 = EToE[k1,f1]
            f2 = EToF[k1,f1]

            # find volume node numbers of left and right nodes
            vidM = vmapM[:,f1,k1]
            vidP = vmapM[:,f2,k2]

            x1  = x[vidM]
            x2  = x[vidP]

            # Compute distance matrix
            D = sum(abs2, x1 -x2)
            if D < 1e-5
                vmapP[:,f1,k1] = vidP
            end
        end
    end

    # Find where vmapP == vmapM
    mapB = findall(iszero, vmapP - vmapM)
    vmapB = vmapM[mapB]

    # Create specific left (inflow) and right (outflow) maps
    mapI = 1
    mapO = K*Nfaces
    vmapI = 1
    vmapO = K*Np;

    return vmapP, vmapM, mapB, vmapB, mapI, mapO, vmapI, vmapO
end

# Compute LIFT
function Lift1D(Np, Nfaces, Nfp, V)
    print(Np,", ", Nfaces, ", ", Nfp)
    Emat = zeros(Np, Nfaces*Nfp)

    # Define Emat
    Emat[1,1] = 1.0
    Emat[Np,2] = 1.0

    # inv(mass matrix)*\s_n (L_i,L_j)_{edge_n}
    V*(V'*Emat);
end

struct Quadrature1D
    # locations in reference cell [-1,1]
    locations::Array{Float64, 1}
    weights::Array{Float64, 1}
end

"Compute LGL quadrature points on reference cell [-1,1]"
function gaussLobatto1D(npoints)
    r, w = gausslobatto(npoints)
    Quadrature1D(r, w)
end

function equidistant1D(npoints)
    r = range(-1, stop=1, length=npoints)
    w = fill(2.0 / npoints, npoints)
    Quadrature1D(r, w)
end

"Describes the basis functions and quadrature nodes on a reference cell.
The matrices derived on this cell can be applied to all cells
by including the geometrical scaling factors."
struct ReferenceCell
    basis
    quadrature::Quadrature1D

    V    # Vandermonde matrix on reference cell
    invV # Inverse Vandermonde matrix
    Dr   # Differentiation matrix

    function ReferenceCell(basis, quadrature)
        V = zeros(length(quadrature.locations), basis.order)
        Dr = zeros(length(quadrature.locations), basis.order)

        # Calculate basis function at each quadrature location
        # (Vandermonde matrix)
        #   V[i,j] = P_j(x_i)
        #
        # DVr is the differentiation matrix
        #
        for j=1:basis.order
            V[:,j] = basisFunction(basis, quadrature.locations, j)
            Dr[:,j] = basisGradient(basis, quadrature.locations, j)
        end
        Dr /= V

        new(basis, quadrature,
            V, inv(V), Dr)
    end
end

"Describes a 1D grid"
mutable struct Grid1D

    mesh::Mesh1D # Defines node locations and cells

    npoints::Int # Number of points per cell, including edges (Np)
    npoints_per_face::Int # Number of points on a face
    nfaces::Int # Number of faces per cell

    reference::ReferenceCell

    LIFT # Surface integral terms

    J    # Jacobian [point, cell]
    rx   # Inverse Jacobian

    x::Array{Float64, 2} # Physical coordinates of the grid points [point,cell]

    Fscale

    nx  # Outward pointing normal at element faces

    vmapP # Map of external points for each cell face
    vmapM # Map of internal points for each cell face

    mapB
    vmapB

    mapI
    mapO
    vmapI
    vmapO

    function Grid1D(mesh::Mesh1D, basis, quadrature::Quadrature1D)

        # Create an instance. Not yet initialised
        # This is to avoid a long constructor parameter list
        g = new()

        g.mesh = mesh

        g.npoints = length(quadrature.locations)
        g.nfaces = 2 # Number of faces per cell
        g.npoints_per_face = 1  # Number of points per face

        g.reference = ReferenceCell(basis, quadrature)

        N = length(quadrature.locations)-1

        g.LIFT = Lift1D(g.npoints, g.nfaces, g.npoints_per_face, g.reference.V)

        # Coordinates of all the vertices
        va = mesh.element_vertices[:,1]'
        vb = mesh.element_vertices[:,2]'
        g.x = ones(N+1,1)*mesh.vertex_coordinates[va] .+ 0.5*(quadrature.locations  .+ 1)*(mesh.vertex_coordinates[vb] - g.mesh.vertex_coordinates[va])

        # Geometric factors
        g.J = g.reference.Dr * g.x
        g.rx = 1.0 ./ g.J

        fmask1 = 1
        fmask2 = N+1
        Fmask = [fmask1;fmask2]
        Fx = g.x[Fmask,:]  # Locations of cell edges for each cell

        g.Fscale = 1.0 ./ g.J[Fmask,:]

        # Compute outward pointing normals
        g.nx = zeros(g.nfaces, mesh.nelements)
        g.nx[1,:] .= -1
        g.nx[2,:] .= 1

        # Build global connectivity arrays
        EToE, EToF = Connect1D(mesh.element_vertices)

        g.vmapP, g.vmapM, g.mapB, g.vmapB, g.mapI, g.mapO, g.vmapI, g.vmapO = BuildMaps1D(mesh.nelements, 2, 1, g.npoints, Fmask, EToE, EToF, g.x)

        return g
    end
end


text:
using BenchmarkTools
import PolyLog

function benchmark_and_print(f, z, str)
    println("----------------------------")
    println("Benchmarking ", str)
    println("----------------------------")
    b = @benchmark $f($z)
    io = IOBuffer()
    show(io, "text/plain", b)
    println(String(take!(io)))
end

benchmark_and_print(PolyLog.li2, 1.1        , "li2::Float64")
benchmark_and_print(PolyLog.li2, 1.1 + 0.5im, "li2::ComplexF64")
benchmark_and_print(PolyLog.li3, 1.1 + 0.5im, "li3::ComplexF64")
benchmark_and_print(PolyLog.li4, 1.1 + 0.5im, "li4::ComplexF64")
benchmark_and_print(PolyLog.li5, 1.1 + 0.5im, "li5::ComplexF64")
benchmark_and_print(PolyLog.li6, 1.1 + 0.5im, "li6::ComplexF64")


text:
#=

al cubic residual helmholtz models are from:

https://doi.org/10.6028/jres.121.011

Helmholtz Energy Transformations of Common Cubic Equations of State for Use with Pure Fluids and Mixtures

=#
function pressure_impl(mt::SingleVT,model::CubicModel,v,t)
    a,b,p = cubic_abp(mt,model,v,t)
    return p
end

function pressure_impl(mt::MultiVT,model::CubicModel,v,t,x)
    a,b,p = cubic_abp(mt,model,v,t,x)
    return p
end

"""
    cubic_mixing_rule(op, x, p,A=nothing)

returns an efficient implementation of:
`sum(x[i]*x[j]*op(i,j) for i = 1:n , j = 1:n)`

"""
function cubic_mixing_rule(op, x,A=nothing)
    if isnothing(A)
        N = length(x)
        @boundscheck checkbounds(x, N)
        @inbounds begin
            res1 = zero(eltype(x))
            for i = 1:N
                res1 += op(i,i) * x[i]^2
                for j = 1:i - 1
                    res1 += 2 * x[i] * x[j] * op(i, j)
                end
            end
        end
        return res1
    else
        N = length(x)
        @boundscheck checkbounds(x, N)
        @inbounds begin
            res1 = zero(eltype(x))
            for i = 1:N
                res1 += op(i,i) * x[i]^2
                for j = 1:i - 1
                    res1 += 2 * x[i] * x[j] * op(i, j) * A[i,j]
                end
            end
        end
        return res1
    end
end

#=
cubic polynomial:

k0 = -B((c1*c2)*B*(B+1)+A)
k1 = c1*c2*B*B-(c1+c2)*B*(B+1) + A
k2 = (c1+c2-1)*B-1
k3 = 1

on vdw:

c1=c2 = 0, c1*c2 = 0, c1+c2 = 0

k0 = -B((0)*B*(B+1)+A)
k1 = 0*B*B-(0)*B*(B+1) + A
k2 = (0-1)*B-1
k3 = 1

k0 = -AB
k1 =  A
k2 = -(B+1)
k3 = 1

on rk:

c1 = 1, c2 = 0, c1*c2 = 0,c1+c2 = 1

k0 = -B(0*B*(B+1)+A)
k1 = 0*B*B-(1)*B*(B+1) + A
k2 = (1-1)*B-1
k3 = 1

k0 = -AB
k1 = -B*(B+1) + A
k2 = -1
k3 = 1

on pr:

c1 = 1+√2, c2 = 1-√2 c1*c2 = -1,c1+c2 = 1

k0 = -B((-1)*B*(B+1)+A)
k1 = -B*B-(1)*B*(B+1) + A
k2 = (1-1)*B-1
k3 = 1

k0 = B*(B*(B+1)-A)
k1 = -2*B*B - B + A
k2 = -1
k3 = 1

k0 = B*(B*(B+1)-A)
k1 = -B*(2*B+1) + A
k2 = -1
k3 = 1

=#


text:
using ConstrainedDynamics
using ConstrainedDynamics: ∂dyn∂vel, discretizestate!
using ForwardDiff
using Rotations
using LinearAlgebra

function dyntestT()
    Δt = 0.01

    x1 = rand(3)
    v1 = rand(3)
    v2 = rand(3)

    q1 = rand(UnitQuaternion)
    ω1 = rand(3)
    ω2 = rand(3)

    origin = Origin{Float64}()
    body1 = Box(1., 1., 1., 1.)
    body1.m = 1.0
    body1.J = diagm(ones(3))

    oc1 = EqualityConstraint(Floating(origin, body1))

    mech = Mechanism(origin, [body1], [oc1])
    discretizestate!(body1,x1,q1,v1,v2,ω1,ω2,Δt)

    res = ForwardDiff.jacobian(dynTvel, [body1.state.vc;body1.state.vsol[2]])
    V2 = res[1:3,4:6]

    n = norm(V2 - ∂dyn∂vel(body1, Δt)[1:3,1:3])

    # display(n)
    return n
end

function dyntestR()
    Δt = 0.01

    x1 = rand(3)
    v1 = rand(3)
    v2 = rand(3)

    q1 = rand(UnitQuaternion)
    ω1 = rand(3)
    ω2 = rand(3)

    origin = Origin{Float64}()
    body1 = Box(1., 1., 1., 1.)
    body1.m = 1.0
    body1.J = diagm(ones(3))

    oc1 = EqualityConstraint(Floating(origin, body1))

    mech = Mechanism(origin, [body1], [oc1])
    discretizestate!(body1,x1,q1,v1,v2,ω1,ω2,Δt)

    res = ForwardDiff.jacobian(dynRvel, [body1.state.ωc;body1.state.ωsol[2]])
    W2 = res[1:3,4:6]

    n = norm(W2 - ∂dyn∂vel(body1, Δt)[4:6,4:6])

    # display(n)
    return n
end

for i=1:10
    @test isapprox(dyntestT(), 0.0; atol = 1e-8)
    @test isapprox(dyntestR(), 0.0; atol = 1e-8)
end


text:
#default params for the simulations

const params = Dict(
    "gca"      => 1000,     # maximal conductance ca membran channel
    "gkca"     => 50,       # maximal conducatnce ca dependet k channels
    "gk"       => 2700,     # maxiaml conductance of k ?
    "gkatpbar" => 25000,    # maximal conductance of atp dependent k channels

    "vca"      => 25,       # nernst potential ca2+
    "vm"       => -20,      # shapre parameter ca2+ activation
    "sm"       => 12,       # shape parameter ca2+ activation

    "vk"       => -75,      # nernst potential k
    "k_ext"    => 30,       # extracellular k
    "k_i_ext"  => 30,       # initial extracellular k
    "k_int"    => 280,      # intracellular k

    "sn"       => 5,        # shape parameter k
    "taun"     => 20,       # time constant rectifying current
    "nin"      => -16,      # shape parameter rectifying current
    "kd"       => 0.5,      # MgADP^2 factor KATP activation function
    "kdd"      => 17,       # MgADP factor KATP activation function
    "ktd"      => 26,       # ADP3- factor KATP activation function
    "ktt"      => 1,        # ATP4- factor KATP activation function
    "atot"     => 3000,     # apt + adp
    "alpha"    => 5.18e-18, # current to ion flux conversion factor
    "vcyt"     => 1.15e-12, # volume cytosol
    "kpmca"    => 0.2,      # pmca pump rate ca2+ cell membrane
    "kserca"   => 0.4,      # serca pump rate ca2+ er
    "pleak"    => 2.0e-4,   # leak across er membrane
    "kuni"     => 0.4,      # rate of mitochonria uniporters
    "knaca"    => 0.001,    # activity of mitochndrial ca na exchanger NCLX
    "Cm"       => 5300,     # membrane capacitance
    "fca"      => 0.01,     # fraction of ca2+ not bound to buffers
    "sigmam"   => 100,      # ratio of mitochonidria to cytosol volume
    "sigmaer"  => 31,       # ratio of ER to cytosol volume
    "Jgk"      => 0.001,    # glucokinase flux

    "k1"       => 30,       # params for pfk activity
    "k2"       => 1,        # params for pfk activity
    "k3"       => 50000,    # params for pfk activity
    "k4"       => 1000,     # params for pfk activity
    "famp"     => 0.02,     # params for pfk activity
    "fatp"     => 20,       # params for pfk activity
    "ffbp"     => 0.2,      # params for pfk activity
    "fbt"      => 20,       # params for pfk activity
    "fmt"      => 20,       # params for pfk activity
    "kpfk"     => 0.06,     # params for pfk activity

    "vpfk"     => 0.01,     # flux rate of pfk
    "taua"     => 300000,   # d ADP/dt time constant
    "kCaPDH"   => 200,      # factor for pdh activation function
    "vpdh"     => 0.009,    # flux rate of pdh

    "dummy"     => []       # set type to Dict{String, Any} (in julia 1.6.3 setting the global variable type
                            # is not supported. Therefore i use this hack, so i can add a Vector as a dict
                            # entry later on.)
    );


text:
@testset "BRDF" begin

@testset "Utilities" begin
    w = normalize(Vector3(1, 0, 1))

    @test GRAYTR.costheta(w) ≈ sqrt(2)/2
    @test GRAYTR.sintheta(w) ≈ sqrt(2)/2
    @test GRAYTR.sintheta2(w) ≈ 0.5
    @test GRAYTR.cosphi(w) ≈ 1.0
    @test GRAYTR.sinphi(w) ≈ 0.0
end

@testset "Materials" begin

S = Sphere()
R = GRAYTR.Ray(Point3(-1, 1/sqrt(2), 0), Vector3(1, 0, 0))
dg, a, b = GRAYTR.shape_intersect(R, S)
@test dg != nothing

@testset "World to local transformation" begin
    @test typeof(dg.T) == Transformation

    invT = inv(dg.T)

    @test dg.T(dg.n) ≈ Normal3(0.0, 0.0, 1.0)
    @test invT(Normal3(0.0, 0.0, 1.0)) ≈ dg.n

    u = invT(Normal3(0.0, 1.0, 0.0))
    v = invT(Normal3(1.0, 0.0, 0.0))

    @test dot(u, dg.n) ≈ 0.0
    @test dot(v, dg.n) ≈ 0.0
    @test dot(u, v) ≈ 0.0
end

end # testset Materials

spec = GRAYTR.SingleLine(532.0, 1.0)
w0 = normalize(Vector3(0, 0, 1))
w1 = normalize(Vector3(1, 0, 1))
ws = normalize(Vector3(-1, 0, 1))

@testset "Lambert" begin
    R = GRAYTR.Lambert(spec)

    @test GRAYTR.evaluate(R, w0, w1) == GRAYTR.SingleLine(532.0, 1) / π
    @test GRAYTR.rho(R, w0) == GRAYTR.SingleLine(532.0, 1)

    @test_broken GRAYTR.compute_pressure(R, w0, )

end

@testset "Lommel-Seeliger" begin
    P(alpha) = 1.0
    R = GRAYTR.LommelSeeliger(spec, P)

    @test GRAYTR.evaluate(R, w0, w1) == GRAYTR.SingleLine(532.0, 1) * 0.25 / (1.0 + sqrt(2)/2)

end

@testset "Ashikhmin-Shirley" begin
#    R = GRAYTR.AshikhminShirleySingle(spec, 0.5, 100.0)
#    @test GRAYTR.evaluate(R, w0, w1) <: GRAYTR.Spectrum
end

@testset "SpecularDiffuse" begin
    L = GRAYTR.Lambert(spec)
    R = GRAYTR.SpecularDiffuse(spec, 1.0, 0.0)

    @test GRAYTR.evaluate(R, w0, w1) == GRAYTR.SingleLine(532.0, 1) / π
    @test GRAYTR.evaluate(R, w0, w1) == GRAYTR.evaluate(L, w0, w1)

    R = GRAYTR.SpecularDiffuse(spec, 0.5, 0.5)
    @test GRAYTR.evaluate(R, w1, ws) == GRAYTR.SingleLine(532.0, 1) * (0.5 + 0.5/π)

    @test GRAYTR.compute_pressure(R, w0, spec) ≈ Vector3(0.0, 0.0, -(1 + 5/6))
end

end # main BRDF testset

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using Gnuplot
x = y = -5:1:5
z = [x^2 + y^2 for x in x, y in y]
@gp x y z "w image notit" "set auto fix" "set size square"
@gp :- xlab = "x" ylab = "y" "set cblabel 'z'" palette(:dense)
@gp :- title = "Heatmap Discrete"
save(term ="pngcairo size 600,400", output ="heatmap2.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# Only NEW SOCP Constraints in comparison to the original SOCP solver

using SparseArrays

include("../other_utils/feasibleCheck.jl")
include("projections.jl")

# -----------------------
# Simple Second-Order Cone Constraints
# -----------------------
@doc raw"""
    AL_simpleCone(tMax)

Creates a _simple_ Second Order Cone Constraint (e.g. for an SOCP) of the form
```math
||t|| ≤ tMax
```

So, in more standard form, this is becomes:
```math
||t|| - tMax ≤ 0
```

To check constraint satisfaction, use:
[`satisified`](@ref) and [`whichSatisfied`](@ref)

To evaluate the constraint, use:
- [`getRaw`](@ref) to evaluate `||t|| - tMax`.
- [`getNormToProjVals`](@ref) to evaluate the cone using projections.

"""
struct AL_simpleCone <: constraint
    tMax
end

"""
    getRaw(r::AL_simpleCone, t)

Evaluate the constraint without projection.
"""
getRaw(r::AL_simpleCone, t) = norm(t) - r.tMax

"""
    satisfied(r::AL_simpleCone, t)

Check constraint satisfaction. See also [`whichSatisfied`](@ref)
"""
satisfied(r::AL_simpleCone, t) = (getRaw(r, t) ≤ 0)

"""
    whichSatisfied(r::AL_simpleCone, tArr)

Check constraint satisfaction for each element in `tArr`.
"""
whichSatisfied(r::AL_simpleCone, tArr) = [satisfied(r, t) for t in tArr]

"""
    getProjVecs(r::AL_simpleCone, t, filled = true)

Gets the projection vectors for `t` onto the the second order cone defined
by `||t|| - tMax`.

See also [`AL_simpleCone`](@ref)
"""
function getProjVecs(r::AL_simpleCone, t, filled = true)
    return projSecondOrderCone(t, r.tMax, filled)
end

"""
    getNormToProjVals(r::AL_simpleCone, t, λ = 0)

Get the constraint violation for the second order cone constraints.

The parameter `λ` describes the dual of `||t|| - tMax` and is used to determine
if the constraint is active. See [`coneActive`](@ref) for more information.
"""
function getNormToProjVals(r::AL_simpleCone, t, λ = 0)

    active, filled = coneActive(t, r.tMax, λ)
    pv = getProjVecs(r, t, filled)

    return [norm([t; r.tMax] - pv)]
end

"""
    getGradC(r::AL_simpleCone, t)

Calculates the gradient of a constraint.

For `r::AL_simpleCone`, this is the `∇c(x) = x / ||x||`

There is a singularity at the tip.
"""
function getGradC(r::AL_simpleCone, t)
    if satisfied(r, t)
        return zeros(size(t, 1))
    end
    return t / norm(t)
end

"""
    getHessC(r::AL_simpleCone, t)

Calculate the hessian of a constraint.

For `r::AL_simpleCone`, `H(c(x)) = 0`
"""
function getHessC(r::AL_simpleCone, t)
    return spzeros(size(t, 1), size(t, 1))
end

"""
    getHessC_ALTerm(r::AL_simpleCone, t, rho = 1)

This is the hessian of the *constraint term* as it appears in the
augmented lagrangian.

For `r::AL_simpleCone`, `H(ρ c(x)'c(x) + λ c(x)) = ρ s * s' / norm(s)^2`
"""
function getHessC_ALTerm(r::AL_simpleCone, t, rho = 1)
    #=
    The augmented lagrangian constraint term is of the form:
    ρ c(x)'c(x) + λ c(x)

    where c(x) = x / ||x||.

    The Hessian matrix is then
    Htot = ρ (c.H + J.J + H.c) + λ H

    But H = 0 for Affine equalities, so when the constraint is not satisfied,
    Htot = ρ (J.J) = ρ s * s' / norm(s)^2
    When the constraint is satisfied, it is just zero.
    =#

    if satisfied(r, t)
        # Returns a sparse zero matrix of the right dimensions
        return getHessC(r, t)
    end

    return rho * (t * t') / (norm(t)^2)
end

# -----------------------
# Simple Second-Order Cone Constraints for Multiple
# -----------------------
@doc raw"""
    AL_Multiple_simpleCone(
                        sc::AL_simpleCone,
                        nDim::Int64,
                        indicatorList::Tuple{Array{Int64,1},Array{Int64,1}})

Handles Multiple Second Order Cone Constraint (SOCP) of the form
```math
||t_k|| ≤ t_{Max} \quad \text{for} \ t_k ∈ T
```
where T is a list of primal decision vectors.

So, in more standard form, this is becomes:
```math
||t_k|| - t_{Max} ≤ 0 \quad \text{for} \ t_k ∈ T
```

To check constraint satisfaction, use:
[`satisified`](@ref) and [`whichSatisfied`](@ref)

To evaluate the constraint, use:
- [`getRaw`](@ref) to evaluate `||t|| - tMax`.
- [`getNormToProjVals`](@ref) to evaluate the projection to the affine
equalities and (more critically) the second order cone.

"""
struct AL_Multiple_simpleCone <: constraint
    sc::AL_simpleCone
    nDim::Int64
    #Type of findnz(sparseArray)[1] or findall(!iszero, sparseArray)
    indicatorList::Array{Int64,1}
end

function makeAL_Multiple_simpleCone(tMax::Float64, nDim::Int64,
                                    indicatorList::Array{Int64,1})
    return AL_Multiple_simpleCone(AL_simpleCone(tMax), nDim, indicatorList)
end

"""
    parseRelevantSteps(r::AL_Multiple_simpleCone, x)

Run through the indicator list and use the dimensions to extract the
relevant vectors from x
"""
function parseRelevantSteps(r::AL_Multiple_simpleCone, x)
    # This array will store the relevant parts of the x vector.
    steps = []

    for ind in r.indicatorList
        # Find the relevant parts and push them as a group
        push!(steps, x[ind:(ind + r.nDim - 1)])
    end

    return steps
end

"""
    getRaw(r::AL_Multiple_simpleCone, x)

Evaluate each of the simple cone constraints without projection.
"""
function getRaw(r::AL_Multiple_simpleCone, x)
    # Use the parsed steps and calculate the raw value for each set.
    steps = parseRelevantSteps(r, x)
    return [getRaw(r.sc, s) for s in steps]
end

"""
    getProjVecs(r::AL_Multiple_simpleCone, steps, filled)

Get the projection defined by `||t|| - tMax` for each of the constraints.

See also [`AL_Multiple_simpleCone`](@ref)
"""
function getProjVecs(r::AL_Multiple_simpleCone, steps, filled)
    # Then loop through the parsed steps to get the projection.
    pv = []
    for (i, s) in enumerate(steps)
        # Filled is needed to determine which projection and may be different
        # for each of the different SOCP constraints in question.
        push!(pv, getProjVecs(r.sc, s, filled[i]))
    end

    return pv
end

function getProjVecs(r::AL_Multiple_simpleCone, x::Array{Float64, 1}, filled)
    steps = parseRelevantSteps(r, x)
    return getProjVecs(r, steps, filled)
end

function getProjVecs(r::AL_Multiple_simpleCone, x::Array{Int64, 1}, filled)
    steps = parseRelevantSteps(r, x)
    return getProjVecs(r, steps, filled)
end

"""
    getNormToProjVals(r::AL_Multiple_simpleCone, x, λ)

Get the constraint violation for multiple second order cone constraints. The
length of the vector `λ` must match the number of steps.

See also [`AL_Multiple_simpleCone`](@ref) and [`coneActive`](@ref) for more
information.
"""
function getNormToProjVals(r::AL_Multiple_simpleCone, x, λ)
    # Use the parsed steps to separate our each constraint
    steps = parseRelevantSteps(r, x)

    # Determine which of the constraints (if any) require an unfilled cone/
    fillVec = Bool[]
    for (ind, s) in enumerate(steps)
        push!(fillVec, coneActive(s, r.sc.tMax, λ[ind])[2])
    end

    # Using the filled information and the steps in question, get the
    # projection vectors
    pvs = getProjVecs(r, steps, fillVec)

    # Use the projection vectors to calculate the constraints.
    constVio = []
    for (s, p) in zip(steps, pvs)
        push!(constVio, norm([s; r.sc.tMax] - p))
    end

    # If there are NaN values (this happens if you take the norm of a zero
    # vector that is large), use the snippet of code below.
    return replace(constVio, NaN=>0.0)
end

function getNormToProjVals(r::AL_Multiple_simpleCone, x)
    lambdaVec = spzeros(size(r.indicatorList, 1))
    return getNormToProjVals(r, x, lambdaVec)
end

"""
    getGradC(r::AL_Multiple_simpleCone, x)

Parses the relevant `x` vector for each cone constraint. It calculates the
gradient of each and then stiches it into a jacobian of the same dimension
as `x` and number of constraints with the zeros in relevant indicies.

(Result is a sparse array)

See also [`AL_Multiple_simpleCone`](@ref) and [`AL_simpleCone`](@ref)
"""
function getGradC(r::AL_Multiple_simpleCone, x)
    # Use the parsed steps to separate our each constraint
    steps = parseRelevantSteps(r, x)

    # Use the AL_simpleCone to get the gradient for each constraint
    gradList = [getGradC(r.sc, s) for s in steps]

    # Stitch it back together.
    gradStitched = spzeros(size(r.indicatorList, 1), size(x, 1))

    for (ind, pos) in enumerate(r.indicatorList)
        posEnd = pos + r.nDim - 1
        gradStitched[ind, pos:posEnd] = gradList[ind]
    end

    return gradStitched
end

"""
    getHessC(r::AL_Multiple_simpleCone, x)

Calculate the hessian of a constraint.

For `r::AL_Multiple_simpleCone`, `H(c(x)) = 0`

See also [`AL_Multiple_simpleCone`](@ref) and [`AL_simpleCone`](@ref)
"""
function getHessC(r::AL_Multiple_simpleCone, x)
    return spzeros(size(x, 1), size(x, 1))
end

"""
    getHessC_ALTerm(r::AL_Multiple_simpleCone, x, rho = 1)

This is the hessian of the *constraint term* as it appears in the
augmented lagrangian.

For `r::AL_Multiple_simpleCone`,
`H(ρ c(x)'c(x) + λ c(x)) = ρ s * s' / norm(s)^2` for each constraint.

Parses the relevant `x` vector for each cone constraint. It calculates the
hessian of each and then stiches it into a vector of the same dimension
as `x`×`x` with the zeros in relevant indicies.

(Result is a sparse array)

See also [`AL_Multiple_simpleCone`](@ref) and [`AL_simpleCone`](@ref)
"""
function getHessC_ALTerm(r::AL_Multiple_simpleCone, x, rho = 1)
    #=
    The augmented lagrangian constraint term is of the form:
    ρ c(x)'c(x) + λ c(x)

    where c(x) = x / ||x||.

    The Hessian matrix is then
    Htot = ρ (c.H + J.J + H.c) + λ H

    But H = 0 for Affine equalities, so
    Htot = ρ (J.J) = ρ s * s' / norm(s)^2

    For each constraint!
    =#

    # Use the parsed steps to separate our each constraint
    steps = parseRelevantSteps(r, x)

    # Use the AL_simpleCone to get the total hessian for each constraint
    hessList = [Symmetric(getHessC_ALTerm(r.sc, s)) for s in steps]

    # Stitch it back together.
    hessStitched = spzeros(size(x, 1), size(x, 1))

    for (ind, pos) in enumerate(r.indicatorList)
        posEnd = pos + r.nDim - 1
        hessStitched[pos:posEnd, pos:posEnd] = hessList[ind]
    end

    return hessStitched

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Nettle_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Nettle")
JLLWrappers.@generate_main_file("Nettle", UUID("4c82536e-c426-54e4-b420-14f461c4ed8b"))
end  # module Nettle_jll


text:
using DashBootstrapComponents, DashHtmlComponents

fade = html_div([
    dbc_button("Toggle fade", id="fade-button", className="mb-3", n_clicks=0),
    dbc_fade(
        dbc_card(
            dbc_cardbody(html_p("This content fades in and out", className="card-text")),
        ),
        id="fade",
        is_in=true,
        appear=false,
    ),
]);

callback!(
    app,
    Output("fade", "is_in"),
    Input("fade-button", "n_clicks"),
    State("fade", "is_in"),
) do n, is_in
    return n > 0 ? is_in == 0 : true
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module ExportAllTest
  using Test
  include("Bar.jl")
  using .Bar
  using ExportAll
@testset "Export test" begin
  @test foo() == 1
  @test bar() == 2
  @test barBar() == 3
end

end #= ExportAllTest =#


text:
module AdventOfCode2018

include("get_input.jl")

for day in 1:11
	include("Day$day.jl")
end

include("Day16.jl")

end


text:
using Rotations

"""
    box
        particle with contacts at each corner
        orientation representation: modified rodrigues parameters
"""
struct Box{T}
    n::Int
    m::Int
    d::Int

    mass::T
    J::T
    μ::T
    g::T

    r
	n_corners
    corner_offset

    nq
    nu
    nc
    nf
    nb

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s

end

# Dimensions
nq = 6 # configuration dimension
nu = 3 # control dimension
nc = 8 # number of contact points
nf = 4 # number of faces for friction cone pyramid
nb = nc*nf
ns = 1

n = 2*nq
m = nu + nc + nb + nc + nb + 1
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
r = 0.5
c1 = @SVector [r, r, r]
c2 = @SVector [r, r, -r]
c3 = @SVector [r, -r, r]
c4 = @SVector [r, -r, -r]
c5 = @SVector [-r, r, r]
c6 = @SVector [-r, r, -r]
c7 = @SVector [-r, -r, r]
c8 = @SVector [-r, -r, -r]

corner_offset = @SVector [c1, c2, c3, c4, c5, c6, c7, c8]

# Parameters
μ = 1.0  # coefficient of friction
g = 9.81
mass = 1.0   # mass
J = 1.0 / 12.0 * mass * ((2.0 * r)^2 + (2.0 * r)^2)

# Methods
M_func(model::Box, q) = Diagonal(@SVector [model.mass, model.mass, model.mass, model.J, model.J, model.J])

G_func(model::Box, q) = @SVector [0., 0., model.mass * model.g, 0., 0., 0.]

function ϕ_func(model::Box, q)
    p = view(q, 1:3)
    rot = view(q, 4:6)

    R = MRP(rot...)

    @SVector [(p + R * model.corner_offset[1])[3],
              (p + R * model.corner_offset[2])[3],
              (p + R * model.corner_offset[3])[3],
              (p + R * model.corner_offset[4])[3],
              (p + R * model.corner_offset[5])[3],
              (p + R * model.corner_offset[6])[3],
              (p + R * model.corner_offset[7])[3],
              (p + R * model.corner_offset[8])[3]]
end

function B_func(::Box, q)
    rot = view(q, 4:6)
    R = MRP(rot...)
    @SMatrix [0. 0. 0. R[1,1] R[2,1] R[3,1];
              0. 0. 0. R[1,2] R[2,2] R[3,2];
              0. 0. 0. R[1,3] R[2,3] R[3,3]]
end

function N_func(model::Box, q)
    tmp(z) = ϕ_func(model, z)
    ForwardDiff.jacobian(tmp, q)
end

function P_func(model::Box, q)
    map = [1. 0.;
           0. 1.;
           -1. 0.;
           0. -1.]

    function p(x)
        pos = view(x, 1:3)
        rot = view(x, 4:6)

        R = MRP(rot...)

        [map * (pos + R * model.corner_offset[1])[1:2];
         map * (pos + R * model.corner_offset[2])[1:2];
         map * (pos + R * model.corner_offset[3])[1:2];
         map * (pos + R * model.corner_offset[4])[1:2];
         map * (pos + R * model.corner_offset[5])[1:2];
         map * (pos + R * model.corner_offset[6])[1:2];
         map * (pos + R * model.corner_offset[7])[1:2];
         map * (pos + R * model.corner_offset[8])[1:2]]
    end

    ForwardDiff.jacobian(p, q)
end

function friction_cone(model::Box,u)
    λ = u[model.idx_λ]
    b = u[model.idx_b]

    @SVector [model.μ * λ[1] - sum(b[1:4]),
              model.μ * λ[2] - sum(b[5:8]),
              model.μ * λ[3] - sum(b[9:12]),
              model.μ * λ[4] - sum(b[13:16]),
              model.μ * λ[5] - sum(b[17:20]),
              model.μ * λ[6] - sum(b[21:24]),
              model.μ * λ[7] - sum(b[25:28]),
              model.μ * λ[8] - sum(b[29:32])]
end

function maximum_dissipation(model::Box, x⁺, u, h)
    q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

    ψ = u[model.idx_ψ]
    ψ_stack = [ψ[1] * ones(4);
               ψ[2] * ones(4);
               ψ[3] * ones(4);
               ψ[4] * ones(4);
               ψ[5] * ones(4);
               ψ[6] * ones(4);
               ψ[7] * ones(4);
               ψ[8] * ones(4)]

    η = u[model.idx_η]

    P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function fd(model::Box, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    [q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{6}(q2⁺) - SVector{6}(q1))
    - M_func(model, q2⁺) * (SVector{6}(q3) - SVector{6}(q2⁺)))
    + transpose(B_func(model, q3)) * SVector{3}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{8}(λ)
    + transpose(P_func(model, q3)) * SVector{32}(b)
    - h * G_func(model, q2⁺))]
end

model = Box(n, m, d,
			mass, J, μ, g,
            r, 8, corner_offset,
            nq, nu, nc, nf, nb,
            idx_u,
            idx_λ,
            idx_b,
            idx_ψ,
            idx_η,
            idx_s)

function visualize!(vis, model::Box, q;
        Δt = 0.1)

    setobject!(vis["box"], Rect(Vec(-1.0 * model.r,
		-1.0 * model.r,
		-1.0 * model.r),
		Vec(2.0 * model.r, 2.0 * model.r, 2.0 * model.r)))

    for i = 1:model.n_corners
        setobject!(vis["corner$i"], Sphere(Point3f0(0),
            convert(Float32, 0.05)),
            MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))
    end

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        MeshCat.atframe(anim, t) do

            settransform!(vis["box"],
				compose(Translation(q[t][1:3]...), LinearMap(MRP(q[t][4:6]...))))

            for i = 1:model.n_corners
                settransform!(vis["corner$i"],
                    Translation((q[t][1:3] + MRP(q[t][4:6]...) * (corner_offset[i]))...))
            end
        end
    end
    # settransform!(vis["/Cameras/default"], compose(Translation(-1, -1, 0),LinearMap(RotZ(pi/2))))
    MeshCat.setanimation!(vis, anim)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""

# shipley_test

$(SIGNATURES)

Test of all independencies implied by a given DAG

Computes a simultaneous test of all independence relationships implied by a
given Gaussian model defined according to a directed acyclic graph, based on
the sample covariance matrix.

The test statistic is C = -2 sum(ln(p_j)) where p_j are the
p-values of tests of conditional independence in the basis set computed by
basiSet(A). The p-values are independent uniform variables on
(0,1) and the statistic has exactly a chi square distribution on
2k degrees of freedom where k is the number of elements of the
basis set.  Shipley (2002) calls this test Fisher's C test.

### Method
```julia
shipley_test(;
* `d::Dag`                             : Directed acyclic graph
)
```

### Returns
```julia
* `res::NamedTuple`                    : (ctest=..., dof=..., pval=...)
```
where:

  ctest: Test statistic C
  dof:   Degrees of freedom.
  pval:  The P-value of the test, assuming a two-sided alternative.

# Extended help

### Example

### Shipley_test for the mathematics marks data
```julia
using StructuralCausalModels, RData

objs = RData.load(scm_path("..", "data", "marks.rda");
marks_df = objs["marks"]

d = OrderedDict(
  :mechanics => [:vectors, :algebra],
  :vectors => [:algebra],
  :statistics => [:algebra, :analysis],
  :analysis => [:algebra]
);
dag = Dag(d; df=df)
shipley_test(dag)
```

### See also
```julia
?Dag
?basis_set
?pcor_test
```

### Acknowledgements

Original author:                       Giovanni M. Marchetti

Translated to Julia:                   Rob J Goedman

### References

Shipley, B. (2000). A new inferential test for path models based
on directed acyclic graphs. Structural Equation Modeling, 7(2),
206--218.

### Licence

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

The Julia translation is licenced under: MIT.

Part of the api, exported.
"""
function shipley_test(d::DAG)

  @assert isposdef(d.s) "Covariance matrix is not positive definite."
  n = nrow(d.df)

  # see pcor()
  function pval(r, q, n)
    df = n - 2 - q
    tval = r * sqrt(df) / sqrt(1 - r * r)
    2 * cdf(TDist(df), -abs(tval))
  end

  l = basis_set(d)
  k = length(l)
  if k == 0
    @info "Basis set empty, Shipley test not possible."
    return nothing
  end
  p = zeros(k)
  for i in 1:k
    u = l[i]
    println(u)
    r = pcor(d, u)
    q = length(u) - 2
    p[i] = pval(r, q, n)
  end
  ctest = -2 * sum(log.(p))
  df = 2 * k
  pv = 1 - cdf(Chisq(df), ctest)
  (ctest=ctest, df=df, pv=pv,)
end

export
  shipley_test