text:
"""
Rgen is an abstract type whose subtypes represent random
distributions and other array fillers.  Currently implemented subtypes
are listed below.

    * Gaussian(mean, std)
    * Uniform(min, max)
    * Constant(val)
    * Identity(scale)
    * Xavier()
    * Bernoulli(p,scale)
"""
abstract Rgen
type Gaussian <: Rgen; mean; std; Gaussian(mean=0,std=1)=new(mean,std); end
type Uniform  <: Rgen; min; max; Uniform(min=0,max=1)=new(min,max); end
type Constant <: Rgen; val; Constant(val=0)=new(val); end
type Identity <: Rgen; scale; Identity(x=1)=new(x); end
type Bernoulli <: Rgen; p; scale; Bernoulli(p=0.5, s=1)=new(p,s); end
type Xavier <: Rgen; scale; Xavier(scale=1)=new(scale); end # See http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf, Lasagne calls this GlorotUniform, See: http://lasagne.readthedocs.org/en/latest/modules/init.html#lasagne.init.GlorotUniform

function rgen!(r::Rgen, y)
    (isa(r, Constant)  ? fillsync!(y, r.val) :
     isa(r, Uniform)   ? axpb!(rand!(y); a=r.max - r.min, b=r.min) :
     isa(r, Gaussian)  ? axpb!(randn!(y); a=r.std, b=r.mean) :
     isa(r, Identity)  ? scale!(r.scale, copysync!(y, eye(eltype(y), size(y)...))) :
     isa(r, Bernoulli) ? bernoulli!(r.p, r.scale, rand!(y)) :
     isa(r, Xavier)    ? (s = r.scale*xavier(y); axpb!(rand!(y); a=2s, b=-s)) :
     error("Unknown Rgen=$r"))
    return y
end

bernoulli!{T}(p,s,x::Array{T})=(p=T(p);s=T(s);@inbounds for i=1:length(x); x[i] = (rand(T) < p ? s : 0); end; x)
@gpu bernoulli!(p,s,x::CudaArray{Float32})=(ccall((:bernoulli32,libknet), Void, (Cint,Cfloat,Cfloat,Ptr{Cfloat}), length(x), p, s, x); gpusync(); x)
@gpu bernoulli!(p,s,x::CudaArray{Float64})=(ccall((:bernoulli64,libknet), Void, (Cint,Cdouble,Cdouble,Ptr{Cdouble}), length(x), p, s, x); gpusync(); x)

function xavier(w)
     # The old implementation was not right for fully connected layers:
     # (fanin = length(y) / (size(y)[end]); scale = sqrt(3 / fanin); axpb!(rand!(y); a=2*scale, b=-scale)) :
    if ndims(w) < 2
        error("ndims=$(ndims(w)) in xavier")
    elseif ndims(w) == 2
        fanout = size(w,1)
        fanin = size(w,2)
    else
        fanout = size(w, ndims(w)) # Caffe disagrees: http://caffe.berkeleyvision.org/doxygen/classcaffe_1_1XavierFiller.html#details
        fanin = div(length(w), fanout)
    end
    # See: http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf
    sqrt(2 / (fanin + fanout))
end

function Base.isequal(a::Rgen,b::Rgen)
    typeof(a)==typeof(b) || return false
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end


text:
"""
    AbstractMeshingAlgorithm

Abstract type to specify an algorithm for isosurface extraction.
See:
- MarchingCubes
- MarchingTetrahedra
- NaiveSurfaceNets
"""
abstract type AbstractMeshingAlgorithm end

function (::Type{MeshAlgo})(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2, MeshAlgo <: AbstractMeshingAlgorithm}
    if isconcretetype(MeshAlgo)
        return MeshAlgo(iso, eps, reduceverts, insidepositive)
    else
        return MeshAlgo{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
    end
end

function (::Type{MeshAlgo})(iso) where MeshAlgo <: AbstractMeshingAlgorithm
    MeshAlgo(iso=iso)
end

function (::Type{MeshAlgo})(iso,eps) where MeshAlgo <: AbstractMeshingAlgorithm
    MeshAlgo(iso=iso,eps=eps)
end

"""
    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    MarchingCubes(iso)
    MarchingCubes(iso,eps)

Specifies the use of the Marching Cubes algorithm for isosurface extraction.
This algorithm provides a good balance between performance and vertex count.
In contrast to the other algorithms, vertices may be repeated, so mesh size
may be large and it will be difficult to extract topological/connectivity information.

- `iso` (default: 0.0) specifies the iso level to use for surface extraction.
- `eps` (default: 1e-3) is the tolerence around a voxel corner to ensure manifold mesh generation.
- `reduceverts` (default: true) if true will merge vertices within a voxel to reduce mesh size by around 30% and with slight performance improvement.
- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
"""
struct MarchingCubes{T} <: AbstractMeshingAlgorithm
    iso::T
    eps::T
    reduceverts::Bool
    insidepositive::Bool
end

"""
    MarchingTetrahedra(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    MarchingTetrahedra(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    MarchingTetrahedra(iso)
    MarchingTetrahedra(iso,eps)

Specifies the use of the Marching Tetrahedra algorithm for isosurface extraction.
This algorithm has a roughly 2x performance penalty compared to Marching Cubes,
and generates more faces. However, each vertex is guaranteed to not be repeated,
making this algorithm useful for topological analysis.

- `iso` specifies the iso level to use for surface extraction.
- `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
- `reduceverts` (default: true) if false, vertices will not be unique and have repeated face references.
- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
"""
struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
    iso::T
    eps::T
    reduceverts::Bool
    insidepositive::Bool
end

"""
    NaiveSurfaceNets(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    NaiveSurfaceNets(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
    NaiveSurfaceNets(iso)
    NaiveSurfaceNets(iso,eps)

Specifies the use of the Naive Surface Nets algorithm for isosurface extraction.
This algorithm has a slight performance advantage in some cases over Marching Cubes.
Each vertex is guaranteed to not be repeated (useful for topological analysis),
however the algorithm does not guarantee accuracy and generates quad faces.

- `iso` specifies the iso level to use for surface extraction.
- `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
- `reduceverts` reserved for future use.
- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
"""
struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
    iso::T
    eps::T
    reduceverts::Bool
    insidepositive::Bool
end

#
# Helper functions
#

default_face_length(::Union{MarchingCubes,MarchingTetrahedra}) = 3
default_face_length(::NaiveSurfaceNets) = 4

ext: jl


text:
"""
    init_weights(WeightsStruct, modelname, n[, pbc=true])
Generates single site and pair interactions Hamiltonians. Interacts with
`WEIGHTS` type, modelname passed as string.
"""
function init_weights(WeightsStruct, modelname, n, pbc=true)
    if modelname == "afh"
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n,n) ; w_yy = zeros(n,n); w_zz = zeros(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] = 1; w_yy[i, j] = 1; w_zz[i, j] = 1
        end
        w_xx[1, n] = 1; w_yy[1, n] = 1; w_zz[1, n] = 1
        w_xx += w_xx'; w_yy += w_yy'; w_zz += w_zz'

    elseif modelname == "U_afh"  # Marshall-Peierls Heisenberg model
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n,n) ; w_yy = zeros(n,n); w_zz = zeros(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] = -1; w_yy[i, j] = -1; w_zz[i, j] = 1
        end
        w_xx[1, n] = -1; w_yy[1, n] = -1; w_zz[1, n] = 1
        w_xx += w_xx'; w_yy += w_yy'; w_zz += w_zz'

    elseif modelname == "constr_x"
        w_x = -0.5 * ones(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n, n) ; w_yy = zeros(n,n) ; w_zz = zeros(n,n)
        for i=1:n -1
            j = i + 1
            w_zz[i, j] = 0
        end
        w_zz[1, n] = 0
        w_zz += w_zz'

    elseif modelname == "tfi_h1/2"
        w_x = -(1/2)*ones(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n, n) ; w_yy = zeros(n,n) ; w_zz = zeros(n,n)
        for i=1:n -1
            j = i + 1
            w_zz[i, j] = -1
        end
        w_zz[1, n] = - 1
        w_zz += w_zz'

    elseif modelname == "tfi_h1"
        w_x = -ones(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n, n) ; w_yy = zeros(n,n) ; w_zz = zeros(n,n)
        for i=1:n -1
            j = i + 1
            w_zz[i, j] = -1
        end
        w_zz[1, n] = - 1
        w_zz += w_zz'

    elseif modelname == "tfi_h2"
        w_x = -2 * ones(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n, n) ; w_yy = zeros(n,n) ; w_zz = zeros(n,n)
        for i=1:n -1
            j = i + 1
            w_zz[i, j] = -1
        end
        w_zz[1, n] = - 1
        w_zz += w_zz'

    elseif modelname == "afh_conn"
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = -ones(n,n) ; w_yy = -ones(n,n) ; w_zz = -ones(n,n)
        w_xx -= Diagonal(w_xx); w_yy -= Diagonal(w_yy); w_zz -= Diagonal(w_zz)

    elseif modelname == "afh_ext_y"
        w_x = zeros(n); w_y = -ones(n); w_z = zeros(n)
        w_xx = zeros(n,n) ; w_yy = zeros(n,n); w_zz = zeros(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] = -1; w_yy[i, j] = -1; w_zz[i, j] = -1
        end
        w_xx[1, n] = -1; w_yy[1, n] = -1; w_zz[1, n] = -1
        w_xx += w_xx'; w_yy += w_yy'; w_zz += w_zz'

    elseif modelname == "afh_xyz"
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = zeros(n,n) ; w_yy = zeros(n,n); w_zz = zeros(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] = -0.2; w_yy[i, j] = -0.6; w_zz[i, j] = -0.4
        end
        w_xx[1, n] = -0.2; w_yy[1, n] = -0.6; w_zz[1, n] = -0.4
        w_xx += w_xx'; w_yy += w_yy'; w_zz += w_zz'

    elseif modelname == "afh_xyz_ext"
        w_x = rand(Normal(0, 1/sqrt(n)), n); w_y = rand(Normal(0, 1/sqrt(n)), n); w_z = rand(Normal(0, 1/sqrt(n)), n)
        w_xx = zeros(n,n) ; w_yy = zeros(n,n); w_zz = zeros(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] = -0.2; w_yy[i, j] = -0.6; w_zz[i, j] = -0.4
        end
        w_xx[1, n] = -0.2; w_yy[1, n] = -0.6; w_zz[1, n] = -0.4
        w_xx += w_xx'; w_yy += w_yy'; w_zz += w_zz'

    elseif modelname == "afh_noisy"
        w_x = 0.01 * rand(n); w_y = 0.01 * rand(n); w_z = 0.01 * rand(n)
        w_xx = 0.01 * rand(n,n) ; w_yy = 0.01 * rand(n,n); w_zz = 0.01 * rand(n,n)
        for i=1:n - 1
            j=i + 1
            w_xx[i, j] -= 1; w_yy[i, j] -= 1; w_zz[i, j] -= 1
        end
        w_xx[1, n] -= 1; w_yy[1, n] -= 1; w_zz[1, n] -= 1
        w_xx = (w_xx + transpose(w_xx)) / 2; w_yy = (w_yy + transpose(w_yy)) / 2; w_zz = (w_zz + transpose(w_zz)) / 2
        w_xx -= Diagonal(w_xx); w_yy -= Diagonal(w_yy); w_zz -= Diagonal(w_zz)

    elseif modelname == "random_normal"
        w_x = rand(Normal(0,1), n); w_y = rand(Normal(0, 1), n); w_z = rand(Normal(0, 1), n)
        w_xx = rand(Normal(0, 1/sqrt(n)), n, n); w_yy = rand(Normal(0, 1/sqrt(n)), n, n); w_zz = rand(Normal(0, 1/sqrt(n)), n, n)
        w_xx = (w_xx + transpose(w_xx)) / 2; w_yy = (w_yy + transpose(w_yy)) / 2; w_zz = (w_zz + transpose(w_zz)) / 2
        w_xx -= Diagonal(w_xx); w_yy -= Diagonal(w_yy); w_zz -= Diagonal(w_zz)

    elseif modelname == "random_normal_ext0_symm"
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = w_yy = w_zz = rand(Normal(0, 1e-5), n, n)
        w_xx = (w_xx + transpose(w_xx)) / 2; w_yy = (w_yy + transpose(w_yy)) / 2; w_zz = (w_zz + transpose(w_zz)) / 2

    elseif modelname == "random"
        w_x = rand(n); w_y = rand(n); w_z = rand(n)
        w_xx = rand(n, n); w_yy = rand(n, n); w_zz = rand(n, n)
        w_xx = (w_xx + transpose(w_xx)) / 2; w_yy = (w_yy + transpose(w_yy)) / 2; w_zz = (w_zz + transpose(w_zz)) / 2
        w_xx -= Diagonal(w_xx); w_yy -= Diagonal(w_yy); w_zz -= Diagonal(w_zz)

    elseif modelname == "random_ext0"
        w_x = zeros(n); w_y = zeros(n); w_z = zeros(n)
        w_xx = rand(n, n); w_yy = rand(n, n); w_zz = rand(n, n)
        w_xx = (w_xx + transpose(w_xx)) / 2; w_yy = (w_yy + transpose(w_yy)) / 2; w_zz = (w_zz + transpose(w_zz)) / 2
        w_xx -= Diagonal(w_xx); w_yy -= Diagonal(w_yy); w_zz -= Diagonal(w_zz)

    else
        error("no model with that name")
    end

    if ~pbc && n>2
        w_xx[1, n] = w_yy[1, n] = w_zz[1, n] = 0
        w_xx[n, 1] = w_yy[n, 1] = w_zz[n, 1] = 0
    end

    weights = WeightsStruct(w_x, w_y, w_z, w_xx, w_yy, w_zz)
    return weights
end

"""
    Base.+(W1::WEIGHTS, W2::WEIGHTS)

Define additions of weight-structs for QBM learning rule.
"""
function Base.:+(W1::WEIGHTS, W2::WEIGHTS)
    SumW = init_weights(WEIGHTS, "zeros", size(W1.w_x, 1))
    for k in fieldnames(WEIGHTS)
        W1_tmp = getfield(W1, k)
        W2_tmp = getfield(W2, k)
        setfield!(SumW, k, W1_tmp + W2_tmp)
    end
    return SumW
end

"""
    Base.-(W1::WEIGHTS, W2::WEIGHTS)

Define additions of weight-structs for QBM learning rule.
"""
function Base.:-(W1::WEIGHTS, W2::WEIGHTS)
    subtrW = init_weights(WEIGHTS, "zeros", size(W1.w_x, 1))
    for k in fieldnames(WEIGHTS)
        W1_tmp = getfield(W1, k)
        W2_tmp = getfield(W2, k)
        setfield!(subtrW, k, W1_tmp - W2_tmp)
    end
    return subtrW
end

"""
    Base.*(x::Float64, W2::WEIGHTS)

Define multipliciation of weight-structs by a scalar for QBM learning rule.
"""
function Base.:*(x::Float64, W::WEIGHTS)
    timesW = init_weights(WEIGHTS, "zeros", size(W.w_x, 1))
    for k in fieldnames(WEIGHTS)
        timesW_tmp = x * getfield(W, k)
        setfield!(timesW, k, timesW_tmp)
    end
    return timesW
end


text:
include("geometryhelper.jl")
include("icosahedronQuadrature.jl")
include("octaQuadrature.jl")
include("standardQuadrature.jl")
include("dim2Quadrature.jl")
using FastGaussQuadrature

export Quadrature
export rotateandinterpolate!
export rotateOnly!

mutable struct Quadrature

	norder::Int64
	nquadpoints::Int64
	pointsxyz::Array{Float64,2}
	pointsmuphi::Array{Float64,2}
    weights::Array{Float64,1}

    muphiOLD::Array{Float64,2}
    xyzOLD::Array{Float64,2}

    triangulation::Array{Int64,2}

    interpolationids::Array{Int64,2}
	interpolationweights::Array{Float64,2}
	neighbors::Array{Int64,2}

    rotationmatrix::Array{Float64,2}

	quadtype::String;

	function Quadrature()
		return Quadrature(5,1)
	end

    function Quadrature(N::Int64,quadtype=1::Int)
		n = N
		quadTypeString = "quad"

		if quadtype==1 # standard quad
			xyz,weights = createPointStandards(n)
			triangulation = zeros(Int,1,1)
			neighbors = zeros(Int64,length(weights),6)
			nquadpoints = size(xyz)[1]
			quadTypeString = "Tensorized"

		elseif quadtype == 2 # wn octa
			xyz,triangulation =  createPointsOcta(n)
			nquadpoints = size(xyz)[1]
			weights = createWeights(nquadpoints,xyz,triangulation)
			neighbors = zeros(Int64,length(weights),6)
			for i=1:length(weights)
				tmp = (1 .+ rem.(findall(triangulation[:].==i) .- 1,size(triangulation,1)))
				neighbors[i,1:length(tmp)] = tmp;
			end
			quadTypeString = "Octaeder"
		elseif quadtype == 3 # wn ico
			xyz,triangulation =  createPointsIco(n)
			nquadpoints = size(xyz)[1]
			weights = createWeights(nquadpoints,xyz,triangulation)
			neighbors = zeros(Int64,length(weights),6)
			for i=1:length(weights)
				tmp = (1 .+ rem.(findall(triangulation[:] .==i) .- 1,size(triangulation,1)))
				neighbors[i,1:length(tmp)] = tmp;
			end
			quadTypeString = "Icosaeder"
		elseif quadtype == 4 # 2d
			xyz,weights = createPoint2d(n)
			triangulation = zeros(Int,1,1)
			neighbors = zeros(Int64,length(weights),6)
			nquadpoints = size(xyz)[1]
			quadTypeString = "test_2d"
		elseif quadtype == 5 # from file
			println("Please enter the quadrature file to be loaded and press enter.")
			println("(empty field + enter = quadrature.txt)")
			file = readline()
			if isempty(file)
				file = "quadrature.txt"
			end
			println("Trying to read $file.")
			M = readdlm(file,',')
			xyz = M[:,1:3];
			weights = M[:,4];
			triangulation = zeros(Int,1,1)
			neighbors = zeros(Int64,length(weights),6)
			nquadpoints = size(xyz)[1]
			quadTypeString = "fromfile"
		end

		muphi = zeros(Float64,nquadpoints,2)
		xyz2muphi!(xyz,muphi)

		interpolationids = zeros(Int64,nquadpoints,3)
		interpolationweights = zeros(Float64,nquadpoints,3)

		rotationmatrix = zeros(3,3);
		rotationmatrix[1,1] = 1.0
		rotationmatrix[2,2] = 1.0
		rotationmatrix[3,3] = 1.0

		xyzOLD = deepcopy(xyz)
		muphiOLD = deepcopy(muphi)

		new(n,nquadpoints,xyz,muphi,weights,muphiOLD,xyzOLD,triangulation,
			interpolationids,interpolationweights,neighbors,rotationmatrix,quadTypeString)
    end
end # struct

function createWeights(n::Int64,xyz::Array{Float64,2},triangles::Array{Int64,2})
	weights = zeros(n)
	#Compute the quadrature weights for each node
	#Therefore we iterate over all triangles and increase the weight
	#of each of the nodes of the triangle by the area of the triangle that
	#is assigned to it

	weights = zeros(n)
	nTriangles = size(triangles)[1]
	xy = zeros(Float64,3)
	yz = zeros(Float64,3)
	zx = zeros(Float64,3)
	mid = zeros(Float64,3)

	for n=1:nTriangles

		# get three nodes of a triangle
		i,j,k = triangles[n,:]

		# get the corresponding points
		x = xyz[i,:]
		y = xyz[j,:]
		z = xyz[k,:]

		# Now get the midpoint of the triangle and the midpoints along the lines
		mid = (x+y+z)/3.;
		xy = (x+y)/2.;
		yz = (y+z)/2.;
		zx = (z+x)/2.;

		# These points still have to be projected onto the sphere
		mid = mid/norm(mid,2);
		xy = xy/norm(xy,2);
		yz = yz/norm(yz,2);
		zx = zx/norm(zx,2);

		#By these four points, plus x,y,z we can span 6 triangles
		#Each of these triangles is assigned to one of the three nodes of the triangle x,y,z.
		#the area = the weight

		# for i
		weights[i] += getarea(x,mid,xy,"sphere")
		weights[i] += getarea(x,mid,zx,"sphere")

		# for j
		weights[j] += getarea(y,mid,xy,"sphere")
		weights[j] += getarea(y,mid,yz,"sphere")

		# for k
		weights[k] += getarea(z,mid,yz,"sphere")
		weights[k] += getarea(z,mid,zx,"sphere")

	end

    return weights
end

function applyinterpolation!(phiNew::Array{Float64,3},phi::Array{Float64,3},Q::Quadrature)
	nq,ni,nj = size(phi)
	#@parallel

	for j=1:nj
			for i=1:ni
				for q=1:nq
					id1 = Q.interpolationids[q,1]
					id2 = Q.interpolationids[q,2]
					id3 = Q.interpolationids[q,3]

					w1 = Q.interpolationweights[q,1]
					w2 = Q.interpolationweights[q,2]
					w3 = Q.interpolationweights[q,3]

					phiNew[q,i,j] = w1*phi[id1,i,j]+w2*phi[id2,i,j]+w3*phi[id3,i,j]

			end
		end
	end

	# fix to be conservative
	for j=1:nj
		for i=1:ni
			newnorm = 0.0;
			oldnorm = 0.0
			for q=1:nq
				newnorm += phiNew[q,i,j]*Q.weights[q]
				oldnorm += phi[q,i,j]*Q.weights[q]
			end
			if oldnorm>0.0
				for q=1:nq
					phiNew[q,i,j] *= (oldnorm/newnorm)
				end
			end
		end
	end
end

#Get the ID of the triangle to which v belongs
#if v belongs to multiple triangles because it is on an edge/vertix this
#returns still just one id as the interpolation would be the same,
#independent of the triangle id returned.
#If no triangle could be assigned, we throw an error
function gettriangleid(v::Array{Float64,1},Q::Quadrature,neighbors::Array{Int64,1})
	# this gives the rows of the triangulation matrix that contain i
	# check only neighbors first, usually this is enough
	for i=1:length(neighbors)
		if neighbors[i]>0 # the poles only have 4, not 6 neighbors
			a = Q.xyzOLD[Q.triangulation[neighbors[i],1],:]
			b = Q.xyzOLD[Q.triangulation[neighbors[i],2],:]
			c = Q.xyzOLD[Q.triangulation[neighbors[i],3],:]

			if isprojectedpointinsidetriangle(v,a,b,c)
				#We still have got to make sure that we don't pick the triangle from the opposite side.
				#Therefore it should be sufficient to check what the distance of the point
				#to one of the vertices of the triangle is.
				if norm(a-v)<1.0
					return neighbors[i];
				end
			end
		end
	end

	# check everything
	for i=1:size(Q.triangulation,1)
		a = Q.xyzOLD[Q.triangulation[i,1],:]
		b = Q.xyzOLD[Q.triangulation[i,2],:]
		c = Q.xyzOLD[Q.triangulation[i,3],:]

		if isprojectedpointinsidetriangle(v,a,b,c)
			#We still have got to make sure that we don't pick the triangle from the opposite side.
			#Therefore it should be sufficient to check what the distance of the point
			#to one of the vertices of the triangle is.
			if norm(a-v)<1.0
				return i;
			end
		end
	end
	println(v)
	println(norm(v)-1.0)
	error("We should not reach this part!")
end

function prepareinterpolation!(Q::Quadrature)
	# find a traingle id
	for i=1:Q.nquadpoints
		triangleID = gettriangleid(Q.pointsxyz[i,:],Q,Q.neighbors[i,:])
		Q.interpolationids[i,:] = Q.triangulation[triangleID,:]
		v = Q.pointsxyz[i,:]
		a = Q.xyzOLD[Q.interpolationids[i,1],:]
		b = Q.xyzOLD[Q.interpolationids[i,2],:]
		c = Q.xyzOLD[Q.interpolationids[i,3],:]

		Q.interpolationweights[i,:] = computeinterpolationweights(v,a,b,c)';
	end
end

function rotate!(Q::Quadrature)
	for i=1:size(Q.pointsxyz,1)
		Q.pointsxyz[i,:] = (Q.rotationmatrix*(Q.pointsxyz[i,:]))';
		Q.pointsxyz[i,:] = Q.pointsxyz[i,:]/norm(Q.pointsxyz[i,:])
	end
	xyz2muphi!(Q.pointsxyz,Q.pointsmuphi);

end

function rotateandinterpolate!(alpha::Float64,beta::Float64,gamma::Float64,
	Q::Quadrature,phiNew::Array{Float64,3},phi::Array{Float64,3},direction::String)

	getrotationmatrix!(Q.rotationmatrix,alpha,beta,gamma,direction)
	rotate!(Q)
	prepareinterpolation!(Q)
	applyinterpolation!(phiNew,phi,Q)
end

function rotateOnly!(u::Array{Float64,1},theta::Float64,
	Q::Quadrature,direction::String)
	getrotationmatrix!(Q.rotationmatrix,u,theta,direction)
	rotate!(Q)
	Q.xyzOLD = deepcopy(Q.pointsxyz)
	Q.muphiOLD = deepcopy(Q.pointsmuphi)

end

function rotateandinterpolate!(u::Array{Float64,1},theta::Float64,
	Q::Quadrature,phiNew::Array{Float64,3},phi::Array{Float64,3},direction::String)
	getrotationmatrix!(Q.rotationmatrix,u,theta,direction)
	rotate!(Q)
	prepareinterpolation!(Q)
	applyinterpolation!(phiNew,phi,Q)
	Q.xyzOLD = deepcopy(Q.pointsxyz)
	Q.muphiOLD = deepcopy(Q.pointsmuphi)

end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
### Effective sample size (ESS)

ess{N<:Real}(mcvariance::N, iidvariance::N, len::Integer) = len*iidvariance/mcvariance

ess(v::AbstractArray; vtype::Symbol=:imse, args...) = ess(mcvar(v; vtype=vtype, args...), mcvar_iid(v), length(v))

ess(v::AbstractArray, region; vtype::Symbol=:imse, args...) = mapslices(x -> ess(x; vtype=vtype, args...), v, region)

ess(s::VariableNState{Univariate}; vtype::Symbol=:imse, args...) = ess(s.value; vtype=vtype, args...)

ess(s::VariableNState{Multivariate}, i::Integer; vtype::Symbol=:imse, args...) = ess(s.value[i, :]; vtype=vtype, args...)

ess(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; vtype::Symbol=:imse, args...) =
  eltype(s)[ess(s, i; vtype=vtype, args...) for i in r]


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{PSY.RegulationDevice{T}, U},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection, U <: DeviceLimitedRegulation}
    devices = get_available_components(get_component_type(model), sys)
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        AreaBalancePowerModel,
    )

    add_variables!(container, DeltaActivePowerUpVariable, devices, U())
    add_variables!(container, DeltaActivePowerDownVariable, devices, U())
    add_variables!(container, AdditionalDeltaActivePowerUpVariable, devices, U())
    add_variables!(container, AdditionalDeltaActivePowerDownVariable, devices, U())
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{PSY.RegulationDevice{T}, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection}
    devices = get_available_components(get_component_type(model), sys)

    add_constraints!(
        container,
        RegulationLimitsConstraint,
        DeltaActivePowerUpVariable,
        devices,
        model,
        AreaBalancePowerModel,
    )

    add_constraints!(container, RampLimitConstraint, devices, model, AreaBalancePowerModel)
    add_constraints!(
        container,
        ParticipationAssignmentConstraint,
        devices,
        model,
        AreaBalancePowerModel,
    )
    objective_function!(container, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{PSY.RegulationDevice{T}, U},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection, U <: ReserveLimitedRegulation}
    devices = get_available_components(get_component_type(model), sys)
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        AreaBalancePowerModel,
    )

    add_variables!(container, DeltaActivePowerUpVariable, devices, U())
    add_variables!(container, DeltaActivePowerDownVariable, devices, U())
    add_variables!(container, AdditionalDeltaActivePowerUpVariable, devices, U())
    add_variables!(container, AdditionalDeltaActivePowerDownVariable, devices, U())
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{PSY.RegulationDevice{T}, ReserveLimitedRegulation},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection}
    devices = get_available_components(get_component_type(model), sys)

    add_constraints!(
        container,
        RegulationLimitsConstraint,
        DeltaActivePowerUpVariable,
        devices,
        model,
        AreaBalancePowerModel,
    )

    add_constraints!(container, ParticipationAssignmentConstraint, devices, model, S)
    objective_function!(container, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{PSY.RegulationDevice{T}, FixedOutput},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection}
    devices = get_available_components(get_component_type(model), sys)
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{PSY.RegulationDevice{T}, FixedOutput},
    ::Type{AreaBalancePowerModel},
) where {T <: PSY.StaticInjection}
    # There is no-op under FixedOutput formulation
    return
end


text:
using SunAsAStar
using Test
using Dates

@testset "SunAsAStar.jl" begin
    # Write your tests here.
    using SunAsAStar

    @testset "Common" begin
        @test_nowarn get_obs_loc(:WIYN)
        @test_nowarn get_obs_loc(:NEID)
        @test_nowarn get_obs_loc(:DCT)
        @test_nowarn get_obs_loc(:EXPRES)
        @test_nowarn get_obs_loc(:LAPALMA)
        @test_nowarn get_obs_loc(:HARPSN)
        @test_nowarn get_lst(Dates.datetime2julian(DateTime(2021,1,1)), obs=:WIYN )
        @test_nowarn get_lst([Dates.datetime2julian(DateTime(2020,1,1)), Dates.datetime2julian(DateTime(2021,1,1))], obs=:WIYN )
    end

    t = Dates.datetime2julian(DateTime(2021,1,1))
    @testset "Solar rotation" begin
        @test_nowarn SolarRotation.calc_Δfwhm_solar_rotation(t, obs=:WIYN)

    end
    @testset "Differential extinction" begin
        @test_nowarn DifferentialExtinction.calc_Δv_diff_extinction(t, obs=:WIYN)
        tlist = collect(t .+ range(0.0,stop=1.0,step=0.05))
        # TODO: Improve agreement.  Likely due to astropy breaking refraction=true option
        @test begin
            output = DifferentialExtinction.compare_Horizons_Astropy(tlist, :WIYN, max_airmass=4.0)
            maximum(abs.(output.delta_vr_old .- output.delta_vr_new)) < 0.01
        end
    end

    include("test_solar_rotation.jl")
end


text:
SetCardType = CardType{:set}()
const SetCard = NastranCard{SetCardType}

type AEFACT <: SetCard
    id::Int64
    factors::Vector{Float64}
end
function convert(::Type{AEFACT},card)
    id = card[2]::Int64
    i_br = 0
    for (i,f) in enumerate(card[3:end])
        if f == ""
            if i_br == 0
                i_br = i + 1
            end
        elseif i_br != 0
            throw(ExceptionError("Embedded blanks not allowed in AEFACT"))
        elseif typeof(f) != Float64
            throw(ExceptionError("AEFACT factors must be floats"))
        end
    end
    if i_br == 0
        factors = card[3:end]
    else
        factors = card[3:i_br]
    end
    AEFACT(id,factors)
end

type SET1 <: SetCard
    id::Int64
end
function convert(::Type{SET1},card)
    SET1(card[2]::Int64)
end

type AELIST <: SetCard
    id::Int64
end
function convert(::Type{AELIST},card)
    AELIST(card[2]::Int64)
end


text:
using DifferentialEquations
using Setfield

# Do the simulation taking the object state
function stepForward(
    objects:: Vector{<:Object}, # array of all objects
    fluid :: FluidParams,
    dt :: Float64, # time in each step
    n :: Int,  # number of steps
    maxPosChange :: Float64, # adapt dt so that position change doesn't exceed this in single time step
) :: Vector{<:Object}
    for j in 1:n
        remainingDt = dt

        while remainingDt > 0
            velocity( d ) = sum( map( s -> velocityField(d, s, fluid), objects ) )
            # velocityGradient( d ) = sum( map( s -> velocityGradientField(d, s, fluid), objects ) )
            velocityCurl( d ) = sum( map( s -> velocityCurlField(d, s, fluid), objects ) )

            vObjects = map( o -> velocity(o.location), objects ) # velocities that objects experience
            # vgObjects = map( o -> velocityGradient(o.location), objects ) # velocity gradients that objects experience
            vcObjects = map( o -> velocityCurl(o.location), objects ) # velocity curls that objects experience

            normVector( s ) = norm( [ s.theta, s.phi ] )
            maxV = maximum( map( normVector, vObjects ) )
            dt = min( remainingDt, maxPosChange / maxV )

            objects = map( (v,vv,vcv) -> stepForward(v, vv, vcv, dt), objects, vObjects, vcObjects )
            remainingDt -= dt
        end
    end

    return( objects )
end

# Do the simulation taking the object state
# NOTE: currently not using this, this is still incomplete
# TODO: add dynamics for dipole orientation
function simulate(
    objects :: Vector{<:Object}, # array of all objects
    fluid :: FluidParams,
    endTime :: Float64, # time in each step
    saveAtTime :: Float64,  # number of steps
    maxPosChange :: Float64, # adapt dt so that position change doesn't exceed this in single time step
) :: Tuple{Vector{Float64},<:Vector{<:Object}}
    n = length(objects) # number of particles

    function newLocation( o, loc )
        no = @set o.location = loc
    end

    locsToCoords(l) = [ SphericalCoord(l[2*i-1], l[2*i]) for i in 1:n ]

    function dynamics!(du,u,p,t)
        @assert n == round(Int, length(u)/2) # number of particles

        nObjects = map( newLocation, objects, locsToCoords(u) )
        velocity( d ) = sum( map( s -> velocityField(d, s, fluid), nObjects ) )

        vObjects = map( o -> velocity(o.location), nObjects ) # velocities that object experience

        newDU = reduce(
            vcat,
            map(
                (v,o) -> [v.theta, csc(o.location.theta)*v.phi],
                vObjects,
                nObjects
            )
        )
        for i in 1:(2*n)
            du[i] = newDU[i]
        end
    end

    u0 = reduce( vcat, map( o -> [ o.location.theta, o.location.phi], objects ) )
    tSpan = (0., endTime)
    prob = ODEProblem( dynamics!, u0, tSpan )
    sol = solve( prob, AutoVern7(Rodas5()), saveat=saveAtTime, reltol=maxPosChange, abstol=maxPosChange )

    simObjects = map(
        locs -> map( newLocation, objects, locsToCoords(locs) ),
        sol.u
    )

    return( sol.t, simObjects )
end


text:
abstract type Prior end

struct DefaultPrior <: Prior end

"""
    CustomPrior(predictors, intercept, auxiliary)

struct to hold information regarding user-specified custom priors.

# Usage

The `CustomPrior` struct has 3 fields:

1. `predictors`: the β coefficients.
2. `intercept`: the α intercept.
3. `auxiliary`: an auxiliary parameter.

In robust models, e.g. Linear Regression with Student-t likelihood or Count Regression
with Negative Binomial likelihood, often there is an extra auxiliary parameter that is
needed to parametrize to model to overcome under- or over-dispersion. If you are specifying
a custom prior for one of these type of models, then you should also specify a prior for
the auxiliary parameter.

Non-robust models do not need an auxiliary parameter and you can pass `nothing` as the
auxiliary argument.
"""
struct CustomPrior <: Prior
    predictors
    intercept
    auxiliary
end

function CustomPrior()
    throw(ArgumentError("CustomPrior has no default parameters, please use DefaultPrior"))
    return nothing
end


text:
@recipe function f(p::EmpiricalPattern; variables=[:rho, :vp], grad=:rainbow)
    @assert variables ⊆ names(p) "Variable name(s) not found"
    depth = cumsum(first(reals(p, :dh)))
    l = labels(p)
    units = ["kg/m³", "ft/s"]
    layout --> (1, length(variables) + 1)
    legend --> false
    yflip --> true
    size --> (800,400)
    @series begin
        subplot := 1
        seriestype --> :bar
        linewidth --> 0
        ylims --> (0.5, length(l) + 0.5)
        bar_width --> 1.
        seriescolor --> grad
        clims --> (1, 14)
        fill_z --> l
        orientation --> :horizontal
        yticks --> 1:length(l)
        xticks --> 1:14
        ylabel --> "Layer Index"
        xlabel --> "Flag"

        l
    end

    for (i, var) in enumerate(variables)
        realizations = reals(p, var)

        ylims --> (0, depth[end])
        ylabel --> "Depth, ft"
        for r in realizations
            @series begin
                subplot := i + 1
                linecolor --> :grey
                r, depth
            end
        end

        @series begin
            subplot := i + 1
            linecolor --> :red
            title --> string(var)
            xlabel --> units[i]
            mean(realizations), depth
        end

    end
end


text:
include("D:/Box Sync/git/Juliog/JuliogParameterize.jl")
include("D:/Box Sync/git/Juliog/JuliogSanitize.jl")
#include("D:/Box Sync/git/Juliog/JuliogSanitizeif.jl")
include("D:/Box Sync/git/Juliog/JuliogSanitizemacro.jl")
include("D:/Box Sync/git/Juliog/JULIOGtoVerilog.jl")
#include("D:/Box Sync/git/Juliog/JULIOGtoJulia.jl")
include("D:/Box Sync/git/Juliog/Juliog.jl")
include("D:/Box Sync/git/Juliog/mixer.jl")

#loadJULIOGexpr(mixer)
#@block MULT_3M5A "M" (DO, DI_X, DI_Y ; WL_X=10, WL_Y=10)
#@block SINCOS "SC" (COS, SIN, X)
#@block MIXER "mixer" (DO, DI, THETA, ENAB, CLK, RST ; WL=10, PL=2)

println("At end of include")

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
import YAML
config = YAML.load(open("config.yaml"))

# Points at which to create images
parameters = readdlm("parameters.dat")

npars = size(parameters, 1)

parnames = config["emulator"]["parameters"]

# Go through each of these and plot images, saving with the filename `image_001.out` etc,
# where each index corresponds to the line number in the parameters file.

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using HDF5

grd = config["grid"]
grid = Grid(grd)

# Load the parameters file from the configuration setting, and make this an object.
pars = convert_dict(config["parameters"], config["model"])
vel = pars.vel # [km/s]
npix = config["npix"] # number of pixels
species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

# Load the visibilities and sort
dv = DataVis(config["data_file"], true)[5]
# Conjugation is necessary for the SMA and ALMA
visibilities.conj!(dv)

qq = get_qq(dv)
# Find the indices that sort these in increasing order
ind = sortperm(qq)
qq = qq[ind]
uu = dv.uu[ind]
vv = dv.vv[ind]
nvis = length(qq)

VV = Array(Complex128, (nvis, npars))

# Iterate over rows
for i in 1:size(parameters,1)
    row = parameters[i,:]
    println("processing ", row)

    # Frequency will always be the first parameter
    nu, incl = row
    # Update the parameters
    pars.incl = incl

    # Read in the image
    image = @sprintf("image_%3.3i.out", i)
    im = imread("images/" * image)
    skim = imToSky(im, pars.dpc)
    corrfun!(skim) # alpha = 1.0

    # Determine dRA and dDEC from the image and distance
    dRA = abs(skim.ra[2] - skim.ra[1])/2. # [arcsec] the half-size of a pixel

    # FFT the image channel
    vis_fft = transform(skim)

    # Interpolate the `vis_fft` to the same uu,vv locations as our dataset
    # In the same order as the sorted qq vector
    for j=1:nvis
        VV[j,i] = interpolate_uv(uu[j], vv[j], vis_fft)
    end

end

# write to a file

fid = h5open("grid_VV.hdf5", "w")
fid["uu"] = uu
fid["vv"] = vv
fid["real"] = real(VV)
fid["imag"] = imag(VV)
close(fid)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
function close(ws::Worksheet)
    if ws.handle != C_NULL
        xls_close_WS(ws.handle)
        ws.handle = C_NULL
    end
end

function Worksheet(wb::Workbook, sheet_index::Integer)
    check_valid_sheetindex(wb, sheet_index)
    handle = xls_getWorkSheet(wb.handle, sheet_index - 1)
    if handle == C_NULL
        error("Couldn't open Worksheet $sheet_index.")
    end

    expect(xls_parseWorkSheet(handle), "Failed parsing sheet $sheet_index")

    # parse c struct xlsWorkSheet
    xlsws = unsafe_load(handle)

    return Worksheet(wb, sheet_index, handle, xlsws.rows)
end

@inline last_row_index(ws::Worksheet) = ws.rows.lastrow + 1
@inline last_column_index(ws::Worksheet) = ws.rows.lastcol + 1

const MAX_WORKSHEET_ROW_INDEX = typemax(UInt16)

# See implementation of `xlsRow *xls_row(xlsWorkSheet* pWS, WORD cellRow)`
@inline is_valid_worksheet_row(ws::Worksheet, row::Integer) = (
    (row <= MAX_WORKSHEET_ROW_INDEX)
    && (0 < row <= last_row_index(ws))
    && ws.rows.row != C_NULL
)

@inline is_valid_worksheet_column(ws::Worksheet, column::Integer) = 0 < column <= last_column_index(ws)

@inline check_valid_worksheet_column(ws::Worksheet, column::Integer) = @assert is_valid_worksheet_column(ws, column) "Worksheet column out of bounds: $column."
@inline check_valid_worksheet_row(ws::Worksheet, row::Integer) = @assert is_valid_worksheet_row(ws, row) "Worksheet Row out of bounds: $row."
Base.size(ws::Worksheet) = (last_row_index(ws), last_column_index(ws))

function WorksheetRow(ws::Worksheet, row::Integer)
    check_valid_worksheet_row(ws, row)

    # adds to buffer if not present
    if row ∉ keys(ws.worksheet_rows)
        row_data = unsafe_load(ws.rows.row, row)
        worksheet_row = WorksheetRow(ws, row_data)
        ws.worksheet_rows[row] = worksheet_row
    end

    return ws.worksheet_rows[row]
end

function celldata(ws::Worksheet, row::Integer, col::Integer)::st_cell_data
    check_valid_worksheet_column(ws, col)
    wsrow = WorksheetRow(ws, row)

    if col ∉ keys(wsrow.cell_data)
        cell_data_ptr = wsrow.row_data.cells.cell
        cell_data = unsafe_load(cell_data_ptr, col)
        wsrow.cell_data[col] = cell_data
    end

    return wsrow.cell_data[col]
end

function Base.getindex(ws::Worksheet, row::Integer, column::Integer)
    cell = celldata(ws, row, column)
    cell_record = XLSRecord(cell.id)

    if cell_record == XLS_RECORD_NUMBER || cell_record == XLS_RECORD_RK
        return cell.d
    elseif cell_record == XLS_RECORD_BLANK
        return missing
    elseif cell_record == XLS_RECORD_LABELSST
        return unsafe_string(cell.str)
    else
        error("Unsupported Record: $cell_record.")
    end
end

sheetindex(ws::Worksheet) = ws.sheet_index
sheetname(ws::Worksheet) = sheetname(ws.parent, sheetindex(ws))


text:
using Test
using NQCBase
using PyCall

@testset "Single frame" begin
    atoms = Atoms([:H, :C, :O, :N])
    cell = PeriodicCell([10 0 0; 0 10 0; 0 0 10])
    R = rand(3, 4) .* 10

    ase_atoms = convert_to_ase_atoms(atoms, R, cell)
    new_atoms, new_R, new_cell = convert_from_ase_atoms(ase_atoms)

    @test new_atoms == atoms
    @test new_cell.vectors ≈ cell.vectors
    @test new_cell.inverse ≈ cell.inverse
    @test new_cell.periodicity ≈ cell.periodicity
    @test new_R ≈ R
end

@testset "Multiple frames" begin
    atoms = Atoms([:H, :C, :O, :N])
    cell = PeriodicCell([10 0 0; 0 10 0; 0 0 10])
    R = [rand(3, 4) .* 10 for _=1:100]

    ase_atoms = convert_to_ase_atoms(atoms, R, cell)
    out = convert_from_ase_atoms.(ase_atoms)
    @test out isa Vector{<:Tuple}
end


text:
using LinearAlgebra
using StaticArrays
using Rotations
using Cartesian
using Test

@testset "Cartesian.jl" begin
    @show @test Vector3()===SVector{3}([0.0, 0.0, 0.0])
    @show @test Vector3(1,2,3)===SVector{3}([1.0, 2.0, 3.0])
    @show @test Vector3([1,2,3])===SVector{3}([1.0, 2.0, 3.0])

    @show @test Matrix3()===SMatrix{3,3}([0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])
    @show @test Matrix3(0*I)===SMatrix{3,3}([0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])
    @show @test Matrix3(1*I)=== SMatrix{3,3}([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])
    @show @test Matrix3(1,2,3,4,5,6,7,8,9) === SMatrix{3,3}([1.0 2.0 3.0; 4.0  5.0  6.0; 7.0  8.0  9.0])

    @show @test ô === Vector3()
    @show @test î === Vector3(1,0,0)
    @show @test ĵ === Vector3(0,1,0)
    @show @test k̂ === Vector3(0,0,1)

    @show @test Ô === Matrix3(0*I)
    @show @test Î === Matrix3(1*I)

    @show r = Vector3(1,2,3)
    @show @test r.X == 1 && r.Y == 2 && r.Z == 3
    @show @test r.SumSquares == 14
    @show @test r.Magnitude == sqrt(14)
    @show n = normalize(r)
    @show @test n == SVector{3}([1.0,2.0,3.0]/sqrt(14))
    @show @test n.Magnitude == 1.0

    @show g = Vector3(-2,3,2)

    @show @test dot(r,g)==10 && r⋅g == 10
    @show @test cross(r,g)===Vector3(-5,-8,7) && r×g === Vector3(-5,-8,7)

    @show @test cross(r)*g === cross(r,g)
    @show @test cross(r) === Matrix3(0,-3,2, 3,0,-1, -2,1,0) && ×(r) === cross(r)
    @show @test cross2(r) === cross(r)*cross(r) && ×(r)*×(r) === cross2(r)
    @show @test -×(r)*×(g) === Matrix3(12, 4, 6, -3, 4, -9, -2, -4, 4)

    @show A = Matrix3(5,-1,2,-1,5,3,-1,0,3)
    @show @test A.Row1 === Vector3(A[1,:])
    @show @test A.Row2 === Vector3(A[2,:])
    @show @test A.Row3 === Vector3(A[3,:])
    @show @test A.Column1 === Vector3(A[:,1])
    @show @test A.Column2 === Vector3(A[:,2])
    @show @test A.Column3 === Vector3(A[:,3])
    @show @test A.Trace == 13.0
    @show @test A.Determinat == 85.0

    @show @test isapprox(85*(A\r) , Vector3(-18, -17, 79))
    @show @test isapprox( 85*inv(A), Matrix3(15,3,-13, 0, 17, -17, 5, 1, 24))

    @show @test isapprox(angle(r, g), 0.86557323388)
    @show @test isapprox(slerp(r,g, 1.2098312542), Vector3([-2.5118334954,2.9375160105,1.5629943723]))

end


text:
abstract type ActionSelector end

rand_action(gp, actions) = rand(actions)

function build_gp(action_dims, belief_dims, mean_type, mean_params, kernel_type, kernel_params, log_noise, k_neighborhood)
    mean = isnothing(mean_params) ? mean_type() : mean_type(mean_params...)
    total_dims = action_dims + belief_dims
    kernel = isnothing(kernel_params) ? kernel_type() : kernel_type(kernel_params...)
    x = zeros(Float64, total_dims, 0)
    y = zeros(Float64, 0)
    GPLA(x, y, k_neighborhood, action_dims, belief_dims, mean, kernel, log_noise)
end

"""
Bayesian Optimization Action Selector
"""
mutable struct BOActionSelector{M, K, IA, GP, D} <: ActionSelector
    action_dims::Int64
    belief_dims::Int64
    discrete_actions::Bool
    k_neighbors::Int64
    mean_type::M
    mean_params::Union{Array, Nothing}
    kernel_type::K
    kernel_params::Union{Array, Nothing}
    belief_λ::Float64 # λ < 0 Use separate GPs
    log_noise::Float64
    buffer_size::Int64
    acquisition_function::Union{AbstractAcquisition, Function}
    lower_bounds::Union{Function, AbstractArray, Nothing}
    upper_bounds::Union{Function, AbstractArray, Nothing}
    optim::Union{NLopt.Opt, Nothing}
    initial_action::IA
    gp::GP
    b_a_dict::Dict
    obs_buffer::Tuple{Matrix{Float64}, Vector{Float64}} # Number of observations to store
    function BOActionSelector{M, K, IA, GP, D}(
                            gp::GP,
                            action_dims::Int64,
                            belief_dims::Int64,
                            discrete_actions::Bool;
                            k_neighbors::Int64 = 10,
                            mean_type::M = MeanConst,
                            mean_params::Union{Array, Nothing} = [0.],
                            kernel_type::K = SEIso,
                            kernel_params::Union{Array, Nothing} = [0. 0.],
                            belief_λ::Float64 = 1.0,
                            log_noise::Float64 = -3.0,
                            buffer_size::Int64 = 0,
                            acquisition_function = ExpectedImprovement(),
                            lower_bounds = nothing,
                            upper_bounds = nothing,
                            optim = nothing,
                            initial_action::IA = rand_action
                            ) where {M, K, IA, GP, D}
        obs_buffer = (zeros(Float64, action_dims + belief_dims, 0), zeros(Float64, 0),)
        if discrete_actions
            return new(action_dims, belief_dims, discrete_actions, k_neighbors, mean_type, mean_params, kernel_type, kernel_params, belief_λ, log_noise, buffer_size, acquisition_function, nothing, nothing, nothing, initial_action, gp, Dict(), obs_buffer)
        else
            if isnothing(lower_bounds) || isnothing(upper_bounds)
                throw(ArgumentError("Lower and upper bounds must be provided for continuous action spaces"))
            end
            if isnothing(optim)
                optim = NLopt.Opt(:LD_LBFGS, action_dims)
                setproperty!(optim, :maxeval, 10)
                NLopt.lower_bounds!(optim, lower_bounds)
                NLopt.upper_bounds!(optim, upper_bounds)
            end
            return new(action_dims, belief_dims, discrete_actions, k_neighbors, mean_type, mean_params, kernel_type, kernel_params, belief_λ, log_noise, buffer_size, acquisition_function, lower_bounds, upper_bounds, optim, initial_action, gp, Dict(), obs_buffer)
        end
    end
end

function BOActionSelector(action_dims::Int64,
                        belief_dims::Int64,
                        discrete_actions::Bool;
                        k_neighbors::Int64 = 10,
                        mean_type::Any = MeanConst,
                        mean_params::Union{Array, Nothing} = [0.],
                        kernel_type::Any = SEIso,
                        kernel_params::Union{Array, Nothing} = [0. 0.],
                        belief_λ::Float64 = 1.0,
                        log_noise::Float64 = -3.0,
                        buffer_size::Int64 = 0,
                        acquisition_function = ExpectedImprovement(),
                        lower_bounds = nothing,
                        upper_bounds = nothing,
                        optim = nothing,
                        initial_action::Any = rand_action,
                        )
    if belief_λ < 0
        gp = Dict{Int64, GPE}()
    else
        gp = build_gp(action_dims, belief_dims, mean_type, mean_params, kernel_type, kernel_params, log_noise, k_neighbors)
    end
    BOActionSelector{typeof(mean_type), typeof(kernel_type), typeof(initial_action), typeof(gp), discrete_actions}(gp,
                                                                                action_dims,
                                                                                belief_dims,
                                                                                discrete_actions,
                                                                                k_neighbors=k_neighbors,
                                                                                mean_type=mean_type,
                                                                                mean_params=mean_params,
                                                                                kernel_type=kernel_type,
                                                                                kernel_params=kernel_params,
                                                                                belief_λ=belief_λ,
                                                                                log_noise=log_noise,
                                                                                buffer_size=buffer_size,
                                                                                acquisition_function=acquisition_function,
                                                                                lower_bounds=lower_bounds,
                                                                                upper_bounds=upper_bounds,
                                                                                optim=optim,
                                                                                initial_action=initial_action
                                                                                )
end


text:
"""
Authors : Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020
"""

include("./helpers/generate.jl")

using .Generator

T = 12
M = 5
R = 6
Emax = 5
f = [[10, 30, 60, 90],[10, 30, 60, 90, 120, 150],[10, 30, 60, 90, 120, 150, 180, 210]]
e = [[8, 7, 6, 5],[8, 7, 6, 5, 4, 3],[8, 7, 6, 5, 4, 3 ,2 ,1]]
h = 1
p = [[0, 0, 0, 0],[0, 0, 0, 0, 0, 0],[0, 0, 0, 0, 0, 0, 0, 0]]
dmin = 20
dmax = 70
# d = Array{Int,1}(undef, T)
d=[48, 34, 35, 69, 43, 70, 45, 62, 41, 25, 37, 46];

# for i in 1:T
# 	d[i] = rand(dmin:dmax)
# end

for r in 1:T
	generate(T, 4, r, Emax, f[1], e[1], h, p[1], d,"instances/var_M/set_4/inst_$(10+r).dat")
end

for r in 1:T
	generate(T, 6, r, Emax, f[2], e[2], h, p[2], d,"instances/var_M/set_6/inst_$(10+r).dat")
end

for r in 1:T
	generate(T, 8, r, Emax, f[3], e[3], h, p[3], d,"instances/var_M/set_8/inst_$(10+r).dat")
end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
using ReservoirComputing

const approx_res_size = 12
const alpha = 1.0
const beta = 1*10^(-5)
const extended_states = false
const nla_type = NLADefault()
const activation = tanh

const input_weight = 0.1
const cyrcle_weight = 0.99
const jump_weight = 0.1
const jumps = 12

const memory_size = 16

input = ones(Float64, 10, 2)
output = ones(Float64, 10, 2)

W = CRJ(approx_res_size, cyrcle_weight, jump_weight, jumps)
W_in = irrational_sign_input(approx_res_size, size(input, 2), input_weight)

rmm = RMM(W, input, output, W_in, memory_size, beta,
        extended_states = extended_states,
        nla_type = nla_type,
        activation = activation,
        alpha = alpha)

@test isequal(approx_res_size, rmm.res_size)
@test isequal(rmm.in_size, size(input, 2))
@test isequal(rmm.out_size, size(output, 2))
@test isequal(input, rmm.train_data)
@test isequal(alpha, rmm.alpha)
@test isequal(nla_type, rmm.nla_type)
@test isequal(activation, rmm.activation)
@test size(rmm.W) == (rmm.res_size, rmm.res_size)
@test size(rmm.W_in) == (rmm.res_size, rmm.in_size)
@test isequal(rmm.extended_states, extended_states)
@test size(rmm.W_out) == (rmm.out_size, rmm.res_size+rmm.in_size)
#test direct predict
rmmout = RMMdirect_predict(rmm, input)
@test size(rmmout) == size(input)


text:
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using Endowment, DataFrames, Distributed, Finometrics, CodecZlib, CSV, GZip#, CSVFiles#, CSV

const WORKERS=1

const CLEAN_WORKERS=false

#set up for parallel operations
initializePar(WORKERS, reset=CLEAN_WORKERS)

#a list of high-level functions to run at the script level
function endowmentTasks()::Nothing
  #downloadIRSPDFs()

  #metadata::Vector{MetaDictionary} = preProcessData(refreshData=false,
  #  validateDownloads=false,  refreshMetadata = false, reserializeData = false,
  #  refreshSummary = true)

  #NCCSData(metadata, refreshProcessed=true, refreshNCCSDataCSV=false, mode=:parallel)
  #NOTE: Persistance plots= return persistance plots
  makeDescriptive(refreshFilter = false, refreshDescriptiveFigures = false,
    refreshDescriptiveTables = false, doOneOffTasks = false, refreshβ=false,
    refreshSP500=false, refreshPersistancePlots=true)

  return nothing
end

@time endowmentTasks()

#cleanup the extra processes if desired
if CLEAN_WORKERS
  cleanup()
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module TestVariableContainer

using Test
import VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_empty()
    a = MOI.Utilities.VariablesContainer([0x0000, 0x0000], [1, 2], [3, 4])
    MOI.empty!(a)
    @test MOI.is_empty(a)
    return
end

function test_resize()
    a = MOI.Utilities.VariablesContainer([0x0000, 0x0000], [1, 2], [3, 4])
    @test length(a.set_mask) == 2
    @test length(a.lower) == 2
    @test length(a.upper) == 2
    resize!(a, 4)
    @test length(a.set_mask) == 4
    @test length(a.lower) == 4
    @test length(a.upper) == 4
    return
end

function test_add_variable()
    a = MOI.Utilities.VariablesContainer{Int}()
    MOI.add_variable(a)
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0000], [0], [0])
    a = MOI.Utilities.VariablesContainer{Float64}()
    MOI.add_variable(a)
    @test a ==
          MOI.Utilities.VariablesContainer{Float64}([0x0000], [-Inf], [Inf])
    return
end

function test__flag_to_set_type()
    T = Int
    @test_throws AssertionError MOI.Utilities._flag_to_set_type(0x0011, T)
    @test MOI.Utilities._flag_to_set_type(0x0010, T) == MOI.Integer
    @test MOI.Utilities._flag_to_set_type(0x0020, T) == MOI.ZeroOne
    return
end

function test_add_constraint()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0000], Int[0], Int[0])
    MOI.add_constraint(a, x, MOI.GreaterThan(3))
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0002], [3], [0])
    MOI.add_constraint(a, x, MOI.LessThan(4))
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0006], [3], [4])
    return
end

function test_add_constraint_LowerBoundAlreadySet()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    MOI.add_constraint(a, x, MOI.GreaterThan(3))
    @test_throws(
        MOI.LowerBoundAlreadySet{MOI.GreaterThan{Int},MOI.GreaterThan{Int}},
        MOI.add_constraint(a, x, MOI.GreaterThan(3)),
    )
    return
end

function test_add_constraint_UpperBoundAlreadySet()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    MOI.add_constraint(a, x, MOI.LessThan(3))
    @test_throws(
        MOI.UpperBoundAlreadySet{MOI.LessThan{Int},MOI.LessThan{Int}},
        MOI.add_constraint(a, x, MOI.LessThan(3)),
    )
    return
end

function test_delete_constraint_LessThan()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    c = MOI.add_constraint(a, x, MOI.LessThan(3))
    @test MOI.is_valid(a, c)
    MOI.delete(a, c)
    @test !MOI.is_valid(a, c)
    return
end

function test_delete_variable()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0000], Int[0], Int[0])
    c = MOI.add_constraint(a, x, MOI.LessThan(3))
    MOI.delete(a, x)
    @test a == MOI.Utilities.VariablesContainer{Int}([0x8000], Int[0], Int[3])
    return
end

function test_delete_constraint_GreaterThan()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    c = MOI.add_constraint(a, x, MOI.GreaterThan(3))
    @test MOI.is_valid(a, c)
    MOI.delete(a, c)
    @test !MOI.is_valid(a, c)
    return
end

function test_set_ConstraintSet()
    a = MOI.Utilities.VariablesContainer{Int}()
    x = MOI.add_variable(a)
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0000], Int[0], Int[0])
    c = MOI.add_constraint(a, x, MOI.GreaterThan(3))
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0002], Int[3], Int[0])
    MOI.set(a, MOI.ConstraintSet(), c, MOI.GreaterThan(2))
    @test a == MOI.Utilities.VariablesContainer{Int}([0x0002], Int[2], Int[0])
    return
end

function test_NumberOfConstraints()
    b = MOI.Utilities.VariablesContainer(
        [0x0008, 0x0002, 0x0004, 0x0006],
        [1.0, 3.0, -Inf, -1.0],
        [2.0, Inf, 4.0, 1.0],
    )
    get(S) = MOI.get(b, MOI.NumberOfConstraints{MOI.VariableIndex,S}())
    @test get(MOI.ZeroOne) == 0
    @test get(MOI.GreaterThan{Float64}) == 2
    @test get(MOI.LessThan{Float64}) == 2
    @test get(MOI.Interval{Float64}) == 1
    return
end

function test_ListOfConstraintTypesPresent()
    b = MOI.Utilities.VariablesContainer(
        [0x0008, 0x0002, 0x0004, 0x0006],
        [1.0, 3.0, -Inf, -1.0],
        [2.0, Inf, 4.0, 1.0],
    )
    @test MOI.get(b, MOI.ListOfConstraintTypesPresent()) == [
        (MOI.VariableIndex, MOI.GreaterThan{Float64}),
        (MOI.VariableIndex, MOI.LessThan{Float64}),
        (MOI.VariableIndex, MOI.Interval{Float64}),
    ]
    return
end

function test_ListOfConstraintIndices()
    b = MOI.Utilities.VariablesContainer{Float64}(
        [0x0008, 0x0002, 0x0004, 0x0006],
        [1.0, 3.0, -Inf, -1.0],
        [2.0, Inf, 4.0, 1.0],
    )
    get(S) = MOI.get(b, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
    @test get(MOI.ZeroOne) == []
    @test get(MOI.GreaterThan{Float64}) == [
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(2),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(4),
    ]
    @test get(MOI.LessThan{Float64}) ==
          MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}.([3, 4])
    @test get(MOI.Interval{Float64}) ==
          [MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}}(1)]
    return
end

###
### Hyperrectangle
###

function test_Hyperrectangle_equal()
    a = MOI.Utilities.Hyperrectangle([1, 2], [3, 4])
    b = MOI.Utilities.Hyperrectangle([1.0, 2.0], [3.0, 4.0])
    c = MOI.Utilities.Hyperrectangle([1.0, 3.0], [3.0, 4.0])
    d = MOI.Utilities.Hyperrectangle([1.0, 2.0], [3.0, 5.0])
    @test a == a
    @test a == b
    @test a != c
    @test a != d
    return
end

function test_Hyperrectangle_empty()
    a = MOI.Utilities.Hyperrectangle([1, 2], [3, 4])
    empty!(a)
    @test a.lower == Int[]
    @test a.upper == Int[]
    return
end

function test_Hyperrectangle_resize()
    a = MOI.Utilities.Hyperrectangle([1, 2], [3, 4])
    @test length(a.lower) == 2
    @test length(a.upper) == 2
    resize!(a, 4)
    @test length(a.lower) == 4
    @test length(a.upper) == 4
    return
end

function test_load_constants()
    a = MOI.Utilities.Hyperrectangle([-Inf, -Inf, -Inf], [Inf, Inf, Inf])
    MOI.Utilities.load_constants(a, 0, MOI.Interval(1.0, 2.0))
    MOI.Utilities.load_constants(a, 1, MOI.GreaterThan(3.0))
    MOI.Utilities.load_constants(a, 2, MOI.LessThan(4.0))
    @test a == MOI.Utilities.Hyperrectangle([1.0, 3.0, -Inf], [2.0, Inf, 4.0])
    return
end

function test_function_constants()
    a = MOI.Utilities.Hyperrectangle([-Inf, -Inf, -Inf], [Inf, Inf, Inf])
    MOI.Utilities.function_constants(a, 0) == 0.0
    return
end

function test_set_from_constants()
    a = MOI.Utilities.Hyperrectangle([1.0, 3.0, -Inf], [2.0, Inf, 4.0])
    @test MOI.Utilities.set_from_constants(a, MOI.Interval{Float64}, 1) ==
          MOI.Interval(1.0, 2.0)
    @test MOI.Utilities.set_from_constants(a, MOI.GreaterThan{Float64}, 2) ==
          MOI.GreaterThan(3.0)
    @test MOI.Utilities.set_from_constants(a, MOI.LessThan{Float64}, 3) ==
          MOI.LessThan(4.0)
    @test MOI.Utilities.set_from_constants(a, MOI.ZeroOne, 2) == MOI.ZeroOne()
    return
end

function _test_bound_vectors(::Type{T}, nolb, noub) where {T}
    variable_bounds = MOI.Utilities.VariablesContainer{T}()
    @test variable_bounds.lower == T[]
    @test variable_bounds.upper == T[]
    x = MOI.add_variable(variable_bounds)
    @test variable_bounds.lower == [nolb]
    @test variable_bounds.upper == [noub]
    ux = MOI.add_constraint(variable_bounds, x, MOI.LessThan(T(1)))
    @test variable_bounds.lower == [nolb]
    @test variable_bounds.upper == [T(1)]
    y = MOI.add_variable(variable_bounds)
    @test variable_bounds.lower == [nolb, nolb]
    @test variable_bounds.upper == [T(1), noub]
    cy = MOI.add_constraint(variable_bounds, y, MOI.Interval(T(2), T(3)))
    @test variable_bounds.lower == [nolb, T(2)]
    @test variable_bounds.upper == [T(1), T(3)]
    lx = MOI.add_constraint(variable_bounds, x, MOI.GreaterThan(T(0)))
    @test variable_bounds.lower == [T(0), T(2)]
    @test variable_bounds.upper == [T(1), T(3)]
    MOI.delete(variable_bounds, lx)
    @test variable_bounds.lower == [nolb, T(2)]
    @test variable_bounds.upper == [T(1), T(3)]
    MOI.delete(variable_bounds, ux)
    @test variable_bounds.lower == [nolb, T(2)]
    @test variable_bounds.upper == [noub, T(3)]
    cx = MOI.add_constraint(variable_bounds, x, MOI.Semicontinuous(T(3), T(4)))
    @test variable_bounds.lower == [T(3), T(2)]
    @test variable_bounds.upper == [T(4), T(3)]
    MOI.delete(variable_bounds, cy)
    @test variable_bounds.lower == [T(3), nolb]
    @test variable_bounds.upper == [T(4), noub]
    sy = MOI.add_constraint(variable_bounds, y, MOI.Semiinteger(T(-2), T(-1)))
    @test variable_bounds.lower == [T(3), T(-2)]
    @test variable_bounds.upper == [T(4), T(-1)]
    MOI.delete(variable_bounds, sy)
    @test variable_bounds.lower == [T(3), nolb]
    @test variable_bounds.upper == [T(4), noub]
    ey = MOI.add_constraint(variable_bounds, y, MOI.EqualTo(T(-3)))
    @test variable_bounds.lower == [T(3), T(-3)]
    @test variable_bounds.upper == [T(4), T(-3)]
    MOI.delete(variable_bounds, ey)
    @test variable_bounds.lower == [T(3), nolb]
    @test variable_bounds.upper == [T(4), noub]
    MOI.delete(variable_bounds, cx)
    @test variable_bounds.lower == [nolb, nolb]
    @test variable_bounds.upper == [noub, noub]
end

function test_bound_vectors()
    _test_bound_vectors(Int, 0, 0)
    _test_bound_vectors(Float64, -Inf, Inf)
    return
end

function test_ListOfConstraintTypesPresent_2()
    for set in (
        MOI.EqualTo(1.0),
        MOI.GreaterThan(1.0),
        MOI.LessThan(1.0),
        MOI.Interval(1.0, 2.0),
        MOI.Semicontinuous(1.0, 2.0),
        MOI.Semiinteger(1.0, 2.0),
        MOI.Integer(),
        MOI.ZeroOne(),
    )
        model = MOI.Utilities.VariablesContainer{Float64}()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, x, set)
        @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) ==
              [(MOI.VariableIndex, typeof(set))]
    end
    return
end

end  # module

TestVariableContainer.runtests()


text:
import Base: show

"""

# CmdStanModels

Base model specialized in:

### Method
```julia
*  SampleModel                         : StanSample.jl
*  OptimizeModel                       : StanOptimize.jl
*  VariationalModel                     : StanVariational.jl

Not yet done:

*  DiagnoseModel                       : StanDiagnose.jl
*  Generate_QuamtitiesModel            : StanGenerate_Quantities.jl
```
"""
abstract type CmdStanModels end

"""

# RandomSeed

Random number generator seed value

### Method
```julia
RandomSeed(;seed=-1)
```
### Optional arguments
```julia
* `seed::Int`           : Starting seed value
```
"""
mutable struct RandomSeed
  seed::Int64
end
RandomSeed(;seed::Number=-1) = RandomSeed(seed)

"""

# Init

Default bound for parameter initial value interval (if not found in init file)

### Method
```julia
Init(;bound=2)
```
### Optional arguments
```julia
* `bound::Number`           : Set interval to [-bound, bound]
```
"""
mutable struct Init
  bound::Int64
end
Init(;bound::Int64=2) = Init(bound)

"""

# Output

Default settings for output portion in cmdstan command

### Method
```julia
Output(;file="", diagnostics_file="", refresh=100)
```
### Optional arguments
```julia
* `fikle::AbstractString`              : File used to store draws
* `diagnostics_fikle::AbstractString`  : File used to store diagnostics
* `refresh::Int64`                     : Output refresh rate
```
"""
mutable struct Output
  file::String
  diagnostic_file::String
  refresh::Int64
end
Output(;file="", diagnostic_file="", refresh=100) =
  Output(file, diagnostic_file, refresh)


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

using MortarContact2DAD, Test
import FEMBase.Test: @test_resource

X = Dict(1 => [-2.0, 0.0], 2 => [ 0.0, 0.0], 3 => [ 2.0, 0.0],
         4 => [-2.0, 0.0], 5 => [ 0.0, 0.0], 6 => [-2.0, -2.0],
         7 => [0.0, -2.0], 8 => [2.0, -2.0], 9 => [-2.0, 2.0],
         10 => [0.0, 2.0])

 D = [0.5, 0.0]

for j in [4, 5, 9, 10]
    X[j] += D
end

contact = Problem(Contact2DAD, "contact", 2, "displacement")
element7 = Element(Seg2, [1, 2])
element8 = Element(Seg2, [2, 3])
element9 = Element(Seg2, [4, 5])
contact_slave_elements = [element7, element8]
contact_master_elements = [element9]
add_slave_elements!(contact, contact_slave_elements)
add_master_elements!(contact, contact_master_elements)
all_elements = [element7, element8, element9]
update!(all_elements, "geometry", X)

function load1(filename)
    prefix = first(splitext(filename))
    objs = ("K", "C1", "C2", "D", "f", "g", "u", "la")
    data = Dict(obj => readdlm("$(prefix)_$(obj).dat") for obj in objs)
    return data
end

i1 = load1(@test_resource("iter_1.jld"))
i2 = load1(@test_resource("iter_2.jld"))
i3 = load1(@test_resource("iter_3.jld"))

function to_dict(u, ndofs, nnodes)
    return Dict(j => [u[ndofs*(j-1)+i] for i=1:ndofs] for j=1:nnodes)
end

ndofs = 20
contact.properties.iteration = 1
for (data, time) in zip([i1, i2, i3], [0.0, 1.0, 2.0])
    println("Testing data for iteration / time $time ")
    empty!(contact.assembly)
    contact.assembly.u = vec(data["u"])
    contact.assembly.la = vec(data["la"])
    update!(contact, "displacement", time => to_dict(data["u"], 2, 5))
    update!(contact, "lambda", time => to_dict(data["la"], 2, 5))
    assemble!(contact, time)
    @test isapprox(Matrix(sparse(contact.assembly.K, ndofs, ndofs)), data["K"])
    @test isapprox(Matrix(sparse(contact.assembly.C1, ndofs, ndofs)), data["C1"])
    @test isapprox(Matrix(sparse(contact.assembly.C2, ndofs, ndofs)), data["C2"])
    @test isapprox(Matrix(sparse(contact.assembly.D, ndofs, ndofs)), data["D"])
    @test isapprox(Vector(sparse(contact.assembly.f, ndofs, 1)[:]), data["f"]; atol=1.0e-9)
    @test isapprox(Vector(sparse(contact.assembly.g, ndofs, 1)[:]), data["g"])
end


text:
#=
fld-read.jl
Jeff Fessler and David Hong
=#

export fld_open, fld_header, fld_read

"""
    head, is_external_file, fid = fld_open(file ; dir::String="", chat=false)

Read header data from AVS format `.fld` file.
Leaves file open for more reading.

in
- `file::String` file name, usually ending in `.fld`

option
- `dir::String` prepend file name with this directory; default ""
- `chat::Bool` verbose? default: `false`

out
- `head::String` array of header information
- `is_external_file::Bool` true if AVS external format
- `fid::IOstream`
"""
function fld_open(
    file::AbstractString ;
    dir::String = "",
    chat::Bool = false,
)

    file = joinpath(dir, file)
    fid = open(file)

    formfeed = '\f' # form feed
    newline = '\n'
    is_external_file = false

    # read header until we find the end of file or the 1st form feed
    header = ""
    while true
        # end of file means external file (or error)
        if eof(fid)
            if occursin("file=", header)
                is_external_file = true
                break # end of header file!
            else
                chat && @info(header)
                throw("end of file before form feeds?")
            end
        end

        inchar = read(fid, Char) # read one character

        # form feed means embedded file
        if inchar == formfeed
            eof(fid) && throw("end of file before 2nd form feed?")
            inchar = read(fid, Char)
            (inchar != formfeed) && throw("not two form feeds?")
            chat && @info("embedded data file")
            break
        end

        # otherwise append this character to header string
        header *= inchar
    end

    header = string_to_array(header) # convert to array
    chat && @info(header)

    return header, is_external_file, fid
end

"""
    head = fld_header(file::String ; dir::String="", chat=false)

Read header data from AVS format `.fld` file, then close file.

in
- `file::String` file name, usually ending in `.fld`

option
- `dir::String` prepend file name with this directory; default ""
- `chat::Bool` verbose? default: `false`

out
- `head::String` array of header information
"""
function fld_header(file::AbstractString ; kwargs...)
    header, _, fid = fld_open(file ; kwargs...)
    close(fid)
    return header
end

# todo:
# + [ ] 'raw'    0|1    1: return raw data class (default), 0: return doubles
# + [ ] 'slice' int    specify which slice to read from 3D file. (0 ... nz-1)
# + [ ] 'chat'    0|1    enable verbosity
# + [ ] 'dim_only' 0|1    returns dims. data and coord are equal to [].
# + [ ] 'coord' 0|1    returns coordinates too (default: 0)
# + [ ] 'coord_format'    default: 'n'; see fopen() for machine format options
#    (needed for some UM .fld files with 'vaxd' coordinates)
# + [ ] multi files
# + [ ] short datatype

"""
    data = fld_read(file::String ; dir::String="", chat=false)

Read data from AVS format `.fld` file

in
- `file` file name, usually ending in `.fld`

option
- `dir::String` prepend file name with this directory; default ""
- `chat::Bool` verbose?

out
- `data` Array (1D - 5D) in the data type of the file itself
"""
function fld_read(
    file::AbstractString ;
    dir::String = "",
    chat::Bool = false,
)

    file = joinpath(dir, file)

    header, is_external_file, fid = fld_open(file ; chat)
    chat && @info("is_external_file = $is_external_file")

    # parse header to determine data dimensions and type
    ndim = arg_get(header, "ndim")
    dims = [arg_get(header, "dim$ii") for ii in 1:ndim]
    fieldtype = arg_get(header, "field", false)
    datatype = arg_get(header, "data", false)
    (arg_get(header, "veclen") != 1) && throw("only veclen=1 done")
    chat && @info("ndim=$ndim")
    chat && @info("dims=$dims")

    # external file (binary data in another file)
    _skip = 0
    if is_external_file
        close(fid)
        extfile = arg_get(header, "file", false)

        filetype = arg_get(header, "filetype", false)
        chat && @info("Current file = '$file', External file = '$extfile', type='$filetype'")

        _skip = occursin("skip=",prod(header)) ?

        arg_get([prod(header)], "skip", false) : 0

        if filetype == "ascii"
            tmp = dirname(file)
            extfile = joinpath(tmp, extfile)
            chat && @info("extfile = $extfile")
            isfile(extfile) || throw("no ascii file $extfile")
            format, _, _ = datatype_fld_to_mat(datatype)
            return fld_read_ascii(extfile, (dims...,), format)

        elseif filetype != "multi"
            if !isfile(extfile)
                fdir = file
                slash = findlast(isequal('/'), fdir) # todo: windows?
                isnothing(slash) && throw("cannot find external file $extfile")
                fdir = fdir[1:slash]
                extfile = fdir * extfile # add directory
                !isfile(extfile) && throw("no external ref file $extfile")
            end

        else
            throw("multi not supported yet") # todo
        end

    else
        filetype = ""
        extfile = ""
    end

    # finally, read the binary data
    format, endian, bytes = datatype_fld_to_mat(datatype)

    # single file reading
    data = fld_read_single(file, fid, dims, datatype, fieldtype,
        is_external_file, extfile, format, endian, bytes, _skip)

    close(fid)

    return data
end

# todo: currently supports only one entry per line (see fld_read.m)
function fld_read_ascii(extfile::AbstractString, dims::Dims, datatype::DataType)
    data = zeros(datatype, dims)
    open(extfile, "r") do fid
        for i in 1:length(data)
            tmp = readline(fid)
            data[i] = parse(Float64, tmp)
        end
    end
    return data
end

function fld_read_single(
    file, fid, dims, datatype, fieldtype,
    is_external_file, extfile, format::DataType, endian, bytes, _skip,
)

    # reopen file to same position, with appropriate endian too.
    if is_external_file
        fid = open(extfile)
    end

    skip(fid,_skip)

    rdims = dims # from handling slice option

    # read binary data and reshape appropriately
    data = Array{format}(undef, rdims...)
    try
        read!(fid, data)
    catch
        @info("rdims=$rdims")
        throw("file count != data count")
    end

    return endian === :le ? htol.(data) :
        endian === :be ? hton.(data) :
        throw("bug $endian")
end

"""
    header = string_to_array(header_lines)

Convert long string with embedded newlines into string array.
"""
function string_to_array(header_lines::String)
    newline = '\n'

    # ensure there is a newline at end, since dumb editors can forget...
    if header_lines[end] != newline
        header_lines *= newline
    end

    ends = findall(isequal(newline), header_lines)
    (length(ends) <= 0) && throw("no newlines?")

    header = split(header_lines, newline, keepempty=false)

    # strip comments (lines that begin with #)
    header = filter(line -> line[1] != '#', header)

    return header
end

"""
    arg_get(head, name, toint)

Parse an argument from header, of the `name=value` form
"""
function arg_get(head::Array{<:AbstractString}, name::String, toint::Bool=true)
    for line in head
        start = findfirst(name * '=',line)
        if !isnothing(start)
            !isnothing(findnext(name*'=',line,start[end]+1)) && throw("bug: multiples?")
            line = line[(start[end]+1):end]
            arg = split(line)[1]
            toint && (arg = parse(Int,arg))
            return arg
        end
    end
    throw("could not find $name in header")
end

"""
    format, endian, bytes = datatype_fld_to_mat(datatype)

Determine data format from `.fld` header datatype.
"""
function datatype_fld_to_mat(datatype::AbstractString)

    dict = Dict([
        ("byte", (UInt8, :be, 1)), # :be irrelevant
        ("short_be", (UInt16, :be, 2)),
        ("short_sun", (UInt16, :be, 2)),
        ("xdr_short", (UInt16, :be, 2)),
        ("short_le", (UInt16, :le, 2)),
        ("int", (Int32, "", 4)), # native int - not portable
        ("int_le", (Int32, :le, 4)),
        ("int_be", (Int32, :be, 4)),
        ("xdr_int", (Int32, :be, 4)),
        ("float", (Float32, "", 4)), # native float - not portable
        ("float_le", (Float32, :le, 4)),
        ("float_be", (Float32, :be, 4)),
        ("xdr_float", (Float32, :be, 4)),
        ("double", (Float64, "", 8)), # native 64oat - not portable
        ("double_le", (Float64, :le, 8)),
        ("double_be", (Float64, :be, 8)),
        ("xdr_double", (Float64, :be, 8)),
    ])

    return dict[datatype] # format, endian, bytes
end

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO