text:
# simulator_quantum_2d_fractal.jl

"""
	Quantum2dFractalSimulator

# arguments

* `lattice::Quantum2dFractalLattice`
* `dimM::Int` the maximum dimension for the tensors
* `inititeration::Int` initial iteration (default=0)

# to do:
1. bug fixing - DONE
2. merge inititer to this

"""
type Quantum2dFractalSimulator{T} <: Quantum2dSimulator{T}
	dimM::Int
	wholeiteration::Int
	inititeration::Int
	countiteration::Int
	trottercount::Int
	normalizationFactor::Array{T,2}
	coefficients::Array{T,2}
	lattice::Quantum2dFractalLattice{T}
	infoCoefficients::Array{AbstractString,1}

	function Quantum2dFractalSimulator{T}(lattice::Quantum2dFractalLattice{T}, dimM::Int, inititeration::Int)
		this = new{T}()
		this.dimM = dimM
		this.inititeration = inititeration
		wholeiteration = getTrotteriteration(lattice) * 2 + inititeration
		this.wholeiteration = wholeiteration
		this.lattice = lattice
		this.infoCoefficients = ["t", "px", "py", "q",
					  "c", "ly", "ey",
					  "lc", "cl", "ec",
					  "lccl", "eccl",
					  "lx", "ex",
					  "lccll", "eccll"]
		numberOfCoefficients = getLengthOfNormList(this)
		this.normalizationFactor = ones(T,wholeiteration+1,numberOfCoefficients)
		this.coefficients = zeros(T,wholeiteration+1,numberOfCoefficients)
		initializeCoefficients!(this)
		initializeCount!(this)
		return this
	end
end

Quantum2dFractalSimulator{T}(lattice::Quantum2dFractalLattice{T}, dimM::Int, inititeration::Int) =
Quantum2dFractalSimulator{T}(lattice,dimM,inititeration)
Quantum2dFractalSimulator{T}(lattice::Quantum2dFractalLattice{T}, dimM::Int) =
Quantum2dFractalSimulator{T}(lattice,dimM, 0)

#--- refactor
function getLengthOfNormList(simulator::Quantum2dFractalSimulator)
	return length(simulator.infoCoefficients)
end

function getListOfNormFactors(simulator::Quantum2dFractalSimulator)
	return simulator.infoCoefficients
end

function getIndexOf(simulator::Quantum2dFractalSimulator, normname::AbstractString)
	listOfNormFactors = getListOfNormFactors(simulator)
	indexof = 0
	for (i,a) in enumerate(listOfNormFactors)
		if a == normname
			indexof = i
		end
	end
	indexof == 0 && error("input for getIndexOf is wrong")
	return indexof
end

#---
# run the simulator

function (simulator::Quantum2dFractalSimulator)(;printlog="none")
	initializeCount!(simulator) # set to be zero
	initializeCoefficients!(simulator)
	countUp!(simulator)
	normalizeTensor!(simulator)

	if printlog in ["coef", "norm"]
		printLog(simulator, printlog="label")
		printLog(simulator, printlog=printlog)
	end

	if getInititeration(simulator) >= 1
		for i = 1:getInititeration(simulator)
			countUp!(simulator, "trotter")
			renormalizeTrotter!(simulator, getDimM(simulator))
			updateCoefficients!(simulator, "trotter")
			printLog(simulator, printlog=printlog)
		end
	end

	while true
		countUp!(simulator)
		renormalizeSpace!(simulator, getDimM(simulator))
 		updateCoefficients!(simulator)
		printLog(simulator, printlog=printlog)

 		countUp!(simulator, "trotter")
 		renormalizeTrotter!(simulator, getDimM(simulator))
 		updateCoefficients!(simulator,"trotter")
		printLog(simulator, printlog=printlog)
 		if getCount(simulator) > getWholeiteration(simulator)
 			break
 		end
	end
	freeenergy = getFreeEnergy(simulator) ###### count the trotter into
	magnetization = getExpectationValue(simulator)
	return freeenergy, magnetization
end

#---
# functions for simulate

# functions about free energy
function getFreeEnergy{T}(simulator::Quantum2dFractalSimulator{T})
	numberofsites = sum(getCoefficient(simulator,0,1)[1:4])

	tensorT = getTensorT(simulator)[1]
	termcorrection = log(traceTensorTPeriodic(tensorT))

	freeenergy = - getFirstTerm(simulator) * (
											termFromNormFactor(simulator) +
										    termcorrection
										   ) / numberofsites
# 	println(getFirstTerm(simulator))
# 	println(termFromNormFactor(simulator))
# 	println(termcorrection)
# 	println(numberofsites)
	return freeenergy
end

function termFromNormFactor{T}(simulator::Quantum2dFractalSimulator{T})
	termNorm = zero(T)
	for i in 1:getCount(simulator)
		for j in 1:getLengthOfNormList(simulator)
			termNorm += getCoefficient(simulator, j, i) * log(getNorm(simulator, j, i))
		end
	end
	return termNorm
end

# coefficients

function initializeCoefficients!{T}(simulator::Quantum2dFractalSimulator{T})
	wholeiteration = getWholeiteration(simulator)
	setCoefficient!(simulator, 1.0, getIndexOf(simulator, "t"), wholeiteration + 1)
end

function updateCoefficients!{T}(simulator::Quantum2dFractalSimulator{T})
	iteration = getWholeiteration(simulator) +2 -getCount(simulator)
	# depends on the order of the storing of the coefficients!
	currentCoef = getCoefficient(simulator, 0, iteration+1)[1:4]
	evolveMatrix = # this is 16 x 4 matrix
	#	t	px	py	q
	[	4	4	4	4; #t
  		4	2	4	2; #px
		4	4	2	2; #py
		0	4	0	4; #q
		2	2	2	2; #c
		2	1	2	1; #ly
		0	1	0	1; #ey
		2	1	2	1; #lc
		0	0	0	0; #cl
		0	1	0	1; #ec
		1	0	1	0; #lccl
		0	1	0	1; #eccl
		2	2	1	1; #lx
		0	0	1	1; #ex
		1	0	1	0; #lccll
		0	1	0	1] #eccll
	nextCoef = evolveMatrix * currentCoef
	setCoefficients!(simulator, nextCoef, iteration)
end

# move them to 'simulator_fractal.jl' when it works
function setCoefficients!{T}(simulator::Quantum2dFractalSimulator{T}, nextCoef::Array{T,1}, iteration::Int; lengthNextCoef::Int = getLengthOfNormList(simulator))
	length(nextCoef) == lengthNextCoef || error("input length does not match")
	for i in 1:lengthNextCoef
		setCoefficient!(simulator, nextCoef[i], i, iteration)
	end
end

##### input: simulator
#=
from the input get
iteration: the current step
currentCoef: the current coefficients in vector (1-dim array) form
	(but what we need is only the four values,,, so, actually
	it does not work directly from the getCoef function. Either
		-del- change the getCoef function
		-change the results by limiting
		WE DO DEAL WITH THE RESULTS (but still check the func getCoef)
-
store a 2-dim array to get the way to evolve the coefficients
-
output a one-dim array with more values than the currentCoef
	-del- Do not, yet, set the values to the simulator?
Or, yes, do store the new values to the simulator?
WE DO SET THE VALUES TO THE SIMULATOR
=#

function updateCoefficients!{T}(simulator::Quantum2dFractalSimulator{T}, trotter::AbstractString)
	if trotter == "trotter"
		iteration = getWholeiteration(simulator) +2 -getCount(simulator)
		currentCoef = getCoefficient(simulator, 0, iteration+1)[1:4]
		newCoef = currentCoef .* 2.0
		setCoefficients!(simulator, newCoef, iteration, lengthNextCoef = 4)
	end
end

# normalization

"""
	normalizeTensor!{T}(simulator::Quantum2dFractalSimulator{T})
normalize T,Px,Py,Q tensors
set the normed tensors back to the lattice in simulator
set the Norm to normalizationFactor
"""
function normalizeTensor!{T}(simulator::Quantum2dFractalSimulator{T})
	tensorT, tensorTtilde = getTensorT(simulator)
	newTensorT, normT = normalizeTensor(tensorT)
	newTensorTtilde = tensorTtilde ./ normT
	setTensorT!(simulator, newTensorT)
	setTensorT!(simulator, newTensorTtilde; tilde = true)

	newTensorPx, normPx = normalizeTensor(getTensorP(simulator, 1))
	setTensorP!(simulator, 1, newTensorPx)
	newTensorPy, normPy = normalizeTensor(getTensorP(simulator, 2))
	setTensorP!(simulator, 2, newTensorPy)

	newTensorQ, normQ = normalizeTensor(getTensorQ(simulator))
	setTensorQ!(simulator, newTensorQ)

	setNorm!(simulator, normT, getIndexOf(simulator, "t"))
	setNorm!(simulator, normPx, getIndexOf(simulator, "px"))
	setNorm!(simulator, normPy, getIndexOf(simulator, "py"))
	setNorm!(simulator, normQ, getIndexOf(simulator, "q"))
end

#---
# functions for initialize

#---
# test function

function testQuantum2dFractalSimulator()
	# parameters:
	externalfield = 0.1
	temperature = 0.0
	trotterparameter = 0.01
	trotteriteration = 30
	dimM = 3
	# build simulator
	spininfo = SpinInfo("quantum_ising_2","asym",externalfield,temperature)
	trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
	spinmodel = QuantumIsingModel(spininfo, trotterinfo)
	lattice = buildLattice("quantum_2d_fractal_1", spinmodel)
	simulator = Quantum2dFractalSimulator(lattice, dimM)
	return simulator
end


text:
using DXFutils
using Test

@testset "groupcode" begin
    for (key, value) in DXFutils.group_code_registry
        @test key == groupcode(value)
        @test key == groupcode(one(value))
    end
end

@testset "DXF Groups" begin
    @test DXFutils.group_code_registry[0] == EntityType
    @test groupcode(EntityType) == 0
    grp = DXFutils.EntityType("SECTION")
    @test groupcode(grp) == 0
    @test grp.value == "SECTION"
end

example_file = "c:/Users/Mark Nahabedian/crafts/crafts/TrapezoidalLapJoint/drawings/strut-side-view.dxf"

@testset "read DXF file" begin
    groups = read_dxf_file(example_file)
    @test all([g isa DXFGroup for g in groups])
    open("DXF-groups.txt", "w") do io
        DXFutils.show_group_code_registry(io)
    end
end

@testset "parse empty document" begin
    parser = Parser(Vector{DXFGroup}([ EntityType("EOF") ]))
    debugparser(parser; rt=true) do
        parse(parser)
        if (length(parser.pending) == 1
            && parser.pending[1] isa DXFDocument
            && length(parser.pending[1].contents) == 2)
            @test true
        else
            showstate(parser)
            @test false
        end
    end
end

@testset "parser reduce to point" begin
    parser = Parser(Vector{DXFGroup}([
        EntityType("SECTION"),
        Name("HEADER"),
        HeaderVariableName("foo"),
        Point10_X(1.0),
        Point20_Y(2.0),
        Point30_Z(3.0) ]))
    debugparser(parser; rt=true) do
        # @test_throws IncompleteDXFInput parse(parser)
        try
            parse(parser)
            @test false     # Did not throw any exception.
        catch e
            if e isa IncompleteDXFInput
                @test true  # expected error was thrown
            else
                rethrow(e)
            end
        end
        if (isa(parser.pending[end], HeaderVariable) &&
            isa(parser.pending[end].value, DXFPoint))
            @test true
        else
            showstate(parser)
            @test false
        end
    end
end

@testset "parser sections" begin
    parser = Parser(Vector{DXFGroup}([
        EntityType("SECTION"),
        Name("HEADER"),
        EntityType("ENDSEC"),
        EntityType("SECTION"),
        Name("CLASSES"),
        EntityType("ENDSEC"),
        EntityType("SECTION"),
        Name("TABLES"),
        EntityType("ENDSEC"),
        EntityType("SECTION"),
        Name("BLOCKS"),
        EntityType("ENDSEC"),
        EntityType("SECTION"),
        Name("ENTITIES"),
        EntityType("ENDSEC"),
        EntityType("SECTION"),
        Name("OBJECTS"),
        EntityType("ENDSEC"),
        EntityType("EOF") ]))
    debugparser(parser; rt=true) do
        parse(parser)
        if last(parser.pending) isa DXFDocument
            @test true
        else
            showstate(parser)
            @test false
        end
    end
end

@testset "parse sample file" begin
    parser = Parser(example_file)
    debugparser(parser; rt=true) do
        parse(parser)
        @test parser isa Parser
        @test length(parser.pending) == 1
        @test parser.pending[1] isa DXFDocument
    end
end

"""
# Hand testing

begin
  parser = Parser(example_file)
  showstate(parse(parser))
  doc = parser.pending[1]
  for s in doc
    println(summary(s))
  end
  println("\n***** Points")
  pointpaths = dxffind(x -> x isa DXFPoint, doc)
  for pp in pointpaths
    println(reverse([summary(x) for x in pp]))
  end
  println("\n***** Lines")
  linepaths = dxffind(doc) do x
    isa(x, DXFutils.DXFEntity) && x.contents[1].value == "LINE"
  end
  for lp in linepaths
    println(reverse([summary(x) for x in lp]))
  end
end

"""

ext: jl


text:
using SLEEF
using Base.Test

@testset "Basic Tests" begin
	include("basic_tests.jl")
end


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test

@testset "GNL" begin
    rmin = 0.0
    rmax = 2000.0

    kmax1 = 0.0614
    amodes = SFB.AnlmModes(kmax1, rmin, rmax, cache=false)
    cmodes = SFB.ClnnModes(amodes, Δnmax=0)
    knl1 = amodes.knl[isfinite.(amodes.knl)]
    lkk1 = SFB.getlkk(cmodes)
    @test all(@. knl1 < kmax1)
    @test all(@. lkk1[2,:] < kmax1)
    @test all(@. lkk1[3,:] < kmax1)

    kmax2 = 0.1
    amodes = SFB.AnlmModes(kmax2, rmin, rmax, cache=false)
    cmodes = SFB.ClnnModes(amodes, Δnmax=0)
    knl2 = amodes.knl[isfinite.(amodes.knl)]
    lkk2 = SFB.getlkk(cmodes)
    @test all(@. knl2 < kmax2)
    @test all(@. lkk2[2,:] < kmax2)
    @test all(@. lkk2[3,:] < kmax2)

    s = @. knl2 < kmax1
    @show length(knl1) length(knl2[s])
end

# vim: set sw=4 et sts=4 :


text:
module IntCodeMachine
export IntCode, execute!, reset!

struct IntCode
    memory::Vector{Int}
    IntCode(memory::Vector{Int}) = new(copy(memory))
end
IntCode(str::AbstractString) = IntCode(parse.(Int, split(str, ',')))

function Base.rpad(mem::IntCode, n::Integer, p::Int = 0)
    append!(mem.memory, repeat([p], max(0, n - length(mem.memory))))
    mem
end

Base.checkbounds(::Type{Bool}, c::IntCode, i::Integer) = checkbounds(Bool, c.memory, i + 1)

Base.getindex(mem::IntCode, i::Integer) = rpad(mem, i).memory[i+1]
function Base.getindex(mem::IntCode, i::Integer, mode::Mode)
    val = mem[i]
    if mode == Position
        val
    elseif mode == Immediate
        mem[val]
    elseif mode == Relative
        @error "UNFIJNASDOFHSA"
    end
end
Base.getindex(mem::IntCode, r::UnitRange{Int}) = view(mem.memory, r .+ 1)
function Base.getindex(mem::IntCode, r::UnitRange{Int}, ms::Vector{Mode})
    [mem[index, mode] for (index, immediate) in zip(r, ms)]
end

Base.setindex!(mem::IntCode, value, i::Integer) = rpad(mem, i).memory[i+1] = value

reset!(mem::IntCode, str::AbstractString) = reset!(mem, parse.(Int, split(str, ",")))
function reset!(mem::IntCode, xs::Vector{Int})
    append!(empty!(mem.memory), copy(xs))
    mem
end

@enum Mode Position Immediate Relative

function Base.iterate(mem::IntCode, pc::Int = 0)
    if !checkbounds(Bool, mem, pc)
        return nothing
    end

    opcode = mem[pc]
    E, D, A, B, C = digits(opcode, pad = 5)

    # The instruction the opcode represents
    instr = 10D + E

    # The access modes, `0` == position , `1` == immediate
    modes = Mode[A, B, C]
    A, B, C = modes

    if instr in 1:2
        a, b = mem[pc+1:pc+2, modes]
        c = mem[pc+3]
        if instr == 1 # ADD
            @debug "$(pc) $(opcode) ADD: mem[$c] = $a + $b"
            mem[c] = a + b
        elseif instr == 2 # MUL
            @debug "$(pc) $(opcode) MUL: mem[$c] = $a * $b"
            mem[c] = a * b
        end
        (mem, pc + 4)

    elseif instr in 3:4
        if instr == 3 # INPUT
            a = mem[pc+1]
            input = parse(Int, readline())
            @debug "$(pc) $(opcode) INPUT: mem[$a] = $(input)"
            mem[a] = input
        elseif instr == 4 # OUTPUT
            a = mem[pc+1, A]
            @debug "$(pc) $(opcode) OUTPUT: $a"
            println(a)
        end
        (mem, pc + 2)

    elseif instr in 5:6
        a, b = mem[pc+1:pc+2, modes]
        @debug if instr == 5
            "$(pc) $(opcode) JUMP-IF-TRUE: if !iszero($(a)) pc = $(b) end"
        elseif instr == 6
            "$(pc) $(opcode) JUMP-IF-FALSE: if iszero($(a)) pc = $(b) end"
        end
        (mem, (instr == 5 && !iszero(a)) || (instr == 6 && iszero(a)) ? b : pc + 3)

    elseif instr in 7:8
        a, b = mem[pc+1:pc+2, modes]
        mem[mem[pc+3]] = (instr == 7 && a < b) || (instr == 8 && a == b)
        (mem, pc + 4)

    elseif instr == 99
        @debug "$(pc) $(opcode) END"
        nothing

    else
        throw("Unimplemented instruction: $(instr) at $(pc)")
    end
end

Base.IteratorSize(::Type{IntCode}) = Base.SizeUnknown()

function execute!(code::IntCode)
    foreach(code) do _
    end
    code
end

load(filename::AbstractString) = IntCode(readchomp(filename))

end # module

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
directories = Dict(
    "free_convection"          => "2DaySuite/three_layer_constant_fluxes_hr48_Qu0.0e+00_Qb1.2e-07_f1.0e-04_Nh256_Nz128_free_convection_statistics.jld2",
    "strong_wind"              => "2DaySuite/three_layer_constant_fluxes_hr48_Qu1.0e-03_Qb0.0e+00_f1.0e-04_Nh256_Nz128_strong_wind_statistics.jld2",
    "strong_wind_no_coriolis"  => "2DaySuite/three_layer_constant_fluxes_hr48_Qu2.0e-04_Qb0.0e+00_f0.0e+00_Nh256_Nz128_strong_wind_no_rotation_statistics.jld2",
    "weak_wind_strong_cooling" => "2DaySuite/three_layer_constant_fluxes_hr48_Qu3.0e-04_Qb1.0e-07_f1.0e-04_Nh256_Nz128_weak_wind_strong_cooling_statistics.jld2",
    "strong_wind_weak_cooling" => "2DaySuite/three_layer_constant_fluxes_hr48_Qu8.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_strong_wind_weak_cooling_statistics.jld2",
    "strong_wind_weak_heating" => "2DaySuite/three_layer_constant_fluxes_hr48_Qu1.0e-03_Qb-4.0e-08_f1.0e-04_Nh256_Nz128_strong_wind_weak_heating_statistics.jld2",
    "-1e-3"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-9e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu9.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-8e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu8.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-7e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu7.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-6e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu6.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-5e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-4e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-3e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "-2e-4"                          => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2",
    "cooling_6e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb6.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "cooling_5e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb5.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "cooling_4e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb4.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "cooling_3e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb3.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "cooling_2e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb2.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "cooling_1e-8"                   => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb1.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "heating_-3e-8"                  => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2",
    "wind_-5e-4_cooling_4e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb4.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-1e-3_cooling_4e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb4.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-2e-4_cooling_1e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-1e-3_cooling_2e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb2.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-5e-4_cooling_1e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-2e-4_cooling_5e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb5.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-5e-4_cooling_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-2e-4_cooling_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-1e-3_cooling_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-1e-3_heating_-4e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb-4.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-1e-3_heating_-1e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-1e-3_heating_-3e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-5e-4_heating_-5e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-5.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-5e-4_heating_-3e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-5e-4_heating_-1e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-2e-4_heating_-5e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-5.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-2e-4_heating_-3e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-2e-4_heating_-1e-8"       => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",

    #anything above this line is old and should not be used
    "wind_-5e-4_cooling_3e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-5e-4_cooling_2e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-5e-4_cooling_1e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-3.5e-4_cooling_3e-8_new"  => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-3.5e-4_cooling_2e-8_new"  => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-3.5e-4_cooling_1e-8_new"  => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-2e-4_cooling_3e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-2e-4_cooling_2e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-2e-4_cooling_1e-8_new"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_WC_8days_new_statistics.jld2",
    "wind_-5e-4_heating_-3e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-5e-4_heating_-2e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-2.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-5e-4_heating_-1e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-3.5e-4_heating_-3e-8_new" => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-3.5e-4_heating_-2e-8_new" => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb-2.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-3.5e-4_heating_-1e-8_new" => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-2e-4_heating_-3e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-3.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-2e-4_heating_-2e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-2.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",
    "wind_-2e-4_heating_-1e-8_new"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb-1.0e-08_f1.0e-04_Nh256_Nz128_WH_8days_new_statistics.jld2",

    "wind_-5e-4_diurnal_5e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb5.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-5e-4_diurnal_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-5e-4_diurnal_1e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-3.5e-4_diurnal_5e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb5.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3.5e-4_diurnal_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3.5e-4_diurnal_1e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-2e-4_diurnal_5e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb5.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-2e-4_diurnal_3e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-2e-4_diurnal_1e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb1.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-2e-4_diurnal_2e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-2e-4_diurnal_3.5e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb3.5e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3.5e-4_diurnal_2e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3.5e-4_diurnal_3.5e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb3.5e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-5e-4_diurnal_2e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-5e-4_diurnal_3.5e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb3.5e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "cooling_5e-8_new"             => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb5.0e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_4.5e-8_new"           => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb4.5e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_4e-8_new"             => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb4.0e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_3.5e-8_new"           => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb3.5e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_3e-8_new"             => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb3.0e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_2.5e-8_new"           => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb2.5e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_2e-8_new"             => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb2.0e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_1.5e-8_new"           => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb1.5e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",
    "cooling_1e-8_new"             => "Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb1.0e-08_f1.0e-04_Nh256_Nz128_FC_8days_statistics.jld2",

    "wind_-5e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-4.5e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-4e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-3.5e-4_new"              => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.5e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-3e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-2.5e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.5e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",
    "wind_-2e-4_new"                => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",

    "wind_-4.5e-4_cooling_2.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb2.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-2.5e-4_cooling_1.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.5e-04_Qb1.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-4.5e-4_cooling_1.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb1.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-2.5e-4_cooling_2.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.5e-04_Qb2.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",

    "wind_-4.5e-4_heating_-2.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb-2.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-2.5e-4_heating_-1.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.5e-04_Qb-1.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-4.5e-4_heating_-1.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb-1.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-2.5e-4_heating_-2.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.5e-04_Qb-2.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",

    "wind_-4.5e-4_diurnal_4e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb4.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-4.5e-4_diurnal_2e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu4.5e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3e-4_diurnal_4e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.0e-04_Qb4.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-3e-4_diurnal_2e-8"      => "Data/three_layer_constant_fluxes_linear_hr192_Qu3.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-2e-4_diurnal_2e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb2.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-5.5e-4_diurnal_5.5e-8"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.5e-04_Qb5.5e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-1.5e-4_diurnal_5.5e-8"    => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.5e-04_Qb5.5e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",
    "wind_-2e-4_diurnal_4e-8"        => "Data/three_layer_constant_fluxes_linear_hr192_Qu2.0e-04_Qb4.0e-08_f1.0e-04_Nh256_Nz128_diurnal_8days_statistics.jld2",

    "wind_-5.5e-4_new"              => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.5e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_WM_8days_new_statistics.jld2",

    "wind_-5.5e-4_heating_-3.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.5e-04_Qb-3.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",
    "wind_-1.5e-4_heating_-3.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.5e-04_Qb-3.5e-08_f1.0e-04_Nh256_Nz128_WH_8days_statistics.jld2",

    "wind_-5.5e-4_cooling_3.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu5.5e-04_Qb3.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
    "wind_-1.5e-4_cooling_3.5e-8"   => "Data/three_layer_constant_fluxes_linear_hr192_Qu1.5e-04_Qb3.5e-08_f1.0e-04_Nh256_Nz128_WC_8days_statistics.jld2",
)

function diurnal_fluxes(train_files, constants)
    α, g = constants.α, constants.g
    output = Array{Function}(undef,length(train_files))

    @inline wT_flux(Qᵇ, t)::Float32 = Qᵇ * sin(2π / (24 * 60 ^ 2) * t) / (α * g)

    for (i, train_file) in enumerate(train_files)
        if occursin("diurnal_5.5e-8", train_file)
            Q = 5.5f-8
        elseif occursin("diurnal_5e-8", train_file)
            Q = 5f-8
        elseif occursin("diurnal_4e-8", train_file)
            Q = 4f-8
        elseif occursin("diurnal_3.5e-8", train_file)
            Q = 3.5f-8
        elseif occursin("diurnal_3e-8", train_file)
            Q = 3f-8
        elseif occursin("diurnal_2e-8", train_file)
            Q = 2f-8
        elseif occursin("diurnal_1e-8", train_file)
            Q = 1f-8
        end
        output[i] = t -> wT_flux(Q, t)
    end
    return output
end

function read_les_output(filename::String)
    filename = joinpath(pwd(), directories[filename])
    return ReadJLD2_LESbraryData(filename)
end

"""
# Description
Takes NzxNt arrays of profiles for variables u, v, T and returns
Nzx(Nt-1) arrays of the profile evolution for u, v, T, u'w', v'w', and w'T'
the horizontally averaged flux for variable V.

# Arguments
Unscaled u, v, T, z, t, and f
"""
function reconstruct_flux_profiles(u, v, T, νₑ_∂z_u, νₑ_∂z_v, κₑ_∂z_T, zF, t, f)

    Δz = diff(zF)
    Δt = diff(t, dims=1)'

    Nz, Nt = size(T)

    ∂t(A) = (A[:,2:Nt] .- A[:,1:Nt-1]) ./ Δt
    dudt = ∂t(u) # Nz x (Nt-1) array of approximate dUdt values
    dvdt = ∂t(v) # Nz x (Nt-1) array of approximate dVdt values
    dTdt = ∂t(T) # Nz x (Nt-1) array of approximate dTdt values

    ∂z(A) = (A[1:Nz,:] .- A[2:Nz+1,:]) ./ Δz
    νₑ_∂²z_u = ∂z(νₑ_∂z_u)
    νₑ_∂²z_v = ∂z(νₑ_∂z_v)
    κₑ_∂²z_T = ∂z(κₑ_∂z_T)

    # remove the last timestep from the variables that were not differentiated w.r.t t
    u = u[:,1:Nt-1]
    v = v[:,1:Nt-1]
    T = T[:,1:Nt-1]
    νₑ_∂²z_u = νₑ_∂²z_u[:,1:Nt-1]
    νₑ_∂²z_v = νₑ_∂²z_v[:,1:Nt-1]
    κₑ_∂²z_T = κₑ_∂²z_T[:,1:Nt-1]

    """ evaluates wϕ = ∫ ∂z(wϕ) dz """
    function wϕ(∂z_wϕ)
        ans = zeros(Nz+1, Nt-1) # one fewer column than T
        for i in 1:Nt-1, h in 1:Nz-1
            ans[h+1, i] = ans[h, i] + Δz[h] * ∂z_wϕ[h, i]
        end
        return ans
    end

    # duw_dz = -dudt .+ f*v .+ νₑ_∂²z_u
    # dvw_dz = -dvdt .- f*u .+ νₑ_∂²z_v
    # dwT_dz = -dTdt .+ κₑ_∂²z_T

    # Without subgrid fluxes:
    duw_dz = -dudt .+ f*v
    dvw_dz = -dvdt .- f*u
    dwT_dz = -dTdt

    # u, v, T, uw, vw, wT, t
    return (u, v, T, wϕ(duw_dz), wϕ(dvw_dz), wϕ(dwT_dz), t[1:Nt-1])
end

struct FluxData{Z, C, S, U, T} # for each of uw, vw, and wT
                z :: Z # z vector for the variable
           coarse :: C # Nz x Nt array of unscaled profiles
           scaled :: S # Nz x Nt array of scaled profiles
       unscale_fn :: U # function to unscaled profile vectors with
       uvT_scaled :: S # unsubsampled uvT profiles
    training_data :: T # subsampled (uvT, scaled) pairs
end

struct uvTData{Z, C, S, U} # for each of u, v, T
             z :: Z # z vector for the variable
        coarse :: C # Nz x Nt array of unscaled profiles
        scaled :: S # Nz x Nt array of scaled profiles
    unscale_fn :: U # function to unscaled profile vectors with
end

struct ProfileData{Σ, U, V, Θ, UW, VW, WT, T, D}
    grid_points :: Int
   uvT_unscaled :: Σ  # 3Nz x Nt array
     uvT_scaled :: Σ  # 3Nz x Nt array
              u :: U
              v :: V
              T :: Θ
             uw :: UW
             vw :: VW
             wT :: WT
              t :: T  # timeseries Vector
       scalings :: D  # Dict mapping names (e.g. "uw") to the AbstractFeatureScaling object associated with it.
end

"""
    data(filenames; animate=false, scale_type=MinMaxScaling, animate_dir="Output", override_scalings=nothing, reconstruct_fluxes=false)

# Arguments
- filenames                "free_convection"
- animate                  Whether to save an animation of all the original profiles
- animate_dir              Directory to save the animation files to
- scale_type               ZeroMeanUnitVarianceScaling or MinMaxScaling
- override_scalings::Dict  For if you want the testing simulation data to be scaled in the same way as the training data.
                           Set to 𝒟train.scalings to use the scalings from 𝒟train.
"""
function data(filenames; animate=false, scale_type=MinMaxScaling, animate_dir="Output",
                override_scalings=nothing, reconstruct_fluxes=false, subsample_frequency=1,enforce_surface_fluxes=false)

    filenames isa String && (filenames = [filenames])

    # Harvest data from Oceananigans simulation output files.
    all_les = Dict()

    for file in filenames
        all_les[file] = read_les_output(file)
    end

    get_array(f) = cat((f(les) for (file, les) in all_les)..., dims=2)

    u  = get_array(les -> les.U)
    v  = get_array(les -> les.V)
    T  = get_array(les -> les.T)
    # νₑ_∂z_u = get_array(les -> les.νₑ_∂z_u)
    # νₑ_∂z_v = get_array(les -> les.νₑ_∂z_v)
    # κₑ_∂z_T = get_array(les -> les.κₑ_∂z_T)
    t  = cat((les.t for (file, les) in all_les)..., dims=1)

    function enforce_top_surface_flux!(A, flux)
        A[end,:] .= flux
        return A
    end

    vw = get_array(les -> les.wv)
    if enforce_surface_fluxes
        uw = get_array(les -> enforce_top_surface_flux!(les.wu, les.u_top))
        wT = get_array(les -> enforce_top_surface_flux!(les.wT, les.θ_top))
    else
        uw = get_array(les -> les.wu)
        wT = get_array(les -> les.wT)
    end

    first = all_les[filenames[1]]
    zC = first.zC
    zF = first.zF

    NzC = length(zC)
    NzF = length(zF)
    Nt = sum([length(les.t) for (file, les) in all_les])

    if reconstruct_fluxes
        Nt -= length(all_les)
        u = zeros(NzC, Nt)
        v = zeros(NzC, Nt)
        T = zeros(NzC, Nt)
        uw = zeros(NzF, Nt)
        vw = zeros(NzF, Nt)
        wT = zeros(NzF, Nt)
        t = zeros(Nt)
        x=1
        for (file,les) in all_les
            Nt = length(les.t)-1
            cols = x:x+Nt-1 # time indices corresponding to the current simulation
            # @. u[:,cols], v[:,cols], T[:,cols], uw[:,cols], vw[:,cols], wT[:,cols], t[cols] = reconstruct_flux_profiles(les.U, les.V, les.T, les.νₑ_∂z_u, les.νₑ_∂z_v, les.κₑ_∂z_T, zF, les.t, les.f⁰)
            a = reconstruct_flux_profiles(les.U, les.V, les.T, les.νₑ_∂z_u, les.νₑ_∂z_v, les.κₑ_∂z_T, zF, les.t, les.f⁰)
            u[:,cols]  .= a[1]
            v[:,cols]  .= a[2]
            T[:,cols]  .= a[3]
            uw[:,cols] .= a[4]
            vw[:,cols] .= a[5]
            wT[:,cols] .= a[6]
            t[cols]    .= a[7]
            if enforce_surface_fluxes
                uw[end,cols] .= les.u_top
                wT[end,cols] .= les.θ_top
            end
            x += Nt
        end
    end

    if animate
        animate_gif([uw], zF, t, "uw", directory=animate_dir)
        animate_gif([vw], zF, t, "vw", directory=animate_dir)
        animate_gif([wT], zF, t, "wT", directory=animate_dir)
        animate_gif([u],  zC, t, "u",  directory=animate_dir)
        animate_gif([v],  zC, t, "v",  directory=animate_dir)
        animate_gif([T],  zC, t, "T",  directory=animate_dir)
    end

    function coarsify_cell(x)
        output = zeros(typeof(x[1]), 32, size(x,2))
        for i in 1:size(x,2)
            col = @view output[:,i]
            col .= coarse_grain(x[:,i], 32, Center)
        end
        return output
    end

    function coarsify_face(x)
        output = zeros(typeof(x[1]), 33, size(x,2))
        for i in 1:size(x,2)
            col = @view output[:,i]
            col .= coarse_grain_linear_interpolation(x[:,i], 33, Face)
            # col .= coarse_grain(x[:,i], 33, Face)
        end
        return output
    end
    # coarsify_cell(x) = cat((coarse_grain(x[:,i], 32, Center) for i in 1:size(x,2))..., dims=2)
    # coarsify_face(x) = cat((coarse_grain_linear_interpolation(x[:,i], 33, Face) for i in 1:size(x,2))..., dims=2)
    # coarsify_face(x) = cat((coarse_grain(x[:,i], 33, Face) for i in 1:size(x,2))..., dims=2)

    u_coarse  = coarsify_cell(u)
    v_coarse  = coarsify_cell(v)
    T_coarse  = coarsify_cell(T)
    uw_coarse = coarsify_face(uw)
    vw_coarse = coarsify_face(vw)
    wT_coarse = coarsify_face(wT)
    # νₑ_∂z_u   = coarsify_face(νₑ_∂z_u)
    # νₑ_∂z_v   = coarsify_face(νₑ_∂z_v)
    # κₑ_∂z_T   = coarsify_face(κₑ_∂z_T)

    zC_coarse = coarse_grain(zC, 32, Center)
    zF_coarse = coarse_grain_linear_interpolation(zF, 33, Face)
    # zF_coarse = coarse_grain(zF, 33, Face)

    function get_scaling(name, coarse)
        if isnothing(override_scalings)
            # set the scaling according to the data (for training simulations)
            return scale_type(coarse)
        else
            # for if you want the testing simulation data to be scaled in the same way as the training data
            return override_scalings[name]
        end
    end

    all_scalings=Dict()

    for (name, coarse) in [("u", u_coarse), ("v", v_coarse), ("T", T_coarse),
                           ("uw", uw_coarse), ("vw", vw_coarse), ("wT", wT_coarse)]
        all_scalings[name] = get_scaling(name, coarse)
    end

    training_set = collect(1:subsample_frequency:length(t))

    function get_scaled(name, coarse)
        scaling = get_scaling(name, coarse)
        scaled = all_scalings[name].(coarse)
        unscale_fn = Base.inv(scaling)
        return (scaled, unscale_fn)
    end

    function get_uvTData(name, coarse, z)
        scaled, unscale_fn = get_scaled(name, coarse)
        return uvTData(z, coarse, scaled, unscale_fn)
    end

    function get_FluxData(name, coarse, z)
        scaled, unscale_fn = get_scaled(name, coarse)
        training_data = [(uvT_scaled[:,i], scaled[:,i]) for i in training_set] # (predictor, target) pairs
        return FluxData(z, coarse, scaled, unscale_fn, uvT_scaled, training_data)
    end

    u = get_uvTData("u", u_coarse, zC_coarse)
    v = get_uvTData("v", v_coarse, zC_coarse)
    T = get_uvTData("T", T_coarse, zC_coarse)
    uvT_unscaled = cat(u.coarse, v.coarse, T.coarse, dims=1)
    uvT_scaled = cat(u.scaled, v.scaled, T.scaled, dims=1)

    uw = get_FluxData("uw", uw_coarse, zF_coarse)
    vw = get_FluxData("vw", vw_coarse, zF_coarse)
    wT = get_FluxData("wT", wT_coarse, zF_coarse)

    return ProfileData(33, uvT_unscaled, uvT_scaled, u, v, T, uw, vw, wT, t, all_scalings)
end


text:
"""
Genie Web Framework を使ったサーバ
"""

using Genie
import Genie.Router: route
import Genie.Renderer.Json: json

# route: / => {message: "Hello, Genie!"}
route("/") do
    (:message => "Hello, Genie!") |> json
end

# Genie設定: http://localhost:8080
Genie.config.run_as_server = true
Genie.config.server_port = 8080
Genie.config.server_host = "0.0.0.0"
# Nginxのように静的ファイルを配信: ./*
Genie.config.server_handle_static_files = true
Genie.config.server_document_root = "./"

# Genie起動
Genie.startup()

max_stars_repo_name: amenoyoya/wsgi-api-server

max_issues_repo_name: amenoyoya/wsgi-api-server

max_forks_repo_name: amenoyoya/wsgi-api-server


text:
# Randomized nonnegative matrix factorization

function rnmf(X::Array{T,2},
              k::Int,
              p::Int=5,
              q::Int=1,
              max_iter::Int=1000,
              tol::Float64=1e-5)
              where T <: Number

    m, n = size(X)
    P = rand(n, k + p)
    Y = X * P
    sub_iterations!(X, Y, q)
    Q = qr(Y).Q
    B = Q' * X

    W = rand(m, k)
    Wr = rand(k + p, k)
    H = rand(k, n)
    res = 1

    for i = 1:max_iter
        HALS!()
        if hasConverged()
            break
    end

    return W, H
end

function HALS!(W, Wr, H, B, k::Int)
    R = B' * Wr
    S = W'W

    updateH!(H, R, S, k)

    T = B * H'
    V = H * H'

    updateWr!(Wr, Q, T, V, Q, k)
end

function updateH!(H, R, S, k)
    for j = 1:k
        H[j, :] += R[:, j] - H'S[:, j] ./ S[j, j]
        H[j, H[j, :] < 0] = 0
    end
end

function updateWr!(Wr, Q, T, V, Q, k)
    for j = 1:k
        Wr[]
    end
end

function hasConverged()
end


text:
# getindex/view/setindex! ======================================================

#### getindex/view ####

# Integer returns a single value, but not for view
@propagate_inbounds Base.getindex(A::AbstractDimArray, i1::Integer, i2::Integer, I::Integer...) =
    Base.getindex(parent(A), i1, i2, I...)
# No indices. These just prevent stack overflows
@propagate_inbounds Base.getindex(A::AbstractDimArray) = Base.getindex(parent(A))
@propagate_inbounds Base.view(A::AbstractDimArray) = rebuild(A, Base.view(parent(A)), ())

const SelectorOrStandard = Union{SelectorOrInterval,StandardIndices}

for f in (:getindex, :view, :dotview)
    @eval begin
        #### Array getindex/view ###
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::Integer) = Base.$f(parent(A), i)
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::CartesianIndex) = Base.$f(parent(A), i)
        # Linear indexing forwards to the parent array
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::Union{Colon,AbstractVector{<:Integer}}) =
            Base.$f(parent(A), i)
        @propagate_inbounds Base.$f(A::AbstractDimArray, i::AbstractArray{<:Bool}) =
            Base.$f(parent(A), i)
        # Except 1D DimArrays
        @propagate_inbounds Base.$f(A::AbstractDimArray{<:Any,1}, i::Union{Colon,AbstractVector{<:Integer}}) =
            rebuildsliced(Base.$f, A, Base.$f(parent(A), i), (i,))
        # Selector/Interval indexing
        @propagate_inbounds Base.$f(A::AbstractDimArray, i1::SelectorOrStandard, I::SelectorOrStandard...) =
            Base.$f(A, dims2indices(A, (i1, I...))...)
        @propagate_inbounds Base.$f(A::AbstractDimArray, extent::Extents.Extent) =
            Base.$f(A, dims2indices(A, extent)...)
        # Dimension indexing. Allows indexing with A[somedim=At(25.0)] for Dim{:somedim}
        @propagate_inbounds Base.$f(A::AbstractDimArray, args::Dimension...; kw...) =
            Base.$f(A, dims2indices(A, (args..., kwdims(values(kw))...))...)
        # Everything else - such as custom indexing types from other packages
        # @propagate_inbounds Base.$f(A::AbstractDimArray, i1, I...) = Base.$f(parent(A), i1, I...)
    end
    # Standard indices
    if f == :view
        @eval @propagate_inbounds function Base.$f(A::AbstractDimArray, i1::StandardIndices, i2::StandardIndices, I::StandardIndices...)
            I = _unwrap_cartesian(i1, i2, I...)
            x = Base.$f(parent(A), I...)
            rebuildsliced(Base.$f, A, x, I)
        end
    else
        @eval @propagate_inbounds function Base.$f(A::AbstractDimArray, i1::StandardIndices, i2::StandardIndices, I::StandardIndices...)
            I = _unwrap_cartesian(i1, i2, I...)
            x = Base.$f(parent(A), I...)
            all(i -> i isa Integer, I) ? x : rebuildsliced(Base.$f, A, x, I)
        end
    end
end

@inline _unwrap_cartesian(i1::CartesianIndex, I...) = (Tuple(i1)..., _unwrap_cartesian(I...)...)
@inline _unwrap_cartesian(i1, I...) = (i1, _unwrap_cartesian(I...)...)
@inline _unwrap_cartesian() = ()

#### setindex ####

@propagate_inbounds Base.setindex!(A::AbstractDimArray, x) = setindex!(parent(A), x)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, args::Dimension...; kw...) =
    setindex!(A, x, dims2indices(A, (args..., kwdims(values(kw))...))...)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, i, I...) =
    setindex!(A, x, dims2indices(A, (i, I...))...)
@propagate_inbounds Base.setindex!(A::AbstractDimArray, x, i1::StandardIndices, I::StandardIndices...) =
    setindex!(parent(A), x, i1, I...)


text:
# Note: These methods need a more helpfull error message than a `NoMethodError`,
#       when the deprecation is removed
import Base.@deprecate
@deprecate (+)(A::AbstractDataArray{Bool},x::Bool)              A .+ x
@deprecate (+)(x::Bool,A::AbstractDataArray{Bool})              x .+ A
@deprecate (-)(A::AbstractDataArray{Bool},x::Bool)              A .- x
@deprecate (-)(x::Bool,A::AbstractDataArray{Bool})              x .- A
@deprecate (+)(A::AbstractDataArray,x::Union(NAtype,Number))    A .+ x
@deprecate (+)(x::Union(NAtype,Number),A::AbstractDataArray)    x .+ A
@deprecate (-)(A::AbstractDataArray,x::Union(NAtype,Number))    A .- x
@deprecate (-)(x::Union(NAtype,Number),A::AbstractDataArray)    x .- A
@deprecate (/)(x::Union(NAtype,Number),A::AbstractDataArray)    x ./ A


text:
function romberg(f::Function, a::Number, b::Number, nl::Int)
  #= Numerical Integration of functions using Romberg's method.
  Keyword arguments:
  f -- Anonymous functon to integrate
  a -- Lower limit of integration
  b -- Upper limit of integration
  nl -- The level of integration to go to.

  Returns:
  The result of the integration of f between a & b
  =#

  intTable = zeros(nl, nl); # Empty matrix to use as div diff table

  # Calculate the first column of the div diff table using the trapezium rule.
  h = (b - a);
  for k = 1:nl
    intTable[k,1] = 0.5 * h * (f(a) + 2*sum(f([a+h:h:b-h])) + f(b));
    h /= 2;
  end

  # Calculate the remaining elements using Romberg's method.
  for k = 2:nl, j = nl-k+1:-1:1
    intTable[j, k] = ((4^(k-1))intTable[j+1, k-1] - intTable[j, k-1])/((4^(k-1)) - 1);
  end

  return intTable[1, end];
end

function romberg(intEstimates::Array{Float64,1})
  #= Estiamte integrals from sets of data using the Romberg Method.

  Keyword arguments:
  intEstimates -- An array of initial estimates of integrals

  Returns:
    A tuple containing the integral and the relative error in the estimate
  =#
  intTable = zeros(length(intEstimates), length(intEstimates));
  relativeErrorTable = zeros(size(intEstimates));

  intTable[:, 1] = intEstimates;

  for k = 2:size(intTable)[1]
    for j = size(intTable)[1]-k+1:-1:1
      intTable[j, k] = ((4^(k-1))intTable[j+1, k-1] - intTable[j, k-1])/((4^(k-1)) - 1);
    end
  end

  for k = 2:size(intTable)[1]
    relativeErrorTable[k, 1] = 100 * abs((intTable[1, k] - intTable[2, k-1]) / intTable[1, k]);
  end

  return (intTable[1, end], relativeErrorTable[1, end]);
end

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_name: alexjohnj/numerical-methods


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames

ui = Builder(filename=(@__DIR__) * "/ui.glade")
showall(ui["win"])

gettext(textview::GtkTextView) = get_gtk_property(textview, :buffer, GtkTextBuffer) |> x -> get_gtk_property(x, :text, AbstractString)

function eqn2val()
	s = replace(gettext(ui["eqn"]), "\n" => ";")
	val = eval(Meta.parse(s))

	set_gtk_property!(ui["resultent"], :text, val)

	return nothing
end

function reset_button()
	set_gtk_property!(ui["resultent"], :text, "")

	return nothing
end

signal_connect(x -> eqn2val(), ui["cal"], "clicked")
signal_connect(x -> reset_button(), ui["reset"], "clicked")

if !isinteractive()
	c = Condition()
	signal_connect(ui["win"], :destroy) do widget
		notify(c)
	end
	wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_tso_limitation" begin

    @testset verbose=true "test_tso_limitation_after_dp" begin

        #=
        TS: [11h]
        S: [S1]
                            bus 1
                            |
        (limitable) wind_1_1|load
        Pmin=0, Pmax=100    | S1: 30
        Csta=0, Cprop=1     | S2: 35
        DP => 9h30          |
        S1: 55              |
        S2: 60              |
        =#
        @testset "tso_limits_limitables" begin
            TS = [DateTime("2015-01-01T11:00:00")]
            ech = DateTime("2015-01-01T10:00:00")
            network = PSCOPF.Networks.Network()
            # Buses
            PSCOPF.Networks.add_new_bus!(network, "bus_1")
            # Limitables
            PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                    0., 100.,
                                                    0., 1.,
                                                    Dates.Second(90*60), Dates.Second(90*60))
            # Uncertainties
            uncertainties = PSCOPF.Uncertainties()
            # initial generators state : need to pay starting cost at TS[1]
            generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()
            #ManagementMode
            mode = PSCOPF.ManagementMode("test_mode", Dates.Minute(5))

            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 30.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 35.)

            # firmness
            firmness = PSCOPF.Firmness(
                    SortedDict(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                    )

            context = PSCOPF.PSCOPFContext(network, TS, mode,
                                            generators_init_state,
                                            uncertainties, nothing)
            tso = PSCOPF.TSOOutFO()

            @test firmness == PSCOPF.compute_firmness(tso,
                                                    ech, DateTime("2015-01-01T10:40:00"),
                                                    TS, context)

            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
            # Limit
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S1"])
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S2"])

            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"])
            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"])

            # active limits ?
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"]) # now 0, should be 1 ?
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"]) # now 0, should be 35 ?
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"]) # now 0, should be 1 ?
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"]) # now 0, should be 35 ?
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
            @test 5. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"])
        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1
                            |
        (limitable) wind_1_1|load
        Pmin=0, Pmax=100    | S1: 30
        Csta=0, Cprop=1     | S2: 60
        S1: 55            |
        S2: 60            |
        =#
        @testset "tso_cannot_choose_injection_level_1" begin
            TS = [DateTime("2015-01-01T11:00:00")]
            ech = DateTime("2015-01-01T10:00:00")
            network = PSCOPF.Networks.Network()
            # Buses
            PSCOPF.Networks.add_new_bus!(network, "bus_1")
            # Limitables
            PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                    0., 100.,
                                                    0., 1.,
                                                    Dates.Second(90*60), Dates.Second(90*60))
            # Uncertainties
            uncertainties = PSCOPF.Uncertainties()
            # initial generators state : need to pay starting cost at TS[1]
            generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()
            #ManagementMode
            mode = PSCOPF.ManagementMode("test_mode", Dates.Minute(5))

            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 30.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

            # firmness
            firmness = PSCOPF.Firmness(
                    SortedDict(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                    )

            context = PSCOPF.PSCOPFContext(network, TS, mode,
                                            generators_init_state,
                                            uncertainties, nothing)
            tso = PSCOPF.TSOOutFO()

            @test firmness == PSCOPF.compute_firmness(tso,
                                                    ech, DateTime("2015-01-01T10:40:00"),
                                                    TS, context)

            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
            # Limit
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S1"])
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S2"])

            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"])
            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"])

            # active limits ?
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"])
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"])
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"])
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"])
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
            @test 30. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"])
        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1
                            |
        (limitable) wind_1_1|load
        Pmin=0, Pmax=100    | S1: 30
        Csta=0, Cprop=1     | S2: 60
        S1: 70            |
        S2: 60            |
        uncertainty(S1) > limit
        =#
        @testset "tso_cannot_choose_injection_level_2" begin
            TS = [DateTime("2015-01-01T11:00:00")]
            ech = DateTime("2015-01-01T10:00:00")
            network = PSCOPF.Networks.Network()
            # Buses
            PSCOPF.Networks.add_new_bus!(network, "bus_1")
            # Limitables
            PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                    0., 100.,
                                                    0., 1.,
                                                    Dates.Second(90*60), Dates.Second(90*60))
            # Uncertainties
            uncertainties = PSCOPF.Uncertainties()
            # initial generators state : need to pay starting cost at TS[1]
            generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()
            #ManagementMode
            mode = PSCOPF.ManagementMode("test_mode", Dates.Minute(5))

            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 70.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 30.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

            # firmness
            firmness = PSCOPF.Firmness(
                    SortedDict(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                    )

            context = PSCOPF.PSCOPFContext(network, TS, mode,
                                            generators_init_state,
                                            uncertainties, nothing)
            tso = PSCOPF.TSOOutFO()

            @test firmness == PSCOPF.compute_firmness(tso,
                                                    ech, DateTime("2015-01-01T10:40:00"),
                                                    TS, context)

            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
            # Limit
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S1"])
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S2"])

            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"])
            @test 30. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"])

            # active limits ?
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"])
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"])
            @test 1 ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"])
            @test 30. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"])
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
            @test 30. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"])
        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1
                            |
        (limitable) wind_1_1|load
        Pmin=0, Pmax=100    | S1: 30
        Csta=0, Cprop=1     | S2: 60
        S1: 30            |
        S2: 60            |
        =#
        @testset "tso_no_need_for_limit" begin
            TS = [DateTime("2015-01-01T11:00:00")]
            ech = DateTime("2015-01-01T10:00:00")
            network = PSCOPF.Networks.Network()
            # Buses
            PSCOPF.Networks.add_new_bus!(network, "bus_1")
            # Limitables
            PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                    0., 100.,
                                                    0., 1.,
                                                    Dates.Second(90*60), Dates.Second(90*60))
            # Uncertainties
            uncertainties = PSCOPF.Uncertainties()
            # initial generators state : need to pay starting cost at TS[1]
            generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()
            #ManagementMode
            mode = PSCOPF.ManagementMode("test_mode", Dates.Minute(5))

            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 30.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 30.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)
            PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

            # firmness
            firmness = PSCOPF.Firmness(
                    SortedDict(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                    )

            context = PSCOPF.PSCOPFContext(network, TS, mode,
                                            generators_init_state,
                                            uncertainties, nothing)
            tso = PSCOPF.TSOOutFO()

            @test firmness == PSCOPF.compute_firmness(tso,
                                                    ech, DateTime("2015-01-01T10:40:00"),
                                                    TS, context)

            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
            # Limit
            @test 30. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S1"])
            @test 60. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1],"S2"])

            @test ( 60. -1e-09 < value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]) < 100 +1e-09 ) # betwen 60 and 100
            @test ( 60. -1e-09 < value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) < 100 +1e-09 ) # betwen 60 and 100

            # active limits ?
            @test value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
            @test value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
            @test value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"]) < 1e-09
            @test value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"]) < 1e-09
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"]) < 1e-09
        end

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
Base.@kwdef struct Case
    # We create a struct to hold case data in a type-declared format
    # `Base.@kwdef` allows us to instantiate this via keywords

    branchid::Array{Int64,1}
    branch_from::Array{Int64,1}
    branch_to::Array{Int64,1}
    branch_reactance::Array{Float64,1}
    branch_rating::Array{Float64,1}

    dclineid::Array{Int64,1}
    dcline_from::Array{Int64,1}
    dcline_to::Array{Int64,1}
    dcline_pmin::Array{Float64,1}
    dcline_pmax::Array{Float64,1}

    busid::Array{Int64,1}
    bus_demand::Array{Float64,1}
    bus_zone::Array{Int64,1}

    genid::Array{Int64,1}
    genfuel::Array{String,1}
    gen_bus::Array{Int64,1}
    gen_status::BitArray{1}
    gen_pmax::Array{Float64,1}
    gen_pmin::Array{Float64,1}
    gen_ramp30::Array{Float64,1}

    gencost::Array{Float64,2}
    gencost_orig::Array{Float64,2}

    demand::DataFrames.DataFrame
    hydro::DataFrames.DataFrame
    wind::DataFrames.DataFrame
    solar::DataFrames.DataFrame
end

Base.@kwdef struct Storage
    enabled::Bool
    gen::Array{Float64,2}
    sd_table::DataFrames.DataFrame
end

Base.@kwdef struct DemandFlexibility
    flex_amt_up::Union{DataFrames.DataFrame,Nothing}
    flex_amt_dn::Union{DataFrames.DataFrame,Nothing}
    cost_dn::Union{DataFrames.DataFrame,Nothing}
    cost_up::Union{DataFrames.DataFrame,Nothing}
    duration::Int64
    enabled::Bool
    interval_balance::Bool
    rolling_balance::Bool
    input_granularity::String
end

Base.@kwdef struct Results
    # We create a struct to hold case results in a type-declared format
    pg::Array{Float64,2}
    pf::Array{Float64,2}
    lmp::Array{Float64,2}
    congu::Array{Float64,2}
    congl::Array{Float64,2}
    pf_dcline::Array{Float64,2}
    storage_pg::Array{Float64,2}
    storage_e::Array{Float64,2}
    load_shed::Array{Float64,2}
    load_shift_up::Array{Float64,2}
    load_shift_dn::Array{Float64,2}
    f::Float64
    status::String
end

Base.@kwdef struct Sets
    # Branch and branch subsets
    num_branch::Int64
    num_branch_ac::Int64
    branch_idx::UnitRange{Int64}
    noninf_branch_idx::Array{Int64,1}
    branch_to_idx::Array{Int64,1}
    branch_from_idx::Array{Int64,1}
    # Bus & bus subsets
    num_bus::Int64
    num_load_bus::Int64
    bus_idx::UnitRange{Int64}
    load_bus_idx::Array{Int64,1}
    bus_id2idx::Dict{Int64,Int64}
    load_bus_map::SparseMatrixCSC{Int64,Int64}
    # Gen & gen sub-sets
    num_gen::Int64
    num_wind::Int64
    num_solar::Int64
    num_hydro::Int64
    gen_idx::UnitRange{Int64}
    gen_wind_idx::Array{Int64,1}
    gen_solar_idx::Array{Int64,1}
    gen_hydro_idx::Array{Int64,1}
    renewable_idx::Array{Int64,1}
    noninf_pmax::Array{Int64,1}
    # Segments
    num_segments::Int64
    segment_idx::UnitRange{Int64}
    ## Demand Flexibility
    flexible_bus_idx::Union{Array{Int64,1},Nothing}
    num_flexible_bus::Int64
    flexible_load_bus_map::Union{SparseMatrixCSC,Nothing}
    # Storage
    num_storage::Int64
    storage_idx::Union{UnitRange{Int64},Nothing}
end


text:
# Inverse Gamma distribution

immutable InverseGamma <: ContinuousUnivariateDistribution
    shape::Float64 # location
    scale::Float64 # scale
    function InverseGamma(sh::Real, sc::Real)
    	sh > zero(sh) && sc > zero(sc) || error("Both shape and scale must be positive")
    	new(float64(sh), float64(sc))
    end
end

@distr_support InverseGamma 0.0 Inf

_inv(d::InverseGamma) = Gamma(d.shape, 1.0 / d.scale)

scale(d::InverseGamma) = d.scale
rate(d::InverseGamma) = 1.0 / d.scale

mean(d::InverseGamma) = d.shape > 1.0 ? d.scale / (d.shape - 1.0) : Inf
var(d::InverseGamma) = d.shape > 2.0 ? abs2(d.scale) / (abs2(d.shape - 1.0) * (d.shape - 2.0)) : Inf
skewness(d::InverseGamma) = d.shape > 3.0 ? (4.0 * sqrt(d.shape - 2.0)) / (d.shape - 3.0) : NaN
kurtosis(d::InverseGamma) = d.shape > 4.0 ? (30.0 * d.shape - 66.0) / ((d.shape - 3.0) * (d.shape - 4.0)) : NaN

mode(d::InverseGamma) = d.scale / (d.shape + 1.0)

function entropy(d::InverseGamma)
    a = d.shape
    b = d.scale
    a + log(b) + lgamma(a) - (1.0 + a) * digamma(a)
end

pdf(d::InverseGamma, x::Float64) = exp(logpdf(d, x))

function logpdf(d::InverseGamma, x::Float64)
    a = d.shape
    b = d.scale
    a * log(b) - lgamma(a) - (a + 1.0) * log(x) - b / x
end

cdf(d::InverseGamma, x::Float64) = ccdf(_inv(d), 1.0 / x)
ccdf(d::InverseGamma, x::Float64) = cdf(_inv(d), 1.0 / x)
logcdf(d::InverseGamma, x::Float64) = logccdf(_inv(d), 1.0 / x)
logccdf(d::InverseGamma, x::Float64) = logcdf(_inv(d), 1.0 / x)

quantile(d::InverseGamma, p::Float64) = 1.0 / cquantile(_inv(d), p)
cquantile(d::InverseGamma, p::Float64) = 1.0 / quantile(_inv(d), p)
invlogcdf(d::InverseGamma, p::Float64) = 1.0 / invlogccdf(_inv(d), p)
invlogccdf(d::InverseGamma, p::Float64) = 1.0 / invlogcdf(_inv(d), p)

function mgf(d::InverseGamma, t::Real)
    a = d.shape
    b = d.scale
    t == zero(t) ? one(float(t)) : 2.0*(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*b*t))
end

function cf(d::InverseGamma, t::Real)
    a = d.shape
    b = d.scale
    t == zero(t) ? complex(one(float(t))) : 2.0*(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*im*b*t))
end

rand(d::InverseGamma) = 1.0 / rand(_inv(d))

function rand!(d::InverseGamma, A::Array{Float64})
    rand!(_inv(d), A)
    for i in 1:length(A)
    	A[i] = 1.0 / A[i]
    end
    A
end


text:
resolutions(::Type{EarthEnv{HabitatHeterogeneity}}) = ("1km", "5km", "25km")
defres(::Type{EarthEnv{HabitatHeterogeneity}}) = "25km"
layers(::Type{EarthEnv{HabitatHeterogeneity}}) = values(heterogeneity_lookup)
layerkeys(::Type{EarthEnv{HabitatHeterogeneity}}) = keys(heterogeneity_lookup)
layerkeys(T::Type{EarthEnv{HabitatHeterogeneity}}, layers) = map(l -> layerkeys(T, l), layers)
function layerkeys(::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol)
    Symbol(lowercase(string(layer)))
end

# Yes, they randomly chose cases
const heterogeneity_lookup = (
    cv = :cv,
    evenness = :evenness,
    range = :range,
    shannon = :shannon,
    simpson = :simpson,
    std = :std,
    contrast = :Contrast,
    correlation = :Correlation,
    dissimilarity = :Dissimilarity,
    entropy = :Entropy,
    homogeneity = :Homogeneity,
    maximum = :Maximum,
    uniformity = :Uniformity,
    variance = :Variance,
)

heterogeneity_layer(x::Symbol) = heterogeneity_lookup[Symbol(lowercase(string(x)))]

"""
    getraster(source::Type{EarthEnv{HabitatHeterogeneity}}, [layer]; res="25km")

Download [`EarthEnv`](@ref) habitat heterogeneity data.

# Arguments
- `layer`: `Symbol` or `Tuple` of `Symbol` from `$(layers(EarthEnv{HabitatHeterogeneity}))`.
    Without a `layer` argument, all layers will be downloaded, and a `NamedTuple` of paths returned.

# Keywords
- `res`: `String` chosen from `$(resolutions(EarthEnv{HabitatHeterogeneity}))`, defaulting to "25km".

Returns the filepath/s of the downloaded or pre-existing files.
"""
function getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layers::Union{Tuple,Symbol};
    res::String=defres(T)
)
    _getraster(T, layers, res)
end

function _getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layers::Tuple, res::String)
    return _map_layers(T, layers, res)
end
function _getraster(T::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol, res::String)
    layer = heterogeneity_layer(layer)
    _check_layer(T, layer)
    _check_res(T, res)
    path = rasterpath(T, layer; res)
    url = rasterurl(T, layer; res)
    return _maybe_download(url, path)
end

function rastername(::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res::String=defres(T))
    "$(layer)_$(res).tif"
end
function rasterpath(::Type{EarthEnv{HabitatHeterogeneity}})
    joinpath(rasterpath(EarthEnv), "HabitatHeterogeneity")
end
function rasterpath(T::Type{<:EarthEnv{HabitatHeterogeneity}}, layer::Symbol; res::String=defres(T))
    joinpath(rasterpath(T), string(res), rastername(T, layer; res))
end
function rasterurl(::Type{EarthEnv{HabitatHeterogeneity}})
    joinpath(rasterurl(EarthEnv), "habitat_heterogeneity")
end
function rasterurl(T::Type{EarthEnv{HabitatHeterogeneity}}, layer; res::String=defres(T))
    prec = _getprecision(layer, res)
    layerpath = "$res/$(layer)_01_05_$(res)_$prec.tif"
    joinpath(rasterurl(T), layerpath)
end

# See http://www.earthenv.org/texture
function _getprecision(layer, res)
    if ((res in ("1km", "5km")) && (layer == :Correlation))
        "int16"
    elseif ((res == "5km") && (layer == :cv)) ||
           ((res == "25km") && (layer == :Entropy)) ||
           layer in (:Contrast, :Dissimilarity, :Variance)
        "uint32"
    else
        "uint16"
    end
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
This is the restart version of Coordinate Linear Variance Reduction (Lazy Update Version).
"""
function clvr_lazy_restart_x_y(
    problem::StandardLinearProgram,
    exitcriterion::ExitCriterion;
    γ=1.0, σ=0.0, R=10, blocksize=10, restartfreq=Inf,
    io=nothing)

    @info("Running clvr_lazy_restart_x_y with")
    @info("blocksize = $(blocksize)")
    @info("γ = $(γ)")
    @info("σ = $(σ)")
    @info("R = $(R)")
    @info("restartfreq = $(restartfreq)")
    if !isnothing(io)
        flush(io)
    end

    # Algorithm 2 from the paper

    A_T, b, c = problem.A_T, problem.b, problem.c
    prox = problem.prox

    d, n = size(A_T)
    x0 = zeros(d)
    y0 = zeros(n)

    _time1 = time()
    # Precomputing blocks, nzrows, sliced_A_T
    blocks, C = compute_nzrows_for_blocks(A_T, blocksize)
    sliced_A_Ts = Array{SparseMatrixCSC{Float64, Int}}([])
    for j in 1:length(C)
        push!(sliced_A_Ts, A_T[C[j], blocks[j]])
    end
    _time2 = time()
    @info ("Initialization time = ", _time2 - _time1)

    ##### Start of clvr_lazy_restart_x_y

    m = length(blocks)

    # Log initial measure
    starttime = time()
    results = Results()
    init_fvaluegap, init_metricLP = compute_fvaluegap_metricLP(x0, y0, problem)
    logresult!(results, 1, 0.0, init_fvaluegap, init_metricLP)

    outer_k = 1
    exitflag = false

    while !exitflag
        # Init of CLVR Lazy
        a = 1 / (R * m)
        pre_a = a
        idx_seq = 1:m

        x = deepcopy(x0)
        y = deepcopy(y0)
        x_tilde = zero(x0)
        y_tilde = zero(y0)
        z = A_T * y + c
        q = a * deepcopy(z)
        θ_x = ones(Int, length(x0))
        θ_y = ones(Int, length(y0))

        k = 1
        restartflag = false
        while !exitflag && !restartflag
            # Line 4
            j = rand(idx_seq)

            # Slice of variables based on nzrowsC[j]
            z_sliced = z[C[j]]
            q_sliced = q[C[j]]
            Adelta_sliced = a * ((k-1) .- θ_x[C[j]])

            # Line 5
            q_hat = q_sliced + Adelta_sliced .* z_sliced

            # Line 6
            x_hat = prox(x0[C[j]] - 1/γ * q_hat, 1/γ * a * k)

            # Line 7 & 12
            sliced_A_T = sliced_A_Ts[j]
            Delta_y = γ * m * a * ((x_hat' * sliced_A_T)' - b[blocks[j]])

            y_tilde[blocks[j]] = y_tilde[blocks[j]] + a * (k .- θ_y[blocks[j]]) .* y[blocks[j]]+ (m-1) * a * Delta_y[:]

            y[blocks[j]] = y[blocks[j]] + Delta_y

            # Line 10
            Delta_Delta_y = sliced_A_T * Delta_y
            q[C[j]] = q_sliced + a * ((k+1) .- θ_x[C[j]]) .* z_sliced + (m+1) * a * Delta_Delta_y

            # Line 11
            x_tilde[C[j]] = x_tilde[C[j]] + a * (k .- θ_x[C[j]]) .* x[C[j]]

            # Line 9
            z[C[j]] = z[C[j]] + Delta_Delta_y

            # Line 13
            x[C[j]] =  prox(x0[C[j]] - 1/γ * q[C[j]], 1/γ * (k+1) * a)

            # Line 14
            θ_x[C[j]] .= k
            θ_y[blocks[j]] .= k

            # Logging and checking exit condition
            # set restartflag when reached some measure
            if outer_k % (exitcriterion.loggingfreq * m) == 0
                x_tilde_tmp = x_tilde[:] + a * ((k+1) .- θ_x[:]) .* x[:]
                y_tilde_tmp = y_tilde[:] + a * ((k+1) .- θ_y[:]) .* y[:]
                x_out = x_tilde_tmp / (a * k)
                y_out = y_tilde_tmp / (a * k)

                # Progress measures
                fvaluegap, metricLP = compute_fvaluegap_metricLP(x_out, y_out, problem)

                elapsedtime = time() - starttime
                @info "elapsedtime: $elapsedtime"
                @info "outer_k: $(outer_k), fvaluegap: $(fvaluegap), metricLP: $(metricLP)"
                elapsedtime = time() - starttime
                logresult!(results, outer_k, elapsedtime, fvaluegap, metricLP)
                if !isnothing(io)
                    flush(io)
                end

                exitflag = checkexitcondition(exitcriterion, outer_k, elapsedtime, metricLP)
                if exitflag
                    break
                end

                if k >= restartfreq * m || (restartfreq == Inf && metricLP <= 0.5 * init_metricLP)
                    @info "<===== RESTARTING"
                    @info "k ÷ m: $(k ÷ m)"
                    @info "elapsedtime: $elapsedtime"
                    @info "outer_k: $(outer_k), fvaluegap: $(fvaluegap), metricLP: $(metricLP)"
                    if !isnothing(io)
                        flush(io)
                    end

                    x0, y0 = deepcopy(x_out), deepcopy(y_out)
                    init_fvaluegap = fvaluegap
                    init_metricLP = metricLP
                    restartflag = true
                    break
                end
            end

            k += 1
            outer_k += 1
        end
    end

    return results
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
module ArensdorfOrbit
 #=
 In my ongoing project to get up to speed in Julia, I'm now
 taking a look at "arensdorf orbits", that were designed
 by Richard Arensdorf while he was figuring out how to
 put men on the moon.

 It's a simple diferential equation, and it gives nice plots
 so it's perfect to exercise two parts of my Julia-Fu

 https://www.johndcook.com/blog/2020/02/08/arenstorf-orbit/

=#

  using DifferentialEquations, ForwardDiff, LinearAlgebra
  using OrdinaryDiffEq
  using Plots

  # Bootstrapping using Lorentz attractor.  This is a first order
  # differential equation, but many of the same
  # techniques are used, so it's a good starting point
  function lorenz!(du,u,p,t)
     x, y, z  = u

     du[1] = 10.0*(y-x)
     du[2] = x*(28.0-z) - y
     du[3] = x*y - (8/3)*z
  end

  # ... works abslutely awsomely!
  function plot_lorenz()
    u0 = [1.0;0.0;0.0]
    tspan = (0.0,100.0)
    prob = ODEProblem(lorenz!,u0,tspan)
    sol = solve(prob)
    plot(sol,vars=(1,2,3))
  end

   # .... work nicely too.
   function plot_hh()
     function HH_acceleration(dv,v,u,p,t)
       x,y  = u
       dx,dy = dv
       dv[1] = -x - 2x*y
       dv[2] = y^2 - y -x^2
     end
     tspan = (0.0,100.0)
     initial_positions = [0.0,0.1]
     initial_velocities = [0.5,0.0]
     prob = SecondOrderODEProblem(HH_acceleration,initial_velocities,initial_positions,tspan)
     sol = solve(prob, KahanLi8(), Dt=1/10);
     plot(sol,vars=(1,2,3))
   end

   #  This also works now
   function plot_arensdorf()
     μ   = 0.012277471
     μ′  = 1 - μ

     function arensdorf_orbit(du, u, z, p, t)
         x, y  = z
	 x′,y′ = u

	 D(z) = ((x + z)^2  + y^2)^(3/2)
	 D1 = D( μ)
	 D2 = D(-μ′)

	 x′′ = x + 2y′  - (μ′*(x + μ)/D1) - (μ*(x - μ′)/D2)
	 y′′ = y - 2x′  - (μ′*y/D1)       - (μ*y/D2)

	 du[1] = x′′
	 du[2] = y′′
     end

     # Sensitivity analysis must be added to this code

     initial_positions =  [0.994,  0]
     initial_velocities = [0.0,   -2.001585106]
     tspan = (0.0, 10.0)
     prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
     sol = solve(prob, KahanLi8(), dt=1/100);
     plot(sol,vars=(1,2))
   end
end


text:
@testset "CF DiskArray" begin
a = _DiskArray(zeros(Int8,10,10), chunksize = (2,5))
a2 = CFDiskArray(a,Dict("missing_value"=>-128, "add_offset"=>10, "scale_factor"=>1/128))

@test eltype(a2) == Union{Float64,Missing}

a2[1,:] .= 9.5
a2[2,:] .= 10.5
a2[3,:] .= missing
@test all(isequal(-64), a.parent[1,:])
@test all(isequal(64), a.parent[2,:])
@test all(isequal(-128), a.parent[3,:])
@test all(iszero, a.parent[4:end,:])
@test a2[1,:] == fill(9.5,10)
@test a2[2,:] == fill(10.5,10)
@test all(ismissing,a2[3,:])
end


text:
include("gross_economy.jl")
include("emissions.jl")

function run_my_model()

    my_model = Model()

    setindex(my_model, :time, [2015:5:2110])
    setindex(my_model, :regions, ["Region1", "Region2", "Region3"])  #Note that the regions of your model must be specified here

    addcomponent(my_model, grosseconomy)
    addcomponent(my_model, emissions)

    setparameter(my_model, :grosseconomy, :l, l)
    setparameter(my_model, :grosseconomy, :tfp, tfp)
    setparameter(my_model, :grosseconomy, :s, s)
    setparameter(my_model, :grosseconomy, :depk,depk)
    setparameter(my_model, :grosseconomy, :k0, k0)
    setparameter(my_model, :grosseconomy, :share, 0.3)

    #set parameters for emissions component
    setparameter(my_model, :emissions, :sigma, sigma)
    connectparameter(my_model, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    run(my_model)
    return(my_model)

end


text:
module TestAPI

using Matte

const title = "Test API"

ui = function()
    sidebar_layout((
    h1("hi there"),
    slider("slider","Slider",-100,11)
    )
    ,"")
end

module Server

function foo(a)
    a + 1
end

end

end

run_app(TestAPI, async = true)

# Test that we can actually fetch the javascript that powers Matte
res = HTTP.request("GET", "http://127.0.0.1:8000/bundle.min.js")
@test res.status == 200

# Test that submitting inputs to the backend gives the right result
HTTP.WebSockets.open("ws://127.0.0.1:8000") do ws
    payload = Dict(
        :channel => "matte",
        :message => "api",
        :payload => Dict(
            :id => "foo",
            :session_id => "test_session",
            :input => Dict(
                :a => 1
            )
        )
    )
    expected = Dict(
        "id" => "foo",
        "value" => 2
    )
    write(ws, JSON.json(payload))
    @test JSON.parse(String(readavailable(ws))) == expected
end

stop_app()


text:
T = Complex{Float64}
N = 51
expect_output_type = Vector{Complex{Float64}}

# Test remove_response function
@testset "remove response" begin
    data = ones(T, N)
    freq_response = ones(T, N)
    remove_response!(data, freq_response)
    @test size(data, 1) == N
    @test typeof(data) == expect_output_type
end

# remove_response with large samples
@testset "remove response with large samples" begin
    samples = 90000
    data = ones(T, samples)
    freq_response = ones(T, samples)
    remove_response!(data, freq_response)
    @test size(data, 1) == samples
    @test eltype(data) == T
end


text:
import Git: git

function current_branch_name(dir=".")
    cd(dir) do
        read(`$(git()) rev-parse --abbrev-ref HEAD`, String) |> strip
    end
end

function fetch(dir=".")::Bool

    branch_name = current_branch_name(dir)

    cd(dir) do
        run(`$(git()) fetch origin $(branch_name) --quiet`)

        local_hash = read(`$(git()) rev-parse HEAD`, String)
        remote_hash = read(`$(git()) rev-parse \@\{u\}`, String)

        return local_hash == remote_hash
    end

end

function pullhard(dir=".")

    branch_name = current_branch_name(dir)

    cd(dir) do
        run(`$(git()) reset --hard origin/$(branch_name)`)
    end

end

function fetch_pull(dir=","; pull_sleep=1)
    in_sync = fetch(dir)

    if !in_sync
        pullhard(dir)
        sleep(pull_sleep)
    end
end

function poll_pull_loop(dir="."; interval=5)
    while true
        fetch_pull(dir)
        sleep(interval)
    end
end


text:
__precompile__()

module QuantumInformation

# package code goes here
using QuantumOptics
import QuantumOptics: sigmam, sigmap, sigmax, sigmay, sigmaz
include("export.jl")

QuantumOptics.set_printing(standard_order=true)
# The reasons why the order of tensor product is inverted.
# https://github.com/qojulia/QuantumOptics.jl/blob/395f7077528d36bded0d5ef0652de19b5a9263db/src/printing.jl#L14-L24
# https://qojulia.org/documentation/quantumobjects/operators.html#tensor_order-1

include("basic.jl")
include("gates.jl")
include("printing.jl")

# Alternative Interface (`using QuantumInformation.ShortNames`)
module ShortNames
       using ..QuantumInformation
       include("export.jl")
       include("shortnames.jl")
end # ShortNames
using .ShortNames
end # module


text:
n=4; T = 1

# Cooling Algorithm Test
P = rand(2^n,2^n)
test_cooling_matrix = P*diagm(collect(1:2^n))*inv(P)
eigen = eigfact(test_cooling_matrix)
val,idx = findmin(eigen[:values])
maxEigen = maximum(eigen[:values])
cooling_state = BellState(n)

while !cooling!(cooling_state, DemonCooling(test_cooling_matrix; k=1, tol=0))
    cooling_state = BellState(n)
end

test_cooling_state = BellState(n)

for i=1:1
test_cooling_state.coeffs = normalize((I-im*expm(-im*test_cooling_matrix*π/(2*maxEigen)))*test_cooling_state.coeffs)
end
@test round(abs(dot(coeff(test_cooling_state),eigen[:vectors][:,idx])),5) == round(abs(dot(coeff(cooling_state),eigen[:vectors][:,idx])),5)

# @show round(norm(coeff(aqc.state)[ans[1]+1])^2,3)
# @show abs(dot(coeff(aqc.state),eigen[:vectors][:,idx]))
# @show abs(dot(coeff(cooled_aqc.state),eigen[:vectors][:,idx]))

# show(cooled_aqc)
# print(cooled_aqc)

# test_cooling_state = BellState(n)

# for i=1:200
# test_cooling_state.coeffs = normalize(test_cooling_state.coeffs-im*expm(-im*test_cooling_matrix*0.522697)*test_cooling_state.coeffs)
# end
# @show abs(dot(coeff(test_cooling_state),eigen[:vectors][:,idx]))

# test_cooling_state = BellState(n)

# for i=1:200
# test_cooling_state.coeffs = normalize(test_cooling_state.coeffs-im*exp(im*1.5368)*expm(-im*test_cooling_matrix*2.06334)*test_cooling_state.coeffs)
# end
# @show abs(dot(coeff(test_cooling_state),eigen[:vectors][:,idx]))

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
export COFFSymbols, COFFSymtabEntry, COFFSymbolRef

# We need to import stuff from Base in order to do our tricksy Symbols stuffage
import Base: iteratorsize, SizeUnknown

"""
    COFFSymbols

COFF symbol table, contains the list of symbols defined within the object file.

Note that because COFF Symbols are variable-length, we store a table of offsets
at which the (non-auxilliary) symbols can be found.
"""
struct COFFSymbols{H<:COFFHandle} <: Symbols{H}
    handle::H
    symbol_offsets::Vector{UInt64}
end

"""
    scan_symbol_offsets(oh::COFFHandle)

Find the offsets for each Symbol within a COFFHandle, skipping over auxilliary
symbols as necessary.
"""
function scan_symbol_offsets(oh::H) where {H <: COFFHandle}
    # num_syms represents the number of total symbols, but we're only
    # interested in the non-auxiliary symbols
    num_syms = num_symbols(header(oh))
    offsets = UInt64[]

    curr_idx = 1
    idx = 1
    while idx < num_syms
        # Load in the next symbol and store its offset
        seek(oh, symtab_entry_offset(oh) + symtab_entry_size(oh)*(idx-1))

        sym = unpack(oh, COFFSymtabEntry{H})
        push!(offsets, idx)

        # Move our offset forward by 1 (because we just read in a SymtabEntry)
        # and again as many times as needed to skip over the auxilliary symbols
        idx += 1 + sym.NumberOfAuxSymbols
        curr_idx += 1
    end

    return offsets
end

function Symbols(oh::H) where {H <: COFFHandle}
    symbol_offsets = scan_symbol_offsets(oh)
    return COFFSymbols(oh, symbol_offsets)
end

handle(syms::COFFSymbols) = syms.handle
# function next(syms::COFFSymbols, idx)
#     # We skip over auxiliary symbols here
#     next_idx = idx + deref(syms[idx]).NumberOfAuxSymbols + 1
#     return (syms[idx], next_idx)
# end
endof(syms::COFFSymbols) = num_symbols(header(handle(syms)))

# We override `iteratorsize` so that list comprehensions and collect() calls on
# our Symbols don't have a bunch of #undef entries at the end of the array
# because it tries to pre-allocate an array of the proper size.
# iteratorsize(::Type{H}) where {H <: COFFSymbols} = SizeUnknown()

@io struct COFFSymtabEntry{H <: COFFHandle} <: SymtabEntry{H}
    Name::fixed_string{UInt64}
    Value::UInt32
    SectionNumber::Int16
    Type::UInt16
    StorageClass::UInt8
    NumberOfAuxSymbols::UInt8
end align_packed

function symbol_name(sym::COFFSymtabEntry)
    # COFF Symbols set the first four bytes of the name to zero to signify a
    # name that must be looked up in the string table
    if (sym.Name.data & 0xffffffff) == 0
        return string("strtab@", (sym.Name.data >> 32))
    end
    return unsafe_string(sym.Name)
end

symbol_value(sym::COFFSymtabEntry) = sym.Value
symbol_section(sym::COFFSymtabEntry) = sym.SectionNumber
symbol_type(sym::COFFSymtabEntry) = sym.Type
function isundef(sym::COFFSymtabEntry)
    sym.StorageClass in (
        IMAGE_SYM_CLASS_EXTERNAL_DEF,
        IMAGE_SYM_CLASS_UNDEFINED_LABEL,
        IMAGE_SYM_CLASS_UNDEFINED_STATIC) ||
    (sym.StorageClass == IMAGE_SYM_CLASS_EXTERNAL &&
     sym.SectionNumber == IMAGE_SYM_CLASS_NULL)
end
isglobal(sym::COFFSymtabEntry) = isundef(sym) || sym.StorageClass == IMAGE_SYM_CLASS_EXTERNAL
islocal(sym::COFFSymtabEntry) = !isglobal(sym)
isweak(sym::COFFSymtabEntry) = sym.StorageClass == IMAGE_SYM_CLASS_WEAK_EXTERNAL

"""
    COFFSymbolRef

Contains a reference to an `COFFSymtabEntry`, as well as an `COFFSymbols`, etc...
"""
struct COFFSymbolRef{H<:COFFHandle} <: SymbolRef{H}
    syms::COFFSymbols{H}
    entry::COFFSymtabEntry{H}
    idx::UInt32
end
function SymbolRef(syms::COFFSymbols, entry::COFFSymtabEntry, idx)
    return COFFSymbolRef(syms, entry, UInt32(idx))
end

deref(sym::COFFSymbolRef) = sym.entry
Symbols(sym::COFFSymbolRef) = sym.syms
symbol_number(sym::COFFSymbolRef) = sym.idx
@derefmethod symbol_type(sym::COFFSymbolRef)
function symbol_name(sym::COFFSymbolRef)
    # COFF Symbols set the first four bytes of the name to zero to signify a
    # name that must be looked up in the string table
    name = deref(sym).Name
    if (name.data & 0xffffffff) == 0
        # If it's zeroes all the way down, just return an empty string
        if name.data == 0
            return ""
        end
        return strtab_lookup(StrTab(handle(sym)), (name.data >> 32))
    end
    return unsafe_string(name)
end

# Symbol printing stuff
"""
    symbol_type_string(s::COFFSymtabEntry)

Return the given `COFFSymtabEntry`'s section type as a string.
"""
function symbol_type_string(sym::COFFSymtabEntry)
    global IMAGE_SYM_TYPE

    # fast-track a zero symbol type
    if sym.Type == 0
        return ""
    end

    # We have to construct a compound type string here
    type_string = ""

    # We start with the "derived" type, e.g. pointer, array, etc...
    dtype = (sym.Type & 0xf0) >> 4
    if haskey(IMAGE_SYM_DTYPE, dtype)
        type_string = string("(", IMAGE_SYM_DTYPE[dtype], ")")
    end

    # Then, move on to the "base" type, e.g. void, char, short...
    btype = sym.Type & 0x0f
    if haskey(IMAGE_SYM_TYPE, btype)
        type_string = "$(type_string) $(IMAGE_SYM_TYPE[btype])"
    end

    # If we didn't get ANY info, complain
    if isempty(type_string)
        string("Unknown Symbol Type (0x", hex(sym.Type), ")")
    end

    return type_string
end
@derefmethod symbol_type_string(s::COFFSymbolRef)


text:
function loaddata(filename::AbstractString)
    numstr, instrstr = split(readchomp(filename), r"\r?\n\r?\n")
    nums = [CartesianIndex(reverse(parse.(Int, t) .+ 1)...) for t in split.(split(numstr, r"\r?\n"), ',')]
    maxi = maximum(nums)
    matrix = zeros(Bool, (maxi[1] + iseven(maxi[1]), maxi[2]))
    matrix[nums] .= true
    matrix, split(instrstr, r"\r?\n")
end

if basename(pwd()) == "aoc"
    cd("2021/13")
end

atlas, instrs = loaddata("example.txt")

function foldalong(atlas, instr::AbstractString)
    d, v = match(r"fold along ([xy])=(\d+)", instr)
    n = parse(Int, v) + 1
    println("Folding $(size(atlas)) at $d=$n")
    if d == "y"
        foldhorizontal(atlas, n)
    else
        foldvertical(atlas, n)
    end
end

function foldhorizontal(atlas::AbstractMatrix, y::Integer)
    top = atlas[begin:y-1, :]
    bottom = atlas[y+1:end, :]
    top .|= reverse(bottom, dims = 1)
end

function foldvertical(atlas::AbstractMatrix, x::Integer)
    left = atlas[:, begin:x-1]
    right = atlas[:, x+1:end]
    left .|= reverse(right, dims = 2)
end

part1(filename::AbstractString) = part1(loaddata(filename)...)
function part1(atlas::AbstractMatrix, instrs::Vector)
    count(foldalong(atlas, first(instrs)))
end

part2(filename::AbstractString) = part2(loaddata(filename)...)
function part2(atlas::AbstractMatrix, instrs)
    for instr in instrs
        atlas = foldalong(atlas, instr)
    end
    atlas
end


text:
function EnumeratedRealDistribution(arg0::Vector{jdouble})
    return EnumeratedRealDistribution((Vector{jdouble},), arg0)
end

function EnumeratedRealDistribution(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return EnumeratedRealDistribution((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function cumulative_probability(obj::EnumeratedRealDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::EnumeratedRealDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_numerical_mean(obj::EnumeratedRealDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::EnumeratedRealDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_pmf(obj::EnumeratedRealDistribution)
    return jcall(obj, "getPmf", List, ())
end

function get_support_lower_bound(obj::EnumeratedRealDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::EnumeratedRealDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::EnumeratedRealDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::EnumeratedRealDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function probability(obj::EnumeratedRealDistribution, arg0::jdouble)
    return jcall(obj, "probability", jdouble, (jdouble,), arg0)
end