text:
module InfiniteOpt

# Import the necessary packages.
import JuMP
import MathOptInterface
import Distributions
import Random
import MutableArithmetics
const _MA = MutableArithmetics
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const JuMPC = JuMP.Containers

# Import all of the datatpyes, methods, macros, and definitions.
include("datatypes.jl")
include("infinite_sets.jl")
include("parameters.jl")
include("variables.jl")
include("reduced_variables.jl")
include("expressions.jl")
include("mutable_arithmetics.jl")
include("measures.jl")
include("operators.jl")
include("constraints.jl")
include("objective.jl")
include("measure_expansions.jl")
include("macros.jl")
include("optimize.jl")
include("results.jl")
include("show.jl")
include("utilities.jl")

include("TranscriptionOpt/TranscriptionOpt.jl")
using .TranscriptionOpt

include("MeasureEvalMethods/MeasureEvalMethods.jl")
using .MeasureEvalMethods: generate_measure_data, mc_sampling, gauss_legendre,
gauss_hermite, gauss_laguerre, infinite_transform, generate_supports_and_coeffs,
trapezoid

# Export the measure eval functions and constants
export generate_measure_data, mc_sampling, gauss_legendre, gauss_hermite,
gauss_laguerre, trapezoid, infinite_transform, measure_dispatch,
generate_supports_and_coeffs

# Export model object datatype
export InfiniteModel

# Export macros
export @infinite_variable, @point_variable, @hold_variable, @infinite_parameter,
@BDconstraint, @finite_parameter, @set_parameter_bounds, @add_parameter_bounds

# Export constants
export Infinite, Hold, Point, Parameter, sampling, quadrature

# Export infinite parameter datatypes
export InfOptParameter, ParameterRef

# Export infinite parameter set types
export AbstractInfiniteSet, IntervalSet, DistributionSet

# Export infinite parameter methods
export build_parameter, add_parameter, infinite_set, set_infinite_set,
num_parameters, parameter_by_name, all_parameters, num_supports, has_supports,
set_supports, add_supports, delete_supports, supports, used_by_constraint,
used_by_measure, used_by_variable, is_used, group_id, is_independent,
set_independent, unset_independent, is_finite_parameter, supports_in_set,
generate_and_add_supports!, generate_support_values

# Export variable datatypes
export InfOptVariable, InfiniteVariable, PointVariable, HoldVariable,
GeneralVariableRef, InfiniteVariableRef, MeasureFiniteVariableRef,
FiniteVariableRef, PointVariableRef, HoldVariableRef, InfOptVariableRef,
ReducedInfiniteVariableRef, AbstractReducedInfo, ReducedInfiniteInfo,
ParameterBounds

# Export variable methods
export used_by_objective, infinite_variable_ref, parameter_refs,
set_parameter_refs, add_parameter_ref, used_by_point_variable, parameter_values,
eval_supports, used_by_reduced_variable, has_parameter_bounds, parameter_bounds,
set_parameter_bounds, add_parameter_bound, delete_parameter_bound,
delete_parameter_bounds

# Export expression datatypes
export InfiniteExpr, ParameterExpr, MeasureExpr

# Export constraint datatypes
export GeneralConstraintRef, InfiniteConstraintRef, MeasureConstraintRef,
FiniteConstraintRef, BoundedScalarConstraint

# Export measure datatypes
export AbstractMeasureData, DiscreteMeasureData, MultiDiscreteMeasureData,
Measure, MeasureRef

# Export measure methods
export add_measure, measure, measure_function, measure_data, expand,
expand_all_measures!, expect, support_sum, measure_name, measure_data_in_hold_bounds,
make_point_variable_ref, make_reduced_variable_ref, expand_measure, integral,
set_integral_defaults, integral_defaults

# Export transcription datatypes
export TranscriptionData, TranscriptionModel

# Export transcription methods
export is_transcription_model, transcription_data, transcription_variable,
transcription_constraint, transcription_model

# Export optimize methods
export optimizer_model, set_optimizer_model, optimizer_model_ready,
set_optimizer_model_ready, build_optimizer_model!, optimizer_model_key,
optimizer_model_variable, optimizer_model_constraint, variable_supports,
constraint_supports, constraint_parameter_refs, add_infinite_model_optimizer,
clear_optimizer_model_build!

# Export result query methods
export map_value, map_optimizer_index, map_dual, map_shadow_price,
map_lp_rhs_perturbation_range, map_lp_objective_perturbation_range

# Export functions that will be included in future JuMP releases

# Export printing methods
export bound_string

# Export support generation and fill-in methods
export fill_in_supports!, generate_supports

end # module


text:
using UHDBindings
import UHDBindings.LibUHD as LibUHD

import UHDBindings.initRxUHD
import UHDBindings.@assert_uhd

# ----------------------------------------------------
# --- Tools
# ----------------------------------------------------
"""
Init an empty string of size n, filled with space. Usefull to have container to get string from UHD.
"""
initEmptyString(n) = String(ones(UInt8,n)*UInt8(32))
"""
Retrict the big string used as log container to its usefull part
"""
truncate(s::String) = s[1 : findfirst("\0",s)[1] - 1]
"""
Interrupt the script just at right place
"""
stop() = error("Decide to stop the script. Be sure to release all ressources")
# ----------------------------------------------------
# --- Parameters
# ----------------------------------------------------
carrierFreq = 868e6
samplingRate= 2e6
gain        = 25
# uhdArgs     = "num_recv_frames=1000"
uhdArgs     = ""
nbAntennaRx = 2

# ----------------------------------------------------
# --- Address extraction
# ----------------------------------------------------
# allR  = uhd_find_devices()
# extractAddr(x) = x[1][findfirst("_addr=",x[1])[end] .+ (1:13)]
# addr  = extractAddr(allR)

# ----------------------------------------------------
# --- Create radio
# ----------------------------------------------------
# --- Open radio, in MIMO mode
# Note the warning about the streamer. It can not be set in MIMO mode
# The configuration will be updated later, and most of the parameters overwritten
# To ensure we have what we want, we need to bypass the streamer, and deactivate all antennas
radio = openUHD(carrierFreq,samplingRate,gain;nbAntennaRx=0,nbAntennaTx=0,bypassStreamer = true)

# ----------------------------------------------------
# --- Set up Radio
# ----------------------------------------------------
# --- Channels  & Antenna
channelIndexes = [0,1]
antennas       = ["TX/RX","RX2"]
# --- Streamer arguments
a1			   =  Base.unsafe_convert(Cstring,"fc32");
a2			   =  Base.unsafe_convert(Cstring,"sc16");
a3			   =  Base.unsafe_convert(Cstring,uhdArgs);
channel        = pointer(channelIndexes)
uhdArgs_0	   = LibUHD.uhd_stream_args_t(a1,a2,a3,channel,nbAntennaRx);
# --- RF Configuration
for (c,currChan) in enumerate(channelIndexes[1:nbAntennaRx])
    # Set the carrier frequencies
    updateCarrierFreq!(radio.rx,carrierFreq,currChan)
    # Set the sampling rate
    updateSamplingRate!(radio.rx,samplingRate,currChan)
    # Set the gains
    updateGain!(radio.rx,gain,currChan)
    # Set the antennas
    LibUHD.uhd_usrp_set_rx_antenna(radio.rx.uhd.pointerUSRP,antennas[c],currChan)
end
# --- Subdev with 2 boards
pointerSubDev = Ref{LibUHD.uhd_subdev_spec_handle}()
@assert_uhd LibUHD.uhd_subdev_spec_make(pointerSubDev,"A:0 A:1")
@assert_uhd LibUHD.uhd_usrp_set_rx_subdev_spec(radio.rx.uhd.pointerUSRP,pointerSubDev[],0)
LibUHD.uhd_subdev_spec_free(pointerSubDev)
# --- Internal streamer and buffer config
pointerArgs	  = Ref{LibUHD.uhd_stream_args_t}(uhdArgs_0);
pointerSamples = Ref{Csize_t}(0);
@assert_uhd LibUHD.uhd_usrp_get_rx_stream(radio.rx.uhd.pointerUSRP,pointerArgs,radio.rx.uhd.pointerStreamer)
@assert_uhd LibUHD.uhd_rx_streamer_max_num_samps(radio.rx.uhd.pointerStreamer,pointerSamples)
println("Internal buffer size is $(pointerSamples[])")
radio.rx.packetSize= pointerSamples[]
# --- Clock source
@assert_uhd LibUHD.uhd_usrp_set_clock_source(radio.rx.uhd.pointerUSRP,"internal",0)
@assert_uhd LibUHD.uhd_usrp_set_time_now(radio.rx.uhd.pointerUSRP,0,0,0)

# ----------------------------------------------------
# --- Set up streamer
# ----------------------------------------------------
streamCmd	= UHDBindings.uhd_stream_cmd_t(UHDBindings.UHD_STREAM_MODE_NUM_SAMPS_AND_MORE,radio.rx.packetSize*10_000,false,1,0.5);
# streamCmd	= UHDBindings.uhd_stream_cmd_t(UHDBindings.UHD_STREAM_MODE_NUM_SAMPS_AND_MORE,rx.packetSize,true,0,0);
pointerCmd	= Ref{UHDBindings.uhd_stream_cmd_t}(streamCmd);
LibUHD.uhd_rx_streamer_issue_stream_cmd(radio.rx.uhd.pointerStreamer,pointerCmd)
sleep(0.1)
# ----------------------------------------------------
# --- Julia buffers
# ----------------------------------------------------
# Define an array to get all the buffers from all the channels
nbSamples             = getBufferSize(radio)
sig                   = [zeros(Complex{Cfloat},radio.rx.packetSize) for n ∈ 1:nbAntennaRx]
listBuffer            = [pointer(sig[n],1) for n ∈ 1 : nbAntennaRx]
# ptr                   = Ref(Ptr{Cvoid}(listBuffer[1]))
ptr                   = listBuffer
pointerCounterSamples = Ref{Csize_t}(0);

# ----------------------------------------------------
# --- Receive data
# ----------------------------------------------------
pointerCounterSamples = Ref{Csize_t}(0);
LibUHD.uhd_rx_streamer_recv(radio.rx.uhd.pointerStreamer,ptr,nbSamples,radio.rx.uhd.addressMD,1.6,true,pointerCounterSamples)
nbReceivedSamples = pointerCounterSamples[]
println("Receive $nbReceivedSamples samples")

# ----------------------------------------------------
# --- Second call
# ----------------------------------------------------
accum = 0
# UHDBindings.restartStreamer(rx)
for iN = 1 : 1 : 10_000
    pointerCounterSamples = Ref{Csize_t}(0);
    LibUHD.uhd_rx_streamer_recv(radio.rx.uhd.pointerStreamer,ptr,nbSamples,radio.rx.uhd.addressMD,0,true,pointerCounterSamples)
    nbReceivedSamples = pointerCounterSamples[]
    global accum += nbReceivedSamples
end
println("End of transmission, received $accum samples")

# ----------------------------------------------------
# --- Logs
# ----------------------------------------------------
n = 100
container = initEmptyString(n)
# If we get 0, we should infvestigate metada structrue
LibUHD.uhd_rx_metadata_strerror(radio.rx.uhd.addressMD[],container,n)
@show truncate(container)

container = initEmptyString(n)
LibUHD.uhd_rx_streamer_last_error(radio.rx.uhd.pointerStreamer,container,n)
@show truncate(container)

# ----------------------------------------------------
# --- Close
# ----------------------------------------------------
close(radio)

ext: jl

max_stars_repo_name: RGerzaguet/UHDBindings

max_issues_repo_name: RGerzaguet/UHDBindings


text:
#
# Chunked Data Transfer
#

immutable ChunkedStream
    io::IO
end
function write_chunked(stream, arg)
    write(stream,string(hex(sizeof(arg)),CRLF))
    write(stream,arg)
    write(stream,string(CRLF))
end

write(io::ChunkedStream, arg) = write_chunked(io.io, arg)

#
# File uploads
#
# Upload a file using multipart form upload. `file` may be one of:
#
#   - An IO object whose contents will be uploaded
#   - A string or Array to be sent
#
# Note that when passing an IO object, the IO object may not otherwise be modified
# until the request completes. Optionally you may set `close` to true to have Requests
# automatically close your file when it's done with it.
#

immutable FileParam
    file::Union{IO,Base.File,AbstractString,Vector{UInt8}}     # The file
    # The content type (default: "", which is interpreted as text/plain serverside)
    ContentType::Compat.UTF8String
    name::Compat.UTF8String                     # The fieldname (in a form)
    filename::Compat.UTF8String                              # The filename (of the actual file)
    # Whether or not to close the file when the request is done
    close::Bool

    function FileParam(str::Union{AbstractString,Vector{UInt8}},ContentType="",name="",filename="")
        new(str,ContentType,name,filename,false)
    end

    function FileParam(io::IO,ContentType="",name="",filename="",close::Bool=false)
        new(io,ContentType,name,filename,close)
    end

    function FileParam(io::Base.File,ContentType="",name="",filename="",close::Bool=false)
        if !isopen(io)
            close = true
        end
        new(io,ContentType,name,filename,close)
    end
end

# Determine whether or not we need to use
datasize(::IO) = -1
datasize(f::Union{AbstractString,Array{UInt8}}) = sizeof(f)
datasize(f::File) = filesize(f)
datasize(f::IOBuffer) = nb_available(f)
function datasize(io::IOStream)
    iofd = fd(io)
    # If this IOStream is not backed by a file, we can't find the filesize
    if iofd == -1
        return -1
    else
        return filesize(iofd) - position(io)
    end
end

const multipart_mime = "multipart/form-data; boundary="
const part_mime = "Content-Disposition: form-data"
const name_file = "; name=\""
const filename_file = "; filename=\""
const ContentType_header = "Content-Type: "

function write_part_header(stream,file::FileParam,boundary)
    buf = IOBuffer()
    write(buf,"--",boundary,CRLF)
    write(buf,part_mime)
    !isempty(file.name) && write(buf,name_file,file.name,'\"')
    !isempty(file.filename) && write(buf,filename_file,file.filename,'\"')
    write(buf,CRLF)
    !isempty(file.ContentType) && write(buf,ContentType_header,file.ContentType,CRLF)
    write(buf,CRLF)
    write(stream,take!(buf))
end

# Write a file by reading it in 1MB chunks (unless we know its size and it's smaller than that)
function write_file(stream,file::IO,datasize,doclose)
    datasize == datasize == -1 : 2^20 : min(2^20,datasize)
    x = Array(UInt8,datasize)
    while !eof(file)
        nread = readbytes!(file,x)
        if nread == 2^20
            write(stream,x)
        else
            write(stream,sub(x,1:nread))
        end
    end
    doclose && close(file)
end

# Write a file by mmaping it
function write_file(stream,file::IOStream,datasize,doclose)
    @assert datasize != -1
    write(stream, Mmap.mmap(file, Vector{UInt8}, datasize, position(file)))
    doclose && close(file)
end

# Write data already in memory
function write_file(stream,file::Union{AbstractString,Array{UInt8}},datasize,doclose)
    @assert datasize != -1
    write(stream,file)
    doclose && close(file)
end

function write_file(stream,file::IOBuffer,datasize,doclose)
    @assert datasize != -1
if VERSION < v"0.5.0-dev+4817"
    write(stream,sub(file.data,(position(file)+1):(position(file)+nb_available(file))))
else
    write(stream,Base.view(file.data,(position(file)+1):(position(file)+nb_available(file))))
end
    doclose && close(file)
end

function partheadersize(file,datasize,boundary)
    totalsize = 0
    # Chucksize =
    #   +  "--" (2) + boundary (sizeof(boundary)) + "\r\n" (2)
    totalsize += (2 + sizeof(boundary) ) + 2
    #   + multipart_mime + optional names + "\r\n"(2)
    totalsize += sizeof(multipart_mime)
    if !isempty(file.name)
        # +1 for "\""
        totalsize += sizeof(name_file) + sizeof(file.name) + 1
    end
    if !isempty(file.filename)
        # +1 for "\""
        totalsize += sizeof(filename_file) + sizeof(file.filename) + 1
    end
    totalsize += 2
    if !isempty(file.ContentType)
        # +2 for "\r\n"
        totalsize += sizeof(ContentType_header) + sizeof(file.ContentType) + 2
    end
    # "\r\n" + The actual data + "\r\n" (2)
    totalsize += 2 + datasize + 2
    totalsize
end

choose_boundary() = hex(rand(UInt128))

function send_multipart(stream, settings, files)
    chunked, boundary, datasizes = settings.chunked, settings.boundary, settings.datasizes
    if chunked
        begin
            for i = 1:length(files)
                file = files[i]
                if datasizes[i] != -1
                    # Make this all one chunk
                    #write(stream,hex(datasizes[i]+partheadersize(file,0,boundary)),CRLF)
                    write_part_header(ChunkedStream(stream),file, boundary)
                    write_file(ChunkedStream(stream),file.file,datasizes[i],file.close)
                    # File CRLF
                    write(ChunkedStream(stream),CRLF)
                    # Chunk CRLF
                    #write(stream,CRLF)
                else
                    phs = partheadersize(file,0,boundary)-2
                    # Make the part header one chunk
                    write(stream,hex(phs),CRLF)
                    write_part_header(stream,file, boundary)
                    # Chunk CRLF
                    write(stream,CRLF)
                    # Write the rest as a chunk
                    write_file(ChunkedStream(stream),file.file,datasizes[i],file.close)
                    # This sucks, I'm making and extra chunk just for CRLF, but
                    # so be it for now
                    write(stream,"1\r\n\r\n\r\n")
                end
            end
            write(ChunkedStream(stream),"--$boundary--")
            write(stream,string(hex(0),CRLF,CRLF))
        end
    else
        begin
            for i = 1:length(files)
                file = files[i]
                write_part_header(stream,file, boundary)
                write_file(stream,file.file,datasizes[i],file.close)
                write(stream,CRLF)
            end
            write(stream, "--$boundary--", CRLF)
        end
    end

end

immutable MultipartSettings
    datasizes::Vector{Int}
    boundary::Compat.UTF8String
    chunked::Bool
end

function prepare_multipart_request!(request, files)
    local boundary
    headers = request.headers
    if !haskey(headers,"Content-Type")
        boundary = choose_boundary()
        headers["Content-Type"] = multipart_mime*boundary
    else

        if headers["Content-Type"][1:sizeof(multipart_mime)] != multipart_mime
            error("Cannot extract boundary from MIME type")
        end
        boundary = headers["Content-Type"][(sizeof(multipart_mime)+1):end]
    end

    chunked = false
    if haskey(headers,"Transfer-Encoding")
        if headers["Transfer-Encoding"] != "chunked"
            error("Unrecognized Transfer-Encoding")
        end
        chunked = true
    end

    datasizes = Array(Int,length(files))

    # Try to determine final size of the request. If this fails,
    # we fall back to chunked transfer
    totalsize = 0
    for i = 1:length(files)
        file = files[i]
        size = datasizes[i] = datasize(file.file)
        if size == -1
            if !chunked
                error("""Tried to pass in an IO object that is not of fixed size.\n
                         This is only support with the chunked Transfer-Encoding.\n
                         Please verify however that the server you are connecting to\n
                         supports chunked transfer encoding as support for this feature\n
                         is broken in a large number of servers.\n""")
            end
            # don't break because we'll still use the datasize later if
            # available to optimize chunked transfer
        end
        totalsize += partheadersize(file,size,boundary)
    end
    # "--" (2) + boundary (sizeof(boundary)) + "--" (2) + CRLF (2)
    totalsize += 2 + sizeof(boundary) + 2 + 2

    if chunked
        headers["Transfer-Encoding"] = "chunked"
    else
        headers["Content-Length"] = dec(totalsize)
    end

    MultipartSettings(datasizes, boundary, chunked)
end

function send_multipart(uri, headers, files, verb, timeout, tls_conf)
    req, datasizes, boundary, chunked = prepare_multipart_send(uri,headers,files,verb)
    stream = open_stream(uri,req,tls_conf)
    do_multipart_send(stream,files,datasizes, boundary, chunked)
    process_response(stream, timeout), req
end


text:
import PushInterface: LIBUSB
import .LIBUSB: libusb_init, libusb_set_debug, libusb_open_device_with_vid_pid, libusb_get_device, libusb_get_device_descriptor, libusb_kernel_driver_active, libusb_detach_kernel_driver, libusb_claim_interface, libusb_error_name, libusb_alloc_transfer, libusb_fill_bulk_transfer, libusb_submit_transfer, libusb_cancel_transfer, libusb_handle_events_completed, libusb_free_transfer, libusb_release_interface, libusb_close
import .LIBUSB: LibusbContext, LibusbDeviceDescriptor, LibusbTransfer # Ptr
import .LIBUSB: LIBUSB_LOG_LEVEL_DEBUG, LIBUSB_SUCCESS

const ABLETON_VENDOR_ID = 0x2982
const PUSH2_PRODUCT_ID  = 0x1967

ctx = Ptr{LibusbContext}(C_NULL)
init_result = libusb_init(Ptr{Ptr{LibusbContext}}(ctx))
libusb_set_debug(ctx, LIBUSB_LOG_LEVEL_DEBUG)
device_handle = libusb_open_device_with_vid_pid(ctx, ABLETON_VENDOR_ID, PUSH2_PRODUCT_ID)
if C_NULL == device_handle
    @error :libusb_open_device_with_vid_pid "could not found ableton push 2"
    exit(1)
end
device = libusb_get_device(device_handle)
descriptor = Ref{Ptr{LibusbDeviceDescriptor}}()
libusb_get_device_descriptor(device, descriptor)
# @info :descriptor_bDeviceClass descriptor[].bDeviceClass # 0 LIBUSB_CLASS_PER_INTERFACE
# @info :descriptor_idVendor     descriptor[].idVendor
# @info :descriptor_idProduct    descriptor[].idProduct

interface_number = Cint(0)

# 0 no kernel driver is active
# 1 if a kernel driver is active
kernel_driver_is_active = libusb_kernel_driver_active(device_handle, interface_number) == 1
if kernel_driver_is_active
    detach_result = libusb_detach_kernel_driver(device_handle, interface_number)
    @info :detach_result detach_result
end

claim_result = libusb_claim_interface(device_handle, interface_number)
if claim_result != LIBUSB_SUCCESS
    @info :claim_result libusb_error_name(claim_result)
end

iso_packets = Cint(0)
frame_header_transfer = libusb_alloc_transfer(iso_packets)
C_NULL == frame_header_transfer && @error "could not allocate frame header transfer handle"

frame_header = [
  0xff, 0xcc, 0xaa, 0x88,
  0x00, 0x00, 0x00, 0x00,
  0x00, 0x00, 0x00, 0x00,
  0x00, 0x00, 0x00, 0x00]

const PUSH2_BULK_EP_OUT = Cuchar(0x01)
const TRANSFER_TIMEOUT  = Cuint(1000) # milliseconds

function libusb_transfer_cb(transfer::Ptr{LibusbTransfer})
    t = unsafe_load(transfer)
    @info :libusb_transfer_cb t.status t.callback t.user_data
end
on_frame_header_transfer_finished = @cfunction libusb_transfer_cb Cvoid (Ptr{LibusbTransfer},)

libusb_fill_bulk_transfer(
    frame_header_transfer,
    device_handle,
    PUSH2_BULK_EP_OUT,
    pointer(frame_header),
    length(frame_header),
    on_frame_header_transfer_finished,
    C_NULL,
    TRANSFER_TIMEOUT)

transfer_result = libusb_submit_transfer(frame_header_transfer)
if transfer_result != LIBUSB_SUCCESS
    @info :transfer_result transfer_result
end

libusb_cancel_transfer(frame_header_transfer)
completed = Ptr{Cint}(C_NULL)
events_completed_result = libusb_handle_events_completed(ctx, completed)
if events_completed_result != LIBUSB_SUCCESS
    @info :events_completed_result events_completed_result completed
end
libusb_free_transfer(frame_header_transfer)
libusb_release_interface(device_handle, interface_number)
libusb_close(device_handle)


text:
module ExamplePlots

using Reexport
@reexport using Plots
import DataFrames, RDatasets

export
  test_examples

include("example_generation.jl")

end # module


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _planck_wave{T<:AbstractFloat}(wavelength::T, temperature::T)
    c1  = 3.741771790075259e-16 # = 2*pi*h*c*c
    c2  = 1.43877735382772e-2   # = h*c/k
    return c1/(wavelength^5*(expm1(c2/wavelength/temperature)))
end

"""
    planck_wave(wavelength, temperature) -> black_body_flux

### Purpose ###

Calculate the flux of a black body per unit wavelength.

### Explanation ###

Return the spectral radiance of a black body per unit wavelength using [Planck's
law](https://en.wikipedia.org/wiki/Planck%27s_law)

\$\$ B_\\lambda(\\lambda, T) =\\frac{2hc^2}{\\lambda^5}\\frac{1}{e^{\\frac{hc}{\\lambda k_\\mathrm{B}T}} - 1} \$\$

### Arguments ###

* `wavelength`: wavelength at which the flux is to be calculated, in meters.
* `temperature`: the equilibrium temperature of the black body, in Kelvin.

Both arguments can be either scalar or arrays of the same length.

### Output ###

The spectral radiance of the black body, in units of W/(sr·m³).

### Example ###

Plot the spectrum of a black body in \$[0, 3]\$ µm at \$5000\$ K.  Use
[PyPlot.jl](https://github.com/stevengj/PyPlot.jl) for plotting.

``` julia
using PyPlot
wavelength=linspace(0, 3e-6, 1000);
temperature=ones(wavelength)*5000;
flux=planck_wave(wavelength, temperature);
plot(wavelength, flux)
```

### Notes ###

`planck_freq` calculates the flux of a black body per unit frequency.

Code of this function is based on IDL Astronomy User's Library.
"""
planck_wave(w::Real, t::Real) = _planck_wave(promote(float(w), float(t))...)

@vectorize_2arg Real planck_wave


text:
function MakeP_bar(A,b,C,d,XLB,XUB,y;Cuts = [])
    m,n = size(A);
    nObj = size(C)[1];
    MyModel = Model();
    @variable(MyModel, XLB[i] <= X[i=1:n] <= XUB[i])
    @variable(MyModel,Z)
    @objective(MyModel,Min,Z)
    @constraintref Ucons[1:m];
    @constraintref LAMBDAcons[1:nObj];
    for i = 1:m
        Ucons[i] = @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    for j = 1:nObj
        LAMBDAcons[j] = @constraint(MyModel,sum(C[j,i]*X[i] for i = 1:n) + d[j] - Z <= y[j] )
    end
    for c = 1:length(Cuts)
        @constraint(MyModel,sum(Cuts[c].Coeff[i]*(sum(C[i,j]*X[j] for j = 1:n )+d[i])  for i = 1:nObj) >= Cuts[c].RHS)
    end
    return MyModel,Ucons,LAMBDAcons;

end


text:
should_precompile = true

# Don't edit the following! Instead change the script for `snoop_bot`.
ismultios = true
ismultiversion = true
# precompile_enclosure
@static if !should_precompile
    # nothing
elseif !ismultios && !ismultiversion
    @static if isfile(
        joinpath(@__DIR__, "../deps/SnoopCompile/precompile/precompile_AcuteML.jl"),
    )
        include("../deps/SnoopCompile/precompile/precompile_AcuteML.jl")
        _precompile_()
    end
else
    @static if Sys.islinux()
        @static if v"1.2.0-DEV" <= VERSION <= v"1.2.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.2/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.2/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.3/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.3/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.4/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.4/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.5.0-DEV" <= VERSION <= v"1.5.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.5/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.5/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.6.0-DEV" <= VERSION <= v"1.6.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.6/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.6/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.7.0-DEV" <= VERSION <= v"1.7.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.7/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.7/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.8.0-DEV" <= VERSION <= v"1.8.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/linux/1.8/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/linux/1.8/precompile_AcuteML.jl")
                _precompile_()
            end
        else
        end

    elseif Sys.iswindows()
        @static if v"1.2.0-DEV" <= VERSION <= v"1.2.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.2/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.2/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.3/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.3/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.4/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.4/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.5.0-DEV" <= VERSION <= v"1.5.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.5/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.5/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.6.0-DEV" <= VERSION <= v"1.6.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.6/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.6/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.7.0-DEV" <= VERSION <= v"1.7.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.7/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.7/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.8.0-DEV" <= VERSION <= v"1.8.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/windows/1.8/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/windows/1.8/precompile_AcuteML.jl")
                _precompile_()
            end
        else
        end

    elseif Sys.isapple()
        @static if v"1.2.0-DEV" <= VERSION <= v"1.2.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.2/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.2/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.3/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.3/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.4/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.4/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.5.0-DEV" <= VERSION <= v"1.5.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.5/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.5/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.6.0-DEV" <= VERSION <= v"1.6.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.6/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.6/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.7.0-DEV" <= VERSION <= v"1.7.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.7/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.7/precompile_AcuteML.jl")
                _precompile_()
            end
        elseif v"1.8.0-DEV" <= VERSION <= v"1.8.9"
            @static if isfile(
                joinpath(
                    @__DIR__,
                    "../deps/SnoopCompile/precompile/apple/1.8/precompile_AcuteML.jl",
                ),
            )
                include("../deps/SnoopCompile/precompile/apple/1.8/precompile_AcuteML.jl")
                _precompile_()
            end
        else
        end

    else
    end

end # precompile_enclosure

max_stars_repo_name: aminya/AML

max_forks_repo_name: aminya/AML


text:
# This file updates `docs/src/registry.md` (it does not build documentation in general)
module RegistryMD

using FileIO: FileIO, LOAD, SAVE, OSX, OS, DataFormat, @format_str, MimeWriter
using UUIDs
using Pkg

const ctx = Pkg.Types.Context()
try
    ctx.registries
catch
    error("this file must be run on Julia 1.6 or higher")
end

function pkg_url(uuid)
    for reg in ctx.registries
        entry = get(reg.pkgs, uuid, nothing)
        entry === nothing && continue
        info = Pkg.Registry.init_package_info!(entry)
        base, ext = splitext(info.repo)
        return ext == ".git" ? base : info.repo
    end
end
library2string(x::Pair{<:Any,UUID}) = "[$(x.first)]($(pkg_url(x.second)))"
library2string(x::Module) = string(nameof(x))

extension2string(x) = join(map(string, x), ", ")
extension2string(x::AbstractString) = x

os2string(x::Vector) = isempty(x) ? "all platforms " : join(map(os2string, x), ", ")
os2string(os::Type{O}) where {O <: OS} = "$(O.name.name)"

magic2string(x::Function) = "has detection function"
magic2string(x::Tuple) = isempty(x) ? "only extension" : string(x)
magic2string(x::AbstractString) = '"' * escape_string(x, ['_']) * '"'
magic2string(x) = string(x)
function loadsave2string(load_save_libraries)
    io = IOBuffer()
    loader_str, saver_str = " ", " "
    for (i, predicates) in enumerate(load_save_libraries)
        library = popfirst!(predicates)

        os, loadsave = FileIO.split_predicates(predicates)
        if isempty(loadsave)
            print(io, "loads and saves on all platforms with ", library2string(library), " ")
        elseif (LOAD in loadsave)
            print(io, "loads with ", library2string(library), " on ", os2string(os), " ")
        elseif (SAVE in loadsave)
            print(io, "saves with ", library2string(library), " on ", os2string(os), " ")
        end
        i < length(load_save_libraries) && print(io, '\n')
    end
    split(String(take!(io)), '\n')
end

fs = open(joinpath(pkgdir(FileIO), "docs", "src", "registry.md"), "w")

    println(fs, "# Registry table")
    println(fs)
    println(fs, "The following formats are registered with FileIO:")
    println(fs)

    function add_format(::Type{DataFormat{Sym}}, magic, extension, io_libs...) where Sym
        liblinks = loadsave2string(io_libs)
        for (i, lib) in enumerate(liblinks)
            if i == 1
                println(fs, "| $(Sym) | $(extension2string(extension)) | $lib | $(magic2string(magic)) |")
            else
                println(fs, "| | | $lib | |")
            end
        end
    end

    function add_format(fmt::Type{DataFormat{sym}}, magic::Union{Tuple,AbstractVector,String}, extension) where sym
        println(sym)
    end

    # for multiple magic bytes
    function add_format(sym::Symbol, magics::Vector{Vector{UInt8}}, extension)
        println(sym)
    end

    # For when "magic" is supplied as a function (see the HDF5 example in
    # registry.jl)
    function add_format(fmt::Type{DataFormat{sym}}, magic, extension) where sym
        println(sym)
    end

    println(fs, """
    | Format Name | extensions | IO library | detection or magic number |
    | ----------- | ---------- | ---------- | ---------- |""")
    include(joinpath(pkgdir(FileIO), "src", "registry.jl"))

close(fs)

end  # module RegistryMD


text:
include("setup.jl")

@testset "JET.jl" begin
    @testset "virtualprocess.jl" begin
        include("test_virtualprocess.jl")
    end

    @testset "abstractinterpretation.jl" begin
        include("test_abstractinterpretation.jl")
    end

    @testset "tfuncs.jl" begin
        include("test_tfuncs.jl")
    end

    @testset "jetcache.jl" begin
        include("test_jetcache.jl")
    end

    @testset "reports.jl" begin
        include("test_reports.jl")
    end

    # tests with Windows-paths is just an hell
    @static Sys.iswindows() || @testset "print.jl" begin
        include("test_print.jl")
    end

    @testset "misc" begin
        include("test_misc.jl")
    end

    @testset "performance" begin
        include("test_performance.jl")
    end

    @testset "interface.jl" begin
        include("test_interfaces.jl")
    end
end


text:
"""
Shared code for file io.  Currently, a convenience wrapper for FITSIO.

Author: Eric Ford
Created: August 2020
"""

""" `read_header( fits_file )`
Read header from FITS file and return Dict with contents.
Optional inputs:
- hdu: Specifies which HDU to read from the FITS file.  (Default: 1)
"""
function read_header(f::FITS; header_idx::Integer = 1)
    @assert 1<=header_idx<=length(f)
    #@assert read_key(f[header_idx],"SKY-OBJ")[1] == "Solar"
    hdr = FITSIO.read_header(f[header_idx])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
end

""" `read_header( filename )`
Read header from FITS file and return Dict with contents.
Optional inputs:
- hdu: Specifies which HDU to read from the FITS file.  (Default: 1)
"""
function read_header(fn::String; header_idx::Integer = 1)
    #println("# Reading: ",fn, " hdu= ",header_idx)
    f = FITS(fn)
    read_header(f, header_idx=header_idx)
end

""" `read_fits_header( filename )`
Read header from FITS file and return Dict with contents.
Optional inputs:
- hdu: Specifies which HDU to read from the FITS file.  (Default: 1)
"""
function read_fits_header(fn::String; header_idx::Integer = 1)
    read_header(fn,header_idx=header_idx)
end

""" `read_metradata_from_fits( filename, fields )`
Read metadata in FITS header for specified keys and return data as a Dict.  `fields` can be an array of symbols or strings.
Optional inputs:
- Passing both fields (an array of symbols) and fields_str (an array of strings) as named parameters allows for differences in string used in FITS file and Symbol used in resulting Dict.
- hdu: Specifies which HDU to read from the FITS file.  (Default: 1)
"""
function read_metradata_from_fits
end

function read_metadata_from_fits(file::Union{String,FITS}, fields::Array{Symbol,1} ; hdu::Union{Integer,String} = 1)
    fields_str=string.(fields)
    read_metadata_from_fits(file,hdu=hdu,fields=fields,fields_str=fields_str)
end

function read_metadata_from_fits(file::Union{String,FITS}, fields_str::AbstractArray{AS,1} ; hdu::Union{Integer,String} = 1)  where { AS<:AbstractString }
    fields = map(f->Symbol(f),fields_str)
    read_metadata_from_fits(file,hdu=hdu,fields=fields,fields_str=fields_str)
end

function read_metadata_from_fits(f::FITS; fields::Array{Symbol,1}, fields_str::AbstractArray{AS,1}, hdu::Union{Integer,String} = 1 )  where { AS<:AbstractString }
    @assert length(fields) == length(fields_str)
    @assert typeof(hdu) == String || 1 <= hdu <= length(f)
    hdr = FITSIO.read_header(f[hdu])
    #values = Vector{Any}(undef,length(fields))
    values_mask = falses(length(fields))
    for (i,field) in enumerate(fields_str)
        #println("  ",field,": ",typeof(hdr[field]))
        #@assert findfirst(isequal(field),keys(hdr)) != nothing
        idx = findfirst(isequal(field),keys(hdr))
        if idx != nothing
            #values[i] = hdr[fields_str[idx]]
            values_mask[i] = true
        end
    end
    values = map(s->hdr[s],fields_str[values_mask])

    df = Dict{Symbol,Any}(zip(fields[values_mask],values))
    return df
end

function read_metadata_from_fits(fn::String; fields::Array{Symbol,1}, fields_str::AbstractArray{AS,1}, hdu::Union{Integer,String} = 1 )  where { AS<:AbstractString }
    f = FITS(fn)
    read_metadata_from_fits(f; fields=fields, fields_str=fields_str, hdu=hdu )
end

function check_metadata_fields_expected_present(df::Dict,  fields::Array{Symbol,1} ) #, fields_str::AbstractArray{AS,1} )  where { AS<:AbstractString }
    # Check that header has all expected fields
    for field in fields
        if !haskey(df, field)
            println("# metadata lack entry:  ",field)
            return false
        end
    end
    return true
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using PyPlot

f(x) = -sin(x[1]^2/2 - x[2]^2/4 + 3) * cos(2x[1] + 1 - exp(x[2]))

function fplot(f)
    x = range(-2.5, stop=2.5, length=200)
    y = range( -1.5, stop=1.5, length=100)
    u = [f([x,y]) for y in y, x in x]

    gcf().set_size_inches(12,5)
    contour(x, y, u, 25, colors="k", linewidths=0.5)
    contourf(x, y, u, 25)
    axis("equal")
    colorbar()
    return
end

fplot(f)

function df(x)
    a1 = x[1]^2/2 - x[2]^2/4 + 3
    a2 = 2x[1] + 1 - exp(x[2])
    b1 = cos(a1)*cos(a2)
    b2 = sin(a1)*sin(a2)
    return -[x[1]*b1 - 2b2, -x[2]/2*b1 + exp(x[2])*b2]
end

function gradient_descent(f, df, x0, α=0.1; tol=1e-4, maxiter=500)
    x = x0
    xs = [x0]
    for i = 1:maxiter
        gradient = df(x)
        if sqrt(sum(gradient.^2)) < tol
            break
        end
        x -= α*gradient
        push!(xs, x)
    end
    return xs
end

function plot_path(xs)
    plot(first.(xs), last.(xs), "w.", markersize=6)
    plot(first.(xs), last.(xs), "r", linewidth=1)
end

x0s = [[-2,.5], [0,.5], [2.2,-0.5]];

fplot(f)
title("Gradient descent, fixed \$\\alpha\$")
for x0 in x0s
    xs = gradient_descent(f, df, x0, 0.3)
    plot_path(xs)
    println("Path length = $(length(xs)), ||gradient|| = $(sqrt(sum(df(xs[end]).^2)))")
end

function line_search(f, direction, x, αmin=1/2^20, αmax=2^20)
    α = αmin
    fold = f(x)
    while true
        if α ≥ αmax
            return α
        end
        fnew = f(x - α*direction)
        if fnew ≥ fold
            return α/2
        else
            fold = fnew
        end
        α *= 2
    end
end

function gradient_descent_line_search(f, df, x0; tol=1e-4, maxiter=500)
    x = x0
    xs = [x0]
    for i = 1:maxiter
        gradient = df(x)
        if sqrt(sum(gradient.^2)) < tol
            break
        end
        α = line_search(f, gradient, x)
        x -= α*gradient
        push!(xs, x)
    end
    return xs
end

fplot(f)
title("Gradient descent, line searches for \$\\alpha\$")
for x0 in x0s
    xs = gradient_descent_line_search(f, df, x0)
    plot_path(xs)
    println("Path length = $(length(xs)), ||gradient|| = $(sqrt(sum(df(xs[end]).^2)))")
end

function finite_difference_gradient(f, x, ϵ=1e-5)
    n = length(x)
    df = zeros(n)
    for k = 1:n
        x1 = copy(x)
        x1[k] += ϵ
        fP = f(x1)
        x1[k] -= 2ϵ
        fM = f(x1)
        df[k] = (fP - fM) / 2ϵ
    end
    return df
end

fplot(f)
title("Gradient descent, line searches for \$\\alpha\$, numerical gradients")
for x0 in x0s
    num_df(x) = finite_difference_gradient(f, x)
    xs = gradient_descent_line_search(f, num_df, x0)
    plot_path(xs)
    println("Path length = $(length(xs)), ||gradient|| = $(sqrt(sum(df(xs[end]).^2)))")
end

function finite_difference_hessian(f, x, ϵ=1e-5)
    n = length(x)
    ddf = zeros(n,n)
    f0 = f(x)
    for i = 1:n
        x1 = copy(x)
        x1[i] += ϵ
        fP = f(x1)
        x1[i] -= 2ϵ
        fM = f(x1)
        ddf[i,i] = (fP - 2*f0 + fM) / ϵ^2
    end
    for i = 1:n, j = i+1:n
        x1 = copy(x)
        x1[i] += ϵ
        x1[j] += ϵ
        fPP = f(x1)
        x1[i] -= 2ϵ
        fMP = f(x1)
        x1[j] -= 2ϵ
        fMM = f(x1)
        x1[i] += 2ϵ
        fPM = f(x1)
        ddf[i,j] = (fPP - fMP - fPM + fMM) / 4ϵ^2
        ddf[j,i] = ddf[i,j]
    end
    return ddf
end

fplot(f)
title("Newton's method, line searches for \$\\alpha\$, numerical gradients and Hessians")
x0s = [[-2,0.0], [0,.5], [2.,0.]];
for x0 in x0s
    search_dir(x) = finite_difference_hessian(f,x) \ finite_difference_gradient(f, x)
    xs = gradient_descent_line_search(f, search_dir, x0)
    plot_path(xs)
    println("Path length = $(length(xs)), ||gradient|| = $(sqrt(sum(df(xs[end]).^2)))")
end


text:
"""Test
Test to launch in complex-modeler/src_test folder for example D:\repo\complex-modeler\src_test>julia test_feasibility_GOC.jl
# Arguments
- instances path
- instance folder
- scenario folder
- epsilon
Test to check feasibility of solution point given by GOC of VOLTM, Genbounds, Smax and coupling constraints
Print the names of the constraints violated along with a message

"""
function read_arguments(ARGS)
    if length(ARGS)!=4
        error("# Arguments
        - instances path
        - instance folder
        - scenario folder
        - epsilon")
    else
        instance_path = joinpath(pwd(),ARGS[1],ARGS[2],ARGS[3])
        epsilon = float(ARGS[4])
        return instance_path, epsilon
    end

end

instance_path, epsilon = read_arguments(ARGS)

ROOT = pwd()
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

violations = test_feasibility_GOC(instance_path,epsilon)
for (scenario, dict) in violations
    for (ctrname, msg) in dict
        @printf("%30s %50s\n", ctrname, msg)
    end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
"""
Returns variable cost bids time-series data.

# Arguments
- `ts::IS.TimeSeriesData`:TimeSeriesData
- `component::Component`: Component
- `start_time::Union{Nothing, Dates.DateTime} = nothing`: Time when the time-series data starts
- `len::Union{Nothing, Int} = nothing`: Length of the time-series to be returned
"""
function get_variable_cost(
    ts::IS.TimeSeriesData,
    component::Component,
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    if start_time === nothing
        start_time = IS.get_initial_timestamp(ts)
    end
    data = IS.get_time_series_array(component, ts, start_time, len = len)
    time_stamps = TimeSeries.timestamp(data)
    return TimeSeries.TimeArray(time_stamps, map(VariableCost, TimeSeries.values(data)))
end

"""
Returns variable cost bids time-series data for  MarketBidCost.

# Arguments
- `device::StaticInjection`: Static injection device
- `cost::MarketBidCost`: Operations Cost
- `start_time::Union{Nothing, Dates.DateTime} = nothing`: Time when the time-series data starts
- `len::Union{Nothing, Int} = nothing`: Length of the time-series to be returned
"""
function get_variable_cost(
    device::StaticInjection,
    cost::OperationalCost;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    time_series_key = get_variable(cost)
    if isnothing(time_series_key)
        error(
            "Cost component has a `nothing` stored in field `variable`, Please use `set_variable_cost!` to add variable cost forecast.",
        )
    end
    raw_data = IS.get_time_series_by_key(
        time_series_key,
        device;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, device, start_time, len)
    return cost
end

"""
Returns variable cost time-series data for a ReserveDemandCurve.

# Arguments
- `service::ReserveDemandCurve`: ReserveDemandCurve
- `start_time::Union{Nothing, Dates.DateTime} = nothing`: Time when the time-series data starts
- `len::Union{Nothing, Int} = nothing`: Length of the time-series to be returned
"""
function get_variable_cost(
    service::ReserveDemandCurve;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    time_series_key = get_variable(service)
    if isnothing(time_series_key)
        error(
            "Cost component has a `nothing` stored in field `variable`, Please use `set_variable_cost!` to add variable cost forecast.",
        )
    end
    raw_data = IS.get_time_series_by_key(
        time_series_key,
        service;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, service, start_time, len)
    return cost
end

"""
Returns service bids time-series data for a device that has MarketBidCost.

# Arguments
- `sys::System`: PowerSystem System
- `cost::MarketBidCost`: Operations Cost
- `service::Service`: Service
- `start_time::Union{Nothing, Dates.DateTime} = nothing`: Time when the time-series data starts
- `len::Union{Nothing, Int} = nothing`: Length of the time-series to be returned
"""
function get_services_bid(
    device::StaticInjection,
    cost::MarketBidCost,
    service::Service;
    start_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Nothing, Int} = nothing,
)
    variable_ts_key = get_variable(cost)
    service_ts_key = IS.TimeSeriesKey(variable_ts_key.time_series_type, get_name(service))
    raw_data = IS.get_time_series_by_key(
        service_ts_key,
        device;
        start_time = start_time,
        len = len,
        count = 1,
    )
    cost = get_variable_cost(raw_data, device, start_time, len)
    return cost
end

"""
Adds energy market bids time-series to the MarketBidCost.

# Arguments
- `sys::System`: PowerSystem System
- `component::StaticInjection`: Static injection device
- `time_series_data::IS.TimeSeriesData`: TimeSeriesData
"""
function set_variable_cost!(
    sys::System,
    component::StaticInjection,
    time_series_data::IS.TimeSeriesData,
)
    add_time_series!(sys, component, time_series_data)
    key = IS.TimeSeriesKey(time_series_data)
    market_bid_cost = get_operation_cost(component)
    set_variable!(market_bid_cost, key)
    return
end

"""
Adds energy market bids time-series to the ReserveDemandCurve.

# Arguments
- `sys::System`: PowerSystem System
- `component::StaticInjection`: Static injection device
- `time_series_data::IS.TimeSeriesData`: TimeSeriesData
"""
function set_variable_cost!(
    sys::System,
    component::ReserveDemandCurve,
    time_series_data::IS.TimeSeriesData,
)
    add_time_series!(sys, component, time_series_data)
    key = IS.TimeSeriesKey(time_series_data)
    set_variable!(component, key)
    return
end

"""
Adds service bids time-series data to the MarketBidCost.

# Arguments
- `sys::System`: PowerSystem System
- `component::StaticInjection`: Static injection device
- `service::Service,`: Service for which the device is eligible to contribute
- `time_series_data::IS.TimeSeriesData`: TimeSeriesData
"""
function set_service_bid!(
    sys::System,
    component::StaticInjection,
    service::Service,
    time_series_data::IS.TimeSeriesData,
)
    if get_name(time_series_data) != get_name(service)
        error(
            "Name provided in the TimeSeries Data $(get_name(time_series_data)), doesn't match the Service $(get_name(service)).",
        )
    end
    verify_device_eligibility(sys, component, service)
    add_time_series!(sys, component, time_series_data)
    ancillary_services = get_ancillary_services(get_operation_cost(component))
    push!(ancillary_services, service)
    return
end

"""
Validates if a device is eligible to contribute to a service.

# Arguments
- `sys::System`: PowerSystem System
- `component::StaticInjection`: Static injection device
- `service::Service,`: Service for which the device is eligible to contribute
"""
function verify_device_eligibility(
    sys::System,
    component::StaticInjection,
    service::Service,
)
    contributing_devices = get_contributing_devices(sys, service)
    if !in(component, contributing_devices)
        error(
            "Device $(get_name(component)) isn't eligible to contribute to service $(get_name(service)).",
        )
    end
    return
end


text:
@testset "Evolutionary Objective" begin

    v = Evolutionary.default_values(ones(3))
    @test eltype(v) == Float64
    @test all(isnan.(v))

    v = Evolutionary.default_values(fill(true, 3))
    @test eltype(v) == Bool
    @test iszero(v)

    v = Evolutionary.default_values(BitVector(ones(Bool, 3)))
    @test eltype(v) == Bool
    @test iszero(v)

    v = Evolutionary.default_values(ones(Int,2, 3))
    @test eltype(v) == Int
    @test iszero(v)

    f1(x) = sum(x)
    objfun = EvolutionaryObjective(f1, [1; 2; 3; 4])
    @test typeof(objfun.F) == Int
    @test eltype(objfun.x_f) == Int
    f2(x) = sum(x)/length(x)
    objfun = EvolutionaryObjective(f2, [1; 2; 3; 4])
    @test typeof(objfun.F) == Float64
    @test eltype(objfun.x_f) == Int
    objfun = EvolutionaryObjective(f1, BitVector(ones(10)))
    @test typeof(objfun.F) == Int
    @test eltype(objfun.x_f) == Bool
    objfun = EvolutionaryObjective(f2, BitVector(ones(10)))
    @test typeof(objfun.F) == Float64
    @test eltype(objfun.x_f) == Bool
    f3(expr::Expr) = Evolutionary.Expression(expr).([1; 2; 3; 4]) |> sum
    objfun = EvolutionaryObjective(f3, Expr(:call, *, :x, :x))
    @test typeof(objfun.F) == Int
    @test typeof(objfun.x_f) == Expr

    f4(x) = x[1]+x[2]
    x = zeros(2)
    objfun = EvolutionaryObjective(f4, x)
    @test objfun.F isa Float64
    @test objfun.x_f isa Vector{Float64}
    @test !ismultiobjective(objfun)
    @test objfun.F === 0.0
    @test all(isnan.(objfun.x_f))
    v = value(objfun, ones(2))
    @test v == 2.0
    @test v isa Float64
    @test f_calls(objfun) == 1
    @test objfun.F === 0.0
    @test all(isnan.(objfun.x_f))
    @test_throws MethodError value(objfun, ones(Int, 2))
    v = value!(objfun, ones(2))
    @test v == 2.0
    @test v isa Float64
    @test f_calls(objfun) == 2
    @test objfun.F === 2.0
    @test all(isnan.(objfun.x_f))
    v = value!!(objfun, zeros(2))
    @test v === 0.0
    @test v isa Float64
    @test f_calls(objfun) == 3
    @test objfun.F === 0.0
    @test all(iszero.(objfun.x_f))

    function f5(F, x)
        F[1] = x[1]+1
        F[2] = x[2]+2
        F
    end
    objfun = EvolutionaryObjective(f5, x, copy(x))
    @test objfun.F   isa Vector{Float64}
    @test objfun.x_f isa Vector{Float64}
    @test ismultiobjective(objfun)
    @test_throws MethodError value(objfun, ones(Int, 2))
    @test iszero(objfun.F)
    @test all(isnan.(objfun.x_f))
    v = value(objfun, ones(2))
    @test eltype(v) == Float64
    @test v == [2.0, 3.0]
    @test f_calls(objfun) == 1
    @test iszero(objfun.F)
    @test all(isnan.(objfun.x_f))
    y = zeros(2)
    v = value(objfun, y, ones(2))
    @test v == [2.0, 3.0]
    @test v == y
    @test f_calls(objfun) == 2
    @test iszero(objfun.F)
    @test all(isnan.(objfun.x_f))
    y = zeros(2)
    v = value!(objfun, y, ones(2))
    @test v == [2.0, 3.0]
    @test v == y
    @test f_calls(objfun) == 3
    @test objfun.F == v
    @test all(isnan.(objfun.x_f))
    y = zeros(2)
    v = value!!(objfun, y, zeros(2))
    @test v == [1.0, 2.0]
    @test v == y
    @test f_calls(objfun) == 4
    @test objfun.F == v
    @test iszero(objfun.x_f)

    f6(x) = [x[1], x[2]+1]
    x = zeros(2)
    objfun = EvolutionaryObjective(f6, x)
    v = value(objfun, ones(2))
    @test v == [1.0, 2.0]
    @test iszero(objfun.F)
    @test all(isnan.(objfun.x_f))
    v = value!(objfun, ones(2))
    @test v == [1.0, 2.0]
    @test objfun.F == v
    @test all(isnan.(objfun.x_f))
    v = value!!(objfun, zeros(2))
    @test objfun.F == v
    @test iszero(objfun.x_f)

    X = [ones(2) for i in 1:10]
    F = zeros(10)
    objfun = EvolutionaryObjective(f4, x)
    V = value!(objfun, F, X);
    @test F == 2ones(10)

    objfun = EvolutionaryObjective(f4, x; eval=:thread)
    V = value!(objfun, F, X);
    @test F == 2ones(10)

    objfun = EvolutionaryObjective(f5, x, zeros(2))
    F = zeros(2,10)
    V = value!(objfun, F, X)
    @test F == [2ones(10) 3ones(10)]'

    objfun = EvolutionaryObjective(f5, x, zeros(2); eval=:thread)
    F = zeros(2,10)
    V = value!(objfun, F, X)
    @test F == [2ones(10) 3ones(10)]'

end


text:
using Matter

import Matter.Particle, Matter.Vector

position = Vector(0., 0.)
velocity = Vector(0., 0.)
gravity = Vector(0., -10.)
damping = 0.995
mass = 2.

p = Particle(position, velocity, gravity, damping, 1. / mass)

for step in 1:10
    dt = 0.033
    Matter.integrate!(p, dt)
    print(p.position)
end


text:
struct Bisection <: AbstractBracketingAlgorithm
  exact_left::Bool
  exact_right::Bool
end

function Bisection(;exact_left=false, exact_right=false)
  Bisection(exact_left, exact_right)
end

struct BisectionCache{uType}
  state::Int
  left::uType
  right::uType
end

function alg_cache(alg::Bisection, left, right, p, ::Val{true})
  BisectionCache(0, left, right)
end

function alg_cache(alg::Bisection, left, right, p, ::Val{false})
  BisectionCache(0, left, right)
end

function perform_step(solver::BracketingImmutableSolver, alg::Bisection, cache)
  @unpack f, p, left, right, fl, fr, cache = solver

  if cache.state == 0
    fzero = zero(fl)
    fl * fr > fzero && error("Bracket became non-containing in between iterations. This could mean that "
    * "input function crosses the x axis multiple times. Bisection is not the right method to solve this.")

    mid = (left + right) / 2

    if left == mid || right == mid
      @set! solver.force_stop = true
      @set! solver.retcode = FLOATING_POINT_LIMIT
      return solver
    end

    fm = f(mid, p)

    if iszero(fm)
      if alg.exact_left
        @set! cache.state = 1
        @set! cache.right = mid
        @set! cache.left = mid
        @set! solver.cache = cache
      elseif alg.exact_right
        @set! solver.left = prevfloat_tdir(mid, left, right)
        solver = sync_residuals!(solver)
        @set! cache.state = 2
        @set! cache.left = mid
        @set! solver.cache = cache
      else
        @set! solver.left = prevfloat_tdir(mid, left, right)
        @set! solver.right = nextfloat_tdir(mid, left, right)
        solver = sync_residuals!(solver)
        @set! solver.force_stop = true
        return solver
      end
    else
      if sign(fm) == sign(fl)
        @set! solver.left = mid
        @set! solver.fl = fm
      else
        @set! solver.right = mid
        @set! solver.fr = fm
      end
    end
  elseif cache.state == 1
    mid = (left + cache.right) / 2

    if cache.right == mid || left == mid
      if alg.exact_right
        @set! cache.state = 2
        @set! solver.cache = cache
        return solver
      else
        @set! solver.right = nextfloat_tdir(mid, left, right)
        solver = sync_residuals!(solver)
        @set! solver.force_stop = true
        return solver
      end
    end

    fm = f(mid, p)

    if iszero(fm)
      @set! cache.right = mid
      @set! solver.cache = cache
    else
      @set! solver.left = mid
      @set! solver.fl = fm
    end
  else
    mid = (cache.left + right) / 2

    if right == mid || cache.left == mid
      @set! solver.force_stop = true
      return solver
    end

    fm = f(mid, p)

    if iszero(fm)
      @set! cache.left = mid
      @set! solver.cache = cache
    else
      @set! solver.right = mid
      @set! solver.fr = fm
    end
  end
  solver
end


text:
module Phylogenetics
	## Dependencies
	using Base.Intrinsics
	importall Base

	## Exported Methods and Types
	export Phylogeny,
		Phylo,
		Clado,
		ReducedTopology,
		readtree,
		writetree,
		getroot,
		getkids,
		@tr_str

	## Load Package Files
	include(Pkg.dir("Phylogenetics", "src", "typedefs.jl"))
	include(Pkg.dir("Phylogenetics", "src", "treeio.jl"))

end


text:
export EncodingOp, lrEncodingOp, fourierEncodingOp2d, fourierEncodingOp3d,
       encodingOps2d_simple, encodingOps3d_simple, encodingOps2d_parallel,
       encodingOps3d_parallel, encodingOp2d_multiEcho, encodingOp3d_multiEcho,
       encodingOp2d_multiEcho_parallel, encodingOp3d_multiEcho_parallel

"""
    encodingOps2d_simple(acqData::AcquisitionData, shape::NTuple{2,Int64}
                              ; kargs...)

generates an Array of LinearOperators which describe 2d signal encoding of the individual
contrasts in an MRI acquisition (for a given slice).

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{2,Int64}`              - size of image to be encoded/reconstructed
"""
function encodingOps2d_simple(acqData::AcquisitionData, shape::NTuple{2,Int64}; kargs...)
  numContr = numContrasts(acqData)
  tr = [trajectory(acqData,i) for i=1:numContr]
  idx = acqData.subsampleIndices
  return [fourierEncodingOp2d(shape, tr[i], "fast",
                 subsampleIdx=idx[i]; kargs...) for i=1:numContr]
end

"""
    encodingOps3d_simple(acqData::AcquisitionData, shape::NTuple{2,Int64}; kargs...)

generates an Array of LinearOperators which describe 3d signal encoding of the individual
contrasts&coils in an MRI acquisition (for a given slice).
Arguments are the same as in the 2d case, with the exception that shape is of type `NTuple{3,Int64}`
and the considered slice is not specified.
"""
function encodingOps3d_simple(acqData::AcquisitionData, shape::NTuple{3,Int64}; kargs...)
  numContr = numContrasts(acqData)
  tr = [trajectory(acqData,i) for i=1:numContr]
  idx = acqData.subsampleIndices
  return [fourierEncodingOp3d(shape, tr[i], "fast", subsampleIdx=idx[i]; kargs...) for i=1:numContr]
end

"""
    encodingOps2d_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                              , senseMaps::Array{ComplexF64}
                              ; kargs...)

generates an Array of LinearOperators which describe 2d signal encoding of the individual
contrasts in an MRI acquisition (for a given slice). The different coils are taken into account
in terms of their sensitivities

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{2,Int64}`              - size of image to be encoded/reconstructed
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
"""
function encodingOps2d_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                                , senseMaps::Array{ComplexF64}
                                ; slice=1, kargs...)

  numContr, numChan = numContrasts(acqData), numChannels(acqData)
  # fourier operators
  ft = encodingOps2d_simple(acqData, shape; slice=slice, kargs...)
  S = SensitivityOp(reshape(senseMaps[:,:,slice,:],:,numChan),1)
  Op = [ diagOp(ft[i], numChan) ∘ S for i=1:numContr]

  return Op
end

"""
    encodingOps3d_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                              , senseMaps::Array{ComplexF64}
                              ; kargs...)

generates an Array of LinearOperators which describe 3d signal encoding of the individual
contrasts in an MRI acquisition (for a given slice). The different coils are taken into account
in terms of their sensitivities
Arguments are the same as in the 2d case, with the exception that shape is of type `NTuple{3,Int64}`
and the considered slice is not specified.
"""
function encodingOps3d_parallel(acqData::AcquisitionData, shape::NTuple{3,Int64}
                                , senseMaps::Array{ComplexF64}
                                ; kargs...)

  numContr, numChan = numContrasts(acqData), numChannels(acqData)
  # fourier operators
  ft = encodingOps3d_simple(acqData, shape; kargs...)
  S = SensitivityOp(reshape(senseMaps,:,numChan),1)
  return [ diagOp(ft[i], numChan) ∘ S for i=1:numContr]
end

"""
    encodingOp2d_multiEcho(acqData::AcquisitionData, shape::NTuple{2,Int64}; kargs...)

generates a LinearOperator which describe combined 2d signal encoding of all
the contrasts in an MRI acquisition (for a given slice).

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{2,Int64}`              - size of image to be encoded/reconstructed
"""
function encodingOp2d_multiEcho(acqData::AcquisitionData, shape::NTuple{2,Int64}
                                ; kargs...)
  # fourier operators
  ft = encodingOps2d_simple(acqData, shape; kargs...)
  return diagOp(ft...)
end

"""
    encodingOp3d_multiEcho(acqData::AcquisitionData, shape::NTuple{2,Int64}; kargs...)

generates a LinearOperator which describe combined 3d signal encoding of all
the contrasts in an MRI acquisition (for a given slice).
Arguments are the same as in the 2d case, with the exception that shape is of type `NTuple{3,Int64}`
and the considered slice is not specified.
"""
function encodingOp3d_multiEcho(acqData::AcquisitionData, shape::NTuple{3,Int64}; kargs...)
  # fourier operators
  ft = encodingOps3d_simple(acqData, shape; kargs...)
  return diagOp(ft...)
end

"""
    encodingOp2d_multiEcho_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                                          , senseMaps::Array{ComplexF64}
                                          ; kargs...)

generates a LinearOperator which describe combined 3d signal encoding of all
the contrasts in an MRI acquisition (for a given slice). The different coils are taken into account
in terms of their sensitivities

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{2,Int64}`              - size of image to be encoded/reconstructed
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
"""
function encodingOp2d_multiEcho_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                                          , senseMaps::Array{ComplexF64}
                                          ; slice::Int64=1, kargs...)

  numChan = numChannels(acqData)
  # fourier operators
  ft = encodingOps2d_simple(acqData, shape; kargs...)
  S = SensitivityOp(reshape(senseMaps[:,:,slice,:],:,numChan),numContrasts(acqData))
  return diagOp(ft, numChan) ∘ S
end

"""
    encodingOp3d_multiEcho_parallel(acqData::AcquisitionData, shape::NTuple{2,Int64}
                                          , senseMaps::Array{ComplexF64}
                                          ; kargs...)

generates a LinearOperator which describe combined 3d signal encoding of all
the contrasts in an MRI acquisition (for a given slice). The different coils are taken into account
in terms of their sensitivities
Arguments are the same as in the 2d case, with the exception that shape is of type `NTuple{3,Int64}`
and the considered slice is not specified.
"""
function encodingOp3d_multiEcho_parallel(acqData::AcquisitionData, shape::NTuple{3,Int64}
                                          , senseMaps::Array{ComplexF64}
                                          ; kargs...)

  numChan = numChannels(acqData)
  # fourier operators
  ft = encodingOps3d_simple(acqData, shape; kargs...)
  S = SensitivityOp(reshape(senseMaps,:,numChan),1)
  return diagOp(ft, numChan) ∘ S
end

###################################
# Encoding with low rank projection
###################################
function lrEncodingOp(acqData::AcquisitionData, shape, params::Dict; numContr::Int64=1, parallel::Bool=false)

  numChan = numChannels(acqData)
  # low rank operator
  N=prod(shape)
  subspace = get( params, :phi, Matrix{Float64}(I,N,N) )
  K = size(subspace,2)
  Φ = MapSliceOp(subspace[:,:,1],2,(N, K, numChan), (N, numContr, numChan))

  # Fourier Operator
  tr = trajectory(acqData,1)
  ft = fourierEncodingOp2d(shape, tr, "fast"; params...)

  # coil sensitivities in case of SENSE-reconstruction
  if parallel
    S = SensitivityOp(params[:senseMaps][:,:],K)
    E = diagOp( [ft for i=1:numChan*K]... )*S
  else
    E = diagOp( [ft for i=1:K]... )
  end

  #  TODO: sampling operator in case of undersampled cartesian acquisitions
  # vectorize sampling pattern
  if get(params,:sampling, nothing) == "binary"
    M = SamplingOp(Array{Bool}(acqData.subsampleIndices))
  else
    # sampling idx for all contrasts but one coil
    subIdx = hcat(acqData.subsampleIndices...)
    M = SamplingOp( hcat([subIdx for c=1:numChan]...), (N, numContr, numChan))
  end

  return M ∘ (Φ ∘ E)

end

"""
return 2d Fourier encoding operator (either Explicit or NFFT)
  * `opname` : "explicit" or "fast"
  * `slice` : slice to which the operator will be applied
  echoImage : calculate signal evolution relative to the echo time
"""
function fourierEncodingOp2d(shape::NTuple{2,Int64}, tr::Trajectory, opName::String;
          subsampleIdx::Vector{Int64}=Int64[], slice::Int64=1, correctionMap=[],
          echoImage::Bool=true,  kargs...)
  # Fourier transformations
  if opName=="explicit"
    @debug "ExplicitOp"
    ftOp = ExplicitOp(shape,tr,correctionMap[:,:,slice], echoImage=echoImage)
  elseif opName=="fast"
    @debug "NFFT-based Op"
    if !isempty(correctionMap) && correctionMap!=zeros(ComplexF64,size(correctionMap))
      ftOp = FieldmapNFFTOp(shape, tr, correctionMap[:,:,slice], echoImage=echoImage; kargs...)
    elseif isCartesian(tr)
      @debug "FFTOp"
      ftOp = FFTOp(ComplexF64, shape; unitary=false)
    else
      ftOp = NFFTOp(shape, tr; kargs...)
    end
  else
    @error "opName $(opName) is not known"
  end

  # subsampling
  if !isempty(subsampleIdx) && (subsampleIdx != collect(1:size(tr,2))) && isCartesian(tr)
    S = SamplingOp(subsampleIdx,(tr.numSamplingPerProfile,tr.numProfiles))
    return S ∘ ftOp
  else
    return ftOp
  end
end

"""
return 3d Fourier encoding operator (either Explicit, FFT or NFFT)
  opname : "explicit", "fft" or "fast"
  echoImage : calculate signal evolution relative to the echo time
"""
function fourierEncodingOp3d(shape::NTuple{3,Int64}, tr::Trajectory, opName::String;
          subsampleIdx::Vector{Int64}=Int64[], correctionMap=ComplexF64[], echoImage::Bool=true,
          kargs...)

  if opName=="explicit"
    ftOp = ExplicitOp(shape,tr,correctionMap,echoImage=echoImage)
  elseif opName=="fast"
    if !isempty(correctionMap) && correctionMap!=zeros(ComplexF64,size(correctionMap))
      ftOp = FieldmapNFFTOp(shape, tr, correctionMap, echoImage=echoImage; kargs...)
    elseif isCartesian(tr)
      #ftOp = sqrt(prod(shape))*FFTOp(ComplexF64, shape)
      ftOp = FFTOp(ComplexF64, shape; unitary=false)
    else
      ftOp = NFFTOp(shape, tr; kargs...)
    end
  else
    error("opName $(opName) is not known")
  end

  # subsampling
  if !isempty(subsampleIdx) && (subsampleIdx != collect(1:size(tr,2))) && isCartesian(tr)
    S = SamplingOp(subsampleIdx,shape)
    return S ∘ ftOp
  else
    return ftOp
  end
end


text:
module PlainMT

using Metatheory
using Metatheory.Library
using Metatheory.EGraphs

@metatheory_init()

static = @theory begin
    a + a => 2 * a
    a * (b * c) => (a * b) * c
end

dump(static[2])

dynamic = @theory begin
    a::Number * b::Number |> a * b
end

t = static ∪ dynamic

ex = Expr(:call, :+, 2, 2, Expr(:call, :*, 2, QuoteNode(:q)))
G = EGraph(ex)
saturate!(G, t)
ex = Metatheory.EGraphs.extract!(G, astsize)
display(ex)

end # module


text:
# include("EBNF.jl")

immutable Node
    name::AbstractString
    value::AbstractString
    first::Int
    last::Int
    children::Array #::Array{Node}
    ruleType::Type
    sym::Any

    Node(node::Node) =
        new(node.name, node.value, node.first, node.last, node.children, node.ruleType, node.sym)

    Node(name::AbstractString, value::AbstractString,
         first::Int, last::Int, children::Array, ruleType::Type) =
        new(name, value, first, last, children, ruleType,
            length(name)==0 ? nothing : symbol(name))
end

Node(name::AbstractString, value::AbstractString, first::Int, last::Int, typ) =
    Node(name, value, first, last, [], typ)

show{T}(io::IO, val::T, indent) = println(io, "$val ($(typeof(val)))")

function show(io::IO, node::Node, indent)
  println(io, "node($(node.name)) {$(displayValue(node.value, node.ruleType))$(node.ruleType)}")
  if isa(node.children, Array)
    for (i, child) in enumerate(node.children)
      print(io, "  "^indent)
      print(io, "$i: ")
      show(io, child, indent+1)
    end
  else
    print(io, "  "^(indent+1))
    show(io, node.children, indent+1)
  end
end

show(io::IO, node::Node) = show(io, node, 0)


text:
struct Language end
struct Algorithm end
struct Implementation end
struct Problem end
struct Experiment end
struct Run end

function New!(db, ::Type{Language}; name = missing, version = missing, commit = "")
    @check_not_missing "Language" (name, "name") (version, "version")
    @check_type String "Language" (name, "name") (version, "version") (commit, "commit")
    id = newid(db.language_versions_table)
    add_row!(db.language_versions_table, @NT(ID = id, Name = name, Version = version, Commit = commit))

    return
end

function New!(db, ::Type{Algorithm}; name = missing, description = "")
    @check_not_missing "Algorithm" (name, "name")
    @check_type String "Algorithm" (name, "name")
    id = newid(db.algorithms_table)
    add_row!(db.algorithms_table, @NT(ID = id, Name = name, Description = description))
    return
end

function New!(db, ::Type{Implementation}; algorithm = missing, language = missing, repository = missing, commit = missing, script = missing, func = missing, parameters = (), parameter_descriptions = (), results = missing)
    @check_not_missing "Implementation's" (algorithm, "algorithm") (language, "programming language") (repository, "repository") (script, "script") (func, "function") (results, "results")
    @check_type String "Implementation's" (repository, "repository") (script, "script")
    @check_type U(String) "Implementation's" (commit, "commit")
    @check_type Symbol "Implementation's" (func, "function")
    @check_type Tuple "Implementation's" (results, "results") (parameters, "parameters")
    length(results) == 0 || eltype(results) <: String || (warn("Results names should be of type String."); return )
    length(parameters) > 0 || eltype(parameters) <: Symbol || (warn("Parameters should be of type Symbol."); return )
    @check_valid_id db.language_versions_table ID language "Language"
    @check_valid_id db.algorithms_table ID algorithm "Algorithm"

    if commit isa Missing
        commit = HEAD(repository)
    end
    impl_id = newid(db.implementations_table)
    add_row!(db.implementations_table, @NT(ID = impl_id, Algorithm_ID = algorithm, Language_ID = language, Repository = repository, Commit = commit, Script = script, Function = func))

    for (i,p) in enumerate(parameters)
        param_id = newid(db.implementation_parameters_table)
        add_row!(db.implementation_parameters_table, @NT(ID = param_id, Implementation_ID = impl_id, Parameter_symbol = p, Parameter_description = parameter_descriptions == () ? "" : parameter_descriptions[i]))
    end

    for r in results
        result_id = newid(db.implementation_results_table)
        add_row!(db.implementation_results_table, @NT(ID = result_id, Implementation_ID = impl_id, Result_name = r))
    end

    return
end

function New!(db, ::Type{Problem}; description = "", kwargs...)
    problemid = newid(db.problems_table)
    add_row!(db.problems_table, @NT(ID = problemid, Description = description))
    for (k, v) in kwargs
        featureid = newid(db.problem_features_table)
        add_row!(db.problem_features_table, @NT(ID = featureid, Problem_ID = problemid, Feature_symbol = k, Feature_value = v))
    end
    return
end

function New!(db, ::Type{Experiment}; directory = pwd(), execute = true)
    id = newid(db.experiments_table)
    if !ispath(directory)
        mkdir(directory)
    end
    add_row!(db.experiments_table, @NT(ID = id, Directory = directory, Execute = execute))
    return
end

function New!(db, ::Type{Run}; experiment = missing, implementation = missing, parameters = missing, problems = missing)
    @check_not_missing "Run's" (experiment, "experiment") (implementation, "implementation") (parameters, "parameters")
    @check_type Union{NamedTuples.NamedTuple, Vector{<: NamedTuples.NamedTuple}} "Run's" (parameters, "parameters")
    @check_valid_id db.experiments_table ID experiment "Experiment"
    @check_valid_id db.implementations_table ID implementation "Implementation"

    if parameters isa Array
        param = parameters[1]
    else
        param = parameters
    end
    param_ids = Dict{Symbol, Int}()
    t = db.implementation_parameters_table
    for k in keys(param)
        x = @from i in t begin
            @where i.Implementation_ID == implementation && i.Parameter_symbol == k
            @select i.ID
            @collect
        end
        param_ids[k] = x[1]
    end
    if problems isa Missing
        all_problem_ids = columns(db.problems_table, :ID)
    else
        all_problem_ids = problems
    end
    for problem in all_problem_ids
        if typeof(parameters) <: AbstractVector
            for param in parameters
                runid = newid(db.runs_table)
                add_row!(db.runs_table, @NT(ID = runid, Experiment_ID = experiment, Problem_ID = problem, Implementation_ID = implementation, Executing = false, Executed = false))
                for (k,v) in zip(keys(param), values(param))
                    valueid = newid(db.run_implementation_parameter_table)
                    add_row!(db.run_implementation_parameter_table, @NT(ID = valueid, Implementation_parameter_ID = param_ids[k], Run_ID = runid, Value = v))
                end
            end
        else
            runid = newid(db.runs_table)
            add_row!(db.runs_table, @NT(ID = runid, Experiment_ID = experiment, Problem_ID = problem, Implementation_ID = implementation, Executing = false, Executed = false))
            for (k,v) in zip(keys(parameters), values(parameters))
                valueid = newid(db.run_implementation_parameter_table)
                add_row!(db.run_implementation_parameter_table, @NT(ID = valueid, Implementation_parameter_ID = param_ids[k], Run_ID = runid, Value = v))
            end
        end
    end

    return
end


text:
###############################################################################
##### Tree data structures for SoftBART
###############################################################################

abstract type Node end

mutable struct Branch <: Node
  var::Int
  cut::Float64
  left::Node
  right::Node
end

mutable struct Leaf <: Node
  μ::Float64
end

mutable struct Tree
  root::Node
  λ::Float64
end

###############################################################################
##### Tree utility functions
###############################################################################

function leafnodes(node::Node)
  isa(node, Leaf) ? [node] :
    reduce(vcat, [leafnodes(node.left), leafnodes(node.right)])
end

function Base.parent(node::Node, tree::Tree)
  parent(node, tree.root)
end

function Base.parent(node::Node, child_node::Branch)
  if (child_node.left == node) || (child_node.right == node)
    return child_node
  else
    isa(parent(node, child_node.left), Nothing) ?
      parent(node, child_node.right) : parent(node, child_node.left)
  end
end

function Base.parent(node::Node, child_node::Leaf)
  nothing
end

function onlyparents(tree::Tree)
  branches = Branch[]
  if isa(tree.root, Leaf)
    return [tree.root]
  else
    onlyparents(tree.root, branches)
  end
end

function onlyparents(branch::Branch, branches::Vector{Branch})
  if isa(branch.left, Leaf) & isa(branch.right, Leaf)
    push!(branches, branch)
  else
    onlyparents(branch.left, branches)
    onlyparents(branch.right, branches)
  end
  return branches
end

onlyparents(leaf::Leaf, branches::Vector{Branch}) = nothing

function depth(node::Node, tree::Tree)
  tree.root == node ? 0 : 1 + depth(parent(node, tree), tree)
end

function depth(tree)
  maximum([depth(leaf, tree) for leaf in leafnodes(tree.root)])
end

function isleft(node::Node, tree::Tree)
  parentnode = parent(node, tree)
  parentnode.left == node ? true : false
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "normaliz"
version = v"3.8.4"

# Collection of sources required to build normaliz
sources = [
    "https://github.com/Normaliz/Normaliz/releases/download/v$version/normaliz-$version.tar.gz" =>
    "80d21ebaf1a2d472ccdc1e1b2e42b4d71f45f3b8df4d7195ff83edf38f8945c8",

]

# Bash recipe for building across all platforms
script = raw"""
cd normaliz-3.8.4
# avoid libtool problems ....
rm /workspace/destdir/lib/libgmpxx.la
# add missing header
sed -i -e 's#^nobase_include_HEADERS = #nobase_include_HEADERS = libnormaliz/output.h #' source/Makefile.am source/Makefile.in
./configure --prefix=$prefix --host=$target --with-gmp=$prefix --disable-scip
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# windows build would require MPIR instead of GMP for 'long long'
platforms = filter(x->!isa(x,Windows),supported_platforms())

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libnormaliz", :libnormaliz)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl",
    "https://github.com/benlorenz/boostBuilder/releases/download/v1.70.0/build_boost.v1.70.0.jl"
]

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

max_stars_repo_name: benlorenz/normalizBuilder

max_issues_repo_name: benlorenz/normalizBuilder

max_forks_repo_name: benlorenz/normalizBuilder


text:
module A2R2
export GameEnv, GameSpec
include("game.jl")
end


text:
FILES = [
    Dict{String, Any}("sha" => "d3053f59200f5009b74bfed9a3f17d2d06783209", "path" => "secretsmanager-2017-10-17.examples.json"),
    Dict{String, Any}("sha" => "d3053f59200f5009b74bfed9a3f17d2d06783209", "path" => "secretsmanager-2021-01-01.examples.json"),
    Dict{String, Any}("sha" => "40526c64ec05df4f77499c2e0ec544e8818f3015", "path" => "secretsmanager-2017-10-17.min.json"),
    Dict{String, Any}("sha" => "0141ce42b03c3617179452e94af490043cecfdd4", "path" => "secretsmanager-2017-10-17.normal.json"),
    Dict{String, Any}("sha" => "23589621f3cb14d6911f04859f90c1af61bd8573", "path" => "secretsmanager-2017-10-17.paginators.json"),
    Dict{String, Any}("sha" => "0ea7e3b0bbe917eb027880396ac01509becd1fa0", "path" => "securityhub-2018-10-26.examples.json"),
    Dict{String, Any}("sha" => "cb76fc3ad2cce30155cfdf6f2c0c313922391d79", "path" => "securityhub-2018-10-26.min.json"),
    Dict{String, Any}("sha" => "8d144a13919a8c8ca017e073b2eea707d48be5e3", "path" => "securityhub-2018-10-26.normal.json"),
    Dict{String, Any}("sha" => "7726e5357a8eedc8b9f6871702eabd7a8c5d512a", "path" => "securityhub-2018-10-26.paginators.json"),
]

_get_service_and_version_patch = @patch function MockAWS._get_service_and_version(
    filename::String,
    files::Vector,
    auth::GitHub.OAuth2
)
    return split(FILES[1]["path"], '-')
end


text:
using GeneralizedSampling
using Base.Test

println("Testing collection of Freq2Wave type...")

#=
Test collecting of Freq2Wave:
- The Fourier transform of the non-boundary part should be the same in a Freq2NoBoundaryWave and Freq2BoundaryWave
- TODO: Test the boundary parts
=#

# ------------------------------------------------------------
# 1D

begin
	J = 3
	M = 2^(J+1)

	samples = grid(M, 0.5)

	# No boundary
	T = Freq2Wave(samples, "haar", J)
	A = collect(T)

	# With boundary (using the same internal Fourier transform as above)
	TMP = Freq2Wave(samples, "db2", J)
	Tb = GeneralizedSampling.Freq2BoundaryWave1D(T.internal, TMP.weights, TMP.J, TMP.wavename, TMP.NFFT, TMP.left, TMP.right)
	Ab = collect(Tb)

	vm = van_moment(Tb)
	# The internal parts are the same
	N = size(T,2)
	@test_approx_eq A[:,vm+1:N-vm] Ab[:,vm+1:N-vm]
end

# ------------------------------------------------------------
# 2D

begin
	J = 3
	M = 2^(J+1)

	samples = grid((M,M), 0.5)

	# No boundary
	T = Freq2Wave(samples, "haar", J)
	A = collect(T)

	# With boundary (using the same internal Fourier transform as above)
	TMP = Freq2Wave(samples, "db2", J)
	Tb = GeneralizedSampling.Freq2BoundaryWave2D(T.internal, TMP.weights, TMP.J, TMP.wavename, TMP.NFFT, TMP.left, TMP.right)
	Ab = collect(Tb)

	# Select the indices of the internal scaling functions
	vm = van_moment(Tb)
	J = zeros(Bool, wsize(T))
	S = split(J,vm)
	for idx in eachindex(S.II)
		S.II[idx] = true
	end

	#@show norm( A[:,vec(J)] - Ab[:,vec(J)] )
	@test_approx_eq A[:,vec(J)] Ab[:,vec(J)]
end


text:
DIR = joinpath("tmp","BSplineManifold")
rm(DIR, force=true, recursive=true)
mkpath(DIR)

@testset "BSplineManifold" begin
    Random.seed!(42)

    @testset "1d2d" begin
        p = 3 # degree of polynomial
        k = KnotVector(1:12) # knot vector
        P = BSplineSpace{p}(k) # B-spline space
        rand_a = [randn(2) for i in 1:dim(P)] / 2
        a = [SVector(2i-6, 0) for i in 1:dim(P)] + rand_a # random generated control points
        M = BSplineManifold(a, (P,)) # Define B-spline manifold

        @testset "luxor-svg" begin
            path = joinpath(DIR, "1d2d.svg")
            save_svg(path,M)
            @test isfile(path)
        end
        @testset "luxor-png" begin
            path = joinpath(DIR, "1d2d.png")
            save_png(path,M)
            @test isfile(path)
        end
    end

    @testset "2d2d" begin
        p = 3 # degree of polynomial
        k = KnotVector(1:12) # knot vector
        P = BSplineSpace{p}(k) # B-spline space
        rand_a = [randn(2) for i in 1:dim(P), j in 1:dim(P)] / 2
        a = [SVector(2i-6, 2j-6) for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
        M = BSplineManifold(a, (P, P)) # Define B-spline manifold

        @testset "luxor-svg" begin
            path = joinpath(DIR, "2d2d.svg")
            save_svg(path,M)
            @test isfile(path)
        end
        @testset "luxor-png" begin
            path = joinpath(DIR, "2d2d.png")
            save_png(path,M)
            @test isfile(path)
        end
        @testset "luxor-color-a-png" begin
            color(u1,u2) = rand(RGB)
            path = joinpath(DIR, "2d2d_color-a.png")
            save_png(path, M, color)
            @test isfile(path)
        end
        @testset "luxor-color-b-png" begin
            n1, n2 = dim.(bsplinespaces(M))
            colors = rand(RGB,n1,n2)
            path = joinpath(DIR, "2d2d_color-b.png")
            save_png(path, M, colors)
            @test isfile(path)
        end
    end

    @testset "1d3d" begin
        p = 4 # degree of polynomial
        k = KnotVector(rand(12)) # knot vector
        P = BSplineSpace{p}(k) # B-spline space
        rand_a = [randn(3) for i in 1:dim(P)]
        a = [SVector(2i-6, 2i-6, 0) for i in 1:dim(P)] + rand_a # random generated control points
        M = BSplineManifold(a, (P,)) # Define B-spline manifold

        @testset "povray" begin
            path = joinpath(DIR, "1d3d.inc")
            save_pov(path,M)
            @test isfile(path)
        end
    end

    @testset "2d3d" begin
        p = 3 # degree of polynomial
        k = KnotVector(1:12) # knot vector
        P = BSplineSpace{p}(k) # B-spline space
        rand_a = [randn(3) for i in 1:dim(P), j in 1:dim(P)] / 2
        a = [SVector(2i-6, 2j-6, 0) for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
        M = BSplineManifold(a, (P, P)) # Define B-spline manifold

        @testset "povray" begin
            path = joinpath(DIR, "2d3d.inc")
            save_pov(path, M , maincolor=RGBA(0,1,1,0.95))
            @test isfile(path)
        end
    end

    @testset "3d3d" begin
        p = 3 # degree of polynomial
        k = KnotVector(1:12) # knot vector
        P = BSplineSpace{p}(k) # B-spline space
        rand_a = [randn(3) for i in 1:dim(P), j in 1:dim(P), k in 1:dim(P)] / 4
        a = [SVector(2i-6, 2j-6, 2k-6) for i in 1:dim(P), j in 1:dim(P), k in 1:dim(P)] + rand_a # random generated control points
        M = BSplineManifold(a, (P, P, P)) # Define B-spline manifold

        @testset "povray" begin
            path = joinpath(DIR, "3d3d.inc")
            save_pov(path, M, maincolor=RGBA(0.5,0.1,1,0.8))
            @test isfile(path)
        end
    end
end