text:
using Simplices
using Simplices: SimplexSplitting

"""
    intersectiontest(dim::Int, n_reps::Int, intersection_type::String, tolerance::Float64)
Generate a `n_reps` pairs of `dim`-dimensional simplices that intersect in the way specified by `intersecton_type`. Valid intersection types are `"nontrivial"` and `"sharingvertices"`. Checks if the total volume spanned by each pair of simplices matches that obtained when accounting for the shared volume.
"""
function intersectiontest(dim::Int, n_reps::Int, intersection_type::String; tolerance = 1/10^12)

    discrepancies = zeros(Float64, n_reps)
    elapsed = zeros(Float64, n_reps)
    for i = 1:n_reps
        s₁, s₂ = intersecting_simplices(dim = dim, intersection_type = intersection_type)
        before = time_ns()/10^9
        intvol = simplexintersection(copy(transpose(s₁)), copy(transpose(s₂)), tolerance = tolerance, what = "volume")
        after =  time_ns()/10^9
        elapsed[i] = after - before
        vol₁ = volume(s₁)
        vol₂ = volume(s₂)
        total_vol = vol₁ + vol₂
        discrepancies[i] = total_vol - ((vol₁ - intvol) + (vol₂ - intvol) + intvol)
        if discrepancies[i] > 1/10^5
            @show vol₁, vol₂, discrepancies[i]
        end
    end
    return discrepancies, elapsed
end

"""
    intersectiontest(dim::Int, n_reps::Int, intersection_type::String, tolerance::Float64)
Generate a `n_reps` pairs of `dim`-dimensional simplices that intersect in the way specified by `intersecton_type`. Valid intersection types are `"nontrivial"` and `"sharingvertices"`. Checks if the total volume spanned by each pair of simplices matches that obtained when accounting for the shared volume.
"""
function intersectiontest(dim::Int, n_reps::Int, intersection_type::String)

    discrepancies = zeros(Float64, n_reps)
    for i = 1:n_reps
        s₁, s₂ = intersecting_simplices(dim = dim, intersection_type = intersection_type)
        intvol = simplexintersection(copy(transpose(s₁)), copy(transpose(s₂)))
        vol₁ = volume(s₁)
        vol₂ = volume(s₂)
        total_vol = vol₁ + vol₂
        discrepancies[i] = total_vol - (intvol + (vol₁ - intvol) + (vol₂ - intvol) )
        if discrepancies[i] > 1/10^5
            @show s₁, s₂
        end
    end
    return discrepancies
end

# Splits a canonical simplex of dimension E with size division factor k. Inside the simplex,
# it generates a random subsimplex. Compute its volume by computing the intersection of
# that simplex with the simplices forming the splitting of the original simplex. Compare
# the volume obtained by the simplex intersection routine to the analytical volume.
# N is the number of times to perform the test, and tolerance is the usual tolerance in the
# simplex intersection routine.
function nd_Test(k, E, N; tolerance = 1/10^12, plot = false)

    # Define vertices of canonical simplex
    canonical_simplex_vertices = zeros(E + 1, E)
    canonical_simplex_vertices[2:(E+1), :] = Matrix(1.0I, E, E)
    simplex_indices = zeros(Int, 1, E + 1)
    simplex_indices[1, :] = round.(Int, collect(1:E+1))

    refined = refine_triangulation(canonical_simplex_vertices, simplex_indices, [1], k)
    triang_vertices, triang_simplex_indices = refined[1], refined[2]

    differences = zeros(Float64, N)

    # Repeat the test N times
    for i = 1:N
        # Build convex expansion coefficients for the random simplex. We need these in a
        # form that guarantees that the random simplex lies within the canonical simplex.
        beta = rand(E + 1, E + 1)
        beta = beta .* repeat(1 ./ sum(beta, dims = 2), 1, E + 1)
        # Ensure that we have convex expansions

        # Linear combination of the original vertices and the convex expansion coefficients.
        # Creates a matrix containing the vertices of the random simplex (which now is
        # guaranteed to lie _within_ the original simplex.
        random_simplex = beta * canonical_simplex_vertices
        random_simplex_orientation = det(hcat(ones(E + 1), random_simplex))
        random_simplex_volume = abs(random_simplex_orientation)
        random_simplex = copy(transpose(random_simplex))

        # Intersection between each of the subsimplices with the random simplex
        intersecting_volumes = zeros(Float64, size(triang_simplex_indices, 1))

        for j = 1:size(triang_simplex_indices, 1)
            # Get the subsimplex vertices
            subsimplex = copy(transpose(triang_vertices[triang_simplex_indices[j, :], :]))
            intvol = simplexintersection(random_simplex, subsimplex, tol = tolerance)
            intersecting_volumes[j] = intvol
        end

        # Compute the discrepancies
        numeric_volume = sum(intersecting_volumes)
        analytic_volume = random_simplex_volume
        differences[i] = abs((numeric_volume - analytic_volume)/analytic_volume)

    end
    return differences
end

function test_nontrivial(dim,N)
   for i = 1:N
        S1,S2 = nontrivially_intersecting_simplices(dim)
        simplexintersection(copy(transpose(S1)), copy(transpose(S2)))
    end

    return true
end

function test_sharing(dim,N)
    for i = 1:N
        S1,S2 = simplices_sharing_vertices(dim)
        simplexintersection(copy(transpose(S1)), copy(transpose(S2)))
    end
    return true
end

@testset "Nontrivial intersection" begin
    n = 10
    @testset "E = $E" for E in 3:5
        @test test_nontrivial(E, n) == true
    end
end

@testset "Sharing a vertex" begin
    n = 10
    @testset "E = $E" for E in 3:5
        @test test_sharing(E, n) == true
    end
end

@test nd_Test(2, 3, 1)[1] < 1e-8
@test nd_Test(2, 4, 1)[1] < 1e-8
@test nd_Test(2, 5, 1)[1] < 1e-8

@testset "nD discrepancy: strictly contained" begin
    # Trigger precompilation on simplest example
    k = 2; reps = 1
    function init()
        for E = 3
            discrepancies = nd_Test(k, E, reps)
            @test maximum(discrepancies) < 1e-8
        end
    end

    init()

    reps = 10
    # Time more involved examples
    println("nD discrepancy test with ", reps, " reps:")
    @testset "E = $E" for E in 3:5
        # Trigger once for precompilation
        @testset "k = $k" for k in 2:3
            t1 = time_ns()/10^9
            discrepancies = nd_Test(k, E, reps)
            t2 = time_ns()/10^9
            @test maximum(discrepancies) < 1e-8
        end
    end
end


text:
import ElasticFusion
import Surfels
import Pangolin
import SLAMData
import CameraModel

function fuse_RGBD_to_3D_model(
    rgbd_reader::T,
    id2Twc, # traj
    cam::CameraModel.RgbdCamParams,
    output_ply_fpath::String
) where T <: SLAMData.RgbdReader

    # hidden_window = Surfels.create_hidden_gl_window(cam.width, cam.height)
    window_name = "EFusion Window"
    Pangolin.create_window(window_name, cam.width, cam.height)

    ef = ElasticFusion.Fusion(
        w=cam.width,
        h=cam.height,
        fx=cam.fx,
        fy=cam.fy,
        cx=cam.cx,
        cy=cam.cy,
        timeDelta = 200,
        countThresh = 35000,
        errThresh = 5e-05,
        covThresh = 1e-05,
        closeLoops = false,
        iclnuim = false,
        reloc = false,
        photoThresh = 115,
        confidence = 10,
        depthCut = 5,
        icpThresh = 10,
        fastOdom = false,
        fernThresh = 0.3095,
        so3 = true,
        frameToFrameRGB = false,
        fileName = "fsd_2021sep"
    )

    frame_ids = sort!(collect(keys(id2Twc)))
    @assert frame_ids[1] == 0

    for f_id = frame_ids

        # TODO: how to get timestamp from id
        time_stamp = Int64(f_id)

        # pose = convert(Matrix{Float32}, inv(all_Twc0[f_id]))
        pose = convert(Matrix{Float32}, id2Twc[f_id])

        img, depth16 = SLAMData.readRGBD(rgbd_reader, id=time_stamp)

        # note depth_factor only support 1000.0 for now
        @assert cam.depth_factor == 1000.0

        ElasticFusion.process_frame!(
            ef, img, depth16, time_stamp, pose=pose, depth_factor=1000.0
        )

        # function process_frame!(
        #     ef::ElasticFusion,
        #     image::Matrix{RGB{Normed{UInt8,8}}},
        #     depth::Matrix{UInt16},
        #     timestamp::Int64,
        #     pose::Matrix{Float32};
        #     weightMultiplier::Union{Float32, Float64}
        #     bootstrap::Bool
        #     )

        # Pangolin.activate(surfel_viewer.d_cam, surfel_viewer.s_cam)
        # glClearColor(0.05 , 0.05 , 0.3 , 0.0);
        # glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        # # end
        # mvp = Pangolin.get_projection_view_mat(surfel_viewer.s_cam)
        # # mvp32 = convert(Matrix{Float32}, mvp)

        # ElasticFusion.render(
        #     ef,
        #     mvp,
        #     threshold=10,#::Union{Float64, Float32},10
        #     drawUnstable=true,#//bool
        #     drawNormals=false,#, //bool
        #     drawColors=true,#,  //bool
        #     drawPoints=true,#,  //bool
        #     drawWindow=false,#,  //bool
        #     drawTimes=false,#,   //bool
        #     time=time_stamp,
        #     timeDelta=200
        # )
        # Pangolin.finish_frame()
    end

    ElasticFusion.save_ply(ef, output_ply_fpath)
    finalize(ef)
    Pangolin.destroy_window(window_name)
end

ext: jl


text:
@testset "Datastructures" begin
    # Test the types and the sizes of the moments
    # For 64 bit floating numbers, one dimensional
    @testset "One dimensional moments" begin
        n = 10
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            moments = JuSwalbe.Macroquant{Vector{type},Vector{type}}(zeros(type,n),zeros(type,n),zeros(type,n),zeros(type,n))
            @test isa(moments, JuSwalbe.Macroquant)
            test_dict = struct2dict(moments)
            for key in keys(test_dict)
                @test isa(test_dict[key], Array{type,1})
                @test size(test_dict[key], 1) == n
            end
        end
    end
    @testset "One dimensional forces" begin
        n = 10
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            forces = JuSwalbe.Forces{Vector{type}}(zeros(type,n),zeros(type,n),zeros(type,n),zeros(type,n))
            @test isa(forces, JuSwalbe.Forces)
            test_dict = struct2dict(forces)
            for key in keys(test_dict)
                @test isa(test_dict[key], Array{type,1})
                @test size(test_dict[key], 1) == n
            end
        end
    end
    # For 32 bit floating numbers, one dimensional
    @testset "Two dimensional moments" begin
        m,n = 10,8
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            twovec = JuSwalbe.Twovector{Matrix{type}}(zeros(type, (m,n)), zeros(type, (m,n)))
            moments = JuSwalbe.Macroquant{Matrix{type},JuSwalbe.Twovector{Matrix{type}}}(zeros(type, (m,n)),twovec ,zeros(type, (m,n)),zeros(type, (m,n)))
            @test isa(moments, JuSwalbe.Macroquant)
            test_dict = struct2dict(moments)
            for key in keys(test_dict)
                if key == :velocity
                    @test isa(test_dict[key], JuSwalbe.Twovector)
                    @test isa(test_dict[key].x, Matrix{type})
                    @test isa(test_dict[key].y, Matrix{type})
                    @test size(test_dict[key].x) == (m,n)
                    @test size(test_dict[key].y) == (m,n)
                else
                    @test isa(test_dict[key], Array{type,2})
                    @test size(test_dict[key], 1) == m
                    @test size(test_dict[key], 2) == n
                end
            end
        end
    end
    @testset "Two dimensional forces" begin
        n, m = 10, 6
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            twovec = JuSwalbe.Twovector{Matrix{type}}(zeros(type, (n,m)),zeros(type, (n,m)))
            forces = JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{type}}}(twovec, twovec, twovec, twovec)
            @test isa(forces, JuSwalbe.Forces)
            test_dict = struct2dict(forces)
            for key in keys(test_dict)
                @test isa(test_dict[key], JuSwalbe.Twovector)
                @test isa(test_dict[key].x, Matrix{type})
                @test isa(test_dict[key].y, Matrix{type})
                @test size(test_dict[key].x) == (n,m)
                @test size(test_dict[key].y) == (n,m)
            end
        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
function waveAnim(path::String, path1::String; itype="wfd", cpt="Vz",cmap="PuOr",
                  pclip=0.98, vlim=1.0, aspect="auto", interpolation="hanning", ox=0.0,
                  oy=0.0, dx=10.0, dy=10.0, dt1=0.002, wbox=8, hbox=8, spd=40)
    fid = open(path1, "r")
    nz = read(fid, Int32); nx = read(fid, Int32); nz = Int64(nz); nx = Int64(nx);
    ext= read(fid, Int32); iflag = read(fid, Int32)
    dt = read(fid, Float32);
    inc = round(Int64, dt1/dt)
    fig = plt.figure(figsize=(wbox, hbox))
    ims = PyCall.PyObject[]
    if itype == "wfd"
       sptSize = sizeof(Float32)*nz*nx*5
       nt = floor(Int64, (filesize(fid)-sizeof(Float32)*5)/sptSize)
       for it = 1 : inc : nt
           wfd = ReadWfd(path1, it)
           if cpt == "Vx"
              d = reshape(wfd.Vx, nz, nx)
           elseif cpt == "Vz"
              d = reshape(wfd.Vz, nz, nx)
           elseif cpt == "Txx"
              d = reshape(wfd.Txx, nz, nx)
           elseif cpt == "Tzz"
              d = reshape(wfd.Tzz, nz, nx)
           elseif cpt == "Txz"
              d = reshape(wfd.Txz, nz, nx)
           end
           if pclip < 1.0
              vlim = quantile(vec(abs(d)), pclip)
           end
           im = plt.imshow(d, cmap=cmap, vmin=-vlim, vmax=vlim, extent=[ox, ox+(size(d,2)-1)*dx, oy+(size(d,1)-1)*dy,oy], aspect=aspect, interpolation=interpolation)
           push!(ims, PyCall.PyObject[im])
       end
    elseif itype == "spt"
       if iflag == 1
          Nz = nz + 2*ext
          upper = ext
       elseif iflag == 2
          Nz = nz +   ext
          upper = Int32(0)
       end
       Nx  = nx + 2*ext
       sptSize = sizeof(Float32)*Nz*Nx*10
       nt = floor(Int64, (filesize(fid)-sizeof(Float32)*5)/sptSize)
       for it = 1 : inc : nt
           spt = ReadSnapShot(path1, it)
           if cpt == "Vx"
              tmp = reshape(spt.Vxx+spt.Vxz, Nz, Nx)
              d = tmp[upper+1:upper+nz, ext+1:nx+ext]
           elseif cpt == "Vz"
              tmp = reshape(spt.Vzx+spt.Vzz, Nz, Nx)
              d = tmp[upper+1:upper+nz, ext+1:nx+ext]
           elseif cpt == "Txx"
              tmp = reshape(spt.Txxx+spt.Txxz, Nz, Nx)
              d = tmp[upper+1:upper+nz, ext+1:nx+ext]
           elseif cpt == "Tzz"
              tmp = reshape(spt.Tzzx+spt.Tzzz, Nz, Nx)
              d = tmp[upper+1:upper+nz, ext+1:nx+ext]
           elseif cpt == "Txz"
              tmp = reshape(spt.Txzx+spt.Txzz, Nz, Nx)
              d = tmp[upper+1:upper+nz, ext+1:nx+ext]
           end
           if pclip < 1.0
              vlim = quantile(vec(abs(d)), pclip)
           end
           im = plt.imshow(d, cmap=cmap, vmin=-vlim, vmax=vlim, extent=[ox, ox+(size(d,2)-1)*dx, oy+(size(d,1)-1)*dy,oy], aspect=aspect, interpolation=interpolation)
           push!(ims, PyCall.PyObject[im])
       end
    end
    ani = anim.ArtistAnimation(fig, ims, interval=spd, blit=true)
    path = join([path "_" cpt ".mp4"])
    ani[:save](path)
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
@doc raw"""
    PolynomialKernel([a=1 [,c=1 [,d=3]]])

The polynomial kernel is a Mercer kernel given by:

```
    κ(x,y) = (a⋅xᵀy + c)ᵈ   α > 0, c ≧ 0, d ∈ ℤ⁺
```

# Examples

```jldoctest; setup = :(using MLKernels)
julia> PolynomialKernel(2.0f0)
PolynomialKernel{Float32}(2.0,1.0,3)

julia> PolynomialKernel(2.0f0, 2.0)
PolynomialKernel{Float64}(2.0,2.0,3)

julia> PolynomialKernel(2.0f0, 2.0, 2)
PolynomialKernel{Float64}(2.0,2.0,2)
```
"""
struct PolynomialKernel{T<:AbstractFloat} <: MercerKernel{T}
    a::T
    c::T
    d::T
    function PolynomialKernel{T}(
            a::Real=T(1),
            c::Real=T(1),
            d::Real=T(3)
        ) where {T<:AbstractFloat}
        @check_args(PolynomialKernel, a, a >  zero(a), "a > 0")
        @check_args(PolynomialKernel, c, c >= zero(c), "c ≧ 0")
        @check_args(PolynomialKernel, d, d >= one(d) && d == trunc(d), "d ∈ ℤ₊")
        return new{T}(a, c, d)
    end
end

function PolynomialKernel(
        a::T₁=1.0,
        c::T₂=T₁(1),
        d::T₃=convert(promote_float(T₁,T₂), 3)
    ) where {T₁<:Real,T₂<:Real,T₃<:Real}
    T = promote_float(T₁,T₂,T₃)
    return PolynomialKernel{T}(a, c, d)
end

@inline basefunction(::PolynomialKernel) = ScalarProduct()

@inline function kappa(κ::PolynomialKernel{T}, xᵀy::T) where {T}
    return (κ.a*xᵀy + κ.c)^(κ.d)
end

function convert(::Type{K}, κ::PolynomialKernel) where {K>:PolynomialKernel{T}} where T
    return PolynomialKernel{T}(κ.a, κ.c, κ.d)
end


text:
module ASDF

using PyCall
using URIParser

const asdf = PyNULL()
function __init__()
    # copy!(asdf, pyimport_conda("asdf", "asdf", "conda-forge"))
    # copy!(asdf, pyimport_conda("asdf", "asdf", "astropy"))
    copy!(asdf, pyimport("asdf"))
end

################################################################################

struct Empty end
const empty = Empty()

const Maybe{T} = Union{Nothing, T}
nothing2tuple(x) = x
nothing2tuple(::Nothing) = ()

################################################################################

const tag2asdftype = Dict{String, Type}()

function makeASDFType(pyobj::PyObject)
    tag = try
        pyobj.yaml_tag
    catch
        # Convert to a nice Julia type if possible
	converted = convert(PyAny, pyobj)
	return (typeof(converted) == Dict{Any,Any} ? Tree(converted) : converted)
    end
    type_ = get(tag2asdftype, tag, empty)
    if type_ === empty
        @show tag
        @assert false           # for debugging
        return pyobj
    end
    type_(pyobj)::ASDFType
end

################################################################################

struct File
    pyobj::PyObject
end

function File(dict::Dict)       # for convenience
    File(asdf.AsdfFile(dict))
end

function open(filename::AbstractString)::File
    File(asdf.open(filename))
end

function close(file::File)::Nothing
    file.pyobj.close()
end

function write_to(file::File, filename::AbstractString)::Nothing
    file.pyobj.write_to(filename)
end

function file_format_version(file::File)::VersionNumber
    VersionNumber(file.pyobj.file_format_version)
end

function version(file::File)::VersionNumber
    VersionNumber(file.pyobj.version)
end

function Base.VersionNumber(obj::PyObject)
    VersionNumber(
        obj.major, obj.minor, obj.patch, obj.prerelease, obj.build)
end
function Base.VersionNumber(obj::Dict)
    VersionNumber(
        obj["major"], obj["minor"], obj["patch"],
        nothing2tuple(obj["prerelease"]), nothing2tuple(obj["build"]))
end

function comments(file::File)::Vector{String}
    file.pyobj.comments
end

function tree(file::File)::Tree
    makeASDFType(file.pyobj."tree"::PyObject)::Tree
end

################################################################################

# id: "http://stsci.edu/schemas/asdf/core/asdf-1.1.0"
# tag: "tag:stsci.edu:asdf/core/asdf-1.1.0"
# title: Top-level schema for every ASDF file.
# additionalProperties: true

"""Top-level schema for every ASDF file"""
struct Tree
    pyobj::PyObject
end
tag2asdftype["tag:stsci.edu:asdf/core/asdf-1.0.0"] = Tree
tag2asdftype["tag:stsci.edu:asdf/core/asdf-1.1.0"] = Tree
additionalProperties(::Tree) = ()

# This constructor must come after the type Tree has been defined
# TODO: Reverse order of types in this file
function File(tree::Tree)
    File(asdf.AsdfFile(tree.pyobj))
end

function asdf_library(tree::Tree)::Maybe{Software}
    software = get(tree.pyobj, PyObject, "asdf_library", empty)
    if software === empty
        return nothing
    end
    software::PyObject
    makeASDFType(software)::Software
end

# History is apparently not supported by the Python ASDF library

################################################################################

# id: "http://stsci.edu/schemas/asdf/core/software-1.0.0"
# tag: "tag:stsci.edu:schemas/asdf/core/software-1.0.0"
# title: Describes a software package.
# required: [name, version]
# additionalProperties: true

"""Describes a software package"""
struct Software
    pyobj::PyObject
end
tag2asdftype["tag:stsci.edu:asdf/core/software-1.0.0"] = Software
additionalProperties(::Software) = ()

function name(software::Software)
    software.pyobj.get("name")::String
end

function author(software::Software)
    software.pyobj.get("author")::Maybe{String}
end

function homepage(software::Software)::Maybe{URI}
    homepage = software.pyobj.get("homepage", empty)
    if homepage === empty
        return nothing
    end
    URI(homepage)::URI
end

function version(software::Software)::Union{VersionNumber, String}
    version = software.pyobj.get("version")::String
    try
        VersionNumber(version)
    catch
        version
    end
end

################################################################################

# id: "http://stsci.edu/schemas/asdf/core/ndarray-1.0.0"
# tag: "tag:stsci.edu:asdf/core/ndarray-1.0.0"
# title: An *n*-dimensional array.

abstract type Datatype end

@enum ScalarType begin
    int8
    uint8
    int16
    uint16
    int32
    uint32
    int64
    uint64
    float32
    float64
    complex64
    complex128
    bool8
    ascii
    ucs4
end
const string2scalartype = Dict{String, ScalarType}(
    "int8"       => int8,
    "uint8"      => uint8,
    "int16"      => int16,
    "uint16"     => uint16,
    "int32"      => int32,
    "uint32"     => uint32,
    "int64"      => int64,
    "uint64"     => uint64,
    "float32"    => float32,
    "float64"    => float64,
    "complex64"  => complex64,
    "complex128" => complex128,
    "bool8"      => bool8,
    "ascii"      => ascii,
    "ucs4"       => ucs4)
const scalartype2type = Dict{ScalarType, Type}(
    int8       => Int8,
    uint8      => UInt8,
    int16      => Int16,
    uint16     => UInt16,
    int32      => Int32,
    uint32     => UInt32,
    int64      => Int64,
    uint64     => UInt64,
    float32    => Float32,
    float64    => Float64,
    complex64  => ComplexF32,
    complex128 => ComplexF64,
    bool8      => Bool,
    ascii      => String,
    ucs4       => String)

struct ScalarDatatype <: Datatype
    type_::ScalarType
    length::Int                 # only for ascii and ucs4; else -1
    function ScalarDatatype(type_::ScalarType)
        @assert type_ in [
            int8, uint8, int16, uint16, int32, uint32, int64, uint64,
            float32, float64, complex64, complex128,
            bool8]
        new(type_, -1)
    end
    function ScalarDatatype(type_::ScalarType, length::Int)
        @assert type_ in [ascii, ucs4]
        @assert length >= 0
        new(type_, length)
    end
end

function ScalarDatatype(type_::String, length...)
    ScalarDatatype(string2scalartype[type_], length...)
end

function julia_type(scalardatatype::ScalarDatatype)
    scalartype2type[scalardatatype.type_]
end

@enum Byteorder big little

struct Field
    name::Maybe{String}
    datatype::Datatype
    byteorder::Maybe{Byteorder}
    shape::Maybe{Vector{Int}}
end

struct DatatypeList <: Datatype
    types::Vector{Field}
end

function Datatype(dtype::PyObject)
    if dtype.names !== nothing
        # fields = []
        # for name in dtype.names:
        #     field = dtype.fields[name][0]
        #     d = {}
        #     d['name'] = name
        #     field_dtype, byteorder = numpy_dtype_to_asdf_datatype(field)
        #     d['datatype'] = field_dtype
        #     if include_byteorder:
        #         d['byteorder'] = byteorder
        #     if field.shape:
        #         d['shape'] = list(field.shape)
        #     fields.append(d)
        # return fields, numpy_byteorder_to_asdf_byteorder(dtype.byteorder)
        @assert false

    elseif dtype.subdtype !== nothing
        # return numpy_dtype_to_asdf_datatype(dtype.subdtype[0])
        @assert false

    elseif dtype.name in keys(string2scalartype)
        return ScalarDatatype(dtype.name)

    elseif dtype.name == "bool"
        return ScalarDatatype(bool8)

    elseif startswith(dtype.name, "string") ||
            startswith(dtype.name, "bytes")
        return ScalarDatatype(ascii, Int(dtype.itemsize))

    elseif startswith(dtype.name, "unicode") ||
            startswith(dtype.name, "str")
        return ScalarDatatype(ucs4, Int(dtype.itemsize) ÷ 4)

    end
    @assert false
end

const PyArrayTypes = Union{
    Bool,
    Int8, Int16, Int32, Int64, Int128,
    UInt8, UInt16, UInt32, UInt64, UInt128,
    Float16, Float32, Float64,
    ComplexF16, ComplexF32, ComplexF64}

"""An *n*-dimensional array"""
struct NDArray{T, D, Repr <: Union{PyArray, PyObject}} <: DenseArray{T, D}
    # An NDArray either uses a PyArray (if possible) or a PyObject
    # as internal representation
    repr::Repr
    function NDArray{T, D, PyArray{T, D}}(pyobj::PyObject) where {T, D}
        new{T, D, PyArray{T, D}}(PyArray(pyobj))
    end
    function NDArray{T, D, PyObject}(pyobj::PyObject) where {T, D}
        new{T, D, PyObject}(pyobj)
    end
    function NDArray{T, D}(pyobj::PyObject) where {T, D}
        isefficient = T <: PyArrayTypes
        if isefficient
            # Python 3 requires array-like Python object to support a
            # "buffer protocol" as defined in PEP 3118. Apparently,
            # asdf.py's "NDArrayType" does not, and the call to
            # "PyArray_Info" fails.
            isefficient = PyCall.isbuftype(pyobj)
        end
        if isefficient
            info = PyArray_Info(pyobj)
            isefficient = info.native # byteorder is native
        end
        if isefficient
            # There is an efficient PyArray implementation
            NDArray{T, D, PyArray{T, D}}(pyobj)
        else
            # Fallback for other (e.g. string) types
            NDArray{T, D, PyObject}(pyobj)
        end
    end
    function NDArray(pyobj::PyObject)
        # Determine element type and rank
        T = julia_type(Datatype(pyobj.dtype))
        D = length(pyobj.shape)
        NDArray{T, D}(pyobj)
    end
end
tag2asdftype["tag:stsci.edu:asdf/core/ndarray-1.0.0"] = NDArray

const FastNDArray = NDArray{T, D, PyArray{T, D}} where {T <: PyArrayTypes, D}
const SlowNDArray = NDArray{T, D, PyObject} where {T, D}
pyarr(arr::FastNDArray{T, D}) where {T, D} = arr.repr::PyArray{T, D}
pyobj(arr::SlowNDArray{T, D}) where {T, D} = arr.repr::PyObject

isefficient(::Type{NDArray{T, D, PyArray{T, D}}}) where {T, D} = true
isefficient(::Type{NDArray{T, D, PyObject}}) where {T, D} = false
isefficient(::T) where {T <: NDArray} = isefficient(T)

# PyArray-based implementation
Base.axes(arr::FastNDArray) = axes(pyarr(arr))
Base.eachindex(ind::IndexCartesian, arr::FastNDArray) =
    eachindex(ind, pyarr(arr))
Base.eachindex(ind::IndexLinear, arr::FastNDArray) = eachindex(ind, pyarr(arr))
Base.getindex(arr::FastNDArray, i) = pyarr(arr)[i]
Base.length(arr::FastNDArray) = length(pyarr(arr))
Base.ndims(arr::FastNDArray) = ndims(pyarr(arr))
Base.size(arr::FastNDArray) = size(pyarr(arr))
Base.strides(arr::FastNDArray) = strides(pyarr(arr))

Base.IteratorSize(::Type{<:NDArray{T, D, PyArray{T, D}}}) where {T, D} =
    HasShape{D}()
Base.iterate(arr::FastNDArray, state...) = iterate(pyarr(arr), state...)

# PyObject-based implementation
function Base.axes(arr::SlowNDArray{T, D}) where {T, D}
    map(sz -> Base.Slice(0:sz-1),
        size(arr))::NTuple{D, Base.Slice{UnitRange{Int}}}
end
function Base.eachindex(::IndexCartesian, arr::SlowNDArray)
    CartesianIndices(axes(arr))
end
function Base.eachindex(::IndexLinear, arr::SlowNDArray)
    axes(arr, 1)
end
function Base.getindex(arr::SlowNDArray{T, D}, i::NTuple{D, Int}) where {T, D}
    @boundscheck @assert all(checkindex(Bool, axes(arr)[d], i[d]) for d in 1:D)
    pycall(pyobj(arr).__getitem__, T, i)::T
end
function Base.getindex(arr::SlowNDArray{T, D}, i::NTuple{D, I}) where {T, D, I}
    arr[NTuple{D, Int}(i)]
end
function Base.getindex(arr::SlowNDArray{T, D}, i::CartesianIndex{D}) where {
        T, D}
    arr[Tuple(i)]
end
function Base.getindex(arr::SlowNDArray, i...)
    arr[i]
end
function Base.ndims(arr::SlowNDArray{T, D}) where {T, D}
    D
end
function Base.size(arr::SlowNDArray{T, D}) where {T, D}
    NTuple{D, Int}(pyobj(arr).shape::NTuple{D, Int64})
end
function Base.strides(arr::SlowNDArray)
    Base.size_to_strides(1, reverse(size(arr))...)
end

Base.IteratorSize(::Type{<:SlowNDArray{T, D}}) where {T, D} = HasShape{D}()
function Base.iterate(arr::SlowNDArray, state...)
    iter = eachindex(arr)
    res = iterate(iter, state...)
    if res === nothing
        return nothing
    end
    idx, state = res
    arr[idx], state
end

################################################################################

const ASDFType = Union{Tree, Software, NDArray}

function Base.getindex(obj::ASDFType, key::String)
    additionalProperties(obj)   # check
    makeASDFType(get(obj.pyobj, PyObject, key))
end
function Base.get(obj::ASDFType, key::String, default)
    additionalProperties(obj)   # check
    value = get(obj.pyobj, PyObject, key, empty)
    if value === empty
        return default
    end
    makeASDFType(value)::ASDFType
end
function Base.setindex!(obj::ASDFType, value::ASDFType, key::String)
    additionalProperties(obj)   # check
    obj.pyobj[key] = value.pyobj
end
function Base.delete!(obj::ASDFType, key::String)
    additionalProperties(obj)   # check
    delete!(obj.pyobj, key)
    obj
end

function Base.length(obj::ASDFType)
    Int(obj.pyobj.__len__())
end
function Base.keys(obj::ASDFType)
    iter = obj.pyobj.keys().__iter__()
    keys = String[]
    while true
        try
            key = iter.__next__()
            push!(keys, key)
        catch
            # TODO: Check for Python StopIteration exception
            break
        end
    end
    Set(keys)
end
function Base.iterate(obj::ASDFType)
    Base.iterate(obj, obj.pyobj.__iter__())
end
function Base.iterate(obj::ASDFType, iter)
    try
        iter.__next__(), iter
    catch
        # TODO: Check for Python StopIteration exception
        nothing
    end
end

end


text:
__precompile__(true)
module Magpie

using SIUnits.ShortUnits

# dimensions dicts
include("dimensions/ball_bearings.jl")

end # module


text:
using Test

using LinearAlgebra
using SymbolicWedderburn
using SymbolicWedderburn.FiniteFields
using LinearAlgebra

using GroupsCore
using PermutationGroups
using Cyclotomics
using SparseArrays

include("smallgroups.jl")

@testset "Characters" begin
    import SymbolicWedderburn.Characters
    include("gf.jl")
    include("eigenspacedecomposition.jl")
    include("ccmatrix.jl")
    include("dixon.jl")
    include("characters.jl")
end
include("projections.jl")
include("sa_basis.jl")
include("action_permutation.jl")
include("action_linear.jl")
include("action_dihedral.jl")

if VERSION >= v"1.7.0" && !haskey(ENV, "CI")
    @testset "Examples" begin
        include("../examples/run_examples.jl")
    end
end


text:
println("---")

using KShiftsClustering, FunctionalData

d = 10
srand(0)
data = @p rand Float32 d 100_000
centers = @p rand Float32 d 100

labels = zeros(Int,len(data))
kshiftslabels!(labels, take(data,100), take(centers,2))
@time kshiftslabels!(labels, data, centers)
Profile.clear_malloc_data()
@time kshiftslabels!(labels, data, centers)
@time kshiftslabels!(labels, data, centers)
@time kshiftslabels!(labels, data, centers)
@time kshiftslabels!(labels, data, centers)


text:
const SlopeDefaultParams = Dict(
    :iterations => 10000,
    :verbosity  => 0, # Print intermediate info if > 0
    :optimIter  => 1,
    :gradIter   => 20,
    :tolInfeas  => 1e-6,
    :tolRelGap  => 1e-6,
    :xInit      => Float64[]
)

# Make sure the SLOPE code is downloaded and installncvreg R library is loaded
load_SLOPE_R_lib()

type SlopeRegression <: RegressionLearner
    params::Dict{Symbol, Any}
    SlopeRegression(; params...) = new(Dict(params))
end

function learn(r::SlopeRegression, dataset::DataSet)
    X = indepvars_as_row_major_matrix(dataset)
    y = depvar_as_vector(dataset)

    # Create params for this specific call based on default ones plus the ones for this
    # regression datum.
    default_params = copy(SlopeDefaultParams)
    params = merge(default_params, r.params)

    # Create lambda that needs to be as long as the num of coefficients we want to find
    numcoefs = size(X, 2)
    lambda = (1.0 / numcoefs) * collect(numcoefs:-1:1)

    # Now call Adlas and copy back the result as a Dict
    retval = R"Adlas($X, $y, $lambda, $params)"
    rcopy(retval)
end

function coefficients(r::SlopeRegression, dataset::DataSet; includeIntercept = true)
    coeffs = learn(r, dataset)[:x]
    includeIntercept ? coeffs : coeffs[2:end]
end


text:
using ViewerGL
GL = ViewerGL

points = rand(100,2)
GL.VIEW([
  GL.GLHull2d(points),
  GL.GLPoints(points),
  GL.GLFrame2
]);


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

makedocs(;
    modules=[ArDCA],
    authors="Andrea Pagnani, Jeanne Trinquier, Guido Uguzzoni, Martin Weigt, Francesco Zamponi",
    repo="https://github.com/pagnani/ArDCA/blob/{commit}{path}#L{line}",
    sitename="ArDCA",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pagnani.github.io/ArDCA",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/pagnani/ArDCA",
)

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
module Day08

import ..data_dir # from parent module
input = readlines(joinpath(data_dir, "day08"))

export part1, part2

"""

  0:      1:      2:      3:      4:
 aaaa    ....    aaaa    aaaa    ....
b    c  .    c  .    c  .    c  b    c
b    c  .    c  .    c  .    c  b    c
 ....    ....    dddd    dddd    dddd
e    f  .    f  e    .  .    f  .    f
e    f  .    f  e    .  .    f  .    f
 gggg    ....    gggg    gggg    ....

  5:      6:      7:      8:      9:
 aaaa    aaaa    aaaa    aaaa    aaaa
b    .  b    .  .    c  b    c  b    c
b    .  b    .  .    c  b    c  b    c
 dddd    dddd    ....    dddd    dddd
.    f  e    f  .    f  e    f  .    f
.    f  e    f  .    f  e    f  .    f
 gggg    gggg    ....    gggg    gggg

"""

function split_display(display::AbstractString)
    raw_patterns, raw_output = split(display, " | ")
    patterns = [Set(string.(split(p, ""))) for p in split(raw_patterns, " ")]
    outputs = [Set(string.(split(p, ""))) for p in split(raw_output, " ")]
    patterns, outputs
end

function patterns_of_length(patterns::Vector{Set{String}}, l::Int64)
    filter(x -> length(x) == l, patterns)
end

function identify_configuration(patterns::Vector{Set{String}})
    # these are uniquely identifiable from the number of segments
    one = first(patterns_of_length(patterns, 2))
    four = first(patterns_of_length(patterns, 4))
    seven = first(patterns_of_length(patterns, 3))
    eight = first(patterns_of_length(patterns, 7))

    # There are three 6-segment elements: six, nine, zero
    patterns_of_length_6 = patterns_of_length(patterns, 6)
    # four is a subset of nine but not six or zero
    nine = filter(s -> four ⊆ s, patterns_of_length_6) |> first
    # six is the only 6-segment digit that is not a superset of one
    six = filter(s -> one ⊈ s, patterns_of_length_6) |> first
    # and by process of elimination
    zero = filter(s -> s ∉ [six, nine], patterns_of_length_6) |> first

    # The three remaining elements are of 5 segments
    patterns_of_length_5 = patterns_of_length(patterns, 5)
    # three overlaps with one
    three = filter(s -> one ⊆ s, patterns_of_length_5) |> first
    # five is a subset of six
    five = filter(s -> s ⊆ six, patterns_of_length_5) |> first
    # and the last is two
    two = filter(s -> s ∉ [three, five], patterns_of_length_5) |> first

    Dict(
        zero => 0,
        one => 1,
        two => 2,
        three => 3,
        four => 4,
        five => 5,
        six => 6,
        seven => 7,
        eight => 8,
        nine => 9
    )
end

function identify_digits_in_output(display::AbstractString)
    patterns, outputs = split_display(display)
    configuration = identify_configuration(patterns)
    [configuration[output] for output in outputs]
end

function part1(input = input)
    digits_in_outputs = identify_digits_in_output.(input)
    target_digits = [1, 4, 7, 8]
    counts = [length(filter(x -> x ∈ target_digits, digits)) for digits in digits_in_outputs]
    sum(counts)
end

function glue_digits(digits::Vector{Int64})
    n_digits = length(digits)
    sum([digits[i] * 10^(n_digits - i) for i = 1:n_digits])
end

function part2(input = input)
    digits_in_outputs = identify_digits_in_output.(input)
    as_numbers = glue_digits.(digits_in_outputs)
    sum(as_numbers)
end

end #module


text:
# # Hyperelasticity
#
# **Keywords**: *hyperelasticity*, *finite strain*, *large deformations*, *Newton's method*,
# *conjugate gradient*, *automatic differentiation*
#
# ![hyperelasticity.png](hyperelasticity.png)
#
# *Figure 1*: Cube loaded in torsion modeled with a hyperelastic material model and
# finite strain.
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`hyperelasticity.ipynb`](@__NBVIEWER_ROOT_URL__/examples/hyperelasticity.ipynb).
#-
# ## Introduction
#
# In this example we will solve a problem in a finite strain setting using an
# hyperelastic material model. In order to compute the stress we will use automatic
# differentiation, to solve the non-linear system we use Newton's
# method, and for solving the Newton increment we use conjugate gradients.
#
# The weak form is expressed in terms of the first Piola-Kirchoff stress ``\mathbf{P}``
# as follows: Find ``u \in \mathbb{U}`` such that
#
# ```math
# \int_{\Omega} \nabla_{\mathbf{X}} \delta \mathbf{u} : \mathbf{P}(\mathbf{u})\ \mathrm{d}\Omega =
# \int_{\Omega} \delta \mathbf{u} \cdot \mathbf{b}\ \mathrm{d}\Omega + \int_{\Gamma^\mathrm{N}}
# \delta \mathbf{u} \cdot \mathbf{t}\ \mathrm{d}\Gamma
# \quad \forall \delta \mathbf{u} \in \mathbb{U}^0,
# ```
#
# where ``\mathbf{u}`` is the unknown displacement field, ``\mathbf{b}`` is the body force acting
# on the reference domain, ``\mathbf{t}`` is the traction acting on the Neumann part of the reference
# domain's boundary, and where ``\mathbb{U}`` and ``\mathbb{U}^0`` are suitable trial and test sets.
# ``\Omega`` denotes the reference domain, which is also called reference or material domain.
# Gradients are defined with respect to the reference domain, here denoted with  an ``X``.
# Formally this is expressed as ``(\nabla_{\mathbf{X}} \bullet)_{ij} := \frac{\partial(\bullet)_i}{\partial X_j}``.
# Note that for large deformation problems it is also possibile that gradients and integrals
# are defined on the deformed domain, which is also called the current or spatial domain, depending
# on the specific formulation.
#

using Ferrite, Tensors, TimerOutputs, ProgressMeter
import KrylovMethods, IterativeSolvers

# ## Hyperelastic material model
#
# The stress can be derived from an energy potential, defined in
# terms of the right Cauchy-Green tensor ``\mathbf{C} = \mathbf{F}^{\mathrm{T}} \mathbf{F}``,
# where ``\mathbf{F} = \mathbf{I} + \nabla_{\mathbf{X}} \mathbf{u}`` is the deformation gradient.
# We shall use a neo-Hookean model, where the potential can be written as
#
# ```math
# \Psi(\mathbf{C}) = \frac{\mu}{2} (I_1 - 3) - \mu \ln(J) + \frac{\lambda}{2} \ln(J)^2,
# ```
#
# where ``I_1 = \mathrm{tr}(\mathbf{C})`` is the first invariant, ``J = \sqrt{\det(\mathbf{C})}``
# and ``\mu`` and ``\lambda`` material parameters.
# From the potential we obtain the second Piola-Kirchoff stress ``\mathbf{S}`` as
#
# ```math
# \mathbf{S} = 2 \frac{\partial \Psi}{\partial \mathbf{C}},
# ```
#
# and the tangent of ``\mathbf{S}`` as
#
# ```math
# \frac{\partial \mathbf{S}}{\partial \mathbf{C}} = 4 \frac{\partial \Psi}{\partial \mathbf{C}}.
# ```
#
# We can implement the material model as follows, where we utilize automatic differentiation
# for the stress and the tangent, and thus only define the potential:

struct NeoHooke
    μ::Float64
    λ::Float64
end

function Ψ(C, mp::NeoHooke)
    μ = mp.μ
    λ = mp.λ
    Ic = tr(C)
    J = sqrt(det(C))
    return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2
end

function constitutive_driver(C, mp::NeoHooke)
    ## Compute all derivatives in one function call
    ∂²Ψ∂C², ∂Ψ∂C = Tensors.hessian(y -> Ψ(y, mp), C, :all)
    S = 2.0 * ∂Ψ∂C
    ∂S∂C = 4.0 * ∂²Ψ∂C²
    return S, ∂S∂C
end;

# Finally, for the finite element problem we need ``\mathbf{P}`` and
# ``\frac{\partial \mathbf{P}}{\partial \mathbf{F}}``, which can be
# obtained by using the following relations:
#
# ```math
# \begin{align*}
# \mathbf{P} &= \mathbf{F} \cdot \mathbf{S},\\
# \frac{\partial \mathbf{P}}{\partial \mathbf{F}} &= [\mathbf{F} \bar{\otimes} \mathbf{I}] :
# \frac{\partial \mathbf{S}}{\partial \mathbf{C}} : [\mathbf{F}^\mathrm{T} \bar{\otimes} \mathbf{I}]
# + \mathbf{I} \bar{\otimes} \mathbf{S}.
# \end{align*}
# ```
#
# ## Newton's Method
#
# As mentioned above, to deal with the non-linear weak form we first linearize
# the problem such that we can apply Newton's method, and then apply the FEM to
# discretize the problem. Skipping a detailed derivation, Newton's method can
# be expressed as:
# Given some initial guess ``\mathbf{u}^0``, find a sequence ``\mathbf{u}^{k}`` by iterating
#
# ```math
# \mathbf{u}^{k+1} = \mathbf{u}^{k} - \Delta \mathbf{u}^{k}
# ```
#
# until some termination condition has been met. Therein we determine ``\Delta \mathbf{u}^{k}``
#
# ```math
# \mathbf{K}(\mathbf{u}^{k}) \Delta \mathbf{u}^{k} = \mathbf{g}(\mathbf{u}^{k})
# ```
#
# where $\mathbf{g}$ is the global residual and $\mathbf{K} = \frac{\partial g}{\partial u}$
# the Jacobi matrix, such that
#
# ```math
# K_{ij} = \int_{\Omega} \nabla_{\mathbf{X}} \delta \mathbf{u}_{i} : \frac{\partial \mathbf{P}}{\partial \mathbf{F}} : \nabla_{\mathbf{X}} \delta \mathbf{u}_{j} \, \mathrm{d} \Omega
# ```
# and
#
# ```math
# g_{i} = \int_{\Omega} \nabla_{\mathbf{X}} \delta \mathbf{u}_{i} : \mathbf{P} - \delta \mathbf{u}_{i} \cdot \mathbf{b} \, \mathrm{d} \Omega
# ```
#
# A detailed derivation can be found in every continuum mechanics book, which has a
# chapter about finite elasticity theory. We used "Nonlinear solid mechanics: a continuum
# approach for engineering science." by Gerhard Holzapfel (chapter 8) as a reference.
#
# ## Finite element assembly
#
# The element routine for assembling the residual and tangent stiffness is implemented
# as usual, with loops over quadrature points and shape functions:

function assemble_element!(ke, ge, cell, cv, fv, mp, ue)
    ## Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    b = Vec{3}((0.0, -0.5, 0.0)) # Body force
    t = Vec{3}((0.1, 0.0, 0.0)) # Traction
    ndofs = getnbasefunctions(cv)

    for qp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, qp)
        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇u = function_gradient(cv, qp, ue)
        F = one(∇u) + ∇u
        C = tdot(F)
        ## Compute stress and tangent
        S, ∂S∂C = constitutive_driver(C, mp)
        P = F ⋅ S
        I = one(S)
        ∂P∂F = otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I) + otimesu(I, S)

        ## Loop over test functions
        for i in 1:ndofs
            ## Test function + gradient
            δui = shape_value(cv, qp, i)
            ∇δui = shape_gradient(cv, qp, i)
            ## Add contribution to the residual from this test function
            ge[i] += ( ∇δui ⊡ P - δui ⋅ b ) * dΩ

            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
            for j in 1:ndofs
                ∇δuj = shape_gradient(cv, qp, j)
                ## Add contribution to the tangent
                ke[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
            end
        end
    end

    ## Surface integral for the traction
    for face in 1:nfaces(cell)
        if onboundary(cell, face)
            reinit!(fv, cell, face)
            for q_point in 1:getnquadpoints(fv)
                dΓ = getdetJdV(fv, q_point)
                for i in 1:ndofs
                    δui = shape_value(fv, q_point, i)
                    ge[i] -= (δui ⋅ t) * dΓ
                end
            end
        end
    end
end;

# Assembling global residual and tangent
function assemble_global!(K, f, dh, cv, fv, mp, u)
    n = ndofs_per_cell(dh)
    ke = zeros(n, n)
    ge = zeros(n)

    ## start_assemble resets K and f
    assembler = start_assemble(K, f)

    ## Loop over all cells in the grid
    @timeit "assemble" for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        ue = u[global_dofs] # element dofs
        @timeit "element assemble" assemble_element!(ke, ge, cell, cv, fv, mp, ue)
        assemble!(assembler, global_dofs, ge, ke)
    end
end;

# Define a main function, with a loop for Newton iterations

function solve()
    reset_timer!()

    ## Generate a grid
    N = 10
    L = 1.0
    left = zero(Vec{3})
    right = L * ones(Vec{3})
    grid = generate_grid(Tetrahedron, (N, N, N), left, right)

    ## Material parameters
    E = 10.0
    ν = 0.3
    μ = E / (2(1 + ν))
    λ = (E * ν) / ((1 + ν) * (1 - 2ν))
    mp = NeoHooke(μ, λ)

    ## Finite element base
    ip = Lagrange{3, RefTetrahedron, 1}()
    qr = QuadratureRule{3, RefTetrahedron}(1)
    qr_face = QuadratureRule{2, RefTetrahedron}(1)
    cv = CellVectorValues(qr, ip)
    fv = FaceVectorValues(qr_face, ip)

    ## DofHandler
    dh = DofHandler(grid)
    push!(dh, :u, 3) # Add a displacement field
    close!(dh)

    function rotation(X, t, θ = deg2rad(60.0))
        x, y, z = X
        return t * Vec{3}(
            (0.0,
            L/2 - y + (y-L/2)*cos(θ) - (z-L/2)*sin(θ),
            L/2 - z + (y-L/2)*sin(θ) + (z-L/2)*cos(θ)
            ))
    end

    dbcs = ConstraintHandler(dh)
    ## Add a homogenous boundary condition on the "clamped" edge
    dbc = Dirichlet(:u, getfaceset(grid, "right"), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3])
    add!(dbcs, dbc)
    dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> rotation(x, t), [1, 2, 3])
    add!(dbcs, dbc)
    close!(dbcs)
    t = 0.5
    Ferrite.update!(dbcs, t)

    ## Pre-allocation of vectors for the solution and Newton increments
    _ndofs = ndofs(dh)
    un = zeros(_ndofs) # previous solution vector
    u  = zeros(_ndofs)
    Δu = zeros(_ndofs)
    ΔΔu = zeros(_ndofs)
    apply!(un, dbcs)

    ## Create sparse matrix and residual vector
    K = create_sparsity_pattern(dh)
    g = zeros(_ndofs)

    ## Perform Newton iterations
    newton_itr = -1
    NEWTON_TOL = 1e-8
    prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:")

    while true; newton_itr += 1
        u .= un .+ Δu # Current guess
        assemble_global!(K, g, dh, cv, fv, mp, u)
        normg = norm(g[Ferrite.free_dofs(dbcs)])
        apply_zero!(K, g, dbcs)
        ProgressMeter.update!(prog, normg; showvalues = [(:iter, newton_itr)])

        if normg < NEWTON_TOL
            break
        elseif newton_itr > 30
            error("Reached maximum Newton iterations, aborting")
        end

        ## Compute increment using cg! from IterativeSolvers.jl
        @timeit "linear solve (KrylovMethods.cg)" ΔΔu′, flag, relres, iter, resvec = KrylovMethods.cg(K, g; maxIter = 1000)
        @assert flag == 0
        @timeit "linear solve (IterativeSolvers.cg!)" IterativeSolvers.cg!(ΔΔu, K, g; maxiter=1000)

        apply_zero!(ΔΔu, dbcs)
        Δu .-= ΔΔu
    end

    ## Save the solution
    @timeit "export" begin
        vtk_grid("hyperelasticity", dh) do vtkfile
            vtk_point_data(vtkfile, dh, u)
        end
    end

    print_timer(title = "Analysis with $(getncells(grid)) elements", linechars = :ascii)
    return u
end

# Run the simulation

u = solve();

## test the result                #src
using Test                        #src
@test norm(u) ≈ 4.865189736192834 #src

#md # ## Plain program
#md #
#md # Here follows a version of the program without any comments.
#md # The file is also available here: [`hyperelasticity.jl`](hyperelasticity.jl).
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using Distributions

@testset "gains" begin
    _,vis, amp, lcamp, cphase = load_data()

    st = scantable(vis)
    stamp = scantable(amp)
    stcp = scantable(cphase)
    stlca = scantable(lcamp)

    gcache = GainCache(st)

    θ = (f1 = 1.0,
         σ1 = 1.2e-10,
         τ1 = 0.25,
         ξ1 = -0.42,
         f2 = 0.5,
         σ2 = 8.0e-12,
         τ2 = 0.4,
         ξ2 = 0.0,
         x = 4.0e-11,
         y = 1.5e-10
        )

    m = test_model(θ)

    tel = stations(vis)
    gamp_prior = NamedTuple{Tuple(tel)}(ntuple(_->LogNormal(0.0, 0.1), length(tel)))
    gph_prior = NamedTuple{Tuple(tel)}(ntuple(_->Normal(0.0, 0.1), length(tel)))

    gamp = GainPrior(gamp_prior, st)
    gpha = GainPrior(gph_prior, st)

    ga = fill(1.0, size(rand(gamp))...)
    gm = GainModel(gcache, ga, m)

    ac = arrayconfig(vis)
    cpac = arrayconfig(cphase)
    lcac = arrayconfig(lcamp)
    @test visibilities(gm, ac) ≈ visibilities(m, ac)
    @test amplitudes(gm, ac) ≈ amplitudes(m, ac)
    @test closure_phases(gm, cpac) ≈ closure_phases(m, cpac)
    @test logclosure_amplitudes(gm, lcac) ≈ logclosure_amplitudes(m, lcac)

    rga = rand(gamp)
    rph = rand(gpha)
    @inferred logdensityof(gamp, rga)
    @inferred logdensityof(gpha, rph)
end


text:
# Autogenerated wrapper script for GTK4_jll for x86_64-unknown-freebsd
export libgtk4

using Glib_jll
using Graphene_jll
using Cairo_jll
using Pango_jll
using FriBidi_jll
using FreeType2_jll
using gdk_pixbuf_jll
using Libepoxy_jll
using HarfBuzz_jll
using xkbcommon_jll
using iso_codes_jll
using Wayland_jll
using Wayland_protocols_jll
using Xorg_libXrandr_jll
using Xorg_libX11_jll
using Xorg_libXrender_jll
using Xorg_libXi_jll
using Xorg_libXext_jll
using Xorg_libXcursor_jll
using Xorg_libXdamage_jll
using Xorg_libXfixes_jll
using Xorg_libXinerama_jll
using Fontconfig_jll
JLLWrappers.@generate_wrapper_header("GTK4")
JLLWrappers.@declare_library_product(libgtk4, "libgtk-4.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, Graphene_jll, Cairo_jll, Pango_jll, FriBidi_jll, FreeType2_jll, gdk_pixbuf_jll, Libepoxy_jll, HarfBuzz_jll, xkbcommon_jll, iso_codes_jll, Wayland_jll, Wayland_protocols_jll, Xorg_libXrandr_jll, Xorg_libX11_jll, Xorg_libXrender_jll, Xorg_libXi_jll, Xorg_libXext_jll, Xorg_libXcursor_jll, Xorg_libXdamage_jll, Xorg_libXfixes_jll, Xorg_libXinerama_jll, Fontconfig_jll)
    JLLWrappers.@init_library_product(
        libgtk4,
        "lib/libgtk-4.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import ViscousFlow: Layers

Δx = 0.02
xlim = (-5.98,5.98)
ylim = (-5.98,5.98)
g = Fields.PhysicalGrid(xlim,ylim,Δx)

w = Nodes(Dual,size(g))
dq = Edges(Dual,w)
n = 150
radius = 1.0
body = Circle(radius,n)
X = VectorData(collect(body))
f = VectorData(X)
ϕ = ScalarData(f)
fill!(ϕ,1.0)

regop = Regularize(X,Δx,I0=origin(g),issymmetric=true,ddftype=Fields.Yang3)
Hv, Ev = RegularizationMatrix(regop,f,dq)
Hs, Es = RegularizationMatrix(regop,ϕ,w)

@testset "Double Layer" begin

  dlayer = DoubleLayer(body,Hv,weight=1/Δx)
  ϕ .= rand(length(ϕ))

  @test abs(sum(dlayer(ϕ))) < 100*eps(1.0)

  dlayer2 = DoubleLayer(body,regop,w)

  @test dlayer2(ϕ) == dlayer(ϕ)

end

@testset "Single Layer" begin

  slayer = SingleLayer(body,Hs)
  ϕ .= 1.0

  @test abs(sum(slayer(ϕ))-2π*radius) < 2e-3

  slayer2 = SingleLayer(body,regop,w)

  @test slayer2(ϕ) == slayer(ϕ)

end

@testset "Masks" begin

  inner = Mask(body,regop,w)

  fill!(w,1)

  @test abs(sum(inner(w))*cellsize(g)^2 - π*radius^2) < 1e-3

  outer = ComplementaryMask(inner)

end


text:
using Minesweeper

import_gui()

# simple model that selects random location and flags it with probablity θflag
# and clears it with probability 1-θflag
struct Model
    θflag::Float64
end

function select_random_cell(game)
    unselected = findall(x->!x.revealed && !x.flagged, game.cells)
    return rand(unselected)
end

function flag!(game, choice)
    game.cells[choice].flagged = true
end

# simulate game
function run!(game; realtime=false, gui=nothing, show_gui=false, pause=.2)
    detonated = false
    while !game.mine_detonated && !game_over(game)
        choice = select_random_cell(game)
        if rand() <= model.θflag
            flag!(game, choice)
        else
            select_cell!(game, choice)
        end
        show_gui ? update_reveal!(game, gui) : nothing
        realtime ? sleep(pause) : nothing
        game.trials += 1
    end
    compute_score!(game)
    return nothing
end

# run with GUI
model = Model(.5)
game = Game(dims=(10,10), n_mines=15)
gui = generate_gui(game)

start_new_game!(gui, game)
run!(game; gui, realtime=true, show_gui=true)
println(game.score)

#Run without GUI
model = Model(.5)
game = Game(dims=(10,10), n_mines=15)
run!(game)
println(game.score)


text:
"""
    The function that compute the cost of a single output value.
"""
single_cost(netout, expout) = ((netout - expout)^2)/2;
"""
    The prime derivative of the single_cost function.
"""
single_cost_prime(netout, expout) = (netout - expout);

function cost(net::Network, testdatum::TestDatum)
    @assert outlength(net) == outlength(testdatum);
    @assert inlength(net) == inlength(testdatum);
    return sum(single_cost.(output(net),output(testdatum)));
end

function cost!(net::Network, testdatum::TestDatum)
    feedforward!(net,input(testdatum));
    return cost(net,testdatum);
end

function cost!(net::Network, testdata::Vector{TestDatum})
    costsum = 0.0;
    for test in testdata
        costsum += cost!(net,test);
    end
    return costsum/length(testdata);
end


text:
using KernelDensityDistributionEsimates
using Test

@testset "KernelDensityDistributionEsimates.jl" begin
    # Write your own tests here.
end


text:
"""
    CPInterval <: ConfidenceInterval

Naive Clopper-Pearson confidence interval using default ordering.
"""
struct CPInterval{TD<:Design,TR<:Real} <: ConfidenceInterval

  design::TD
  confidence::TR

  function CPInterval{TD,TR}(
    design::TD, confidence::TR
  ) where {TD<:Design,TR<:Real}

    @checkprob confidence
    new(design, confidence)

  end # inner constructor

end # CPInterval

"""
    CPInterval{T<:Real}(
      design::Design;
      confidence::T = .9
    )

Naive Clopper-Pearson confidence interval using default ordering.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | binary two-stage design |
| confidence   | confidence level of the interval |
"""
function CPInterval(
  design::TD; confidence::Real = .9
) where {TD<:Design}

  return CPInterval{TD,typeof(confidence)}(design, confidence)

end # CPInterval

Base.show(io::IO, ci::CPInterval) = print("CPInterval")

function limits(ci::CPInterval, x1::T, x2::T) where {T<:Integer}

    d              = design(ci)
    ispossible(d, x1, x2) ? nothing : error("(x1, x2) not compatible with underlying design")
    alpha::Float64 = 1 - confidence(ci)
    n1             = interimsamplesize(d)
    nmax::Int64    = maximum(samplesize(d))
    x              = x1 + x2
    n              = samplesize(d, x1)
    if x == 0
        limits = [0.0; (1.0 - (alpha/2)^(1/n))]
    elseif x == n
        limits = [(alpha/2)^(1/n); 1.0]
    else
        limits = [quantile(Distributions.Beta(x, n - x + 1), alpha/2); quantile(Distributions.Beta(x + 1, n - x), 1 - alpha/2)]
    end
    return limits

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
struct IdealGas
    AdiabaticIndex
end

function calcSoundSpeed(eos::IdealGas,
                        d::Float64,
                        p::Float64)
    g = eos.AdiabaticIndex
    return sqrt(g*p/d)
end

function calcSpecificThermalEnergy(eos::IdealGas,
                                   d::Float64,
                                   p::Float64)
    g = eos.AdiabaticIndex
    return p/d/(g-1)
end

function calcPressure(eos::IdealGas,
                      d::Float64,
                      e::Float64)
    g = eos.AdiabaticIndex
    return (g-1)*e*d
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
"""
$(TYPEDEF)

Abstract representation of a finite horizon game.
"""
abstract type AbstractGame end

"""
    $(FUNCTIONNAME)(g::AbstractGame)

Returns the dynamics of the game, a `ControlSystem`.
"""
function dynamics end

"""
    $(FUNCTIONNAME)(g::AbstractGame)

Returns the cost representation for the game (a vector of PlayerCost) for each
player.
"""
function player_costs end

# TODO add documentation
function uindex end
function horizon end

n_players(g::AbstractGame) = length(uindex(g))

"""
    $(TYPEDSIGNATURES)

Returns the type of the strategy that is a solution to this game.
"""
# TODO: I would really prefer if we did not have to use this!
function strategytype(g::AbstractGame)
    elt = AffineStrategy{n_states(g), n_controls(g),
                         SMatrix{n_controls(g), n_states(g), Float64,
                                 n_controls(g)*n_states(g)},
                         SVector{n_controls(g), Float64}}

    return SizedArray{Tuple{horizon(g)}, elt, 1, 1}
end

"""
    $(TYPEDSIGNATURES)

Returns the type of the state that is used in this game.
"""
statetype(g::AbstractGame) = SVector{n_states(g), Float64}

"""
    $(FUNCTIONNAME)(g::AbstractGame)

Returns the lq approximation of the game.
"""
function lq_approximation end

# delegate some function calls to the dynamics
n_states(g::AbstractGame) = n_states(dynamics(g))
n_controls(g::AbstractGame) = n_controls(dynamics(g))
xindex(g::AbstractGame) = xindex(dynamics(g))
xyindex(g::AbstractGame) = xyindex(dynamics(g))
samplingtime(g::AbstractGame) = samplingtime(dynamics(g))

# additional convenience methods
time_disc2cont(g::AbstractGame, k::Int, t0::Float64=0.) = (t0 +
                                                           (k-1)*samplingtime(g))
pindex(g) = @S 1:n_players(g)


text:
using HeartRateVariability
using Test

n=HeartRateVariability.infile("e1304.txt")
td=HeartRateVariability.time_domain(n)
fd=HeartRateVariability.frequency(n)
nl=HeartRateVariability.nonlinear(n)
g=HeartRateVariability.geometric(n)

@testset "HeartRateVariability.jl" begin

    @testset "HeartRateVariability.infile" begin
        @test HeartRateVariability.infile("e1304","atr")==n
    end

    @testset "HeartRateVariability.time_domain" begin
        @test td.mean≈917.24 atol=0.1
        @test td.sdnn≈137.19 atol=0.1
        @test td.rmssd≈27.85 atol=0.1
        @test td.sdsd≈20.35 atol=0.1
        @test td.nn50≈342 atol=1
        @test td.pnn50≈4.41 atol=0.1
        @test td.nn20≈2831 atol=1
        @test td.pnn20≈36.53 atol=0.1
        @test td.rRR≈2.67 atol=0.1
    end

    @testset "HeartRateVariability.frequency" begin
        @test fd.vlf≈1317.96 atol=0.01*fd.vlf
        @test fd.lf≈90.36 atol=0.01*fd.lf
        @test fd.hf≈176.05 atol=0.01*fd.hf
        @test fd.lfhf_ratio≈0.51 atol=0.01*fd.lfhf_ratio
        @test fd.tp≈1584.35 atol=0.01*fd.tp
    end

    @testset "HeartRateVariability.nonlinear" begin
        @test nl.apen≈2.16 atol=0.1
        @test nl.sampen≈2.16 atol=0.1
        @test nl.hurst≈0.37 atol=0.1
        @test nl.renyi0≈-6.82 atol=0.1
        @test nl.renyi1≈-6.83 atol=0.1
        @test nl.renyi2≈-6.84 atol=0.1

        #testing if get_rs from module Nonlinear returns 0 when S or R is 0
        @test HeartRateVariability.Nonlinear.get_rs(ones(100))==0

        #testing if warning is thwown
        @test_logs (:warn,"To obtain a valid value for the hurst coefficient, the length of the data series must be greater than or equal to 100.") HeartRateVariability.nonlinear([1.0,2.0,1.0,2.0,1.0,2.0,1.0,2.0,1.0,2.0])
    end

    @testset "HeartRateVariability.geometric" begin
        @test g.poincare!=nothing
        @test g.recurrence!=nothing
        @test_throws ErrorException HeartRateVariability.geometric(n,"error")
    end
end


text:
############################################################
# FactCheck.jl
# A testing framework for Julia
# http://github.com/JuliaLang/FactCheck.jl
# MIT Licensed
############################################################

module TestFactCheck

using FactCheck
using Test

############################################################
# Before we excerse the other various parts of FactCheck,
# check we actually catch and report errors correctly. This
# also allows us to test printing code for the Failure and
# Error cases, which wouldn't be tested otherwise.
print("Testing Result counting and printing, not actual errors!\n")
facts("Test error pathways") do
    a_success = @fact 1 --> 1 "I will never be seen"
    println(a_success)
    a_failure = @fact 1 --> 2 "one doesn't equal two!"
    a_error   = @fact sqrt(-1) --> im "domains are tricky"
    a_pending = @pending not_really_pending() "sorta pending"
    println(a_pending)
end
# stats = getstats()
# FactCheck.clear_results()
# @test stats["nSuccesses"] == 1
# @test stats["nFailures"] == 1
# @test stats["nErrors"] == 1
# @test stats["nPending"] == 1
# @test stats["nNonSuccessful"] == 2
# print("Done, begin actual FactCheck tests\n")

############################################################
# Begin actual tests
mutable struct Foo a end
mutable struct Bar a end
mutable struct Baz end
mutable struct Bazz a end
import Base.(==)
==(x::Foo, y::Foo) = x.a == y.a

struct MyError <: Exception
end

module MyModule
    struct MyError <: Exception
    end
end

facts("Testing core functionality") do
    @fact 1 --> 1
    @fact 2*2 --> 4
    @fact uppercase("foo") --> "FOO"
    @fact_throws 2^-1
    @fact_throws 2^-1 "a domain error happend"
    @fact_throws DomainError 2^-1
    @fact_throws DomainError 2^-1 "a domain error happened"
    @fact_throws MyError throw(MyError())
    @fact_throws MyError throw(MyError()) "my error happend"
    @fact_throws MyModule.MyError throw(MyModule.MyError())
    @fact_throws MyModule.MyError throw(MyModule.MyError()) "my error happend"
    @fact 2*[1,2,3] --> [2,4,6]
    @fact Foo(1) --> Foo(1)
    if VERSION >= v"0.5-"
        hmm = function()
            2
            3
            FactCheck.getline()
        end
        @fact hmm() --> 72
    end
end

facts("Testing invalid @fact_throws macro") do
    @fact_throws ArgumentError eval(:(@fact_throws "this needs to be an expression"))
    @fact_throws ArgumentError eval(:(@fact_throws "wrong type" :wrong_type))
    @fact_throws ArgumentError eval(:(@fact_throws "wrong type" error("this is an error")))
end

facts("Testing 'context'") do
    # FactCheck.LEVEL starts from 0
    @fact FactCheck.LEVEL --> 0

    context("context will increase LEVEL and set contexts") do
        @fact FactCheck.LEVEL --> 1
        @fact FactCheck.contexts[end] --> "context will increase LEVEL and set contexts"
    end

    @fact FactCheck.LEVEL --> 0

    # context called without 'desc' will still increase LEVEL
    context() do
        @fact FactCheck.LEVEL --> 1
    end

    context("nested context") do
        @fact FactCheck.LEVEL --> 1
        @fact FactCheck.contexts[end] --> "nested context"

        context("inner") do
            @fact FactCheck.LEVEL --> 2
            @fact FactCheck.contexts[end] --> "inner"
        end
    end

    facts("'facts' doesn't increase LEVEL") do
        @fact FactCheck.LEVEL --> 0
    end

    context("will execute the function which is passed to the 'context'") do
        executed = false
        f() = (executed = true)

        @fact executed --> false
        context(f, "Run f")
        @fact executed --> true
    end

    context("indent by current LEVEL") do
        original_STDOUT = stdout
        (out_read, out_write) = redirect_stdout()
        system_output = @async readstring(out_read)

        context("intended") do
            close(out_write)
            close(out_read)

            redirect_stdout(original_STDOUT)
            # current LEVEL is 2
            expected_str = string('\t'^2,"> intended\n")
            @fact wait(system_output) --> expected_str
        end
    end
end

facts("FactCheck assertion helper functions") do

    context("`not` works for values and functions") do
        notone = not(1)
        @fact notone(2) --> true
        @fact notone(1) --> false
        @fact 2 --> not(1)

        noteven = not(iseven)
        @fact noteven(3) --> true
        @fact noteven(2) --> false
        @fact not(iseven)(2) --> false
        @fact 3 --> not(iseven)
    end

    context("`exactly` can be used to check object equality") do
        @fact exactly(exactly)(exactly) --> true

        x() = ()
        @fact exactly(x)(x) --> true

        # types with no fields return a singleton object when instantiated
        @fact exactly(Baz())(Baz()) --> true

        @fact exactly(Bazz(1))(Bazz(1)) --> false
    end

    context("`roughly` compares numbers... roughly") do
        @fact 2.4999999999999 --> roughly(2.5)
        @fact 9.5 --> roughly(10; atol=1.0)
        @fact 10.5 --> roughly(10; atol=1.0)
        @fact 10.5 --> roughly(10, 1.0)
    end

    context("`roughly` compares matrixes... roughly") do
        X = [1.1 1.2; 2.1 2.2]
        Y = X + [+1e-8 -1e-8; -1e-8 +1e-8]
        Z = [1 1; 2 2]
        @fact X --> roughly(Y)
        @fact X --> roughly(Z; atol=0.4)
        @fact X --> roughly(Z, 0.4)
    end

    context("`anyof` compares with all arguments") do
        @fact 2+2 --> anyof(4, :four, "four")
        @fact 5   --> not(anyof(:five, "five"))
    end

    context("`anyof` works for functions") do
        @fact 5 --> anyof(5.1, roughly(5.1,0.01), roughly(4.9,0.2))
    end

    context("less_than") do
        @fact 1 --> less_than(2)
    end

    context("less_than_or_equal") do
        @fact 1 --> less_than_or_equal(2)
        @fact 1 --> less_than_or_equal(1)
    end

    context("greater_than") do
        @fact 2 --> greater_than(1)
    end

    context("greater_than_or_equal") do
        @fact 2 --> greater_than_or_equal(1)
        @fact 2 --> greater_than_or_equal(2)
    end
end

# exitstatus()

end # module


text:
function sucr_model(day::Int64, df_gens::DataFrame,
    df_peakers::DataFrame, df_storage::DataFrame, fuelprice::Matrix{Float64},
    fuelprice_peakers::Matrix{Float64},
	map_gens::Matrix{Int64}, map_storage::Matrix{Int64},
	sucr_results::UC_Results, sucr_prepoc::DemandWAPreprocGens,
	lb_MDT::Matrix{Int64}, lb_MUT::Matrix{Int64},
	lb_MDT_peaker::Matrix{Int64}, lb_MUT_peaker::Matrix{Int64})

	t1_SUCRmodel = time_ns()

	uc_gens = GensResults();
	uc_peakers = PeakersResults(
		genout_block = zeros(Float64, PEAKERS,BLOCKS,24-INIT_HR_SUCR+INIT_HR_FUCR));
	uc_storg = StorageResults();
	uc_results = UC_Results(uc_gens, uc_peakers, uc_storg);

	SUCRmodel = direct_model(CPLEX.Optimizer())
	set_optimizer_attribute(SUCRmodel, "CPX_PARAM_EPGAP", SOLVER_EPGAP)

	# Declaring the decision variables for conventional generators
	@variable(SUCRmodel, sucrm_genOnOff[1:GENS, 0:HRS_SUCR], Bin) #Bin
	@variable(SUCRmodel, sucrm_genStartUp[1:GENS, 1:HRS_SUCR], Bin) # startup variable
	@variable(SUCRmodel, sucrm_genShutDown[1:GENS, 1:HRS_SUCR], Bin) # shutdown variable
	@variable(SUCRmodel, sucrm_genOut[1:GENS, 0:HRS_SUCR]>=0) # Generator's output schedule
	@variable(SUCRmodel, sucrm_genOut_Block[1:GENS, 1:BLOCKS, 1:HRS_SUCR]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)
	@variable(SUCRmodel, sucrm_genResUp[1:GENS, 1:HRS_SUCR]>=0) # Generators' up reserve schedule
	@variable(SUCRmodel, sucrm_genResNonSpin[1:GENS, 1:HRS_SUCR]>=0) # Generators' up reserve schedule
	@variable(SUCRmodel, sucrm_genResDn[1:GENS, 1:HRS_SUCR]>=0) # Generator's down rserve schedule
	@variable(SUCRmodel, sucrm_TotGenVioP[g=1:GENS, 1:HRS_SUCR]>=0)
	@variable(SUCRmodel, sucrm_TotGenVioN[g=1:GENS, 1:HRS_SUCR]>=0)
	@variable(SUCRmodel, sucrm_MaxGenVioP[g=1:GENS, 1:HRS_SUCR]>=0)
	@variable(SUCRmodel, sucrm_MinGenVioP[g=1:GENS, 1:HRS_SUCR]>=0)

	# Declaring the decision variables for peakers
	@variable(SUCRmodel, 0<=sucrm_peakerOnOff[1:PEAKERS, 0:HRS_SUCR]<=1) #Bin, relaxed
	@variable(SUCRmodel, 0<=sucrm_peakerStartUp[1:PEAKERS, 1:HRS_SUCR]<=1) # startup variable
	@variable(SUCRmodel, 0<=sucrm_peakerShutDown[1:PEAKERS, 1:HRS_SUCR]<=1) # shutdown variable
	@variable(SUCRmodel, sucrm_peakerOut[1:PEAKERS, 0:HRS_SUCR]>=0) # Generator's output schedule
	@variable(SUCRmodel, sucrm_peakerOut_Block[1:PEAKERS, 1:BLOCKS, 1:HRS_SUCR]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)
	@variable(SUCRmodel, sucrm_peakerResUp[1:PEAKERS, 1:HRS_SUCR]>=0) # Generators' up reserve schedule
	@variable(SUCRmodel, sucrm_peakerResNonSpin[1:PEAKERS, 1:HRS_SUCR]>=0) # Generators' up reserve schedule
	@variable(SUCRmodel, sucrm_peakerResDn[1:PEAKERS, 1:HRS_SUCR]>=0) # Generator's down rserve schedule

	# Declaring decision variables for storage Units
	@variable(SUCRmodel, sucrm_storgChrg[1:STORG_UNITS, 1:HRS_SUCR], Bin) #Bin variable equal to 1 if unit runs in the charging mode
	@variable(SUCRmodel, sucrm_storgDisc[1:STORG_UNITS, 1:HRS_SUCR], Bin) #Bin variable equal to 1 if unit runs in the discharging mode
	@variable(SUCRmodel, sucrm_storgIdle[1:STORG_UNITS, 1:HRS_SUCR], Bin) ##Bin variable equal to 1 if unit runs in the idle mode
	@variable(SUCRmodel, sucrm_storgChrgPwr[1:STORG_UNITS, 0:HRS_SUCR] >= 0) #Chargung power
	@variable(SUCRmodel, sucrm_storgDiscPwr[1:STORG_UNITS, 0:HRS_SUCR] >= 0) # Discharging Power
	@variable(SUCRmodel, sucrm_storgSOC[1:STORG_UNITS, 0:HRS_SUCR] >= 0) # state of charge (stored energy level for storage unit at time t)
	@variable(SUCRmodel, sucrm_storgResUp[1:STORG_UNITS, 0:HRS_SUCR] >= 0) # Scheduled up reserve
	@variable(SUCRmodel, sucrm_storgResDn[1:STORG_UNITS, 0:HRS_SUCR] >= 0) # Scheduled down reserve

	# declaring decision variables for renewable generation
	@variable(SUCRmodel, sucrm_solarG[1:N_ZONES, 1:HRS_SUCR] >= 0) # solar energy schedules
	@variable(SUCRmodel, sucrm_windG[1:N_ZONES, 1:HRS_SUCR] >= 0) # wind energy schedules
	@variable(SUCRmodel, sucrm_hydroG[1:N_ZONES, 1:HRS_SUCR] >= 0) # hydro energy schedules
	@variable(SUCRmodel, sucrm_solarGSpil[1:N_ZONES, 1:HRS_SUCR] >= 0) # solar energy schedules
	@variable(SUCRmodel, sucrm_windGSpil[1:N_ZONES, 1:HRS_SUCR] >= 0) # wind energy schedules
	@variable(SUCRmodel, sucrm_hydroGSpil[1:N_ZONES, 1:HRS_SUCR] >= 0) # hydro energy schedules

	# Declaring decision variables for hourly dispatched and curtailed demand
	@variable(SUCRmodel, sucrm_Demand[1:N_ZONES, 1:HRS_SUCR]>=0) # Hourly scheduled demand
	@variable(SUCRmodel, sucrm_Demand_Curt[1:N_ZONES, 1:HRS_SUCR]>=0) # Hourly schedule demand

	@variable(SUCRmodel, sucrm_OverGen[1:N_ZONES, 1:HRS_SUCR]>=0) # Hourly schedule demand

	# Declaring variables for transmission system
	@variable(SUCRmodel, sucrm_voltAngle[1:N_ZONES, 1:HRS_SUCR]) #voltage angle at zone/bus n in t//
	@variable(SUCRmodel, sucrm_powerFlow[1:N_ZONES, 1:M_ZONES, 1:HRS_SUCR]) #transmission Flow from zone n to zone m//

	# Defining the objective function that minimizes the overal cost of supplying electricity (=total variable, no-load, startup, and shutdown costs)
	@objective(SUCRmodel, Min, sum(sum(df_gens.IHRC_B1_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,1,t]
			+df_gens.IHRC_B2_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,2,t]
			+df_gens.IHRC_B3_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,3,t]
			+df_gens.IHRC_B4_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,4,t]
			+df_gens.IHRC_B5_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,5,t]
			+df_gens.IHRC_B6_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,6,t]
			+df_gens.IHRC_B7_HR[g]*fuelprice[g,day]*sucrm_genOut_Block[g,7,t]
			+df_gens.NoLoadHR[g]*fuelprice[g,day]*sucrm_genOnOff[g,t] +
			((df_gens.HotStartU_FixedCost[g]+(df_gens.HotStartU_HeatRate[g]*fuelprice[g,day]))*sucrm_genStartUp[g,t])
			+df_gens.ShutdownCost[g]*sucrm_genShutDown[g, t]
			+(sucrm_TotGenVioP[g,t]*VIOLATION_PENALTY)+(sucrm_TotGenVioN[g,t]*VIOLATION_PENALTY)+
			(sucrm_MaxGenVioP[g,t]*VIOLATION_PENALTY)+(sucrm_MinGenVioP[g,t]*VIOLATION_PENALTY) for g in 1:GENS)
			+sum(df_peakers.IHRC_B1_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,1,t]
			+df_peakers.IHRC_B2_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,2,t]
			+df_peakers.IHRC_B3_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,3,t]
			+df_peakers.IHRC_B4_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,4,t]
			+df_peakers.IHRC_B5_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,5,t]
			+df_peakers.IHRC_B6_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,6,t]
			+df_peakers.IHRC_B7_HR[k]*fuelprice_peakers[k,day]*sucrm_peakerOut_Block[k,7,t]
			+df_peakers.NoLoadHR[k]*fuelprice_peakers[k,day]*sucrm_peakerOnOff[k,t] +
			((df_peakers.HotStartU_FixedCost[k]+(df_peakers.HotStartU_HeatRate[k]*fuelprice_peakers[k,day]))*sucrm_peakerStartUp[k,t])
			+df_peakers.ShutdownCost[k]*sucrm_peakerShutDown[k, t] for k in 1:PEAKERS) for t in 1:HRS_SUCR)
		   +sum(sum((sucrm_Demand_Curt[n,t]*LOAD_SHED_PENALTY)+(sucrm_OverGen[n,t]*OVERGEN_PENALTY) for n=1:N_ZONES) for t=1:HRS_SUCR))

	#Initialization of commitment and dispatch variables for slow-start generators at t=0 (representing the last hour of previous scheduling horizon day=day-1 and t=24)
	@constraint(SUCRmodel, conInitGenOnOff[g=1:GENS], sucrm_genOnOff[g,0]==sucr_results.gens.onoff_init[g]) # initial generation level for generator g at t=0
	@constraint(SUCRmodel, conInitGenOut[g=1:GENS], sucrm_genOut[g,0]==sucr_results.gens.power_out_init[g]) # initial on/off status for generators g at t=0

	#Initialization of commitment and dispatch variables for peakers at t=0 (representing the last hour of previous scheduling horizon day=day-1 and t=24)
	@constraint(SUCRmodel, conInitGenOnOff_Peaker[k=1:PEAKERS], sucrm_peakerOnOff[k,0]==sucr_results.peakers.onoff_init[k]) # initial generation level for peaker k at t=0
	@constraint(SUCRmodel, conInitGenOut_Peaker[k=1:PEAKERS], sucrm_peakerOut[k,0]==sucr_results.peakers.power_out_init[k]) # initial on/off status for peaker k at t=0
	#Initialization of commitment and dispatch variables for storage units at t=0 (representing the last hour of previous scheduling horizon day=day-1 and t=24)
	@constraint(SUCRmodel, conInitSOC[p=1:STORG_UNITS], sucrm_storgSOC[p,0]==sucr_results.storg.soc_init[p]) # SOC for storage unit p at t=0

	# Baseload Operation of nuclear units
	@constraint(SUCRmodel, conNuckBaseLoad[t=1:HRS_SUCR, g=1:GENS], sucrm_genOnOff[g,t]>=df_gens.Nuclear[g]) #
	@constraint(SUCRmodel, conNuclearTotGenZone[t=1:HRS_SUCR, n=1:N_ZONES], sum((sucrm_genOut[g,t]*map_gens[g,n]*df_gens.Nuclear[g]) for g=1:GENS) -sucr_prepoc.nuclear_wa[t,n] ==0)

	#Limits on generation of cogen units
	@constraint(SUCRmodel, conCoGenBaseLoad[t=1:HRS_SUCR, g=1:GENS], sucrm_genOnOff[g,t]>=df_gens.Cogen[g]) #
	@constraint(SUCRmodel, conCoGenTotGenZone[t=1:HRS_SUCR, n=1:N_ZONES], sum((sucrm_genOut[g,t]*map_gens[g,n]*df_gens.Cogen[g]) for g=1:GENS) -sucr_prepoc.cogen_wa[t,n] ==0)

	# Constraints representing technical limits of conventional generators
	#Status transition trajectory of
	@constraint(SUCRmodel, conStartUpAndDn[t=1:HRS_SUCR, g=1:GENS], (sucrm_genOnOff[g,t] - sucrm_genOnOff[g,t-1] - sucrm_genStartUp[g,t] + sucrm_genShutDown[g,t])==0)
	# Max Power generation limit in Block 1
	@constraint(SUCRmodel, conMaxPowBlock1[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,1,t] <= df_gens.IHRC_B1_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 2
	@constraint(SUCRmodel, conMaxPowBlock2[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,2,t] <= df_gens.IHRC_B2_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 3
	@constraint(SUCRmodel, conMaxPowBlock3[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,3,t] <= df_gens.IHRC_B3_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 4
	@constraint(SUCRmodel, conMaxPowBlock4[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,4,t] <= df_gens.IHRC_B4_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 5
	@constraint(SUCRmodel, conMaxPowBlock5[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,5,t] <= df_gens.IHRC_B5_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 6
	@constraint(SUCRmodel, conMaxPowBlock6[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,6,t] <= df_gens.IHRC_B6_Q[g]*sucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 7
	@constraint(SUCRmodel, conMaxPowBlock7[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut_Block[g,7,t] <= df_gens.IHRC_B7_Q[g]*sucrm_genOnOff[g,t] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(SUCRmodel, conTotalGen[t=1:HRS_SUCR, g=1:GENS],  sum(sucrm_genOut_Block[g,b,t] for b=1:BLOCKS) + sucrm_TotGenVioP[g,t] - sucrm_TotGenVioN[g,t] >=sucrm_genOut[g,t])

	#Max power generation limit
	@constraint(SUCRmodel, conMaxPow[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut[g,t]+sucrm_genResUp[g,t] - sucrm_MaxGenVioP[g,t]  <= df_gens.MaxPowerOut[g]*sucrm_genOnOff[g,t] )
	# Min power generation limit
	@constraint(SUCRmodel, conMinPow[t=1:HRS_SUCR, g=1:GENS],  sucrm_genOut[g,t]-sucrm_genResDn[g,t] + sucrm_MinGenVioP[g,t] >= df_gens.MinPowerOut[g]*sucrm_genOnOff[g,t] )
	# Up reserve provision limit
	@constraint(SUCRmodel, conMaxResUp[t=1:HRS_SUCR, g=1:GENS], sucrm_genResUp[g,t] <= df_gens.SpinningRes_Limit[g]*sucrm_genOnOff[g,t] )
	# Non-Spinning Reserve Limit
	#@constraint(SUCRmodel, conMaxNonSpinResUp[t=1:HRS_SUCR, g=1:GENS], sucrm_genResNonSpin[g,t] <= (df_gens.NonSpinningRes_Limit[g]*(1-sucrm_genOnOff[g,t])) )
	@constraint(SUCRmodel, conMaxNonSpinResUp[t=1:HRS_SUCR, g=1:GENS], sucrm_genResNonSpin[g,t] <= 0 )
	#Down reserve provision limit
	@constraint(SUCRmodel, conMaxResDown[t=1:HRS_SUCR, g=1:GENS],  sucrm_genResDn[g,t] <= df_gens.SpinningRes_Limit[g]*sucrm_genOnOff[g,t] )
	#Up ramp rate limit
	@constraint(SUCRmodel, conRampRateUp[t=1:HRS_SUCR, g=1:GENS], (sucrm_genOut[g,t] - sucrm_genOut[g,t-1] <=(df_gens.RampUpLimit[g]*sucrm_genOnOff[g, t-1]) + (df_gens.RampStartUpLimit[g]*sucrm_genStartUp[g,t])))
	# Down ramp rate limit
	@constraint(SUCRmodel, conRampRateDown[t=1:HRS_SUCR, g=1:GENS], (sucrm_genOut[g,t-1] - sucrm_genOut[g,t] <=(df_gens.RampDownLimit[g]*sucrm_genOnOff[g,t]) + (df_gens.RampShutDownLimit[g]*sucrm_genShutDown[g,t])))
	# Min Up Time limit with alternative formulation
	@constraint(SUCRmodel, conUpTime[t=1:HRS_SUCR, g=1:GENS], (sum(sucrm_genStartUp[g,k] for k=lb_MUT[g,t]:t)<=sucrm_genOnOff[g,t]))
	# Min down Time limit with alternative formulation
	@constraint(SUCRmodel, conDownTime[t=1:HRS_SUCR, g=1:GENS], (1-sum(sucrm_genShutDown[g,i] for i=lb_MDT[g,t]:t)>=sucrm_genOnOff[g,t]))

	# Constraints representing technical limits of Peakers
	#Status transition trajectory of peakers
	@constraint(SUCRmodel, conStartUpAndDn_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], (sucrm_peakerOnOff[k,t] - sucrm_peakerOnOff[k,t-1] - sucrm_peakerStartUp[k,t] + sucrm_peakerShutDown[k,t])==0)
	# Max Power generation limit in Block 1
	@constraint(SUCRmodel, conMaxPowBlock1_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,1,t] <= df_peakers.IHRC_B1_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 2
	@constraint(SUCRmodel, conMaxPowBlock2_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,2,t] <= df_peakers.IHRC_B2_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 3
	@constraint(SUCRmodel, conMaxPowBlock3_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,3,t] <= df_peakers.IHRC_B3_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 4
	@constraint(SUCRmodel, conMaxPowBlock4_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,4,t] <= df_peakers.IHRC_B4_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 5
	@constraint(SUCRmodel, conMaxPowBlock5_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,5,t] <= df_peakers.IHRC_B5_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 6
	@constraint(SUCRmodel, conMaxPowBlock6_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,6,t] <= df_peakers.IHRC_B6_Q[k]*sucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 7
	@constraint(SUCRmodel, conMaxPowBlock7_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut_Block[k,7,t] <= df_peakers.IHRC_B7_Q[k]*sucrm_peakerOnOff[k,t] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(SUCRmodel, conTotalGen_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sum(sucrm_peakerOut_Block[k,b,t] for b=1:BLOCKS)>=sucrm_peakerOut[k,t])
	#Max power generation limit
	@constraint(SUCRmodel, conMaxPow_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut[k,t]+sucrm_peakerResUp[k,t] <= df_peakers.MaxPowerOut[k]*sucrm_peakerOnOff[k,t] )
	# Min power generation limit
	@constraint(SUCRmodel, conMinPow_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerOut[k,t]-sucrm_peakerResDn[k,t] >= df_peakers.MinPowerOut[k]*sucrm_peakerOnOff[k,t] )
	# Up reserve provision limit
	@constraint(SUCRmodel, conMaxResUp_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], sucrm_peakerResUp[k,t] <= df_peakers.SpinningRes_Limit[k]*sucrm_peakerOnOff[k,t] )
	# Non-Spinning Reserve Limit
	@constraint(SUCRmodel, conMaxNonSpinResUp_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], sucrm_peakerResNonSpin[k,t] <= (df_peakers.NonSpinningRes_Limit[k]*(1-sucrm_peakerOnOff[k,t])))
	#Down reserve provision limit
	@constraint(SUCRmodel, conMaxResDown_Peaker[t=1:HRS_SUCR, k=1:PEAKERS],  sucrm_peakerResDn[k,t] <= df_peakers.SpinningRes_Limit[k]*sucrm_peakerOnOff[k,t] )
	#Up ramp rate limit
	@constraint(SUCRmodel, conRampRateUp_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], (sucrm_peakerOut[k,t] - sucrm_peakerOut[k,t-1] <=(df_peakers.RampUpLimit[k]*sucrm_peakerOnOff[k, t-1]) + (df_peakers.RampStartUpLimit[k]*sucrm_peakerStartUp[k,t])))
	# Down ramp rate limit
	@constraint(SUCRmodel, conRampRateDown_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], (sucrm_peakerOut[k,t-1] - sucrm_peakerOut[k,t] <=(df_peakers.RampDownLimit[k]*sucrm_peakerOnOff[k,t]) + (df_peakers.RampShutDownLimit[k]*sucrm_peakerShutDown[k,t])))
	# Min Up Time limit with alternative formulation
	@constraint(SUCRmodel, conUpTime_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], (sum(sucrm_peakerStartUp[k,r] for r=lb_MUT_peaker[k,t]:t)<=sucrm_peakerOnOff[k,t]))
	# Min down Time limit with alternative formulation
	@constraint(SUCRmodel, conDownTime_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], (1-sum(sucrm_peakerShutDown[k,s] for s=lb_MDT_peaker[k,t]:t)>=sucrm_peakerOnOff[k,t]))

	# Renewable generation constraints
	@constraint(SUCRmodel, conSolarLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_solarG[n, t] + sucrm_solarGSpil[n,t]<=sucr_prepoc.solar_wa[t,n])
	@constraint(SUCRmodel, conWindLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_windG[n, t] + sucrm_windGSpil[n,t]<=sucr_prepoc.wind_wa[t,n])
	@constraint(SUCRmodel, conHydroLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_hydroG[n, t] + sucrm_hydroGSpil[n,t]<=sucr_prepoc.hydro_wa[t,n])

	# Constraints representing technical characteristics of storage units
	# status transition of storage units between charging, discharging, and idle modes
	@constraint(SUCRmodel, conStorgStatusTransition[t=1:HRS_SUCR, p=1:STORG_UNITS], (sucrm_storgChrg[p,t]+sucrm_storgDisc[p,t]+sucrm_storgIdle[p,t])==1)
	# charging power limit
	@constraint(SUCRmodel, conStrgChargPowerLimit[t=1:HRS_SUCR, p=1:STORG_UNITS], (sucrm_storgChrgPwr[p,t] - sucrm_storgResDn[p,t])<=df_storage.Power[p]*sucrm_storgChrg[p,t])
	# Discharging power limit
	@constraint(SUCRmodel, conStrgDisChgPowerLimit[t=1:HRS_SUCR, p=1:STORG_UNITS], (sucrm_storgDiscPwr[p,t] + sucrm_storgResUp[p,t])<=df_storage.Power[p]*sucrm_storgDisc[p,t])
	# Down reserve provision limit
	@constraint(SUCRmodel, conStrgDownResrvMax[t=1:HRS_SUCR, p=1:STORG_UNITS], sucrm_storgResDn[p,t]<=df_storage.Power[p]*sucrm_storgChrg[p,t])
	# Up reserve provision limit`
	@constraint(SUCRmodel, conStrgUpResrvMax[t=1:HRS_SUCR, p=1:STORG_UNITS], sucrm_storgResUp[p,t]<=df_storage.Power[p]*sucrm_storgDisc[p,t])
	# State of charge at t
	@constraint(SUCRmodel, conStorgSOC[t=1:HRS_SUCR, p=1:STORG_UNITS], sucrm_storgSOC[p,t]==sucrm_storgSOC[p,t-1]-(sucrm_storgDiscPwr[p,t]/df_storage.TripEfficDown[p])+(sucrm_storgChrgPwr[p,t]*df_storage.TripEfficUp[p])-(sucrm_storgSOC[p,t]*df_storage.SelfDischarge[p]))
	# minimum energy limit
	@constraint(SUCRmodel, conMinEnrgStorgLimi[t=1:HRS_SUCR, p=1:STORG_UNITS], sucrm_storgSOC[p,t]-(sucrm_storgResUp[p,t]/df_storage.TripEfficDown[p])+(sucrm_storgResDn[p,t]/df_storage.TripEfficUp[p])>=0)
	# Maximum energy limit
	@constraint(SUCRmodel, conMaxEnrgStorgLimi[t=1:HRS_SUCR, p=1:STORG_UNITS], sucrm_storgSOC[p,t]-(sucrm_storgResUp[p,t]/df_storage.TripEfficDown[p])+(sucrm_storgResDn[p,t]/df_storage.TripEfficUp[p])<=(df_storage.Power[p]/df_storage.PowerToEnergRatio[p]))
	# Constraints representing transmission grid capacity constraints
	# DC Power Flow Calculation
	#@constraint(SUCRmodel, conDCPowerFlowPos[t=1:HRS_SUCR, n=1:N_ZONES, m=1:N_ZONES], sucrm_powerFlow[n,m,t]-(TranS[n,m]*(sucrm_voltAngle[n,t]-sucrm_voltAngle[m,t])) ==0)
	@constraint(SUCRmodel, conDCPowerFlowNeg[t=1:HRS_SUCR, n=1:N_ZONES, m=1:N_ZONES], sucrm_powerFlow[n,m,t]+sucrm_powerFlow[m,n,t]==0)
	# Tranmission flow bounds (from n to m and from m to n)
	@constraint(SUCRmodel, conPosFlowLimit[t=1:HRS_SUCR, n=1:N_ZONES, m=1:N_ZONES], sucrm_powerFlow[n,m,t]<=TranC[n,m])
	@constraint(SUCRmodel, conNegFlowLimit[t=1:HRS_SUCR, n=1:N_ZONES, m=1:N_ZONES], sucrm_powerFlow[n,m,t]>=-TranC[n,m])
	# Voltage Angle bounds and reference point
	#@constraint(SUCRmodel, conVoltAnglUB[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_voltAngle[n,t]<=π)
	#@constraint(SUCRmodel, conVoltAnglLB[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_voltAngle[n,t]>=-π)
	#@constraint(SUCRmodel, conVoltAngRef[t=1:HRS_SUCR], sucrm_voltAngle[1,t]==0)

	# Demand-side Constraints
	@constraint(SUCRmodel, conDemandLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_Demand[n,t]+ sucrm_Demand_Curt[n,t] == sucr_prepoc.wk_ahead[t,n])

	# Demand Curtailment and wind generation limits
	@constraint(SUCRmodel, conDemandCurtLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_Demand_Curt[n,t] <= LOAD_SHED_MAX);
	@constraint(SUCRmodel, conOverGenLimit[t=1:HRS_SUCR, n=1:N_ZONES], sucrm_OverGen[n,t] <= OVERGEN_MAX);

	# System-wide Constraints
	#nodal balance constraint
	@constraint(SUCRmodel, conNodBalanc[t=1:HRS_SUCR, n=1:N_ZONES], sum((sucrm_genOut[g,t]*map_gens[g,n]) for g=1:GENS) + sum((sucrm_storgDiscPwr[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) - sum((sucrm_storgChrgPwr[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) +sucrm_solarG[n, t] +sucrm_windG[n, t] +sucrm_hydroG[n, t] - sucrm_OverGen[n,t] - sucrm_Demand[n,t] == sum(sucrm_powerFlow[n,m,t] for m=1:M_ZONES))

	# Minimum zonal up reserve requirement, if there are more than two zones, we should  define reserve regions for DEC and DEP
	#@constraint(SUCRmodel, conMinUpReserveReq[t=1:HRS_SUCR, n=1:N_ZONES], sum((sucrm_genResUp[g,t]*map_gens[g,n]) for g=1:GENS) + sum((sucrm_storgResUp[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) >= Reserve_Req_Up[n] )
	@constraint(SUCRmodel, conMinUpReserveReq[t=1:HRS_SUCR], sum((sucrm_genResUp[g,t]+sucrm_genResNonSpin[g,t]) for g=1:GENS) + sum((sucrm_storgResUp[p,t]) for p=1:STORG_UNITS) >= sum(Reserve_Req_Up[n] for n=1:N_ZONES))

	# Minimum down reserve requirement
	#    @constraint(SUCRmodel, conMinDnReserveReq[t=1:HRS_SUCR], sum(genResDn[g,t] for g=1:GENS) + sum(storgResDn[p,t] for p=1:STORG_UNITS) >= Reserve_Req_Dn[t] )

	t2_SUCRmodel = time_ns()
	time_SUCRmodel = (t2_SUCRmodel -t1_SUCRmodel)/1.0e9;
	@info "SUCRmodel for day: $day setup executed in (s): $time_SUCRmodel";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("SUCRmodel", time_SUCRmodel, "day: $day",
					"", "", "Model Setup"), ',')
	end; # close file

	# solve the Second WA-UC model (SUCR)
	JuMP.optimize!(SUCRmodel)

	open(".//outputs//objective_values_v76.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("SUCRmodel", "day: $day",
					"", "", JuMP.objective_value(SUCRmodel)), ',')
	end;

	# Pricing general results in the terminal window
	println("Objective value: ", JuMP.objective_value(SUCRmodel))
	println("------------------------------------")
	println("------- SURC OBJECTIVE VALUE -------")
	println("Objective value for day ", day, " is ", JuMP.objective_value(SUCRmodel))
	println("------------------------------------")
	println("------- SURC PRIMAL STATUS -------")
	println(primal_status(SUCRmodel))
	println("------------------------------------")
	println("------- SURC DUAL STATUS -------")
	println(JuMP.dual_status(SUCRmodel))
	println("Day: ", day, " solved")
	println("---------------------------")
	println("SUCRmodel Number of variables: ", JuMP.num_variables(SUCRmodel))
	@info "SUCRmodel Number of variables: " JuMP.num_variables(SUCRmodel)

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("SUCRmodel", JuMP.num_variables(SUCRmodel), "day: $day",
					"", "", "Variables"), ',')
	end;

	@debug "SUCRmodel for day: $day optimized executed in (s):  $(solve_time(SUCRmodel))";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("SUCRmodel", solve_time(SUCRmodel), "day: $day",
					"", "", "Model Optimization"), ',')
	end; # close file

	## Write the optimal outcomes into spreadsheets
	#TODO: Later we need to include a variable for day so the cell
	# number in which the results are printed is updated accordingly

	# Write the conventional generators' schedules
	t1_write_SUCRmodel_results = time_ns()
	open(".//outputs//SUCR_GenOutputs.csv", FILE_ACCESS_APPEND) do io
		for t in 1:HRS_SUCR, g=1:GENS
			writedlm(io, hcat(day, t+INIT_HR_SUCR, g, df_gens.UNIT_ID[g],
				df_gens.MinPowerOut[g], df_gens.MaxPowerOut[g],
				JuMP.value.(sucrm_genOut[g,t]), JuMP.value.(sucrm_genOnOff[g,t]),
				JuMP.value.(sucrm_genShutDown[g,t]), JuMP.value.(sucrm_genStartUp[g,t]),
				JuMP.value.(sucrm_genResUp[g,t]), JuMP.value.(sucrm_genResNonSpin[g,t]),
				JuMP.value.(sucrm_genResDn[g,t]), JuMP.value.(sucrm_TotGenVioP[g,t]),
				JuMP.value.(sucrm_TotGenVioN[g,t]), JuMP.value.(sucrm_MaxGenVioP[g,t]),
				JuMP.value.(sucrm_MinGenVioP[g,t])), ',')
		end # ends the loop
	end; # close file

	open(".//outputs//SUCR_PeakerOutputs.csv", FILE_ACCESS_APPEND) do io
		for t in 1:HRS_SUCR, k=1:PEAKERS
			writedlm(io, hcat(day, t+INIT_HR_SUCR, k, df_peakers.UNIT_ID[k],
				df_peakers.MinPowerOut[k], df_peakers.MaxPowerOut[k],
				JuMP.value.(sucrm_peakerOut[k,t]), JuMP.value.(sucrm_peakerOnOff[k,t]),
				JuMP.value.(sucrm_peakerShutDown[k,t]), JuMP.value.(sucrm_peakerStartUp[k,t]),
				JuMP.value.(sucrm_peakerResUp[k,t]), JuMP.value.(sucrm_peakerResNonSpin[k,t]),
				JuMP.value.(sucrm_peakerResDn[k,t])), ',')
		end # ends the loop
	end; # close file

	# Writing storage units' optimal schedules in CSV file
	open(".//outputs//SUCR_StorageOutputs.csv", FILE_ACCESS_APPEND) do io
		for t in 1:HRS_SUCR, p=1:STORG_UNITS
			writedlm(io, hcat(day, t+INIT_HR_SUCR, p, df_storage.Power[p],
				df_storage.Power[p]/df_storage.PowerToEnergRatio[p],
				JuMP.value.(sucrm_storgChrg[p,t]), JuMP.value.(sucrm_storgDisc[p,t]),
				JuMP.value.(sucrm_storgIdle[p,t]), JuMP.value.(sucrm_storgChrgPwr[p,t]),
				JuMP.value.(sucrm_storgDiscPwr[p,t]), JuMP.value.(sucrm_storgSOC[p,t]),
				JuMP.value.(sucrm_storgResUp[p,t]), JuMP.value.(sucrm_storgResDn[p,t]) ), ',')
		end # ends the loop
	end; # close file

	# Writing the transmission flow schedules into spreadsheets
	open(".//outputs//SUCR_TranFlowOutputs.csv", FILE_ACCESS_APPEND) do io
		for t in 1:HRS_SUCR, n=1:N_ZONES, m=1:M_ZONES
			writedlm(io, hcat(day, t+INIT_HR_SUCR, n, m,
				JuMP.value.(sucrm_powerFlow[n,m,t]), TranC[n,m] ), ',')
		end # ends the loop
	end; # close file

	# Writing the curtilment, overgeneration, and spillage outcomes in CSV file
	open(".//outputs//SUCR_Curtail.csv", FILE_ACCESS_APPEND) do io
		for t in 1:HRS_SUCR, n=1:N_ZONES
		   writedlm(io, hcat(day, t+INIT_HR_SUCR, n,
			   JuMP.value.(sucrm_OverGen[n,t]), JuMP.value.(sucrm_Demand_Curt[n,t]),
			   JuMP.value.(sucrm_windGSpil[n,t]), JuMP.value.(sucrm_solarGSpil[n,t]),
			   JuMP.value.(sucrm_hydroGSpil[n,t])), ',')
		end # ends the loop
	end; # close file

	t2_write_SUCRmodel_results = time_ns()
	time_write_SUCRmodel_results = (t2_write_SUCRmodel_results -t1_write_SUCRmodel_results)/1.0e9;
	@info "Write SUCRmodel results for day $day: $time_write_SUCRmodel_results executed in (s)";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("SUCRmodel", time_write_SUCRmodel_results, "day: $day",
					"", "", "Write CSV files"), ',')
	end; #close file

	# Create and save the following parameters, which are transefered to BAUC2
	for h=1:24-INIT_HR_SUCR+INIT_HR_FUCR
		for g=1:GENS
			uc_results.gens.onoff[g,h]= round(JuMP.value.(sucrm_genOnOff[g,h]));
			uc_results.gens.power_out[g,h]= JuMP.value.(sucrm_genOut[g,h]);
			uc_results.gens.startup[g,h]= round(JuMP.value.(sucrm_genStartUp[g,h]));
			uc_results.gens.shutdown[g,h]=round(JuMP.value.(sucrm_genShutDown[g,h]));
			for b=1:BLOCKS
				uc_results.gens.genout_block[g,b,h]=JuMP.value.(sucrm_genOut_Block[g,b,h]);
			end
		end
		for k=1:PEAKERS
			uc_results.peakers.onoff[k,h]= round(JuMP.value.(sucrm_peakerOnOff[k,h]));
			uc_results.peakers.power_out[k,h]=JuMP.value.(sucrm_peakerOut[k,h]);
			uc_results.peakers.startup[k,h]= round(JuMP.value.(sucrm_peakerStartUp[k,h]));
			uc_results.peakers.shutdown[k,h]= round(JuMP.value.(sucrm_peakerShutDown[k,h]));
			for b=1:BLOCKS
				uc_results.peakers.genout_block[k,b,h]=JuMP.value.(sucrm_peakerOut_Block[k,b,h]);
			end
		end
		for p=1:STORG_UNITS
			uc_results.storg.chrg[p,h]=round(JuMP.value.(sucrm_storgChrg[p,h]));
			uc_results.storg.disc[p,h]=round(JuMP.value.(sucrm_storgDisc[p,h]));
			uc_results.storg.idle[p,h]=round(JuMP.value.(sucrm_storgIdle[p,h]));
			uc_results.storg.chrgpwr[p,h]=JuMP.value.(sucrm_storgChrgPwr[p,h]);
			uc_results.storg.discpwr[p,h]=JuMP.value.(sucrm_storgDiscPwr[p,h]);
			uc_results.storg.soc[p,h]=JuMP.value.(sucrm_storgSOC[p,h]);
		end
	end

	return uc_results;

end; #end function

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
@testset "ALWB" begin
    using RasterDataSources: rastername, rasterurl, Values, Deciles

    alwb_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "ALWB")
    @test rasterpath(ALWB) == alwb_path
    @test rasterpath(ALWB{Values,Year}) == joinpath(alwb_path, "values", "year")
    @test rastername(ALWB{Values,Year}, :ss_pct; date=Date(2001, 1)) == "ss_pct.nc"
    @test rastername(ALWB{Values,Month}, :ss_pct; date=Date(2001, 1)) == "ss_pct.nc"
    @test rastername(ALWB{Values,Day}, :ss_pct; date=Date(2001, 1)) == "ss_pct_2001.nc"
    @test rasterpath(ALWB{Values,Day}, :ss_pct; date=Date(2001, 1)) ==
        joinpath(alwb_path, "values", "day", "ss_pct_2001.nc")

    @test rasterurl(ALWB{Values,Year}, :ss_pct; date=Date(2001, 1)) ==
        URI(scheme="http", host="www.bom.gov.au", path="/jsp/awra/thredds/fileServer/AWRACMS/values/year/ss_pct.nc")
    @test rasterurl(ALWB{Values,Day}, :ss_pct; date=Date(2001, 1)) ==
        URI(scheme="http", host="www.bom.gov.au", path="/jsp/awra/thredds/fileServer/AWRACMS/values/day/ss_pct_2001.nc")

    date = DateTime(2018, 01, 01), DateTime(2018, 03, 02)
    getraster(ALWB{Values,Day}, :ss_pct; date)
    @test isfile(joinpath(alwb_path, "values", "day", "ss_pct_2018.nc"))
    getraster(ALWB{Deciles,Month}, :s0_pct; date)
    @test isfile(joinpath(alwb_path, "deciles", "month", "s0_pct.nc"))

    getraster(ALWB{Values,Day}, :ma_wet; date)
    @test isfile(joinpath(alwb_path, "values", "day", "ma_wet_2018.nc"))
    getraster(ALWB{Values,Month}, :etot; date)
    @test isfile(joinpath(alwb_path, "values", "month", "etot.nc"))
    getraster(ALWB{Values,Year}, :asce_pet; date)
    @test isfile(joinpath(alwb_path, "values", "year", "asce_pet.nc"))
    getraster(ALWB{Values,Year}, :dd; date)
    @test isfile(joinpath(alwb_path, "values", "year", "dd.nc"))
end