text:
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots

# external velocities
Vr(t, x) = 2.
Vl(t, x) = -2.
# interactions
W_attr(t, r) = 5 * log(abs(r) + 1)
W_rep(t, r) = -2 * log(abs(r) + 1)
W(t, r) = 2 * (exp(abs(r)/4) + exp(-2abs(r)))
# mobilities
mobρ(ρ, σ) = max(2 - ρ - 0.5σ, 0)
mobσ(ρ, σ) = max(2 - σ - 0.5ρ, 0)

imodel = IntegratedModel(
    (Vr, Vl),
    ((W, W_rep),
     (W_rep, W)),
    (mobρ, mobσ))

n = 20
x0 = ArrayPartition(
    vcat(range(-2., -1.5, length=n), range(-1., -.5, length=n)),
    vcat(range(.5, 1.5, length=2n)))

tspan = (0., 1.2)

prob = ODEProblem(velocities!, x0, tspan, imodel)

sol = solve(prob, BS5(), reltol=1e-6, abstol=1e-6)

plot(legend=false)
plot!(sol, vars=1:2n, color=:blue)
plot!(sol, vars=2n+1:4n, color=:red)
plot!(title="2-species integrated scheme", xlabel="time", ylabel="position")

savefig("two-species.png")

n = 60
x0 = ArrayPartition(
    vcat(range(-2., -1.5, length=n), range(-1., -.5, length=n)),
    vcat(range(.5, 1.5, length=2n)))

prob = ODEProblem(velocities!, x0, tspan, imodel)

sol = solve(prob, BS5(), reltol=1e-7, abstol=1e-7)

anim = @animate for t in tspan[1]:1/48:tspan[2]
    p = plot(title="2-species integrated scheme", xlabel="position", ylabel="density",
        legend=false, xlims=(-4,4), ylims=(0,1))
    plot_density!(p, sol(t).x[1], color=:blue)
    plot_density!(p, sol(t).x[2], color=:red)
end

gif(anim, "two-species.gif")


text:
push!(LOAD_PATH,"../src")

using Test
using TypesParticles
using TypesBasis

"""
		run_typesparticles_test()

This is a function to test the user types in TypesParticles.jl. In asserts basics
construction of the user data type for hydrogen.

""" function run_typesparticles_test()
	 @testset "Types" begin
	 #H2 structure
	 x0,y0,z0 = 0.00e0,0.00e0,0.00e0;
	 x1,y1,z1 = 0.77e0,0.00e0,0.00e0;
	 Z0,Z1 = 1,1;
	 H1 = Atom(x0,y0,z0,Z0);
	 H2 = Atom(x1,y1,z1,Z1);

	 cell = [1.00e0 0.00e0 0.00e0;
	      	 0.00e0 1.00e0 0.00e0;
		 0.00e0 0.00e0 1.00e0];

		@test begin
	       	H1.x == x0;
	       	H1.y == y0;
	       	H1.z == z0;
	       	H1.Z == Z0;
	       	end

		@test begin
	       	atomlist = [H1,H2];
	       	system = SystemOfAtoms(2,atomlist,cell);
	       	system.natoms == 2;
	       	typeof(system.atoms) == typeof(atomlist);
		system.atoms[1] == H1;
		system.atoms[2] == H2;
	       	system.cell == cell;
	       	end

	 end #testset
end #run_typeparticles_test

"""
		run_typesbasis_test()

This is a function to test the user types in TypesBasis.jl. In asserts basic
construction of the user data type for H 6-31G (Gaussian) basis set.

""" function run_typesbasis_test()
	 @testset "Basis" begin
	 #Hydrogen 6-31G
	 basistype="6-31G";
	 nprimatives=4;
	 coefficients = [ 0.03349460;
	 	          0.23472695;
			  0.81375733;
			  1.00000000];
	 alphas = [18.7311370;
	 	   2.8253937;
		   0.6401217;
		   0.1612778];
         norms = (2.00e0.*alphas./pi).^(3/4);
	 @test begin
	      basisset = Gaussian(basistype,nprimatives,
		       coefficients,alphas);

	      basisset.info == basistype;
	      basisset.nprims == nprimatives;
	      basisset.coefs == coefficients;
	      basisset.alphas == alphas;
	 end

	 @test begin
	       x0,y0,z0 = 0.00e0, 0.00e0, 0.00e0;
	       Z0 = 1;
	       H1 = Atom(x0,y0,z0,Z0);

	       basisset = Gaussian(basistype,nprimatives,
		       coefficients,alphas);
	       gaussorb1 = GaussOrbitals(x0,y0,z0,basisset);
	       gaussorb2 = GaussOrbitals(H1,basisset);

	       gaussorb1.ax == x0;
	       gaussorb2.ax == x0;
	       gaussorb1.ay == y0;
	       gaussorb2.ay == y0;
	       gaussorb1.az == z0;
	       gaussorb2.az == z0;

	       gaussorb1.coefs == coefficients;
	       gaussorb2.coefs == coefficients;

	       gaussorb1.alphas == alphas;
	       gaussorb2.alphas == alphas;

	       gaussorb1.norms == norms;
	       gaussorb2.norms == norms;
	   end
	   end #testset
end #run_typesbasis_test()

ext: jl


text:
export linear_regression

import ..LinearRegression

using mlpack._Internal.io

import mlpack_jll
const linear_regressionLibrary = mlpack_jll.libmlpack_julia_linear_regression

# Call the C binding of the mlpack linear_regression binding.
function linear_regression_mlpackMain()
  success = ccall((:linear_regression, linear_regressionLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module linear_regression_internal
  import ..linear_regressionLibrary

import ...LinearRegression

# Get the value of a model pointer parameter of type LinearRegression.
function IOGetParamLinearRegression(paramName::String)::LinearRegression
  LinearRegression(ccall((:IO_GetParamLinearRegressionPtr, linear_regressionLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type LinearRegression.
function IOSetParamLinearRegression(paramName::String, model::LinearRegression)
  ccall((:IO_SetParamLinearRegressionPtr, linear_regressionLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeLinearRegression(stream::IO, model::LinearRegression)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeLinearRegressionPtr, linear_regressionLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeLinearRegression(stream::IO)::LinearRegression
  buffer = read(stream)
  LinearRegression(ccall((:DeserializeLinearRegressionPtr, linear_regressionLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    linear_regression(; [input_model, lambda, test, training, training_responses, verbose])

An implementation of simple linear regression and simple ridge regression using
ordinary least squares. This solves the problem

  y = X * b + e

where X (specified by `training`) and y (specified either as the last column of
the input matrix `training` or via the `training_responses` parameter) are known
and b is the desired variable.  If the covariance matrix (X'X) is not
invertible, or if the solution is overdetermined, then specify a Tikhonov
regularization constant (with `lambda`) greater than 0, which will regularize
the covariance matrix to make it invertible.  The calculated b may be saved with
the `output_predictions` output parameter.

Optionally, the calculated value of b is used to predict the responses for
another matrix X' (specified by the `test` parameter):

   y' = X' * b

and the predicted responses y' may be saved with the `output_predictions` output
parameter.  This type of regression is related to least-angle regression, which
mlpack implements as the 'lars' program.

For example, to run a linear regression on the dataset `X` with responses `y`,
saving the trained model to `lr_model`, the following command could be used:

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> y = CSV.read("y.csv")
julia> lr_model, _ = linear_regression(training=X,
            training_responses=y)
```

Then, to use `lr_model` to predict responses for a test set `X_test`, saving the
predictions to `X_test_responses`, the following command could be used:

```julia
julia> using CSV
julia> X_test = CSV.read("X_test.csv")
julia> _, X_test_responses = linear_regression(input_model=lr_model,
            test=X_test)
```

# Arguments

 - `input_model::LinearRegression`: Existing LinearRegression model to
      use.
 - `lambda::Float64`: Tikhonov regularization for ridge regression.  If 0,
      the method reduces to linear regression.  Default value `0`.

 - `test::Array{Float64, 2}`: Matrix containing X' (test regressors).
 - `training::Array{Float64, 2}`: Matrix containing training set X
      (regressors).
 - `training_responses::Array{Float64, 1}`: Optional vector containing y
      (responses). If not given, the responses are assumed to be the last row of
      the input file.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::LinearRegression`: Output LinearRegression model.
 - `output_predictions::Array{Float64, 1}`: If --test_file is specified,
      this matrix is where the predicted responses will be saved.

"""
function linear_regression(;
                           input_model::Union{LinearRegression, Missing} = missing,
                           lambda::Union{Float64, Missing} = missing,
                           test = missing,
                           training = missing,
                           training_responses = missing,
                           verbose::Union{Bool, Missing} = missing,
                           points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, linear_regressionLibrary), Nothing, ());

  IORestoreSettings("Simple Linear Regression and Prediction")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    linear_regression_internal.IOSetParamLinearRegression("input_model", convert(LinearRegression, input_model))
  end
  if !ismissing(lambda)
    IOSetParam("lambda", convert(Float64, lambda))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if !ismissing(training_responses)
    IOSetParamRow("training_responses", training_responses)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("output_predictions")
  # Call the program.
  linear_regression_mlpackMain()

  return linear_regression_internal.IOGetParamLinearRegression("output_model"),
         IOGetParamRow("output_predictions")
end


text:
using Base: @deprecate, @deprecate_binding, depwarn

# Deprecate in Julia 0.6 cycle
function Base.isnan(da::DataArray)
    depwarn("vectorized method isnan(da) is deprecated, use isnan.(da) instead", :isnan)
    return isnan.(da)
end

@deprecate isna(x::AbstractArray) ismissing.(x)
@deprecate anyna(x) any(ismissing, x)
@deprecate allna(x) all(ismissing, x)

function reldiff(v::Vector{T}) where T
    depwarn("reldiff is deprecated.", :reldiff)
    n = length(v)
    res = Array(T, n - 1)
    for i in 2:n
        res[i - 1] = v[i] / v[i - 1]
    end
    return res
end

function percent_change(v::Vector{T}) where T
    depwarn("percent_change is deprecated.", :percent_change)
    n = length(v)
    res = Array(T, n - 1)
    for i in 2:n
        res[i - 1] = (v[i] - v[i - 1]) / v[i - 1]
    end
    return res
end

mutable struct xtab{T}
    vals::Array{T}
    counts::Vector{Int}
end

function xtab(x::AbstractArray{T}) where T
    depwarn("xtab is deprecated; use the FreqTables package instead.", :xtab)
    d = Dict{T, Int}()
    for el in x
        d[el] = get(d, el, 0) + 1
    end
    kk = sort(keys(d))
    cc = Array(Int, length(kk))
    for i in 1:length(kk)
        cc[i] = d[kk[i]]
    end
    return xtab(kk, cc)
end

function xtabs(x::AbstractArray{T}) where T
    depwarn("xtabs is deprecated; use the FreqTables package instead.", :xtabs)
    d = Dict{T, Int}()
    for el in x
        d[el] = get(d, el, 0) + 1
    end
    return d
end

# Implicitly vectorized bitwise operations

import Base: &, |, xor

for f in [:(&), :(|), :(xor)]
    for T in [:(BitArray), :(AbstractRange{<:Integer}), :(AbstractArray{<:Integer}), :(Integer)]
        @eval begin
            @deprecate ($f)(a::DataArray{<:Integer}, b::$T) ($f).(a, b)
            @deprecate ($f)(a::$T, b::DataArray{<:Integer}) ($f).(a, b)
        end
    end
    @eval begin
        @deprecate ($f)(a::DataArray{<:Integer}, b::DataArray{<:Integer}) ($f).(a, b)
    end
end

@deprecate_binding NAtype Missing
@deprecate_binding NA missing
@deprecate_binding NAException MissingException
@deprecate isna ismissing
@deprecate dropna(x) collect(skipmissing(x))
@deprecate padna padmissing
@deprecate each_failna Missings.fail
@deprecate each_dropna skipmissing
@deprecate each_replacena Missings.replace
@deprecate_binding EachFailNA DataArrays.EachFailMissing
@deprecate_binding EachDropNA DataArrays.EachDropMissing
@deprecate_binding EachReplaceNA DataArrays.EachReplaceMissing
import SpecialFunctions: digamma, erf, erfc
@deprecate digamma(x::Missing) ismissing(x) ? missing : SpecialFunctions.digamma(x)
@deprecate erf(x::Missing) ismissing(x) ? missing : SpecialFunctions.erf(x)
@deprecate erfc(x::Missing) ismissing(x) ? missing : SpecialFunctions.erfc(x)


text:
open("sample.txt") do f
    for line in readlines(f)[1:5]
        println("line: $line")
    end
end


text:
using SeaPearl
using SeaPearlExtras
using ReinforcementLearning
const RL = ReinforcementLearning
using Flux
using GeometricFlux
using JSON
using BSON: @load, @save
using Random
using Dates

# -------------------
# Generator
# -------------------
n_city = 21
grid_size = 100
max_tw_gap = 0
max_tw = 100
tsptw_generator = SeaPearl.TsptwGenerator(n_city, grid_size, max_tw_gap, max_tw, true)
SR = SeaPearl.TsptwStateRepresentation{SeaPearl.TsptwFeaturization, SeaPearl.TsptwTrajectoryState}

# -------------------
# Internal variables
# -------------------
numInFeatures=SeaPearl.feature_length(SR)

# -------------------
# Experience variables
# -------------------
nbEpisodes = 500
evalFreq = 200
nbInstances = 10
nbRandomHeuristics = 1

# -------------------
# Agent definition
# -------------------
include("agents.jl")

# -------------------
# Value Heuristic definition
# -------------------
learnedHeuristic = SeaPearl.LearnedHeuristic{SR, SeaPearl.TsptwReward, SeaPearl.VariableOutput}(agent)
include("nearest_heuristic.jl")
nearest_heuristic = SeaPearl.BasicHeuristic(select_nearest_neighbor) # Basic value-selection heuristic

function select_random_value(x::SeaPearl.IntVar; cpmodel=nothing)
    selected_number = rand(1:length(x.domain))
    i = 1
    for value in x.domain
        if i == selected_number
            return value
        end
        i += 1
    end
    @assert false "This should not happen"
end

randomHeuristics = []
for i in 1:nbRandomHeuristics
    push!(randomHeuristics, SeaPearl.BasicHeuristic(select_random_value))
end

valueSelectionArray = [learnedHeuristic, nearest_heuristic]
append!(valueSelectionArray, randomHeuristics)

# -------------------
# Variable Heuristic definition
# -------------------
struct TsptwVariableSelection{TakeObjective} <: SeaPearl.AbstractVariableSelection{TakeObjective} end
TsptwVariableSelection(;take_objective=false) = TsptwVariableSelection{take_objective}()
function (::TsptwVariableSelection{false})(cpmodel::SeaPearl.CPModel; rng=nothing)
    for i in 1:length(keys(cpmodel.variables))
        if haskey(cpmodel.variables, "a_"*string(i)) && !SeaPearl.isbound(cpmodel.variables["a_"*string(i)])
            return cpmodel.variables["a_"*string(i)]
        end
    end
end
variableSelection = TsptwVariableSelection()

# -------------------
# -------------------
# Core function
# -------------------
# -------------------
function trytrain(nbEpisodes::Int)
    experienceTime = now()
    dir = mkdir(string("exp_",Base.replace("$(round(experienceTime, Dates.Second(3)))",":"=>"-")))
    expParameters = Dict(
        :nbEpisodes => nbEpisodes,
        :evalFreq => evalFreq,
        :nbInstances => nbInstances
    )
    open(dir*"/params.json", "w") do file
        JSON.print(file, expParameters)
    end

    metricsArray, eval_metricsArray = SeaPearl.train!(
    valueSelectionArray=valueSelectionArray,
    generator=tsptw_generator,
    nbEpisodes=nbEpisodes,
    strategy=SeaPearl.DFSearch(),
    variableHeuristic=variableSelection,
    out_solver=true,
    verbose = false,
    evaluator=SeaPearl.SameInstancesEvaluator(valueSelectionArray,tsptw_generator; evalFreq = evalFreq, nbInstances = nbInstances)
)

    trained_weights = params(agent.policy.learner.approximator.model)
    @save dir*"/model_weights_tsptw"*string(n_city)*".bson" trained_weights

    SeaPearlExtras.storedata(metricsArray[1]; filename=dir*"/tsptw_$(n_city)_training")
    SeaPearlExtras.storedata(eval_metricsArray[:,1]; filename=dir*"/tsptw_$(n_city)_trained")
    SeaPearlExtras.storedata(eval_metricsArray[:,2]; filename=dir*"/tsptw_$(n_city)_nearest")
    for i = 1:nbRandomHeuristics
        SeaPearlExtras.storedata(eval_metricsArray[:,i+2]; filename=dir*"/tsptw_$(n_city)_random$(i)")
    end

    return metricsArray, eval_metricsArray
end

metricsArray, eval_metricsArray = trytrain(nbEpisodes)
nothing

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using ASTInterpreter

immutable DummyState; end
Base.LineEdit.transition(s::DummyState, _) = nothing
dummy_state(interp) = ASTInterpreter.InterpreterState(interp, interp, 1, DummyState(), nothing, nothing)

# Steps through the whole expression using `s`
function step_through(interp)
    state = dummy_state(interp)
    while !(isa(interp.next_expr[2], Expr) && interp.next_expr[2].head == :return)
        ASTInterpreter.execute_command(state, interp, Val{:s}(), "s")
    end
    return interp.retval
end

@assert step_through(ASTInterpreter.enter_call_expr(nothing, :($(+)(1,2.5)))) == 3.5
@assert step_through(ASTInterpreter.enter_call_expr(nothing, :($(sin)(1)))) == sin(1)

# Step into generated functions
@generated function generatedfoo(T)
    :(return $T)
end
callgenerated() = generatedfoo(1)
interp = ASTInterpreter.enter_call_expr(nothing, :($(callgenerated)()))
state = dummy_state(interp)

# Step into the generated function itself
ASTInterpreter.execute_command(state, state.top_interp, Val{:sg}(), "sg")

# Should now be in generated function
interp = state.interp
ASTInterpreter.execute_command(state, state.top_interp, Val{:finish}(), "finish")
@assert interp.retval == :(return $(Int64))

# Now finish the regular function
interp = state.interp
ASTInterpreter.execute_command(state, state.top_interp, Val{:finish}(), "finish")
@assert interp.retval == Int64

# Optional arguments
function optional(n = sin(1))
    x = asin(n)
    cos(x)
end

interp = ASTInterpreter.enter_call_expr(nothing, :($(optional)()))
state = dummy_state(interp)
# First call steps in
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == nothing
# cos(1.0)
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == nothing
# return
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == cos(1.0)

# Macros
macro insert_some_calls()
    esc(quote
        x = sin(b)
        y = asin(x)
        z = sin(y)
    end)
end

# Work around the fact that we can't detect macro expansions if the macro
# is defined in the same file
include_string("""
function test_macro()
    a = sin(5)
    b = asin(a)
    @insert_some_calls
    z
end
""","file.jl")

interp = ASTInterpreter.enter_call_expr(nothing, :($(test_macro)()))
state = dummy_state(interp)
# a = sin(5)
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == nothing
# b = asin(5)
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == nothing
# @insert_some_calls
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == nothing
# return z
ASTInterpreter.execute_command(state, state.top_interp, Val{:n}(), "n")
@assert interp.retval == sin(5)


text:
function zeck(n)
    n <= 0 && return 0
    fib = [2,1]; while fib[1] < n unshift!(fib,sum(fib[1:2])) end
    dig = Int[]; for f in fib f <= n ? (push!(dig,1); n = n-f;) : push!(dig,0) end
    return dig[1] == 0 ? dig[2:end] : dig
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##### Steady State Reflectance from a Slab #####
# Taken from Eqn. 45 from Contini 1997

function DA_reflslabCW(ρ::Float64, β::Array{Float64,1},ndet::Float64, nmed::Float64, s::Float64)
	#s is slab thickness

	n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	z3m = Array{Float64}(undef, length(xs))
	z4m = Array{Float64}(undef, length(xs))

    Rt1 = Array{Float64}(undef, length(xs))
    Rt2 = Array{Float64}(undef, length(xs))

	if n == 1.0
		A= 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

    z3m = Float64[-2m.*s .- 4m.*ze .- zs for m in xs]
	z4m = Float64[-2m.*s .- (4m .- 2).*ze .+ zs for m in xs]

    R1 = @. (z3m*(ρ^2 + z3m^2)^(-3/2))*(1+(μa*(ρ^2 +z3m^2))^(1/2))*exp(-(μa*(ρ^2 +z3m^2)/D)^(1/2))
    R2 = @. (z4m*(ρ^2 + z4m^2)^(-3/2))*(1+(μa*(ρ^2 +z4m^2))^(1/2))*exp(-(μa*(ρ^2 +z4m^2)/D)^(1/2))

    Rt = -1/(4π)*(sum(R1) - sum(R2))

	return Rt
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
#This program is to test if only one node exists.
using Plasmo
using Gurobi
using PlasmoTest

graph = OptiGraph()
optimizer = Gurobi.Optimizer

#Add only one node to a ModelGraph
n1 = @optinode(graph)

#For this node
@variable(n1,0 <= x <= 2)
@variable(n1,2 <= y <= 4)
@constraint(n1,x + y >= 5)

#Graph objective
@objective(graph, Min, n1[:y])

#Use Ipopt to solve
optimize!(graph, optimizer)

println("y = ",value(n1,n1[:y]))

max_stars_repo_path: Julia studying/PlasmoTest/test/1Node.jl

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_path: Julia studying/PlasmoTest/test/1Node.jl

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_path: Julia studying/PlasmoTest/test/1Node.jl

max_forks_repo_name: BoyangSheng/Skill-studying


text:
module Simple

include("../src/Problox.jl")
using .Problox

coin = term("coin")
heads = term("heads")
tails = term("tails")
win = term("win")

# Little generator - generates worlds :)
@problox function generator(p, q)
    C = variable(:C);
    coin(:c1);
    coin(:c2);
    (p :: heads(C), q :: tails(C)) :- coin(C);
    win << heads(C);
    query(win);
end

net = generator(0.3, 0.5)
ret = evaluate(net)
println(ret)

end # module


text:
import Base: ∈

export AbstractZonotope,
       genmat,
       generators,
       ngens,
       order,
       togrep

"""
    AbstractZonotope{N<:Real} <: AbstractCentrallySymmetricPolytope{N}

Abstract type for zonotopic sets.

### Notes

Mathematically, a zonotope is defined as the set

```math
Z = \\left\\{ c + ∑_{i=1}^p ξ_i g_i,~~ ξ_i \\in [-1, 1]~~ ∀ i = 1,…, p \\right\\},
```
where ``c \\in \\mathbb{R}^n`` is its center and ``\\{g_i\\}_{i=1}^p``,
``g_i \\in \\mathbb{R}^n``, is the set of generators.
This characterization defines a zonotope as the finite Minkowski sum of line
segments.
Zonotopes can be equivalently described as the image of a unit infinity-norm
ball in ``\\mathbb{R}^n`` by an affine transformation.

See [`Zonotope`](@ref) for a standard implementation of this interface.

Every concrete `AbstractZonotope` must define the following functions:
- `genmat(::AbstractZonotope{N})` -- return the generator matrix
- `generators(::AbstractZonotope{N})` -- return an iterator over the generators

Since the functions `genmat` and `generators` can be defined in terms of each
other, it is sufficient to only genuinely implement one of them and let the
implementation of the other function call the fallback implementation
`genmat_fallback` resp. `generators_fallback`.

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractZonotope)
3-element Array{Any,1}:
 AbstractHyperrectangle
 LineSegment
 Zonotope
```
"""
abstract type AbstractZonotope{N<:Real} <: AbstractCentrallySymmetricPolytope{N} end

isconvextype(::Type{<:AbstractZonotope}) = true

# --- fallback implementations of AbstractZonotope functions ---

"""
    genmat_fallback(Z::AbstractZonotope{N}) where {N<:Real}

Fallback definition of `genmat` for zonotopic sets.

### Input

- `Z`     -- zonotopic set
- `gens`  -- (optional; default: `generators(Z)`) iterator over generators
- `ngens` -- (optional; default: `nothing`) number of generators or `nothing` if
             unknown

### Output

A matrix where each column represents one generator of `Z`.

### Notes

Passing the number of generators is much more efficient as otherwise the
generators have to be obtained from the iterator (`gens`) and stored in an
intermediate vector until the final result matrix can be allocated.
"""
function genmat_fallback(Z::AbstractZonotope{N};
                         gens=generators(Z),
                         ngens=nothing) where {N<:Real}
    if isempty(gens)
        return Matrix{N}(undef, dim(Z), 0)
    elseif ngens == nothing
        return _genmat_fallback_generic(Z, gens)
    else
        return _genmat_fallback_ngens(Z, gens, ngens)
    end
end

function _genmat_fallback_generic(Z::AbstractZonotope{N}, gens) where {N<:Real}
    Gv = Vector{Vector{N}}()
    @inbounds for (i, g) in enumerate(gens)
        push!(Gv, g)
    end
    G = Matrix{N}(undef, dim(Z), length(Gv))
    @inbounds for (i, g) in enumerate(Gv)
        G[:, i] = g
    end
    return G
end

function _genmat_fallback_ngens(Z::AbstractZonotope{N}, gens, ngens) where {N<:Real}
    G = Matrix{N}(undef, dim(Z), ngens)
    @inbounds for (i, g) in enumerate(gens)
        G[:, i] = g
    end
    return G
end

# iterator that wraps the generator matrix
struct FallbackGeneratorIterator{M<:AbstractMatrix}
    G::M
    n_plus_one::Int

    FallbackGeneratorIterator(G::M) where {M<:AbstractMatrix} =
        new{M}(G, size(G, 2) + 1)
end

Base.length(it::FallbackGeneratorIterator) = it.n_plus_one - 1

Base.eltype(::Type{<:FallbackGeneratorIterator{<:AbstractMatrix{N}}}) where {N} =
    AbstractVector{N}

Base.eltype(::Type{<:FallbackGeneratorIterator{<:Matrix{N}}}) where {N} =
    Vector{N}

Base.eltype(::Type{<:FallbackGeneratorIterator{<:SparseMatrixCSC{N}}}) where {N} =
    SparseVector{N}

function Base.iterate(it::FallbackGeneratorIterator, state::Int=1)
    if state == it.n_plus_one
        return nothing
    end
    g = it.G[:, state]
    state += 1
    return (g, state)
end

"""
    generators_fallback(Z::AbstractZonotope{N}) where {N<:Real}

Fallback definition of `generators` for zonotopic sets.

### Input

- `Z` -- zonotopic set

### Output

An iterator over the generators of `Z`.
"""
function generators_fallback(Z::AbstractZonotope{N}) where {N<:Real}
    return FallbackGeneratorIterator(genmat(Z))
end

"""
    togrep(Z::AbstractZonotope)

Return a generator representation of a zonotopic set.

### Input

- `Z` -- zonotopic set

### Output

The same set in generator representation.
This fallback implementation returns a `Zonotope`; however, more specific
implementations may return other generator representations.
"""
function togrep(Z::AbstractZonotope)
    return convert(Zonotope, Z)
end

# --- common AbstractZonotope functions ---

"""
    ngens(Z::AbstractZonotope)

Return the number of generators of a zonotopic set.

### Input

- `Z` -- zonotopic set

### Output

An integer representing the number of generators.
"""
function ngens(Z::AbstractZonotope)
    return length(generators(Z))
end

"""
    order(Z::AbstractZonotope)

Return the order of a zonotope.

### Input

- `Z` -- zonotope

### Output

A rational number representing the order of the zonotope.

### Notes

The order of a zonotope is defined as the quotient of its number of generators
and its dimension.
"""
function order(Z::AbstractZonotope)
    return ngens(Z) // dim(Z)
end

# --- LazySet interface functions ---

"""
    ρ(d::AbstractVector{N}, Z::AbstractZonotope{N}) where {N<:Real}

Return the support function of a zonotopic set in a given direction.

### Input

- `d` -- direction
- `Z` -- zonotopic set

### Output

The support function of the zonotopic set in the given direction.

### Algorithm

The support value is ``cᵀ d + ‖Gᵀ d‖₁`` where ``c`` is the center and ``G`` is
the generator matrix of `Z`.

"""
function ρ(d::AbstractVector{N}, Z::AbstractZonotope{N}) where {N<:Real}
    return dot(center(Z), d) + sum(abs.(transpose(genmat(Z)) * d))
end

"""
    σ(d::AbstractVector{N}, Z::AbstractZonotope{N}) where {N<:Real}

Return the support vector of a zonotopic set in a given direction.

### Input

- `d` -- direction
- `Z` -- zonotopic set

### Output

A support vector in the given direction.
If the direction has norm zero, the vertex with ``ξ_i = 1 \\ \\ ∀ i = 1,…, p``
is returned.
"""
function σ(d::AbstractVector{N}, Z::AbstractZonotope{N}) where {N<:Real}
    G = genmat(Z)
    return center(Z) .+ G * sign_cadlag.(_At_mul_B(G, d))
end

"""
    ∈(x::AbstractVector{N}, Z::AbstractZonotope{N};
      solver=default_lp_solver(N)) where {N<:Real}

Check whether a given point is contained in a zonotopic set.

### Input

- `x`      -- point/vector
- `Z`      -- zonotopic set
- `solver` -- (optional, default: `default_lp_solver(N)`) the backend used to
              solve the linear program

### Output

`true` iff ``x ∈ Z``.

### Examples

```jldoctest
julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1]);

julia> [1.0, 0.2] ∈ Z
false
julia> [1.0, 0.1] ∈ Z
true
```

### Algorithm

The membership problem is computed by stating and solving the following linear
program.
Let ``p`` and ``n`` be the number of generators and ambient dimension,
respectively.
We consider the minimization of ``x_0`` in the ``p+1``-dimensional space of
elements ``(x_0, ξ_1, …, ξ_p)`` constrained to ``0 ≤ x_0 ≤ ∞``,
``ξ_i ∈ [-1, 1]`` for all ``i = 1, …, p``, and such that ``x-c = Gξ`` holds.
If a feasible solution exists, the optimal value ``x_0 = 0`` is achieved.
"""
function ∈(x::AbstractVector{N}, Z::AbstractZonotope{N};
           solver=default_lp_solver(N)) where {N<:Real}
    @assert length(x) == dim(Z)

    p, n = ngens(Z), dim(Z)
    # (n+1) x (p+1) matrix with block-diagonal blocks 1 and genmat(Z)
    A = [[one(N); zeros(N, p)]'; [zeros(N, n) genmat(Z)]]
    b = [zero(N); (x - center(Z))]
    lbounds = [zero(N); fill(-one(N), p)]
    ubounds = [N(Inf); ones(N, p)]
    sense = ['>'; fill('=', n)]
    obj = [one(N); zeros(N, p)]

    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
    return (lp.status == :Optimal) # Infeasible or Unbounded => false
end

"""
    linear_map(M::AbstractMatrix{N}, Z::AbstractZonotope{N}) where {N<:Real}

Concrete linear map of a zonotopic set.

### Input

- `M` -- matrix
- `Z` -- zonotopic set

### Output

The zonotope obtained by applying the linear map to the center and generators
of ``Z``.
"""
function linear_map(M::AbstractMatrix{N}, Z::AbstractZonotope{N}
                   ) where {N<:Real}
    @assert dim(Z) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
                                 "applied to a set of dimension $(dim(Z))"

    c = M * center(Z)
    gi = M * genmat(Z)
    return Zonotope(c, gi)
end

"""
    translate(Z::AbstractZonotope{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:Real}

Translate (i.e., shift) a zonotope by a given vector.

### Input

- `Z`     -- zonotope
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated zonotope.

### Notes

The generator matrix is shared with the original zonotope if `share == true`.

### Algorithm

We add the vector to the center of the zonotope.
"""
function translate(Z::AbstractZonotope{N}, v::AbstractVector{N};
                   share::Bool=false) where {N<:Real}
    @assert length(v) == dim(Z) "cannot translate a $(dim(Z))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    c = center(Z) + v
    G = share ? genmat(Z) : copy(genmat(Z))
    return Zonotope(c, G)
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(Z::AbstractZonotope{N}; [apply_convex_hull]::Bool=true
                 ) where {N<:Real}

Return the vertices of a zonotopic set.

### Input

- `Z`                 -- zonotopic set
- `apply_convex_hull` -- (optional, default: `true`) if `true`, post-process the
                         computation with the convex hull of the points

### Output

List of vertices as a vector of vectors.

### Algorithm

If the zonotopic set has ``p`` generators, each vertex is the result of summing
the center with some linear combination of generators, where the combination
factors are ``ξ_i ∈ \\{-1, 1\\}``.

There are at most ``2^p`` distinct vertices. Use the flag `apply_convex_hull` to
control whether a convex hull algorithm is applied to the vertices computed by
this method; otherwise, redundant vertices may be present.
"""
function vertices_list(Z::AbstractZonotope{N};
                       apply_convex_hull::Bool=true) where {N<:Real}
    p = ngens(Z)
    if p == 0
        return [center(Z)]
    end

    vlist = Vector{Vector{N}}()
    sizehint!(vlist, 2^p)
    G = genmat(Z)

    for ξi in Iterators.product([[1, -1] for i = 1:p]...)
        push!(vlist, center(Z) .+ G * collect(ξi))
    end

    return apply_convex_hull ? convex_hull!(vlist) : vlist
end

"""
    constraints_list(P::AbstractZonotope{N}) where {N<:Real}

Return the list of constraints defining a zonotopic set.

### Input

- `Z` -- zonotopic set

### Output

The list of constraints of the zonotopic set.

### Algorithm

This is the (inefficient) fallback implementation for rational numbers.
It first computes the vertices and then converts the corresponding polytope
to constraint representation.
"""
function constraints_list(Z::AbstractZonotope{N}) where {N<:Real}
    return constraints_list(VPolytope(vertices_list(Z)))
end

"""
    constraints_list(Z::AbstractZonotope{N}; check_full_rank::Bool=true
                    ) where {N<:AbstractFloat}

Return the list of constraints defining a zonotopic set.

### Input

- `Z`               -- zonotopic set
- `check_full_rank` -- (optional; default: `true`) flag for checking whether the
                       generator matrix has full rank

### Output

The list of constraints of the zonotopic set.

### Notes

The algorithm assumes that no generator is redundant.
The result has ``2 \\binom{p}{n-1}`` (with ``p`` being the number of generators
and ``n`` being the ambient dimension) constraints, which is optimal under this
assumption.

If ``p < n`` or the generator matrix is not full rank, we fall back to the
(slower) computation based on the vertex representation.

### Algorithm

We follow the algorithm presented in *Althoff, Stursberg, Buss: Computing
Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes.
2009.*

The one-dimensional case is not covered by that algorithm; we manually handle
this case, assuming that there is only one generator.
"""
function constraints_list(Z::AbstractZonotope{N}; check_full_rank::Bool=true
                         ) where {N<:AbstractFloat}
    G = genmat(Z)
    p = ngens(Z)
    n = dim(Z)

    # use fallback implementation if order < 1 or matrix is not full rank
    if p < n || (check_full_rank && rank(G) < n)
        return invoke(constraints_list, Tuple{AbstractZonotope{<:Real}}, Z)
    end

    # special handling of 1D case
    if n == 1
        if p > 1
            error("1D-zonotope constraints currently only support a single " *
                  "generator")
        end

        c = center(Z)[1]
        g = G[:, 1][1]
        constraints = [LinearConstraint([N(1)], c + g),
                       LinearConstraint([N(-1)], g - c)]
        return constraints
    end

    c = center(Z)
    m = binomial(p, n - 1)
    constraints = Vector{LinearConstraint{N, Vector{N}}}()
    sizehint!(constraints, 2m)
    for columns in StrictlyIncreasingIndices(p, n-1)
        c⁺ = cross_product(view(G, :, columns))
        iszero(c⁺) && continue
        normalize!(c⁺, 2)

        Δd = sum(abs.(transpose(G) * c⁺))

        d⁺ = dot(c⁺, c) + Δd
        c⁻ = -c⁺
        d⁻ = -d⁺ + 2 * Δd  # identical to dot(c⁻, c) + Δd

        push!(constraints, LinearConstraint(c⁺, d⁺))
        push!(constraints, LinearConstraint(c⁻, d⁻))
    end
    return constraints
end


text:
#######################
## Prepare libraries ##
#######################

# export JULIA_NUM_THREADS=4

BASE_FOLDER = pwd()
(basename(BASE_FOLDER) == "examples") && (BASE_FOLDER = dirname(BASE_FOLDER))
using Pkg
Pkg.activate(BASE_FOLDER)
# Install all required packages if not present
isfile(joinpath(BASE_FOLDER,"Manifest.toml")) || Pkg.instantiate()

using Combinatorics, Distributions, LightGraphs, LinearAlgebra, QuantEcon
using Roots, SymPy, Manifolds
using Nash

## Create game from provided payoff matrices
# Player 1 payoff matrix is
    # 1  0
    # 0  1
# Player 2 payoff matrix is
    # 1  0
    # 1  1
generate_game([1 0; 0 1], [1 0; 1 1])

## Generate random game for 2 players.
# First player has 2 actions and second player has 3 actions
# Payoffs are sampled from standard normal Distributions

random_2players_game(Normal(0,1),2,3)

## Generate random game for 3 players.
# First player has 2 actions, second player has 3 actions, third player has 3 actions
# Payoffs are sampled from Binomial distribution with 5 trials and p = 0.5

random_nplayers_game(Binomial(5,0.5),[2,3,3])

## Get payoffs for each player from given mixed strategy profiles
# Player 1 plays first and second action with p = 0.5
# Player 2 plays first action with p = 0.75 and second action with p = 0.25

game = random_nplayers_game(Binomial(5,0.5),[2,2]);
s = [[0.5,0.5],[0.75,0.25]]
s[1]'*game["player2"]*s[2]
get_payoff(game,s)

## Return best reply of first (k = 1) player
# with perturbation epsil=0.5
best_reply(generate_game([1 0; 0 1], [1 0; 0 1]), [[1.0, 0.], [0, 1]], 1, 0.5)
# no perturbation, with convex hull
a2 = best_reply(generate_game(Matrix(I,3,3), Matrix(I,3,3)), [[1/2,1/2,0],[1/3,1/3,1/3]], 1, return_val = "chull")
# plot best responses
plot_br(a2)

# Check if passed strategy is best response for each player and is it a Nash equilibrium
is_nash_q(generate_game([1 0 ; 0 1], [1 0; 0 1]), [[1, 0], [1, 0]])

# Iteratively find best reply with perturbation (epsil=1/3) and limit of 10 iterations
game = generate_game(Matrix(I, 3, 3), Matrix(I, 3, 3))
s = [[1/2,1/2, 0], [1/3,1/3,1/3]]
game_history = iterate_best_reply(game, s, 1/3, 10)

# Transform game into Markov Chain
game = generate_game(Matrix(1I, 3, 3), Matrix(I, 3, 3))
s = [[1, 0, 0],[0, 1, 0]]
mc = game2markov(game, s)
simulate(mc, 9, init = 3)
simulate(mc, 9) # random initial condition

is_irreducible(mc) # can we reach any state from any other state?
communication_classes(mc)
period(mc)
is_aperiodic(mc)
stationary_distributions(mc)

plot_markov(10, mc)

# Generate random symmetric game
game = random_symmetric_2players_game(Binomial(10,1/2),2)

# Find necessary values to find symmetric nash equilibrium
x = symbols("x", real=true)
find_symmetric_nash_equilibrium_2players_game(game, [x, 1-x])

# Create equations for replicator
create_replicator_eqs(game, [x, 1-x])

# Generate graph of relations that need to hold to make a game symmetric
actions_no_vector = [2, 2, 2];
graph = create_symmetries_graph(actions_no_vector)

# Check single relation
check_equality_condition(graph,((2,2,1),1), ((2,1,2),1))

# Find all needed relations
find_all_equalities(graph)

# counting vFunction - three players coordination game
pay = reshape([1 0 0 0 0 0 0 1], (2,2,2))
g1 = generate_game(pay, pay, pay)
vFunction(g1, [[1, 0], [1, 0], [1, 0]])
vFunction(g1, [[1, 0], [1, 0], [0, 1]])

# Creates cartersian product of simplices
s1 = [[1, 0], [1, 0], [0, 1]]
v1 = cart_prod_simplices(s1)

s2 = [[1, 0], [1, 0], [1, 2, 3]]
v2 = cart_prod_simplices(s2)


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

struct RichardsonIterator{T, VT, KKT, LinSolver} <: AbstractIterator
    linear_solver::LinSolver
    kkt::KKT
    residual::VT
    max_iter::Int
    tol::T
    acceptable_tol::T
    logger::Logger
end
function RichardsonIterator(
    linear_solver::AbstractLinearSolver,
    kkt::AbstractKKTSystem,
    res::AbstractVector;
    max_iter=10, tol=1e-10, acceptable_tol=1e-5, logger=Logger(),
)
    return RichardsonIterator(
        linear_solver, kkt, res, max_iter, tol, acceptable_tol, logger,
    )
end

# Solve reduced KKT system. Require only the primal/dual values.
function solve_refine!(
    x::AbstractKKTVector{T, VT},
    solver::RichardsonIterator{T, VT, KKT},
    b::AbstractKKTVector{T, VT},
) where {T, VT, KKT<:AbstractReducedKKTSystem}
    solve_refine!(primal_dual(x), solver, primal_dual(b))
end

# Solve unreduced KKT system. Require UnreducedKKTVector as inputs.
function solve_refine!(
    x::UnreducedKKTVector{T, VT},
    solver::RichardsonIterator{T, VT, KKT},
    b::UnreducedKKTVector{T, VT},
) where {T, VT, KKT<:AbstractUnreducedKKTSystem}
    solve_refine!(full(x), solver, full(b))
end

function solve_refine!(
    x::AbstractVector{T},
    solver::RichardsonIterator{T},
    b::AbstractVector{T},
) where T
    @debug(solver.logger, "Iterative solver initiated")

    ε = solver.residual
    norm_b = norm(b, Inf)

    fill!(x, zero(T))
    fill!(ε, zero(T))

    ε = solver.residual
    axpy!(-1, b, ε)
    norm_res = norm(ε, Inf)
    residual_ratio = norm_res / (one(T) + norm_b)

    iter = 0
    residual_ratio_old = Inf
    noprogress = 0

    while true
        mod(iter, 10)==0 &&
            @debug(solver.logger,"iter ||res||")
        @debug(solver.logger, @sprintf("%4i %6.2e", iter, residual_ratio))
        iter += 1
        if (iter > solver.max_iter) || (residual_ratio < solver.tol)
            break
        end

        solve!(solver.linear_solver, ε)
        axpy!(-1, ε, x)
        mul!(ε, solver.kkt, x)
        axpy!(-1, b, ε)
        norm_res = norm(ε, Inf)

        residual_ratio_old = residual_ratio
        residual_ratio = norm_res / (one(T)+norm_b)
    end

    @debug(solver.logger, @sprintf(
        "Iterative solver terminated with %4i refinement steps and residual = %6.2e",
        iter, residual_ratio),
    )

    if residual_ratio < solver.acceptable_tol
        return :Solved
    else
        return :Failed
    end
end


text:
function calc_norm(array :: Matrix{Float64})
  sqrt(sumabs2(array) / size(array, 1))
end

function calc_norm(vector :: Array{Float64, 1})
  sqrt(sumabs2(vector))
end

function restriction(vector)
  vector / (calc_norm(vector) + 0.001)
end

function init(a...)
    restriction(randn(a...))
end

const DefaultRate = 0.01

immutable Class
    class_name :: String
    type_len :: Int64
end

function Base.show(io :: IO, m :: Class)
    print(io, string("Class(", m.class_name, ")"))
end

immutable Object
    obj_name :: String
    class :: Class
    value :: Array{Float64, 1}
end

Object(obj_name :: String, class :: Class) = Object(obj_name, class, init(class.type_len))

function Base.show(io :: IO, m :: Object)
    print(io, string("Obj(", m.obj_name, "){", m.class.class_name, "}"))
end

type DFunction
    func_name :: String
    in_classes :: Array{Class, 1}
    out_class :: Class
    f_matrix :: Matrix{Float64}
    f_transformation :: Transformation
    learning_rate :: Float64
end

function Base.show(io :: IO, m :: DFunction)
    in_names = join([c.class_name for c in m.in_classes], ",")
    print(io, string("Func(", m.func_name, "){", in_names, "->", m.out_class.class_name, "}"))
end

function LFunction(func_name, in_classes :: Array{Class, 1}, out_class :: Class, learning_rate = DefaultRate)
    in_len = sum([class.type_len for class in in_classes])
    out_len = out_class.type_len
    DFunction(func_name, in_classes, out_class, init(out_len, in_len), Identity, learning_rate)
end

function SFunction(func_name, in_classes :: Array{Class, 1}, out_class :: Class, learning_rate = DefaultRate)
    in_len = sum([class.type_len for class in in_classes])
    out_len = out_class.type_len
    DFunction(func_name, in_classes, out_class, init(out_len, in_len), Sigmoid, learning_rate)
end

function DFunction(func_name, in_classes :: Array{Class, 1}, out_class :: Class, f_transformation = Sigmoid, learning_rate = DefaultRate)
    in_len = sum([class.type_len for class in in_classes])
    out_len = out_class.type_len
    DFunction(func_name, in_classes, out_class, init(out_len, in_len), f_transformation, learning_rate)
end

function CoFunction(func_name, in_class :: Class, f_transformation = Sigmoid)
    DFunction(func_name, [in_class], in_class, diagm(ones(in_class.type_len)), f_transformation, 0.0)
end

function apply(func :: DFunction, objs :: Array{Object, 1})
    inputs = vcat([obj.value for obj in objs]...)
    outputs = eval(func.f_transformation).(func.f_matrix * inputs)
    Object("untitled", func.out_class, outputs)
end

function apply!(func :: DFunction, in_objs :: Array{Object, 1}, out_obj :: Object)
    inputs = vcat([obj.value for obj in in_objs]...)
    out_obj.value .= eval(func.f_transformation).(func.f_matrix * inputs)
end

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
using JuLIP: stress
using JuLIP.Potentials

println("-------------------------------------------------")
println("   Variable Cell Test")
println("-------------------------------------------------")
calc = lennardjones(r0=rnn("Al"))
at = bulk("Al") * 2      # cubic=true,
set_pbc!(at, true)
set_calculator!(at, calc)

# perturb the cell shape
set_defm!(at, defm(at) + 0.1*rand(JMatF), updatepositions=true)
# check that under this deformation the atom positions are still
# in a homogeneous lattice (i.e. they are deformed correctly with the cell)
@test maxnorm(forces(at)) < 1e-12
# now perturb the atom positions as well
rattle!(at, 0.1)

# set the constraint >>> this means the deformed cell defines F0 and X0
set_constraint!(at, ExpVariableCell(at))

# check that energy, forces, stress, dofs, gradient evaluate at all
print("check that energy, forces, stress, dofs, gradient evaluate ... ")
energy(at)
forces(at)
gradient(at)
stress(at)
println("ok.")
@test true

print("Check that setting and getting dofs is consistent ... ")
x = dofs(at)
@assert length(x) == 3 * length(at) + 6
y = x + 0.01 * rand(size(x))
set_dofs!(at, y)
z = dofs(at)
@test vecnorm(y - z, Inf) < 1e-12
println("ok")

# perform the finite-difference test
set_dofs!(at, x)
JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1)

println("-------------------------------------------------")
println("Test optimisation with ExpVariableCell")
# # start with a clean `at`
at = bulk("Si") * 2   # cubic=true,
set_pbc!(at, true)
set_calculator!(at, StillingerWeber())
set_constraint!(at, ExpVariableCell(at))
#
println("For the initial state, stress is far from 0:")
@show vecnorm(virial(at), Inf)
JuLIP.Solve.minimise!(at, verbose=2)
# dt = 1e-6
# for n = 1:100
#    x = dofs(at)
#    ∇E = gradient(at)
#    @printf(" %3d  |  %4.2e \n", n, vecnorm(∇E, Inf))
#    x -= dt * ∇E
#    set_dofs!(at, x)
# end
println("After optimisation, stress is 0:")
@show vecnorm(virial(at), Inf)
# @test vecnorm(virial(at), Inf) < 1e-4

println("Test for comparison with the standard VariableCell")
at = bulk("Si") * 2   # cubic=true,
set_pbc!(at, true)
set_calculator!(at, StillingerWeber())
set_constraint!(at, VariableCell(at))
#
println("For the initial state, stress is far from 0:")
@show vecnorm(virial(at), Inf)
JuLIP.Solve.minimise!(at, verbose=2)

# println("-------------------------------------------------")
# println("And now with pressure . . .")
# println("-------------------------------------------------")
# set_constraint!(at, ExpVariableCell(at, pressure=10.0123))
# JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1)
# at = Atoms("Al", pbc=(true,true,true)) * 2   # cubic=true,
# set_calculator!(at, calc)
# set_constraint!(at, VariableCell(at, pressure=0.01))
# JuLIP.Solve.minimise!(at)
# @show vecnorm(stress(at), Inf)
# @test vecnorm(gradient(at), Inf) < 1e-4


text:
function is_jll_name(name::AbstractString)::Bool
    return endswith(name, "_jll")
end

function _get_all_dependencies_nonrecursive(working_directory::AbstractString,
                                            pkg,
                                            version)
    all_dependencies = Vector{String}(undef, 0)
    deps = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Deps.toml"))
    for version_range in keys(deps)
        if version in Pkg.Types.VersionRange(version_range)
            for name in keys(deps[version_range])
                push!(all_dependencies, name)
            end
        end
    end
    unique!(all_dependencies)
    return all_dependencies
end

function meets_allowed_jll_nonrecursive_dependencies(working_directory::AbstractString,
                                                     pkg,
                                                     version)
    # If you are a JLL package, you are only allowed to have five kinds of dependencies:
    # 1. Pkg
    # 2. Libdl
    # 3. Artifacts
    # 4. JLLWrappers
    # 5. other JLL packages
    all_dependencies = _get_all_dependencies_nonrecursive(working_directory,
                                                          pkg,
                                                          version)
    for dep in all_dependencies
        if !((dep == "Pkg") | (dep == "Libdl") | (dep == "Artifacts") | (dep == "JLLWrappers") | (is_jll_name(dep)))
            return false, "JLL packages are only allowed to depend on Pkg, Libdl, Artifacts, JLLWrappers and other JLL packages"
        end
    end
    return true, ""
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

type OSMattributes
    oneway::Bool
    oneway_override::Bool
    oneway_reverse::Bool
    visible::Bool
    lanes::Int

    name::UTF8String
    class::UTF8String
    detail::UTF8String
    cycleway::UTF8String
    sidewalk::UTF8String
    bicycle::UTF8String

    # XML elements
    element::Symbol # :None, :Node, :Way, :Tag[, :Relation]
    parent::Symbol # :Building, :Feature, :Highway
    way_nodes::Vector{Int} # for buildings and highways

    id::Int # Uninitialized
    lat::Float64 # Uninitialized
    lon::Float64 # Uninitialized

    OSMattributes() = new(false,false,false,false,1,
                          "","","","","","",:None,:None,[])
end

type OSMdata
    nodes::Dict{Int,LLA}
    highways::Dict{Int,Highway}
    buildings::Dict{Int,Building}
    features::Dict{Int,Feature}
    attr::OSMattributes
    OSMdata() = new(Dict(),Dict(),Dict(),Dict(),OSMattributes())
end

function reset_attributes!(osm::OSMattributes)
    osm.oneway = osm.oneway_override = osm.oneway_reverse = osm.visible = false
    osm.lanes = 1
    osm.name = osm.class = osm.detail = osm.cycleway = osm.sidewalk = osm.bicycle = ""
    osm.element = osm.parent = :None
    empty!(osm.way_nodes)
end

### PARSE XML ELEMENTS ###

function parse_node(attr::OSMattributes, attrs_in::Dict{@compat(AbstractString),@compat(AbstractString)})
    attr.visible = true
    attr.element = :Node
    if haskey(attrs_in, "id")
        attr.id = @compat( parse(Int,attrs_in["id"]) )
        attr.lat = float(attrs_in["lat"])
        attr.lon = float(attrs_in["lon"])
    end
end

function parse_way(attr::OSMattributes, attrs_in::Dict{@compat(AbstractString),@compat(AbstractString)})
    attr.visible = true
    attr.element = :Way
    if haskey(attrs_in, "id")
        attr.id = @compat( parse(Int,attrs_in["id"]) )
    end
end

function parse_nd(attr::OSMattributes, attrs_in::Dict{@compat(AbstractString),@compat(AbstractString)})
    if haskey(attrs_in, "ref")
        push!(attr.way_nodes, @compat( parse(Int64,attrs_in["ref"]) ) )
    end
end

function parse_tag(attr::OSMattributes, attrs_in::Dict{@compat(AbstractString),@compat(AbstractString)})
    if haskey(attrs_in, "k") && haskey(attrs_in, "v")
        k, v = attrs_in["k"], attrs_in["v"]
        if k == "name"
            if isempty(attr.name)
                attr.name = v # applicable to roads (highways), buildings, features
            end
        elseif attr.element == :Way
            if k == "building"
                parse_building(attr, v)
            else
                parse_highway(attr, k, v) # for other highway tags
            end
        elseif attr.element == :Node
            if haskey(FEATURE_CLASSES, k)
                parse_feature(attr, k, v)
            end
        end
    else
        # Nothing to be done here?
    end
end

### PARSE OSM ENTITIES ###

function parse_highway(attr::OSMattributes, k::@compat(AbstractString), v::@compat(AbstractString))
    if k == "highway"
        attr.class = v
        if v == "services" # Highways marked "services" are not traversable
            attr.visible = false
            return
        end
        if v == "motorway" || v == "motorway_link"
            attr.oneway = true # motorways default to oneway
        end
    elseif k == "oneway"
        if v == "-1"
            attr.oneway = true
            attr.oneway_reverse = true
        elseif v == "false" || v == "no" || v == "0"
            attr.oneway = false
            attr.oneway_override = true
        elseif v == "true" || v == "yes" || v == "1"
            attr.oneway = true
        end
    elseif k == "junction" && v == "roundabout"
        attr.oneway = true
    elseif k == "cycleway"
        attr.cycleway = v
    elseif k == "sidewalk"
        attr.sidewalk = v
    elseif k == "bicycle"
        attr.bicycle = v
    elseif k == "lanes" && length(v)==1 && '1' <= v[1] <= '9'
        attr.lanes = @compat parse(Int,v)
    else
        return
    end
    attr.parent = :Highway
end

function parse_building(attr::OSMattributes, v::@compat(AbstractString))
    attr.parent = :Building
    if isempty(attr.class)
        attr.class = v
    end
end

function parse_feature(attr::OSMattributes, k::@compat(AbstractString), v::@compat(AbstractString))
    attr.parent = :Feature
    attr.class = k
    attr.detail = v
end

### LibExpat.XPStreamHandlers ###

function parseElement(handler::LibExpat.XPStreamHandler, name::@compat(AbstractString), attrs_in::Dict{@compat(AbstractString),@compat(AbstractString)})
    attr = handler.data.attr::OSMattributes
    if attr.visible
        if name == "nd"
            parse_nd(attr, attrs_in)
        elseif name == "tag"
            parse_tag(attr, attrs_in)
        end
    elseif !(haskey(attrs_in, "visible") && attrs_in["visible"] == "false")
        if name == "node"
            parse_node(attr, attrs_in)
        elseif name == "way"
            parse_way(attr, attrs_in)
        end
    end # no work done for "relations" yet
end

function collectValues(handler::LibExpat.XPStreamHandler, name::@compat(AbstractString))
    # println(typeof(name))
    osm = handler.data::OSMdata
    attr = osm.attr::OSMattributes
    if name == "node"
        osm.nodes[attr.id] = LLA(attr.lat, attr.lon)
        if attr.parent == :Feature
            osm.features[attr.id] = Feature(attr.class, attr.detail, attr.name)
        end
    elseif name == "way"
        if attr.parent == :Building
            osm.buildings[attr.id] = Building(attr.class, attr.name, copy(attr.way_nodes))
        elseif attr.parent == :Highway
            if attr.oneway_reverse
                reverse!(attr.way_nodes)
            end
            osm.highways[attr.id] = Highway(attr.class, attr.lanes,
                                            (attr.oneway && !attr.oneway_override),
                                            attr.sidewalk, attr.cycleway, attr.bicycle,
                                            attr.name, copy(attr.way_nodes))
        end
    else # :Tag or :Nd (don't reset values!)
        return
    end
    reset_attributes!(osm.attr)
end

### Parse the data from an openStreetMap XML file ###
function parseMapXML(filename::@compat(AbstractString))

    # Parse the file
    street_map = LightXML.parse_file(filename)

    if LightXML.name(LightXML.root(street_map)) != "osm"
        throw(ArgumentError("Not an OpenStreetMap datafile."))
    end

    return street_map
end

function getOSMData(filename::@compat(AbstractString); args...)
    osm = OSMdata()

    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parseElement
    callbacks.end_element = collectValues

    LibExpat.parsefile(filename, callbacks, data=osm; args...)
    osm.nodes, osm.highways, osm.buildings, osm.features
end


text:
using JuliaGridSynth
using Test
using TestSetExtensions

@testset "JuliaGridSynth.jl" begin
    # Write your tests here.

end

max_stars_repo_name: isolatedinformation/JuliaGridSynth

max_issues_repo_name: isolatedinformation/JuliaGridSynth

max_forks_repo_name: isolatedinformation/JuliaGridSynth


text:
# # Implementations and proof of concept
#
# In this examples I want to show some of the implementations and how they are used
# as a proof of concept. This means that we will use the "low-level API", i.e. calling
# the methods directly instead of the `optimize` interface.
#
# Further, we wish to show that the implementations can at least solve some of the
# most common benchmark optimization problems.
#
# Before we start, I will define a seed and an RNG to enable reproducibility of the
# results presented here.

using Random

RANDOM_SEED = 458012;
rng = MersenneTwister(RANDOM_SEED);

#
# ## Nonlinear ``d``-dimensional global optimization problem
#
# Using `BioMetaheuristics.jl` is fairly straightforward, we will start by defining
# an d-dimensional nonlinear function to minimize,
# in this case we will use a popular function, the
# [Griewank function](http://mathworld.wolfram.com/GriewankFunction.html)
# defined as
#
# ```math
# f(\mathbf{x}) = \sum_{i=1}^d \frac{x_i^2}{4000} - \prod_{i=1}^d \cos{\left(
# \frac{x_i}{\sqrt{i}}\right)} + 1
# ```
#
# where ``d`` is the dimension of the problem. It's mostly evaluated within the
# boundaries ``-100 \leq x_i \leq 100``, and it has a **minimum** at ``\mathbf
# {x^*} = (0, \cdots, 0)``, and it evaluates to ``f(\mathbf{x^*}) = 0``.
#
# We define the function in `Julia` like this

function griewank(x)
    first_term = sum(x.^2) / 4000
    ## This variable will hold the result of the product,
    ## the second term in the function definition from above
    second_term = 1.0
    for (idx, val) in enumerate(x)
        second_term *= cos(val / sqrt(idx))
    end

    return first_term - second_term + 1.0
end

# Now, we wish to find the minimum of this function, and fortunately we know the
# true value so we can compare it later, we can use some of the implementations
# from `BioMetaheuristics.jl`, for example, [`PSO`](@ref).
# In this script we have chosen 30 particles within the population, `d` is equal
# to 20, next we define the boundaries and finally we declare that the algorithm
# will run for 20000 maximum iterations until it stops, having _converged_.

using BioMetaheuristics

val = PSO(
    griewank,
    Population(35, 10, -600.0, 600.0, rng),
    20_000,
    rng
)
println(val)

# Within a certain tolerance of about ``\epsilon = 1 \times 10^{-6}`` we have found
# the _global_ minimum of the function. We can actually check the value with the
# evaluation, notice that it actually returns `0`, as expected.

griewank(val.x)

# ## Nonlinear 2-dimensional global optimization problem
#
# Let us now tackle one of the most common optimization problems, which is
# finding the minimum of the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function)
# which is a non-convex function, meaning that is does not have just one minimum
# or stationary point, it has several, so it is a difficult problem for classical
# optimization algorithms. In this example we will try to solve it using the
# [`SimulatedAnnealing`](@ref) implementation from `BioMetaheuristics.jl`.
#
# First, we define the Rosenbrock function in `Julia`

rosenbrock2d(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2;

# We will apply the _Simulated Annealing_ algorithm to find the global optimum
val = SimulatedAnnealing(
    rosenbrock2d, -5.0, 5.0, 2, rng; low_temp=10_000
)
println(val)

# Again, within a certain tolerance we find the expected result which is
#
# ```math
# \mathbf{x}^{*} = (1, 1)
# ```
#
# and if we account for rounding errors and floating-point arithmetic, we
# can safely take this result as the best.


text:
@mix @columns struct Catabolism{MoMoD}
    k::MoMoD  | 0.2 | mol*mol^-1*d^-1 | Beta(2.0, 2.0)  | [0.0,1.0]   | _  | "Reserve turnover rate"
end

@mix @columns struct CatabolismCN{MoMoD}
    kC::MoMoD | 0.2 | mol*mol^-1*d^-1 | Beta(2.0, 2.0)  | [0.0,1.0]   | _  | "C-reserve turnover rate"
    kN::MoMoD | 0.2 | mol*mol^-1*d^-1 | Beta(2.0, 2.0)  | [0.0,1.0]   | _  | "N-reserve turnover rate"
end

abstract type AbstractCatabolism end
abstract type AbstractCatabolismCN <: AbstractCatabolism end
abstract type AbstractCatabolismCNE <: AbstractCatabolism end
@Catabolism struct Catabolism{} <: AbstractCatabolism end
@Catabolism struct CatabolismCNshared{} <: AbstractCatabolismCN end
@CatabolismCN struct CatabolismCN{} <: AbstractCatabolismCN end
@Catabolism @CatabolismCN struct CatabolismCNE{} <: AbstractCatabolismCNE end

# κEC::F             | 0.3             | _               | Beta(2.0, 2.0)  | [0.0,1.0]  | _ | "Non-processed C-reserve returned to C-reserve"
# κEN::F             | 0.3             | _               | Beta(2.0, 2.0)  | [0.0,1.0]  | _ | "Non-processed N-reserve returned to N-reserve"

kC(p::CatabolismCNshared) = p.k
kN(p::CatabolismCNshared) = p.k
kC(p::AbstractCatabolism) = p.kC
kN(p::AbstractCatabolism) = p.kN
kE(p::AbstractCatabolism) = p.k

"""
    catabolism!(o, u, t::Number)
Catabolism for E, C and N, or C, N and E reserves.
Does not alter flux in J - operates only on J1 (intermediate storage)
"""
catabolism!(o, u) = catabolism!(catabolism_pars(o), o, u)
catabolism!(p::AbstractCatabolismCNE, o, u) = begin
    v, J, J1 = vars(o), flux(o), flux1(o)
    turnover = (kC(p), kN(p), kE(p)) .* tempcorrection(v) .* shape(v)
    reserve = (u.C, u.N, u.E)
    rel_reserve = reserve ./ u.V
    corr_j_E_mai = j_E_mai(o) * tempcorrection(v)

    r = calc_rate(rate_formula(o), su_pars(o), rel_reserve, turnover, corr_j_E_mai, y_E_EC(o), y_E_EN(o), y_V_E(o), κsoma(o))
    set_rate!(o, r)
    r < zero(r) && return false

    J1[:C,:ctb], J1[:N,:ctb], J1[:EE,:ctb] = non_growth_flux.(reserve, turnover, r)
    J1[:C,:rej], J1[:N,:rej], J1[:CN,:ctb] = stoich_merge(su_pars(o), J1[:C,:ctb], J1[:N,:ctb], y_E_EC(o), y_E_EN(o))

    J1[:E,:ctb] = J1[:EE,:ctb] + J1[:CN,:ctb] # Total catabolic flux
    set_θE!(o, J1[:EE,:ctb]/J1[:E,:ctb]) # Proportion of general reserve flux in total catabolic flux
    true
end
catabolism!(p::AbstractCatabolismCN, o, u) = begin
    v, J, J1 = vars(o), flux(o), flux1(o)
    turnover = (kC(p), kN(p)) .* tempcorrection(v) .* shape(v)
    reserve = (u.C, u.N)
    rel_reserve = reserve ./ u.V
    corr_j_E_mai = j_E_mai(o) * tempcorrection(v)

    r, found = calc_rate(rate_formula(o), su_pars(o), rel_reserve, turnover, corr_j_E_mai, y_E_EC(o), y_E_EN(o), y_V_E(o), κsoma(o))
    if !found
        println("Root for rate not found at t - $(tstep(o))")
        return false # dead
    elseif r < zero(r)
        println("Rate is less than zero at t - $(tstep(o))")
        return false # dead
    else
        set_rate!(o, r)
    end

    J1[:C,:ctb], J1[:N,:ctb] = non_growth_flux.(reserve, turnover, r)
    J1[:C,:rej], J1[:N,:rej], J1[:E,:ctb] = stoich_merge(su_pars(o), J1[:C,:ctb], J1[:N,:ctb], y_E_EC(o), y_E_EN(o))
    true
end


text:
path = dirname(@__FILE__)
fred = CSV.read("$path/../reference/fred_160812.csv")
custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :n_anticipated_shocks     => Setting(:n_anticipated_shocks, 6))

m = Model990(custom_settings = custom_settings, testing = true)
m <= Setting(:rate_expectations_source, :ois)
@testset "Check FRED data is properly loaded" begin
    @test @test_matrix_approx_eq Matrix(fred[:,2:end]) Matrix(load_fred_data(m, end_date = date_mainsample_end(m),
                                                                             verbose = :none)[:,2:end])
end


text:
import Celeste: detect_sources
import Celeste.Coordinates: match_coordinates
using Celeste.SDSSIO

@testset "detection" begin

# Compare the initial catalog produced by detect_sources() versus the
# SDSS primary catalog.
@testset "detect_sources() vs SDSS catalog" begin
    images = SampleData.get_sdss_images(4263, 5, 119)
    sdss_catalog = SampleData.get_sdss_catalog(4263, 5, 119)

    catalog, _ = detect_sources(images)

    ra = [ce.pos[1] for ce in catalog]
    dec = [ce.pos[2] for ce in catalog]
    sdss_ra = [ce.pos[1] for ce in sdss_catalog]
    sdss_dec = [ce.pos[2] for ce in sdss_catalog]
    idx, dists = match_coordinates(ra, dec, sdss_ra, sdss_dec)

    # Test that there are a bunch of coordinates that match within 0.5 arcsec
    # (This is a basic sanity check, not a very strict test.)
    @test sum(dists .< 0.5/3600.) > 600
end

end


text:
# Digraph
mutable struct Digraph{T} <: AbstractNexus{T}
    nn::Int
    na::Int
    fadj::Dict{T,Dict{T,Dict{Int,Nothing}}}
    badj::Dict{T,Dict{T,Dict{Int,Nothing}}}
end

is_directed(::Digraph{T}) where {T} = true
Digraph{T}() where {T} = begin
    fadj = Dict{T,Dict{T,Dict{Int,Nothing}}}()
    badj = Dict{T,Dict{T,Dict{Int,Nothing}}}()
    return Digraph(0, 0, fadj, badj)
end

badj(g::Digraph) = g.badj

function add_arc!(G::Digraph{T}, u::T, v::T, key::Union{Int,Nothing} = nothing) where {T}
    haskey(G.fadj, u) || begin
        G.fadj[u] = Dict()
        G.badj[u] = Dict()
        G.nn += 1
    end
    haskey(G.badj, v) || begin
        G.badj[v] = Dict()
        G.fadj[v] = Dict()
        G.nn += 1
    end
    isnothing(key) && (key = new_arc_key(G, u, v))

    get!(G.fadj[u], v, Dict())
    get!(G.badj[v], u, Dict())

    haskey(G.fadj[u][v], key) || (G.na += 1)
    get!(G.fadj[u][v], key, nothing)
    get!(G.badj[v][u], key, nothing)
    return key
end

function new_arc_key(g::Digraph{T}, u::T, v::T) where {T}
    haskey(g.fadj[u], v) ? begin
        keylist = keys(g.fadj[u][v])
    end : return 1
    key = length(keylist)
    while key in keylist
        key += 1
    end # find lowest integer not used
    return key
end

function rem_arc!(G::Digraph{T}, u::T, v::T, key::Union{Int,Nothing} = nothing) where {T}
    haskey(G.fadj[u], v) ? d = G.fadj[u][v] : return false # no v on u
    isnothing(key) ? pop!(d) : haskey(d, key) ? delete!(d, key) : return false # rem key
    length(d) == 0 && begin
        delete!(G.fadj[u], v)
        delete!(G.badj[v], u)
    end
    G.na -= 1
    return true
end

function add_node!(g::Digraph{T}, u::T) where {T}
    !haskey(g.fadj, u) ? begin # doesn't have node, create new dict
        g.fadj[u] = Dict()
        g.badj[u] = Dict()
        g.nn += 1
        return u
    end : return nothing
end

function rem_node!(g::Digraph{T}, u::T) where {T}
    haskey(g.fadj, u) || return false

    for v in keys(g.fadj[u])
        g.na -= length(g.badj[v][u])
        delete!(g.badj[v], u)
    end
    delete!(g.fadj, u)

    for v in keys(g.badj[u])
        g.na -= length(g.fadj[v][u])
        delete!(g.fadj[v], u)
    end
    delete!(g.badj, u)

    g.nn -= 1
    return true
end


text:
const multiples_of_coefficients = (

# 1-10
    (  0,  0,  0,  0,  1),
    (  0,  0,  0,  0,  2),
    ( -2,  0,  2,  0,  1),
    (  2,  0, -2,  0,  0),
    ( -2,  0,  2,  0,  2),
    (  1, -1,  0, -1,  0),
    (  0, -2,  2, -2,  1),
    (  2,  0, -2,  0,  1),
    (  0,  0,  2, -2,  2),
    (  0,  1,  0,  0,  0),

    # 11-20
    (  0,  1,  2, -2,  2),
    (  0, -1,  2, -2,  2),
    (  0,  0,  2, -2,  1),
    (  2,  0,  0, -2,  0),
    (  0,  0,  2, -2,  0),
    (  0,  2,  0,  0,  0),
    (  0,  1,  0,  0,  1),
    (  0,  2,  2, -2,  2),
    (  0, -1,  0,  0,  1),
    ( -2,  0,  0,  2,  1),

    # 21-30
    (  0, -1,  2, -2,  1),
    (  2,  0,  0, -2,  1),
    (  0,  1,  2, -2,  1),
    (  1,  0,  0, -1,  0),
    (  2,  1,  0, -2,  0),
    (  0,  0, -2,  2,  1),
    (  0,  1, -2,  2,  0),
    (  0,  1,  0,  0,  2),
    ( -1,  0,  0,  1,  1),
    (  0,  1,  2, -2,  0),

    # 31-40
    (  0,  0,  2,  0,  2),
    (  1,  0,  0,  0,  0),
    (  0,  0,  2,  0,  1),
    (  1,  0,  2,  0,  2),
    (  1,  0,  0, -2,  0),
    ( -1,  0,  2,  0,  2),
    (  0,  0,  0,  2,  0),
    (  1,  0,  0,  0,  1),
    ( -1,  0,  0,  0,  1),
    ( -1,  0,  2,  2,  2),

    # 41-50
    (  1,  0,  2,  0,  1),
    (  0,  0,  2,  2,  2),
    (  2,  0,  0,  0,  0),
    (  1,  0,  2, -2,  2),
    (  2,  0,  2,  0,  2),
    (  0,  0,  2,  0,  0),
    ( -1,  0,  2,  0,  1),
    ( -1,  0,  0,  2,  1),
    (  1,  0,  0, -2,  1),
    ( -1,  0,  2,  2,  1),

    # 51-60
    (  1,  1,  0, -2,  0),
    (  0,  1,  2,  0,  2),
    (  0, -1,  2,  0,  2),
    (  1,  0,  2,  2,  2),
    (  1,  0,  0,  2,  0),
    (  2,  0,  2, -2,  2),
    (  0,  0,  0,  2,  1),
    (  0,  0,  2,  2,  1),
    (  1,  0,  2, -2,  1),
    (  0,  0,  0, -2,  1),

    # 61-70
    (  1, -1,  0,  0,  0),
    (  2,  0,  2,  0,  1),
    (  0,  1,  0, -2,  0),
    (  1,  0, -2,  0,  0),
    (  0,  0,  0,  1,  0),
    (  1,  1,  0,  0,  0),
    (  1,  0,  2,  0,  0),
    (  1, -1,  2,  0,  2),
    ( -1, -1,  2,  2,  2),
    ( -2,  0,  0,  0,  1),

    # 71-80
    (  3,  0,  2,  0,  2),
    (  0, -1,  2,  2,  2),
    (  1,  1,  2,  0,  2),
    ( -1,  0,  2, -2,  1),
    (  2,  0,  0,  0,  1),
    (  1,  0,  0,  0,  2),
    (  3,  0,  0,  0,  0),
    (  0,  0,  2,  1,  2),
    ( -1,  0,  0,  0,  2),
    (  1,  0,  0, -4,  0),

    # 81-90
    ( -2,  0,  2,  2,  2),
    ( -1,  0,  2,  4,  2),
    (  2,  0,  0, -4,  0),
    (  1,  1,  2, -2,  2),
    (  1,  0,  2,  2,  1),
    ( -2,  0,  2,  4,  2),
    ( -1,  0,  4,  0,  2),
    (  1, -1,  0, -2,  0),
    (  2,  0,  2, -2,  1),
    (  2,  0,  2,  2,  2),

    # 91-100
    (  1,  0,  0,  2,  1),
    (  0,  0,  4, -2,  2),
    (  3,  0,  2, -2,  2),
    (  1,  0,  2, -2,  0),
    (  0,  1,  2,  0,  1),
    ( -1, -1,  0,  2,  1),
    (  0,  0, -2,  0,  1),
    (  0,  0,  2, -1,  2),
    (  0,  1,  0,  2,  0),
    (  1,  0, -2, -2,  0),

    # 101-106
    (  0, -1,  2,  0,  1),
    (  1,  1,  0, -2,  1),
    (  1,  0, -2,  2,  0),
    (  2,  0,  0,  2,  0),
    (  0,  0,  2,  4,  2),
    (  0,  1,  0,  1,  0)
    )

const multiples_of_arguments = (

    (-171996.0, -174.2,  92025.0,    8.9 ),
    (   2062.0,    0.2,   -895.0,    0.5 ),
    (     46.0,    0.0,    -24.0,    0.0 ),
    (     11.0,    0.0,      0.0,    0.0 ),
    (     -3.0,    0.0,      1.0,    0.0 ),
    (     -3.0,    0.0,      0.0,    0.0 ),
    (     -2.0,    0.0,      1.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    ( -13187.0,   -1.6,   5736.0,   -3.1 ),
    (   1426.0,   -3.4,     54.0,   -0.1 ),
    (   -517.0,    1.2,    224.0,   -0.6 ),
    (    217.0,   -0.5,    -95.0,    0.3 ),
    (    129.0,    0.1,    -70.0,    0.0 ),
    (     48.0,    0.0,      1.0,    0.0 ),
    (    -22.0,    0.0,      0.0,    0.0 ),
    (     17.0,   -0.1,      0.0,    0.0 ),
    (    -15.0,    0.0,      9.0,    0.0 ),
    (    -16.0,    0.1,      7.0,    0.0 ),
    (    -12.0,    0.0,      6.0,    0.0 ),
    (     -6.0,    0.0,      3.0,    0.0 ),
    (     -5.0,    0.0,      3.0,    0.0 ),
    (      4.0,    0.0,     -2.0,    0.0 ),
    (      4.0,    0.0,     -2.0,    0.0 ),
    (     -4.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (  -2274.0,   -0.2,    977.0,   -0.5 ),
    (    712.0,    0.1,     -7.0,    0.0 ),
    (   -386.0,   -0.4,    200.0,    0.0 ),
    (   -301.0,    0.0,    129.0,   -0.1 ),
    (   -158.0,    0.0,     -1.0,    0.0 ),
    (    123.0,    0.0,    -53.0,    0.0 ),
    (     63.0,    0.0,     -2.0,    0.0 ),
    (     63.0,    0.1,    -33.0,    0.0 ),
    (    -58.0,   -0.1,     32.0,    0.0 ),
    (    -59.0,    0.0,     26.0,    0.0 ),
    (    -51.0,    0.0,     27.0,    0.0 ),
    (    -38.0,    0.0,     16.0,    0.0 ),
    (     29.0,    0.0,     -1.0,    0.0 ),
    (     29.0,    0.0,    -12.0,    0.0 ),
    (    -31.0,    0.0,     13.0,    0.0 ),
    (     26.0,    0.0,     -1.0,    0.0 ),
    (     21.0,    0.0,    -10.0,    0.0 ),
    (     16.0,    0.0,     -8.0,    0.0 ),
    (    -13.0,    0.0,      7.0,    0.0 ),
    (    -10.0,    0.0,      5.0,    0.0 ),
    (     -7.0,    0.0,      0.0,    0.0 ),
    (      7.0,    0.0,     -3.0,    0.0 ),
    (     -7.0,    0.0,      3.0,    0.0 ),
    (     -8.0,    0.0,      3.0,    0.0 ),
    (      6.0,    0.0,      0.0,    0.0 ),
    (      6.0,    0.0,     -3.0,    0.0 ),
    (     -6.0,    0.0,      3.0,    0.0 ),
    (     -7.0,    0.0,      3.0,    0.0 ),
    (      6.0,    0.0,     -3.0,    0.0 ),
    (     -5.0,    0.0,      3.0,    0.0 ),
    (      5.0,    0.0,      0.0,    0.0 ),
    (     -5.0,    0.0,      3.0,    0.0 ),
    (     -4.0,    0.0,      0.0,    0.0 ),
    (      4.0,    0.0,      0.0,    0.0 ),
    (     -4.0,    0.0,      0.0,    0.0 ),
    (     -3.0,    0.0,      0.0,    0.0 ),
    (      3.0,    0.0,      0.0,    0.0 ),
    (     -3.0,    0.0,      1.0,    0.0 ),
    (     -3.0,    0.0,      1.0,    0.0 ),
    (     -2.0,    0.0,      1.0,    0.0 ),
    (     -3.0,    0.0,      1.0,    0.0 ),
    (     -3.0,    0.0,      1.0,    0.0 ),
    (      2.0,    0.0,     -1.0,    0.0 ),
    (     -2.0,    0.0,      1.0,    0.0 ),
    (      2.0,    0.0,     -1.0,    0.0 ),
    (     -2.0,    0.0,      1.0,    0.0 ),
    (      2.0,    0.0,      0.0,    0.0 ),
    (      2.0,    0.0,     -1.0,    0.0 ),
    (      1.0,    0.0,     -1.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,     -1.0,    0.0 ),
    (     -2.0,    0.0,      1.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,     -1.0,    0.0 ),
    (     -1.0,    0.0,      1.0,    0.0 ),
    (     -1.0,    0.0,      1.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,     -1.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 ),
    (     -1.0,    0.0,      0.0,    0.0 ),
    (      1.0,    0.0,      0.0,    0.0 )
)


text:
#!/usr/bin/env julia
#
# @license Apache-2.0
#
# Copyright (c) 2018 The Stdlib Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import JSON

"""
    gen( x, filepath )

Generate fixture data and write to file.

# Arguments

* `x`: domain
* `name::AbstractString`: filepath of the output file

# Examples

``` julia
julia> x = linspace( -0.5, 1.5, 2001 );
julia> gen( x, \"./data.json\" );
```
"""
function gen( x, filepath )
	y = Array( Float64, length( x ) );
	for i in eachindex(x)
		y[i] = erfcinv( x[i] );
	end

	data = Dict([
		("x", x),
		("expected", y)
	]);

	outfile = open( filepath, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# Get the filename:
file = @__FILE__;

# Extract the directory in which this file resides:
dir = dirname( file );

# 0.5 <= x <= 1.5 (linear regime)
x = linspace( 0.5, 1.5, 3000 );
out = joinpath( dir, "x_0.5_1.5.json" );
gen( x, out );

# 0.25 < 2-x < 0.5
x = linspace( 1.5, 1.75, 500 );
out = joinpath( dir, "x_1.5_1.75.json" );
gen( x, out );

# 0.25 < x < 0.5
x = linspace( 0.25, 0.5, 500 );
out = joinpath( dir, "x_0.25_0.5.json" );
gen( x, out );

# 0.25 < 2-x < 0.00012340980408664937
x = linspace( 1.75, 1.9998, 500 );
out = joinpath( dir, "x_1.75_1.9998.json" );
gen( x, out );

# 0.00012340980408664937 < x < 0.25
x = linspace( 0.0002, 0.25, 500 );
out = joinpath( dir, "x_0.0002_0.25.json" );
gen( x, out );

# 0.00012340980408664937 < 2-x < 2.220446049250313e-16
x = linspace( 1.9998, 1.9999999999999998, 500 );
out = joinpath( dir, "x_1.9998_1.9999..8.json" );
gen( x, out );

# 2.220446049250313e-16 < 2-x < 0
x = linspace( 1.9999999999999998, 2, 500 );
out = joinpath( dir, "x_1.9999..8_2.json" );
gen( x, out );

max_stars_repo_name: andbmme/stdlib

max_issues_repo_name: andbmme/stdlib

max_forks_repo_name: andbmme/stdlib


text:
using LinearAlgebra: eigvals,norm2,×,⋅
"""
    ke(I::CartesianIndex,u,U=0)

Compute ½|u-U|² at center of cell `I` where `U` can be used
to subtract a background flow.
"""
ke(I::CartesianIndex{m},u,U=zeros(m)) where m = 0.125sum(@inbounds(abs2(u[I,i]+u[I+δ(i,I),i]-2U[i])) for i ∈ 1:m)
"""
    ∂(i,j,I,u)

Compute ∂uᵢ/∂xⱼ at center of cell `I`. Cross terms are computed
less accurately than inline terms because of the staggered grid.
"""
@fastmath @inline ∂(i,j,I,u) = (i==j ? ∂(i,I,u) :
    @inbounds 0.25*(u[I+δ(j,I),i]+u[I+δ(j,I)+δ(i,I),i]
                   -u[I-δ(j,I),i]-u[I-δ(j,I)+δ(i,I),i]))
"""
    λ₂(I::CartesianIndex{3},u)

λ₂ is a deformation tensor metric to identify vortex cores.
See https://en.wikipedia.org/wiki/Lambda2_method and
Jeong, J., & Hussain, F. doi:10.1017/S0022112095000462
"""
function λ₂(I::CartesianIndex{3},u)
    J = [∂(i,j,I,u) for i ∈ 1:3, j ∈ 1:3]
    S,Ω = (J+J')/2,(J-J')/2
    eigvals(S^2+Ω^2)[2]
end

function permute(f,i)
    j,k = i%3+1,(i+1)%3+1
    f(j,k)-f(k,j)
end
"""
    curl(i,I,u)

Compute component `i` of ∇×u at the __edge__ of cell `I`.
For example `curl(3,CartesianIndex(2,2,2),u)` will compute
`ω₃(x=1.5,y=1.5,z=2)` as this edge produces the highest
accuracy for this mix of cross derivatives on a staggered grid.
"""
curl(i,I,u) = permute((j,k)->∂(j,CI(I,k),u), i)
"""
    ω(I::CartesianIndex{3},u)

Compute 3-vector ω=∇×u at the center of cell `I`.
"""
ω(I::CartesianIndex{3},u) = [permute((j,k)->∂(k,j,I,u), i) for i ∈ 1:3]
"""
    ω_mag(I::CartesianIndex{3},u)

Compute |ω| at the center of cell `I`.
"""
ω_mag(I::CartesianIndex{3},u) = norm2(ω(I,u))
"""
    ω_θ(I::CartesianIndex{3},z,center,u)

Compute ω⋅θ at the center of cell `I` where θ is the azimuth
direction around vector `z` passing through `center`.
"""
function ω_θ(I::CartesianIndex{3},z,center,u)
    θ = z × (collect(I.I)-center)
    n = norm2(θ)
    n<=eps(n) ? 0. : θ ⋅ ω(I,u) / n
end


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

"""
    Beam - Euler-Bernoulli beam for 3d problems

# Load types for distributed load
- P1 and P1 for loads in beam normal 1 and normal 2 direction
- PX, PY, PZ for loads in global coordinate system
"""
struct Beam <: FieldProblem end

FEMBase.get_unknown_field_name(::Problem{Beam}) = "displacement"

function FEMBase.get_integration_points(::Problem{Beam}, ::Element{Seg2})
    return FEMQuad.get_quadrature_points(Val{:GLSEG3})
end

"""
    get_rotation_matrix(X1, X2, n_1)

Given element coordinates ``X_1``, ``X_2``  and the first beam section axis
``n1``, return rotation matrix ``R``, which can be used to define the
orientation of the beam cross-section.

The orientation of a beam cross-section is defined in terms of local,
right-handed ``(t, n_1, n_2)`` axis system, where ``t`` is the tangent to the
axis of the element, positive in the direction from the first to the second node
of the element, and ``n_1`` and ``n_2`` are basis vectors that define the local
1- and 2-directions of the cross section. ``n_1`` is referred to as the first
beam section axis, and ``n_2`` is referred to as the normal to the beam.

Notice that rotation matrix is made of three row vectors which must be linearly
independent. That is, if tangent vector ``t`` is parallel to ``n_1``, rotation
matrix cannot be defined because ``t × n_2 = 0``.

Connection to 2d-beams: for beams in a plane the `n_1`-direction is always
``(0.0, 0.0, -1.0)``, which is normal to the plane in which the motion occurs.
therefore, planar beams can bend only about the first beam section axis.
"""
function get_rotation_matrix(X1, X2, n1)
    t = X2-X1
    L = norm(t)
    a = L/2.0
    Pe = [2.0*a 0.0 0.0; 0.0 0.0 -2.0*a; 0.0 1.0 0.0]
    n2 = cross(t, n1)
    if iszero(n2)
        error("get_rotation_matrix(X1, X2, n1) = get_rotation_matrix($X1, $X2, $n1) ",
              "is failing because the first beam section axis n1 is parallel ",
              "to the tangent vector t. To fix this issue and construct an ",
              "orthonormal set of linearly independent basis vectors, choose ",
              "n1 such that cross-product t × n1 gives a non-zero answer. Now, ",
              "cross($t, $n1) = $n2.")
    end
    Pg = [t n1 n2]
    T = Pe*inv(Pg)
    return T
end

function FEMBase.assemble_elements!(problem::Problem{Beam}, assembly::Assembly,
                                    elements::Vector{Element{Seg2}}, time::Float64)

    B = zeros(4, 12)
    N = zeros(6, 12)
    D = zeros(4, 4)
    Ke = zeros(12, 12)
    Me = zeros(12, 12)
    fe = zeros(12)
    b_loc = zeros(6)
    b_glob = zeros(6)
    Rho = zeros(6,6)

    for element in elements
        fill!(Ke, 0.0)
        fill!(Me, 0.0)
        fill!(fe, 0.0)

        for (w, xi) in get_integration_points(problem, element)

            xi1 = xi[1]
            X1, X2 = element("geometry", time)
            L = norm(X2-X1)

            if haskey(element, "orientation")
                T = element("orientation", xi1, time)
            elseif haskey(element, "normal")
                n1 = element("normal", xi1, time)
                T = get_rotation_matrix(X1, X2, n1)
            end

            Z = zeros(3,3)
            Rd = [T Z Z Z; Z T Z Z; Z Z T Z; Z Z Z T]
            Rf = [T Z; Z T]
            detJ = L/2.0
            s = 2.0/L

            fill!(B, 0.0)
            fill!(b_loc, 0.0)
            fill!(b_glob, 0.0)
            fill!(N, 0.0)
            fill!(D, 0.0)
            fill!(Rho, 0.0)

            N1 = -xi1/2 + 1/2
            N2 =  xi1/2 + 1/2
            dN1 = -1/2*s
            dN2 =  1/2*s

            M1 = xi1^3/4 - 3*xi1/4 + 1/2
            M2 = -xi1^3/4 + 3*xi1/4 + 1/2
            dM1 =  (3*xi1^2/4 - 3/4)*s
            dM2 = (-3*xi1^2/4 + 3/4)*s
            d2M1 = 3*xi1/2*s^2
            d2M2 = -3*xi1/2*s^2

            L1 = (xi1^3/4 - xi1^2/4 - xi1/4 + 1/4)/s
            L2 = (xi1^3/4 + xi1^2/4 - xi1/4 - 1/4)/s
            dL1 = (3*xi1^2/4 - xi1/2 - 1/4)
            dL2 = (3*xi1^2/4 + xi1/2 - 1/4)
            d2L1 = (3*xi1 - 1)/2*s
            d2L2 = (3*xi1 + 1)/2*s

            # Shape functions N

            # for node 1
            N[1,1] =  N1
            N[2,2] =  M1
            N[6,2] =  dM1
            N[3,3] =  M1
            N[5,3] =  dM1
            N[4,4] =  N1
            N[3,5] = -L1
            N[5,5] = -dL1
            N[2,6] =  L1
            N[6,6] =  dL1

            # for node 2
            N[1,7]  =  N2
            N[2,8]  =  M2
            N[6,8]  =  dM2
            N[3,9]  =  M2
            N[5,9]  = -dM2
            N[4,10] =  N2
            N[3,11] = -L2
            N[5,11] = -dL2
            N[2,12] =  L2
            N[6,12] =  dL2

            # Kinematic matrix B

            # for node 1
            B[1,1] = dN1
            B[2,2] = -d2M1
            B[2,6] = -d2L1
            B[3,3] = -d2M1
            B[3,5] = d2L1
            B[4,4] = dN1

            # for node 2
            B[1,7] = dN2
            B[2,8] = -d2M2
            B[2,12] = -d2L2
            B[3,9] = -d2M2
            B[3,11] = d2L2
            B[4,10] = dN2

            # Material matrix D

            E = element("youngs modulus", xi1, time)
            A = element("cross-section area", xi1, time)
            I1 = element("torsional moment of inertia 1", xi1, time)
            I2 = element("torsional moment of inertia 2", xi1, time)
            G = element("shear modulus", xi1, time)
            J = element("polar moment of inertia", xi1, time)
            D[1,1] = E*A
            D[2,2] = E*I1
            D[3,3] = E*I2
            D[4,4] = G*J

            # Assemble stiffness matrix

            Ke += w * Rd'*(B'*D*B)*Rd * detJ

            # Assemble mass matrix (if density is defined)

            if haskey(element, "density")
                rho = element("density", xi1, time)
                Rho[1,1] = rho*A
                Rho[2,2] = rho*A
                Rho[3,3] = rho*A
                Rho[4,4] = rho*J
                #Rho[5,5] = I1 # Looks these have only a very small influence to
                #Rho[6,6] = I2 # result and are neglected in ABAQUS
                Me += w * Rd'*(N'*Rho*N)*Rd * detJ
            end

            # Assemble distributed loads / moments
            if haskey(element, "distributed load 1")
                b_loc[3] = element("distributed load 1", xi1, time)
            end
            if haskey(element, "distributed load 2")
                b_loc[2] = element("distributed load 2", xi1, time)
            end
            if haskey(element, "distributed load x")
                b_glob[1] = element("distributed load x", xi1, time)
            end
            if haskey(element, "distributed load y")
                b_glob[2] = element("distributed load y", xi1, time)
            end
            if haskey(element, "distributed load z")
                b_glob[3] = element("distributed load z", xi1, time)
            end

            fe += w * Rd'*N'*Rf*b_loc * detJ
            fe += w * N'*b_glob * detJ
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.M, gdofs, gdofs, Me)
        add!(assembly.f, gdofs, fe)
    end

    return nothing

end

function FEMBase.assemble_elements!(problem::Problem{Beam}, assembly::Assembly,
                                    elements::Vector{Element{Poi1}}, time::Float64)

    for element in elements
        gdofs = get_gdofs(problem, element)
        for i=1:3
            if haskey(element, "point force $i")
                P = element("point force $i", time)
                add!(assembly.f, gdofs[i], P)
            end
            if haskey(element, "point moment $i")
                P = element("point moment $i", time)
                add!(assembly.f, gdofs[3+i], P)
            end
            if haskey(element, "fixed displacement $i")
                g = element("fixed displacement $i", time)
                add!(assembly.C1, gdofs[i], gdofs[i], 1.0)
                add!(assembly.C2, gdofs[i], gdofs[i], 1.0)
                add!(assembly.g, gdofs[i], g)
            end
            if haskey(element, "fixed rotation $i")
                g = element("fixed rotation $i", time)
                add!(assembly.C1, gdofs[3+i], gdofs[3+i], 1.0)
                add!(assembly.C2, gdofs[3+i], gdofs[3+i], 1.0)
                add!(assembly.g, gdofs[3+i], g)
            end
        end
    end

    return nothing

end


text:
"""
    Δθ = angle_in_plane(vec, mat)
    Δθ = angle_in_plane(vec1, vec2, mat)
    Δθ = angle_in_plane(vec, orbit)
    Δθ = angle_in_plane(vec, vec, orbit)
    Δθ = angle_in_plane(vec, orbit, t)
    Δθ = angle_in_plane(orbit1, orbit2, t)

Computes the angle between two vectors after projecting them to the orbital plane defined by an orbit's
basis matrix `mat`.

If only one vector is provided, the angle between the vector `vec` and the orbit's periapsis
is found.

If a vector, orbit, and time `t` are provided, the in-plane angle between the vector and the orbit's
position at `t` is found.

If two orbits and a time `t` are provide, the in-plane angle between the two orbits' positions at time
`t` is found.
"""
function angle_in_plane(vec::AbstractVector{<:Real}, mat::Union{Rotation{3,Float64},AbstractMatrix}=MRP(1I))
    planevec = mat\vec
    return atan(planevec[2],planevec[1])
end

angle_in_plane(vec::AbstractVector{<:Real}, vecref::AbstractVector{<:Real}, mat::Union{Rotation{3,Float64},AbstractMatrix}=MRP(1I)
    ) = wrap_angle(angle_in_plane(vec,mat) - angle_in_plane(vecref,mat))
angle_in_plane(vec::AbstractVector{<:Real}, vecref::AbstractVector{<:Real}, xdir::AbstractVector{<:Real}, ydir::AbstractVector{<:Real}
    ) = angle_in_plane(vec, vecref, basis_rotation(xdir,ydir))

# angle between the vector `vec` projected to the orbital plane and the periapsis of `orb`
angle_in_plane(vec::AbstractVector{<:Real}, orb::Orbit
    ) = angle_in_plane(vec, orb.basis)

# angle between the vectors `vec` and `vecref` after projecting both to the orbital plane of `orb`
angle_in_plane(vec::AbstractVector{<:Real}, vecref::AbstractVector{<:Real}, orb::Orbit
    ) = angle_in_plane(vec, vecref, orb.basis)

# angle between the vector `vec` projected to the orbital plane and the orbit's position at time `t`
angle_in_plane(vec::AbstractVector{<:Real}, orb::Orbit, t::Real
    ) = wrap_angle(angle_in_plane(vec, time_state_vector(t, orb)[1], orb))

# angle between the orbital positions of `orb` from `orbref` at time `t`, after projecting the position `orb` to the orbital plane of `orbref`
angle_in_plane(orb::Orbit, orbref::Orbit, t::Real
    ) = wrap_angle(angle_in_plane(time_state_vector(t, orb)[1], orbref, t))


text:
using DepositionMTK
using Test
using StaticArrays

# include("../src/DepositionMTK.jl")

# Results from Wesely (1989) table 3; updated to values from Walmsley (1996) table 1.

const SO2 = SA_F32[
	130.0 140 160 380 1000 100 1200
	1400 1400 1400 1400 1500 100 1300
	1100 1100 1100 1100 1200 90 1000
	1000 1000 1000 1000 1100 1100 1100
	270 290 330 620 1100 90 1000]

const O3 = SA_F32[
	100.0 110 130 320 960 960 580
	430 470 520 710 1300 950 580
	390 420 460 610 960 770 510
	560 620 710 1100 3200 3200 3200
	180 200 230 440 950 820 530]

const NO2 = SA_F32[
	120.0 130 160 480 2900 2700 2300
	1900 1900 1900 2000 2700 2500 2200
	1700 1700 1800 1900 2400 2300 2000
	3900 4000 4100 4500 9999 9999 9999
	270 290 350 850 2500 2300 2000]

const H2O2 = SA_F32[
	90.0 90 110 250 640 90 80
	400 430 480 650 1100 90 90
	370 390 430 550 840 90 80
	400 430 470 620 1000 1000 1000
	160 170 200 370 750 90 80]

const ALD = SA_F32[
	330.0 340 370 800 9999 9999 9999
	9999 9999 9999 9999 9999 9999 9999
	9999 9999 9999 9999 9999 9999 9999
	9999 9999 9999 9999 9999 9999 9999
	520 550 630 1700 9999 9999 9999]

const HCHO = SA_F32[
	100.0 110 140 450 6700 1400 1400
	8700 8700 8700 8700 8700 1400 1400
	8300 8300 8300 8300 8400 1400 1400
	2900 2900 2900 2900 2900 2900 2900
	250 270 340 1000 7500 1400 1400]

const OP = SA_F32[
	120.0 130 160 480 2800 2500 2200
	1900 1900 1900 2000 2700 2400 2000
	1700 1700 1800 1800 2400 2100 1900
	3700 3700 3800 4200 8600 8600 8600
	270 290 350 850 2500 2200 1900]

const PAA = SA_F32[
	150.0 160 200 580 2800 2400 2000
	1900 1900 1900 2000 2700 2200 1900
	1700 1700 1700 1800 2400 2000 1800
	3400 3400 3500 3800 7200 7200 7200
	330 350 420 960 2400 2100 1800]

const ORA = SA_F32[
	30.0 30 30 40 50 10 10
	140 140 150 170 190 10 10
	130 140 140 160 180 10 10
	310 340 390 550 910 910 910
	60 60 70 80 90 10 10]

const NH3 = SA_F32[
	80.0 80 100 320 2700 430 430
	3400 3400 3400 3400 3400 440 440
	3000 3000 3000 3000 3100 430 430
	1500 1500 1500 1500 1500 1500 1500
	180 200 240 680 2800 430 430]
const PAN = SA_F32[
	190.0 210 250 700 2900 2700 2300
	1900 1900 1900 2000 2700 2500 2200
	1700 1700 1800 1900 2400 2300 2000
	3900 4000 4100 4500 9999 9999 9999
	410 430 510 1100 2500 2300 2000]
const HNO2 = SA_F32[
	110.0 120 140 330 950 90 90
	1000 1000 1000 1100 1400 90 90
	860 860 870 910 1100 90 90
	820 830 830 870 1000 1000 1000
	220 240 280 530 1000 90 90]

function different(a::Float64, b::Float32) where T<:AbstractFloat
	c = abs(a - b)
	return c/b > 0.1 && c >= 11.0
end

function TestWesely()
	iLandUse = 4                       			# deciduous forest
	Ts = [25.0, 10, 2, 0, 10]        			# Surface Temperature [C]
	Garr = [800.0, 500, 300, 100, 0] 			# Solar radiation [W m-2]
	θ = 0.0                                  	# Slope [radians]

	polNames = [
		"SO2", "O3", "NO2", "H2O2",
		"ALD", "HCHO", "OP", "PAA",
		"ORA", "NH3", "PAN", "HNO2"]

	testData = [SO2, O3, NO2, H2O2, ALD, HCHO, OP, PAA, ORA, NH3, PAN, HNO2]
	gasData = [
		So2Data, O3Data, No2Data,
		H2o2Data, AldData, HchoData,
		OpData, PaaData, OraData,
		Nh3Data, PanData, Hno2Data]

	for i in 1:12
		pol = polNames[i]
		polData = testData[i]
		isSO2, isO3 = false, false
		if pol == "SO2"
			isSO2 = true
		end
		if pol == "O3"
			isO3 = true
		end
		for iSeason in 1:5
			for ig in 1:5
				G = Garr[ig]
				r_c = SurfaceResistance(gasData[i], G, Ts[iSeason], θ,
					iSeason, iLandUse, false, false, isSO2, isO3)
				if different(r_c, polData[iSeason, ig])
					println(pol, iSeason, G, r_c, polData[iSeason, ig])
					return false
				end
			end
			r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ,
				iSeason, iLandUse, false, true, isSO2, isO3) # dew
			if different(r_c, polData[iSeason, 6])
				println(pol, iSeason, "dew", r_c, polData[iSeason, 6])
				return false;
			end
			r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ,
				iSeason, iLandUse, true, false, isSO2, isO3) # rain
			if different(r_c, polData[iSeason, 7])
				println(pol, iSeason, "rain", r_c, polData[iSeason, 7])
				return false;
			end
		end
	end
	return true
end

@testset "DepositionMTK.jl" begin
    @test TestWesely() == true
    @test r_s(1.0, 1.0, 1, 1, true) ≈ 3.0772306344553016e30
    @test r_dc(1.0, 1.0) ≈ 9.181727363545544
    @test r_mx(1.0, 1.0) ≈ 0.009999966666777778
    @test r_smx(1.0, 1.0, 1.0) ≈ 2.0
    @test r_lux(1.0, 1.0, 1, 1, true, false, true, false) ≈ 50
    @test r_clx(1.0, 1.0, 1, 1) ≈ 9.999899563027895e24
    @test r_gsx(1.0, 1.0, 1, 1) ≈ 299.9977500168749
    @test SurfaceResistance(So2Data, 1.0, 1.0, 1.0, 1, 1, true, true, true, false) ≈ 45.45454545454546
end