text:
@testset "BorgMOEA" begin
    @testset "Schaffer1" begin
        res = bboptimize(BlackBoxOptim.Schaffer1Family; Method=:borg_moea,
                         FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true),
                         SearchRange=(-10.0, 10.0), NumDimensions=2, ϵ=0.01,
                         MaxSteps=5000, TraceMode=:silent)
        @test BlackBoxOptim.IGD(BlackBoxOptim.Schaffer1Family.opt_value, pareto_frontier(res),
                                fitness_scheme(res), Val{length(best_candidate(res))}) < 0.05
    end
    @testset "CEC09_UP8" begin
        res = bboptimize(BlackBoxOptim.CEC09_Unconstrained_Set[8]; Method=:borg_moea,
                         NumDimensions=4, ϵ=0.025,
                         MaxSteps=50000, TraceMode=:silent)
        @test BlackBoxOptim.IGD(BlackBoxOptim.CEC09_Unconstrained_Set[8].opt_value, pareto_frontier(res),
                                fitness_scheme(res), Val{length(best_candidate(res))}) < 0.13
    end
end


text:
## --- tc1/tc1.jl

   @test find_tc1_crust(33,-100) == 40.0
   @test find_tc1_crust([33,-30,70],[-100,20,-130]) == [40.0, 38.0, 46.0]

   @test find_tc1_lith(33,-100) == 132.0
   @test find_tc1_lith([33,-30,70],[-100,20,-130]) == [132.0, 123.0, 167.0]

   @test find_tc1_age(33,-100) == (1400.0, 1100.0, 1700.0)
   @test find_tc1_age([33,-30,70],[-100,20,-130]) ==
      ([1400.0, 1400.0, 2100.0], [1100.0, 1100.0, 1700.0], [1700.0, 1700.0, 2500.0])

## --- PartitionCoefficients/PartitionCoefficients.jl

   @test Float16.(claiborne_zircon_kd.(["Hf","Th","U","Y","Nb","Nd","Sm","Tb","Eu","Dy","Gd","Ho","Er","Yb","Tm","Lu"], 800)) ==
      Float16[1.303e3, 6.46, 36.97, 33.47, 0.2556, 0.02368, 0.4363, 10.695, 0.8784, 19.56, 3.293, 40.34, 62.88, 101.75, 94.94, 126.4]

## --- Elevation.jl

   lat = [43.7022,-26.2041,-19.5723,-34.9285,46.4908]
   lon = [-72.2896,28.0473,65.7550,138.6007,9.8355]
   @test find_geolcont(lat, lon) == [3, 1, 7, 5, 2]

## ---

ext: jl


text:
"""
This device only has one channel so channel specific functions
and optional arguments are not available.

# Available functions
- [`enable_output(obj::Instr{VersatilePower})`](@ref)
- [`disable_output(obj::Instr{VersatilePower})`](@ref)
- [`set_voltage(obj::Instr{VersatilePower})`](@ref)
- [`get_voltage(obj::Instr{VersatilePower})`](@ref)
- [`set_current_limit(obj::Instr{VersatilePower}, current)`](@ref)
- [`get_current_limit(obj::Instr{VersatilePower})`](@ref)

# Helpers
- `remote_mode()`: sets the device to remote mode. Automatically called on initialize
- `local_mode()`: sets the device to local mode. Automatically called on terminate

This instrument has a remote and local mode. Some commands do not
work while the device is in local mode, thus when initializing this
device `remote_mode` is called automatically and the device is always in
remote mode by default.

These functions should not be directly needed but if
for some reason you need to switch modes while using the device
you can use `local_mode` to turn the device back to
local mode.

"""
struct VersatilePower <: PowerSupply end

"""
    enable_output(obj::Instr{VersatilePower})

This will enable an output on a device.

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply
"""
enable_output(obj::Instr{VersatilePower}) = write(obj, "OUTPUT ON")

"""
    disable_output(obj::Instr{VersatilePower})

This will disable an output on a device.

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply
"""
disable_output(obj::Instr{VersatilePower}) = write(obj, "OUTPUT OFF")

"""
    get_output(obj::Instr{VersatilePower})

This will return the state of an output on a device.

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply

Returns:
  String: {"OFF"|"ON"}
"""
get_output(obj::Instr{VersatilePower}) = query(obj, "OUTPUT?")

"""
    set_voltage(obj::Instr{VersatilePower}, num::Voltage)

This will change the voltage output of a device.

Supported Instruments:
   - Power supply

Returns:
  Nothing
"""
set_voltage(obj::Instr{VersatilePower}, num::Voltage) = write(obj, "VOLTAGE $(raw(num))")

"""
    get_voltage(obj::Instr{VersatilePower})

This will return the voltage of a device

Supported Instruments:
   - Power supply

Returns:
  Voltage
"""
get_voltage(obj::Instr{VersatilePower}) = f_query(obj, "VOLTAGE?") * V

"""
    set_current_limit(obj::Instr{VersatilePower}, num::Current)

This will change the current limit of a device

Supported Instruments:
   - Power supply

Returns:
  Nothing
"""
set_current_limit(obj::Instr{VersatilePower}, num::Current) = write(obj, "CURRENT $(raw(num))")

"""
    get_current_limit(obj::Instr{VersatilePower})

This will return the current limit of a device.

Supported Instruments:
   - Power supply

Returns:
  Current Limit
"""
get_current_limit(obj::Instr{VersatilePower}) = query(obj, "CURRENT?")

"""
    remote_mode(obj::Instr{VersatilePower})

Set device to remote mode
"""
remote_mode(obj::Instr{VersatilePower}) = write(obj, "SYSTEM:MODE REMOTE")

"""
    local_mode(obj::Instr{VersatilePower})

Set device to remote mode
"""
local_mode(obj::Instr{VersatilePower}) =   write(obj, "SYSTEM:MODE LOCAL")


text:
#####
##### Dynamic launch configuration
#####

using CUDAdrv

function launch_config(grid, dims)
    return function (kernel)
        fun = kernel.fun
        config = launch_configuration(fun)

        # Adapt the suggested config from 1D to the requested grid dimensions
        if dims == :xyz
            t = floor(Int, cbrt(config.threads))
            threads = [t, t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Ny, grid.Nz] ./ t)
        elseif dims == :xy
            t = floor(Int, sqrt(config.threads))
            threads = [t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Ny] ./ t)
        elseif dims == :xz
            t = floor(Int, sqrt(config.threads))
            threads = [t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Nz] ./ t)
        else
            error("Unsupported launch configuration: $dims")
        end

        return (threads=Tuple(threads), blocks=Tuple(blocks))
    end
end


text:
function compute(n::Float64)::Int
    function is_prime_fast(x::Int)::Bool
        if x == 3
            return true
        elseif !(x % 6 == 1 || x % 6 == 5)
            return false
        end
        for i in 6:12:trunc(Int, √x)
            if x % (i - 1) == 0 || x % (i + 1) == 0 || x % (i + 7) == 0
                return false
            end
        end
        return true
    end
    j = 3
    primes = 0
    while true
        k = (j - 3)j + 3
        primes += is_prime_fast(k) + is_prime_fast(k + j - 1) + is_prime_fast(k + 2j - 2)
        if primes / (2j - 1) < n
            return j
        end
        j += 2
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export logforward!

"""
`logforward!(logalpha,Nt,a,A,y)`

In-place version of [`logforward`](@ref).
"""
function logforward!(alpha::AbstractArray{T},Nt::Integer,
                     a::AbstractVector{T},
                     A::AbstractMatrix{T},
                     y::AbstractArray{T}) where {T <: AbstractFloat}
  Ns = size(A,1)
  fill!(alpha,-T(Inf))
  for j = 1:Ns
    alpha[1,j] = a[j] + y[1,j]
  end

  for t = 2:Nt
    for j = 1:Ns
      for i = 1:Ns
        aij = A[i,j]
        if aij != -T(Inf)
          alpha[t,j] = logadd(alpha[t,j],alpha[t-1,i]+aij+y[t,j])
        end
      end
    end
  end

  logML = -T(Inf)
  for i = 1:Ns
    logML = logadd(logML,alpha[Nt,i])
  end
  return alpha, logML/Nt
end

# forward with constrained path
function logforward!(alpha::AbstractArray{T},
                     Nt::Integer,
                     t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                     A::AbstractMatrix{T},
                     y::AbstractArray{T}) where {D<:Integer, T <: AbstractFloat}
  Ns = size(A,1)

  fill!(alpha,-T(Inf))
  idx0 = unique(t2IJ[1][1])
  for j in idx0
    alpha[1,j] = y[1,j] - log(length(idx0))
  end

  for t = 2:Nt
    I,J = t2IJ[t-1]
    for k in eachindex(I)
      i,j = I[k],J[k]
      alpha[t,j] = logadd(alpha[t,j],alpha[t-1,i]+A[i,j]+y[t,j])
    end
  end

  logML = -T(Inf)
  for i = 1:Ns
    logML = logadd(logML,alpha[Nt,i])
  end
  return alpha, logML/Nt
end

export logforward

"""
`logalpha, logML = logforward(Nt,a,A,y)`

Computes the scaled forward log-probabilities (``\\hat{\\alpha}``) of an Hidden Markov Model (HMM) of `Ns` states with log-transition matrix `A` (must be a `Ns` × `Ns` matrix), initial log-probabilities `a` (must be a `Ns`-long vector) and observation log-likelihoods `y` (must be of size `Nt2` × `Ns` with `Nt2` ≥ `Nt`).

Returns:
* `logalpha` a matrix of size `(Nt,Ns)` containing the forward log-probabilities
* `logML` the log likelihood of the observation normalized by the sequence length, i.e. ``\\log(P(\\mathbf{X}))/N_t``.

`logalpha, logML = logforward(Nt,t2tr,A,y)`

Returns the forward log-probabilities with a constrained path (see [`forward`](@ref)).
"""
function logforward(Nt::Integer,a,A::AbstractMatrix{T},y) where {T <: AbstractFloat}
  Ns = size(A,1)
  alpha = -T(Inf)*ones(T,Nt,Ns)
  logforward!(alpha,Nt,a,A,y)
end

# y is a tensor Nt2 × Ns × Nb
function logforward(Nts::Array{D},
                    as::AbstractMatrix{T},
                    As::AbstractArray{T,3},
                    y::AbstractArray{T,3}) where {D <: Integer,
                                          T <: AbstractFloat}
  Ns = size(As,1)
  Nb = size(y,3) # batch number
  @assert Nb == length(Nts) == size(as,2) == size(As,3)
  alphas = [-T(Inf)*ones(T,Nt,Ns) for Nt in Nts]
  logMLs = zeros(T,Nb)
  for i = 1:Nb
    _, logML = logforward!(alphas[i],Nts[i],view(as,:,i),view(As,:,:,i),view(y,:,:,i))
    logMLs[i] = logML
  end
  return alphas, logMLs
end

# y is a tensor Nt2 × Ns × Nb
function logforward(Nts::Array{D},
                    as::AbstractVector,
                    A::AbstractMatrix{T},
                    y::Array{T,3}) where {D <: Integer,
                                          T <: AbstractFloat}
  Ns = size(A,1)
  Nb = size(y,3) # batch number
  @assert Nb == length(Nts) == length(as)
  alphas = [-T(Inf)*ones(T,Nt,Ns) for Nt in Nts]
  logMLs = zeros(T,Nb)
  for i = 1:Nb
    _, logML = logforward!(alphas[i],Nts[i],as[i],A,view(y,:,:,i))
    logMLs[i] = logML
  end
  return alphas, logMLs
end


text:
# parameters that govern counterfactual

## All parameters remain unchanged.


text:
#
#   Thomson Problem on Oblique(n,m)
#
using Manopt, Manifolds, ManifoldsBase, LinearAlgebra, Random
Random.seed!(42)
n = 50
m = 20
M = Oblique(n, m)

function F(X::Array{Float64,2})
    f = 0
    for i in 1:m
        for j in 1:m
            if i != j
                f = f + 1 / (norm(X[:, i] - X[:, j])^2)
            end
        end
    end
    return f
end

function ∇F(X::Array{Float64,2})
    g = zeros(n, m)
    Id = Matrix(I, n, n)
    for i in 1:m
        f = zeros(n, 1)
        for j in 1:m
            if i != j
                f = f + 1 / (1.0 - X[:, i]' * X[:, j]) * X[:, j]
            end
        end
        g[:, i] = (Id - X[:, i] * X[:, i]') * f
    end
    return g
end

x = random_point(M)

@time quasi_Newton(
    M,
    F,
    ∇F,
    x;
    memory_size=100,
    vector_transport_method=PowerVectorTransport(ParallelTransport()),
    stopping_criterion=StopWhenGradientNormLess(norm(M, x, ∇F(x)) * 10^(-6)),
    debug=[:Iteration, " ", :Cost, "\n", 1, :Stop],
)


text:
using DashBootstrapComponents, DashHtmlComponents

input_group = html_div([
    dbc_inputgroup(
        [dbc_inputgroupaddon("Large", addon_type="prepend"), dbc_input()],
        size="lg",
    ),
    html_br(),
    dbc_inputgroup([dbc_inputgroupaddon("Default", addon_type="prepend"), dbc_input()]),
    html_br(),
    dbc_inputgroup(
        [dbc_inputgroupaddon("Small", addon_type="prepend"), dbc_input()],
        size="sm",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""
Module Email

Provides functions for operations in an Email client.

# Functions

###### Open inbox
- [`inbox`](@ref)

To see a description of a function type `?<functionname>`.

See also: [`Internet`](@ref)
"""
module Email

import DefaultApplication
import ..JustSayIt: @voiceargs

## CONSTANTS

const EMAIL_ENGINE = "https://mail.google.com"

## FUNCTIONS

@doc """
    email `inbox` | `outbox`

Manage e-mails, performing one of the following actions:
- open `inbox`
- open `outbox`
"""
email
@enum Action inbox outbox
@voiceargs action=>(valid_input_auto=true) function email(action::Action)
    if     (action == inbox)  open_inbox()
    elseif (action == outbox) open_outbox()
    else                      @info "unknown action"  #NOTE: this should never happen.
    end
end

@doc "Open E-mail inbox."
open_inbox() = DefaultApplication.open(EMAIL_ENGINE)

@doc "Open E-mail outbox."
open_outbox() = DefaultApplication.open(EMAIL_ENGINE * "/" * "#sent")

end # module Email


text:
require("ZMQ")
require("JSON")

module Mongrel2
using ZMQ
import JSON

export M2Request,M2Connection,
		parse_request, is_disconnected, connect, reply, reply_http, run_server, disconnect_client

type M2Request
	sender_id
	connection_id
	path
	headers
	body
	json_data

	function M2Request(s,c,p,h,b)
		if h["METHOD"] == "JSON"
			new(s,c,p,h,b,JSON.parse(b))
		else
			new(s,c,p,h,b,Dict())
		end
	end
end

type M2Connection
	ctx::ZMQContext
	reqs::ZMQSocket
	resp::ZMQSocket
	sub_addr::String
	pub_addr::String

end

function parse_netstring (str::String)
	 i = search(str, ':')
	 s = str[i+1:end]
	 len = int(str[1:i-1])
	 if s[len+1] != ','; error ("Netstring does not end with comma : $str"); end
	 if len==0
	 	return "", ""
	 else
		return s[1:len], s[len+2:end]
	end
end

function parse_request(msg::ZMQMessage)
	s=ASCIIString[msg]
	#Uncomment for debug:
	#print ("Recd message: $s \n\n");flush(stdout_stream)

	r = split(s, ' ')
	sender = r[1]
	connection_id = r[2]
	path = r[3]
	rest = r[4]
	if length(r) > 4
		for i = r[5:end]
			rest = string(rest, " ", i)
		end
	end

	headers, head_rest = parse_netstring(rest)

	(body, ) = parse_netstring(head_rest)

	return M2Request(sender, connection_id, path, JSON.parse(headers), body)
end

function is_disconnected(req::M2Request)
	if req.headers["METHOD"] == "JSON"
		return req.json_data["type"] == "disconnect"
	else
		return false
	end

end

function connect(sender_id::String, sub_addr::String, pub_addr::String)
	ctx = ZMQContext(1)
	reqs = ZMQSocket(ctx, ZMQ_UPSTREAM)
	ZMQ.connect(reqs, sub_addr)
	resp = ZMQSocket(ctx, ZMQ_PUB)
	ZMQ.connect(resp, pub_addr)
	ZMQ.set_identity(resp, sender_id)
	return M2Connection(ctx, reqs, resp, sub_addr, pub_addr)
end

#internal function
function send_resp(conn::M2Connection, uuid, conn_id, msg)
	str = "$uuid $(length(conn_id)):$(conn_id), $msg"
	ZMQ.send(conn.resp, ZMQMessage(str))
	#Uncomment for debug:
	#print("Sent Reply to $conn_id: $str \n\n");flush(stdout_stream)
end

recv(conn::M2Connection) = parse_request(ZMQ.recv(conn.reqs))::M2Request

#This shouldnt be required  except for HTTP 1.0 clients, or websockets.
#Normally, Let the browser manage keepalive.
disconnect_client(conn::M2Connection, req::M2Request) = send_resp(conn, req.sender_id, req.connection_id, "")

reply(conn::M2Connection, req::M2Request, msg::String) = send_resp(conn, req.sender_id, req.connection_id, msg)
reply_http (conn::M2Connection, req::M2Request, body, code, headers::Dict{String, String}) =
								reply(conn, req, http_response(body, code, headers))
reply_http (conn::M2Connection, req::M2Request, body, headers) = reply_http(conn, req, body, 200, headers)
reply_http (conn::M2Connection, req::M2Request, body) = reply_http(conn, req, body, 200, Dict{String, String}())

function http_response(body, code, headers::Dict{String, String})
	headers["Content-Length"] = string(length(body))
	headers_s = ""
	for (k, v) = headers
		headers_s = string(headers_s, "$(k): $(v)\r\n")
	end
	return "HTTP/1.1 $code $(StatusMessage[int(code)])\r\n$(headers_s)\r\n\r\n$(body)"
end

function run_server(sender_id, sub_addr, pub_addr)
	conn::M2Connection = connect(sender_id, sub_addr, pub_addr)

	function runner()
		while true
			request::M2Request = recv(conn)
			produce((conn, request))
		end
	end
	print("Julia Mongrel2 hander started, connecting back on [$sub_addr] and [$pub_addr] \n")
	#flush(stdout_stream)
	return Task(runner)
end

const StatusMessage = Dict{Int, String}()
StatusMessage[100] = "Continue"
StatusMessage[101] = "Switching Protocols"
StatusMessage[200] = "OK"
StatusMessage[201] = "Created"
StatusMessage[202] = "Accepted"
StatusMessage[203] = "Non-Authoritative Information"
StatusMessage[204] = "No Content"
StatusMessage[205] = "Reset Content"
StatusMessage[206] = "Partial Content"
StatusMessage[300] = "Multiple Choices"
StatusMessage[301] = "Moved Permanently"
StatusMessage[302] = "Found"
StatusMessage[303] = "See Other"
StatusMessage[304] = "Not Modified"
StatusMessage[305] = "Use Proxy"
StatusMessage[307] = "Temporary Redirect"
StatusMessage[400] = "Bad Request"
StatusMessage[401] = "Unauthorized"
StatusMessage[402] = "Payment Required"
StatusMessage[403] = "Forbidden"
StatusMessage[404] = "Not Found"
StatusMessage[405] = "Method Not Allowed"
StatusMessage[406] = "Not Acceptable"
StatusMessage[407] = "Proxy Authentication Required"
StatusMessage[408] = "Request Timeout"
StatusMessage[409] = "Conflict"
StatusMessage[410] = "Gone"
StatusMessage[411] = "Length Required"
StatusMessage[412] = "Precondition Failed"
StatusMessage[413] = "Request Entity Too Large"
StatusMessage[414] = "Request-URI Too Large"
StatusMessage[415] = "Unsupported Media Type"
StatusMessage[416] = "Request Range Not Satisfiable"
StatusMessage[417] = "Expectation Failed"
StatusMessage[500] = "Internal Server Error"
StatusMessage[501] = "Not Implemented"
StatusMessage[502] = "Bad Gateway"
StatusMessage[503] = "Service Unavailable"
StatusMessage[504] = "Gateway Timeout"
StatusMessage[505] = "HTTP Version Not Supported"

end


text:
# adapted from plots because its not in Recipes Base
function partialcircle(start_θ, end_θ, n = 20, r=1)
    [r*cos(u) for u in linspace(start_θ, end_θ, n)],
    [r*sin(u) for u in linspace(start_θ, end_θ, n)]
end

function make_circle(cx, cy, r)
    x, y = partialcircle(0, 2pi, 20, r)
    x .+= cx
    y .+= cy
    collect(zip(x, y))
end

@recipe function plot(env::CartpoleEnvironment; force=0.)
    sim = env.sim
    xlims --> (sim.xlimits[1], sim.xlimits[2])
    ylims --> sim.pole_length .* (-1.2,1.2)
    ratio --> true
    xlabel --> "x [meter]"
    ylabel --> "y [meter]"
    title --> @sprintf("Cartpole at t=%0.2f s", sim.time)
    min_force, max_force = extrema(env.actions)
    cart_pos = configuration(sim.state)[1]
    cart_offset = -0.04
    cart_height = 0.1
    cart_width = 0.6
    arrow_height = cart_height / 2
    arrow_width_plus = 0.15 / max_force * force
    arrow_width_minus = 0.15 / min_force * force
    pole_width = sim.pole_com == sim.pole_length ? 0.005 : 0.02
    wheel_radius = 0.07
    wheel_offset = 0.2
    wheel_turns = (cart_pos+10) / (2pi * wheel_radius)
    wheel_angle = -2pi*(wheel_turns - floor(wheel_turns))
    actual_xlims = plotattributes[:xlims]
    @series begin # floor
        seriescolor := :grey
        linewidth := 1
        label := ""
        [actual_xlims...], cart_offset .- cart_height .- wheel_radius .+ zeros(2)
    end
    if force > 0
        @series begin
            seriescolor := :red
            seriestype := :shape
            alpha := 0.2
            label := ""
            (cart_pos .+ cart_width./2 .+ [0, 0, arrow_width_plus, arrow_width_plus*1.2, arrow_width_plus],
             cart_offset .- cart_height./2 .+ arrow_height./2 .* [-1, 1, 1, 0, -1])
        end
    elseif force < 0
        @series begin
            seriescolor := :red
            seriestype := :shape
            alpha := 0.2
            label := ""
            (cart_pos .- cart_width./2 .- [0, 0, arrow_width_minus, arrow_width_minus*1.2, arrow_width_minus],
             cart_offset .- cart_height./2 .+ arrow_height./2 .* [-1, 1, 1, 0, -1])
        end
    end
    @series begin # triangle
        seriescolor := :lightgrey
        seriestype := :shape
        label := ""
        cart_pos .+ [-.1, 0., .1], cart_offset .+ [0., 0.1, 0.]
    end
    @series begin # big cart
        seriescolor := :lightgrey
        seriestype := :shape
        label := ""
        cart_pos .+ cart_width .* [-1, 1, 1, -1] ./ 2, cart_offset .+ [-cart_height, -cart_height, 0., 0.]
    end
    @series begin # wheel left
        seriestype := :shape
        seriescolor := :grey
        label := ""
        make_circle(cart_pos - wheel_offset, cart_offset - cart_height, wheel_radius)
    end
    @series begin # wheel right
        seriestype := :shape
        seriescolor := :grey
        label := ""
        make_circle(cart_pos + wheel_offset, cart_offset - cart_height, wheel_radius)
    end
    wheel_x = [0.7wheel_radius*cos(wheel_angle)]
    wheel_y = [0.7wheel_radius*sin(wheel_angle)]
    @series begin # wheel left scatter
        seriescolor := :black
        seriestype := :scatter
        markersize := 1
        label := ""
        wheel_x .+ cart_pos .- wheel_offset, wheel_y .+ cart_offset .- cart_height
    end
    @series begin # wheel right scatter
        seriescolor := :black
        seriestype := :scatter
        markersize := 1
        label := ""
        wheel_x .+ cart_pos .+ wheel_offset, wheel_y .+ cart_offset .- cart_height
    end
    p_bot_l = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector(-pole_width,0,0.)),
            default_frame(bodies(sim.mechansim)[1])).v
    p_bot_r = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector(pole_width,0,0.)),
            default_frame(bodies(sim.mechansim)[1])).v
    p_top_l = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector{3}(-pole_width,0,sim.pole_length)),
            default_frame(bodies(sim.mechansim)[1])).v
    p_top_r = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector{3}(pole_width,0,sim.pole_length)),
            default_frame(bodies(sim.mechansim)[1])).v
    @series begin # pole itself
        seriescolor := :grey
        seriestype := :shape
        label := ""
        [p_bot_l[1],p_bot_r[1],p_top_r[1],p_top_l[1]], [p_bot_l[3],p_bot_r[3],p_top_r[3],p_top_l[3]]
    end
    if sim.pole_com == sim.pole_length
        p_bot_c = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector(0,0,0.)),
            default_frame(bodies(sim.mechansim)[1])).v
        p_top_c = transform(sim.state,
            Point3D(frame_after(joints(sim.mechansim)[3]), SVector(0.,0,sim.pole_length)),
            default_frame(bodies(sim.mechansim)[1])).v
        @series begin # sphere on top
            seriestype := :shape
            seriescolor := :grey
            label := ""
            make_circle(p_top_c[1], p_top_c[3], 0.07)
        end
    end
    @series begin # sphere on bottom
        seriestype := :shape
        seriescolor := :lightgrey
        label := ""
        make_circle(cart_pos, 0., 0.03)
    end
    if get(plotattributes, :legend, :best) != :none
        @series begin # x line
            seriescolor := :grey
            linewidth := 1
            linestyle := :dash
            label := ""
            [cart_pos, cart_pos], [cart_offset .- cart_height .- wheel_radius .- 0.05, sim.pole_length]
        end
        @series begin # x text
            seriescolor := :grey
            label := ""
            series_annotations := [@sprintf("x = %0.2f", cart_pos), ""]
            [cart_pos], [cart_offset .- cart_height .- wheel_radius .- 0.15]
        end
    end
    nothing
end


text:
using Random
using ForwardDiff
using StatsBase
using GalacticOptim
using Flux
using Optim
using Plots

Random.seed!(0)

# taken from DiffEqFlux.jl
abstract type FastLayer <: Function end

struct FastDense{F,P} <: FastLayer
    in::Int
    out::Int
    σ::F
    initial_params::P
    function FastDense(in::Integer, out::Integer, σ=identity; initW=Flux.glorot_uniform, initb=Flux.zeros32)
        initial_params() = vcat(vec(initW(out,in)), initb(out))
        new{typeof(σ),typeof(initial_params)}(in,out,σ,initial_params)
    end
end

initial_params(f::FastDense) = f.initial_params()
paramlength(f::FastDense) = f.out*(f.in+1)

function (f::FastDense)(x,θ)
    W = reshape(θ[1:(f.out*f.in)], f.out, f.in)
    b = θ[(f.out*f.in+1):end]
    f.σ.(W*x .+ b)
end

struct FastChain{T<:Tuple} <: FastLayer
    layers::T
end
FastChain(ls...) = FastChain(ls)

applychain(::Tuple{}, x, p) = x
function applychain(fs::Tuple, x, p)
    applychain(Base.tail(fs), first(fs)(x,p[1:paramlength(first(fs))]), p[(paramlength(first(fs))+1):end])
end
(c::FastChain)(x,p) = applychain(c.layers, x, p)
paramlength(c::FastChain) = sum(paramlength(x) for x in c.layers)
initial_params(c::FastChain) = vcat(initial_params.(c.layers)...)

##########  PINN  ##############################################################

function lotkavolterra(u::AbstractVector,θ)
    α, β, γ, δ = θ
    u₁, u₂ = u

    du₁ = α*u₁ - β*u₁*u₂
    du₂ = δ*u₁*u₂ - γ*u₂

    [du₁, du₂]
end
#lotkavolterra(u::AbstractMatrix,θ) = mapreduce(ui->lotkavolterra(ui,θ), hcat, eachcol(u))
function lotkavolterra(u::AbstractMatrix,θ)
    α, β, γ, δ = θ
    u₁, u₂ = u[1,:], u[2,:]

    du₁ = α.*u₁ - β.*u₁.*u₂
    du₂ = δ.*u₁.*u₂ - γ.*u₂

    [du₁ du₂]'
end

function forward_derivative(u::FastLayer,x::AbstractVector,ε::AbstractVector,θ)
    J = ForwardDiff.jacobian(i->u(i,θ),x)
    J * ε
end
function forward_derivative(u::FastLayer,xs::AbstractMatrix,ε::AbstractVector,θ)
    mapreduce(x->forward_derivative(u,x,ε,θ), hcat, eachcol(xs))
end
function forward_derivative(u::FastLayer,x::AbstractArray,v::Int,θ)
    ε = Flux.onehot(v,1:size(x,1))
    forward_derivative(u,x,ε,θ)
end
function get_ε(dim, der_num)
    epsilon = cbrt(eps(Float32))
    e = map(1:dim) do i
        # NOTE: needed to prevent Zygote mutation error
        i==der_num ? epsilon : 0f0
    end
end
function numeric_derivative(u::FastLayer,x,ε::Vector,order::Int,θ::Vector)
    _epsilon = 1 / (2*cbrt(eps(eltype(θ))))
    if order > 1
        a = numeric_derivative(u,x .+ ε, ε, order-1, θ)
        b = numeric_derivative(u,x .- ε, ε, order-1, θ)
        (a .- b) .* _epsilon
    else
        (u(x .+ ε,θ) .- u(x .- ε,θ)) .* _epsilon
    end
end

struct PINN{U<:FastLayer,L,R}
    u::U
    lhs::L
    rhs::R
end

function internal_loss(m::PINN, x::AbstractArray, θ::AbstractVector)
    u = m.u
    mean(abs2, m.lhs(u,x,θ) - m.rhs(u,x,θ))
end
testloss(θ) = internal(pinn, rand(1,10), θ)

function data_loss(m::PINN, x::AbstractArray, y::AbstractArray, θ::AbstractVector)
    mean(abs2, m.u(x,θ) - y)
end

hdim = 10
m = FastChain(
    FastDense(1,hdim,tanh),
    FastDense(hdim,hdim,tanh),
    FastDense(hdim,2))
θ = initial_params(m)

lhs(u,x,θ) = forward_derivative(u,x,1,θ)
lotka_params = [0.1f0, 0.2f0, 0.3f0, 0.2f0]
rhs(u,x,θ) = lotkavolterra(u(x,θ), lotka_params)
pinn = PINN(m,lhs,rhs)

using JLD2
data = load("../../lecture_12/lotkadata.jld2")
dx = 100
de = 0.5f0
xs = Float32.(reshape(data["t"],1,:)[:,1:dx:round(Int,end*de)])
ys = Float32.(data["u"][:,1:dx:round(Int,end*de)])
#xs = Float32.(reshape([data["t"][1]],1,1))
#ys = Float32.(reshape(data["u"][:,1],:,1))

_xs = reshape(collect(0.0f0:2.0f0:100.0f0*de),1,:)
function plotprogress(θ)
    t = vec(xs)
    p1 = scatter(t, ys[1,:], label="Data x", ms=7)
    scatter!(p1, t, ys[2,:], label="Data y", ms=7)
    t = data["t"]
    plot!(p1, t, data["u"][1,:], label=false, lw=3, c=:gray)
    plot!(p1, t, data["u"][2,:], label=false, lw=3, c=:gray)

    us = pinn.u(reshape(t,1,:),θ)
    plot!(p1, t, us[1,:], label="PINN x", lw=3, c=1)
    plot!(p1, t, us[2,:], label="PINN y", lw=3, c=2)
    plot!(p1, xlim=(_xs[1],_xs[end]))
end

function loss(θ,p=nothing)
    l1 = internal_loss(pinn, _xs, θ)
    l2 = data_loss(pinn, xs, ys, θ)
    l1+l2, (l1,l2)
end

function callback(θ,l,ls)
    (internal, dataloss) = ls
    plotprogress(θ) |> display
    @info l internal dataloss
    false
end

display(loss(θ))
ForwardDiff.gradient(first ∘ loss, θ)

func = OptimizationFunction(loss, GalacticOptim.AutoForwardDiff())
prob = OptimizationProblem(func, θ)
opt  = LBFGS()

@time res = solve(prob, opt, maxiters=10_000, cb=Flux.throttle(callback,1), progress=true)
p = plotprogress(res.minimizer)


text:
export IntersectionGraph

AnySet = Union{Set,BitSet}

"""
`IntersectionGraph(setlist)` creates an intersection graph on vertex set
`1:n` from a list of `n` sets (all of type `Set` or all of type
`IntSet`).
"""
function IntersectionGraph(slist::Vector{T}) where {T<:AnySet}
    n = length(slist)
    f = Dict{Int,T}()
    for k = 1:n
        f[k] = slist[k]
    end
    return IntersectionGraph(f)
end

"""
`IntersectionGraph(f)` creates an intersection graph from the
dictionary `f`. The keys are the vertices and the values are the sets
(all of type `Set` or all of type `IntSet`).
"""
function IntersectionGraph(f::Dict{S,T}) where {S,T<:AnySet}
    G = SimpleGraph{S}()
    for v in keys(f)
        add!(G, v)
    end

    vtcs = vlist(G)
    n = length(vtcs)

    for i = 1:n-1
        u = vtcs[i]
        A = f[u]
        for j = i+1:n
            v = vtcs[j]
            B = f[v]

            if length(intersect(A, B)) > 0
                add!(G, u, v)
            end
        end
    end
    cache_save(G, :IntersectionRepresentation, f)
    cache_save(G, :name, "Intersection graph")
    return G
end


text:
#     The forecast_ functions create the forecast dictionary, which contains time-
#     varying parameters.

function forecast_base!(fileNm::String, p::Dict; f::Dict=Dict(),
    startDate::Date=Date(), endDate::Date=Date(), startT::Int=1,
    T::Int=-1, history::Dict=Dict())
    return forecast_base!(CSV.read(fileNm), p; f=f, startDate=startDate,
        endDate=endDate, startT=startT, T=T, history=history)
end

function forecast_base!(data::DataFrame, p::Dict; f::Dict=Dict(),
    startDate::Date=Date(), endDate::Date=Date(), startT::Int=1,
    T::Int=-1, history::Dict=Dict())
    # in the optimization, we are in local time
    data[:Date] += Dates.Hour(p[:utc_offset])

    if startDate != Date()
        startT = findDate(startDate, Date.(data[:Date]))
    end

    if endDate != Date()
        endT = findDate(endDate, Date.(data[:Date]))
        T = endT - startT + 2;
    end

    # Time
    if T == -1
        T = 1 + size(data)[1] - startT+1;
    end
    endT = startT + T - 2;
    f[:startT] = startT; # remember where I started reading from.
    f[:endT] = endT; # remember where I stopped reading from.
    f[:T] = T;
    f[:timestamp] = data[:Date][startT:endT]

    # number of months
    f[:M] = getMonthNum(T-1, f);
    f[:loadKWH] = data[:Campus_KWH][startT:endT]; # kWh

    f[:hourCA] = mod.(Dates.hour.(f[:timestamp])-7, 24)
    return f;
end

function forecast_water!(f::Dict, fileNm::String, p::Dict; history::Dict=Dict())
    return forecast_water!(f, readtable(fileNm), p; history)
end

function forecast_water!(f::Dict, data::DataFrame, p::Dict;
    history::Dict=Dict())

    startT = f[:startT]
    endT = f[:endT]

    if haskey(data, :chillerPE)
        f[:chillerPE] = data[:chillerPE][startT:endT] # kW/ton
    else
        f[:chillerPE] = ones(endT-startT+1) * p[:chillerPE]
    end
    f[:loadHW] = data[:HW_loads][startT:endT] # tons
    f[:loadCW] = data[:CW_loads][startT:endT] # mmbtu
    if haskey(history, :hwTank)
        f[:hwTankInit] = history[:hwTank][end]
    else
        f[:hwTankInit] = p[:hwTankInit]
    end
    if haskey(history, :cwTank)
        f[:cwTankInit] = history[:cwTank][end]
    else
        f[:cwTankInit] = p[:cwTankInit]
    end
    f[:cwTankEnd] = f[:cwTankInit]/2
    f[:hwTankEnd] = f[:hwTankInit]/2
    return f;
end

function forecast_storage!(f::Dict, fileNm::String, p::Dict;
    history::Dict=Dict())
    return forecast_storage!(f, readtable(fileNm), p; history=history)
end

function forecast_storage!(f::Dict, data::DataFrame, p::Dict;
    history::Dict=Dict())

    startT = f[:startT]
    endT = f[:endT]

    if haskey(history, :elecStorageInit)
        f[:elecStorageInit] = history[:elecStorageInit][end]
    else
        f[:elecStorageInit] = p[:elecStorageInit]
    end
    return f;
end

function forecast_prices!(f::Dict, fileNm::String, p::Dict;
    history::Dict=Dict())
    return forecast_prices!(f, readtable(fileNm), p; history=history)
end

function forecast_prices!(f::Dict, data::DataFrame, p::Dict;
    history::Dict=Dict())

    startT = f[:startT]
    endT = f[:endT]
    # add a fixed cost here?
    f[:elecPrice] = data[:Electricity_prices][startT:endT] # $/kW
    f[:gasPrice] = p[:gasPrice]
    #f[:demandCharge] = p[:demandCharge] * ones(f[:M])
    # use demandCharge data to populate f[:demandCharge]
    simMonths = unique(DateTime.(Dates.Year.(f[:timestamp]),
                                 Dates.Month.(f[:timestamp])))
    f[:demandCharge] = [p[:demandCharge][m] for m in simMonths]
    return f
end

function forecast_carbon!(f::Dict, fileNm::String, p::Dict;
    history::Dict=Dict())
    return forecast_carbon!(f, readtable(fileNm), p, history)
end

function forecast_carbon!(f::Dict, data::DataFrame, p::Dict;
    history::Dict=Dict())

    startT = f[:startT]
    endT = f[:endT]

    f[:carbonIntensityGrid] = data[:carbon][startT:endT]*p[:carbon_units]
    f[:carbonIntensityGas] = p[:carbonIntensityGas]
    f[:carbonPrice] = p[:carbonPrice]
    return f
end

function forecast_CBP!(f::Dict, fileNm::String, p::Dict; history::Dict=Dict())
    if !haskey(f, :eGrid)
        error("There should be :eGrid in f")
    end
    f[:event_ts] # contains a list of timestamps at which we are doing demand response

    f[:ts_to_index] #maps ts->integer index
    return f;
end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
module LinearSolversTests

using Test
using Gridap.Algebra
using Gridap.Algebra: NonlinearOperatorMock

using LinearAlgebra
using SparseArrays

diff1(M) = [ [1.0 zeros(1,M-1)]; diagm(1=>ones(M-1)) - Matrix(1.0I,M,M) ]

sdiff1(M) = sparse(diff1(M))

function Laplacian(Nx, Ny, Lx, Ly)
   dx = Lx / (Nx+1)
   dy = Ly / (Ny+1)
   Dx = sdiff1(Nx) / dx
   Dy = sdiff1(Ny) / dy
   Ax = Dx' * Dx
   Ay = Dy' * Dy
   return kron(sparse(1.0I,Ny,Ny), Ax) + kron(Ay, sparse(1.0I,Nx,Nx))
end

n = 10
A = Laplacian(n,n,1,1)
x = rand(n^2)
b = A*x

op = AffineOperator(A,b)
test_nonlinear_operator(op,x,zeros(size(x)),jac=A)

ls = LUSolver()
test_linear_solver(ls,A,b,x)

n = 10
A = Laplacian(n,n,1,1)
x = rand(n^2)
b = A*x

ls = BackslashSolver()
test_linear_solver(ls,A,b,x)

nls = NLSolver(show_trace=false,method=:newton)
x0 = zeros(length(x))
op = AffineOperator(A,b)
@test jacobian(op,x) === op.matrix

cache = solve!(x0,nls,op)
test_nonlinear_solver(nls,op,x0,x)

x1 = copy(x0)
cache = solve!(x1,nls,op,cache)
@test all(x1 .≈ x0)

op = AffineOperator(A,2*b)
solve!(x1,nls,op,cache)
@test all(x1 .≈ 2*x0)

op = AffineOperator(A,b)
x0 = copy(x)
cache = solve!(x0,ls,op,nothing,true)
@test all(x .≈ x0)
cache = solve!(x0,ls,op,nothing)
@test all(x .≈ x0)
cache = solve!(x0,ls,op)
@test all(x .≈ x0)

x0 = copy(x)
cache = solve!(x0,ls,op,cache)
@test all(x .≈ x0)
cache = solve!(x0,ls,op,cache,false)
@test all(x .≈ x0)

op2 = AffineOperator(A,2*b)
x0 = copy(x)
cache = solve!(x0,ls,op2,cache)
@test all(2*x .≈ x0)
cache = solve!(x0,ls,op2,cache,false)
@test all(2*x .≈ x0)
cache = solve!(x0,ls,op2,nothing,false)
@test all(2*x .≈ x0)

nlop = NonlinearOperatorMock()
x = [1.0, 3.0]
b = [0.0,0.0]
A = [-1.0 0.0; 0.0 1.0]
test_nonlinear_operator(nlop,x,b,jac=A)
x0 = copy(x)
nl_cache = solve!(x0,ls,nlop,nothing)
nl_cache = solve!(x0,ls,nlop,nl_cache)

end # module


text: using DataFrames


text:
function initialize!(integrator, cache::ROCK2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ROCK2ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, fp1, fp2, recf = cache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + dt * integrator.eigen_est)/0.811) + 1))
  if mdeg >= 200
    mdeg = 200
  end
  cache.mdeg = max(mdeg, 3) - 2
  cache.mdeg != cache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[cache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  gprev2 = copy(uprev)
  gprev = uprev + temp1 * fsalfirst
  ms[cache.mdeg] < 2 && ( u = gprev )
  # for the second to the ms[cache.mdeg] th stages
  for i in 2:ms[cache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = -1 - κ
    dtμ = dt*μ
    ci1 = dtμ - ν * ci2 - κ * ci3
    u = dtμ * u - ν * gprev - κ * gprev2
    i < ms[cache.mdeg] && (gprev2 = gprev; gprev = u)
    ci3 = ci2
    ci2 = ci1
  end # end if
  # two-stage finishing procedure.
  temp1 = dt * fp1[cache.mdeg]
  temp2 = dt * fp2[cache.mdeg]
  gprev2 = f(u, p, ci1)
  gprev = u + temp1 * gprev2
  ci1 += temp1
  u = f(gprev, p, ci1)
  temp3 = temp2 * (u - gprev2)
  u = gprev + temp1 * u + temp3
  # error estimate
  if integrator.opts.adaptive
    atmp = calculate_residuals(temp3, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.u = u
end

function initialize!(integrator, cache::ROCK2Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator, cache::ROCK2Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack k, k2, tmp, gprev2, gprev, atmp = cache
  @unpack ms, fp1, fp2, recf = cache.constantcache
  ccache = cache.constantcache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + dt * integrator.eigen_est)/0.811) + 1))
  if mdeg >= 200
    mdeg = 200
  end
  ccache.mdeg = max(mdeg, 3) - 2
  ccache.mdeg != ccache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[ccache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  @. gprev2 = uprev
  @. gprev = uprev + temp1 * fsalfirst
  ms[ccache.mdeg] < 2 && ( @. u = gprev )
  # for the second to the ms[ccache.mdeg] th stages
  for i in 2:ms[ccache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = κ - 1
    temp1 = dt * μ
    temp2 = 1 + κ
    temp3 = -κ
    ci1 = temp1 + temp2 * ci2 + temp3 * ci3
    @. u = temp1 * u + temp2 * gprev + temp3 * gprev2
    i < ms[ccache.mdeg] && (gprev2 .= gprev; gprev .= u)
    ci3 = ci2
    ci2 = ci1
  end # end if
  # two-stage finishing procedure.
  temp1 = dt * fp1[ccache.mdeg]
  temp2 = dt * fp2[ccache.mdeg]
  f(k, u, p, ci1)
  @. gprev = u + temp1 * k
  ci1 += temp1
  f(k2, gprev, p, ci1)
  @. tmp = temp2 * (k2 - k)
  @. u = gprev + temp1 * k2 + tmp
  # error estimate
  if integrator.opts.adaptive
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end


text:
@AbstractPlotting.recipe(Bezier) do scene
    merge(
        default_theme(scene, Lines),
        Attributes(
            npoints = 30,
            colorrange = AbstractPlotting.automatic
        )
    )
end

AbstractPlotting.conversion_trait(::Type{<: Bezier}) = AbstractPlotting.PointBased()

function calculated_attributes!(::Type{<: Bezier}, plot)
    color_and_colormap!(plot)
    pos = plot[1][]
    # extend one color per linesegment to be one (the same) color per vertex
    # taken from @edljk  in PR #77
    if haskey(plot, :color) && isa(plot[:color][], AbstractVector) && iseven(length(pos)) && (length(pos) ÷ 2) == length(plot[:color][])
        plot[:color] = lift(plot[:color]) do cols
            map(i-> cols[(i + 1) ÷ 2], 1:(length(cols) * 2))
        end
    end
end

# used in the pipeline too (for poly)
function from_nansep_vec(v::Vector{T}) where T
    idxs = findall(isnan, v)

    if isempty(idxs)
        return [v]
    end
    vs = Vector{Vector{T}}(undef, length(idxs))
    prev = 1
    num = 1
    for i in idxs
        vs[num] = v[prev:i-1]

        prev = i + 1
        num += 1
    end

    return vs
end

function bezier_value(pts::AbstractVector, t::Real)
    val = 0.0
    n = length(pts) - 1
    for (i, p) in enumerate(pts)
        val += p * binomial(n, i - 1) * (1 - t)^(n - i + 1) * t^(i - 1)
    end
    val
end

function to_bezier(p::Vector{Point2f0}, npoints::Int)
    curves = Point2f0[]

    rawvecs = [getindex.(p, n) for n in 1:2]

    for rng in from_nansep_vec(p)
        ts = LinRange(0, 1, npoints)
        xs = map(t -> bezier_value(getindex.(rng, 1), t), ts)
        ys = map(t -> bezier_value(getindex.(rng, 2), t), ts)

        append!(curves, Point2f0.(xs, ys))

        push!(curves, Point2f0(NaN))
    end

    return curves
end

function AbstractPlotting.plot!(plot::Bezier)
    positions = plot[1]

    @extract plot (npoints,)

    curves = lift(to_bezier, positions, npoints)

    lines!(
        plot,
        curves;
        linestyle = plot.linestyle,
        linewidth = plot.linewidth,
        color = plot.color,
        colormap = plot.colormap,
        colorrange = plot.colorrange
    )
end


text:
module PhysCon

using Parameters
#-----------------------------------------------------------------------
# DESCRIPTION:
# Calculate saturation vapor pressure and latent heat of vaporization
#
export physcon
@with_kw struct phys{T<:Number}
    grav::T = 9.80665;               # Gravitational acceleration (m/s2)
    tfrz::T = 273.15;                # Freezing point of water (K)
    σ::T = 5.67e-08;                 # Stefan-Boltzmann constant (W/m2/K4)
    mmdry::T = 28.97 / 1000;         # Molecular mass of dry air (kg/mol)
    mmh2o::T = 18.02 / 1000;         # Molecular mass of water (kg/mol)
    Cpd::T = 1005;                   # Specific heat of dry air at constant pressure (J/kg/K)
    Cpw::T = 1846;                   # Specific heat of water vapor at constant pressure (J/kg/K)
    Rgas::T = 8.31446;               # Universal gas constant (J/K/mol)
    visc_0::T = 13.3e-06;            # Kinematic viscosity at 0C and 1013.25 hPa (m2/s)
    Dh_0::T = 18.9e-06;              # Molecular diffusivity (heat) at 0C and 1013.25 hPa (m2/s)
    Dv_0::T = 21.8e-06;              # Molecular diffusivity (H2O) at 0C and 1013.25 hPa (m2/s)
    Dc_0::T = 13.8e-06;              # Molecular diffusivity (CO2) at 0C and 1013.25 hPa (m2/s)
end

physcon = phys{Float32}()

end

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
#= Generates gradient of the objective function with ReverseDiff
using a representative input x0
=#
function generate_gradient(f, x0)
	cfg_f = ReverseDiff.GradientConfig(x0)

	function grad!(g, x)
		ReverseDiff.gradient!(g, f, x, cfg_f)
	end

	return grad!
end

#= Generates gradient of the objective function with ReverseDiff
using a representative input x0 and compiles the tape.

This is useful if there are no branches in the code. If there
are branches, then compiling the tape could result in incorrect
gradients for certain inputs.
=#
function generate_gradient_tape(f, x0)
	cfg_f = ReverseDiff.GradientConfig(x0)
	tape_f = ReverseDiff.GradientTape(f, x0, cfg_f)
	ctape_f = ReverseDiff.compile(tape_f)

	function grad!(g, x)
		ReverseDiff.gradient!(g, ctape_f, x)
	end

	return grad!
end

#= Generates Jacobian of the constraints with Reversediff
using a representative input x0 and output size m
=#
function generate_jacobian(c!, x0, m::Int)
	tmp_m = similar(x0, m)
	cfg_c = ReverseDiff.JacobianConfig(tmp_m, x0)

	function jac!(Jc, cval, x)
		ReverseDiff.jacobian!(Jc, c!, cval, x, cfg_c)
	end

	return jac!
end

#= Generates Jacobian of the constraints with Reversediff
using a representative input x0 and output size m and compiles
a tape.

This is useful if there are no branches in the code. If there
are branches, then compiling the tape could result in incorrect
gradients for certain inputs.
=#
function generate_jacobian_tape(c!, x0, m::Int)
	tmp_m = similar(x0, m)
	cfg_c = ReverseDiff.JacobianConfig(tmp_m, x0)
	tape_c = ReverseDiff.JacobianTape(c!, tmp_m, x0, cfg_c)
	ctape_c = ReverseDiff.compile(tape_c)

	function jac!(Jc, cval, x)
		ReverseDiff.jacobian!(Jc, ctape_c, x)
		ReverseDiff.extract_result_value!(cval, ReverseDiff.output_hook(ctape_c))
	end

	return jac!
end

#= Generates function to calculate the application of the Hessian
of the Lagrangian on a vector with a representative input of x0
=#
function generate_hess_lag_vec(grad_dual!, jac_dual!, x0, m::Int)
	n = length(x0)

	duals_n = zeros(Dual{nothing,Float64,1}, n)
	duals_m = zeros(Dual{nothing,Float64,1}, m)
	duals_jac = zeros(Dual{nothing,Float64,1}, m, n)
	x_dual = zeros(Dual{nothing,Float64,1}, n)

	function hess_lag_vec!(dest, src, x, λ)
		a = Dual{nothing}(0.0, 1.0)
		x_dual .= x .+ a.*src

		# calculate dual objective function gradient
		grad_dual!(duals_n, x_dual)

		# extract Hessian-vector products
		@inbounds for i = 1:n
			dest[i] = duals_n[i].partials[1]
		end

		if m == 0
			return dest
		end

		# calculate dual constraint Jacobian
		jac_dual!(duals_jac, duals_m, x_dual)

		@inbounds for j = 1:n
			dest[j] += sum(duals_jac[i,j].partials[1]*λ[i] for i in 1:m)
		end

		return dest
	end

	return hess_lag_vec!
end


text:
# ---
# title: 691. Stickers to Spell Word
# id: problem691
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/stickers-to-spell-word/description/>
# hidden: true
# ---
#
# We are given N different types of stickers. Each sticker has a lowercase
# English word on it.
#
# You would like to spell out the given `target` string by cutting individual
# letters from your collection of stickers and rearranging them.
#
# You can use each sticker more than once if you want, and you have infinite
# quantities of each sticker.
#
# What is the minimum number of stickers that you need to spell out the
# `target`? If the task is impossible, return -1.
#
# **Example 1:**
#
# Input:
#
#
#
#     ["with", "example", "science"], "thehat"
#
#
# Output:
#
#
#
#     3
#
#
# Explanation:
#
#
#
#     We can use 2 "with" stickers, and 1 "example" sticker.
#     After cutting and rearrange the letters of those stickers, we can form the target "thehat".
#     Also, this is the minimum number of stickers necessary to form the target string.
#
#
# **Example 2:**
#
# Input:
#
#
#
#     ["notice", "possible"], "basicbasic"
#
#
# Output:
#
#
#
#     -1
#
#
# Explanation:
#
#
#
#     We can't form the target "basicbasic" from cutting letters from the given stickers.
#
#
# **Note:**
#
# * `stickers` has length in the range `[1, 50]`.
# * `stickers` consists of lowercase English words (without apostrophes).
# * `target` has length in the range `[1, 15]`, and consists of lowercase English letters.
# * In all test cases, all words were chosen _randomly_ from the 1000 most common US English words, and the target was chosen as a concatenation of two random words.
# * The time limit may be more challenging than usual. It is expected that a 50 sticker test case can be solved within 35ms on average.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    Options

Various options for NNFEM simulation.

- `save_history::Int64`

  - 0: no histroy is saved
  - 1: save stress, displacement, acceleration history
  - 2: in addition to 1, save internal and external force
"""
mutable struct Options
    save_history::Int64
end

function Options()
    Options(2)
end

options = Options()


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.reduce
# Date: 12/15/2014

module EncounterDBN

export
    AbstractEncounterDBN,
    AddObserver,

    getInitialState,
    initialize,
    update,
    get,

    CorrAEMDBN,
    StarDBN,
    SideOnDBN,
    PairwiseCorrAEMDBN

import CommonInterfaces: initialize, update

using AbstractEncounterDBNImpl
import AbstractEncounterDBNInterfaces: getInitialState, get

using CorrAEMDBNImpl
using StarDBNImpl
#using PairwiseCorrAEMDBNImpl
using SideOnDBNImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
using Base: @info, printstyled, stderr, stdout
using Suppressor
using Test: @testset, @test, @test_throws

@testset "Suppressor" begin

# everything that prints to stdout and stderr should be prefixed with an
# incrementing number so we can check that the output is as expected

@testset "stdout capture" begin
    output = @capture_out begin
        println("CAPTURED STDOUT")
        println(stderr, "01 PRINTED STDERR")
    end
    @test output == "CAPTURED STDOUT\n"
    println("02 PRINTED STDOUT")
end

@testset "stderr capture" begin
    output = @capture_err begin
        println("03 PRINTED STDOUT")
        println(stderr, "CAPTURED STDERR")
    end
    @test output == "CAPTURED STDERR\n"
    println(stderr, "04 PRINTED STDERR")
end

# we're assuming the global context here has color enabled
@testset "disabling color" begin
    printstyled("05 PRINTED GREEN STDOUT\n", color=:green)
    printstyled(stderr, "06 PRINTED GREEN STDERR\n", color=:green)

    @color_output false begin
        printstyled("07 PRINTED NORMAL STDOUT\n", color=:green)
        printstyled(stderr, "08 PRINTED NORMAL STDERR\n", color=:green)
    end

    printstyled("09 PRINTED GREEN STDOUT\n", color=:green)
    printstyled(stderr, "10 PRINTED GREEN STDERR\n", color=:green)
end

@testset "enabling color" begin
    @color_output false begin
        printstyled("11 PRINTED NORMAL STDOUT\n", color=:green)
        printstyled(stderr, "12 PRINTED NORMAL STDERR\n", color=:green)

        @color_output true begin
            printstyled("13 PRINTED GREEN STDOUT\n", color=:green)
            printstyled(stderr, "14 PRINTED GREEN STDERR\n", color=:green)
        end

        printstyled("15 PRINTED NORMAL STDOUT\n", color=:green)
        printstyled(stderr, "16 PRINTED NORMAL STDERR\n", color=:green)
    end
end

@testset "stdout suppression" begin
    @test @suppress_out begin
        println("SUPPRESSED STDOUT")
        println(stderr, "17 PRINTED STDERR")
        42
    end == 42
    println("18 PRINTED STDOUT")
end

@testset "stderr suppression" begin
    @test @suppress_err begin
        println("19 PRINTED STDOUT")
        println(stderr, "SUPPRESSED STDERR")
        42
    end == 42
    println(stderr, "20 PRINTED STDERR")
end

@testset "stderr and stdout suppression" begin
    @test @suppress begin
        println("SUPPRESSED STDOUT")
        println(stderr, "SUPPRESSED STDERR")
        42
    end == 42
    println("21 PRINTED STDOUT")
    println(stderr, "22 PRINTED STDERR")
end

# make sure that things still work after an exception is thrown
@testset "exception cleanup" begin
    try
        @capture_out throw(ErrorException(""))
    catch
    end
    println("23 PRINTED STDOUT")
    println(stderr, "24 PRINTED STDERR")

    try
        @capture_err throw(ErrorException(""))
    catch
    end
    println("25 PRINTED STDOUT")
    println(stderr, "26 PRINTED STDERR")

    try
        @suppress throw(ErrorException(""))
    catch
    end
    println("27 PRINTED STDOUT")
    println(stderr, "28 PRINTED STDERR")

    try
        @suppress_out throw(ErrorException(""))
    catch
    end
    println("29 PRINTED STDOUT")
    println(stderr, "30 PRINTED STDERR")

    try
        @suppress_err throw(ErrorException(""))
    catch
    end
    println("31 PRINTED STDOUT")
    println(stderr, "32 PRINTED STDERR")
end

@test_throws ErrorException @suppress begin
    println("SUPPRESSED STDOUT")
    println(stderr, "SUPPRESSED STDERR")
    error("errors would normally get printed but are caught here by @test_throws")
end

@testset "logging capture" begin
    output = @capture_err @info "CAPTURED LOGINFO"
    # 0.6.2 output:
    if isdefined(Base, :CoreLogging)
        @test output == "[ Info: CAPTURED LOGINFO\n"
    else
        @test output == "\e[1m\e[36mInfo: \e[39m\e[22m\e[36mCAPTURED LOGINFO\n\e[39m"
    end
    @info "33 PRINTED LOGINFO"
end

@testset "logging suppression" begin
    @suppress_err @info "SUPPRESSED LOGINFO"
    @info "34 PRINTED LOGINFO"

    @suppress @info "SUPPRESSED LOGINFO"
    @info "35 PRINTED LOGINFO"
end

@testset "color output exception handling" begin
    @color_output true begin
        try
            @color_output false begin
                throw(Exception())
            end
        catch
        end
        printstyled("36 PRINTED GREEN STDOUT\n", color=:green)
        printstyled("37 PRINTED GREEN STDERR\n", color=:green)
    end
end

end # @testset "Suppressor"


text:
using Test, LinearMaps, LinearAlgebra

@testset "identity/scaling map" begin
    @test_throws ArgumentError LinearMaps.UniformScalingMap(true, -1)
    m = 5
    A = 2 * rand(ComplexF64, (m, m)) .- 1
    B = rand(size(A)...)
    M = @inferred 1 * LinearMap(A)
    N = @inferred LinearMap(B)
    LC = @inferred M + N
    v = rand(ComplexF64, m)
    w = similar(v)
    Id = @inferred LinearMap(I, m)
    @test occursin("$m×$m LinearMaps.UniformScalingMap{Bool}", sprint((t, s) -> show(t, "text/plain", s), Id))
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, m, 2m)
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, (m, 2m))
    @test size(Id) == (m, m)
    @test @inferred isreal(Id)
    @test @inferred issymmetric(Id)
    @test @inferred ishermitian(Id)
    @test @inferred isposdef(Id)
    @test Id * v == v
    @test (2 * M' + 3 * I) * v == 2 * A'v + 3v
    @test (3 * I + 2 * M') * v ≈ 3v + 2 * A'v rtol=2eps()
    @test (2 * M' - 3 * I) * v == 2 * A'v - 3v
    @test (3 * I - 2 * M') * v ≈ 3v - 2 * A'v rtol=2eps()
    @test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v
    @test LinearMap(2 * M' + 0I)' * v ≈ (2 * A')' * v
    for λ in (0, 1, rand()), α in (0, 1, rand()), β in (0, 1, rand()), sz in (m, (m,5))
        Λ = @inferred LinearMap(λ*I, m)
        x = rand(Float64, sz)
        y = ones(Float64, sz)
        mul!(y, Λ, x, α, β)
        @test (@allocated mul!(y, Λ, x, α, β)) == 0
        y = copy(x)
        @inferred mul!(y, Λ, x, α, β)
        @test y ≈ λ * x * α + β * x
    end
    for elty in (Float64, ComplexF64), transform in (identity, transpose, adjoint)
        λ = rand(elty)
        x = rand(m)
        J = @inferred LinearMap(LinearMaps.UniformScalingMap(λ, m))
        @test transform(J) * x == transform(λ) * x
        J = @inferred LinearMap(λ*I, m)
        @test (λ * J) * x == (J * λ) * x == (λ * λ) * x
    end
    X = rand(m, m); Y = similar(X)
    @test mul!(Y, Id, X) == X
    @test Id*X isa LinearMap
    @test X*Id isa LinearMap
    @test Matrix(Id*X) == X
    @test Matrix(X*Id) == X
    @test Matrix(Id + X) == I + X
    @test Matrix(X + Id) == X + I
end


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Benders Decomposition

# **Originally Contributed by**: Shuvomoy Das Gupta

# This notebook describes how to implement Benders decomposition in JuMP, which
# is a large scale optimization scheme.

# We only discuss the classical approach (using loops) here. The approach using
# lazy constraints is showed in the correponding tutorial.

# To illustrate an implementation of the Benders decomposition in JuMP, we apply
# it to the following general mixed integer problem:

# ```math
# \begin{aligned}
# & \text{maximize} \quad &&c_1^T x+c_2^T v \\
# & \text{subject to} \quad &&A_1 x+ A_2 v \preceq b \\
# & &&x \succeq 0, x \in \mathbb{Z}^n \\
# & &&v \succeq 0, v \in \mathbb{R}^p \\
# \end{aligned}
# ```

# where $b \in \mathbb{R}^m$, $A_1 \in \mathbb{R}^{m \times n}$,
# $A_2 \in \mathbb{R}^{m \times p}$ and $\mathbb{Z}$ is the set of integers.
#
# Here the symbol $\succeq$ ($\preceq$) stands for element-wise greater (less)
# than or equal to. Any mixed integer programming problem can be written in the
# form above.

# We want to write the Benders decomposition algorithm for the problem above.
# Consider the polyhedron $\{u \in \mathbb{R}^m| A_2^T u \succeq 0, u \succeq 0\}$.
# Assume the set of vertices and extreme rays of the polyhedron is denoted by
# $P$ and $Q$ respectively.
#
# Assume on the $k$th iteration the subset of vertices of the polyhedron
# mentioned is denoted by $T(k)$ and the subset of extreme rays are denoted by
# $Q(k)$, which will be generated by the Benders decomposition problem below.

# ### Benders decomposition algorithm

# **Step 1 (Initialization)**

# We start with $T(1)=Q(1)=\emptyset$.
# Let $f_m^{(1)}$ be arbitrarily large and $x^{(1)}$ be any non-negative integer
# vector and go to Step 2.

# **Step 2 (Solving the master problem)**

# Solve the master problem, $f_\text{m}^{(k)}$ =

# ```math
# \begin{aligned}
# &\text{maximize} &&\quad t \\
# &\text{subject to} \quad &&\forall \bar{u} \in T(k) \qquad t + (A_1^T \bar{u} - c_1)^T x \leq b^T \bar{u} \\
# & && \forall \bar{y} \in Q(k) \qquad (A_1 ^T \bar{y})^T x \leq b^T \bar{y} \\
# & && \qquad \qquad \qquad \; x \succeq 0, x \in \mathbb{Z}^n
# \end{aligned}
# ```

# Let the maximizer corresponding to the objective value $f_\text{m}^{(k)}$ be
# denoted by $x^{(k)}$. Now there are three possibilities:

# - If $f_\text{m}^{(k)}=-\infty$, i.e., the master problem is infeasible,
#   then the original proble is infeasible and sadly, we are done.
#
# - If $f_\text{m}^{(k)}=\infty$, i.e. the master problem is unbounded above,
#   then we take $f_\text{m}^{(k)}$ to be arbitrarily large and $x^{(k)}$ to be
#   a corresponding feasible solution. Go to Step 3.
#
# - If $f_\text{m}^{(k)}$ is finite, then we collect $t^{(k)}$ and $x^{(k)}$ and
#   go to Step 3.

# **Step 3 (Solving the subproblem and add Benders cut when needed)**

# Solve the subproblem, $f_s(x^{(k)})$ =

# ```math
# \begin{aligned}
#   c_1^T x^{(k)} + & \text{minimize} &&  (b-A_1 x^{(k)})^T u \\
#   & \text{subject to} && A_2^T u \succeq c_2 \\
#   & && u \succeq 0, u \in \mathbb{R}^m
# \end{aligned}
# ```

# Let the minimizer corresponding to the objective value $f_s(x^{(k)})$ be
# denoted by $u^{(k)}$. There are three possibilities:

# - If $f_s(x^{(k)}) = \infty$, the original problem is either infeasible or
#   unbounded. We quit from Benders algorithm and use special purpose algorithm
#   to find a feasible solution if there exists one.
#
# - If $f_s(x^{(k)}) = - \infty$, we arrive at an extreme ray $y^{(k)}$.
#   We add the Benders cut corresponding to this extreme ray $(A_1 ^T y^{(k)})^T x \leq b^T y^{(k)}$
#   to the master problem, i.e., $Q(k+1):= Q(k) \cup \{y^{(k)}\}$. Take $k:=k+1$
#   and go to Step 3.
#
# - If $f_s(x^{(k)})$ is finite, then
#
#  * If $f_s(x^{(k)})=f_m^{(k)}$ we arrive at the optimal solution.
#    The optimum objective value of the original problem is $f_s(x^{(k)})=f_m^{(k)}$,
#    an optimal $x$ is $x^{(k)}$ and an optimal $v$ is the dual values for the
#    second constraints of the subproblem. We are happily done!
#
#  * If $f_s(x^{(k)}) < f_m^{(k)}$ we get an suboptimal vertex $u^{(k)}$. We add
#    the corresponding Benders cut $u_0 + (A_1^T u^{(k)} - c_1)^T x \leq b^T u^{(k)}$
#    to the master problem, i.e., $T(k+1) := T(k) \cup \{u^{(k)}\}$. Take $k:=k+1$
#    and go to Step 3.

# For a more general approach to Bender's Decomposition you can have a look at
# [Mathieu Besançon's blog](https://matbesancon.github.io/post/2019-05-08-simple-benders/).

# ## Data for the problem

# The input data is from page 139, Integer programming by Garfinkel and Nemhauser[[1]](#c1).

c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [1 -3;
     -1 -3]
A2 = [1 -2;
     -1 -1]

M = 1000;

# ## How to implement the Benders decomposition algorithm in JuMP

# There are two ways we can implement Benders decomposition in JuMP:
# - *Classical approach:* Adding the Benders cuts in a loop,  and
# - *Modern approach:* Adding the Benders cuts as lazy constraints.

# The classical approach might be inferior to the modern one, as the solver
# - might revisit previously eliminated solution, and
# - might discard the optimal solution to the original problem in favor of a
#   better but ultimately infeasible solution to the relaxed one.

# For more details on the comparison between the two approaches, see
# [Paul Rubin's blog on Benders Decomposition](http://orinanobworld.blogspot.ca/2011/10/benders-decomposition-then-and-now.html).

# ## Classical Approach: Adding the Benders Cuts in a Loop

# Let's describe the master problem first. Note that there are no constraints,
# which we will added later using Benders decomposition.

# ### Loading the necessary packages

using JuMP
import GLPK
import Test

# ### Master Problem Description

master_problem_model = Model(GLPK.Optimizer);

# ### Variable Definition

@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)

# ### Objective Setting

@objective(master_problem_model, Max, t)
global iter_num = 1

print(master_problem_model)

# Here is the loop that checks the status of the master problem and the
# subproblem and then adds necessary Benders cuts accordingly.

iter_num = 1

while true
    println("\n-----------------------")
    println("Iteration number = ", iter_num)
    println("-----------------------\n")
    println("The current master problem is")
    print(master_problem_model)

    optimize!(master_problem_model)

    t_status = termination_status(master_problem_model)
    p_status = primal_status(master_problem_model)

    if p_status == MOI.INFEASIBLE_POINT
        println("The problem is infeasible :-(")
        break
    end

    (fm_current, x_current) = if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
        (M, M * ones(dim_x))
    elseif p_status == MOI.FEASIBLE_POINT
        (value(t), value.(x))
    else
        error("Unexpected status: $((t_status, p_status))")
    end

    println("Status of the master problem is ", t_status,
            "\nwith fm_current = ", fm_current,
            "\nx_current = ", x_current)

    sub_problem_model = Model(GLPK.Optimizer)

    c_sub = b - A1 * x_current

    local u = @variable(sub_problem_model, u[1:dim_u] >= 0)

    @constraint(
        sub_problem_model,
        constr_ref_subproblem[j = 1:size(A2, 2)],
        A2[:, j]' * u >= c2[j],
    )
    ## The second argument of @constraint macro,
    ## constr_ref_subproblem[j=1:size(A2,2)] means that the j-th constraint is
    ## referenced by constr_ref_subproblem[j].

    @objective(sub_problem_model, Min, c1' * x_current + c_sub' * u)

    print("\nThe current subproblem model is \n")
    print(sub_problem_model)

    optimize!(sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = b' * u_current

    println("Status of the subproblem is ", t_status_sub,
        "\nwith fs_x_current = ", fs_x_current,
        "\nand fm_current = ", fm_current)

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current == fm_current # we are done
        Test.@test value(t) ≈ -4 #hide
        println("\n################################################")
        println("Optimal solution of the original problem found")
        println("The optimal objective value t is ", fm_current)
        println("The optimal x is ", x_current)
                println("The optimal v is ", dual.(constr_ref_subproblem))
        println("################################################\n")
        break
    end

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv = A1' * u_current - c1
        @constraint(master_problem_model, t + cv' * x <= γ)
        println("t + ", cv, "ᵀ x <= ", γ)
    end

    if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1'* u_current
        @constraint(master_problem_model, ce' * x <= γ)
        println(ce, "ᵀ x <= ", γ)
    end

    global iter_num += 1
end

# ### References

# ```@raw html
# <a id='c1'></a>
# ```

# 1. Garfinkel, R. & Nemhauser, G. L. Integer programming. (Wiley, 1972).

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# joins(x) = length(x) > 0 ? join(x, " ") : x
#
# joinsample = (joins ∘ sample)

lastcommand() = "Last command: \"$(replset.commands[end])\" "

Base.nothing(; verbose = true) = verbose && println(lastcommand() * " - " * reportfeasible())

markdownescape(x) = replace(x, "|"=>"\\|") |> (x -> replace(x, "#"=>"\\#"))

printmarkdown(x) = show(stdout, MIME("text/plain"), Markdown.parse(x))

reportfeasible(logicset::LogicalCombo) = "Feasible Outcomes: $(nfeasible(logicset)) \t $(join(sample(logicset),", "))"

reportfeasible() = "Feasible Outcomes: $(nfeasible(replset)) \t $(join(sample(replset),", "))"

tounicode(x) = replace(x, r"\bin\b"=>"∈")


text:
@testset "eBPF" begin

include("definitions/bpf.jl")

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

@testset "No-op" begin
    kernel() = 0

    output = bpf_code_native(kernel, ())[1]
    @test occursin("\tr0 = 0\n\texit", output)
end
@testset "Return argument" begin
    kernel(x) = x

    output = bpf_code_native(kernel, (UInt64,); strip=true)[1]
    @test occursin("\tr0 = r1\n\texit", output)
end
@testset "Addition" begin
    kernel(x) = x+1

    output = bpf_code_native(kernel, (UInt64,); strip=true)[1]
    @test occursin("\tr0 = r1\n\tr0 += 1\n\texit", output)
end
@testset "Errors" begin
    kernel(x) = fakefunc(x)

    @test_throws GPUCompiler.KernelError bpf_code_native(kernel, (UInt64,); strip=true)[1]
end
@testset "Function Pointers" begin
    @testset "valid" begin
        goodcall(x) = Base.llvmcall("%2 = call i64 inttoptr (i64 3 to i64 (i64)*)(i64 %0)\nret i64 %2", Int, Tuple{Int}, x)
        kernel(x) = goodcall(x)

        output = bpf_code_native(kernel, (Int,); strip=true)[1]
        @test occursin("\tcall 3\n\texit", output)
    end
    @testset "invalid" begin
        badcall(x) = Base.llvmcall("%2 = call i64 inttoptr (i64 3000 to i64 (i64)*)(i64 %0)\nret i64 %2", Int, Tuple{Int}, x)
        kernel(x) = badcall(x)

        @test_throws GPUCompiler.InvalidIRError bpf_code_native(kernel, (Int,); strip=true)[1]
    end
end

end


text:
module CUDNN

using CUDAapi
using CUDAapi: libraryPropertyType

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

using ..CuArrays
using ..CuArrays: @argout, @workspace
import ..CuArrays.unsafe_free!

import NNlib

const libcudnn = Ref("libcudnn")

# core library
include("libcudnn_common.jl")
include("error.jl")
include("libcudnn.jl")

# low-level wrappers
include("util.jl")
include("base.jl")
include("tensor.jl")
include("conv.jl")
include("pooling.jl")
include("activation.jl")
include("filter.jl")
include("softmax.jl")
include("batchnorm.jl")
include("dropout.jl")
include("rnn.jl")

# high-level integrations
include("nnlib.jl")

include("compat.jl")

const created_handles = IdDict{CuContext,cudnnHandle_t}()
const active_handles = Vector{Union{Nothing,cudnnHandle_t}}()

function handle()
    tid = Threads.threadid()
    if @inbounds active_handles[tid] === nothing
        ctx = context()
        active_handles[tid] = get!(created_handles, ctx) do
            handle = cudnnCreate()
            atexit(()->CUDAdrv.isvalid(ctx) && cudnnDestroy(handle))
            handle
        end
    end
    @inbounds active_handles[tid]
end

function __init__()
    resize!(active_handles, Threads.nthreads())
    fill!(active_handles, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        # we don't eagerly initialize handles, but do so lazily when requested
        active_handles[tid] = nothing
    end
end

end


text:
export moveispseudo,
    pseudocaptures, pseudochecks, pseudoislegal, pseudomoves, pseudoquiets, pseudoevasions

function is_attacked_after_move(b::Board, s::Square, m::Move)::Bool
    f = from(m)
    t = to(m)
    us = sidetomove(b)
    theirs = pieces(b, -us)

    # Remove captured piece, if any:
    if pieceon(b, t) != EMPTY
        theirs -= t
    elseif moveisep(b, m)
        theirs -= t - (us == WHITE ? DELTA_N : DELTA_S)
    end

    # Test for attacks by non-sliding pieces:
    if !isempty(pawnattacks(us, s) ∩ pawns(b) ∩ theirs)
        return true
    end
    if !isempty(knightattacks(s) ∩ knights(b) ∩ theirs)
        return true
    end
    if !isempty(kingattacks(s) ∩ kings(b) ∩ theirs)
        return true
    end

    # Test for attacks by sliding pieces:
    ours = pieces(b, us) - f + t
    occ = ours ∪ theirs

    if !isempty(bishopattacks(occ, s) ∩ bishoplike(b) ∩ theirs)
        return true
    end
    if !isempty(rookattacks(occ, s) ∩ rooklike(b) ∩ theirs)
        return true
    end

    # Square not attacked!
    return false
end

"""
    moveispseudo(b::Board, m::Move)

Test whether the move `m` is pseudo-legal.

This means that the move `m` is legal except for possibly leaving the king in
check.

# Examples
```julia-repl
julia> b = @startboard e4 c5 Nf3 d6 Bb5;

julia> moveispseudo(b, Move(SQ_G8, SQ_F6))
true

julia> moveispseudo(b, Move(SQ_B8, SQ_C6))
true

julia> moveispseudo(b, Move(SQ_G8, SQ_G6))
false
```
"""
function moveispseudo(b::Board, m::Move)::Bool
    f = from(m)
    t = to(m)
    us = sidetomove(b)
    them = -us

    # The piece on the source square must belong to the current side to move:
    pcolor(pieceon(b, f)) == us || return false

    # The piece on the destination square cannot belong to the current side to
    # move. TODO: Chess960 castling
    pcolor(pieceon(b, t)) ≠ us || return false

    p = pieceon(b, f)
    Δ = t - f

    if ptype(p) == PAWN
        # If the destination square is on the 1st or 8th rank, the move must be
        # a promotion:
        if t ∈ SS_RANK_1 ∪ SS_RANK_8 && !ispromotion(m)
            return false
        end

        # Is this a normal pawn push?
        push = us == WHITE ? DELTA_N : DELTA_S
        if Δ == push
            return pieceon(b, t) == EMPTY
        end

        # Is this a double pawn push?
        if Δ == 2 * push && f ∈ SS_RANK_2 ∪ SS_RANK_7
            return pieceon(b, t) == EMPTY && pieceon(b, f + push) == EMPTY
        end

        # Is this a capture?
        if t ∈ pawnattacks(us, f)
            return pcolor(pieceon(b, t)) == them || t == epsquare(b)
        end

        # Neither a pawn push nor a pawn capture -- cannot be a legal move:
        return false
    else # Non-pawn move
        # Non-pawn  moves cannot be promotions:
        ispromotion(m) && return false

        pt = ptype(p)

        if pt == KNIGHT
            return t ∈ knightattacks(f)
        elseif pt == BISHOP
            return t ∈ bishopattacks(b, f)
        elseif pt == ROOK
            return t ∈ rookattacks(b, f)
        elseif pt == QUEEN
            return t ∈ queenattacks(b, f)
        elseif pt == KING
            if t ∈ kingattacks(f)
                return true
            elseif Δ == 2 * DELTA_E
                # Kingside castling
                return cancastlekingside(b, us) &&
                       !ischeck(b) &&
                       pieceon(b, f + DELTA_E) == EMPTY &&
                       pieceon(b, t) == EMPTY
            elseif Δ == 2 * DELTA_W
                # Queenside castling
                return cancastlequeenside(b, us) &&
                       !ischeck(b) &&
                       pieceon(b, f + DELTA_W) == EMPTY &&
                       pieceon(b, t) == EMPTY &&
                       pieceon(b, f + 3 * DELTA_W) == EMPTY
            else
                return false
            end
        else
            return false
        end
    end
end

"""
    pseudoislegal(b::Board, m::Move)::Bool

Tests whether the pseudo-legal move `m` is legal.

# Examples
```julia-repl
julia> b = @startboard e4 c5 Nf3 d6 Bb5;

julia> pseudoislegal(b, Move(SQ_C8, SQ_D7))
true

julia> pseudoislegal(b, Move(SQ_G8, SQ_F6))
false
```
"""
function pseudoislegal(b::Board, m::Move)::Bool
    f = from(m)
    t = to(m)
    us = sidetomove(b)
    them = -us
    ksq = kingsquare(b, us)
    check = ischeck(b)

    # If we're not in check, any non-king, non-en-passant move of a non-pinned
    # piece is legal:
    if !check && f ∉ pinned(b) && f ≠ ksq && !moveisep(b, m)
        return true
    end

    # Castling moves are legal if we're not in check and both the destination
    # square and the square we're passing over are not under attack by the
    # opponent:
    if moveisshortcastle(b, m)
        return !check && !isattacked(b, f + DELTA_E, them) && !isattacked(b, t, them)
    end
    if moveislongcastle(b, m)
        return !check && !isattacked(b, f + DELTA_W, them) && !isattacked(b, t, them)
    end

    # For king moves, make sure the king is not under attack on the new square:
    if f == ksq
        return !is_attacked_after_move(b, t, m)
    end

    # When in check or when making an en passant capture, check in a slow way
    # that the move does not leave the king in check:
    if check || moveisep(b, m)
        return !is_attacked_after_move(b, ksq, m)
    end

    # If we are here, we are not in check, and we're moving a pinned piece. The
    # move is legal if it's moving along the ray towards/away from the friendly
    # king.
    return t ∈ squaresbetween(ksq, f) || f ∈ squaresbetween(ksq, t)
end

function addmove!(list::MoveList, f::Square, t::Square)
    push!(list, Move(f, t))
end

function addpromotions!(list::MoveList, f::Square, t::Square)
    for p ∈ [QUEEN, ROOK, BISHOP, KNIGHT]
        push!(list, Move(f, t, p))
    end
end

function addpawnmoves!(list::MoveList, target::SquareSet, Δ::SquareDelta)
    for t ∈ target
        addmove!(list, t - Δ, t)
    end
end

function addpawnpromotions!(list::MoveList, target::SquareSet, Δ::SquareDelta)
    for t ∈ target
        addpromotions!(list, t - Δ, t)
    end
end

function knight_pseudomoves(b::Board, list::MoveList, src::SquareSet, target::SquareSet)
    for f ∈ src
        for t ∈ knightattacks(f) ∩ target
            addmove!(list, f, t)
        end
    end
end

function bishop_pseudomoves(b::Board, list::MoveList, src::SquareSet, target::SquareSet)
    for f ∈ src
        for t ∈ bishopattacks(b, f) ∩ target
            addmove!(list, f, t)
        end
    end
end

function rook_pseudomoves(b::Board, list::MoveList, src::SquareSet, target::SquareSet)
    for f ∈ src
        for t ∈ rookattacks(b, f) ∩ target
            addmove!(list, f, t)
        end
    end
end

function queen_pseudomoves(b::Board, list::MoveList, src::SquareSet, target::SquareSet)
    for f ∈ src
        for t ∈ queenattacks(b, f) ∩ target
            addmove!(list, f, t)
        end
    end
end

function king_pseudomoves(b::Board, list::MoveList, src::SquareSet, target::SquareSet)
    for f ∈ src
        for t ∈ kingattacks(f) ∩ target
            addmove!(list, f, t)
        end
    end
end

function pawn_pseudocaptures(b::Board, list::MoveList)
    us = sidetomove(b)
    them = -us
    ps = pawns(b, us)

    if us == WHITE
        # Promotions (captures and non-captures)
        source1 = ps ∩ SS_RANK_7
        addpawnpromotions!(list, pawnshift_ne(source1) ∩ pieces(b, them), DELTA_NE)
        addpawnpromotions!(list, pawnshift_nw(source1) ∩ pieces(b, them), DELTA_NW)
        addpawnpromotions!(list, pawnshift_n(source1) ∩ emptysquares(b), DELTA_N)

        # Non-promotion captures
        source2 = ps ∩ -SS_RANK_7
        addpawnmoves!(list, pawnshift_ne(source2) ∩ pieces(b, them), DELTA_NE)
        addpawnmoves!(list, pawnshift_nw(source2) ∩ pieces(b, them), DELTA_NW)
    else
        # Promotions (captures and non-captures)
        source1 = ps ∩ SS_RANK_2
        addpawnpromotions!(list, pawnshift_se(source1) ∩ pieces(b, them), DELTA_SE)
        addpawnpromotions!(list, pawnshift_sw(source1) ∩ pieces(b, them), DELTA_SW)
        addpawnpromotions!(list, pawnshift_s(source1) ∩ emptysquares(b), DELTA_S)

        # Non-promotion captures
        source2 = ps ∩ -SS_RANK_2
        addpawnmoves!(list, pawnshift_se(source2) ∩ pieces(b, them), DELTA_SE)
        addpawnmoves!(list, pawnshift_sw(source2) ∩ pieces(b, them), DELTA_SW)
    end
end

function pawn_pseudopushes(b::Board, list::MoveList, target::SquareSet)
    us = sidetomove(b)
    ps = pawns(b, us)
    if us == WHITE
        source1 = ps ∩ -SS_RANK_7
        source2 = ps ∩ SS_RANK_2
        target1 = pawnshift_n(source1) ∩ emptysquares(b) ∩ target
        target2 =
            pawnshift_n(pawnshift_n(source2) ∩ emptysquares(b)) ∩ emptysquares(b) ∩ target
        addpawnmoves!(list, target1, DELTA_N)
        addpawnmoves!(list, target2, 2 * DELTA_N)
    else
        source1 = ps ∩ -SS_RANK_2
        source2 = ps ∩ SS_RANK_7
        target1 = pawnshift_s(source1) ∩ emptysquares(b) ∩ target
        target2 =
            pawnshift_s(pawnshift_s(source2) ∩ emptysquares(b)) ∩ emptysquares(b) ∩ target
        addpawnmoves!(list, target1, DELTA_S)
        addpawnmoves!(list, target2, 2 * DELTA_S)
    end
end

function pawn_promotion_pseudopushes(b::Board, list::MoveList, target::SquareSet)
    us = sidetomove(b)
    ps = pawns(b, us)
    if us == WHITE
        source = ps ∩ SS_RANK_7
        target = pawnshift_n(source) ∩ emptysquares(b) ∩ target
        addpawnpromotions!(list, target, DELTA_N)
    else
        source = ps ∩ SS_RANK_2
        target = pawnshift_s(source) ∩ emptysquares(b) ∩ target
        addpawnpromotions!(list, target, DELTA_N)
    end
end

function pseudo_ep_captures(b::Board, list::MoveList)
    us = sidetomove(b)
    if epsquare(b) != SQ_NONE
        for f ∈ pawnattacks(-us, epsquare(b)) ∩ pawns(b, us)
            addmove!(list, f, epsquare(b))
        end
    end
end

function pseudo_castles(b::Board, list::MoveList)
    us = sidetomove(b)
    ksq = kingsquare(b, us)
    if cancastlekingside(b, us)
        if isempty(squaresbetween(ksq, ksq + 3 * DELTA_E) ∩ occupiedsquares(b))
            addmove!(list, ksq, ksq + 2 * DELTA_E)
        end
    end
    if cancastlequeenside(b, us)
        if isempty(squaresbetween(ksq, ksq + 4 * DELTA_W) ∩ occupiedsquares(b))
            addmove!(list, ksq, ksq + 2 * DELTA_W)
        end
    end
end

function non_pawn_moves_to_target(b::Board, list::MoveList, target::SquareSet)
    us = sidetomove(b)
    knight_pseudomoves(b, list, knights(b, us) ∩ -pinned(b), target)
    bishop_pseudomoves(b, list, bishops(b, us), target)
    rook_pseudomoves(b, list, rooks(b, us), target)
    queen_pseudomoves(b, list, queens(b, us), target)
    king_pseudomoves(b, list, kings(b, us), target)
end

function pawn_discovered_checks(b::Board, list::MoveList, src::SquareSet)
    us = sidetomove(b)
    source = src ∩ pawns(b, us) ∩ -filesquares(file(kingsquare(b, -us)))
    push = us == WHITE ? DELTA_N : DELTA_S
    for f ∈ source
        if f + push ∉ (SS_RANK_1 ∪ SS_RANK_8) && pieceon(b, f + push) == EMPTY
            addmove!(list, f, f + push)
            if f ∈ SS_RANK_2 ∪ SS_RANK_7 && pieceon(b, f + 2 * push) == EMPTY
                addmove!(list, f, f + 2 * push)
            end
        end
    end
end

function pawn_plain_checks(b::Board, list::MoveList, src::SquareSet)
    us = sidetomove(b)
    ksq = kingsquare(b, -us)
    source = src ∩ pawns(b, us) ∩ adjacentfilesquares(ksq)
    push = us == WHITE ? DELTA_N : DELTA_S

    for f ∈ source
        if f + push ∉ (SS_RANK_1 ∪ SS_RANK_8) && pieceon(b, f + push) == EMPTY
            if ksq ∈ pawnattacks(us, f + push)
                addmove!(list, f, f + push)
            end
            if f ∈ (SS_RANK_2 ∪ SS_RANK_7) && pieceon(b, f + 2 * push) == EMPTY
                if ksq ∈ pawnattacks(us, f + 2 * push)
                    addmove!(list, f, f + 2 * push)
                end
            end
        end
    end
end

function discovered_checks(b::Board, list::MoveList, source::SquareSet)
    us = sidetomove(b)
    them = -us
    ksq = kingsquare(b, them)
    target = emptysquares(b)

    knight_pseudomoves(b, list, knights(b, us) ∩ source, target)
    bishop_pseudomoves(b, list, bishops(b, us) ∩ source, target)
    rook_pseudomoves(b, list, rooks(b, us) ∩ source, target)
    pawn_discovered_checks(b, list, source)
    for f ∈ kings(b, us) ∩ source
        for t ∈ kingattacks(f) ∩ target
            if t ∉ squaresbetween(f, ksq) && f ∉ squaresbetween(t, ksq)
                addmove!(list, f, t)
            end
        end
    end
end

function plain_checks(b::Board, list::MoveList, source::SquareSet)
    us = sidetomove(b)
    them = -us
    ksq = kingsquare(b, them)
    target = emptysquares(b)
    b_check_sqs = bishopattacks(b, ksq) ∩ target
    r_check_sqs = rookattacks(b, ksq) ∩ target
    q_check_sqs = b_check_sqs ∪ r_check_sqs

    knight_pseudomoves(b, list, knights(b, us) ∩ source, knightattacks(ksq) ∩ target)
    bishop_pseudomoves(b, list, bishops(b, us) ∩ source, b_check_sqs)
    rook_pseudomoves(b, list, rooks(b, us) ∩ source, r_check_sqs)
    queen_pseudomoves(b, list, queens(b, us) ∩ source, q_check_sqs)
    pawn_plain_checks(b, list, source)
end

"""
    pseudocaptures(b::Board, list::MoveList)

Generates all pseudo-legal captures and promotions.
"""
function pseudocaptures(b::Board, list::MoveList)
    target = pieces(b, -sidetomove(b))
    pawn_pseudocaptures(b, list)
    non_pawn_moves_to_target(b, list, target)
    pseudo_ep_captures(b, list)
end

"""
    pseudoquiets(b::Board, list::MoveList)

Generates all pseudo-legal quiet moves.

Captures and promotions will not be included.
"""
function pseudoquiets(b::Board, list::MoveList)
    target = emptysquares(b)
    pawn_pseudopushes(b, list, target)
    non_pawn_moves_to_target(b, list, target)
    pseudo_castles(b, list)
end

"""
    pseudomoves(b::Board, list::MoveList)

Generates all pseudo-legal moves.
"""
function pseudomoves(b::Board, list::MoveList)
    pseudocaptures(b, list)
    pseudoquiets(b, list)
end

"""
    pseudoevasions(b::Board, list::MoveList)

Generates pseudo-legal moves when in check.

Most obviously illegal moves are excluded.
"""
function pseudoevasions(b::Board, list::MoveList)
    us = sidetomove(b)

    # Generate all pseudo-legal moves for the king.
    king_pseudomoves(b, list, kings(b, us), -pieces(b, us))

    # Other moves are only possible if this is not a double check:
    if issingleton(b.checkers)
        ksq = kingsquare(b, us)
        chsq = first(b.checkers)
        target = squaresbetween(ksq, chsq) + chsq

        # Blocking pawn pushes:
        pawn_pseudopushes(b, list, target)
        pawn_promotion_pseudopushes(b, list, target)

        # Pawn captures of checking piece:
        for f ∈ pawnattacks(-us, chsq) ∩ pawns(b, us) ∩ -pinned(b)
            if chsq ∈ (SS_RANK_8 ∪ SS_RANK_1)
                addpromotions!(list, f, chsq)
            else
                addmove!(list, f, chsq)
            end
        end
        if chsq == epsquare(b) - (us == WHITE ? DELTA_N : DELTA_S)
            pseudo_ep_captures(b, list)
        end

        # Remaining pieces: Generate moves to target for non-pinned pieces.
        knight_pseudomoves(b, list, knights(b, us) ∩ -pinned(b), target)
        bishop_pseudomoves(b, list, bishops(b, us) ∩ -pinned(b), target)
        rook_pseudomoves(b, list, rooks(b, us) ∩ -pinned(b), target)
        queen_pseudomoves(b, list, queens(b, us) ∩ -pinned(b), target)
    end
end

"""
    pseudochecks(b::Board, list::MoveList)

Generates all pseudo-legal non-capturing checking moves.
"""
function pseudochecks(b::Board, list::MoveList)
    dc = discovered_check_candidates(b)
    discovered_checks(b, list, dc)
    plain_checks(b, list, pieces(b, sidetomove(b)) ∩ -dc)
end


text:
using ArgParse, DelimitedFiles, JLD

function readScenarioNum()
    s = ArgParseSettings(description = "arg parse here for parallel experimentation")
    @add_arg_table! s begin
        "--topology", "-t"  # topology name
        "--topoindex", "-i"  # index of IP topology
        "--cutoff", "-c"  # failure scenario cutoff
        "--scenarioID", "-o"  # generate new failure scenarios from weibull
    end
    parsed_args = parse_args(s) # the result is a Dict{String,Any}
    topology = parsed_args["topology"]
    topology_index = parsed_args["topoindex"]
    scenario_cutoff = parse(Float64, parsed_args["cutoff"])
    scenario_id = parsed_args["scenarioID"]

    failure_num = 0
    failureFileName =  "./data/topology/$(topology)/IP_topo_$(topology_index)/$(scenario_cutoff)_ip_scenarios_$(scenario_id).jld"
    if isfile(failureFileName)
        data = load(failureFileName)
        failure_num = length(data["IPScenarios"]["prob"])
    else
        failure_num = 0
    end
    println(failure_num)
end

readScenarioNum()

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
using QContinuousMeasurement
using Test

@testset "css" begin

    @testset "LocalDephasing" begin
        params = ModelParameters(Nj=2, kind=0, kcoll=0)
        model = LocalDephasingModel(params)
        state = blockdiag_css(params.Nj)

        @test state isa BlockDiagonalState
        @test coherentspinstate(model).ρ == state.ρ
    end

    @testset "CollectiveDephasing" begin
        params = ModelParameters(Nj=2, kind=0, kcoll=0)
        model = CollectiveDephasingModel(params)
        state = fixedj_css(params.Nj)

        @test state isa FixedjState
        @test coherentspinstate(model).ρ == state.ρ
    end
end

@testset "SpinSqueezed" begin
    params = ModelParameters(Nj=2, kind=0, kcoll=0)
    model = LocalDephasingModel(params)
    state = blockdiag_sss(params.Nj, 0.0)
    css_state = blockdiag_css(params.Nj)

    @test state.ρ ≈ css_state.ρ
    @test squeezing_param(state) ≈ 1.0
    @test state isa BlockDiagonalState

    sq_state = blockdiag_sss(params.Nj, 0.1)
    @test squeezing_param(sq_state) > 1.0

    params = ModelParameters(Nj=2, kind=0, kcoll=0)
    model = CollectiveDephasingModel(params)
    state = fixedj_sss(params.Nj, 0.0)
    css_state = fixedj_css(params.Nj)

    @test state.ρ ≈ css_state.ρ
    @test squeezing_param(state) ≈ 1.0
    @test state isa FixedjState

    sq_state_f = fixedj_sss(params.Nj, 0.1)
    @test squeezing_param(sq_state_f) > 1.0
    @test squeezing_param(sq_state) ≈ squeezing_param(sq_state_f)
end


text:
# OpenCL Memory Object

abstract CLMemObject <: CLObject

#This should be implemented by all subtypes
# type CLMemType <: CLMemObject
#     valid::Bool
#     id::CL_mem
#     ...
# end

Base.pointer(mem::CLMemObject) = mem.id

Base.sizeof(mem::CLMemObject) = begin
    val = Csize_t[0,]
    @check api.clGetMemObjectInfo(mem.id, CL_MEM_SIZE, sizeof(Csize_t),
                                  val, C_NULL)
    return val[1]
end

function _finalize(mem::CLMemObject)
    if !mem.valid
        throw(CLMemoryError("attempted to double free mem object $mem"))
    end
    if mem.id != C_NULL
        @check_release api.clReleaseMemObject(mem.id)
        mem.id = C_NULL
    end
    mem.valid = false
end

context(mem::CLMemObject) = begin
    param = Array(CL_context, 1)
    @check api.clGetMemObjectInfo(mem.id, CL_MEM_CONTEXT,
                                  sizeof(Csize_t), param, C_NULL)
    return Context(param[1], retain=true)
end

let mem_type(m::CLMemObject) = begin
        result = Array(CL_mem_object_type, 1)
        @check api.clGetMemObjectInfo(m.id, CL_MEM_TYPE,
                        sizeof(CL_mem_object_type), result, C_NULL)
        return result[1]
    end

    mem_flags(m::CLMemObject) = begin
        result = Array(CL_mem_flags)
        @check api.clGetMemObjectInfo(m.id, CL_MEM_FLAGS,
                        sizeof(CL_mem_flags), result, C_NULL)
        mf = result[1]
        flags = Symbol[]
        if (mf & CL_MEM_READ_WRITE) != 0
            push!(flags, :rw)
        end
        if (mf & CL_MEM_WRITE_ONLY) != 0
            push!(flags, :w)
        end
        if (mf & CL_MEM_READ_ONLY) != 0
            push!(flags, :r)
        end
        if (mf & CL_MEM_USE_HOST_PTR) != 0
            push!(flags, :use)
        end
        if (mf & CL_MEM_ALLOC_HOST_PTR) != 0
            push!(flags, :alloc)
        end
        if (mf & CL_MEM_COPY_HOST_PTR) != 0
            push!(flags, :copy)
        end
        return tuple(flags...)
    end

    size(m::CLMemObject) = begin
        result = Array(Csize_t, 1)
        @check api.clGetMemObjectInfo(m.id, CL_MEM_SIZE,
                        sizeof(Csize_t), result, C_NULL)
        return result[1]
    end

    reference_count(m::CLMemObject) = begin
        result = Array(CL_uint, 1)
        @check api.clGetMemObjectInfo(m.id, CL_MEM_REFERENCE_COUNT,
                        sizeof(CL_uint), result, C_NULL)
        return result[1]
    end

    map_count(m::CLMemObject) = begin
        result = Array(CL_uint, 1)
        @check api.clGetMemObjectInfo(m.id, CL_MEM_MAP_COUNT,
                        sizeof(CL_uint), result, C_NULL)
        return result[1]
    end

    info_map = Dict{Symbol, Function}(
        :mem_type => mem_type,
        :mem_flags => mem_flags,
        :size => size,
        :reference_count => reference_count,
        :map_count => map_count
    )

    function info(mem::CLMemObject, minfo::Symbol)
        try
            func = info_map[minfo]
            func(mem)
        catch err
            if isa(err, KeyError)
                throw(ArgumentError("OpenCL.MemObject has no info for: $minfo"))
            else
                throw(err)
            end
        end
    end
end

#TODO: enqueue_migrate_mem_objects(queue, mem_objects, flags=0, wait_for=None)
#TODO: enqueue_migrate_mem_objects_ext(queue, mem_objects, flags=0, wait_for=None)


text:
using StackOverflow, HTTP

"""
    getsuggested_edits(sort::string = "creation", order::string = "desc")

    Prints 30 suggested edits in descending(default)/ascending order of their creation_date(default)/approval_date/rejection_date in the systems.
"""

function getsuggested_edits(sort::String = "creation", order::String = "desc" )

    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/suggested-edits?order=$(order)&sort=$(sort)&site=stackoverflow")

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj

        if occursin("items", k)
            for item in v
                for (key, value) in item
                    println("$(key) : $(value)")
                end
                println("\n")
            end
        end

    end

end