text:
export PiecewiseSegment

struct PiecewiseSegment{T} <: Domain{T}
    points::Vector{T}
    PiecewiseSegment{T}(d::Vector{T}) where {T} = new{T}(d)
end
PiecewiseSegment(d::AbstractVector) = PiecewiseSegment{eltype(d)}(collect(d))
PiecewiseSegment(d...) = PiecewiseSegment(collect(mapreduce(eltype,promote_type,d),d))

function PiecewiseSegment(pcsin::AbstractVector{IT}) where IT<:IntervalOrSegment
    pcs=collect(pcsin)
    p=∂(pop!(pcs))
    successful=true
    while successful
        successful=false
        for k=1:length(pcs)
            if leftendpoint(pcs[k]) == last(p)
                push!(p,rightendpoint(pcs[k]))
                deleteat!(pcs,k)
                successful=true
                break
            end
        end
    end
    @assert isempty(pcs)
    PiecewiseSegment(p)
end

==(a::PiecewiseSegment,b::PiecewiseSegment) = a.points==b.points

indomain(x, d::PiecewiseSegment) = any(in.(x, components(d)))

canonicaldomain(d::PiecewiseSegment)=d
ncomponents(d::PiecewiseSegment)=length(d.points)-1
component(d::PiecewiseSegment,j::Integer) = Segment(d.points[j],d.points[j+1])
components(d::PiecewiseSegment{T}) where {T} = Segment{T}[component(d,k) for k=1:ncomponents(d)]

for OP in (:arclength,:complexlength)
    @eval $OP(d::PiecewiseSegment) = mapreduce($OP,+,components(d))
end

isperiodic(_) = false
isperiodic(d::PiecewiseSegment) = first(d.points) == last(d.points)

reverseorientation(d::PiecewiseSegment) = PiecewiseSegment(reverse(d.points))

isambiguous(d::PiecewiseSegment) = isempty(d.points)
convert(::Type{PiecewiseSegment{T}},::AnyDomain) where {T<:Number} = PiecewiseSegment{T}([])
convert(::Type{IT},::AnyDomain) where {IT<:PiecewiseSegment}=PiecewiseSegment(Float64[])

function points(d::PiecewiseSegment,n)
   k=div(n,ncomponents(d))
    r=n-ncomponents(d)*k

    float(eltype(d))[vcat([points(component(d,j),k+1) for j=1:r]...);
        vcat([points(component(d,j),k) for j=r+1:ncomponents(d)]...)]
end

rand(d::PiecewiseSegment) = rand(d[rand(1:ncomponents(d))])
checkpoints(d::PiecewiseSegment{T}) where {T} = mapreduce(checkpoints,union,components(d))

leftendpoint(d::PiecewiseSegment) = first(d.points)
rightendpoint(d::PiecewiseSegment) = last(d.points)

# Comparison with UnionDomain
for OP in (:(isapprox),:(==))
    @eval begin
        $OP(a::PiecewiseSegment,b::UnionDomain) = $OP(UnionDomain(components(a)),b)
        $OP(b::UnionDomain,a::PiecewiseSegment) = $OP(UnionDomain(components(a)),b)
    end
end

function union(S::PiecewiseSegment{<:Real}, D::IntervalOrSegment{<:Real})
    isempty(D) && return S
    a,b = endpoints(D)
    (a ∈ S || b ∈ S) && return PiecewiseSegment(sort!(union(S.points, a, b)))
    UnionDomain(S, D)
end
union(D::IntervalOrSegment{<:Real}, S::PiecewiseSegment{<:Real}) = union(S,D)

for OP in (:minimum, :maximum)
    @eval $OP(d::PiecewiseSegment) = $OP(d.points)
end


text:
## IntSet

# Construction, collect
data_in = (1,5,100)
s = IntSet(data_in)
data_out = collect(s)
@test all(map(d->in(d,data_out), data_in))
@test length(data_out) == length(data_in)

# eltype, similar
@test eltype(IntSet()) ===  Int
@test isequal(similar(IntSet([1,2,3])), IntSet())

# show
@test sprint(show, IntSet()) == "IntSet([])"
@test sprint(show, IntSet([1,2,3])) == "IntSet([1, 2, 3])"
@test contains(sprint(show, complement(IntSet())), "...,")

s = IntSet([0,1,10,20,200,300,1000,10000,10002])
@test last(s) == 10002
@test first(s) == 0
@test length(s) == 9
@test pop!(s) == 10002
@test length(s) == 8
@test shift!(s) == 0
@test length(s) == 7
@test !in(0,s)
@test !in(10002,s)
@test in(10000,s)
@test_throws ArgumentError first(IntSet())
@test_throws ArgumentError last(IntSet())
t = copy(s)
sizehint!(t, 20000) #check that hash does not depend on size of internal Array{UInt32, 1}
@test hash(s) == hash(t)
@test hash(complement(s)) == hash(complement(t))

@test setdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3])
@test symdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) == IntSet([1, 3, 5, 6])

s2 = IntSet([1, 2, 3, 4])
setdiff!(s2, IntSet([2, 4, 5, 6]))

@test s2 == IntSet([1, 3])

# issue #7851
@test_throws ArgumentError IntSet(-1)
@test !(-1 in IntSet(0:10))

# # issue #8570
# This requires 2^29 bytes of storage, which is too much for a simple test
# s = IntSet(2^32)
# @test length(s) == 1
# for b in s; b; end

# Copy, copy!, similar
s1 = IntSet([1,2,3])
s2 = similar(s1)
copy!(s2, s1)
s3 = copy(s2)
@test s3 == s2 == s1
@test collect(s3) == collect(s2) == [1,2,3]

c1 = complement!(IntSet())
pop!(c1, 1)
pop!(c1, 2)
pop!(c1, 3)
c2 = similar(c1)
copy!(c2, c1)
c3 = copy(c2)
c4 = complement(s1)
@test c1 == c2 == c3 == c4
@test c4 === sizehint!(c4, 100)
@test c1 == c4
@test last(c1) == typemax(Int)-1
@test last(complement(IntSet())) == typemax(Int)-1

# Push, union
s1 = IntSet()
@test_throws ArgumentError push!(s1, -1)
push!(s1, 1, 10, 100, 1000)
@test collect(s1) == [1, 10, 100, 1000]
push!(s1, 606)
@test collect(s1) == [1, 10, 100, 606, 1000]
s2 = IntSet()
@test s2 === union!(s2, s1)
s3 = IntSet([1, 10, 100])
union!(s3, [1, 606, 1000])
s4 = union(IntSet([1, 100, 1000]), IntSet([10, 100, 606]))
@test s1 == s2 == s3 == s4

c1 = complement(s1)
@test !(1 in c1)
push!(c1, 1)
@test 1 in c1
push!(c1, 10, 100, 10)
@test collect(complement(c1)) == [606, 1000]
c2 = complement(IntSet([606, 1000]))
@test c2 === union!(c2, c1)
c3 = union!(complement(IntSet([10, 606, 1000])), complement(IntSet([1, 606, 1000, 2000])))
@test c3 == union(complement(IntSet([1, 606, 1000, 2000])), complement(IntSet([10, 606, 1000])))
c4 = union!(complement(IntSet([1, 10, 606, 1000, 1001])), complement(IntSet([606, 1000])))
@test c4 == union(complement(IntSet([606, 1000])), complement(IntSet([1, 10, 606, 1000, 1001])))
c5 = union!(complement(IntSet([10, 606, 1000])), IntSet([0, 10, 20, 30]))
@test c5 == union(IntSet([0, 10, 20, 30]), complement(IntSet([10, 606, 1000])))
c6 = union!(complement(IntSet([10, 606, 1000])), IntSet([10, 4000, 3]))
@test c6 == union(IntSet([10, 4000, 3]), complement(IntSet([10, 606, 1000])))
@test c1 == c2 == c3 == c4 == c5 == c6

## pop!, delete!
s = IntSet(1:2:10)
@test pop!(s, 1) == 1
@test !(1 in s)
@test_throws KeyError pop!(s, 1)
@test_throws ArgumentError pop!(s, -1)
@test_throws ArgumentError pop!(s, -1, 1)
@test_throws ArgumentError pop!(()->throw(ErrorException()), s, -1)
@test pop!(s, 1, 0) == 0
@test s === delete!(s, 1)
for i in s; pop!(s, i); end
@test isempty(s)
x = 0
@test 1 == pop!(()->(global x; x+=1), s, 100)
@test x == 1
push!(s, 100)
@test pop!(()->throw(ErrorException()), s, 100) == 100
push!(s, 1:2:10...)
@test pop!(s) == 9
@test pop!(s) == 7
@test shift!(s) == 1
@test shift!(s) == 3
@test collect(s) == [5]
empty!(s)
@test isempty(s)

c = complement(IntSet())
@test pop!(c, 1) == 1
@test !(1 in c)
@test_throws KeyError pop!(c, 1)
@test_throws ArgumentError pop!(c, -1)
@test_throws ArgumentError pop!(c, -1, 1)
@test_throws ArgumentError pop!(()->throw(ErrorException()), c, -1)
@test pop!(c, 1, 0) == 0
@test c === delete!(c, 1)
@test shift!(c) == 0
@test shift!(c) == 2
@test_throws ArgumentError pop!(c)
@test collect(complement(c)) == [0,1,2]
@test empty!(c) == IntSet()

## Intersect
@test isempty(intersect(IntSet()))
@test isempty(intersect(IntSet(1:10), IntSet()))
@test isempty(intersect(IntSet(), IntSet(1:10)))
@test isempty(intersect(IntSet(), complement(IntSet())))
@test isempty(intersect(IntSet(), complement(IntSet(1:10))))
@test isempty(intersect(complement(IntSet()), IntSet()))
@test isempty(intersect(complement(IntSet(1:10)), IntSet()))

@test intersect(IntSet([1,2,3])) == IntSet([1,2,3])
@test intersect(complement!(IntSet()), IntSet(1)) ==
      intersect(IntSet(1), complement!(IntSet())) == IntSet(1)

@test intersect(IntSet(0:7), IntSet(3:10)) ==
      intersect(IntSet(3:10), IntSet(0:7)) == IntSet(3:7)
@test intersect(complement(IntSet([0:2; 11:16])), IntSet(0:7)) ==
      intersect(IntSet(0:7), complement(IntSet([0:2; 11:16]))) == IntSet(3:7)

@test intersect(complement(IntSet(5:12)), complement(IntSet(7:10))) ==
      intersect(complement(IntSet(7:10)), complement(IntSet(5:12))) == complement(IntSet(5:12))

@test intersect(IntSet(0:10), IntSet(1:4), 0:5, [1,2,10]) == IntSet(1:2)

## Setdiff
s1 = IntSet(1:100)
setdiff!(s1, IntSet(1:2:100))
s2 = setdiff(IntSet(1:100), IntSet(1:2:100))
@test s1 == s2 == IntSet(2:2:100)
@test collect(s1) == collect(2:2:100)

s1 = IntSet(1:10)
s2 = complement(IntSet(3:5))
@test setdiff(s1, s2) == setdiff(s1, [0:2; 6:100]) == IntSet(3:5)
@test isempty(setdiff(complement(IntSet()), complement(IntSet())))
@test setdiff(complement(IntSet()), complement(IntSet(3:5))) == IntSet(3:5)
@test setdiff(complement(IntSet(1:5)), complement(IntSet(3:10))) == IntSet([6, 7, 8, 9, 10])
@test setdiff(complement(IntSet(2:2:10)), IntSet(1:5)) == complement(IntSet([1:5; 6:2:10]))
@test setdiff!(complement(IntSet()), complement(IntSet())) == IntSet()
@test setdiff!(complement(IntSet(0:2:10)), complement(IntSet(0:10))) == IntSet(1:2:9)

## Symdiff
@test symdiff(IntSet([1, 2, 3, 4]), IntSet([2, 4, 5, 6])) ==
      symdiff(IntSet([2, 4, 5, 6]), IntSet([1, 2, 3, 4])) ==
      symdiff(IntSet([1, 2, 3, 4]), [2, 4, 5, 6]) ==
      symdiff(IntSet([2, 4, 5, 6]), [1, 2, 3, 4]) == IntSet([1, 3, 5, 6])
@test symdiff(complement(IntSet()), IntSet(2:3)) ==
      symdiff(IntSet(2:3), complement(IntSet())) ==
      symdiff(complement(IntSet()), 2:3) ==complement(IntSet(2:3))

@test symdiff(complement(IntSet(2:7)), IntSet(5:10)) ==
      symdiff(IntSet(5:10), complement(IntSet(2:7))) == complement(IntSet([2:4; 8:10]))
@test symdiff(complement(IntSet(3:7)), complement(IntSet(5:10))) == IntSet([3:4; 8:10])

## Subsets, equality
@test IntSet(2:2:10) < IntSet(1:10)
@test !(IntSet(2:2:10) < IntSet(2:2:10))
@test IntSet(2:2:10) <= IntSet(2:10)
@test IntSet(2:2:10) <= IntSet(2:2:10)
@test IntSet(1) < complement!(IntSet())
@test IntSet(1) <= complement!(IntSet())
@test !(IntSet(1) < complement!(IntSet(1)))

# Test logic against Set
p = IntSet([0,1,4,5])
resize!(p.bits, 6)
q = IntSet([0,2,4,6])
resize!(q.bits, 8)
p′ = complement(p)
q′ = complement(q)
function collect10(itr)
    r = eltype(itr)[]
    for i in itr
        i > 10 && break
        push!(r, i)
    end
    r
end
a = Set(p)
b = Set(q)
a′ = Set(collect10(p′))
b′ = Set(collect10(q′))
for f in (union, intersect, setdiff, symdiff)
    @test collect(f(p, p)) == sort(collect(f(a, a)))
    @test collect(f(q, q)) == sort(collect(f(b, b)))
    @test collect(f(p, q)) == sort(collect(f(a, b)))
    @test collect(f(q, p)) == sort(collect(f(b, a)))

    @test collect10(f(p′, p)) == sort(collect(f(a′, a)))
    @test collect10(f(q′, q)) == sort(collect(f(b′, b)))
    @test collect10(f(p′, q)) == sort(collect(f(a′, b)))
    @test collect10(f(q′, p)) == sort(collect(f(b′, a)))

    @test collect10(f(p, p′)) == sort(collect(f(a, a′)))
    @test collect10(f(q, q′)) == sort(collect(f(b, b′)))
    @test collect10(f(p, q′)) == sort(collect(f(a, b′)))
    @test collect10(f(q, p′)) == sort(collect(f(b, a′)))

    @test collect10(f(p′, p′)) == sort(collect(f(a′, a′)))
    @test collect10(f(q′, q′)) == sort(collect(f(b′, b′)))
    @test collect10(f(p′, q′)) == sort(collect(f(a′, b′)))
    @test collect10(f(q′, p′)) == sort(collect(f(b′, a′)))
end

## Other
s = IntSet()
push!(s, 0, 2, 100)
@test 0 in s
@test !(1 in s)
@test 2 in s
@test 100 in s
@test !(101 in s)
@test !(1000 in s)
@test first(s) == 0
@test last(s) == 100
@test s == IntSet([0, 2, 100])
push!(s, 1000)
@test [i for i in s] == [0, 2, 100, 1000]
@test pop!(s) == 1000
@test s == IntSet([0, 2, 100])
@test hash(s) == hash(IntSet([0, 2, 100]))

b = 0:1000
s = IntSet(b)
@test collect(s) == collect(b)
@test length(s) == length(b)
@test pop!(s, 100) == 100
@test collect(s) == [0:99; 101:1000]
@test_throws KeyError pop!(s, 100)
@test pop!(s, 100, 0) == 0
@test pop!(s, 99, 0) == 99
@test pop!(()->1, s, 99) == 1
@test pop!(()->1, s, 98) == 98

show(IOBuffer(), IntSet())
show(IOBuffer(), complement(IntSet()))

ext: jl


text:
"""
    methoddataset(data, method, context)

Transform data container `data` of samples into a data container of `(x, y)`-pairs.
Maps `encode(method, context, sample)` over the observations in `data`.
"""
@with_kw struct MethodDataset{M<:LearningMethod}
    data
    method::M
    context::Context
end

LearnBase.nobs(ds::MethodDataset) = nobs(ds.data)

function LearnBase.getobs(ds::MethodDataset, idx)
    return encode(ds.method, ds.context, getobs(ds.data, idx))
end

function LearnBase.getobs!(buf, ds::MethodDataset, idx)
    return encode!(buf, ds.method, ds.context, getobs(ds.data, idx))
end

"""
    methoddataset(data, method, context)

Transform data container `data` of samples into a data container of `(x, y)`-pairs.
Maps `encode(method, context, sample)` over the observations in `data`.
"""
const methoddataset = MethodDataset

"""
    methoddataloaders(data, method)
    methoddataloaders(traindata, validdata, method[, batchsize; shuffle = true, dlkwargs...])

Create training and validation `DataLoader`s from two data containers `(traindata, valdata)`.
If only one container `data` is passed, splits it into two with `pctgvalid`% of the data
going into the validation split.

## Keyword arguments

- `batchsize = 16`
- `shuffle = true`: Whether to shuffle the training data container
- `validbsfactor`: Factor to multiply batchsize for validation data loader with (validation
    batches can be larger since no GPU memory is needed for the backward pass)

All remaining keyword arguments are passed to [`DataLoader`](#).
"""
function methoddataloaders(
        traindata,
        validdata,
        method::LearningMethod,
        batchsize = 16;
        shuffle = true,
        validbsfactor = 2,
        kwargs...)
    traindata = shuffle ? shuffleobs(traindata) : traindata
    return (
        DataLoader(methoddataset(traindata, method, Training()), batchsize; kwargs...),
        DataLoader(methoddataset(validdata, method, Validation()), validbsfactor * batchsize; kwargs...),
    )
end

function methoddataloaders(
        data,
        method::LearningMethod,
        batchsize = 16;
        pctgval = 0.2,
        shuffle = true,
        kwargs...)
    traindata, validdata = splitobs(shuffleobs(data), at = 1-pctgval)
    methoddataloaders(traindata, validdata, method, batchsize; shuffle = false, kwargs...)
end


text:
function fucr_model(day::Int64, df_gens::DataFrame, df_peakers::DataFrame,
	df_storage::DataFrame, fuelprice::Matrix{Float64},
	fuelprice_peakers::Matrix{Float64},  map_gens::Matrix{Int64},
	map_peakers::Matrix{Int64}, map_storage::Matrix{Int64},
	fucr_results::UC_Results, fucr_prepoc::DemandWAPreprocGens,
	lb_MDT::Matrix{Int64}, lb_MUT::Matrix{Int64},
	lb_MDT_peaker::Matrix{Int64}, lb_MUT_peaker::Matrix{Int64})

	t1_FUCRmodel = time_ns()

	uc_gens = GensResults();
	uc_peakers = PeakersResults();
	uc_storg = StorageResults();
	uc_results = UC_Results(uc_gens, uc_peakers, uc_storg);

	FUCRmodel = direct_model(CPLEX.Optimizer())
	# Enable Benders strategy
	#MOI.set(FUCRmodel, MOI.RawParameter("CPXPARAM_Benders_Strategy"), 3)
	set_optimizer_attribute(FUCRmodel, "CPX_PARAM_EPGAP", SOLVER_EPGAP)

  	# Declaring the decision variables for conventional generators
	@variable(FUCRmodel, fucrm_genOnOff[1:GENS, 0:HRS_FUCR], Bin) #Bin
	@variable(FUCRmodel, fucrm_genStartUp[1:GENS, 1:HRS_FUCR], Bin) # startup variable
	@variable(FUCRmodel, fucrm_genShutDown[1:GENS, 1:HRS_FUCR], Bin) # shutdown variable
	@variable(FUCRmodel, fucrm_genOut[1:GENS, 0:HRS_FUCR]>=0) # Generator's output schedule
	@variable(FUCRmodel, fucrm_genOut_Block[1:GENS, 1:BLOCKS, 1:HRS_FUCR]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)
	@variable(FUCRmodel, fucrm_genResUp[1:GENS, 1:HRS_FUCR]>=0) # Generators' up reserve schedule
	@variable(FUCRmodel, fucrm_genResNonSpin[1:GENS, 1:HRS_FUCR]>=0) # Scheduled up reserve on offline fast-start peakers
	@variable(FUCRmodel, fucrm_genResDn[1:GENS, 1:HRS_FUCR]>=0) # Generator's down rserve schedule
	@variable(FUCRmodel, fucrm_TotGenVioP[g=1:GENS, 1:HRS_FUCR]>=0)
	@variable(FUCRmodel, fucrm_TotGenVioN[g=1:GENS, 1:HRS_FUCR]>=0)
	@variable(FUCRmodel, fucrm_MaxGenVioP[g=1:GENS, 1:HRS_FUCR]>=0)
	@variable(FUCRmodel, fucrm_MinGenVioP[g=1:GENS, 1:HRS_FUCR]>=0)

  	# Declaring the decision variables for peaker units
	@variable(FUCRmodel, 0<=fucrm_peakerOnOff[1:PEAKERS, 0:HRS_FUCR]<=1)
	@variable(FUCRmodel, 0<=fucrm_peakerStartUp[1:PEAKERS, 1:HRS_FUCR]<=1)
	@variable(FUCRmodel, 0<=fucrm_peakerShutDown[1:PEAKERS, 1:HRS_FUCR]<=1)
	@variable(FUCRmodel, fucrm_peakerOut[1:PEAKERS, 0:HRS_FUCR]>=0) # Generator's output schedule
	@variable(FUCRmodel, fucrm_peakerOut_Block[1:PEAKERS, 1:BLOCKS, 1:HRS_FUCR]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)
	@variable(FUCRmodel, fucrm_peakerResUp[1:PEAKERS, 1:HRS_FUCR]>=0) # Generators' up reserve schedule
	@variable(FUCRmodel, fucrm_peakerResNonSpin[1:PEAKERS, 1:HRS_FUCR]>=0) # Scheduled up reserve on offline fast-start peakers
	@variable(FUCRmodel, fucrm_peakerResDn[1:PEAKERS, 1:HRS_FUCR]>=0) # Generator's down rserve schedule

    # declaring decision variables for storage Units
	@variable(FUCRmodel, fucrm_storgChrg[1:STORG_UNITS, 1:HRS_FUCR], Bin) #Bin variable equal to 1 if unit runs in the charging mode
	@variable(FUCRmodel, fucrm_storgDisc[1:STORG_UNITS, 1:HRS_FUCR], Bin) #Bin variable equal to 1 if unit runs in the discharging mode
	@variable(FUCRmodel, fucrm_storgIdle[1:STORG_UNITS, 1:HRS_FUCR], Bin) ##Bin variable equal to 1 if unit runs in the idle mode
	@variable(FUCRmodel, fucrm_storgChrgPwr[1:STORG_UNITS, 0:HRS_FUCR]>=0) #Chargung power
	@variable(FUCRmodel, fucrm_storgDiscPwr[1:STORG_UNITS, 0:HRS_FUCR]>=0) # Discharging Power
	@variable(FUCRmodel, fucrm_storgSOC[1:STORG_UNITS, 0:HRS_FUCR]>=0) # state of charge (stored energy level for storage unit at time t)
	@variable(FUCRmodel, fucrm_storgResUp[1:STORG_UNITS, 0:HRS_FUCR]>=0) # Scheduled up reserve
	@variable(FUCRmodel, fucrm_storgResDn[1:STORG_UNITS, 0:HRS_FUCR]>=0) # Scheduled down reserve

	# declaring decision variables for renewable generation
	@variable(FUCRmodel, fucrm_solarG[1:N_ZONES, 1:HRS_FUCR]>=0) # solar energy schedules
	@variable(FUCRmodel, fucrm_windG[1:N_ZONES, 1:HRS_FUCR]>=0) # wind energy schedules
	@variable(FUCRmodel, fucrm_hydroG[1:N_ZONES, 1:HRS_FUCR]>=0) # hydro energy schedules
	@variable(FUCRmodel, fucrm_solarGSpil[1:N_ZONES, 1:HRS_FUCR]>=0) # solar energy schedules
	@variable(FUCRmodel, fucrm_windGSpil[1:N_ZONES, 1:HRS_FUCR]>=0) # wind energy schedules
	@variable(FUCRmodel, fucrm_hydroGSpil[1:N_ZONES, 1:HRS_FUCR]>=0) # hydro energy schedules

	# Declaring decision variables for hourly dispatched and curtailed demand
	@variable(FUCRmodel, fucrm_Demand[1:N_ZONES, 1:HRS_FUCR]>=0) # Hourly scheduled demand
	@variable(FUCRmodel, fucrm_Demand_Curt[1:N_ZONES, 1:HRS_FUCR]>=0) # Hourly schedule demand

	# Declaring variables for transmission system
	@variable(FUCRmodel, fucrm_voltAngle[1:N_ZONES, 1:HRS_FUCR]) #voltage angle at zone/bus n in t//
	@variable(FUCRmodel, fucrm_powerFlow[1:N_ZONES, 1:M_ZONES, 1:HRS_FUCR]) #transmission Flow from zone n to zone m//

	# Declaring over and undergeneration decision variable
	@variable(FUCRmodel, fucrm_OverGen[1:N_ZONES, 1:HRS_FUCR]>=0) #overgeneration at zone n and time t//

    # Defining the objective function that minimizes the overal cost of supplying electricity (=total variable, no-load, startup, and shutdown costs)
	@objective(FUCRmodel, Min, sum(sum(df_gens.IHRC_B1_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,1,t]
		 +df_gens.IHRC_B2_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,2,t]
	     +df_gens.IHRC_B3_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,3,t]
         +df_gens.IHRC_B4_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,4,t]
         +df_gens.IHRC_B5_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,5,t]
         +df_gens.IHRC_B6_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,6,t]
         +df_gens.IHRC_B7_HR[g]*fuelprice[g,day]*fucrm_genOut_Block[g,7,t]
         +df_gens.NoLoadHR[g]*fuelprice[g,day]*fucrm_genOnOff[g,t]
         +((df_gens.HotStartU_FixedCost[g]
         +(df_gens.HotStartU_HeatRate[g]*fuelprice[g,day]))*fucrm_genStartUp[g,t])
         +df_gens.ShutdownCost[g]*fucrm_genShutDown[g, t]
         +(fucrm_TotGenVioP[g,t]*VIOLATION_PENALTY)
         +(fucrm_TotGenVioN[g,t]*VIOLATION_PENALTY)
         +(fucrm_MaxGenVioP[g,t]*VIOLATION_PENALTY)
         +(fucrm_MinGenVioP[g,t]*VIOLATION_PENALTY) for g in 1:GENS)
         +sum(df_peakers.IHRC_B1_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,1,t]
         +df_peakers.IHRC_B2_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,2,t]
         +df_peakers.IHRC_B3_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,3,t]
         +df_peakers.IHRC_B4_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,4,t]
         +df_peakers.IHRC_B5_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,5,t]
         +df_peakers.IHRC_B6_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,6,t]
         +df_peakers.IHRC_B7_HR[k]*fuelprice_peakers[k,day]*fucrm_peakerOut_Block[k,7,t]
         +df_peakers.NoLoadHR[k]*fuelprice_peakers[k,day]*fucrm_peakerOnOff[k,t]
         +((df_peakers.HotStartU_FixedCost[k]
         +(df_peakers.HotStartU_HeatRate[k]*fuelprice_peakers[k,day]))*fucrm_peakerStartUp[k,t])
         +df_peakers.ShutdownCost[k]*fucrm_peakerShutDown[k, t] for k in 1:PEAKERS) for t in 1:HRS_FUCR)
         +sum(sum((fucrm_Demand_Curt[n,t]*LOAD_SHED_PENALTY)
         +(fucrm_OverGen[n,t]*OVERGEN_PENALTY) for n=1:N_ZONES) for t=1:HRS_FUCR))

	#Initialization of commitment and dispatch variables for convnentioal generatoes at t=0 (representing the last hour of previous scheduling horizon day=day-1 and t=24)
	@constraint(FUCRmodel, conInitGenOnOff[g=1:GENS], fucrm_genOnOff[g,0]==fucr_results.gens.onoff_init[g]) # initial generation level for generator g at t=0
	@constraint(FUCRmodel, conInitGenOut[g=1:GENS], fucrm_genOut[g,0]==fucr_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(FUCRmodel, conInitGenOnOff_Peakers[k=1:PEAKERS], fucrm_peakerOnOff[k,0]==fucr_results.peakers.onoff_init[k]) # initial generation level for peaker k at t=0
	@constraint(FUCRmodel, conInitGenOut_Peakers[k=1:PEAKERS], fucrm_peakerOut[k,0]==fucr_results.peakers.power_out_init[k]) # initial on/off status for peaker k at t=0
	#Initialization of SOC variables for storage units at t=0 (representing the last hour of previous scheduling horizon day=day-1 and t=24)
	@constraint(FUCRmodel, conInitSOC[p=1:STORG_UNITS], fucrm_storgSOC[p,0]==fucr_results.storg.soc_init[p]) # SOC for storage unit p at t=0

	#Base-Load Operation of nuclear Generators
	@constraint(FUCRmodel, conNuckBaseLoad[t=1:HRS_FUCR, g=1:GENS], fucrm_genOnOff[g,t]>=df_gens.Nuclear[g]) #
	@constraint(FUCRmodel, conNuclearTotGenZone[t=1:HRS_FUCR, n=1:N_ZONES], sum((fucrm_genOut[g,t]*map_gens[g,n]*df_gens.Nuclear[g]) for g=1:GENS) -fucr_prepoc.nuclear_wa[t,n] ==0)

	#Limits on generation of cogen units
	@constraint(FUCRmodel, conCoGenBaseLoad[t=1:HRS_FUCR, g=1:GENS], fucrm_genOnOff[g,t]>=df_gens.Cogen[g]) #
	@constraint(FUCRmodel, conCoGenTotGenZone[t=1:HRS_FUCR, n=1:N_ZONES], sum((fucrm_genOut[g,t]*map_gens[g,n]*df_gens.Cogen[g]) for g=1:GENS) -fucr_prepoc.cogen_wa[t,n] ==0)

	# Constraints representing technical limits of conventional generators
	#Status transition trajectory of
	@constraint(FUCRmodel, conStartUpAndDn[t=1:HRS_FUCR, g=1:GENS], (fucrm_genOnOff[g,t] - fucrm_genOnOff[g,t-1] - fucrm_genStartUp[g,t] + fucrm_genShutDown[g,t])==0)
	# Max Power generation limit in Block 1
	@constraint(FUCRmodel, conMaxPowBlock1[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,1,t] <= df_gens.IHRC_B1_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 2
	@constraint(FUCRmodel, conMaxPowBlock2[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,2,t] <= df_gens.IHRC_B2_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 3
	@constraint(FUCRmodel, conMaxPowBlock3[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,3,t] <= df_gens.IHRC_B3_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 4
	@constraint(FUCRmodel, conMaxPowBlock4[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,4,t] <= df_gens.IHRC_B4_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 5
	@constraint(FUCRmodel, conMaxPowBlock5[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,5,t] <= df_gens.IHRC_B5_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 6
	@constraint(FUCRmodel, conMaxPowBlock6[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,6,t] <= df_gens.IHRC_B6_Q[g]*fucrm_genOnOff[g,t] )
	# Max Power generation limit in Block 7
	@constraint(FUCRmodel, conMaxPowBlock7[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut_Block[g,7,t] <= df_gens.IHRC_B7_Q[g]*fucrm_genOnOff[g,t] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(FUCRmodel, conTotalGen[t=1:HRS_FUCR, g=1:GENS],  sum(fucrm_genOut_Block[g,b,t] for b=1:BLOCKS) + fucrm_TotGenVioP[g,t] - fucrm_TotGenVioN[g,t] ==fucrm_genOut[g,t])
	#Max power generation limit
	@constraint(FUCRmodel, conMaxPow[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut[g,t]+fucrm_genResUp[g,t] - fucrm_MaxGenVioP[g,t] <= df_gens.MaxPowerOut[g]*fucrm_genOnOff[g,t] )
	# Min power generation limit
	@constraint(FUCRmodel, conMinPow[t=1:HRS_FUCR, g=1:GENS],  fucrm_genOut[g,t]-fucrm_genResDn[g,t] + fucrm_MinGenVioP[g,t] >= df_gens.MinPowerOut[g]*fucrm_genOnOff[g,t] )
	# Up reserve provision limit
	@constraint(FUCRmodel, conMaxResUp[t=1:HRS_FUCR, g=1:GENS], fucrm_genResUp[g,t] <= df_gens.SpinningRes_Limit[g]*fucrm_genOnOff[g,t] )
	# Non-Spinning Reserve Limit
	#    @constraint(FUCRmodel, conMaxNonSpinResUp[t=1:HRS_SUCR, g=1:GENS], fucrm_genResNonSpin[g,t] <= (df_gens.NonSpinningRes_Limit[g]*(1-fucrm_genOnOff[g,t])*df_gens.FastStart[g]))
	@constraint(FUCRmodel, conMaxNonSpinResUp[t=1:HRS_SUCR, g=1:GENS], fucrm_genResNonSpin[g,t] <= 0)
	#Down reserve provision limit
	@constraint(FUCRmodel, conMaxResDown[t=1:HRS_FUCR, g=1:GENS],  fucrm_genResDn[g,t] <= df_gens.SpinningRes_Limit[g]*fucrm_genOnOff[g,t] )
	#Up ramp rate limit
	@constraint(FUCRmodel, conRampRateUp[t=1:HRS_FUCR, g=1:GENS], (fucrm_genOut[g,t] - fucrm_genOut[g,t-1] <=(df_gens.RampUpLimit[g]*fucrm_genOnOff[g, t-1]) + (df_gens.RampStartUpLimit[g]*fucrm_genStartUp[g,t])))
	# Down ramp rate limit
	@constraint(FUCRmodel, conRampRateDown[t=1:HRS_FUCR, g=1:GENS], (fucrm_genOut[g,t-1] - fucrm_genOut[g,t] <=(df_gens.RampDownLimit[g]*fucrm_genOnOff[g,t]) + (df_gens.RampShutDownLimit[g]*fucrm_genShutDown[g,t])))
	# Min Up Time limit with alternative formulation
	@constraint(FUCRmodel, conUpTime[t=1:HRS_FUCR, g=1:GENS], (sum(fucrm_genStartUp[g,r] for r=lb_MUT[g,t]:t)<=fucrm_genOnOff[g,t]))
	# Min down Time limit with alternative formulation
	@constraint(FUCRmodel, conDownTime[t=1:HRS_FUCR, g=1:GENS], (1-sum(fucrm_genShutDown[g,s] for s=lb_MDT[g,t]:t)>=fucrm_genOnOff[g,t]))

	# Peaker Units' constraints
	#Status transition trajectory of
	@constraint(FUCRmodel, conStartUpAndDn_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], (fucrm_peakerOnOff[k,t] - fucrm_peakerOnOff[k,t-1] - fucrm_peakerStartUp[k,t] + fucrm_peakerShutDown[k,t])==0)
	# Max Power generation limit in Block 1
	@constraint(FUCRmodel, conMaxPowBlock1_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,1,t] <= df_peakers.IHRC_B1_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 2
	@constraint(FUCRmodel, conMaxPowBlock2_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,2,t] <= df_peakers.IHRC_B2_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 3
	@constraint(FUCRmodel, conMaxPowBlock3_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,3,t] <= df_peakers.IHRC_B3_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 4
	@constraint(FUCRmodel, conMaxPowBlock4_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,4,t] <= df_peakers.IHRC_B4_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 5
	@constraint(FUCRmodel, conMaxPowBlock5_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,5,t] <= df_peakers.IHRC_B5_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 6
	@constraint(FUCRmodel, conMaxPowBlock6_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,6,t] <= df_peakers.IHRC_B6_Q[k]*fucrm_peakerOnOff[k,t] )
	# Max Power generation limit in Block 7
	@constraint(FUCRmodel, conMaxPowBlock7_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut_Block[k,7,t] <= df_peakers.IHRC_B7_Q[k]*fucrm_peakerOnOff[k,t] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(FUCRmodel, conTotalGen_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  sum(fucrm_peakerOut_Block[k,b,t] for b=1:BLOCKS)>=fucrm_peakerOut[k,t])
	#Max power generation limit
	@constraint(FUCRmodel, conMaxPow_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut[k,t]+fucrm_peakerResUp[k,t] <= df_peakers.MaxPowerOut[k]*fucrm_peakerOnOff[k,t] )
	# Min power generation limit
	@constraint(FUCRmodel, conMinPow_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerOut[k,t]-fucrm_peakerResDn[k,t] >= df_peakers.MinPowerOut[k]*fucrm_peakerOnOff[k,t] )
	# Up reserve provision limit
	@constraint(FUCRmodel, conMaxResUp_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], fucrm_peakerResUp[k,t] <= df_peakers.SpinningRes_Limit[k]*fucrm_peakerOnOff[k,t] )
	# Non-Spinning Reserve Limit
	@constraint(FUCRmodel, conMaxNonSpinResUp_Peaker[t=1:HRS_SUCR, k=1:PEAKERS], fucrm_peakerResNonSpin[k,t] <= (df_peakers.NonSpinningRes_Limit[k]*(1-fucrm_peakerOnOff[k,t])))
	#Down reserve provision limit
	@constraint(FUCRmodel, conMaxResDown_Peaker[t=1:HRS_FUCR, k=1:PEAKERS],  fucrm_peakerResDn[k,t] <= df_peakers.SpinningRes_Limit[k]*fucrm_peakerOnOff[k,t] )
	#Up ramp rate limit
	@constraint(FUCRmodel, conRampRateUp_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], (fucrm_peakerOut[k,t] - fucrm_peakerOut[k,t-1] <=(df_peakers.RampUpLimit[k]*fucrm_peakerOnOff[k, t-1]) + (df_peakers.RampStartUpLimit[k]*fucrm_peakerStartUp[k,t])))
	# Down ramp rate limit
	@constraint(FUCRmodel, conRampRateDown_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], (fucrm_peakerOut[k,t-1] - fucrm_peakerOut[k,t] <=(df_peakers.RampDownLimit[k]*fucrm_peakerOnOff[k,t]) + (df_peakers.RampShutDownLimit[k]*fucrm_peakerShutDown[k,t])))
	# Min Up Time limit with alternative formulation
	@constraint(FUCRmodel, conUpTime_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], (sum(fucrm_peakerStartUp[k,r] for r=lb_MUT_peaker[k,t]:t)<=fucrm_peakerOnOff[k,t]))
	# Min down Time limit with alternative formulation
	@constraint(FUCRmodel, conDownTime_Peaker[t=1:HRS_FUCR, k=1:PEAKERS], (1-sum(fucrm_peakerShutDown[k,s] for s=lb_MDT_peaker[k,t]:t)>=fucrm_peakerOnOff[k,t]))

	# Renewable generation constraints
	@constraint(FUCRmodel, conSolarLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_solarG[n, t] + fucrm_solarGSpil[n,t]<=fucr_prepoc.solar_wa[t,n])
	@constraint(FUCRmodel, conWindLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_windG[n, t] + fucrm_windGSpil[n,t]<=fucr_prepoc.wind_wa[t,n])
	@constraint(FUCRmodel, conHydroLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_hydroG[n, t] + fucrm_hydroGSpil[n,t]<=fucr_prepoc.hydro_wa[t,n])

	#=
	@constraint(FUCRmodel, conSolarLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_solarG[n, t] + fucrm_solarGSpil[n,t]<=0)
	@constraint(FUCRmodel, conWindLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_windG[n, t] + fucrm_windGSpil[n,t]<=0)
	@constraint(FUCRmodel, conHydroLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_hydroG[n, t] + fucrm_hydroGSpil[n,t]<=0)
	=#

	# Constraints representing technical characteristics of storage units
	# status transition of storage units between charging, discharging, and idle modes
	@constraint(FUCRmodel, conStorgStatusTransition[t=1:HRS_FUCR, p=1:STORG_UNITS], (fucrm_storgChrg[p,t]+fucrm_storgDisc[p,t]+fucrm_storgIdle[p,t])==1)
	# charging power limit
	@constraint(FUCRmodel, conStrgChargPowerLimit[t=1:HRS_FUCR, p=1:STORG_UNITS], (fucrm_storgChrgPwr[p,t] - fucrm_storgResDn[p,t])<=df_storage.Power[p]*fucrm_storgChrg[p,t])
	# Discharging power limit
	@constraint(FUCRmodel, conStrgDisChgPowerLimit[t=1:HRS_FUCR, p=1:STORG_UNITS], (fucrm_storgDiscPwr[p,t] + fucrm_storgResUp[p,t])<=df_storage.Power[p]*fucrm_storgDisc[p,t])
	# Down reserve provision limit
	@constraint(FUCRmodel, conStrgDownResrvMax[t=1:HRS_FUCR, p=1:STORG_UNITS], fucrm_storgResDn[p,t]<=df_storage.Power[p]*fucrm_storgChrg[p,t])
	# Up reserve provision limit`
	@constraint(FUCRmodel, conStrgUpResrvMax[t=1:HRS_FUCR, p=1:STORG_UNITS], fucrm_storgResUp[p,t]<=df_storage.Power[p]*fucrm_storgDisc[p,t])
	# State of charge at t
	@constraint(FUCRmodel, conStorgSOC[t=1:HRS_FUCR, p=1:STORG_UNITS], fucrm_storgSOC[p,t]==fucrm_storgSOC[p,t-1]-(fucrm_storgDiscPwr[p,t]/df_storage.TripEfficDown[p])+(fucrm_storgChrgPwr[p,t]*df_storage.TripEfficUp[p])-(fucrm_storgSOC[p,t]*df_storage.SelfDischarge[p]))
	# minimum energy limit
	@constraint(FUCRmodel, conMinEnrgStorgLimi[t=1:HRS_FUCR, p=1:STORG_UNITS], fucrm_storgSOC[p,t]-(fucrm_storgResUp[p,t]/df_storage.TripEfficDown[p])+(fucrm_storgResDn[p,t]/df_storage.TripEfficUp[p])>=0)
	# Maximum energy limit
	@constraint(FUCRmodel, conMaxEnrgStorgLimi[t=1:HRS_FUCR, p=1:STORG_UNITS], fucrm_storgSOC[p,t]-(fucrm_storgResUp[p,t]/df_storage.TripEfficDown[p])+(fucrm_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(FUCRmodel, conDCPowerFlowPos[t=1:HRS_FUCR, n=1:N_ZONES, m=1:N_ZONES], fucrm_powerFlow[n,m,t]-(TranS[n,m]*(fucrm_voltAngle[n,t]-fucrm_voltAngle[m,t])) ==0)
	@constraint(FUCRmodel, conDCPowerFlowNeg[t=1:HRS_FUCR, n=1:N_ZONES, m=1:N_ZONES], fucrm_powerFlow[n,m,t]+fucrm_powerFlow[m,n,t]==0)
	# Tranmission flow bounds (from n to m and from m to n)
	@constraint(FUCRmodel, conPosFlowLimit[t=1:HRS_FUCR, n=1:N_ZONES, m=1:N_ZONES], fucrm_powerFlow[n,m,t]<=TranC[n,m])
	@constraint(FUCRmodel, conNegFlowLimit[t=1:HRS_FUCR, n=1:N_ZONES, m=1:N_ZONES], fucrm_powerFlow[n,m,t]>=-TranC[n,m])
	# Voltage Angle bounds and reference point
	#@constraint(FUCRmodel, conVoltAnglUB[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_voltAngle[n,t]<=π)
	#@constraint(FUCRmodel, conVoltAnglLB[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_voltAngle[n,t]>=-π)
	#@constraint(FUCRmodel, conVoltAngRef[t=1:HRS_FUCR], fucrm_voltAngle[1,t]==0)

	# Demand-side Constraints
	@constraint(FUCRmodel, conDemandLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_Demand[n,t]+ fucrm_Demand_Curt[n,t] == fucr_prepoc.wk_ahead[t,n])

	# Demand Curtailment and wind generation limits
	@constraint(FUCRmodel, conDemandCurtLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_Demand_Curt[n,t] <= LOAD_SHED_MAX);
	@constraint(FUCRmodel, conOverGenLimit[t=1:HRS_FUCR, n=1:N_ZONES], fucrm_OverGen[n,t] <= OVERGEN_MAX);

	# System-wide Constraints
	#nodal balance constraint
	@constraint(FUCRmodel, conNodBalanc[t=1:HRS_FUCR, n=1:N_ZONES], sum((fucrm_genOut[g,t]*map_gens[g,n]) for g=1:GENS) +sum((fucrm_peakerOut[k,t]*map_peakers[k,n]) for k=1:PEAKERS)  + sum((fucrm_storgDiscPwr[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) - sum((fucrm_storgChrgPwr[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) +fucrm_solarG[n, t] +fucrm_windG[n, t] +fucrm_hydroG[n, t] - fucrm_Demand[n,t] - fucrm_OverGen[n,t]== sum(fucrm_powerFlow[n,m,t] for m=1:M_ZONES))

	#@constraint(FUCRmodel, conNodBalanc[t=1:HRS_FUCR], sum(fucrm_genOut[g,t] for g=1:GENS) + sum((fucrm_storgDiscPwr[p,t]) for p=1:STORG_UNITS) - sum((fucrm_storgChrgPwr[p,t]) for p=1:STORG_UNITS) +sum(fucrm_solarG[n, t] for n=1:N_ZONES) + sum(fucrm_windG[n, t] for n=1:N_ZONES)+ sum(fucrm_hydroG[n, t] for n=1:N_ZONES) - sum(fucr_prepoc.wk_ahead[t,n] for n=1:N_ZONES) == 0)

	# @constraint(FUCRmodel, conNodBalanc[t=1:HRS_FUCR], sum((fucrm_genOut[g,t]) for g=1:GENS) + sum((fucrm_storgDiscPwr[p,t]) for p=1:STORG_UNITS) - sum((fucrm_storgChrgPwr[p,t]) for p=1:STORG_UNITS) +sum((fucrm_solarG[n, t]) for n=1:N_ZONES) +sum((fucrm_windG[n, t]) for n=1:N_ZONES) +sum((fucrm_hydroG[n, t]) for n=1:N_ZONES) - sum((fucr_prepoc.wk_ahead[t,n]) for n=1:N_ZONES) == 0)
	# Minimum zonal up reserve requirement, if there are more than two zones, we should  define reserve regions for DEC and DEP
	#@constraint(FUCRmodel, conMinUpReserveReq[t=1:HRS_FUCR, n=1:N_ZONES], sum((fucrm_genResUp[g,t]*map_gens[g,n]) for g=1:GENS) + sum((fucrm_storgResUp[p,t]*map_storage[p,n]) for p=1:STORG_UNITS) >= Reserve_Req_Up[n] )
	#@constraint(FUCRmodel, conMinUpReserveReq[t=1:HRS_FUCR], sum((fucrm_genResUp[g,t]+fucrm_genResNonSpin[g,t]) for g=1:GENS) + sum((fucrm_storgResUp[p,t]) for p=1:STORG_UNITS) >= sum(Reserve_Req_Up[n] for n=1:N_ZONES))
	@constraint(FUCRmodel, conMinUpReserveReq[t=1:HRS_FUCR], sum((fucrm_genResUp[g,t]) for g=1:GENS) + sum((fucrm_peakerResUp[k,t]+fucrm_peakerResNonSpin[k,t]) for k=1:PEAKERS)+ sum((fucrm_storgResUp[p,t]) for p=1:STORG_UNITS) >= sum(Reserve_Req_Up[n] for n=1:N_ZONES))

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

	t2_FUCRmodel = time_ns()
	time_FUCRmodel = (t2_FUCRmodel -t1_FUCRmodel)/1.0e9;
	@info "FUCRmodel for day: $day setup executed in (s): $time_FUCRmodel";

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

	# solve the First WAUC model (FUCR)
	JuMP.optimize!(FUCRmodel)

	# Pricing general results in the terminal window
	println("Objective value: ", JuMP.objective_value(FUCRmodel))

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

	println("------------------------------------")
	println("------- FUCR OBJECTIVE VALUE -------")
	println("Objective value for day ", day, " is ", JuMP.objective_value(FUCRmodel))
	println("------------------------------------")
	println("-------FUCR PRIMAL STATUS -------")
	println(primal_status(FUCRmodel))
	println("------------------------------------")
	println("------- FUCR DUAL STATUS -------")
	println(JuMP.dual_status(FUCRmodel))
	println("Day: ", day, " solved")
	println("---------------------------")
	println("FUCRmodel Number of variables: ", JuMP.num_variables(FUCRmodel))
	@info "FUCRmodel Number of variables: " JuMP.num_variables(FUCRmodel)

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

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

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

	# Write the conventional generators' schedules in CSV file
	t1_write_FUCRmodel_results = time_ns()
	open(".//outputs//FUCR_GenOutputs.csv", FILE_ACCESS_APPEND) do io
	  for t in 1:HRS_FUCR, g=1:GENS
	      writedlm(io, hcat(day, t+INIT_HR_FUCR, g, df_gens.UNIT_NAME[g],
	          df_gens.MinPowerOut[g], df_gens.MaxPowerOut[g],
	          JuMP.value.(fucrm_genOut[g,t]), JuMP.value.(fucrm_genOnOff[g,t]),
	          JuMP.value.(fucrm_genShutDown[g,t]), JuMP.value.(fucrm_genStartUp[g,t]),
	          JuMP.value.(fucrm_genResUp[g,t]), JuMP.value.(fucrm_genResNonSpin[g,t]),
	          JuMP.value.(fucrm_genResDn[g,t]), JuMP.value.(fucrm_TotGenVioP[g,t]),
	          JuMP.value.(fucrm_TotGenVioN[g,t]), JuMP.value.(fucrm_MaxGenVioP[g,t]),
	          JuMP.value.(fucrm_MinGenVioP[g,t]) ), ',')
	  end # ends the loop
	end; # closes file
	# Write the peakers' schedules in CSV file
	open(".//outputs//FUCR_PeakerOutputs.csv", FILE_ACCESS_APPEND) do io
	  for t in 1:HRS_FUCR, k=1:PEAKERS
	      writedlm(io, hcat(day, t+INIT_HR_FUCR, k, df_peakers.UNIT_NAME[k],
	         df_peakers.MinPowerOut[k], df_peakers.MaxPowerOut[k],
	         JuMP.value.(fucrm_peakerOut[k,t]), JuMP.value.(fucrm_peakerOnOff[k,t]),
	         JuMP.value.(fucrm_peakerShutDown[k,t]), JuMP.value.(fucrm_peakerStartUp[k,t]),
	         JuMP.value.(fucrm_peakerResUp[k,t]), JuMP.value.(fucrm_peakerResNonSpin[k,t]),
	         JuMP.value.(fucrm_peakerResDn[k,t]) ), ',')
	  end # ends the loop
	end; # closes file

	# Writing storage units' optimal schedules in CSV file
	open(".//outputs//FUCR_StorageOutputs.csv", FILE_ACCESS_APPEND) do io
	   for t in 1:HRS_FUCR, p=1:STORG_UNITS
	      writedlm(io, hcat(day, t+INIT_HR_FUCR, p, df_storage.Name[p],
	          df_storage.Power[p], df_storage.Power[p]/df_storage.PowerToEnergRatio[p],
	          JuMP.value.(fucrm_storgChrg[p,t]), JuMP.value.(fucrm_storgDisc[p,t]),
	          JuMP.value.(fucrm_storgIdle[p,t]), JuMP.value.(fucrm_storgChrgPwr[p,t]),
	          JuMP.value.(fucrm_storgDiscPwr[p,t]), JuMP.value.(fucrm_storgSOC[p,t]),
	          JuMP.value.(fucrm_storgResUp[p,t]), JuMP.value.(fucrm_storgResDn[p,t]) ), ',')
	   end # ends the loop
	end; # closes file

	# Writing the transmission line flows in CSV file
	open(".//outputs//FUCR_TranFlowOutputs.csv", FILE_ACCESS_APPEND) do io
	  for t in 1:HRS_FUCR, n=1:N_ZONES, m=1:M_ZONES
	     writedlm(io, hcat(day, t+INIT_HR_FUCR, n,
	         JuMP.value.(fucrm_powerFlow[n,m,t]), TranC[n,m] ), ',')
	  end # ends the loop
	end; # closes file

	# Writing the curtilment, overgeneration, and spillage outcomes in CSV file
	open(".//outputs//FUCR_Curtail.csv", FILE_ACCESS_APPEND) do io
	  for t in 1:HRS_FUCR, n=1:N_ZONES
	     writedlm(io, hcat(day, t+INIT_HR_FUCR, n,
	         JuMP.value.(fucrm_OverGen[n,t]), JuMP.value.(fucrm_Demand_Curt[n,t]),
	         JuMP.value.(fucrm_windGSpil[n,t]), JuMP.value.(fucrm_solarGSpil[n,t]),
	         JuMP.value.(fucrm_hydroGSpil[n,t])), ',')
	  end # ends the loop
	end; # closes file

	t2_write_FUCRmodel_results = time_ns()

	time_write_FUCRmodel_results = (t2_write_FUCRmodel_results -t1_write_FUCRmodel_results)/1.0e9;
	@info "Write FUCRmodel results for day $day: $time_write_FUCRmodel_results executed in (s)";

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

	t1_FUCRtoBUCR1_data_hand = time_ns()

	for h=1:INIT_HR_SUCR-INIT_HR_FUCR
	  for g=1:GENS
		  uc_results.gens.onoff[g,h]= round(JuMP.value.(fucrm_genOnOff[g,h]));
		  uc_results.gens.power_out[g,h]= JuMP.value.(fucrm_genOut[g,h]);
		  uc_results.gens.startup[g,h]= round(JuMP.value.(fucrm_genStartUp[g,h]));
		  uc_results.gens.shutdown[g,h]= round(JuMP.value.(fucrm_genShutDown[g,h]));
		  for b=1:BLOCKS
			  uc_results.gens.genout_block[g,b,h]=JuMP.value.(fucrm_genOut_Block[g,b,h]);
		  end
	  end
	  for k=1:PEAKERS
		  uc_results.peakers.onoff[k,h]= round(JuMP.value.(fucrm_peakerOnOff[k,h]));
		  uc_results.peakers.power_out[k,h]=JuMP.value.(fucrm_peakerOut[k,h]);
		  uc_results.peakers.startup[k,h]= round(JuMP.value.(fucrm_peakerStartUp[k,h]));
		  uc_results.peakers.shutdown[k,h]= round(JuMP.value.(fucrm_peakerShutDown[k,h]));
		  for b=1:BLOCKS
			  uc_results.peakers.genout_block[k,b,h]=JuMP.value.(fucrm_peakerOut_Block[k,b,h]);
		  end
	  end
	  for p=1:STORG_UNITS
		  uc_results.storg.chrg[p,h]=round(JuMP.value.(fucrm_storgChrg[p,h]));
		  uc_results.storg.disc[p,h]=round(JuMP.value.(fucrm_storgDisc[p,h]));
		  uc_results.storg.idle[p,h]=round(JuMP.value.(fucrm_storgIdle[p,h]));
		  uc_results.storg.chrgpwr[p,h]=JuMP.value.(fucrm_storgChrgPwr[p,h]);
		  uc_results.storg.discpwr[p,h]=JuMP.value.(fucrm_storgDiscPwr[p,h]);
		  uc_results.storg.soc[p,h]=JuMP.value.(fucrm_storgSOC[p,h]);
	  end
	end

	t2_FUCRtoBUCR1_data_hand = time_ns();

	time_FUCRtoBUCR1_data_hand = (t2_FUCRtoBUCR1_data_hand -t1_FUCRtoBUCR1_data_hand)/1.0e9;
	@info "FUCRtoBUCR1 data handling for day $day executed in (s): $time_FUCRtoBUCR1_data_hand";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
		  writedlm(io, hcat("FUCRmodel", time_FUCRtoBUCR1_data_hand, "day: $day",
				  " ", "Pre-processing variables", "Data Manipulation"), ',')
	end; #closes file

	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:
using Gaugefields
using LinearAlgebra

function MDtest!(gauge_action,U,Dim,nn)
    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients.
    Uold = similar(U)
    dSdU = similar(U)

    substitute_U!(Uold,U)
    MDsteps = 10
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    numtrj = 100
    for itrj = 1:numtrj
        accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,nn,dSdU)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("$itrj plaq_t = $plaq_t")
        println("acceptance ratio ",numaccepted/itrj)
    end
end

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC #evaluate_GaugeAction(gauge_action,U) = tr(evaluate_GaugeAction_untraced(gauge_action,U))
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,nn,dSdU)

    Δτ = 1/MDsteps
    gauss_distribution!(p)

    Uout,Uout_multi,_ = calc_smearedU(U,nn)
    Sold = calc_action(gauge_action,Uout,p)

    substitute_U!(Uold,U)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action,dSdU,nn)

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end

    Uout,Uout_multi,_ = calc_smearedU(U,nn)
    Snew = calc_action(gauge_action,Uout,p)

    println("Sold = $Sold, Snew = $Snew")
    println("Snew - Sold = $(Snew-Sold)")
    ratio = min(1,exp(Snew-Sold))
    if rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action,dSdU,nn) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    factor =  -ϵ*Δτ/(NC)
    temps = get_temporary_gaugefields(gauge_action)
    Uout,Uout_multi,_ = calc_smearedU(U,nn)

    for μ=1:Dim
        calc_dSdUμ!(dSdU[μ],gauge_action,μ,Uout)
    end

    dSdUbare = back_prop(dSdU,nn,Uout_multi,U)

    for μ=1:Dim
        mul!(temps[1],U[μ],dSdUbare[μ]) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function test1()
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 1
    Dim = 4
    NC = 3

    U  =Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 5.7/2
    push!(gauge_action,β,plaqloop)

    show(gauge_action)

    L = [NX,NY,NZ,NT]
    nn = CovNeuralnet()
    ρ = [0.1]
    layername = ["plaquette"]
    st = STOUT_Layer(layername,ρ,L)
    push!(nn,st)
    #push!(nn,st)

    MDtest!(gauge_action,U,Dim,nn)

end

test1()


text:
using DataFileNames
import DataFileNames: parse_arg, _set_default_SEPS!
# using Random
using Test

@testset "DataFileNames.jl" begin

    # ------------------------------------------
    # Use default SEPS (mainly because ext sepataror)
    _set_default_SEPS!()

    # ------------------------------------------
    @info("Testing 'dfname' produces 'isvalid_dfname' names")

    for fn in [
        dfname("dat"),
        dfname("dat", "test"),
        dfname("dat", ".png"),
        dfname("dat", ".gz.tar"),
        dfname((;A = rand())),
        dfname((;A = rand()), "png"),
        dfname("dat", "dyn", (;ϵ = rand(), B = "main")),
        dfname("dat", (;A = :v1), "png"),
        dfname(["dir1", "dir2"], "dat", (;A = :v1), "png")
    ]
        isvalid = isvalid_dfname(fn)
        @info("Testing", fn, isvalid)
        @test isvalid
    end
    @test dfname() == ""
    @test dfname("") == ""
    println()

    # ------------------------------------------
    #  Test dfname
    @info("Testing dfname")
    @test dfname("dat") == "dat"
    @test dfname("dat") == dfname(:dat)

    @test dfname("dat", ".ext") != dfname("dat", "ext")
    @test dfname("dat", ".") == "dat."

    @test dfname("", "dat") == "dat"
    @test dfname("dat", "") == "dat"

    @test dfname("dat", (;A = 1)) == dfname("dat", :A => 1)
    @test dfname("dat", (;A = 1)) == dfname("dat", Dict(:A => 1))
    @test dfname("dat", "", (;A = 1)) == dfname("dat", (;A = 1))
    @test dfname("dat", (;A = 1), "") == dfname("dat", (;A = 1))
    @test dfname("dat", (;A = 1), ".ext") == dfname("dat", (;A = 1), "ext")

    @test basename(dfname(["dir"], "dat")) == "dat"
    @test dfname(["dir"], "dat") == joinpath("dir", "dat")
    @test dfname(["dir1", "dir2"], "dat") == joinpath("dir1", "dir2", "dat")
    @test dfname(["dir1"], ["dir2"], "dat") == joinpath("dir1", "dir2", "dat")
    println()

    # ------------------------------------------
    # Custom Type
    @info("Custom type")
    struct Foo
        f::Float64
        i::Int
        s::String
    end
    foo = Foo(1.0, 1, "hi")

    DataFileNames.parse_arg(f::Foo) = (;f.f, f.i, f.s)
    @test dfname(foo) == dfname(parse_arg(foo))
    @test dfname("bla", foo) == dfname("bla", parse_arg(foo))
    println()

    # ------------------------------------------
    # Parsing
    @info("Testing parse_dfname")
    let
        head = "dat"
        ext = ".jls"
        fname = dfname(head, ext)
        par_head, par_params, par_ext = parse_dfname(fname)
        @test head in par_head
        @test ext == par_ext
    end

    for it in 1:100
        head = [rand(-5:10), "Hello", rand()]
        params = (;S = "string", F = rand(), I = rand(-10:10))
        ext = ".ext"
        fname = dfname(head..., params, ext)
        par_head, par_params, par_ext = parse_dfname(fname)

        _isapprox(s1::String, s2::String) = (s1 == s2)
        _isapprox(v1, v2; atol = 1e-2) = isapprox(v1, v2; atol)

        @test length(head) == length(par_head)
        @test length(params) == length(par_params)

        @test all(_isapprox.(head, par_head))
        for (k, param) in pairs(params)
            kstr = string(k)
            @test haskey(par_params, kstr)
            par_param = par_params[kstr]
            @test all(_isapprox.(param, par_param))
        end
        @test length(par_ext) == length(ext) ?
            par_ext == ext : endswith(par_ext, ext)

        @test dfname(head..., params, ext) == dfname(par_head..., par_params, par_ext)
    end
    println()

    # ------------------------------------------
    # Try parse
    @test tryparse_dfname("Invalid=name") == nothing

    # ------------------------------------------
    # Create file
    @info("Testing creating a file")
    let
        file = dfname([tempdir()], "dat", (;I = rand(Int), F = rand()), "jls")
        @assert !isfile(file)
        @show file
        touch(file)
        @test isfile(file)
    end

end


text:
#===============================
<-----------Lathe.jl----------->
Programmed by Emmett Boudreau
    <emmett@emmettboudreau.com>
        <http://emmettboudreau.com>
MIT General Open Source License
    (V 3.0.0)
        Free for Modification and
        Redistribution
Thank you for your forks!
<-----------Lathe.jl----------->
38d8eb38-e7b1-11e9-0012-376b6c802672
#[deps]
DataFrames.jl
Random.jl
================================#
module Lathe
using DataFrames
using Random
#================
Stats
    Module
================#
module stats
#<----Mean---->
function mean(array)
    observations = length(array)
    average = sum(array)/observations
    return(average)
end
#<----Mode---->
function mode(array)
    m = findmax(array)
    return(m)
end
#<----Variance---->
function variance(array)
    me = mean(array)
    sq = sum(array) - me
    squared_mean = sq ^ 2
    return(squared_mean)
end
#<----Confidence Intervals---->
function confiints(data, confidence=.95)
    mean = mean(data)
    std = std(data)
    stderr = standarderror(data)
#    interval = stderr * scs.t.ppf((1 + confidence) / 2.0, n-1)
#    return (mean-interval, mean+interval)
end
#<----Standard Error---->
function standarderror(data)
    std = std(data)
    sample = length(data)
    ste = (std/sqrt(sample))
    return(ste)
end
#<----Standard Deviation---->
function std(array3)
    m = mean(array3)
    [i = (i-m) ^ 2 for i in array3]
    m = mean(array3)
    m = sqrt(m)
    return(m)
end
#<---- Correlation Coefficient --->
function correlationcoeff(x,y)
    n = length(x)
    yl = length(y)
    if n != yl
        throw(ArgumentError("The array shape does not match!"))
    end
#    sx = std(x)
#    sy = std(y)
#    x̄ = mean(x)
#    ȳ = mean(x)
#    [i = (i-x̄) / sx for i in x]
#    [i = (i-ȳ) / sy for i in y]
#    n1 = n-1
#    mult = x .* y
#    sq = sum(mult)
#    corrcoff = sq / n1
    xy = x .* y
    sx = sum(x)
    sy = sum(y)
    sxy = sum(xy)
    x2 = x .^ 2
    y2 = y .^ 2
    sx2 = sum(x2)
    sy2 = sum(y2)
    r = ((n*sxy) - (sx * sy)) / (sqrt((((n*sx2)-(sx^2)) * ((n*sy2)-(sy^2)))))
    return(r)
end
#<----Z score---->
function z(array)
    x̄ = mean(array)
    σ = std(array)
    return map(x -> (x - x̄) / σ, array)
end
#<----Quartiles---->
# - First
function firstquar(array)
    m = median(array)
    q15 = array / m
    q1 = array / m
    return(q)
end
# - Second(median)
function secondquar(array)
    m = median(array)
    return(m)
end
# - Third
function thirdquar(array)
    q = median(array)
    q = q * 1.5
end
# <---- Rank ---->
function getranks(array,rev = false)
    sortedar = sort(array,rev=rev)
    num = 1
    list = []
    for i in sortedar
        append!(list,i)
        num = num + 1
    end
    return(list)
end
#-------Inferential-----------__________
#<----Inferential Summary---->
function inf_sum(data,grdata)
    #Doing our calculations
    t = independent_t(data,grdata)
    f = f_test(data,grdata)
#    low,high = confiints(data)
    var = variance(data)
    grvar = variance(grdata)
    avg = mean(data)
    gravg = mean(grdata)
    sampstd = std(data)
    grstd = std(grdata)
    #Printing them out
    println("================")
    println("     Lathe.stats Inferential Summary")
    println("     _______________________________")
    println("N: ",length(grdata))
    println("x̅: ",avg)
    println("μ: ",gravg)
    println("s: ",sampstd)
    println("σ: ",grstd)
    println("var(x): ",var)
    println("σ2: ",grvar)
#    println("Low Confidence interval: ",low)
#    println("High Confidence interval: ",high)
    println("α ",t)
    println("Fp: ",f)
    println("================")
end
#<----T Test---->
# - Independent
function independent_t(sample,general)
    sampmean = mean(sample)
    genmean = mean(general)
    samples = length(sample)
    m = genmean
    [i = (i-m) ^ 2 for i in general]
    m = mean(general)
    m = sqrt(m)
    std = m
    t = (sampmean - genmean) / (std / sqrt(samples))
    return(t)
end
# - Paired
function paired_t(var1,var2)
    d = var1 .- var2
    d̄ = mean(x)
end
#<---- Correlations ---->
# - Spearman
function spearman(var1,var2)
    rgX = getranks(var1)
    rgY = getranks(var2)
    ρ = rgX*rgY / (std(rgX)*std(rgY))
    return(ρ)
end
# - Pearson
function pearson(x,y)
    sx = std(x)
    sy = std(y)
    x̄ = mean(x)
    ȳ = mean(x)
    [i = (i-x̄) / sx for i in x]
    [i = (i-ȳ) / sy for i in y]
    n1 = n-1
    mult = x .* y
    sq = sum(mult)
    corrcoff = sq / n1
    return(corrcoff)
end
# <---- Chi Distribution --->
function chidist(x,e)
# it is tough to calculate -> is it really needed?
    # A little less necessary, as its certainly not the most useful,
    # But this stats library could serve as a foundation for models that
    # Utilize Chi-Distributions, and although I wouldn't say having
    # A function to do so is urgent, It definitely would be cool,
    # Rather than an end user having to add another package just
    # To do one or two things, if you know what I mean.
    # But certainly there are other more important things to get through
    # Before 1.0, and I wouldn't consider any of these statistics incredibly
    # Necessary, but the template is there for what I want to include,
    # So people adding the module now can kindof know what to expect.
    # So hopefully that answers your question!
end
#<---- Chi-Square ---->
function chisq(var1,var2)
    chistat(obs, exp) = (obs - exp)^2/exp
    return chistat.(x, e) |> sum
end
#<---- ANOVA ---->
function anova(var1,var2)

end
#<---- Wilcoxon ---->
# - Wilcoxon Rank-Sum Test
function wilcoxrs(var1,var2)

end
function wilcoxsr(var1,var2)

end
#<---- Binomial Distribution ---->
function binomialdist(positives,size)
    # p = n! / x!(n-x!)*π^x*(1-π)^N-x
    n = size
    x = positives
    factn = factorial(n)
    factx = factorial(x)
    nx = factn / (factx * (n-x))
    return(nx)
end
#<---- Sign Test ---->
function sign(var1,var2)
    sets = var1 .- var2
    positives = []
    negatives = []
    zeros = []
    for i in sets
        if i == 0
            zeros.append(i)
        elseif i > 0
            positives.append(i)
        elseif i < 0
            negatives.append(i)
        end
    end
    totalpos = length(positives)
    totallen = length(sets)
    ans = binomialdist(positives,totallen)
    return(ans)
end
#<---- F-Test---->
function f_test(sample,general)
    totvariance = variance(general)
    sampvar = variance(sample)
    f =  sampvar / totvariance
    return(f)
end
#-------Bayesian--------------___________
#<----Bayes Theorem---->
#P = prob, A = prior, B = Evidence,
function bay_ther(p,a,b)
    psterior = (p*(b|a) * p*(a)) / (p*b)
    return(psterior)
end
function cond_prob(p,a,b)
    psterior = bay_ther(p,a,b)
    cond = p*(a|b)
    return(cond)
end

#---------------------------
end
#================
Model
    Validation
        Module
================#
module validate
#-------Model Metrics--------____________
using Lathe
## <---- Mean Absolute Error ---->
function mae(actual,pred)
    l = length(actual)
    lp = length(pred)
    if l != lp
        throw(ArgumentError("The array shape does not match!"))
    end
    result = actual-pred
    maeunf = Lathe.stats.mean(result)
    if maeunf < 0
        maeunf = (maeunf - maeunf) - maeunf
    end
    return(maeunf)
end
# <---- R Squared ---->
function r2(actual,pred)
    l = length(actual)
    lp = length(pred)
    if l != lp
        throw(ArgumentError("The array shape does not match!"))
    end
    r = Lathe.stats.correlationcoeff(actual,pred)
    rsq = r^2
    rsq = rsq * 100
    return(rsq)
end
function binomialdistribution(actual,pred)
    # p = n! / x!(n-x!)*π^x*(1-π)^N-x
    Lathe.stats.binomialdist(pos,neg,tot)
end
# --- Get Permutation ---
function getPermutation(model)

end
#--------------------------------------------
# End
end
#================
Preprocessing
     Module
================#
module preprocess
using Random
using Lathe
#===============
Generalized
    Data
        Processing
===============#
# Train-Test-Split-----
function TrainTestSplit(df,at = 0.75)
    sample = randsubseq(1:size(df,1), at)
    trainingset = df[sample, :]
    notsample = [i for i in 1:size(df,1) if isempty(searchsorted(sample, i))]
    testset = df[notsample, :]
    return(trainingset,testset)
end
# Array-Split ----------
function ArraySplit(data, at = 0.7)
    n = length(data)
    idx = Random.shuffle(1:n)
    train_idx = view(idx, 1:floor(Int, at*n))
    test_idx = view(idx, (floor(Int, at*n)+1):n)
    data[train_idx,:], data[test_idx,:]
    return(test_idx,train_idx)
end
# Sort-Split -------------
function SortSplit(data, at = 0.25, rev=false)
  n = length(data)
  sort!(data, rev=rev)  # Sort in-place
  train_idx = view(data, 1:floor(Int, at*n))
  test_idx = view(data, (floor(Int, at*n)+1):n)
  return(test_idx,train_idx)
end
# Unshuffled Split ----
function Uniform_Split(data, at = 0.7)
    n = length(data)
    idx = data
    train_idx = view(idx, 1:floor(Int, at*n))
    test_idx = view(idx, (floor(Int, at*n)+1):n)
    data[train_idx,:], data[test_idx,:]
    return(test_idx,train_idx)
end
#=======
Numerical
    Scaling
=======#
# ---- Rescalar (Standard Deviation) ---
function Rescalar(array)
    v = []
    min = minimum(array)
    max = maximum(array)
    for i in array
        x = (i-min) / (max - min)
        append!(v,x)
    end
    return(v)
end
# ---- Arbitrary Rescalar ----
function ArbitraryRescale(array)
    v = []
    a = minimum(array)
    b = maximum(array)
    for i in array
        x = a + ((i-a*i)*(b-a)) / (b-a)
        append!(v,x)
    end
    return(v)
end
# ---- Mean Normalization ----
function MeanNormalization(array)
    avg = Lathe.stats.mean(array)
    v = []
    a = minimum(array)
    b = maximum(array)
    for i in array
        m = (i-avg) / (b-a)
        append!(v,m)
    end
end
# ---- Z Normalization ----
function StandardScalar(array)
    q = Lathe.stats.std(array)
    avg = Lathe.stats.mean(array)
    v = []
    for i in array
        y = (i-avg) / q
        append!(v,y)
    end
    return(v)
end
# ---- Unit L-Scale normalize ----
function Unit_LScale(array)

end
#==========
Categorical
    Encoding
==========#
# <---- One Hot Encoder ---->
function OneHotEncode(array::Number)
    flatarr = Iterators.flatten(array)
    len = size(flatarr, 2)
    poslen = size(unique(flatarr), 2)
    out = Array{Number}(undef, len, poslen)
    for i in 1:len
        el = flatarr[i]
        idx = findall(x -> x == el, flatarr |> unique)[1][2]
        out[i, idx] = 1
    end
    return(out)
end
#-----------------------------
end
#================
Predictive
    Learning
        Models
================#
module models
#==
Baseline
    Model
==#
using Lathe
using Random
#Show models shows all the models that are stable
#And ready for use in the library
function showmodels()
    println("    Lathe.JL    ")
    println("________________")
    println("Current")
    println("    Usable")
    println("       Models")
    println("================")
    println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Continuous Models")
    println("-----------------")
    println("_____Linear_____")
    println("meanBaseline(y)")
    println("LinearLeastSquare(x,y,Type)")
    println("LinearRegression(x,y)")
    println("-----------------")
    println("___Non-Linear___")
    println("FourSquare(x,y)")
    println("ExponentialScalar(x,y)")
end
#Takes model, and X to predict, and returns a y prediction
function predict(m,x)
    if typeof(m) == FourSquare
        y_pred = pred_foursquare(m,x)
    end
    if typeof(m) == majBaseline
        y_pred = pred_catbaseline(m,x)
    end
    if typeof(m) == RegressionTree
        y_pred = pred_regressiontree(m,x)
    end
    if typeof(m) == LinearRegression
        y_pred = pred_LinearRegression(m,x)
    end
    if typeof(m) == meanBaseline
        y_pred = pred_meanbaseline(m,x)
    end
    if typeof(m) == RidgeRegression
        y_pred = pred_ridgeregression(m,x)
    end
    if typeof(m) == LinearLeastSquare
        y_pred = pred_linearleastsquare(m,x)
    end
    if typeof(m) == LogisticRegression
        y_pred = pred_logisticregression(m,x)
    end
    if typeof(m) == ExponentialScalar
        y_pred = pred_exponentialscalar(m,x)
    end
    if typeof(m) == MultipleLinearRegression
        y_pred = pred_multiplelinearregression(m,x)
    end
    return(y_pred)
end
# The help function:
function help(args)
    if typeof(args) == LogisticRegression
        println("Logistic")
    elseif typeof(args) == meanBaseline
        println("MeanBaseline")
    end
end
#======================================================================
=======================================================================
            CONTINUOS MODELS               CONTINUOS MODELS
            CONTINUOS MODELS               CONTINUOS MODELS
======================================================================
======================================================================#
#==
Mean
    Baseline
==#
 # Model Type
mutable struct meanBaseline
    y
end
#----  Callback
function pred_meanbaseline(m,xt)
    e = []
    m = Lathe.stats.mean(m.y)
    print("-Lathe.models Mean Baseline-")
    print("mean: ",m)
    for i in xt
        append!(e,m)
    end
    return(e)
end
#==
Multi-
    Gap
 - A quad range predictor, on steroids. -
==#
# Model Type
mutable struct RegressionTree
    x
    y
    n_divisions
    divisionsize
end
#----  Callback
function pred_regressiontree(m,xt)
    # x = q1(r(floor:q1)) |x2 = q2(r(q1:μ)) |x3 = q3(r(q2:q3)) |x4 q4(r(q3:cieling))
    # y' = q1(x * (a / x)) | μ(x * (a / x2)) | q3(x * (a / x3) | q4(x * (a / x4))
    # Original 4 quartile math ^^
        x = m.x
        y = m.y
        xtcopy = xt
        divs = m.n_divisions
        size = m.divisionsize
        # Go ahead and throw an error for the wrong input shape:
        xlength = length(x)
        ylength = length(y)
        if xlength != ylength
            throw(ArgumentError("The array shape does not match!"))
        end
        # Now we also need an error for when the total output of the
        #    division size and n divisions is > 100 percent
        divisions = size * divs
        if divisions != 1
            throw(ArgumentError("Invalid hyperparameters!: divisions * number of
            divisions must be = to 100 percent!"))
        end
        # Empty list
        e = []
        while divs > 0
            predictorx,x = Lathe.preprocess.SortSplit(x,size)
            predictory,y = Lathe.preprocess.SortSplit(y,size)
            predictorxt,xtcopy = Lathe.preprocess.SortSplit(xtcopy,size)
            currentrange = (minimum(predictorxt):maximum(predictorxt))
            linregmod = LinearRegression(predictorx,predictory)
            # Recursion replacement method:
            [predict(LinearRegression(predictorx,
            predictory),x) for x in currentrange]
            divs = divs - 1
        end
        return(xt)
end
#==
Four
    Square
==#
# Model Type
mutable struct FourSquare
    x
    y
end
#----  Callback
function pred_foursquare(m,xt)
    # x = q1(r(floor:q1)) |x2 = q2(r(q1:μ)) |x3 = q3(r(q2:q3)) |x4 q4(r(q3:cieling))
    # y' = q1(x * (a / x)) | μ(x * (a / x2)) | q3(x * (a / x3) | q4(x * (a / x4))
        x = m.x
        y = m.y
        # Go ahead and throw an error for the wrong input shape:
        xlength = length(x)
        ylength = length(y)
        if xlength != ylength
            throw(ArgumentError("The array shape does not match!"))
        end
        # Our empty Y prediction list==
        e = []
        # Quad Splitting the data ---->
        # Split the Y
        y2,range1 = Lathe.preprocess.SortSplit(y)
        y3,range2 = Lathe.preprocess.SortSplit(y2)
        y4,range3 = Lathe.preprocess.SortSplit(y3)
        y5,range4 = Lathe.preprocess.SortSplit(y4)
        yrange5 = y5
        # Split the x train
        x1,xrange1 = Lathe.preprocess.SortSplit(x)
        x2,xrange2 = Lathe.preprocess.SortSplit(x1)
        x3,xrange3 = Lathe.preprocess.SortSplit(x2)
        x4,xrange4 = Lathe.preprocess.SortSplit(x3)
        xrange5 = y5
        # Fitting the 4 linear regression models ---->
        regone = LinearLeastSquare(xrange1,range1)
        regtwo = LinearLeastSquare(xrange2,range2)
        regthree = LinearLeastSquare(xrange3,range3)
        regfour = LinearLeastSquare(xrange4,range4)
        regfive = LinearLeastSquare(xrange5,yrange5)
        # Split the train Data
        xt1,xtrange1 = Lathe.preprocess.SortSplit(xt)
        xt2,xtrange2 = Lathe.preprocess.SortSplit(xt1)
        xt3,xtrange3 = Lathe.preprocess.SortSplit(xt2)
        xt4,xtrange4 = Lathe.preprocess.SortSplit(xt3)
        xtrange5 = xt4
        # Get min-max
        xtrange1min = minimum(xtrange1)
        xtrange1max = maximum(xtrange1)
        xtrange2min = minimum(xtrange2)
        xtrange2max = maximum(xtrange2)
        xtrange3min = minimum(xtrange3)
        xtrange3max = maximum(xtrange3)
        xtrange4min = minimum(xtrange4)
        xtrange4max = maximum(xtrange4)
        xtrange5min = minimum(xtrange5)
        # Ranges for ifs
        condrange1 = (xtrange1min:xtrange1max)
        condrange2 = (xtrange2min:xtrange2max)
        condrange3 = (xtrange3min:xtrange3max)
        condrange4 = (xtrange4min:xtrange4max)
        # This for loop is where the dimension's are actually used:
        for i in xt
            if i in condrange1
                ypred = predict(regone,i)
            elseif i in condrange2
                ypred = predict(regtwo,i)
            elseif i in condrange3
                ypred = predict(regthree,i)
            elseif i in condrange4
                ypred = predict(regfour,i)
            else
                ypred = predict(regfive,i)
            end

            append!(e,ypred)
        end
        return(e)
end
#==
Isotonic
    Regression
==#
mutable struct IsotonicRegression
    x
    y
end
function pred_isotonicregression(m,xt)
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end
end
#==
Multiple
    Linear
        Regression
==#
mutable struct MultipleLinearRegression
    x
    y
end
function pred_multiplelinearregression(m,xt)
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end
    if length(m.x) != length(xt)
        throw(ArgumentError("Bad Feature Shape |
        Training Features are not equal!",))
    end
    y_pred = []
    for z in xt
        predavg = []
        for i in matrice
            m = LinearRegression(i,y)
            pred = predict(m,z)
            append!(predavg,pred)
        end
        mn = Lathe.stats.mean(predavg)
        append!(y_pred,mn)
        return(y_pred)
    end
end
#==
Linear
    Regression
==#
mutable struct LinearRegression
    x
    y
end
#----  Callback
function pred_LinearRegression(m,xt)
    # a = ((∑y)(∑x^2)-(∑x)(∑xy)) / (n(∑x^2) - (∑x)^2)
    # b = (x(∑xy) - (∑x)(∑y)) / n(∑x^2) - (∑x)^2
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end
    # Get our x and y as easier variables
    x = m.x
    y = m.y
    # Get our Summatations:
    Σx = sum(x)
    Σy = sum(y)
    # dot x and y
    xy = x .* y
    # ∑dot x and y
    Σxy = sum(xy)
    # dotsquare x
    x2 = x .^ 2
    # ∑ dotsquare x
    Σx2 = sum(x2)
    # n = sample size
    n = length(x)
    # Calculate a
    a = (((Σy) * (Σx2)) - ((Σx * (Σxy)))) / ((n * (Σx2))-(Σx^2))
    # Calculate b
    b = ((n*(Σxy)) - (Σx * Σy)) / ((n * (Σx2)) - (Σx ^ 2))
    [i = a+(b*i) for i in xt]
    return(xt)
end
#==
Linear
    Least
     Square
==#
mutable struct LinearLeastSquare
    x
    y
    Type
end
function pred_linearleastsquare(m,xt)
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end
    if m.Type == :REG
        x = m.x
        y = m.y
        # Summatation of x*y
        xy = x .* y
        sxy = sum(xy)
        # N
        n = length(x)
        # Summatation of x^2
        x2 = x .^ 2
        sx2 = sum(x2)
        # Summatation of x and y
        sx = sum(x)
        sy = sum(y)
        # Calculate the slope:
        slope = ((n*sxy) - (sx * sy)) / ((n * sx2) - (sx)^2)
        # Calculate the y intercept
        b = (sy - (slope*sx)) / n
        # Empty prediction list:
        y_pred = []
        for i in xt
            pred = (slope*i)+b
            append!(y_pred,pred)
        end
    end
    if m.Type == :OLS

    end
    if m.Type == :WLS

    end
    if m.Type == :GLS

    end
    if m.Type == :GRG

    end
    return(y_pred)
end
#==
Ridge
    Regression
==#
mutable struct RidgeRegression
    x
    y
end
function pred_ridgeregression(m,xt)
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end
end
#==
Logistic
    Regression
==#
mutable struct LogisticRegression
    x
    y
end
function pred_logisticregression(m,xt)
    if length(m.x) != length(m.y)
        throw(ArgumentError("The array shape does not match!"))
    end

end
#==
Linear
    Scalar
==#
mutable struct ExponentialScalar
    x
    y
end
function pred_exponentialscalar(m,xt)
    x = m.x
    y = m.y
    at = 0.25
    xdiv1,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv2,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv3,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv4,x = Lathe.preprocess.SortSplit(x,.05)
    ydiv1,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv2,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv3,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv4,y = Lathe.preprocess.SortSplit(y,.05)
    scalarlist1 = ydiv1 ./ xdiv1
    scalarlist2 = ydiv2 ./ xdiv2
    scalarlist3 = ydiv3 ./ xdiv3
    scalarlist4 = ydiv3 ./ xdiv3
    xdiv1,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv2,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv3,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv4,x = Lathe.preprocess.SortSplit(x,.05)
    ydiv1,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv2,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv3,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv4,y = Lathe.preprocess.SortSplit(y,.05)
    scalarlist6 = ydiv1 ./ xdiv1
    scalarlist7 = ydiv2 ./ xdiv2
    scalarlist8 = ydiv3 ./ xdiv3
    scalarlist9 = ydiv3 ./ xdiv3
    xdiv1,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv2,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv3,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv4,x = Lathe.preprocess.SortSplit(x,.05)
    ydiv1,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv2,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv3,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv4,y = Lathe.preprocess.SortSplit(y,.05)
    scalarlist10 = ydiv1 ./ xdiv1
    scalarlist11 = ydiv2 ./ xdiv2
    scalarlist12 = ydiv3 ./ xdiv3
    scalarlist13 = ydiv3 ./ xdiv3
    xdiv1,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv2,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv3,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv4,x = Lathe.preprocess.SortSplit(x,.05)
    ydiv1,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv2,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv3,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv4,y = Lathe.preprocess.SortSplit(y,.05)
    scalarlist14 = ydiv1 ./ xdiv1
    scalarlist15 = ydiv2 ./ xdiv2
    scalarlist16 = ydiv3 ./ xdiv3
    scalarlist17 = ydiv3 ./ xdiv3
    xdiv1,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv2,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv3,x = Lathe.preprocess.SortSplit(x,.05)
    xdiv4,x = Lathe.preprocess.SortSplit(x,.05)
    ydiv1,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv2,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv3,y = Lathe.preprocess.SortSplit(y,.05)
    ydiv4,y = Lathe.preprocess.SortSplit(y,.05)
    scalarlist18 = ydiv1 ./ xdiv1
    scalarlist19 = ydiv2 ./ xdiv2
    scalarlist20 = y ./ x
    # Now we sortsplit the x train
    xtdiv1,xt2 = Lathe.preprocess.SortSplit(xt,.05)
    xtdiv2,xt2 = Lathe.preprocess.SortSplit(xt2,.05)
    xtdiv3,xt2 = Lathe.preprocess.SortSplit(xt2,.05)
    xtdiv4,null = Lathe.preprocess.SortSplit(xt2,.05)
    range1 = minimum(xtdiv1):maximum(xtdiv1)
    range2 = minimum(xtdiv2):maximum(xtdiv2)
    range3 = minimum(xtdiv3):maximum(xtdiv3)
    range4 = minimum(xtdiv4):maximum(xtdiv4)
    range5 = minimum(null):maximum(null)
    returnlist = []
    for i in xt
        if i in range1
            res = i * rand(scalarlist1)
            append!(returnlist,res)
        elseif i in range2
            predlist = []
            res = i * rand(scalarlist2)
            append!(returnlist,res)

        elseif i in range3
            predlist = []
            res = i * rand(scalarlist3)
            append!(returnlist,res)
        elseif i in range4
            predlist = []
            res = i * rand(scalarlist4)
            append!(returnlist,res)
        else
            predlist = []
            res = i * rand(scalarlist20)
            append!(returnlist,res)
        end
    end
    return(returnlist)
end
#======================================================================
=======================================================================
            CATEGORICAL MODELS             CATEGORICAL MODELS
            CATEGORICAL MODELS             CATEGORICAL MODELS
======================================================================
======================================================================#
#==
Majority
    Class
        Baseline
==#
# Model Type
mutable struct majBaseline
    y
end
#----  Callback
function pred_catbaseline(m,xt)
    y = m.y
    e = []
    mode = Lathe.stats.mode(xt)
    for i in xt
        append!(e,i)
    end

end
#
#----------------------------------------------
end
#================
Pipeline
    Module
================#
module Pipelines

# Note to future self, or other programmer:
# It is not necessary to store these as constructors!
# They can just be strings, and use the model's X and Y!
using Lathe
mutable struct Pipeline
    model
    categoricalenc
    contenc
    imputer
end
mutable struct fitpipeline
    pipeline
    catx
    conx
    y
end
function pipe_predict(fitpipeline,xtcats,xtcons)
    """ Takes a fit pipeline, and an X and predicts. """
    fitpipeline.conx = fitpipeline.pipeline.contenc(x)
    ypr = Lathe.models.predict(fitpipeline.pipeline.model(fitpipeline.conx,
    fitpipeline.y),xt)
    return(ypr)
end
function serialize(fitpipeline,filename)
    """Outputs pipeline as sav file."""
end
#----------------------------------------------
end
#==
This is the end of the main
module, nothing is to be written
beyond here
==#
end


text:
abstract type AbstractBinningAlgorithm end

struct Bin{T<:Real}
    """Number of samples."""
    nsamples::Ref{Int}
    """Sum of predictions."""
    prediction_sum::Vector{T}
    """Counts of labels."""
    label_counts::Vector{Int}

    function Bin{T}(nsamples::Ref{Int}, prediction_sum::Vector{T},
                    label_counts::Vector{Int}) where T
        nsamples[] ≥ 0 || throw(ArgumentError("the number of samples must be non-negative"))
        nclasses = length(prediction_sum)
        nclasses > 1 || throw(ArgumentError("the number of classes must be greater than 1"))
        nclasses == length(label_counts) ||
            throw(DimensionMismatch("the number of predicted classes has to be equal to the number of classes"))

        new{T}(nsamples, prediction_sum, label_counts)
    end
end

Bin(nsamples::Ref{Int}, prediction_sum::Vector{<:Real}, label_counts::Vector{Int}) =
    Bin{eltype(prediction_sum)}(nsamples, prediction_sum, label_counts)

"""
    Bin(predictions, labels)

Create bin of `predictions` and corresponding `labels`.
"""
function Bin(predictions::AbstractMatrix{<:Real}, labels::AbstractVector{<:Integer})
    # check arguments
    nclasses, nsamples = check_size(predictions, labels)

    # compute sum of predictions
    prediction_sum = vec(sum(predictions, dims=2))

    # compute counts of labels
    label_counts = counts(labels, 1:nclasses)

    Bin(Ref(nsamples), prediction_sum, label_counts)
end

"""
    Bin(T, m::Int)

Create empty container that keeps track of the sum of predictions of type `T` and the
counts of labels of a model with `m` classes.
"""
Bin(::Type{T}, m::Int) where T = Bin(Ref(0), zeros(T, m), zeros(Int, m))

"""
    adddata!(bin::Bin, predictions, labels)

Update running statistics of the `bin` by integrating additional pair(s) of `prediction`s
and `label`s.
"""
function adddata!(bin::Bin, prediction::AbstractVector{<:Real}, label::Integer)
    @unpack nsamples, prediction_sum, label_counts = bin

    # check arguments
    length(prediction) == length(prediction_sum) ||
        throw(DimensionMismatch("the number of predicted classes is incorrect"))

    # update number of samples
    nsamples[] += 1

    # update sum of predictions
    prediction_sum .+= prediction

    # update counts of labels
    if 1 ≤ label ≤ length(label_counts)
        @inbounds label_counts[label] += 1
    end

    nothing
end

function adddata!(bin::Bin, predictions::AbstractMatrix{<:Real},
                  labels::AbstractVector{<:Integer})
    @unpack nsamples, prediction_sum, label_counts = bin

    # check arguments
    addnclasses, addnsamples = check_size(predictions, labels)
    addnclasses == length(prediction_sum) ||
        throw(DimensionMismatch("the number of predicted classes is incorrect"))

    # update number of samples
    nsamples[] += addnsamples

    # update sum of predictions
    for p in eachcol(predictions)
        prediction_sum .+= p
    end

    # update counts of labels
    addcounts!(label_counts, labels, 1:addnclasses)

    nothing
end

# distance computations

"""
    bin_eval(bin::Bin, distance)

Evaluate `distance` between the average prediction and the distribution of labels in the
`bin`, multiplied by the number of samples.

The multiplication with the number of samples simplifies the computation of the expected
miscalibration error (ECE)[@ref].
"""
function bin_eval(bin::Bin, distance::Union{TotalVariation,Cityblock,Euclidean,SqEuclidean})
    @unpack nsamples, prediction_sum, label_counts = bin

    # compute output type
    T = typeof(bin_eval_end(distance,
                            zero(result_type(distance, prediction_sum, label_counts)), 1))

    # if bin is empty return 0
    iszero(nsamples[]) && return zero(T)

    # compute distance between sum of predictions and counts of labels
    d = evaluate(distance, prediction_sum, label_counts)

    # perform normalization (e.g., in the case of squared Euclidean distance)
    bin_eval_end(distance, d, nsamples[])::T
end

# normalization is required since the distance is computed between sum of predictions
# and counts of labels instead of between the mean of predictions and distribution of labels
# thus in the case of the squared Euclidean distance without normalization, the result
# would be multiplied by the squared number of samples instead of only the number of
# samples.
bin_eval_end(::Union{TotalVariation,Cityblock,Euclidean}, d::Real, ::Int) = d
bin_eval_end(::SqEuclidean, d::Real, nsamples::Int) = d * inv(nsamples)


text:
function convertvideo(file; pathout::String=splitdir(file)[1], filenameout::String=splitext(splitdir(file)[2])[1])
    file = abspath(file)
    pathout = abspath(pathout)
    filenameout = filenameout*".mp4"
    println("Creating " * filenameout * ".")
    MeshCat.convert_frames_to_video(file, joinpath(pathout, filenameout), overwrite=true)
    return
end


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

name = "HYPRE"
version = v"2.23.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/hypre-space/hypre/archive/refs/tags/v$(version).tar.gz", "8a9f9fb6f65531b77e4c319bf35bfc9d34bf529c36afe08837f56b635ac052e2")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/hypre-*
cd src/

mkdir build
cd build/

CMAKE_FLAGS=()

#help find mpi on mingw, subset of https://github.com/JuliaPackaging/Yggdrasil/blob/b4fdb545c3954cff218051d7520c7418991d3416/T/TauDEM/build_tarballs.jl#L28-L53
if [[ "$target" == x86_64-w64-mingw32 ]]; then
    CMAKE_FLAGS+=(
        -DMPI_HOME=${prefix}
        -DMPI_GUESS_LIBRARY_NAME=MSMPI
    )
    if [[ "${target}" == x86_64-* ]]; then
        for lang in C CXX; do
            CMAKE_FLAGS+=(-DMPI_${lang}_LIBRARIES=msmpi64)
        done
    fi
fi

cmake .. \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DHYPRE_ENABLE_SHARED=ON \
-DHYPRE_ENABLE_HYPRE_BLAS=ON \
-DHYPRE_ENABLE_HYPRE_LAPACK=ON \
-DHYPRE_ENABLE_CUDA_STREAMS=OFF \
-DHYPRE_ENABLE_CUSPARSE=OFF \
-DHYPRE_ENABLE_CURAND=OFF \
"${CMAKE_FLAGS[@]}"

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS_jll", uuid="4536629a-c528-5b80-bd46-f80d51c5b363"))
    Dependency(PackageSpec(name="LAPACK_jll", uuid="51474c39-65e3-53ba-86ba-03b1b862ec14"))
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4"))
    Dependency(PackageSpec(name="MicrosoftMPI_jll", uuid="9237b28f-5490-5468-be7b-bb81f5f5e6cf"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    LazyString <: AbstractString

A lazy representation of string interpolation. This is useful when a string
needs to be constructed in a context where performing the actual interpolation
and string construction is unnecessary or undesirable (e.g. in error paths
of functions).

This type is designed to be cheap to construct at runtime, trying to offload
as much work as possible to either the macro or later printing operations.

!!! compat "Julia 1.8"
    `LazyString` requires Julia 1.8 or later.
"""
mutable struct LazyString <: AbstractString
    parts::Tuple
    # Created on first access
    str::String
    LazyString(args...) = new(args)
end

"""
    lazy"str"

Create a [`LazyString`](@ref) using regular string interpolation syntax.
Note that interpolations are *evaluated* at LazyString construction time,
but *printing* is delayed until the first access to the string.

!!! compat "Julia 1.8"
    `lazy"str"` requires Julia 1.8 or later.
"""
macro lazy_str(text)
    parts = Any[]
    lastidx = idx = 1
    while (idx = findnext('$', text, idx)) !== nothing
        lastidx < idx && push!(parts, text[lastidx:idx-1])
        idx += 1
        expr, idx = Meta.parseatom(text, idx; filename=string(__source__.file))
        push!(parts, esc(expr))
        lastidx = idx
    end
    lastidx <= lastindex(text) && push!(parts, text[lastidx:end])
    :(LazyString($(parts...)))
end

function String(l::LazyString)
    if !isdefined(l, :str)
        l.str = sprint() do io
            for p in l.parts
                print(io, p)
            end
        end
    end
    return l.str
end

hash(s::LazyString, h::UInt64) = hash(String(s), h)
lastindex(s::LazyString) = lastindex(String(s))
iterate(s::LazyString) = iterate(String(s))
iterate(s::LazyString, i::Integer) = iterate(String(s), i)
isequal(a::LazyString, b::LazyString) = isequal(String(a), String(b))
==(a::LazyString, b::LazyString) = (String(a) == String(b))
ncodeunits(s::LazyString) = ncodeunits(String(s))
codeunit(s::LazyString) = codeunit(String(s))
codeunit(s::LazyString, i::Integer) = codeunit(String(s), i)
isvalid(s::LazyString, i::Integer) = isvalid(String(s), i)

max_stars_repo_name: lawless-m/julia

max_issues_repo_name: lawless-m/julia

max_forks_repo_name: lawless-m/julia


text:
using FileRules
using Test

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


text:
# Testing part 1
using Test
using ReferenceTests
using ParallelStencil
using ParallelStencil.FiniteDifferences3D

if !ParallelStencil.is_initialized()
    @init_parallel_stencil(Threads, Float64, 3)
end

if !@isdefined USE_GPU
    const USE_GPU = false
end
@test @isdefined USE_GPU
include("../scripts-part1/part1_kernel_programming.jl") # modify to include the correct script
include("../scripts-part1/part1_array_programming.jl") # modify to include the correct script

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(d1[k], d2[k], atol=1e-5) for k in keys(d1) ])

MPI.Init()

Xc_g, H_g = diffusion_3D_array_programming(nx=32, ny=32, nz=32, init_and_finalize_MPI=false, verbose=false)
inds = Int.(ceil.(LinRange(1, length(Xc_g), 12)))
d_array = Dict(:X=> Xc_g[inds], :H=>H_g[inds, inds, 15])

Xc_g, H_g, _ = diffusion_3D_kernel_programming(nx=32, ny=32, nz=32, use_shared_memory=true, init_and_finalize_MPI=false, verbose=false)
inds = Int.(ceil.(LinRange(1, length(Xc_g), 12)))
d_kernel_shared_memory = Dict(:X=> Xc_g[inds], :H=>H_g[inds, inds, 15])

Xc_g, H_g, _ = diffusion_3D_kernel_programming(nx=32, ny=32, nz=32, use_shared_memory=false, init_and_finalize_MPI=false, verbose=false)
inds = Int.(ceil.(LinRange(1, length(Xc_g), 12)))
d_kernel = Dict(:X=> Xc_g[inds], :H=>H_g[inds, inds, 15])

if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    @testset "Ref-file" begin
        @test_reference "reftest-files/test_1.bson" d_array by=comp
        @test_reference "reftest-files/test_1.bson" d_kernel by=comp
        @test_reference "reftest-files/test_1.bson" d_kernel_shared_memory by=comp
    end
end

MPI.Finalize()
@ParallelStencil.reset_parallel_stencil()


text:
Pkg.update()
Pkg.add("DataFrames")
Pkg.add("SQLite")
Pkg.add("Gtk")
Pkg.update()

max_stars_repo_name: DarienFlamont/Shogi

max_issues_repo_name: DarienFlamont/Shogi

max_forks_repo_name: DarienFlamont/Shogi


text:
@deprecate autodiff(mode, block) markdiff(block)
@deprecate autodiff(mode) markdiff
@deprecate opdiff faithful_grad
@deprecate statdiff faithful_grad
@deprecate QFTCircuit(n) qft_circuit(n)

struct QFTBlock{N} <: PrimitiveBlock{N}
    @deprecate QFTBlock{n}() where n qft(n)
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
# Code for 1D example with highly oscillating coefficient,
# using AMEn solver via Python interface (see TTPy.jl)
module MultiscaleAMEn

using TT, TTFEM, TTAlg, TTPy

export amenOscillatory1D

function amenOscillatory1D(L::Integer, ε::AbstractFloat, K::Number)
    Dinv = ttdiagm(add(2., ttones(L), 1., ttcos(L, K*π, .5)))
    Drhs = ttones(L)
    d, _ = stsolve(L, 1, x->add(1., matvec(Dinv,x), -1., Drhs), 3., 1., .5*ɛ, .75)
    D = ttdiagm(d)
    C = ttBPXdnC(L, 1)
    F = ttBPXdnQLaplace(L,1,1)
    p = size(F[L],1)
    q = size(F[L],4)
    T = reshape(F[L],(p*4,q))
    S = reshape(F[L+1],(q,1))
    F[L] = reshape(T*S,(p,2,2,1))
    pop!(F)
    Ft = transpose(F)
    f = ttrhsdn(L)
    g = matvec(C, f)
    svdtrunc!(g, svd!(g), .5*ɛ)
    B = matmat(Ft, matmat(D, F))
    u = amen(B, g, g, ε)
    Cu = matvec(C, u)
    Fu = matvec(F, u)
    return u, Cu, Fu
end

end