text:
"""
    MeshExportModule

Module for export of meshes and data defined on meshes.
"""
module MeshExportModule

module VTK
################################################################################
# VTK export
################################################################################
using ...FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ...FESetModule: AbstractFESet, FESetP1, FESetL2, FESetT3, FESetQ4, FESetT4, FESetH8, FESetQ8, FESetL3, FESetT6, FESetT10, FESetH20, connasarray
import ...FENodeSetModule: FENodeSet
import Base.close

using Printf
import LinearAlgebra: norm, cross

const P1=1
const L2=3
const T3=5
const Q4=9
const T4=10
const H8=12
const L3=21
const T6=22
const Q8=23
const T10=24
const H20=25

VTKtypemap = Dict{DataType, Int}(FESetP1=>P1, FESetL2=>L2, FESetT3=>T3,
    FESetQ4=>Q4, FESetT4=>T4, FESetH8=>H8, FESetQ8=>Q8, FESetL3=>L3, FESetT6=>T6,
    FESetT10=>T10, FESetH20=>H20)

numnodesmap = Dict{Int, Int}(P1=>1, L2=>2, T3=>3,
    Q4=>4, T4=>4, H8=>8, Q8=>8, L3=>3, T6=>6,
    T10=>10, H20=>20)

"""
    vtkexportmesh(theFile::String, fens::FENodeSet, fes::T; opts...) where {T<:AbstractFESet}

Export mesh to a VTK 1.0 file as an unstructured grid.

`opts` = keyword argument list, where
`scalars` = array of tuples, (name, data)
`vectors` = array of tuples, (name, data)

For the `scalars`: If `data` is a vector, that data is exported as a single field.
On the other hand, if it is an 2d array, each column is exported  as a separate field.

"""
function vtkexportmesh(theFile::String, fens::FENodeSet, fes::T;
    opts...) where {T<:AbstractFESet}
    Cell_type = get(()->error("Unknown VTK type!"), VTKtypemap, typeof(fes));
    return vtkexportmesh(theFile, connasarray(fes), fens.xyz, Cell_type; opts...)
end

"""
    vtkexportmesh(theFile::String, Connectivity, Points, Cell_type;
        vectors=nothing, scalars=nothing)

Export mesh to a VTK 1.0 file as an unstructured grid.

`opts` = keyword argument list, where
`scalars` = array of tuples, (name, data)
`vectors` = array of tuples, (name, data)

For the `scalars`: If `data` is a vector, that data is exported as a single field.
On the other hand, if it is an 2d array, each column is exported  as a separate field.
"""
function vtkexportmesh(theFile::String, Connectivity, Points, Cell_type;
    vectors=nothing, scalars=nothing)
    X = Points
    if size(Points, 2) < 3
        X = zeros(size(Points,1),3)
        for j  = 1:size(Points,1)
            X[j,1:size(Points,2)] =  Points[j,:]
        end
    end
    numnodes = get(()->error("Wrong number of connected nodes!"), numnodesmap, Cell_type);
    if typeof(Connectivity[1]) <: Tuple # Vector of connectivity tuples: convert to an array
        c = fill(0, length(Connectivity), length(Connectivity[1]))
        for i = 1:length(Connectivity)
            c[i, :] = [Connectivity[i]...]
        end
        Connectivity = c
    end
    @assert numnodes == size(Connectivity,2)

    fid=open(theFile,"w");
    if (fid == -1)
        error(["Could not open " * theFile])
        return nothing
    end
    print(fid,"# vtk DataFile Version 1.0\n");
    print(fid,"FinEtools Export\n");
    print(fid,"ASCII\n");
    print(fid,"\n");
    print(fid,"DATASET UNSTRUCTURED_GRID\n");
    print(fid,"POINTS ", size(X,1), " double\n");
    for i= 1:size(X, 1)
        for j= 1:size(X,2)-1
            print(fid,X[i,j]," ");
        end
        print(fid,X[i,end],"\n");
    end
    print(fid,"\n");
    print(fid,"\n");

    print(fid,"CELLS ",size(Connectivity,1)," ",(size(Connectivity,1)*(size(Connectivity,2)+1)),"\n");
    for i= 1:size(Connectivity, 1)
        print(fid,size(Connectivity,2)," ");
        for j= 1:size(Connectivity,2)-1
            print(fid,Connectivity[i,j]-1," ");
        end
        print(fid,Connectivity[i,end]-1,"\n");
    end
    print(fid,"\n");
    print(fid,"\n");
    print(fid,"CELL_TYPES ",size(Connectivity,1),"\n");
    for i= 1:size(Connectivity,1)
        print(fid,Cell_type,"\n");
    end
    print(fid,"\n");
    print(fid,"\n");

    did_point_data = false

    # First try to write point data
    if (scalars != nothing)
        did_point_data = false
        for  sx = 1:length(scalars)
            name, data = scalars[sx]
            if (size(data, 1) == size(Points, 1)) # point data
                if (!did_point_data)
                    print(fid,"POINT_DATA ",size(data, 1),"\n");
                    did_point_data = true
                end
                if size(data, 2) > 1 # there are multiple scalar fields here
                    for i = 1:size(data, 2)
                        print(fid,"SCALARS ", name * "$(i)"," double\n");
                        print(fid,"LOOKUP_TABLE default\n");
                        for j= 1:size(data,  1)
                            print(fid,data[j, i],"\n");
                        end
                    end
                else  # there's just one scalar field here
                    print(fid,"SCALARS ", name," double\n");
                    print(fid,"LOOKUP_TABLE default\n");
                    for j= 1:size(data,  1)
                        print(fid,data[j],"\n");
                    end
                end
                print(fid,"\n");
            elseif (size(data, 1) == size(Connectivity, 1)) # cell data
                # Handled below, in the section for cell data
            end
        end
    end

    print(fid,"\n");

    if vectors != nothing
        for  vx = 1:length(vectors)
            name, data = vectors[vx]
            if (!did_point_data)
                print(fid,"POINT_DATA ",size(data, 1),"\n");
                did_point_data = true
            end
            print(fid,"VECTORS ", name," double\n");
            #print(fid,"LOOKUP_TABLE default\n");
            if size(data, 2) < 3
                X = zeros(size(data,1),3)
                for j  = 1:size(data,1)
                    X[j, 1:size(data,2)] =  data[j,:]
                end
            else
                X = data
            end
            for j= 1:size(X,1)
                k=1;
                print(fid,X[j,k]);
                for k = 2:size(X,2)
                    print(fid," ", X[j,k]);
                end
                print(fid,"\n");
            end
            print(fid,"\n");
        end
    end
    print(fid,"\n");

    # Are there any cell data?  If so, write out.
    if (scalars != nothing)
        did_cell_data = false
        for  sx = 1:length(scalars)
            name, data = scalars[sx]
            if (size(data, 1) == size(Points, 1)) # point data
                # Handled above for the case of point data
            elseif (size(data, 1) == size(Connectivity, 1)) # cell data
                if (!did_cell_data)
                    print(fid,"CELL_DATA ",size(data, 1),"\n");
                    did_cell_data = true
                end
                if size(data, 2) > 1 # there are multiple scalar fields here
                    for i = 1:size(data, 2)
                        print(fid,"SCALARS ", name * "$(i)"," double\n");
                        print(fid,"LOOKUP_TABLE default\n");
                        for j= 1:size(data,  1)
                            print(fid,data[j, i],"\n");
                        end
                    end
                else  # there's just one scalar field here
                    print(fid,"SCALARS ", name," double\n");
                    print(fid,"LOOKUP_TABLE default\n");
                    for j= 1:size(data,  1)
                        print(fid,data[j],"\n");
                    end
                end
                print(fid,"\n");
            end
        end
    end

    fid=close(fid);
    return true
end

"""
    vtkexportvectors(theFile::String, Points, vectors)

Export vector data to a VTK 1.0 file.

- `Points` = array of collection of coordinates (tuples or vectors),
- `vectors` = array of tuples, `(name, data)`, where `name` is a string, and
  `data` is array of collection of coordinates (tuples or vectors).

# Example
```
Points = [(1.0, 3.0), (0.0, -1.0)]
vectors = [("v", [(-1.0, -2.0), (1.0, 1.0)])]
vtkexportvectors("theFile.VTK", Points, vectors)
```
will produce file with
```
# vtk DataFile Version 1.0
FinEtools Export
ASCII

DATASET UNSTRUCTURED_GRID
POINTS 2 double
1.0 3.0 0.0
0.0 -1.0 0.0

POINT_DATA 2
VECTORS v double
-1.0 -2.0 0.0
1.0 1.0 0.0
```

!!! note

The filter "Glyph" must be used Paraview. Also in the drop-down
menu "Glyph mode" select "all points".
"""
function vtkexportvectors(theFile::String, Points, vectors)
    fid=open(theFile,"w");
    if (fid == -1)
        error(["Could not open " * theFile])
        return nothing
    end
    print(fid,"# vtk DataFile Version 1.0\n");
    print(fid,"FinEtools Export\n");
    print(fid,"ASCII\n");
    print(fid,"\n");
    print(fid,"DATASET UNSTRUCTURED_GRID\n");
    print(fid,"POINTS ", length(Points), " double\n");
    pt = fill(0.0, 3)
    for i = 1:length(Points)
        for k = 1:length(Points[i])
            pt[k] = Points[i][k]
        end
        for j= 1:length(pt)-1
            print(fid, pt[j]," ");
        end
        print(fid,pt[end],"\n");
    end
    print(fid,"\n");
    print(fid,"\n");

    did_point_data = false

    if vectors != nothing
        for  vx = 1:length(vectors)
            name, data = vectors[vx]
            if (!did_point_data)
                print(fid,"POINT_DATA ",length(data),"\n");
                did_point_data = true
            end
            print(fid,"VECTORS ", name," double\n");
            #print(fid,"LOOKUP_TABLE default\n");
            for i= 1:length(data)
                for k = 1:length(data[i])
                    pt[k] = data[i][k]
                end
                for j= 1:length(pt)-1
                   print(fid, pt[j]," ");
                end
                print(fid,pt[end],"\n");
            end
            print(fid,"\n");
        end
    end
    print(fid,"\n");

    fid=close(fid);
    return true
end

end # VTK

module Abaqus
################################################################################
# Abaqus export
################################################################################
using ...FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ...FESetModule: AbstractFESet, FESetP1, FESetL2, FESetT3, FESetQ4, FESetT4, FESetH8, FESetQ8, FESetL3, FESetT6, FESetT10, FESetH20, connasarray
import ...FENodeSetModule: FENodeSet
import Base.close
using LinearAlgebra

"""
    AbaqusExporter

Export mesh to Abaqus.
"""
mutable struct AbaqusExporter
  filename::AbstractString
  ios::IO
  element_range::Tuple{Int64, Int64}
  function AbaqusExporter(filename::AbstractString)
    if match(r".*\.inp$", filename) == nothing
      filename = filename * ".inp"
    end
    ios = open(filename,"w+")
    return new(deepcopy(filename), ios, (typemax(Int64), 0))
  end
end

"""
    COMMENT(self::AbaqusExporter, Text::AbstractString)

Write out the `**` comment option.
"""
function COMMENT(self::AbaqusExporter, Text::AbstractString)
  println(self.ios, "**" * Text)
end

"""
    HEADING(self::AbaqusExporter, Text::AbstractString)

Write out the `*HEADING` option.
"""
function HEADING(self::AbaqusExporter, Text::AbstractString)
  println(self.ios, "*HEADING");
  println(self.ios, Text)
end

"""
    PART(self::AbaqusExporter, NAME::AbstractString)

Write out the `*PART` option.
"""
function PART(self::AbaqusExporter, NAME::AbstractString)
  println(self.ios, "*PART, NAME=" * NAME)
end

"""
    END_PART(self::AbaqusExporter)

Write out the `*END PART` option.
"""
function END_PART(self::AbaqusExporter)
  println(self.ios, "*END PART")
end

"""
    ASSEMBLY(self::AbaqusExporter, NAME::AbstractString)

Write out the `*ASSEMBLY` option.
"""
function ASSEMBLY(self::AbaqusExporter, NAME::AbstractString)
  println(self.ios, "*ASSEMBLY, NAME=" * NAME)
end

"""
    END_ASSEMBLY(self::AbaqusExporter)

Write out the `*END ASSEMBLY` option.
"""
function END_ASSEMBLY(self::AbaqusExporter)
  println(self.ios, "*END ASSEMBLY")
end

"""
    INSTANCE(self::AbaqusExporter, NAME::AbstractString, PART::AbstractString)

Write out the `*INSTANCE` option.
"""
function INSTANCE(self::AbaqusExporter, NAME::AbstractString, PART::AbstractString)
  println(self.ios, "*INSTANCE, NAME=" * NAME * ", PART=" * PART)
end

"""
    END_INSTANCE(self::AbaqusExporter)

Write out the `*END INSTANCE` option.
"""
function END_INSTANCE(self::AbaqusExporter)
  println(self.ios, "*END INSTANCE")
end

"""
    NODE(self::AbaqusExporter, xyz::AbstractArray{T, 2}) where {T}

Write out the `*NODE` option.

`xyz`=array of node coordinates
"""
function NODE(self::AbaqusExporter, xyz::AbstractArray{T, 2}) where {T}
  println(self.ios, "*NODE")
  for j = 1:size(xyz,1)
    print(self.ios, "$j,$(xyz[j, 1])")
    for ixxxx = 2:1:size(xyz,2)
      print(self.ios, ",$(xyz[j, ixxxx])")
    end
    print(self.ios, "\n")
  end
end

"""
    ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
      start::Integer, conn::AbstractArray{T, 2}) where {T<:Integer}

Write out the `*ELEMENT` option.

`TYPE`= element type code,
`ELSET`= element set to which the elements belong,
`start`= start the element numbering at this integer,
`conn`= connectivity array for the elements, one row per element
"""
function ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
    start::Integer, conn::AbstractArray{T, 2}) where {T<:Integer}
    # Check that start is valid
    @assert start > 0 "The  starting element number must be > 0"
    # Check that the current element number range is disjoint from the
    # range of the elements on input: they must not overlap
    @assert (self.element_range[2] < start) ||
            (self.element_range[1] > start+size(conn, 1)-1) "Elements  must be given unique numbers"
    # Update the element range
    self.element_range = (min(self.element_range[1], start),
                            max(self.element_range[2], start+size(conn, 1)-1))
    println(self.ios, "*ELEMENT, TYPE =" * TYPE * ", ELSET=" * ELSET)
    for j=1:size(conn,1)
        print(self.ios, "$(j+start-1),")
        for ixxxx = 1:size(conn,2)-1
        if ixxxx>15
            print(self.ios, "\n")
        end
        print(self.ios, "$(conn[j,ixxxx]),")
        end
        print(self.ios, "$(conn[j,size(conn,2)])\n")
    end
end

function ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
    conn::AbstractArray{T, 2}) where {T<:Integer}
    start = self.element_range[2] + 1
    ELEMENT(self, TYPE, ELSET, start, conn)
end

function ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
    start::Integer, conn::C) where {C}
    @assert typeof(conn[1])<:Tuple
    c = fill(0, length(conn), length(conn[1]))
    for i = 1:length(conn)
        c[i, :] = [conn[i]...]
    end
    ELEMENT(self, TYPE, ELSET, start, c)
end

function ELEMENT(self::AbaqusExporter, TYPE::AbstractString, ELSET::AbstractString,
    conn::C) where {C}
    start = self.element_range[2] + 1
    ELEMENT(self, TYPE, ELSET, start, c)
end

"""
    NSET_NSET(self::AbaqusExporter, NSET::AbstractString,
      n::AbstractVector{T}) where {T<:Integer}

Write out the `*NSET` option.

`NSET` = name of the set,
`n` = array of the node numbers
"""
function NSET_NSET(self::AbaqusExporter, NSET::AbstractString, n::AbstractVector{T}) where {T<:Integer}
  println(self.ios, "*NSET, NSET=" * NSET)
  for j=1:length(n)-1
    println(self.ios, "$(n[j]),")
  end
  println(self.ios, "$(n[length(n)])")
end

"""
    ELSET_ELSET(self::AbaqusExporter, ELSET::AbstractString, n::AbstractArray{T, 1}) where {T<:Integer}

Write out the `*ELSET` option.

`ELSET` = name of the set,
`n` = array of the node numbers
"""
function ELSET_ELSET(self::AbaqusExporter, ELSET::AbstractString, n::AbstractArray{T, 1}) where {T<:Integer}
  println(self.ios, "*ELSET, ELSET=" * ELSET)
  for j=1:length(n)-1
    println(self.ios, "$(n[j]),")
  end
  println(self.ios, "$(n[length(n)])")
end

"""
    ORIENTATION(self::AbaqusExporter, ORIENTATION::AbstractString,
      a::AbstractArray{T,1}, b::AbstractArray{T,1})

Write out the `*ORIENTATION` option.

Invoke at level: Part,  Part instance,  Assembly
"""
function ORIENTATION(self::AbaqusExporter, ORIENTATION::AbstractString,
  a::AbstractArray{T,1}, b::AbstractArray{T,1}) where {T<:Real}
  println(self.ios, "*ORIENTATION,NAME=" * ORIENTATION)
  println(self.ios, "$(a[1]),$(a[2]),$(a[3]),$(b[1]),$(b[2]),$(b[3])")
  println(self.ios, "1,0.0")
end

"""
    MATERIAL(self::AbaqusExporter, MATERIAL::AbstractString)

Write out the `*MATERIAL` option.
"""
function MATERIAL(self::AbaqusExporter, MATERIAL::AbstractString)
  println(self.ios, "*MATERIAL,NAME=" * MATERIAL);
end

"""
    ELASTIC(self::AbaqusExporter, E1::F, E2::F, E3::F, nu12::F, nu13::F, nu23::F,
      G12::F, G13::F, G23::F) where {F}

Write out the `*ELASTIC,TYPE=ENGINEERING CONSTANTS` option.
"""
function ELASTIC(self::AbaqusExporter, E1::F, E2::F, E3::F, nu12::F, nu13::F, nu23::F,
  G12::F, G13::F, G23::F) where {F}
  println(self.ios, "*ELASTIC,TYPE=ENGINEERING CONSTANTS");
  println(self.ios, "$E1,$E2,$E3,$nu12,$nu13,$nu23,$G12,$G13,");
  println(self.ios, "$G23");
end

"""
    ELASTIC(self::AbaqusExporter, E::F, nu::F) where {F}

Write out the `*ELASTIC,TYPE=ISOTROPIC` option.
"""
function ELASTIC(self::AbaqusExporter, E::F, nu::F) where {F}
  println(self.ios, "*ELASTIC,TYPE=ISOTROPIC ");
  println(self.ios, "$E,$nu");
end

"""
    EXPANSION(self::AbaqusExporter, CTE::F) where {F}

Write out the `*EXPANSION` option.
"""
function EXPANSION(self::AbaqusExporter, CTE::F) where {F}
  println(self.ios, "*EXPANSION");
  println(self.ios, "$CTE,");
end

"""
    DENSITY(self::AbaqusExporter, rho)

Write out the `*DENSITY` option.
"""
function DENSITY(self::AbaqusExporter, rho)
  println(self.ios, "*DENSITY ");
  println(self.ios, "$rho");
end

"""
    SECTION_CONTROLS(self::AbaqusExporter, NAME::AbstractString,
      OPTIONAL::AbstractString)

Write out the `*SECTION CONTROLS` option.

`OPTIONAL` = string, for instance
                     HOURGLASS=ENHANCED
"""
function SECTION_CONTROLS(self::AbaqusExporter, NAME::AbstractString,
  OPTIONAL::AbstractString)
  println(self.ios, "*SECTION CONTROLS, NAME=" * NAME * "," * OPTIONAL);
end

"""
    SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
      ORIENTATION::AbstractString, ELSET::AbstractString,
      CONTROLS::AbstractString)

Write out the `*SOLID SECTION` option.

Level: Part,  Part instance

"""
function SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString,
  CONTROLS::AbstractString)
  println(self.ios, "*SOLID SECTION,MATERIAL=" * MATERIAL *
      ",ORIENTATION =" * ORIENTATION * ",ELSET=" * ELSET *
      ", CONTROLS  =" * CONTROLS);
end

"""
    SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
      ORIENTATION::AbstractString, ELSET::AbstractString)

Write out the `*SOLID SECTION` option.

Level: Part,  Part instance

"""
function SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString)
  println(self.ios, "*SOLID SECTION,MATERIAL=" * MATERIAL *
      ",ORIENTATION =" * ORIENTATION * ",ELSET=" * ELSET);
end

"""
    SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
      ORIENTATION::AbstractString, ELSET::AbstractString)

Write out the `*SOLID SECTION` option.

Level: Part,  Part instance

"""
function SOLID_SECTION(self::AbaqusExporter, MATERIAL::AbstractString,
  ORIENTATION::AbstractString, ELSET::AbstractString, thickness::F) where {F}
  SOLID_SECTION(self, MATERIAL, ORIENTATION, ELSET)
  println(self.ios, "$(thickness),");
end

# """
#     HOURGLASS(self::AbaqusExporter, KIND::AbstractString, VALUE::F) where {F}
#
# Write out the `*HOURGLASS` option.
#
# ```
# Example:
# *SOLID SECTION,ELSET=SOLID3,MATERIAL=MAT,CONTROL=Ac
# *HOURGLASS STIFFNESS
# 5.E8
# ```
# """
# function HOURGLASS(self::AbaqusExporter, KIND::AbstractString, VALUE::F) where {F}
#   println(self.ios, "*HOURGLASS " * KIND);
#   println(self.ios, "$( VALUE )");
# end

"""
    SURFACE_SECTION(self::AbaqusExporter, ELSET::AbstractString)

Write out the `*SURFACE SECTION` option.
"""
function SURFACE_SECTION(self::AbaqusExporter, ELSET::AbstractString)
  println(self.ios, "*SURFACE SECTION, ELSET=" * ELSET);
end

"""
    STEP_PERTURBATION_STATIC(self::AbaqusExporter)

Write out the `*STEP,PERTURBATION` option for linear static analysis.
"""
function STEP_PERTURBATION_STATIC(self::AbaqusExporter)
  println(self.ios, "*STEP,PERTURBATION");
  println(self.ios, "*STATIC");
end

"""
    STEP_PERTURBATION_BUCKLE(self::AbaqusExporter, neigv::Integer)

Write out the `*STEP,PERTURBATION` option for linear buckling analysis.
"""
function STEP_PERTURBATION_BUCKLE(self::AbaqusExporter, neigv::Integer)
  println(self.ios, "*STEP, name=Buckling, nlgeom=NO, perturbation");
  println(self.ios, "*BUCKLE");
  println(self.ios, "$(neigv), , , , ");
end

"""
    STEP_FREQUENCY(self::AbaqusExporter, Nmodes::Integer)

Write out the `*STEP,FREQUENCY` option.
"""
function STEP_FREQUENCY(self::AbaqusExporter, Nmodes::Integer)
  println(self.ios, "*STEP");
  println(self.ios, "*FREQUENCY, EIGENSOLVER=LANCZOS");
  println(self.ios, "$(Nmodes), , ,-1.E6 \n");
end

"""
    BOUNDARY(self::AbaqusExporter, mesh, nodes, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}

Write out the `*BOUNDARY` option.

The boundary condition is applied to the nodes specified in
the array `nodes`, in the mesh (or node set) `meshornset`.

`meshornset` = mesh or node set (default is empty)
`nodes` = array of node numbers, the node numbers are attached to the mesh label,
`is_fixed`= array of Boolean flags (true means fixed, or prescribed),  one row per node,
`fixed_value`=array of displacements to which the corresponding displacement components is fixed

# Example

```
BOUNDARY(AE, "ASSEM1.INSTNC1", 1:4, fill(true, 4, 1), reshape([uy(fens.xyz[i, :]...) for i in 1:4], 4, 1))
```
"""
function BOUNDARY(self::AbaqusExporter, meshornset, nodes, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}
	println(self.ios, "*BOUNDARY");
	if meshornset == ""
		meshlabel = ""
	else
		meshlabel = meshornset * "."
	end
	for j in nodes
		for k=1:size(is_fixed,2)
			#<node number>, <first dof>, <last dof>, <magnitude of displacement>
			if is_fixed[j,k]
				println(self.ios, "$(meshlabel)$(j),$k,$k,$(fixed_value[j,k])");
			end
		end
	end
end

function BOUNDARY(self::AbaqusExporter, nodes, is_fixed::AbstractVector{B},  fixed_value::AbstractVector{F}) where {B, F}
	BOUNDARY(self, nodes, reshape(is_fixed, length(is_fixed), 1), reshape(fixed_value, length(fixed_value), 1))
end

"""
    BOUNDARY(self::AbaqusExporter, meshornset, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}

Write out the `*BOUNDARY` option.

`meshornset` = name of a mesh or a node set,
`is_fixed`= array of Boolean flags (true means fixed, or prescribed),  one row per node, as many columns as there are degrees of freedom per node,
`fixed_value`=array of displacements to which the corresponding displacement components is fixed, as many columns as there are degrees of freedom per node
"""
function BOUNDARY(self::AbaqusExporter, meshornset, is_fixed::AbstractArray{B,2}, fixed_value::AbstractArray{F,2}) where {B, F}
	BOUNDARY(self, meshornset, 1:size(is_fixed, 1), is_fixed, fixed_value)
end

"""
    BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer,  fixed_value::FFlt)

Write out the `*BOUNDARY` option.

`NSET` = name of a node set,
`is_fixed`= array of Boolean flags (true means fixed, or prescribed),  one row per node,
`fixed_value`=array of displacements to which the corresponding displacement components is fixed
"""
function BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer,  fixed_value::FFlt)
	println(self.ios, "*BOUNDARY");
	println(self.ios, NSET * ",$dof,$dof,$(fixed_value)");
end

"""
    BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer)

Write out the `*BOUNDARY` option to fix displacements at zero for a node set.

Invoke at Level: Model,  Step

`NSET`= node set,
`dof`=Degree of freedom, 1, 2, 3
"""
function BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer)
	BOUNDARY(self, NSET, dof,  0.0)
end

"""
    BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer,
      value::F) where {F}

Write out the `*BOUNDARY` option to fix displacements at nonzero value for a
node set.

`NSET`= node set,
`dof`=Degree of freedom, 1, 2, 3
"""
function BOUNDARY(self::AbaqusExporter, NSET::AbstractString, dof::Integer, value::F) where {F}
	println(self.ios, "*BOUNDARY,TYPE=DISPLACEMENT");
	println(self.ios, NSET * ",$dof,$dof,$value");
end

"""
    DLOAD(self::AbaqusExporter, ELSET::AbstractString,
      traction::AbstractVector{F}) where {F}

Write out the `*DLOAD` option.
"""
function DLOAD(self::AbaqusExporter, ELSET::AbstractString, traction::AbstractVector{F}) where {F}
	println(self.ios, "*DLOAD, follower=NO");
	nt = norm(traction)
	print(self.ios, ELSET * ",TRVEC," * "$(nt),")
	println(self.ios, "$(traction[1]/nt),$(traction[2]/nt),$(traction[3]/nt)");
end

"""
    CLOAD(self::AbaqusExporter, NSET::AbstractString, dof::Integer,
      magnitude::F) where {F}

Write out the `*CLOAD` option.

NSET=Node set
dof= 1, 2, 3,
magnitude= signed multiplier
"""
function CLOAD(self::AbaqusExporter, NSET::AbstractString, dof::Integer,
  magnitude::F) where {F}
  println(self.ios, "*CLOAD");
  println(self.ios, NSET * ",$dof,$magnitude");
end

"""
    CLOAD(self::AbaqusExporter, nodenumber::Integer, dof::Integer,
      magnitude::F) where {F}

Write out the `*CLOAD` option.

nodenumber=Number of node
dof= 1, 2, 3,
magnitude= signed multiplier
"""
function CLOAD(self::AbaqusExporter, nodenumber::Integer, dof::Integer,
  magnitude::F) where {F}
  println(self.ios, "*CLOAD");
  println(self.ios, "$(nodenumber),$dof,$magnitude");
end

"""
    TEMPERATURE(self::AbaqusExporter, nlist::AbstractArray{I, 1},
      tlist::AbstractArray{F, 1}) where {I, F}

Write out the `*TEMPERATURE` option.
"""
function TEMPERATURE(self::AbaqusExporter, Classifier::AbstractString,
  nlist::AbstractArray{I, 1},
  tlist::AbstractArray{F, 1}) where {I, F}
  @assert length(nlist) == length(tlist)
  println(self.ios, "*TEMPERATURE");
  for ixxxx = 1:length(nlist)
    println(self.ios, Classifier * "$(nlist[ixxxx]),$(tlist[ixxxx])");
  end
end

"""
    NODE_PRINT(self::AbaqusExporter, NSET::AbstractString)

Write out the `*NODE PRINT` option.
"""
function NODE_PRINT(self::AbaqusExporter, NSET::AbstractString; KEYS = "U")
  println(self.ios, "*NODE PRINT, NSET=" * NSET);
  println(self.ios, KEYS);
end

"""
    EL_PRINT(self::AbaqusExporter, ELSET::AbstractString, KEYS::AbstractString)

Write out the `*EL PRINT` option.
"""
function EL_PRINT(self::AbaqusExporter, ELSET::AbstractString, KEYS::AbstractString)
  println(self.ios, "*EL PRINT, ELSET=" * ELSET * ", POSITION=INTEGRATION POINTS, SUMMARY= YES");
  println(self.ios, "");
  println(self.ios, KEYS);
end

"""
    ENERGY_PRINT(self::AbaqusExporter)

Write out the `*ENERGY PRINT` option.
"""
function ENERGY_PRINT(self::AbaqusExporter)
  println(self.ios, "*ENERGY PRINT");
end

"""
    END_STEP(self::AbaqusExporter)

Write out the `*END STEP` option.
"""
function END_STEP(self::AbaqusExporter)
  println(self.ios, "*END STEP");
end

"""
    close(self::AbaqusExporter)

Close  the stream opened by the exporter.
"""
function close(self::AbaqusExporter)
  close(self.ios)
end

end # Abaqus

module CSV
################################################################################
# CSV export
################################################################################

"""
    savecsv(name::String; kwargs...)

Save arrays as a CSV file.
"""
function savecsv(name::String; kwargs...)
    colnames = Symbol[]
    columns = Any[]
    for apair in pairs(kwargs)
        sy, val = apair
        push!(colnames, sy)
        push!(columns, val)
    end
    ncol = length(colnames)
    nrow = length(columns[1])
    for j = 1:ncol
        @assert length(columns[j]) == nrow "Columns must have the same number of rows"
    end
    if !occursin(r"^.*\.csv$", name) && !occursin(r"^.*\.CSV$", name)
        name = name * ".csv"
    end
    open(name, "w") do fid
        for (j, c) in enumerate(colnames)
            print(fid, "\"", c, "\"")
            j < ncol && print(fid, ",")
        end
        print(fid, "\n")
        for r = 1:nrow
            for j = 1:ncol
                print(fid, columns[j][r])
                j < ncol && print(fid, ",")
            end
            print(fid, "\n")
        end
        print(fid, "\n")
    end
    return true
end

end # CSV

module NASTRAN
################################################################################
# NASTRAN export
################################################################################
using ...FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ...FESetModule: AbstractFESet, FESetP1, FESetL2, FESetT3, FESetQ4, FESetT4, FESetH8, FESetQ8, FESetL3, FESetT6, FESetT10, FESetH20, connasarray
import ...FENodeSetModule: FENodeSet
import Base.close
using Printf

"""
NASTRANExporter

Export mesh to Abaqus.
"""
mutable struct NASTRANExporter
    filename::AbstractString
    ios::IO
    element_range::Tuple{Int64, Int64}
    function NASTRANExporter(filename::AbstractString)
        if match(r".*\.nas$", filename) == nothing
            filename = filename * ".nas"
        end
        ios = open(filename,"w+")
        return new(deepcopy(filename), ios, (typemax(Int64), 0))
    end
end

NASTRANtypemap = Dict{DataType, AbstractString}(FESetT4=>"CTETRA", FESetT10=>"CTETRA")

"""
    CEND(self::NASTRANExporter)

Terminate the Executive Control section.
"""
function CEND(self::NASTRANExporter)
    println(self.ios, "CEND")
end

"""
    BEGIN_BULK(self::NASTRANExporter)

Terminate the Case Control section by starting the bulk section.
"""
function BEGIN_BULK(self::NASTRANExporter)
    println(self.ios, "BEGIN BULK")
end

"""
    ENDDATA(self::NASTRANExporter)

Terminate the bulk section.
"""
function ENDDATA(self::NASTRANExporter)
    println(self.ios, "ENDDATA")
end

"""
    GRID(self::NASTRANExporter, n::Int, xyz::Vector{FFlt})

Write a grid-point statement.
"""
function GRID(self::NASTRANExporter, n::Int, xyz::Vector{FFlt})
    @printf self.ios "GRID,%d,,%g,%g,%g\n" n xyz[1] xyz[2] xyz[3]
end

"""
    PSOLID(self::NASTRANExporter, pid::Int, mid::Int)

Write solid-property statement.
"""
function PSOLID(self::NASTRANExporter, pid::Int, mid::Int)
    @printf self.ios "PSOLID,%d,%d\n" pid mid
end

"""
    MAT1(self::NASTRANExporter, mid::Int, E::FFlt, G::FFlt, nu::FFlt, rho::FFlt, A::FFlt, TREF::FFlt, GE::FFlt)

Write a statement for an isotropic elastic material.
"""
function MAT1(self::NASTRANExporter, mid::Int, E::FFlt, G::FFlt, nu::FFlt, rho::FFlt, A::FFlt, TREF::FFlt, GE::FFlt)
    @printf self.ios "MAT1,%d,%g,%g,%g,%g,%g,%g,%g\n" mid E G nu rho A TREF GE
end

"""
    MAT1(self::NASTRANExporter, mid::Int, E::FFlt, nu::FFlt, rho::FFlt = 0.0, A::FFlt = 0.0, TREF::FFlt = 0.0)

Write a statement for an isotropic elastic material.
"""
function MAT1(self::NASTRANExporter, mid::Int, E::FFlt, nu::FFlt, rho::FFlt = 0.0, A::FFlt = 0.0, TREF::FFlt = 0.0, GE::FFlt = 0.0)
    @printf self.ios "MAT1,%d,%g,,%g,%g,%g,%g,%g\n" mid E nu rho A TREF GE
end

function CTETRA(self::NASTRANExporter, eid::Int, pid::Int, conn::Vector{Int})
    nc = length(conn)
    @printf self.ios "CTETRA,%d,%d" eid pid
    for k = 1:nc
        @printf self.ios ",%d" conn[k]
        if k == 6
            @printf self.ios "\n"
        end
    end
    @printf self.ios "\n"
end

"""
    close(self::NASTRANExporter)

Close  the stream opened by the exporter.
"""
function close(self::NASTRANExporter)
    close(self.ios)
end

end # NASTRAN

module STL
################################################################################
# STL export
################################################################################
using ...FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
using Printf
import LinearAlgebra: norm, cross
import Base.close

"""
STLExporter

Export surface mesh as STL file.
"""
mutable struct STLExporter
    filename::AbstractString
    ios::IO
    element_range::Tuple{Int64, Int64}
    function STLExporter(filename::AbstractString)
        if match(r".*\.stl$", filename) == nothing
            filename = filename * ".stl"
        end
        ios = open(filename,"w+")
        return new(deepcopy(filename), ios, (typemax(Int64), 0))
    end
end

function solid(self::STLExporter, name::AbstractString = "thesolid")
    @printf self.ios "solid %s\n" name
end

function facet(self::STLExporter, v1::Vector{FFlt}, v2::Vector{FFlt}, v3::Vector{FFlt})
    V = v2 - v1
    W = v3 - v1
    normal = cross(V, W)
    normal = normal / norm(normal)
    @printf self.ios "facet normal %f %f %f\n" normal[1] normal[2] normal[3]
    @printf self.ios "    outer loop\n"
    @printf self.ios "        vertex %f %f %f\n" v1[1] v1[2] v1[3]
    @printf self.ios "        vertex %f %f %f\n" v2[1] v2[2] v2[3]
    @printf self.ios "        vertex %f %f %f\n" v3[1] v3[2] v3[3]
    @printf self.ios "    endloop\n"
    @printf self.ios "endfacet\n"
end

function endsolid(self::STLExporter, name::AbstractString = "thesolid")
    @printf self.ios "endsolid %s\n" name
end

"""
    close(self::STLExporter)

Close  the stream opened by the exporter.
"""
function close(self::STLExporter)
    close(self.ios)
end

end # STL

module H2Lib
################################################################################
# H2Lib triangular-surface export
################################################################################
using Printf
import Base.BitSet

mutable struct _HypFace
    n::Int # serial number of the hyperface
    hf::Array{Int} # numbers of vertices defining the hyperface
    a::Int # anchor node number
    o::BitSet # numbers of the other nodes on the hyperface
end

mutable struct _HypFaceContainer
    nv::Int # total number of vertices per hyperface
    d::Dict{Int, Array{_HypFace}} # dictionary
    n::Int # total number of hyperfaces in the container
end

makehypfacedict(nv) = _HypFaceContainer(nv, Dict{Int, Array{_HypFace}}(), 0);

function hypfacen(container, hypf)
    h = sort([i for i in hypf])
    anchor = h[1]; other = BitSet(h[2:end]);
    ca = get(container.d, anchor, _HypFace[]);
    for k = 1:length(ca)
        if (ca[k].a == anchor) && (ca[k].o == other)
            return ca[k].n # do have this hyper face: here is its serial number
        end
    end
    return 0 # do not have this hyper face
end

function inserthypface!(container, hypf)
    if hypfacen(container, hypf) == 0
        @assert length(hypf) == container.nv "Hyperface does not have the right number of vertices"
        h = sort([i for i in hypf])
        anchor = h[1]; other = BitSet(h[2:end]);
        ca = get(container.d, anchor, _HypFace[]);
        container.n = container.n + 1
        push!(ca, _HypFace(container.n, hypf, anchor, other));
        container.d[anchor] = ca
    end
    return container
end

function hypfacedicttoarray(container)
    a = fill(0, container.n, container.nv)
    for ca in values(container.d)
        for i = 1:length(ca)
            a[ca[i].n, :] = ca[i].hf
        end
    end
    return a
end

function h2libexporttri(theFile::String, Connectivity, Points)
    @assert size(Connectivity, 2) == 3 "Only triangles accepted"
    adjust = 1 # the library is C-based and expects indexes to be zero-based

    ed = makehypfacedict(2)
    for i = 1:size(Connectivity, 1)
        inserthypface!(ed, Connectivity[i, [1, 2]])
        inserthypface!(ed, Connectivity[i, [2, 3]])
        inserthypface!(ed, Connectivity[i, [3, 1]])
    end
    ea = hypfacedicttoarray(ed)

    fid=open(theFile,"w");
    if (fid == -1)
        error(["Could not open " * theFile])
        return nothing
    end

    print(fid, size(Points, 1), " ", size(ea, 1), " ", size(Connectivity, 1), "\n");

    for i= 1:size(Points, 1)
        for j= 1:size(Points,2)-1
            print(fid,Points[i,j]," ");
        end
        print(fid,Points[i,end],"\n");
    end

    for i= 1:size(ea, 1)
        for j= 1:size(ea,2)-1
            print(fid, ea[i,j]-adjust, " ");
        end
        print(fid, ea[i,end]-adjust, "\n");
    end

    for i= 1:size(Connectivity, 1)
        s = [hypfacen(ed, Connectivity[i, [2, 3]]) hypfacen(ed, Connectivity[i, [3, 1]]) hypfacen(ed, Connectivity[i, [1, 2]])]
        for j= 1:size(Connectivity,2)
            print(fid, Connectivity[i,j]-adjust, " ");
        end
        for j= 1:length(s)-1
            print(fid, s[j]-adjust, " ");
        end
        print(fid, s[end]-adjust, "\n");
    end

    fid=close(fid);
    return true
end

end # H2Lib

module MESH
################################################################################
# MESH export
################################################################################
using DelimitedFiles
using ...FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ...FESetModule: AbstractFESet, FESetP1, FESetL2, FESetT3, FESetQ4, FESetT4, FESetH8, FESetQ8, FESetL3, FESetT6, FESetT10, FESetH20, connasarray
import ...FENodeSetModule: FENodeSet

MESHtypemap = Dict{DataType, String}(FESetP1=>"P1", FESetL2=>"L2", FESetT3=>"T3",
    FESetQ4=>"Q4", FESetT4=>"T4", FESetH8=>"H8", FESetQ8=>"Q8", FESetL3=>"L3", FESetT6=>"T6",
    FESetT10=>"T10", FESetH20=>"H20")

function write_MESH(meshfile::String, fens::FENodeSet, fes::T) where {T<:AbstractFESet}
    meshfilebase, ext = splitext(meshfile)
    dinfo = [
    meshfilebase * "-xyz.dat",
    MESHtypemap[typeof(fes)],
    meshfilebase * "-conn.dat"]
    open(dinfo[1], "w") do file
        writedlm(file, fens.xyz, ' ')
    end
    open(dinfo[3], "w") do file
        writedlm(file, connasarray(fes), ' ')
    end
    open(meshfilebase * ".mesh", "w") do file
        writedlm(file, dinfo)
    end
    return true
end

end

end


text:
################################################################################
#  [2] MAIN FUNCTIONS
################################################################################
module GSComp
using Random, DataFrames, XLSX, LinearAlgebra, Statistics, StatsBase, Distributions
include("mainsup.jl")
function ComParison(dataset,p,h,Lτ,Uτ,nmodls, cut)
    quint = [0.16 0.50 0.84]
    nc = nfields(dataset)
    m  = size(dataset[1])[2]
    ΔIRF = Array{Float32,3}(undef,m,h,nc*nmodls)
    ΔIRFcoun = Array{Float32,3}(undef,2*m,h,nc)
    j = 1
    for y in dataset
        range = nmodls*(j-1)+1:nmodls*j
        cout = GScomparison(y,p,h,Lτ, Uτ, nmodls)
        ΔIRF[:, :, range] = cout
        Δmn = dropdims(mean(cout, dims=3),dims=3)
        Δmd = dropdims(median(cout, dims=3),dims=3)
        ΔIRFcoun[:,:,j] = [Δmn; Δmd]
        j+=1
    end
    g1 = ΔIRF[:, :, 1:cut*nmodls]
    g2 = ΔIRF[:, :, cut*nmodls+1:end]
    Δmn1 = dropdims(mean(g1, dims=3),dims=3)
    Δmd1 = dropdims(median(g1, dims=3),dims=3)
    Δqt1 = Qntls(g1, nmodls, quint, m, h)
    Δmn2 = dropdims(mean(g2, dims=3),dims=3)
    Δmd2 = dropdims(median(g2, dims=3),dims=3)
    Δqt2 = Qntls(g2, (nc-cut)*nmodls, quint, m, h)
    ecx = (mean = Δmn1, median = Δmd1, bands = Δqt1)
    dcx = (mean = Δmn2, median = Δmd2, bands = Δqt2)
    return ecx, dcx, ΔIRFcoun
end

function GScomparison(y, p, h, Lτ, Uτ, nmodls)
    nx = 3
    varGS = 1:3
    pos_tot = 2
    # --------------------------------------------------------
    # [2.1.1] VAR in reduce form
    m = size(y)[2]
    Y, X = VARData(y, p)
    B, Φ = OLSbetas(Y, X, m, p, Xblock = true, Nx = nx)
    # --------------------------------------------------------
    # [2.1.2] Residual of the reduce form
    E = Y - X * B
    SS = (E' * E)
    # --------------------------------------------------------
    # [2.1.3] Global Shock identification
    ΔIRF = Array{Float64,3}(undef, m, h, nmodls)
    for i = 1:nmodls
        Φ, Γ, Γcom = GSComSim(Y, X, B, SS, p, Lτ, Uτ, pos_tot, nx, varGS)
        #Φ, Γ, ξ_tot
        IRF1, f1 = irf_fevd(Φ, Γ, h, m)
        IRF2, f2 = irf_fevd(Φ, Γcom, h, m)
        ΔIRF[:, :, i] = IRF1[:, 1, :] - IRF2[:, 1, :]
    end
    return ΔIRF
end
function GSComSim(Y, X, B, SS, p, Lτ, Uτ, pos_tot, nx, varGS)
    # --------------------------------------------------------
    # [2.2.1] Drawing βₙ
    Xblock = true
    T, m = size(Y)
    σdist = InverseWishart(T, SS)
    b = vec(B)
    rc, cc = size(B)
    Σ = rand(σdist)
    XX = Symmetric(inv(X'X), :L)
    kr = kron(Σ, XX)
    disti = MvNormal(b, kr)
    β, Φ = BetaDraw(disti, rc, cc, m, p)
    # --------------------------------------------------------
    # [2.2.2] Cholesky Identification
    E1 = Y - X * β
    Σ1 = (E1' * E1) / (T - size(X, 2))
    C1 = cholesky(Σ1).L
    # --------------------------------------------------------
    # [2.2.3] Impulse-Response Function :
    aux_IRF = Array{Float64,3}(undef, m * p, m * p, Uτ)
    IRF = Array{Float64,3}(undef, m, m, Uτ)
    aux_IRF[:, :, 1] = Φ
    IRF[:, :, 1] = aux_IRF[1:m, 1:m, 1] * C1
    for i = 2:Uτ
        aux_IRF[:, :, i] = Φ * aux_IRF[:, :, i-1]
        IRF[:, :, i] = aux_IRF[1:m, 1:m, i] * C1
    end
    # --------------------------------------------------------
    # [2.2.4] New Augmented ToT as in Zeev
    Λtot = Λmatrix(IRF, pos_tot, m; Lτ = Lτ, Uτ = Uτ, Xblock = true, nx = nx)
    ξcom = eigen(Λtot).vectors[:, end:-1:1]
    ξcom = sign.(diag(ξcom))' .* ξcom
    Γcom = [C1[:, 1:nx] * ξcom C1[:, nx+1:end]]
    Γcom = sign.(diag(Γcom))' .* Γcom
    # --------------------------------------------------------
    # [2.2.4] Global Shocks Identification
    # [Step 1] Matrix Λ and their weights
    k = length(varGS)
    n = (Xblock * nx) + (1 - Xblock) * m
    Λ = zeros(n, n, k)
    λ = zeros(k)
    for i = 1:k
        Λ[:, :, i] = Λmatrix(IRF, varGS[i], m; Lτ = Lτ, Uτ = Uτ, Xblock = true, nx = nx)
        λ[i] = tr(Λ[:, :, i])
    end
    λ = prod(λ) ./ λ
    Ξ = dropdims(sum(reshape(λ, 1, 1, :) .* Λ, dims = 3), dims = 3)
    ξ = eigen(Ξ).vectors[:, end:-1:1]
    Γ = [C1[:, 1:nx] * ξ C1[:, nx+1:end]]
    Γ = sign.(diag(Γ))' .* Γ  # Normalization to improve identification
    return Φ, Γ, Γcom # Identification, New-augmented
end
end

ext: jl

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
"""
Create an enum and export the enum type and all values.
"""
macro _exported_enum(name, args...)
    return esc(quote
        @enum($name, $(args...))
        export $name
        $([:(export $arg) for arg in args]...)
    end)
end

"""
Create and export a GoogleSheets enum.  The enum values are prefixed to avoid collisions.

Also creates a function which returns the GoogleSheets string value.
"""
macro _gsheet_enum(type_name, prefix, args...)
    values = [ Symbol(prefix,arg) for arg in args ]
    map = zip(values,String.(args))

    return esc(quote
        @_exported_enum($type_name, $(values...))
        function gsheet_string(x::$(type_name))::String
            d = Dict(
                $([:($k => $v) for (k,v) in map]...)
            )

            return d[x]
        end
    end)
end

"""
An authorization scope for accessing Google resources.

- `AUTH_SCOPE_READONLY`:  Allows read-only access to the user's sheets and their properties.
- `AUTH_SCOPE_READWRITE`: Allows read/write access to the user's sheets and their properties.
"""
@_gsheet_enum AuthScope AUTH_SCOPE_ READONLY READWRITE

"""
A Google Sheets value type.

- `VALUE_TYPE_MIN`: minimum value in a range
- `VALUE_TYPE_MAX`: maximum value in a range
- `VALUE_TYPE_NUMBER`: exact number
- `VALUE_TYPE_PERCENT`: percentage of a range
- `VALUE_TYPE_PERCENTILE`: percentile of a range
"""
@_gsheet_enum ValueType VALUE_TYPE_ MIN MAX NUMBER PERCENT PERCENTILE

"""
A Google Sheets condition type.

- `CONDITION_TYPE_NUMBER_GREATER`: value must be greater than the condition value.
- `CONDITION_TYPE_NUMBER_GREATER_THAN_EQ`: value must be greater than or equal to the condition value.
- `CONDITION_TYPE_NUMBER_LESS`: value must be less than the condition value.
- `CONDITION_TYPE_NUMBER_LESS_THAN_EQ`: value must be less than or equal to the condition value.
- `CONDITION_TYPE_NUMBER_EQ`: value must be equal to the condition value.
- `CONDITION_TYPE_NUMBER_NOT_EQ`: value must not be equal to the condition value.
- `CONDITION_TYPE_NUMBER_BETWEEN`: value must be between the two condition values.
- `CONDITION_TYPE_NUMBER_NOT_BETWEEN`: value must not be between the two condition values.
- `CONDITION_TYPE_TEXT_CONTAINS`: value must contain the condition value.
- `CONDITION_TYPE_TEXT_NOT_CONTAINS`: value must not contain the condition value.
- `CONDITION_TYPE_TEXT_STARTS_WITH`: value must start with the condition value.
- `CONDITION_TYPE_TEXT_ENDS_WITH`: value must end with the condition value.
- `CONDITION_TYPE_TEXT_EQ`: value must equal the condition value.
- `CONDITION_TYPE_TEXT_NOT_EQ`: value must not equal the condition value.
- `CONDITION_TYPE_TEXT_IS_EMAIL`: value must be a valid email address.
- `CONDITION_TYPE_TEXT_IS_URL`: value must be a valid URL.
- `CONDITION_TYPE_DATE_EQ`: value must have the same date as the condition value.
- `CONDITION_TYPE_DATE_NOT_EQ`: value must not have the same date as the condition value.
- `CONDITION_TYPE_DATE_BEFORE`: value's date must be before the date of the condition value.
- `CONDITION_TYPE_DATE_AFTER`: value's date must be before the date of the condition value.
- `CONDITION_TYPE_DATE_ON_OR_BEFORE`: value's date must be on or before the date of the condition value.
- `CONDITION_TYPE_DATE_ON_OR_AFTER`: value's date must be on or after the date of the condition value.
- `CONDITION_TYPE_DATE_BETWEEN`: value's date must be between the dates of the condition values.
- `CONDITION_TYPE_DATE_NOT_BETWEEN`: value's date must not be between the dates of the condition values.
- `CONDITION_TYPE_DATE_IS_VALID`: value must be a date.
- `CONDITION_TYPE_ONE_OF_LIST`: value must be present in the condition values.
- `CONDITION_TYPE_BLANK`: value must be empty.
- `CONDITION_TYPE_NOT_BLANK`: value must not be empty.
- `CONDITION_TYPE_CUSTOM_FORMULA`: condition formula must evaluate to true.
- `CONDITION_TYPE_BOOLEAN`: value must be TRUE/FALSE.
- `CONDITION_TYPE_ONE_OF_RANGE`: value is present in the condition value's cell range.
"""
@_gsheet_enum ConditionType CONDITION_TYPE_ NUMBER_GREATER NUMBER_GREATER_THAN_EQ NUMBER_LESS NUMBER_LESS_THAN_EQ NUMBER_EQ NUMBER_NOT_EQ NUMBER_BETWEEN NUMBER_NOT_BETWEEN TEXT_CONTAINS TEXT_NOT_CONTAINS TEXT_STARTS_WITH TEXT_ENDS_WITH TEXT_EQ TEXT_NOT_EQ TEXT_IS_EMAIL TEXT_IS_URL DATE_EQ DATE_NOT_EQ DATE_BEFORE DATE_AFTER DATE_ON_OR_BEFORE DATE_ON_OR_AFTER DATE_BETWEEN DATE_NOT_BETWEEN DATE_IS_VALID ONE_OF_LIST BLANK NOT_BLANK CUSTOM_FORMULA BOOLEAN ONE_OF_RANGE

"""
A Google Sheets number format type.

- `NUMBER_FORMAT_TYPE_TEXT`: text formatting, e.g. `1000.12`
- `NUMBER_FORMAT_TYPE_NUMBER`: number formatting, e.g. `1,000.12`
- `NUMBER_FORMAT_TYPE_PERCENT`: percent formatting, e.g. `10.12%`
- `NUMBER_FORMAT_TYPE_CURRENCY`: currency formatting, e.g. `\u0024 1,000.12`
- `NUMBER_FORMAT_TYPE_DATE`: date formatting, e.g. `9/26/2008`
- `NUMBER_FORMAT_TYPE_TIME`: time formatting, e.g. `3:59:00 PM`
- `NUMBER_FORMAT_TYPE_DATE_TIME`: date+time formatting, e.g. `9/26/08 15:59:00`
- `NUMBER_FORMAT_TYPE_SCIENTIFIC`: scientific number formatting, e.g. `1.01E+03`
"""
@_gsheet_enum NumberFormatType NUMBER_FORMAT_TYPE_ TEXT NUMBER PERCENT CURRENCY DATE TIME DATE_TIME SCIENTIFIC


text:
struct SRS{R<:AbstractRNG} <: Imputor
    rng::R
end

# Docstring below uses julia-repl cause the rng may give different result on different
# versions of julia.
"""
    SRS(; rng=Random.GLOBAL_RNG)

Simple Random Sampling (SRS) imputation is a method for imputing both continuous and categorical
variables. Furthermore, it completes imputation while preserving the distributional
properties of the variables (e.g., mean, standard deviation).

The basic idea is that for a given variable, `x`, with missing data, we randomly draw
from the observed values of `x` to impute the missing elements. Since the random draws
from `x` for imputation are done in proportion to the frequency distribution of the values
in `x`, the univariate distributional properties are generally not impacted; this is true
for both categorical and continuous data.

# Keyword Arguments
* `rng::AbstractRNG`: A random number generator to use for observation selection

# Example
```julia-repl
julia> using Random; using Impute: SRS, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, SRS(; rng=MersenneTwister(1234)); dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  1.0  2.0  5.0
 1.1  2.2  3.3  3.3  5.5
```
"""
SRS(; rng=Random.GLOBAL_RNG) = SRS(rng)

function _impute!(data::AbstractArray{Union{T, Missing}}, imp::SRS) where T
    obs_values = collect(skipmissing(data))
    if !isempty(obs_values)
        for i in eachindex(data)
            if ismissing(data[i])
                data[i] = rand(imp.rng, obs_values)
            end
        end
    end

    return data
end


text:
using ReinforcementLearningZoo
using Test
using ReinforcementLearningBase
using ReinforcementLearningCore
using ReinforcementLearningEnvironments
using Flux
using Statistics
using Random
using StableRNGs
using OpenSpiel

function get_optimal_kuhn_policy(α = 0.2)
    TabularRandomPolicy(
        table = Dict(
            "0" => [1 - α, α],
            "0pb" => [1.0, 0.0],
            "1" => [1.0, 0.0],
            "1pb" => [2.0 / 3.0 - α, 1.0 / 3.0 + α],
            "2" => [1 - 3 * α, 3 * α],
            "2pb" => [0.0, 1.0],
            "0p" => [2.0 / 3.0, 1.0 / 3.0],
            "0b" => [1.0, 0.0],
            "1p" => [1.0, 0.0],
            "1b" => [2.0 / 3.0, 1.0 / 3.0],
            "2p" => [0.0, 1.0],
            "2b" => [0.0, 1.0],
        ),
    )
end

@testset "ReinforcementLearningZoo.jl" begin
    include("cfr/cfr.jl")
end


text:
using Test
using Random
using ClimaCoupler, Dates, Unitful

@testset "Coupler Interface" begin
    Random.seed!(26)

    coupler = CouplerState()

    data = rand(10, 10)
    date = DateTime(2021)
    coupler_add_field!(coupler, :test1, data, nothing, date, u"kg")
    coupler_add_field!(coupler, :test2, data, nothing, date, u"km/hr")
    coupler_add_field!(coupler, :test3, data, nothing, date, u"°C")

    @show coupler

    @testset "coupler_get" begin
        @test data === coupler_get(coupler, :test1, nothing, date, u"kg")
        # unit conversion
        @test data .* 1000 == coupler_get(coupler, :test1, nothing, date, u"g")
        @test data .* (5 / 18) == coupler_get(coupler, :test2, nothing, date, u"m/s")
        @test ustrip.(u"°F", data * u"°C") == coupler_get(coupler, :test3, nothing, date, u"°F")

        # key not in coupler dict
        @test_throws KeyError coupler_get(coupler, :idontexist, nothing, date, u"kg")
        # retreival at inconsistent datetime
        @test_throws ErrorException coupler_get(coupler, :test1, nothing, DateTime(2022), u"kg")
        # incompatible units
        @test_throws Unitful.DimensionError coupler_get(coupler, :test1, nothing, date, u"°C")
    end

    @testset "coupler_put!" begin
        newdata = rand(10, 10)
        newdate = DateTime(2022)
        coupler_put!(coupler, :test1, newdata, nothing, newdate, u"kg")

        @test newdata == coupler_get(coupler, :test1, nothing, newdate, u"kg")
        # coupler_put! is in-place; original data array has been modified
        @test data === coupler_get(coupler, :test1, nothing, newdate, u"kg")
        # unit conversion; data is converted from g -> kg when put in coupler
        coupler_put!(coupler, :test1, newdata, nothing, newdate, u"g")
        @test newdata ≈ 1000 * data

        # coupler_put! must be to a previously add_fielded field
        @test_throws KeyError coupler_put!(coupler, :idontexist, newdata, nothing, newdate, u"kg")
        # incoming data must match dimensions of add_fielded field
        @test_throws DimensionMismatch coupler_put!(coupler, :test1, rand(10, 5), nothing, newdate, u"kg")
        # incompatible units
        @test_throws Unitful.DimensionError coupler_put!(coupler, :test1, newdata, nothing, newdate, u"J/m")
        # coupler_put! updates coupler timestamp
        @test_throws ErrorException coupler_get(coupler, :test1, nothing, date, u"kg")
    end

    @testset "grid interpolation" begin
        # placeholder
    end
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine

max_forks_repo_name: CliMA/CouplerMachine


text:
# Autogenerated wrapper script for LibGit2_jll for armv6l-linux-musleabihf
export libgit2

using MbedTLS_jll
using LibSSH2_jll
JLLWrappers.@generate_wrapper_header("LibGit2")
JLLWrappers.@declare_library_product(libgit2, "libgit2.so.1.1")
function __init__()
    JLLWrappers.@generate_init_header(MbedTLS_jll, LibSSH2_jll)
    JLLWrappers.@init_library_product(
        libgit2,
        "lib/libgit2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module AlgebraHelpersTest

using Test
using ReactiveMP
using Random

using LinearAlgebra

@testset "Helpers" begin

    @testset "mirrorlog" begin
        rng  = MersenneTwister(1234)
        vals = rand(rng, 10)
        @test ReactiveMP.mirrorlog.(vals) == map(x -> log(1 - x), vals)
    end

    @testset "xtlog" begin
        rng  = MersenneTwister(1234)
        vals = rand(rng, 10)
        @test ReactiveMP.xtlog.(vals) == map(x -> x * log(x), vals)
    end

    @testset "negate_inplace!" begin
        rng = MersenneTwister(1234)

        for size in 2:4
            A = rand(rng, size, size)
            B = copy(A)
            @test B == A
            @test ReactiveMP.negate_inplace!(B) == -A
            @test B == -A

            a = rand(rng)
            b = copy(a)
            @test ReactiveMP.negate_inplace!(b) == -a
            @test b == a # numbers stay unchanged
        end
    end

    @testset "negate_inplace!" begin
        rng = MersenneTwister(1234)

        for size in 2:4
            A = rand(rng, size, size)
            B = copy(A)
            c = randn(rng)
            @test B == A
            @test ReactiveMP.mul_inplace!(c, B) == c * A
            @test B == c * A

            a = rand(rng)
            b = copy(a)
            c = randn(rng)
            @test ReactiveMP.mul_inplace!(c, b) == c * a
            @test b == a # numbers stay unchanged
        end
    end

    @testset "rank1update" begin
        rng = MersenneTwister(1234)

        for size in 2:4, T1 in (Float32, Float64), T2 in (Float32, Float64), T3 in (Float32, Float64)
            A = rand(rng, T1, size, size)
            x = rand(rng, T2, size)
            @test ReactiveMP.rank1update(A, x) ≈ (A + x*x')
            y = rand(rng, T3, size)
            @test ReactiveMP.rank1update(A, x, y) ≈ (A + x*y')
        end

    end

    @testset "mul_trace" begin
        rng = MersenneTwister(1234)

        for size in 2:4, T1 in (Float32, Float64), T2 in (Float32, Float64)
            A = rand(rng, T1, size, size)
            B = rand(rng, T2, size, size)

            # See: https://github.com/JuliaSIMD/LoopVectorization.jl/issues/377
            # Remove if statement once fixed
            if T1 !== Float32 || T2 !== Float32
                @test ReactiveMP.mul_trace(A, B) ≈ tr(A * B)
            end

            a = rand(rng, T1)
            b = rand(rng, T2)
            @test ReactiveMP.mul_trace(a, b) ≈ a * b
        end

    end

end

end


text:
# Autogenerated wrapper script for gperftools_jll for aarch64-apple-darwin
export libprofiler, libtcmalloc, libtcmalloc_and_profiler, libtcmalloc_debug, libtcmalloc_minimal, libtcmalloc_minimal_debug

using LibUnwind_jll
JLLWrappers.@generate_wrapper_header("gperftools")
JLLWrappers.@declare_library_product(libprofiler, "@rpath/libprofiler.0.dylib")
JLLWrappers.@declare_library_product(libtcmalloc, "@rpath/libtcmalloc.4.dylib")
JLLWrappers.@declare_library_product(libtcmalloc_and_profiler, "@rpath/libtcmalloc_and_profiler.4.dylib")
JLLWrappers.@declare_library_product(libtcmalloc_debug, "@rpath/libtcmalloc_debug.4.dylib")
JLLWrappers.@declare_library_product(libtcmalloc_minimal, "@rpath/libtcmalloc_minimal.4.dylib")
JLLWrappers.@declare_library_product(libtcmalloc_minimal_debug, "@rpath/libtcmalloc_minimal_debug.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header(LibUnwind_jll)
    JLLWrappers.@init_library_product(
        libprofiler,
        "lib/libprofiler.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc,
        "lib/libtcmalloc.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_and_profiler,
        "lib/libtcmalloc_and_profiler.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_debug,
        "lib/libtcmalloc_debug.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_minimal,
        "lib/libtcmalloc_minimal.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libtcmalloc_minimal_debug,
        "lib/libtcmalloc_minimal_debug.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ITK
using Test
using TestImages

# Test if we can load cxx
@testset "Sanity Check" begin
    @test ITK.verifycxx(0) == 0
    @test ITK.verifycxx(100) == 100
end

# Simple registration check on still image
@testset "Registration Control Check" begin
    testimage_path = string(joinpath(dirname(dirname(pathof(TestImages)))), "/images")
    lighthouse_path = joinpath(testimage_path, "lighthouse.png")

    testreg1 = ITK.registerframe(lighthouse_path, lighthouse_path, "NA", false, "Amoeba")
    testreg2 = ITK.MeanSquaresGradientRegistration(lighthouse_path, lighthouse_path, "NA", false, 5.0, 0.005, 30, 0.5)

    @test (testreg1[1], testreg1[2]) == (0, 0)
    @test (testreg2[1], testreg2[2]) == (0, 0)
end

# TODO registration test case on moving images


text:
using CSFML.LibCSFML
using ColorTypes
using ModernGL

# fully discrete, stochastic agent-based simulation for the SIR infection model.
# https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology

struct Coordinate
	x::Int64
	y::Int64
end

function Base.:+(a::Coordinate, b::Coordinate)
	Coordinate(a.x+b.x,a.y+b.y)
end

possible_moves = [
	Coordinate( 1, 0),
	Coordinate( 0, 1),
 	Coordinate(-1, 0),
 	Coordinate( 0,-1),
]

function collide_boundary(c::Coordinate, L::Number)
	cx = c.x
	cy = c.y

	if cx < -L
		cx = -L
	end
	if cy < -L
			cy = -L
	end
	if cx > L
		cx = L
	end
	if cy > L
		cy = L
	end
	Coordinate(cx,cy)
end

@enum InfectionStatus S I R

mutable struct Agent
	position::Coordinate
	status::InfectionStatus
	num_infected::Int64
end

function initialize(N::Number, L::Number)
	agents=[]
	for _ in 1:N
		c = Coordinate(rand(-L:L), rand(-L:L) )
		a = Agent(c,S::InfectionStatus,0)
		push!(agents,a)
	end

	agents[rand(1:N)].status = I::InfectionStatus
	return agents
end

function getColor(s::InfectionStatus)
    color_S = RGB(0.4, 0.5, 1.0)
    color_I = RGB(1.0, 0.4 ,0.4 )
    color_R = RGB(0.4, 1.0, 0.4)

    if s == S
        color_S
	elseif s == I
		color_I
	else
		color_R
    end
end

position(a::Agent) = a.position
getColor(a::Agent) = getColor(a.status)

abstract type AbstractInfection end

struct CollisionInfectionRecovery <: AbstractInfection
	p_infection::Float64
	p_recovery::Float64
end

function interact!(agent::Agent, source::Agent, infection::CollisionInfectionRecovery)
	if agent.position == source.position
		if source.status == I && agent.status == S
		# Chance of infection
			if rand() < infection.p_infection
				agent.status = I
				source.num_infected += 1
			end
		end
		if agent.position == source.position
			if agent.status == I
			# Chance of recovery
				if rand() < infection.p_recovery
					agent.status = R
				end
			end
		end
	end
end

function step!(agents::Vector, L::Number, infection::AbstractInfection)
	a = rand(agents)
	a.position = collide_boundary(a.position+rand(possible_moves), L)
	#interact
    for other in agents
		if other !== a
			interact!(a,other,infection)
		end
	end
	agents
end

function sweep!(agents::Vector, L::Number, infection::AbstractInfection)
    for _ in 1:length(agents)
        step!(agents, L, infection)
    end
	agents
end

#===================================================================================================

                        ▒█████        ██▓      ▒█████        ██▓      ▒█████
                        ▒██▒  ██▒     ▓██▒     ▒██▒  ██▒     ▓██▒     ▒██▒  ██▒
                        ▒██░  ██▒     ▒██▒     ▒██░  ██▒     ▒██▒     ▒██░  ██▒
                        ▒██   ██░     ░██░     ▒██   ██░     ░██░     ▒██   ██░
                    ██▓ ░ ████▓▒░ ██▓ ░██░ ██▓ ░ ████▓▒░ ██▓ ░██░ ██▓ ░ ████▓▒░ ██▓
                    ▒▓▒ ░ ▒░▒░▒░  ▒▓▒ ░▓   ▒▓▒ ░ ▒░▒░▒░  ▒▓▒ ░▓   ▒▓▒ ░ ▒░▒░▒░  ▒▓▒
                    ░▒    ░ ▒ ▒░  ░▒   ▒ ░ ░▒    ░ ▒ ▒░  ░▒   ▒ ░ ░▒    ░ ▒ ▒░  ░▒
                    ░   ░ ░ ░ ▒   ░    ▒ ░ ░   ░ ░ ░ ▒   ░    ▒ ░ ░   ░ ░ ░ ▒   ░
                    ░      ░ ░    ░   ░    ░      ░ ░    ░   ░    ░      ░ ░    ░
                    ░             ░        ░             ░        ░             ░
===================================================================================================#

# init simulation

number_of_agents = 600
L=100

agents = initialize(number_of_agents, L)
pandemic = CollisionInfectionRecovery( 0.3, 0.01 )

# init rendering engine

target_frequency = 60.0

mode = sfVideoMode(1000, 1000, 32)
window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL
sfWindow_setVerticalSyncEnabled(window,true)
sfWindow_setActive(window, sfTrue)

font = sfFont_createFromFile(joinpath("/julia/playgrounds/resources/NotoSans-Black.ttf"))
@assert font != C_NULL

clock = sfClock_create()

# express time in seconds.
function get_time(clock ::Ptr{Nothing})
    sfTime_asSeconds(sfClock_getElapsedTime(clock))
end
function restart(clock ::Ptr{Nothing})
     sfTime_asSeconds(sfClock_restart(clock))
end

function drawSquare(window, square, x, y, color::RGB)
        function f2i(f)
            round(f*255)
        end
        function to_sfColor(color::RGB)
            sfColor_fromRGB(f2i(color.r) ,f2i(color.g) , f2i(color.b) )
        end

        sfRectangleShape_setPosition(square, sfVector2f(x, y))
        sfRectangleShape_setFillColor(square, to_sfColor(color) )

        sfRenderWindow_drawRectangleShape(window, square, C_NULL)
end

# main loop

running = true
t0 = 0
event_ref = Ref{sfEvent}()
while (running)

    while Bool(sfRenderWindow_isOpen(window)) && running
        frame_timestamp = get_time(clock)
        frame_time = frame_timestamp - t0
        t0 = frame_timestamp
        actual_frequency = 1 / frame_time

        # process events
        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            if event_ref[].type == sfEvtClosed
                sfRenderWindow_close(window)
                running = false
                println("Render window closed.")
            elseif event_ref[].type == sfEvtResized
                glViewport(9, 9, event_ref[].size.width, event_ref[].size.height)
                println("Viewport updated.")
            elseif event_ref[].type == sfEvtKeyPressed
                if event_ref[].key.code == sfKeyQ
                    running = false
                    println("Quitting application.")
                end
            end
        end
        time_event_processing = get_time(clock) - t0
        #println("Event processing took $time_event_processing seconds.")

        # Render scene

        sfRenderWindow_clear(window, sfColor_fromRGBA(24,20,18,255))

        s = 1000/ (2*L)
        square = sfRectangleShape_create()
        sfRectangleShape_setSize(square, sfVector2f(s,s) )
        sfRectangleShape_setFillColor(square,sfColor_fromRGBA(255,155,105,255) )

        for a in agents
            x = (a.position.x + L )* 1000/(2L)
            y = (a.position.y + L )* 1000/(2L)
            drawSquare(window, square, x, y, getColor(a) )
        end

        sweep!(agents, L, pandemic)

        time_render = get_time(clock) - time_event_processing - t0

        sf_text = sfText_create()
        sfText_setPosition(sf_text, sfVector2f(4,1000-20))
        sfText_setString(sf_text, "frame render time : $time_render seconds.")
        sfText_setFont(sf_text, font)
        sfText_setCharacterSize(sf_text, 14)
        sfRenderWindow_drawText(window, sf_text, C_NULL)

        sf_text2 = sfText_create()
        sfText_setPosition(sf_text2, sfVector2f(4,1000-40))
        sfText_setString(sf_text2, "fps : $actual_frequency")
        sfText_setFont(sf_text2, font)
        sfText_setCharacterSize(sf_text2, 14)
        sfRenderWindow_drawText(window, sf_text2, C_NULL)

        sfRenderWindow_display(window)
        time_to_sleep = 1 / target_frequency - time_render - time_event_processing
        if time_to_sleep > 0.005
            sleep(time_to_sleep-0.005)
        end
    end

end

sfFont_destroy(font)
sfRenderWindow_destroy(window)

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
x = [1.0, -1.0, 1.0, 0.0, 0.0, 0.0]
function optimize!(mock)
    MOIU.mock_optimize!(mock, [zeros(MOI.get(mock, MOI.NumberOfVariables()) - 12); x; x],
        (MOI.VectorAffineFunction{Float64},MOI.Zeros) => [zeros(16)])
end
for mock in mocks(optimize!)
    Tests.rearrangement_test(mock, config)
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
gradient_descent(x::ANode, g::ANode, alpha::Real=1e-3) = axpy!(x, g, -alpha)

function momentum(x::ANode, g::ANode, α::Real=1e-3, μ::Real=0.9)
    v = ANode(zeros(eltype(arg(x)), size(x))) # this is really v/α
    v1 = μ*v + g
    v2 = store!(v, v1)
    axpy!(x, v2, -α)
end

function adadelta(w::ANode, g::ANode, ϵ::Real=1e-6, ρ::Real=0.95)
  Eg² = ANode("E[g²]", zeros(eltype(arg(w)), size(w)))
  EΔ² = ANode("E[Δ²]", zeros(eltype(arg(w)), size(w)))

  u = store!(Eg², ρ*Eg² + (1-ρ)*(g.^2))
  Δ = g.*√((EΔ² + ϵ)./(u + ϵ))
  order(store!(EΔ², ρ*EΔ² + (1-ρ)*(Δ.^2)), axpy!(w, Δ, -1))
end

function adadelta_fast(w::ANode, g::ANode, ϵ::Real=1e-6, ρ::Real=0.95)
  Eg² = ANode("E[g²]", zeros(eltype(arg(w)), size(w)))
  EΔ² = ANode("E[Δ²]", zeros(eltype(arg(w)), size(w)))
  ANode(:adadelta, (w, g, Eg², EΔ²), size(w), (ϵ, ρ))
end

mutates(::ANode{:adadelta}) = true
function do_forward!(n::ANode{:adadelta}, w, g, Eg², EΔ²)
  ϵ = arg(n)[1]
  ρ = arg(n)[2]
  adadelta_update(Csize_t(length(w)), w, g, convert(eltype(w), ϵ), convert(eltype(w), ρ), Eg², EΔ²)
end

@cuda_kernel adadelta_update(n::Csize_t, v::&T, g::&T, e::T, r::T, Eg2::&T, ED2::&T) [(T, Cfloat, Cdouble)] """
  int i = blockDim.x*blockIdx.x + threadIdx.x;
  if (i < n) {
    Eg2[i] = r*Eg2[i] + (1-r)*(g[i]*g[i]);
    T D = g[i] * sqrt((ED2[i] + e)/(Eg2[i] + e));
    v[i] -= D;
    g[i] = 0;
    ED2[i] = r*ED2[i] + (1-r)*(D*D);
  }
"""

function adadelta_update{T}(n::Csize_t, v::Array{T}, g::Array{T}, ϵ::T, ρ::T, Eg²::Array{T}, EΔ²::Array{T})
  @assert n == length(v) == length(g) == length(g) == length(Eg²) == length(EΔ²)
  for i in 1:n
    Eg²[i] = ρ*Eg²[i] + (1-ρ)g[i].^2
    Δ = g[i] * √((EΔ²[i] + ϵ)/(Eg²[i] + ϵ))
    v[i] -= Δ
    EΔ²[i] = ρ*EΔ²[i] + (1-ρ)Δ.^2
  end
end

function adam(w::ANode, g::ANode, α::Real=1e-3, β1::Real=0.9, β2::Real=0.999, ϵ::Real=1e-8)
  t = ANode("t", ones(eltype(arg(w)), 1))
  m = ANode("m", zeros(eltype(arg(w)), size(w)))
  v = ANode("v", zeros(eltype(arg(w)), size(w)))

  αt = α*sqrt(1-β2.^t)./(1-β1.^t)

  sm = store!(m, β1*m + (1-β1)*g)
  sv = store!(v, β2*v + (1-β2)*(g.^2))

  order(store!(t,t+1), store!(w, w - αt.*(sm./(sqrt(sv)+ϵ))))
end

function adam_fast(w::ANode, g::ANode, α::Real=1e-3, β1::Real=0.9, β2::Real=0.999, ϵ::Real=1e-8)
  t = ANode("t", ones(eltype(arg(w)), 1))
  m = ANode("m", zeros(eltype(arg(w)), size(w)))
  v = ANode("v", zeros(eltype(arg(w)), size(w)))
  ANode(:adam, (w, g, t, m, v), size(w), (α, β1, β2, ϵ))
end

mutates(::ANode{:adam}) = true
function do_forward!(n::ANode{:adam}, w, g, t, m, v)
  α,β1,β2,ϵ = arg(n)
  adam_update(Csize_t(length(w)), w, g, convert(eltype(w), α), convert(eltype(w), β1), convert(eltype(w), β2), convert(eltype(w), ϵ), t, m, v)
end

@cuda_kernel adam_update(n::Csize_t, val::&T, g::&T, a::T, b1::T, b2::T, e::T, t::&T, m::&T, v::&T) [(T, Cfloat, Cdouble)] """
  int i = blockDim.x*blockIdx.x + threadIdx.x;
  if (i < n) {
    T gi = g[i];
    m[i] = b1*m[i] + (1-b1)*gi;
    v[i] = b2*v[i] + (1-b2)*(gi*gi);
    T at = a*sqrt(1-pow(b2,t[0]))/(1-pow(b1,t[0]));
    val[i] -= at*m[i]/(sqrt(v[i])+e);
    g[i] = 0;
  }
  if (i == 0) t[0] += 1;
"""

function adam_update{T}(n::Csize_t, val::Array{T}, g::Array{T}, α::T, β1::T, β2::T, ϵ::T, t::Array{T}, m::Array{T}, v::Array{T})
  @assert n == length(val) == length(g) == length(m) == length(v)
  αt = α*sqrt(1-β2^t[1])/(1-β1^t[1])
  t[1] += 1
  @inbounds for i in 1:n
    m[i] = β1*m[i] + (1-β1)*g[i]
    v[i] = β2*v[i] + (1-β2)*(g[i].^2)
    val[i] -= αt*m[i]/(sqrt(v[i])+ϵ)
  end
end

function rmsprop(w::ANode, g::ANode, α::Real=1e-4, ρ::Real=0.95, μ::Real=0.9, ϵ::Real=1e-6)
  Eg² = ANode("E[g²]", zeros(eltype(arg(w)), size(w)))
  u = store!(Eg², ρ*Eg² + (1-ρ)*(g.^2))
  w1 = store!(w, μ*w - α*g./√(u + ϵ))
  axpy!(w, w1, 1)
end

function rmsprop_fast(w::ANode, g::ANode, α::Real=1e-4, μ::Real=0.9, ρ::Real=0.95, ϵ::Real=1e-6)
  Eg² = ANode("E[g²]", zeros(eltype(arg(w)), size(w)))
  ANode(:rmsprop, (w, g, Eg²), size(w), (α, ρ, μ, ϵ))
end

mutates(::ANode{:rmsprop}) = true
function do_forward!(n::ANode{:rmsprop}, w, g, Eg²)
  α = arg(n)[1]
  ρ = arg(n)[2]
  μ = arg(n)[3]
  ϵ = arg(n)[4]
  rmsprop_update(Csize_t(length(w)), w, g, convert(eltype(w), α), convert(eltype(w), ρ), convert(eltype(w), μ), convert(eltype(w), ϵ), Eg²)
end

@cuda_kernel rmsprop_update(n::Csize_t, v::&T, g::&T, a::T, r::T, m::T, e::T, Eg2::&T) [(T, Cfloat, Cdouble)] """
  int i = blockDim.x*blockIdx.x + threadIdx.x;
  if (i < n) {
    Eg2[i] = r*Eg2[i] + (1-r)*(g[i]*g[i]);
    T D = m*v[i] - a * g[i] /sqrt(Eg2[i] + e);
    v[i] += D;
    g[i] = 0;
  }
"""

function rmsprop_update{T}(n::Csize_t, w::Array{T}, g::Array{T}, α::T, ρ::T, μ::T, ϵ::T, Eg²::Array{T})
  @assert n == length(w) == length(g) == length(Eg²)
  @inbounds for i in 1:n
    Eg²[i] = ρ*Eg²[i] + (1-ρ)*g[i]^2
    Δ = μ*w[i] - α * g[i]/sqrt(Eg²[i] + ϵ)
    w[i] = Δ
  end
end


text:
using D4M

n = 2 .^ [6 6 7 8 9 10 11 12 13 14]

K = 8

MaxGB = 2
MaxGF = 4 * 1.6

m = K .* n
assoc_gbytes = zeros(1,length(n))
assoc_flops = zeros(1,length(n))
assoc_gflops = zeros(1,length(n))
assoc_time = zeros(1,length(n))

for i = 1:length(n)

    ii = round(Int, floor(rand(m[i]) .* n[i]) +1)
    jj = round(Int, floor(rand(m[i]) .* n[i]) +1)
    A = Assoc(ii,jj,1.0)

    ii = round(Int, floor(rand(m[i]) .* n[i]) +1)
    jj = round(Int, floor(rand(m[i]) .* n[i]) +1)
    B = Assoc(ii,jj,1.0)

    tic()
        C = CatValMul(A,B)
    assoc_time[i] = toq()
    assoc_flops[i] = 2*sum(C)
    ii, jj, vv = find(C)
    assoc_gbytes[i] = assoc_gbytes[i] + (length(ii) + length(jj)) + 8 .* m[i] ./ 1e9
    assoc_gflops[i] = assoc_flops[i] ./ assoc_time[i] ./ 1e9
    print("Time: ", assoc_time[i])
    print(", GFlops: ", assoc_gflops[i])
    println(", GBytes: ", assoc_gbytes[i])
end


text:
using BayesianDataFusion
using Base.Test

using DataFrames

A = randn(30,2);
B = randn(40,2);

X = Float64[ sum(A[i,:].*B[j,:]) for i in 1:size(A,1), j in 1:size(B,1) ]

df = DataFrame(A=Int64[], B=Int64[], v=Float64[])
for i=1:size(A,1), j=1:size(B,1)
  push!(df, Any[i, j, X[i,j]])
end

## adding features
feat = randn(size(A,1), 2)

rd = RelationData(df)
rd.entities[1].F = feat
rd.entities[1].lambda_beta_sample = true

assignToTest!(rd.relations[1], 10)

result = macau(rd, burnin=50, psamples=10, num_latent=2, verbose=false)


text:
#
# WIP: work in progress on deepvalidate. This code is not expected to work at this time
#

#
# Julia implementation of LSTM/RNN classifier for classifying program
# execution traces.
#

#
# Loop to check if packages are installed and install them if not, then
# import.
#

# packages = ['Flux','DelimitedFiles','MLDataPattern']
# for package in packages
#     haskey(Pkg.installed(),package) || Pkg.add(package)

using DelimitedFiles
using Flux
using Flux: onehot, throttle, crossentropy, onehotbatch, params, shuffle
using MLDataPattern: stratifiedobs
using Base.Iterators: partition

include("../utils.jl")

#
# Set up inputs for model
#

# Read lines from traces.dat text in to arrays of characters
# Convert to onehot matrices

cd(@__DIR__)

text, alphabet, N = get_data("traces.csv")
stop = onehot('\n', alphabet);

# Partition into subsequences to input to our model

seq_len = 50

Xs = [collect(partition(t,seq_len)) for t in text];
Ys = readdlm("y_results.csv");

dataset = [(onehotbatch(x, alphabet, '\n'), onehot(Ys[i], unique(Ys)))
           for i in 1:length(Ys) for x in Xs[i]] |> shuffle

# Pad sequences to equal lengths

Xs = [hcat(x,repeat(stop,1,seq_len-size(x)[2])) for x in first.(dataset)]
Ys = last.(dataset)

# There are 972,290 items in our data. We use a train:test split of 90:10

(Xtrain, Ytrain), (Xtest, Ytest) = stratifiedobs((Xs, Ys), p=0.9)

train = [(Xtrain[i], Ytrain[i]) for i in 1:length(Ytrain)];
test = [(Xtest[i], Ytest[i]) for i in 1:length(Ytest)];

scanner = Chain(Dense(length(alphabet), seq_len, σ), LSTM(seq_len, seq_len))
encoder = Dense(seq_len, 2)

function model(x)
  state = scanner.([x])[end]
  Flux.reset!(scanner)
  softmax(encoder(state))
end

loss(tup) = crossentropy(mod(tup[1]), tup[2])
accuracy(tup) = mean(argmax(m(tup[1])) .== argmax(tup[2]))

testacc() = mean(accuracy(t) for t in test)
testloss() = mean(loss(t) for t in test)

opt = ADAM(0.01)
ps = params(mod)
evalcb = () -> @show testloss(), testacc()

Flux.train!(loss, ps, train, opt, cb = throttle(evalcb, 10))


text:
function hmm_forward_alg(prior::Vector{Float64},
        emission_dists::AbstractArray{Float64,2}, transition_dists::AbstractArray{Float64,2},
        emissions::Vector{Int})
    marg_lik = 1.
    alpha = prior # p(z_1)
    for i=2:length(emissions)

        # p(z_{i-1} , y_{i-1} | y_{1:i-2}) for each z_{i-1}
        prev_posterior = alpha .* emission_dists[emissions[i-1], :]

        # p(y_{i-1} | y_{1:i-2})
        denom = sum(prev_posterior)

        # p(z_{i-1} | y_{1:i-1})
        prev_posterior = prev_posterior / denom

        # p(z_i | y_{1:i-1})
        alpha = transition_dists * prev_posterior

        # p(y_{1:i-1})
        marg_lik *= denom
    end
    prev_posterior = alpha .* emission_dists[emissions[end], :]
    denom = sum(prev_posterior)
    marg_lik *= denom
    marg_lik
end

@testset "hmm forward alg" begin

    # test the hmm_forward_alg on a hand-calculated example
    prior = [0.4, 0.6]
    emission_dists = [0.1 0.9; 0.7 0.3]'
    transition_dists = [0.5 0.5; 0.2 0.8]'
    obs = [2, 1]
    expected_marg_lik = 0.
    # z = [1, 1]
    expected_marg_lik += prior[1] * transition_dists[1, 1] * emission_dists[obs[1], 1] * emission_dists[obs[2], 1]
    # z = [1, 2]
    expected_marg_lik += prior[1] * transition_dists[2, 1] * emission_dists[obs[1], 1] * emission_dists[obs[2], 2]
    # z = [2, 1]
    expected_marg_lik += prior[2] * transition_dists[1, 2] * emission_dists[obs[1], 2] * emission_dists[obs[2], 1]
    # z = [2, 2]
    expected_marg_lik += prior[2] * transition_dists[2, 2] * emission_dists[obs[1], 2] * emission_dists[obs[2], 2]
    actual_marg_lik = hmm_forward_alg(prior, emission_dists, transition_dists, obs)
    @test isapprox(actual_marg_lik, expected_marg_lik)

end

@testset "particle filtering" begin

    prior = [0.2, 0.3, 0.5]

    emission_dists = [
        0.1 0.2 0.7;
        0.2 0.7 0.1;
        0.7 0.2 0.1
    ]'

    transition_dists = [
        0.4 0.4 0.2;
        0.2 0.3 0.5;
        0.9 0.05 0.05
    ]'

    @gen function kernel(t::Int, prev_z::Int, params::Nothing)
        z = @trace(categorical(transition_dists[:,prev_z]), :z)
        @trace(categorical(emission_dists[:,z]), :x)
        return z
    end

    chain = Unfold(kernel)

    @gen function model(num_steps::Int)
        z_init = @trace(categorical(prior), :z_init)
        @trace(categorical(emission_dists[:,z_init]), :x_init)
        @trace(chain(num_steps-1, z_init, nothing), :chain)
    end

    num_steps = 4
    obs_x = [1, 1, 2, 3]

    # latents:
    # :z_init
    # :chain => 1 => :z
    # :chain => 2 => :z
    # :chain => 3 => :z

    # observations :
    # :x_init
    # :chain => 1 => :x
    # :chain => 2 => :x
    # :chain => 3 => :x

@testset "custom proposal" begin

    Random.seed!(0)
    num_particles = 10000
    ess_threshold = 10000 # make sure we exercise resampling

    # initialize particle filter

    @gen function init_proposal(x::Int)
        dist = prior .* emission_dists[x,:]
        @trace(categorical(dist ./ sum(dist)), :z_init)
    end

    init_proposal_args = (obs_x[1],)
    init_observations = choicemap((:x_init, obs_x[1]))

    state = initialize_particle_filter(model, (1,), init_observations,
        init_proposal, init_proposal_args, num_particles)

    # do particle filter steps

    @gen function step_proposal(prev_trace, T::Int, x::Int)
        @assert T > 1
        choices = get_choices(prev_trace)
        if T > 2
            prev_z = choices[:chain => (T-2) => :z]
        else
            prev_z = choices[:z_init]
        end
        dist = transition_dists[:,prev_z] .* emission_dists[x,:]
        @trace(categorical(dist ./ sum(dist)), :chain => T-1 => :z)
    end

    argdiffs = (UnknownChange(),) # the length may change
    for T=2:length(obs_x)
        maybe_resample!(state, ess_threshold=ess_threshold)
        new_args = (T,)
        observations = choicemap((:chain => (T-1) => :x, obs_x[T]))
        proposal_args = (T, obs_x[T])
        particle_filter_step!(state, new_args, argdiffs, observations,
            step_proposal, proposal_args)
    end

    # check log marginal likelihood estimate
    expected_log_ml = log(hmm_forward_alg(prior, emission_dists, transition_dists, obs_x))
    actual_log_ml_est = log_ml_estimate(state)
    @test isapprox(expected_log_ml, actual_log_ml_est, atol=0.01)
end

@testset "default proposal" begin

    Random.seed!(0)
    num_particles = 10000
    ess_threshold = 10000 # make sure we exercise resampling

    # initialize the particle filter
    init_observations = choicemap((:x_init, obs_x[1]))
    state = initialize_particle_filter(model, (1,), init_observations, num_particles)

    # do steps
    argdiffs = (UnknownChange(),) # the length may change
    for T=2:length(obs_x)
        maybe_resample!(state, ess_threshold=ess_threshold)
        new_args = (T,)
        observations = choicemap((:chain => (T-1) => :x, obs_x[T]))
        particle_filter_step!(state, new_args, argdiffs, observations)
    end

    # check log marginal likelihood estimate
    expected_log_ml = log(hmm_forward_alg(prior, emission_dists, transition_dists, obs_x))
    actual_log_ml_est = log_ml_estimate(state)
    @test isapprox(expected_log_ml, actual_log_ml_est, atol=0.01)
end

end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen