text:
# build test case for using Pose2Pose2 to landmarks -- to help AprilTag

using RoME  #, Distributions
using RoMEPlotting, Gadfly

# const IIF = IncrementalInference
# const TU = TransformUtils

# true landmark locations
global gtl = Dict{Symbol, Vector{Float64}}()
gtl[:l1] = [1.0;-1.0;0]
gtl[:l2] = [2.0;-1.0;0]
gtl[:l3] = [3.0;-1.0;0]
gtl[:l4] = [4.0;-1.0;0]

gtl[:l11] = [1.0;1.0;0]
gtl[:l12] = [2.0;1.0;0]
gtl[:l13] = [3.0;1.0;0]
gtl[:l14] = [4.0;1.0;0]

# true pose locations
global gtp = Dict{Symbol,Vector{Float64}}()
gtp[:x0] = [0.0;0;0]
gtp[:x1] = [0.5;0;0]
gtp[:x2] = [1.0;0;0]
gtp[:x3] = [2.0;0;0]
gtp[:x4] = [2.5;0;0]
gtp[:x5] = [3.0;0;0]
gtp[:x6] = [3.0;0;0]

# calculate true measurements
global gtpt = Dict{Symbol, Vector{Float64}}()
for (ps,gp) in gtp, (ls,gl) in gtl
  gtpt[Symbol("$(ps)$(ls)")] = se2vee(SE2(gp) \ SE2(gl))
end

# Pose and Velocity noise model
global xposesig, yposesig = 0.5, 0.1
global thetasig = 0.3
global xvelsig, yvelsig = 0.2, 0.05

# Figure export folder
global currdirtime = now()
global imgdir = joinpath(ENV["HOME"], "Pictures", "testimgs", "$(currdirtime)")
mkdir(imgdir)

# construct the factor graph

global N=100
global fg = initfg()

addVariable!(fg, :x0, DynPose2; nanosecondtime=0)
addFactor!(fg, [:x0], DynPose2VelocityPrior(MvNormal(zeros(3),Matrix(Diagonal([0.01;0.01;0.001].^2))),
                                            MvNormal(zeros(2),Matrix(Diagonal([0.05;0.05].^2)))) )

IIF.doautoinit!(fg, :x0)

addVariable!(fg, :l1, Pose2)
addFactor!(fg, [:x0;:l1], DynPose2Pose2(MvNormal(gtpt[:x0l1],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

addVariable!(fg, :l11, Pose2)
addFactor!(fg, [:x0;:l11], DynPose2Pose2(MvNormal(gtpt[:x0l11],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

global pts = IIF.approxConv(fg, :x0l1f1, :l1, N=N)
# setVal!(fg, :l1, pts)

global pts = IIF.approxConv(fg, :x0l1f1, :x0, N=N)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTreeR!(fg, tree, N=N)

global psid = 0
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x1
addVariable!(fg, :x1, DynPose2; nanosecondtime=1000_000_000)
addFactor!(fg, [:x0;:x1],VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                          MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

#
# getVal(fg, :x0)
# getVal(fg, :x1)
# getSample(IIF.getfnctype(getVert(fg, :x0x1f1, nt=:fnc)),2)
# IIF.approxConv(fg, :x0x1f1, :x1)

addFactor!(fg, [:x1;:l1], DynPose2Pose2(MvNormal(gtpt[:x1l1],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addFactor!(fg, [:x1;:l11], DynPose2Pose2(MvNormal(gtpt[:x1l11],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l2, Pose2)
addFactor!(fg, [:x1;:l2], DynPose2Pose2(MvNormal(gtpt[:x1l2],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l12, Pose2)
addFactor!(fg, [:x1;:l12], DynPose2Pose2(MvNormal(gtpt[:x1l12],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTreeR!(fg, tree, N=N)

global psid = 1
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x2
addVariable!(fg, :x2, DynPose2; nanosecondtime=2000_000_000)
addFactor!(fg, [:x1;:x2], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                           MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

addFactor!(fg, [:x2;:l2], DynPose2Pose2(MvNormal(gtpt[:x2l2],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addFactor!(fg, [:x2;:l12], DynPose2Pose2(MvNormal(gtpt[:x2l12],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l3, Pose2)
addFactor!(fg, [:x2;:l3], DynPose2Pose2(MvNormal(gtpt[:x2l3],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 2
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# gtp[:x3] = [2.0;0;0]
addVariable!(fg, :x3, DynPose2; nanosecondtime=3000_000_000)
addFactor!(fg, [:x2;:x3], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                           MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

addFactor!(fg, [:x3;:l3], DynPose2Pose2(MvNormal(gtpt[:x3l3],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l13, Pose2)
addFactor!(fg, [:x3;:l13], DynPose2Pose2(MvNormal(gtpt[:x3l13],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 3
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x4
# gtp[:x4] = [2.5;0;0]
addVariable!(fg, :x4, DynPose2; nanosecondtime=4000_000_000)
addFactor!(fg, [:x3;:x4], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                           MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

addFactor!(fg, [:x4;:l3], DynPose2Pose2(MvNormal(gtpt[:x4l3],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addFactor!(fg, [:x4;:l13], DynPose2Pose2(MvNormal(gtpt[:x4l13],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l4, Pose2)
addFactor!(fg, [:x4;:l4], DynPose2Pose2(MvNormal(gtpt[:x4l4],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addVariable!(fg, :l14, Pose2)
addFactor!(fg, [:x4;:l14], DynPose2Pose2(MvNormal(gtpt[:x4l14],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 4
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x5
# gtp[:x5] = [3.0;0;0]
addVariable!(fg, :x5, DynPose2; nanosecondtime=5000_000_000)
addFactor!(fg, [:x4;:x5], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                           MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

addFactor!(fg, [:x5;:l4], DynPose2Pose2(MvNormal(gtpt[:x5l4],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
addFactor!(fg, [:x5;:l14], DynPose2Pose2(MvNormal(gtpt[:x5l14],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 5
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x6
# gtp[:x6] = [3.0;0;0]
addVariable!(fg, :x6, DynPose2; nanosecondtime=6000_000_000)
addFactor!(fg, [:x5;:x6], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([xposesig;yposesig;thetasig].^2))),
                                           MvNormal(zeros(2),Matrix(Diagonal([xvelsig;yvelsig].^2))) ))

addFactor!(fg, [:x6;:l4], DynPose2Pose2(MvNormal(gtpt[:x6l4],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))
# addFactor!(fg, [:x6;:l14], DynPose2Pose2(MvNormal(gtpt[:x6l14],Matrix(Diagonal([0.1;0.1;0.01].^2))) ))

# writeGraphPdf(fg)

tree = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 6
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

global coord = Coord.Cartesian(xmin=-2.0,xmax = 2.0)

global pl = plotPose2Vels(fg, :x0, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x0.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x1, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x1.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x2, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x2.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x3, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x3.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x4, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x4.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x5, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x5.png"),20cm, 10cm),pl)
global pl = plotPose2Vels(fg, :x6, coord=coord)
Gadfly.draw(PNG(joinpath(imgdir,"vel_x6.png"),20cm, 10cm),pl)

#


text:
using Snopt
using DelimitedFiles
using PyPlot
import ForwardDiff

# set up boundary constraint wrapper function
function boundary_wrapper(x)
    # include relevant globals
    global boundary_center
    global boundary_radius

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return boundary distances
    return ff.circle_boundary(boundary_center, boundary_radius, turbine_x, turbine_y)
end

# set up spacing constraint wrapper function
function spacing_wrapper(x)
    # include relevant globals
    global rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return spacing distances
    return 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)
end

# set up objective wrapper function
function aep_wrapper(x)
    # include relevant globals
    global turbine_z
    global rotor_diameter
    global hub_height
    global turbine_yaw
    global ct_models
    global generator_efficiency
    global cut_in_speed
    global cut_out_speed
    global rated_speed
    global rated_power
    global windresource
    global power_models
    global model_set
    global rotor_points_y
    global rotor_points_z
    global obj_scale

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_models, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z)

    # return the objective as an array
    return [AEP]
end

# set up optimization problem wrapper function
function wind_farm_opt(g, df, dg, x, deriv)

    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x)
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    dcdx = [ds_dx; db_dx]
    dg[:] = c[:]
    g[:] = c[:]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]
    dAEP_dx = -ForwardDiff.jacobian(aep_wrapper,x)
    df = dAEP_dx
    # set fail flag to false
    fail = false

    # return objective, constraint, and jacobian values
    # return AEP, c, dAEP_dx, dcdx, fail
    return AEP, fail
end

# import model set with wind farm and related details
cd("../")
include("../model_sets/model_set_6.jl")
cd("examples/")
# scale objective to be between 0 and 1
obj_scale = 1E-11

# set wind farm boundary parameters
boundary_center = [0.0,0.0]
boundary_radius = 300.0

# set globals for use in wrapper functions
global model_set
global rotor_points_y
global rotor_points_z
global turbine_z
global ambient_ti
global rotor_diameter
global boundary_center
global boundary_radius
global obj_scale

# initialize design variable array
x = [copy(turbine_x);copy(turbine_y)]

# report initial objective value
println("starting objective value: ", aep_wrapper(x)[1])

# add initial turbine location to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# set general lower and upper bounds for design variables
lb = zeros(length(x)) .- boundary_radius
ub = zeros(length(x)) .+ boundary_radius
lg = -Inf*ones(Int(nturbines + (nturbines-1)*nturbines/2))
ug = zeros(Int(nturbines + (nturbines-1)*nturbines/2))

# set up options for SNOPT
options = Dict{String, Any}()
options["Derivative option"] = 1
options["Verify level"] = 3
options["Major optimality tolerance"] = 1e-5
options["Major iteration limit"] = 1e6
options["Summary file"] = "summary.out"
options["Print file"] = "print.out"

# run and time optimization
t1 = time()
xopt, fopt, info, out = snopta(wind_farm_opt, x, lb, ub, lg, ug, [], [], options)
t2 = time()
clkt = t2-t2

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("end objective value: ", aep_wrapper(xopt))

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# add wind farm boundary to plot
plt.gcf().gca().add_artist(plt.Circle((boundary_center[1],boundary_center[2]), boundary_radius, fill=false,color="C2"))

# set up and show plot
axis("square")
xlim(-boundary_radius-200,boundary_radius+200)
ylim(-boundary_radius-200,boundary_radius+200)
plt.show()

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DomainView(domain, locations)

Return a view of `domain` at `locations`.

### Notes

This type implements the `AbstractDomain` interface.
"""
struct DomainView{T,N,
                  D<:AbstractDomain{T,N},
                  I<:AbstractVector{Int}} <: AbstractDomain{T,N}
  domain::D
  locations::I
end

npoints(view::DomainView) = length(view.locations)

coordinates!(buff::AbstractVector, view::DomainView, location::Int) =
  coordinates!(buff, view.domain, view.locations[location])

# ------------
# IO methods
# ------------
function Base.show(io::IO, view::DomainView{T,N,D,I}) where {T,N,
                                                             D<:AbstractDomain{T,N},
                                                             I<:AbstractVector{Int}}
  npts = npoints(view)
  print(io, "$npts DomainView{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", view::DomainView)
  println(io, view)
  Base.print_array(io, coordinates(view))
end


text:
# ---------------------------------------------------------------------------- #
#
#   genSTLmodel.jl
#
#   Write STL file for the model, based on just the mesh
#
#   intrico
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function genSTLmodel( mesh::Union{divido.Mesh3D,MeshF3D}, flnameBase::String )

    if length(flnameBase) > 4 && flnameBase[end-3:end] == ".stl"
        flnameBase = flnameBase[1:end-4]
    end

    nbound = length( mesh.boundtags )

    for bb in 1:nbound

        nface = count( mesh.fb[:,2] == bb )

        # open STL file
        fid = open( string(flnameBase, "_", mesh.boundtags[bb], ".stl" ), "w" )

        # write header
        for jj in 1:80
            write(fid,' ')
        end
        write(fid,UInt32( nface ))

        # generate all faces (and write)
        for jj in 1:size( mesh.fb, 1 )

            if mesh.fb[jj,2] != bb
                continue
            end

            indf = mesh.fb[jj,1]

            vert = SVector{3, SVector{3,Float64} }( SVector{3}( mesh.p[ mesh.f[indf,1], : ] ),
                                                    SVector{3}( mesh.p[ mesh.f[indf,2], : ] ),
                                                    SVector{3}( mesh.p[ mesh.f[indf,3], : ] ) )

            writeFacetSTLB( vert, fid )

        end

        # close STL
        close( fid )

    end

end


text:
module Integrators

    using CompactBasisFunctions
    using Documenter: @doc
    using LinearAlgebra
    using OffsetArrays
    using QuadratureRules
    using RungeKutta
    using SimpleSolvers

    using GeometricBase
    using GeometricBase.Config
    using GeometricBase.Utils
    using GeometricEquations
    using GeometricEquations: _get_v, _get_f, _get_v̄, _get_f̄

    using ..Discontinuities
    using ..Solutions

    import RungeKutta: nstages

    export Extrapolation,
           EulerExtrapolation, EulerExtrapolationODE,
           MidpointExtrapolation, MidpointExtrapolationODE, MidpointExtrapolationIODE, MidpointExtrapolationPODE,
           HermiteExtrapolation

    include("integrators/extrapolation/extrapolation.jl")
    include("integrators/extrapolation/aitken_neville.jl")
    include("integrators/extrapolation/euler.jl")
    include("integrators/extrapolation/hermite.jl")
    include("integrators/extrapolation/midpoint.jl")

    export InitialGuess, InitialGuessODE, InitialGuessIODE, InitialGuessPODE,
           initialize!

    include("integrators/initial_guess/initial_guess_ode.jl")
    include("integrators/initial_guess/initial_guess_iode.jl")
    include("integrators/initial_guess/initial_guess_pode.jl")

    export AbstractCoefficients, AbstractTableau

    include("integrators/abstract_coefficients.jl")
    include("integrators/abstract_tableau.jl")

    export Integrator, DeterministicIntegrator, StochasticIntegrator
    export ODEIntegrator, DAEIntegrator, SDEIntegrator,
           PODEIntegrator, PDAEIntegrator, PSDEIntegrator,
           IODEIntegrator, IDAEIntegrator,
           HODEIntegrator, HDAEIntegrator,
           LODEIntegrator, LDAEIntegrator,
           SPSDEIntegrator

    export IntegratorCache, IntegratorConstructor
    export integrate, integrate!, integrate_step!, equation, timestep
    export function_stages!#, NonlinearFunctionParameters

    include("integrators/abstract_integrator.jl")
    include("integrators/integrator_cache.jl")
    include("integrators/integrators_common.jl")

    export IntegratorExplicitEuler

    include("integrators/various/integrators_explicit_euler.jl")

    export AbstractIntegratorRK, AbstractIntegratorIRK, AbstractIntegratorPRK, IntegratorRK

    export get_symplectic_conjugate_coefficients, symplecticize,
           check_symplecticity, symplecticity_conditions,
           check_symmetry, compute_symplecticity_error

    include("integrators/rk/abstract_integrator_rk.jl")
    include("integrators/rk/tableaus.jl")

    export IntegratorERK
    export IntegratorDIRK
    export IntegratorFIRK
    export IntegratorFIRKimplicit, IntegratorMidpointImplicit, IntegratorSRKimplicit

    include("integrators/rk/integrators_erk.jl")
    include("integrators/rk/integrators_dirk.jl")
    include("integrators/rk/integrators_firk.jl")
    include("integrators/rk/integrators_firk_implicit.jl")
    include("integrators/rk/integrators_midpoint_implicit.jl")
    include("integrators/rk/integrators_srk_implicit.jl")

    export IntegratorEPRK
    export IntegratorIPRK
    export IntegratorPRKimplicit
    export IntegratorFLRK

    include("integrators/rk/integrators_eprk.jl")
    include("integrators/rk/integrators_iprk.jl")
    include("integrators/rk/integrators_prk_implicit.jl")
    include("integrators/rk/integrators_flrk.jl")

    export IntegratorPGLRK, CoefficientsPGLRK

    include("integrators/rk/integrators_pglrk.jl")

    export IntegratorSplitting,
           IntegratorComposition,
           IntegratorExactODE,
           AbstractTableauSplitting,
           TableauSplitting,
           TableauSplittingGS,
           TableauSplittingNS,
           TableauSplittingSS

    include("integrators/splitting/integrators_exact_ode.jl")
    include("integrators/splitting/splitting_tableau.jl")
    include("integrators/splitting/integrators_composition.jl")
    include("integrators/splitting/integrators_splitting.jl")

    export IntegratorCGVI, IntegratorDGVI, IntegratorDGVIEXP,
           IntegratorDGVIPI, IntegratorDGVIP0, IntegratorDGVIP1

    include("integrators/cgvi/integrators_cgvi.jl")
    include("integrators/dgvi/integrators_dgvi.jl")
    include("integrators/dgvi/integrators_dgvi_experimental.jl")
    include("integrators/dgvi/integrators_dgvi_path_integral.jl")
    include("integrators/dgvi/integrators_dgvi_projection_initial.jl")
    include("integrators/dgvi/integrators_dgvi_projection_final.jl")

    export IntegratorDVIA, IntegratorDVIB,
           IntegratorCMDVI, IntegratorCTDVI

    include("integrators/dvi/integrators_dvi_a.jl")
    include("integrators/dvi/integrators_dvi_b.jl")
    include("integrators/dvi/integrators_cmdvi.jl")
    include("integrators/dvi/integrators_ctdvi.jl")

    include("integrators/SPARK.jl")
    include("integrators/VPRK.jl")

    include("integrators/integrators.jl")

    function __init__()
        default_params = (
            (:ig_extrapolation, HermiteExtrapolation),
            (:ig_extrapolation_stages, 5),
            (:int_show_progress_nmin,  1000),
        )

        for param in default_params
            add_config(param...)
        end
    end

end


text:
module SimpleBirthDeathProcess
  using LinearAlgebra
  using SpecialFunctions

  export
  loglik,
  mle,
  observation_continuous_time,
  observation_discrete_time_even,
  observation_discrete_time_uneven,
  rand_continuous,
  rand_discrete,
  trans_prob

  include("common/common.jl")
  include("simple/simple.jl")
end


text:
module GenieDevTools

using Genie
using Genie.Renderers.Json, Genie.Renderers.Html

using Revise
using Dates
using RemoteREPL

import Stipple, Stipple.Pages

const defaultroute = "/_devtools_"
const _HOOKS = Function[]

include("route_helpers.jl")
import .RouteHelpers as RH

function hook!(f::Function)
  push!(_HOOKS, f)
end

function runhooks()
  for f in _HOOKS
    Base.invokelatest(f)
  end
end

function register_routes(defaultroute = defaultroute) :: Nothing
  RH.errors(defaultroute)
  RH.dir(defaultroute)
  RH.edit(defaultroute)
  RH.save(defaultroute)
  RH.exec(defaultroute)
  RH.id(defaultroute)
  RH.log(defaultroute)
  RH.exit(defaultroute)
  RH.up(defaultroute)
  RH.down(defaultroute)
  RH.pages(defaultroute)
  RH.startrepl(defaultroute)

  nothing
end

end


text:
module Spiral
using Random

export load_data

function load_data(;seed=1984)
    Random.seed!(seed)  # 関数直後の!は関数が破壊的変更を行うことを示す
    @show seed
    N = 100             # クラスごとのサンプルサイズ
    DIM = 2             # データの要素数
    CLS_DIM = 3         # クラス数

    x = zeros(N * CLS_DIM, DIM)
    t = zeros(Integer, N * CLS_DIM, CLS_DIM)

    for j = 0:CLS_DIM - 1 # Juliaではインデックスが1始まりのため
        for i = 0:N - 1
            rate = i / N
            radius = 1.0 * rate
            θ = j * 4.0 + 4.0 * rate + Random.randn() * 0.2

            ix = N * j + i
            x[ix + 1, :] = [radius * sin(θ) radius * cos(θ)]

            t[ix + 1, j + 1] = 1
        end
    end

    return x, t
end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export simplegrid

"""
    simplegrid(;kwargs...)
    simplegrid(children::Any;kwargs...)
    simplegrid(children_maker::Function;kwargs...)

A SimpleGrid component.
Responsive grid where each item takes equal amount of space. For more information, see: https://mantine.dev/core/simple-grid/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): <Col /> components only
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `breakpoints` (optional): Breakpoints data to change items per row and spacing based on max-width. breakpoints has the following type: lists containing elements 'maxWidth', 'cols', 'spacing'.
Those elements have the following types:
  - `maxWidth` (Real; required)
  - `cols` (Real; required)
  - `spacing` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional)
- `className` (String; optional): Often used with CSS to style elements with common properties
- `cols` (Real; required): Default amount of columns, used when none of breakpoints can be applied
- `spacing` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Default spacing between columns, used when none of breakpoints can be applied
- `style` (Dict; optional): Inline style override
"""
function simplegrid(; kwargs...)
        available_props = Symbol[:children, :id, :breakpoints, :className, :cols, :spacing, :style]
        wild_props = Symbol[]
        return Component("simplegrid", "SimpleGrid", "dash_mantine_components", available_props, wild_props; kwargs...)
end

simplegrid(children::Any; kwargs...) = simplegrid(;kwargs..., children = children)
simplegrid(children_maker::Function; kwargs...) = simplegrid(children_maker(); kwargs...)


text:
using ConditionalDists
using Flux

batch = 10000
xlength = 300
zlength = 150

f = Dense(zlength,xlength*2)
p = CMeanVarGaussian{DiagVar}(f)

z = randn(Float32, zlength, batch)
x = randn(Float32, xlength, batch)

@time logpdf(p, x, z)
@time logpdf(p, x, z)
@time logpdf(p, x, z)


text:
using GenericSchur: trsylvester!, adjtrsylvester!
using GenericSVD

Random.seed!(1234)

function checkgsyl(m=4,n=2,Ty=Float64)
    rctr = Ty(0.5) * ((Ty <: Complex) ? (1+1im) : 1)
    A = triu(rand(Ty,m,m) .- rctr) + 2I
    B = triu(rand(Ty,n,n) .- rctr) + I
    X = rand(Ty,m,n)

    # separate spectra so problem is not too hard
    D = triu(rand(Ty,m,m) .- rctr) + 4I
    E = triu(rand(Ty,n,n) .- rctr) + I
    Y = rand(Ty,m,n)
    C = A*X-Y*B
    C2 = copy(C)
    F = D*X-Y*E
    F2 = copy(F)

    X2,Y2,scale = trsylvester!(A,B,C2,D,E,F2)
    fwd2a = norm(scale*X-X2) / norm(X)
    fwd2b = norm(scale*Y-Y2) / norm(Y)
    bwd2a = norm(A*X2-Y2*B-scale*C)/((norm(A)+norm(B))*norm(X))
    bwd2b = norm(D*X2-Y2*E-scale*F)/((norm(D)+norm(E))*norm(X))
    if verbosity[] > 0
        println("our    scale $scale error norms fwd: ",fwd2a," ",fwd2b)
        println("                                bwd: ", bwd2a," ",bwd2b)
    end
    tol = 2
    stol = tol * max(m,n) * eps(real(Ty))
    # forward errors for Sylvester are fiddly
    if max(m,n) < 10
        @test fwd2a < sqrt(stol)
        @test fwd2b < sqrt(stol)
    end
    @test bwd2a < tol
    @test bwd2b < tol

end

function checkgasyl(m=4,n=2,Ty=Float64)
    A = triu(rand(Ty,m,m) .- 0.5) + 2I
    B = triu(rand(Ty,n,n) .- 0.5) + I
    X = rand(Ty,m,n)

    # separate spectra so problem is not too hard
    D = triu(rand(Ty,m,m) .- 0.5) + 4I
    E = triu(rand(Ty,n,n) .- 0.5) + I
    Y = rand(Ty,m,n)
    C = A'*X + D'*Y
    C2 = copy(C)
    F = -X*B' - Y*E'
    F2 = copy(F)

    X2, Y2, scale = adjtrsylvester!(A,B,C2,D,E,F2)
    fwd2a = norm(scale*X-X2) / norm(X)
    fwd2b = norm(scale*Y-Y2) / norm(Y)
    bwd2a = norm(A'*X2+D'*Y2-scale*C)/((norm(A)+norm(B))*norm(X))
    bwd2b = norm(X2*B'+Y2*E'+scale*F)/((norm(D)+norm(E))*norm(X))
    if verbosity[] > 0
        println("our    scale $scale error norms fwd: ",fwd2a," ",fwd2b)
        println("                                bwd: ", bwd2a," ",bwd2b)
    end

    tol = 2
    stol = tol * max(m,n) * eps(real(Ty))
    if max(m,n) < 10
        @test fwd2a < sqrt(stol)
        @test fwd2b < sqrt(stol)
    end
    @test bwd2a < tol
    @test bwd2b < tol

end

for Ty in [ComplexF64, Float64, Complex{BigFloat}]
    @testset "gen. sylvester $Ty" begin
        checkgsyl(2,4,Ty)
        checkgsyl(4,2,Ty)
        checkgsyl(24,24,Ty)
    end
    @testset "gen. adj. sylvester $Ty" begin
        checkgasyl(2,4,Ty)
        checkgasyl(4,2,Ty)
        checkgasyl(24,24,Ty)
    end
end

function checkord(A::Matrix{Ty}, B::Matrix{Ty}, tol=10) where {Ty<:Complex}
    n = size(A,1)
    ulp = eps(real(Ty))

    GS = schur(A,B)
    S2 = copy(GS.S)
    Q2 = copy(GS.Q)
    T2 = copy(GS.T)
    Z2 = copy(GS.Z)
    select = fill(true,n)
    for i=1:(n>>1)
        select[i] = false
    end
    # we extend the LinearAlgebra function, but want to test
    # our version with BLAS types:
    S2,T2,a,b,Q2,Z2 = invoke(LinearAlgebra._ordschur!,
                             Tuple{StridedMatrix{T}, StridedMatrix{T},
                                   StridedMatrix{T}, StridedMatrix{T},
                                   Union{Vector{Bool},BitVector}}
                             where T <: Complex{Tr} where Tr <: AbstractFloat,
                           S2,T2,Q2,Z2,select)

    # usual tests for Schur
    @test norm(tril(S2,-1)) / (n * norm(A) * ulp) < tol
    @test norm(tril(T2,-1)) / (n * norm(A) * ulp) < tol
    # @test all(tril(S.S,-1) .== 0)
    # @test all(tril(S.T,-1) .== 0)
    # test 2: norm(A - S.Q * S.T * S.Z') / (n * norm(A) * ulp) < tol
    @test norm(A - Q2 * S2 * Z2') / (n * norm(A) * ulp) < tol
    @test norm(B - Q2 * T2 * Z2') / (n * norm(B) * ulp) < tol
    # test 3: S.Z is unitary: norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - Z2 * Z2') / (n * ulp) < tol
    @test norm(I - Q2 * Q2') / (n * ulp) < tol

    # make sure we got the ones we asked for
    # WARNING: this test may fail if selected ev are ill-conditioned
    nwanted = count(select)
    wwanted = [GS.values[j] for j=1:n if select[j]]
    errs = zeros(nwanted)
    for i=1:nwanted
        w = a[i] / b[i]
        t,j = findmin(abs.(wwanted .- w))
        errs[i] =  t / (ulp + abs(w))
        wwanted = deleteat!(wwanted, j)
    end
    @test all(errs .< tol)
end

for Ty in [ComplexF64, Complex{BigFloat}]
    @testset "gen. ordschur $Ty" begin
        n = 16
        A = rand(Ty,n,n)
        rctr = Ty(0.5) * ((Ty <: Complex) ? (1+1im) : 1)
        B = rand(Ty,n,n) .- rctr + I
        checkord(A,B)
    end
end

# WARNING: KPTest was adapted from a more thorough validation study.
# For interesting arguments the tests don't pass as written, but results
# agree with LAPACK. With benign arguments, however, it is a good
# test of functionality.

# Test with CHK4 matrices from Kågström and Poromaa, LAWN 87, 1994.
function KPTest(α,β,x,y,dtype,Ty=ComplexF64)
    Db = I
    if dtype == 1
        Da = diagm(0=>Ty.([1,2,3,4,5])) + α * I
        atrue = [1,2,3,4,5] .+ α
    else
        Da = diagm(0=>Ty.([1,1,1,1+α,1+α]), 1=>Ty.([-1,0,0,1+β]),
                   -1=>Ty.([1,0,0,-1-β]))
        atrue = [1-1im,1+1im,1,1+α-(1+β)*im,1+α+(1+β)*im]
    end
    YH = diagm(0=>ones(Ty,5), 1=>Ty.([0,-y,0,0]), 2=>Ty.([-y,y,0]),
               3=>Ty.([y,-y]), 4=>Ty.([-y]))
    X = diagm(0=>ones(Ty,5), 1=>Ty.([0,x,0,0]), 2=>Ty.([-x,-x,0]),
              3=>Ty.([-x,-x]), 4=>Ty.([x]))
    A = YH \ (Da / X)
    B = YH \ (Db / X)
    S  = schur(A,B)
    aest = S.α ./ S.β
    for j=1:5
        λtrue = atrue[j]
        _,idx = findmin(abs.(aest .- λtrue))
        sel = falses(5)
        sel[idx] = true
        if idx == 1
            S2 = S
        else
            S2 = ordschur(S,sel)
        end
        pl, pr = eigvalscond(S2, 1)
        sep =  subspacesep(S2, 1)

        # known values from construction
        if dtype == 1 || j==3
            y = YH[j,:] # a column vector!
            x = X[:,j]
        else
            if j==1
                y = YH[1,:] + im * YH[2,:]
                x = X[:,1] + im * X[:,2]
            elseif j==2
                y = YH[1,:] - im * YH[2,:]
                x = X[:,1] - im * X[:,2]
            elseif j==4
                y = YH[4,:] - im * YH[5,:]
                x = X[:,4] - im * X[:,5]
            elseif j==5
                y = YH[4,:] + im * YH[5,:]
                x = X[:,4] + im * X[:,5]
            end
        end

        rcond = (sqrt(abs(dot(y,(A * x)))^2 + abs(dot(y,(B * x)))^2)
                 / (norm(x)*norm(y'))) / sqrt(3.0)

        n = 5
        Inm1 = Matrix{Ty}(I,n-1,n-1)
        Zl = hcat(vcat(λtrue * Inm1, Inm1),
                  vcat(-S2.S[2:n,2:n], -S2.T[2:n,2:n]))
        difl = minimum(svdvals(Zl))
        if verbosity[] > 0
            println("λ true $λtrue error $(λtrue - aest[idx])")
            println("ours   pl $(pl) pr $(pr) sep $(sep)")
            println("true rcond $rcond fudged $(rcond)  Difl $(difl)")
        end

        @test abs(log10(pl / rcond)) < 2
        @test abs(log10(pr / rcond)) < 2
        @test abs(log10(sep[1]/difl)) < 2
        @test abs(log10(sep[2]/difl)) < 2

    end
end

for Ty in [ComplexF64, Complex{BigFloat}]
    @testset "gen. condition $Ty" begin
        α = 1.0
        β = 1.0
        x = 1.0
        y = 1.0
        dtype = 2
        KPTest(α,β,x,y,dtype,Ty)
    end
end


text:
include(homedir() * "/.julia/dev/WordleSolver/extras/worst_words.jl")
reps = 10
T = run(reps)


text:
include("fenkar.jl")
using StaticArrays, DifferentialEquations, ModelingToolkit, Plots

#= 	fenkar.jl has an in-place FenKar function (fenkar!), but it is slower than the
	out-of-place version using StaticArrays (fenkar), go figure.
		Basic version is  ~33712x faster than real time (Tsit5())
		SArray version is ~48685x faster than real time (Tsit5())
	So we use the SArray one (fenkar), unless managing static arrays gets annoying,
	which they often do, then just use fenkar!.
=#

tspan = Fenkar.tspan
u0 = rand(Float64,3)
p  = [ 22.0, 333.0,   40.0, 10.0, 65.0, 1000.0, 0.12, 12.5, 25.0, 10.0, 0.13, 0.025,  0.85]

prob = ODEProblem(Fenkar.fenkar!, u0, tspan, p)

sol = solve(prob, Tsit5())
plot(sol)

using OrdinaryDiffEq, ForwardDiff, MinimallyDisruptiveCurves, Statistics, Plots, LinearAlgebra, LaTeXStrings

nom_prob = prob
nom_sol = sol
t = collect(range(tspan[1], tspan[2]; length=1025))

## Model features of interest are mean prey population, and max predator population (over time)
function features(p)
  prob = remake(nom_prob; p=p)
  sol  = solve(prob, Tsit5(); saveat = t)
  return [maximum(sol[1,:]), sum(abs.(sol[1,:])), sum(sol[1,:].^2)]
  #return sol[1,1:2:100]
end

nom_features = features(p)

## loss function, we can take as l2 difference of features vs nominal features
function loss(p)
  prob = remake(nom_prob; p=p)
  p_features = features(p)
  loss = sum(abs2, p_features - nom_features)
  return loss
end

## gradient of loss function
function lossgrad(p,g)
  g[:] = ForwardDiff.gradient(p) do p
    loss(p)
  end
  return loss(p)
end

## package the loss and gradient into a DiffCost structure
cost = DiffCost(loss, lossgrad)

#=
We evaluate the hessian once only, at p.
Why? to find locally insensitive directions of parameter perturbation
The small eigenvalues of the Hessian are one easy way of defining these directions
=#
#which md curve to plot
#which_dir = 2
hess0 = ForwardDiff.hessian(loss,p)
ev(i) = eigen(hess0).vectors[:,i]
which_dir = 6 #argmin(abs.(eigvals(hess0)))

## Now we set up a minimally disruptive curve, with nominal parameters p and initial direction ev(1)
init_dir = ev(which_dir); momentum = 1.; span = (-100.,100.)
curve_prob = MDCProblem(cost, p, init_dir, momentum, span)
@time mdc = evolve(curve_prob, Tsit5)

function sol_at_p(p)
  prob = remake(nom_prob; p=p)
  sol = solve(prob, Tsit5())
end

cost_vec = [mdc.cost(el) for el in eachcol(trajectory(mdc))]

#tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi
p1 = plot(mdc; pnames=["tsi","tv1m","tv2m","tvp","twm","twp","td","to","tr","xk","uc","uv","ucsi"])
p2 = plot(distances(mdc), log.(cost_vec), ylabel = "log(cost)", xlabel = "distance", title = "cost over MD curve");
mdc_plot = plot(p1,p2, layout=(2,1), size = (800,800))

nominal_trajectory = plot(sol_at_p(mdc(0.)[:states]), label = ["u" "v" "w"])
perturbed_trajectory = plot(sol_at_p(mdc(75.)[:states]), label = ["u" "v" "w"])
traj_comparison = plot(nominal_trajectory, perturbed_trajectory, layout = (2,1), size = (800,800), xlabel = "time", ylabel = "state variable")

plot(mdc_plot, traj_comparison, layout=(1,2), size=(1200,900))
savefig("./fig.pdf")


text:
using DeferredAcceptance

#=  Most school-choice algorithms start with each side ranking the other
    in terms of preferability. First we input the student preferences.
    Each column is a student, and each row is a school. So, the 2 in cell
    (3, 4) means that student 4 has named school 3 as her 2nd choice. =#
students = [3 3 4 3 4 3 3 3 3 4;
            4 4 3 4 3 4 4 4 4 3;
            2 1 2 2 2 1 2 2 2 2;
            1 2 1 1 1 2 1 1 1 1]

#=  Now we input the school preferences over the students. Notice that the
    shape is transposed: each column is a school, and row a student. This
    format improves computation speed. Also notice that now we have ties:
    for example, school 1 has four students tied for 1st place. You can
    also write the first column as [1, 1, 2, 2, 4, 1, 2, 3, 3, 1]. Basically,
    they just need to be positive integers in descending order of preference. =#
schools = [1 5 7 5;
           1 1 1 1;
           5 1 1 1;
           5 1 1 1;
           10 9 7 8;
           1 5 4 5;
           5 1 4 1;
           8 5 7 8;
           8 9 7 8;
           1 5 4 5]

# Finally we have the capacities, or number of students each school can accept.
capacities = [3, 2, 2, 3]

#=  Let's use the DA algorithm to find a stable match. The DA algorithm requires
    that both students and schools have strict preferences, so we need to run
    some random tiebreaking to make the school preferences strict. We will use
    STB, which is usually the best choice.        =#
schools_tiebroken = STB(schools)

# Now we run DA.
assn, ranks = DA(students, schools_tiebroken, capacities; verbose=true)

# Make sure the match is stable.
@assert isstable(students, schools, capacities, assn)

#=  The first output is the assignment. Student 1 goes to school 1, student 2
    goes to school 3, etc. If the market had too many students, we would have
    some students assigned to school 5, which means no school accepted them. =#
println(assn)
# [1, 3, 4, 4, 2, 3, 4, 1, 1, 2]

#=  The second output is the rank each student gave to their assigned school.
    So student 1 got her 3rd choice, student 2 got her 2nd, etc.    =#
println(ranks)
# [3, 1, 1, 1, 3, 1, 1, 3, 3, 3]

# A common measure of student disutility
println(sum(ranks))
# 22

#=  Since the tiebreaking mechanism involves randomness, you might have
    found a slightly different assignment.

    Now let's try something different. Suppose the schools don't have any
    preferences over the students at all (as in New Orleans). Then we can
    start with an arbitrary assignment and exchange Pareto-improving pairs
    to get a better one. This is called top trading cycles and is a classic
    solution to the "housing market" problem. In school choice, it produces
    an assignment that tends to maximize student welfare.       =#
assn, ranks = TTC_match(students, capacities; verbose=true)

# Same output format
println(assn)
# [2, 3, 4, 4, 2, 1, 4, 1, 1, 3]
println(ranks)
# [4, 1, 1, 1, 3, 3, 1, 3, 3, 2]

# Student disutility
println(sum(ranks))
# 21

# Looks like TTC did a little better, but what are the tradeoffs?
isstable(students, schools, capacities, assn; verbose=true) # Returns false
# Student feas. :  true
# School feas.  :  true
# Stability     :  false

#=  Let's try something different: the nonatomic DA model. Under this model, we
    have a measure (usually a percentage) of students associated with each preference
    list. Instead of admitting individual students, stable matches can be characterized
    by score cutoffs at each school.            =#

# Possible preference lists: Most students have (1, 2, 3, 4), but some prefer 2 to 1.
students = [1 2;
            2 1;
            3 3;
            4 4]

# Percentage of students with each preference list.
students_dist = [0.75, 0.25]

#=  School capacities. Schools have no preference list; assume student preferability is
    independently distributed in each group.        =#
capacities = [0.2, 0.2, 0.2, 0.2]

# Run nonatomic DA.
assn, rdist, cutoffs = DA_nonatomic(students, students_dist, nothing, capacities;
                                    verbose=true, return_cutoffs=true)

# Equivalently,
# cutoffs = DA_nonatomic_lite(students, students_dist, capacities)

# As expected, school 1 is most selective.
println(cutoffs)
# [0.7873499783936067, 0.7620499351812947, 0.666666666665277, 0.49999999999882583]

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
struct Hydro end

const curtailconstraints = PowerSimulations.powerconstraints

function constructdevice!(category::Type{Hydro}, transmission::Type{CopperPlatePowerModel}, m::JuMP.AbstractModel, devices_netinjection::T, sys::PowerSystems.PowerSystem, constraints::Array{<:Function}=[powerconstraints]) where T <: JumpExpressionMatrix

    devices = [d for d in sys.generators.hydro if (d.available == true && !isa(d, HydroFix))]

    if !isempty(devices)

        pre, devices_netinjection = activepowervariables(m, devices_netinjection, devices, sys.time_periods)

        for c in constraints

            m = c(m, devices, sys.time_periods)

        end

    else

        @warn "Hydro dispatch variables created without constraints, the problem might be unbounded"

    end

    return m, devices_netinjection

end


text:
module Scale

using Color
using Compose
using DataArrays
using DataStructures
using Gadfly

import Gadfly: element_aesthetics, isconcrete, concrete_length,
               nonzero_length, formatter, setfield!, set

include("color.jl")

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#   datas: Zero or more data objects. (Yes, I know "datas" is not a real word.)
#
# Returns:
#   nothing
#
function apply_scales(scales,
                      aess::Vector{Gadfly.Aesthetics},
                      datas::Gadfly.Data...)
    for scale in scales
        apply_scale(scale, aess, datas...)
    end

    for (aes, data) in zip(aess, datas)
        aes.titles = data.titles
    end
end

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   datas: Zero or more data objects.
#
# Returns:
#   A vector of Aesthetics of the same length as datas containing scaled data.
#
function apply_scales(scales, datas::Gadfly.Data...)
    aess = Gadfly.Aesthetics[Gadfly.Aesthetics() for _ in datas]
    apply_scales(scales, aess, datas...)
    aess
end

# Transformations on continuous scales
immutable ContinuousScaleTransform
    f::Function     # transform function
    finv::Function  # f's inverse

    # A function taking one or more values and returning an array of
    # strings.
    label::Function
end

function identity_formatter(xs::AbstractArray; format=:auto)
    fmt = formatter(xs, fmt=format)
    [fmt(x) for x in xs]
end

const identity_transform =
    ContinuousScaleTransform(identity, identity, identity_formatter)

function log10_formatter(xs::AbstractArray; format=:plain)
    fmt = formatter(xs, fmt=format)
    [@sprintf("10%s", fmt(x)) for x in xs]
end

const log10_transform =
    ContinuousScaleTransform(log10, x -> 10^x, log10_formatter)

function log2_formatter(xs::AbstractArray; format=:plain)
    fmt = formatter(xs, fmt=format)
    [@sprintf("2%s", fmt(x)) for x in xs]
end

const log2_transform =
    ContinuousScaleTransform(log2, x -> 2^x, log2_formatter)

function ln_formatter(xs::AbstractArray; format=:plain)
    fmt = formatter(xs, fmt=format)
    [@sprintf("e%s", fmt(x)) for x in xs]
end

const ln_transform =
    ContinuousScaleTransform(log, exp, ln_formatter)

function asinh_formatter(xs::AbstractArray; format=:plain)
    fmt = formatter(xs, fmt=format)
    [@sprintf("asinh(%s)", fmt(x)) for x in xs]
end

const asinh_transform =
    ContinuousScaleTransform(asinh, sinh, asinh_formatter)

function sqrt_formatter(xs::AbstractArray; format=:plain)
    fmt = formatter(xs, fmt=format)
    [@sprintf("√%s", fmt(x)) for x in xs]
end

const sqrt_transform = ContinuousScaleTransform(sqrt, x -> x^2, sqrt_formatter)

# Continuous scale maps data on a continuous scale simple by calling
# `convert(Float64, ...)`.
immutable ContinuousScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
    trans::ContinuousScaleTransform
    minvalue
    maxvalue
    format

    function ContinuousScale(vars::Vector{Symbol},
                             trans::ContinuousScaleTransform;
                             minvalue=nothing, maxvalue=nothing,
                             format=nothing)
        new(vars, trans, minvalue, maxvalue, format)
    end
end

function make_labeler(scale::ContinuousScale)
    if scale.format == nothing
        scale.trans.label
    else
        function f(xs)
            scale.trans.label(xs, format=scale.format)
        end
    end
end

const x_vars = [:x, :xmin, :xmax, :xintercept]
const y_vars = [:y, :ymin, :ymax, :yintercept, :middle,
                :upper_fence, :lower_fence, :upper_hinge, :lower_hinge]

function continuous_scale_partial(vars::Vector{Symbol},
                                  trans::ContinuousScaleTransform)
    function f(;minvalue=nothing, maxvalue=nothing, format=nothing)
        ContinuousScale(vars, trans, minvalue=minvalue, maxvalue=maxvalue,
                        format=format)
    end
end

# Commonly used scales.
const x_continuous = continuous_scale_partial(x_vars, identity_transform)
const y_continuous = continuous_scale_partial(y_vars, identity_transform)
const x_log10      = continuous_scale_partial(x_vars, log10_transform)
const y_log10      = continuous_scale_partial(y_vars, log10_transform)
const x_log2       = continuous_scale_partial(x_vars, log2_transform)
const y_log2       = continuous_scale_partial(y_vars, log2_transform)
const x_log        = continuous_scale_partial(x_vars, ln_transform)
const y_log        = continuous_scale_partial(y_vars, ln_transform)
const x_asinh      = continuous_scale_partial(x_vars, asinh_transform)
const y_asinh      = continuous_scale_partial(y_vars, asinh_transform)
const x_sqrt       = continuous_scale_partial(x_vars, sqrt_transform)
const y_sqrt       = continuous_scale_partial(y_vars, sqrt_transform)

const size_continuous = continuous_scale_partial([:size], identity_transform)

function element_aesthetics(scale::ContinuousScale)
    return scale.vars
end

# Apply a continuous scale.
#
# Args:
#   scale: A continuos scale.
#   datas: Zero or more data objects.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#
# Return:
#   nothing
#
function apply_scale(scale::ContinuousScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            if getfield(data, var) === nothing
                continue
            end

            T = Any
            for (i, d) in enumerate(getfield(data, var))
                if isconcrete(d)
                    T = typeof(scale.trans.f(d))
                    break
                end
            end

            ds = DataArray(T, length(getfield(data, var)))
            for (i, d) in enumerate(getfield(data, var))
                if isconcrete(d)
                    ds[i] = scale.trans.f(d)
                    i += 1
                else
                    ds[i] = d
                end
            end

            setfield!(aes, var, ds)

            if var in x_vars
                label_var = :x_label
            elseif var in y_vars
                label_var = :y_label
            else
                label_var = symbol(@sprintf("%s_label", string(var)))
            end

            if in(label_var, set(names(aes)))
                setfield!(aes, label_var, make_labeler(scale))
            end
        end

        if scale.minvalue != nothing
            if scale.vars === x_vars
                aes.xviewmin = scale.trans.f(scale.minvalue)
            elseif scale.vars === y_vars
                aes.yviewmin = scale.trans.f(scale.minvalue)
            end
        end

        if scale.maxvalue != nothing
            if scale.vars === x_vars
                aes.xviewmax = scale.trans.f(scale.maxvalue)
            elseif scale.vars === y_vars
                aes.yviewmax = scale.trans.f(scale.maxvalue)
            end
        end
    end
end

# Reorder the levels of a pooled data array
function reorder_levels(da::PooledDataArray, order::AbstractVector)
    level_values = levels(da)
    if length(order) != length(level_values)
        error("Discrete scale order is not of the same length as the data's levels.")
    end
    permute!(level_values, order)
    return PooledDataArray(da, level_values)
end

function discretize_make_pda(values::Vector, levels=nothing)
    if levels == nothing
        return PooledDataArray(values)
    else
        return PooledDataArray(convert(Vector{eltype(levels)}, values), levels)
    end
end

function discretize_make_pda(values::DataArray, levels=nothing)
    if levels == nothing
        return PooledDataArray(values)
    else
        return PooledDataArray(convert(DataArray{eltype(levels)}, values), levels)
    end
end

function discretize_make_pda(vaues::Range1, levels=nothing)
    if levels == nothing
        return PooledDataArray(collect(values))
    else
        return PooledDataArray(collect(values), levels)
    end
end

function discretize_make_pda(values::PooledDataArray, levels=nothing)
    if levels == nothing
        return values
    else
        return PooledDataArray(values, convert(Vector{eltype(values)}, levels))
    end
end

function discretize(values, levels=nothing, order=nothing,
                    preserve_order=true)
    if levels == nothing
        if preserve_order
            levels = OrderedSet()
            for value in values
                push!(levels, value)
            end
            da = discretize_make_pda(values, collect(levels))
        else
            da = discretize_make_pda(values)
        end
    else
        da = discretize_make_pda(values, levels)
    end

    if order != nothing
        return reorder_levels(da, order)
    else
        return da
    end
end

immutable DiscreteScaleTransform
    f::Function
end

immutable DiscreteScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}

    # If non-nothing, give values for the scale. Order will be respected and
    # anything in the data that's not represented in values will be set to NA.
    levels::Union(Nothing, AbstractVector)

    # If non-nothing, a permutation of the pool of values.
    order::Union(Nothing, AbstractVector)

    # If true, order levels as they appear in the data.
    preserve_order::Bool

    function DiscreteScale(vals::Vector{Symbol}; levels=nothing, order=nothing,
                           preserve_order::Bool=true)
        new(vals, levels, order, preserve_order)
    end
end

const discrete = DiscreteScale

element_aesthetics(scale::DiscreteScale) = scale.vars

function x_discrete(; levels=nothing, order=nothing, preserve_order=true)
    return DiscreteScale(x_vars, levels=levels, order=order,
                         preserve_order=preserve_order)
end

function y_discrete(; levels=nothing, order=nothing, preserve_order=true)
    return DiscreteScale(y_vars, levels=levels, order=order,
                         preserve_order=preserve_order)
end

function apply_scale(scale::DiscreteScale, aess::Vector{Gadfly.Aesthetics},
                     datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            label_var = symbol(@sprintf("%s_label", string(var)))

            if getfield(data, var) === nothing
                continue
            end

            disc_data = discretize(getfield(data, var), scale.levels, scale.order)
            setfield!(aes, var, PooledDataArray(int64(disc_data.refs)))

            # The leveler for discrete scales is a closure over the discretized data.
            function labeler(xs)
                lvls = levels(disc_data)
                vals = {1 <= x <= length(lvls) ? lvls[x] : "" for x in xs}
                if all([isa(val, FloatingPoint) for val in vals])
                    format = formatter(vals)
                    [format(val) for val in vals]
                else
                    [string(val) for val in vals]
                end
            end

            if in(label_var, set(names(aes)))
                setfield!(aes, label_var, labeler)
            end
        end
    end
end

immutable DiscreteColorScale <: Gadfly.ScaleElement
    f::Function # A function f(n) that produces a vector of n colors.

    # If non-nothing, give values for the scale. Order will be respected and
    # anything in the data that's not represented in values will be set to NA.
    levels::Union(Nothing, AbstractVector)

    # If non-nothing, a permutation of the pool of values.
    order::Union(Nothing, AbstractVector)

    # If true, order levels as they appear in the data
    preserve_order::Bool

    function DiscreteColorScale(f::Function; levels=nothing, order=nothing,
                                preserve_order=true)
        new(f, levels, order, preserve_order)
    end
end

function element_aesthetics(scale::DiscreteColorScale)
    [:color]
end

# Common discrete color scales
function discrete_color_hue(; levels=nothing, order=nothing)
    DiscreteColorScale(
        h -> convert(Vector{ColorValue},
             distinguishable_colors(h, ColorValue[LCHab(70, 60, 240)],
                                    lchoices=Float64[65, 70, 75, 80],
                                    cchoices=Float64[0, 50, 60, 70],
                                    hchoices=linspace(0, 330, 24),
                                    transform=c -> deuteranopic(c, 0.5))),
        levels=levels, order=order)
end

const discrete_color = discrete_color_hue

function discrete_color_manual(colors...; levels=nothing, order=nothing)
    cs = ColorValue[color(c) for c in colors]
    function f(n)
        distinguishable_colors(n, cs)
    end
    DiscreteColorScale(f, levels=levels, order=order)
end

function apply_scale(scale::DiscreteColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    levelset = OrderedSet()
    for (aes, data) in zip(aess, datas)
        if data.color === nothing
            continue
        end
        for d in data.color
            if !isna(d)
                push!(levelset, d)
            end
        end
    end

    if scale.levels == nothing
        scale_levels = [levelset...]
        if !scale.preserve_order
            sort!(scale_levels)
        end
    else
        scale_levels = scale.levels
    end
    if scale.order != nothing
        permute!(scale_levels, scale.order)
    end
    colors = convert(Vector{ColorValue}, scale.f(length(scale_levels)))

    color_map = {color => string(label)
                 for (label, color) in zip(scale_levels, colors)}
    function labeler(xs)
        [color_map[x] for x in xs]
    end

    for (aes, data) in zip(aess, datas)
        if data.color === nothing
            continue
        end
        ds = discretize(data.color, scale_levels)
        colorvals = Array(ColorValue, nonzero_length(ds.refs))
        i = 1
        for k in ds.refs
            if k != 0
                colorvals[i] = colors[k]
                i += 1
            end
        end

        colored_ds = PooledDataArray(colorvals, colors)
        aes.color = colored_ds

        aes.color_label = labeler
        aes.color_key_colors = colors
    end
end

immutable ContinuousColorScale <: Gadfly.ScaleElement
    # A function of the form f(p) where 0 <= p <= 1, that returns a color.
    f::Function

    minvalue
    maxvalue

    function ContinuousColorScale(f::Function; minvalue=nothing, maxvalue=nothing)
        new(f, minvalue, maxvalue)
    end
end

element_aesthetics(::ContinuousColorScale) = [:color]

function continuous_color_gradient(;minvalue=nothing, maxvalue=nothing)
    ContinuousColorScale(
        lab_gradient(LCHab(20, 44, 262), LCHab(100, 44, 262)),
        minvalue=minvalue, maxvalue=maxvalue)
end

const continuous_color = continuous_color_gradient

function apply_scale(scale::ContinuousColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    cmin = Inf
    cmax = -Inf
    for data in datas
        if data.color === nothing
            continue
        end

        for c in data.color
            if c === NA
                continue
            end

            c = convert(Float64, c)
            if c < cmin
                cmin = c
            end

            if c > cmax
                cmax = c
            end
        end
    end

    if cmin == Inf || cmax == -Inf
        return nothing
    end

    if scale.minvalue != nothing
        cmin = scale.minvalue
    end

    if scale.maxvalue  != nothing
        cmax = scale.maxvalue
    end

    cmin, cmax = promote(cmin, cmax)

    ticks, viewmin, viewmax = Gadfly.optimize_ticks(cmin, cmax)
    if ticks[1] == 0 && cmin >= 1
        ticks[1] = 1
    end

    cmin = ticks[1]
    cmax = ticks[end]
    cspan = cmax != cmin ? cmax - cmin : 1

    for (aes, data) in zip(aess, datas)
        if data.color === nothing
            continue
        end

        nas = [c === NA for c in data.color]
        cs = Array(ColorValue, length(data.color))
        for (i, c) in enumerate(data.color)
            if c === NA
                continue
            end
            cs[i] = scale.f((convert(Float64, c) - cmin) / cspan)
        end

        aes.color = DataArray(cs, nas)

        color_key_colors = Array(ColorValue, 0)
        color_key_labels = Array(String, 0)

        num_steps = 1
        tick_labels = identity_formatter(ticks)
        for (i, j, label) in zip(ticks, ticks[2:end], tick_labels)
            r = (i - cmin) / cspan
            push!(color_key_colors, scale.f(r))
            push!(color_key_labels, label)

            for step in 1:num_steps
                k = i + (j - i) * (step / (1 + num_steps))
                r = (k - cmin) / cspan
                push!(color_key_colors, scale.f(r))
                push!(color_key_labels, "")
            end
        end
        push!(color_key_colors, scale.f((ticks[end] - cmin) / cspan))
        push!(color_key_labels, tick_labels[end])

        color_key_color_label =
            [c => l for (c, l) in zip(color_key_colors, color_key_labels)]
        function labeler(xs)
            [get(color_key_color_label, x, "") for x in xs]
        end

        aes.color_label = labeler
        aes.color_key_colors = color_key_colors
        reverse!(aes.color_key_colors)
        aes.color_key_continuous = true
    end
end

# Label scale is always discrete, hence we call it 'label' rather
# 'label_discrete'.
immutable LabelScale <: Gadfly.ScaleElement
end

function apply_scale(scale::LabelScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        if data.label === nothing
            continue
        end

        aes.label = discretize(data.label)
    end
end

element_aesthetics(::LabelScale) = [:label]

const label = LabelScale

# Scale applied to grouping aesthetics.
immutable GroupingScale <: Gadfly.ScaleElement
    var::Symbol
end

function xgroup(; levels=nothing, order=nothing)
    return DiscreteScale([:xgroup], levels=levels, order=order)
end

function ygroup(; levels=nothing, order=nothing)
    return DiscreteScale([:ygroup], levels=levels, order=order)
end

# Catchall scale for when no transformation of the data is necessary
immutable IdentityScale <: Gadfly.ScaleElement
    var::Symbol
end

function element_aesthetics(scale::IdentityScale)
    return [scale.var]
end

function apply_scale(scale::IdentityScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip (aess, datas)
        if getfield(data, scale.var) === nothing
            continue
        end

        setfield!(aes, scale.var, getfield(data, scale.var))
    end
end

function func()
    return IdentityScale(:func)
end

end # module Scale


text:
@testset "Util" begin
    @testset "Perpendicular vector" begin
        for i in 1 : 100
            vec = randn(SVector{3, Float64})
            perp = Rotations.perpendicular_vector(vec)
            @test norm(perp) >= maximum(abs.(vec))
            @test isapprox(dot(vec, perp), 0.; atol = 1e-10)
        end
        let vec = randn(SVector{3, Float64})
            allocs = @allocated Rotations.perpendicular_vector(vec)
            @test allocs == 0
        end
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

ioslength(io::IOBuffer) = (io.seekable ? io.size : nb_available(io))

let io = IOBuffer()
@test eof(io)
@test_throws EOFError read(io,UInt8)
@test write(io,"abc") === 3
@test isreadable(io)
@test iswritable(io)
@test isopen(io)
@test ioslength(io) == 3
@test position(io) == 3
@test eof(io)
seek(io, 0)
@test read(io,UInt8) == convert(UInt8, 'a')
a = Array{UInt8}(2)
@test read!(io, a) == a
@test a == UInt8['b','c']
@test String(io) == "abc"
seek(io, 1)
truncate(io, 2)
@test position(io) == 1
@test !eof(io)
seekend(io)
@test position(io) == 2
truncate(io, 0)
@test position(io) == 0
truncate(io, 10)
@test position(io) == 0
@test all(io.data .== 0)
@test write(io,Int16[1,2,3,4,5,6]) === 12
seek(io, 2)
truncate(io, 10)
@test ioslength(io) == 10
io.readable = false
@test_throws ArgumentError read!(io,UInt8[0])
truncate(io, 0)
@test write(io,"boston\ncambridge\n") > 0
@test String(take!(io)) == "boston\ncambridge\n"
@test String(take!(io)) == ""
@test write(io, Complex{Float64}(0)) === 16
@test write(io, Rational{Int64}(1//2)) === 16
close(io)
@test_throws ArgumentError write(io,UInt8[0])
@test_throws ArgumentError seek(io,0)
@test eof(io)
end

let io = IOBuffer("hamster\nguinea pig\nturtle")
@test position(io) == 0
@test readline(io) == "hamster"
@test readstring(io) == "guinea pig\nturtle"
@test_throws EOFError read(io,UInt8)
seek(io,0)
@test read(io,UInt8) == convert(UInt8, 'h')
@test_throws ArgumentError truncate(io,0)
@test_throws ArgumentError write(io,UInt8(0))
@test_throws ArgumentError write(io,UInt8[0])
@test String(take!(io)) == "hamster\nguinea pig\nturtle"
@test String(take!(io)) == "hamster\nguinea pig\nturtle" #should be unchanged
close(io)
end

let io = PipeBuffer()
@test_throws EOFError read(io,UInt8)
@test write(io,"pancakes\nwaffles\nblueberries\n") > 0
@test position(io) == 0
@test readline(io) == "pancakes"
Base.compact(io)
@test readline(io) == "waffles"
@test write(io,"whipped cream\n") > 0
@test readline(io) == "blueberries"
@test_throws ArgumentError seek(io,0)
@test_throws ArgumentError truncate(io,0)
@test readline(io) == "whipped cream"
@test write(io,"pancakes\nwaffles\nblueberries\n") > 0
@test readlines(io) == String["pancakes", "waffles", "blueberries"]
write(io,"\n\r\n\n\r \n") > 0
@test readlines(io, chomp=false) == String["\n", "\r\n", "\n", "\r \n"]
write(io,"\n\r\n\n\r \n") > 0
@test readlines(io, chomp=true) == String["", "", "", "\r "]
@test write(io,"α\nβ\nγ\nδ") > 0
@test readlines(io, chomp=false) == String["α\n","β\n","γ\n","δ"]
@test write(io,"α\nβ\nγ\nδ") > 0
@test readlines(io, chomp=true) == String["α", "β", "γ", "δ"]
@test readlines(IOBuffer(""), chomp=false) == []
@test readlines(IOBuffer(""), chomp=true) == []
@test readlines(IOBuffer("first\nsecond"), chomp=false) == String["first\n", "second"]
@test readlines(IOBuffer("first\nsecond"), chomp=true) == String["first", "second"]
Base.compact(io)
@test position(io) == 0
@test ioslength(io) == 0
Base.ensureroom(io,50)
@test position(io) == 0
@test ioslength(io) == 0
@test length(io.data) == 50
Base.ensureroom(io,10)
@test ioslength(io) == 0
@test length(io.data) == 50
io.maxsize = 75
Base.ensureroom(io,100)
@test ioslength(io) == 0
@test length(io.data) == 75
seekend(io)
@test ioslength(io) == 0
@test position(io) == 0
write(io,zeros(UInt8,200))
@test ioslength(io) == 75
@test length(io.data) == 75
write(io,1)
@test ioslength(io) == 75
@test length(io.data) == 75
write(io,[1,2,3])
@test ioslength(io) == 75
@test length(io.data) == 75
skip(io,1)
@test write(io,UInt8(104)) === 1
skip(io,3)
@test write(io,b"apples") === 3
skip(io,71)
@test write(io,'y') === 1
@test readstring(io) == "happy"
@test eof(io)
write(io,zeros(UInt8,73))
write(io,'a')
write(io,'b')
write(io,'c')
write(io,'d')
write(io,'e')
@test ioslength(io) == 75
@test length(io.data) == 75
@test position(io) == 0
skip(io,72)
@test String(take!(io)) == "\0ab"
@test String(take!(io)) == ""

# issues 4021
print(io, true)
close(io)
end

# issue 5453
let io=IOBuffer("abcdef")
a = Array{UInt8}(1024)
@test_throws EOFError read!(io,a)
@test eof(io)
end

@test isempty(readlines(IOBuffer(), chomp=false))

# issue #8193
let io=IOBuffer("asdf")
    @test position(skip(io, -1)) == 0
    @test position(skip(io, 6)) == 4
    @test position(seek(io, -1)) == 0
    @test position(seek(io, 6)) == 4
end

# issue #10658
let io=IOBuffer("hello")
    @test position(skip(io, 4)) == 4
    @test position(skip(io, 10)) == 5
    @test position(skip(io, -2)) == 3
    @test position(skip(io, -3)) == 0
    @test position(skip(io, -3)) == 0
end

# pr #11554
let io=IOBuffer(SubString("***αhelloworldω***",4,16)), io2 = IOBuffer(b"goodnightmoon", true, true)
    @test read(io, Char) == 'α'
    @test_throws ArgumentError write(io,"!")
    @test_throws ArgumentError write(io,'β')
    a = Array{UInt8}(10)
    @test read!(io, a) === a
    @test String(a) == "helloworld"
    @test read(io, Char) == 'ω'
    @test_throws EOFError read(io,UInt8)
    skip(io, -3)
    @test readstring(io) == "dω"
    @test String(io) == "αhelloworldω"
    @test_throws ArgumentError write(io,"!")
    @test take!(io) == b"αhelloworldω"
    seek(io, 2)
    seekend(io2)
    write(io2, io)
    @test readstring(io) == ""
    @test readstring(io2) == ""
    @test String(take!(io)) == "αhelloworldω"
    seek(io2, 0)
    truncate(io2, io2.size - 2)
    @test readstring(io2) == "goodnightmoonhelloworld"
    seek(io2, 0)
    write(io2, io2)
    @test readstring(io2) == ""
    @test String(io2) == "goodnightmoonhelloworld"
end

# issue #11917
# (previous tests triggered this sometimes, but this should trigger nearly all the time)
let io = IOBuffer(0)
   write(io, ones(UInt8, 1048577))
end

let bstream = BufferStream()
    @test isopen(bstream)
    @test isreadable(bstream)
    @test iswritable(bstream)
    @test nb_available(bstream) == 0
    @test sprint(show, bstream) == "BufferStream() bytes waiting:$(nb_available(bstream.buffer)), isopen:true"
    a = rand(UInt8,10)
    write(bstream,a)
    @test !eof(bstream)
    @test flush(bstream) === nothing
    b = read(bstream,UInt8)
    @test a[1] == b
    b = read(bstream,UInt8)
    @test a[2] == b
    c = zeros(UInt8,8)
    @test nb_available(bstream) == 8
    @test !eof(bstream)
    read!(bstream,c)
    @test c == a[3:10]
    @test close(bstream) === nothing
    @test eof(bstream)
    @test nb_available(bstream) == 0
end

@test flush(IOBuffer()) === nothing # should be a no-op

# pr #19461
let io = IOBuffer()
    @test Base.buffer_writes(io) === io
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: fcard/julia


text:
# This file was generated, do not modify it. # hide
using DelimitedFiles
using Plots
img = readdlm("/tmp/benchC.txt",',',Int)
M,N = img |> size
p=heatmap(1:N, 1:M, img, aspect_ratio=1, size=(400,400))
savefig(p, joinpath(@OUTPUT, "mandC.png")) #hide


text:
using ClusterTrees
using Test

N = ClusterTrees.SimpleTrees.TreeNode
nodes = N[
    N(8,"a"),
        N(2,"b"),
            N(0,"c"),
            N(0,"d"),
        N(4,"e"),
            N(0,"f"),
            N(2,"g"),
                N(0,"h"),
                N(0,"i")]

tree = ClusterTrees.SimpleTrees.SimpleTree(nodes)
ClusterTrees.print_tree(tree)

@test data(tree, root(tree)) == "a"
@test collect(data(tree,ch) for ch in ClusterTrees.children(tree, root(tree))) == ["b","e"]
@test collect(data(tree,lv) for lv in ClusterTrees.leaves(tree, root(tree))) == ["c","d","f","h","i"]

mtree = ClusterTrees.mutable(tree)
ClusterTrees.print_tree(mtree)

# n1 = ClusterTrees.root(mtree)
# n2 = first(ClusterTrees.children(mtree,n1))
# ClusterTrees.SimpleTrees.insert_child!(mtree, n2, "z")

chd = ClusterTrees.children(mtree, ClusterTrees.root(mtree))
i1 = ClusterTrees.SimpleTrees.start(chd)
n1, i2 = ClusterTrees.SimpleTrees.next(chd, i1)
@test depth(mtree, n1) == 2
@test ClusterTrees.parent(mtree, n1) == root(mtree)
@test depth(mtree, root(mtree)) == 1
n2, i3 = ClusterTrees.SimpleTrees.next(chd, i2)
ClusterTrees.SimpleTrees.insert!(chd, "Q", i3)

ClusterTrees.print_tree(tree)
A = [ClusterTrees.data(mtree, n) for n in collect(ClusterTrees.DepthFirstIterator(mtree, root(mtree)))]
@test A == ["c","d","b","f","h","i","g","e","Q","a"]

const N2 = ClusterTrees.PointerBasedTrees.Node
nodes2 = N2[
    N2("a",2,-1,-1,2),
        N2("b",2,5,1,3),
            N2("c",0,4,2,-1),
            N2("d",0,-1,2,-1),
        N2("e",2,-1,1,6),
            N2("f",0,7,5,-1),
            N2("g",2,-1,5,8),
                N2("h",0,9,7,-1),
                N2("i",0,-1,7,-1)]
tree2 = ClusterTrees.PointerBasedTrees.PointerBasedTree(nodes2,1)

@test data(tree2) == "a"
@test collect(data(tree2,ch) for ch in ClusterTrees.children(tree2)) == ["b","e"]
@test collect(data(tree2,lv) for lv in ClusterTrees.leaves(tree2)) == ["c","d","f","h","i"]

chd = ClusterTrees.children(tree2, ClusterTrees.root(tree2))
i1 = ClusterTrees.SimpleTrees.start(chd)
n1, i2 = ClusterTrees.SimpleTrees.next(chd, i1)
n2, i3 = ClusterTrees.SimpleTrees.next(chd, i2)
insert!(chd, "Q", i3)

ClusterTrees.print_tree(tree)
A = [ClusterTrees.data(mtree, n) for n in collect(ClusterTrees.DepthFirstIterator(mtree, root(mtree)))]
@test A == ["c","d","b","f","h","i","g","e","Q","a"]


text:
digits_prec(prec::Integer) = floor(Int, prec * log(2) / log(10))

function _string(x::MagOrRef)
    Libc.flush_cstdio()
    Base.flush(stdout)
    io = IOStream("arb")
    original_stdout = stdout
    out_rd, out_wr = redirect_stdout()
    try
        ccall(@libarb(mag_print), Cvoid, (Ref{mag_struct},), x)
        Libc.flush_cstdio()
    finally
        redirect_stdout(original_stdout)
        close(out_wr)
    end
    return read(out_rd, String)
end

function Base.show(io::IO, x::MagOrRef)
    if isdefined(Main, :IJulia) && Main.IJulia.inited
        print(io, Float64(x))
    else
        print(io, _string(x))
    end
end

Base.show(io::IO, x::ArfOrRef) = print(io, BigFloat(x))

function Base.show(io::IO, x::ArbOrRef)
    cstr = ccall(
        @libarb(arb_get_str),
        Ptr{UInt8},
        (Ref{arb_struct}, Int, UInt),
        x,
        digits_prec(precision(x)),
        UInt(0),
    )
    print(io, unsafe_string(cstr))
    ccall(@libflint(flint_free), Nothing, (Ptr{UInt8},), cstr)
end

function Base.show(io::IO, x::AcbOrRef)
    show(io, realref(x))
    if !iszero(imagref(x))
        print(io, " + ")
        show(io, imagref(x))
        print(io, "im")
    end
end

function Base.show(io::IO, poly::T) where {T<:Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}}
    if (T == ArbPoly || T == AcbPoly) && iszero(poly)
        print(io, "0")
    end

    # We only want to print up until the last non-zero element, for
    # ArbPoly and AcbPoly this is given by the degree, for ArbSeries
    # and AcbSeries we get this by taking the degree of the underlying
    # polynomial.
    N = degree(cstruct(poly))
    for i = 0:N
        x = poly[i]
        if !iszero(x)
            str =
                ifelse(!isreal(x), "(", "") *
                "$x" *
                ifelse(!isreal(x), ")", "") *
                ifelse(i > 0, "⋅x", "") *
                ifelse(i > 1, "^$i", "") *
                ifelse(i != N, " + ", "")
            print(io, str)
        end
    end

    if T == ArbSeries || T == AcbSeries
        str =
            ifelse(iszero(poly), "", " + ") *
            "𝒪(x" *
            ifelse(degree(poly) == 0, "", "^$(degree(poly) + 1)") *
            ")"
        print(io, str)
    end
end

for ArbT in (Mag, MagRef, Arf, ArfRef, Arb, ArbRef, Acb, AcbRef)
    arbf = Symbol(cprefix(ArbT), :_, :print)
    @eval begin
        function string_nice(
            x::$ArbT,
            digits::Integer = digits_prec(precision(x)),
            flags::UInt = UInt(0),
        )
            Libc.flush_cstdio()
            Base.flush(stdout)
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            try
                printn(x, digits, flags)
                Libc.flush_cstdio()
            finally
                redirect_stdout(original_stdout)
                close(out_wr)
            end
            return read(out_rd, String)
        end
    end
end

for ArbT in (Mag, MagRef, Arf, ArfRef, Arb, ArbRef)
    @eval begin
        function load_string!(x::$ArbT, str::AbstractString)
            res = load!(x, str)
            iszero(res) || throw(ArgumentError("could not load $str as " * $(string(ArbT))))
            return x
        end

        function dump_string(x::$ArbT)
            char_ptr = dump(x)
            str = unsafe_string(char_ptr)
            ccall(@libflint(flint_free), Cvoid, (Cstring,), char_ptr)
            return str
        end
    end
end

for T in [
    :MagLike,
    :ArfLike,
    :ArbLike,
    :AcbLike,
    :ArbVectorLike,
    :AcbVectorLike,
    :ArbMatrixLike,
    :AcbMatrixLike,
]
    @eval Base.show(io::IO, ::Type{$T}) = print(io, $(QuoteNode(T)))
end


text:
using ToyFHE
using OffsetArrays
using Test
using GaloisFields
using Primes
using StructArrays

n = 2048
p1 = nextprime(Int128(2)^50 + 1, 1; interval=2n)
p2 = nextprime(p1+2n, interval=2n)

ℛ = let CT = ToyFHE.CRTEncoded{2, Tuple{GaloisField(p1), GaloisField(p2)}}
    ζ₄₀₉₆ = GaloisFields.minimal_primitive_root(CT, 2n)
    ToyFHE.NegacyclicRing{CT, n}(ζ₄₀₉₆)
end

p3 = nextprime(p2+2n, 1; interval=2n)
p4 = nextprime(p3+2n, 1; interval=2n)
p5 = nextprime(p4+2n, 1; interval=2n)
p6 = nextprime(p5+2n, 1; interval=2n)

ℛbig = let CT = ToyFHE.CRTEncoded{4, Tuple{GaloisField(p3), GaloisField(p4), GaloisField(p5), GaloisField(p6)}}
    ζ₄₀₉₆ = GaloisFields.minimal_primitive_root(CT, 2n)
    ToyFHE.NegacyclicRing{CT, n}(ζ₄₀₉₆)
end

ℛplain = plaintext_space(ℛ, 53)

params = BFVParams(
    ℛ,
    ℛbig,
    ℛplain,
    1,
    3.2,
    div(ToyFHE.NTT.modulus(ToyFHE.NTT.coefftype(ℛ)), ToyFHE.NTT.modulus(ToyFHE.NTT.coefftype(ℛplain)))
)
kp = keygen(params)

plain = zero(plaintext_space(params))
plain[0] = 6

c = encrypt(kp, plain)
@test decrypt(kp, c)[0] == 6

let y = c*c
    @test decrypt(kp, y)[0] == 0x24
end


text:
module GalacticOptim

using Reexport
@reexport using DiffEqBase
using Requires
using DiffResults, ForwardDiff, Zygote, ReverseDiff, Tracker, FiniteDiff
using Optim, Flux
using Logging, ProgressLogging, Printf, ConsoleProgressMonitor, TerminalLoggers, LoggingExtras
using ArrayInterface

using ForwardDiff: DEFAULT_CHUNK_THRESHOLD
import DiffEqBase: OptimizationProblem, OptimizationFunction, AbstractADType

include("solve.jl")
include("function.jl")

export solve

export BBO, CMAEvolutionStrategyOpt

end # module


text:
module Bridges

include("Variable/Variable.jl")
include("Constraint/Constraint.jl")

end


text:
using Documenter, MIToS

include("literate.jl")

makedocs(
    doctest=true,
	format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    ),
    sitename="MIToS",
    authors="Diego Javier Zea",
    modules=[MIToS],
    pages=[
        "Home" => "index.md",
        "Installation.md",
        "Example.md",
        "Modules" => [  "MSA.md",
                        "Information.md",
                        "SIFTS.md",
                        "PDB.md",
                        "Pfam.md"
                          ],
        "Cookbook" => [ "01_Change_B_factors.md",
                        "02_Linking_structural_and_evolutionary_information.md",
                        "03_RMSF.md"],
        "API" => [  "MSA_API.md",
                    "Information_API.md",
                    "SIFTS_API.md",
                    "PDB_API.md",
                    "Pfam_API.md",
                    "Utils_API.md"          ],
        "Scripts.md",
        ]
)

deploydocs(
    repo="github.com/diegozea/MIToS.jl.git",
    target="build",
    deps=nothing,
    make=nothing
)


text:
using Sobol

"""
    Landau( α, kx)

Test structure to initialize a particles distribtion for
Landau damping test case in 1D1V and 1D2V

"""
struct LandauDamping
    alpha::Float64
    kx::Float64
end

"""
    sample!(d, pg)

Sampling from a probability distribution to initialize a Landau damping in
1D2V space.

```math
f_0(x,v,t) = \\frac{n_0}{2π v_{th}^2} ( 1 + \\alpha cos(k_x x))
 exp( - \\frac{v_x^2+v_y^2}{2 v_{th}^2})
```
The newton function solves the equation ``P(x)-r=0`` with Newton’s method
```math
x^{n+1} = x^n – (P(x)-(2\\pi r / k)/f(x)
```
with
```math
P(x) = \\int_0^x (1 + \\alpha cos(k_x y)) dy = x + \\frac{\\alpha}{k_x} sin(k_x x)
```
"""
function sample!(d::LandauDamping, pg::ParticleGroup{1,2})
    alpha, kx = d.alpha, d.kx

    function newton(r)
        x0, x1 = 0.0, 1.0
        r *= 2π / kx
        while (abs(x1 - x0) > 1e-12)
            p = x0 + alpha * sin(kx * x0) / kx
            f = 1 + alpha * cos(kx * x0)
            x0, x1 = x1, x0 - (p - r) / f
        end
        return x1
    end

    s = SobolSeq(2)
    nbpart = pg.n_particles

    for i in 1:nbpart
        v = sqrt(-2 * log((i - 0.5) / nbpart))
        r1, r2 = Sobol.next!(s)
        θ = r1 * 2π
        set_x!(pg, i, newton(r2))
        set_v!(pg, i, [v * cos(θ), v * sin(θ)])
        set_weights!(pg, i, 2 * pi / kx / nbpart)
    end
end


text:
function test_hinge()
    ϑ = rand()
    l = Hinge(ϑ)
    @testset "hinge loss" begin
        @test typeof(l) <: AbstractSurrogate
        @test l.ϑ == ϑ
        @test l.value.([-1/ϑ - 1, 0, -1/ϑ + 1])    ≈ [0, 1, ϑ]
        @test l.gradient.([-1/ϑ - 1, 0, -1/ϑ + 1]) ≈ [0, ϑ, ϑ]
    end
end

function test_quadratic()
    ϑ = rand()
    l = Quadratic(ϑ)
    @testset "quadratic loss" begin
        @test typeof(l) <: AbstractSurrogate
        @test l.ϑ == ϑ
        @test l.value.([-1/ϑ - 1, 0, -1/ϑ + 1])    ≈ [0, 1, ϑ^2]
        @test l.gradient.([-1/ϑ - 1, 0, -1/ϑ + 1]) ≈ [0, 2*ϑ, 2*ϑ^2]
    end
end


text:
@testset "polynomial" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    c = rand(rng)
    @testset "LinearKernel" begin
        k = LinearKernel()
        @test kappa(k, x) ≈ x
        @test k(v1, v2) ≈ dot(v1, v2)
        @test kappa(LinearKernel(), x) == kappa(k, x)
        @test metric(LinearKernel()) == KernelFunctions.DotProduct()
        @test metric(LinearKernel(; c=c)) == KernelFunctions.DotProduct()
        @test repr(k) == "Linear Kernel (c = 0.0)"

        # Errors.
        @test_throws ArgumentError LinearKernel(; c=-0.5)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x -> LinearKernel(; c=x[1]), [c])
        test_params(LinearKernel(; c=c), ([c],))
    end
    @testset "PolynomialKernel" begin
        k = PolynomialKernel()
        @test kappa(k, x) ≈ x^2
        @test k(v1, v2) ≈ dot(v1, v2)^2
        @test kappa(PolynomialKernel(), x) == kappa(k, x)
        @test repr(k) == "Polynomial Kernel (c = 0.0, degree = 2)"

        # Coherence tests.
        @test kappa(PolynomialKernel(; degree=1, c=c), x) ≈ kappa(LinearKernel(; c=c), x)
        @test metric(PolynomialKernel()) == KernelFunctions.DotProduct()
        @test metric(PolynomialKernel(; degree=3)) == KernelFunctions.DotProduct()
        @test metric(PolynomialKernel(; degree=3, c=c)) == KernelFunctions.DotProduct()

        # Deprecations.
        k = @test_deprecated PolynomialKernel(; d=1)
        @test k.degree == 1

        # Errors.
        @test_throws ArgumentError PolynomialKernel(; d=0)
        @test_throws ArgumentError PolynomialKernel(; degree=0)
        @test_throws ArgumentError PolynomialKernel(; c=-0.5)
        @test_throws ErrorException PolynomialKernel(; d=2.5)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x -> PolynomialKernel(; c=x[1]), [c])
        test_params(PolynomialKernel(; c=c), ([c],))
    end
end


text:
"""
Void placeholder communicating that there is no redistribution happening
"""
function noneRedistribution(votingResult, challenger, deposit,dispensation_pct)

end

"""
Redistribution model when only the challenger receives the reward
in the form of the initial applicant deposit and the challenge fee
both being the same amount
"""
function onlyChallengerRewardRedistribution(result, challenger, deposit,dispensation_pct)
    if (!result[1])
        challenger.balance += 2 * deposit
    end
end

"""
Redistribution model when the challenger is rewarded and all who voted
against him are punished according to the dispensation parameter who defines severity
"""
function punishmentRedistribution(result, challenger, deposit, dispensation_pct)
    if (!result[1])
        challenger.balance += (2+dispensation_pct) * deposit
        foreach(agent -> agent.balance -= dispensation_pct*deposit / length(result[2]), result[2])
    end
end

"""
Redistribution model when all of voters who supported him are rewarded
"""
function rewardRedistribution(result, challenger, deposit,dispensation_pct)
    if (!result[1])
        challenger.balance += (1+dispensation_pct)*deposit
        foreach(agent -> agent.balance += (1-dispensation_pct)*deposit / length(result[3]), result[3])
    end
end

"""
Redistribution model when both the challenger is rewarded and voters who supported
him are rewarded and all who voted against him are punished
"""
function punishmentAndRewardRedistribution(result, challenger, deposit,dispensation_pct)
    if (!result[1])
        challenger.balance += 2* deposit
        foreach(agent -> agent.balance -= deposit / length(result[2]), result[2])
        foreach(agent -> agent.balance += deposit / length(result[3]), result[3])
    end
end

"""
Redistribution model when both the challenger is rewarded with his deposit and dispensation percantage and his supporters are rewarded with a deposit
and the remains of the dispensation.
"""

function punishmentAndRewardRedistributionDispensation(result, challenger, deposit,dispensation_pct)
    if (!result[1])
        challenger.balance += (1+dispensation_pct)* deposit
        foreach(agent -> agent.balance -= deposit / length(result[2]), result[2])
        foreach(agent -> agent.balance += (2-dispensation_pct)*deposit / length(result[3]), result[3])
    end
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
using AbstractTrees
import AbstractTrees: children
using CSTParser
using CSTParser: FunctionDef, BinarySyntaxOpCall, EXPR
using Tokenize: Tokens

struct OverlayNode{T}
    parent::Union{OverlayNode, Void}
    buffer::String
    expr::Union{EXPR{T}, T}
    fullspan::UnitRange{Int}
    span::UnitRange{Int}
end

Base.Expr(n::OverlayNode) = Base.Expr(n.expr)
AbstractTrees.parentlinks(::Type{<:OverlayNode}) = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{<:OverlayNode}) = AbstractTrees.ImplicitSiblings()
AbstractTrees.parent(node::OverlayNode) = node.parent
AbstractTrees.isroot(node::OverlayNode) = node.parent == nothing

OverlayNode(expr, buffer) = OverlayNode(nothing, buffer, expr, 1:(expr.fullspan), expr.span)
AbstractTrees.printnode(io::IO, o::OverlayNode{T}) where {T} = print(io, sprint(AbstractTrees.printnode, o.expr), " -- ", o.fullspan, " (", o.span, ")")
Base.show(io::IO, o::OverlayNode) = AbstractTrees.print_tree(io, o)

Base.length(o::OverlayNode) = length(children(o.expr))
Base.endof(o::OverlayNode)  = endof(children(o.expr))

Base.length(::CSTParser.LITERAL) = 0
Base.length(::CSTParser.OPERATOR) = 0
Base.length(::CSTParser.KEYWORD) = 0
Base.length(::CSTParser.PUNCTUATION) = 0
Base.length(::CSTParser.IDENTIFIER) = 0

children(node::OverlayNode) = node

# CSTParser gives argument in filter in the opposite lexical order,
# remove this when fixed in CSTParser
children(filter::EXPR{CSTParser.Filter}) = reverse(filter.args)

function Base.getindex(node::OverlayNode, idx::Integer)
    offset = first(node.fullspan)
    c = children(node.expr)
    for i = 1:idx-1
        offset += c[i].fullspan
    end
    expr = children(node.expr)[idx]
    OverlayNode(node, node.buffer, expr, offset:(offset+expr.fullspan-1), offset - 1 + expr.span)
end
function Base.getindex(node::OverlayNode, range::Range)
    map(x->node[x], range)
end

Base.setindex!(o::OverlayNode, x::OverlayNode, idx) = Base.setindex!(o.expr.args, x.expr, idx)

Base.start(node::OverlayNode) = (start(node.expr), first(node.fullspan))
function Base.next(node::OverlayNode, state::Tuple{Int, Int})
    idx, offset = state
    expr, idx = next(children(node.expr), idx)
    (OverlayNode(node, node.buffer, expr, offset:(offset+expr.fullspan-1), offset - 1 + expr.span), (idx, offset + expr.fullspan))
end
Base.done(node::OverlayNode, state::Tuple{Int, Int}) = done(node.expr, state[1])

function function_def_call(expr)
    if isexpr(expr, FunctionDef)
        expr = children(expr)[2]
    else
        @assert (isexpr(expr, BinarySyntaxOpCall) || isexpr(expr, CSTParser.WhereOpCall)) && isexpr(children(expr)[2], CSTParser.OPERATOR, Tokens.EQ)
    end
    while isexpr(expr, BinarySyntaxOpCall) || isexpr(expr, CSTParser.WhereOpCall)
        expr = children(expr)[1]
    end
    return expr
end


text:
using Test
using CompScienceMeshes

T = Float64
tol = eps(T) * 10^3
const third = one(T)/3

# test creation of a rectangle
m = meshrectangle(1.0,1.0,0.5);
@test numvertices(m) == 9
@test numcells(m) == 8

# test edge detection
edges = skeleton(m, 1)
@test numcells(edges) == 16

# test volume pair construction from a list of edges
vp = cellpairs(m,edges)
internal_edges = count(x->x>0, vp[2,:])
@test internal_edges == 8

# test the relative orientation routine
@test CompScienceMeshes.relorientation(index(1,2,3),index(1,2,3,4)) == -4
@test CompScienceMeshes.relorientation(index(2,3,4),index(1,2,3,4)) == +1
@test CompScienceMeshes.relorientation(index(1,3,2),index(1,2,3,4)) == +4
@test CompScienceMeshes.relorientation(index(2,4,3),index(1,2,3,4)) == -1

o, e1, e2, e3 = euclidianbasis(3)
cell = simplex(o, e1, e2)

A  = volume(cell)
@test abs(A - 0.5) < tol

mp = neighborhood(cell,[third, third]);
r = cartesian(mp)
@test norm(r - point(third, third, zero(T))) < tol

# repeat the test but now on a random cell
randpoint() = point(2*rand(T)-1, 2*rand(T)-1, 2*rand(T)-1)
v1 = randpoint()
v2 = randpoint()
v3 = randpoint()
cell = simplex(v1, v2, v3)

mp = neighborhood(cell,[third, third]);
r = cartesian(mp)

@test norm(third*(v1 + v2 + v3) - r) < tol
@test abs(dot(cross(r-v1, r-v2), r-v3)) < tol

# test that the normal is directed outwards
m = meshcircle(1.0, 2π/51)
#for i in 1:numcells(m)
for cl in cells(m)
    #p = chart(m, cells(m,i))
    v = vertices(m, cl)
    p = simplex(v)
    _c = (p.vertices[1] + p.vertices[2]) / 2
    @test dot(_c, p.normals[1]) > 0
end

## test cellpairs on non-oriented meshes
using CompScienceMeshes
using Test

p1 = point(1,0,0)
p2 = point(0,1,0)
p3 = point(-1,0,0)
p4 = point(0,-1,0)

i1 = index(1,2,3)
i2 = index(1,4,2)

m = Mesh([p1,p2,p3,p4], [i1,i2])
e = skeleton(m,1)

@test numcells(m) == 2
@test numcells(e) == 5

cps = cellpairs(m,e)


text:
using Documenter, Catwalk, Literate

Literate.markdown("docs/src/usage.jl", "docs/src"; documenter=false, execute=false)

makedocs(
    modules = [Catwalk],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Krisztián Schäffer",
    sitename = "Catwalk.jl",
    pages = Any["index.md", "usage.md", "howitworks.md", "tuning.md"],
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/tisztamo/Catwalk.jl.git",
    push_preview = true,
    devbranch = "main",
)


text:
mutable struct Agent
    x::Vector{Float64}; f::Float64; v::Bool; c::Float64
    # x := parameters passed into models
    # f := function-value of fn!(x)
    # v := viability / feasibility
    # c := contravention / violation

    Agent(ND::Int) = new(Vector{Float64}(undef, ND), Inf, false, 0.0) # @code_warntype ✓
end

function Base.isequal(a1::Agent, a2::Agent) # @code_warntype ✓
    # a1, a2 are both feasible
    a1.v && a2.v && return a1.f == a2.f
    # a1, a2 are both infesasible
    !a1.v && !a2.v && return a1.c == a2.c
    return false
end

function Base.isless(a1::Agent, a2::Agent) # @code_warntype ✓
    # a1, a2 are both feasible
    a1.v && a2.v && return a1.f < a2.f
    # a1, a2 are both infesasible
    !a1.v && !a2.v && return a1.c < a2.c
    # if (a1, a2) = (feasible, infeasible), then a1 < a2 is true
    # if (a1, a2) = (infeasible, feasible), then a2 < a1 is false
    return a1.v
end

#### random initialization, @code_warntype ✓
function born!(x::VecO, lb::NTuple{ND}, ub::NTuple{ND}) where ND
    @simd for i in eachindex(x)
        @inbounds x[i] = lb[i] + rand() * (ub[i] - lb[i])
    end
end

#### groups, @code_warntype ✓
function return_agents(ND::Int, NP::Int)
    agents = Vector{Agent}(undef, NP)
    @inbounds for i in eachindex(agents)
        agents[i] = Agent(ND)
    end
    return agents
end

#### subgroups, @code_warntype ✓
return_elites(agents::VecI{Agent}, NR::Int)          = view(agents, 1:NR)
return_throng(agents::VecI{Agent}, NR::Int, NP::Int) = view(agents, NR+1:NP)


text:
function hartree_pot_bc(rho::Array{Float64,2}, Ls::Float64, YukawaK::Float64, epsil0::Float64)
# Calculate the Hartree potential and the energy in periodic boundary potential
# this function is called from the one in the Ham. This is anotated to give the jit
# compiler more information

    Ns_glb = length(rho);
    kx = 2*pi* vcat(collect(0:Ns_glb/2),
                    collect(-Ns_glb/2+1:-1))./Ls ;
    # Use Yukawa rather than the bare Coulomb potential.  The method does
    # not depend critically on the fact that the Coulomb kernel is to be
    # used.
    invkx2 = [0;1 ./(kx[2:end].^2 .+ YukawaK^2)];
    # invkx2 = [1/YukawaK^2;1./(kx(2:end).^2+YukawaK^2)];
    # The 1/YukawaK^2 factor is not important for neutral defect calculation, but
    # might be important for the charged defect calculation.
    Vtemp = 4*pi*invkx2.*fft(rho)  # we can decorate this one
    V  = real(ifft(Vtemp));
 return V / epsil0;
end

############################################################################
##  Trying different optimizations

function hartree_pot_bc_opt(rho::Array{Float64,2}, Ls::Float64,
                            YukawaK::Float64, epsil0::Float64)
# Calculate the Hartree potential and the energy in periodic boundary potential
# this function is called from the one in the Ham. This is anotated to give the jit
# compiler more information
    rhoFourier = rfft(rho)
    inv_yukawa_fourier_mult!(rhoFourier,Ls,YukawaK)
    V  = irfft(rhoFourier, size(rho,1));
    return V / epsil0;
end

@inline function inv_yukawa_fourier_mult!(R::Array{ComplexF64,2}, Ls::Float64,
                                          YukawaK::Float64 )
    c = (2 * pi / Ls)^2
    c1 = 4*pi
    R[1] = 0
    @inbounds @simd for ii = 2:size(R)[1]
        R[ii] = c1*R[ii]/((ii-1)^2*c + YukawaK^2);
    end
end

# using 1D vector in general should produce faster resutls

function hartree_pot_bc_opt_vec(rho::Vector{Float64}, Ls::Float64,
                                YukawaK::Float64, epsil0::Float64)
# Calculate the Hartree potential and the energy in periodic boundary potential
# this function is called from the one in the Ham. This is anotated to give the jit
# compiler more information
    rhoFourier = rfft(rho)
    inv_yukawa_fourier_mult_vec!(rhoFourier,Ls,YukawaK)
    V  = irfft(rhoFourier, size(rho,1))/ epsil0;
    return V;
end

@inline function inv_yukawa_fourier_mult_vec!(R::Vector{ComplexF64}, Ls::Float64,
                                              YukawaK::Float64 )
    c = (2 * pi / Ls)^2
    c1 = 4*pi
    R[1] = 0
    @inbounds @simd for ii = 2:size(R)[1]
        R[ii] = c1*R[ii]/((ii-1)^2*c + YukawaK^2);
    end
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
module TestKalman

using QuantEcon
using Base.Test
using FactCheck

# set up
A = [.95 0; 0. .95]
Q = eye(2) .* 0.5
G = eye(2) .* .5
R = eye(2) .* 0.2
kf = Kalman(A, G, Q, R)

rough_kwargs = Dict(:atol => 1e-2, :rtol => 1e-4)

sig_inf, kal_gain = stationary_values(kf)

# Compute the Kalman gain and sigma infinity according to the
# recursive equations and compare
mat_inv = inv(G * sig_inf * G' + R)

kal_recursion = A * sig_inf * (G') * mat_inv
sig_recursion = A * sig_inf * A' - kal_recursion * G * sig_inf * A' + Q

facts("Testing kalman.jl") do
    # test stationary
    @fact sig_inf --> roughly(sig_recursion; rough_kwargs...)
    @fact kal_gain --> roughly(kal_recursion; rough_kwargs...)

    # test_update_using_stationary
    set_state!(kf, [0.0 0.0]', sig_inf)
    update!(kf, [0.0 0.0]')
    @fact kf.cur_sigma --> roughly(sig_inf; rough_kwargs...)
    @fact kf.cur_x_hat --> roughly([0.0 0.0]'; rough_kwargs...)

    # test update nonstationary
    curr_x, curr_sigma = ones(2, 1), eye(2) .* .75
    y_observed = fill(0.75, 2, 1)
    set_state!(kf, curr_x, curr_sigma)
    update!(kf, y_observed)
    mat_inv = inv(G * curr_sigma * G' + R)
    curr_k = A * curr_sigma * (G') * mat_inv
    new_sigma = A * curr_sigma * A' - curr_k * G * curr_sigma * A' + Q
    new_xhat = A * curr_x + curr_k * (y_observed - G * curr_x)
    @fact kf.cur_sigma --> roughly(new_sigma; rough_kwargs...)
    @fact kf.cur_x_hat --> roughly(new_xhat; rough_kwargs...)

end  # facts
end  # module


text:
####
# NOT FUNCTIONAL
#TODO: Convert this to v2 model - This will not work for v2
####

import ..GPU

import ..YOLO: loadWeights!

export loadWeights!

"""
    flipkernel(x)

Flips given kernel
"""
flipkernel(x) = x[end:-1:1, end:-1:1, :, :]

"""
    updateconv!(c, gamma, mean, variance)

Applies batch-normalization to given convolutional layer
"""
function updateconv!(c, gamma, mean, variance)
    gamma4 = reshape(gamma, 1, 1, 1, :)
    variance4 = reshape(variance, 1, 1, 1, :)
    gamma3 = reshape(gamma, 1, 1, :, 1)
    mean3 = reshape(mean, 1, 1, :, 1)
    variance3 = reshape(variance, 1, 1, :, 1)
    #update c.w
    a = convert(Float32, 0.001)
    c.w = (c.w .* gamma4) ./ sqrt.(variance4 .+ a)
    #update c.b
    c.b = c.b .- (gamma3 .* mean3 ./ sqrt.(variance3 .+ a))
    return c.w, c.b
end

"""
    loadWeights!(model::YOLO.v2_tiny.v2Chain, weights_filepath::String)

Loads layers' weights from given weights file path.
"""
function loadWeights!(model::v2Chain, settings::Settings)
    outlength = settings.cell_bboxes * (5 + settings.num_classes)

    open(settings.weights_filepath, "r") do io
        readconstants!(io)
        #First Conv layer
        loadconv!(model.layers[1], io, 3, 3, 3, 16)
        #Second Conv layer
        loadconv!(model.layers[3], io, 3, 3, 16, 32)
        #Third Conv layer
        loadconv!(model.layers[5], io, 3, 3, 32, 64)
        #4th Conv layer
        loadconv!(model.layers[7], io, 3, 3, 64, 128)
        #5th Conv layer
        loadconv!(model.layers[9], io, 3, 3, 128, 256)
        #6th Conv layer
        loadconv!(model.layers[11], io, 3, 3, 256, 512)
        #YoloPad
        model.layers[13].w[1, 1, 1, 1] = 1
        #7th Conv layer
        loadconv!(model.layers[14], io, 3, 3, 512, 1024)
        #8th Conv layer
        loadconv!(model.layers[15], io, 3, 3, 1024, 1024)
        #last layer
        read!(io, model.layers[16].b)
        toRead = Array{Float32}(undef, 1024 * outlength)
        read!(io, toRead)
        toRead = reshape(toRead, 1, 1, 1024, outlength)
        model.layers[16].w = permutedims(toRead, [2, 1, 3, 4])
        model.layers[16].w = flipkernel(model.layers[16].w)

        if GPU >= 0
            model.layers[16].w = KnetArray(model.layers[16].w)
            model.layers[16].b = KnetArray(model.layers[16].b)
            model.layers[13].w = KnetArray(model.layers[13].w)
        end
    end
end

"""
    loadconv!(c, io, d1, d2, d3, d4)

Loads the io to given convolutional layer and updates it by batch-normalization
"""
function loadconv!(c, io::IOStream, d1, d2, d3, d4)
    read!(io, c.b)
    gamma = Array{Float32}(undef, d4)
    mean = Array{Float32}(undef, d4)
    variance = Array{Float32}(undef, d4)
    read!(io, gamma)
    read!(io, mean)
    read!(io, variance)
    toRead = Array{Float32}(undef, d4 * d3 * d2 * d1)
    read!(io, toRead)
    toRead = reshape(toRead, d1, d2, d3, d4)
    c.w = permutedims(toRead, [2, 1, 3, 4])
    c.w, c.b = updateconv!(c, gamma, mean, variance)
    c.w = flipkernel(c.w)
    if GPU >= 0
        c.w = KnetArray(c.w)
        c.b = KnetArray(c.b)
    end
end

"""
    readconstants!(io)

Read constant and unnecessary numbers from the io stream
"""
function readconstants!(io::IOStream)
    major = read(io, Int32)
    minor = read(io, Int32)
    revision = read(io, Int32)
    iseen = read(io, Int32)
end