text:
export getFEMMatrices2D, getFEMMatrices3D, getFEMsource3D, getFEMsource2D, getFEMrhs3D, getFEMrhs2D

function Phi(i::Int,x::Float64,y::Float64)
    if i==1
        return (1-x)*(1-y)
    elseif i==2
        return  x*(1-y)
    elseif i==3
        return  x*y
    elseif i==4
        return y*(1-x)
    end
end

function Phi(i::Int, x::Float64, y::Float64, z::Float64)
    if i==1
        return   (1-x)*(1-y)*(1-z)
    elseif i==2
        return   x*(1-y)*(1-z)
    elseif i==3
        return   x*y*(1-z)
    elseif i==4
        return   y*(1-x)*(1-z)
	elseif i==5
	    return   (1-x)*(1-y)*z
	elseif i==6
	    return   x*(1-y)*z
	elseif i==7
	    return   x*y*z
	elseif i==8
	    return   y*(1-x)*z
    end
end
function gradPhi(i::Int, x::Float64, y::Float64, z::Float64)
    if i==1
        return [-(1-y)*(1-z) -(1-x)*(1-z) -(1-x)*(1-y)]
    elseif i==2
        return [(1-y)*(1-z) -x*(1-z) -x*(1-y)]
    elseif i==3
        return [y*(1-z) x*(1-z) -x*y]
    elseif i==4
        return [-y*(1-z) (1-x)*(1-z) -y*(1-x)]
    elseif i==5
        return [-(1-y)*z -(1-x)*z  (1-x)*(1-y)]
    elseif i==6
        return [(1-y)*z -x*z  x*(1-y)]
    elseif i==7
        return [y*z x*z x*y]
    elseif i==8
        return [-y*z (1-x)*z y*(1-x)]
    end
end

function gradPhi(i::Int,x::Float64, y::Float64)
    if i==1
        return  [y-1 x-1]
    elseif i==2
        return  [1-y -x]
    elseif i==3
        return  [y x]
    elseif i==4
        return  [-y 1-x]
    end
end

function getFEMMatrices2D(M::RegularMesh;sig=0.01,v=[1.;0])
h      = M.h;
Jac    = Diagonal(h)
invJac = inv(Jac)
detJac = prod(h)

gaussRule=5
wGauss = [0.118463442528094; 0.239314335249683; 0.284444444444445;
0.239314335249683; 0.118463442528094];
xGauss = [0.9530899229693319; 0.7692346550528416; 0.5000000000000000;
0.2307653449471585; 0.0469100770306681];
xGauss = sort(xGauss)

X,Y = meshgrid(xGauss,xGauss);
W   = kron(wGauss,wGauss);
# wGauss*wGauss';
numquadnodes = gaussRule^M.dim;
X = vec(X)
Y = vec(Y)
W = vec(W)

# element stiffness matrix for Laplacian
sm_laplace = zeros(4,4)
for ii=1:4
    for jj=1:4
        for q=1:numquadnodes
            x = X[q];
            y = Y[q];
            sm_laplace[ii,jj] +=  dot(gradPhi(ii,x,y)*invJac,gradPhi(jj,x,y)*invJac)*W[q]*detJac;
        end
    end
end

# element mass matrix
mass = zeros(4,4);
for ii=1:4
    for jj=1:4
        for q=1:numquadnodes
            x = X[q];
            y = Y[q];
            mass[ii,jj] += Phi(ii,x,y)*Phi(jj,x,y)*W[q]*detJac
        end
    end
end

# element mass matrix with constant ansatz functions
mass_const = zeros(4);
for ii=1:4
    for q=1:numquadnodes
        x = X[q];
        y = Y[q];
        mass_const[ii] +=  Phi(ii,x,y)*1.0*W[q]*detJac
    end
end

# element stiffness matrix for advection
sm_advection = zeros(4,4)
for ii=1:4
    for jj=1:4
        for q=1:numquadnodes
            x = X[q]
            y = Y[q]
            sm_advection[ii,jj] +=  Phi(ii,x,y)*dot(v,gradPhi(jj,x,y)*invJac)*W[q]*detJac
        end
    end
end

nn = M.n.+1;

Ism = zeros(Int,0)
Jsm = zeros(Int,0)
Vsm = zeros(Float64,0)

Imass = zeros(Int,0)
Jmass = zeros(Int,0)
Vmass = zeros(Float64,0)

Imass_const = zeros(Int,0)
Jmass_const = zeros(Int,0)
Vmass_const = zeros(Float64,0)

for ii=1:nn[2]-1
    for jj=1:nn[1]-1
        nodes = [jj+nn[1]*(ii-1); jj+nn[1]*(ii-1)+1; jj+nn[1]*ii+1; jj+nn[1]*ii];
        for k=1:4
            for l=1:4
                push!(Imass,nodes[k])
				push!(Jmass,nodes[l])
				push!(Vmass,mass[k,l])

				push!(Ism, nodes[k])
				push!(Jsm, nodes[l])
				vsm = sig*sm_laplace[k,l] + sm_advection[k,l]
				push!(Vsm, vsm)
           end
			push!(Imass_const,nodes[k])
			push!(Jmass_const,(ii-1)*(nn[1]-1)+jj)
			push!(Vmass_const,mass_const[k])
        end
    end
end
SM         = sparse(Ism,Jsm,Vsm, prod(nn), prod(nn))
Mass       = sparse(Imass, Jmass, Vmass, prod(nn), prod(nn))
Mass_const = sparse(Imass_const, Jmass_const, Vmass_const, prod(nn), M.nc)
return Mass, Mass_const, SM
end

function getFEMMatrices3D(M::RegularMesh;sig=0.01,v=[1.;0;0])

h      = M.h;
Jac    = Diagonal(h)
invJac = inv(Jac)
detJac = prod(h)

gaussRule=5
wGauss = [0.118463442528094; 0.239314335249683; 0.284444444444445;
0.239314335249683; 0.118463442528094];
xGauss = [0.9530899229693319; 0.7692346550528416; 0.5000000000000000;
0.2307653449471585; 0.0469100770306681];

X,Y,Z = meshgrid(xGauss,xGauss,xGauss);
W     = kron(wGauss,kron(wGauss,wGauss));
# wGauss*wGauss';
numquadnodes = gaussRule^M.dim;
X = vec(X)
Y = vec(Y)
Z = vec(Z)
W = vec(W)

# element stiffness matrix for Laplacian
sm_laplace = zeros(8,8)
for ii=1:8
    for jj=1:8
        for q=1:numquadnodes
            x = X[q]
            y = Y[q]
			z = Z[q]
			sm_laplace[ii,jj] +=  dot(gradPhi(ii,x,y,z)*invJac,gradPhi(jj,x,y,z)*invJac)*W[q]*detJac;
        end
    end
end

# element mass matrix
mass = zeros(8,8);
for ii=1:8
    for jj=1:8
        for q=1:numquadnodes
            x = X[q]
            y = Y[q]
            z = Z[q]
			mass[ii,jj] += Phi(ii,x,y,z)*Phi(jj,x,y,z)*W[q]*detJac
        end
    end
end

# element mass matrix with constant ansatz functions
mass_const = zeros(8);
for ii=1:8
    for q=1:numquadnodes
        x = X[q]
        y = Y[q]
        z = Z[q]
		mass_const[ii] +=  Phi(ii,x,y,z)*1.0*W[q]*detJac
    end
end

# element stiffness matrix for advection for v = [1;0;0]
sm_advection = zeros(8,8)
for ii=1:8
    for jj=1:8
        for q=1:numquadnodes
            x = X[q]
            y = Y[q]
            z = Z[q]
			sm_advection[ii,jj] += Phi(ii,x,y,z)*dot(v,gradPhi(jj,x,y,z)*invJac)*W[q]*detJac
        end
    end
end

nn         = M.n.+1;

Ism = zeros(Int,0)
Jsm = zeros(Int,0)
Vsm = zeros(Float64,0)

Imass = zeros(Int,0)
Jmass = zeros(Int,0)
Vmass = zeros(Float64,0)

Imass_const = zeros(Int,0)
Jmass_const = zeros(Int,0)
Vmass_const = zeros(Float64,0)

nodes = zeros(Int,8)
incz = prod(nn[1:2])

for kk=1:nn[3]-1
	for ii=1:nn[2]-1
	    for jj=1:nn[1]-1
	 		nodes[1] =  jj+nn[1]*(ii-1)+incz*(kk-1);
	 		nodes[2] =  jj+nn[1]*(ii-1)+1+incz*(kk-1);
			nodes[3] =  jj+nn[1]*ii+1+incz*(kk-1);
			nodes[4] =  jj+nn[1]*ii+incz*(kk-1)
			nodes[5] =  jj+nn[1]*(ii-1)+incz*kk;
			nodes[6] =  jj+nn[1]*(ii-1)+1+incz*kk;
			nodes[7] =  jj+nn[1]*ii+1+incz*kk;
			nodes[8] =  jj+nn[1]*ii+incz*kk;
	       for k=1:8
	            for l=1:8
	                push!(Imass,nodes[k])
					push!(Jmass,nodes[l])
					push!(Vmass,mass[k,l])

					push!(Ism, nodes[k])
					push!(Jsm, nodes[l])
					vsm = sig*sm_laplace[k,l] + sm_advection[k,l]
					push!(Vsm, vsm)
	           end
			end
			for k=1:8
				push!(Imass_const,nodes[k])
				push!(Jmass_const,(ii-1)*(nn[1]-1)+jj+ prod(nn[1:2].-1)*(kk-1))
				push!(Vmass_const,mass_const[k])
	        end
	    end
	end
end
SM         = sparse(Ism,Jsm,Vsm, prod(nn), prod(nn))
Mass       = sparse(Imass, Jmass, Vmass, prod(nn), prod(nn))
Mass_const = sparse(Imass_const, Jmass_const, Vmass_const, prod(nn), M.nc)
return Mass, Mass_const, SM
end

function getFEMsource3D(M::RegularMesh)
	xc  = getCellCenteredGrid(M)
	blob1 = (xc[:,1].>.1) .& (xc[:,1] .< 2.2 ) .& (xc[:,2].>.3) .& (xc[:,2] .< 0.5) .& (xc[:,3] .> .3 ) .& (xc[:,3] .< 1.1);
	f     = 1.0*(blob1 )
	return f
end

function getFEMsource2D(M::RegularMesh)

	xc  = getCellCenteredGrid(M)
	blob1 = (xc[:,1].>.5) .& (xc[:,1] .< .7 ) .& (xc[:,2].>.3) .& (xc[:,2] .< 0.5)
	blob2 = (xc[:,1].>1.0) .& (xc[:,1] .< 1.5 ) .& (xc[:,2].>.6) .& (xc[:,2] .< 0.8)
	f     = 1.0*(blob1 + blob2)
	return f
end

function getFEMrhs2D(M::RegularMesh)
nn  = M.n
hh  = M.h
rhs = zeros(M.nc,1);
for ii=1:nn[2]
	for jj=1:nn[1]
		x      = (jj-0.5)*hh[1];
		y      = (ii-0.5)*hh[2];
		blob1  = (x > 0.5) & (x < 0.7 ) & (y > 0.3) & (y < 0.5)
		blob2  = (x > 1.0) & (x < 1.5 ) & (y > 0.6) & (y < 0.8)
		p      =  (ii-1)*nn[1]+jj
		rhs[p] = 1.0*(blob1 + blob2)

	end
end
return rhs
end

function getFEMrhs3D(M::RegularMesh)
nn   = M.n
hh   = M.h
rhs  = zeros(M.nc,1)
incz = prod(nn[1:2])
for kk=1:nn[3]
	for ii=1:nn[2]
		for jj=1:nn[1]
			x      = (jj-0.5)*hh[1];
			y      = (ii-0.5)*hh[2];
			z      = (kk-0.5)*hh[3];
			blob1  = (x >.1) & (x < 1 ) & (y > 0.3) & (y < 0.5) & (z > 0.3 ) & (z < .5);
			blob2  = (x > 1) & (x < 2 ) & (y > 0.6) & (y < 0.9) & (z > 0.6 ) & (z < .9);
			p      = (ii-1)*nn[1]+jj+(kk-1)*incz
			#rhs[p] = 1.0*(blob1)
			rhs[p] = 1.0*(blob1 + blob2)
		end
	end
end
return rhs
end

#
# xnodes = linspace(0,xlength, nn[1]);
# ynodes = linspace(0,ylength, nn[2]);
# [X,Y] = meshgrid(xnodes, ynodes);
#
# rhs = zeros((nn[1]-1)*(nn[2]-1),1);
# for ii=1:(nn[2]-1)
#     for jj=1:(nn[1]-1)
#         x = (jj+0.5)*hx;
#         y = (ii+0.5)*hy;
#         rhs((ii-1)*(nn[1]-1)+jj) = (x>0.5 && x<0.7 && y>0.3 && y<0.5) || (x>1.0 && x<1.5 && y>0.6 && y<0.8);
#     end
# end
# V = Mass_const*rhs;
# u = SM\V;                         % ... solve SM*u = Mass_const*rhs, where whs are element controls
#                                   % ... u has Nx*nn[2] entries U(i,j) = u((i-1)*Nx+j) for i=1,Nx, j=1,nn[2]
#                                   % ... rhs has (Nx-1)*(Ny-1) entries W(i,j) = rhs((i-1)*(Nx-1)+j) for i=1,Nx-1, j=1,Ny-1
# PrintAMPL = 0;
# if PrintAMPL
#   outFid = fopen('FEM.dat','w');
#   for k=1:Nx*Ny
#     % SM(k,:)*u = Mass_const(k,:)*rhs
#     irows = find(SM(k,:));
#     fprintf(outFid,'FEM_%-6i:',k);
#     for l=1:length(irows)
#       jj = mod(irows(l),Nx) + 1;
#       ii = floor(irows(l) / Nx) + 1;
#       aa = full(SM(k,irows(l)));
#       fprintf(outFid,' + %12.8g * U[%6i,%6i]', aa,ii,jj);
#     end
#     fprintf(outFid,';\n');
#   end
#   fclose(outFid);
# end
#
# u_plot = reshape(u, Nx, Ny)';
# figure(1); clf;
# imagesc(reshape(rhs,Nx-1,Ny-1))
# figure(2); clf;
# surf(X,Y,u_plot)
# %daspect("auto")
# pbaspect([xlength ylength max(u)-min(u)])

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
$ julia -e 'for x in ARGS; println(x); end' foo bar
foo
bar

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Autogenerated wrapper script for SBC_jll for aarch64-linux-gnu
export libsbc, sbcdec, sbcenc, sbcinfo

using libsndfile_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libsbc`
const libsbc_splitpath = ["lib", "libsbc.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libsbc_path = ""

# libsbc-specific global declaration
# This will be filled out by __init__()
libsbc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libsbc = "libsbc.so.1"

# Relative path to `sbcdec`
const sbcdec_splitpath = ["bin", "sbcdec"]

# This will be filled out by __init__() for all products, as it must be done at runtime
sbcdec_path = ""

# sbcdec-specific global declaration
function sbcdec(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(sbcdec_path)
    end
end

# Relative path to `sbcenc`
const sbcenc_splitpath = ["bin", "sbcenc"]

# This will be filled out by __init__() for all products, as it must be done at runtime
sbcenc_path = ""

# sbcenc-specific global declaration
function sbcenc(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(sbcenc_path)
    end
end

# Relative path to `sbcinfo`
const sbcinfo_splitpath = ["bin", "sbcinfo"]

# This will be filled out by __init__() for all products, as it must be done at runtime
sbcinfo_path = ""

# sbcinfo-specific global declaration
function sbcinfo(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(sbcinfo_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"SBC")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (libsndfile_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (libsndfile_jll.LIBPATH_list,))

    global libsbc_path = normpath(joinpath(artifact_dir, libsbc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libsbc_handle = dlopen(libsbc_path)
    push!(LIBPATH_list, dirname(libsbc_path))

    global sbcdec_path = normpath(joinpath(artifact_dir, sbcdec_splitpath...))

    push!(PATH_list, dirname(sbcdec_path))
    global sbcenc_path = normpath(joinpath(artifact_dir, sbcenc_splitpath...))

    push!(PATH_list, dirname(sbcenc_path))
    global sbcinfo_path = normpath(joinpath(artifact_dir, sbcinfo_splitpath...))

    push!(PATH_list, dirname(sbcinfo_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
using Test
using ColorBlendModes
using ColorTypes, FixedPointNumbers

if VERSION >= v"1.6.0-DEV.816" # JuliaLang/julia PR #36962
    @test isempty(detect_ambiguities(ColorBlendModes))
else
    @test isempty(detect_ambiguities(ColorBlendModes, Base, Core))
end

const blend_modes = map(k -> parse(BlendMode, k), ColorBlendModes.mode_keywords)
const separable_modes = filter(m -> m isa ColorBlendModes.SeparableBlendMode, blend_modes)

const Ts = (Float64, Float32, N0f8, N0f16)

# RGB -> Gray conversions are implemented in Colors.jl, and the following are
# slightly different from the implementations of Colors.jl.
function ColorTypes._convert(::Type{Cout}, ::Type{C1}, ::Type{C2},
                             c) where {Cout<:AbstractGray,C1<:AbstractGray,C2<:AbstractRGB}
    Cout(ColorBlendModes.lum100(c) / 100)
end
function ColorTypes._convert(::Type{A}, ::Type{C1}, ::Type{C2},
                             c, alpha=alpha(c)) where {A<:TransparentGray,C1<:AbstractGray,C2<:AbstractRGB}
    A(ColorBlendModes.lum100(color(c)) / 100, alpha)
end

const gray_n0f8 = Gray{N0f8}[0, 0.004, 0.251, 0.498, 0.502, 0.506, 0.753, 0.996, 1]
const rgb1_n0f8 = RGB{N0f8}.(gray_n0f8)
const rgb2_n0f8 = [RGB{N0f8}(gray_n0f8[3i+mod1(i+1,3)], gray_n0f8[3i+mod1(i+2,3)], gray_n0f8[3i+mod1(i+3,3)]) for i = 0:2]

const gray_n0f16 = Gray{N0f16}[0, 2e-5, 0.25, 0.49999, 0.50001, 0.50002, 0.75001, 0.99998, 1]
const rgb1_n0f16 = RGB{N0f16}.(gray_n0f16)
const rgb2_n0f16 = [RGB{N0f16}(gray_n0f16[3i+mod1(i+1,3)], gray_n0f16[3i+mod1(i+2,3)], gray_n0f16[3i+mod1(i+3,3)]) for i = 0:2]

const gray_f32 = Gray{Float32}[0.0f0, Float32(0x1p-126), 0.25f0,
                               prevfloat(0.5f0), 0.5f0, nextfloat(0.5f0),
                               0.75f0, prevfloat(1.0f0), 1.0f0]
const rgb1_f32 = RGB{Float32}.(gray_f32)
const rgb2_f32 = [RGB{Float32}(gray_f32[3i+mod1(i+1,3)], gray_f32[3i+mod1(i+2,3)], gray_f32[3i+mod1(i+3,3)]) for i = 0:2]

const gray_f64 = Gray{Float64}[0.0, 0x1p-1022, 0.25,
                               prevfloat(0.5), 0.5, nextfloat(0.5),
                               0.75, prevfloat(1.0), 1.0]
const rgb1_f64 = RGB{Float64}.(gray_f64)
const rgb2_f64 = [RGB{Float64}(gray_f64[3i+mod1(i+1,3)], gray_f64[3i+mod1(i+2,3)], gray_f64[3i+mod1(i+3,3)]) for i = 0:2]

f64(c::C) where C<:Colorant = convert(base_colorant_type(C){Float64}, c)
f32(c::C) where C<:Colorant = convert(base_colorant_type(C){Float32}, c)
n0f8(c::C) where C<:Colorant = convert(base_colorant_type(C){N0f8}, c)
n0f16(c::C) where C<:Colorant = convert(base_colorant_type(C){N0f16}, c)

@testset "low level optimization" begin
    r_n0f8 = 0N0f8:eps(N0f8):1N0f8
    r_n0f16 = 0N0f16:eps(N0f16):1N0f16
    @testset "mul" begin
        function test_mul(v1::N, v2::N) where {T, N <: Normed{T}}
            rv1, rv2 = reinterpret(v1), reinterpret(v2)
            expected = reinterpret(N, round(T, rv1 / typemax(T) * rv2))
            ColorBlendModes.mul(v1, v2) === expected
        end
        @test all(test_mul(v1, v2) for v1 in r_n0f8,  v2 in r_n0f8)
        @test all(test_mul(v1, v2) for v1 in r_n0f16, v2 in r_n0f16)
    end

    @testset "div01" begin
        function test_div01(v1::N, v2::N) where {T, N <: Normed{T}}
            v2 == zero(N) && return true
            rv1, rv2 = reinterpret(v1), reinterpret(v2)
            expected = reinterpret(N, round(T, min(1.0, rv1 / rv2) * typemax(T)))
            ColorBlendModes.div01(v1, v2) === expected
        end
        @test all(test_div01(v1, v2) for v1 in r_n0f8,  v2 in r_n0f8)
    end
end

@testset "hue operation" begin
    Hue = ColorBlendModes.Hue

    @test ColorBlendModes._w(Hue(30.0), Hue(130.0), 0.0) ≈ 30.0
    @test ColorBlendModes._w(Hue(30.0), Hue(130.0), 0.6) ≈ 90.0
    @test ColorBlendModes._w(Hue(30.0), Hue(130.0), 1.0) ≈ 130.0
    @test ColorBlendModes._w(Hue(130.0), Hue(30.0), 0.6) ≈ 70.0
    @test ColorBlendModes._w(Hue(175.0), Hue(185.0), 0.6) ≈ 181.0
    @test ColorBlendModes._w(Hue(230.0), Hue(30.0), 0.6) ≈ 326.0
    @test ColorBlendModes._w(Hue(5.0), Hue(355.0), 0.6) ≈ 359.0

    @test ColorBlendModes._w(Hue(30.0), 1.0, Hue(130.0), 0.0) ≈ 30.0
    @test ColorBlendModes._w(Hue(30.0), 0.4, Hue(130.0), 0.6) ≈ 90.0
    @test ColorBlendModes._w(Hue(30.0), 0.0, Hue(130.0), 1.0) ≈ 130.0
    @test ColorBlendModes._w(Hue(130.0), 0.4, Hue(30.0), 0.6) ≈ 70.0
    @test ColorBlendModes._w(Hue(175.0), 0.4, Hue(185.0), 0.6) ≈ 181.0
    @test ColorBlendModes._w(Hue(230.0), 0.4, Hue(30.0), 0.6) ≈ 326.0
    @test ColorBlendModes._w(Hue(5.0), 0.4, Hue(355.0), 0.6) ≈ 359.0

    @test ColorBlendModes._w(Hue(30.0, true), Hue(130.0), 0.0) ≈ 130.0
    @test ColorBlendModes._w(Hue(30.0, true), Hue(130.0), 0.6) ≈ 130.0
    @test ColorBlendModes._w(Hue(30.0, true), Hue(130.0), 1.0) ≈ 130.0
    @test ColorBlendModes._w(Hue(30.0), Hue(130.0, true), 0.0) ≈ 30.0
    @test ColorBlendModes._w(Hue(30.0), Hue(130.0, true), 0.6) ≈ 30.0
    @test ColorBlendModes._w(Hue(30.0), Hue(130.0, true), 1.0) ≈ 30.0
end

@testset "blend: RGB" begin
    include("blend_rgb.jl")
end

@testset "blend: Gray (and RGB)" begin
    include("blend_gray.jl")
end

@testset "blend: Lab/Luv" begin
    include("blend_lab.jl")
end

@testset "blend: HSV/HSL/HSI" begin
    include("blend_hsx.jl")
end

@testset "composite" begin
    include("composite.jl")
end

@testset "combination" begin
    include("combination.jl")
end

@testset "keyword" begin
    @test keyword(BlendNormal) == "normal"
    @test keyword(BlendColorDodge) == "color-dodge"

    @test keyword(CompositeSourceOver) == "source-over"
end

@testset "parse" begin
    @test parse(BlendMode, "color-burn") === BlendColorBurn
    @test parse(BlendMode, "Hard-Light") === BlendHardLight
    @test_throws ArgumentError parse(BlendMode, "SoftLight")

    @test parse(CompositeOperation, "source-over") === CompositeSourceOver
    @test parse(CompositeOperation, "Source-Atop") === CompositeSourceAtop
    @test_throws ArgumentError parse(CompositeOperation, "SourceAtop")
end


text:
module MultiFieldTests

using Test

@testset "MultiFieldCellFields" begin include("MultiFieldCellFieldsTests.jl") end

@testset "MultiFieldFESpaces" begin include("MultiFieldFESpacesTests.jl") end

@testset "MultiFieldFEFunctions" begin include("MultiFieldFEFunctionsTests.jl") end

@testset "MultiFieldSparseMatrixAssemblers" begin include("MultiFieldSparseMatrixAssemblersTests.jl") end

@testset "MultiFieldFEOperators" begin include("MultiFieldFEOperatorsTests.jl") end

@testset "MultiFieldFESpacesWithLinearConstraints" begin include("MultiFieldFESpacesWithLinearConstraintsTests.jl") end

end # module


text:
using BED
using Distributions
using GenomicFeatures
using GenomePermutations
using HypothesisTests
using Statistics
using Test

@testset "GenomicPermutations.jl" begin

    # set up intervals to test
    interval1 = GenomicFeatures.Interval("chr1", 1150, 1180)  # test2
    interval2 = GenomicFeatures.Interval("chr1", 250,280) # test3
    interval3 = GenomicFeatures.Interval("chr1", 99,280) # test7
    interval4 = GenomicFeatures.Interval("chr3", 99,280) # test7
    interval5 = GenomicFeatures.Interval("chr2", 2150, 2180)
    interval6 = GenomicFeatures.Interval("chr2", 2350, 2380)
    collection1 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 100),
        GenomicFeatures.Interval("chr1", 200, 300),
        GenomicFeatures.Interval("chr1", 400, 500)
        ]) # test1
    collection2 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 50, 80),
        GenomicFeatures.Interval("chr1", 250, 280)
        ]) # test4
    collection3 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 50, 80),
        GenomicFeatures.Interval("chr1", 2250, 2280)
        ]) # test5
    collection4 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 2150, 2180),
        GenomicFeatures.Interval("chr1", 2250, 2280)
        ]) # test6
    collection5 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 2150, 2180),
        GenomicFeatures.Interval("chr1", 2250, 2280),
        GenomicFeatures.Interval("chr2", 2250, 2280)
        ])
    collection6 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr2", 2150, 2180),
        GenomicFeatures.Interval("chr2", 2250, 2280)
    ])
    collection7 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 0, 1000),
        GenomicFeatures.Interval("chr1", 2000, 3000)
    ])
    collection8 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr2", 2150, 2180),
        GenomicFeatures.Interval("chr2", 2250, 2280),
        interval6
    ])
    collection9 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 100),
        GenomicFeatures.Interval("chr1", 200, 300),
        GenomicFeatures.Interval("chr1", 2400, 2500),
    ])
    collection10 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 100),
        GenomicFeatures.Interval("chr1", 200, 300),
        GenomicFeatures.Interval("chr1", 2400, 2500),
        GenomicFeatures.Interval("chr2", 1, 100),
        GenomicFeatures.Interval("chr2", 200, 300),
        GenomicFeatures.Interval("chr2", 2400, 2500)
    ])
    collection11 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 5),
        GenomicFeatures.Interval("chr1", 21, 25),
        GenomicFeatures.Interval("chr1", 41, 45)
    ])
    collection12 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 5),
        GenomicFeatures.Interval("chr1", 21, 25),
        GenomicFeatures.Interval("chr1", 41, 45),
        GenomicFeatures.Interval("chr2", 1, 5),
        GenomicFeatures.Interval("chr2", 21, 25),
        GenomicFeatures.Interval("chr2", 41, 45),
    ])
    collection13 = GenomicFeatures.IntervalCollection([
            GenomicFeatures.Interval("chr1", 101, 199),
            GenomicFeatures.Interval("chr1", 501, 600),
            GenomicFeatures.Interval("chr1", 601, 701)
    ])
    collection14 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 100),
        GenomicFeatures.Interval("chr1", 200, 300),
        GenomicFeatures.Interval("chr1", 400, 500),
        GenomicFeatures.Interval("chr2", 1, 100),
        GenomicFeatures.Interval("chr2", 200, 300),
        GenomicFeatures.Interval("chr2", 400, 500)
    ])
    collection15 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 50, 80),
        GenomicFeatures.Interval("chr1", 250, 280),
        GenomicFeatures.Interval("chr2", 50, 80),
        GenomicFeatures.Interval("chr2", 250, 280)
    ])
    collection16 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 101, 199),
        GenomicFeatures.Interval("chr1", 501, 600),
        GenomicFeatures.Interval("chr1", 601, 701),
        GenomicFeatures.Interval("chr2", 101, 199),
        GenomicFeatures.Interval("chr2", 501, 600),
        GenomicFeatures.Interval("chr2", 601, 701)
    ])
    collection17 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr2", 50, 80),
        GenomicFeatures.Interval("chr2", 250, 780)
    ])

    regions1 = collection7
    regions2 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1", 0, 1000),
        GenomicFeatures.Interval("chr1", 2000, 3000),
        GenomicFeatures.Interval("chr2", 0, 1000),
        GenomicFeatures.Interval("chr2", 2000, 3000)
    ])
    regions3 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1",500,650),
        GenomicFeatures.Interval("chr1", 1000, 1101)
    ])
    regions4 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chr1",500,650),
        GenomicFeatures.Interval("chr1", 1000, 1101),
        GenomicFeatures.Interval("chr2",500,650),
        GenomicFeatures.Interval("chr2", 1000, 1101)
    ])
    regions5 = GenomicFeatures.IntervalCollection([
        GenomicFeatures.Interval("chrM", 0, 1000),
        GenomicFeatures.Interval("chrM", 2000, 3000),
    ])
    regions6 = IntervalCollection([GenomicFeatures.Interval("chr1", 1, 500000)])
    regions7 = IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 500000),
        GenomicFeatures.Interval("chr2", 1, 500000)
        ])
    regions8 = IntervalCollection([GenomicFeatures.Interval("chrM", 1, 500000)])
    regions9 = IntervalCollection([GenomicFeatures.Interval("chr1", 1, 701)])
    regions10 = IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 701),
        GenomicFeatures.Interval("chr2", 1, 701)
    ])
    regions11 = IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 701),
        GenomicFeatures.Interval("chrM", 1, 701)
    ])
    regions12 =  IntervalCollection([
        GenomicFeatures.Interval("chr1", 1, 500000),
        GenomicFeatures.Interval("chr2", 1, 500000),
        GenomicFeatures.Interval("chrM", 1, 500000)
        ])

    # no need to test isoverlapping(::GenomicFeatures.Interval, ::GenomicFeatures.Interval)
    @testset "overlaps and isin functions" begin
        @test ([false, true, true, false] == map(x->anyoverlapping(x,collection1), [interval1, interval2, interval3, interval4]))
        @test GenomePermutations.alloverlapping(collection1, collection1)
        @test ([true, false, false]  == map(x->alloverlapping(x,collection1), [collection2, collection3, collection4]))

        @test isin(interval1,interval1) && isin(interval2,interval3) && !(isin(interval1,interval2)) && !(isin(interval3,interval4))
        @test isin(interval2, collection1) && !(isin(interval1, collection1))
        @test isin(collection1, collection1) && isin(collection2, collection1) && !(isin(collection3, collection2)) && !(isin(collection6, collection5))

        # NEED TO ADD A TEST FOR WHEN INTERVALS ARE contained in a collection --> need to differentiatie in/issubset and contain
        # test subsetting functions are equivalent
        @test collection4 == GenomePermutations.vec_getcollection(collection5, "chr1") == GenomePermutations.iter_getcollection(collection5, "chr1")

        # test overlap counts
        @test countoverlapping(collection1, collection1) == collection1.length
        @test countoverlapping(collection2, collection1) == 2
        @test countoverlapping(collection6, collection1) == countoverlapping(collection1, collection6) == 0
        @test countoverlapping(collection5, collection6) == 1
    end

    @testset "distance functions" begin
        # for single intervals
        @test ismissing(dist(interval3, interval4)) && dist(interval1, interval2) == 870 && dist(interval2, interval3) == 0

        # for an interval and a collection
        @test [20,150,96.66666666666667,120] == map( x->dist(interval2 ,collection1, x), [minimum, maximum, mean, median])
        @test [970, 1070, 1020.0, 1020.0 ] == map( x->dist(interval1 ,collection5, x), [minimum, maximum, mean, median])
        #@test_logs (:warn, r"could not calculcate distance") match_mode=:any map( x->dist(interval1 ,collection5, x), [minimum, maximum, mean, median])

        # for 2 collections
        @test dist(collection5, collection5) == [0, 0, 0]
        @test  [[0, 170], 0, 170, 85.0, 85.0] == map( x->dist(collection2 ,collection3, x), [x -> x, minimum, maximum, mean, median])
        @test  [[2050, 1850, 1650], 1650, 2050, 1850.0, 1850.0] == map( x->dist(collection1 ,collection5, x), [x -> x, minimum, maximum, mean, median])
        #@test_logs (:warn, r"could not calculcate distance") match_mode=:any map( x->dist(collection1 ,collection5, x), [x -> x, minimum, maximum, mean, median])
        @test [[1650, 1750], 1650, 1750, 1700.0, 1700.0] ==  map( x->dist(collection5 ,collection1, x), [x -> x, minimum, maximum, mean, median]) # why does this fail
    end

    # test distribution generation
    @testset "generate distributions" begin
        @test generatedistribution(collection1) ==
            MixtureModel([DiscreteUniform(1,100), DiscreteUniform(200, 300), DiscreteUniform(400, 500)], Distributions.Categorical([99/299, 100/299, 100/299]))
        @test typeof(generatedistribution(collection1)) <: Distributions.Distribution
    end

    dist2 = generatedistribution(collection2)
    dist3 = generatedistribution(collection3)
    dist6 = generatedistribution(collection6)
    dist7 = generatedistribution(collection7)
    dist8 = generatedistribution(collection8)

    @testset "private interval randomisation" begin
        @test_throws ErrorException GenomePermutations._randominterval(interval1, dist7, GenomicFeatures.IntervalCollection{nothing}(); allow_overlap = true)
        @test isin(GenomePermutations._randominterval(interval2, dist7, collection7; allow_overlap = true), collection7)
        @test isin(GenomePermutations._randominterval(interval2, dist7, collection7; allow_overlap = true), collection7)
        @test_throws ErrorException GenomePermutations._randominterval(interval5, dist6, collection6; collection = collection6, allow_overlap = false)  # must throw exception
        @test isin(GenomePermutations._randominterval(interval5, dist6, collection8; collection = collection6, allow_overlap = true), collection8)
        @test GenomePermutations._randominterval(interval5, dist8, collection8; collection = collection6, allow_overlap = false) == interval6 # can only be 1 place
        @test_throws ErrorException GenomePermutations._randominterval(interval4, dist6, collection1; allow_overlap = false) # all stupidly mismatched. Indeed it fails gracefully
        @testset "fail behaviour" begin
            # test default behaviour
            @test_throws ErrorException  GenomePermutations._randominterval(interval3, dist2, collection2; allow_overlap = false)
            # test :orig
            @test GenomePermutations._randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :orig) == interval3
            # test :skip
            @test isnothing(GenomePermutations._randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :skip))
            # test :throw
            @test_throws ErrorException   GenomePermutations._randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :throw)
        end
    end

    @testset "public interval randomisation" begin
        # now we get prettier errors if sequences mismatch
        @test_throws KeyError GenomePermutations.randominterval(interval1, dist7, GenomicFeatures.IntervalCollection{nothing}(); allow_overlap = true)
        @test_nowarn GenomePermutations.randominterval(interval1, dist7, collection7; allow_overlap = true)
        @test isin(GenomePermutations.randominterval(interval2, dist7, collection7; allow_overlap = true, onfail = :orig), collection7)
        @testset "public fail behaviour" begin
            # test default
            @test_throws ErrorException  GenomePermutations.randominterval(interval3, dist2, collection2; allow_overlap = false)
            # test :orig
            @test GenomePermutations.randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :orig) == interval3
            # test skip
            @test isnothing(GenomePermutations.randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :skip))
            # test :throw
            @test_throws ErrorException   GenomePermutations.randominterval(interval3, dist2, collection2; allow_overlap = false, onfail = :throw)
        end
    end

    @testset "region randomisation" begin
        result = randomiseregions(collection1, dist7, collection7)
        @test isin(result, collection7) && result != collection1 # check perfect case
        # check we fail if we must overlap
        @test_throws ErrorException randomiseregions(collection8, dist6, collection6)
        # check if allowing overlaps rescues test 2
        @test_nowarn randomiseregions(collection8, dist6, collection6; allow_overlap = true)
        # check we converge to the only available solution
        @test randomiseregions(collection2, dist3, collection3) == collection3
        # check we converge to self
        @test randomiseregions(collection2, dist2, collection2; allow_overlap = false, max_tries=100000) == collection2
        #= This mess should erorr in an undefinend way. God only knows what kind of exception
        but it should do any error and refuse to go through an infinite loop =#
        @test_throws ErrorException randomiseregions(collection3, dist2, collection6)
        @testset "check error handling" begin
            @test_throws ErrorException randomiseregions(collection17, dist6, collection6)
            @test_throws ErrorException randomiseregions(collection17, dist6, collection6, onfail = :throw)
            @test_throws ErrorException randomiseregions(collection17, dist6, collection6)
            @test randomiseregions(collection17, dist6, collection6; onfail = :skip).length == 1
        end
    end

    @testset "genome randomisation" begin
        # test 1 sequence
        @test isin(randomisegenome(collection9, regions1), regions1)
        @test randomisegenome(collection9, regions1) != collection9
        # test 2 sequences
        @test isin(randomisegenome(collection10, regions2), regions2)
        @test randomisegenome(collection10, regions2) != collection10
        #= test over itself - note we need tiny interval to make sure we
        always draw the only possible solution within a reasonable number of tries =#
        @test randomisegenome(collection11, collection11; max_tries = 10000) == collection11
        @test randomisegenome(collection12, collection12; max_tries = 10000) == collection12
        # test a case that only works with no overlaps
        @test_throws Core.Exception randomisegenome(collection9, regions3)
        @test isin(randomisegenome(collection9, regions3; allow_overlap=true), regions3)
        @test_throws Core.Exception randomisegenome(collection10, regions4)
        @test isin(randomisegenome(collection10, regions4; allow_overlap=true), regions4)
        # test an impossible case
        @test_logs (:warn, r"is not found in the target regions")  randomisegenome(collection10, regions1) # skips wrong sequences
        # all sequences are wrong - should return empty
        @test randomisegenome(collection9, regions5) == GenomicFeatures.IntervalCollection{Nothing}()

        @testset "onfail behaviour"  begin
            # test we successfully throw if we fail to randomise an interval
            @test_throws ErrorException randomisegenome(collection9, regions3; onfail = :throw) == collection9
            # test a case where we pass the original interval on fail
            @testset ":orig" begin
                scollection1 = GenomicFeatures.IntervalCollection([
                    GenomicFeatures.Interval("chr1", 1, 2),
                    GenomicFeatures.Interval("chr1",3, 4),
                    GenomicFeatures.Interval("chr1", 200, 500),
                    ])
                sregions1 = GenomicFeatures.IntervalCollection([
                    GenomicFeatures.Interval("chr1", 1, 4),
                    GenomicFeatures.Interval("chr1", 300, 400),
                    ])
                # test a case where we pass the original interval on fail and allow overlaps
                res = randomisegenome(scollection1, sregions1; onfail = :orig, allow_overlap = true)
                @test isin(GenomicFeatures.Interval("chr1", 200, 500), res)
                # really odd corner case
                @test_broken randomisegenome(scollection1, sregions1; onfail = :orig, allow_overlap = false) == scollection1
            end
        end
    end

    @testset "simple P structure (constructor)" begin
        sp = SimplePTest(100, 0.1, :Less)
        @test sp.iterations == 100
        @test sp.p_val == 0.1
        @test sp.alternative == :Less
        @test sp.min_p_val == 0.009900990099009901
        @test pvalue(sp) == 0.1
    end

    @testset "PermTestResult structure (constructor)" begin
        @testset "single value" begin
            res = PermTestResult(3, 10, [10, 10, 10], SimplePTest(3, 1, :Less) , "test1", "test2", overlappermtest ,"guess and check")
            @test typeof(res) == PermTestResult{Int, Int, SimplePTest}
            @test length(res.ran) == res.iterations
            @test pvalue(res) == 1.0
        end
        @testset "array" begin
            res = GenomePermutations.PermTestResult(3, [1, 2, 3], [[1, 2, 3],[1, 2, 3],[1, 2, 3]], SimplePTest(3, 1, :Less) , "test1", "test2", overlappermtest ,"guess and check")
            @test typeof(res) == PermTestResult{Vector{Int}, Vector{Int}, SimplePTest}
            @test length(res.ran) == res.iterations
            @test pvalue(res) == 1.0
        end
    end

    @testset "simple P test" begin
        @test GenomePermutations.simpleP(80, [80, 80, 80, 80, 80], 5; alternative = "Auto") == GenomePermutations.SimplePTest(5, 1, :Less)
        @test GenomePermutations.simpleP(80, [400, 400, 400, 400, 400, 400, 400, 400, 400, 400], 10; alternative = "Auto") == GenomePermutations.SimplePTest(10, 0.09090909090909091, :Less)
        @test GenomePermutations.simpleP(80, [78, 78, 78, 99, 78], 5; alternative = "Auto") == GenomePermutations.simpleP(80, [2, 2, 2, 99, 2], 5; alternative = "Auto") ==GenomePermutations.SimplePTest(5, 0.3333333333333333,:More)
    end

    @testset "overlap permutation test" begin
        @testset "single sequence" begin
            # 1 Overlap with itself should be as significant as it can
            @test overlappermtest(collection1,  collection1, regions6, 100).test.p_val ==
                GenomePermutations.simpleP(1,repeat([1], 100),100).min_p_val
            # 2 Less overlaps than expected
            test = overlappermtest(collection1, collection13, regions9, 100)
            @test test.test.p_val < 0.01 && test.test.alternative == :Less
            # 3 More overlaps than expected
            test = overlappermtest(collection1, collection2, regions6, 100)
            @test test.test.p_val < 0.01 && test.test.alternative == :More
            # 4 sequenecs to randomise do not belong to regions
            @test_throws ErrorException overlappermtest(collection1,  collection2, regions8, 100)
            # 5 Regions to randomise have an extra sequence
            @test_nowarn overlappermtest(collection1,  collection2, regions11, 100)
        end
        @testset "multiple sequences" begin
            # 1 Overlap with itself should be as significant as it can
            @test overlappermtest(collection10,  collection10, regions7, 100).test.p_val ==
                GenomePermutations.simpleP(1,repeat([1], 100),100).min_p_val
            # 2 Less overlaps than expected
            test = overlappermtest(collection14,  collection16, regions10, 100)
            @test test.test.p_val < 0.01 && test.test.alternative == :Less
            # 3 More overlaps than expected
            test = overlappermtest(collection14,  collection15, regions7, 100)
            @test test.test.p_val < 0.01 && test.test.alternative == :More
            # 4 All sequences wrong
            @test_throws ErrorException overlappermtest(collection14,  collection15, regions8, 100)
            # 5 One of 2 sequences are wrong
            @test_throws ErrorException overlappermtest(collection14,  collection15, regions11, 100)
            # 6 there is an extra chromosome - This should work fine
            @test_nowarn overlappermtest(collection14,  collection15, regions12, 100)
        end
    end

    @testset "Random Genome Generation" begin
        @test_nowarn GenomePermutations._randomgenome(regions2, 10, Distributions.Categorical(2), Distributions.Binomial(10),["chr1","chr2"])
        @test isin(GenomePermutations._randomgenome(regions2, 10, Distributions.Categorical(2), Distributions.Binomial(10),["chr1","chr2"]), regions2)
        @test_throws KeyError GenomePermutations._randomgenome(regions1, 10, Distributions.Categorical(2), Distributions.Binomial(10),["chr1","chr2","chrZ"])
        @test_throws KeyError GenomePermutations._randomgenome(regions1, 10, Distributions.Categorical(1), Distributions.Binomial(10),["chr1","chr2","chrZ"])
    end

    @testset "Generic permutation test" begin
        # define some supplementary regions as needed
        scollection1 =  GenomicFeatures.IntervalCollection([
            GenomicFeatures.Interval("chr1", 2, 101),
            GenomicFeatures.Interval("chr1", 201, 301),
            GenomicFeatures.Interval("chr1", 402, 503)])
        scollection2 =  GenomicFeatures.IntervalCollection([
            GenomicFeatures.Interval("chr1", 2, 101),
            GenomicFeatures.Interval("chr1", 201, 301),
            GenomicFeatures.Interval("chr1", 402, 503),
            GenomicFeatures.Interval("chr2", 2, 101),
            GenomicFeatures.Interval("chr2", 201, 301),
            GenomicFeatures.Interval("chr2", 402, 503)
            ])
        scollection3 =  GenomicFeatures.IntervalCollection([
            GenomicFeatures.Interval("chr1",  499400, 499500),
            GenomicFeatures.Interval("chr1",  499700, 499800),
            GenomicFeatures.Interval("chr1", 499900, 500000)])
        scollection4 =  GenomicFeatures.IntervalCollection([
            GenomicFeatures.Interval("chr1",  499400, 499500),
            GenomicFeatures.Interval("chr1",  499700, 499800),
            GenomicFeatures.Interval("chr1", 499900, 500000),
            GenomicFeatures.Interval("chr2",  499400, 499500),
            GenomicFeatures.Interval("chr2",  499700, 499800),
            GenomicFeatures.Interval("chr2", 499900, 500000)])
        scollection5 = GenomicFeatures.IntervalCollection([
                GenomicFeatures.Interval("chr1", 1, 100),
                GenomicFeatures.Interval("chr1", 200, 300),
                GenomicFeatures.Interval("chr1", 400, 500),
                GenomicFeatures.Interval("chr2", 1, 100),
                GenomicFeatures.Interval("chr2", 200, 300),
                GenomicFeatures.Interval("chr2", 400, 500)
                ])

        @testset "default distance" begin
            @testset "Single sequence" begin
                # 1 distance to iteslf should be as significant as it gets
                test =  permtest(collection1, collection1, regions6, 100)
                @test pvalue(test.test) < 0.01
                # 2 closer than expected
                test = permtest(collection1, scollection1, regions6, 100)
                @test pvalue(test.test) < 0.01
                # 3 more distant than expected
                test = permtest(collection1, scollection3, regions6, 100)
                @test pvalue(test.test) < 0.01
                # 4 sequenecs to randomise do not belong to regions
                @test_throws ErrorException permtest(collection1,  collection2, regions8, 100)
                # 5 Regions to randomise have an extra sequence
                @test_nowarn permtest(collection1,  collection2, regions11, 100)
            end
            @testset "Mutliple sequences" begin
                # 1 distance to iteslf should be as significant as it gets
                test =  permtest(collection10, collection10, regions7, 100)
                @test pvalue(test.test) < 0.01
                # 2 closer than expected
                test = permtest(collection10, scollection2, regions7, 100)
                @test pvalue(test.test) < 0.01
               # 3 more distant than expected
                test = permtest(scollection5, scollection4, regions7, 100)
                @test pvalue(test.test) < 0.01
                # Need a test for the case where we expect a p of 1.
                # 4 All sequences wrong
                @test_throws ErrorException permtest(collection14,  collection15, regions8, 100)
                # 5 One of 2 sequences are wrong
                @test_throws ErrorException permtest(collection14,  collection15, regions11, 100)
                # 6 there is an extra chromosome - This should work fine
                @test_nowarn permtest(collection14,  collection15, regions12, 100)
            end

        end

        @testset "Large Collection test" begin
            lr1 = IntervalCollection(sort!([Interval(i, 1, 100000) for i in "chr" .* map(string, collect(1:22))]))
            lr2 = IntervalCollection(sort!([Interval(i, 1, 100000) for i in "chr" .* map(string, collect(1:22))]))
            lc1 = GenomePermutations._randomgenome(lr1, 400, Distributions.Categorical(22), Distributions.Binomial(500),"chr" .* map(string, collect(1:22)))
            lc2 = GenomePermutations._randomgenome(lr2, 400, Distributions.Categorical(22), Distributions.Binomial(500),"chr" .* map(string, collect(1:22)))
            @test isin(lc1, lr1) &&  isin(lc2, lr2)
            @test_nowarn permtest(lc1,  lc2, lr1, 1000)
        end
    end

end


text:
using Meshes
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
isappveyor = "APPVEYOR" ∈ keys(ENV)
isCI = istravis || isappveyor
datadir = joinpath(@__DIR__,"data")

# list of tests
testfiles = [
  "points.jl",
  "angles.jl",
  "geometries.jl",
  "polytopes.jl",
  "primitives.jl",
  "mesh.jl",
  "sampling.jl",
  "discretization.jl",
  "boundingboxes.jl"
]

# --------------------------------
# RUN TESTS WITH SINGLE PRECISION
# --------------------------------
T = Float32
P2, P3 = Point{2,T}, Point{3,T}
V2, V3 = Vec{2,T}, Vec{3,T}
@testset "Meshes.jl ($T)" begin
  for testfile in testfiles
    include(testfile)
  end
end

# --------------------------------
# RUN TESTS WITH DOUBLE PRECISION
# --------------------------------
T = Float64
P2, P3 = Point{2,T}, Point{3,T}
V2, V3 = Vec{2,T}, Vec{3,T}
@testset "Meshes.jl ($T)" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
import Mongoc
using Dates
using NamedTupleTools
using DataStructures

"""
Get the types from a union type.
https://stackoverflow.com/a/43305886/1133481
"""
union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type) = (x,)

struct TypeValueMismatch <: Exception
    msg::AbstractString
end

struct FieldsMismatch <: Exception
    msg::AbstractString
end

field_names_as_string(type) = map(fnSymbol -> string(fnSymbol), fieldnames(type))

"""
Convert the given bson document or value to the specified type, recursively converting
arrays, tuples, dicts and sub-documents according to the same structure produced by
`Mongoc.BSON(value::Any)`
"""
function as_type(type, value)
    valueType = typeof(value)

    # println("as_type   $(type)   $(value)   valueType: $(valueType)")

    if valueType <: type
        return value
    end

    # Attempt to `convert` the value.
    try
        converted = convert(type, value)
        # println("converted: $(converted)")
        return converted
    catch
        # println("could not convert")
    end

    if type <: Tuple
        vals = []
        for (index, val) in enumerate(value)
            elementType = eltype(type)
            convertedVal = as_type(elementType, val)
            push!(vals, (index, convertedVal))
        end
        sort!(vals, by = x -> x[1])
        vals = map(x -> x[2], vals)
        return (vals...,)
    elseif type == DataType
        return eval(Meta.parse(value))
    elseif type <: Dict
        kt = keytype(type)
        vt = valtype(type)

        dkeys = []
        for k in keys(value)
            push!(dkeys, as_type(kt, k))
        end

        dvalues = []
        for v in Base.values(value)
            push!(dvalues, as_type(vt, v))
        end

        return Dict(zip(dkeys, dvalues))
    elseif type <: NamedTuple
        symbols = fieldnames(type)
        vals = []
        for (s, t) in zip(symbols, type.types)
            push!(vals, as_type(t, value[string(s)]))
        end
        return namedtuple(symbols, vals)
    elseif type <: Number && valueType <: AbstractString
        return parse(type, value)
    elseif typeof(type) <: Union
        # Try to convert the value to each type in the union type.
        for ut in union_types(type)
            try
                return as_type(ut, value)
            catch
            end
        end
        throw(TypeValueMismatch("could not convert $(value) to any of the types in $(type)"))
    end

    typeFieldNames = field_names_as_string(type)

    if valueType <: AbstractDict
        valueFieldNames = [ keys(value)... ]
        # println("keys from dict/bson: $(valueFieldNames)")
    else
        valueFieldNames = field_names_as_string(valueType)
        # println("keys from other: $(valueFieldNames)")
    end

    if isempty(typeFieldNames) || isempty(valueFieldNames)
        throw(TypeValueMismatch("invalid value for type $(type): $(value)"))
    end

    if Set(typeFieldNames) != Set(valueFieldNames)
        throw(FieldsMismatch("value and type have mismatched fields: type fields = $(typeFieldNames); value fields = $(valueFieldNames)"))
    end

    values = []
    for (fname, ftype) in zip(typeFieldNames, type.types)
        push!(values, as_type(ftype, value[fname]))
    end
    return type(values...)
end

max_stars_repo_name: OliverColeman/BSONify

max_forks_repo_name: OliverColeman/BSONify


text:
export LPLoss, LPPenalty,
        L1Loss, L1Penalty,
        L2Loss, L2Penalty,
        LogisticLoss, MultinomialLoss

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# No Loss / No Penalty
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

(l::NoLoss)(a::AVR, b::AVR) = 0.0
(p::NoPenalty)(θ::AVR)      = 0.0

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# LP Losses and Penalties
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

"""
$TYPEDEF

Scaled L-p loss of the residual.

``L(x, y) = ||x-y||_{p}^{p}/p``

The scaling simplifies expressions in the common L2 case.
"""
struct LPLoss{p} <: AtomicLoss where p <: Real end

"""
$TYPEDEF

Scaled L-p norm of the parameter vector.

``P(θ) = ||θ||_{p}^{p}/p``

The scaling simplifies expressions in the common L2 case.
"""
struct LPPenalty{p} <: AtomicPenalty where p <: Real end

# Useful Shortcuts
const L1Loss    = LPLoss{1}
const L1Penalty = LPPenalty{1}
const L2Loss    = LPLoss{2}
const L2Penalty = LPPenalty{2}
const LPCost{p} = Union{LPLoss{p},LPPenalty{p}}

const L1R = ScaledPenalty{L1Penalty}
const L2R = Union{NoPenalty,ScaledPenalty{L2Penalty}}
const ENR = Union{L1R,CompositePenalty}

"""
$SIGNATURES

Return the `p` in an `LPCost{p}`.
"""
getp(lpc::LPCost{p}) where p = p

"""
$SIGNATURES

Compute the lp norm to the `p`-th power of a vector given `p` scaled by `p`.
"""
function lp(v::AbstractVector{<:Real}, p)
    p == Inf && return maximum(v)
    p == 1   && return sum(abs.(v))
    p == 2   && return sum(abs2.(v)) / 2
    p  > 0   && return sum(abs.(v).^p) / p
    throw(DomainError("[lp] `p` has to be greater than 0"))
end

(l::LPLoss)(a::AVR, b::AVR) = lp(a .- b, getp(l))
(l::LPLoss)(r::AVR)         = lp(r, getp(l))
(p::LPPenalty)(θ::AVR)      = lp(θ, getp(p))

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Logistic loss
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

"""
$TYPEDEF

``L(x, y) = -∑logσ(xᵢyᵢ)``

where `logσ` is the log of the sigmoid function; `yᵢ ∈ {±1}`. In a logistic regression `x`
corresponds to `Xθ` where `X` is the design matrix and `θ` the vector of parameters.
See [`logsigmoid`](@ref).
"""
struct LogisticLoss <: AtomicLoss end

(::LogisticLoss)(x::AVR, y::AVR) = -sum(logsigmoid.(x .* y))

"""
$TYPEDEF

``L(P, y) = ∑log Zᵢ - ∑∑ 1(yᵢ=j)Pᵢⱼ``

where `P` is a matrix where each row contains class probabilities, `yᵢ ∈ {1, 2, ..., K}`
corresponding to column indices and,

``Zᵢ = ∑ exp(Pᵢ)``

In a multinomial regression, `P` corresponds to `XW` where `X` is the design matrix and `W` the
matrix of size `p * K` where each column corresponds to the parameters corresponding to that class.
"""
struct MultinomialLoss <: AtomicLoss end

(::MultinomialLoss)(P::Matrix{<:Real}, y::Vector{Int}) = begin
    L = 0.0
    @inbounds for i in eachindex(y)
        Pᵢ = P[i, :]
        m  = maximum(Pᵢ)
        sᵢ = sum(exp.(Pᵢ .- m)) # avoid overflow
        logZᵢ = log(sᵢ) + m
        L += logZᵢ - P[i, y[i]]
    end
    return L
end
(l::MultinomialLoss)(y::Vector{Int}, P::Matrix{<:Real}) = l(P, y)


text:
using AutoGrad: Param, @diff, value, grad
using Knet: relu, sigm, tanh, softmax

in_size       = 12
hidden_size   = 20
out_size      = 12

memory_size   = 128
location_size = 32

hm_readers    = 2
hm_writers    = 1

Layer_Type    = :FeedForward

normalize(arr) = # why cant u normalize a 2d array LinearAlgebra
begin
    norm  = sqrt(sum(arr .^2))
    norm != 0 ? arr./norm : arr
end

mutable struct Recurrent
    state::Param
    wf1::Param
    wf2::Param
    bf::Param
    wk1::Param
    wk2::Param
    bk::Param
    wi::Param
    bi::Param

Recurrent(in_size,layer_size) = new(
    Param(zeros(1,layer_size)),
    Param(randn(in_size,layer_size)),
    Param(randn(layer_size,layer_size)),
    Param(zeros(1,layer_size)),
    Param(randn(in_size,layer_size)),
    Param(randn(layer_size,layer_size)),
    Param(zeros(1,layer_size)),
    Param(randn(in_size,layer_size)),
    Param(zeros(1,layer_size)),
)
end

(layer::Recurrent)(in) =
begin
    focus  = sigm.(in * layer.wf1 + layer.state * layer.wf2 + layer.bf)
    keep   = sigm.(in * layer.wk1 + layer.state * layer.wk2 + layer.bk)
    interm = tanh.(in * layer.wi  + layer.state .* focus    + layer.bi)

layer.state = Param(keep .* interm + (1 .- keep) .* layer.state)
end

mutable struct FeedForward
    w::Param
    b::Param

FeedForward(in_size, layer_size) = new(
    Param(randn(in_size, layer_size)),
    Param(zeros(1, layer_size)),
)
end

(layer::FeedForward)(in) =
begin

tanh.(in * layer.w + layer.b)
end

Layer = eval(Layer_Type)

mutable struct Processor
    layer_in::Layer
    layer_out::Layer
    read_key::Layer
    write_key::Layer
    memory_creator::Layer

Processor(in_size, hidden_size, out_size, location_size, hm_readers) = new(
    Layer(in_size+location_size*hm_readers, hidden_size),
    Layer(hidden_size, out_size),
    Layer(hidden_size, location_size),
    Layer(hidden_size, location_size),
    Layer(hidden_size, location_size),
)
end

(processor::Processor)(input, memory_attended) =
begin
    hidden = processor.layer_in(hcat(input, memory_attended))
    output = processor.layer_out(hidden)
    read_key = processor.read_key(hidden)
    write_key = processor.write_key(hidden)
    new_data = processor.memory_creator(hidden)

output, read_key, write_key, new_data
end

mutable struct Reader
    location_attend::Layer

Reader(location_size) = new(
    Layer(location_size*2, 1),
)
end

(reader::Reader)(read_key, memory) =
begin
    memory_attentions = softmax([reader.location_attend(hcat(read_key, location))[end] for location in memory])
    memory_attended = sum([location .* attention for (location, attention) in zip(memory, memory_attentions)])

memory_attended
end

mutable struct Writer
    data_interpret::Layer
    location_attend::Layer
    location_free::Layer
    location_alloc::Layer

Writer(location_size) = new(
    Layer(location_size*2, location_size),
    Layer(location_size*2, 1),
    Layer(location_size*2, location_size),
    Layer(location_size*2, location_size),
)
end

(writer::Writer)(memory, write_key, new_data) =
begin
    new_content = writer.data_interpret(hcat(write_key, new_data))
    memory_attentions = softmax([writer.location_attend(hcat(write_key, location))[end] for location in memory])
    free_attentions = [writer.location_free(hcat(write_key, location)) for location in memory]
    alloc_attentions = [writer.location_alloc(hcat(write_key, location)) for location in memory]
    new_memory = [(1 .- attention) .* location + (attention) .* (location .* free + new_content .* alloc) for (location, attention, free, alloc) in zip(memory, memory_attentions, free_attentions, alloc_attentions)]
    # [location .* (1 .- location_attention) .* (1 .- free_attention) .+ new_data .* location_attention .* alloc_attention for (location, attention, free, alloc) in zip(memory, memory_attentions, free_attentions, alloc_attentions)]

new_memory
end

mutable struct Model
    processor::Processor
    readers::Array{Reader}
    writers::Array{Writer}

Model(in_size, hidden_size, out_size, memory_size, location_size, hm_readers, hm_writers) = new(
    Processor(in_size, hidden_size, out_size, location_size, hm_readers),
    [Reader(location_size) for _ in 1:hm_readers],
    [Writer(location_size) for _ in 1:hm_writers],
)
end

(model::Model)(input, memory, read_key) =
begin
    memory_attended = hcat([reader(read_key, memory) for reader in model.readers]...)
    output, read_key, write_key, new_data = model.processor(input, memory_attended)
    new_memory = sum([writer(memory, write_key, new_data) for writer in model.writers]) ./ hm_writers

output, new_memory, read_key
end

propogate_timeseries(sequence, model; memory=nothing, read_key=nothing) =
begin
    memory == nothing ? memory = [randn(1, location_size) for _ in 1:memory_size] : ()
    read_key == nothing ? read_key = randn(1, location_size) : ()

    response = []
    for timestep in sequence
        output, memory, read_key = model(timestep, memory, read_key)
        push!(response, output)
    end

    if Layer_Type == "Recurrent"
        for component_name in fieldnames(typeof(model))
            component = getfield(model, component_name)
            for layer_name in fieldnames(typeof(component))
                layer = getfield(component, layer_name)
                layer.state = Param(zeros(size(layer.state)))
            end
        end
    end

response
end

##TESTS

model = Model(in_size, hidden_size, out_size, memory_size, location_size, hm_readers, hm_writers)
memory = [randn(1, location_size) for _ in 1:memory_size]

sample_input = randn(1, in_size)
sample_output = randn(1, out_size)

output, new_memory = model(sample_input, memory, randn(1,location_size))

##TESTS2

sample_timeserie = [randn(1, in_size) for _ in 1:5]

response = propogate_timeseries(sample_timeserie, model, memory=memory)


text:
using PackageCompiler

try
    using NSGAII
catch
    include("install.jl")
end
#using Pkg
#Pkg.clone("https://github.com/JuliaCN/Py2Jl.jl")

#using Py2Jl

include("local_iz_neuron.jl")
using NSGAII
include("../src/units.jl")

using .HHNSGA

function fdecode(x)
   decoded = NSGAII.decode(x, HHNSGA.bc)
   return decoded
end

using Plots
unicodeplots()

function plot_pop(P)
    P = filter(x -> x.rank == 1, P)
    plot(map(x -> x.y[1], P), map(x -> x.y[2], P))#, "bo", markersize = 1)
end
#@bp

#HHNSGA.init_function()
#@bp
@time repop = NSGAII.nsga(4,4,HHNSGA.z,HHNSGA.init_function)#, fplot = plot_pop)
@time best = sort(repop, by = ind -> ind.y[1])[end];
@time worst = sort(repop, by = ind -> ind.y[1])[1];

@time decoded = fdecode(best.pheno)
println("z = $(x1.y)")
py"""
from izhi import IZModel
model = IZModel()
"""
py"model.attrs" = decoded

py"""
rt = test['Rheobase test']
rheo = rt.generate_prediction(model)
"""

vm = py"model.get_membrane_potential()"
Pkg.add("Plotly")
using Plotly

trace1 = [
  "x" => vm,
  "y" => vm.times,
  "mode" => "lines+markers",
  "type" => "scatter"
]
data = [trace1]
response = Plotly.plot(data, ["filename" => "line-scatter", "fileopt" => "overwrite"])
plot_url = response["url"]

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
using Clang
using Test

include("cbasic.jl")
include("functions.jl")
include("wrap_c.jl")
include("wrap_struct_typedef.jl")


text:
export DummySurveillanceUnit, DummySurveillanceUnitParams

Base.@kwdef struct DummySurveillanceUnitParams <: DeviceParams

end
DummySurveillanceUnitParams(dict::Dict) = params_from_dict(DummySurveillanceUnitParams, dict)

Base.@kwdef mutable struct DummySurveillanceUnit <: SurveillanceUnit
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::DummySurveillanceUnitParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  acPowerEnabled::Bool = false
end

function init(su::DummySurveillanceUnit)
  @debug "Initializing dummy surveillance unit with ID `$(su.deviceID)`."

  su.present = true
end

neededDependencies(::DummySurveillanceUnit) = []
optionalDependencies(::DummySurveillanceUnit) = []

Base.close(su::DummySurveillanceUnit) = nothing

getTemperatures(su::DummySurveillanceUnit) = 30.0.*ones(4) .+ 1.0.*randn(4)
getACStatus(su::DummySurveillanceUnit) = su.acPowerEnabled

function enableACPower(su::DummySurveillanceUnit)
  @debug "Enable AC Power"
  su.acPowerEnabled = true
end

function disableACPower(su::DummySurveillanceUnit)
  @debug "Disable AC Power"
  su.acPowerEnabled = false
end

resetDAQ(su::DummySurveillanceUnit) = nothing
hasResetDAQ(su::DummySurveillanceUnit) = false


text:
# This file was generated, do not modify it. # hide
using Soss
using MeasureTheory
using TupleVectors

m = @model n begin
    α ~ Normal()
    β ~ Normal()
    x ~ Normal() |> iid(n)
    σ ~ Exponential(λ=1)
    y ~ For(x) do xj
        Normal(α + β * xj, σ)
    end
    return y
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
#=
Author: Andrew Ning

A general blade element momentum (BEM) method for propellers/fans and turbines.

Some unique features:
- a simple yet very robust solution method
- allows for non-ideal conditions like reversed flow, or free rotation
- allows arbitrary inflow
- convenience methods for common wind turbine inflow scenarios

=#

module CCBlade

import Dierckx  # cubic b-spline for airfoil data
import Roots  # solve residual equation
import Parameters: @unpack

export Section, Rotor, Inflow, Outputs
export af_from_file, af_from_data
export simpleinflow, windturbineinflow, windturbineinflow_az
export solve, loads, effectivewake, thrusttorque, thrusttorque_azavg, nondim

# --------- structs -------------

"""
    Section(r, chord, theta, af)

Define section properties at a given radial location on the rotor

**Arguments**
- `r::Float64`: radial location (Rhub < r < Rtip)
- `chord::Float64`: local chord length
- `theta::Float64`: twist angle (radians)
- `af`: a function of the form: cl, cd = af(alpha, Re, Mach)
"""
struct Section{TF, TAF}

    r::TF
    chord::TF
    theta::TF  # includes pitch
    af::TAF

end

"""
    Rotor(Rhub, Rtip, B, turbine, precone=0.0)

Define geometry common to the entire rotor.

**Arguments**
- `Rhub::Float64`: hub radius (along blade length)
- `Rtip::Float64`: tip radius (along blade length)
- `B::Int64`: number of blades
- `turbine::Bool`: true if turbine, false if propeller
- `precone::Float64`: precone angle
"""
struct Rotor{TF, TI, TB}

    Rhub::TF
    Rtip::TF
    B::TI
    turbine::TB
    precone::TF

end

# convenience constructor for no precone
Rotor(Rhub, Rtip, B, turbine) = Rotor(Rhub, Rtip, B, turbine, 0.0)

# make rotor broadcastable as a single entity
Base.Broadcast.broadcastable(r::Rotor) = Ref(r)

"""
    Inflow(Vx, Vy, rho, mu=1.0, asound=1.0)

Operation point for a rotor.  The x direction is the axial direction, and y direction is the tangential direction in the rotor plane.  See Documentation for more detail on coordinate systems.

**Arguments**
- `Vx::Float64`: velocity in x-direction
- `Vy::Float64`: velocity in y-direction
- `rho::Float64`: fluid density
- `mu::Float64`: fluid dynamic viscosity (unused if Re not included in airfoil data)
- `asound::Float64`: fluid speed of sound (unused if Mach not included in airfoil data)
"""
struct Inflow{TF}
    Vx::TF
    Vy::TF
    rho::TF
    mu::TF
    asound::TF
end

# convenience constructor when Re and Mach are not used.
Inflow(Vx, Vy, rho) = Inflow(Vx, Vy, rho, 1.0, 1.0)

"""
    Outputs(Np, Tp, a, ap, u, v, phi, W, cl, cd, F)

Outputs from the BEM solver.

**Arguments**
- `Np::Float64`: normal force per unit length
- `Tp::Float64`: tangential force per unit length
- `a::Float64`: axial induction factor
- `ap::Float64`: tangential induction factor
- `u::Float64`: axial induced velocity
- `v::Float64`: tangential induced velocity
- `phi::Float64`: inflow angle
- `W::Float64`: inflow velocity
- `cl::Float64`: lift coefficient
- `cd::Float64`: drag coefficient
- `F::Float64`: hub/tip loss correction
"""
struct Outputs{TF}

    Np::TF
    Tp::TF
    a::TF
    ap::TF
    u::TF
    v::TF
    phi::TF
    W::TF
    cl::TF
    cd::TF
    F::TF

end

# constructor for case with no solution found
Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

# -------------------------------

# ----------- airfoil ---------------

"""
    af_from_file(filename)

Read an airfoil file.
Currently only reads one Reynolds number.
Additional data like cm is optional but will be ignored.
alpha should be in degrees

format:

header\n
alpha1 cl1 cd1\n
alpha2 cl2 cd2\n
alpha3 cl3 cd3\n
...

Returns a function of the form `cl, cd = func(alpha, Re, M)` although Re and M are currently ignored.
"""
function af_from_file(filename)

    alpha = Float64[]
    cl = Float64[]
    cd = Float64[]

    open(filename) do f

        # skip header
        readline(f)

        for line in eachline(f)
            parts = split(line)
            push!(alpha, parse(Float64, parts[1]))
            push!(cl, parse(Float64, parts[2]))
            push!(cd, parse(Float64, parts[3]))
        end

    end

    return af_from_data(alpha*pi/180.0, cl, cd)
end

"""
    af_from_data(alpha, cl, cd)

Create an AirfoilData object directly from alpha, cl, and cd arrays.
alpha should be in radians.

`af_from_file` calls this function indirectly.  Uses a cubic B-spline
(if the order of the data permits it).  A small amount of smoothing of
lift and drag coefficients is also applied to aid performance
for gradient-based optimization.

Returns a function of the form `cl, cd = func(alpha, Re, M)` although Re and M are currently ignored.
"""
function af_from_data(alpha, cl, cd, spl_k=3)

    # # TODO: update once smoothing is implemented: https://github.com/JuliaMath/Interpolations.jl/issues/254
    # afcl = CubicSplineInterpolation(alpha, cl)
    # afcd = CubicSplineInterpolation(alpha, cd)
    # af = AirfoilData(afcl, afcd)

    k = min(length(alpha)-1, spl_k)  # can't use cubic spline is number of entries in alpha is small

    # 1D interpolations for now.  ignoring Re dependence (which is very minor)
    afcl_1d = Dierckx.Spline1D(alpha, cl; k=k, s=0.1)
    afcd_1d = Dierckx.Spline1D(alpha, cd; k=k, s=0.001)

    afeval(alpha, Re, M) = afcl_1d(alpha), afcd_1d(alpha)  # ignore Re, M

    return afeval
end

# ---------------------------------

# ------------ BEM core ------------------

"""
(private) residual function
"""
function residual(phi, section, inflow, rotor)

    # unpack inputs
    @unpack r, chord, theta, af = section
    @unpack Vx, Vy, rho, mu, asound = inflow
    @unpack Rhub, Rtip, B, turbine = rotor

    # check if turbine or propeller and change input sign if necessary
    swapsign = turbine ? 1 : -1
    theta *= swapsign
    Vx *= swapsign

    # constants
    sigma_p = B*chord/(2.0*pi*r)
    sphi = sin(phi)
    cphi = cos(phi)

    # angle of attack
    alpha = phi - theta

    # Reynolds number
    W0 = sqrt(Vx^2 + Vy^2)  # ignoring induction, which is generally a very minor error and only affects Reynolds number
    Re = rho * W0 * chord / mu

    # Mach number
    Mach = W0/asound  # also ignoring induction

    # airfoil cl/cd
    cl, cd = af(alpha, Re, Mach)

    # resolve into normal and tangential forces
    cn = cl*cphi + cd*sphi
    ct = cl*sphi - cd*cphi

    # Prandtl's tip and hub loss factor
    factortip = B/2.0*(Rtip - r)/(r*abs(sphi))
    Ftip = 2.0/pi*acos(exp(-factortip))
    factorhub = B/2.0*(r - Rhub)/(Rhub*abs(sphi))
    Fhub = 2.0/pi*acos(exp(-factorhub))
    F = Ftip * Fhub

    # sec parameters
    k = cn*sigma_p/(4.0*F*sphi*sphi)
    kp = ct*sigma_p/(4.0*F*sphi*cphi)

    # parameters used in Vx=0 and Vy=0 cases
    k0 = cn*sigma_p/(4.0*F*sphi*cphi)
    k0p = ct*sigma_p/(4.0*F*sphi*sphi)

    # --- solve for induced velocities ------
    if isapprox(Vx, 0.0, atol=1e-6)

        u = sign(phi)*k0*Vy
        v = 0.0
        a = 0.0
        ap = 0.0

    elseif isapprox(Vy, 0.0, atol=1e-6)

        u = 0.0
        v = k0p*abs(Vx)
        a = 0.0
        ap = 0.0

    else

        if phi < 0
            k = -k
        end

        if isapprox(k, -1.0, atol=1e-6)  # state corresopnds to Vx=0, return any nonzero residual
            return 1.0, Outputs()
        end

        if k <= 2.0/3  # momentum region
            a = k/(1 + k)

        else  # empirical region
            g1 = 2.0*F*k - (10.0/9-F)
            g2 = 2.0*F*k - (4.0/3-F)*F
            g3 = 2.0*F*k - (25.0/9-2*F)

            if isapprox(g3, 0.0, atol=1e-6)  # avoid singularity
                a = 1.0 - 1.0/(2.0*sqrt(g2))
            else
                a = (g1 - sqrt(g2)) / g3
            end
        end

        u = a * Vx

        # -------- tangential induction ----------
        if Vx < 0
            kp = -kp
        end

        if isapprox(kp, 1.0, atol=1e-6)  # state corresopnds to Vy=0, return any nonzero residual
            return 1.0, Outputs()
        end

        ap = kp/(1 - kp)
        v = ap * Vy

    end

    # ------- residual function -------------
    # R = sin(phi)*(Vy + v) - cos(phi)*(Vx - u)
    # R = sin(phi)/(1 - a) - Vx/Vy*cos(phi)/(1 + ap)
    R = sin(phi)/(Vx - u) - cos(phi)/(Vy + v)

    # ------- loads ---------
    W2 = (Vx - u)^2 + (Vy + v)^2
    Np = cn*0.5*rho*W2*chord
    Tp = ct*0.5*rho*W2*chord

    # ---- swap output signs as needed -----
    Tp *= swapsign
    v *= swapsign

    return R, Outputs(Np, Tp, a, ap, u, v, phi, sqrt(W2), cl, cd, F)  # multiply by F because a and ap as currently used are only correct in combination with the loads.  If you want a wake model then you need to add the hub/tip loss factors separately.

end

"""
(private) Find a bracket for the root closest to xmin by subdividing
interval (xmin, xmax) into n intervals.

Returns found, xl, xu.
If found = true a bracket was found between (xl, xu)
"""
function firstbracket(f, xmin, xmax, n, backwardsearch=false)

    xvec = range(xmin, xmax, length=n)
    if backwardsearch  # start from xmax and work backwards
        xvec = reverse(xvec)
    end

    fprev = f(xvec[1])
    for i = 2:n
        fnext = f(xvec[i])
        if fprev*fnext < 0  # bracket found
            if backwardsearch
                return true, xvec[i], xvec[i-1]
            else
                return true, xvec[i-1], xvec[i]
            end
        end
        fprev = fnext
    end

    return false, 0.0, 0.0

end

"""
    solve(rotor, section, inflow)

Solve the BEM equations for one section, with given inflow conditions, and rotor properties.
If multiple sections are to be solved (typical usage) then one can use broadcasting:
`solve.(sections, inflows, rotor)` where sections and inflows are arrays.

**Arguments**
- `rotor::Rotor`: rotor properties
- `section::Section`: section properties
- `inflow::Inflow`: inflow conditions

**Returns**
- `outputs::Outputs`: BEM output data including loads, induction factors, etc.
"""
function solve(rotor, section, inflow)

    # error handling
    if typeof(section) <: Array
        error("You passed in an array for section, but this funciton does not accept an array.\nProbably you intended to use broadcasting (notice the dot): solve.(sections, inflows, rotor)")
    end

    # parameters
    npts = 20  # number of discretization points to find bracket in residual solve

    # unpack
    swapsign = rotor.turbine ? 1 : -1
    Vx = inflow.Vx * swapsign  # TODO: shouldn't be doing the sign swap in two different places
    Vy = inflow.Vy
    theta = section.theta * swapsign

    # ---- determine quadrants based on case -----
    Vx_is_zero = isapprox(Vx, 0.0, atol=1e-6)
    Vy_is_zero = isapprox(Vy, 0.0, atol=1e-6)

    # quadrants
    epsilon = 1e-6
    q1 = [epsilon, pi/2]
    q2 = [-pi/2, -epsilon]
    q3 = [pi/2, pi-epsilon]
    q4 = [-pi+epsilon, -pi/2]

    if Vx_is_zero && Vy_is_zero
        return Outputs()

    elseif Vx_is_zero

        startfrom90 = false  # start bracket search from 90 deg instead of 0 deg.

        if Vy > 0 && theta > 0
            order = (q1, q2)
        elseif Vy > 0 && theta < 0
            order = (q2, q1)
        elseif Vy < 0 && theta > 0
            order = (q3, q4)
        else  # Vy < 0 && theta < 0
            order = (q4, q3)
        end

    elseif Vy_is_zero

        startfrom90 = true  # start bracket search from 90 deg

        if Vx > 0 && abs(theta) < pi/2
            order = (q1, q3)
        elseif Vx < 0 && abs(theta) < pi/2
            order = (q2, q4)
        elseif Vx > 0 && abs(theta) > pi/2
            order = (q3, q1)
        else  # Vx < 0 && abs(theta) > pi/2
            order = (q4, q2)
        end

    else  # normal case

        startfrom90 = false

        if Vx > 0 && Vy > 0
            order = (q1, q2, q3, q4)
        elseif Vx < 0 && Vy > 0
            order = (q2, q1, q4, q3)
        elseif Vx > 0 && Vy < 0
            order = (q3, q4, q1, q2)
        else  # Vx < 0 && Vy < 0
            order = (q4, q3, q2, q1)
        end

    end

    # ----- solve residual function ------

    # wrapper to residual function to accomodate format required by fzero
    function R(phi)
        zero, _ = residual(phi, section, inflow, rotor)
        return zero
    end

    success = false
    for j = 1:length(order)  # quadrant orders.  In most cases it should find root in first quadrant searched.
        phimin, phimax = order[j]

        # check to see if it would be faster to reverse the bracket search direction
        backwardsearch = false
        if !startfrom90
            if phimin == -pi/2 || phimax == -pi/2  # q2 or q4
                backwardsearch = true
            end
        else
            if phimax == pi/2  # q1
                backwardsearch = true
            end
        end

        # find bracket
        success, phiL, phiU = firstbracket(R, phimin, phimax, npts, backwardsearch)

        # once bracket is found, solve root finding problem and compute loads
        if success

            phistar = Roots.fzero(R, phiL, phiU)
            _, outputs = residual(phistar, section, inflow, rotor)

            return outputs
        end
    end

    # it shouldn't get to this point.  if it does it means no solution was found
    # it will return empty outputs
    # alternatively, one could increase npts and try again

    return Outputs()
end

# ------------ inflow ------------------

"""
    simpleinflow(Vinf, Omega, r, rho, mu=1.0, asound=1.0, precone=0.0)

Uniform inflow through rotor.  Returns an Inflow object.

**Arguments**
- `Vinf::Float64`: freestream speed (m/s)
- `Omega::Float64`: rotation speed (rad/s)
- `r::Float64`: radial location where inflow is computed (m)
- `precone::Float64`: precone angle (rad)
- `rho::Float64`: air density (kg/m^3)
"""
function simpleinflow(Vinf, Omega, r, rho, mu=1.0, asound=1.0, precone=0.0)

    # error handling
    if typeof(r) <: Array
        error("You passed in an array for r, but this function does not accept an array.\nProbably you intended to use broadcasting")
    end

    Vx = Vinf * cos(precone)
    Vy = Omega * r * cos(precone)

    return Inflow(Vx, Vy, rho, mu, asound)

end

"""
    windturbineinflow(Vinf, Omega, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho)

Compute relative wind velocity components along blade accounting for inflow conditions
and orientation of turbine.  See Documentation for angle definitions.

**Arguments**
- `Vhub::Float64`: freestream speed at hub (m/s)
- `Omega::Float64`: rotation speed (rad/s)
- `r::Array{Float64, 1}`: radial locations where inflow is computed (m)
- `precone::Float64`: precone angle (rad)
- `yaw::Float64`: yaw angle (rad)
- `tilt::Float64`: tilt angle (rad)
- `azimuth::Float64`: azimuth angle (rad)
- `hubHt::Float64`: hub height (m) - used for shear
- `shearExp::Float64`: power law shear exponent
- `rho::Float64`: air density (kg/m^3)
"""
function windturbineinflow(Vhub, Omega, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho, mu=1.0, asound=1.0)

    sy = sin(yaw)
    cy = cos(yaw)
    st = sin(tilt)
    ct = cos(tilt)
    sa = sin(azimuth)
    ca = cos(azimuth)
    sc = sin(precone)
    cc = cos(precone)

    # coordinate in azimuthal coordinate system
    x_az = -r*sin(precone)
    z_az = r*cos(precone)
    y_az = 0.0  # could omit (the more general case allows for presweep so this is nonzero)

    # get section heights in wind-aligned coordinate system
    heightFromHub = (y_az*sa + z_az*ca)*ct - x_az*st

    # velocity with shear
    V = Vhub*(1 + heightFromHub/hubHt)^shearExp

    # transform wind to blade c.s.
    Vwind_x = V * ((cy*st*ca + sy*sa)*sc + cy*ct*cc)
    Vwind_y = V * (cy*st*sa - sy*ca)

    # wind from rotation to blade c.s.
    Vrot_x = -Omega*y_az*sc
    Vrot_y = Omega*z_az

    # total velocity
    Vx = Vwind_x + Vrot_x
    Vy = Vwind_y + Vrot_y

    # operating point
    return Inflow(Vx, Vy, rho, mu, asound)

end

function windturbineinflow_az(Vhub, Omega, r, precone, yaw, tilt, azimuth_array, hubHt, shearExp, rho, mu=1.0, asound=1.0)

    naz = length(azimuth_array)
    azinflows = Array{Array{Inflow}}(undef, naz)
    for i = 1:naz
        azinflows[i] = windturbineinflow.(Vhub, Omega, r, precone, yaw, tilt, azimuth_array[i], hubHt, shearExp, rho)
    end

    return azinflows
end

# -------------------------------------

# -------- convenience methods ------------

"""
    loads(outputs)

Extract arrays for Np and Tp from the outputs.

This does not do any calculations.  It is simply a syntax shorthand as loads are usually what we mostly care about from the outputs.
"""
function loads(outputs)
    Np = getfield.(outputs, :Np)
    Tp = getfield.(outputs, :Tp)

    return Np, Tp
end

"""
    effectivewake(outputs)

Computes rotor wake velocities.

Note that the BEM methodology applies hub/tip losses to the loads rather than to the velocities.
This is the most common way to implement a BEM, but it means that the raw velocities are misleading
as they do not contain any hub/tip loss corrections.
To fix this we compute the effective hub/tip losses that would produce the same thrust/torque.
In other words:\n
CT = 4 a (1 - a) F = 4 a G (1 - a G)\n
This is solved for G, then multiplied against the wake velocities.
"""
function effectivewake(outputs)

    u = getfield.(outputs, :u)
    v = getfield.(outputs, :v)
    F = getfield.(outputs, :F)
    a = getfield.(outputs, :a)

    # the effective "F" if it were multiplied against the velocities instead of the forces
    G = (1.0 .- sqrt.(1.0 .- 4*a.*(1.0 .- a).*F))./(2*a)

    ueff = u.*G
    veff = v.*G

    return ueff, veff

end

"""
    thrusttorque(rotor, sections, outputs)

integrate the thrust/torque across the blade,
including 0 loads at hub/tip, using a trapezoidal rule.

**Arguments**
- `rotor::Rotor`: rotor object
- `sections::Array{Section, 1}`: section data along blade
- `outputs::Array{Outputs, 1}`: output data along blade

**Returns**
- `T::Float64`: thrust (along x-dir see Documentation)
- `Q::Float64`: torque (along x-dir see Documentation)
"""
function thrusttorque(rotor, sections, outputs)

    Np, Tp = loads(outputs)
    r = getfield.(sections, :r)

    # add hub/tip for complete integration.  loads go to zero at hub/tip.
    rfull = [rotor.Rhub; r; rotor.Rtip]
    Npfull = [0.0; Np; 0.0]
    Tpfull = [0.0; Tp; 0.0]

    # integrate Thrust and Torque (trapezoidal)
    thrust = Npfull*cos(rotor.precone)
    torque = Tpfull.*rfull*cos(rotor.precone)

    T = rotor.B * trapz(rfull, thrust)
    Q = rotor.B * trapz(rfull, torque)

    return T, Q
end

"""
ttodo
"""
function thrusttorque_azavg(rotor, sections, azinflows)

    T = 0.0
    Q = 0.0
    n = length(azinflows)

    for i = 1:n
        outputs = solve.(rotor, sections, azinflows[i])
        Tsub, Qsub = thrusttorque(rotor, sections, outputs)
        T += Tsub / n
        Q += Qsub / n
    end

    return T, Q
end

"""
(private) trapezoidal integration
"""
function trapz(x, y)  # integrate y w.r.t. x

    integral = 0.0
    for i = 1:length(x)-1
        integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
    end
    return integral
end

"""
    nondim(T, Q, Vhub, Omega, rho, rotor)

Nondimensionalize the outputs.

**Arguments**
- `T::Float64`: thrust (N)
- `Q::Float64`: torque (N-m)
- `Vhub::Float64`: hub speed used in turbine normalization (m/s)
- `Omega::Float64`: rotation speed used in propeller normalization (rad/s)
- `rho::Float64`: air density (kg/m^3)
- `rotor::Rotor`: rotor object

**Returns**

if windturbine
- `CP::Float64`: power coefficient
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient

if propeller
- `eff::Float64`: efficiency
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient
"""
function nondim(T, Q, Vhub, Omega, rho, rotor)

    P = Q * Omega
    Rp = rotor.Rtip*cos(rotor.precone)

    if rotor.turbine  # wind turbine normalizations

        q = 0.5 * rho * Vhub^2
        A = pi * Rp^2

        CP = P / (q * A * Vhub)
        CT = T / (q * A)
        CQ = Q / (q * Rp * A)

        return CP, CT, CQ

    else  # propeller

        n = Omega/(2*pi)
        Dp = 2*Rp

        if T < 0
            eff = 0.0  # creating drag not thrust
        else
            eff = T*Vhub/P
        end
        CT = T / (rho * n^2 * Dp^4)
        CQ = Q / (rho * n^2 * Dp^5)

        return eff, CT, CQ

    # elseif rotortype == "helicopter"

    #     A = pi * Rp^2

    #     CT = T / (rho * A * (Omega*Rp)^2)
    #     CP = P / (rho * A * (Omega*Rp)^3)
    #     FM = CT^(3/2)/(sqrt(2)*CP)

    #     return FM, CT, CP
    end

end

end  # module


text:
using Distances
import Distances: Metric, result_type, eval_reduce, eval_end, eval_op, eval_start, evaluate, parameters
const MinkowskiMetric = Union{Euclidean,Chebyshev,Cityblock,Minkowski,WeightedEuclidean,WeightedCityblock,WeightedMinkowski}
struct KDNode{T}
    lo::T           # The low boundary for the hyper rect in this dimension
    hi::T           # The high boundary for the hyper rect in this dimension
    split_val::T    # The value the hyper rectangle was split at
    split_dim::Int  # The dimension the hyper rectangle was split at
end

struct KDTree{V <: AbstractVector,M <: MinkowskiMetric,T}
    data::Vector{V}
    indices::Vector{Int}
    metric::M
    nodes::Vector{KDNode{T}}
    reordered::Bool
end

a = KDTree()


text:
"infinite vector with non-trivial values within a,b. Outside this range it simply repeats a default value"
struct RangeInfiniteVector{T} <: AbstractInfiniteVector{T}
    support::UnitRange{Int}
    default::T
    values::Vector{T}
    function RangeInfiniteVector{T}(support::UnitRange{Int}, default::T, values::Vector{T}) where {T}
        length(values) == length(support) || throw(ArgumentError("values and support must have the same length"))
        if support.start > support.stop
            support = 1:0  # standardize empty support
        end
        new(support, default, values)
    end
end
RangeInfiniteVector(support::UnitRange{Int}, default::T, values::Vector{T}) where {T} = RangeInfiniteVector{T}(support, default, values)
RangeInfiniteVector(support::UnitRange{Int}, default::T) where {T} = RangeInfiniteVector(support, default, fill(default, length(support)))
RangeInfiniteVector(a::Int, default::T, values::Vector{T}) where {T} = RangeInfiniteVector(a : a + length(values) - 1, default, values)
RangeInfiniteVector(a::Int, values::Vector{T}) where {T} = RangeInfiniteVector(a : a + length(values) - 1, zero(T), values)

EmptyRangeInfiniteVector(default::T) where {T} = RangeInfiniteVector(1:0, default)

"Resizes an infinite vector to a new support"
function resize(x::RangeInfiniteVector{T}, support::UnitRange{Int})::RangeInfiniteVector{T} where {T}
    RangeInfiniteVector{T}(support, x[support], x.default)
end

Base.getindex(x::RangeInfiniteVector{T}, i::Int) where {T} = i ∈ x.support ? x.values[i - x.support.start + 1] : x.default
Base.getindex(x::RangeInfiniteVector{T}, I) where {T} = [x[i] for i in I]

function Base.setindex!(x::RangeInfiniteVector{T}, v::T, i::Int) where {T}
    if i ∈ x.support
        x.values[i - x.support.start + 1] = v
    else
        throw(BoundsError())
    end
end

function Base.:+(x::RangeInfiniteVector, y::RangeInfiniteVector)
    a = min(x.support.start, y.support.start)
    b = max(x.support.stop, y.support.stop)
    values = [x[i] + y[i] for i = a:b]
    default = x.default + y.default
    return RangeInfiniteVector(a:b, default, values)
end

Base.:-(x::RangeInfiniteVector) = RangeInfiniteVector(x.support, -x.default, -x.values)
Base.:*(x::RangeInfiniteVector, y::Number) = RangeInfiniteVector(x.support, x.default * y, x.values * y)
Base.:*(x::Number, y::RangeInfiniteVector) = y * x
Base.:-(x::RangeInfiniteVector, y::RangeInfiniteVector) = x + (-y)
Base.:/(x::RangeInfiniteVector, y::Number) = RangeInfiniteVector(x.support, x.default / y, x.values / y)

maximum(x::RangeInfiniteVector) = max(x.default, maximum(x.values))
minimum(x::RangeInfiniteVector) = min(x.default, minimum(x.values))
extrema(x::RangeInfiniteVector) = extrema([x.default; x.values])