text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots

f(t, x) = exp(-x^2/(2t))/√(2π*t)
u(t, x) = f(t+0.1, x+2) + 2f(t+0.2, x-2)
t = range(0, 3; length=101)
x = range(-5, 5; length=101)
T = u.(t', x)

surface(t, x, T; colorbar=false,
    xlabel="t", ylabel="x", zlabel="T",
    camera=(30, 70), color=reverse(cgrad(:RdYlBu_11)))

# %%
#using LinearAlgebra, SparseArrays, Plots, DifferentialEquations
#using Surrogates
using Plots
#gr()
L=20
N=100
k=50
x=range(0, stop=100, length=100);
dx=x[2]-x[1];
dt=0.5*dx^2/(2*k);
t=0:dt:10
# print("the vlues of t",t)
#T=zeros(length(t),length(x))
T=zeros(length(x),length(t))
s=size(T)
print("size of matrix",s)
#T[1,2]=T[1,1]=100 #complete 1st row
#T[end,end]=150  #last row of a matrix
@. T[:, 1] = exp(-(x-50)^2/(2*5^2))
for j in 1:1960
    for i=2:99
        T[i,j+1]=T[i,j]+k*dt/(dx^2)*(T[i+1,j]-2*T[i,j] + T[i-1,j])
    end
end
surface(t, x, T; color=:gist_earth)

# %%
#using LinearAlgebra, SparseArrays, Plots, DifferentialEquations
#using Surrogates
using Plots
#gr()
L=20
N=100
k=50
x=range(0, stop=100, length=100);
dx=x[2]-x[1];
dt=0.5*dx^2/(2*k);
t=0:dt:10
# print("the vlues of t",t)
#T=zeros(length(t),length(x))
T=zeros(length(x),length(t)) # t <--> x
s=size(T)
print("size of matrix",s)
#T[1,2]=T[1,1]=100 #complete 1st row
#T[end,end]=150  #last row of a matrix
@. T[:, 1] = exp(-(x-50)^2/(2*5^2)) # initial values
#for i in 2:1960
for j in 1:1960 # i -> j
#    for j=1:99
    for i=2:99 # j -> i, 1:99 -> 2:99
        T[i,j+1]=T[i,j]+k*dt/(dx^2)*(T[i+1,j]-2*T[i,j] + T[i-1,j])
    end
end
#surface(x,t,T)
surface(t,x,T) # x <--> t

# %%
T

# %%
using Plots

k = 50
x = range(0, 100; length=100);
dx = step(x)
dt = 0.5*dx^2/(2k)
t = 0:dt:10
T = zeros(length(x), length(t))
U = zeros(length(x), length(t))

@. T[:, 1] = exp(-(x-50)^2/(2*5^2))

for i in eachindex(x)
    U[i, 1] = exp(-(x[i]-50)^2/(2*5^2))
end

@show T ≈ U

plot(x, T[:, 1]; label="T[:,1]", lw=2)
plot!(x, U[:, 1]; label="U[:,1]", lw=2, ls=:dash)

# %%
@show T[:, 1];

# %%
plot(T[:, 1])

# %%
?@.

# %%
@show VERSION
using Pkg
Pkg.status("Plots")
Pkg.status("PlotUtils"; mode=PKGMODE_MANIFEST)
Pkg.status("GR"; mode=PKGMODE_MANIFEST)

# %%

max_stars_repo_path: 0015/Surrogates/3D plot example - heat eq.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# Bogumił Kamiński, 2019-03-25

using DataFrames, CSV

df = CSV.read("iris.csv")
describe(df)
df2 = stack(df)
CSV.write("iris2.csv", df2)

ext: jl

max_stars_repo_path: 8. IO/ioexample_3.jl

max_issues_repo_path: 8. IO/ioexample_3.jl

max_forks_repo_path: 8. IO/ioexample_3.jl


text:
function getDriverConfigDescriptor()

    return [
            ConfigEntry(
                :casename,
                :required,
                [String,],
            ),

            ConfigEntry(
                :caseroot,
                :required,
                [String,],
            ),

            ConfigEntry(
                :caserun,
                :required,
                [String,],
            ),

            ConfigEntry(
                :archive_root,
                :required,
                [String,],
            ),

            ConfigEntry(
                :archive_list,
                :optional,
                [String,],
                "archive_list.txt",
            ),
   ]
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
abstract type AbstractMonteCarloConfiguration end

mutable struct MonteCarloConfiguration{num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod, baseMode <: BaseMode, rngType <: Random.AbstractRNG} <: AbstractMonteCarloConfiguration
    Nsim::num1
    Nstep::num2
    monteCarloMethod::abstractMonteCarloMethod
    parallelMode::baseMode
    seed::Int64
    offset::Union{Int64, Nothing}
    rng::rngType
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, seed::Number) where {num1 <: Integer, num2 <: Integer}
        return MonteCarloConfiguration(Nsim, Nstep, StandardMC(), SerialMode(), Int64(seed), MersenneTwister())
    end
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::abstractMonteCarloMethod = StandardMC(), seed::Number = 0) where {num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod}
        return MonteCarloConfiguration(Nsim, Nstep, monteCarloMethod, SerialMode(), Int64(seed), MersenneTwister())
    end
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, parallelMethod::baseMode, seed::Number = 0, rng::rngType_ = MersenneTwister()) where {num1 <: Integer, num2 <: Integer, baseMode <: BaseMode, rngType_ <: Random.AbstractRNG}
        return MonteCarloConfiguration(Nsim, Nstep, StandardMC(), parallelMethod, Int64(seed), rng)
    end
    #Most General, no default argument, offset is controllable from outside
    function GeneralMonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::abstractMonteCarloMethod, parallelMethod::baseMode, seed::Number, rng::rngType_) where {num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod, baseMode <: BaseMode, rngType_ <: Random.AbstractRNG}
        if Nsim <= zero(num1)
            error("Number of Simulations must be positive")
        elseif Nstep <= zero(num2)
            error("Number of Steps must be positive")
        else
            return new{num1, num2, abstractMonteCarloMethod, baseMode, rngType_}(Nsim, Nstep, monteCarloMethod, parallelMethod, Int64(seed), 0, rng)
        end
    end
    function GeneralMonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::abstractMonteCarloMethod, parallelMethod::baseMode, seed::Number, rng::MersenneTwister, offset::Number) where {num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod, baseMode <: BaseMode}
        if Nsim <= zero(num1)
            error("Number of Simulations must be positive")
        elseif Nstep <= zero(num2)
            error("Number of Steps must be positive")
        else
            return new{num1, num2, abstractMonteCarloMethod, baseMode, MersenneTwister}(Nsim, Nstep, monteCarloMethod, parallelMethod, Int64(seed), offset, rng)
        end
    end
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::abstractMonteCarloMethod, parallelMethod::baseMode, seed::Number, rng::MersenneTwister, offset::Number) where {num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod, baseMode <: BaseMode}
        if Nsim <= zero(num1)
            error("Number of Simulations must be positive")
        elseif Nstep <= zero(num2)
            error("Number of Steps must be positive")
        else
            return GeneralMonteCarloConfiguration(Nsim, Nstep, monteCarloMethod, parallelMethod, Int64(seed), offset, rng)
        end
    end
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::abstractMonteCarloMethod, parallelMethod::baseMode, seed::Number = 0, rng::rngType_ = MersenneTwister()) where {num1 <: Integer, num2 <: Integer, abstractMonteCarloMethod <: AbstractMonteCarloMethod, baseMode <: BaseMode, rngType_ <: Random.AbstractRNG}
        return GeneralMonteCarloConfiguration(Nsim, Nstep, monteCarloMethod, parallelMethod, Int64(seed), rng)
    end
    function MonteCarloConfiguration(Nsim::num1, Nstep::num2, monteCarloMethod::AntitheticMC, parallelMethod::baseMode, seed::Number = 0, rng::rngType_ = MersenneTwister()) where {num1 <: Integer, num2 <: Integer, baseMode <: BaseMode, rngType_ <: Random.AbstractRNG}
        if div(Nsim, 2) * 2 != Nsim
            error("Antithetic support only even number of simulations")
        else
            return GeneralMonteCarloConfiguration(Nsim, Nstep, monteCarloMethod, parallelMethod, Int64(seed), rng)
        end
    end
end

export MonteCarloConfiguration;


text:
using Base.Iterators
using PredictiveControl
using Test

import Base: iterate

struct TestIter end

mutable struct TestIterState
    iter::UInt
end

function iterate( iter::TestIter )
    state = TestIterState( 0 )
    return (state, state)
end

function iterate( iter::TestIter, state::TestIterState )
    state.iter = state.iter + 1
    return (state, state)
end

struct TestStoppingCondition <: PredictiveControl.FGM.AbstractStopCondition end

function (t::TestStoppingCondition)( state::TestIterState )
    return state.iter > 4
end

iter = TestIter()
iter = PredictiveControl.FGM.halt( iter, TestStoppingCondition() )
iter = take( iter, 10 )
iter = enumerate( iter )

x = nothing

for y in iter
    global x = y
end

(_, finalstate) = x

# It should have ended after 5 iterations now
@test finalstate.iter == 5


text:
module sep_Compute
using MosekTools ; using JuMP
using LinearAlgebra ; const la = LinearAlgebra

srcDir = dirname(dirname(@__FILE__))*"\\src\\"
include(srcDir*"Model\\C_sep_Model.jl")
include(srcDir*"Model\\R_sep_Model.jl")
using .C_sep_Model
using .R_sep_Model

export Computeξₜˢᵉᵖ,
       get_sol_vals,
       get_sol_min_eigval

"""
Inpute: Jump model of separability problem.
Output: MOSEK solution response.
"""
function Computeξₜˢᵉᵖ(model)
    JuMP.set_optimizer(model, Mosek.Optimizer)
    # optimize
    optimize!(model)

    # output results
    println("Primal: ", primal_status(model))
    println("Dual: ", dual_status(model))
    println("Objective: ", objective_value(model))
    return model
end
Computeξₜˢᵉᵖ(ρ,d,t;con_list ="S1G",noBlock = false,isRe = false) = isRe ?
    Computeξₜˢᵉᵖ(R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list)) :
    Computeξₜˢᵉᵖ(C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list,noBlock))

get_sol_vals(arr) = JuMP.value.(arr)
get_sol_min_eigval(arr) = la.eigmin(get_sol_vals(arr))

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
function generateArray(a::Int64, b::Int64, n::Int64=30)::Array{Int64}
    X = [ ]
    for _ in 1:n
        push!(X, Int64(rand(a:b)))
    end
    return X
end

function funk(x)
    mergeSortStep(x,1,length(x))
    return x
end

function mergeSortStep(x,p,r)
    if p < r
        q = Int(floor((p+r)/2))
        mergeSortStep(x,p,q)
        mergeSortStep(x,q+1,r)
        mergeParts(x,p,q,r)
    end
end

function mergeParts(x,p,q,r)
    L = copy(x[p:q])
    R = copy(x[q+1:r])
    push!(L,typemin(Int64))
    push!(R,typemin(Int64))
    i = j = 1
    for k = p:r
        if L[i] >= R[j]
            x[k] = L[i]
            i +=1
        else
            x[k] = R[j]
            j +=1
        end
    end

end

function funk1(x)
    for i in 1:length(x)-1
        minI = i
        for j in i+1:length(x)
            if x[j] < x[minI]
                minI = j
            end
        end
        if i != minI
            x[i], x[minI] = x[minI], x[i]
        end
    end
    return x
end

function funk2(x)
    for i in 2:length(x)
        key = x[i]
        j = i - 1
        while j>0 && x[j] < key
            x[j+1] = x[j]
            j -=1
        end
        x[j+1] = key
    end
    return x
end

function funk3(x)
    quickSortStep(x,1,length(x))
    return x
end

function quickSortStep(x,p,r)
    if p < r
        q = partition(x,p,r)
        quickSortStep(x,p,q-1)
        quickSortStep(x,q+1,r)
    end
end

function partition(x,p,r)
    pivot = x[r]
    i = p - 1
    for j in p:r-1
        if x[j] <= pivot
            i += 1
            x[i], x[j] = x[j], x[i]
        end
    end
    x[i+1], x[r] = x[r], x[i+1]
    return i+1
end

funk3(generateArray(1,10,10))


text:
println("testing geometry")
a = [0.0 0.0;
     0.0 10.0;
     10.0 10.0;
     10.0 0.0;
     0.0 0.0];
@test pt_in_poly([5,5], a)
@test pt_in_poly([0,0], a)
@test !pt_in_poly([-1,0], a)
@test pt_in_poly([1,1], a)
@test !pt_in_poly([11,1], a)
@test pt_in_poly([10,1], a)
@test pt_in_poly([9.99,1], a)

b = [5.0 0.0;
     9.0 10.0;
     3.0 11.0;
     4.0 5.0];
@test pt_in_poly([5,1], b)
@test !pt_in_poly([0,0], b)
@test !pt_in_poly([-1,0], b)
@test !pt_in_poly([1,1], b)
@test !pt_in_poly([11,1], b)
@test !pt_in_poly([10.01,1], b)
@test !pt_in_poly([9.99,1], b)
@test pt_in_poly([4,5], b)
@test pt_in_poly([4,6], b)

@test poly_intersects(a, b)
c = [0.0 0.0;
     3.0 0.0;
     3.0 4.0]
d = [3.0 0.0;
     3.0 4.0;
     6.0 0.0]
@test poly_intersects(c, d)
@test poly_intersects(d, c)
@test !poly_intersects(b, c)
@test !poly_intersects(c, b)
@test poly_intersects(c, a)
@test poly_intersects(b, d)
@test poly_intersects(a, d)

@test pt_on_line_segment([0,0], ([-1,0], [1,0]))
@test pt_on_line_segment([0,0], ([-1,0], [1,0.01]))
@test pt_on_line_segment([0,0], ([-1,0], [1.0001,0]))
@test pt_on_line_segment([0,0], ([-1,0], [1.001,0]))
@test pt_on_line_segment([0,0], ([-1,0], [1.01,0]))
@test pt_on_line_segment([0,0], ([-1,0], [1.1,0]))
@test !pt_on_line_segment([0,0], ([-1,0], [1.1,0.1]))

@test pt_on_line_segment([0,0], ([0,-1], [0,1]))
@test pt_on_line_segment([0,0], ([0,-1], [0.01,1]))
@test pt_on_line_segment([0,0], ([0,-1], [0,1.0001]))
@test pt_on_line_segment([0,0], ([0,-1], [0,1.001]))
@test pt_on_line_segment([0,0], ([0,-1], [0,1.01]))
@test pt_on_line_segment([0,0], ([0,-1], [0,1.1]))
@test !pt_on_line_segment([0,0], ([0,-1], [0.1,1.1]))

@test pt_on_line_segment([1,1], ([0,0], [10,10]))

e = [3.0 3.0;
     6.0 3.0;
     6.0 6.0;
     3.0 6.0;
     3.0 3.0];
@test poly_intersects(a, e)
@test poly_intersects(e, a)

a = [3,4]
b = [4,5]
@test cross2(a, b) == -1
@test cross2(b, a) == 1

a = [0,0,2,2]
b = [0,2,2,0]
@test find_line_segment_intersection(a,b) == [1,1]
@test find_line_intersection(a,b) == [1,1]
b = [0,4,2,2]
@test find_line_segment_intersection(a,b) == [2,2]
@test find_line_intersection(a,b) == [2,2]
b = [0,10,2,8]
@test find_line_segment_intersection(a,b) == nothing
@test find_line_intersection(a,b) == [5,5]
a = [1,0,1,2]
b = [0,2,2,0]
@test find_line_segment_intersection(a,b) == [1,1]
@test find_line_intersection(a,b) == [1,1]
a = [1,0,1,2]
b = [1,2,1,1]
@test find_line_segment_intersection(a,b) == [1,1]
@test find_line_intersection(a,b) == [1,1]

a = [0,0,3,0]
b = [4,5]
@test pt_on_left(b, a) == true
b = [1000,1000]
@test pt_on_left(b, a) == true
b = [-1,-1]
@test pt_on_left(b, a) == false
a = [3,0,0,0]
@test pt_on_left(b, a) == true
b = [4,5]
@test pt_on_left(b, a) == false

subject = [0.0 0.0
           10.0 0.0
           10.0 10.0
           0.0 10.0
           0.0 0.0]
clip = [-1.0 3.0
      11.0 3.0
      11.0 8.0
      -1.0 8.0
      -1.0 3.0]
output = [10.0 3.0
      10.0 8.0
      0.0 8.0
      0.0 3.0
      10.0 3.0]
@test clip_polygon(subject, clip) ≈ output

subject = [0.0 0.0
           10.0 0.0
           10.0 10.0
           0.0 10.0
           0.0 0.0]
clip = [3.0 3.0
      8.0 3.0
      8.0 8.0
      3.0 8.0
      3.0 3.0]
output = [3.0 3.0
      8.0 3.0
      8.0 8.0
      3.0 8.0
      3.0 3.0]
@test clip_polygon(subject, clip) ≈ output

a = [0.0 0.0
    4.0 0.0
    4.0 3.0
    0.0 0.0]
@test poly_area(a) ≈ 6
a = [0.0 0.0
    -4.0 0.0
    -4.0 -3.0
    0.0 0.0]
@test poly_area(a) ≈ 6


text:
print("What is your name? ")
name = chomp(readline())
println("Hello, $name, nice to meet you!")


text:
import AutoML
using Test

@testset "config" begin
    inputs = [("feature1", "String"),
              ("feature2", "Text"),
              ("feature3", "Int"),
              ("feature4", "Image"),
              ("feature5", "Binary Category")]
    outputs = [("label1", "Category"),
               ("label2", "Float")]
    @test typeof(AutoML.Config(inputs, outputs)) == AutoML.Config
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
# This test file test the auxillary functions in FractionalDiffEq.jl

using FractionalDiffEq
using Test

@testset "Test eliminator" begin
    @test isapprox(eliminator(3, 2), [1.0 0.0 0.0; 0.0 0.0 1.0]; atol=1e-5)
end

@testset "Test Riesz Matrix" begin
    @test isapprox(RieszMatrix(0.5, 3, 0.01), [-5.0 4.375 -0.3125; 4.375 -5.0 4.375; -0.3125 4.375 -5.0]; atol=0.00001)
end

@testset "Test ω function" begin
    @test isapprox(omega(3, 0.5), [1.0; -0.5; -0.125; -0.0625]; atol=1e-5)
    @test isapprox(omega(5, 0.5), [1.0; -0.5; -0.125; -0.0625; -0.0390625; -0.02734375]; atol=1e-5)
end

@testset "Test distributed order utils function" begin
    @test isapprox(DOB(x->6*x*(1-x), [0, 1], 0.1, 1, 0.1), [3.547014602981364]; atol=1e-5)
    @test isapprox(DOF(x->6*x*(1-x), [0, 1], 0.1, 1, 0.1), [-3.547014602981364]; atol=1e-5)
    @test isapprox(DORANORT(x->6*x*(1-x), [0, 1], 0.1, 1, 0.1), [2.120468497843435]; atol=1e-5)
end

#FIXME: Add array type tests
@testset "Test Mittag Leffler function" begin
    @test isapprox(mittleff(1, 2, 1), exp(1)-1; atol=1e-5)
    @test isapprox(mittleff(2, 2, 1), sinh(1); atol=1e-5)
    @test isapprox(mittleff(1, 1, 0), 1; atol=1e-5)
    @test isapprox(mittleff(1, 2, [1, 2, 3]), [1.718281828459045; 3.194528049465325; 6.361845641062556]; atol=1e-2)

    @test isapprox(mittleffderiv(1, 1, 1), 2.718281828459045; atol=1e-5)
    @test isapprox(mittleffderiv(1, 1), 2.718281828459045; atol=1e-5)

    @test isapprox(mittleffderiv(1, 1, 1/3), 1.39561242508609; atol=1e-5)
    @test isapprox(mittleffderiv(2, 2, 1), 0.18393972058572117; atol=1e-5)

    @test isapprox(mittleff(1, 2, 1, 1), 1.7182818284590438; atol=1e-5)
    @test isapprox(mittleff(2, 2, 1, 1), 1.1752011936438016; atol=1e-3)
    @test isapprox(mittleff(2.1, 2, 1, 1), 1.1527981788744044; atol=1e-4)

    # Also tests for Complex z:

    m = mittleff(1, 2, im)
    @test isapprox(real(m), 0.8414709848078965; atol=1e-5)
    @test isapprox(imag(m), 0.4596976941318611; atol=1e-5)
end

@testset "Test meshgrid" begin
    a, b = meshgrid(0:1, 0:1)
    isapprox(a, [0 1; 0 1]; atol=1e-2)
    isapprox(b, [0 0; 1 1]; atol=1e-2)
end

@testset "Test isFunction" begin
    @test isFunction(x->x)==true
    @test isFunction("Hello")==false
end

@testset "Test singleterm show method" begin
    # SingleTermFODEProblem
    singlefun(x, y) = 1-y
    singletermprob = SingleTermFODEProblem(singlefun, 1.8, 0, 5)
    @test_nowarn show(singletermprob)
end

@testset "Test Multiterms show method" begin
    # MultiTermsFODEProblem
    rightfun(x, y) = 172/125*cos(4/5*x)
    multitermsprob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0 0 0 0 0 0], (0, 10))
    multitermssol = solve(multitermsprob, 0.5, PIEX())

    @test_nowarn show(multitermsprob)
    @test_nowarn show(multitermssol)
end

@testset "Test FODESystem show method" begin
    # FODESystem
    h=0.5; tspan=(0, 1)
    alpha = [0.99, 0.99, 0.99]
    x0 = [1, 0, 1]
    function testf!(du, u, p, t)
        a, b, c = 10, 28, 8/3
        du[1] = a*(u[2]-u[1])
        du[2] = u[1]*(b-u[3])-u[2]
        du[3] = u[1]*u[2]-c*u[3]
    end
    fodesystemprob = FODESystem(testf!, alpha, x0, tspan)
    fodesystemsol = solve(fodesystemprob, h, GL())

    @test_nowarn show(fodesystemprob)
    @test_nowarn show(fodesystemsol)
end

@testset "Test FDDEProblem show method" begin
    # FDDEProblem
    function ϕ(x)
        if x == 0
            return 19.00001
        else
            return 19.0
        end
    end
    function delayf(t, y, ϕ)
        return 3.5*y*(1-ϕ/19)
    end
    fddeprob = FDDEProblem(delayf, ϕ, 0.97, 0.8, (0, 2))
    fddesol = solve(fddeprob, 0.5, DelayPI())

    @test_nowarn show(fddeprob)
    @test_nowarn show(fddesol)
end

@testset "Test FDDESystem show method" begin
    function EnzymeKinetics!(dy, y, ϕ, t)
        dy[1] = 10.5-y[1]/(1+0.0005*ϕ[4]^3)
        dy[2] = y[1]/(1+0.0005*ϕ[4]^3)-y[2]
        dy[3] = y[2]-y[3]
        dy[4] = y[3]-0.5*y[4]
    end
    q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
    fddesystemprob = FDDESystem(EnzymeKinetics!, q, α, 4, 6)

    @test_nowarn show(fddesystemprob)
end

@testset "Test DODEProblem show method" begin
    h = 0.5; t = collect(0:h:1)
    dodefun(t)=0
    dodeprob = DODEProblem([1, 0.1], [x->6*x*(1-x), 0], (0, 1), dodefun, 1, t)
    dodesol = solve(dodeprob, h, DOMatrixDiscrete())

    @test_nowarn show(dodeprob)
    @test_nowarn show(dodesol)
end

@testset "Test FractionalDifferenceProblem" begin
    differencefun(x) = 0.5*x+1
    α=0.5;x0=1;
    T=1; h=0.1
    differenceprob = FractionalDifferenceProblem(differencefun, α, x0)
    differencesol = solve(differenceprob, T, h, PECEDifference())

    @test_nowarn show(differenceprob)
    @test_nowarn show(differencesol)
end

@testset "Test FIEProblem show method" begin
    fietspan = LinRange(-1, 1, 5)
    fieprob = FIEProblem([1, 1], [1, 0.5], 1, fietspan)
    fiesol = solve(fieprob, 20, SpectralUltraspherical())

    @test_nowarn show(fieprob)
    @test_nowarn show(fiesol)
end


text:
include("kern_proc.jl")

d = 10
n = [50, 100, 500]

λ = rand(d)

cov_procs = Array(Proc, 0)
push!(cov_procs, KernelTest(SEArd(λ, 1.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(Mat12Ard(λ, 1.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(Mat32Ard(λ, 1.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(Mat52Ard(λ, 1.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(RQArd(λ, 1.0, 1.0), d, GaussianProcesses.cov))
cov_table = run(cov_procs, n)

grad_stack_procs = Array(Proc, 0)
push!(grad_stack_procs, KernelTest(SEArd(λ, 1.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(Mat12Ard(λ, 1.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(Mat32Ard(λ, 1.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(Mat52Ard(λ, 1.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(RQArd(λ, 1.0, 1.0), d, GaussianProcesses.grad_stack))
grad_stack_table = run(grad_stack_procs, n)

println("\ncov calculations")
show(cov_table)

println("\grad_stack calculations")
show(grad_stack_table)


text:
BENCHMARK_ENV = abspath((@__DIR__) * "/..")
BENCHMARK_RESULTS_DIR = (@__DIR__) * "/results"
PACKAGE_DIR = BENCHMARK_ENV

using Pkg
Pkg.activate(BENCHMARK_ENV)

using DimensionalData, BenchmarkTools, Unitful, SparseArrays, Dates, Statistics, LibGit2

if !isdir(BENCHMARK_RESULTS_DIR)
    mkdir(BENCHMARK_RESULTS_DIR)
end

# Utils

"""Generate a path for benchmarking results using current repository status"""
function gen_results_pth()
    repo = LibGit2.GitRepo(PACKAGE_DIR)
    head = LibGit2.head(repo)
    branchname = replace(LibGit2.shortname(head), "/" => "\\/")  # Normalize paths
    hash_str = string(LibGit2.GitHash(LibGit2.peel(LibGit2.GitCommit, head)))
    basename = "$(BENCHMARK_RESULTS_DIR)/$(branchname)-$(hash_str)"
    if LibGit2.isdirty(repo)
        basename *= "-dirty"
    end
    if isfile(basename * ".json")
        i = 1
        while isfile(basename * "-$i.json")
            i += 1
        end
        basename *= "-$i"
    end
    basename *= ".json"
end

#= Benchmarks

Test how much the recalculation of coordinates and dim types
costs over standard getindex/view.

Indexing with Y(1) has no overhead at all, but ranges
have an overhead for slicing the dimensions. =#

const suite = BenchmarkGroup()

suite["view"] = BenchmarkGroup(["view"])
g = DimArray(rand(100, 50), (X(51:150), Y(-40:9)))

vi1(g) = view(data(g), 1, 2)
vd1(g) = view(g, X(1), Y(2))
vi2(g) = view(data(g), :, :)
vd2(g) = view(g, X(:), Y(:))
vi3(g) = view(data(g), 10:40, 1:20)
vd3(g) = view(g, X(10:40), Y(1:20))

suite["view"]["Parent indices with Number"] = @benchmarkable vi1($g)
suite["view"]["Dims with Number"] = @benchmarkable vd1($g)
suite["view"]["Parent indices with Colon"] = @benchmarkable vi2($g);
suite["view"]["Dims with Colon"] = @benchmarkable vd2($g);
suite["view"]["Parent indices with UnitRange"] = @benchmarkable vi3($g);
suite["view"]["Dims with UnitRange"] = @benchmarkable vd3($g);

suite["getindex"] = BenchmarkGroup()
i1(g) = data(g)[10, 20]
d1(g) = g[Y(10), X(20)]
i2(g) = data(g)[:, :]
d2(g) = g[Y(:), X(:)]
i3(g) = data(g)[1:20, 10:40]
d3(g) = g[Y(1:20), X(10:40)]

suite["getindex"]["Parent indices with Number"] = @benchmarkable i1($g)
suite["getindex"]["Dims with Number"] = @benchmarkable d1($g)
suite["getindex"]["Parent indices with Colon"] = @benchmarkable i2($g)
suite["getindex"]["Dims with Colon"] = @benchmarkable d2($g)
suite["getindex"]["Parent indices with UnitRange"] = @benchmarkable i3($g)
suite["getindex"]["Dims with UnitRange"] = @benchmarkable d3($g);

a = rand(5, 4, 3);
da = DimArray(a, (Y((1u"m", 5u"m")), X(1:4), Ti(1:3)))
dimz = dims(da)

if VERSION > v"1.1-"
    suite["eachslice"] = BenchmarkGroup()
    suite["eachslice"]["array_intdim"] = @benchmarkable (()->eachslice($a; dims = 2))();
    suite["eachslice"]["dimarray_intdim"] = @benchmarkable (()->eachslice($da; dims = 2))();
    suite["eachslice"]["dimarray_dim"] = @benchmarkable (()->eachslice($da; dims = Y()))();
    suite["eachslice_to_vector"] = BenchmarkGroup()
    suite["eachslice_to_vector"]["array_intdim"] = @benchmarkable [slice for slice in eachslice($a; dims = 2)];
    suite["eachslice_to_vector"]["dimarray_intdim"] = @benchmarkable [slice for slice in eachslice($da; dims = 2)];
    suite["eachslice_to_vector"]["dimarray_dim"] = @benchmarkable [slice for slice in eachslice($da; dims = X())];
    # @test [slice for slice in eachslice(da; dims=1)] == [slice for slice in eachslice(da; dims=Y)]
end

suite["mean"] = BenchmarkGroup()
suite["mean"]["array_intdim"] = @benchmarkable mean($a; dims = 2);
suite["mean"]["dimarray_intdim"] = @benchmarkable mean($da; dims = 2);
suite["mean"]["dimarray_dim"] = @benchmarkable mean($da; dims = X());
suite["permutedims"] = BenchmarkGroup()
suite["permutedims"]["array_intdim"] = @benchmarkable permutedims($a, (2, 1, 3))
suite["permutedims"]["dimarray_intdim"] = @benchmarkable permutedims($da, (2, 1, 3))
suite["permutedims"]["dimarray_dim"] = @benchmarkable permutedims($da, (Y(), X(), Ti()))
suite["reverse"] = BenchmarkGroup()
suite["reverse"]["array_intdim"] = @benchmarkable reverse($a; dims = 1)
suite["reverse"]["dimarray_intdim"] = @benchmarkable reverse($da; dims = 1)
suite["reverse"]["dimarray_dim"] = @benchmarkable reverse($da; dims = Y())

# Sparse (and similar specialised arrays)

@dim Var "Variable"
@dim Obs "Observation"

sparse_a = sprand(1000, 1000, 0.1)
sparse_d = DimArray(sparse_a, (Var <| 1:1000, Obs <| 1:1000))

suite["sparse"] = BenchmarkGroup()
# Benchmarks
suite["sparse"]["mean with dims arge: regular sparse"] = @benchmarkable mean($sparse_a, dims = $1)
suite["sparse"]["mean with dims arge: dims sparse"] = @benchmarkable mean($sparse_d, dims = $(Var()))

suite["sparse"]["mean: regular sparse"] = @benchmarkable mean($sparse_a)
suite["sparse"]["mean: dims sparse"] = @benchmarkable mean($sparse_d)

suite["sparse"]["copy: regular sparse"] = @benchmarkable copy($sparse_a)
suite["sparse"]["copy: dims sparse"] = @benchmarkable copy($sparse_d)

suite["sparse"]["reduce: regular sparse"] = @benchmarkable reduce(+, $sparse_a)
suite["sparse"]["reduce: dims sparse"] = @benchmarkable reduce(+, $sparse_a)

suite["sparse"]["map: regular sparse"] = @benchmarkable map(sin, $sparse_a)
suite["sparse"]["map: dims sparse"] = @benchmarkable map(sin, $sparse_a)

paramspath = joinpath(dirname(@__FILE__), "params.json")
resultpath = gen_results_pth()

if isfile(paramspath)
    println("Found tuning info at $(paramspath)")
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    println("Running tuning...")
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
    println("Tuning results cached to $(paramspath)")
end

results = run(suite, verbose = true)
println("Writing results to: $(resultpath)")
BenchmarkTools.save(resultpath, results)


text:
@testset "intrinsics" begin

@testset "math" begin
    a = ones(Float32,1)
    a .* Float32(3.14)
    bufferA = MtlArray(a)
    vecA = unsafe_wrap(Vector{Float32}, bufferA.buffer, 1)

    function intr_test(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] = cos(buf[idx])
        return nothing
    end
    @metal intr_test(bufferA.buffer)
    synchronize()
    @test vecA ≈ cos.(a)

    function intr_test2(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] = Metal.rsqrt(buf[idx])
        return nothing
    end
    @metal intr_test2(bufferA.buffer)
end

@testset "sync" begin
    function sync_test_kernel(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] += UInt8(1)
        return nothing
    end
    buf = MtlArray{UInt8,1}(undef, tuple(1024); storage=Shared)
    vec = unsafe_wrap(Vector{UInt8}, buf.buffer, (1024))
    @metal threads=1024 sync_test_kernel(buf)
    synchronize()
    @test all(vec .== UInt8(1))

    function barrier_test_kernel(buf)
        idx = thread_position_in_grid_1d()
        if thread_position_in_threadgroup_1d() != UInt32(1)
            for i in range(1,threads_per_threadgroup_1d())
                buf[idx] += UInt32(i)
            end
            buf[idx] = 1
        end

        threadgroup_barrier()

        if thread_position_in_threadgroup_1d() == UInt32(1)
            for i in range(1,threads_per_threadgroup_1d())
                buf[idx] += buf[idx+i-1]
            end
        end
        return nothing
    end

    buf = MtlArray{Int,1}(undef, tuple(1024); storage=Shared)
    vec = unsafe_wrap(Vector{Int}, buf.buffer, (1024))
    @metal threads=1024 barrier_test_kernel(buf)
    synchronize()
    @test vec[1] == 992

    # TODO: simdgroup barrier test
end

end


text:
include("../roblib.jl")

function mark_chess_and_return(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    direction       = Right
    put_marker_flag = true

    if (sum(num_steps[1:2:end]) % 2 != 0)
        put_marker_flag = false
    end

    while (isborder(r, Up) == false) || (isborder(r, Right) == false)
        if put_marker_flag
            putmarker!(r)
            put_marker_flag = false
        else
            put_marker_flag = true
        end
        if isborder(r, direction)
            direction = inverse(direction)
            move!(r, Up)
        else move!(r, direction)
        end
    end

    through_rectangles_into_angle(r, (Down, Left))
    movements!(r, (Up, Right), reverse(num_steps))
end


text:
@testset "solve" begin
    @testset "Invalid input" begin
        @polyvar x y z
        @test_throws ArgumentError solve([x-2y+2, 0])
        @test_throws ArgumentError solve([x-2, y-2], [x-2, y-2,y+2], [[2, -3]])
        @test_throws ArgumentError solve([x-2z, y^2+3z, z^3+x^3], homvar=z)

        # constant term
        @test_throws ArgumentError solve([subs(x+2, x=>2), y^2+3x])
        @polyvar u v
        @test_throws ArgumentError solve([subs(x+2, x=>2), y^2+3x] ∘ [u, v])
        # overdetermined and abstract system
        @test_throws ArgumentError solve(FPSystem([x-2z, y^2+3z^2, z^3+x^3, z+x]))
        @test_throws ArgumentError solve(FPSystem([x-2z, y^2+3z^2, z^3+x^3, z+x]), homvar=4)

        @test_throws ArgumentError solve([x-2z, y^2+3z^2, z^3+x^3])

        @test_throws ArgumentError solve([x-2z, y^2+3z], homvar=z)
        @test_throws ArgumentError solve([x - 1], patch=OrthogonalPatch())
        # invalid kwargs
        @test_throws ArgumentError solve(equations(cyclic(5)), def=0.4, abc=23)

        # test numerical homogeneous check fails
        @polyvar x y z
        G = FPSystem([x-2z, y^2+3z])
        @test_throws ArgumentError solve(G, homvar=3)
    end

    @testset "solve" begin
        @polyvar x
        @test nfinite(solve([x - 1], threading=false, save_all_paths=true)) == 1
        F = equations(katsura(5))
        @test nfinite(solve(F, threading=false)) == 32
        @test nfinite(solve(F, system=SPSystem, threading=false)) == 32
        @test nfinite(solve(F, system=FPSystem, threading=false)) == 32

        # Simple step size
        F = equations(katsura(5))
        @test nfinite(solve(F, simple_step_size_alg=true, threading=false)) == 32

        # scaling
        F = equations(katsura(5))
        @test nfinite(solve(F, system_scaling=nothing, threading=false)) == 32
        @test nfinite(solve(F, system_scaling=:equations, threading=false)) == 32
        @test nfinite(solve(F, system_scaling=:equations_and_variables, threading=false)) == 32
        @test_throws ArgumentError solve(F, system_scaling=:lalala)

        @test nfinite(solve(F, homotopy=StraightLineHomotopy)) == 32
        result = solve(F, predictor=Euler(), homotopy=StraightLineHomotopy)
        @test nresults(result) == 32
        @test nfinite(solve(F, accuracy=1e-5)) == 32

        result = solve(F)
        @test nfinite(result) == 32
        @test string(result) isa String

        F_hom = homogenize(F)
        @test nfinite(solve(F_hom, threading=false, patch=RandomPatch())) == 32
        @test nfinite(solve(F_hom, threading=false, patch=EmbeddingPatch())) == 32
        @test nfinite(solve(F_hom, threading=false, patch=OrthogonalPatch())) == 32

        @polyvar w
        F = equations(cyclic(5))
        result = solve(homogenize(F, w), threading=false, homvar=w)
        @test nfinite(result) == 70
        @test result isa Result

        @polyvar x y z
        G = [x-2, y+3]
        F = [x+2, y-3]
        @test nfinite(solve(G, F, [[2, -3]])) == 1
        @test nfinite(solve(G, F, [[2+0.0im, -3.0+0im]])) == 1
        @test nfinite(solve(G, F, [@SVector [2, -3]])) == 1

        F = FPSystem(homogenize(equations(cyclic(5))))
        result = solve(F, homvar=6, save_all_paths=true)
        @test nfinite(result) == 70

        F = equations(katsura(5))
        result = solve(homogenize(F), threading=false)
        @test result isa Result{<:ProjectiveVectors.PVector}
        @test is_projective(first(result))
        @test nnonsingular(result) == 32

        prob, startsolutions = problem_startsolutions(TargetSystemInput(F))
        result = solve(prob.homotopy, map(s -> HC.embed(prob, s), startsolutions))
        @test result isa Result{<:Vector}
        @test nnonsingular(result) == 32

        prob, startsolutions = problem_startsolutions(TargetSystemInput(F); affine_tracking=false)
        proj_starts = HC.embed.(prob, startsolutions)
        result = solve(prob.homotopy, proj_starts, homvar=homvars(prob)[1])
        @test result isa Result{<:Vector}
        @test nnonsingular(result) == 32
        @test nfinite(result) == 32

        prob, startsolutions = problem_startsolutions(TargetSystemInput(F); projective_tracking=true)
        proj_starts = HC.embed.(prob, startsolutions)
        result = solve(prob.homotopy, proj_starts)
        @test result isa Result{<:ProjectiveVectors.PVector}

        # Composition
        @polyvar a b c x y z u v
        e = [u + 1, v - 2]
        f = [a * b - 2, a*c- 1]
        g = [x+y, y + 3, x + 2]
        res = solve(e ∘ f ∘ g)
        @test nnonsingular(res) == 2
        @test nnonsingular(solve(e ∘ f ∘ g, system_scaling=nothing)) == 2
        @test nnonsingular(solve(e ∘ f ∘ g, system_scaling=:equations_and_variables)) == 2

        res = solve(e ∘ f ∘ g, system=SPSystem)
        @test nnonsingular(res) == 2

        # Affine SPSystem input
        @test nsolutions(solve(SPSystem(equations(katsura(5))); threading=false)) == 32
    end

    @testset "Random seed" begin
        R = solve(equations(katsura(5)), seed=1234)
        @test seed(R) == 1234

        @polyvar x y

        G = [x-2, y+3]
        F = [x+2, y-3]
        @test seed(solve(G, F, [[2, -3]], seed=222)) == 222
    end

    @testset "Singular solutions" begin
        @polyvar x y z
        z = 1
        F = [x^2 + 2*y^2 + 2*im*y*z, (18 + 3*im)*x*y + 7*im*y^2 - (3 - 18*im)*x*z - 14*y*z - 7*im*z^2]
        result = solve(F)
        @test multiplicity(first(nonsingular(result))) == 1
        @test multiplicity(first(singular(result))) == 3
        @test nsingular(result) == 1
        @test nsingular(result; counting_multiplicities=true) == 3
        @test length(singular(result)) == 1
        @test all(r -> multiplicity(r) == 3, singular(result; multiple_results = true))
    end

    @testset "Path jumping" begin
        result = solve(equations(katsura(5)); accuracy=5e-4, refinement_accuracy=1e-8,
            seed=39813, threading=false, affine_tracking=true)
        @test nreal(result) == 16
    end

    @testset "Affine vs projective" begin
        @polyvar x y z
        F = solve([x-2y, y-2z])
        @test length(F[1].solution) == 3
        @test F[1].solution isa ProjectiveVectors.PVector

        G = solve([x-2, y-2])
        @test length(G[1].solution) == 2
        @test G[1].solution isa Vector
    end

    @testset "Parameter Homotopies" begin
        @polyvar x a y b
        F = [x^2-a, x*y-a+b]
        p = [a, b]
        S = solve(F, [[1.0, 1.0 + 0.0*im]], parameters=p, p₁=[1, 0], p₀=[2, 4])

        @test S[1].solution ≈ [complex(√2), -complex(√2)]
        @test nfinite(S) == 1

        @polyvar x a y b z
        F = [x^2-a*z^2, x*y-(a-b)*z^2]
        S = solve(F, [[1.0, 1.0 + 0.0*im, 1.0]];
                parameters=[a, b],
                start_parameters=[1, 0], target_parameters=[2, 4])
        @test S isa Result{<:ProjectiveVectors.PVector}
        @test ProjectiveVectors.affine_chart(solution(S[1])) ≈ [√2, -√2]
        @test nnonsingular(S) == 1

        S2 = solve(F, [[1.0, 1.0 + 0.0*im, 1.0]];
                parameters=[a, b], p₁=[1, 0], p₀=[2, 4],
                homvar=z, threading=false)
        @test S2 isa Result{<:Vector}
        @test solution(S2[1]) ≈ [√2, -√2]
        @test nfinite(S2) == 1

        γ₁, γ₀ =randn(ComplexF64, 2)
        S2 = solve(F, [[1.0, 1.0 + 0.0*im, 1.0]];
                parameters=[a, b], p₁=[1, 0], p₀=[2, 4],
                γ₁=γ₁, γ₀=γ₀, homvar=z)
        @test solution(S2[1]) ≈ [√2, -√2]
        @test nfinite(S2) == 1

        γ₁, γ₀ = randn(ComplexF64, 2)
        # test SVector input
        S2 = solve(F, [@SVector [1.0, 1.0 + 0.0*im, 1.0]];
                    parameters=[a, b], p₁=[1, 0], p₀=[2, 4],
                    γ₁=γ₁, γ₀=γ₀, homvar=z)
        @test solution(S2[1]) ≈ [√2, -√2]
        @test nfinite(S2) == 1
    end

    @testset "Keep invalid start values" begin
        @polyvar x[1:2] a[1:2]
        F = [x[1]^2-a[1], x[1]*x[2]-a[1]+a[2]]
        startsolutions = [[100, 100]]
        p₁ = [1, 1]
        p₀ = [3im, 0.5+2im]
        res = solve(F, startsolutions; parameters=a, start_parameters=p₁, target_parameters=p₀)
        @test nfailed(res) == 1
        @test res[1].return_code == :terminated_invalid_startvalue
    end

    @testset "ParameterHomotopy with Composition" begin
        @polyvar p q a b c x y z u v
        f = [a * b - 2, a*c- 1]
        g = [x + y, y + 3, x + 2]
        res = solve(f ∘ g; system=SPSystem, threading=false)
        @test nnonsingular(res) == 2
        # parameters at the end
        f2 = [a * b - q, a*c- p]
        g = [x + y, y + 3, x + 2]
        r = solve(f2 ∘ g, solutions(res);
                parameters=[p, q], p₁=[1, 2], p₀=[2, 3],
                threading=false)
        @test nnonsingular(r) == 2

        # parameters at the beginning
        f = [a * b - 2, a*c- 1]
        g2 = [x + y, y + u, x + v]
        r = solve(f ∘ g2, solutions(res);
                    parameters=[u, v], p₁=[3, 2],
                    p₀=[-2, 3], threading=false)
        @test nnonsingular(r) == 2

        # parameter in the middle
        e = [u + 1, v - 2]
        res2 = solve(e ∘ f ∘ g; system=SPSystem)
        f2 = [a * b - q, a * c - p]
        r = solve(e ∘ f2 ∘ g, solutions(res2);
                    parameters=[p, q], p₁=[1, 2],
                    p₀=[2, 3], threading=false)
        @test nnonsingular(r) == 2
    end

    @testset "Coefficient Homotopy" begin
        @polyvar x a y b
        E = [[2 1 0; 0 0 0], [1 0; 1 0]]
        start = [[1.0+0im, -3.0, 2.0], [2.0+0im, -2.0]]
        target = [randn(ComplexF64, 3), randn(ComplexF64, 2)]
        H = CoefficientHomotopy(E, start, target)
        result = solve(H, [[1, 1]])
        @test is_success(result[1])
    end

    @testset "Overdetermined" begin
        Random.seed!(1234567)
        @polyvar x y z w
        a = [0.713865+0.345317im, 0.705182+0.455495im, 0.9815+0.922608im, 0.337617+0.508932im]

        f = [x*z-y^2, y*w-z^2, x*w-y*z]
        L₁ = [1, -1, 1, -1] ⋅ [x, y, z, w]
        L₂ = rand(ComplexF64, 4) ⋅ [x, y, z, w]
        S = solve([f; L₁], [f; L₂], [[1, 1, 1, 1]])
        @test nnonsingular(S) == 1

        f = [x*z-y^2, y-z^2, x-y*z]
        L₁ = [1, -1, 1, -1] ⋅ [x, y, z, 1]
        L₂ = rand(ComplexF64, 4) ⋅ [x, y, z, 1]
        S = solve([f; L₁], [f; L₂], [[1, 1, 1]])
        @test nfinite(S) == 1

        @polyvar x y z
        p₁ = (y - x^2) * (x^2+y^2+z^2 - 1) * (x - 0.5)
        p₂ = (z - x^3) * (x^2+y^2+z^2 - 1) * (y - 0.5)
        p₃ = (y - x^2) * (z - x^3) * (x^2+y^2+z^2 - 1) * (z - 0.5)
        L₁ = [rand(ComplexF64, 4) ⋅ [x, y, z, 1], rand(ComplexF64, 4) ⋅ [x, y, z, 1]]
        L₂ = [rand(ComplexF64, 4) ⋅ [x, y, z, 1], rand(ComplexF64, 4) ⋅ [x, y, z, 1]]

        s = first(solutions(solve([x^2+y^2+z^2 - 1; L₁])))
        tracker, starts = pathtracker_startsolutions(
                [[p₁, p₂, p₃]; L₁], [[p₁, p₂, p₃]; L₂], s;
                affine_tracking=false)
        @test is_success(track(tracker, s))

        s = first(solutions(solve([x^2+y^2+z^2 - 1; L₁])))
        tracker, starts = pathtracker_startsolutions(
                [[p₁, p₂, p₃]; L₁], [[p₁, p₂, p₃]; L₂], s;
                projective_tracking=true)
        @test is_success(track(tracker, s))
    end

    @testset "MultiHomogeneous" begin
        @polyvar x y u v

        f = [x*y - 6, x^2 - 5]
        @test bezout_number(f, variable_groups=[(x,), (y,)]) == 2
        @test bezout_number(f, variable_groups=((x,), (y,))) == 2
        @test bezout_number(f, variable_groups=[[x], [y]]) == 2
        S = solve(f, variable_groups=[(x,), (y,)])
        @test nnonsingular(S) == 2
        S = solve(f, variable_groups=[[x], [y]])
        @test nnonsingular(S) == 2
        g = [x*y - 6u*v, x^2 - u^2]
        S = solve(g, variable_groups=[(x,u), (y,v)], homvars=(u,v))
        @test nnonsingular(S) == 2

        # The 6-R inverse problem

        ## initialize the variables
        @polyvar z[1:6,1:3]
        F = let
            p = [1, 1, 0]
            α = randn(5)
            a = randn(9)
            ## define the system of polynomials
            f = [z[i,:] ⋅ z[i,:] for i=2:5]
            g = [z[i,:] ⋅ z[i+1,:] for i=1:5]
            h = sum(a[i] .* (z[i,:] × z[i+1,:]) for i=1:3) +
                sum(a[i+4] .* z[i,:] for i=2:5)
            F′ = [f .- 1; g .- cos.(α); h .- p]
            ## assign values to z₁ and z₆
            [subs(f, z[1,:] => [1, 0, 0], z[6,:] => [1,0,0]) for f in F′]
        end
        group1 = [[z[2,:]; z[4,:]], [z[3,:]; z[5,:]]]
        group2 = [z[2,:], z[4,:], z[3,:], z[5,:]]

        @test bezout_number(F) == 1024
        @test bezout_number(F; variable_groups=group1) == 320
        @test bezout_number(F; variable_groups=group2) == 576

        @test nnonsingular(solve(F; seed=1234, show_progress=false)) == 16
        @test nnonsingular(solve(F; variable_groups=group1, seed=1234, show_progress=false)) == 16
    end

    @testset "Singular1" begin
        @polyvar x z y
        # This has two roots of multiplicity 6 at the hyperplane z=0
        F = [0.75x^4 + 1.5x^2*y^2-2.5x^2*z^2+0.75*y^4 - 2.5y^2*z^2 + 0.75z^4; 10x^2*z + 10*y^2*z-6z^3]
        res = solve(F; threading=false)
        @test nsingular(res) == 2
        @test nsingular(res; counting_multiplicities=true) == 12
        @test multiplicity.(results(res)) == [6,6]

        F_affine = subs.(F, Ref(y => 1))
        res_affine = solve(F_affine; threading=false)
        @test nsingular(res_affine) == 2
        @test nsingular(res_affine; counting_multiplicities=true) == 12
        @test multiplicity.(results(res_affine)) == [6,6]
        @test_nowarn sprint(show, res_affine)

        res_affine = solve(F_affine; threading=true)
        @test nsingular(res_affine) == 2
        @test nsingular(res_affine; counting_multiplicities=true) == 12
        @test multiplicity.(results(res_affine)) == [6,6]
        @test_nowarn sprint(show, res_affine)

        res_affine2 = solve(F_affine; affine_tracking=false, threading=false)
        @test nsingular(res_affine2) == 2
        @test nsingular(res_affine2; counting_multiplicities=true) == 12
        @test multiplicity.(results(res_affine2)) == [6,6]
        @test_nowarn sprint(show, res_affine2)

        # change multiplicity tolerance to only have one
        multiplicities!(res_affine2; tol=10)
        @test nsingular(res_affine2) == 1
        @test multiplicity.(results(res_affine2)) == [12]
    end

    @testset "Parameter promotion to Float64" begin
        d = 3
        k = 1

        R = 1.0
        r = 0.5

        @polyvar lam[1:k] p[1:6]
        F = [(R^2 - r^2 + (p[1] + p[4] * lam[1])^2 + (p[2] + p[5] * lam[1])^2 + (p[3] + p[6] * lam[1])^2)^2 - 4 * R^2 * ((p[1] + p[4] * lam[1])^2 + (p[2] + p[5] * lam[1])^2)]
        #Randomly choose a start system
        # with only Complex{Float32 parameters
        p0 = randn(Complex{Float32}, 6)
        F0 = subs(F, p => p0)
        res0 = solve(F0)
        S0 = solutions(res0)

        #Construct the tracker
        tracker = pathtracker(F; parameters=p, generic_parameters=p0)

        @test last(tracker.core_tracker.homotopy.homotopy.p) isa Vector{ComplexF64}
        @test first(tracker.core_tracker.homotopy.homotopy.p) isa Vector{ComplexF64}

        #Let's solve another system using PathTracking
        p_current = randn(6)
        result = track(tracker, S0[1]; target_parameters=p_current, accuracy=1e-9)
        sol = solution(result)
        F_SP = SPSystem(F; parameters=p)
        @test residual(result) ≈ euclidean_norm(F_SP(sol, p_current)) atol=1e-13
    end
end


text:
export time_modeling

# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to devito
function time_modeling(model_full::Model, srcGeometry, srcData, recGeometry, recData, dm, op::Char, mode::Int64, options)
    # Load full geometry for out-of-core geometry containers
    recGeometry = Geometry(recGeometry)
    srcGeometry = Geometry(srcGeometry)

    # Reutrn directly for J*0
    if op=='J' && mode == 1
        if norm(dm) == 0 && options.return_array == false
            return judiVector(recGeometry, zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
        elseif norm(dm) == 0 && options.return_array == true
            return vec(zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
        end
    end

    # limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        model, dm = limit_model_to_receiver_area(srcGeometry, recGeometry, model, options.buffer_size; pert=dm)
    else
        model = model_full
    end

    # Set up Python model structure
    modelPy = devito_model(model, options; dm=dm)

    # Load shot record if stored on disk
    typeof(recData) == SegyIO.SeisCon && (recData = convert(Array{Float32,2}, recData[1].data))

    # Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
    recGeometry, recData = remove_out_of_bounds_receivers(recGeometry, recData, model)

    # Devito interface
    argout = devito_interface(modelPy, model, srcGeometry, srcData, recGeometry, recData, dm, options)

    # Extend gradient back to original model size
    if op=='J' && mode==-1 && options.limit_m==true
        argout = extend_gradient(model_full, model, argout)
    end

    return argout
end

# Function instance without options
time_modeling(model::Model, srcGeometry::Geometry, srcData, recGeometry::Geometry, recData, perturbation, srcnum::Int64, op::Char, mode::Int64) =
    time_modeling(model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum, op, mode, Options())

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using Tokens
using Test

function testvector()
end


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

@doc """
$(convertext(:angle,"angle(U,S)"))

Extent of one-dimensional angle or `angle` (rad), unit conversion factor.

```Julia
julia> angle(CGS,Metric) # rad⋅rad⁻¹
$(angle(CGS,Metric))
```
""" angle, A

@doc """
$(convertext(:solidangle,"angle(U,S)^2"))

Extent of two-dimensional angle or `solidangle` (rad²), unit conversion factor.

```Julia
julia> solidangle(CGS,Metric) # rad²⋅rad⁻²
$(solidangle(CGS,Metric))
```
""" solidangle

# spacetime

@doc """
$(convertext(:area,"length(U,S)^2"))

Extent of two-dimensional shape or `area` (m²), unit conversion factor.

```Julia
julia> area(CGS,Metric) # m²⋅cm⁻²
$(area(CGS,Metric))

julia> area(English,Metric) # m²⋅ft⁻²
$(area(English,Metric))

julia> area(Survey,English) # ft²⋅ftUS⁻²
$(area(Survey,English))
```
""" area

@doc """
$(convertext(:volume,"length(U,S)^3"))

Extent of three-dimensional shape or `volume` (m³), unit conversion factor.

```Julia
julia> volume(CGS,Metric) # m³⋅cm⁻³
$(volume(CGS,Metric))

julia> volume(English,Metric) # m³⋅ft⁻³
$(volume(English,Metric))

julia> volume(Survey,English) # ft³⋅ftUS⁻³
$(volume(Survey,English))
```
""" volume

@doc """
$(convertext(:wavenumber,"1/length(U,S)"))

Number of occurences per unit of space (m⁻¹), unit conversion factor.

```Julia
julia> wavenumber(CGS,Metric) # cm⋅m⁻¹
$(wavenumber(CGS,Metric))

julia> wavenumber(English,Metric) # ft⋅m⁻¹
$(wavenumber(English,Metric))
```
""" wavenumber

@doc """
$(convertext(:angularwavenumber,"angle(U,S)/length(U,S)"))

Number of occurences per unit of space (m⁻¹), unit conversion factor.

```Julia
julia> angularwavenumber(CGS,Metric) # cm⋅m⁻¹
$(angularwavenumber(CGS,Metric))

julia> angularwavenumber(English,Metric) # ft⋅m⁻¹
$(angularwavenumber(English,Metric))
```
""" angularwavenumber

@doc """
$(convertext(:fuelefficiency,"1/area(U,S)"))

Distance per volume or `fuelefficiency` (m⋅m⁻³, m⁻²), unit conversion factor.

```Julia
julia> fuelefficiency(CGS,Metric) # cm²⋅m⁻²
$(fuelefficiency(CGS,Metric))

julia> fuelefficiency(English,Metric) # ft²⋅m⁻²
$(fuelefficiency(English,Metric))
```
""" fuelefficiency

@doc """
$(convertext(:numberdensity,"1/volume(U,S)"))

Number per `volume` or `numberdensity` (m⁻³ or ft⁻³), unit conversion factor.

```Julia
julia> numberdensity(CGS,Metric) # cm³⋅m⁻³
$(numberdensity(CGS,Metric))

julia> numberdensity(English,Metric) # ft³⋅m⁻³
$(numberdensity(English,Metric))
```
""" numberdensity

@doc """
$(convertext(:frequency,"1/time(U,S)"))

Number of occurences per unit of time (Hz or s⁻¹), unit conversion factor.

```Julia
julia> frequency(IAU,Metric) day⋅s⁻¹
$(frequency(IAU,Metric))
```
""" frequency

@doc """
$(convertext(:angularfrequency,"angle(U,S)/time(U,S)"))

Circular radian frequency (rad⋅Hz or rad⋅s⁻¹), unit conversion factor.

```Julia
julia> angularfrequency(IAU,Metric) day⋅s⁻¹
$(frequency(IAU,Metric))
```
""" angularfrequency

@doc """
$(convertext(:frequencydrift,"1/time(U,S)^2"))

Drift of `frequency` per `time` or `frequencydrift` (Hz⋅s⁻¹, s⁻²), unit conversion factor.
```Julia
julia> frequencydrift(IAU,Metric) day²⋅Hz⋅s⁻¹
$(frequencydrift(IAU,Metric))
```
""" frequencydrift

@doc """
$(convertext(:speed,"lightspeed(S)/lightspeed(U)"))

Velocity or `length` per `time` or `speed` (m⋅s⁻¹), unit conversion factor.

```Julia
julia> speed(CGS,Metric) # m⋅cm⁻¹
$(speed(CGS,Metric))

julia> speed(IAU,Metric) # m⋅day⋅s⁻¹⋅au⁻¹
$(speed(IAU,Metric))

julia> speed(English,Metric) # m⋅ft⁻¹
$(speed(English,Metric))

julia> speed(Survey,English) # ft⋅ftUS⁻¹
$(speed(Survey,English))
```
""" speed

@doc """
$(convertext(:acceleration,"speed(U,S)/time(U,S)"))

Specific force or `speed` per `time` or `acceleration` (m⋅s⁻²), unit conversion factor.

```Julia
julia> acceleration(CGS,Metric) # m⋅s⁻¹⋅gal⁻¹
$(acceleration(CGS,Metric))

julia> acceleration(IAU,Metric) # m⋅day²⋅s⁻²⋅au⁻¹
$(acceleration(IAU,Metric))

julia> acceleration(English,Metric) # m⋅ft⁻¹
$(acceleration(English,Metric))

julia> acceleration(Survey,English) # ft⋅ftUS⁻¹
$(acceleration(Survey,English))
```
""" acceleration

@doc """
$(convertext(:jerk,"speed(U,S)/time(U,S)^2"))

Jolt or `acceleration` per `time` or `jerk` (m⋅s⁻³), unit conversion factor.

```Julia
julia> jerk(CGS,Metric) # m⋅cm⁻¹
$(jerk(CGS,Metric))

julia> jerk(IAU,Metric) # m⋅day³⋅s⁻³⋅au⁻¹
$(jerk(IAU,Metric))

julia> jerk(English,Metric) # m⋅ft⁻¹
$(jerk(English,Metric))

julia> jerk(Survey,English) # ft⋅ftUS⁻¹
$(jerk(Survey,English))
```
""" jerk

@doc """
$(convertext(:snap,"speed(U,S)/time(U,S)^3"))

Jounce or `jerk` per `time` or `snap` (m⋅s⁻⁴), unit conversion factor.

```Julia
julia> snap(CGS,Metric) # m⋅cm⁻¹
$(snap(CGS,Metric))

julia> snap(IAU,Metric) # m⋅day⁴⋅s⁻⁴⋅au⁻¹
$(snap(IAU,Metric))

julia> snap(English,Metric) # m⋅ft⁻¹
$(snap(English,Metric))

julia> snap(Survey,English) # ft⋅ftUS⁻¹
$(snap(Survey,English))
```
""" snap

@doc """
$(convertext(:crackle,"speed(U,S)/time(U,S)^4"))

A `snap` per `time` or `crackle` (m⋅s⁻⁵), unit conversion factor.

```Julia
julia> crackle(CGS,Metric) # m⋅cm⁻¹
$(crackle(CGS,Metric))

julia> crackle(IAU,Metric) # m⋅day⁵⋅s⁻⁵⋅au⁻¹
$(crackle(IAU,Metric))

julia> crackle(English,Metric) # m⋅ft⁻¹
$(crackle(English,Metric))

julia> crackle(Survey,English) # ft⋅ftUS⁻¹
$(crackle(Survey,English))
```
""" crackle

@doc """
$(convertext(:pop,"speed(U,S)/time(U,S)^5"))

A `crackle` per `time` or `pop` (m⋅s⁻⁶), unit conversion factor.

```Julia
julia> pop(CGS,Metric) # m⋅cm⁻¹
$(pop(CGS,Metric))

julia> pop(IAU,Metric) # m⋅day⁶⋅s⁻⁶⋅au⁻¹
$(pop(IAU,Metric))

julia> pop(English,Metric) # m⋅ft⁻¹
$(pop(English,Metric))

julia> pop(Survey,English) # ft⋅ftUS⁻¹
$(pop(Survey,English))
```
""" pop

@doc """
$(convertext(:volumeflow,"area(U,S)*speed(U,S)"))

Volumetric flow rate or `volumeflow` (m³⋅s⁻¹), unit conversion factor.

```Julia
julia> volumeflow(CGS,Metric) # m³⋅cm⁻³
$(volume(CGS,Metric))

julia> volumeflow(English,Metric) # m³⋅ft⁻³
$(volume(English,Metric))

julia> volumeflow(Survey,English) # ft³⋅ftUS⁻³
$(volume(Survey,English))
```
""" volumeflow

@doc """
$(convertext(:specificenergy,"speed(U,S)^2"))

Massic energy or `energy` per `mass` or `specificenergy` (J⋅kg⁻¹), unit conversion factor.

```Julia
julia> specificenergy(CGS,Metric) # m²⋅cm⁻²
$(specificenergy(CGS,Metric))

julia> specificenergy(IAU,Metric) # m²⋅day²⋅s⁻²⋅au⁻²
$(specificenergy(IAU,Metric))

julia> specificenergy(English,Metric) # m²⋅ft⁻²
$(specificenergy(English,Metric))

julia> specificenergy(Survey,English) # ft²⋅ftUS⁻²
$(specificenergy(Survey,English))
```
""" specificenergy

# kinematic

@doc """
$(convertext(:mass,"electronmass(S)/electronmass(U)"))

Inertal `mass` or matter quantity or resistance to aceleration (kg), unit conversion factor.

```Julia
julia> mass(CGS,Metric) # kg⋅g⁻¹
$(mass(CGS,Metric))

julia> mass(CODATA,Metric) # kg⋅kg⁻¹
$(mass(CODATA,Metric))

julia> mass(Conventional,Metric) # kg⋅kg⁻¹
$(mass(Conventional,Metric))

julia> mass(English,Metric) # kg⋅slug⁻¹
$(mass(English,Metric))

julia> mass(IAU,Metric) # kg⋅m⊙⁻¹
$(mass(IAU,Metric))

julia> mass(PlanckGauss,Metric) # kg⋅mP⁻¹
$(mass(PlanckGauss,Metric))
```
""" mass, M

@doc """
$(convertext(:inertia,"mass(U,S)/gravity(U,S)"))

Inertal `mass` or matter quantity or resistance to aceleration (kg), unit conversion factor.

```Julia
julia> inertia(CGS,Metric) # kg⋅g⁻¹
$(inertia(CGS,Metric))

julia> inertia(CODATA,Metric) # kg⋅kg⁻¹
$(inertia(CODATA,Metric))

julia> inertia(Conventional,Metric) # kg⋅kg⁻¹
$(inertia(Conventional,Metric))

julia> inertia(English,Metric) # kg⋅slug⁻¹
$(inertia(English,Metric))

julia> inertia(IAU,Metric) # kg⋅m⊙⁻¹
$(inertia(IAU,Metric))

julia> inertia(PlanckGauss,Metric) # kg⋅mP⁻¹
$(inertia(PlanckGauss,Metric))
```
""" inertia

@doc """
$(convertext(:energy,"mass(U,S)*specificenergy(U,S)"))

Work or heat or `energy` (J, N⋅m, kg⋅m²⋅s⁻²), unit conversion factor.

```Julia
julia> energy(CGS,Metric) # J⋅erg⁻¹
$(energy(CGS,Metric))

julia> energy(CGS,English) # ft⋅lb⋅erg⁻¹
$(energy(CGS,English))

julia> energy(English,Metric) # J⋅ft⁻¹⋅lb⁻¹
$(energy(English,Metric))

julia> 0.001/3600 # J⋅kW⁻¹⋅h⁻¹
$(0.001/3600)

julia> 1/elementarycharge(SI2019) # J⋅eV⁻¹
$(inv(elementarycharge(SI2019)))
```
""" energy
#julia> inv(kilocalorie(SI2019) # J⋅kcalₜₕ⁻¹
#$(inv(kilocalorie(SI2019)))
#julia> 1/BTUJ # J⋅BTU⁻¹
#$(1/1055.036345118633)

@doc """
$(convertext(:power,"energy(U,S)/time(U,S))"))

Radiant flux or `power` or `energy` per `time` (W, J⋅s⁻¹, kg⋅m²⋅s⁻³), unit conversion factor.

```Julia
julia> power(CGS,Metric) # W⋅s⋅erg⁻¹
$(power(CGS,Metric))

julia> power(English,Metric) # W⋅s⋅ft⁻¹⋅lb⁻¹
$(power(English,Metric))
```
""" power

@doc """
$(convertext(:force,"inertia(U,S)*acceleration(U,S)"))

Weight or force or `inertia` times `acceleration` (N, kg⋅m⋅s⁻²), unit conversion factor.

```Julia
julia> force(CGS,Metric) # N⋅dyn⁻¹
$(force(CGS,Metric))

julia> force(CGS,English) # lb⋅dyn⁻¹
$(force(CGS,English))

julia> force(English,Metric) # N⋅lb⁻¹
$(force(English,Metric))

julia> force(FPS,Metric) # pdl⋅N⁻¹
$(force(FPS,Metric))

julia> force(MetricEngineering,Metric) # kp⋅N⁻¹
$(force(MetricEngineering,Metric))
```
""" force, F

@doc """
$(convertext(:specificforce,"acceleration(U,S)/gravity(U,S)"))

Weight or `force` per `mass` or `gforce` (N/kg, m⋅s⁻²), unit conversion factor.
```Julia
julia> specificforce(CGS,Metric)
$(specificforce(CGS,Metric))

julia> specificforce(MetricEngineering,Metric)
$(specificforce(MetricEngineering,Metric))

julia> specificforce(English,Metric)
$(specificforce(English,Metric))
```
""" specificforce

@doc """
$(convertext(:gravityforce,"acceleration(U,S)/specificforce(U,S)"))

Reference `acceleration` per `specificforce` (𝟏, F⁻¹MLT⁻²), unit conversion factor.
```Julia
julia> gravityforce(Metric,CGS)
$(gravityforce(Metric,CGS))

julia> gravityforce(Metric,MetricEngineering)
$(gravityforce(Metric,MetricEngineering))

julia> gravityforce(Metric,English)
$(gravityforce(Metric,English))
```
""" gravityforce

@doc """
$(convertext(:pressure,"force(U,S)/area(U,S)"))

Pressure or stress or `force` per `area` (Pa, N⋅m⁻², kg⋅m⁻¹⋅s⁻²), unit conversion factor.

```Julia
julia> pressure(CGS,Metric) # Pa⋅Ba⁻¹
$(pressure(CGS,Metric))

julia> 1/atm # Pa⋅atm⁻¹
$(inv(atm))

julia> pressure(English,Metric) # Pa⋅ft²⋅lb⁻¹
$(pressure(English,Metric))

julia> pressure(Metric,IPS) # psi⋅Pa⁻¹
$(pressure(Metric,IPS))
```
""" pressure

# mechanical

@doc """
$(convertext(:impulse,"force(U,S)*time(U,S)"))

Linear `impulse` or `force` times `time` (N⋅s, kg⋅m⋅s⁻¹), unit conversion factor.

```Julia
julia> impulse(CGS,Metric) # N⋅dyn⁻¹
$(impulse(CGS,Metric))

julia> impulse(CGS,English) # lb⋅dyn⁻¹
$(impulse(CGS,English))

julia> impulse(English,Metric) # N⋅lb⁻¹
$(impulse(English,Metric))
```
""" impulse

@doc """
$(convertext(:momentum,"mass(U,S)*speed(U,S)"))

Linear `momentum` or `mass` times `speed` (N⋅s, kg⋅m⋅s⁻¹), unit conversion factor.

```Julia
julia> momentum(CGS,Metric) # N⋅dyn⁻¹
$(momentum(CGS,Metric))

julia> momentum(CGS,English) # lb⋅dyn⁻¹
$(momentum(CGS,English))

julia> momentum(British,Metric) # N⋅lb⁻¹
$(momentum(British,Metric))
```
""" momentum

@doc """
$(convertext(:angularmomentum,"momentum(U,S)*lengt(U,S)*angle(U,S)"))

Rotational momentum or `angularmomentum` (N⋅m⋅s, kg⋅m²⋅s⁻¹), unit conversion factor.

```Julia
julia> momentum(CGS,Metric) # N⋅m⋅dyn⁻¹⋅cm⁻¹
$(momentum(CGS,Metric))

julia> momentum(CGS,English) # lb⋅ft⋅dyn⁻¹⋅cm⁻¹
$(momentum(CGS,English))

julia> momentum(British,Metric) # N⋅m⋅lb⁻¹⋅ft⁻¹
$(momentum(British,Metric))
```
""" angularmomentum

@doc """
$(convertext(:yank,"mass(U,S)*jerk(U,S)"))

Rate of change of `force` or `yank` (N⋅s⁻¹, kg⋅m⋅s⁻³), unit conversion factor.

```Julia
julia> yank(CGS,Metric) # N⋅dyn⁻¹
$(yank(CGS,Metric))

julia> yank(CGS,English) # lb⋅dyn⁻¹
$(yank(CGS,English))

julia> yank(British,Metric) # N⋅lb⁻¹⋅
$(yank(British,Metric))
```
""" yank

@doc """
$(convertext(:areadensity,"mass(U,S)/area(U,S)"))

Surface or `areadensity` or `mass` per `area` (kg⋅m⁻²), unit conversion factor.

```Julia
julia> areadensity(CGS,Metric) # kg⋅cm²⋅g⁻¹⋅m⁻²
$(areadensity(CGS,Metric))

julia> areadensity(CGS,English) # lb⋅cm²⋅g⁻¹⋅ft⁻²
$(areadensity(CGS,English))

julia> areadensity(English,Metric) # kg⋅ft²⋅lb⁻¹⋅m⁻²
$(areadensity(English,Metric))

julia> areadensity(British,Metric) # kg⋅ft²⋅slug⁻¹⋅m⁻²
$(areadensity(British,Metric))
```
""" areadensity

@doc """
$(convertext(:density,"mass(U,S)/volume(U,S)"))

Specific mass or `mass` per `volume` or `density` (kg⋅m⁻³), unit conversion factor.

```Julia
julia> density(CGS,Metric) # kg⋅cm³⋅g⁻¹⋅m⁻³
$(density(CGS,Metric))

julia> density(CGS,Brtish) # slug⋅cm³⋅g⁻¹⋅ft⁻³
$(density(CGS,British))

julia> density(English,Metric) # kg⋅ft³⋅lb⁻¹⋅m⁻³
$(density(English,Metric))
```
""" density

@doc """
$(convertext(:specificweight,"force(U,S)/volume(U,S)"))

Specific weight or `force` per `volume` (N⋅m⁻³ or lb⋅ft⁻³), unit conversion factor.
```Julia
julia> specificweight(CGS,Metric) # N⋅cm³⋅dyn⁻¹⋅m⁻³
$(specificweight(CGS,Metric))

julia> specificweight(CGS,Brtish) # lb⋅cm³⋅dyn⁻¹⋅ft⁻³
$(specificweight(CGS,British))

julia> specificweight(English,Metric) # N⋅ft³⋅lb⁻¹⋅m⁻³
$(specificweight(English,Metric))
```
""" specificweight

@doc """
$(convertext(:specificvolume,"volume(U,S)/mass(U,S)"))

Reciprocal `density` or `volume` per `mass` or `specificvolume` (m³⋅kg), unit conversion factor.

```Julia
julia> specificvolume(CGS,Metric) # g⋅m³⋅kg⁻¹⋅cm⁻³
$(specificvolume(CGS,Metric))

julia> specificvolume(CGS,British) # kg⋅ft³⋅slug⁻¹⋅cm⁻³
$(specificvolume(CGS,British))

julia> specificvolume(English,Metric) # lb⋅m³⋅kg⁻¹⋅ft⁻³
$(specificvolume(English,Metric))
```
""" specificvolume

@doc """
$(convertext(:action,"energy(U,S)*time(U,S)"))

Integrated `momentum` over `length` or `action` (J⋅s, N⋅m⋅s), unit conversion factor.

```Julia
julia> action(CGS,Metric) # J⋅erg⁻¹
$(action(CGS,Metric))

julia> action(CGS,English) # ft⋅lb⋅erg⁻¹
$(action(CGS,English))

julia> action(English,Metric) # J⋅ft⁻¹⋅lb⁻¹
$(action(English,Metric))
```
""" action

#=@doc """
$(convertext(:stiffness,"force(U,S)/length(U,S)"))

Amount of `force` per `length` or `stiffness` (N⋅m⁻¹, J⋅m⁻², kg⋅s⁻²), unit conversion factor.

```Julia
julia> stiffness(CGS,Metric) # kg⋅g⁻¹
$(stiffness(CGS,Metric))

julia> stiffness(CGS,English) # lb⋅g⁻¹
$(stiffness(CGS,English))

julia> stiffness(English,Metric) # kg⋅lb⁻¹
$(stiffness(English,Metric))
```
""" stiffness=#

@doc """
$(convertext(:intensity,"power(U,S)/area(U,S)"))

Heat flux density or irradiance or `power` per `area` (W⋅m⁻², kg⋅s⁻³), unit conversion factor.

```Julia
julia> intensity(CGS,Metric) # kg⋅g⁻¹
$(intensity(CGS,Metric))

julia> intensity(CGS,English) # lb⋅g⁻¹
$(intensity(CGS,English))

julia> intensity(English,Metric) # kg⋅lb⁻¹
$(intensity(English,Metric))
```
""" intensity

@doc """
$(convertext(:diffusivity,"speed(U,S)*length(U,S)"))

Thermal `diffusivity` or kinematic viscostiy (m²⋅s⁻¹), unit conversion factor.

```Julia
julia> diffusivity(CGS,Metric) # m²⋅cm⁻²
$(diffusivity(CGS,Metric))

julia> diffusivity(English,Metric) # m²⋅ft⁻²
$(diffusivity(English,Metric))

julia> diffusivity(Survey,English) # ft²⋅ftUS⁻²
$(diffusivity(Survey,English))
```
""" diffusivity

@doc """
$(convertext(:viscosity,"inertia(U,S)/length(U,S)/time(U,S)"))

Resistance to deformation or `viscosity` (Pa⋅s, kg⋅m⁻¹⋅s⁻¹), unit conversion factor.

```Julia
julia> viscosity(CGS,Metric) # Pa⋅Ba⁻¹
$(viscosity(CGS,Metric))

julia> viscosity(English,Metric) # Pa⋅ft²⋅lb⁻¹
$(viscosity(English,Metric))

julia> viscosity(British,Metric) # Pa⋅ft²⋅lb⁻¹
$(viscosity(British,Metric))
```
""" viscosity

@doc """
$(convertext(:lineardensity,"mass(U,S)/length(U,S)"))

Amount of `lineardensity` or `mass` per `length` (kg⋅m⁻¹), unit conversion factor.

```Julia
julia> lineardensity(CGS,Metric) # kg⋅cm¹⋅g⁻¹⋅m⁻¹
$(lineardensity(CGS,Metric))

julia> lineardensity(CGS,British) # slug⋅cm¹⋅g⁻¹⋅ft⁻¹
$(lineardensity(CGS,British))

julia> lineardensity(English,Metric) # kg⋅ft¹⋅lb⁻¹⋅m⁻¹
$(lineardensity(English,Metric))
```
""" lineardensity

@doc """
$(convertext(:massflow,"mass(U,S)/time(U,S)"))

Rate of `massflow` or `mass` per `time` (kg⋅s⁻¹), unit conversion factor.

```Julia
julia> massflow(CGS,Metric) # kg⋅g⁻¹
$(massflow(CGS,Metric))

julia> massflow(CODATA,Metric) # kg⋅kg⁻¹
$(massflow(CODATA,Metric))

julia> massflow(Conventional,Metric) # kg⋅kg⁻¹
$(massflow(Conventional,Metric))

julia> massflow(English,Metric) # kg⋅slug⁻¹
$(massflow(English,Metric))
```
""" massflow

@doc """
$(convertext(:spectralflux,"power(U,S)/length(U,S)"))

Spectral power or `power` per wave `length` (W⋅m⁻¹), unit conversion factor.

```Julia
julia> spectralflux(CGS,Metric) # kg⋅m⋅g⁻¹⋅cm⁻¹
$(spectralflux(CGS,Metric))

julia> spectralflux(CGS,English) # lb⋅ft⋅g⁻¹⋅cm⁻¹
$(spectralflux(CGS,English))

julia> spectralflux(English,Metric) # kg⋅m⋅lb⁻¹⋅ft⁻¹
$(spectralflux(English,Metric))
```
""" spectralflux

@doc """
$(convertext(:powerdensity,"power(U,S)/volume(U,S)"))

Spectral irradiance (volume) or `powerdensity` (W⋅m⁻³), unit conversion factor.

```Julia
julia> powerdensity(CGS,Metric) # kg⋅cm⋅g⁻¹⋅m⁻¹
$(powerdensity(CGS,Metric))

julia> powerdensity(CGS,English) # lb⋅cm⋅g⁻¹⋅ft⁻¹
$(powerdensity(CGS,English))

julia> powerdensity(English,Metric) # kg⋅ft⋅lb⁻¹⋅m⁻¹
$(powerdensity(English,Metric))
```
""" powerdensity

@doc """
$(convertext(:compressibility,"1/pressure(U,S)"))

Relative volume change or `compressibility` (Pa⁻¹), unit conversion factor.

```Julia
julia> compressibility(CGS,Metric) # Ba⋅Pa⁻¹
$(compressibility(CGS,Metric))

julia> compressibility(English,Metric) # lb⋅ft⁻²⋅Pa⁻¹
$(compressibility(English,Metric))

julia> compressibility(Metric,IPS) # Pa⋅psi⁻¹
$(compressibility(Metric,IPS))
```
""" compressibility

@doc """
$(convertext(:fluence,"energy(U,S)/area(U,S"))

Radiant exposure or `force` per `length` or stiffness (N⋅m⁻¹, J⋅m⁻²), unit conversion factor.

```Julia
julia> fluence(CGS,Metric) # kg⋅g⁻¹
$(mass(CGS,Metric))

julia> fluence(CGS,English) # lb⋅g⁻¹
$(fluence(CGS,English))

julia> fluence(CODATA,Metric) # kg⋅kg⁻¹
$(mass(CODATA,Metric))

julia> fluence(Conventional,Metric) # kg⋅kg⁻¹
$(fluence(Conventional,Metric))

julia> fluence(English,Metric) # kg⋅lb⁻¹
$(fluence(English,Metric))
```
""" fluence

@doc """
$(convertext(:rotationalinertia,"mass(U,S)*area(U,S)"))

Moment of inertia or `rotationalinertia` (kg⋅m²), unit conversion factor.

```Julia
julia> rotationalinertia(CGS,Metric) # kg⋅m²⋅g⁻¹⋅cm⁻²
$(rotationalinertia(CGS,Metric))

julia> rotationalinertia(CGS,British) # slug⋅ft²⋅g⁻¹⋅cm⁻²
$(rotationalinertia(CGS,British))

julia> rotationalinertia(English,Metric) # kg⋅m²⋅lb⁻¹⋅ft⁻²
$(rotationalinertia(English,Metric))
```
""" rotationalinertia

# acoustic

@doc """
$(convertext(:soundexposure,"pressure(U,S)^2*time(U,S)"))

Square of `pressure` by `time` or `soundexposure` (Pa²⋅s, N²⋅m⁻⁴), unit conversion factor.

```Julia
julia> soundexposure(CGS,Metric) # Pa²⋅Ba⁻²
$(soundexposure(CGS,Metric))

julia> soundexposure(English,Metric) # Pa²⋅ft⁴⋅lb⁻²
$(soundexposure(English,Metric))
```
""" soundexposure

@doc """
$(convertext(:specificimpedance,"pressure(U,S)/speed(U,S)"))

Characteristic specific acoustic impedance (Rayl, Pa⋅s⋅m⁻¹), unit conversion factor.

```Julia
julia> specificimpedance(CGS,Metric) # Pa⋅cm⋅m⁻¹⋅Ba⁻¹
$(specificimpedance(CGS,Metric))

julia> specificimpedance(English,Metric) # Pa⋅ft³⋅m⁻¹⋅lb⁻¹
$(specificimpedance(English,Metric))
```
""" specificimpedance

@doc """
$(convertext(:impedance,"specificimpedance(U,S)/area(U,S)"))

Acoustic `impedance` (Rayl⋅m⁻², Pa⋅s⋅m⁻³, kg⋅s⁻¹⋅m⁻⁴), unit conversion factor.

```Julia
julia> impedance(CGS,Metric) # Pa⋅cm³⋅m⁻³⋅Ba⁻¹
$(impedance(CGS,Metric))

julia> impedance(English,Metric) # Pa⋅ft⁵⋅m⁻³⋅lb⁻¹
$(impedance(English,Metric))
```
""" impedance(U::UnitSystem,S::UnitSystem)

@doc """
$(convertext(:admittance,"area(U,S)/specificimpedance(U,S)"))

Acoustic `admittance` (m²⋅Rayl⁻¹, m³⋅s⁻¹⋅Pa⁻¹, m⁴⋅s⋅kg⁻¹), unit conversion factor.

```Julia
julia> admittance(CGS,Metric) # Ba⋅m³⋅cm⁻³⋅Pa⁻¹
$(admittance(CGS,Metric))

julia> admittance(English,Metric) # lb⋅m³⋅ft⁻⁵⋅Pa⁻¹
$(admittance(English,Metric))
```
""" admittance

@doc """
$(convertext(:compliance,"time(U,S)^2/mass(U,S)"))

Acoustic `compliance` is reciprocal of `fluence` (m⋅N⁻¹, m³⋅Pa⁻¹), unit conversion factor.

```Julia
julia> compliance(CGS,Metric) # kg⋅g⁻¹
$(compliance(CGS,Metric))

julia> compliance(CGS,English) # slug⋅g⁻¹
$(compliance(CGS,English))

julia> compliance(English,Metric) # kg⋅lb⁻¹
$(compliance(English,Metric))
```
""" compliance

@doc """
$(convertext(:inertance,"mass(U,S)/length(U,S)^4"))

Acoustic mass or `inertance` (kg⋅m⁴, Pa⋅s²⋅m⁻³), unit conversion factor.

```Julia
julia> inertance(CGS,Metric) # kg⋅cm⁴⋅g⁻¹⋅m⁻⁴
$(inertance(CGS,Metric))

julia> inertance(CGS,English) # slug⋅cm⁴⋅g⁻¹⋅ft⁻⁴
$(inertance(CGS,English))

julia> inertance(English,Metric) # kg⋅ft⁴⋅lb⁻¹⋅m⁻⁴
$(inertance(English,Metric))
```
""" inertance


text:
using Reactive
using StaticArrays,Rotations
const PM = StaticArrays.SArray{Tuple{4,4},Float32,2,16}
(::Type{PM})(R, T) = begin
    s_R = convert(SArray{Tuple{3,3},Float32,2,9},R)
    s_T = convert(SVector{3,Float32},T)
    PM(s_R,s_T)
end
(::Type{PM})(R::SArray{Tuple{3,3},Float32,2,9}, T::SVector{3,Float32}) = begin
    T0 = zero(Float32)
    T1 = one(Float32)
    @SMatrix [R[1,1] R[1,2] R[1,3] T[1];
              R[2,1] R[2,2] R[2,3] T[2];
              R[3,1] R[3,2] R[3,3] T[3];
              T0     T0     T0      T1];
end
export PM

const SPM = Signal{PM}
export SPM

using GeometryTypes
import Base.*
*(M::PM,p::Union{Point3,Vec3,SVector{3,Float32}}) = begin
    p = M * Vec4(p[1],p[2],p[3],1)
    Point3(p[1]/p[4],p[2]/p[4],p[3]/p[4])
end

import Base.convert
convert(::Type{SPM},x::Matrix) = Signal(PM(x))
convert(::Type{Signal{T}},x::T) where T = Signal(x)

#TODO deprecate this in vafor of vec
decompose_PM(trans::PM) = begin
    R = SMatrix{3,3}(trans[1:3,1:3])
    R = RotXYZ(R)
    angles = SVector{3}(R.theta1,R.theta2,R.theta3)
    T = SVector{3}(trans[1:3,4])
    (angles,T)
end

#TODO ue vec instead pf decompose
function lerp_PM(start_trans::PM,end_trans::PM,alpha)
    (o_angles,o_T) = decompose_PM(start_trans)
    (n_angles,n_T) = decompose_PM(end_trans)
    angles = alpha * n_angles + (1-alpha) * o_angles
    T = alpha * n_T + (1-alpha) * o_T
    return PM(Array(RotXYZ(angles...)),T)
end
export lerp_PM

import Base.vec
vec(trans::PM) = begin
    R = RotXYZ(trans[1:3,1:3])
    angles = SVector{3}(R.theta1,R.theta2,R.theta3)
    T = trans[1:3,4]
    [angles...,T...]
end

(::Type{PM})(V::AbstractVector) = begin
    if length(V) != 6
        error("6 elements required to construct PM(Projection Matrix) from vector")
    end
    PM(RotXYZ(V[1],V[2],V[3]),V[4:6])
end


text:
using Gadfly, DataArrays, RDatasets

if VERSION < v"0.4-dev"
    using Dates
else
    using Base.Dates
end

economics = dataset("HistData", "Prostitutes")
dates = Date[Date(d) for d in economics[:Date]]
economics[:Date] = dates

p = plot(economics, x=:Date, y=:Count, Geom.line)


text:
using Genie.Renderer, Genie.Requests
using Random
import Base.Filesystem: mktemp

greeting = "Welcome"
name = "Genie"

function htmlviewfile_withvars()
  raw"

$(@vars(:greeting))

This is a $(@vars(:name)) test

  "
end

function htmltemplatefile_withvars()
  raw"
  <!DOCTYPE HTML>

    <title>$(@vars(:name)) test</title>

    <% @yield %>

    Just a footer

  "
end

viewfile = mktemp()
write(viewfile[2], htmlviewfile_withvars())
close(viewfile[2])

templatefile = mktemp()
write(templatefile[2], htmltemplatefile_withvars())
close(templatefile[2])

@testset "HTML rendering with view files" begin
  r = Requests.HTTP.Response()

  @testset "HTML rendering with view file no layout with vars" begin
    r = html(Renderer.FilePath(viewfile[1]), greeting = greeting, name = Genie)

    @test String(r.body) == "
$greeting

This is a $name test

"
  end;

  @testset "HTML rendering with view file and layout with vars" begin
    r = html(Renderer.FilePath(viewfile[1]), layout = Renderer.FilePath(templatefile[1]), greeting = greeting, name = Genie)

    @test String(r.body) == "<title>$name test</title>

$greeting

This is a $name test

\n
Just a footer"
  end;

  @test r.status == 200
  @test r.headers[1]["Content-Type"] == "text/html; charset=utf-8"

  @testset "HTML rendering with view file no layout with vars custom headers" begin
    r = html(Renderer.FilePath(viewfile[1]), headers = Renderer.HTTPHeaders("Cache-Control" => "no-cache"), greeting = greeting, name = Genie)

    @test String(r.body) == "
$greeting

This is a $name test

"
    @test r.headers[1]["Cache-Control"] == "no-cache"
  end;

  @testset "HTML rendering with view file and layout with vars custom headers" begin
    r = html(Renderer.FilePath(viewfile[1]), layout = Renderer.FilePath(templatefile[1]), headers = Renderer.HTTPHeaders("Cache-Control" => "no-cache"), greeting = greeting, name = Genie)

    @test String(r.body) == "<title>$name test</title>

$greeting

This is a $name test

\n
Just a footer"
    @test r.headers[1]["Cache-Control"] == "no-cache"
  end;

  @testset "HTML rendering with view file no layout with vars custom headers custom status" begin
    r = html(Renderer.FilePath(viewfile[1]), headers = Renderer.HTTPHeaders("Cache-Control" => "no-cache"), status = 500, greeting = greeting, name = Genie)

    @test String(r.body) == "
$greeting

This is a $name test

"
    @test r.headers[1]["Cache-Control"] == "no-cache"
    @test r.status == 500
  end;

  @testset "HTML rendering with view file and layout with vars custom headers custom status" begin
    r = html(Renderer.FilePath(viewfile[1]), layout = Renderer.FilePath(templatefile[1]),
              headers = Renderer.HTTPHeaders("Cache-Control" => "no-cache"), status = 404, greeting = greeting, name = Genie)

    @test String(r.body) == "<title>$name test</title>

$greeting

This is a $name test

\n
Just a footer"
    @test r.headers[1]["Cache-Control"] == "no-cache"
    @test r.status == 404
  end;
end;


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

name = "SoftPosit"
version = v"0.4.1"

# Collection of sources required to build SoftPosit
sources = [
    "https://gitlab.com/cerlane/SoftPosit/-/archive/$version/SoftPosit-$version.tar.gz" =>
    "13f7360c5b91ad3704f66537a754ba3748a764e1291eaa33940866ca37c7dbf5",
    "./patches"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SoftPosit-*/build/Linux-x86_64-GCC/
patch -p1 < $WORKSPACE/srcdir/makefile.patch
if [[ "${target}" == *darwin* ]]; then
    make SLIB=".dylib" julia
elif [[ "${target}" == *mingw* ]]; then
    make SLIB=".dll" julia
else
    make julia
fi
mkdir -p $prefix/lib
cp softposit.* $prefix/lib/.
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [

]

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


text:
using NIDAQ, Base.Test

info("NIDAQ: Checking installation - no measurement hardware is needed for this section")
@test typeof(getproperties()) == Dict{ASCIIString,Tuple{Any,Bool}}
@test haskey(getproperties(), "NIDAQMajorVersion")
d = devices()

if length(d) == 0
    info("NIDAQ: no measurement HW found, skipping rest of the test suite")
    device_available = false
else
    d=d[1]
    info("NIDAQ: found at least one measurement HW: $(d)")
    @test typeof(getproperties(d)) == Dict{ASCIIString,Tuple{Any,Bool}}
    pr = getproperties(d)
    device_available = true

    @test typeof(getproperties(d)) == Dict{ASCIIString,Tuple{Any,Bool}}
    @test typeof(analog_input_channels()) == Vector{ASCIIString}
    @test typeof(analog_output_channels()) == Vector{ASCIIString}
    @test typeof(digital_input_channels()) == Vector{ASCIIString}
    @test typeof(digital_output_channels()) == Vector{ASCIIString}
    @test typeof(counter_input_channels()) == Vector{ASCIIString}
    @test typeof(counter_output_channels()) == Vector{ASCIIString}
    @test typeof(analog_input_ranges()) == Matrix{Float64}
    @test typeof(analog_output_ranges()) == Matrix{Float64}
end

if !(device_available && length(pr["AIPhysicalChans"][1]) > 0)
    info("NIDAQ: measurement HW does not support any AIPhysicalChans, skipping part of the suite")
else
    t = analog_input(d*"/ai0")
    @test typeof(getproperties(t)) == Dict{ASCIIString,Tuple{Any,Bool}}
    @test typeof(getproperties(t,d*"/ai0")) == Dict{ASCIIString,Tuple{Any,Bool}}
    @test typeof(t) == NIDAQ.AITask
    @test start(t) == nothing
    @test length(read(t, Float64, 3)) == 3
    @test stop(t) == nothing
    @test analog_input(t, d*"/ai1") == nothing
    @test start(t) == nothing
    @test length(read(t, UInt32, 6)) == 12
    @test stop(t) == nothing
    @test clear(t) == nothing
end

if !(device_available && length(pr["AOPhysicalChans"][1]) > 0)
    info("NIDAQ: measurement HW does not support any AOPhysicalChans, skipping part of the suite")
else
    t = analog_output(d*"/ao0")
    @test typeof(t) == NIDAQ.AOTask
    @test start(t) == nothing
    @test write(t, rand(3)) == 3
    @test stop(t) == nothing
    @test analog_output(t, d*"/ao1") == nothing
    @test start(t) == nothing
    @test write(t, rand(UInt32,6,2)) == 6
    @test stop(t) == nothing
    @test clear(t) == nothing
end

if !(device_available && length(pr["DILines"][1]) > 0)
    info("NIDAQ: measurement HW does not support DILines, skipping part of the suite")
else
    t = digital_input(d*"/Port0/Line0")
    @test typeof(t) == NIDAQ.DITask
    @test start(t) == nothing
    @test length(read(t, UInt32, 3)) == 3
    @test stop(t) == nothing
    @test digital_input(t, d*"/Port0/Line1") == nothing
    @test start(t) == nothing
    @test length(read(t, UInt8, 6)) == 12
    @test stop(t) == nothing
    @test clear(t) == nothing
end
if !(device_available && length(pr["DOLines"][1]) > 0)
    info("NIDAQ: measurement HW does not support DILines, skipping part of the suite")
else
    t = digital_output(d*"/Port0/Line0")
    @test typeof(t) == NIDAQ.DOTask
    @test start(t) == nothing
    @test write(t, round(UInt32, [1,0,1,0,1,0])) == 6
    @test stop(t) == nothing
    @test digital_output(t, d*"/Port0/Line1") == nothing
    @test start(t) == nothing
    @test write(t, round(UInt8, [1 0; 0 0; 1 0; 0 1; 1 1; 0 1])) == 6
    @test stop(t) == nothing
    @test clear(t) == nothing
end

if !(device_available && length(pr["COPhysicalChans"][1]) > 0)
    info("NIDAQ: measurement HW does not support COPhysicalChans, skipping part of the suite")
else
    t = generate_pulses(d*"/ctr0")
    @test typeof(t) == NIDAQ.COTask
    @test NIDAQ.CfgImplicitTiming(t.th, NIDAQ.Val_FiniteSamps, UInt64(10)) == 0
    @test start(t) == nothing
    @test stop(t) == nothing
    @test clear(t) == nothing
end


text:
using Test, QuHamiltonian
using Lattices

@testset "TFIM on Chain Lattice" begin

J = 0.1
h = 0.1

@nearest i, j
@vertex l
h = J * @sum(Z[i]Z[j]) - h * @sum(X[l])
lattice = Chain(4)
H = h(lattice)

end

@testset "TFIM on Square Lattice" begin

J = 0.1
h = 0.1

@nearest i, j
H = J * @sum(Z[i]Z[j]) - h * @sum(X[i])
lattice = Sqaure(4, 4)
H(lattice)

end

@testset "J1J2 on Square Lattice" begin

@neighbor 1 i, j
@neighbor 2 k, l

J1, J2 = 0.1, 0.2

H = J1 * @sum(S[i]S[j]) - J2 * @sum(S[k]S[l])
lattice = Square(4, 4)
H(lattice)

end

@testset "Heisenberg Model with rand correlation" begin

lattice = Chain(4)

@nearrest i, j
J = i, j -> rand()
H = @sum J(i, j)S[i]S[j]

lattice = Square(4, 4)
H(lattice)

end

@testset "Heisenberg Model with given correlation" begin

lattice = Chain(4)
@nearrest i, j
A = size(lattice, i, j) |> zeros
A[1, 3] = 3
J = i, j -> A[i, j]

H = @sum J(i, j)Z[i]Z[j]
H(lattice)

end

@testset "Toric Code" begin

@surround v
@face p

# A_v = @prod X[v]
# B_p = @prod Z[p]

H = -J * @sum(@prod X[v]) - J * @sum(@prod B[p])

lattice = Square(8, 8)
H(lattice)

end

@testset "SYK Model" begin

# NOTE:
# url: https://nationalmaglab.org/images/news_events/searchable_docs/winterschool/2018/Theory_2018_sachdev_syk.pdf

@region :unknow i, j, k, l

U = @parameter i, j, k, l->randn()
H = 1/(2N^(3/2)) * @sum(U[i, j, k, l]C[i]C[j]C[k]C[l]) - μ @sum(C[i] * conj(C)[j])
H(N)

end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
using BinDeps
using Compat

@BinDeps.setup

# https://support.hdfgroup.org/HDF5/ lists "Current Releases"
# make sure we have one of those
const MINVERSION = v"1.8.0"
function h5_get_libversion(lib, handle)
    sym = Libdl.dlsym_e(handle, "H5get_libversion")
    majnum, minnum, relnum = Ref{Cuint}(), Ref{Cuint}(), Ref{Cuint}()
    status = ccall(sym, Cint, (Ptr{Cuint}, Ptr{Cuint}, Ptr{Cuint}), majnum, minnum, relnum)
    status < 0 && error("Error getting HDF5 library version")
    VersionNumber(majnum[], minnum[], relnum[])
end
compatible_version(lib, handle) = h5_get_libversion(lib, handle) >= MINVERSION

hdf5 = library_dependency("libhdf5",
    aliases = ["libhdf5", "libhdf5_serial", "libhdf5_serial.so.10", "libhdf5_openmpi", "libhdf5_mpich"],
    validate=compatible_version)

provides(AptGet, "hdf5-tools", hdf5, os=:Linux)
provides(Pacman, "hdf5", hdf5, os=:Linux)
provides(Yum, "hdf5", hdf5, os=:Linux)
provides(BSDPkg, "hdf5", hdf5, os=:FreeBSD)

if Compat.Sys.iswindows()
    using WinRPM
    provides(WinRPM.RPM, "hdf5", hdf5, os=:Windows)
end

if Compat.Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "hdf5", hdf5, os=:Darwin)
end

provides(Sources, URI("https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/src/hdf5-1.10.0-patch1.tar.gz"), hdf5)
provides(BuildProcess, Autotools(libtarget=joinpath("src", "libhdf5.la")), hdf5)

@BinDeps.install Dict(:libhdf5 => :libhdf5)


text:
import Observables
import Makie

q1 = Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
scene = Makie.Scene()
r = 1 + rand()
R = 10 + rand()
segments = rand(5:10)
color = Makie.RGBAf(rand(4)...)
transparency = false
torus = Torus(q1,
              scene,
              r = r,
              R = R,
              segments = segments,
              color = color,
              transparency = transparency)
q2 = q1 * Biquaternion(Quaternion(rand(4)), ℝ³(rand(3)))
value1 = getsurface(torus.observable, segments)
update(torus, q2)
value2 = getsurface(torus.observable, segments)

@test isapprox(torus.q, q2)
@test isapprox(value1, value2) == false

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module DefModelDatasets
using DataDeps

function __init__()
    register(DataDep(
        "Definition",
        """
        Dataset: WordNet
        Reference: https://arxiv.org/abs/1612.00394
        Source: https://github.com/northanapon/torch-defseq
        Website: https://github.com/DefinitionModeling/torch-defseq
        """,
        "https://github.com/DefinitionModeling/torch-defseq/archive/master.tar.gz",
        "94d4a4257038ef0c9865d7eb6564134629d115efaa80d01a19004d17f99eb038",
        post_fetch_method=(file -> unpack(file))
    ))

    register(DataDep(
        "Learning",
        """
        Datasets: Oxford, Urban Dictionary, Wikipedia, WordNet
        Reference: https://aclanthology.org/N19-1350/
        Source: https://github.com/shonosuke/ishiwatari-naacl2019
        Website: https://github.com/DefinitionModeling/ishiwatari-naacl2019
        """,
        "http://www.tkl.iis.u-tokyo.ac.jp/~ishiwatari/naacl_data.zip",
        "d8775ea04757c9281f13e9a45ef5234eec5f9eea91b5411ea0781a93e3a134b2",
        post_fetch_method=(file -> unpack(file))
    ))

    register(DataDep(
        "Evaluating",
        """
        Datasets: WordNet, Wiktionary, OmegaWiki
        Website: https://github.com/DefinitionModeling/Multi-sense-Multi-lingual-Definition-Modeling
        """,
        "https://github.com/DefinitionModeling/Multi-sense-Multi-lingual-Definition-Modeling/archive/master.tar.gz",
        "b95d1cd007eba7b4e734dc93692cab517a13ef86402dcb9bde9aefb76d041a98",
        post_fetch_method=(file -> unpack(file))
    ))

    register(DataDep(
        "Generationary",
        """
        Dataset: Hei++, SamplEval
        Reference: https://aclanthology.org/2020.emnlp-main.585/
        Source: https://sapienzanlp.github.io/generationary-web/
        """,
        "https://sapienzanlp.github.io/generationary-web/res/HEI++1.0_SamplEval.zip",
        "e5b1335d92bc9c4a3659d962147f72c547125dac7333b9266ffa07d9da24f450",
        post_fetch_method=(file -> unpack(file))
    ))
end

end


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

"""
```
call(gd::Guard, [lk::Link,] f::Function, args...)
```
Call a guard actor `gd` to execute `f(var, args...)` on
its guarded variable `var` and to respond with a
deep copy of the result.

# Arguments
- `gd::Guard`: link to the `:guard` actor,
- `lk::Link`: if a link `lk` is provided send the response
    to that, else do a synchronous `request`,
- `f`: function taking the guarded variable `var`
    as first argument,
- `args...`: further arguments to `f`.
"""
Actors.call(gd::Guard, lk::Link, f::F, args...) where F<:Function = call(gd.link, lk, f, args...)
Actors.call(gd::Guard, f::F, args...) where F<:Function = call(gd.link, f, args...)

"""
```
call(gd::Guard [, lk::Link])
```
Call a guard actor `gd` to respond with a deep copy of
the guarded variable `var`.

# Arguments
# Arguments
- `gd::Guard`: link to the `:guard` actor,
- `lk::Link`: if a link `lk` is provided send the response
    to that, else do a synchronous `request`,
"""
Actors.call(gd::Guard, lk::Link) = call(gd.link, lk)
Actors.call(gd::Guard) = call(gd.link)

"""
    cast(gd::Guard, f, args...)

Cast a message to a guard actor `gd` to execute `f(var, args...)` on
its guarded variable var.

# Arguments
- `gd::Guard`: link to the `:guard` actor,
- `f`: callable object taking the guarded variable `var`
    as first argument,
- `args...`: further arguments to `f`.
"""
Actors.cast(gd::Guard, f, args...) = cast(gd.link, f, args...)

"""
    update!(gd::Guard, var)

Update a guarded variable represented by `gd` with `var`.

**Note:** `var` must be of the same type as the guarded
variable!
"""
function Actors.update!(gd::Guard, var)
    @assert guardtype(gd) == typeof(var) "var must have the same type as the guarded variable"
    init!(gd.link, id, var)
    return deepcopy(var)
end

"""
    @grd f(gd, args...)

Execute a function `f` on a guarded variable `var`
represented by an actor link `gd` and return a deep copy
of the result or return a deep copy of `var`.  This is
a wrapper to [`call`](@ref).

    @grd gd

Return a deep copy of the guarded variable.

# Parameters
- `gd::Guard`: a link to the `:guard` actor,
- `f`: callable object taking the guarded variable `var`
    as first argument,
- `args...`:  further arguments to `f`.
"""
macro grd(expr)
    if expr isa Symbol
        esc(:(call($expr)))
    elseif expr isa Expr && expr.head == :call
        f = expr.args[1]
        gd = expr.args[2]
        args = expr.args[3:end]
        esc(:(call($gd, $f, $(args...))))
    else
        "@grd: not a symbol or a call!"
    end
end


text:
using Dates
using HTTP
using HDF5
using LinearAlgebra
using MPIFiles
using Statistics
using Test
using UUIDs
using Unitful

if !isdir("data")
  @info "download data.zip"
  HTTP.open("GET", "http://media.tuhh.de/ibi/MPIFiles/data.zip") do http
    open("data.zip", "w") do file
        write(file, http)
    end
  end
  @info "extracting data"
  run(`unzip -oq data.zip`)
  rm("data.zip")
end

include("Positions.jl")
include("General.jl")
include("MDFv1.jl")
include("MultiMPIFile.jl")
include("Reco.jl")
include("IMT.jl")
include("TransferFunction.jl")

@info "The unit tests are done!"


text:
module QuantilePriors

using ApproxFun
using Statistics

export QuantilePrior,
    mean,
    quantile

"""
    const S = Chebyshev(0..1)

All functions in this package use a Chebyshev basis on the [0, 1] interval.
"""
const S = Chebyshev(0..1)
"""
    const QP_INT = Integral(S)

The integral operator for the Chebyshev basis on the [0, 1] interval. Used for
calculating the mean of a prior.
"""
const QP_INT = Integral(S)

"""
    QuantilePrior{T} where T<:Fun
    QuantilePrior(p, q)

A quantile (inverse CDF) function that can be used as a prior, inspired by [Ben
Goodrich's StanCon 2020 talk](https://youtu.be/_wfZSvasLFk).

Usually specified via a vector of probabilities `p` and corresponding quantiles
`q`. These are interpolated using Chebyshev polynomials from the ApproxFun
package. This must result in a monotone function to represent a valid
distribution (this is tested when the function is constructed).

This would be used in e.g. MCMC settings by sampling from [0, 1] and then using
`quantile` to recover a parameter value that can be used. In this case the
density function of the prior is not needed; it is included implicitly in the
posterior density due to the change of variables.

Note that the resulting prior has a compact support, which may or may not be a
problem.
"""
struct QuantilePrior{T}
    f::T

    function QuantilePrior(f::T) where T <: Fun
        domain(f) == 0..1 || throw(DomainError("Domain must be [0, 1]"))
        ismonotoneinc(f) || throw("Not monotone; specify additional quantiles")
        new{T}(f)
    end
end

function QuantilePrior(p, q)
    all(0 .≤ p .≤ 1) || throw(DomainError("All `p` must be in [0, 1]"))
    V = vandermonde(p)
    f = Fun(S, V \ q)
    QuantilePrior(f)
end

"""
    vandermonde(p)

Construct a Vandermonde matrix at grid locations `p`, as described in ApproxFun
FAQ.
"""
function vandermonde(p)
    T = eltype(p)
    # Following ApproxFun FAQ
    n = length(p)
    V = Array{T}(undef, n, n)
    for k in 1:n
        V[:, k] = Fun(S, [zeros(T, k - 1); 1]).(p)
    end
    V
end

"""
    ismonotoneinc(f::Fun, x)

Test that the derivatives of `f` at `x` are all positive.
"""
function ismonotoneinc(f::Fun, x=range(0, 1, length=1025))
    df = f'.(x)
    all(df .> 0)
end

# Needed for broadcasting
Statistics.quantile(qp::QuantilePrior, p) = qp.f(p)
# Statistics.mean(qp::QuantilePrior) = (QP_INT * qp.f)(1)
Statistics.mean(qp::QuantilePrior) = stop("Not implemented (correctly)")

end


text:
using ArgoData, MeshArrays, Test

@testset "ArgoData.jl" begin

    files_list=GDAC.files_list()
    GDAC.download_file(files_list[10000,:])

    ftp="ftp://usgodae.org/pub/outgoing/argo/dac/"
    fil=GDAC.download_file(files_list[10000,:],"meta",ftp)

    @test isfile(fil)

    tmp=GDAC.grey_list()
    @test isa(tmp,GDAC.DataFrame)

    fil=joinpath(dirname(pathof(ArgoData)),"..","examples","ArgoToMITprof.jl")
    include(fil)

    @test isa(output_file,String)

    mp=MITprofStandard(output_file)

    MITprof.write(output_file*".tmp1",mp);
    MITprof.write(output_file*".tmp2",[mp,mp]);

    @test isa(mp,MITprofStandard)

    pth=dirname(output_file)
    nt,np,nz,cost=AnalysisMethods.cost_functions(pth,"prof_S")
    @test isapprox(cost[1],1.365848840650727)

    γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
    Γ=GridLoad(γ)
    df=AnalysisMethods.profile_positions(pth,Γ)
    @test isapprox(maximum(df.lat),-39.894)

end


text:
using CounterfactualExplanations
using CounterfactualExplanations.Models
using CounterfactualExplanations.Counterfactuals
using Random, LinearAlgebra
Random.seed!(1234)

@testset "Generic" begin

    # Data:
    using CounterfactualExplanations.Data
    using Random
    Random.seed!(1234)
    N = 25
    xs, ys = Data.toy_data_linear(N)
    X = hcat(xs...)
    counterfactual_data = CounterfactualData(X,ys')

    # Generator:
    generator = GenericGenerator()

    @testset "LogisticModel" begin
        # Model
        using CounterfactualExplanations.Models: LogisticModel, probs
        # Logit model:
        w = [1.0 1.0] # true coefficients
        b = 0
        M = LogisticModel(w, [b])

        # Randomly selected factual:
        Random.seed!(123)
        x = select_factual(counterfactual_data,rand(1:size(X)[2]))
        y = round(probs(M, x)[1])

        @testset "Predetermined outputs" begin
            γ = 0.9
            target = round(probs(M, x)[1])==0 ? 1 : 0
            counterfactual = generate_counterfactual(x, target, counterfactual_data, M, generator; γ=γ)
            @test counterfactual.target == target
            @test counterfactual.x == x
        end

        @testset "Convergence" begin

            # Already in target and exceeding threshold probability:
            γ = probs(M, x)[1]
            target = round(γ)
            counterfactual = generate_counterfactual(x, target, counterfactual_data, M, generator)
            @test length(path(counterfactual))==1
            @test counterfactual.x == counterfactual.x′
            @test converged(counterfactual) == true

            # Threshold reached if converged:
            γ = 0.9
            target = round(probs(M, x)[1])==0 ? 1 : 0
            T = 1000
            counterfactual = generate_counterfactual(x, target, counterfactual_data, M, generator; γ=γ, T=T)
            import CounterfactualExplanations.Counterfactuals: counterfactual_probability
            @test !converged(counterfactual) || counterfactual_probability(counterfactual)[1] >= γ # either not converged or threshold reached
            @test !converged(counterfactual) || length(path(counterfactual)) <= T

        end
    end

end


text:
### order matters!!!
# if Experiment.jl loads Neurons.jl the latter one has to be included before

include("MyUtils.jl")
include("Behavior.jl")
include("Neurons.jl")
include("Experiment.jl")
include("MyDSP.jl")
include("MyLinearAlgebra.jl")

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
using Random, Parameters

using FundStabABM.Types
using FundStabABM.Params
using FundStabABM.Func

params = Params.default()
@unpack bigt, bigm, bign, bigk, marketstartval = params
@testset "Market Type" begin

    market = Types.MarketIndex(vcat(
        marketstartval,
        zeros(bigt-1)))

    @test isa(market, Types.MarketIndex)

    @test typeof(market.value) == Array{Float64, 1}

    @test size(market.value) == (bigt, )

    @test market.value[1] == marketstartval
    @test all(market.value .>= 0)
end

@testset "Stocks Type" begin

    stocks = Types.Equity(
        zeros(bigm, bigt),
        zeros(bigm),
        zeros(bigm),
        zeros(bigm),
        zeros(bigm, bigt))

    @test isa(stocks, Types.Equity)

    @test typeof(stocks.value) == Array{Float64, 2}
    @test typeof(stocks.vol) == Array{Float64, 1}
    @test typeof(stocks.beta) == Array{Float64, 1}
    @test typeof(stocks.impact) == Array{Float64, 1}
    @test typeof(stocks.volume) == Array{Float64, 2}

    @test size(stocks.value) == (bigm, bigt)
    @test size(stocks.beta) == (bigm,)
    @test size(stocks.vol) == (bigm,)
    @test size(stocks.impact) == (bigm,)
    @test size(stocks.volume) == (bigm, bigt)

    # NOTE: No test of stock.beta now
    @test all(stocks.value .>= 0)
    @test all(stocks.vol .>= 0)
    @test all(stocks.impact .>= 0)
    @test all(stocks.volume .>= 0)
end

@testset "Fund Type" begin

    funds = Types.EquityFund(
    zeros(bigk, bigm),
    zeros(bigk, bign),
    zeros(bigk, bigt))

    @test isa(funds, Types.EquityFund)

    @test typeof(funds.holdings) == Array{Float64, 2}
    @test typeof(funds.stakes) == Array{Float64, 2}
    @test typeof(funds.value) == Array{Float64, 2}

    @test all(funds.holdings .>= 0)
    @test all(funds.stakes .>= 0)
    @test all(funds.value .>= 0)

    @test size(funds.holdings) == (bigk, bigm)
    @test size(funds.stakes) == (bigk, bign)
    @test size(funds.value) == (bigk, bigt)
end

@testset "Investor Type" begin

    investors = Types.RetailInvestor(
    zeros(bign, bigk + 1),
    zeros(bign),
    zeros(bign))

    @test isa(investors, Types.RetailInvestor)

    @test typeof(investors.assets) == Array{Float64, 2}
    @test typeof(investors.horizon) == Array{Int64, 1}
    @test typeof(investors.threshold) == Array{Float64, 1}

    @test size(investors.assets) == (bign, bigk + 1,)
    @test size(investors.horizon) == (bign,)
    @test size(investors.threshold) == (bign,)

    @test all(investors.assets .>= 0)
    @test all(investors.horizon .>= 0)
    @test all(1 .>= investors.threshold .>= -1)
end

@testset "Order Types" begin

    buyorder = Types.BuyMarketOrder(
    Array{Float64}(undef, 0, bigm),
    Array{Int64}(undef, 0))

    sellorder = Types.SellMarketOrder(
    Array{Float64}(undef, 0, bigm),
    Array{Int64}(undef, 0))

    @test isa(buyorder, Types.BuyMarketOrder)
    @test isa(sellorder, Types.SellMarketOrder)

    @test typeof(buyorder.values) == Array{Float64, 2}
    @test typeof(buyorder.funds) == Array{Int64, 1}
    @test typeof(sellorder.values) == Array{Float64, 2}
    @test typeof(sellorder.investors) == Array{Int64, 1}

    @test size(buyorder.values) == (0, bigm)
    @test size(buyorder.funds) == (0, )
    @test size(sellorder.values) == (0, bigm)
    @test size(sellorder.investors) == (0, )

    @test all(buyorder.values .>= 0)
    @test all(buyorder.funds .>= 0)
    @test all(sellorder.values .<= 0)
    @test all(sellorder.investors .>= 0)
end

#  Market makes perhaps not necessary
# @testset "Market Maker Type" begin
#
#     mktmaker = Types.MarketMaker(zeros(Params.bign, Params.bigm))
#
#     @test isa(mktmaker, Types.MarketMaker)
#
#     @test typeof(mktmaker.orderbook) == Array{Float64, 2}
#
#     @test size(mktmaker.orderbook) == (Params.bign, Params.bigm)
#
#end