text:
function evalPsiFromBEM(BEMesh,quadRule,f,q,xr,xz; derivative=0)
#=
Given a BE contour BEMesh, a quadrature rule on this contour, the Dirichlet P1 dof and the Neumann P0 dof,
evaluate psi at any point (xr,xz) outside the contour resulting from BEM
q is a vector of size n edges
f is a vector of size n nodes (the BEMesh contour is closed)
=#

	#constant
	mu0=4*pi*10^-7; mu0inv=1/mu0

	# dimension
	n=size(BEMesh.edges,1)

	#quadrature rule
	xGauss=quadRule.x
	wGauss=quadRule.w
	nQuad=size(quadRule.w,1)

	#in order to remap the quadrature we need the midpoints and length of edges
	e=collect(1:n)
	edgeToNode=BEMesh.edges
	h=BEMesh.length

    P=0.5*(BEMesh.coordinates[edgeToNode[:,1],:]+BEMesh.coordinates[edgeToNode[:,2],:])
	#@show P

	#normal
	normal=BEMesh.normalsIn

	#quadrature points on all edges
	PG=zeros(size(P,1),2,nQuad)
	for i in 1:nQuad
		PG[:,:,i]=P + xGauss[i]*0.5*(BEMesh.coordinates[edgeToNode[:,2],:]-BEMesh.coordinates[edgeToNode[:,1],:])
	end

	#evaluate Green function at all y=Gauss points PG and x=(xr,xz)
	#need to broadcast on the number of edges

	DE=zeros(n,2); DI=zeros(n,2); DJ=zeros(n,2)
	NE=zeros(n); NI=zeros(n); NJ=zeros(n)

	dxrDE=zeros(n,2);
	dxzDE=zeros(n,2);
	dxrNE=zeros(n)
	dxzNE=zeros(n)

	for i in 1:nQuad
		#Evaluate Green function and first derivative (in y=B)
		G = @. Green(xr,xz,PG[:,1,i],PG[:,2,i],derivative=derivative+1)
		#G is a vector of named tuples, we need the broadcasted getindex function to access a value (with .() or @. see below)
		#@show G

		#Green function value
		G_PXPG = getindex.(G,:G_xy)
		#@show G_PXPG

		#Neumann trace
		dnG_PXPG = @.getindex(G,:dyrG_xy) .* normal[:,1] +  @.getindex(G,:dyzG_xy) .* normal[:,2]
		#@show dnG_P1GP2G

		# Matrices. Matrix.vector product to be summed afterwards for integral values
		# for psi(x)
		NE=NE+0.5*wGauss[i]*h.*G_PXPG

		#DE=DE-0.5*wGauss[i]*[ h.*dnG_PXPG*(1-(0.5+xGauss[i]*0.5))  h.*dnG_PXPG*((0.5+xGauss[i]*0.5)) ]
		DE=DE-0.5*wGauss[i]*mu0inv*[ h.*(1.0./PG[:,1,i]).*dnG_PXPG*(1-(0.5+xGauss[i]*0.5))  h.*(1.0./PG[:,1,i]).*dnG_PXPG*((0.5+xGauss[i]*0.5)) ]

		# for d_xr_psi(x) and d_xz_psi(x)
		if(derivative>0)
			dxrG= @.getindex(G,:dxrG_xy);
			dxzG= @.getindex(G,:dxzG_xy);
			dxrNE=dxrNE+0.5*wGauss[i]*h.*dxrG
			dxzNE=dxzNE+0.5*wGauss[i]*h.*dxzG

			dxrdyn = @.getindex(G,:dxrdyrG_xy) .* normal[:,1] + @.getindex(G,:dxrdyzG_xy) .* normal[:,2]
			dxzdyn = @.getindex(G,:dxzdyrG_xy) .* normal[:,1] + @.getindex(G,:dxzdyzG_xy) .* normal[:,2]

			dxrDE=dxrDE-0.5*wGauss[i]*mu0inv*[ h.*(1.0./PG[:,1,i]).*dxrdyn*(1-(0.5+xGauss[i]*0.5))  h.*(1.0./PG[:,1,i]).*dxrdyn*((0.5+xGauss[i]*0.5)) ]
			dxzDE=dxzDE-0.5*wGauss[i]*mu0inv*[ h.*(1.0./PG[:,1,i]).*dxzdyn*(1-(0.5+xGauss[i]*0.5))  h.*(1.0./PG[:,1,i]).*dxzdyn*((0.5+xGauss[i]*0.5)) ]
		end

	end

	NI=e; NJ=e
	DI=[e e]; DJ=edgeToNode

	DI=reshape(DI,size(DI,1)*size(DI,2))
	DJ=reshape(DJ,size(DJ,1)*size(DJ,2))
	DE=reshape(DE,size(DE,1)*size(DE,2))

	Dspmat=sparse(DI,DJ,DE,n,n)
	Nspmat=sparse(NI,NJ,NE,n,n)
	Ival=Dspmat*f+Nspmat*q
	Ival=sum(Ival)

	if(derivative==0)
		return (psi=Ival,)
	end

	if(derivative>0)
		dxrDE=reshape(dxrDE,size(dxrDE,1)*size(dxrDE,2))
		dxzDE=reshape(dxzDE,size(dxzDE,1)*size(dxzDE,2))
		dxrDspmat=sparse(DI,DJ,dxrDE,n,n)
		dxzDspmat=sparse(DI,DJ,dxzDE,n,n)
		dxrNspmat=sparse(NI,NJ,dxrNE,n,n)
		dxzNspmat=sparse(NI,NJ,dxzNE,n,n)

		dxrIval=dxrDspmat*f+dxrNspmat*q
		dxrIval=sum(dxrIval)

		dxzIval=dxzDspmat*f+dxzNspmat*q
		dxzIval=sum(dxzIval)

		return (psi=Ival, drpsi=dxrIval, dzpsi=dxzIval)
	end

end


text:
@testset "deprecations.jl" begin
    @test (@test_deprecated VarName(:x)) == VarName{:x}()
    @test (@test_deprecated VarName(:x, ((1,), (:, 2)))) == VarName{:x}(((1,), (:, 2)))
end

ext: jl


text:
@reexport module ExampleSystems

using Distributions
using StatsBase
using DynamicalSystemsBase, DelayEmbeddings

include("discretemaps/ar1.jl")
include("discretemaps/ar1_bidir.jl")
include("discretemaps/anishchenko1.jl")
include("discretemaps/henon2.jl")
include("discretemaps/henontriple.jl")
include("discretemaps/ikeda.jl")
include("discretemaps/linearmap.jl")
include("discretemaps/nonlinear3D_linear_and_nonlinear_coupling.jl")
include("discretemaps/nontrivial_pegiun.jl")
include("discretemaps/logistic2_unidir.jl")
include("discretemaps/logistic2_bidir.jl")
include("discretemaps/logistic3.jl")
include("discretemaps/logistic4.jl")
include("discretemaps/ulammap.jl")
include("discretemaps/var1.jl")
include("discretemaps/verdes.jl")

include("continuous_systems/chuacircuits_driven.jl")
include("continuous_systems/chuacircuit_nscroll_sine.jl")
include("continuous_systems/hindmarsh_rose.jl")
include("continuous_systems/mediated_link.jl")
include("continuous_systems/lorenz_lorenz_bidir.jl")
include("continuous_systems/lorenz_lorenz_lorenz_bidir_forced.jl")
include("continuous_systems/lorenz_lorenz_lorenz_transitive.jl")
include("continuous_systems/repressilator.jl")
include("continuous_systems/rossler_rossler_bidir.jl")
include("continuous_systems/rossler_rossler_rossler_bidir_forced.jl")
include("continuous_systems/rossler_lorenz.jl")

include("noise.jl")

end


text:
Base.@kwdef mutable struct Parameters
    mode::String = "Normal"      # If Debug it will allow printing some useful information including collecting values for analysis parameters.
    algorithm::String = "SQP-TR" # SQP-TR: sequential quadratic programming with trust region

    # Defines the external solver for suproblems
    external_optimizer::Union{Nothing,DataType,MOI.OptimizerWithAttributes,Function} =
        nothing

    # Whether to use approximation hessian (limited-memory), exact, or none
    hessian_type::String = "none"

    # flags
    OutputFlag::Int = 1   # 0 supresses all outputs except warnings and errors
    StatisticsFlag::Int = 0   # 0 supresses collection of statistics parameters

    # Algorithmic parameters
    tol_direction::Float64 = 1.e-6  # tolerance for the norm of direction
    tol_residual::Float64 = 1.e-6   # tolerance for Kuhn-Tucker residual
    tol_infeas::Float64 = 1.e-6     # tolerance for constraint violation
    max_iter::Int = 1000            # Defines the maximum number of iterations
    time_limit::Float64 = Inf       # Defines the time limit for the solver. (This hasn't been implemented yet)
    init_mu::Float64 = 1.e+4        # initial mu value
    max_mu::Float64 = 1.e+10        # maximum mu value allowed
    rho::Float64 = 0.8              # parameter in (0,1) used for updating merit function penalty
    eta::Float64 = 0.4              # descent step test parameter defined in (0,0.5)
    tau::Float64 = 0.9              # line search step decrease parameter defined in (0,1)
    min_alpha::Float64 = 1.e-6      # minimum step size
    tr_size::Float64 = 10.          # trust region size
    use_soc::Bool = false           # use second-order correction
end

function get_parameter(params::Parameters, pname::String)
    return getfield(params, Symbol(pname))
end

function set_parameter(params::Parameters, pname::String, val)
    setfield!(params, Symbol(pname), val)
    return nothing
end


text:
module _TestIndirectImportsPing

using IndirectImports
@indirect import _TestIndirectImportsPong: pong

@indirect function ping end
@indirect pong(x::Integer) = x < 2 ? x : ping(x - 1) + pong(x - 2)

end # module


text:
module Writers

export AbstractWriter, NetCDFWriter, JLD2Writer, write_data, full_name

abstract type AbstractWriter end

"""
    write_data(
        writer,
        filename,
        dims,
        varvals,
    )

Writes the specified dimension names, dimensions, axes, variable names
and variable values to a file. Specialized by every `Writer` subtype.

# Arguments:
# - `writer`: instance of a subtype of `AbstractWriter`.
# - `filename`: into which to write data (without extension).
# - `dims`: Dict of dimension name to axis.
# - `varvals`: Dict of variable name to array of values.
"""
function write_data end

"""
    full_name(
        writer,
        filename
    )

Appends the appropriate (based on `writer`) extension to the specified
filename.
"""
function full_name(writer::AbstractWriter, filename::String)
    if writer isa NetCDFWriter
        return string(filename, ".nc")
    elseif writer isa JLD2Writer
        return string(filename, ".jld2")
    else
        error("Unrecognized Writer $(writer)")
    end
end

include("netcdf_writer.jl")
include("jld2_writer.jl")

end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
export pressure!, pressure,
       pressure_FD!, pressure_FD,
       pressure_AD,
       measure_state

function measure_state(state, t, config::VortexConfig; withfreestream::Bool=false)
    if withfreestream == false
        return pressure(config.ss, state_to_lagrange(state, config), t)
    else
        freestream = Freestream(config.U)
        return pressure(config.ss, state_to_lagrange(state, config), freestream, t)
    end
end

function pressure!(press, targetvels, sourcevels, target, source, t)
    source = deepcopy(source)

    reset_velocity!(sourcevels)
    reset_velocity!(targetvels)

    # Compute the self-induced velocity of the system
    self_induce_velocity!(sourcevels, source, t)

    # Compute the induced velocity on the target elements
    induce_velocity!(targetvels, target, source, t)

    #Only the vortices contribute to the unsteady term
    press .= -real.(convective_complexpotential(target, source, sourcevels)) -0.5*abs2.(targetvels)

    return press
end

pressure(target, source, t) = pressure!(zeros(Float64, length(target)),
                                                        allocate_velocity(target),
                                                        allocate_velocity(source),
                                                        target, source, t)

# Version of the pressure calculation with freestream
function pressure!(press, targetvels, sourcevels, target, source, freestream, t)
    source = deepcopy(source)

    reset_velocity!(sourcevels)
    reset_velocity!(targetvels)

    # Compute the self-induced velocity of the system
    self_induce_velocity!(sourcevels, source, t)
    induce_velocity!(sourcevels, source, freestream, t)
    # Compute the induced velocity on the target elements
    induce_velocity!(targetvels, target, (source, freestream), t)

    #Only the vortices contribute to the unsteady term
    press .= -real.(convective_complexpotential(target, source, sourcevels))
    press .+= -0.5*abs2.(targetvels)

    return press
end

pressure(target, source, freestream, t) = pressure!(zeros(Float64, length(target)),
                                                        allocate_velocity(target),
                                                        allocate_velocity(source),
                                                        target, source, freestream, t)

function pressure_FD!(press, targetvels, targetϕ, sourcevels, target, source, t, Δt)
    source = deepcopy(source)

    reset_velocity!(sourcevels)
    reset_velocity!(targetvels)

    # Compute the self-induced velocity of the system
    self_induce_velocity!(sourcevels, source, t)

    # Compute the induced velocity on the target elements
    induce_velocity!(targetvels, target, source, t)

    targetϕ .= real.(complexpotential(target, source))

    # advective term
    fill!(press, 0.0)
    press .= -0.5*abs2.(targetvels)

    # unsteady term
    advect!(source, source, sourcevels, Δt)
    press .+= (targetϕ - real.(complexpotential(target, source)))/Δt

    return press
end

pressure_FD(target, source, t, Δt) = pressure_FD!(zeros(Float64, length(target)),
                                                        allocate_velocity(target),
                                                        zeros(Float64, length(target)),
                                                        allocate_velocity(source),
                                                        target, source, t, Δt)

function pressure_AD(target, source, t)

    source = deepcopy(source)

    sourcevels = zeros(Complex{Elements.property_type(eltype(source))},length(source))
    self_induce_velocity!(sourcevels, source, t)

    targetvels = zeros(Complex{Elements.property_type(eltype(source))},length(target))
    induce_velocity!(targetvels, target, source, t)

    press = zeros(Complex{Elements.property_type(eltype(source))}, length(target))

    # Unsteady term (only the vortices contribute)
    convective_complexpotential!(press, target, source, sourcevels)
    press .= -real(press)

    # Convective term
    press .-= 0.5*abs2.(targetvels)

    return press
end

function pressure_AD(target, source, freestream, t)

    source = deepcopy(source)

    sourcevels = zeros(Complex{Elements.property_type(eltype(source))},length(source))
    self_induce_velocity!(sourcevels, source, t)
    induce_velocity!(sourcevels, source, freestream, t)

    targetvels = zeros(Complex{Elements.property_type(eltype(source))},length(target))
    induce_velocity!(targetvels, target, (source, freestream), t)

    press = zeros(Complex{Elements.property_type(eltype(source))}, length(target))

    # Unsteady term (only the vortices contribute)
    convective_complexpotential!(press, target, source, sourcevels)
    press .= -real(press)

    # Convective term
    press .-= 0.5*abs2.(targetvels)

    return press
end


text:
include("kiwiComposer.jl")
import kiwiComposer

const doc = """
kiwiComposer

Usage:
    compose.jl FILENAME [-o OUTPUTNAME]
"""
const version = v"0.1"

using DocOpt

arguments = docopt(doc, version=version)
print(arguments)

targetFile = open(arguments["FILENAME"])
setting = readline(targetFile)
sheet = readstring(targetFile)
close(targetFile)

unit = 1.0
eval(parse(setting))
@show setting
print(sheet)

if arguments["-o"]
    kiwiComposer.composer(sheet,arguments["OUTPUTNAME"],unit=unit)
else
    name = match(r"(?<name>.*?)\.\w?", arguments["FILENAME"])[:name]
    kiwiComposer.composer(sheet,name*".wav",unit=unit)
end

max_stars_repo_name: ZaneMuir/kiwi-composer

max_issues_repo_name: ZaneMuir/kiwi-composer

max_forks_repo_name: ZaneMuir/kiwi-composer


text:
const type_prior = ones(4) / 4
const type_transition =
    [0.1 * ones(2, 2) 0.4 * ones(2, 2);
     0.4 * ones(2, 2) 0.1 * ones(2, 2);
     type_prior']
const directions = [0 -1; 0 1; -1 0; 1 0]

generate_points(start, speed, dir, n) =
    start .+ speed .* (dir * collect(1:n)')

@gen function stroke_step(t::Int, prev_state::Tuple)
    prev_stroke, prev_points = prev_state
    n_points ~ shifted_poisson(10, 1)
    speed ~ gamma(2.0, 0.05)
    disconnect ~ bernoulli(0.1)
    if disconnect
        stroke ~ categorical(type_prior)
        start_x ~ uniform(0.0, 5.0)
        start_y ~ uniform(0.0, 5.0)
        start = [start_x, start_y]
    else
        stroke ~ categorical(type_transition[prev_stroke,:])
        start = prev_points[:, end]
    end
    points = generate_points(start, speed, directions[stroke, :], n_points)
    return (stroke, points)
end

@gen (static) function get_points(state::Tuple)
    return collect(eachcol(state[2]))
end

@gen function observe_point(point::AbstractVector)
    return pt ~ broadcasted_normal(point, [0.05, 0.05])
end

@gen (static) function curve_model(T::Int)
    n_strokes ~ poisson(10)
    init_x ~ uniform(0.0, 5.0)
    init_y ~ uniform(0.0, 5.0)
    states = {:strokes} ~ Unfold(stroke_step)(n_strokes, (5, [init_x, init_y]))
    points_per_stroke ~ Map(get_points)(states)
    points = reduce(vcat, points_per_stroke)
    points_observed = points[1:min(T,length(points))]
    obs ~ Map(observe_point)(points_observed)
    return obs
end

@load_generated_functions()

export curve_model


text:
module Risbo

using ArgCheck

using WignerFunctions: Index, initial_wigner
using OffsetArrays: OffsetArray

const Cache = Dict{Index, T} where T <: AbstractFloat

nosym(β::AbstractFloat, index::Index) = nosym(β, index, Cache{typeof(β)}())
function nosym(β::AbstractFloat)
    cache = Cache{typeof(β)}()
    (l::Int64, m₁::Int64, m₂::Int64) -> nosym(β, Index((l, m₁, m₂)), cache)
end
nosym(β::Integer) = nosym(float(β))

function nosym(β::AbstractFloat, index::Index, cache::Cache)
    @argcheck 0 <= index.l
    @argcheck abs(index.m₁) <= index.l
    @argcheck abs(index.m₂) <= index.l

    index ∈ keys(cache) && return cache[index]

    if index.l <= 1
        result = initial_wigner(convert(valtype(cache), β), index)
        cache[index] = result
    else
        nosym_recurrence(convert(valtype(cache), β), index.l, cache)
        result = cache[index]
    end

    return result
end

function nosym_recurrence(β::AbstractFloat, l::Integer, cache::Cache)
    coshb = -cos(β / 2);
    sinhb = sin(β / 2);
    T = promote_type(typeof(β), valtype(cache))

    intermediate = fill!(OffsetArray{T}(undef, -l:l,  -l:l), 0)

    for k ∈ -(l - 1):(l - 1), i ∈ -(l - 1):(l - 1)
        rec = nosym(β, Index((l - 1, k, i)), cache) / (2l - 1)
        intermediate[i, k] += √((l - i) * (l - k)) * rec * coshb
        intermediate[i + 1, k] -= √((l + i) * (l - k)) * rec * sinhb
        intermediate[i, k + 1] += √((l - i) * (l + k)) * rec * sinhb
        intermediate[i + 1, k + 1] += √((l + i) * (l + k)) * rec * coshb
    end

    for k ∈ -l:l, i ∈ -l:l
        cache[Index((l, k, i))] = 0
    end

    for k ∈ -l:(l - 1), i ∈ -l:(l - 1)
        rec = intermediate[i + 1, k + 1] / 2l
        cache[Index((l, k, i))] += √((l - i) * (l - k)) * rec * coshb
        cache[Index((l, k, i + 1))] -= √((l + i + 1) * (l - k)) * rec * sinhb
        cache[Index((l, k + 1, i))] += √((l - i) * (l + k + 1)) * rec * sinhb
        cache[Index((l, k + 1, i + 1))] += √((l + i + 1) * (l + k + 1)) * rec * coshb
    end
end

half(β::AbstractFloat, index::Index) = half(β, index, Cache{typeof(β)}())
function half(β::AbstractFloat)
    cache = Cache{typeof(β)}()
    (l::Int64, m₁::Int64, m₂::Int64) -> half(β, Index((l, m₁, m₂)), cache)
end
half(β::Integer) = half(float(β))

function half(β::AbstractFloat, index::Index, cache::Cache)
    @argcheck 0 <= index.l
    @argcheck abs(index.m₁) <= index.l
    @argcheck abs(index.m₂) <= index.l

    if index.m₂ <= 0
        index ∈ keys(cache) && return cache[index]

        if index.l <= 1
            cache[index] = initial_wigner(convert(valtype(cache), β), index)
        else
            half_recurrence(convert(valtype(cache), β), index.l, cache)
            cache[index]
        end
    elseif (index.m₁ + index.m₂) % 2 == 0
        half(β, Index((index.l, -index.m₁, -index.m₂)), cache)
    else
        -half(β, Index((index.l, -index.m₁, -index.m₂)), cache)
    end
end

function half_recurrence(β::AbstractFloat, l::Integer, cache::Cache)
    coshb = -cos(β / 2);
    sinhb = sin(β / 2);
    T = promote_type(typeof(β), valtype(cache))

    intermediate = fill!(OffsetArray(T, -l:l,  -l:l), 0)

    for k ∈ -(l - 1):(l - 1), i ∈ -(l - 1):1
        rec = half(β, Index((l - 1, k, i)), cache) / (2l - 1)
        intermediate[i, k] += √((l - i) * (l - k)) * rec * coshb
        intermediate[i + 1, k] -= √((l + i) * (l - k)) * rec * sinhb
        intermediate[i, k + 1] += √((l - i) * (l + k)) * rec * sinhb
        intermediate[i + 1, k + 1] += √((l + i) * (l + k)) * rec * coshb
    end

    for k ∈ -l:l, i ∈ -l:0
        cache[Index((l, k, i))] = 0
    end

    for k ∈ -l:(l - 1), i ∈ -l:0
        rec = intermediate[i + 1, k + 1] / 2l
        cache[Index((l, k, i))] += √((l - i) * (l - k)) * rec * coshb
        cache[Index((l, k, i + 1))] -= √((l + i + 1) * (l - k)) * rec * sinhb
        cache[Index((l, k + 1, i))] += √((l - i) * (l + k + 1)) * rec * sinhb
        cache[Index((l, k + 1, i + 1))] += √((l + i + 1) * (l + k + 1)) * rec * coshb
    end
end

end


text:
module ProcessingStd

import GLFW
using ModernGL
using GLUtil
using GLText
using GLWindow
using Events
using Color
using Tau

include("constants.jl")

export animate, coordSystem
export size, Height, Width, displayHeight, displayWidth
export focused, cursor, noCursor
# export frameCount, frameRate
# export createShape, loadShape
export arc, ellipse, line, point, quad, rect, triangle
export bezier
# export bezierDetail, bezierPoint, bezierTangent, curve, curveDetail, curvePoint, curveTangent, curveTightness
export strokeWeight, strokeCap, strokeJoin, noSmooth, smooth
export ellipseMode, rectMode
# export beginShape, endShape, vertex
# export bezierVertex, curveVertex, quadraticVertex
# export shape, shapeMode
export mouseButton, mouseClicked, mouseDragged, mouseMoved, mousePressed
export mouseReleased, mouseWheel, mouseX, mouseY, pmouseX, pmouseY
export background, fill, noFill, colorMode, noStroke, stroke
# export applyMatrix
export popMatrix, printMatrix, pushMatrix, resetMatrix, rotate
export scale, shearX, shearY, translate
export alpha, blue, brightness, color, green, hue, lerpColor, red, saturation
# export createImage
export image, loadImage, noTint
# export imageMode, tint, requestImage
# export texture, textureMode, textureWrap
# export blend, copy, filter, get, loadPixels, set, updatePixels
# export blendMode
# export createGraphics
export createFont, text, textFont
# export loadFont
export textSize, textWidth
# export textAlign, textLeading, textMode
# export textAscent, textDescent

# state variables

# bgCol = RGB(0.94, 0.92, 0.9) #pleasant default background color from Jeff Bezanson's Fractals iJulia notebook (http://nbviewer.ipython.org/url/beowulf.csail.mit.edu/18.337/fractals.ipynb)
# fillStuff = true #should drawing elements be filled or not?
# fillCol = RGB(0.7,0.7,0.7) #color to be used for filling elements
# strokeStuff = true #should drawing elements be stroked or not?
# strokeCol = RGB(0,0,0) #color to be used for stroking elements
# h = 275 #height of display window
# w = 275 #width of display window
# left = -1 #scaling of left-hand of x-axis in plotting coordinate system
# right = 1 #scaling of right-hand of x-axis in plotting coordinate system
# top = 1 #scaling of top-hand of x-axis in plotting coordinate system
# bottom = -1 #scaling of bottom-hand of x-axis in plotting coordinate system
# cMode = "RGB" #what color space are we using for colors?
# title = "Processing.jl" #window title
# ellipseMode = CENTER #specify ellipses center, width, and height when drawing it by default
# rectMode = CORNER #specficy rects upper-left corner coordinate and its width and height when drawing it by default
# cmX
# cmY
# pmX
# pmY
# mouse1Pressed
# mouse2Pressed
# mouse3Pressed
# mouse1Released
# mouse2Released
# mouse3Released
# mouse1Dragged

# state structure

type stateStruct
    bgCol::ColorValue
    fillStuff::Bool
    fillCol::ColorValue
    strokeStuff::Bool
    strokeCol::ColorValue
    h::Int
    w::Int
    left::Float32
    right::Float32
    top::Float32
    bottom::Float32
    cMode::String
    title::String
    tintStuff::Bool
    ellipseMode::Int
    rectMode::Int
    cmX::Int
    cmY::Int
    pmX::Int
    pmY::Int
    mouse1Pressed::Bool
    mouse2Pressed::Bool
    mouse3Pressed::Bool
    mouse1Released::Bool
    mouse2Released::Bool
    mouse3Released::Bool
    mouse1Dragged::Bool
end

# initialize state structure and open drawing window
state = stateStruct(RGB(0.94, 0.92, 0.9), true, RGB(0.7,0.7,0.7), true, RGB(0,0,0), 275, 275, -1., 1., 1., -1., "RGB", "Processing.jl", false, CENTER, CORNER, 0, 0, 0, 0, false, false, false, true, true, true, false)
export state

function size(w, h)
    state.w = w
    state.h = h

    GLFW.Init()

    @osx_only begin
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
    end

    state.window = GLFW.CreateWindow(state.w, state.h, "Processing.jl")
    GLFW.MakeContextCurrent(state.window)
end

function size()
    GLFW.Init()

    @osx_only begin
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
    end

    state.window = GLFW.CreateWindow(state.w, state.h, "Processing.jl")
    GLFW.MakeContextCurrent(state.window)
end

# special Processing.jl animate() command for smoother animations
function animate()
    GLFW.SwapBuffers(state.window)
    GLFW.PollEvents()
end

# allow user to control coordinate system
function coordSystem(left, right, top, bottom)
    Base.Graphics.set_coords(cr, 0, 0, state.w, state.h, left, right, top, bottom)
end

# Environment

# exported environment variables
displayHeight = tcl("winfo", "screenwidth", win)
displayWidth = tcl("winfo", "screenheight", win)

function cursor(cursorType)
    if cursorType == ARROW

    elseif cursorType == CROSS

    elseif cursorType == HAND

    elseif cursorType == MOVE

    elseif cursorType == TEXT

    elseif cursorType == WAIT

    end
end

function focused()
    if
        return false
    else
        return true
    end
end

#frameCount
#frameRate

function Height()
	(pixelwidth::Float64, pixelheight::Float64) = GLFW.GetFramebufferSize(window)
    return pixelheight
end

function noCursor()

end

function Width()
	(pixelwidth::Float64, pixelheight::Float64) = GLFW.GetFramebufferSize(window)
    return pixelwidth
end

# Shape

#createShape
#loadShape

## 2D Primitives

function arc(xcent, ycent, ellipseW, ellipseH, angle1, angle2, mode)
    save(cr)
    move_to(cr, xcent, ycent)
    scale(cr, ellipseW, ellipseH)
    new_sub_path(cr)
    if mode == OPEN
        arc(cr, 0, 0, 1, angle1, angle2)
    elseif mode == PIE
    end
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
    if state.fillStuff
        set_source(cr, state.fillCol)
        fill(cr)
    end
    if mode == CHORD
        move_to(cos(angle1), sin(angle1))
        line_to(cos(angle2), sin(angle2))
    end
    restore(cr)
end

function ellipse(xcent, ycent, ellipseW, ellipseH)
    save(cr)
    if state.ellipseMode == RADIUS
        move_to(cr, xcent, ycent)
        scale(cr, ellipseW/2, ellipseH/2)
    elseif state.ellipseMode == CENTER
        move_to(cr, xcent, ycent)
        scale(cr, ellipseW, ellipseH)
    elseif state.ellipseMode == CORNER
        move_to(cr, xcent+ellipseW/2, ycent+ellipseH/2)
        scale(cr, ellipseW, ellipseH)
    elseif state.ellipseMode == CORNERS
        move_to(cr, xcent, ycent)
        scale(cr, ellipseW, ellipseH)
    end
    new_sub_path(cr)
    arc(cr, 0, 0, 1, 0, 2*pi)
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
    if state.fillStuff
        set_source(cr, state.fillCol)
        fill(cr)
    end
    restore(cr)
end

function line(x1, y1, x2, y2)
    move_to(cr,x1,y1)
    line_to(cr,x2,y2)
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke(cr)
    end
end

function point(x, y)
    move_to(cr,x,y)
    dx, dy = device_to_user_distance!(cr,[1., 0.])
    rectangle(cr,x,y,dx,dx)
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke(cr)
    end
end

function quad(x1, y1, x2, y2, x3, y3, x4, y4)
    move_to(cr, x1, y1)
    line_to(cr, x2, y2)
    line_to(cr, x3, y3)
    line_to(cr, x4, y4)
    close_path(cr)
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
    if state.fillStuff
        set_source(cr, state.fillCol)
        fill(cr)
    end
end

function rect(xtopleft, ytopleft, width, height)
    if state.rectMode == CORNER
        rectangle(cr, xtopleft, ytopleft, width, height)
    elseif state.rectMode == CORNERS # in this case, width and height are
                                     # reinterpreted as (x,y) coords of
                                     # bottom-right corner
        rectangle(cr, xtopleft, ytopleft, width-xtopleft, height-ytopleft)
    elseif state.rectMode == CENTER
        rectangle(cr, xtopleft-width/2, ytopleft-height/2, width, height)
    elseif state.rectMode == RADIUS
        rectangle(cr, xtopleft-width/2, ytopleft-height/2, width/2, height/2)
    end
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
    if state.fillStuff
        set_source(cr, state.fillCol)
        fill(cr)
    end
end

function triangle(x1,y1,x2,y2,x3,y3)
    move_to(cr, x1, y1)
    line_to(cr, x2, y2)
    line_to(cr, x3, y3)
    close_path(cr)
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
    if state.fillStuff
        set_source(cr, state.fillCol)
        fill(cr)
    end
end

## Curves

function bezier(x1, y1, x2, y2, x3, y3, x4, y4)
    move_to(cr, x1, y1);
    curve_to(cr, x2, y2, x3, y3, x4, y4);
    if state.strokeStuff
        set_source(cr, state.strokeCol)
        stroke_preserve(cr)
    end
end

#bezierDetail
#bezierPoint
#bezierTangent
#curve
#curveDetail
#curvePoint
#curveTangent
#curveTightness

## Attributes

function ellipseMode(eMode)
    state.ellipseMode = eMode
end

function noSmooth()

end

function rectMode(rMode)
    state.rectMode = rMode
end

function smooth()

end

function strokeCap(capType)
    if capType == ROUND

    elseif capType == SQUARE

    elseif capType == PROJECT

    end
end

function strokeJoin(joinType)
    if joinType == MITER

    elseif joinType == BEVEL

    elseif joinType == ROUND

    end
end

function strokeWeight(newWeight)

end

## Vertex

# function beginShape()
#
# end

#bezierVertex
#curveVertex

# function endShape()
#
# end

#quadraticVertex

# function vertex(x,y)
#
# end

## Loading & Displaying

#shape
#shapeMode

## Input

# Mouse

function mouseButton()
    if state.mouse1Pressed
        return LEFT
    elseif state.mouse2Pressed
        return CENTER
    elseif state.mouse3Pressed
        return RIGHT
    end
end

function mouseClicked()

end

function mouseDragged()

end

function mouseMoved()

end

function mousePressed()
    return pressed = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_LEFT) == GLFW.PRESS
end

function mouseReleased()

end

function mouseWheel()

end

function mouseX()
    (mx, my) = GLFW.GetCursorPos(window)
    return mx
end

function mouseY()
    (mx, my) = GLFW.GetCursorPos(window)
    return my
end

function pmouseX()
    return state.pmX
end

function pmouseY()
    return state.pmY
end

## Transform

#applyMatrix()

function popMatrix()

end

function printMatrix()

end

function pushMatrix()

end

function resetMatrix()

end

function rotate(ang)

end

function scale(sx, sy)

end

#shearX()
#shearY()

function translate(x, y)

end

# Color

## Setting

function background(r, g, b, a)
    state.bgCol = RGB(r, g, b)
    glClearColor(r, g, b, a)
    glClear(GL_COLOR_BUFFER_BIT)
end

function colorMode(mode::String)
    state.cMode = mode
end

function fill(r, g, b, a)
    if state.fillStuff == false
        state.fillStuff = true
    end
    if state.cMode == "RGB"
        state.fillCol = RGB(r, g, b)
    else
    end
end

function noFill()
    state.fillStuff = false
end

function noStroke()
    state.strokeStuff = false
end

function stroke(r, g, b, a)
    if state.strokeStuff == false
        state.strokeStuff = true
    end
    if state.cMode == "RGB"
        state.strokeCol = RGB(r, g, b)
    else
    end
end

## Creating & Reading

function alpha(c::ColorValue)
    return c.a
end

function blue(c::ColorValue)
    return c.b
end

function brightness(c::ColorValue)
    hsv = convert(HSV, c)
    return hsv.v
end

function color(r, g, b)
    return RGB(r, g, b)
end

function green(c::ColorValue)
    return c.g
end

function hue(c::ColorValue)
    hsv = convert(HSV, c)
    return hsv.h
end

function lerpColor(c1::ColorValue, c2::ColorValue, amt::Float32)
    return weighted_color_mean(amt, c1, c2)
end

function red(c::ColorValue)
    return c.r
end

function saturation(c::ColorValue)
    hsv = convert(HSV, c)
    return hsv.s
end

# Image

#function createImage(x, y, colorSpace::String)
#
#end

## Loading & Displaying

function image(img, x, y, w, h)

end

#imageMode

function loadImage(fileName::String)
    return
end

function noTint()
    state.tintStuff = false
end

#requestImage

#function tint()
#
#end

## Textures

#texture
#textureMode
#textureWrap

## Pixels

#blend
#copy
#filter
#get
#loadPixels
#set
#updatePixels

# Rendering

# function blendMode()
#
# end

#createGraphics

# Typography

## Loading & Displaying

function createFont(fontName::String, fontSize::Float32)

end

#loadFont

function text(str::String, x, y)

end

function textFont(fontName::String)

end

## Attributes

#function textAlign()
#
#end

#textLeading
#textMode

function textSize(fontSize)

end

function textWidth(str::String)

end

## Metrics

#function textAscent(str::String)
#    extents = scaled_font_extents(cr, str)
#    return extents[1]
#end

#function textDescent(str::String)
#    extents = scaled_font_extents(cr, str)
#    return extents[2]
#end

end # module Processing


text:
using Dates
using OpenSSL

function create_self_signed_certificate()
    # Create a root certificate.
    x509_certificate = X509Certificate()

    evp_pkey_ca = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey_ca

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:TRUE"))
    add_extension(x509_certificate, X509Extension("keyUsage", "keyCertSign"))

    sign_certificate(x509_certificate, evp_pkey_ca)

    root_certificate = x509_certificate

    # Create a certificate sign request.
    x509_request = X509Request()

    evp_pkey = EvpPKey(rsa_generate_key())

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_request.subject_name = x509_name

    x509_exts = StackOf{X509Extension}()

    ext = X509Extension("subjectAltName", "DNS:localhost")
    push!(x509_exts, ext)
    add_extensions(x509_request, x509_exts)
    finalize(ext)

    finalize(x509_exts)

    sign_request(x509_request, evp_pkey)

    # Create a certificate.
    x509_certificate = X509Certificate()
    x509_certificate.version = 2

    # Set issuer and subject name of the cert from the req and CA.
    x509_certificate.subject_name = x509_request.subject_name
    x509_certificate.issuer_name = root_certificate.subject_name

    x509_exts = x509_request.extensions

    ext = pop!(x509_exts)

    add_extension(x509_certificate, ext)
    add_extension(x509_certificate, X509Extension("keyUsage", "digitalSignature, nonRepudiation, keyEncipherment"))
    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:FALSE"))

    # Set public key
    x509_certificate.public_key = x509_request.public_key

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    sign_certificate(x509_certificate, evp_pkey_ca)

    # Store the certificate.
    private_key_io = IOBuffer()
    write(private_key_io, evp_pkey)

    public_key_io = IOBuffer()
    write(public_key_io, x509_certificate)

    private_key_pem = String(take!(private_key_io))
    public_key_pem = String(take!(public_key_io))

    return private_key_pem, public_key_pem
end

"""
    Creates a self signed certificate.
"""
function create_self_signed_root_certificate()
    x509_certificate = X509Certificate()

    evp_pkey = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:TRUE"))
    add_extension(x509_certificate, X509Extension("keyUsage", "keyCertSign"))

    sign_certificate(x509_certificate, evp_pkey)

    private_key_io = IOBuffer()
    write(private_key_io, evp_pkey)

    public_key_io = IOBuffer()
    write(public_key_io, x509_certificate)

    private_key_pem = String(take!(private_key_io))
    public_key_pem = String(take!(public_key_io))

    return private_key_pem, public_key_pem
end


text:
using GrammarExpts, ExprSearch, ACASXProblem, Gallium
using ExprSearch.GP

problem = ACASXClustering("dasc")
grammar = get_grammar(problem)


text:
mutable struct Tree <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    status::Bool  # true is green and false is burning
end

function model_initiation(; f = 0.02, d = 0.8, p = 0.01, griddims=(100,100), seed = 111)
    Random.seed!(seed)
    space = GridSpace(griddims, moore = true)
    properties = Dict(:f => f, :d => d, :p => p)
    forest = AgentBasedModel(Tree, space; properties = properties)

    ## create and add trees to each node with probability d,
    ## which determines the density of the forest
    for node in nodes(forest)
        if rand() ≤ forest.d
            add_agent!(node, forest, true)
        end
    end
    return forest
end

function forest_step!(forest)
  for node in nodes(forest, by = :random)
    nc = get_node_contents(node, forest)
    ## the cell is empty, maybe a tree grows here
    if length(nc) == 0
        rand() ≤ forest.p && add_agent!(node, forest, true)
    else
      tree = forest[nc[1]] # by definition only 1 agent per node
      if tree.status == false  # if it is has been burning, remove it.
        kill_agent!(tree, forest)
      else
        if rand() ≤ forest.f  # the tree ignites spontaneously
          tree.status = false
        else  # if any neighbor is on fire, set this tree on fire too
          for cell in node_neighbors(node, forest)
            neighbors = get_node_contents(cell, forest)
            length(neighbors) == 0 && continue
            if any(n -> !forest.agents[n].status, neighbors)
              tree.status = false
              break
            end
          end
        end
      end
    end
  end
end


text:
function correctorfluxes!(Fbarforward::Vector{Float64},Fbarbackward::Vector{Float64},
    Qforward::Vector{Float64},Aforward::Vector{Float64},Qbackward::Vector{Float64},
    Abackward::Vector{Float64},beta::Float64,rho::Float64,
    acolscor::Vector{Int8},qcolscor::Vector{Int8})
    Fbarforward[acolscor] .= Qforward;
    Fbarforward[qcolscor] .= Qforward.^2./Aforward .+ (1./3).*beta./rho.*
        Aforward.^1.5;
        # system.branches.Fbarforward[i][system.solverparams.acolscor] .= (
        #     system.branches.Qforward[i]);
        # system.branches.Fbarforward[i][system.solverparams.qcolscor] .= (
        #     system.branches.Qforward[i].^2./system.branches.Aforward[i] .+
        #     (1./3).*system.branches.beta[i][end]./system.solverparams.rho.*
        #     system.branches.Aforward[i].^1.5);
    Fbarbackward[acolscor] .= Qbackward;
    Fbarbackward[qcolscor] .= Qbackward.^2./Abackward .+ (1./3).*beta./rho.*
        Abackward.^1.5;
        # system.branches.Fbarbackward[i][system.solverparams.acolscor] .= (
        #     system.branches.Qbackward[i]);
        # system.branches.Fbarbackward[i][system.solverparams.qcolscor] .= (
        #     system.branches.Qbackward[i].^2./system.branches.Abackward[i] .+
        #     (1./3).*system.branches.beta[i][end]./system.solverparams.rho.*
        #     system.branches.Abackward[i].^1.5);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
grid_size = 25
c_vals = linspace(2, 12, grid_size)
w_bar_vals = similar(c_vals)

mcm = McCallModel()

for (i, c) in enumerate(c_vals)
    mcm.c = c
    w_bar = compute_reservation_wage(mcm)
    w_bar_vals[i] = w_bar
end

plot(c_vals,
     w_bar_vals,
     lw=2,
     α=0.7,
     xlabel="unemployment compensation",
     ylabel="reservation wage",
     label=L"$\bar w$ as a function of $c$")

max_stars_repo_name: yywingliang/lecture-source-jl

max_issues_repo_name: yywingliang/lecture-source-jl

max_forks_repo_name: yywingliang/lecture-source-jl


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "Nettle"
version = v"3.4.1"
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/nettle/nettle-$(version).tar.gz",
                  "f941cf1535cd5d1819be5ccae5babef01f6db611f9b5a777bae9c7604b8a92ad"),
]

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

# Force c99 mode
export CFLAGS="${CFLAGS} -std=c99"

update_configure_scripts
./configure --prefix=$prefix --host=$target --with-include-path=$prefix/include
make -j${nproc}
make install
"""

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

# Disable FreeBSD for now, because hogweed needs alloca()?
platforms = [p for p in platforms if !(typeof(p) <: FreeBSD)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libnettle", :libnettle),
    LibraryProduct("libhogweed", :libhogweed),
    ExecutableProduct("nettle-hash", :nettle_hash)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "GMP_jll",
]

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


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = y = -10:0.11:10
    y1d =  sin.(x) ./ x
    # 3D heatmap
    sinc2d(x,y) = sin.(sqrt.(x.^2 + y.^2))./sqrt.(x.^2+y.^2)
    fxy = [sinc2d(x,y) for x in x, y in y]

    fig = Figure(resolution = (700, 450))
    ax1 = Axis(fig, xlabel = "x", ylabel = "f(x)", xgridvisible = true,
            ygridvisible = true)
    ax2 = Axis(fig, bbox = BBox(140, 260, 260, 350), xticklabelsize = 12,
          yticklabelsize = 12, showgrid=false,
          title = "inset  at (140, 260, 260, 350)")
    line1 = lines!(ax1, x, y1d, color = :red)
    # inset
    hmap = heatmap!(ax2, x, y,fxy, colormap = :Spectral_11)
    limits!(ax2, -10,10,-10,10)
    hidespines!(ax2)
    ax2.yticks = [-10,0,10]
    ax2.xticks = [-10,0,10]

    leg = Legend(fig, [line1], ["sinc(x)"], halign = :right, valign = :top,
        tellheight = false, tellwidth = false,)
    cbar = Colorbar(fig, hmap, label ="sinc(x,y)", labelpadding = 5,
        tellheight = false, tellwidth = false, ticklabelsize = 12,
        width = 10, height = Relative(1.5/4), halign = :right, valign = :center)

    fig[1, 1] = ax1
    fig[1, 1] = leg
    fig[1, 1] = cbar
    fig
    #save("LineWithInsetHeatmap.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "LineWithInsetHeatmap.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
" The following initializes the master problem and runs a loop iterating between adding scenarios
and computing a solution for the current subset of scenarios. It is described in Algorithms 1 and 3 of the paper.
  * It is assumed throughout that the number of uncertain constraints is equal to data.ncons.
  * It is assumed throught that f(ξ,x) = ξᵀCx h(x).
  * Uncertain constraints are assumed to be linear functions g_ℓ(ξ,x) = ξᵀ(A_ℓ)x."

function scenario_generation(heuristic,warmup)
  # cc and aa record the uncertainty vectors ξ generated so far and multiplies them by the matrices C and A_ℓ.
  # Specifically:
  # cc[j,:] = Cᵀξ^j for each j in [m]
  # aa[ℓ][j,:] = Aᵀ_ℓ ξ^j for eacsh j in [m] and ℓ in [L]
  aa, cc = populate_scenarios()
  tΞ = [] #set of scenarios generated so far
  incumbent = []
  if !heuristic
    # If the algorithm is called in exact mode then a heuristic will be run first to obtain an upper bound
    # This upper bound serves to filter the set of solutions to be generated X(UB)
    if heuristicmode == "DUAL"
      @timeit to "heuristic" heurvalue,incumbent = heuristic_dualization()
    elseif heuristicmode == "HEURCG"
      @timeit to "heuristic" heurvalue,incumbent = scenario_generation(true,warmup)
    end
    # Current upper bound is set to the best solution value returned by the heuristic
    ub = heurvalue
    data.heuristic_ub = ub
    !warmup && @info "exact mode"
    # All the solutions in the set X(UB) are generated by using a problem-specific enumeration algorithm
    @timeit to "gen_sol" solutions = generate_solutions(ub)
  else
    # If the algoritm is called in heuristic mode then it is initalized simply with a set of solutions
    # Below this is done simply by starting with the optimal solution using the data of the first scenario and its neighbors
    !warmup && @info "heuristic mode"

    if data.constraint_uncertainty
      a_vec = aa[:][1,:]
    else
      a_vec = []
    end
    # A deterministic optimal solution is obtained using the data realization of first scenario
    detsol = find_new_solution(cc[1,:],a_vec,Inf,[],[],false,false,TIME_LIM)
    # The incumbent is initialized by creating k copies of this deterministic solution
    incumbent = Array{Int64}(undef, k, n)
    for k in K incumbent[k,:] = round_solution(detsol) end
    # The set of solutions is enriched by adding some neighbors of the incumbent
    @timeit to "gen_sol" solutions = generate_neighbors(incumbent,aa)
    # An initial (poor) upper bound is obtained by computing the robust cost of detsol
    ub = compute_basic_ub(detsol)
  end
  # In some pathological cases where the number of enumerated solutions is < k we add copies of the first solution
  for kk in (size(solutions,1)+1):k push!(solutions,solutions[1]) end
  # S is the number of solutions generated
  S = 1:size(solutions,1)
  # M is the number of scenarios added so far
  M = 1:size(cc,1)

  # Now we proceed to the iterations of the algorithm
  # solutions containts all solutions from X(UB) in the exact mode and a subset of these solutions in the heuristic mode

  # indicator new_scenario_added is used as a stopping criteria of convergence
  new_scenario_added = true

  # First initialize the structures that will be used throughout the algorithm
  # x_sol_array contains all solution vectors
  # sol_items contains the indices of non-zero variables in a solution vector
  x_sol_array = Vector{Vector{Float64}}()
  sol_items = Vector{Vector{Int}}()
  # Populate x_sol_array and sol_items using the enumerated solutions
  fill_solutions_bit(solutions, x_sol_array, sol_items)
  # The following vector contains h(x) for each x in X
  hh = [h(sol_items[s]) for s in S]
  # Starting lb: min_x∈X ξ^Tx for some ξ ∈ Ξ
  # In the following the initial scenario is used as ξ
  lb = minimum(sum(x_sol_array[s][l]*cc[1,l]*hh[s] for l in N) for s in S)

  # d records cost of each solution in each scenario, specifically
  # d[s,j] = ξ_jᵀCx_s h(x_s) for s in [r] and j in [m]
  d = Matrix{Float64}(undef,size(solutions,1),size(cc,1))
  # a_cov records constraint satisfaction of each solution in each scenario, specifically
  # (a_cov[ℓ][s,j] = 1) ⇔ (ξ_jᵀ(A_ℓ)x_s h(x_s) ≤ b_ℓ), for s in [r] and j in [m]
  a_cov = Vector{Matrix{Bool}}()
  if data.constraint_uncertainty
    for cons in 1:data.ncons
      push!(a_cov,Matrix{Bool}(undef,size(solutions,1),size(cc,1)))
    end
  end
  # initialize the number of iterations
  NUMITER = 0

  # calculate the initial gap
  gap = 100*round((ub - lb)/abs(ub), digits = 4)

  if time()-STARTTIME < TIME_LIM
    println("IT LB  UB  gap time")
    while new_scenario_added
      ξ = []
      NUMITER += 1
      if heuristic
        # In the heuristic mode d and a_cov are recalculated from scratch in each iteration
        # as this is not very costly with a small set of solutions
        @timeit to "re-calculate d" d = calculate_d(x_sol_array,cc,S,M,hh)
        if data.constraint_uncertainty
          @timeit to "re-calculate a_cov" a_cov = calculate_a(x_sol_array,aa,S,M)
        end
      else
        # In the exact mode d and a_cov are updated each time a new scenario is introduced
        # This permits reusing all d and a_cov values for the previously added scenarios
        @timeit to "update d" d = update_d(x_sol_array,cc,d,S,M,NUMITER,hh)
        if data.constraint_uncertainty
          @timeit to "update a_cov" a_cov = update_a(x_sol_array,aa,a_cov,S,M,NUMITER)
        end
      end
      # The master problem is solved using a binary search algorithm
      # It returns the current optimal solution and its value
      if !(data.constraint_uncertainty)
        @timeit to "binary search" x̃, ṽ_down, ṽ_up = binary_search(x_sol_array,sol_items,lb,ub,min(gap,1),d,a_cov,cc,δ)
      else
        @timeit to "binary search" x̃, ṽ_down, ṽ_up = binary_search(x_sol_array,sol_items,lb,ub,ϵ,d,a_cov,cc,δ)
      end

      # The optimal value of the master problem provides a lower bound
      lb = ṽ_down
      # If the master problem is solved succesfully then we proceed to the separation
      # If not then the binary search algorithm was stopped due to time limit restrictions
      # In this case it returns an empty solution permitting the algorithm to skip the separation step and break out of the while loop
      if length(x̃) > 0
        # Separation problem returns an identified scenario along with the separation objective value
        @timeit to "separation" ξ, sep_value = build_and_solve_separation(x̃, ṽ_up)
        # Because of the tolerance garantee when building the covcering matrix, it may happen that a given scenario be generated more than once
        # If that happens, we consider instead that the solution did not return any scenario.
        if data.constraint_uncertainty
          @timeit to "scenario generated?" begin
            if length(ξ) > 0 && length(tΞ) > 0
              for ξ′ in tΞ
                if norm(ξ-ξ′) < ϵ*min(norm(ξ),norm(ξ′))
                  @warn "scenario already generated"
                  ξ = []
                  sep_value = 0
                  break
                end
              end
            end
          end
        end
        length(ξ) > 0 && push!(tΞ,ξ)

        if !data.constraint_uncertainty
          # If there is no constraint uncertainty then upper bound updates are done whenever the separation value
          # improves upon the current bound. The upper bound is obtained simply as the separation value.
          # As the adversarial function is continuous this is a valid upper bound.
          if sep_value < ub  && sep_value > -Inf
            ub = sep_value
            incumbent = x̃
          end
        end
      end
      # Update the optimality gap
      gap = 100*round((ub - lb)/abs(ub), digits = 4)

      # If a scenario is returned from the separation problem and the convergence is not yet achieved
      # then add the realization of this scenario to vectors cc and aa.
      # Otherwise the indicator new_scenario_added is turned false to break out of the while loop.
      # The convergence is tested using the relative gap between ub and lb in the case where there is only objective uncertainty
      # and using the value returned by the function S(x̃,ṽ) in the case of contraint uncertainty
      if !(data.constraint_uncertainty)
        if length(ξ) == 0 || gap ≤ ϵ
          new_scenario_added = false
        end
      else
        if length(ξ) == 0 || sep_value ≤ ϵ
          new_scenario_added = false
        end
      end

      if new_scenario_added
        # We compute Cᵀξ and Aᵀ_ℓ ξ for the scenario ξ just obtained
        # println("Scenario added: $ξ")
        new_c, new_a = add_scenario(ξ)
        # M is increased by 1
        M = 1:(length(M)+1)
        # cc and aa are augmented with the new realization
        cc = [cc; new_c']
        if data.constraint_uncertainty
          for cons in 1:data.ncons
            aa[cons] = [aa[cons]; new_a[cons]']
          end
        end
      else
        # The current solution leads to S(x̃,ṽ) so the ub on its optimal value is an ub for the true problem
        ub = ṽ_up
        gap = 100*round((ub - lb)/abs(ub), digits = 4)
      end

      if heuristic && new_scenario_added
        # HEURISTIC MODE: the new scenario offers an opportunity to enrich the set of solutions
        if data.constraint_uncertainty
          a_vec = aa[:][end,:]
        else
          a_vec = []
        end
        # A good solution for the newly added scenario and its neighbors are added to the set of enumerated solutions
        @timeit to "relax and fix" new_solution = relax_and_fix(cc[end,:],a_vec,ub)

        if length(new_solution)>0
          @timeit to "generate neighbors" solutions_to_add = generate_neighbors(reshape(round_solution(new_solution),length(new_solution),1)',aa)
          for sol in solutions_to_add
            if !(sol in solutions)
              push!(solutions,sol)
            end
          end
        end
        # The following updates data-structures used throughout with the newly discovered solutions
        x_sol_array = Vector{Vector{Float64}}()
        sol_items = Vector{Vector{Int}}()
        fill_solutions_bit(solutions, x_sol_array, sol_items)
        S = 1:size(solutions,1)
        hh = [ h(sol_items[s]) for s in S ]
      end
      TotalTime = round(time()-STARTTIME, digits = 2)
      LowerBound = round(lb, digits = 4)
      println("$NUMITER $LowerBound $ub $gap $TotalTime")
      TotalTime > TIME_LIM && break
    end
  end
  if heuristic
    return ub,incumbent
  else
    return lb, ub, incumbent, gap, length(M), NUMITER, length(S)
  end
end

#-----------------------------------------------------------------------------------

"The following computes an optimal solution for the current set of scenarios using
a binary search algorithm based on the optimal solution cost, in line with the radius
formulation used for the p-center problem. It is described in Algorithm 2 of the paper."

function binary_search(x_sol_array,sol_items,lb,ub,ϵ_BS,d,a_cov,cc,δ)

  " For a given coverage matrix identify the non-dominated solutions."
  function calculate_nondom_solutions(cov,last,M,S)
    map = [BitSet([j for j in M if cov[s, j]]) for s in S]
    discarded = falses(length(S))
    for s in 1:(length(S)+1-k)
      if length(S)-length(discarded)==k break end
      if discarded[s] continue end
      for ss in 1:(length(S)+1-k)
        if s == ss continue
        elseif discarded[ss] continue
        elseif issubset(map[ss], map[s])
            discarded[ss] = true
            break
        elseif issubset(map[s], map[ss])
            discarded[s] = true
            break
        end
      end
    end
    SS = [s for s in S if !discarded[s]]
    return SS
  end

  " For a given coverage matrix identify the non-dominated scenarios."
  function calculate_nondom_scenarios(cov,M,S)
    map = [BitSet([s for s in S if cov[s, j]]) for j in M]
    discarded = falses(length(M))
    for m in 1:(length(M)+1-k)
        if discarded[m] continue end
        for mm in 1:(length(M)+1-k)
          if m == mm continue
          elseif discarded[mm] continue
          elseif issubset(map[mm],map[m])
            discarded[m] = true
            break
          elseif issubset(map[m],map[mm])
            discarded[mm] = true
            break
          end
        end
    end
    MM = [m for m in M if !discarded[m]]
    return MM
  end

  " For a given radius, test whether the associated covering problem is feasible using
  a MILP formulation with slack variables."
  function build_and_solve_covering_MIP(SS,MM,cov)
    time_remaining = TIME_LIM - (time()-STARTTIME)
    cover = create_model(time_remaining)
    @variable(cover, y[SS], Bin)
    @variable(cover, slack[MM] ≥ 0)
    @constraint(cover, covering[m in MM], sum(y[s] for s in SS if cov[s,m]) >= 1-slack[m]) # covering constraints
    @constraint(cover, lessthank, sum(y[s] for s in SS) == k) # cardinality constraint
    @objective(cover, Min, sum(slack[m] for m in MM))
    optimize!(cover)
    sol_used = []
    if termination_status(cover) == MOI.OPTIMAL
      for s in SS
        if value(y[s]) > 0.9
          push!(sol_used,s)
        end
      end
      return objective_value(cover),sol_used
    else
      return Inf,sol_used
    end
  end

  " Translate the binary search solution to a solution vector x."
  function compute_x_values(S, sol_items, solutions_used)::Array{Int64,2}
      x = zeros(k, data.n)
      k_ = 0
      for s in solutions_used
          k_ += 1
          for l in sol_items[S[s]]
              x[k_ ,l] = 1
          end
      end
      return x
  end

  S = 1:length(x_sol_array)
  M = 1:size(d,2)
  # Indicator for verifying if this is the last iteration of the algorithm
  # When last=true a feasible solution corresponding to the ub is calculated
  last = false
  # Count the number of iterations of the binary search algorithm
  ITER_BS = 0
  # Initialize the optimal solution vector
  x = []
  # The while loop is executed at least for the first and the last iteration
  # and until the convergence criteria on the relative gap is met
  while (ub - lb)/abs(ub) > ϵ_BS/100 || last || ITER_BS == 0
    ITER_BS += 1
    time()-STARTTIME > TIME_LIM && break

    if last
      # In the last iteraion r is set large enough to find a feasible solution
      r = ub
    else
      # Otherwise binary search classical update of r
      r = (lb+ub)/2
    end
    # A coverage matrix is built based on whether a solution covers a scenario with a cost <= r
    # For Smuggler, increase the precision to ϵ
    @timeit to "build cov" cov = [d[s,m] ≤ r+ϵ*abs(r) for s in S, m in M]
    if data.constraint_uncertainty
      for cons in 1:data.ncons
        # The objective coverage matrix is multiplied with the constraint coverage matrix to filter out infeasible solutions
        # (this might be sped-up by avoiding repeating 0 elements and using sparse representations)
        @timeit to "update cov" cov = cov .& a_cov[cons]
      end
    end
    infeasible = false
    # Non-dominated solutions and scenarios are filtered based on the coverage matrix
    @timeit to "dominance solutions" SS = calculate_nondom_solutions(cov,last,M,S)
    @timeit to "dominance scenarios" MM = calculate_nondom_scenarios(cov,M,S)
    # A MIP determining whether it is possible to cover all scenarios by choosing K solutions is solved
    # (notice only non-dominated scenarios and solutions are sent to this MIP)
    @timeit to "solve covering problem" objval,solutions_selected = build_and_solve_covering_MIP(SS,MM,cov)
    # If its optimal value is >= 1 then the MIP is infeasible
    infeasible = objval > 0.9
    if !infeasible
      # If the covering MIP is feasible then a feasible solution with value r can be found and therefore
      # the upper bound is updated
      ub = r
      if last
        # If this is the last iteration of the algorithm a feasible solution is returned
        x = compute_x_values(S, sol_items, solutions_selected)
        break
      end
    else
      # If the covering MIP is infeasible then the lower bound can be increased to r
      if last && objval < Inf
        @error "Shouldn't be here"
        exit()
      elseif last && objval >= Inf
        return [], lb, ub
      else
        lb = r
      end
    end
    # When convergence is achieved last iteration indicator is turned true and one last iteration of BS is done
    if (ub - lb)/abs(ub) ≤ ϵ_BS
      last = true
    end
  end
  # The algorithm returns with an optimal solution and a lower bound on the optimal master value
  # Note that if the algorithm is broken due to time limit restrictions then we return x=[]
  return x, lb, ub
end

#-----------------------------------------------------------------------------------

function update_d(x_sol_array,cc,d,S,M,NUMITER,hh)
  if NUMITER == 1
    for s in S, m in M
      d[s,m] = hh[s]*sum(x_sol_array[s][l]*cc[m,l] for l in N)
    end
  else
    new_costs = [hh[s]*sum(x_sol_array[s][l]*cc[length(M),l] for l in N) for s in S]
    d = [d new_costs]
  end
  return d
end

#-----------------------------------------------------------------------------------

function update_a(x_sol_array,aa,a_cov,S,M,NUMITER)
  if NUMITER == 1
    for cons in data.ncons, s in S, m in M
      a_cov[cons][s,m] = sum(x_sol_array[s][l]*aa[cons][m,l] for l in N) ≤ data.B[cons]
    end
  else
    for cons in data.ncons
      new_cov = [sum(x_sol_array[s][l]*aa[cons][length(M),l] for l in N) ≤ data.B[cons] for s in S]
      a_cov[cons] = [a_cov[cons] new_cov]
    end
  end
  return a_cov
end

#-----------------------------------------------------------------------------------

function calculate_d(x_sol_array,cc,S,M,hh)
  d = Matrix{Float64}(undef,length(S),length(M))
  for s in S, m in M
    d[s,m] = hh[s]*sum(x_sol_array[s][l]*cc[m,l] for l in N)
  end
  return d
end

#-----------------------------------------------------------------------------------

function calculate_a(x_sol_array,aa,S,M)
  a_cov = Vector{Matrix{Bool}}()
  for cons in 1:data.ncons
    push!(a_cov,Matrix{Bool}(undef,length(S),length(M)))
    for s in S, m in M
      a_cov[cons][s,m] = sum(x_sol_array[s][l]*aa[cons][m,l] for l in N) ≤ data.B[cons]
    end
  end
  return a_cov
end

#-----------------------------------------------------------------------------------

function find_new_solution(c_vec,a_vec,UB,fixedzero,fixedone,relax,bound,TIME)
  sol_model = create_model(TIME)
  if relax
    @variable(sol_model, 0<=x[i in N]<=1, container=Array)
  else
    @variable(sol_model, x[i in N], Bin, container=Array)
  end
  add_constraints_X(sol_model,x)
  if data.constraint_uncertainty
    @constraint(sol_model, [cons in 1:data.ncons], sum(a_vec[cons][i]*x[i] for i in N) ≤ data.B[cons] )
  end
  @constraint(sol_model,[i in fixedzero],x[i]==0)
  @constraint(sol_model,[i in fixedone],x[i]==1)
  if bound @constraint(sol_model,sum(x[i]*c_vec[i] for i in N) ≤ UB) end
  @objective(sol_model,Min,sum(x[i]*c_vec[i] for i in N))
  optimize!(sol_model)
  if termination_status(sol_model) == MOI.OPTIMAL
    return value.(x)
  else
    return []
  end
end

#-----------------------------------------------------------------------------------

function relax_and_fix(cc,aa,UB)
  fixedzero = []
  fixedone = []
  time_remaining = TIME_LIM-(time()-STARTTIME)
  @timeit to "relaxation" solution = find_new_solution(cc,aa,UB,fixedzero,fixedone,true,true,time_remaining)
  if size(solution,1) > 0
    for i in N
      if solution[i]>=1-ϵ
        push!(fixedone,i)
      elseif solution[i]<=ϵ
        push!(fixedzero,i)
      end
    end
    @timeit to "integer" solution = find_new_solution(cc,aa,UB,fixedzero,fixedone,false,false,time_remaining)
    if size(solution,1) > 0
      return solution
    else
      return []
    end
  else
    return []
  end
end

#-----------------------------------------------------------------------------------|
" Traslate the enumerated solutions (vector of bool) to solution vectors in 0/1 (x_sol_array). For each solution we also
store the vector of indices of the non-zero variables (sol_items)."
function fill_solutions_bit(solutions, x_sol_array, sol_items)
    S = 1:size(solutions,1)
    for s in S
        sol = zeros(data.n)
        for i in N
            if solutions[s][i]
              sol[i] = 1
            end
        end
        push!(x_sol_array, sol)
        items = Vector{Int}()
        for l in N
            if sol[l] == 1
                push!(items, l)
            end
        end
        push!(sol_items, items)
    end
end

#-----------------------------------------------------------------------------------

" Round the variable values that pass the solver tolerance checks."
function round_solution(sol)
  for i=1:length(sol)
    if sol[i] <= ϵ
      global sol[i] = 0
    elseif sol[i] >= 1-ϵ
      global sol[i]= 1
    end
  end
  return sol
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
"Fits a straight line through a set of points, `y = a₁ + a₂ * x`"
linear_fit(x, y) = hcat(fill!(similar(x),1), x) \ y

"Fits a log function through a set of points: `y = a₁+ a₂*log(x)`"
log_fit(x, y) = linear_fit(log.(x), y)

"Fits a power law through a set of points: `y = a₁*x^a₂`"
function power_fit(x, y)
    fit = linear_fit(log.(x), log.(y))
    [exp(fit[1]), fit[2]]
end

"Fits an `exp` through a set of points: `y = a₁*exp(a₂*x)`"
function exp_fit(x, y)
    fit = linear_fit(x, log.(y))
    [exp(fit[1]), fit[2]]
end

"""
Fits a polynomial of degree `n` through a set of points.

Simple algorithm, doesn't use orthogonal polynomials or any such thing
and therefore unconditioned matrices are possible. Use it only for low
degree polynomials.

This function returns a the coefficients of the polynomial.
"""
function poly_fit(x, y, n)

    nx = length(x)
    A = zeros(eltype(x), nx, n+1)
    A[:,1] .= 1.0
    for i in 1:n
        for k in 1:nx
            A[k,i+1] = A[k,i] * x[k]
        end
    end
    A\y
end

"""
High Level interface for fitting straight lines
"""
struct LinearFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    LinearFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
LinearFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = LinearFit{T}(linear_fit(x, y))

"High Level interface for fitting log laws"
struct LogFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    LogFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
LogFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = LogFit{T}(log_fit(x, y))

"High Level interface for fitting power laws"
struct PowerFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    PowerFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
PowerFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = PowerFit{T}(power_fit(x, y))

"High Level interface for fitting exp laws"
struct ExpFit{T<:Number} <: LeastSquares
    coefs::Array{T,1}
    ExpFit{T}(coefs) where {T<:Number} = new(copy(coefs))
end
ExpFit(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Number} = ExpFit{T}(exp_fit(x, y))

"""
# Generic interface for curve fitting.

The same function `curve_fit` can be used to fit the data depending on fit type,
shich is specified in the first parameter. This function returns an object that
can be used to estimate the value of the fitting model using function `apply_fit`.

## A few examples:

 * `f = curve_fit(LinearFit, x, y)`
 * `f = curve_fit(Poly, x, y, n)`
"""
curve_fit(::Type{T}, x, y) where {T<:LeastSquares} = T(x, y)
curve_fit(::Type{T}, x, y, args...) where {T<:LeastSquares} = T(x, y, args...)
curve_fit(::Type{Poly}, x, y, n=1) = Poly(poly_fit(x, y, n))

(f::LinearFit)(x) = f.coefs[1] + f.coefs[2] *x
(f::PowerFit)(x) = f.coefs[1] * x ^ f.coefs[2]
(f::LogFit)(x) = f.coefs[1] + f.coefs[2] * log(x)
(f::ExpFit)(x) = f.coefs[1] * exp(f.coefs[2] * x)


text:
"""
    enforce_no_flow_through!(p::Plate, motion, elements, t)

Update the plate, `p`, to enforce the no-flow-through condition given ambient vortex elements, `elements`, and while moving with kinematics specified by `motion`.

# Example

```jldoctest
julia> plate = Plate(128, 2.0, 0.0, π/3)
Plate: N = 128, L = 2.0, c = 0.0 + 0.0im, α = 60.0ᵒ
       LESP = 0.0, TESP = 0.0

julia> motion = Plates.RigidBodyMotion(1.0, 0.0);

julia> point = Vortex.Point(0.0 + 2im, 1.0);

julia> Plates.enforce_no_flow_through!(plate, motion, point, 0.0)

julia> plate
Plate: N = 128, L = 2.0, c = 0.0 + 0.0im, α = 60.0ᵒ
       LESP = 1.27, TESP = -1.93
```
"""
function enforce_no_flow_through!(p::Plate, ṗ, elements, t)
    @get p (L, C, α, dchebt!)
    @get ṗ (ċ, α̇)

    fill!(C, zero(ComplexF64))
    induce_velocity!(C, p, elements, t)

    n̂ = exp(-im*α)
    rmul!(C, n̂)

    dchebt! * C

    p.Γ = -circulation(elements)
    p.B₀ = normal(ċ, α)
    p.B₁ = 0.5α̇*L

    nothing
end

"""
    influence_on_plate!(∂A::Vector{ComplexF64}, plate::Plate, v, t)

Compute in-place the change in plate's Chebyshev coefficients `∂A` by a vortex element `v`
"""
function influence_on_plate!(∂A::Vector{ComplexF64}, plate::Plate, v, t)
    fill!(∂A, zero(ComplexF64))
    induce_velocity!(∂A, plate, v, t)
    rmul!(∂A, exp(-im*plate.α))
    plate.dchebt! * ∂A
    nothing
end

function influence_on_plate(plate::Plate, v, t)
    ∂A = Vector{ComplexF64}(undef, plate.N)
    influence_on_plate!(∂A, plate, v, t)
    return ∂A
end

"""
    vorticity_flux(p::Plate, v₁, v₂,
                   lesp = 0.0, tesp = 0.0,
                   ∂C₁ = Vector{ComplexF64}(undef, plate.N),
                   ∂C₂ = Vector{ComplexF64}(undef, plate.N))

Return strengths of new vortex elements that satisfies edge suction parameters.
For a given edge, if the current suction parameter is less than the criticial suction parameter, then no vorticity is released.  If it is higher, however, vorticity will be released so that the suction parameter equals the critical value.

# Arguments

- `p`: the plate
- `v₁, v₂`: the vortex elements (with unit circulation) that the vorticity flux is going into
- `lesp`, `tesp`: the critical leading and trailing edge suction parameters we want to enforce.  By default, both parameters are set to 0.0 to enforce the Kutta condition on both edges.  We can disable vortex shedding from an edge by setting the its critical suction parameter to `Inf`

# Returns

- `Γ₁, Γ₂`: the strengths that the vortex element should have in order to satisfy the edge suction parameters
- `∂C₁, ∂C₂`: Chebyshev coefficients of the normal velocity induced by the vortex elements
  Instead of running `enforce_bc!` with the new vortex elements, we can use this matrix to directly update the Chebyshev coefficients associated with the bound vortex sheet without recomputing all the velocities.

# Example

Enforcing the trailing edge Kutta condition with an point vortex at negative infinity:

```jldoctest
julia> plate = Plate(128, 2.0, 0.0, π/6)
Plate: N = 128, L = 2.0, c = 0.0 + 0.0im, α = 30.0ᵒ
       LESP = 0.0, TESP = 0.0

julia> motion = Plates.RigidBodyMotion(1.0, 0.0);

julia> Plates.enforce_no_flow_through!(plate, motion, (), 0.0)

julia> point = Vortex.Point(-Inf, 1.0);

julia> _, Γ, _, _ = Plates.vorticity_flux(plate, (), point, 0.0, Inf);

julia> Γ # should equal -πULsin(α) = -π
-3.1415926535897927
```
"""
function vorticity_flux(plate::Plate, v₁, v₂, t, lesp = 0.0, tesp = 0.0,
                        ∂C₁ = Vector{ComplexF64}(undef,plate.N),
                        ∂C₂ = Vector{ComplexF64}(undef,plate.N))

    @get plate (N, α, B₀, B₁, L, A, Γ)

    b₊ = +2(A[0] - B₀) + (A[1] - B₁) + 2Γ/(π*L)
    b₋ = -2(A[0] - B₀) + (A[1] - B₁) + 2Γ/(π*L)

    influence_on_plate!(∂C₁, plate, v₁, t)
    influence_on_plate!(∂C₂, plate, v₂, t)

    Γ₁ = circulation(v₁)
    Γ₂ = circulation(v₂)

    A₁₊ =  2imag(∂C₁[1]) + imag(∂C₁[2]) - 2Γ₁/(π*L)
    A₂₊ =  2imag(∂C₂[1]) + imag(∂C₂[2]) - 2Γ₂/(π*L)
    A₁₋ = -2imag(∂C₁[1]) + imag(∂C₁[2]) - 2Γ₁/(π*L)
    A₂₋ = -2imag(∂C₂[1]) + imag(∂C₂[2]) - 2Γ₂/(π*L)

    if (abs2(lesp) > abs2(b₊)) && (abs2(tesp) ≤ abs2(b₋))
        K₁, K₂ = 0.0, (sign(b₋)*tesp - b₋)/A₂₋
    elseif (abs2(lesp) ≤ abs2(b₊)) && (abs2(tesp) > abs2(b₋))
        K₁, K₂ = (sign(b₊)*lesp - b₊)/A₁₊, 0.0
    elseif (abs2(lesp) > abs2(b₊)) && (abs2(tesp) > abs2(b₋))
        K₁ = K₂ = 0.0
    else
        b₊ = sign(b₊)*lesp - b₊
        b₋ = sign(b₋)*tesp - b₋

        detA = A₁₊*A₂₋ - A₂₊*A₁₋

        @assert (detA != 0) "Cannot enforce suction parameters"

        K₁ = (A₂₋*b₊ - A₂₊*b₋)/detA
        K₂ = (A₁₊*b₋ - A₁₋*b₊)/detA
    end

    return K₁*Γ₁, K₂*Γ₂, rmul!(∂C₁, K₁), rmul!(∂C₂, K₂)
end

"""
    vorticity_flux!(p::Plate, v₁, v₂,
                    lesp = 0.0, tesp = 0.0,
                    ∂C₁ = Vector{ComplexF64}(undef,plate.N),
                    ∂C₂ = Vector{ComplexF64}(undef,plate.N))

In-place version of [`vorticity_flux`](@ref), except instead of just
returning the possible changes in plate Chebyshev coefficients, we
modify `plate.C` with those changes so that no-flow-through is
enforced in the presence of `v₁` and `v₂` with strengths that satisfy
the suction parameters.
"""
function vorticity_flux!(plate::Plate, v₁, v₂, t, lesp = 0.0, tesp = 0.0,
                        ∂C₁ = Vector{ComplexF64}(undef, plate.N),
                        ∂C₂ = Vector{ComplexF64}(undef, plate.N))
    Γ₁, Γ₂, _, _ = vorticity_flux(plate, v₁, v₂, t, lesp, tesp, ∂C₁, ∂C₂)
    @. plate.C += ∂C₁ + ∂C₂
    plate.Γ -= Γ₁ + Γ₂

    return Γ₁, Γ₂, ∂C₁, ∂C₂
end

function suction_parameters(plate)
    @get plate (N, α, B₀, B₁, L, A, Γ)

    b₊ = +2(A[0] - B₀) + (A[1] - B₁) + 2Γ/(π*L)
    b₋ = -2(A[0] - B₀) + (A[1] - B₁) + 2Γ/(π*L)
    return b₊, b₋
end


text:
@reaction_func new_hill(x, v, k, n) = v*x^n/(k^n+x^n)
@reaction_func new_poly(x) = 3x^2+1
@reaction_func new_exp(x) = exp(x)

network1 = @reaction_network rn begin
    hill(X,4,3,2), X + Y --> Z1
    mm(X,4.0,3), X + Y --> Z2
    3X^2+1, X + Y --> Z3
    exp(Y), X + Y --> Z4
    hillR(X,2.0,3.0,2), X + Y --> Z5
    mmR(X,4.0,1), X + Y --> Z6
end
network2 = @reaction_network rn begin
    new_hill(X,4,3,2), X + Y --> Z1
    4.0*X/(X+3), X + Y --> Z2
    new_poly(X), X + Y --> Z3
    new_exp(Y), X + Y --> Z4
    2.0*(3.0^2)/(3.0^2+X^2), X + Y --> Z5
    4.0*1/(X+1), X + Y --> Z6
end

for i = 1:100
    u = 5*rand(length(network1.syms))
    du1 = 3*rand(length(network1.syms)); du2 = du1;
    du1g = 2.5*rand(length(network1.syms),length(network1.reactions)); du2g = du1g;
    t = 9*rand()
    p = []

    @test network1.f(du1,u,p,t) == network2.f(du2,u,p,t)
    @test network1.g(du1g,u,p,t) == network2.g(du2g,u,p,t)
    @test network1.jumps[1].rate(u,p,t) == network2.jumps[1].rate(u,p,t)
    @test network1.jumps[2].rate(u,p,t) == network2.jumps[2].rate(u,p,t)
    @test network1.jumps[3].rate(u,p,t) == network2.jumps[3].rate(u,p,t)
end

model1 = @reaction_network rn begin
    (5,5000), X ↔ 0
end
model2 = @reaction_network rn η begin
    (5,5000), X ↔ 0
end
model3 = @reaction_network rn η begin
    (1,5000), X ↔ 0
end
function tmp_std(sol)
    vect = Vector{Float64}(undef,length(sol.u))
    for i = 1:length(sol.u)
        vect[i] = sol.u[i][1]
    end
    return std(vect)
end
prob1 = SDEProblem(model1,[1000.],(0.,200.))
prob2 = SDEProblem(model2,[1000.],(0.,200.), 0.1)
prob3 = SDEProblem(model3,[1000.],(0.,200.),10)
sol1 = solve(prob1, ImplicitEM())
sol2 = solve(prob2, ImplicitEM())
sol3 = solve(prob3,  ImplicitEM())
@test tmp_std(sol2) < tmp_std(sol1) < tmp_std(sol3)


text:
################################################################################
# NewtonCore
################################################################################

mutable struct NewtonCore{Vr,Vrt,SMj,SVhi,SVvi,SAr,SAj,Sd}
	res::Vr               # residual vector
	res_tmp::Vrt          # holds a temporary copy of the residual vector
	jac::SMj              # residual sparse jacobian
	probsize::ProblemSize # size of the problem
	horiz_inds::SVhi      # indices for each variable in an horizontal block
	verti_inds::SVvi      # indices for each variable in an vertical block
	res_sub::SAr          # Residual views dictionary
	jac_sub::SAj          # Jacobian views dictionary
	dyn::Sd               # Dictionary of the Jacobian dynamics indices
end

function NewtonCore(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S
	res = zeros(S)
	res_tmp = deepcopy(res)
	jac = spzeros(S,S)
	verti_inds = vertical_indices(probsize)
	horiz_inds = horizontal_indices(probsize)
	res_sub = residual_views(res, probsize, verti_inds)
	jac_sub = jacobian_views(jac, probsize, verti_inds, horiz_inds)
	dyn = dynamics_indices(probsize)
	TYPE = typeof.((res, res_tmp, jac, horiz_inds, verti_inds, res_sub, jac_sub, dyn))
	return NewtonCore{TYPE...}(res, res_tmp, jac, probsize,
		horiz_inds, verti_inds, res_sub, jac_sub, dyn)
end

################################################################################
# Indices
################################################################################

function vertical_indices(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	p = probsize.p
	mi = probsize.mi
	verti_inds = Dict()
	off = 0
	for i = 1:p
		for k = 1:N-1
			stamp = stampify(:opt, i, :x, 1, k+1)
			verti_inds[stamp] = SVector{n,Int}(off .+ (1:n))
			off += n
			stamp = stampify(:opt, i, :u, i, k)
			verti_inds[stamp] = SVector{mi[i],Int}(off .+ (1:mi[i]))
			off += mi[i]
		end
	end
	for k = 1:N-1
		stamp = stampify(:dyn, 1, :x, 1, k)
		verti_inds[stamp] = SVector{n,Int}(off .+ (1:n))
		off += n
	end
	return verti_inds
end

function horizontal_indices(probsize::ProblemSize)
	N = probsize.N
	n = probsize.n
	p = probsize.p
	mi = probsize.mi

	horiz_inds = Dict()
	off = 0
	for k = 1:N-1
		stamp = stampify(:x, 1, k+1)
		horiz_inds[stamp] = SVector{n,Int}(off .+ (1:n))
		off += n
		for i = 1:p
			stamp = stampify(:u, i, k)
			horiz_inds[stamp] = SVector{mi[i],Int}(off .+ (1:mi[i]))
			off += mi[i]
		end
		for i = 1:p
			stamp = stampify(:λ, i, k)
			horiz_inds[stamp] = SVector{n,Int}(off .+ (1:n))
			off += n
		end
	end
	return horiz_inds
end

function idx(core::NewtonCore, stamp::Stamp)
	vstamp = stampify(stamp.prob, stamp.i0, stamp.n1, stamp.i1, stamp.v1)
	hstamp = stampify(stamp.n2, stamp.i2, stamp.v2)
	verti = vertical_idx(core, vstamp)
	horiz = horizontal_idx(core, hstamp)
	return (verti, horiz)
end

function idx(verti_inds::Dict, horiz_inds::Dict, stamp::Stamp)
	vstamp = stampify(stamp.prob, stamp.i0, stamp.n1, stamp.i1, stamp.v1)
	hstamp = stampify(stamp.n2, stamp.i2, stamp.v2)
	verti = verti_inds[vstamp]
	horiz = horiz_inds[hstamp]
	return (verti, horiz)
end

function horizontal_idx(core::NewtonCore, stamp::HStamp)
	return core.horiz_inds[stamp]
end

function horizontal_idx(horiz_inds::Dict, stamp::HStamp)
	return horiz_inds[stamp]
end

function vertical_idx(core::NewtonCore, stamp::VStamp)
	return core.verti_inds[stamp]
end

function vertical_idx(verti_inds::Dict, stamp::VStamp)
	return verti_inds[stamp]
end

################################################################################
# SubArrays
################################################################################

function residual_views(res::SV, probsize::ProblemSize, verti_inds::Dict) where {SV}
	N = probsize.N
	p = probsize.p

	# Residual Views
	res_sub = Dict{VStamp,SubArray}()
	for prob ∈ (:dyn, :opt)
		for i0 = 1:p
			for n1 in (:x,:u)
				for i1 = 1:p
					for v1 = 1:N
						stamp = stampify(prob, i0, n1, i1, v1)
						if valid(stamp,N,p)
							res_sub[stamp] = view(res, vertical_idx(verti_inds, stamp))
						end
					end
				end
			end
		end
	end
	return res_sub
end

function jacobian_views(jac::SM, probsize::ProblemSize, verti_inds::Dict, horiz_inds::Dict) where {SM}
	N = probsize.N
	p = probsize.p

	# Jacobian Views
	jac_sub = Dict{Stamp,SubArray}()
	for prob ∈ (:dyn, :opt)
		for i0 = 1:p
			for n1 in (:x,:u)
				for i1 = 1:p
					for v1 = 1:N
						for n2 in (:x, :u, :λ)
							for i2 = 1:p
								for v2 = 1:N
									stamp = stampify(prob, i0, n1, i1, v1, n2, i2, v2)
									if valid(stamp,N,p)
										jac_sub[stamp] = view(jac, idx(verti_inds, horiz_inds, stamp)...)
									end
								end
							end
						end
					end
				end
			end
		end
	end
	return jac_sub
end

################################################################################
# Dynamics Indices
################################################################################

function dynamics_indices(probsize::ProblemSize)
	n = probsize.n
	p = probsize.p
	mi = probsize.mi
	pu = probsize.pu
	dyn = Dict()
	dyn[:x] = Dict()
	dyn[:u] = Dict()
	dyn[:x][1] = SVector{n,Int}(1:n)
	for i = 1:p
		dyn[:u][i] = SVector{mi[i],Int}(n .+ pu[i])
	end
	return dyn
end

################################################################################
# Player specific mask
################################################################################

function vertical_mask(core::NewtonCore, i::Int, splitted_state::Bool=false)
	vertical_mask(core.probsize, core.verti_inds, i, splitted_state=splitted_state)
end

function vertical_mask(probsize::ProblemSize, verti_inds::Dict, i::Int; splitted_state::Bool=false)
	N = probsize.N
	p = probsize.p
	pi = splitted_state ? probsize.pz[i] : Vector(1:probsize.n)

	msk = Vector{Int}()
	stamp = VStamp()

	# Select the optimality indices
	prob = :opt
	i0 = i
	for n1 in (:x,:u)
		for i1 = 1:p
			for v1 = 1:N
				stampify!(stamp, prob, i0, n1, i1, v1)
				if valid(stamp,N,p)
					ind = Vector(vertical_idx(verti_inds, stamp))
					(n1 == :x) && (ind = ind[pi])
					push!(msk, ind...)
				end
			end
		end
	end

	# Select the dynamics indices
	prob = :dyn
	i0 = 1
	n1 = :x
	i1 = 1
	for v1 = 1:N
		stampify!(stamp, prob, i0, n1, i1, v1)
		if valid(stamp,N,p)
			ind = Vector(vertical_idx(verti_inds, stamp)[pi])
			push!(msk, ind...)
		end
	end
	return msk
end

function horizontal_mask(core::NewtonCore, i::Int, splitted_state::Bool=false)
	horizontal_mask(core.probsize, core.horiz_inds, i, splitted_state=splitted_state)
end

function horizontal_mask(probsize::ProblemSize, horiz_inds::Dict, i::Int; splitted_state::Bool=false)
	N = probsize.N
	p = probsize.p
	pi = splitted_state ? probsize.pz[i] : Vector(1:probsize.n)

	msk = Vector{Int}()
	stamp = HStamp()

	# Select the state indices
	n2 = :x
	i2 = 1
	for v2 = 1:N
		stampify!(stamp, n2, i2, v2)
		if valid(stamp,N,p)
			ind = Vector(horizontal_idx(horiz_inds, stamp)[pi])
			push!(msk, ind...)
		end
	end

	# Select the control indices
	n2 = :u
	i2 = i
	for v2 = 1:N
		stampify!(stamp, n2, i2, v2)
		if valid(stamp,N,p)
			ind = Vector(horizontal_idx(horiz_inds, stamp))
			push!(msk, ind...)
		end
	end

	# Select the multipliers indices
	n2 = :λ
	i2 = i
	for v2 = 1:N
		stampify!(stamp, n2, i2, v2)
		if valid(stamp,N,p)
			ind = Vector(horizontal_idx(horiz_inds, stamp)[pi])
			push!(msk, ind...)
		end
	end
	return msk
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LibUV_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LibUV")
JLLWrappers.@generate_main_file("LibUV", UUID("183b4373-6708-53ba-ad28-60e28bb38547"))
end  # module LibUV_jll


text:
using Images, ImageView, Statistics

# function to return magnitude of image elements
function brightness(img_element::AbstractRGB)
    return mean((img_element.r + img_element.g + img_element.b))
end

function find_energy(img)
    energy_x = imfilter(brightness.(img), Kernel.sobel()[2])
    energy_y = imfilter(brightness.(img), Kernel.sobel()[1])
    return sqrt.(energy_x.^2 + energy_y.^2)
end

function find_energy_map(energy)
    energy_map = zeros(size(energy))
    energy_map[end,:] = energy[end,:]
    next_elements = zeros(Int, size(energy))

    for i = size(energy)[1]-1:-1:1
        for j = 1:size(energy)[2]
            left = max(1, j-1)
            right = min(j+1, size(energy)[2])

            local_energy. next_element = findmin(energy_map[i+1, left:right])

            energy_map[i,j] = local_energy + energy[i,j]

            next_elements[i,j] = next_element - 2

            if left == 1
                next_elements[i,j] += 1
            end
        end
    end
    return energy_map, next_elements
end

function find_seam_at(next_elements, element)
    seam = zeros(Int, size(next_elements)[1])
    seam[1] = element

    for i = 2:length(seam)
        seam[1] = seam[i-1] + next_elements[i, seam[i-1]]
    end
    return seam
end

# function to create seams and return seam of minimum energy
function find_seam(energy)
    energy_map, next_elements = find_energy_map(energy)
    _, min_element = findmin(energy_map[1,:])
    return find_seam_at(next_elements, min_element)
end

#function to remove seam
function remove_seam(img, seam)
    im_res = (size(img)[i], size(img)[2]-1)

    # preallocate image
    new_img = Array{RGB}(undef, img_res)

    for i = 1:length(seam)
        if seam[i] > 1 && seam[i] < size(img)[2]
            new_img[i, :] .= vcat(img[i, 1:seam[i]-1],
                                  img[i, seam[i]+1:end])
        elseif seam[i] == 1
            new_img[i, :] .= img[i,2:end]
        elseif seam[i] == size(img)[2]
            new_img[i, :] .= img[i, 1:end-1]
        end
    end

    return new_img
end

function seam_carver(img, res)
    if res < 0 || res > size(img)[2]
        error("resolution not acceptable")
    end

    for i = (1:size(img)[2] - res)
        energy = find_energy(img)
        seam = find_seam(energy)
        img = remove_seam(img, seam)
    end
    return img
end

max_stars_repo_name: prtpydv/seam-carver

max_issues_repo_name: prtpydv/seam-carver

max_forks_repo_name: prtpydv/seam-carver