text:
module WindowManager

import GLFW
import GLFW.Window

WINDOW_HANDLER_ID_COUNTER = 0
lastKey=Number(0)
KEYS=Dict{Number,Number}()

type WindowHandler
	id						::Symbol
	name					::AbstractString
	pos						::Tuple{Number,Number}
	size					::Tuple{Number,Number}
	lastSize			::Tuple{Number,Number}
	listenList		::Dict{Symbol, Function}
	events				::Array{Array{Any,1},1}
	lastPos	      ::AbstractVector
	lastCursorPos	::AbstractVector
	nativeWindow	::Window
	fullScreen		::Bool

	function WindowHandler(name::AbstractString, size::Tuple{Number,Number})
		global WINDOW_HANDLER_ID_COUNTER

		if WINDOW_HANDLER_ID_COUNTER == 0
			GLFW.Init()
		end

    new(
      Symbol("display"*string(WINDOW_HANDLER_ID_COUNTER+=1)),
			name,
			(0,0),
			size,
			size,
			Dict(),
			[],[],[],
			Window(GLFW.WindowHandle(C_NULL)),
			false
		)
	end
end

WINDOW_TO_HANDLER_DICT = Dict{Window, WindowHandler}()

windowHandler = nothing

getWindowHandler() = windowHandler
setWindowHandler(this::WindowHandler) = (global windowHandler = this)

terminate() = GLFW.Terminate()

function setListener(this::WindowHandler, key::Symbol, listener::Function)
	this.listenList[key] = listener
end

################################################################################
# Window Events

function OnWindowResize(window::Window, width::Number, height::Number)
	h = WINDOW_TO_HANDLER_DICT[window]
	h.size = (width,height)
	OnEvent(h, :OnWindowResize, [width,height])
end

OnWindowIconify(window::Window, iconified::Number) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnWindowIconify, iconified)
OnWindowClose(window::Window) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnWindowClose)
OnWindowFocus(window::Window, focused::Number) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnWindowFocus, focused)
OnWindowRefresh(window::Window) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnWindowRefresh)
OnFramebufferResize(window::Window, width::Number, height::Number) =  OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnFramebufferResize, [width,height])
OnCursorEnter(window::Window, entered::Number) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnCursorEnter, entered)
OnDroppedFiles(window::Window, files::AbstractArray) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnDroppedFiles, files)

function OnWindowPos(window::Window, x::Number, y::Number)
  h = WINDOW_TO_HANDLER_DICT[window]
  t = [-x,y]
  OnEvent(h, :OnWindowPos, t)
  l = h.lastPos; if l == [] l = t end
  n = l - t
  h.lastPos = t
  OnEvent(h, :OnWindowMove, n)
end

function OnCursorPos(window::Window, x::Number, y::Number)
  h = WINDOW_TO_HANDLER_DICT[window]
  t = [x,y]
  OnEvent(h, :OnMousePos, t)
  l = h.lastCursorPos; if l == [] l = t end
  n = l - t
  h.lastCursorPos = t
  OnEvent(h, :OnMouseMove, n)
end

# Input Events
OnScroll(window::Window, x::Number, y::Number) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnScroll, [x,y])
OnCharMods(window::Window, code::Char, mods::Number) = OnEvent(WINDOW_TO_HANDLER_DICT[window], :OnCharMods, code, mods)

function OnUpdateEvents()
  GLFW.PollEvents()  # Poll for and process events
	for (window,h) in WINDOW_TO_HANDLER_DICT
		for e in h.events
			if length(e) > 5 && e[1] == :OnKey
				key=e[2]
				if haskey(KEYS, key) e[6]=KEYS[key] end
			end
			OnEvent(h, e...)
		end
		if length(h.events) > 0	h.events	= [] end
	end
end

function OnKey(window::Window, key::Number, scancode::Number, action::Number, mods::Number)
	h = WINDOW_TO_HANDLER_DICT[window]
	global lastKey = key
	push!(h.events, [ :OnKey, key, scancode, action, mods, Number(0) ]) # unicode = 0
end

function OnMouseKey(window::Window, key::Number, action::Number, mods::Number)
  h = WINDOW_TO_HANDLER_DICT[window]
  OnEvent(h, :OnMouseKey, key, action, mods)
end

function OnUnicodeChar(window::Window, unicode::Char)
	h = WINDOW_TO_HANDLER_DICT[window]
	if lastKey > 0 KEYS[lastKey]=Number(unicode) end
	global lastKey=0
end

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

function ShowError(debugging::Bool)
  @static if is_apple() return end
  GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT, convert(Cint, debugging))
end

function OnEvent(this::WindowHandler, eventName::Symbol, args...)
 #a = values(d) # convert Dict Values to Array
  for (k,listener) in this.listenList
		listener(eventName, args...)
    #s = listener.storage
		#if haskey(s.events, eventName) s.events[eventName](args...) end
  end
end

title(this::WindowHandler, name::String) = GLFW.SetWindowTitle(this.nativeWindow, name)

function cursor(this::WindowHandler, mode::Symbol)
	GLFW.SetInputMode(this.nativeWindow, GLFW.CURSOR, eval(:(GLFW.$mode)))
	#GLFW.STICKY_KEYS or GLFW.STICKY_MOUSE_BUTTONS
end

function fullscreen(this::WindowHandler, full::Bool)
	monitor = GLFW.GetPrimaryMonitor()
	mode = GLFW.GetVideoMode(monitor)
	pos = GLFW.GetMonitorPos(monitor)

	if !full
		pos = this.pos
		this.size = this.lastSize
		mode = GLFW.VidMode(this.size[1], this.size[2], mode.redbits, mode.greenbits, mode.bluebits, mode.refreshrate)
		monitor = GLFW.Monitor(C_NULL)
	else
		this.lastSize = this.size
		this.pos = GLFW.GetWindowPos(this.nativeWindow)
		this.size = GLFW.GetWindowSize(this.nativeWindow)
	end

	this.fullScreen = full
	SetWindowMonitor(this.nativeWindow, monitor, pos[1], pos[2], mode.width, mode.height, mode.refreshrate)
end

# GLFW.SetWindowMonitor
function SetWindowMonitor(window::Window, monitor::GLFW.Monitor, xpos, ypos, width, height, refreshRate)
    ccall((:glfwSetWindowMonitor, GLFW.lib), Void, (GLFW.WindowHandle, GLFW.Monitor, Cint, Cint, Cint, Cint, Cint), window, monitor, xpos, ypos, width, height, refreshRate)
end

# Create a window and its OpenGL context
function open(this::WindowHandler, windowhints=[])
    for (id,val) in windowhints
        GLFW.WindowHint(eval(:(GLFW.$id)),val)
    end

    w = this.size[1]; h = this.size[2]

    window = GLFW.CreateWindow(w, h, this.name)
    this.nativeWindow = window

		this.pos = GLFW.GetWindowPos(this.nativeWindow)

    WINDOW_TO_HANDLER_DICT[window] = this

  	GLFW.MakeContextCurrent(window)

    GLFW.SwapInterval(0)

    GLFW.ShowWindow(window)

    #GLFW.SetWindowSize(window, w, h) # Seems to be necessary to guarantee that window > 0

    # Window Callbacks
    GLFW.SetWindowIconifyCallback(window, OnWindowIconify)
    GLFW.SetWindowSizeCallback(window, OnWindowResize)
    GLFW.SetWindowCloseCallback(window, OnWindowClose)
    GLFW.SetWindowFocusCallback(window, OnWindowFocus)
    GLFW.SetWindowPosCallback(window, OnWindowPos)
    GLFW.SetWindowRefreshCallback(window, OnWindowRefresh)
    GLFW.SetFramebufferSizeCallback(window, OnFramebufferResize)
    GLFW.SetCursorPosCallback(window, OnCursorPos)
    GLFW.SetCursorEnterCallback(window, OnCursorEnter)
    GLFW.SetDropCallback(window, OnDroppedFiles)

	   # Input Callbacks
    GLFW.SetKeyCallback(window, OnKey)
    GLFW.SetMouseButtonCallback(window, OnMouseKey)
    GLFW.SetScrollCallback(window, OnScroll)
    GLFW.SetCharCallback(window, OnUnicodeChar)
    GLFW.SetCharModsCallback(window, OnCharMods)
end

isClosing(this::Window) = GLFW.WindowShouldClose(this.nativeWindow)
isClosing(this::WindowHandler) = GLFW.WindowShouldClose(this.nativeWindow)

function create(name::AbstractString, size::Tuple{Number,Number})
	h=WindowHandler(name,size)
	setWindowHandler(h)
	h
end

function close(this::WindowHandler)
	global WINDOW_HANDLER_ID_COUNTER
	if WINDOW_HANDLER_ID_COUNTER == 1 GLFW.Terminate() end
	if WINDOW_HANDLER_ID_COUNTER > 0 WINDOW_HANDLER_ID_COUNTER -= 1 end
end

function swap(this::WindowHandler)
  GLFW.SwapBuffers(this.nativeWindow)   # Swap front and back buffers
end

function update(this::WindowHandler)
  swap(this)
	OnUpdateEvents()
end

""" Loop until the user closes the window """
function loop(this::WindowHandler, repeat::Function)
	while !isClosing(this)
		repeat()
		update(this) # update window
	end
	close(this)
end

end # WindowManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# <?xml version="1.0" encoding="utf-8" standalone="yes" ?>
# <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9">
# <url>
#      <loc>http://www.example.com/</loc>
#      ?<lastmod>2005-01-01</lastmod>
#      ?<changefreq>monthly</changefreq>
#      ?<priority>0.8</priority>
# </url>
# </urlset>
#

struct SMOpts
    lastmod::Date        # 2005-01-01 -- format YYYY-MM-DD
    changefreq::String   # one of...
    priority::Float64    #
    function SMOpts(l, c, p)
        c = lowercase(c)
        p = isa(p, String) ? parse(Float64, p) : float(p)
        allowedf = ("always", "hourly", "daily", "weekly", "monthly",
        "yearly", "never")
        assertf = """
        Given change frequency on $(FD_ENV[:SOURCE]) is invalid according to
        sitemap specifications, expected one of $allowedf.
        """
        assertp = """
        Given priority on $(FD_ENV[:SOURCE]) is invalid according to sitemap
        specifications, expected a floating point number between 0 and 1.
        """
        @assert c in allowedf assertf
        @assert 0 <= p <= 1 assertp
        return new(l, c, p)
    end
end

const SITEMAP_DICT = LittleDict{String,SMOpts}()

"""
$SIGNATURES

Add an entry to `SITEMAP_DICT`.
"""
function add_sitemap_item(; html=false)
    loc = url_curpage()
    locvar(:sitemap_exclude)::Bool && return nothing
    if !html
        lastmod = locvar(:fd_mtime_raw)
        changefreq = locvar(:sitemap_changefreq)
        priority = locvar(:sitemap_priority)
    else
        # use default which can be overwritten in a {{sitemap_opts ...}}
        fp = joinpath(path(:folder), locvar(:fd_rpath)::String)
        lastmod = Date(unix2datetime(stat(fp).mtime))
        changefreq = "monthly"
        priority = 0.5
    end
    res = SITEMAP_DICT[loc] = SMOpts(lastmod, changefreq, priority)
    return res
end

"""
$SIGNATURES

Generate a `sitemap.xml`, if one already exists, it will be replaced.
"""
function sitemap_generator()
    dst = joinpath(path(:site), "sitemap.xml")
    isfile(dst) && rm(dst)
    io = IOBuffer()
    println(io, """
        <?xml version="1.0" encoding="utf-8" standalone="yes" ?>
        <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9">
        """)
    base_url = globvar(:website_url)
    for (k, v) in SITEMAP_DICT
        key = joinpath(escapeuri.(split(k, '/'))...)
        loc = "<loc>$(joinpath(base_url, key))</loc>"
        lastmod = "<lastmod>$(v.lastmod)</lastmod>"
        changefreq = "<changefreq>$(v.changefreq)</changefreq>"
        priority = "<priority>$(v.priority)</priority>"
        write(io, """
            <url>
                $loc
                $lastmod
                $changefreq
                $priority
            </url>
            """)
    end
    println(io, "</urlset>")
    write(dst, take!(io))
    return nothing
end

ext: jl


text:
function decompose_vector(v::Vector3, k::Vector3)
    v_parallel = (v ⋅ k) * k
    v_perpendicular = v - v_parallel
    v_parallel, v_perpendicular
end

function rotate_around(v::Vector3, k::Vector3, φ::Real)
    k = normalize(k)
    v_parallel, v_perpendicular = decompose_vector(v, k)
    ω = k × v_perpendicular
    v_parallel + v_perpendicular * cos(φ) + ω * sin(φ)
end

@_with_default_rng function random_3d_direction()
    θ = 2π * rand(rng)
    z = 2 * rand(rng) - 1
    r = √(1 - z ^ 2)
    Vector3(r * cos(θ), r * sin(θ), z)
end

"""
    random_direction_on_plane([rng], normal_vector)
    random_direction_on_plane([rng], u, v)

Generate a random direction on a plane defined by a normal vector or two vectors on the plane.
"""
random_direction_on_plane(rng::AbstractRNG, u::Vector3, v::Vector3) = let φ = 2π * rand(rng)
    w = cos(φ) * u + sin(φ) * v
end
@_with_default_rng function random_direction_on_plane(normal_vector::Vector3)
    a, b, c = normal_vector
    u = normalize(Vector3(b - c, -a + c, a - b))
    v = normalize(normal_vector × u)
    random_direction_on_plane(rng, u, v)
end


text:
using CodeCosts
using CodeCosts: costs
using InteractiveUtils
using Test

if get(ENV, "TRAVIS", "false") == "true"
    # force :color => true
    Base.get(::Base.PipeEndpoint, key::Symbol, default) = key === :color ? true : default
    # avoid the `light` colors
    ENV["JULIA_ERROR_COLOR"] = :red
end

@inline func_000(x::T) where T = (x + 1/3, x - oneunit(T), x / 3.0, x * 5.0)

func_001(x) = :(if $x > 0; true; else false; end)

buf = IOBuffer()

v"1.6.0-DEV.0" <= VERSION < v"1.7.0-DEV.0" && @testset "Basic v1.6" begin
    @testset "func_000" begin
        res = @code_costs func_000(1.0f0)
        show(res)
        println()
        @test costs(res) == [1, 1, 1, 1, 20, 1, 4, 0, 0]

        print(buf, res.summary)
        @test String(take!(buf)) == """
                                    CodeCostsSummary(
                                         zero:  2|
                                        cheap:  5| 11111
                                       middle:  4| 4===
                                    expensive: 20| 20==================
                                        total: 29| 100 (default threshold)
                                    )"""

        res = @code_costs debuginfo=:source func_000(1.0)
        show(res)
        println()
        @test costs(res) == [1, 1, 20, 4, 0, 0]
    end

    @testset "func_001" begin
        res = @code_costs func_001(1)
        show(res)
        println()

        print(buf, res)
        @test occursin(r"\n   0 └──      return %\d\n", String(take!(buf)))
        @test costs(res) == [100, 100, 100, 100, 0]
    end
end

v"1.5.0-DEV.0" <= VERSION < v"1.6.0-DEV.0" && @testset "Basic v1.5" begin
    @testset "func_000" begin
        res = @code_costs func_000(1.0f0)
        show(res)
        println()
        @test costs(res) == [1, 1, 1, 1, 20, 1, 4, 0, 0]

        print(buf, res.summary)
        @test String(take!(buf)) == """
                                    CodeCostsSummary(
                                         zero:  2|
                                        cheap:  5| 11111
                                       middle:  4| 4===
                                    expensive: 20| 20==================
                                        total: 29| 100 (default threshold)
                                    )"""

        res = @code_costs debuginfo=:source func_000(1.0)
        show(res)
        println()
        @test costs(res) == [1, 1, 20, 4, 0, 0]
    end

    @testset "func_001" begin
        res = @code_costs func_001(1)
        show(res)
        println()

        print(buf, res)
        @test occursin(r"\n   0 └──      return %\d\n", String(take!(buf)))
        @test costs(res) == [100, 100, 100, 100, 0]
    end
end

v"1.4" <= VERSION < v"1.5.0-DEV.0" && @testset "Basic v1.4" begin
    @testset "func_000" begin
        res = @code_costs func_000(1.0f0)
        show(res)
        println()
        @test costs(res) == [1, 1, 1, 1, 20, 1, 4, 0, 0]

        print(buf, res.summary)
        @test String(take!(buf)) == """
                                    CodeCostsSummary(
                                         zero:  2|
                                        cheap:  5| 11111
                                       middle:  4| 4===
                                    expensive: 20| 20==================
                                        total: 29| 100 (default threshold)
                                    )"""

        res = @code_costs debuginfo=:source func_000(1.0)
        show(res)
        println()
        @test costs(res) == [1, 1, 20, 4, 0, 0]
    end

    @testset "func_001" begin
        res = @code_costs func_001(1)
        show(res)
        println()

        print(buf, res)
        @test occursin(r"\n   0 └──      return %\d\n", String(take!(buf)))
        @test costs(res) == [100, 100, 100, 100, 0]
    end
end

v"1.0" <= VERSION < v"1.1" && @testset "Basic v1.0" begin
    @testset "func_000" begin
        res = @code_costs func_000(1.0f0)
        show(res)
        println()
        @test costs(res) == [1, 1, 1, 1, 20, 1, 4, 0, 0]

        print(buf, res.summary)
        @test String(take!(buf)) == """
                                    CodeCostsSummary(
                                         zero:  2|
                                        cheap:  5| 11111
                                       middle:  4| 4===
                                    expensive: 20| 20==================
                                        total: 29| 100 (default threshold)
                                    )"""
        res = @code_costs debuginfo=:source func_000(1.0)
        show(res)
        println()
        @test costs(res) == [1, 1, 20, 4, 0, 0]
    end

    @testset "func_001" begin
        res = @code_costs func_001(1)
        show(res)
        println()

        print(buf, res)
        @test occursin(r"\n   0    └──      return %\d\n", String(take!(buf)))
        @test costs(res) == [100, 100, 100, 100, 0]
    end
end


text:
using NCDatasets
using OrderedCollections

struct NetCDFWriter <: AbstractWriter end

function write_data(nc::NetCDFWriter, filename, dims, varvals, simtime)
    Dataset(full_name(nc, filename), "c") do ds
        ds.dim["time"] = 1
        for (dn, (dv, da)) in dims
            ds.dim[dn] = length(dv)
        end
        defVar(
            ds,
            "time",
            [simtime],
            ("time",),
            attrib = OrderedDict(
                "units" => "seconds",
                "long_name" => "simulation time",
            ),
        )
        for (dn, (dv, da)) in dims
            defVar(ds, dn, dv, (dn,), attrib = da)
        end
        for (vn, (vd, vv, va)) in varvals
            defVar(ds, vn, vv, vd, attrib = va)
        end
    end
    return nothing
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
@testset "Default post correlation Filter" begin
    post_corr_filter = Tracking.get_default_post_corr_filter(EarlyPromptLateCorrelator())
    @test @inferred(post_corr_filter(1.0 + 0.0im)) == 1.0 + 0.0im

    post_corr_filter = Tracking.get_default_post_corr_filter(
        EarlyPromptLateCorrelator(NumAnts(2))
    )
    @test @inferred(post_corr_filter(SVector(1.0 + 0.0im, 2.0 + 0.0im))) == 1.0 + 0.0im
end

@testset "Resize downconverted signal ($type) for multiple ants" for type in (Int16, Int32, Int64, Float32, Float64)
    downconverted_signal_temp = Tracking.DownconvertedSignalCPU(NumAnts(4))
    signal = StructArray(ones(Complex{type}, 2500, 4))
    downconverted_signal = Tracking.resize!(downconverted_signal_temp, 2500, signal)
    if type == Float64
        @test size(downconverted_signal.downconverted_signal_f32) == (0, 4)
        @test size(downconverted_signal.downconverted_signal_f64) == (2500, 4)
    else
        @test size(downconverted_signal.downconverted_signal_f32) == (2500, 4)
        @test size(downconverted_signal.downconverted_signal_f64) == (0, 4)
    end
end

@testset "Resize downconverted signal ($type) for single ant" for type in (Int16, Int32, Int64, Float32, Float64)
    downconverted_signal_temp = Tracking.DownconvertedSignalCPU(NumAnts(1))
    signal = StructArray(ones(Complex{type}, 2500))
    downconverted_signal = Tracking.resize!(downconverted_signal_temp, 2500, signal)
    if type == Float64
        @test size(downconverted_signal.downconverted_signal_f32) == (0,)
        @test size(downconverted_signal.downconverted_signal_f64) == (2500,)
    else
        @test size(downconverted_signal.downconverted_signal_f32) == (2500,)
        @test size(downconverted_signal.downconverted_signal_f64) == (0,)

    end
end

@testset "Integration time" begin
    gpsl1 = GPSL1()
    galileo_e1b = GalileoE1B()
    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found)
    @test time == 1ms

    max_integration_time = 2ms
    bit_found = true
    time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found)
    @test time == 2ms

    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        galileo_e1b,
        max_integration_time,
        bit_found
    )
    @test time == 2ms

    max_integration_time = 10ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        galileo_e1b,
        max_integration_time,
        bit_found
    )
    @test time == 4ms
end

@testset "Number of chips to integrate" begin
    gpsl1 = GPSL1()
    max_integration_time = 1ms
    code_phase = 200
    bit_found = true
    chips = @inferred Tracking.get_num_chips_to_integrate(
        gpsl1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 1023 - 200

    code_phase = 1200
    chips = @inferred Tracking.get_num_chips_to_integrate(
        gpsl1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 2046 - 1200
end

@testset "Number of samples to integrate" begin
    gpsl1 = GPSL1()
    max_integration_time = 1ms
    sampling_frequency = 4e6Hz
    current_code_doppler = 0Hz
    current_code_phase = 0
    bit_found = true
    samples = @inferred Tracking.get_num_samples_left_to_integrate(
        gpsl1,
        max_integration_time,
        sampling_frequency,
        current_code_doppler,
        current_code_phase,
        bit_found
    )
    @test samples == 4000
end

@testset "Carrier phase delta" begin
    intermediate_frequency = 100Hz
    carrier_doppler = 50Hz
    frequency = @inferred Tracking.get_current_carrier_frequency(
        intermediate_frequency,
        carrier_doppler
    )
    @test frequency == intermediate_frequency + carrier_doppler
end

@testset "Code phase delta" begin
    gpsl1 = GPSL1()
    code_doppler = 1Hz
    frequency = @inferred Tracking.get_current_code_frequency(gpsl1, code_doppler)
    @test frequency == 1023e3Hz + 1Hz
end

@testset "Doppler aiding" begin
    gpsl1 = GPSL1()
    init_carrier_doppler = 10Hz
    init_code_doppler = 1Hz
    carrier_freq_update = 2Hz
    code_freq_update = -0.5Hz
    velocity_aiding = 3Hz

    carrier_freq, code_freq = @inferred Tracking.aid_dopplers(
        gpsl1,
        init_carrier_doppler,
        init_code_doppler,
        carrier_freq_update,
        code_freq_update,
        velocity_aiding
    )

    @test carrier_freq == 10Hz + 2Hz + 3Hz
    @test code_freq == 1Hz + (2Hz + 3Hz) / 1540 - 0.5Hz
end

@testset "Tracking with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64)
    gpsl1 = GPSL1()
    carrier_doppler = 200Hz
    start_code_phase = 100
    code_frequency = carrier_doppler / 1540 + 1023e3Hz
    sampling_frequency = 4e6Hz
    prn = 1
    range = 0:3999
    start_carrier_phase = π / 2
    state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase)

    signal_temp = cis.(
            2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase
        ) .*
        get_code.(
            gpsl1,
            code_frequency .* range ./ sampling_frequency .+ start_code_phase,
            prn
        )
    scaling = 512
    signal = type <: Integer ? complex.(floor.(type, real.(signal_temp) * scaling), floor.(type, imag.(signal_temp) * scaling)) : Complex{type}.(signal_temp)

    track_result = @inferred track(signal, state, sampling_frequency)

    iterations = 2000
    code_phases = zeros(iterations)
    carrier_phases = zeros(iterations)
    tracked_code_phases = zeros(iterations)
    tracked_carrier_phases = zeros(iterations)
    tracked_code_dopplers = zeros(iterations)
    tracked_carrier_dopplers = zeros(iterations)
    tracked_prompts = zeros(ComplexF64, iterations)
    for i = 1:iterations
        carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency +
            start_carrier_phase + π) - π
        code_phase = mod(
            code_frequency * 4000 * i / sampling_frequency + start_code_phase,
            1023
        )
        signal = cis.(
                2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase
            ) .*
            get_code.(
                gpsl1,
                code_frequency .* range ./ sampling_frequency .+ code_phase,
                prn
            )
        track_result = @inferred track(
            signal,
            get_state(track_result),
            sampling_frequency
        )
        comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) /
            sampling_frequency + start_carrier_phase + π) - π
        comp_code_phase = mod(
            code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase,
            1023
        )
        tracked_code_phases[i] = get_code_phase(track_result)
        tracked_carrier_phases[i] = get_carrier_phase(track_result)
        tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz
        tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz
        tracked_prompts[i] = get_prompt(track_result)
        code_phases[i] = comp_code_phase
        carrier_phases[i] = comp_carrier_phase
    end
    @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4
    @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2

#    using PyPlot
#    pygui(true)
#    figure("carrier_phases")
#    plot(tracked_carrier_phases)
#    plot(carrier_phases)
#    grid(true)
#    figure("Code Phases")
#    plot(300 * (tracked_code_phases .- code_phases))
#    figure("Carrier Doppler")
#    plot(tracked_carrier_dopplers)
#    figure("Code Doppler")
#    plot(tracked_code_dopplers)
#    figure("Prompt")
#    plot(real.(tracked_prompts))
#    plot(imag.(tracked_prompts))
end

@testset "Track multiple signals with signal of type $type" for type in (Int16, Int32, Int64, Float32, Float64)
    gpsl1 = GPSL1()
    carrier_doppler = 200Hz
    start_code_phase = 100
    code_frequency = carrier_doppler / 1540 + 1023e3Hz
    sampling_frequency = 4e6Hz
    prn = 1
    range = 0:3999
    start_carrier_phase = π / 2
    state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_ants = NumAnts(3))

    signal = cis.(
            2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase
        ) .*
        get_code.(
            gpsl1,
            code_frequency .* range ./ sampling_frequency .+ start_code_phase,
            prn
        )
    signal_mat_temp = repeat(signal, outer = (1,3))
    scaling = 512
    signal_mat = type <: Integer ? complex.(floor.(type, real.(signal_mat_temp) * scaling), floor.(type, imag.(signal_mat_temp) * scaling)) : Complex{type}.(signal_mat_temp)

    @test_throws ArgumentError track(signal_mat', state, sampling_frequency)

    track_result = @inferred track(signal_mat, state, sampling_frequency)

    iterations = 2000
    code_phases = zeros(iterations)
    carrier_phases = zeros(iterations)
    tracked_code_phases = zeros(iterations)
    tracked_carrier_phases = zeros(iterations)
    tracked_code_dopplers = zeros(iterations)
    tracked_carrier_dopplers = zeros(iterations)
    for i = 1:iterations
        carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency +
            start_carrier_phase + π) - π
        code_phase = mod(
            code_frequency * 4000 * i / sampling_frequency + start_code_phase,
            1023
        )
        signal = cis.(
                2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase
            ) .*
            get_code.(
                gpsl1,
                code_frequency .* range ./ sampling_frequency .+ code_phase,
                prn
            )
        signal_mat = repeat(signal, outer = (1,3))
        track_result = @inferred track(
            signal_mat,
            get_state(track_result),
            sampling_frequency
        )
        comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) /
            sampling_frequency + start_carrier_phase + π) - π
        comp_code_phase = mod(
            code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase,
            1023
        )
        tracked_code_phases[i] = get_code_phase(track_result)
        tracked_carrier_phases[i] = get_carrier_phase(track_result)
        tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz
        tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz
        code_phases[i] = comp_code_phase
        carrier_phases[i] = comp_carrier_phase
    end
    @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4
    @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2

#    using PyPlot
#    pygui(true)
#    figure("carrier_phases")
#    plot(tracked_carrier_phases)
#    plot(carrier_phases)
#    grid(true)
#    figure("Code Phases")
#    plot(300 * (tracked_code_phases .- code_phases))
#    figure("Carrier Doppler")
#    plot(tracked_carrier_dopplers)
#    figure("Code Doppler")
#    plot(tracked_code_dopplers)
end


text:
if !("." in LOAD_PATH) # for easier local testing
    push!(LOAD_PATH, ".")
end
using Test
import TurbulenceConvection
const TC = TurbulenceConvection
const tc_dir = dirname(dirname(pathof(TC)))

include(joinpath(tc_dir, "driver", "main.jl"))
include(joinpath(tc_dir, "driver", "generate_namelist.jl"))
include(joinpath(tc_dir, "post_processing", "compute_mse.jl"))
include(joinpath(tc_dir, "post_processing", "mse_tables.jl"))
import .NameList
best_mse = all_best_mse["Bomex"]

@testset "Unit tests" begin
    # TODO: add unit tests.
end

@testset "Bomex" begin

    case_name = "Bomex"
    println("Running $case_name...")
    namelist = NameList.default_namelist(case_name)
    namelist["meta"]["uuid"] = "01"
    ds_tc_filename, return_code = main(namelist)

    computed_mse = compute_mse_wrapper(
        case_name,
        best_mse,
        ds_tc_filename;
        plot_comparison = true,
        t_start = 4 * 3600,
        t_stop = 6 * 3600,
    )
    nothing
end

@testset "ClimaCore extensions" begin
    include("clima_core_extensions.jl")
end


text:
"""
    `Body2D(X,Y)`

[insert description here]

# ARGUMENTS
* `X::Array{Float64}`           : List of X body points
* `Y::Array{Float64}`           : List of Y body points

# OUTPUTS
* `body::Array{LineSegment}`      : Array of Panel2D objects representing a body
"""
function Body2D(X::Array{Float64},
                Y::Array{Float64})

    # Iterate around body to create geometry
    NPTS = length(X)
    NPAN = NPTS - 1
    body = [LineSegment([X[n], Y[n]], [X[n+1], Y[n+1]]) for n = 1:NPTS-1]
    # body = Array{LineSegment,NPAN}
    # for n=1:NPTS-1
    #     r1 = [X[n];Y[n]]
    #     r2 = [X[n+1];Y[n+1]]
    #     body[n] = LineSegment(r1,r2)
    # end
    return body
end


text:
using CuArrays, CUDAnative, GPUArrays
using NiLang, NiLang.AD

const RotGates = Union{Val{:Rz}, Val{:Rx}, Val{:Ry}}

@i @inline function instruct!(state::CuVector, gate::RotGates, loc::Int, theta::Real)
    XY ← GPUArrays.thread_blocks_heuristic(length(state))
    mask ← 1<<(loc-1)
    @cuda threads=tget(XY,1) blocks=tget(XY,2) rot_kernel(gate, state, mask, theta)
end

@i @inline function rot_kernel(gate::Val{:Rz}, state, mask, θ)
    @invcheckoff b ← (blockIdx().x-1) * blockDim().x + threadIdx().x
    @invcheckoff if (b < length(state) && b & mask == 0, ~)
        ROT_INSTRUCT(gate, state[b+1], state[b⊻mask+1], θ)
    end
end

@i @inline function ROT_INSTRUCT(gate::Val{:Rz}, a::T, b, θ) where T
    # make sure `invcheck` is turned off!
    @routine @invcheckoff begin
        @zeros T anc1 anc2 anc3 anc4
        anc1 += θ*(0.5im)
        anc2 += CUDAnative.exp(anc1)
    end
    anc3 += a * anc2'
    anc4 += b * anc2
    NiLang.SWAP(a, anc3)
    NiLang.SWAP(b, anc4)
    anc3 -= a / anc2'
    anc4 -= b / anc2
    ~@routine
end

v = randn(ComplexF64, 128) |> CuArray
v1 = instruct!(copy(v), Val(:Rz), 3, 0.5)[1]
# we can not obtain the gradient for the race condition.

# TODO: Rx and Ry gates, not finished!
@i @inline function ROT_INSTRUCT(gate::Val{:Rx}, a, b, θ)
    ROT_INSTRUCT(Val(:Rz), a, b, π/2)
    ROT_INSTRUCT(Val(:Ry), a, b, θ)
    ROT_INSTRUCT(Val(:Rz), a, b, -π/2)
end

@i @inline function ROT_INSTRUCT(gate::Val{:Ry}, a, b, θ)
    divint(θ, 2)
    ROT(a, b, θ)
    mulint(θ, 2)
end


text:
xr = 56:76
yr = -162:-134
# xr = 20:30
# yr = -10:-5

function test(sx, sy, n)
  y_end = (2 * sy - n + 1) * n / 2
  if (n <= sx)
      x_end = (2 * sx - n + 1) * n / 2
  else
      x_end = (sx + 1) * sx / 2
  end
  if ((x_end in xr) && (y_end in yr))
    return true
  end
  return false
end

# Coarse but safe ranges for x/y velocities
max_sx = maximum(xr)
max_sy = max(abs(minimum(yr)), abs(maximum(yr)))
max_n = max_sy * 2 + 1

answers = 0
for sy in (-max_sy):max_sy
  for sx in 0:max_sx
    for n in 0:max_n
      if (test(sx, sy, n))
        # println(sx, ",", sy)
        global answers = answers + 1
        break
      end
    end
  end
end

println("answers found ", answers)


text:
module Day3

export day3

function check_slope(terrain, a, b)
    counter = 0
    n = length(terrain[1])
    for (x,row) in enumerate(terrain)
        (q,r) = divrem(a*(x-1), b)
        if r == 0
            y = q+1
            d = mod(y,n)
            y = d == 0 ? n : d
            if row[y] == '#'
                counter += 1
            end
        end
    end
   counter
end

function day3(::Val{:a}, filename)
    terrain = readlines(filename)
    check_slope(terrain, 3,1)
end

function day3(::Val{:b}, filename)
    terrain = readlines(filename)
    slopes = [(1,1), (3,1), (5,1), (7,1), (1,2)]
    result = UInt(1)
    for (a,b) in slopes
        result *= check_slope(terrain, a,b)
    end
    result
end

end


text:
using LinearAlgebra
using Polynomials
include("Modeling.jl")

mutable struct Controller
    C::Matrix{Float64}
    a::Matrix{Float64}
    A::Matrix{Float64}
    α::Matrix{Float64}
    K::Matrix{Float64}
    kr::Float64

    function get_C(ss::StateSpace)
	C = zeros(size(ss.A)...)
	for i in 1:length(ss.B)
	    C[:,i] = ss.A ^ (i - 1) * ss.B
	end
	if det(C) == 0
	    throw("NOT CONTROLLABLE")
	end
	return C
    end

    function get_a(A::Matrix{Float64})
	return Matrix(reverse(coeffs(fromroots(A)))[2:end]')
    end

    function get_A(a::Matrix{Float64})
	A = Matrix(1.0I, length(a), length(a))
	for i in 1:(length(a) - 1)
	    for j in (i + 1):length(a)
		A[i, j] = a[j - i]
	    end
	end
	return A
    end

    function Controller(ss::StateSpace)
	this = new()
	this.C = get_C(ss)
	this.a = get_a(ss.A)
	this.A = get_A(this.a)
	return this
    end

    function Controller(ss::StateSpace, poles::Vector)
	this = Controller(ss)
	set_poles!(this, ss, poles)
	return this
    end
end

function set_poles!(c::Controller, ss::StateSpace, poles::Vector{Float64})
    c.α = Matrix(reverse(coeffs(fromroots(poles)))[2:end]')
    c.K = (c.α - c.a) * inv(c.A) * inv(c.C)
    c.kr = -1 / (ss.C * inv(ss.A - ss.B * c.K) * ss.B)[1]
end

function set_poles!(c::Controller, ss::StateSpace, poles::Vector{ComplexF64})
    Δ = coeffs(fromroots(poles))
    if all(isreal, Δ)
	c.α = Matrix(reverse(Δ)[2:end]')
	c.K = (c.α - c.a) * inv(c.A) * inv(c.C)
	c.kr = -1 / (ss.C * inv(ss.A - ss.B * c.K) * ss.B)[1]
    else
	throw("COMPLEX POLES MUST BE CONJUGATE PAIRS")
    end
end

function get_u(c::Controller, x::Vector{Float64}, r::Float64)
    return c.kr * r - (c.K * x)[1]
end


text:
module TestInference

include("preamble.jl")
using Transducers: _nonidtype, DefaultInitOf, OptInit

collections = [
    1:1,
    1:0,
    Base.Generator(identity, 1:1),
    Base.Generator(identity, 1:0),
]

constant(x) = (_...) -> x

@testset "_nonidtype" begin
    # This is required for "foldl without init" tests below work:
    @test _nonidtype(Float64) === Float64
    @test _nonidtype(Union{DefaultInitOf{typeof(+)}, Float64}) === Float64
end

@testset "foldl" begin
    @testset for xs in collections
        @test_inferred foldl(+, Map(exp), xs; init=0.0)
        @test_inferred foldl(+, Map(exp), xs; init=OptInit)
        @test_inferred foldl(+, Map(exp) |> Filter(x -> x > 0), xs; init=0.0)
        @test_inferred foldl(+, Map(exp) |> Filter(x -> x > 0), xs; init=OptInit)
    end
    @testset for xs in [
        1:1,
        Base.Generator(identity, 1:1),
    ]
        @test_inferred foldl(+, Map(exp), xs)
        @test_inferred foldl(+, Map(exp) |> Filter(x -> x > 0), xs)
        @test_inferred foldl(+, Map(exp) |> Scan(+), xs)
        @test_inferred foldl(+, TakeLast(1), xs)
        @test_inferred foldl(+, PartitionBy(identity) |> Map(first), xs)
        @test_inferred foldl(+, Unique(), xs)
        @test_inferred foldl(right, ZipSource(Filter(isodd) |> Map(inc)), xs)

        # Nested stateful transducers.  (The ones with `right` and
        # `Map(x -> x::Int)` actually succeeded in some REPL
        # sessions...  Maybe some caching problem?)
        @test_broken_inferred foldl(right, Scan(+) |> Scan(+), xs)
        @test_broken_inferred foldl(*, Scan(+) |> Scan(+), xs)
        @test_broken_inferred foldl(*, Scan(+) |> Scan(+), xs; init=1)
        @test_broken_inferred foldl(*, Scan(+) |> Map(x -> x::Int) |> Scan(+), xs)
    end
    @test_inferred foldl(+, Cat(), [[1]])

    @testset "type constructor" begin
        xf = Map(Float64)
        @test_inferred foldl(+, xf, Int[1])

        xf = MapSplat(ComplexF64)
        @test_inferred foldl(+, xf, Tuple{Int, Int}[(1, 1)])
    end
end

@testset "collect" begin
    @testset for xs in collections
        @test_inferred collect(Map(exp), xs)
        @test_inferred collect(Map(exp) |> Filter(x -> x > 0), xs)
        @test_inferred collect(
            Enumerate() |> Map(last) |> Map(exp) |> Filter(x -> x > 0),
            xs,
        )
    end
end

hasnothing(xs) = foreach(Map(identity), xs) do x
    if x === nothing
        return reduced(Val(true))
    end
end |> ifunreduced() do _
    Val(false)
end

@testset "foreach" begin
    @testset for xs in collections
        @test_inferred foreach(constant(nothing), Map(exp), xs)
        @test_inferred foreach(constant(nothing), Map(exp) |> Filter(x -> x > 0),
                               xs)
    end
    @test (@inferred hasnothing((1, 2, 3))) === Val(false)
    @test (@inferred hasnothing((1, nothing, 3))) === Val(true)
end

@testset "eduction" begin
    @testset for xs in collections
        @test_inferred eduction(Map(exp), xs)
        @test_inferred eduction(Map(exp) |> Filter(x -> x > 0), xs)
    end
end

end  # module


text:
@testset "posterior" begin

end


text:
#' ## Генерация априорной информации

using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;

u_l(t) = -8 #+ cos(2*π * t);
u_r(t) =  4 #+ (1 + sin(2*π * t));
qf(x) = sin(3 * π * x);         # Коэффициент линейного усиления, который в обратной
                                # задаче необходимо определить, но при генерации априорной
                                # информации мы задаем некоторый коэффициент, который,
                                # собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.36;                # Область по T
N, M = 100, 80;                 # Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =     qf.(Xₙ);               # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ);               # Начальные условия
nothing #hide

u, XX, TP   = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);                              # Нам вернут решение, сетку, f_1^s
#
ϕl          = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                               # Левый вырожденный корень
ϕr          = phidetermination(qₙ, urₘ,Xₙ, N, Tₘ, M, reverseX = true);               # Нужно подать инвертированную сетку
ϕr          = reverse(ϕr, dims=1);                                                   # А после — инвертировать решение по X
ϕ           = Φ(ϕl, ϕr, N, M);                                                       # Серединный корень
f1_data     = f1(ϕ, u, Xₙ, N, M);                                                    # Положение переходного слоя
f2_data     = f2(f1_data, u, Xₙ, N, M);                                              # Значение функции на переходном слое
nothing #hide
#########################################################################################

q₀ = [ 0.0 * x for x in qₙ];
q₀ = q_guess;

S = 500;
β = 0.01;
qs, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data, S = S, β = β, w = 0.00001)

frames = [1:20; 21:div(S,50):S];    # Первые двадцать без пропусков
frames = collect(1:div(S, 50):S);   # С пропусками, чтобы всего было 50 кадров
make_minimzation_gif(Js, Qs, qₙ, Xₙ, name = "Minimization_uniform.gif", frames_to_write = frames, β = β)


text:
"""
    get_state_series(
        res::SimulationResults,
        ref::Tuple{String, Symbol}
    )
    end

Function to obtain series of states out of DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `ref:Tuple{String, Symbol}` : Tuple used to identify the dynamic device, via its name, as a `String`, and the associated state as a `Symbol`.
"""
function get_state_series(res::SimulationResults, ref::Tuple{String, Symbol})
    if allunique(res.solution.t)
        return res.solution.t, post_proc_state_series(res, ref)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], post_proc_state_series(res, ref)[ix]
    end
    return
end

"""
    get_voltage_magnitude_series(
        res::SimulationResults,
        bus_number::Int
    )

Function to obtain the voltage magnitude series out of the DAE Solution.

# Arguments:

- `res::SimulationResults` : Simulation Results object that contains the solution
- `bus_number::Int` : Bus number identifier
"""
function get_voltage_magnitude_series(res::SimulationResults, bus_number::Int)
    n_buses = get_bus_count(res)
    bus_ix = get(get_bus_lookup(res), bus_number, 0)
    V_R, V_I = post_proc_voltage_series(res.solution, bus_ix, n_buses)
    if allunique(res.solution.t)
        return res.solution.t, sqrt.(V_R .^ 2 + V_I .^ 2)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], sqrt.(V_R[ix] .^ 2 + V_I[ix] .^ 2)
    end
end

"""
    get_voltage_angle_series(
        res::SimulationResults,
        bus_number::Int
    )

Function to obtain the voltage angle series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `bus_number::Int` : Bus number identifier
"""
function get_voltage_angle_series(res::SimulationResults, bus_number::Int)
    n_buses = get_bus_count(res)
    bus_ix = get(get_bus_lookup(res), bus_number, 0)
    V_R, V_I = post_proc_voltage_series(res.solution, bus_ix, n_buses)
    if allunique(res.solution.t)
        return res.solution.t, atan.(V_I ./ V_R)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], atan.(V_I[ix] ./ V_R[ix])
    end
end

"""
    get_real_current_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the real current time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_real_current_series(res::SimulationResults, name::String)
    if allunique(res.solution.t)
        return res.solution.t, post_proc_real_current_series(res, name)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], post_proc_real_current_series(res, name)[ix]
    end
end

"""
    get_imaginary_current_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the imaginary current time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_imaginary_current_series(res::SimulationResults, name::String)
    if allunique(res.solution.t)
        return res.solution.t, post_proc_imaginary_current_series(res, name)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], post_proc_imaginary_current_series(res, name)[ix]
    end
end

"""
    get_activepower_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the active power output time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_activepower_series(res::SimulationResults, name::String)
    if allunique(res.solution.t)
        return res.solution.t, post_proc_activepower_series(res, name)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], post_proc_activepower_series(res, name)[ix]
    end
end

"""
    get_reactivepower_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the reactive power output time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_reactivepower_series(res::SimulationResults, name::String)
    if allunique(res.solution.t)
        return res.solution.t, post_proc_reactivepower_series(res, name)
    else
        @debug "found repeated time steps removing repetitions"
        ix = unique(i -> res.solution.t[i], eachindex(res.solution.t))
        return res.solution.t[ix], post_proc_reactivepower_series(res, name)[ix]
    end
end

"""
    show_states_initial_value(res::SimulationResults)

Function to print initial states.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
"""
function show_states_initial_value(res::SimulationResults)
    bus_size = get_bus_count(res)
    system = get_system(res)
    x0_init = res.solution.u[1]
    println("Voltage Variables")
    println("====================")
    buses_sorted =
        sort(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
    for bus in buses_sorted
        name = PSY.get_name(bus)
        println(name)
        println("====================")
        bus_n = PSY.get_number(bus)
        bus_ix = get_bus_lookup(res)[bus_n]
        V_R = x0_init[bus_ix]
        V_I = x0_init[bus_ix + bus_size]
        Vm = sqrt(V_R^2 + V_I^2)
        θ = angle(V_R + V_I * 1im)
        print("Vm ", round(Vm, digits = 4), "\n")
        print("θ ", round(θ, digits = 4), "\n")
        println("====================")
    end
    println("====================")
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        println("Differential States")
        println(name)
        println("====================")
        global_index = get_global_index(res)[name]
        for s in states
            print(s, " ", round(x0_init[global_index[s]], digits = 4), "\n")
        end
        println("====================")
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        println("====================")
        println("Line Current States")
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            printed_name = "Line " * name
            println("====================")
            println(printed_name)
            global_index = get_global_index(res)[name]
            for (i, s) in enumerate(states)
                print(s, " ", round(x0_init[global_index[s]], digits = 5), "\n")
            end
            println("====================")
        end
    end
    return
end

"""
    show_states_initial_value(sim::Simulation)

Function to print initial states.

# Arguments

- `sim::Simulation` : Simulation object that contains the initial condition
"""
function show_states_initial_value(sim::Simulation)
    inputs = sim.inputs
    bus_size = get_bus_count(inputs)
    system = get_system(sim)
    global_state_map = make_global_state_map(inputs)
    println("Voltage Variables")
    println("====================")
    buses_sorted =
        sort(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
    for bus in buses_sorted
        name = PSY.get_name(bus)
        println(name)
        println("====================")
        bus_n = PSY.get_number(bus)
        bus_ix = PSID.get_lookup(sim.inputs)[bus_n]
        V_R = sim.x0_init[bus_ix]
        V_I = sim.x0_init[bus_ix + bus_size]
        Vm = sqrt(V_R^2 + V_I^2)
        θ = angle(V_R + V_I * 1im)
        print("Vm ", round(Vm, digits = 4), "\n")
        print("θ ", round(θ, digits = 4), "\n")
        println("====================")
    end
    println("====================")
    for device in get_dynamic_injectors(inputs)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        println("Differential States")
        println(name)
        println("====================")
        global_index = global_state_map[name]
        for s in states
            print(s, " ", round(sim.x0_init[global_index[s]], digits = 4), "\n")
        end
        println("====================")
    end
    dyn_branches = get_dynamic_branches(inputs)
    if !isempty(dyn_branches)
        println("====================")
        println("Line Current States")
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            printed_name = "Line " * name
            println("====================")
            println(printed_name)
            global_index = global_state_map[name]
            x0_br = Dict{Symbol, Float64}()
            for (i, s) in enumerate(states)
                print(s, " ", round(sim.x0_init[global_index[s]], digits = 5), "\n")
            end
            println("====================")
        end
    end
    return
end


text:
# generate canonical hexagonal graph

"""
    $SIGNATURES

Generate a canonical factor graph: driving in a hexagonal circular pattern with one landmark.

Notes
- 7 Poses, :x0-:x6, Pose2,
- 1 Landmark, :l1, Point2,
- 6 Odometry, :x0x1f1, etc., Pose2Pose2 (Gaussian)
- 2 Sightings, :x0l1f1, :x6l1f1, RangeBearing (Gaussian)

Example
```julia
using RoME

fg = generateCanonicalFG_Hexagonal()
drawGraph(fg, show=true)
```

Related

[`generateCanonicalFG_Circle`](@ref), [`generateCanonicalFG_Kaess`](@ref), [`generateCanonicalFG_TwoPoseOdo`](@ref), [`generateCanonicalFG_Boxes2D!`](@ref)
"""
function generateCanonicalFG_Hexagonal(;fg::AbstractDFG=initfg(),
                                        N::Int=100,
                                        autoinit::Union{Bool, Nothing}=nothing,
                                        graphinit::Bool=true  )
  #
  getSolverParams(fg).N = N
  graphinit = if autoinit === nothing
    graphinit
  else
    @warn "autoinit is deprecated, use graphinit instead"
    autoinit
  end
  return generateCanonicalFG_Circle(6, graphinit=graphinit, landmark=true, loopClosure=true; fg=fg)
end


text:
@testset "Newton Trust Region" begin
@testset "Subproblems I" begin
    # verify that solve_tr_subproblem! finds the minimum
    n = 2
    gr = [-0.74637,0.52388]
    H = [0.945787 -3.07884; -3.07884 -1.27762]

    s = zeros(n)
    m, interior = Optim.solve_tr_subproblem!(gr, H, 1., s, max_iters=100)

    for j in 1:10
        bad_s = rand(n)
        bad_s ./= norm(bad_s)  # boundary
        model(s2) = (gr' * s2)[] + .5 * (s2' * H * s2)[]
        @test model(s) <= model(bad_s) + 1e-8
    end
end

@testset "Subproblems II" begin
    # random Hessians--verify that solve_tr_subproblem! finds the minimum
    for i in 1:10000
        n = rand(1:10)
        gr = randn(n)
        H = randn(n, n)
        H += H'

        s = zeros(n)
        m, interior = Optim.solve_tr_subproblem!(gr, H, 1., s, max_iters=100)

        model(s2) = (gr' * s2)[] + .5 * (s2' * H * s2)[]
        @test model(s) <= model(zeros(n)) + 1e-8  # origin

        for j in 1:10
            bad_s = rand(n)
            bad_s ./= norm(bad_s)  # boundary
            @test model(s) <= model(bad_s) + 1e-8
            bad_s .*= rand()  # interior
            @test model(s) <= model(bad_s) + 1e-8
        end
    end
end

@testset "Test problems" begin
    #######################################
    # First test the subproblem.
    Random.seed!(42)
    n = 5
    H = rand(n, n)
    H = H' * H + 4 * I
    H_eig = eigen(H)
    U = H_eig.vectors

    gr = zeros(n)
    gr[1] = 1.
    s = zeros(Float64, n)

    true_s = -H \ gr
    s_norm2 = dot(true_s, true_s)
    true_m = dot(true_s, gr) + 0.5 * dot(true_s, H * true_s)

    # An interior solution
    delta = sqrt(s_norm2) + 1.0
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @test interior
    @test !hard_case
    @test reached_solution
    @test abs(m - true_m) < 1e-12
    @test norm(s - true_s) < 1e-12
    @test abs(lambda) < 1e-12

    # A boundary solution
    delta = 0.5 * sqrt(s_norm2)
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @test !interior
    @test !hard_case
    @test reached_solution
    @test m > true_m
    @test abs(norm(s) - delta) < 1e-12
    @test lambda > 0

    # A "hard case" where the gradient is orthogonal to the lowest eigenvector

    # Test the checking
    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., 2., 3.], [0., 1., 1.])
    @test hard_case
    @test lambda_1_multiplicity == 1

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., 3.], [0., 0., 1.])
    @test hard_case
    @test lambda_1_multiplicity == 2

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., -1.], [0., 0., 0.])
    @test hard_case
    @test lambda_1_multiplicity == 3

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([1., 2., 3.], [0., 1., 1.])
    @test !hard_case

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., -1.], [0., 0., 1.])
    @test !hard_case

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., 2., 3.], [1., 1., 1.])
    @test !hard_case

    # Now check an actual had case problem
    L = fill(0.1, n)
    L[1] = -1.
    H = U * Matrix(Diagonal(L)) * U'
    H = 0.5 * (H' + H)
    @test issymmetric(H)
    gr = U[:,2][:]
    @test abs(dot(gr, U[:,1][:])) < 1e-12
    true_s = -H \ gr
    s_norm2 = dot(true_s, true_s)
    true_m = dot(true_s, gr) + 0.5 * dot(true_s, H * true_s)

    delta = 0.5 * sqrt(s_norm2)
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @test !interior
    @test hard_case
    @test reached_solution
    @test abs(lambda + L[1]) < 1e-4
    @test abs(norm(s) - delta) < 1e-12

    #######################################
    # Next, test on actual optimization problems.

    function f(x::Vector)
        (x[1] - 5.0)^4
    end

    function g!(storage::Vector, x::Vector)
        storage[1] = 4.0 * (x[1] - 5.0)^3
    end

    function h!(storage::Matrix, x::Vector)
        storage[1, 1] = 12.0 * (x[1] - 5.0)^2
    end

    d = TwiceDifferentiable(f, g!, h!, [0.0,])

    options = Optim.Options(store_trace = false, show_trace = false,
                            extended_trace = true)
    results = Optim.optimize(d, [0.0], NewtonTrustRegion(), options)
    @test_throws ErrorException Optim.x_trace(results)
    @test length(results.trace) == 0
    @test results.g_converged
    @test norm(Optim.minimizer(results) - [5.0]) < 0.01
    @test summary(results) == "Newton's Method (Trust Region)"

    eta = 0.9

    function f_2(x::Vector)
        0.5 * (x[1]^2 + eta * x[2]^2)
    end

    function g!_2(storage::Vector, x::Vector)
        storage[1] = x[1]
        storage[2] = eta * x[2]
    end

    function h!_2(storage::Matrix, x::Vector)
        storage[1, 1] = 1.0
        storage[1, 2] = 0.0
        storage[2, 1] = 0.0
        storage[2, 2] = eta
    end

    d = TwiceDifferentiable(f_2, g!_2, h!_2, Float64[127, 921])

    results = Optim.optimize(d, Float64[127, 921], NewtonTrustRegion())
    @test results.g_converged
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    # Test Optim.newton for all twice differentiable functions in
    # MultivariateProblems.UnconstrainedProblems.examples
    @testset "Optim problems" begin
        run_optim_tests(NewtonTrustRegion(); skip = ("Trigonometric", ),
                        show_name = debug_printing)
    end
end

@testset "PR #341" begin
    # verify that no PosDef exception is thrown
    Optim.solve_tr_subproblem!([0, 1.], [-1000 0; 0. -999], 1e-2, ones(2))
end

end


text:
@block AnshulSinghvi ["geo", "geomakie"] begin

    @cell "Geographical axes" [geoaxis, projection] begin
        using GeoMakie

        projections = [
            Projection("+proj=longlat") => "Standard lon-lat",
            Projection("+proj=robin")   => "Robinson",
            Projection("+proj=moll")    => "Molleweide",
            Projection("+proj=vandg")   => "van der Grinten I",
            Projection("+proj=wintri")  => "Winkel Tripel",
        ]

        geoaxis_scenes = Scene[]

        for (projection, titlestr) in projections
            push!(geoaxis_scenes, title(geoaxis(-180, 180, -90, 90; show_axis = false, crs = (dest = projection,), scale_plot = false), titlestr))
        end

        vbox(geoaxis_scenes...; parent = Scene(resolution = (1440, 200)))

    end

    @cell "Earth and Coastlines" [earth, coastlines, geoaxis] begin
        using GeoMakie

        coastlines(; crs = (dest = WinkelTripel(),), show_axis = false, scale_plot = false)
        earth!(; crs = (dest = WinkelTripel(),))
        geoaxis!(-180, 180, -90, 90; crs = (dest = WinkelTripel(),))
    end

    @cell "Air Particulates" [record, animation] begin
        using GeoMakie
        using GeoMakie: ImageMagick, Glob

        source = LonLat()
        dest = WinkelTripel()

        imgdir = observations("rgb/MYDAL2_M_AER_RA")

        imgpaths = sort(Glob.glob("*.PNG", imgdir)) # change this to your requirements

        img = ImageMagick.load(imgpaths[1])

        re = GeoMakie.date_regex("MYDAL2_M_AER_RA", "PNG")
        titletext = Node(join(match(re, basename(imgpaths[1])).captures, '-'))

        lons = LinRange(-179.5, 179.5, size(img)[2])
        lats = LinRange(89.5, -89.5, size(img)[1])

        xs = [lon for lat in lats, lon in lons]
        ys = [lat for lat in lats, lon in lons]

        Proj4.transform!(source, dest, vec(xs), vec(ys))

        scene = surface(xs, ys, zeros(size(xs)); color = img, shading = false, show_axis = false, scale_plot = false)
        surfplot = scene.plots[end]

        geoaxis!(scene, -180, 180, -90, 90; crs = (src = source, dest = dest,));

        titletext = Node("07/2016")

        fullsc = title(scene, titletext; fontsize = 40);

        record(fullsc, @replace_with_a_path(mp4), imgpaths; framerate = 10) do img
            year, month = match(re, img).captures
            surfplot.color = ImageMagick.load(img)
            titletext[] = "$month/$year"
        end

    end

end


text:
#!usr/bin/env julia
#testing version: 0.6.0

importall Base

struct Physical{T, m, kg, s, A, K, mol, cd}
    value :: T
end

Physical(x, m, kg, s, A, K, mol, cd) = Physical{typeof(x), m, kg, s, A, K, mol, cd}(x)
Physical(x; m=0, kg=0, s=0, A=0, K=0, mol=0, cd=0) = Physical{typeof(x), m, kg, s, A, K, mol, cd}(x)

@inline character_lift(x) = Dict(
    '0' => '⁰',
    '1' => '¹',
    '2' => '²',
    '3' => '³',
    '4' => '⁴',
    '5' => '⁵',
    '6' => '⁶',
    '7' => '⁷',
    '8' => '⁸',
    '9' => '⁹',
)[x]

function show(io::IO, x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd}
    print(io, "$(x.value)")
    for (u, v) in [(kg, :kg), (A, :A), (m, :m), (K, :K), (mol, :mol), (cd, :cd), (s, :s)]
        if u != 0
            print(io, "⋅", v)
            if u < 0
                print(io, '⁻')
            end
            if u != 1
                print(io, map(character_lift, collect("$(abs(u))")))
            end
        end
    end
end

#addition
(+)(x::Physical{T1, m, kg, s, A, K, mol, cd}, y::Physical{T2, m, kg, s, A, K, mol, cd}) where {T1, T2, m, kg, s, A, K, mol, cd} =
    Physical(x.value + y.value, m, kg, s, A, K, mol, cd)

#negate
(-)(x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd} =
    Physical(x.value - y.value, m, kg, s, A, K, mol, cd)

#substration
(-)(x::Physical{T1, m, kg, s, A, K, mol, cd}, y::Physical{T2, m, kg, s, A, K, mol, cd}) where {T1, T2, m, kg, s, A, K, mol, cd} =
    Physical(x.value - y.value, m, kg, s, A, K, mol, cd)

#multiplication with pure number
(*)(k, x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd} =
    Physical{T, m, kg, s, A, K, mol, cd}(k * x.value)

(*)(x::Physical{T, m, kg, s, A, K, mol, cd}, k) where {T, m, kg, s, A, K, mol, cd} =
    Physical{T, m, kg, s, A, K, mol, cd}(x.value * k)

#general multiplication
(*)(x::Physical{T1, m1, kg1, s1, A1, K1, mol1, cd1}, y::Physical{T2, m2, kg2, s2, A2, K2, mol2, cd2}) where {T1, T2, m1, kg1, s1, A1, K1, mol1, cd1, m2, kg2, s2, A2, K2, mol2, cd2} =
    Physical(x.value * y.value, m1+m2, kg1+kg2, s1+s2, A1+A2, K1+K2, mol1+mol2, cd1+cd2)

#division with a pure number
(/)(k, x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd} =
    Physical{T, -m, -kg, -s, -A, -K, -mol, -cd}(k / x.value)

(/)(x::Physical{T, m, kg, s, A, K, mol, cd}, k) where {T, m, kg, s, A, K, mol, cd} =
    Physical{T, m, kg, s, A, K, mol, cd}(x.value * k)

#general division
(/)(x::Physical{T1, m1, kg1, s1, A1, K1, mol1, cd1}, y::Physical{T2, m2, kg2, s2, A2, K2, mol2, cd2}) where {T1, T2, m1, kg1, s1, A1, K1, mol1, cd1, m2, kg2, s2, A2, K2, mol2, cd2} =
    Physical(x.value / y.value, m1-m2, kg1-kg2, s1-s2, A1-A2, K1-K2, mol1-mol2, cd1-cd2)

#power and square root
(^)(x::Physical{T, m, kg, s, A, K, mol, cd}, k) where {T, m, kg, s, A, K, mol, cd} =
    Physical(x.value ^ k, m*k, kg*k, s*k, A*k, K*k, mol*k, cd*k)

sqrt(x::Physical{T, m, kg, s, A, K, mol, cd}) where {T, m, kg, s, A, K, mol, cd} = let k = 1//2
    Physical(sqrt(x.value), m*k, kg*k, s*k, A*k, K*k, mol*k, cd*k)
end

#exponent and logrithm
#warning: all units are dropped. Probably I should come up with a better solution later
(^)(x, y::Physical) = x.value ^ k
(^)(x::Physical, y::Physical) = x ^ y.value

log(x::Physical) = log(x.value)
log2(x::Physical) = log2(x.value)
log10(x::Physical) = log10(x.value)
log(b, m) = log(m) - log(b)

#handling ambigious
(^)(x::Physical{T, m, kg, s, A, K, mol, cd}, k::Integer) where {T, m, kg, s, A, K, mol, cd} =
    Physical(x.value ^ k, m*k, kg*k, s*k, A*k, K*k, mol*k, cd*k)

(+)(x::Physical, y::Physical) = begin
    warn("Type mismatch while adding")
    show(x)
    print(" and ")
    show(y)
    x.value + y.value
end

meter = Physical(1, m=1)
kilogram = Physical(1, kg=1)
second = Physical(1, s=1)
ampere = Physical(1, A=1)
kelvin = Physical(1, K=1)
mole = Physical(1, mol=1)
candela = Physical(1, cd=1)
hertz = 1/second
newton = kilogram*meter/second^2
pascal = newton/meter^2
joule = newton*meter
watt = joule/second
coulomb = ampere * second
volt = watt/ampere
farad = coulumb/volt
ohm = volt/ampere
weber = volt*second
tesla = weber/meter^2
henry = weber/ampere
lumen = candela
lux = lumen/meter^2

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
using OSF
using Test

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


text:
# ORDER BY clause.

mutable struct OrderClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    by::Vector{SQLClause}

    OrderClause(;
                over = nothing,
                by = SQLClause[]) =
        new(over, by)
end

OrderClause(by...; over = nothing) =
    OrderClause(over = over, by = SQLClause[by...])

"""
    ORDER(; over = nothing, by = [])
    ORDER(by...; over = nothing)

A `ORDER BY` clause.

# Examples

```jldoctest
julia> c = FROM(:person) |>
           ORDER(:year_of_birth) |>
           SELECT(:person_id);

julia> print(render(c))
SELECT "person_id"
FROM "person"
ORDER BY "year_of_birth"
```
"""
ORDER(args...; kws...) =
    OrderClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(ORDER), pats::Vector{Any}) =
    dissect(scr, OrderClause, pats)

function PrettyPrinting.quoteof(c::OrderClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(ORDER))
    append!(ex.args, quoteof(c.by, ctx))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::OrderClause, c′) =
    OrderClause(over = rebase(c.over, c′), by = c.by)


text:
# STOCHASTIC APPROXIMATION ALGORITHM FOR SOLVING DCOPF WITH RESERVE SATURATION

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi

# generate and store Gurobi environment once and for all
const gurobi_env = Gurobi.Env()

tic() # overall time

#*========= MODEL INFORMATION ==========
const model = "6bus"
const modelFile = model * ".jl"
include(modelFile)
#*======================================

#*========= INCLUDE SUPPORTING FILES ==============
include("utilities.jl")
include("modelFns.jl")
include("solveRecourse.jl")
include("evaluateSoln.jl")
include("projection.jl")
include("stochgrad.jl")
include("stepLength.jl")
include("initialguess_stoch.jl")

#*============ LINE PENALTY ==========
# set line penalty coefficients to get comparable results
const gamma_line_val = logspace(1.0,5.0,17)
# finer discretization for the no wind reserves case to try and get better approximation of EF
if(caseID == "no_wind_reserves")
	gamma_line_val1 = logspace(1.0,3.0,9)
	gamma_line_val2 = logspace(3.25,3.5,6)
	gamma_line_val3 = logspace(3.75,5.0,6)
	gamma_line_val = [gamma_line_val1;gamma_line_val2;gamma_line_val3]
end

# set seed for reproducibility
srand(1135)

# directory name for storing results
const baseDirName = "C:/Users/rkannan/Desktop/DCOPF-reserve-saturation/experiments/" * model * "/" * "stochapprox_" * caseID
mkpath(baseDirName)

#*========= PRINT OPTIONS ====================
const storeResults = true
const printIterationStats = false
const printFinalSolution = false

const initFile = "initial_guess.txt"
const solnFile = "solution.txt"

const costFile = "expected_cost.txt"
const costPowerFile = "expected_cost_power.txt"
const costReservesFile = "expected_cost_reserves.txt"
const sumLineViolFile = "sumLineViol.txt"
const lineViolProbFile = "lineViolProb.txt"

const avgLineViolFile = "avgAbsLineFlowViol.txt"
const timeFile = "solutionTime.txt"

const bestSolnFile = "best_solution.txt"
const p0BestSolnFile = "p0_best_solution.txt"
const rpBestSolnFile = "rp_best_solution.txt"
const rmBestSolnFile = "rm_best_solution.txt"
const alphaBestSolnFile = "alpha_best_solution.txt"
const bestCostFile = "best_expected_cost.txt"
const bestCostPowerFile = "best_expected_cost_power.txt"
const bestSumLineViolFile = "best_sumLineViol.txt"
const bestLineViolProbFile = "best_lineViolProb.txt"

const bestCostReservesFile = "best_expected_cost_reserves.txt"
const bestCostPosResPenFile = "best_expected_cost_pos_res_pen.txt"
const bestCostNegResPenFile = "best_expected_cost_neg_res_pen.txt"
const powerGenSolnFile = "best_expected_power_gen.txt"
const windSpillSolnFile = "best_expected_wind_spill.txt"

#*========== ALGORITHMIC PARAMETERS ==========
# mini-batch size (m)
const minibatchSize = 20

# max. number of iterations for each replicate
const maxNumIterations = 1000

# minimum and maximum number of replicates (S)
const minNumReplicates = 20
const maxNumReplicates = 1000

# number of trials of the algorithm
const numTrials = 5
if(model == "118bus")
	numTrials = 1
end

# improvement requirement for considering another replicate
const impReqForTerm = 1E-05
const termCheckPeriod = 10

# step size modification parameters
const stepLengthDecrCheckIter = 3
const stepLengthIncrCheckIter = 3
const impReqForIncrStepSize = 1E-03
const impReqForDecrStepSize = 1E-03
const stepLengthIncrFactor = 2.0
const stepLengthDecrFactor = 2.0

#*===========================================

# number of line penalty parameter iterations
const numGammaLineIter = size(gamma_line_val,1)

println("")
println("***** STARTING THE STOCHASTIC APPROXIMATION METHOD FOR MODEL ",model," AND CASE ",caseID)
println("")

# run trials/replications
for trial = 1:numTrials

	println("\n******************************")
	println("******************************")
	println("TRIAL ",trial,"/",numTrials)
	println("******************************")

	tic() # trial time

	# sweep through the line limit penalty parameters
	for iter_line = 1:numGammaLineIter

		tic() # line penalty iteration time

		gamma_line_iter::Float64 = gamma_line_val[iter_line]

		println("\n==============================")
		println("==============================")
		println("LINE PENALTY ITER ",iter_line,"/",numGammaLineIter," with penalty ",gamma_line_iter)
		println("==============================")

		tic() # preprocessing time

		tic()

		#*========== COMPUTE INITIAL GUESS ===========
		cost_init2::Float64 = Inf
		cost_power_init2::Float64 = Inf
		probAnyLineFlowViol_init2::Float64 = 1.0
		sumLineFlowViol_init2::Float64 = Inf
		cost_reserves_init2::Float64 = Inf

		# in subsequent iterations, try and use the best solution from the previous iteration
		if(iter_line > 1)
			# project initial guess onto feasible set
			p0_initguess3, rp_initguess3, rm_initguess3, alpha_initguess3, ~ = projectOntoFeasibleSet(p0_init,rp_init,rm_init,alpha_init)

			cost_init2, cost_power_init2, probAnyLineFlowViol_init2, sumLineFlowViol_init2, cost_reserves_init2, ~, ~, ~, ~, ~, ~, ~ = estimateCostOfSoln(p0_initguess3,rp_initguess3,rm_initguess3,alpha_initguess3,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line,false,false)
		end

		# compute initial guess for the first-stage variables using the GP model
		p0_initguess, rp_initguess, rm_initguess, alpha_initguess, ig_solve_status = getInitialGuess(gamma_line_iter,delta_line)

		if(ig_solve_status != :Optimal && ig_solve_status != :Suboptimal)
			println("--- Gurobi did not solve IG model to optimality! Going to try IPOPT ---")
			p0_initguess, rp_initguess, rm_initguess, alpha_initguess, ig_solve_status = getInitialGuess(gamma_line_iter,delta_line,false)
			if(ig_solve_status != :Optimal && ig_solve_status != :Suboptimal)
				println("--- Ipopt also did not solve IG model to optimality! Going to use a trivial initial guess ---")
				p0_initguess = zeros(Float64,numGen)
				rp_initguess = zeros(Float64,numGen)
				rm_initguess = zeros(Float64,numGen)
				alpha_initguess = ones(Float64,numGen)/numGen
			end
		end

		# project initial guess onto feasible set
		p0_initguess2, rp_initguess2, rm_initguess2, alpha_initguess2, ~ = projectOntoFeasibleSet(p0_initguess,rp_initguess,rm_initguess,alpha_initguess)

		#*========== EVALUATE QUALITY OF INITIAL GUESS ===========

		IGSolnTime::Float64 = toq()

		tic() # cost time

		cost_init::Float64, cost_power_init::Float64, probAnyLineFlowViol_init::Float64, sumLineFlowViol_init::Float64, cost_reserves_init::Float64, ~, ~, ~, ~, ~, ~, ~ = estimateCostOfSoln(p0_initguess2,rp_initguess2,rm_initguess2,alpha_initguess2,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line,false,false)

		if(cost_init < cost_init2)
			for i = 1:numGen
				p0_init[i] = p0_initguess2[i]
				rp_init[i] = rp_initguess2[i]
				rm_init[i] = rm_initguess2[i]
				alpha_init[i] = alpha_initguess2[i]
			end
		else
			for i = 1:numGen
				p0_init[i] = p0_initguess3[i]
				rp_init[i] = rp_initguess3[i]
				rm_init[i] = rm_initguess3[i]
				alpha_init[i] = alpha_initguess3[i]
			end
			cost_init = cost_init2
			cost_power_init = cost_power_init2
			probAnyLineFlowViol_init = probAnyLineFlowViol_init2
			sumLineFlowViol_init = sumLineFlowViol_init2
			cost_reserves_init = cost_reserves_init2
		end

		println("\n*** INITIAL GUESS ***")

		println("Expected cost (incl. all penalties): ",round.(cost_init,displayPrecision))
		println("Expected cost without line penalty costs: ",round.(cost_power_init,displayPrecision))
		println("Expected cost of reserves (incl. reserve penalties): ",round.(cost_reserves_init,displayPrecision))
		println("Expected sum of line flow violations: ",round.(sumLineFlowViol_init,displayPrecision))
		println("Probability that ANY line flow is violated: ",round.(probAnyLineFlowViol_init,displayPrecision),"\n")

		costTime::Float64 = toq()
		println("Time for evaluating initial guess: ",round.(IGSolnTime,2),"  time for evaluating quality: ",round.(costTime,2))

		tic() # step length time

		#*============ ESTIMATE STEP LENGTH ==========

		numSamplesForLipEst::Int64 = 200

		stepLength = estimateStepLength(p0_initguess2,rp_initguess2,rm_initguess2,alpha_initguess2,minibatchSize,numSamplesForLipEst,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line)

		initialStepLength = 1.0*stepLength

		#*============================================

		stepLengthTime = toq()
		preprocessingTime = toq()

		println("Step length: ",stepLength,"  time for estimating it: ",round.(stepLengthTime,2))
		println("Total preprocessing time: ",round.(preprocessingTime,2),"\n")

		#*========= STORE RESULTS ====================
		if(storeResults)

			subDirName = baseDirName * "/" * "/rep" * string(trial) * "/lp" * string(iter_line) * "/"
			mkpath(subDirName)

			# write details to text file, including some key details about the test instance
			details_file = subDirName * model * ".txt"
			open(details_file, "w") do f
				write(f,"Model: $model \n")
				write(f,"caseID: $caseID \n")
				write(f,"Mini-batch size (m): $minibatchSize \n")
				write(f,"Max. number of iterations (N): $maxNumIterations \n")
				write(f,"Min./Max. number of replicates: $minNumReplicates/$maxNumReplicates \n")
				write(f,"improvement requirement for termination: $impReqForTerm \n")
				write(f,"improvement requirement for increasing stepLength: $impReqForIncrStepSize \n")
				write(f,"non-improvement requirement for decreasing stepLength: $impReqForDecrStepSize \n")
				write(f,"check for decreasing step length every $stepLengthDecrCheckIter iterations \n")
				write(f,"check for increasing step length every $stepLengthIncrCheckIter iterations \n")
				write(f,"initial step length: $stepLength \n")
				write(f,"step length increase factor: $stepLengthIncrFactor \n")
				write(f,"step length decrease factor: $stepLengthDecrFactor \n")
				write(f,"number of samples for estimating Lipschitz constant and variance: $numSamplesForLipEst \n")
				write(f,"Preprocessing time: $preprocessingTime \n\n")

				write(f,"======================\n")
				write(f,"Line flow violation penalty parameters \n")
				write(f,"----------------------\n")
				write(f,"Penalty factors (gamma_line): $gamma_line_val \n")
				write(f,"Beyond what fraction to penalize? (delta_line): $delta_line \n\n")

				write(f,"======================\n")
				write(f,"Reserves parameters \n")
				write(f,"----------------------\n")
				write(f,"How much costlier are reserves compared to normal generation? (reservesCostFactor): $reservesCostFactor \n")
				write(f,"What is the penalty factor for exceeding reserves? (gamma_res): $gamma_res \n")
				write(f,"How much costlier are wind reserves compared to the cheapest regular generator? (windReservesCostFactor): $windReservesCostFactor \n")
				write(f,"\n\n")

				write(f,"**********************\n")
				write(f,"SOME PARAMETERS OF THE INPUT MODEL \n\n")
				write(f,"**********************\n")

				write(f,"======================\n")
				write(f,"Generator bounds \n")
				write(f,"----------------------\n")
				write(f,"Regular generator upper limits: $regGenMaxPow \n\n")

				write(f,"======================\n")
				write(f,"Generator cost parameters \n")
				write(f,"----------------------\n")
				write(f,"Regular generator linear cost coefficients: $regGenCost \n")
				write(f,"Regular generator positive reserves cost coefficients: $regGenPosResCost \n")
				write(f,"Regular generator negative reserves cost coefficients: $regGenNegResCost \n")
				write(f,"Wind generator positive reserves cost coefficients: $windGenPosResCost \n")
				write(f,"Wind generator negative reserves cost coefficients: $windGenNegResCost \n\n")

				write(f,"======================\n")
				write(f,"Line flow parameters \n")
				write(f,"----------------------\n")
				write(f,"Line flow factor: $lineFlowFactor \n")
				write(f,"Line flow limits: $flow_max \n\n")

				write(f,"======================\n")
				write(f,"Demand parameters \n")
				write(f,"----------------------\n")
				write(f,"Average loads/demands: $meanDemand \n")
				write(f,"Standard deviation of loads/demands: $demandStdev \n")
				write(f,"Demand covariance matrix: $demandCovMat \n\n")

				write(f,"======================\n")
				write(f,"Wind parameters \n")
				write(f,"----------------------\n")
				write(f,"Wind factor: $windFactor \n")
				write(f,"Average wind output: $windPredMean \n")
				write(f,"Wind prediction stdev factor: $windPredStdevFactor \n")
				write(f,"Standard deviation of wind output: $windPredStdev \n")
				write(f,"Do the wind generators provide reserves?: $windGenRes \n")
			end

			# write data to text file
			init_file = subDirName * initFile
			open(init_file, "w") do f
				write(f,"NOTE: regular generator variables come first \n\n")
				write(f,"Expected cost: $cost_init \n")
				write(f,"Expected cost of generation (without line penalty costs): $cost_power_init \n")
				write(f,"Expected sum of line flow limit violations: $sumLineFlowViol_init \n")
				write(f,"Probability that ANY line limit is violated: $probAnyLineFlowViol_init \n\n")
				write(f,"p0: $p0_init \n")
				write(f,"rp: $rp_init \n")
				write(f,"rm: $rm_init \n")
				write(f,"alpha: $alpha_init \n\n")
				write(f,"Time for computing initial guess: ",IGSolnTime)
				write(f,"Time for evaluating quality of initial guess: ",costTime)
			end

			# write the line penalty coefficients to a file
			gamma_line_file = subDirName * "gamma_line.txt"
			open(gamma_line_file,"w") do f
				for iter_line2 = 1:numGammaLineIter
					write(f,"$(gamma_line_val[iter_line2])\n")
				end
			end

			# write data to text file
			cost_file = subDirName * costFile
			open(cost_file, "w") do f
			end

			cost_power_file = subDirName * costPowerFile
			open(cost_power_file, "w") do f
			end

			cost_reserves_file = subDirName * costReservesFile
			open(cost_reserves_file, "w") do f
			end

			sum_viol_file = subDirName * sumLineViolFile
			open(sum_viol_file, "w") do f
			end

			viol_prob_file = subDirName * lineViolProbFile
			open(viol_prob_file, "w") do f
			end

			best_cost_file = subDirName * bestCostFile
			open(best_cost_file, "w") do f
			end

			# write data to text file
			best_cost_power_file = subDirName * bestCostPowerFile
			open(best_cost_power_file, "w") do f
			end

			best_sum_viol_file = subDirName * bestSumLineViolFile
			open(best_sum_viol_file, "w") do f
			end

			best_viol_prob_file = subDirName * bestLineViolProbFile
			open(best_viol_prob_file, "w") do f
			end

			avg_viol_file = subDirName * avgLineViolFile
			open(avg_viol_file, "w") do f
			end

			soln_file = subDirName * solnFile
			open(soln_file, "w") do f
			end

			best_soln_file = subDirName * bestSolnFile
			open(best_soln_file, "w") do f
			end

			p0_best_soln_file = subDirName * p0BestSolnFile
			open(p0_best_soln_file, "w") do f
			end

			rp_best_soln_file = subDirName * rpBestSolnFile
			open(rp_best_soln_file, "w") do f
			end

			rm_best_soln_file = subDirName * rmBestSolnFile
			open(rm_best_soln_file, "w") do f
			end

			alpha_best_soln_file = subDirName * alphaBestSolnFile
			open(alpha_best_soln_file, "w") do f
			end

			best_cost_reserves_file = subDirName * bestCostReservesFile
			open(best_cost_reserves_file, "w") do f
			end

			best_cost_pos_res_file = subDirName * bestCostPosResPenFile
			open(best_cost_pos_res_file, "w") do f
			end

			best_cost_neg_res_file = subDirName * bestCostNegResPenFile
			open(best_cost_neg_res_file, "w") do f
			end

			power_gen_soln_file = subDirName * powerGenSolnFile
			open(power_gen_soln_file, "w") do f
			end

			wind_spill_soln_file = subDirName * windSpillSolnFile
			open(wind_spill_soln_file, "w") do f
			end

			time_file = subDirName * timeFile
			open(time_file, "w") do f
			end

		end
		#*============================================

		#*========= STORE BEST FOUND SOLUTION ========

		# store the best found solution in terms of expected cost
		# initialize using the initial guess solution
		best_iter::Int64 = 0
		best_cost::Float64 = cost_init
		best_cost_power::Float64 = cost_power_init
		best_lineViolProb::Float64 = probAnyLineFlowViol_init
		best_sumLineViol::Float64 = sumLineFlowViol_init

		p0_best_solution = zeros(Float64,numGen)
		rp_best_solution = zeros(Float64,numGen)
		rm_best_solution = zeros(Float64,numGen)
		alpha_best_solution = zeros(Float64,numGen)

		for i = 1:numGen
			p0_best_solution[i] = p0_init[i]
			rp_best_solution[i] = rp_init[i]
			rm_best_solution[i] = rm_init[i]
			alpha_best_solution[i] = alpha_init[i]
		end

		# store the final iterate and objective for each replicate
		cost_soln = Inf*ones(Float64,maxNumReplicates)
		cost_power_soln = Inf*ones(Float64,maxNumReplicates)
		anyLineViolProb_soln = Inf*ones(Float64,maxNumReplicates)
		sumLineViol_soln = Inf*ones(Float64,maxNumReplicates)

		p0_solution = zeros(Float64,maxNumReplicates,numGen)
		rp_solution = zeros(Float64,maxNumReplicates,numGen)
		rm_solution = zeros(Float64,maxNumReplicates,numGen)
		alpha_solution = zeros(Float64,maxNumReplicates,numGen)

		numStepSizeIncr::Int64 = 0
		numStepSizeDecr::Int64 = 0

		local_best_cost::Float64 = best_cost
		local_best_cost2::Float64 = best_cost
		local_best_cost3::Float64 = best_cost

		actualNumIter::Int64 = 0

		#*========= OPTIMIZATION LOOP ==========
		for iter_S = 1:maxNumReplicates

			tic()
			objTime::Float64 = 0
			gradTime::Float64 = 0
			projTime::Float64 = 0

			if(printIterationStats)
				println("=====================================")
				println("Optimization iteration # ", iter_S)
			end

			# determine actual number of SGD iterations
			uniformRand = rand(1)
			numIterations::Int64 = ceil.(uniformRand[1]*maxNumIterations)

			if(printIterationStats)
				println("Number of iterations: ",numIterations)
			end

			#*========= SET INITIAL GUESS ==========
			p0 = zeros(Float64,numGen)
			rp = zeros(Float64,numGen)
			rm = zeros(Float64,numGen)
			alpha = zeros(Float64,numGen)

			# if first replicate, use user-specified initial guess; otherwise, use solution from previous iteration
			if(iter_S == 1)
				for i = 1:numGen
					p0[i] = p0_init[i]
					rp[i] = rp_init[i]
					rm[i] = rm_init[i]
					alpha[i] = alpha_init[i]
				end
			else
				for i = 1:numGen
					p0[i] = p0_solution[iter_S-1,i]
					rp[i] = rp_solution[iter_S-1,i]
					rm[i] = rm_solution[iter_S-1,i]
					alpha[i] = alpha_solution[iter_S-1,i]
				end
			end
			#*======================================

			#*===========================================================
			#* CORE STOCHASTIC APPROXIMATION ALGORITHM IN THIS BLOCK
			#*===========================================================
			for t = 1:numIterations

				tic()

				#*==== COMPUTE STOCHASTIC GRADIENT =====
				grad_p0, grad_rp, grad_rm, grad_alpha = getStochasticGradient(p0,rp,rm,alpha,minibatchSize,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line)
				#*======================================

				gradTime += toq()

				#*==== SGD STEP =====
				p0 = p0 - stepLength[1]*grad_p0
				rp = rp - stepLength[2]*grad_rp
				rm = rm - stepLength[3]*grad_rm
				alpha = alpha - stepLength[4]*grad_alpha
				#*===================

				tic()

				#*==== PROJECTION STEP =====
				p0_2, rp_2, rm_2, alpha_2, ~ = projectOntoFeasibleSet(p0,rp,rm,alpha)
				#*==========================

				for i = 1:numGen
					p0[i] = p0_2[i]
					rp[i] = rp_2[i]
					rm[i] = rm_2[i]
					alpha[i] = alpha_2[i]
				end

				projTime += toq()

			end # SGD iterations
			#*========= END SGD LOOP ==========

			actualNumIter = iter_S

			tic()

			#*==== ESTIMATE COST AND LINE VIOLATION PROB =====
			cost_soln[iter_S], cost_power_soln[iter_S], anyLineViolProb_soln[iter_S], sumLineViol_soln[iter_S], cost_reserves_soln::Float64, ~, ~, ~, ~, ~, avgLineFlowViol, maxLineFlowViol = estimateCostOfSoln(p0,rp,rm,alpha,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line,false,false)
			#*====================================

			objTime += toq()

			# update solution for this replicate
			for j = 1:numGen
				p0_solution[iter_S,j] = p0[j]
				rp_solution[iter_S,j] = rp[j]
				rm_solution[iter_S,j] = rm[j]
				alpha_solution[iter_S,j] = alpha[j]
			end

			# update best found solution
			if(cost_soln[iter_S] < best_cost)
				best_iter = iter_S
				best_cost = cost_soln[iter_S]
				best_cost_power = cost_power_soln[iter_S]
				best_lineViolProb = anyLineViolProb_soln[iter_S]
				best_sumLineViol = sumLineViol_soln[iter_S]
				for j = 1:numGen
					p0_best_solution[j] = p0[j]
					rp_best_solution[j] = rp[j]
					rm_best_solution[j] = rm[j]
					alpha_best_solution[j] = alpha[j]
				end
			end

			solutionTime = toq()

			#*========= STORE RESULTS ==========
			if(storeResults)

				open(cost_file, "a") do f
					write(f,"$(cost_soln[iter_S]) \n")
				end

				open(cost_power_file, "a") do f
					write(f,"$(cost_power_soln[iter_S]) \n")
				end

				open(cost_reserves_file, "a") do f
					write(f,"$cost_reserves_soln \n")
				end

				open(sum_viol_file, "a") do f
					write(f,"$(sumLineViol_soln[iter_S]) \n")
				end

				open(viol_prob_file, "a") do f
					write(f,"$(anyLineViolProb_soln[iter_S]) \n")
				end

				open(avg_viol_file, "a") do f
					write(f,"ITER NUM. $iter_S \n")
					for k = 1:numLines
						write(f,"($(lines[k][1]),$(lines[k][2])): $(avgLineFlowViol[k])/$(maxLineFlowViol[k])  \n")
					end
					write(f,"\n")
				end

				open(soln_file, "a") do f
					write(f,"ITER NUM. $iter_S \n")
					write(f,"Expected cost: $(cost_soln[iter_S]) \n")
					write(f,"Expected cost without line penalty costs: $(cost_power_soln[iter_S]) \n")
					write(f,"Expected sum of line flow limit violations: $(sumLineViol_soln[iter_S]) \n")
					write(f,"ANY line limit violation prob: $(anyLineViolProb_soln[iter_S]) \n")
					write(f,"p0: $(p0_solution[iter_S,:]) \n")
					write(f,"rp: $(rp_solution[iter_S,:]) \n")
					write(f,"rm: $(rm_solution[iter_S,:]) \n")
					write(f,"alpha: $(alpha_solution[iter_S,:]) \n")
					write(f,"\n\n")
				end

				open(time_file, "a") do f
					write(f,"$solutionTime \n")
				end
			end
			#*==================================

			if(printIterationStats)
				println("EXPECTED VALUES:  COST (incl. pen.): ",round.(cost_soln[iter_S],displayPrecision),"  COST OF POWER ALONE (incl. reserve pen.): ",round.(cost_power_soln[iter_S],displayPrecision),"  SUM OF LINE VIOL: ",round.(sumLineViol_soln[iter_S],displayPrecision),"  PROB ANY LINE LIMIT VIOLATED: ",round.(anyLineViolProb_soln[iter_S],displayPrecision))

				otherTime::Float64 = solutionTime - gradTime - projTime - objTime

				println("SOLN TIME: ",round.(solutionTime,2)," GRAD TIME: ",round.(gradTime,2)," PROJ TIME: ",round.(projTime,2)," OBJ TIME: ", round.(objTime,2)," OTHER TIME: ", round.(otherTime,2))
			end

			#*========= HEURISTICS FOR UPDATING STEP LENGTH ==========

			# if the min. number of replicates are done, check if termination criteria are satisfied
			if(iter_S >= minNumReplicates)
				improvements = zeros(Float64,termCheckPeriod)
				for j = 1:termCheckPeriod
					improvements[j] = local_best_cost - cost_soln[iter_S+1-j]
				end

				localImpReq::Float64 = impReqForTerm*local_best_cost
				max_improvement::Float64 = maximum(improvements)

				if(max_improvement < localImpReq)
					break
				end
			end

			# once every stepLengthDecrCheckIter iterations, see if the step length can be decreased
			if(rem(iter_S,stepLengthDecrCheckIter) == 0)
				improvements = zeros(Float64,stepLengthDecrCheckIter)
				for j = 1:stepLengthDecrCheckIter
					improvements[j] = local_best_cost2 - cost_soln[iter_S+1-j]
				end

				localImpReq3::Float64 = impReqForDecrStepSize*local_best_cost2

				max_improvement = maximum(improvements)
				if(max_improvement < -localImpReq3)
					stepLength /= stepLengthDecrFactor
					numStepSizeDecr += 1
				end
			end

			# once every stepLengthIncrCheckIter iterations, see if the step length can be increased
			if(rem(iter_S,stepLengthIncrCheckIter) == 0)
				improvements = zeros(Float64,stepLengthIncrCheckIter)
				for j = 1:stepLengthIncrCheckIter
					improvements[j] = local_best_cost2 - cost_soln[iter_S+1-j]
				end

				localImpReq2::Float64 = impReqForIncrStepSize*local_best_cost3

				max_improvement = maximum(improvements)
				min_improvement = minimum(improvements)

				if(min_improvement > -localImpReq2)
					if(max_improvement < localImpReq2)
						stepLength *= stepLengthIncrFactor
						numStepSizeIncr += 1
					end
				end
			end

			# update the best found solutions used for termination and step length modification criteria
			if(iter_S > termCheckPeriod)
				if(cost_soln[iter_S-termCheckPeriod] < local_best_cost)
					local_best_cost = cost_soln[iter_S-termCheckPeriod]
				end
			end
			if(iter_S > stepLengthDecrCheckIter)
				if(cost_soln[iter_S-stepLengthDecrCheckIter] < local_best_cost2)
					local_best_cost2 = cost_soln[iter_S-stepLengthDecrCheckIter]
				end
			end
			if(iter_S > stepLengthIncrCheckIter)
				if(cost_soln[iter_S-stepLengthIncrCheckIter] < local_best_cost3)
					local_best_cost3 = cost_soln[iter_S-stepLengthIncrCheckIter]
				end
			end
			#*========================================================

		end # S
		#*========= END OPTIMIZATION LOOP ==========

		println("")
		println("Number of step size incr/decr: ",numStepSizeIncr,"/",numStepSizeDecr)
		println("Step length change factor: ",round.(stepLength./initialStepLength,displayPrecision),"\n")

		# extract the best found solution for the regular and wind generators
		const p0_reg_soln = zeros(Float64,numRegGen)
		const rp_reg_soln = zeros(Float64,numRegGen)
		const rm_reg_soln = zeros(Float64,numRegGen)
		const alpha_reg_soln = zeros(Float64,numRegGen)

		const p0_wind_soln = zeros(Float64, numWindGen)
		const rp_wind_soln = zeros(Float64, numWindGen)
		const rm_wind_soln = zeros(Float64, numWindGen)
		const alpha_wind_soln = zeros(Float64, numWindGen)

		for i = 1:numRegGen
			p0_reg_soln[i] = p0_best_solution[i]
			rp_reg_soln[i] = rp_best_solution[i]
			rm_reg_soln[i] = rm_best_solution[i]
			alpha_reg_soln[i] = alpha_best_solution[i]
		end

		for i = 1:numWindGen
			p0_wind_soln[i] = p0_best_solution[i+numRegGen]
			rp_wind_soln[i] = rp_best_solution[i+numRegGen]
			rm_wind_soln[i] = rm_best_solution[i+numRegGen]
			alpha_wind_soln[i] = alpha_best_solution[i+numRegGen]
		end

		tic()

		best_cost, best_cost_power, best_lineViolProb, best_sumLineViol, best_cost_reserves::Float64, best_cost_pos_res_pen::Float64, best_cost_neg_res_pen::Float64, power_gen_soln, wind_spill_soln, ~, ~, ~ = estimateCostOfSoln(p0_best_solution,rp_best_solution,rm_best_solution,alpha_best_solution,tau_sat,gamma_res,tau_pos,gamma_line_iter,delta_line,true,true)

		finalCostTime::Float64 = toq()

		#*========= PRINT BEST FOUND SOLUTION ==========

		println("\n*********************")
		println("Lowest cost solution corresponding to iteration #",best_iter)
		println("Expected cost: ",round.(best_cost,displayPrecision))
		println("Expected cost without line penalty costs: ",round.(best_cost_power,displayPrecision))
		println("Expected sum of line flow limit violations: ",round.(best_sumLineViol,displayPrecision))
		println("ANY line limit violation prob: ",round.(best_lineViolProb,displayPrecision))

		if(printFinalSolution)
			println("Regular generator variables come first")
			println("p0: ",round.(p0_best_solution,displayPrecision))
			println("rp: ",round.(rp_best_solution,displayPrecision))
			println("rm: ",round.(rm_best_solution,displayPrecision))
			println("alpha: ",round.(alpha_best_solution,displayPrecision),"\n")

			println("\n")
			println("Regular generator variables")
			println("p0_reg: ",p0_reg_soln)
			println("rp_reg: ",rp_reg_soln)
			println("rm_reg: ",rm_reg_soln)
			println("alpha_reg: ",alpha_reg_soln)

			println("\n")
			println("Wind generator variables")
			println("p0_wind: ",p0_wind_soln)
			println("rp_wind: ",rp_wind_soln)
			println("rm_wind: ",rm_wind_soln)
			println("alpha_wind: ",alpha_wind_soln)
		end

		println("")
		println("Final cost time: ",round.(finalCostTime,2))
		println("")

		iterTime = toq()
		println("RAN ",actualNumIter," ITERATIONS WITH TOTAL TIME FOR LINE PENALTY ITERATION: ",round.(iterTime,2))

		#*========= STORE RESULTS ==========
		if(storeResults)

			open(best_soln_file, "a") do f
				write(f,"NOTE: regular generator variables come first \n\n")
				write(f,"Expected cost: $best_cost \n")
				write(f,"Expected cost without line penalty costs: $best_cost_power \n")
				write(f,"Expected sum of line flow limit violations: $best_sumLineViol \n")
				write(f,"Probability that ANY line limit is violated: $best_lineViolProb \n\n")
				write(f,"p0: $p0_best_solution \n")
				write(f,"rp: $rp_best_solution \n")
				write(f,"rm: $rm_best_solution \n")
				write(f,"alpha: $alpha_best_solution \n\n")
				write(f,"Regular generator variables: \n")
				write(f,"p0_reg: $p0_reg_soln \n")
				write(f,"rp_reg: $rp_reg_soln \n")
				write(f,"rm_reg: $rm_reg_soln \n")
				write(f,"alpha_reg: $alpha_reg_soln \n \n")
				write(f,"Wind generator variables: \n")
				write(f,"p0_wind: $p0_wind_soln \n")
				write(f,"rp_wind: $rp_wind_soln \n")
				write(f,"rm_wind: $rm_wind_soln \n")
				write(f,"alpha_wind: $alpha_wind_soln \n")
				write(f,"\n")
			end

			open(p0_best_soln_file, "a") do f
				for i = 1:numGen
					write(f,"$(p0_best_solution[i]) ")
				end
				write(f,"\n")
			end

			open(rp_best_soln_file, "a") do f
				for i = 1:numGen
					write(f,"$(rp_best_solution[i]) ")
				end
				write(f,"\n")
			end

			open(rm_best_soln_file, "a") do f
				for i = 1:numGen
					write(f,"$(rm_best_solution[i]) ")
				end
				write(f,"\n")
			end

			open(alpha_best_soln_file, "a") do f
				for i = 1:numGen
					write(f,"$(alpha_best_solution[i]) ")
				end
				write(f,"\n")
			end

			open(best_cost_file, "a") do f
				write(f,"$best_cost \n")
			end

			open(best_cost_power_file, "a") do f
				write(f,"$best_cost_power \n")
			end

			open(best_cost_reserves_file, "a") do f
				write(f,"$best_cost_reserves \n")
			end

			open(best_cost_pos_res_file, "a") do f
				write(f,"$best_cost_pos_res_pen \n")
			end

			open(best_cost_neg_res_file, "a") do f
				write(f,"$best_cost_neg_res_pen \n")
			end

			open(best_sum_viol_file, "a") do f
				write(f,"$best_sumLineViol \n")
			end

			open(best_viol_prob_file, "a") do f
				write(f,"$best_lineViolProb \n")
			end

			open(power_gen_soln_file, "a") do f
				for i = 1:numGen
					write(f,"$(power_gen_soln[i]) ")
				end
				write(f,"\n")
			end

			open(wind_spill_soln_file, "a") do f
				for i = 1:numWindGen
					write(f,"$(wind_spill_soln[i]) ")
				end
				write(f,"\n")
			end

			open(time_file, "a") do f
				write(f,"\n\n")
				write(f,"TOTAL TIME: $iterTime")
			end
		end
		#*==================================

		# set initial guess for next iteration as best solution of this iteration
		for i = 1:numGen
			p0_init[i] = p0_best_solution[i]
			rp_init[i] = rp_best_solution[i]
			rm_init[i] = rm_best_solution[i]
			alpha_init[i] = alpha_best_solution[i]
		end

	end # line penalties

	trialTime::Float64 = toq()
	println("\nTime for trial #",trial,": ",round.(trialTime,2))

end # trials

toc() # overall time

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation