text:
###################################################
## data/synthetic.jl
##     Creates a synthetic school-bus routing problem
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

"""
    Stores information about one student's walking time to a corner-stop.
"""
struct WalkingToStop
    "Id of the corresponding corner stop"
    cornerStopId::Int
    "Walking distance, in meters"
    distance::Float64
end

"""
    Represents all the data needed for 1 student
"""
struct Student
    "Student id, correspond to the index in the array of students"
    id::Int
    "Unique student identifier, does not change if students are re-ordered"
    originalId::Int
    "Student's house position"
    position::Point
    "id of school the student is assigned to"
    school::Int
    "whether the student is door to door"
    isD2D::Bool
    "Max walking distance, in meters, 0 if d2d"
    maxwalking::Float64
end

"""
    Represents a potential bus-stop (not a door to door one)
"""
struct CornerStop
    "Stop id"
    id::Int
    "Stop location"
    position::Point
end

"""
    Creates a synthetic school bus problem.
"""
function syntheticproblem(;
    nschools::Int=10,
    nyards::Int=1,
    districtsize::Float64 = 5000.,
    dwelltime::Float64 = 15.*60.,
    studentsperschool::Int = 120,
    schoolrange=(1000., 8000.),
    studentspread::Float64=200.,
    studentclustering::Float64=0.85,
    maxwalking::Float64=1610*0.5,
    d2dfraction::Float64=0.
    )
    data = SchoolBusData()

    districtshape = syntheticdistrictshape(districtsize)

    data.schools = synthetic_schools(nschools, districtshape, dwelltime)
    data.yards = synthetic_yards(districtshape, nyards)
    students = synthetic_students(data.schools, districtshape, studentsperschool,
                                  schoolrange, studentspread, studentclustering,
                                  maxwalking, d2dfraction)
    cornerStops = synthetic_cornerstops(students, districtshape, maxwalking)
    data.stops = simplestopassignment(data, students, cornerStops,
                                      euclideanwalkingdistances(data, students, cornerStops))
    data.withBaseData = true

    data.params.bus_capacity                    = 65
    data.params.constant_stop_time              = 30.
    data.params.max_time_on_bus                 = 3600.
    data.params.stop_time_per_student           = 5.
    data.params.velocity                        = 35/3.6
    data.params.metric                          = EUCLIDEAN
    return data
end

"""
    Create the synthetic corner stops
    - Easy: generate one corner stop for each student, randomly within a radius around it.
"""
function synthetic_cornerstops(students, districtshape, maxwalking)
    stops = Vector{CornerStop}(length(students))

    for (i,student) in enumerate(students)
        centerposX, centerposY = (student.position.x, student.position.y)

        newpos = Point(Inf, Inf)
        while !isindistrict(districtshape, newpos)
            θ = 2π*rand()
            r = rand()+rand()
            r = maxwalking*min(r,2-r)

            newpos = Point(centerposX + r*cos(θ), centerposY + r*sin(θ))
        end
        position = newpos
        stops[i] = CornerStop(i, position)
    end
    return stops
end

"""
    Create synthetic students
    - So far all corner-stop students
"""
function synthetic_students(schools, districtshape, studentsperschool, schoolrange, studentspread, studentclustering, maxwalking, d2dfraction)
    studentsposition = synthetic_studentspositions(schools, districtshape, studentsperschool, schoolrange, studentspread, studentclustering)
    students = Student[]
    studentid = 1
    for (schoolid, positions) in enumerate(studentsposition), position in positions
        isD2D = rand() < d2dfraction
        push!(students,
        Student(studentid, studentid, position, schoolid, isD2D, maxwalking)
        )
        studentid += 1
    end
    students
end

"""
    Create the location of synthetic students and assign them to schoos
    - Schools have different range
    - students are "clustered"
"""
function synthetic_studentspositions(schools, districtshape, studentsperschool, schoolrange, studentspread, studentclustering)
    # total number of students is fixed
    totalstudents = studentsperschool * length(schools)
    # students are partitioned between schools
    schoolsplits = sort(shuffle(1:(totalstudents-1))[1:(length(schools)-1)])
    schoolsplits = vcat([0], schoolsplits, [totalstudents])
    studentsinschool = [schoolsplits[i+1] - schoolsplits[i] for i in 1:length(schools)]

    # school "range" i.e. typical student distance from school : random
    schoolrange = [schoolrange[1] + (schoolrange[2]-schoolrange[1])*rand() for i in 1:length(schools)]

    # generate the positions
    studentsposition = Vector{Point}[]

    function findfeasiblepositionaroundpoint(position, districtshape, range)
        centerposX, centerposY = position.x, position.y
        newpos = Point(Inf, Inf)
        while !isindistrict(districtshape, newpos)
            newpos = Point(centerposX +range * randn(), centerposY + range* randn())
        end
        return newpos
    end

    for school in schools
        push!(studentsposition, Point[])
        currentcluster = Point[]
        for i in 1:studentsinschool[school.id]
            if isempty(currentcluster)
                push!(currentcluster, findfeasiblepositionaroundpoint(school.position, districtshape, schoolrange[school.id]))
            else
                push!(currentcluster,
                      findfeasiblepositionaroundpoint(rand(currentcluster),
                                                      districtshape, studentspread))
            end
            if rand() > studentclustering
                append!(studentsposition[end],[point for point in currentcluster])
                currentcluster = Point[]
            end
        end
        append!(studentsposition[end],[point for point in currentcluster])
    end
    return studentsposition
end

"""
    Create synthetic yards and buses
    so far, very simple: a large number of full buses in each yard.
"""
function synthetic_yards(districtshape, nyards::Int)
    yards = Vector{Yard}(nyards)
    for i in 1:nyards
        yardposition = randomdistrictposition(districtshape)
        yards[i] = Yard(i, yardposition)
    end
    return yards
end

"""
    Create synthetic schools
"""
function synthetic_schools(nschools::Int, districtshape, dwelltime::Float64)
    # first generate the schools uniformly inside the district shape
    schoolpositions = [randomdistrictposition(districtshape) for i in 1:nschools]

    starttimes = randomstarttimes(nschools)

    schools = School[]
    minStartTime = 7.5 * 3600
    maxStartTime = 9.5 * 3600
    for i in 1:nschools
        push!(schools, School(
            i, i, schoolpositions[i], dwelltime,
            starttimes[i], minStartTime, maxStartTime
        ))
    end
    return schools
end

"""
    returns the shape of the district, a rotated cube centered on 0
"""
function syntheticdistrictshape(districtsize)
    return districtsize
end

"""
    return a uniformly random position in the district
"""
function randomdistrictposition(districtshape)
    x = districtshape*(rand()-0.5)
    y = districtshape*(rand()-0.5)
    return Point((x+y)/sqrt(2), (x-y)/sqrt(2))
end

"""
    returns true if the given position is in the district
"""
function isindistrict(districtshape, position)
    x,y = position.x, position.y
    return (-districtshape/2 <= (x+y)/sqrt(2) <= districtshape/2) && (-districtshape/2 <= (x-y)/sqrt(2) <= districtshape/2)
end

"""
    returns a random set of start times, given different schools
"""
function randomstarttimes(nschools::Int)
    return rand([7.5,8.5,9.5].*3600., nschools)
end

"""
    adds euclidean walking distance to a schoolbus problem
"""
function euclideanwalkingdistances(data::SchoolBusData, students::Vector{Student},
                                   cornerStops::Vector{CornerStop})
    walkingtostops = Vector{WalkingToStop}[WalkingToStop[] for s in students]

    for student in students
        if !student.isD2D
            for stop in cornerStops
                dist = euclideandistance(student.position, stop.position)
                if dist <= student.maxwalking
                    push!(walkingtostops[student.id], WalkingToStop(stop.id, dist))
                end
            end
        end
        sort!(walkingtostops[student.id], by=w->w.distance)
    end
    return walkingtostops
end

"""
    simple stop assignment (no limit on #students per stop)
"""
function simplestopassignment(data::SchoolBusData, students::Vector{Student},
                              cornerStops::Vector{CornerStop},
                              walkingtostops::Vector{Vector{WalkingToStop}};
                              λ::Real=1e4, args...)
    studentsinschool = [Int[] for school in data.schools]
    for student in students
        if !student.isD2D
            push!(studentsinschool[student.school], student.id)
        end
    end
    studentstops = zeros(Int, length(students))
    for school in data.schools
        studentids = studentsinschool[school.id]

        if !isempty(studentids)
            stops  = Set{Int}(stop.cornerStopId for i in studentids for stop in walkingtostops[i])

            model = Model(solver=GurobiSolver(Threads=getthreads(), OutputFlag=0;args...))
            @variable(model, z[i in studentids, stop in walkingtostops[i]], Bin)
            # student i choose stop j

            @variable(model, usedStop[stopId in stops], Bin) # stops used

            @constraint(model, oneStopPerStudent[i in studentids],
                        sum(z[i,stop] for stop in walkingtostops[i]) == 1)

            @constraint(model, stopIsUsed[i in studentids, stop in walkingtostops[i]],
                        z[i,stop] <= usedStop[stop.cornerStopId])

            @objective(model, Min,
                    sum(usedStop[stopid] for stopid in stops) +
                λ * sum(z[i,stop] * stop.distance for i in studentids, stop in walkingtostops[i])
            )

            status = solve(model)

            for i in studentids, stop in walkingtostops[i]
                if getvalue(z[i, stop]) >= 0.5
                    studentstops[i] = stop.cornerStopId
                end
            end
        end
    end
    selectedbusstopsandstudents = Dict{Tuple{Int,Int}, Vector{Int}}()
    # (school, cornerStop) => students
    for (studentid, stopid) in enumerate(studentstops)
        students[studentid].isD2D && continue
        schoolid = students[studentid].school
        if haskey(selectedbusstopsandstudents, (schoolid, stopid))
            push!(selectedbusstopsandstudents[schoolid, stopid], studentid)
        else
            selectedbusstopsandstudents[schoolid, stopid] = [studentid]
        end
    end
    stops = Vector{Stop}[Stop[] for i in eachindex(data.schools)]
    # Will be used to create the unique stop id.
    uniqueid = 1
    for (schoolid,stopid) in sort(collect(keys(selectedbusstopsandstudents)))
        stopstudents = selectedbusstopsandstudents[schoolid, stopid]
        push!(stops[schoolid],
            Stop(length(stops[schoolid])+1, uniqueid, stopid, schoolid,
                 cornerStops[stopid].position, length(stopstudents))
        )
        uniqueid += 1
    end
    for student in students
        if student.isD2D
            push!(stops[student.school],
                Stop(length(stops[student.school])+1, uniqueid, 0, student.school,
                     student.position, 1)
            )
            uniqueid += 1
        end
    end
    return stops
end

max_stars_repo_name: jihoyeo/SchoolBusRouting

max_issues_repo_name: jihoyeo/SchoolBusRouting

max_forks_repo_name: jihoyeo/SchoolBusRouting


text:
#+Tag L:target
# This is an macro example [[target][]]
macro swap(x,y)
   quote
      local tmp = $(esc(x))
      $(esc(x)) = $(esc(y))
      $(esc(y)) = tmp
    end
end

ext: jl


text:
type Racos
    rc::RacosCommon
    function Racos()
      return new(RacosCommon())
    end
end

# racos optimization function
function racos_opt!(racos::Racos, objective::Objective, parameter::Parameter; ub=1)
    rc = racos.rc
    rc_clear!(rc)
    rc.objective = objective
    rc.parameter = parameter
    init_attribute!(rc)
    t = parameter.budget / parameter.negative_size
    time_log1 = now()
    for i in 1:t
        j = 0
        iteration_num = length(rc.negative_data)
        while j < iteration_num
            if rand(rng, Float64) < rc.parameter.probability
                classifier = RacosClassification(rc.objective.dim, rc.positive_data,
                rc.negative_data, ub=ub)
                mixed_classification(classifier)
                solution, distinct_flag = distinct_sample_classifier(rc, classifier, data_num=rc.parameter.train_size)
            else
                solution, distinct_flag = distinct_sample(rc, rc.objective.dim)
            end
            #painc stop
            if isnull(solution)
                return rc.best_solution
            end
            # If the solution had been sampled, skip it
            if !distinct_flag
                continue
            end
            obj_eval(objective, solution)
            push!(rc.data, solution)
            j += 1
        end
        selection!(rc)
        rc.best_solution = rc.positive_data[1]
        # display expected running time
        if i == 4
            time_log2 = now()
            # second
            expected_time = t * (Dates.value(time_log2 - time_log1) / 1000) / 5
            if !isnull(rc.parameter.time_limit)
                expected_time = min(expected_time, rc.parameter.time_limit)
            end
            if expected_time > 5
                zoolog(string("expected remaining running time: ", convert_time(expected_time)))
            end
        end
        # time budget check
        if !isnull(rc.parameter.time_limit)
            if expected_time >= rc.parameter.time_limit
                zoolog("Exceed time limit")
                return rc.best_solution
            end
        end
        # terminal_value check
        if !isnull(rc.parameter.terminal_value)
            if rc.best_solution.value <= rc.parameter.terminal_value
                zoolog("terminal_value function value reached")
                return rc.best_solution
            end
        end
    end
    return rc.best_solution
end


text:
using TermWin

TermWin.initsession()
arr = readdir()
#s = "a test"
v = newTwPopup( rootTwScreen, arr )
v.title = "Input: "
activateTwObj( rootTwScreen )
ret = v.value
#unregisterTwObj( rootTwScreen, v )
TermWin.endsession()
println( "You chose ", string( ret ) )


text:
__precompile__(false)
module Ahorn

# Fixes theme issues on Windows
if Sys.iswindows()
    ENV["GTK_THEME"] = get(ENV, "GTK_THEME", "win32")
    ENV["GTK_CSD"] = get(ENV, "GTK_CSD", "0")
end

using Gtk, Gtk.ShortNames, Gtk.GConstants
using Cairo
using Maple
using Random

function abspath(path::String)
    return normpath(joinpath(@__DIR__, path))
end

macro abs_str(path)
    abspath(path)
end

include("helpers/gtk_helpers.jl")

# Force all widgets to be Left to Right
# Program is not really useable in Right to Left mode, which Gtk respects on RTL languages for some reason
setDefaultDirection(1)

const baseTitle = "Ahorn"
const iconFile = abs"../docs/logo-256.png"
const logoFile = abs"../docs/logo-1024.png"

const windowIcon = missing

const canvas = Canvas(0, 0)

storageDirectory = ""
configFilename = ""
persistenceFilename = ""

config = Dict{String, Any}()
persistence = Dict{String, Any}()

langdata = nothing

include("config.jl")
include("debug.jl")

include("loaded_state.jl")
include("camera.jl")

loadedState = LoadedState("", "")
camera = Camera(0, 0, 4)

# Time Machine workaround
sleep(0)

include("helpers/helpers.jl")
include("helpers/macros.jl")
include("lang.jl")
include("validation_entry.jl")
include("mods.jl")
include("hotkey.jl")
include("cursor.jl")
include("multi_range.jl")
include("color_constants.jl")
include("module_loader.jl")
include("shapes/rectangle.jl")
include("shapes/simple_curve.jl")
include("shapes/line.jl")
include("shapes/circle.jl")
include("shapes/ellipse.jl")
include("helpers/list_view_helper.jl")
include("helpers/ahorn_list_helper.jl")
include("helpers/color_helper.jl")
include("helpers/xna_colors.jl")
include("helpers/form_helper.jl")
include("helpers/processing_helper.jl")
include("helpers/entity_id_helper.jl")
include("windows/settings_window.jl")
include("windows/test_eval_window.jl")
include("materials_filter.jl")
include("helpers/menubar_helper.jl")
include("celeste_render.jl")
include("effects.jl")
include("windows/map_image_dumper_window.jl")
include("windows/sprite_dumper_window.jl")
include("selections.jl")
include("property_menu.jl")
include("event_handler.jl")
include("pre_save_sanitizers.jl")
include("windows/file_dialog_window.jl")
include("windows/room_window.jl")
include("windows/about_window.jl")
include("windows/map_window.jl")
include("windows/styleground_window.jl")
include("windows/metadata_window.jl")
include("windows/update_window.jl")
include("windows/exit_window.jl")
include("history_manager.jl")
include("roomlist.jl")
include("everest_rcon.jl")
include("hotkeys.jl")
include("tools.jl")
include("menubar.jl")
include("backups.jl")
include("loading_spinner.jl")
include("object_tile_names.jl")

# Time Machine workaround
sleep(0)

include("windows/main_window.jl")
include("interactive_workaround.jl")

include("init_ahorn.jl")
include("init_external_modules.jl")
include("init_signals.jl")
include("init_globals.jl")

include("display_main_window.jl")

end


text:
module Rebugger

using UUIDs
using REPL
import REPL.LineEdit, REPL.Terminals
using REPL.LineEdit: buffer, bufend, content, edit_splice!
using REPL.LineEdit: transition, terminal, mode, state

using Revise
using Revise: ExLike, RelocatableExpr, get_signature, funcdef_body, get_def, striplines!
using Revise: printf_maxsize
using HeaderREPLs

const msgs = []  # for debugging. The REPL-magic can sometimes overprint error messages

include("debug.jl")
include("ui.jl")
include("deepcopy.jl")

# Set up keys that enter rebug mode from the regular Julia REPL
# This should be called from your ~/.julia/config/startup.jl file
function repl_init(repl)
    repl.interface = REPL.setup_interface(repl; extra_repl_keymap = get_rebugger_modeswitch_dict())
end

function __init__()
    # Set up the Rebugger REPL mode with all of its key bindings
    repl_inited = isdefined(Base, :active_repl)
    @async begin
        while !isdefined(Base, :active_repl)
            sleep(0.05)
        end
        sleep(0.1) # for extra safety
        # Set up the custom "rebug" REPL
        main_repl = Base.active_repl
        repl = HeaderREPL(main_repl, RebugHeader())
        interface = REPL.setup_interface(repl; extra_repl_keymap=[get_rebugger_modeswitch_dict(), rebugger_keys])
        rebug_prompt_ref[] = interface.modes[end]
        add_keybindings(; override=repl_inited, deprecated_keybindings..., keybindings...)
    end
end

end # module


text:
#Demo 6: Bode Plot
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.EasyPlot
using CMDimData.MDDatasets
using Colors
import Printf: @sprintf

#==Attributes
===============================================================================#
dfltline = cons(:a, line=set(style=:solid, color=:blue, width=3))
markerline = cons(:a, line=set(style=:dash, width=2.5))
markerline_light = cons(:a, line=set(style=:dash, width=2.5, color=RGB24(.4,.4,.4)))
dfltglyph = cons(:a, glyph=set(shape=:square, size=3))

#==Input data
===============================================================================#
j = im
G = 20.0; 𝑓_p1 = 50e6; 𝑓_p2 = .9e9; 𝑓_z = Inf #.8e9
𝑓step = 5e6; 𝑓max = 10e9
𝑓min_disp = 10e6

#==Helper functions
===============================================================================#
phase(x) = rad2deg(angle(x))
function roots(a,b,c)
	r = sqrt(b*b-4*a*c)
	result = (-b .+ [r, -r]) / (2a)
	return (result[1], result[2])
end
xf_o2(f, G, f_z, f_p1, f_p2) = #2nd order transfer function
	G*(1 + j*(f/f_z)) / ( (1 + j*(f/f_p1)) * (1 + j*(f/f_p2)) )
function f_unitygain(f, G, f_z, f_p1, f_p2) #Compute unity gain frequency
	invsq(x) = (result = 1/x; return result*result)
	_z² = invsq(f_z); G² = G*G
	_p1² = invsq(f_p1); _p2² = invsq(f_p2)
	r = roots(_p1²*_p2², _p1²+_p2²-G²*_z², 1-G²)
	return sqrt(max(r...)) #Select sqrt(positive root)
end
f_3dB(f, G, f_z, f_p1, f_p2) = #Compute 3dB frequency
	f_unitygain(f, sqrt(2), f_z, f_p1, f_p2)

#==Calculations
===============================================================================#
𝑓 = DataF1(0:𝑓step:𝑓max)
X = xf_o2(𝑓, G, 𝑓_z, 𝑓_p1, 𝑓_p2)
𝑓0 = f_unitygain(𝑓, G, 𝑓_z, 𝑓_p1, 𝑓_p2)
𝑓BW = f_3dB(𝑓, G, 𝑓_z, 𝑓_p1, 𝑓_p2)
phase0 = phase(xf_o2(𝑓0, G, 𝑓_z, 𝑓_p1, 𝑓_p2))
pmargin = 180+phase0
@show pmargin

#==Generate EasyPlot
===============================================================================#
gridconf = set(vmajor = true, vminor = true, hmajor=false, hminor=false)
plot = cons(:plot, nstrips=2, legend=false, title = "Sample Bode Plot",
	xaxis = set(scale=:log, label="Frequency [Hz]", min=𝑓min_disp), #Avoid issues with log scale
	ystrip1 = set(scale=:dB20, axislabel="Magnitude [dB]", min=-10),
	ystrip2 = set(scale=:lin, axislabel="Phase [°]", min=-180),
	grid1 = gridconf, grid2 = gridconf
)

push!(plot,
	cons(:wfrm, abs(X), dfltline, label="|X|", strip=1),
	cons(:wfrm, phase(X), dfltline, label="∠X", strip=2),
)

𝑓3dBstr = @sprintf("𝑓3dB=%.1f MHz", 𝑓BW/1e6)
𝑓0str = @sprintf("𝑓0=%.1f MHz", 𝑓0/1e6)
pmstr = @sprintf("PM=%.1f°", pmargin)

push!(plot,
	cons(:atext, 𝑓3dBstr, x=𝑓BW, offset=set(x=-3), reloffset=set(y=.5), angle=-90, align=:bc, strip=1),
	cons(:atext, 𝑓0str, x=𝑓0, offset=set(x=-3), reloffset=set(y=.5), angle=-90, align=:bc, strip=1),
	cons(:atext, "0dB", y=0, reloffset=set(x=0.5), offset=set(y=3), align=:bc, strip=1),
	cons(:atext, pmstr, y=(phase0-180)/2, reloffset=set(x=0.5), align=:cc, strip=2),
)

#Add H/V markers to plot:
push!(plot,
	cons(:vmarker, 𝑓0, markerline, strip=0),
	cons(:vmarker, 𝑓BW, markerline, strip=0),
	cons(:hmarker, 0, markerline_light, strip=1),
	cons(:hmarker, phase0, markerline, strip=2),
)

pcoll = push!(cons(:plotcoll, title="Bode Test"), plot)

CMDimData.@includepkg EasyPlotInspect
plotgui = EasyPlot.displaygui(:InspectDR, pcoll)
ploth = 500; plotw = round(Int, ploth*1.6)
EasyPlot._write(:png, "demo6_test.png", plotgui, dim=set(w=plotw, h=ploth))

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
module TestPlanning

include("preamble.jl")

@testset begin
    A1 = ones(0)
    A2 = ones(0, 0)
    A3 = ones(0, 0, 0)
    A4 = ones(0, 0, 0, 0)
    @test plan(
        A1 => (:i,),
        A2 => (:i, :j),
    ).indices === _indices("ji")
    @test plan(
        A1 => (:i,),
        A2 => (:j, :i),
        A2 => (:j, :k),
    ).indices === _indices("ikj")
    @test plan(
        A1 => (:i,),
        A2' => (:j, :i),
        A2 => (:j, :k),
    ).indices === _indices("jki")
    @test plan(
        A1 => (:i,),
        transpose(A2) => (:j, :i),
        A2 => (:j, :k),
    ).indices === _indices("jki")
    @test plan(
        A1 => (:i,),
        A2' => (:j, :i),
        A2' => (:j, :k),
    ).indices === _indices("jik")
    @test plan(
        A1 => (:i,),
        A2' => (:i, :j),
        A2' => (:k, :j),
    ).indices === _indices("ikj")
    @test plan(
        A1 => (:i,),
        A2' => (:i, :j),
        A2 => (:i, :j),
        A3 => (:k, :i, :j),
        A4 => (:k, :j, :l, :m),
    ).indices === _indices("ijlmk")
end

end  # module


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function subpath_clipping(fname)
    width, height = 1000, 1000
    Drawing(width, height, fname)
    origin()
    background("ivory")
    sethue("black")
    image = readpng(dirname(@__FILE__) * "/stackoverflow.png")
    w = image.width
    h = image.height
    # all these graphics will make up the clipping mask:
    # reversepath makes holes in the path
    b = false
    for i in 3:15
        randomhue()
        poly(ngon(0, 0, 50 + (i * 30), i, pi/2, vertices=true), :path, reversepath=b)
        if i < 15
            newsubpath()
        end
        b  = !b
    end
    # clip them
    clippreserve()
    # draw image 'through' clipping mask
    scale(40, 40)
    placeimage(image, Point(-w/2, -h/2), 1)
    clipreset()
    @test finish() == true
end

fname = "subpath-clip-image.pdf"
subpath_clipping(fname)
println("...finished test: output in $(fname)")


text:
using LinearAlgebra

abstract type AbstractInterface end
export AbstractInterface

abstract type AbstractBoundaryInterface <: AbstractInterface end
export AbstractBoundaryInterface

abstract type AbstractInternalInterface <: AbstractInterface end
export AbstractInternalInterface

struct EmptyInterface <: AbstractInterface end
export EmptyInterface

@with_kw struct IdealGasCatalystInterface{T<:AbstractPhase,N<:AbstractPhase,Q<:AbstractReaction} <: AbstractInternalInterface
    gas::T
    catalyst::N
    reactions::Array{Q,1}
end
export IdealGasCatalystInterface

struct Inlet{Q<:Real,S,V<:AbstractArray,U<:Real,X<:Real} <: AbstractBoundaryInterface
    domain::S
    y::V
    F::Function
    T::U
    P::X
    H::Q
end

function Inlet(domain::V,conddict::Dict{String,X},F::Function) where {V,X<:Real,B<:Real}
    y = makespcsvector(domain.phase,conddict)
    T = conddict["T"]
    P = conddict["P"]
    yout = y./sum(y)
    H = dot(getEnthalpy.(getfield.(domain.phase.species,:thermo),T),yout)
    return Inlet(domain,yout,F,T,P,H)
end

export Inlet

struct Outlet{V} <: AbstractBoundaryInterface
    domain::V
    F::Function
end
export Outlet


text:
"""
    SeisFKFilter(in; <keyword arguments>)

Removes energy from seismic data by applying filters to a gather in the Frequency-wavenumber domain

# Arguments
- `in`: 2D input data in TX domain. First dimension is time.

# Keyword arguments
- `dt=0.004`: time sampling interval in secs.
- `dx=10`: space sampling interval in meters.
- `va=-2000,vb=-3000,vc=3000,vd=2000`: apparent velocity corners to filter

# Output
- `out`: Filtered data in TX domain

# Example
```julia
julia> d = SeisLinearEvents(); df = SeisFKFilter(d,va=-4000.,vb=-6000.,vc=6000.,vd=4000.);
```
*Credits: Aaron Stanton,2017*

"""
function SeisFKFilter(d::Matrix{Tv};dt::Tv=0.004,dx::Tv=10.0,va::Tv=-2000.,vb::Tv=-3000.,vc::Tv=3000.,vd::Tv=2000.) where { Tv <:AbstractFloat}

#MDS: This one does not work, revamp asap and add working example.
#Fan function needs to be clearly explain and output

    nt = size(d,1)
    nx = size(d,2)

    pa = 1. /va
    pb = 1. /vb
    pc = 1. /vc
    pd = 1. /vd
    m = fft(d)

    nf = size(m,1)
    dw = 1. /nf/dt
    nw = isodd(nf) ? round(Int,floor(nf/2)) + 1 : round(Int,floor(nf/2))
    nk = size(m,2)
    dk = 1. /nk/dx

    for iw=1:nw
	w = dw*(iw-1)
	for ik=1:nk
	    k = ik < round(Int,floor(nk/2)) + 1 ? dk*(ik-1) : -(dk*nk - dk*(ik-1))
	    p = k/w;
	    if p < pa
		m[iw,ik] *= 0.
	    elseif p >= pa && p < pb
		m[iw,ik] *= (p-pa)/(pb-pa)
	    elseif p >= pb && p <= pc
				m[iw,ik] *= 1.
	    elseif p > pc && p <= pd
		m[iw,ik] *= 1. - (p-pc)/(pd-pc)
	    elseif p > pd
				m[iw,ik] *= 0.
	    end
	end
    end
    # symmetries
    m = ifft(m,2)
    for iw=nw+1:nf
	m[iw,:] = conj(m[nf-iw+2,:])
    end
    d = real(ifft(m,1))

    return d[1:nt,1:nx];
end


text:
export
    WindowFilter,
    SyncPF,
    AsyncPF,
    CoherentPF

"""
    struct WindowFilter <: Filter

    General construction for a filter based on a flexible windowing scheme.

#arguments
    window_creator Defines the method used to create windows
    inference_algorithm Defines the algorithm to use on a window
    postprocess! A postprocessing function, that takes the runtime and does any additional processing needed to carry to the next iteration. Defaults to doing nothing.

"""
mutable struct WindowFilter <: Filter
    window_creator :: WindowCreator
    inference_algorithm :: InstantAlgorithm
    postprocess! :: Function
    latest_window :: Union{InstantRuntime, Nothing}
    WindowFilter(wc, ia, pp) = new(wc, ia, pp, nothing)
    WindowFilter(wc, ia) = new(wc, ia, run -> nothing, nothing)
end

function init_filter(wf::WindowFilter, dynrun::DynamicRuntime)
    ensure_all!(dynrun, current_time(dynrun))
    instrun = initial_instant_runtime(dynrun)
    # We assume no evidence or interventions at time 0
    # TODO: Handle placeholder beliefs
    infer(wf.inference_algorithm, instrun)
    wf.latest_window = instrun
    retrieve_values_from_instant_runtime!(dynrun, instrun)
    _store_beliefs(wf, dynrun, instrun)
end

function _store_beliefs(wf::WindowFilter, dynrun::DynamicRuntime{T}, instrun::InstantRuntime) where T
    dynnet = get_network(dynrun)
    for instvar in get_variables(get_network(instrun))
        instinst = current_instance(instrun, instvar)
        belief = marginal(wf.inference_algorithm, instrun, instinst)
        (dynname, t) = dynamic_name_and_time(instvar, T)
        dyninst = get_instance(dynrun, get_node(dynnet, dynname), t)
        set_value!(dynrun, dyninst, :belief, belief)
    end
end

function filter_step(wf::WindowFilter, dynrun::DynamicRuntime{T}, variables::Vector{Variable}, time::T, evidence::Dict{Symbol, Score}) where T
    dynnet = get_network(dynrun)
    insts = create_window(wf.window_creator, dynrun, variables, time)
    instrun = instant_runtime_from_instances(dynrun, insts)
    # Apply the dynamic evidence to the instant runtime
    instev = Dict{Symbol, Score}()
    for (name, sc) in evidence
        instev[instant_name(name, time)] = sc
    end
    # Apply beliefs in the dynamic network as placeholder beliefs in the instant network.
    placeholder_beliefs = Dict{Symbol,Dist}()
    inst_phs = get_placeholders(get_network(instrun))
    for instnode in inst_phs
        (dynname, t) = dynamic_name_and_time(instnode, T)
        dynnode = get_node(dynnet, dynname)
        dyninst = get_instance(dynrun, dynnode, t)
        belief = get_value(dynrun, dyninst, :belief)
        placeholder_beliefs[get_name(instnode)] = belief
    end
    # TODO: Handle interventions
    instinterv = Dict{Symbol,Dist}()
    infer(wf.inference_algorithm, instrun, instev, instinterv, placeholder_beliefs)
    wf.latest_window = instrun
    retrieve_values_from_instant_runtime!(dynrun, instrun)
    set_time!(dynrun, time)
    _store_beliefs(wf, dynrun, instrun)
end

function answer(::Marginal, ::WindowFilter, dynrun::Runtime, target::VariableInstance)
    return get_value(dynrun, target, :belief)
end

# function answer(query::Query, wf::WindowFilter, dynrun::Runtime, targets::Vector{VariableInstance})
#     # TODO: This code assumes that targets are in the latest_window, which might not be true
#     # for an asynchronous filter. We need to construct an instant window for the targets.
#     # For that, we need to put information from the dynamic window into the instant window,
#     # which is not done yet.
#     instrun = wf.latest_window
#     insttargets = VariableInstance[]
#     for target in targets
#         instname = instant_name(get_name(target), current_time(dynrun))
#         instnode = get_node(get_network(instrun), instname)
#         insttarget = current_instance(instrun, instnode)
#         push!(insttargets, insttarget)
#     end
#     answer(query, wf.inference_algorithm, wf.latest_window, insttargets)
# end


text:
##Test ABC algorithms and related functions.
##This is done by running a simple normal example.
##It is simply tested that these do not produce errors.
using ABCDistances
using Distributions

##Set up abcinput
function sample_sumstats(pars::Array)
    success = true
    stats = [pars[1] + 0.1*randn(1); randn(1)]
    (success, stats)
end

sobs = [0.0,0.0]

abcinput = ABCInput();
abcinput.prior = MvNormal(1, 100.0);
abcinput.sample_sumstats = sample_sumstats;
abcinput.abcdist = WeightedEuclidean(sobs);
abcinput.nsumstats = 2;

##ABC rejection
srand(1)
out1 = abcRejection(abcinput, 1000, 200)
abcRejection(abcinput, 1000, 0.1)

##ABC-PMC (doesn't quite test all options)
srand(1)
out2 = abcPMC3V(abcinput, 200, 1/2, 5000, silent=true)
abcPMC5(abcinput, 200, 1/2, 5000, silent=true)
abcPMC3V(abcinput, 200, 1/2, 5000, store_init=true, silent=true)
abcPMC3V(abcinput, 200, 1/2, 5000, store_init=true, silent=true)

abcPMC3(abcinput, 200, 1/2, 5000, silent=true)
abcPMC3(abcinput, 200, 1/2, 5000, store_init=true, silent=true)
abcPMC3(abcinput, 200, 1/2, 5000, diag_perturb=true, silent=true)

abcPMC4(abcinput, 200, 1/2, 5000, silent=true);
abcPMC4(abcinput, 200, 1/2, 5000, store_init=true, silent=true);

##Mean and variance functions
parameter_means(out1)
parameter_means(out2)
parameter_vars(out1)
parameter_vars(out2)


text:
using Pkg, DiffEqSensitivity, OrdinaryDiffEq, Zygote, ComponentArrays, ReverseDiff
using Parameters
Pkg.status("DiffEqSensitivity")
# DiffEqSensitivity v6.58.0

## Proof that ReverseDiff works with ComponentArrays
ReverseDiff.gradient(x->sum(x),ComponentArray(a=1.,b=7.,c=8.))
# ComponentVector{Float64}(a = 1.0, b = 1.0, c = 1.0)

## Proof that ReverseDiff for ODEs doesn't work with ComponentArrays
u0 = [1.0;1.0]
p = ComponentArray(a=rand(3),b=1.0,_c=rand(2))
function fiip(du,u,p,t)
  @unpack a, b, _c = p
  c,d = _c
  du[1] = dx = sum(a)*u[1] - b*u[1]*u[2]
  du[2] = dy = -c*u[2] + d*u[1]*u[2]
end

function foop(u,p,t)
  @unpack a, b, _c = p
  c,d = _c
  dx = sum(a)*u[1] - b*u[1]*u[2]
  dy = -c*u[2] + d*u[1]*u[2]
  [dx, dy]
end
prob = ODEProblem(foop,u0,(0.0,10.0),p)

prob2 = ODEProblem(fiip,u0,(0.0,10.0),p)

# ForwardDiff works
loss(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=ForwardDiffSensitivity()))
dp = Zygote.gradient(loss,p)[1]
# ComponentVector{Float64}(a = 7.349039781610272, b = -159.31079871982794, c = 74.93924771425637, d = -339.3272371527868)

# ReverseDiff fails
loss2(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=ReverseDiffAdjoint()))
dp2 = Zygote.gradient(loss2,p)[1]

### OOP
# some zeros
loss3(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
dp3 = Zygote.gradient(loss3,p)[1]

#  24.963 ms (200048 allocations: 10.48 MiB)
loss33(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))
@btime dp33 = Zygote.gradient(loss33,p)[1]

#  25.901 ms (199973 allocations: 10.48 MiB)
loss4(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint()))
@btime dp4 = Zygote.gradient(loss4,p)[1]

#  25.885 ms (204438 allocations: 10.60 MiB)
loss5(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=BacksolveAdjoint()))
@btime dp5 = Zygote.gradient(loss5,p)[1]

#  29.091 ms (217545 allocations: 11.31 MiB)
loss55(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP())))
@btime dp55 = Zygote.gradient(loss55,p)[1]

#  69.483 ms (565278 allocations: 29.98 MiB)
loss6(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=QuadratureAdjoint()))
@btime dp6 = Zygote.gradient(loss6,p)[1]

#  73.882 ms (573837 allocations: 30.47 MiB)
loss7(p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=QuadratureAdjoint(autojacvec=ZygoteVJP())))
@btime dp7 = Zygote.gradient(loss7,p)[1]

### IP

# some zeros
loss3(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
dp3 = Zygote.gradient(loss3,p)[1]

#  zeros
loss33(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))
@btime dp33 = Zygote.gradient(loss33,p)[1]

#  some zeros
loss4(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=InterpolatingAdjoint()))
@btime dp4 = Zygote.gradient(loss4,p)[1]

#  some zeros
loss5(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=BacksolveAdjoint()))
@btime dp5 = Zygote.gradient(loss5,p)[1]

#
loss6(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=QuadratureAdjoint()))
@btime dp6 = Zygote.gradient(loss6,p)[1]

loss7(p) = sum(solve(prob2,Tsit5(),u0=u0,p=p,saveat=0.1,sensealg=QuadratureAdjoint(autojacvec=ZygoteVJP())))
@btime dp7 = Zygote.gradient(loss7,p)[1]


text:
using OmegaCore
using OmegaTest
using Distributions
using Test

# function test_logpdf_ucond()
#   a = 1 ~ Normal(0, 1)
#   b = 2 ~ Normal(2, 3)
#   z(ω) = (a(ω), b(ω))
#   ω, ℓ = propose(Random.GLOBAL_RNG, z)
#   @test isinferred(propose, Random.GLOBAL_RNG, z)
#   a_, b_ = z(ω)
#   @test ℓ == logpdf(Normal(0, 1), a_) + logpdf(Normal(2, 3), b_)
# end

# # Test doesn't double count
# function test_logpdf_ucond_nodouble()
#   a = 1 ~ Normal(0, 1)
#   b = 2 ~ Normal(2, 3)
#   z(ω) = (a(ω), b(ω), a(ω))
#   ω, ℓ = propose(Random.GLOBAL_RNG, z)
#   @test isinferred(propose, Random.GLOBAL_RNG, z)
#   a_, b_ = z(ω)
#   @test ℓ == logpdf(Normal(0, 1), a_) + logpdf(Normal(2, 3), b_)
# end

# function test_logpdf_cond()
#   μ = 1 ~ Normal(0, 1)
#   x = 2 ~ Normalₚ(μ, 1)
#   μ_ = 0.1234
#   x_ = -0.54321
#   μc = μ |ᶜ x == x_
#   ω, ℓ = propose(Random.GLOBAL_RNG, μc)
# end

# @testset "Logpdf" begin
#   test_logpdf_ucond()
#   test_logpdf_ucond_nodouble()
# end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
@testset "gen_mask" begin
    p, q = synthetic_gradient(Prism(50))
    mask = gen_mask(p, q, 1)
    @test mask[50, 50, 1] == 0.0
    @test mask[1, 1, 1] == 1.0
end


text:
using SMLMSim
using Test

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


text:
module Seleroute

using JuMP
using LightGraphs, MetaGraphs
using Compose, Colors, GraphPlot, Cairo, Fontconfig

using Printf, Statistics, LinearAlgebra, SparseArrays, DelimitedFiles, Random
import Base: summary, getindex
import LightGraphs: edges, vertices
# import JuMP: copy_model

const MOF = MOI.FileFormats

const AbstractSimpleGraph = LightGraphs.AbstractSimpleGraph

inedges(g, v) = [Edge(iv, v) for iv in inneighbors(g, v)]
outedges(g, v) = [Edge(v, ov) for ov in outneighbors(g, v)]

# https://github.com/JuliaOpt/JuMP.jl/pull/1982
function write_LP(io, model::Model)
    lp_model = MOF.Model(format = MOF.FORMAT_LP)
    MOI.copy_to(lp_model, backend(model))
    MOI.write_to_file(lp_model, io)
end

function write_MOF(io, model::Model)
    mof_model = MOF.Model(format = MOF.FORMAT_MOF)
    MOI.copy_to(mof_model, backend(model))
    MOI.write_to_file(mof_model, io)
end

# https://github.com/JuliaOpt/JuMP.jl/issues/2019
function Base.getindex(reference_map::ReferenceMap, container::Union{Containers.DenseAxisArray, Containers.SparseAxisArray})
    return getindex.(reference_map, container)
end

include("modeltype.jl")
include("data/parameters.jl")
include("data/model.jl")
include("data/solution.jl")
include("data/topology.jl")
include("export.jl")
include("basemodels/helpers.jl")
include("basemodels/knowncapacities.jl")
include("basemodels/knowncapacities_flow.jl")
include("basemodels/knowncapacities_path.jl")
include("basemodels/objectives.jl")
include("compute/compute.jl")
include("compute/metrics_evaluate.jl")
include("compute/certain/certain.jl")
include("compute/certain/certain_cg.jl")
include("compute/certain/mmf.jl")
include("compute/oblivious/iter.jl")
include("compute/oblivious/iter_pathcg.jl")
include("compute/oblivious/rr.jl")
include("compute/oblivious/rr_flow.jl")
include("compute/oblivious/rr_path.jl")
include("compute/oblivious/rr_pathcg.jl") # Must be included after iter_pathcg.jl due to shared helpers.
include("output/plots.jl")
include("output/reports.jl")
include("experiments/random.jl")

# Export all symbols. Code copied from JuMP.
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    @eval export $sym
end

end


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

"""
    FIST

Fast Industrial-Strength Triangulation (FIST) of polygons.

This triangulation method is the method behind the famous Mapbox's
Earcut library. It is based on a ear clipping algorithm adapted
for complex n-gons with holes. It has O(n²) time complexity where
n is the number of vertices. In practice it is very efficient due
to heuristics implemented in the algorithm.

## References

* Held, M. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
* Eder et al. 2018. [Parallelized ear clipping for the triangulation and
  constrained Delaunay triangulation of polygons]
  (https://www.sciencedirect.com/science/article/pii/S092577211830004X)
"""
struct FIST <: DiscretizationMethod end

function discretize(polyarea::PolyArea, ::FIST)
  # build bridges in case the polygonal area has
  # holes, i.e. reduce to a single outer boundary
  𝒫 = polyarea |> unique |> bridge

  # points of resulting mesh
  points = vertices(𝒫)

  # keep track of indices of vertices
  inds = CircularVector(1:nvertices(𝒫))

  # perform ear clipping
  𝒬 = ears(𝒫)
  𝒯 = Connectivity{Triangle,3}[]
  clipped = false
  while nvertices(𝒫) > 3
    if !isempty(𝒬)
      i = pop!(𝒬)
      push!(𝒯, connect((inds[i-1], inds[i], inds[i+1]), Triangle))
      inds = [inds[begin:i-1]; inds[i+1:end]]
      𝒫 = Chain(points[inds])
      clipped = true
    elseif clipped
      𝒬 = ears(𝒫)
      clipped = false
    else
      # recovery process
      @warn "entered in recovery process"
    end
  end
  push!(𝒯, connect((inds[1], inds[2], inds[3]), Triangle))

  UnstructuredMesh(collect(points), 𝒯)
end

# return index of all ears of 𝒫
ears(𝒫) = filter(i -> isear(𝒫, i), 1:nvertices(𝒫))

# tells whether or not vertex i is an ear of 𝒫
# assuming that 𝒫 has counter-clockwise orientation
function isear(𝒫::Chain{Dim,T}, i) where {Dim,T}
  v = vertices(𝒫)

  # CE1.1: classify angle as convex vs. reflex
  α = ∠(v[i-1], v[i], v[i+1]) # oriented angle
  θ = α > 0 ? 2*T(π) - α : -α # inner angle
  isconvex = θ < π

  # CE1.2: check if segment vᵢ-₁ -- vᵢ+₁ intersects 𝒫
  sᵢ = Segment(v[i-1], v[i+1])
  intersects = false
  for j in 1:nvertices(𝒫)
    sⱼ = Segment(v[j], v[j+1])
    I = intersecttype(sᵢ, sⱼ)
    if !(I isa CornerTouchingSegments || I isa NonIntersectingSegments)
      intersects = true
      break
    end
  end

  # CE1.3: check if vᵢ-1 ∈ C(vᵢ, vᵢ+1, vᵢ+2) and vᵢ+1 ∈ C(vᵢ-2, vᵢ-1, vᵢ)
  c1 = sideof(v[i-1], Segment(v[i+1], v[i  ])) != :LEFT
  c2 = sideof(v[i-1], Segment(v[i+1], v[i+2])) != :RIGHT
  c3 = sideof(v[i+1], Segment(v[i-1], v[i-2])) != :LEFT
  c4 = sideof(v[i+1], Segment(v[i-1], v[i  ])) != :RIGHT
  incone = all((c1, c2, c3, c4))

  isconvex && !intersects && incone
end


text:
# Inspired from example1.cpp from SDPA
p = SDPA.SDPAProblem()
mDIM   = 3
nBlock = 1
SDPA.inputConstraintNumber(p, mDIM)
SDPA.inputBlockNumber(p, nBlock)
SDPA.inputBlockSize(p, 1, 2)
SDPA.inputBlockType(p, 1, SDPA.SDP)

SDPA.initializeUpperTriangleSpace(p)

SDPA.inputCVec(p, 1, 48)
SDPA.inputCVec(p, 2, -8)
SDPA.inputCVec(p, 3, 20)

SDPA.inputElement(p, 0, 1, 1, 1, -11., false)
SDPA.inputElement(p, 0, 1, 2, 2,  23., false)

SDPA.inputElement(p, 1, 1, 1, 1,  10., false)
SDPA.inputElement(p, 1, 1, 1, 2,   4., false)

SDPA.inputElement(p, 2, 1, 2, 2,  -8., false)

SDPA.inputElement(p, 3, 1, 1, 2,  -8., false)
SDPA.inputElement(p, 3, 1, 2, 2,  -2., false)

SDPA.initializeUpperTriangle(p, false)
SDPA.initializeSolve(p)

SDPA.solve(p)

@test SDPA.getIteration(p) == 10
@test isapprox(SDPA.getPrimalObj(p), -41.8999961638664)
@test isapprox(SDPA.getDualObj(p), -41.89999999999982)
@test SDPA.getPrimalError(p) < 1e-10
@test SDPA.getDualError(p) < 1e-10

SDPA.terminate(p)


text:
#
# types.jl -
#
# Type definitions for the Julia interface to FLI SDK
#
#------------------------------------------------------------------------------

# Julia type for a C enumeration.
const Cenum = Cint

const FLIDOMAIN_USB = 0x002
const FLIDEVICE_CAMERA = 0x100

struct CallError <: Exception
    code::Cint      # return value 0 success, non-zero on failure
    func::Symbol    # function causing error
end

abstract type CameraStatus end

# structs for acquisition signal
struct Stop<:CameraStatus end
struct Continue<:CameraStatus end

# Julia type alias
const flidev = Clong
const CameraHandle = Ptr{flidev}

"""

`TaoBindings.AbstractHighResolutionTime` is the parent type of time types with
a resolution of one nanosecond, that is [`TaoBindings.TimeSpec`](@ref) and
[`TaoBindings.HighResolutionTime`](@ref).

"""
abstract type AbstractHighResolutionTime end

"""
The structure `TaoBindings.HighResolutionTime` is the Julia equivalent to the
TAO `tao_time_t` structure.  Its members are `sec`, an integer number of
seconds, and `nsec`, an integer number of nanoseconds.

Also see [`TaoBindings.TimeSpec`](@ref).

"""
struct HighResolutionTime <: AbstractHighResolutionTime
    sec::Int64
    nsec::Int64
end

"""
    imgConfigContext
    stores configuration parameters for the hardware to be set prior to the
    next acquisition loop
""" ImageConfigContext

mutable struct ImageConfigContext

    width::Clong
    height::Clong
    offsetX::Clong
    offsetY::Clong
    exposuretime::Cdouble   # msec

    bytePerPixel::Int

    # TODO: add binning??

    function ImageConfigContext()
        max_width = 1600
        max_height = 1200
        return new(max_width, max_height, 0, 0,1,100.0)
    end
end

mutable struct CameraList
    camList::Vector{String}
    domainList::Vector{Clong}
    numCam::Int64

    function CameraList(interface::UInt16, deviceType::UInt16)

        inputDomain = convert(Clong, interface | deviceType)
        maxStrLength::UInt32 = 100

        filename = Vector{UInt8}(undef, 100)
        fill!(filename, 0x00)
        filename_ptr = pointer(filename)

        name  = Vector{UInt8}(undef, 100)
        fill!(name, 0x00)
        name_ptr =  pointer(name)

        domain = Ref{Clong}()

        numCam = 0
        camList = Vector{String}(undef,2)
        domainList = Vector{Clong}(undef,2)

        @checked_call(:FLICreateList,(Clong,), inputDomain)

        @checked_call(:FLIListFirst,(Ptr{Clong}, Ptr{UInt8},UInt32,
                    Ptr{UInt8} ,UInt32,), domain,filename_ptr, maxStrLength,
                    name_ptr, maxStrLength)

        numCam = numCam +1
        jfilename = strip(String(filename), Char(0x00))
        println(jfilename)
        camList[numCam] = jfilename
        domainList[numCam] = domain[]

        err = 0
        # Retrieve the camera list
        while true

            err = ccall((:FLIListNext, lib), Cint,
                        (Ptr{Clong}, Ptr{UInt8}, UInt32, Ptr{UInt8}, UInt32,),
                         domain,filename_ptr, maxStrLength,name_ptr, maxStrLength)

            if err != 0
                break
            end

            numCam = numCam +1
            jfilename = strip(String(filename), Char(0x00))
            println(jfilename)

            camList[numCam] = jfilename
            domainList[numCam] = domain[]

        end

        @checked_call(:FLIDeleteList,(Ptr{Cvoid},), Ref{Cvoid}())

        # Return the instanciated object.
        return new(camList, domainList, numCam)
    end
    nothing

end

mutable struct Camera
    handle::CameraHandle
    domain::Clong
    status::CameraStatus

    # get camera from indexing a camera list
    function Camera(devList::CameraList , i::Integer)
        1 ≤ i ≤ length(devList.numCam) || error("out of bound index in camera list ")
        devName = devList.camList[i]
        domainName = devList.domainList[i]
        ref = Ref{CameraHandle}(0)

        @checked_call(:FLIOpen, (Ptr{CameraHandle}, Cstring, Int64,),
                      ref, devName, domainName)
        return finalizer(_finalize,new(ref[], domainName))

    end
end


text:
export sample

using DataFrames: AbstractDataFrame, nrow
import StatsBase: sample
using StatsBase

function StatsBase.sample(df::AbstractDataFrame, args...; kwargs...)
    rows_sampled = sample(axes(df, 1), args...; kwargs...)
    df[rows_sampled, :]
end

function StatsBase.sample(rng, df::AbstractDataFrame, args...; kwargs...)
    rows_sampled = sample(rng, axes(df, 1), args...; kwargs...)
    df[rows_sampled, :]
end

function StatsBase.sample(df::AbstractDataFrame, frac::T; kwargs...) where T <: Union{AbstractFloat, Rational}
    @assert 0 <= frac <= 1
    n = round(Int, nrow(df)*frac)
    rows_sampled = sample(axes(df, 1), n, ; kwargs...)
    df[rows_sampled, :]
end

function StatsBase.sample(rng, df::AbstractDataFrame, frac::T; kwargs...) where T <: Union{AbstractFloat, Rational}
    @assert 0 <= frac <= 1
    n = round(Int, nrow(df)*frac)
    rows_sampled = sample(rng, axes(df, 1), n; kwargs...)
    df[rows_sampled, :]
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
@testset "Basic Figures" begin
    fig = Figure()
    @test current_figure() === fig

    fig2 = Figure()
    @test fig !== fig2
    @test current_figure() === fig2

    current_figure!(fig)
    @test current_figure() === fig
end

@testset "FigureAxisPlot" begin
    fap = scatter(rand(100, 2))
    @test fap isa Makie.FigureAxisPlot
    fig, ax, p = scatter(rand(100, 2))
    @test fig isa Figure
    @test ax isa Axis
    @test p isa Scatter

    fig2, ax2, p2 = scatter(rand(100, 3))
    @test fig2 isa Figure
    @test ax2 isa LScene # 3d plot
    @test p2 isa Scatter
end

@testset "AxisPlot and Axes" begin
    fig = Figure()
    @test current_axis() === nothing
    @test current_figure() === fig

    gridpos = fig[1, 1]
    @test gridpos isa GridPosition
    ap = scatter(gridpos, rand(100, 2))
    @test ap isa Makie.AxisPlot
    @test current_axis() === ap.axis

    ax2, p2 = scatter(fig[1, 2], rand(100, 2))
    @test ax2 isa Axis
    @test p2 isa Scatter
    @test current_axis() === ax2

    ax3, p3 = scatter(fig[1, 3], rand(100, 3))
    @test ax3 isa LScene
    @test p3 isa Scatter
    @test current_axis() === ax3

    @test ap.axis in fig.content
    @test ax2 in fig.content
    @test ax3 in fig.content

    current_axis!(fig, ax2)
    @test current_axis(fig) === ax2
    @test current_axis() === ax2

    fig2 = Figure()
    @test current_figure() === fig2
    @test current_axis() === nothing
    @test current_axis(fig) === ax2

    # current axis can also switch current figure when called without figure argument
    current_axis!(ax2)
    @test current_axis() === ax2
    @test current_figure() === fig
end

@testset "Deleting from figures" begin
    fig = Figure()
    ax = fig[1, 1] = Axis(fig)
    @test current_axis() === ax
    @test ax in fig.content
    delete!(ax)
    @test !(ax in fig.content)
    @test ax.parent === nothing
    @test current_axis() === nothing
end

@testset "Getting figure content" begin
    fig = Figure()
    ax = fig[1, 1] = Axis(fig)
    @test contents(fig[1, 1], exact = true) == [ax]
    @test contents(fig[1, 1], exact = false) == [ax]
    @test contents(fig[1:2, 1:2], exact = true) == []
    @test contents(fig[1:2, 1:2], exact = false) == [ax]

    @test content(fig[1, 1]) == ax
    @test_throws ErrorException content(fig[2, 2])
    @test_throws ErrorException content(fig[1:2, 1:2])

    label = fig[1, 1] = Label(fig)
    @test contents(fig[1, 1], exact = true) == [ax, label]
    @test contents(fig[1, 1], exact = false) == [ax, label]
    @test contents(fig[1:2, 1:2], exact = true) == []
    @test contents(fig[1:2, 1:2], exact = false) == [ax, label]

    @test_throws ErrorException content(fig[1, 1])

    ax2 = fig[1, 2][1, 1] = Axis(fig)
    @test contents(fig[1, 2][1, 1], exact = true) == [ax2]
    @test contents(fig[1, 2][1, 1], exact = false) == [ax2]
    @test contents(fig[1, 2][1:2, 1:2], exact = true) == []
    @test contents(fig[1, 2][1:2, 1:2], exact = false) == [ax2]

    label2 = fig[1, 2][1, 1] = Label(fig)
    @test contents(fig[1, 2][1, 1], exact = true) == [ax2, label2]
    @test contents(fig[1, 2][1, 1], exact = false) == [ax2, label2]
    @test contents(fig[1, 2][1:2, 1:2], exact = true) == []
    @test contents(fig[1, 2][1:2, 1:2], exact = false) == [ax2, label2]

    @test_throws ErrorException content(fig[1, 2][1, 1])
end

@testset "Nested axis assignment" begin
    fig = Figure()
    @test Axis(fig[1, 1]) isa Axis
    @test Axis(fig[1, 1][2, 3]) isa Axis
    @test Axis(fig[1, 1][2, 3][4, 5]) isa Axis
    @test_throws ErrorException scatter(fig[1, 1])
    @test_throws ErrorException scatter(fig[1, 1][2, 3])
    @test_throws ErrorException scatter(fig[1, 1][2, 3][4, 5])
    @test scatter(fig[1, 2], 1:10) isa Makie.AxisPlot
    @test scatter(fig[1, 1][1, 1], 1:10) isa Makie.AxisPlot
    @test scatter(fig[1, 1][1, 1][1, 1], 1:10) isa Makie.AxisPlot

    fig = Figure()
    fig[1, 1] = GridLayout()
    @test Axis(fig[1, 1][1, 1]) isa Axis
    fig[1, 1] = GridLayout()
    @test_throws ErrorException Axis(fig[1, 1][1, 1])
end

@testset "Not implemented error" begin
    @test_throws ErrorException("Not implemented for scatter. You might want to put:  `using Makie` into your code!") scatter()
end


text:
"""
$(SIGNATURES)

Convert a Franklin-Markdown file read as `mds` into a Franklin HTML string.
Returns the html string as well as a dictionary of page variables.

**Arguments**

* `mds`:         the markdown string to process
* `pre_lxdefs`:  a vector of `LxDef` that are already available.

**Keyword arguments**

* `isrecursive=false`: a bool indicating whether the call is the parent call or
                        a child call
* `isinternal=false`:  a bool indicating whether the call stems from `fd2html`
                        in internal mode
* `isconfig=false`:    a bool indicating whether the file to convert is the
                        configuration file
* `has_mddefs=true`:   a bool indicating whether to look for definitions of
                        page variables
"""
function convert_md(mds::AbstractString,
                    pre_lxdefs=collect(values(GLOBAL_LXDEFS));
                    isrecursive::Bool=false,
                    isinternal::Bool=false,
                    isconfig::Bool=false,
                    has_mddefs::Bool=true,
                    pagevar::Bool=false, # whether it's called from pagevar
                    )::String
    # instantiate page dictionaries
    isrecursive || isinternal || set_page_env()
    # if we're given a substring, force it to a string
    mds = String(mds)

    #
    # Parsing of the markdown string
    # (to find latex command, latex definitions, math envs etc.)
    #

    # ------------------------------------------------------------------------
    #> 1. Tokenize
    tokens = find_tokens(mds, MD_TOKENS, MD_1C_TOKENS)
    # distinguish fnref/fndef
    validate_footnotes!(tokens)
    # ignore header tokens that are not at the start of a line
    validate_headers!(tokens)
    # capture some hrule (see issue #432)
    hrules = find_hrules!(tokens)
    #> 1b. Find indented blocks (ONLY if not recursive to avoid ambiguities!)
    if !isrecursive
        find_indented_blocks!(tokens, mds)
        if has_mddefs
            # look for multi-line definitions, for those that meet the
            # requirements, discard *all* tokens in the span leaving only
            # the opening and closing one. Note: expressions will be parsed
            # twice as a result but really it's negligible time.
            preprocess_candidate_mddefs!(tokens)
        end
        # any remaining indented lines that don't meet the requirements for
        # an indented code block get removed here
        filter_lr_indent!(tokens, mds)
    end

    # ------------------------------------------------------------------------
    #> 2. Open-Close blocks (OCBlocks)
    #>> a. find them
    blocks, tokens  = find_all_ocblocks(tokens, MD_OCB)
    toks_pre_ocb    = deepcopy(tokens) # see find_double_brace_blocks
    #>> a'. find the rest
    blocks2, tokens = find_all_ocblocks(tokens, vcat(MD_OCB2, MD_OCB_MATH))
    append!(blocks, blocks2)
    ranges = deactivate_inner_blocks!(blocks)
    #>> b. merge CODE_BLOCK_IND which are separated by emptyness
    merge_indented_blocks!(blocks, mds)
    #>> b'. only keep indented code blocks which are not contained in larger
    # blocks (#285)
    filter_indented_blocks!(blocks)
    #>> c. now that blocks have been found, line-returns can be dropped
    filter!(τ -> τ.name ∉ L_RETURNS, tokens)
    #>> e. keep track of literal content of possible link definitions to use
    validate_and_store_link_defs!(blocks)

    if globvar("autocode") && any(b -> b.name in CODE_BLOCKS_NAMES, blocks)
        set_var!(LOCAL_VARS, "hascode", true)
    end
    if globvar("automath") && any(b -> b.name in MATH_BLOCKS_NAMES, blocks)
        set_var!(LOCAL_VARS, "hasmath", true)
    end

    # ------------------------------------------------------------------------
    #> 3. LaTeX commands
    #>> a. find "newcommands", update active blocks/braces
    lxdefs, tokens, braces, blocks = find_lxdefs(tokens, blocks)
    #>> b. if any lxdefs are given in the context, merge them. `pastdef` specifies
    # that the definitions appeared "earlier"
    lprelx = length(pre_lxdefs)
    (lprelx > 0) && (lxdefs = cat(pastdef.(pre_lxdefs), lxdefs, dims=1))
    #>> c. find latex commands
    lxcoms, _ = find_lxcoms(tokens, lxdefs, braces)

    #> 3[ex]. find double brace blocks, note we do it on pre_ocb tokens
    # as the step `find_all_ocblocks` possibly found and deactivated {...}.
    dbb = find_double_brace_blocks(toks_pre_ocb)
    deactivate_inner_dbb!(dbb, ranges)

    # ------------------------------------------------------------------------
    #> 4. Page variable definition (mddefs), also if in config, update lxdefs
    if has_mddefs
        process_mddefs(blocks, isconfig, pagevar)
    end

    #> 4.b if config, update global lxdefs as well
    if isconfig
        for lxd ∈ lxdefs
            GLOBAL_LXDEFS[lxd.name] = pastdef(lxd)
        end
        # if it's the config file we don't need to go further
        return ""
    end

    # ------------------------------------------------------------------------
    #> 5. Process special characters and html entities so that they can be
    # injected as they are in the HTML later
    sp_chars = find_special_chars(tokens)

    # ========================================================================
    #
    # Forming of the html string
    #
    # filter out the fnrefs that are left (still active)
    # and add them to the blocks to insert
    fnrefs = filter(τ -> τ.name == :FOOTNOTE_REF, tokens)

    # Discard indented blocks unless locvar("indented_code")
    if !locvar("indented_code")
        filter!(b -> b.name != :CODE_BLOCK_IND, blocks)
    end

    #> 1. Merge all the blocks that will need further processing before
    # insertion
    b2insert = merge_blocks(lxcoms, deactivate_divs(vcat(blocks, dbb)),
                            sp_chars, fnrefs, hrules)

    #> 2. Form intermediate markdown + html
    inter_md, mblocks = form_inter_md(mds, b2insert, lxdefs)
    inter_html = md2html(inter_md; stripp=isrecursive)

    #> 3. Plug resolved blocks in partial html to form the final html
    hstring = convert_inter_html(inter_html, mblocks, lxdefs)

    # final var adjustment, infer title if not given
    if isnothing(locvar("title")) && !isempty(PAGE_HEADERS)
        title = first(values(PAGE_HEADERS))[1]
        set_var!(LOCAL_VARS, "title", title)
        ALL_PAGE_VARS[splitext(locvar("fd_rpath"))[1]]["title"] =
            deepcopy(LOCAL_VARS["title"])
    end

    # Return the string
    return hstring
end

"""
$(SIGNATURES)

Same as `convert_md` except tailored for conversion of the inside of a math
block (no command definitions, restricted tokenisation to latex tokens). The
offset keeps track of where the math block was, which is useful to check
whether any of the latex command used in the block have not yet been defined.

**Arguments**

* `ms`:     the string to convert
* `lxdefs`: existing latex definitions prior to the math block
* `offset`: where the mathblock is with respect to the parent string
"""
function convert_md_math(ms::AS, lxdefs::Vector{LxDef}=Vector{LxDef}(),
                         offset::Int=0)::String
    # if a substring is given, copy it as string
    ms = String(ms)
    #
    # Parsing of the markdown string
    # (to find latex command, latex definitions, math envs etc.)
    #

    #> 1. Tokenize (with restricted set)
    tokens = find_tokens(ms, MD_TOKENS_LX, MD_1C_TOKENS_LX)

    #> 2. Find braces and drop line returns thereafter
    blocks, tokens = find_all_ocblocks(tokens, MD_OCB_ALL, inmath=true)
    braces = filter(β -> β.name == :LXB, blocks)

    #> 3. Find latex commands (indicate we're in a math environment + offset)
    lxcoms, _ = find_lxcoms(tokens, lxdefs,  braces, offset; inmath=true)

    #
    # Forming of the html string
    # (see `form_inter_md`, it's similar but simplified since there are fewer
    # conditions)
    #
    htmls = IOBuffer()

    strlen   = lastindex(ms)
    len_lxc  = length(lxcoms)
    next_lxc = iszero(len_lxc) ? BIG_INT : from(lxcoms[1])

    # counters to keep track of where we are and which command we're looking at
    head, lxc_idx = 1, 1
    while (next_lxc < BIG_INT) && (head < strlen)
        # add anything that may occur before the first command
        (head < next_lxc) && write(htmls, subs(ms, head, prevind(ms, next_lxc)))
        # add the first command after resolving, bool to indicate that we're in
        # a math env
        write(htmls, resolve_lxcom(lxcoms[lxc_idx], lxdefs, inmath=true))
        # move the head to after the lxcom and increment the com counter
        head     = nextind(ms, to(lxcoms[lxc_idx]))
        lxc_idx += 1
        next_lxc = from_ifsmaller(lxcoms, lxc_idx, len_lxc)
    end
    # add anything after the last command
    (head <= strlen) && write(htmls, subs(ms, head, strlen))
    return String(take!(htmls))
end

"""
INSERT

String that is plugged as a placeholder of blocks that need further processing.
The spaces allow to handle overzealous inclusion of `<p>...</p>` from the base
Markdown to HTML conversion.
"""
const INSERT     = " ##FDINSERT## "
const INSERT_    = strip(INSERT)
const INSERT_PAT = Regex(INSERT_)
const INSERT_LEN = length(INSERT_)

"""
$SIGNATURES

Form an intermediate MD file where special blocks are replaced by a marker
(`INSERT`) indicating that a piece will need to be plugged in there later.

**Arguments**

* `mds`:    the (sub)string to convert
* `blocks`: vector of blocks
* `lxdefs`: existing latex definitions prior to the math block
"""
function form_inter_md(mds::AS, blocks::Vector{<:AbstractBlock},
                       lxdefs::Vector{LxDef}
                       )::Tuple{String, Vector{AbstractBlock}}
    strlen  = lastindex(mds)
    intermd = IOBuffer()
    # keep track of the matching blocks for each insert
    mblocks = Vector{AbstractBlock}()

    len_b   = length(blocks)
    len_lxd = length(lxdefs)

    # check when the next block is
    next_b = iszero(len_b) ? BIG_INT : from(blocks[1])

    # check when the next lxblock is, extra work because there may be lxdefs
    # passed through in *context* (i.e. that do not appear in mds) therefore
    # search first lxdef actually in mds (nothing if lxdefs is empty)
    first_lxd = findfirst(δ -> (from(δ) > 0), lxdefs)
    next_lxd  = isnothing(first_lxd) ? BIG_INT : from(lxdefs[first_lxd])

    # check what's next: a block or a lxdef
    b_or_lxd = (next_b < next_lxd)
    nxtidx = min(next_b, next_lxd)

    # keep track of a few counters (where we are, which block, which command)
    head, b_idx, lxd_idx = 1, 1, first_lxd

    while (nxtidx < BIG_INT) & (head < strlen)
        # check if there's anything before head and next block and write it
        (head < nxtidx) && write(intermd, subs(mds, head, prevind(mds, nxtidx)))
        # check whether it's a block first or a newcommand first
        if b_or_lxd # it's a block, check if should be written
            β = blocks[b_idx]
            # check whether the block should be skipped
            if isa(β, OCBlock) && β.name ∈ MD_OCB_IGNORE
                head = nextind(mds, to(β))
            else
                if isa(β, OCBlock) && β.name ∈ MD_HEADER
                    # this is a trick to allow whatever follows the title to be
                    # properly parsed by Markdown.parse; it could otherwise
                    # cause issues for instance if a table starts immediately
                    # after the title
                    write(intermd, INSERT * "\n ")
                else
                    write(intermd, INSERT)
                end
                push!(mblocks, β)
                head = nextind(mds, to(blocks[b_idx]))
            end
            b_idx += 1
            next_b = from_ifsmaller(blocks, b_idx, len_b)
        else
            # newcommand or ignore --> skip, increase counters, move head
            head     = nextind(mds, to(lxdefs[lxd_idx]))
            lxd_idx += 1
            next_lxd = from_ifsmaller(lxdefs, lxd_idx, len_lxd)
        end
        # check which block is next
        b_or_lxd = (next_b < next_lxd)
        nxtidx = min(next_b, next_lxd)
    end
    # add whatever is after the last block
    (head <= strlen) && write(intermd, subs(mds, head, strlen))

    # combine everything and return
    return String(take!(intermd)), mblocks
end

"""
$(SIGNATURES)

Take a partial markdown string with the `INSERT` marker and plug in the
appropriately processed block.

**Arguments**

* `ihtml`:  the intermediary html string (with `INSERT`)
* `blocks`: vector of blocks
* `lxdefs`: latex context
"""
function convert_inter_html(ihtml::AS,
                            blocks::Vector{<:AbstractBlock},
                            lxdefs::Vector{LxDef})::String
    # Find the INSERT indicators
    allmatches = collect(eachmatch(INSERT_PAT, ihtml))
    strlen = lastindex(ihtml)

    # write the pieces of the final html in order, gradually processing the
    # blocks to insert
    htmls = IOBuffer()
    head  = 1
    for (i, m) ∈ enumerate(allmatches)
        # two cases can happen based on whitespaces around an insertion that
        # we want to get rid of, potentially both happen simultaneously.
        # 1.
##FDINSERT##...
        # 2. ...##FDINSERT##

        # exceptions,
        # - list items introduce

 and
\n
 which shouldn't remove
        # - end of doc introduces
(\n?) which should not be removed
        δ1, δ2 = 0, 0 # keep track of the offset at the front / back
        # => case 1
        c10 = prevind(ihtml, m.offset, 7) # *li>

        c1a = prevind(ihtml, m.offset, 3) # *p>
        c1b = prevind(ihtml, m.offset)    #
 0) && ihtml[c10:c1b] == "

"
        !(hasli1) && (c1a > 0) && ihtml[c1a:c1b] == "
" && (δ1 = 3)

        # => case 2
        iend = m.offset + INSERT_LEN
        c2a  = nextind(ihtml, iend)
        c2b  = nextind(ihtml, iend, 4)  # </p*
        c20  = nextind(ihtml, iend, 10) #
\n 0 && ihtml[prev] == ' '
            prev = prevind(ihtml, prev)
        end
        (head ≤ prev) && write(htmls, subs(ihtml, head:prev))
        # move head appropriately
        head = iend + δ2
        if head ≤ strlen
            head = ifelse(ihtml[head] in (' ', '>'), nextind(ihtml, head), head)
        end
        # store the resolved block
        resolved = convert_block(blocks[i], lxdefs)
        write(htmls, resolved)
    end
    # store whatever is after the last INSERT if anything
    (head ≤ strlen) && write(htmls, subs(ihtml, head:strlen))
    # return the full string
    return String(take!(htmls))
end


text:
module TestIntrospection

using CUDA
using FoldsCUDA: FailedInference
using Test

function test_failedinference_regular()
    err = FailedInference(nothing, (1,), Union{}, (1,), Union{})
    msg = sprint(showerror, err)
    @test occursin("not return", msg)
    @test !occursin("on the host,", msg)
end

function test_failedinference_ok_on_host()
    err = FailedInference(nothing, (1,), Union{}, (1,), Int)
    msg = sprint(showerror, err)
    @test occursin("not return", msg)
    @test occursin("on the host,", msg)
end

function test_failedinference_invalid_type()
    err = FailedInference(nothing, (1,), Int, (1,), Int)
    msg = sprint(showerror, err)
    @test occursin("invalid type", msg)
end

function test_failedinference_code_typed()
    CUDA.functional() || return
    err = FailedInference(identity, (1,), Int, (1,), Int)
    @test CUDA.code_typed(err) isa Any
end

end  # module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

@given("put attribute type: license, with value type: string") do context
    _put_attribute_to_db(context, "license", VALUE_TYPE.STRING)
end

@given("put attribute type: date, with value type: datetime") do context
    _put_attribute_to_db(context, "date", VALUE_TYPE.DATETIME)
end

@given("put relation type: marriage") do context
    res = g.put(ConceptManager(context[:transaction]), RelationType, "marriage")
end

@given("relation(marriage) set relates role: wife") do context
    res = g.set_relates(context[:transaction],
                    g.Label("","marriage"),
                    "wife")
end

@given("relation(marriage) set relates role: husband") do context
    g.set_relates(context[:transaction],
                    g.Label("","marriage"),
                    "husband")
end

@given("relation(marriage) set owns key type: license") do context
    rel_type = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    lic_typ = g.get(ConceptManager(context[:transaction]), AttributeType, "license")
    set_owns(g.as_remote(rel_type, context[:transaction]),
                        lic_typ,
                        true)
end

@given("relation(marriage) set owns attribute type: date") do context
    rel_type = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    lic_typ = g.get(ConceptManager(context[:transaction]), AttributeType, "date")
    set_owns(g.as_remote(rel_type, context[:transaction]),
                        lic_typ,
                        false)
end

@given("entity(person) set plays role: marriage:wife") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("marriage", "wife"),false))
end

@given("entity(person) set plays role: marriage:husband") do context
    ent_pers = g.get(ConceptManager(context[:transaction]), EntityType, "person")
    g.set_plays(g.as_remote(ent_pers, context[:transaction]), g.RoleType(g.Label("marriage", "husband"),false))
end

@when("\$m = relation(marriage) create new instance with key(license): m") do context
    attr = g.get(context[:cm], AttributeType, "license")
    rel_license = g.put(g.as_remote(attr, context[:transaction]), "m")

    mar_type = g.get(context[:cm], RelationType, "marriage")
    context[:m] = g.create(g.as_remote(mar_type, context[:transaction]))
    g.set_has(context[:transaction], context[:m], rel_license)
end

@then("relation \$m is null: false") do context
    @expect context[:m] !== nothing
end

@then("relation \$m has type: marriage") do context
    @expect context[:m].type.label.name == "marriage"
end

@then("relation(marriage) get instances contain: \$m") do context
    rel_type = g.get(ConceptManager(context[:transaction]), RelationType, "marriage")
    res_inst = g.get_instances(g.as_remote(rel_type, context[:transaction]))

    @expect in(context[:m], res_inst)
end

@when("\$b = entity(person) create new instance with key(username): bob") do context
    attr_user = g.get(context[:cm], AttributeType, "username")
    user_bob = g.put(g.as_remote(attr_user, context[:transaction]), "bob")

    type_pers = g.get(context[:cm], EntityType, "person")
    person = g.create(g.as_remote(type_pers, context[:transaction]))

    g.set_has(context[:transaction], person, user_bob)
    context[:b] = person
end

@when("relation \$m add player for role(wife): \$a") do context
    g.add_player(context[:transaction],
                context[:m],
                g.RoleType(g.Label("marriage","wife"),false),
                context[:a])
end

@when("relation \$m add player for role(husband): \$b") do context
    g.add_player(context[:transaction],
                context[:m],
                g.RoleType(g.Label(context[:m].type.label.name,"husband"),false),
                context[:b])
end

@then("relation \$m get players for role(wife) contain: \$a") do context
    roles = [g.RoleType(g.Label(context[:m].type.label.name, "wife"),false)]
    res = g.get_players(context[:transaction],
                context[:m],
                roles)
    @expect in(context[:a], res)
end

@then("relation \$m get players for role(husband) contain: \$b") do context
    roles = [g.RoleType(g.Label("marriage","husband"),false)]
    res = g.get_players(context[:transaction],
                context[:m],
                roles)
    @expect in(context[:b], res)
end

@then("relation \$m get players contain: \$a") do context
    res = g.get_players(context[:transaction], context[:m])
    @expect in(context[:a], res)
end

@then("relation \$m get players contain: \$b") do context
    res = g.get_players(context[:transaction], context[:m])
    @expect in(context[:b], res)
end

@then("relation \$m get players contain:") do context
    roles_inp = [db[1] for db in context.datatable]
    entities = [db[2] for db in context.datatable]

    context[Symbol("a")] = context[:a]

    for i in 1:length(roles_inp)
        roles = [g.RoleType(g.Label("marriage", roles_inp[i]),false)]
        res = g.get_players(context[:transaction],
                    context[:m],
                    roles)
        @expect in(context[Symbol(string(SubString(entities[i],2,2)))], res)
    end
end

@when("\$m = relation(marriage) get instance with key(license): m") do context
    attr_license  = g.get(context[:cm], AttributeType, "license")
    license = g.get(as_remote(attr_license, context[:transaction]),"m")
    res_owner = g.get_owners(context[:transaction], license)
    context[:m] = isempty(res_owner) ? nothing : res_owner[1]
end

@when("\$b = entity(person) get instance with key(username): bob") do context
    attr  = g.get(context[:cm], AttributeType, "username")
    obj_attr = g.get(as_remote(attr, context[:transaction]),"bob")
    res_owner = g.get_owners(context[:transaction], obj_attr)[1]
    context[:b] = res_owner
end

# Scenario: Role players can get relations
@then("entity \$a get relations(marriage:wife) do not contain: \$m") do context
    rel_alice = g.get_relations(context[:transaction], context[:a])
    @expect !in(context[:m], rel_alice)
end

@then("entity \$b get relations(marriage:husband) do not contain: \$m") do context
    rel_bob= g.get_relations(context[:transaction], context[:b])
    @expect !in(context[:m], rel_bob)
end

@then("entity \$a get relations(marriage:wife) contain: \$m") do context
    rel_alice = g.get_relations(context[:transaction],
                                context[:a],
                                [g.RoleType(g.Label("marriage","wife"),false)])
    @expect in(context[:m], rel_alice)
end

@then("entity \$b get relations(marriage:husband) contain: \$m") do context
    rel_bob = g.get_relations(context[:transaction],
                                context[:b],
                                [g.RoleType(g.Label("marriage","husband"),false)])
    @expect in(context[:m], rel_bob)
end

# Scenario: Role player can be unassigned from relation
@when("relation \$m remove player for role(wife): \$a") do context
    rem_req = g.remove_player(context[:transaction],
                        context[:m],
                        RoleType(g.Label("marriage","wife"),false),
                        context[:a])
end

@then("relation \$m get players for role(wife) do not contain: \$a") do context
    res = g.get_players(context[:transaction], context[:m])
    @expect !in(context[:a], res)
end

@then("relation \$m get players do not contain:") do context
    roles_inp = [db[1] for db in context.datatable]
    entities = [db[2] for db in context.datatable]

    context[Symbol("a")] = context[:a]

    for i in 1:length(roles_inp)
        roles = [g.RoleType(g.Label("marriage", roles_inp[i]),false)]
        res = g.get_players(context[:transaction],
                    context[:m],
                    roles)
        @expect !in(context[Symbol(string(SubString(entities[i],2,2)))], res)
    end
end

# Scenario: Relation without role players get deleted
@then("relation \$m is deleted: true") do context
    res = g.is_deleted(g.as_remote(context[:m], context[:transaction]))
    @expect res
end

@then("entity \$a get relations do not contain: \$m") do  context
    rel_alice = g.get_relations(context[:transaction],
                                context[:a])
    @expect !in(context[:m], rel_alice)
end

@then("relation(marriage) get instances do not contain: \$m") do context
    type_marriage = g.get(g.ConceptManager(context[:transaction]),
                            g.RelationType,
                            "marriage")

    res_marriage = g.get_instances(g.as_remote(type_marriage, context[:transaction]))
    @expect !in(context[:m], res_marriage)
end

@then("relation \$m is null: true") do context
    @expect context[:m] === nothing
end

@then("relation(marriage) get instances is empty") do context
    type_marriage = g.get(g.ConceptManager(context[:transaction]),
    g.RelationType,
    "marriage")
    res_marriage = g.get_instances(g.as_remote(type_marriage, context[:transaction]))

    @expect isempty(res_marriage)
end

# Scenario: Relation with role players can be deleted
@when("delete relation: \$m") do context
    g.delete(g.as_remote(context[:m], context[:transaction]))
end

@then("entity \$b get relations do not contain: \$m") do context
   relations = g.get_relations(context[:transaction], context[:b])
   @expect !in(context[:m], relations)
end

# Scenario: Relation cannot have roleplayers inserted after deletion
@then("relation \$m add player for role(wife): \$a; throws exception") do context
    try
        g.add_player(context[:transaction],
                context[:m],
                g.RoleType(g.Label("marriage","wife"),false),
                context[:res_entity][1].data["x"])
    catch ex
        @expect ex !== nothing
    end
end


text:
using LinearAlgebra
using BlockFactorizations
using BenchmarkTools

d = 512
n = 16
# testing Diagonal BlockFactorization
A = Diagonal([(randn(d, d)) for i in 1:n])
B = BlockFactorization(A)

# x = [randn(d) for _ in 1:n]
# y = [zeros(d) for _ in 1:n]
# @btime A*x; # this doesn't work in 1.6
# @btime mul!(y, A, x);

x = randn(d*n)
y = zeros(d*n)
@btime B*x;
@btime mul!(y, B, x);


text:
module SignalBroadcastingSim

using Test
using Random
using OpenStreetMapX
using LightGraphs
using DataFrames
using Distributions
using HypothesisTests
using SparseArrays
using Serialization
using Parameters
using LinearAlgebra
using DelimitedFiles
using Statistics
using StatsBase
using PyCall
using Serialization
using Distributed

export Simulation
export ModelParams
export Agent
export init!
export step!

export enu_weighted
export VertexMeta
export plot_sim
export simulate_all!
export simulate_steps!

include("agents.jl")
include("visuals.jl")
include("run.jl")

end


text:
module Jags

using Compat, Pkg, Documenter, DelimitedFiles, Unicode, MCMCChains, StatsPlots

#### Includes ####

include("jagsmodel.jl")
include("jagscode.jl")

if !isdefined(Main, :Stanmodel)
  include("utilities.jl")
end

"""The directory which contains the executable `bin/stanc`. Inferred
from `Main.JAGS_HOME` or `ENV["JAGS_HOME"]` when available. Use
`set_jags_home!` to modify."""
JAGS_HOME=""

function __init__()
    global JAGS_HOME = if isdefined(Main, :JAGS_HOME)
        eval(Main, :JAGS_HOME)
    elseif haskey(ENV, "JAGS_HOME")
        ENV["JAGS_HOME"]
    else
        println("Environment variable JAGS_HOME not found. Use set_jags_home!.")
        ""
    end
end

"""Set the path for `Jags`.

Example: `set_jags_home!(homedir() * "/src/src/cmdstan-2.11.0/")`
"""
set_jags_home!(path) = global JAGS_HOME=path

#### Exports ####

export
# From this file
  set_jags_home!,
# From Jags.jl
  JAGS_HOME,

# From jagsmodel.jl
  Jagsmodel,

# From jagscode.jl
  jags

end # module


text:
# Autogenerated wrapper script for libdrm_jll for armv7l-linux-musleabihf
export libdrm

using Xorg_libpciaccess_jll
JLLWrappers.@generate_wrapper_header("libdrm")
JLLWrappers.@declare_library_product(libdrm, "libdrm.so.2")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libpciaccess_jll)
    JLLWrappers.@init_library_product(
        libdrm,
        "lib/libdrm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const GBVecOrMat{T} = Union{GBVector{T}, GBMatrix{T}}
const GBMatOrTranspose{T} = Union{GBMatrix{T}, Transpose{T, GBMatrix{T}}}
const GBVecOrTranspose{T} = Union{GBVector{T}, Transpose{T, GBVector{T}}}
const GBArray{T} = Union{GBVecOrTranspose{T}, GBMatOrTranspose{T}}
const ptrtogbtype = Dict{Ptr, AbstractGBType}()

const GrBOp = Union{
    libgb.GrB_Monoid,
    libgb.GrB_UnaryOp,
    libgb.GrB_Semiring,
    libgb.GrB_BinaryOp,
    libgb.GxB_SelectOp
}

const TypedOp = Union{
    TypedUnaryOperator,
    TypedBinaryOperator,
    TypedMonoid,
    TypedSemiring
}

const MonoidBinaryOrRig = Union{
    TypedMonoid,
    TypedSemiring,
    TypedBinaryOperator,
    AbstractSemiring,
    AbstractBinaryOp,
    AbstractMonoid
}

const OperatorUnion = Union{
    AbstractOp,
    GrBOp
}


text:
export agent2string

##########################################################################################
# ABMPlot recipe
##########################################################################################

"Define ABMPlot plotting function with some attribute defaults."
@recipe(ABMPlot, agent_pos, model) do scene
    Theme(
        # insert InteractiveDynamics theme here?
    )
    Attributes(
        ac = JULIADYNAMICS_COLORS[1],
        as = 10,
        am = :circle,
        scatterkwargs = NamedTuple(),
    )
end

# 2D space
function Makie.plot!(abmplot::ABMPlot{<:Tuple{Vector{Point2f0}, <:Agents.ABM}})
    scatter!(abmplot, abmplot[:agent_pos];
        color=abmplot[:ac], marker=abmplot[:am], markersize=abmplot[:as],
        abmplot[:scatterkwargs]...
    )

    return abmplot
end

# 3D space
function Makie.plot!(abmplot::ABMPlot{<:Tuple{Vector{Point3f0}, <:Agents.ABM}})
    abmplot.am[] == :circle && (abmplot.am = Sphere(Point3f0(0), 1))

    meshscatter!(abmplot, abmplot[:agent_pos];
        color=abmplot[:ac], marker=abmplot[:am], markersize=abmplot[:as],
        abmplot[:scatterkwargs]...
    )

    return abmplot
end

# 2D polygons
function Makie.plot!(abmplot::ABMPlot{<:Tuple{Vector{<:Polygon{2}}, <:Agents.ABM}})
    poly!(abmplot, abmplot[:agent_pos];
        color=abmplot[:ac],
        abmplot[:scatterkwargs]...
    )

    return abmplot
end

##########################################################################################
# Agent inspection on mouse hover
##########################################################################################

# 2D space
function Makie.show_data(inspector::DataInspector,
            plot::ABMPlot{<:Tuple{Vector{Point2f0}, <:Agents.ABM}},
            idx, ::Scatter)
    a = inspector.plot.attributes
    scene = Makie.parent_scene(plot)

    proj_pos = Makie.shift_project(scene, plot, to_ndim(Point3f0, plot[1][][idx], 0))
    Makie.update_tooltip_alignment!(inspector, proj_pos)
    as = plot.as[]

    cursor_pos = (plot[1][][idx].data[1], plot[1][][idx].data[2])
    s = typeof(plot.model[].space)
    if s <: Agents.ContinuousSpace
        cursor_pos = Float64.(cursor_pos)
    elseif s <: Agents.GridSpace
        cursor_pos = Int.(cursor_pos)
    end
    a._display_text[] = agent2string(plot.model[], cursor_pos)
    a._bbox2D[] = FRect2D(proj_pos .- 0.5 .* as .- Vec2f0(5), Vec2f0(as) .+ Vec2f0(10))
    a._px_bbox_visible[] = true
    a._bbox_visible[] = false
    a._visible[] = true

    return true
end

# 3D space
function Makie.show_data(inspector::DataInspector,
            plot::ABMPlot{<:Tuple{Vector{Point3f0}, <:Agents.ABM}},
            idx, ::MeshScatter)
    a = inspector.plot.attributes
    scene = Makie.parent_scene(plot)

    proj_pos = Makie.shift_project(scene, plot, to_ndim(Point3f0, plot[1][][idx], 0))
    Makie.update_tooltip_alignment!(inspector, proj_pos)
    as = plot.as[]

    cursor_pos = (plot[1][][idx].data[1], plot[1][][idx].data[2], plot[1][][idx].data[3])
    s = typeof(plot.model[].space)
    if s <: Agents.ContinuousSpace
        cursor_pos = Float64.(cursor_pos)
    elseif s <: Agents.GridSpace
        cursor_pos = Int.(cursor_pos)
    end
    a._display_text[] = agent2string(plot.model[], cursor_pos)
    a._bbox2D[] = FRect2D(proj_pos .- 0.5 .* as .- Vec2f0(5), Vec2f0(as) .+ Vec2f0(10))
    a._px_bbox_visible[] = true
    a._bbox_visible[] = false
    a._visible[] = true

    return true
end

# 2D polygons
# TODO: Fix this tooltip
function Makie.show_data(inspector::DataInspector,
            plot::ABMPlot{<:Tuple{Vector{<:Polygon{2}}, <:Agents.ABM}},
            idx, ::Makie.Poly)
    a = inspector.plot.attributes
    scene = Makie.parent_scene(plot)

    proj_pos = Makie.shift_project(scene, plot, to_ndim(Point3f0, plot[1][][idx], 0))
    Makie.update_tooltip_alignment!(inspector, proj_pos)
    as = plot.as[]

    cursor_pos = (plot[1][][idx].data[1], plot[1][][idx].data[2])
    s = typeof(plot.model[].space)
    if s <: Agents.ContinuousSpace
        cursor_pos = Float64.(cursor_pos)
    elseif s <: Agents.GridSpace
        cursor_pos = Int.(cursor_pos)
    end
    a._display_text[] = agent2string(plot.model[], cursor_pos)
    a._bbox2D[] = FRect2D(proj_pos .- 0.5 .* as .- Vec2f0(5), Vec2f0(as) .+ Vec2f0(10))
    a._px_bbox_visible[] = true
    a._bbox_visible[] = false
    a._visible[] = true

    return true
end

DiscretePos = Union{NTuple{2, Int}, NTuple{3, Int}}
ContinuousPos = Union{NTuple{2, Float64}, NTuple{3, Float64}}

function agent2string(model::Agents.ABM, cursor_pos::DiscretePos)
    ids = Agents.ids_in_position(cursor_pos, model)
    s = ""

    for id in ids
        s *= agent2string(model[id]) * "\n"
    end

    return s
end

function agent2string(model::Agents.ABM, cursor_pos::ContinuousPos)
    ids = Agents.nearby_ids(cursor_pos, model, 0.01)
    s = ""

    for id in ids
        s *= agent2string(model[id]) * "\n"
    end

    return s
end

"""
    agent2string(agent::A)
Convert agent data into a string which is used to display all agent variables and their
values in the tooltip on mouse hover. Concatenates strings if there are multiple agents
at one position.

Custom tooltips for agents can be implemented by adding a specialised method
for `agent2string`.

Example:

```julia
function InteractiveDynamics.agent2string(agent::SpecialAgent)
    \"\"\"
    ✨ SpecialAgent ✨
    ID = \$(agent.id)
    Main weapon = \$(agent.charisma)
    Side weapon = \$(agent.pistol)
    \"\"\"
end
```
"""
function agent2string(agent::A) where {A<:Agents.AbstractAgent}
    agentstring = "▶ $(nameof(A))\n"

    agentstring *= "id: $(getproperty(agent, :id))\n"

    agent_pos = getproperty(agent, :pos)
    typeof(agent_pos) <: ContinuousPos && (agent_pos = round.(agent_pos, digits=2))
    agentstring *= "pos: $(agent_pos)\n"

    for field in fieldnames(A)[3:end]
        val = getproperty(agent, field)
        V = typeof(val)
        if V <: AbstractFloat
            val = round(val, digits=2)
        elseif V <: AbstractArray{<:AbstractFloat}
            val = round.(val, digits=2)
        elseif V <: Tuple && V <: NTuple{length(val), <:AbstractFloat}
            val = round.(val, digits=2)
        end
        agentstring *= "$(field): $val\n"
    end

    return agentstring
end


text:
package ln

import (
	"math"
	"math/rand"
)

struct Sphere
    Center::Vector
    Radius::Float64
    Box::Box
end

type Sphere struct {
	Center Vector
	Radius float64
	Box    Box
}

func NewSphere(center Vector, radius float64) *Sphere {
	min := Vector{center.X - radius, center.Y - radius, center.Z - radius}
	max := Vector{center.X + radius, center.Y + radius, center.Z + radius}
	box := Box{min, max}
	return &Sphere{center, radius, box}
}

func (s *Sphere) Compile() {
}

func (s *Sphere) BoundingBox() Box {
	return s.Box
}

func (s *Sphere) Contains(v Vector, f float64) bool {
	return v.Sub(s.Center).Length() <= s.Radius+f
}

func (s *Sphere) Intersect(r Ray) Hit {
	radius := s.Radius
	to := r.Origin.Sub(s.Center)
	b := to.Dot(r.Direction)
	c := to.Dot(to) - radius*radius
	d := b*b - c
	if d > 0 {
		d = math.Sqrt(d)
		t1 := -b - d
		if t1 > 1e-2 {
			return Hit{s, t1}
		}
		t2 := -b + d
		if t2 > 1e-2 {
			return Hit{s, t2}
		}
	}
	return NoHit
}

func (s *Sphere) Paths4() Paths {
	var paths Paths
	var seen []Vector
	var radii []float64
	for i := 0; i < 140; i++ {
		var v Vector
		var m float64
		for {
			v = RandomUnitVector()
			m = rand.Float64()*0.25 + 0.05
			ok := true
			for i, other := range seen {
				threshold := m + radii[i] + 0.02
				if other.Sub(v).Length() < threshold {
					ok = false
					break
				}
			}
			if ok {
				seen = append(seen, v)
				radii = append(radii, m)
				break
			}
		}
		p := v.Cross(RandomUnitVector()).Normalize()
		q := p.Cross(v).Normalize()
		n := rand.Intn(4) + 1
		for k := 0; k < n; k++ {
			var path Path
			for j := 0; j <= 360; j += 5 {
				a := Radians(float64(j))
				x := v
				x = x.Add(p.MulScalar(math.Cos(a) * m))
				x = x.Add(q.MulScalar(math.Sin(a) * m))
				x = x.Normalize()
				x = x.MulScalar(s.Radius).Add(s.Center)
				path = append(path, x)
			}
			paths = append(paths, path)
			m *= 0.75
		}
	}
	return paths
}

func (s *Sphere) Paths3() Paths {
	var paths Paths
	for i := 0; i < 20000; i++ {
		v := RandomUnitVector()
		v = v.MulScalar(s.Radius).Add(s.Center)
		paths = append(paths, Path{v, v})
	}
	return paths
}

func (s *Sphere) Paths2() Paths {
	var equator Path
	for lng := 0; lng <= 360; lng++ {
		v := LatLngToXYZ(0, float64(lng), s.Radius)
		equator = append(equator, v)
	}
	var paths Paths
	for i := 0; i < 100; i++ {
		m := Identity()
		for j := 0; j < 3; j++ {
			v := RandomUnitVector()
			m = m.Rotate(v, rand.Float64()*2*math.Pi)
		}
		m = m.Translate(s.Center)
		paths = append(paths, equator.Transform(m))
	}
	return paths
}

func (s *Sphere) Paths() Paths {
	var paths Paths
	n := 10
	o := 10
	for lat := -90 + o; lat <= 90-o; lat += n {
		var path Path
		for lng := 0; lng <= 360; lng++ {
			v := LatLngToXYZ(float64(lat), float64(lng), s.Radius).Add(s.Center)
			path = append(path, v)
		}
		paths = append(paths, path)
	}
	for lng := 0; lng <= 360; lng += n {
		var path Path
		for lat := -90 + o; lat <= 90-o; lat++ {
			v := LatLngToXYZ(float64(lat), float64(lng), s.Radius).Add(s.Center)
			path = append(path, v)
		}
		paths = append(paths, path)
	}
	return paths
}

func LatLngToXYZ(lat, lng, radius float64) Vector {
	lat, lng = Radians(lat), Radians(lng)
	x := radius * math.Cos(lat) * math.Cos(lng)
	y := radius * math.Cos(lat) * math.Sin(lng)
	z := radius * math.Sin(lat)
	return Vector{x, y, z}
}

type OutlineSphere struct {
	Sphere
	Eye Vector
	Up  Vector
}

func NewOutlineSphere(eye, up, center Vector, radius float64) *OutlineSphere {
	sphere := NewSphere(center, radius)
	return &OutlineSphere{*sphere, eye, up}
}

func (s *OutlineSphere) Paths() Paths {
	var path Path
	center := s.Sphere.Center
	radius := s.Sphere.Radius

	hyp := center.Sub(s.Eye).Length()
	opp := radius
	theta := math.Asin(opp / hyp)
	adj := opp / math.Tan(theta)
	d := math.Cos(theta) * adj
	r := math.Sin(theta) * adj

	w := center.Sub(s.Eye).Normalize()
	u := w.Cross(s.Up).Normalize()
	v := w.Cross(u).Normalize()
	c := s.Eye.Add(w.MulScalar(d))
	for i := 0; i <= 360; i++ {
		a := Radians(float64(i))
		p := c
		p = p.Add(u.MulScalar(math.Cos(a) * r))
		p = p.Add(v.MulScalar(math.Sin(a) * r))
		path = append(path, p)
	}
	return Paths{path}
}


text:
using Cones
using Test

println("Running test_utils.jl...")
include("test_utils.jl")

println("Running test_FiniteCone.jl...")
include("test_FiniteCone.jl")

println("Running test_chernikova.jl...")
include("test_chernikova.jl")

println("Running test_chernikova_general.jl...")
include("test_chernikova_general.jl")

println("Running test_PolyhedralCone.jl...")
include("test_PolyhedralCone.jl")

println("Running test_polytope.jl...")
include("test_polytope.jl")


text:
function holtslag(Tv,q,U,V,LH,SH,ustar,ZF)
# a vertical diffusion scheme based on Holtslag and Boville (1993)
# Includes enhenced diffusion below the PBL

#TF-turbulent fluxes, positive upward
#Tv-virtual potential temperature

Rcrit=0.5;
gravity_mks=9.81;
karman=0.4;
a=7.8;
eps=0.1;
rhoa=1.2;
cpa=1005; #J/K/kg
Av=2500000;
lamdac=300;
humid_fac=0.606;

Z=ZF[2:end];
DZ=diff(ZF,dims=1);
ZM=Z-DZ./2;
SZM=length(ZM);

ZD=Z[1]:0.01*(Z[1]):Z[end];

#TF=SH.*(1 .+humid_fac.*q)./rhoa./cpa.+humid_fac.*(Tv[1]./(1 .+humid_fac*q)).*LH./rhoa./Av;
TF=SH./rhoa./cpa.+humid_fac.*Tv[1].*LH./rhoa./Av;

#TF=(SH+LH)./rhoa./cpa;
ustar=abs(ustar)
L=Tv[1].*ustar.^3 ./(karman.*gravity_mks.*TF);
HD=eps.*ZD./L;
H=eps.*ZF[2:end]./L;

wsd=ustar./holtslag_psis.(HD);
Ts=Tv[1].+max.(min.(a.*TF./wsd,3),0);

spl = Spline1D(ZM, U); UD=spl(ZD);
spl = Spline1D(ZM, V);  VD=spl(ZD);
spl = Spline1D(ZM, Tv); TD=spl(ZD);

ha=abs.(Rcrit.*Ts.*(UD.^2+VD.^2)./(gravity_mks.*(TD-Ts)))-ZD;
ha[ha.>0].=Inf;
if all(isinf.(ha))
    if TF>0
        id=SZM;
    else
        id=1;
    end
else
    ~,id=findmin(abs.(ha));
end
id=max(id,1);
h=abs(ZD[id[1]]);
n=(h.-ZF).*((h.-ZF).>0);
n[n.==0].=Inf;
~,i=findmin(n);
frac=(h-ZF[i])./DZ[i];

nsl=(0.1 .*h.-ZF).*((0.1 .*h.-ZF).>0);
nsl[nsl.==0].=Inf;
~,isl=findmin(nsl);

ipbl=[1:i;];
if i==length(ZM)
    npbl=[];
else
    npbl=i+1:length(ZM);
end

kM=zeros(SZM);
kS=zeros(SZM);
wm=zeros(SZM);
ws=((gravity_mks/Tv[1]).*abs(TF).*h).^(1/3);

wm[1:isl]=ustar./holtslag_psim.(H[1:isl]);
wm[isl+1:i].=(ustar.^3 .+0.6.*ws.^3).^(1/3);

kM[ipbl]=karman.*(ZF[ipbl]+DZ[ipbl]/2).*
     wm[i].*(1 .-(ZF[ipbl]+DZ[ipbl]/2)./h).^2;
Pr=(holtslag_psis.(H[isl])./
    holtslag_psim.(H[isl])+
    a.*karman.*0.1.*ws./wm[isl]);
kS[ipbl]=kM[ipbl]./Pr;

S=abs.(diff([0; U])./diff(ZF));
GRi=(gravity_mks./Tv).*diff([Tv[1]; Tv])./diff(ZF)./S.^2;

fm=zeros(SZM);
fm[GRi.>=0]=1 ./(1 .+10 .*GRi[GRi.>=0].*(1 .+8 .*GRi[GRi.>=0]));
fm[GRi.<0]=(1 .-18 .*GRi[GRi.<0]).^(0.5);

lc=1 ./((1 ./(karman.*diff(ZF))).+(1 ./(karman.*lamdac)));

kM[npbl]=lc[npbl].^2 .*S[npbl].*fm[npbl];
kS[npbl]=kM[npbl];

kM[i]=kM[i].*frac+lc[i].^2 .*S[i].*fm[i].*(1-frac);
kS[i]=kS[i].*frac+lc[i].^2 .*S[i].*fm[i].*(1-frac);

return kM,kS
#[i max(kM) max(kS) TF max(Ts) U(1)]

end

function holtslag_psim(xi)
    if xi > 0
        y = 1 + 5 * xi
    else
        y = (1 - 16 * xi)^(-1 / 4)
    end
    return y
end

function holtslag_psis(xi)
    if xi > 0
        y = 1 + 5 * xi
    else
        y = (1 - 16 * xi)^(-1 / 2)
    end
    return y
end


text:
# Should I be using ZygoteRules instead?

Zygote.@adjoint function Zygote.literal_getproperty(x::ComponentArray, ::Val{s}) where s
    function getproperty_adjoint(Δ)
        zero_x = zero(x)
        setproperty!(zero_x, s, Δ)
        return (zero_x, nothing)
    end

    return getproperty(x, s), getproperty_adjoint
end

Zygote.@adjoint getdata(x::ComponentArray) = getdata(x), Δ->ComponentArray(Δ, getaxes(x))

Zygote.@adjoint getaxes(x::ComponentArray) = getaxes(x), Δ->ComponentArray(getdata(x), Δ)

Zygote.@adjoint ComponentArray(data, axes) = ComponentArray(data, axes), Δ->(getdata(Δ), getaxes(Δ))

Zygote.@adjoint function Base.convert(::Type{CA}, x::ComponentArray) where CA<:ComponentArray
    return ComponentArray(getdata(x), getaxes(CA)), Δ->(CA, ComponentArray(getdata(Δ), getaxes(x)))
end

Zygote.refresh()


text:
function  restoring(r,x,br,N2)

#   constructs part of sparse matrix wave operator Lsp for the
#   restoring force. It involves the L_x^2 term in the wave eqn.
#
#
#   input
#   r       - radial grid (includes boundaries)
#   x       - latitude grid x = cos(theta) (excludes poles and equator)
#   br      - dimensionless radial field
#   N2      - dimensionless stratification (N/Omega)^2
#

   nr = length(r);
   nx = length(x);
   dx = x[2]-x[1];
   dx2 = dx * dx;
   rcmb = 1.0;

   # some useful arrays
   s = sqrt.(1 .- x.*x);
   t = s./x;
   b = 0.25 * br / (rcmb^2);   # less useful now that L = Rcmb

   # initialize arrays
   ndim = nx * nr;
   nelements = 3 * ndim;
   row = Array{Int64}(undef,nelements);
   col = Array{Int64}(undef,nelements);
   val = Array{ComplexF64}(undef,nelements);

#println("nr ",nr," nx ",nx)
#
#  fill matrix (including boundaries)
   count = 0;
   for l = 1 : nx
       for k = 1 : nr

           ipos = (l-1) * nr +  k;

           # second derivative in x (central difference)
           if (l>1 && l < nx)

             count = count + 1;
             row[count] = ndim + ipos;
             col[count] = ipos-nr;
             val[count] = t[l] * b[l] * N2[k] * br[l-1] * t[l-1] /dx2;

             count = count + 1;
             row[count] = ndim + ipos;
             col[count] = ipos;
             val[count] = -2*t[l] * b[l] * N2[k] * br[l] * t[l] /dx2;

             count = count + 1;
             row[count] = ndim + ipos;
             col[count] = ipos + nr;
             val[count] = t[l] * b[l] * N2[k] * br[l+1] * t[l+1] /dx2;

           end

           if (l == 1)   # center difference with pole condition; unequal dx)

               count = count+1;
               row[count] = ndim + ipos;
               col[count] = ipos;
               val[count] = -5*t[l]*b[l]*N2[k] * br[l] * t[l] /dx2;

               count = count + 1;
               row[count] = ndim + ipos;
               col[count] = ipos + nr;
               val[count] = (2)*t[l]*b[l]*N2[k] * br[l+1] * t[l+1] /dx2;

               count = count + 1;
               row[count] = ndim + ipos;
               col[count] = ipos + 2*nr;
               val[count] = -(1/5)*t[l]*b[l]*N2[k]* br[l+2] * t[l+2]/dx2;

           end

           if (l == nx)   # backward derivative

               count = count + 1;
               row[count] = ndim + ipos;
               col[count] = ipos;
               val[count] = -5*t[l] * b[l] * N2[k] * br[l] * t[l] /dx2;

               count = count + 1;
               row[count] = ndim + ipos;
               col[count] = ipos-nr;
               val[count] = (2)*t[l]*b[l]*N2[k] * br[l-1] * t[l-1] /dx2;

               count = count + 1;
               row[count] = ndim + ipos;
               col[count] = ipos - 2*nr;
               val[count] = -(1/5)*t[l]*b[l]*N2[k] * br[l-2]* t[l-2]/dx2;

           end

       end
   end

#  add boundary conditions on top and bottom radii
   rsp = sparse(row,col,val,2*ndim,2*ndim);

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves