text:
GeoInterface.coordinates(obj::Point) = Float64[obj.x, obj.y]
GeoInterface.coordinates(obj::PointM) = Float64[obj.x, obj.y]
GeoInterface.coordinates(obj::PointZ) = Float64[obj.x, obj.y, obj.z]

GeoInterface.coordinates(obj::MultiPoint) =
    Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points]
GeoInterface.coordinates(obj::MultiPointM) =
    Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points]
GeoInterface.coordinates(obj::MultiPointZ) =
    Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points]

function GeoInterface.coordinates(obj::Polyline)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[nparts] <= npoints
    push!(obj.parts, npoints)
    coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts]
    pop!(obj.parts)
    coords
end

function GeoInterface.coordinates(obj::PolylineM)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[nparts] <= npoints
    push!(obj.parts, npoints)
    coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts]
    pop!(obj.parts)
    coords
end

function GeoInterface.coordinates(obj::PolylineZ)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[nparts] <= npoints
    push!(obj.parts, npoints)
    coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts]
    pop!(obj.parts)
    coords
end

# Only supports 2D geometries for now
# pt is [x,y] and ring is [[x,y], [x,y],..]
# ported from https://github.com/Turfjs/turf/blob/3d0efd96d59878f82c6d6baf8ed75695e0b6dbc0/packages/turf-inside/index.js#L76-L91
function inring(pt::Vector{Float64}, ring::Vector{Vector{Float64}})
    intersect(i::Vector{Float64}, j::Vector{Float64}) =
        (i[2] >= pt[2]) != (j[2] >= pt[2]) && (pt[1] <= (j[1] - i[1]) *
                                                        (pt[2] - i[2]) /
                                                        (j[2] - i[2]) + i[1])
    isinside = intersect(ring[1], ring[end])
    for k = 2:length(ring)
        isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside
    end
    isinside
end

# ported from https://github.com/Esri/terraformer-arcgis-parser/blob/master/terraformer-arcgis-parser.js#L168-L253
function GeoInterface.coordinates(obj::Polygon)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[end] <= npoints
    coords = Vector{Vector{Vector{Float64}}}[]
    holes = Vector{Vector{Vector{Float64}}}()
    push!(obj.parts, npoints)
    for i = 1:nparts
        ringlength = obj.parts[i+1] - obj.parts[i]
        ringlength < 4 && continue
        test = 0.0
        ring = Vector{Vector{Float64}}()
        for j = (obj.parts[i]+1):(obj.parts[i+1]-1)
            prev = obj.points[j]
            cur = obj.points[j+1]
            test += (cur.x - prev.x) * (cur.y + prev.y)
            push!(ring, GeoInterface.coordinates(prev))
        end
        push!(ring, GeoInterface.coordinates(obj.points[obj.parts[i+1]]))
        @assert length(ring) == ringlength
        if test > 0 # clockwise
            push!(coords, Vector{Vector{Float64}}[ring]) # create new polygon
        else # anti-clockwise
            push!(holes, ring)
        end
    end
    pop!(obj.parts)
    nrings = length(coords)
    for hole in holes
        for i = 1:nrings
            if inring(hole[1], coords[i][1])
                push!(coords[i], hole)
                break
            end
        end
        push!(coords, Vector{Vector{Float64}}[hole]) # hole is not inside any ring; make it a polygon
    end
    coords
end

# copied from GeoInterface.coordinates(obj::Polygon) to match signature here
function GeoInterface.coordinates(obj::PolygonM)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[end] <= npoints
    coords = Vector{Vector{Vector{Float64}}}[]
    holes = Vector{Vector{Vector{Float64}}}()
    push!(obj.parts, npoints)
    for i = 1:nparts
        ringlength = obj.parts[i+1] - obj.parts[i]
        ringlength < 4 && continue
        test = 0.0
        ring = Vector{Vector{Float64}}()
        for j = (obj.parts[i]+1):(obj.parts[i+1]-1)
            prev = obj.points[j]
            cur = obj.points[j+1]
            test += (cur.x - prev.x) * (cur.y + prev.y)
            push!(ring, GeoInterface.coordinates(prev))
        end
        push!(ring, GeoInterface.coordinates(obj.points[obj.parts[i+1]]))
        @assert length(ring) == ringlength
        if test > 0 # clockwise
            push!(coords, Vector{Vector{Float64}}[ring]) # create new polygon
        else # anti-clockwise
            push!(holes, ring)
        end
    end
    pop!(obj.parts)
    nrings = length(coords)
    for hole in holes
        for i = 1:nrings
            if inring(hole[1], coords[i][1])
                push!(coords[i], hole)
                break
            end
        end
        push!(coords, Vector{Vector{Float64}}[hole]) # hole is not inside any ring; make it a polygon
    end
    coords
end

# copied from GeoInterface.coordinates(obj::Polygon) to match signature here
function GeoInterface.coordinates(obj::PolygonZ)
    npoints = length(obj.points)
    nparts = length(obj.parts)
    @assert obj.parts[end] <= npoints
    coords = Vector{Vector{Vector{Float64}}}[]
    holes = Vector{Vector{Vector{Float64}}}()
    push!(obj.parts, npoints)
    for i = 1:nparts
        ringlength = obj.parts[i+1] - obj.parts[i]
        ringlength < 4 && continue
        test = 0.0
        ring = Vector{Vector{Float64}}()
        for j = (obj.parts[i]+1):(obj.parts[i+1]-1)
            prev = obj.points[j]
            cur = obj.points[j+1]
            test += (cur.x - prev.x) * (cur.y + prev.y)
            push!(ring, GeoInterface.coordinates(prev))
        end
        push!(ring, GeoInterface.coordinates(obj.points[obj.parts[i+1]]))
        @assert length(ring) == ringlength
        if test > 0 # clockwise
            push!(coords, Vector{Vector{Float64}}[ring]) # create new polygon
        else # anti-clockwise
            push!(holes, ring)
        end
    end
    pop!(obj.parts)
    nrings = length(coords)
    for hole in holes
        for i = 1:nrings
            if inring(hole[1], coords[i][1])
                push!(coords[i], hole)
                break
            end
        end
        push!(coords, Vector{Vector{Float64}}[hole]) # hole is not inside any ring; make it a polygon
    end
    coords
end

# bbox
const HasMBR = Union{Polyline,PolylineM, PolylineZ, Polygon, PolygonM, PolygonZ,
                     MultiPoint, MultiPointM, MultiPointZ, MultiPatch, Handle}

function GeoInterface.bbox(x::HasMBR)
    rect = x.MBR
    return rect.left, rect.bottom, rect.right, rect.top
end


text:
module MetricsTests

using Photon, Test
using Photon.Layers
using Photon.Metrics
using Photon.Losses
using Photon: history

function simple_conv_model()
    model = Sequential(
        Conv2D(4, 3, :relu),
        Conv2D(8, 3, :relu),
        MaxPool2D(),
        Dense(16, :relu),
        Dense(10)
    )
    return model
end

function getdata(s=28)
    [(
        KorA(randn(Float32,s,s,1,16)),
        KorA(randn(Float32,10,16))
    ) for i=1:10]
end

function test_core()
    model = simple_conv_model()
    workout = Workout(model, MSELoss())

    data = getdata()
    train!(workout, data, epochs=2)
    h = history(workout, :loss)

    @test h isa Tuple
    @test length(h[1]) == length(h[2])
end

function test_metrics()
    y_pred = rand(10,16)
    y_true = rand(0:1,10,16)
    b = BinaryAccuracy()
    acc = b(y_pred, y_true)

    o = OneHotBinaryAccuracy()
    cc = o(y_pred, y_true)

end

function test_algo()
    pred   = [[0.3 0.7]; [0. 1.]; [0.4 0.6]]'
    labels = [[1 0];[0 1];[0 1]]'
    loss = CrossEntropyLoss()
    @assert loss(pred, labels) ≈ 0.5715992760

    pred   = reshape(Array([3, -0.5, 2, 7]), (4,1))
    labels = reshape(Array([2.5, 0.0, 2, 8]), (4,1))
    loss = L1Loss()
    @assert loss(pred,labels) ≈ 0.5

    pred   = reshape(Array([3, -0.5, 2, 7]), (4,1))
    labels = reshape(Array([2.5, 0.0, 2, 8]), (4,1))
    loss = L2Loss()
    @assert loss(pred,labels) ≈ 0.375

    pred   = [[0.3 0.7]; [0. 1.]; [0.4 0.6]]
    labels = [[1 0];[0 1];[0 1]]
    loss = BinaryAccuracy()
    @assert loss(pred, labels) ≈ 0.666666666666

end

@testset "Metrics" begin
    test_core()
    test_metrics()
    test_algo()
end

end

ext: jl


text:
#!/usr/bin/env julia
using Pope, ZMQ, DocOpt

# doc = """
# zmq_listener
# Provide <lo> and <hi> in eV, and point to a <calfile> that has filt_value to energy calibrations.
# This is a sketch of how a popepipe_zmq<->spec interface would work it is not functional
# Usage:
#   zmq_listener.jl <lo> <hi> <calfile>
#
# Options:
#
# """
# arguments = docopt(doc, version=v"0.0.1")
# calfilename = expanduser(arguments["<calfile>"])
# lo = parse(Float64,arguments["<lo>"])
# hi = parse(Float64,arguments["<hi>"])
# This is a sketch of how a popepipe_zmq<->spec interface
# would work
# it is not functional

ctx=Context()
socket = ZMQ.Socket(ctx, ZMQ.SUB)
ZMQ.subscribe(socket)
ZMQ.connect(socket, "tcp://localhost:$(Pope.ZMQ_PORT)")
function recv_multipart(socket)
  out = ZMQ.Message[]
  push!(out, recv(socket))
  while ZMQ.ismore(socket)
    push!(out, recv(socket))
  end
  out
end
recv_dataproduct(socket) = make_dataproduct(recv_multipart(socket))
function make_dataproduct(out)
  channel_string = String(take!(convert(IOStream, out[1])))
  channel_number = parse(Int, channel_string)
  dp = read(seekstart(convert(IOStream,out[2])),Pope.MassCompatibleDataProductFeb2017,1)[1]
  return channel_string, channel_number, dp
end
isheader(out) = startswith(unsafe_string(out[1]),"header")
type ROI
  name::String
  lo::Float64
  hi::Float64
  count::Int
end
inroi(roi::ROI, v) = roi.lo <= v <= roi.hi
Base.push!(roi::ROI, v) = (roi.count+=inroi(roi,v))
reportcounts(roi::ROI) = (println("$(roi.name) $(roi.count)");roi.count=0)
roi1 = ROI("test ROI 5k to 10k",5000,10000,0)
roi2 = ROI("test ROI 0 to ∞",0,Inf,0)
roi3 = ROI("input ROI $lo to $hi",lo,hi,0)
rois = [roi1,roi2,roi3]

"stand in for actual calibration"
applycal(dp,channel_number::Int) = dp.filt_value

"stand in for actual cuts"
iscut(dp,channel_number) = false

endchannel = Channel{Bool}(1)

@schedule while !isready(endchannel)
  out = recv_multipart(socket)
  isheader(out) && continue
  channel_string, channel_number, dp = make_dataproduct(out)
  energy = applycal(dp, channel_number)
  iscut(dp,channel_number) && continue
  for roi in rois push!(roi,energy) end
end

@schedule while !isready(endchannel)
  sleep(1)
  reportcounts.(rois)
end

sleep(10)
out = recv_multipart(socket)
channel_string, channel_number, dp = make_dataproduct(out)
print(dp)


text:
using SubmodularMaximization
using POMDPs
using MCTS

grid_size = 10
iterations = 10
horizon = SubmodularMaximization.default_horizon
max_num_robots = 10
num_partitions = 3

grid = Grid(grid_size, grid_size)
sensor = RangingSensor()

histogram_filter = Filter(grid)

for num_robots = 1:max_num_robots
  println("Running tests for ", num_robots, " robots")

  robot_states = map(x->random_state(grid), 1:num_robots)

  problem = MultiRobotTargetTrackingProblem(grid, sensor, horizon,
                                            [histogram_filter], robot_states)

  println("Solving sequentially")
  solution = solve_sequential(problem)

  println("Solving with ", num_partitions, " partitions")

  println("Serially")
  solve_n_partitions(num_partitions, problem)

  println("Threaded")
  solve_n_partitions(num_partitions, problem, threaded=true)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
struct RtMidiWrapper
    ptr::Ptr{Cvoid}
    data::Ptr{Cvoid}
    ok::Cint
    msg::Ptr{Cchar}
end

const RtMidiPtr = Ptr{RtMidiWrapper}

const RtMidiInPtr = Ptr{RtMidiWrapper}

const RtMidiOutPtr = Ptr{RtMidiWrapper}

@cenum RtMidiApi::UInt32 begin
    RTMIDI_API_UNSPECIFIED = 0
    RTMIDI_API_MACOSX_CORE = 1
    RTMIDI_API_LINUX_ALSA = 2
    RTMIDI_API_UNIX_JACK = 3
    RTMIDI_API_WINDOWS_MM = 4
    RTMIDI_API_RTMIDI_DUMMY = 5
    RTMIDI_API_NUM = 6
end

@cenum RtMidiErrorType::UInt32 begin
    RTMIDI_ERROR_WARNING = 0
    RTMIDI_ERROR_DEBUG_WARNING = 1
    RTMIDI_ERROR_UNSPECIFIED = 2
    RTMIDI_ERROR_NO_DEVICES_FOUND = 3
    RTMIDI_ERROR_INVALID_DEVICE = 4
    RTMIDI_ERROR_MEMORY_ERROR = 5
    RTMIDI_ERROR_INVALID_PARAMETER = 6
    RTMIDI_ERROR_INVALID_USE = 7
    RTMIDI_ERROR_DRIVER_ERROR = 8
    RTMIDI_ERROR_SYSTEM_ERROR = 9
    RTMIDI_ERROR_THREAD_ERROR = 10
end

# typedef void ( * RtMidiCCallback ) ( double timeStamp , const unsigned char * message , size_t messageSize , void * userData )
const RtMidiCCallback = Ptr{Cvoid}

function rtmidi_get_compiled_api(apis, apis_size)
    ccall((:rtmidi_get_compiled_api, rtmidi_jll.rtmidi), Cint, (Ptr{RtMidiApi}, Cuint), apis, apis_size)
end

function rtmidi_api_name(api)
    ccall((:rtmidi_api_name, rtmidi_jll.rtmidi), Ptr{Cchar}, (RtMidiApi,), api)
end

function rtmidi_api_display_name(api)
    ccall((:rtmidi_api_display_name, rtmidi_jll.rtmidi), Ptr{Cchar}, (RtMidiApi,), api)
end

function rtmidi_compiled_api_by_name(name)
    ccall((:rtmidi_compiled_api_by_name, rtmidi_jll.rtmidi), RtMidiApi, (Ptr{Cchar},), name)
end

function rtmidi_error(type, errorString)
    ccall((:rtmidi_error, rtmidi_jll.rtmidi), Cvoid, (RtMidiErrorType, Ptr{Cchar}), type, errorString)
end

function rtmidi_open_port(device, portNumber, portName)
    ccall((:rtmidi_open_port, rtmidi_jll.rtmidi), Cvoid, (RtMidiPtr, Cuint, Ptr{Cchar}), device, portNumber, portName)
end

function rtmidi_open_virtual_port(device, portName)
    ccall((:rtmidi_open_virtual_port, rtmidi_jll.rtmidi), Cvoid, (RtMidiPtr, Ptr{Cchar}), device, portName)
end

function rtmidi_close_port(device)
    ccall((:rtmidi_close_port, rtmidi_jll.rtmidi), Cvoid, (RtMidiPtr,), device)
end

function rtmidi_get_port_count(device)
    ccall((:rtmidi_get_port_count, rtmidi_jll.rtmidi), Cuint, (RtMidiPtr,), device)
end

function rtmidi_get_port_name(device, portNumber)
    ccall((:rtmidi_get_port_name, rtmidi_jll.rtmidi), Ptr{Cchar}, (RtMidiPtr, Cuint), device, portNumber)
end

function rtmidi_in_create_default()
    ccall((:rtmidi_in_create_default, rtmidi_jll.rtmidi), RtMidiInPtr, ())
end

function rtmidi_in_create(api, clientName, queueSizeLimit)
    ccall((:rtmidi_in_create, rtmidi_jll.rtmidi), RtMidiInPtr, (RtMidiApi, Ptr{Cchar}, Cuint), api, clientName, queueSizeLimit)
end

function rtmidi_in_free(device)
    ccall((:rtmidi_in_free, rtmidi_jll.rtmidi), Cvoid, (RtMidiInPtr,), device)
end

function rtmidi_in_get_current_api(device)
    ccall((:rtmidi_in_get_current_api, rtmidi_jll.rtmidi), RtMidiApi, (RtMidiPtr,), device)
end

function rtmidi_in_set_callback(device, callback, userData)
    ccall((:rtmidi_in_set_callback, rtmidi_jll.rtmidi), Cvoid, (RtMidiInPtr, RtMidiCCallback, Ptr{Cvoid}), device, callback, userData)
end

function rtmidi_in_cancel_callback(device)
    ccall((:rtmidi_in_cancel_callback, rtmidi_jll.rtmidi), Cvoid, (RtMidiInPtr,), device)
end

function rtmidi_in_ignore_types(device, midiSysex, midiTime, midiSense)
    ccall((:rtmidi_in_ignore_types, rtmidi_jll.rtmidi), Cvoid, (RtMidiInPtr, Cint, Cint, Cint), device, midiSysex, midiTime, midiSense)
end

function rtmidi_in_get_message(device, message, size)
    ccall((:rtmidi_in_get_message, rtmidi_jll.rtmidi), Cdouble, (RtMidiInPtr, Ptr{Cuchar}, Ptr{Cint}), device, message, size)
end

function rtmidi_out_create_default()
    ccall((:rtmidi_out_create_default, rtmidi_jll.rtmidi), RtMidiOutPtr, ())
end

function rtmidi_out_create(api, clientName)
    ccall((:rtmidi_out_create, rtmidi_jll.rtmidi), RtMidiOutPtr, (RtMidiApi, Ptr{Cchar}), api, clientName)
end

function rtmidi_out_free(device)
    ccall((:rtmidi_out_free, rtmidi_jll.rtmidi), Cvoid, (RtMidiOutPtr,), device)
end

function rtmidi_out_get_current_api(device)
    ccall((:rtmidi_out_get_current_api, rtmidi_jll.rtmidi), RtMidiApi, (RtMidiPtr,), device)
end

function rtmidi_out_send_message(device, message, length)
    ccall((:rtmidi_out_send_message, rtmidi_jll.rtmidi), Cint, (RtMidiOutPtr, Ptr{Cuchar}, Cint), device, message, length)
end


text:
function scale(wp::WaveletParams, s::Int)
    WaveletParams(wp.s * s, wp.N * s, wp.J, wp.L, wp.K, wp.σ)
end

function get_params(x0, wp_x)
    vs_x = v_λ_k_all(x0, wp_x)
    k_x = K_all(x0, wp_x, vs_x)
    w_x = 1 ./ W_all(x0, wp_x, vs_x)

    if maximum(w_x) > 10_000
        w_x = min.(w_x, 10_000)
        @warn "weights are too big"
    end

    (;
        x0,
        wp_x,
        vs_x,
        k_x,
        w_x,
    )
end

function get_colored_params(x0, y0, wp_x, wp_y=wp_x, wp_xy=wp_y; s::Int=2)
    vs_x = v_λ_k_all(x0, wp_x)
    vs_y = v_λ_k_all(y0, wp_y)
    vs_x_cross = v_λ_k_all(x0, wp_xy)
    vs_y_cross = v_λ_k_all(y0, wp_xy)

    k_x = K_all(x0, wp_x, vs_x)
    k_y = K_all(y0, wp_y, vs_y)
    k_y_cross = K_all(x0, y0, wp_xy, vs_x_cross, vs_y_cross)

    # w_x = 1 ./ abs.(k_x)
    # w_y = 1 ./ abs.(k_y)
    # w_y_cross = 1 ./ abs.(k_y_cross)
    w_x = 1 ./ W_all(x0, wp_x, vs_x)
    w_y = 1 ./ W_all(y0, wp_y, vs_y)
    w_y_cross = 1 ./ W_all(x0, y0, wp_xy, vs_x_cross, vs_y_cross)

    if maximum([w_x; w_y; w_y_cross]) > 10_000
        w_x = min.(w_x, 10_000)
        w_y = min.(w_y, 10_000)
        w_y_cross = min.(w_y_cross, 10_000)
        @warn "weights are too big"
    end

    (;
        x0,
        y0,
        wp_x = scale(wp_x, s),
        wp_y = scale(wp_y, s),
        wp_xy = scale(wp_xy, s),
        vs_x,
        vs_y,
        vs_x_cross,
        vs_y_cross,
        k_x,
        k_y,
        k_y_cross,
        w_x,
        w_y,
        w_y_cross
    )
end


text:
struct Pulse{T}
    data #pulse envelope
    FWHM::T
    λ0::T
    τ::Array{T,1}
    P0::T # peak power
    τ0::T
    energy::T
end

function createTimeGrid(time_window; npts = 2^14)
    dt = time_window/npts
    τ = LinRange(-time_window/2, time_window/2, npts)
    τ = collect(τ)
end

function gaussianPulse(;τ= τ, energy=energy, FWHM=FWHM, λ0=λ0)
    τ0 = FWHM/1.664
    P0 = 0.94*energy/FWHM
    res0 = sqrt(P0).*sqrt.((exp.(-(τ.^2.0) ./ (τ0^2.0))))
    res::Pulse = Pulse{Float64}(res0, FWHM, λ0, τ, P0, τ0, energy)
    res
end

function chirpedGaussian(;τ= τ, energy=energy, FWHM=FWHM, λ0=λ0, C=C)
    τ0 = FWHM/1.664
    P0 = 0.94*energy/FWHM
    res0 = sqrt(P0).*sqrt.((exp.(-(1+1im*C)*(τ.^2.0) ./ (τ0^2.0))))
    res::Pulse = Pulse{Float64}(res0, FWHM, λ0, τ, P0, energy)
    res
end

function sechPulse(;τ= τ, energy=energy, FWHM=FWHM, λ0=λ0)
    τ0 = FWHM/1.76
    P0 = energy/(2*τ0)
    res0 = sqrt(P0).*sech.(τ/τ0)
    res::Pulse = Pulse{Float64}(res0, FWHM, λ0, τ, P0, τ0, energy)
    res
end

function chirpSechPulse(;τ= τ, N=N, FWHM=FWHM, λ0=λ0, C=C)
    τ0 = FWHM/1.76
    # P0 = energy/(2*τ0)
    res0 = @. N*sech(τ/τ0)*exp(-1im*C*(τ/τ0)^2/2)
    res::Pulse = Pulse{Float64}(res0, FWHM, λ0, τ, N, τ0, N)
    res
end

function pulseFromFile(file::String; timeFactor=1e-15, npts=2^14, energy=1e-6, phase_sign = 1.0)
    data = readdlm(file, Float64)
    # if length(data[1,:]) < 20
    #     time = data[1,:]
    #     intensity = data[2,:]
    #     phase = data[3,:]
    # else
    #     time = data[:1,]
    #     intensity = data[:2,]
    #     phase = data[:3,]
    # end
    timep = data[:,1]
    intensity = data[:,2]
    phase = data[:,3]
    timep *= timeFactor
    timeGrid = LinRange(-1e-12,1e-12,npts)
    # print(time)
    spl_intensity = Spline1D(timep, intensity)
    spl_phase = Spline1D(timep, phase)
    generated_intensity = spl_intensity(timeGrid)
    generated_phase = spl_phase(timeGrid)
    generated_intensity[generated_intensity.<0] .= 0
    xprofile = sqrt.(generated_intensity) .* exp.(phase_sign*1im .* generated_phase)
    norm = trapz(timeGrid, abs2.(xprofile))
    res = xprofile.*sqrt(energy)./sqrt(norm)
    P0 = maximum(abs2.(res))
    fwhm = FWHM(timeGrid, abs2.(res))
    τ0 = fwhm/1.664
    pulse = Pulse{Float64}(res, fwhm, 800e-9, timeGrid, P0, τ0, energy)
    return pulse
end

function FWHM(pulse::Pulse)
    # for this function to be accurte, the number of points in the time grid must be >= 2^14
    a = (pulse.data).^2
    hm = maximum(a)/2.0
    out = pulse.τ[findlast(x -> x > hm, a)] - pulse.τ[findfirst(x -> x > hm, a)]
    out
end

function FWHM(τ, data)
    # for this function to be accurte, the number of points in the time grid must be >= 2^14
    a = abs2.(data)
    hm = maximum(a)/2.0
    out = τ[findlast(x -> x > hm, a)] - τ[findfirst(x -> x > hm, a)]
    out
end

function Waist(τ, data)
    # for this function to be accurte, the number of points in the time grid must be >= 2^14
    a = abs2.(data)
    hm = maximum(a)/2.718281828459045
    out = τ[findlast(x -> x > hm, a)] - τ[findfirst(x -> x > hm, a)]
    out
end

function getElectricField(pulse::Pulse, fiber) #avergae electric field
    ω0 = 2pi*c/pulse.λ0
    Et = @. pulse.data * exp(-1im * ω0 * pulse.τ)
    Et = real(Et)
    Et = sqrt(2/(ϵ0*c*effectiveArea(fiber.diameter)))*Et
    return Et
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
#=
Example of defining a docstring for a function. The docstring is rendered in
the Julia REPL by typing `?` then the name of the function.
=#

"""
    sayhello(x, y)

Print a hello greeting given a first `x` and last `y` name.

# Inputs

- `x`: the first name.
- `y`: the last name.

### Examples

```
julia> sayhello("Bart", "Simpson")
Saying hello to Bart Simpson.
```

"""
function sayhello(x, y)
    s = "Saying hello to $x $y."
    println(s)
end

a = "Bart"
b = "Simpson"
sayhello(a, b)

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
# 3-step file upload as documenter here: https://canvas.instructure.com/doc/api/file.file_uploads
function Internals.upload_file(endpoint::String, file::String; api::CanvasAPI=getapi(), params::Dict=Dict{String,Any}(), kwargs...)
    params = Dict{String,Any}(
        "name" => basename(file),
        "size" => stat(file).size,
        params...
    )

    # Step 1: Telling Canvas about the file upload and getting a token
    json = request("POST", endpoint; api=api, params=params, kwargs...)

    # Step 2: Upload the file data to the URL given in the previous response
    uri = HTTP.URI(json["upload_url"])
    req = open(file, "r") do io
        body = collect(json["upload_params"])
        push!(body, "file"=>io) # file must be last argument
        form = HTTP.Form(body)
        headers = Internals.canvas_headers(; auth=nothing) # should not be authenticated
        return HTTP.request("POST", uri, headers, form)
    end

    # Step 3: Confirm the upload's success
    if HTTP.isredirect(req)
        # 3XX Redirect
        loc = HTTP.header(req, "Location")
        json = request("GET", loc; api=api) # Should be authenticated
    else # req.status == 201
        # 201 Created
        json = JSON.parse(HTTP.payload(req, String))
    end
    return File(json)
end


text:
p1 = PortHamiltonian.Phs([0], [1 1], [0 0; 0 0], [1])
p2 = PortHamiltonian.Phs([1], [1], [0], [1])

phct = PortHamiltonian.coupled_transformer(p1,[1], p2, [1],[1])

@test (phct.G == ones(2,1))

@test (phct.J == [0. 0.; 0. 1.])

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
###############################################
# Solve laplace equation

# by: Jose Miguel Muñoz
# please contact me if this doesn't works
###############################################

using Plots, Statistics, BenchmarkTools
###############################################
#       BOUNDARY
###############################################
N = 100
a , b = 0,2

Rx = a:(b-a)/(N-1):b
Ry = Rx

boundary_X_Up = map(rx -> 0,Rx)

boundary_X_Down = map(rx -> 0, Rx)

boundary_Y_R = map(rx -> 10,Rx)

boundary_Y_L = map(rx -> 0,Rx)

V = zeros(N, N)

V[1,:] = boundary_X_Up
V[:,1] = boundary_Y_R
V[N,:] = boundary_X_Down
V[:,N] = boundary_Y_L

###############################################
#       SOLVER
###############################################
function solver!(V,n_inter,N=N, tol=1)
    for rep = 1:n_inter

        for row = 2:N-1
            for column = 2:N-1
                V[row, column] = mean([V[row-1, column],V[row+1, column],V[row, column-1],V[row, column+1]])
                #if (row+10)^2+(column+10)^2<1000
                #    V[row, column]=1
                #end
    end end
 end
    return V
end

###############################################
#       Ask for the boundaries
###############################################
@benchmark solver!(V,1000)
V = solver!(V,100000)

plot(Ry, Ry, V, linetype=:contourf,title = "Solucion por Gauss_Seidel")

#gr()
#anim = @animate for n = 1:100:1000
#    heatmap(vss[:,:,n])
#end
#gif(anim)
using Statistics
plot(Ry, Ry, (V-Space*1.01), linetype=:contourf,title = "Solucion por Jacobi")
mean(V-Space)
mean(V.^2-Space.^2)
std(V.^2-Space.^2)

#######################################################
#       Convergence Analysis
#######################################################
listconvGauss=[]
for n_reps=500:100:5000
    V = rand(N, N)

    V[1,:] = boundary_X_Up
    V[:,1] = boundary_Y_R
    V[N,:] = boundary_X_Down
    V[:,N] = boundary_Y_L

    V = solver!(V,n_reps)
    append!(listconvGauss,abs(mean(V-Space)))
end
listconvGauss

alistconvSOR=map(x->abs(x),listconvSOR)
alistconvJacobi=map(x->abs(x),listconvJacobi)
alistconvGauss=map(x->abs(x),listconvGauss)
plot()
plot!(Vector(500:100:5000),alistconvGauss,label="Metodo Gauss_Seidel")
plot!(ps,alistconvSOR[begin+6:end],label="Metodo SOR")
plot!(ps,alistconvJacobi.+randn(),label="Metodo Jacobi")
xlabel!("Numero de iteraciones")
ylabel!("Δφ")
ps
alistconvSOR[begin:end-6]

alistconvJacobi
ps

max_stars_repo_name: munozariasjm/Laplace-and-Julia

max_issues_repo_name: munozariasjm/Laplace-and-Julia

max_forks_repo_name: munozariasjm/Laplace-and-Julia


text:
# Modelo de desaerador (tanque) com 5 entradas e sa죡da l찵quido saturado (x=0)
type Desaerador5
	Desaerador5()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in1}",
				:PosX=>1,
				:PosY=>0.2
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in2}",
				:PosX=>1,
				:PosY=>0.45
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in3}",
				:PosX=>1,
				:PosY=>0.7
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in4}",
				:PosX=>1,
				:PosY=>0.9
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in5}",
				:PosX=>0.5,
				:PosY=>0
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>0,
				:PosY=>0.85
			)),
			[
				:(Fout.F = Fin1.F + Fin2.F + Fin3.F + Fin4.F + Fin5.F),
				:(Fout.F * Fout.H = Fin1.F * Fin1.H + Fin2.F * Fin2.H + Fin3.F * Fin3.H + Fin4.F * Fin4.H + Fin5.F * Fin5.H),
				:([Fout.S,Fout.H] = PP2.propPTl(Fout.P,Fout.T)),
			],
			[
				"","","",
			],
			[:PP2,],
			[:Fin1,:Fin2,:Fin3,:Fin4,:Fin5,:Fout,]
		)
	end
	PP2::DanaPlugin
	Fin1::Corrente
	Fin2::Corrente
	Fin3::Corrente
	Fin4::Corrente
	Fin5::Corrente
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Desaerador5
function setEquationFlow(in::Desaerador5)
	addEquation(1)
	addEquation(2)
	#Fout.T = PP2.Tsat(Fout.P);
	addEquation(3)
end
function atributes(in::Desaerador5,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/desaerador"
	drive!(fields,_)
	return fields
end
Desaerador5(_::Dict{Symbol,Any})=begin
	newModel=Desaerador5()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Contiguous on-device arrays

export oneDeviceArray, oneDeviceVector, oneDeviceMatrix

## construction

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're currently requiring a trailing pointer argument.

struct oneDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::LLVMPtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    oneDeviceArray{T,N,A}(shape::Dims{N}, ptr::LLVMPtr{T,A}) where {T,A,N} = new(shape,ptr)
end

const oneDeviceVector = oneDeviceArray{T,1,A} where {T,A}
const oneDeviceMatrix = oneDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
oneDeviceArray(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A})                where {T,A,N} = oneDeviceArray{T,N,A}(dims, p)
oneDeviceArray(len::Integer,              p::LLVMPtr{T,A})                where {T,A}   = oneDeviceVector{T,A}((len,), p)

# outer constructors, partially parameterized
oneDeviceArray{T}(dims::NTuple{N,<:Integer},   p::LLVMPtr{T,A}) where {T,A,N} = oneDeviceArray{T,N,A}(dims, p)
oneDeviceArray{T}(len::Integer,                p::LLVMPtr{T,A}) where {T,A}   = oneDeviceVector{T,A}((len,), p)
oneDeviceArray{T,N}(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A}) where {T,A,N} = oneDeviceArray{T,N,A}(dims, p)
oneDeviceVector{T}(len::Integer,               p::LLVMPtr{T,A}) where {T,A}   = oneDeviceVector{T,A}((len,), p)

# outer constructors, fully parameterized
oneDeviceArray{T,N,A}(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A}) where {T,A,N} = oneDeviceArray{T,N,A}(Int.(dims), p)
oneDeviceVector{T,A}(len::Integer,               p::LLVMPtr{T,A}) where {T,A}   = oneDeviceVector{T,A}((Int(len),), p)

## getters

Base.pointer(a::oneDeviceArray) = a.ptr
Base.pointer(a::oneDeviceArray, i::Integer) =
    pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:oneDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::oneDeviceArray) = g.shape
Base.length(g::oneDeviceArray) = prod(g.shape)

## conversions

Base.unsafe_convert(::Type{LLVMPtr{T,A}}, a::oneDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

## indexing intrinsics

# NOTE: these intrinsics are now implemented using plain and simple pointer operations;
#       when adding support for isbits union arrays we will need to implement that here.

# TODO: how are allocations aligned by the level zero API? keep track of this
#       because it enables optimizations like Load Store Vectorization
#       (cfr. shared memory and its wider-than-datatype alignment)

@inline function arrayref(A::oneDeviceArray{T}, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_load(pointer(A), index, Val(align))
end

@inline function arrayset(A::oneDeviceArray{T}, x::T, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_store!(pointer(A), x, index, Val(align))
    return A
end

@inline function const_arrayref(A::oneDeviceArray{T}, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_cached_load(pointer(A), index, Val(align))
end

## indexing

Base.@propagate_inbounds Base.getindex(A::oneDeviceArray{T}, i1::Int) where {T} =
    arrayref(A, i1)
Base.@propagate_inbounds Base.setindex!(A::oneDeviceArray{T}, x, i1::Int) where {T} =
    arrayset(A, convert(T,x)::T, i1)

Base.IndexStyle(::Type{<:oneDeviceArray}) = Base.IndexLinear()

## const indexing

"""
    Const(A::oneDeviceArray)

Mark a oneDeviceArray as constant/read-only. The invariant guaranteed is that you will not
modify an oneDeviceArray for the duration of the current kernel.

This API can only be used on devices with compute capability 3.5 or higher.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N,AS} <: DenseArray{T,N}
    a::oneDeviceArray{T,N,AS}
end
Base.Experimental.Const(A::oneDeviceArray) = Const(A)

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
Base.@propagate_inbounds Base.getindex(A::Const, i1::Int) = const_arrayref(A.a, i1)

# deprecated
Base.@propagate_inbounds ldg(A::oneDeviceArray, i1::Integer) = const_arrayref(A, Int(i1))

## other

Base.show(io::IO, a::oneDeviceVector) =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::oneDeviceArray) =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

Base.show(io::IO, mime::MIME"text/plain", a::oneDeviceArray) = show(io, a)

@inline function Base.iterate(A::oneDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end


text:
function _always_assert(condition::Bool)
    if !condition
        throw(ErrorException("assertion failed"))
    end
    return nothing
end

function _default_config()
    return Config(; pkg_server = "https://pkg.julialang.org")
end

function _contents_to_filename(; contents::String)
    filename = joinpath(mktempdir(; cleanup = true), "myfile")
    open(filename, "w") do io
        println(io, contents)
    end
    return filename
end

function project_contents_to_manifest_contents(cfg = _default_config();
                                               project_file_contents::AbstractString)::String
    proj_or_mnfst = Project(FileContents(project_file_contents))
    temporary_environment = PkgServerCacheWarmer.create_new_temporary_environment(cfg, proj_or_mnfst)
    project_directory = temporary_environment.env_with_project["JULIA_PROJECT"]::String

    proj_toml   = joinpath(project_directory, "Project.toml")
    man_toml    = joinpath(project_directory, "Manifest.toml")
    j_proj_toml = joinpath(project_directory, "JuliaProject.toml")
    j_man_toml  = joinpath(project_directory, "JuliaManifest.toml")

    _always_assert(isfile(proj_toml))
    _always_assert(!isfile(man_toml))
    _always_assert(!isfile(j_proj_toml))
    _always_assert(!isfile(j_man_toml))

    commands = PkgServerCacheWarmer.generate_julia_pkg_commands(
        cfg,
        temporary_environment,
        proj_or_mnfst,
    )
    PkgServerCacheWarmer.run_commands(commands)

    _always_assert(isfile(proj_toml))
    _always_assert(isfile(man_toml))
    _always_assert(!isfile(j_proj_toml))
    _always_assert(!isfile(j_man_toml))

    manifest_file_contents = read(man_toml, String)::String
    return manifest_file_contents
end


text:
using XxxxYyyy
using Test
using Documenter

@testset "XxxxYyyy.jl" begin

    @test hello("Julia") == "Hello, Julia"
    @test domath(2.0) ≈ 7.0

    # Documenter test of module documentation
    doctest(XxxxYyyy; manual=false)

end


text:
##################################################
#Experimental API
export Words, Sentences

const tokenizers = [poormans_tokenize, punctuation_space_tokenize,
       penn_tokenize, improved_penn_tokenize, nltk_word_tokenize, tweet_tokenize]
const sentence_splitters = [rulebased_split_sentences]

const Words = tokenize
const Sentences = split_sentences

for rule in [tokenizers; sentence_splitters]
    @eval function Base.split(str::T, ::typeof($(rule))) where T<:AbstractString
        $rule(str)
    end

    @eval function Base.split(str::T, ::typeof($(rule))) where T<:SubString
        $rule(str)
    end
end


text:
########################################################
function realresampleddataforpca()
########################################################

    A = h5read("Data/AGNs.hdf5","simulations/Mrk509/AGNGALresamp")

    x = Matrix(A[1:3, 2, :]')  # select filters out of available six

    Σ =  Matrix(A[1:3, 3, :]')

    return x, Σ

end

########################################################
function realresampleddataforblr()
########################################################

    A = h5read("Data/AGNs.hdf5","simulations/Mrk509/AGNGALresamp")

    x = Matrix(A[1:6, 2, :]')  # select filters out of available six

    srt = sortperm(vec(A[1,1,:])) # sort according to time,
                                  # not necessary, just convenient for plotting

    idx  = [ones(Int, 110)   ; ones(Int, 110)*2 ; ones(Int, 110)*3 ; ones(Int, 110)*4 ; ones(Int, 110)*5 ; ones(Int, 110)*6]
    flux = [vec(x[srt,1])    ; vec(x[srt,2])    ; vec(x[srt,3])    ; vec(x[srt,4])    ; vec(x[srt,5])    ; vec(x[srt,6])]
    Σ    = [vec(A[1, 3, srt]); vec(A[2, 3, srt]); vec(A[3, 3, srt]); vec(A[4, 3, srt]); vec(A[5, 3, srt]); vec(A[6, 3, srt])];

    t    = A[1,1,srt] .- minimum(A[1,1,srt]) # no need to subtract,
                                             # just convenient that times starts at 0

    idx1 = findall(idx .== 1)
    idx2 = findall(idx .== 2)
    idx3 = findall(idx .== 3)
    idx4 = findall(idx .== 4)
    idx5 = findall(idx .== 5)
    idx6 = findall(idx .== 6)

    figure(1)
    cla()
    title("plotting three first LSST filters")
    plot3D(flux[idx1], flux[idx2], flux[idx3], "bo", alpha=0.4)

    figure(2)
    cla()
    plot(t, flux[idx1], ".", label = "1 - u")
    plot(t, flux[idx2], ".", label = "2 - g")
    plot(t, flux[idx3], ".", label = "3 - r")
    plot(t, flux[idx4], ".", label = "4 - i")
    plot(t, flux[idx5], ".", label = "5 - y")
    plot(t, flux[idx6], ".", label = "6 - z")

    legend()

    # return [t;t;t; t;t;t], flux, Σ, idx
    return [t;t;t], flux[1:330], Σ[1:330], idx[1:330]

end


text:
using NLOptControl, VehicleModels, PrettyPlots, JuMP, Plots, Ipopt, DataFrames, Parameters
include("setup.jl")
include("DegeneracyHunter.jl")

function nlp1()
  pa=Vpara(x_min=185.,y_min=-0.1);
  n=NLOpt(); @unpack_Vpara pa
  X0=[  x0_,   y0_, v0_, r0_,   psi0_,   sa0_,  u0_ , ax0_];
  XF=[x_ref, y_ref, NaN, NaN,     NaN,    NaN,   u0_,  NaN];
  XL=[x_min, y_min, NaN, NaN, psi_min, sa_min, u_min, NaN];
  XU=[x_max, y_max, NaN, NaN, psi_max, sa_max, u_max, NaN];
  CL = [sr_min, jx_min]; CU = [sr_max, jx_max];
  n = define(n,stateEquations=ThreeDOFv2,numStates=8,numControls=2,X0=X0,XF=XF,XL=XL,XU=XU,CL=CL,CU=CU,tf_max=40.0)
  n = configure(n,Ni=2,Nck=[10,10];(:integrationMethod => :ps),(:integrationScheme => :lgrExplicit),(:finalTimeDV => true))
  #n = configure(n,N=60;(:integrationMethod => :tm),(:integrationScheme => :bkwEuler),(:finalTimeDV => true))
  #n = configure(n,N=60;(:integrationMethod => :tm),(:integrationScheme => :trapezoidal),(:finalTimeDV => true))
  mdl     = Model(solver = IpoptSolver(print_level=0,
                                       max_iter=500,
                                       tol=1e-1,
                                       dual_inf_tol=1100.,
                                       constr_viol_tol=1e-1,
                                       compl_inf_tol=1.,
                                       acceptable_tol=1e-3,
                                       acceptable_constr_viol_tol=0.01,
                                       acceptable_dual_inf_tol=1e10,
                                       acceptable_compl_inf_tol=0.01,
                                       acceptable_obj_change_tol=1e20,
                                       diverging_iterates_tol=1e20
                                       )
                      )
  names = [:x,:y,:v,:r,:psi,:sa,:ux,:ax];
  descriptions = ["X (m)","Y (m)","Lateral Velocity (m/s)", "Yaw Rate (rad/s)","Yaw Angle (rad)", "Steering Angle (rad)", "Longitudinal Velocity (m/s)", "Longitudinal Acceleration (m/s^2)"];
  stateNames(n,names,descriptions)
  names = [:sr,:jx];
  descriptions = ["Steering Rate (rad/s)","Longitudinal Jerk (m/s^3)"];
  controlNames(n,names,descriptions)
  params = [pa];   # vehicle parameters

  # define tolerances
  #XF_tol=[0.5,0.5,NaN,NaN,NaN,NaN,0.5,NaN];
  #X0_tol=[0.05,0.05,0.05,0.05,0.1,0.001,0.05,0.05];
  #defineTolerances(n;X0_tol=X0_tol,XF_tol=XF_tol);
  n,r = OCPdef(mdl,n,params);
  sr_obj = integrate(mdl,n,r.u[:,1];C=w_sr,(:variable=>:control),(:integrand=>:squared))
  @NLobjective(mdl, Min, n.tf + sr_obj)

  # obstacles
  Q = size(A)[1]; # number of obstacles
  @NLparameter(mdl, a[i=1:Q] == A[i]);
  @NLparameter(mdl, b[i=1:Q] == B[i]);
  @NLparameter(mdl, X_0[i=1:Q] == X0_obs[i]);
  @NLparameter(mdl, Y_0[i=1:Q] == Y0_obs[i]);

  # obstacle postion after the intial postion
  X_obs=@NLexpression(mdl, [j=1:Q,i=1:n.numStatePoints], X_0[j])
  Y_obs=@NLexpression(mdl, [j=1:Q,i=1:n.numStatePoints], Y_0[j])

  # constraint on position
  obs_con=@NLconstraint(mdl, [j=1:Q,i=1:n.numStatePoints-1], 1 <= ((r.x[(i+1),1]-X_obs[j,(i+1)])^2)/((a[j]+sm)^2) + ((r.x[(i+1),2]-Y_obs[j,(i+1)])^2)/((b[j]+sm)^2))
  newConstraint(r,obs_con,:obs_con);
  # TODO  make y0_ and x0_  JuMP parameters
  # ensure that all of the points are in the LiDAR range
  LiDAR_con=@NLconstraint(mdl, [i=1:n.numStatePoints-1], ((r.x[(i+1),1]-x0_)^2+(r.x[(i+1),2]-y0_)^2) <= (L_R + L_Rd)^2); # not constraining the first state
  newConstraint(r,LiDAR_con,:LiDAR_con);
  # ensure that the vehicle is at the edge of the LiDAR range at Tp
  #LiDAR_conF=@NLconstraint(mdl, (L_R - L_Rd)^2 <= ((r.x[n.numStatePoints,1]-x0_)^2+(r.x[n.numStatePoints,2]-y0_)^2) );
  #optimize(mdl,n,r)
  #dfs=Vector{DataFrame}(2);label_string = ["opt.","RK4"];
  #dfs[1]=x_to_DF(;t=(r.t_st),x=r.X[:,1],y=r.X[:,2],v=r.X[:,3],r=r.X[:,4],psi=r.X[:,5],sa=r.X[:,6],ux=r.X[:,7],ax=r.X[:,8],sr=r.U[:,1],jx=r.U[:,2]);
  #all_plots(dfs,["opt. inf"],results_dir,obs_data,s_data,pa)
  return mdl
end

m=nlp1()
tic()
status = solve(m)
tm = toq()
#=
DegeneracyHunter.printInfeasibleEquations(m, 0.1)
DegeneracyHunter.printInfeasibleEquations(mdl, 0.1)

DegeneracyHunter.printInactiveEquations(m)

ps = DegeneracyHunter.assembleProblemStats(m,status,tm)
DegeneracyHunter.printProblemStats(ps)
=#


text:
using FunctionEnvelope
using Test

# Run tests
println("Basic Tests")
@time @test include("basic_tests.jl")
println("Full Tests")
@time @test include("envelope_functions_tests.jl")
println("Basic Tests")
@time @test include("one_function_test.jl")


text:
nanmean(x) = mean(filter(!isnan,x))
nanmean(x,y) = mapslices(nanmean,x,y)
nansem(x) = sem(filter(!isnan,x))
nansem(x,y) = mapslices(nansem,x,y)

function safemean(x)
    if sum(.!(isnan.(x))) < 2
        return NaN
    else
        return nanmean(x)
    end
end

function safesem(x)
    if sum(.!(isnan.(x))) < 4
        return NaN
    else
        return nansem(x)
    end
end

iscolumn(df,x) = x in propertynames(df)

function model_res(m)
    DataFrame(Area = coefnames(m), P = m.pvalues, Z = m.beta ./ m.stderror)
end


text:
using PyCall, JSON

ee = pyimport("ee")
ee.Initialize()

kh = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017").filter(ee.Filter.eq("country_na","Cambodia"));
basins = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7").filterBounds(kh)

output_dir = "../testing"

files = filter(x->endswith(x, ".json"), readdir(output_dir))

features = Array{PyObject}(undef,length(files))

for i in 1:length(files)
    file = files[i]

    code = parse(Int,split(split(file,"_")[end],".")[1])

    f = open(file,"r")
    data_dict = JSON.parse(f)
    close(f)

    features[i] = ee.Feature(basins.filter(ee.Filter.eq("HYBAS_ID",code)).first()).set(data_dict)

end

fc = ee.FeatureCollection(features)

task = ee.batch.Export.table.toAsset(collection=fc, description="fier_testing_results", assetId="users/kelmarkert/fier_testing_resutls")
task.start()


text:
using YaoLocations
using MLStyle
using YaoHIR
using YaoHIR.IntrinsicOperation
using Test

module TestIntrinsic
using YaoHIR: @intrinsic

@intrinsic X
@intrinsic R(theta::T) where {T <: Real}
@intrinsic SWAP

end

@test YaoHIR.routine_name(TestIntrinsic.X) == :X
@test YaoHIR.routine_name(TestIntrinsic.R(1.0)) == :R
@test TestIntrinsic.R(1.0).theta == 1.0

display(X)
display(Rx(1.0))
display(YaoHIR.Operation(X, 2.0))

circ = Chain(
    Gate(X, Locations(1)),
    Core.SSAValue(1),
    Ctrl(Gate(Core.SSAValue(1), Locations(3)), CtrlLocations(2))
)

print(circ)

@test YaoHIR.leaves(circ) == [Gate(X, Locations(1)),
    Core.SSAValue(1),
    Ctrl(Gate(Core.SSAValue(1), Locations(3)), CtrlLocations(2))
]

@testset "test match" begin
    gate = Gate(X, Locations(2))

    @match gate begin
        Gate(op, locs) => begin
            @test op == X
            @test locs == Locations(2)
        end
    end

    ctrl = Ctrl(Gate(X, Locations(2)), CtrlLocations(3))

    @match ctrl begin
        Ctrl(Gate(op, locs), ctrl) => begin
            @test op == X
            @test locs == Locations(2)
            @test ctrl == CtrlLocations(3)
        end
    end
end

@testset "isequal" begin
    circuit1 = Chain(Gate(H, Locations((1, ))), Gate(H, Locations((1, ))))
    circuit2 = Chain(Gate(H, Locations((1, ))), Gate(H, Locations((1, ))))
    @test circuit1 == circuit2
end


text:
import Yao: gatecount

function gatecount!(c::Union{PauliString, Bag, Sequence}, storage::AbstractDict)
    (gatecount!.(c |> subblocks, Ref(storage)); storage)
end

function gatecount!(c::ConditionBlock, storage::AbstractDict)
    k = typeof(c)
    if haskey(storage, k)
        storage[k] += 1
    else
        storage[k] = 1
    end
    storage
end

gatecount!(c::TrivialGate, storage::AbstractDict) = storage

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
module escape_dimer
export escape_exit_num_dimer
using DifferentialEquations
using PolynomialRoots
using Printf
# using MATLAB    H=.17    H=.17

# using ProgressMeter

include("leap_frog_definitions.jl")
# include("Dimer_Eq_of_M.jl")
max_hit_q1 = 11
max_hit_p2 = 100
barrier = 5

H=.20
condition_max_hits_back_forward(u, t, integrator) = u[3] > max_hit_q1 || u[6] > max_hit_p2 #|| maximum([abs(u[1]),abs(u[2]),abs(u[4]),abs(u[5])]) > barrier

affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function condition_hits_PSS_p2(u, t, integrator) # Event when event_f(u,t) == 0
    u[5]
end
function affect_update_iterator_p2!(integrator)
    integrator.u[6] = integrator.u[6] + 1
end

function affect_update_iterator_q1!(integrator)
    q2 = integrator.u[2]
    p2 = integrator.u[5]
    p1 = integrator.u[4]
    P1 = P1_poly(q2, p2, H)
    if !isempty(P1) &&  abs(p1 - P1) < .5
        integrator.u[3] = integrator.u[3] + 1
    end
end

callback_max_hits_bf = DiscreteCallback(condition_max_hits_back_forward, affect_stop!)
callback_hits_PSS_q1_forward = ContinuousCallback(condition_hits_PSS_q1, nothing,  affect_update_iterator_q1!, rootfind = true)
callback_hits_PSS_q1_backward = ContinuousCallback(condition_hits_PSS_q1,  affect_update_iterator_q1!, nothing, rootfind = true)

callback_hits_PSS_p2_forward = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, rootfind = false)
callback_hits_PSS_p2_backward = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, rootfind = false)

cb_forward = CallbackSet(callback_hits_PSS_q1_forward, callback_hits_PSS_p2_forward, callback_max_hits_bf)
cb_backward = CallbackSet(callback_hits_PSS_q1_backward, callback_hits_PSS_p2_backward, callback_max_hits_bf)
# println(Q2)
# println(P2)
function escape_exit_num_dimer(mesh_list, t_end, H)

    Q2 = mesh_list[1]; P2 = mesh_list[2];
    P1 = P1_poly(Q2, P2, H)
    if ~isempty(P1)
        # println(Q1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        # q0[4] = H

        p0[1] = P1
        p0[2] = P2
        p0[3] = 0
        # p0[4] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));

        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol_f = solve(prob, Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15, callback = cb_forward, save_start = false, save_end = true, save_everystep = false)
        # sol_f = solve(prob, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10, callback = cb_forward, save_start = true, save_end = true, save_everystep = true)
        # sol_f=solve(prob, RK4(), maxiters = 1e20,callback=cb_forward,save_start=false,save_end=true,save_everystep=true)
        sol_f=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_forward,save_start=true,save_end=true,save_everystep=true)

        uf = zeros(6)
        uf[1] = sol_f[1,end]
        uf[2] = sol_f[2,end]
        uf[3] = sol_f[3,end]
        uf[4] = sol_f[4,end]
        uf[5] = sol_f[5,end]
        uf[6] = sol_f[6,end]
        q1 = uf[1]
        q2 = uf[2]
        p1 = uf[4]
        p2 = uf[5]

        q = [q1, q2]
        p = [p1, p2]
        forward_hits =  uf[3]
        # println(sol_f.t)

        diff = abs(Hamiltonian_Dimer(q, p, 0) - H)

        # prob_sign = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., 1));
        # sol_f_sign=solve(prob, RK4(),maxiters=1e20,callback=cb_forward,save_start=true,save_end=true,save_everystep=true))
        # sign_sol=sign(sol_f_sign[1,2])
        if length(sol_f[1,:]) > 1
            sign_sol = sign(sol_f[1,2])
        else
            sign_sol = 0
        end

        # if diff<1e-10 && maximum([abs(q1),abs(q2),abs(p1),abs(p2)]) > barrier
        #             return forward_hits
        #         else
        #             return -1
        #         end
        if sign_sol < 0 && diff<1e-3  && maximum([abs(q1),abs(q2),abs(p1),abs(p2)]) > barrier
            return forward_hits
        else
            return -1
        end
        #     prob_b = HamiltonianProblem{true}(Hamiltonian_Dimer_Backwards, q0, p0, (0., t_end));
        #     #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        #     # sol_b=solve(prob_b,  RK4(), maxiters = 1e20,callback=cb_backward,save_start=false,save_end=true,save_everystep=false)
        #     sol_b=solve(prob_b,  RK4(), maxiters = 1e20,reltol=1e-8,abstol=1e-10, callback=cb_backward,save_start=false,save_end=true,save_everystep=false)
        #     # sol_b=solve(prob_b, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_backward,save_start=false,save_end=true,save_everystep=false)

        #     uf = zeros(6)
        #     uf[1] = sol_b[1,end]
        #     uf[2] = sol_b[2,end]
        #     uf[3] = sol_b[3,end]
        #     uf[4] = sol_b[4,end]
        #     uf[5] = sol_b[5,end]
        #     uf[6] = sol_b[6,end]
        #     q1 = uf[1]
        #     q2 = uf[2]
        #     p1 = uf[4]
        #     p2 = uf[5]
        #     q = [q1, q2]
        #     p = [p1, p2]

        #     backward_hits = uf[3]
        #     diff = abs(Hamiltonian_Dimer(q, p, 0) - H)
        #     if diff < 1e-3 && maximum([abs(q1),abs(q2),abs(p1),abs(p2)]) > barrier
        #         return forward_hits+ backward_hits# sol_f.t, sol_b.t
        #     else
        #         return -1
        #     end
        # else
        #     return -1
        # end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    else
        return -2
    end

end
end


text:
function DF1(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = abs(sin(0.5*pi*t))
    H = 0.75*sin(0.5*pi*t)+1.25
    g = 1+sum((x[2:end]-G).^2)
    f1 = x[1]
    f2 = g*(1-(x(1)/g)^H)
    [f1, f2]
end

function DF1_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    G = 0.75*sin(0.5*pi*t)
    H = G+1.25
    f1 = x
    f2 = g*(1-(x/g).^H)
    get_PF([f1, f2], false)
end

function DF2(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = abs(sin(0.5*pi*t))
    r = 1+floor((n-1)*G)
    tmp = setdiff(1:n,r)
    g = 1+sum((x(tmp)-G).^2)
    f1 = x[r]
    f2 = g*(1-(x(r)/g)^0.5)x
    [f1, f2]
end

function DF2_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    g = abs(sin(0.5*pi*t))
    f1 = x
    f2 = g*(1-(x/g).^0.5)
    get_PF([f1, f2], false)
end

function DF3(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = sin(0.5*π*t)
    H = g+1.5
    g = 1+sum((x[2:end]-G-x[1]^H).^2)
    f1 = x[1]
    f2 = g*(1-(x(1)/g)^H)
    [f1, f2]
end

function DF3_PS(tau, taut, nt)
    T0 = 50

    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    g = 1
    H = 50

    G = sin(0.5*pi*t)
    H = G+1.5

    x1 = LinRange(0, 1, 1500)
    x2 = G .+ x1.^H
    hcat(x1, x2)
end

function DF3_PF(tau, taut, nt)
    T0 = 50

    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    g = 1
    H = 50

    x = LinRange(0, 1, 1500)
    G = sin(0.5*pi*t)
    H = G+1.5
    f1 = x
    f2 = g*(1 .- (x/g).^H)
    get_PF([f1, f2], false)
end

function get_PF(f, nondominate)
    n = length(f)
    s = length(f[1])
    h = zeros(s, n)
    for i=1:n
        fi = reshape(f[i], s, 1)
        h[:,i] = fi
    end
    h
end

function DF3_Plot()
    anim = @animate for t in 50:10:440
        ps = MOEA.DF3_PS(t, 10, 20)
        pf = MOEA.DF3_PF(t, 10, 20)
        p = Plots.plot(title="DF3", lw=2, layout=(1,2), size=(850, 400))
        scatter!(p[1],
                 ps[:,1], ps[:,2],
                 seriestype = :scatter,
                 markersize=2,
                 color="green",
                 title = "Pareto optimal solutions (t = $(lpad(t,3,"0")))",
                 legend = false,
                 xlims = (0, 1),
                 ylims = (-1, 2)
                 )
        scatter!(p[2],
                 pf[:,1], pf[:,2],
                 seriestype = :scatter,
                 markersize=2,
                 color="blue",
                 title = "true Pareto front (t = $(lpad(t,3,"0")))",
                 legend = false,
                 xlims = (0, 1),
                 ylims = (0, 1)
                 )
    end
    gif(anim, "DF3.gif", fps = 20, loop = 0)
end

function DF3_PS_A()
    anim = @animate for t in 50:10:440
        ps = MOEA.DF3_PS(t, 10, 10)
        fig = Plots.plot(
            ps[:,1], ps[:,2],
            seriestype = :scatter,
            markersize=2,
            color=:blue,
            title = "Time $(t)",
            legend = false,
        )
        xlims!((0, 1))
        ylims!((0, 2))
    end
    gif(anim, "DF3_PS.gif", fps = 20, loop = 0)
end

function DF3_PF_A()
    anim = @animate for t in 50:10:440
        pf = MOEA.DF3_PF(t, 10, 10)
        fig = Plots.plot(
            pf[:,1], pf[:,2],
            seriestype = :scatter,
            markersize=2,
            color=:blue,
            title = "Time $(t)",
            legend = false,
        )
    end
    gif(anim, "DF3_PF.gif", fps = 20, loop = 0)
end

function DF4(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    a = sin(0.5*pi*t)
    b = 1+abs(cos(0.5*pi*t))
    c = max(abs(a), a+b)
    H = 1.5+a
    g = 0 # 1+sum((x[2:end]-a*(x[1]/c).^(2./[2:n])).^2)
    f1 = g*abs(x(1)-a).^H
    f2 = g*abs(x(1)-a-b).^H
    [f1, f2]
end

function DF4_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    a = sin(0.5*pi*t)
    b = 1+abs(cos(0.5*pi*t))
    x = linspace(a,a+b)
    H = 1.5+a
    f1 = g*abs(x-a).^H
    f2 = g*abs(x-a-b).^H
    get_PF([f1, f2], false)
end

function DF5(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = sin(0.5*pi*t)
    w = floor(10*G)
    g = 1+sum((x[2:end]-G).^2)
    f1 = g*(x(1)+0.02*sin(w*pi*x(1)))
    f2 = g*(1-x(1)+0.02*sin(w*pi*x(1)))
    [f1, f2]
end

function DF5_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = sin(0.5*pi*t)
    w = floor(10*G)
    f1 = g*(x+0.02*sin(w*pi*x))
    f2 = g*(1-x+0.02*sin(w*pi*x))
    get_PF([f1, f2], false)
end

function DF6(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G = sin(0.5*pi*t)
    a = 0.2+2.8*abs(G)
    y = x[2:end]-G
    g = 1+sum((abs(G)*y.^2-10*cos(2*pi*y)+10))
    f1 = g*(x1+0.1*sin(3*pi*x1)).^a
    f2 = g*(1-x1+0.1*sin(3*pi*x1)).^a
    [f1, f2]
end

function DF6_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    a=0.2+2.8*abs(G)
    f1=g*(x+0.1*sin(3*pi*x)).^a
    f2=g*(1-x+0.1*sin(3*pi*x)).^a
    get_PF([f1, f2], false)
end

function DF7(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    a=5*cos(0.5*pi*t)
    tmp=1/(1+exp(a*(x(1)-2.5)))
    g=1+sum((x[2:end]-tmp).^2)
    f1=g*(1+t)/x(1)
    f2=g*x(1)/(1+t)
    [f1, f2]
end

function DF7_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    f1=g*(1+t)./x
    f2=g*x/(1+t)
    get_PF([f1, f2], false)
end

function DF8(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    a=2.25+2*cos(2*pi*t)
    tmp=G*sin(4*pi*x(1))/(1+abs(G))
    g=1+sum((x[2:end]-tmp).^2)
    f1=g*(x(1)+0.1*sin(3*pi*x(1)))
    f2=g*(1-x(1)+0.1*sin(3*pi*x(1))).^a
    [f1, f2]
end

function DF8_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    a=2.25+2*cos(2*pi*t)
    f1=g*(x+0.1*sin(3*pi*x))
    f2=g*(1-x+0.1*sin(3*pi*x)).^a
    get_PF([f1, f2], false)
end

function DF9(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    N=1+floor(10*abs(sin(0.5*pi*t)))
    g=1
    for i=2:n
        tmp=x(i)-cos(4*t+x(1)+x(i-1))
        g=g+tmp.^2
    end
    f1=g*(x(1)+max(0, (0.1+0.5/N)*sin(2*N*pi*x(1))))
    f2=g*(1-x(1)+max(0, (0.1+0.5/N)*sin(2*N*pi*x(1))))
    [f1, f2]
end

function DF9_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x = LinRange(0, 1, 1500)
    N=1+floor(10*abs(sin(0.5*pi*t)));
    f1=g*(x+max(0, (0.1+0.5/N)*sin(2*N*pi*x)));
    f2=g*(1-x+max(0, (0.1+0.5/N)*sin(2*N*pi*x)));
    get_PF([f1, f2], false)
end

function DF10(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    H=2.25+2*cos(0.5*pi*t)
    tmp=sin(4*pi*(x(1)+x(2)))/(1+abs(G))
    g=1+sum((x[3:end]-tmp).^2)
    f1=g*sin(0.5*pi*x(1)).^H
    f2=g*sin(0.5*pi*x(2)).^H.*cos(0.5*pi*x(1)).^H
    f3=g*cos(0.5*pi*x(2)).^H.*cos(0.5*pi*x(1)).^H
    [f1, f2, f3]
end

function DF10_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x1,x2=meshgrid(linspace(0,1,H))
    H=2.25+2*cos(0.5*pi*t)
    f1=g*sin(0.5*pi*x1).^H
    f2=g*sin(0.5*pi*x2).^H.*cos(0.5*pi*x1).^H
    f3=g*cos(0.5*pi*x2).^H.*cos(0.5*pi*x1).^H
    get_PF([f1, f2, f3], false)
end

function DF11(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=abs(sin(0.5*pi*t))
    g=1+G+sum((x[3:end]-0.5*G*x(1)).^2)
    y=pi*G/6+(pi/2-pi*G/3)*x(1:2)
    f1=g*sin(y(1))
    f2=g*sin(y(2))*cos(y(1))
    f3=g*cos(y(2))*cos(y(1))
    [f1, f2]
end

function DF11_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=abs(sin(0.5*pi*t))
    y1=pi*G/6+(pi/2-pi*G/3)*x1
    y2=pi*G/6+(pi/2-pi*G/3)*x2
    f1=g.*sin(y1)
    f2=g.*sin(y2)*cos(y1)
    f3=g.*cos(y2)*cos(y1)
    get_PF([f1, f2, f3], false)
end

function DF12(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    k=10*sin(pi*t)
    tmp1=x[3:end]-sin(t*x(1))
    tmp2=abs(sin(floor(k*(2*x(1:2)-1))*pi/2))
    g=1+sum(tmp1.^2)+prod(tmp2)
    f1=g*cos(0.5*pi*x(2)).*cos(0.5*pi*x(1))
    f2=g*sin(0.5*pi*x(2)).*cos(0.5*pi*x(1))
    f3=g*sin(0.5*pi*x(1))
    [f1, f2, f3]
end

function DF12_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x1,x2=meshgrid(linspace(0,1,H))
    k=10*sin(pi*t)
    tmp2=abs(sin(floor(k*(2*x1-1))*pi/2).*sin(floor(k*(2*x2-1))*pi/2))
    g=1+tmp2
    f1=g.*cos(0.5*pi*x2).*cos(0.5*pi*x1)
    f2=g.*sin(0.5*pi*x2).*cos(0.5*pi*x1)
    f3=g.*sin(0.5*pi*x1)
    get_PF([f1, f2, f3], false)
end

function DF13(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    p=floor(6*G)
    g=1+sum((x[3:end]-G).^2)
    f1=g*cos(0.5*pi*x(1)).^2
    f2=g*cos(0.5*pi*x(2)).^2
    f3=g*sin(0.5*pi*x(1)).^2+sin(0.5*pi*x(1)).*cos(p*pi*x(1)).^2 +
        sin(0.5*pi*x(2)).^2+sin(0.5*pi*x(2)).*cos(p*pi*x(2)).^2
    [f1, f2, f3]
end

function DF13_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    p=floor(6*G)
    f1=g.*cos(0.5*pi*x1).^2
    f2=g.*cos(0.5*pi*x2).^2
    f3=g.*sin(0.5*pi*x1).^2+sin(0.5*pi*x1).*cos(p*pi*x1).^2 +
        sin(0.5*pi*x2).^2+sin(0.5*pi*x2).*cos(p*pi*x2).^2
    get_PF([f1, f2, f3], false)
end

function DF14(x, tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    G=sin(0.5*pi*t)
    g=1+sum((x[3:end]-G).^2)
    y=0.5+G*(x(1)-0.5)
    f1=g*(1-y+0.05*sin(6*pi*y))
    f2=g*(1-x(2)+0.05*sin(6*pi*x(2)))*(y+0.05*sin(6*pi*y))
    f3=g*(x(2)+0.05*sin(6*pi*x(2)))*(y+0.05*sin(6*pi*y))
    [f1, f2, f3]
end

function DF14_PF(tau, taut, nt)
    T0 = 50
    tau_tmp = max(tau+taut-(T0+1), 0)
    t = 1/nt*floor(tau_tmp/taut)

    n = length(x)

    x1,x2=meshgrid(linspace(0,1,H))
    G=sin(0.5*pi*t)
    y=0.5+G*(x1-0.5)
    f1=g.*(1-y+0.05*sin(6*pi*y))
    f2=g.*(1-x2+0.05*sin(6*pi*x2)).*(y+0.05*sin(6*pi*y))
    f3=g.*(x2+0.05*sin(6*pi*x2)).*(y+0.05*sin(6*pi*y))
    get_PF([f1, f2, f3], false)
end


text:
using Distributed

if nworkers()==1
	addprocs(2);
end

@everywhere begin
	using jInv.Mesh
	using jInv.LinearSolvers
	using Multigrid.ParallelJuliaSolver
	using Multigrid.DomainDecomposition
	using Multigrid
	using KrylovMethods

	using SparseArrays
	using LinearAlgebra
	using Test
	using jInv.Mesh;
end
@everywhere begin
	include(string(pwd(),"/DomainDecomposition/DDPoissonFuncs.jl"));
end

println("************************************************* Dirichlet DomainDecomposition for Poisson 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [128,128];
Mr = getRegularMesh(domain,n)
G = getNodalGradientMatrix(Mr);
Ar = G'*G;
Ar = Ar + 1e-5*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

NumCells = [4,4];
overlap = [1,1];
n_tup = tuple(n...);
dim = length(n);

numCores 	= 2;
maxIter     = 10;
relativeTol = 1e-6;

q = randn(prod(n.+1));
x = zeros(size(q));
Ainv = getParallelJuliaSolver(Float64,Int64,numCores=numCores,backend=1);

DDparam = getDomainDecompositionParam(Mr,NumCells,overlap,getNodalIndicesOfCell,Ainv);

println("workers are: ",workers())

setupDDParallel(Ar,DDparam,workers());

Prec = r->solveDDParallel(Ar,r,zeros(eltype(q),size(q)),DDparam,workers(),1,0)[1];
x = zeros(size(q));

x, flag,rnorm,iter = KrylovMethods.fgmres(getAfun(Ar,zeros(eltype(q),size(q)),4),q,10,tol = 1e-10,maxIter = 2,M = Prec, x = x,out=2,flexible=true);

println(norm(Ar'*x - q)/norm(q))

println("************************************************* Neumann DomainDecomposition for Poisson 2D ******************************************************");

Ctor = DomainDecompositionOperatorConstructor(Mr,getSubParams,getLap,getDirichletMassNodalMesh);

setupDDParallel(Ctor,DDparam,workers());

Prec = r->solveDDParallel(Ar,r,zeros(eltype(q),size(q)),DDparam,workers(),1,0)[1];
x = zeros(size(q));

x, flag,rnorm,iter = KrylovMethods.fgmres(getAfun(Ar,zeros(eltype(q),size(q)),4),q,10,tol = 1e-10,maxIter = 2,M = Prec, x = x,out=2,flexible=true);


text:
using Test
prob = Problems.DubinsCar(:three_obstacles)[1]
Z0 = copy(prob.Z)
X0 = deepcopy(states(Z0))
U0 = deepcopy(controls(Z0))

# Solve with iLQR
ilqr = iLQRSolver(Problems.DubinsCar(:three_obstacles)...)
Z0 = copy(ilqr.Z)
X0 = deepcopy(states(Z0))
U0 = deepcopy(controls(Z0))
solve!(ilqr)
@test cost(ilqr) < 8.47
initial_controls!(ilqr, U0)
@test (@allocated solve!(ilqr)) == 0

# Solve with AL-iLQR
al = AugmentedLagrangianSolver(Problems.DubinsCar(:three_obstacles)...)
initial_controls!(al, U0)
solve!(al)
@test max_violation(al) < al.opts.constraint_tolerance
@test cost(al) < 10
initial_controls!(al, U0)
@test (@allocated solve!(al)) == 0

# Solve with ALTRO
altro = ALTROSolver(Problems.DubinsCar(:three_obstacles)...)
altro.opts.opts_pn.verbose = false
initial_controls!(altro, U0)
solve!(altro)
Z_sol = copy(get_trajectory(altro))
X_sol = deepcopy(states(altro))
@test max_violation(altro) < 1e-8
@test cost(altro) < 10
altro.opts.projected_newton = false
initial_controls!(altro, U0)
solve!(altro)
initial_controls!(altro, U0)
@test (@allocated solve!(altro)) == 0

# Ipopt
ipopt = DIRCOLSolver(Problems.DubinsCar(:three_obstacles)..., integration=HermiteSimpson)
ipopt.opts.verbose = false
solve!(ipopt)
@test max_violation(ipopt) < 2e-8
@test cost(ipopt) < 12


text:
# functions related to noncentral T distribution

# R implementations
using .RFunctions:
    ntdistpdf,
    ntdistlogpdf,
    ntdistcdf,
    ntdistccdf,
    ntdistlogcdf,
    ntdistlogccdf,
    ntdistinvcdf,
    ntdistinvccdf,
    ntdistinvlogcdf,
    ntdistinvlogccdf


text:
module Sys

export cleardir, memory

"""
clear a floder
"""
function cleardir(dir::AbstractString)
    for fname in readdir(dir)
        rm(joinpath(dir, fname), recursive=true)
    end
end

end # end of module


text:
const MOIT = MOI.Test

@testset "Linear tests" begin
    MOIT.contlineartest(diff_optimizer(GLPK.Optimizer), MOIT.TestConfig(basis = true), [
        "partial_start",  # see below
        "linear12",       # see below
    ])
    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.partial_start_test(
        model,
        MOIT.TestConfig(basis = true, optimal_status=MOI.LOCALLY_SOLVED, atol=ATOL, rtol=RTOL)
    )

    # This requires an infeasiblity certificate for a variable bound.
    MOIT.linear12test(
        diff_optimizer(GLPK.Optimizer),
        MOIT.TestConfig(infeas_certificates=false)
    )
end

@testset "Convex Quadratic tests" begin
    model = diff_optimizer(OSQP.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.qp1test(model, MOIT.TestConfig(atol=1e-2, rtol=1e-2))
    model = diff_optimizer(OSQP.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.qp2test(model, MOIT.TestConfig(atol=1e-2, rtol=1e-2))
    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.qp3test(
        model,
        MOIT.TestConfig(optimal_status=MOI.LOCALLY_SOLVED, atol=1e-3),
    )
end

@testset "FEASIBILITY_SENSE zeros objective" begin
    model = diff_optimizer(GLPK.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))

    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol=ATOL rtol=RTOL

    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) >= 1.0
end

@testset "ModelLike" begin
    for opt in [GLPK.Optimizer]
        MODEL = diff_optimizer(opt)
        @testset "default_objective_test" begin
            MOIT.default_objective_test(MODEL)
        end
        @testset "default_status_test" begin
            MOIT.default_status_test(MODEL)
        end
        @testset "nametest" begin
            MOIT.nametest(MODEL)
        end
        @testset "validtest" begin
            MOIT.validtest(MODEL)
        end
        @testset "emptytest" begin
            MOIT.emptytest(MODEL)
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(MODEL)
        end
        @testset "copytest" begin
            # Requires VectorOfVariables
            # MOIT.copytest(MODEL, MOIU.CachingOptimizer(
            #     diff_optimizer(GLPK.Optimizer),
            #     GLPK.Optimizer()
            # ))
        end
    end
end

@testset "Unit" begin
    MOIT.unittest(diff_optimizer(GLPK.Optimizer), MOIT.TestConfig(), [
        "number_threads", # might not work on all solvers

        # not testing integer constraints
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3",
        "solve_integer_edge_cases",

        "delete_soc_variables",
        "solve_qcp_edge_cases",  # currently only affine or conic constraints
        "solve_objbound_edge_cases",
        "solve_qp_edge_cases",  # No quadratics
        "solve_qp_zero_offdiag",
        "update_dimension_nonnegative_variables", # TODO: fix this

        # see below
        "solve_duplicate_terms_vector_affine"
    ])
    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.solve_duplicate_terms_obj(model, MOIT.TestConfig())
end

@testset "basic_constraint_tests" begin
    # it contains SOCP constraints
    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    MOIT.basic_constraint_tests(model, MOIT.TestConfig())
end

# TODO: re-organiza conic tests
@testset "contconic.jl tests" begin
    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    # linear tests
    for (name, test) in MOI.Test.lintests
        test(model, MOIT.TestConfig())
    end

    CONFIG_LOW_TOL = MOIT.TestConfig(atol = 1e-3, rtol = 1e-2, duals = false, infeas_certificates = false)
    # SOCP tests
    for (name, test) in MOI.Test.soctests
        test(model, CONFIG_LOW_TOL)
    end
end


text:
using AbstractPlotting, CairoMakie, MakieRecipes; nothing# hide
using MakieRecipes.RecipesPipeline; using MakieRecipes.RecipesPipeline: PlotUtils; using MakieRecipes.Colors: @colorant_str; nothing # hide
# ```julia
# using Makie, MakieRecipes
# ```

# ## The simplest example model
using MakieRecipes.RecipesBase

struct T end

RecipesBase.@recipe function plot(::T, n = 1; customcolor = :green)
    markershape --> :auto        # if markershape is unset, make it :auto
    markercolor :=  customcolor  # force markercolor to be customcolor
    xrotation   --> 45           # if xrotation is unset, make it 45
    zrotation   --> 90           # if zrotation is unset, make it 90
    rand(10,n)                   # return the arguments (input data) for the next recipe
end

recipeplot(T(); seriestype = :path)

AbstractPlotting.save("basic.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](basic.svg)

# ## Testing out series decomposition

sc = Scene()
recipeplot!(sc, rand(10, 2); seriestype = :scatter)
recipeplot!(sc, 1:10, rand(10, 1); seriestype = :path)

AbstractPlotting.save("series.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](series.svg)

# ## Differential Equations

using OrdinaryDiffEq, StochasticDiffEq, DiffEqNoiseProcess

# ### A simple exponential growth model

f(u,p,t) = 1.01.*u
u0 = [1/2, 1]
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

recipeplot(sol)

AbstractPlotting.save("exp.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](exp.svg)

# ### Matrix DiffEq

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

recipeplot(sol)

AbstractPlotting.save("mat.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](mat.svg)

# ### Stochastic DiffEq

f(du,u,p,t) = (du .= u)
g(du,u,p,t) = (du .= u)
u0 = rand(4,2)

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRIW1())

recipeplot(sol)

AbstractPlotting.save("stochastic.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](stochastic.svg)

# ### Animating a differential equation solution

# You can update arguments and attributes arbitrarily, and the recipe pipeline is rerun.

function lorenz(du, u, p, t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

prob = ODEProblem(lorenz, [1., 5., 10.], (0., 100.), (10.0, 28.0, 8/3))

ρ = Node(28.0)

sol = @lift solve(remake(prob; p = (10.0, $ρ, 8/3)), Tsit5())

sc = recipeplot(sol)

ylims!(sc, -30, 70) # avoid jitter when animating

record(sc, "lorenz.gif", LinRange(0, 35, 100)) do ρᵢ
    ρ[] = ρᵢ
end

# ![](lorenz.gif)

# ## Phylogenetic tree

using Phylo

assetpath = joinpath(dirname(pathof(MakieRecipes)), "..", "docs", "src", "assets")
hummers = open(t -> parsenewick(t, NamedPolytomousTree), joinpath(assetpath, "hummingbirds.tree"))

evolve(tree) = Phylo.map_depthfirst((val, node) -> val + randn(), 0., tree, Float64)

trait = evolve(hummers)

scp = recipeplot(
    hummers;
    treetype = :fan,
    line_z = trait,
    linewidth = 5,
    showtips = false,
    cgrad = :RdYlBu,
    seriestype = :path,
    scale_plot = false, # Makie attributes can be used here as well!
    show_axis = false
)

AbstractPlotting.save("phylo.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](phylo.svg)

# ### Animation with different colormaps (changing attributes)

# You can update arguments and attributes arbitrarily, and the recipe pipeline is rerun.

record(scp, "phylo_colormaps.gif", PlotUtils.cgradients(:colorcet), framerate = 3) do cmap
    scp.plots[1].cgrad[] = cmap
end

# ![](phylo_colormaps.gif)

# ## GraphRecipes

using GraphRecipes
using GraphRecipes: GraphPlot

# In tis case, we don't need access to all of the functions which GraphRecipes provides - just the
# `GraphPlot` struct, which is used for recipe dispatch.  The rest is handled by the GraphRecipes internals.

# ### The simplest graph recipe
# Because `userplot` is unsupported, we have to use the low-level interface to user plot recipes,
# by wrapping our arguments in the userplot type.  Arguments must be supplied as a `Tuple` or some other iterable.

recipeplot(
    GraphPlot(
        (
        [0  1  1;
         1  0  1;
         1  1  0],
         )
    );
    show_axis = false,
    scale_plot = false
)

AbstractPlotting.save("simplegraph.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](simplegraph.svg)

# ### Undirected graph with node labels and sizes
n = 15
A = Float64[ rand() < 0.5 ? 0 : rand() for i=1:n, j=1:n]
for i=1:n
    A[i, 1:i-1] = A[1:i-1, i]
    A[i, i] = 0
end

recipeplot(
        GraphPlot((A,));
        markersize = 0.2,
        node_weights = 1:n,
        markercolor = range(colorant"yellow", stop=colorant"red", length=n),
        names = 1:n,
        fontsize = 10,
        linecolor = :darkgrey,
        show_axis = false,
        scale_plot = false
    )

AbstractPlotting.save("undigraph.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](undigraph.svg)

# ### Multigraphs

recipeplot(GraphPlot(([[1,1,2,2],[1,1,1],[1]],)); names="node_".*string.(1:3), nodeshape=:circle, self_edge_size=0.25, show_axis = false, scale_plot = false)

AbstractPlotting.save("digraph.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](digraph.svg)

# ### Julia AST

code = quote
    function mysum(list)
        out = 0
        for value in list
            out += value
        end
        out
    end
end

recipeplot(code; fontsize = 12, shorten = 0.01, axis_buffer = 0.15, nodeshape = :rect)
AbstractPlotting.save("ast.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](ast.svg)

# ### Type tree with GraphRecipes

recipeplot(AbstractFloat; method = :tree, fontsize = 10)
AbstractPlotting.save("typetree.svg", AbstractPlotting.current_scene()); nothing #hide
# ![](typetree.svg)


text:
using DataDeps

register(DataDep(
    "letter-recognition",
    """
    Dataset: letter-recognition
    Author: David J. Slate
    Odesta Corporation; 1890 Maple Ave; Suite 115; Evanston, IL 60201
    Website: https://archive.ics.uci.edu/ml/datasets/Letter+Recognition

    Observations: 20000
    Features:     16
    Classes:      26
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/letter-recognition/letter-recognition.data",
    "2b89f3602cf768d3c8355267d2f13f2417809e101fc2b5ceee10db19a60de6e2",
    post_fetch_method = (path -> begin
        MVDA.process_dataset(path,
            header=false,
            target_index=1,
            feature_indices=2:17,
            ext=".csv.gz",
        )
    end),
))


text:
import FastTransforms
import QuadratureRules: shift_nodes

@testset "$(rpad("Lobatto-Chebyshev",80))" begin

    @test_throws ErrorException LobattoChebyshevQuadrature(1)

    for s in 2:10
        @test LobattoChebyshevQuadrature(s) == ChebyshevQuadrature(s, 2)
        @test LobattoChebyshevQuadrature(s) == ChebyshevQuadrature(Float64, s, Val(2))

        @test sum(weights(LobattoChebyshevQuadrature(s))) ≈ 1

        c = shift_nodes(reverse(FastTransforms.chebyshevpoints(Float64, s, Val(2))))

        @test chebyshev_points(Float64, s, Val(2)) ≈ c
        @test nodes(LobattoChebyshevQuadrature(s)) ≈ c

        # TODO test weights

    end

end


text:
##
# Utility methods
#########################################################
const _cache = Dict()
cache(k, v) = _cache[k] = v
cache_clear(k) = delete!(_cacle, k)
cache_clear() = empty!(_cache)
function cached_get(k, fn)
    v = get(_cache, k, nothing)
    (nothing != v) && return v
    cache(k, fn())
end

openable(path::File) = path
openable(path::String) = File(path)

# cache deserialized objs to help the simple search implementation speed up by preventing repeated loading of indices
as_deserialized(f::File) = as_deserialized(f.path)
function as_deserialized(path::String)
    io = as_io(path)
    obj = deserialize(io)
    close(io)
    obj
end

function preprocess(entity::Union(StringDocument,Corpus))
    #prepare!(entity, strip_corrupt_utf8 | strip_case)
    prepare!(entity, strip_case)
    remove_html_tags!(entity)
    prepare!(entity, strip_whitespace | strip_non_letters | strip_articles | strip_prepositions | strip_pronouns | strip_stopwords)
    stem!(entity)

    if isa(entity, Corpus)
        standardize!(entity, TokenDocument)
    else
        entity = convert(TokenDocument, entity)
    end
    entity
end

as_serialized(obj, f::File) = as_serialized(obj, f.path)
function as_serialized(obj, path::String)
    iob = IOBuffer()
    serialize(iob, obj)

    io = open(path, "w")
    write(io, takebuf_array(iob))
    close(io)
    close(iob)
    path
end

function setup_folders(clean::Bool=false, everywhere::Bool=true)
    if everywhere
        for p in procs()
            remotecall_fetch(p, setup_folders, clean, false)
        end
    else
        println("Setting up work folders...")
        clean && isdir(fs_pfx) && rm(fs_pfx; recursive=true)
        for d in [fs_pfx, part_idx_location, doc_to_id_location, id_to_doc_location, stat_location, docs_location]
            isdir(d) || mkdir(d)
        end
        nothing
    end
end


text:
module IanSmithdistmath

export remap

include("utilities.jl")

end


text:
#

using MultithreadingBenchmarks

function work(i, v, n)
    out = v
    for i in 1:n
        out *= v
    end
    out
end

@noinline mul_barrier(a,b) = a*b
@noinline mul_barrier(a::Vector,b) = a[1]*b
@noinline mul_barrier(a,b::Vector) = a*b[1]

# Allocate in a loop until anything is pushed into `signal` Channel.
function allocate_in_background(signal::Channel, v)
    # Since there is
    ALLOC_EVERY_N = 2
    GC_EVERY_N = 1024
    out = v
    i = 0
    while !isready(signal)
        out_temp = if (i % ALLOC_EVERY_N == 0)
            vi = rand([[v], v])
            typeof(v)(mul_barrier(out, v))
        else
            out * v
        end
        out_temp::typeof(v)
        out = out_temp
        i += 1

        # Manually trigger the GC every so often, since it wasn't running on its own
        if (i % GC_EVERY_N == 0)
            GC.gc()
        end

        # TODO: But if we add this yield-point does it occupy the thread correctly?
        #yield()  # We need this yield-point so this "background" task doesn't keep julia alive.
    end
    return out
end

# ----------------------------
# Start a Task that will take over one thread and allocate garbage in the background
signal = Channel(1)
t = Threads.@spawn allocate_in_background(signal, 2)
# ----------------------------

MultithreadingBenchmarks.measure_work(work)

# ------------------------
put!(signal, 0)
wait(t)
# ------------------------


text:
immutable Node{T}
    state::T
    score::Float64
    prev::Node{T}

    Node(state) = new(state, 0.0)
end

lessthan{T}(x::Node{T}, y::Node{T}) = x.score > y.score

function getseq{T}(finalstate::T)
    seq = T[]
    s = finalstate
    while true
        unshift!(seq, s)
        isdefined(s,:prev) || break
        s = s.prev
    end
    unshift!(seq, s)
    seq
end

"""
* next: state -> score, state
"""
function beamsearch{T}(initstate::T, beamsize::Int)
    chart = Vector{T}[]
    push!(chart, [initstate])

    k = 1
    while k <= length(chart)
        states = chart[k]
        length(states) > beamsize && sort!(states, lt=lessthan)
        for i = 1:min(beamsize, length(states))
            for (s,score) in next(states[i])
                while s.step > length(chart)
                    push!(chart, T[])
                end
                push!(chart[s.step], s)
            end
        end
        k += 1
    end

    state = initstate
    k = 1
    while true
        nexts = next(state)

    end
end

function beamsearch{T}(initstate::T, beamsize::Int)
    chart = Vector{T}[]
    push!(chart, [initstate])

    k = 1
    while k <= length(chart)
        states = chart[k]
        length(states) > beamsize && sort!(states, lt=lessthan)
        for i = 1:min(beamsize, length(states))
            for (s,score) in next(states[i])
                while s.step > length(chart)
                    push!(chart, T[])
                end
                push!(chart[s.step], s)
            end
        end
        k += 1
    end
    sort!(chart[end], lt=lessthan)
    chart
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
#-------------------------------------------------------------------------------------------------
#--
#-- Linear interpolation of Y values for a new set of X values.
#--
function linearInterpolation(x, xvar, yvar)
    l = length(xvar)

    # If l = 1, ratio will be the only one
    if l == 1
        return yvar[1]
    else
        for i in 2:l
            x1 = xvar[i-1]
            x2 = xvar[i  ]

            if x < x2
                deltaY = yvar[i] - yvar[i-1]
                return yvar[i-1] + deltaY * (x - x1) / (x2 - x1)
            end
        end

        # Last possible choice
        return yvar[l]
    end
end


text:
using Test
using StaticArrays
using Polyhedra

function test_planar_square(VT, d)
    ei = ntuple(i -> Polyhedra.basis(VT, d, i), 2)
    so = -ei[1] + ei[2]
    no = -ei[1] - ei[2]
    ne = ei[1] - ei[2]
    se = ei[1] + ei[2]
    expected = [no, so, se, ne]
    for ps in [
        [no, so, ne, se],
        [no, se, so, ne],
        [se, so, no, ne],
        expected
    ]
        hull = Polyhedra.planar_hull(vrep(ps))
        @test collect(points(hull)) == expected
        @test !hasallrays(hull)
        v = vrep(ps) + conichull(se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne, no, so]
        @test collect(rays(hull)) == [Ray(se)]
        @test !haslines(hull)
        v = vrep(ps) + conichull(se, se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne, no, so]
        @test collect(rays(hull)) == [Ray(se)]
        @test !haslines(hull)
        ne = ei[1] - ei[2]
        for v in [vrep(ps) + conichull(se, ne),
                  vrep(ps) + conichull(se, ne, ei[1])]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no, so]
            @test collect(rays(hull)) == [Ray(ne), Ray(se)]
            @test !haslines(hull)
        end
        v = vrep(ps) + conichull(se, ne, -se)
        hull = Polyhedra.planar_hull(v)
        @test collect(points(hull)) == [ne]
        @test collect(rays(hull)) == [Ray(ne)]
        @test collect(lines(hull)) == [Line(se)]
        for v in [vrep(ps) + conichull(se, ne, -ne),
                  vrep(ps) + conichull(se, Line(ne)),
                  vrep(ps) + conichull(ne, se, -ne)]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no]
            @test collect(rays(hull)) == [Ray(se)]
            @test collect(lines(hull)) == [Line(ne)]
        end
        for v in [vrep(ps) + conichull(se, ne, -ne, -se),
                  vrep(ps) + conichull(se, ne, -ei[1])]
            hull = Polyhedra.planar_hull(v)
            @test collect(points(hull)) == [no]
            @test !hasrays(hull)
            @test collect(lines(hull)) == [Line(ne), Line(se)]
        end
    end
end

# Square + redundant points + [+-1.5, -+0.25]
function test_planar_square_with_more()
    expected = [[-1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [1.0, -1.0]]
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0 + 1e-10])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    v0 = convexhull([0.0, 1.0 + 1e-10], [-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    expected = [[-1.0, -1.0], [-1.0, 1.0], [0.0, 1.0 + 1e-4], [1.0, 1.0], [1.0, -1.0]]
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0 + 1e-4])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    expected = [[-1.0, -1.0], [-1.5, 0.25], [-1.0, 1.0], [1.0, 1.0], [1.5, -0.25], [1.0, -1.0]]
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0], [-1.5, 0.25], [1.5, -0.25])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, 1.0], [-0.5, -1.0], [-1.5, 0.25], [0.0, -0.75], [0.0, 0.75], [1.5, -0.25])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.0, -1.0], [0.0, 1.0], [-0.5, -1.0], [-1.5, 0.25], [0.0, -0.75], [0.0, 0.75], [1.5, -0.25])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
    v0 = convexhull([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0], [1.0, 1.0], [0.5, 1.0], [0.0, -1.0], [0.0, 1.0], [-0.5, -1.0], [-1.5, 0.25], [0.0, -0.75], [0.0, 0.75], [1.5, -0.25])
    v1 = Polyhedra.planar_hull(v0)
    @test collect(points(v1)) == expected
end

@testset "Planar" begin
    for T in [Float64, Int]
        for (VT, d) in [(Vector{T}, 2), (SVector{2, T}, StaticArrays.Size(2))]
            test_planar_square(VT, d)
        end
    end
    test_planar_square_with_more()
end


text:
# You might not want to run this file completely, as the Pkg-commands can take a
# long time to complete.

# list all available packages:
Pkg.available()

# install one package (e.g. [Calculus](https://github.com/johnmyleswhite/Calculus.jl)) and all its dependencies:
Pkg.add("Calculus")

# to list all installed packages
Pkg.installed()

# to update all packages to their newest version
Pkg.update()

# to use a package:
using Calculus
# will import all functions of that package into the current namespace, so that
# it is possible to call
derivative(x -> sin(x), 1.0)
# without specifing the package it is included in.

import Calculus
# will enable you to specify which package the function is called from
Calculus.derivative(x -> cos(x), 1.0)

# Using `import` is especially useful if there are conflicts in function/type-names
# between packages.
# Example:
# Winston as well as Gadfly provide a plot() function (see below).
Pkg.add("Winston")
Pkg.add("Gadfly")

# If you were to "import" both of the packages with `using`, there would be a conflict.
# That can be prevented by using `import`, as follows:
import Winston
import Gadfly

Winston.plot(rand(4))
Gadfly.plot(x=[1:10], y=rand(10))

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
"""
Main module for `Helium.jl` -- a binary format writer and reader for matrix.
Three functions are exported from this module for public use:
- [`csv2he`](@ref). Convert a CSV file containing a matrix to the binary Helium format.
- [`writehe`](@ref). Write matrix in Helium format.
- [`readhe`](@ref). Read file in Helium format.

"""

module Helium

using DelimitedFiles

include("he.jl")
export HeAttributes
export HeInfo

include("writehe.jl")
export writehe

include("readhe.jl")
export readhe, getcolnames, getrownames, getsupp, readheader

include("csv2helium.jl")
export csv2he, csv2mat

include("xtract.jl")
export getskipmat, getmat

include("helium2csv.jl")
export he2csv

end # module


text:
import Base.FloatRange

type CoordInterpGrid{T<:FloatingPoint, N, BC<:BoundaryCondition, IT<:InterpType, R} <: AbstractInterpGrid{T,N}
    coord::R
    grid::InterpGrid{T,N,BC,IT}
    function CoordInterpGrid(coord::NTuple{N,Range}, grid::InterpGrid{T,N,BC,IT})
        map(length,coord) == size(grid) || throw(DimensionMismatch("Coordinate lengths do not match grid size."))
        new(coord,grid)
    end
end

function CoordInterpGrid{T,N,BC,IT}(coord::NTuple{N,Range}, grid::InterpGrid{T,N,BC,IT})
    CoordInterpGrid{T,N,BC,IT,typeof(coord)}(coord,grid)
end

function CoordInterpGrid{T<:FloatingPoint,R<:Range}(coord::R, grid::InterpGrid{T,1})
    CoordInterpGrid((coord,),grid)
end

function CoordInterpGrid{N,T<:FloatingPoint}(coord::NTuple{N,Range}, A::Array{T,N}, args...)
    CoordInterpGrid(coord,InterpGrid(A,args...))
end
function CoordInterpGrid{R<:Range,T<:FloatingPoint}(coord::R,A::Array{T,1},args...)
    CoordInterpGrid((coord,),InterpGrid(A,args...))
end

size(C::CoordInterpGrid) = size(C.grid)

coordlookup(r::FloatRange,x::Real) = (r.divisor*x-r.start)/r.step + 1.0
coordlookup(r::StepRange,x::Real) = (x-r.start)/r.step + 1.0
coordlookup(r::UnitRange,x::Real) = x-r.start + 1
coordlookup{N,R<:Range,T<:Real}(r::NTuple{N,R},x::NTuple{N,T}) = map(coordlookup,r,x)

function getindex(C::CoordInterpGrid, x::Real...)
    getindex(C.grid,coordlookup(C.coord,x)...)
end

# gradients
# chain rule: dy/dx = (dy/du) / (dx/du)

coordiscale(r::FloatRange,x::Real) = x*r.divisor/r.step
coordiscale(r::StepRange,x::Real) = x/r.step
coordiscale(r::UnitRange,x::Real) = x

function valgrad{T,BC,IT,R}(C::CoordInterpGrid{T,1,BC,IT,R}, x::Real)
    val, g = valgrad(C.grid,coordlookup(C.coord[1],x))
    val, coordiscale(C.coord[1],g)
end

function valgrad{T,N,BC,IT,R}(g::Vector{T}, C::CoordInterpGrid{T,N,BC,IT,R}, x::Real...)
    val = valgrad(g,C.grid,coordlookup(C.coord,x)...)
    for i=1:length(g)
        g[i] = coordiscale(C.coord[i],g[i])
    end
    val
end