text:
using Test

using FlightMechanicsUtils

using FlightMechanicsSimulator

# Create different dynamic systems
x_sixdofaeroeuler = SixDOFAeroEuler([
    502.0 * FT2M,
    0.2392628,
    0.0005061803,
    1.366289,
    0.05000808,
    0.2340769,
    -0.01499617,
    0.2933811,
    0.06084932,
    0.0 * FT2M,
    0.0 * FT2M,
    0.0 * FT2M,
    64.12363,
])

x_sixdofbodyeuler = SixDOFBodyEuler([
    wind2body(502.0 * FT2M, 0, 0,  0.2392628, 0.0005061803)...,
    1.366289,
    0.05000808,
    0.2340769,
    -0.01499617,
    0.2933811,
    0.06084932,
    0.0 * FT2M,
    0.0 * FT2M,
    0.0 * FT2M,
    64.12363,
])

DSS_ARR = [x_sixdofaeroeuler, x_sixdofbodyeuler]

# Test the interface for each system.
@testset "Interface $(typeof(dss))" for dss in DSS_ARR

    # Check contruction of a DSState type from any other DSState type (including itself)
    @testset "Constructor from $(typeof(dss2))" for dss2 in DSS_ARR
        @test hasmethod(typeof(dss), (typeof(dss2),))
        # TODO: check that DSState are equivalent
    end

    # Contructor for trimmer
    @testset "Constructor for Trimmer" begin
        @test hasmethod(
            typeof(dss),
            (
                FlightMechanicsSimulator.TrimConditions,
                FlightMechanicsSimulator.TrimSolution,
                FlightMechanicsSimulator.Aircraft,
            ),
        )
    end

    methods_scalar_output = [
        get_height,
        get_tas,
        get_α,
        get_β,
        get_engine_power,
    ]

    # Check that methods that return a scalar exist and have right return type
    @testset "$im" for im in methods_scalar_output
        @test hasmethod(im, (typeof(dss),))
        @test isa(im(dss), Number)
    end

    methods_vec3_output = [
        get_earth_position,
        get_tasαβ,
        get_euler_angles,
        get_body_velocity,
        get_horizon_velocity,
        get_ang_vel_body,
        get_euler_angles_rates,
    ]

    # Check that methods that return a len 3 vector exist and have right return type
    @testset "$im" for im in methods_vec3_output
        @test hasmethod(im, (typeof(dss),))
        rv = im(dss)
        @test isa(rv, AbstractArray)
        @test length(rv) == 3
        @test eltype(rv) <: Number
    end

    # Check get_n_states implementation
    @testset "$im" for im in [get_n_states]
        @test hasmethod(im, (typeof(dss),))
        rv = im(dss)
        @test isa(rv, Integer)
        @test length(dss.x) == rv
    end

    # Check get_x implementation
    @testset "$im" for im in [get_x]
        @test hasmethod(im, (typeof(dss),))
        rv = im(dss)
        @test isa(rv, AbstractArray)
        @test length(rv) == get_n_states(dss)
        @test eltype(rv) <: Number
    end

    # Check get_x_names implementation
    @testset "$im" for im in [get_x_names]
        @test hasmethod(im, (typeof(dss),))
        rv = im(dss)
        @test isa(rv, AbstractArray)
        @test length(rv) == get_n_states(dss)
        @test eltype(rv) == Symbol
    end

    # Check state_eqs implementation
    @testset "state_eqs" begin
        @test hasmethod(
            state_eqs,
            (
                typeof(dss),
                Number,
                Number,
                AbstractMatrix,
                AbstractVector,
                AbstractVector,
                AbstractVector,
                Number,
            ),
        )

        time = 0.0
        mass = 10500.0
        inertia = [555.5 0.0 1234.5; 0.0 5674.0 0.0; 987.2 0.0  7895.3]
        forces = [1.0, 2.0, 3.0]
        moments = [4.0, 5.0, 6.0]
        hx = [12.0, 20.0, 30.0]
        pdot = 12.5

        dssd = state_eqs(dss, time, mass, inertia, forces, moments, hx, pdot)

        @test isa(dssd, DSStateDot)
    end

    # Check DSStateDot for each subtype: constructor and mandatory methods.
    @testset "DSStateDot{$(typeof(dss)), $(get_n_states(dss)), $(eltype(dss))}" begin
        # Test constructor
        dssd = DSStateDot(dss, get_x(dss))

        @testset "$im" for im in [get_ds_state]
            @test hasmethod(im, (typeof(dssd),))
            rv = im(dssd)
            @test isa(rv, typeof(dss))
        end

        @testset "$im" for im in [get_xdot]
            @test hasmethod(im, (typeof(dssd),))
            rv = im(dssd)
            @test isa(rv, AbstractArray)
            @test length(rv) == get_n_states(dss)
            @test eltype(rv) <: Number
        end

        @testset "$im" for im in [get_n_states]
            @test hasmethod(im, (typeof(dss),))
            rv = im(dssd)
            @test isa(rv, Integer)
            @test length(dss.x) == rv
        end

        @testset "$im" for im in [get_x]
            @test hasmethod(im, (typeof(dss),))
            rv = im(dssd)
            @test isa(rv, AbstractArray)
            @test length(rv) == get_n_states(dss)
            @test eltype(rv) <: Number
        end

        @testset "$im" for im in [get_x_names]
            @test hasmethod(im, (typeof(dss),))
            rv = im(dssd)
            @test isa(rv, AbstractArray)
            @test length(rv) == get_n_states(dss)
            @test eltype(rv) == Symbol
        end

        methods_scalar_output = [
            get_height,
            get_tas,
            get_α,
            get_β,
            get_engine_power,
            get_tas_dot,
            get_α_dot,
            get_β_dot,
            get_engine_power_dot,
        ]

        @testset "$im" for im in methods_scalar_output
            @test hasmethod(im, (typeof(dssd),))
            @test isa(im(dssd), Number)
        end

        methods_vec3_output = [
            get_earth_position,
            get_tasαβ,
            get_euler_angles,
            get_body_velocity,
            get_horizon_velocity,
            get_ang_vel_body,
            get_euler_angles_rates,
            get_tasαβ_dot,
            get_accel_body,
            get_ang_vel_body,
        ]

        @testset "$im" for im in methods_vec3_output
            @test hasmethod(im, (typeof(dssd),))
            rv = im(dssd)
            @test isa(rv, AbstractArray)
            @test length(rv) == 3
            @test eltype(rv) <: Number
        end

        @testset "$im" for im in [get_trimmer_cost]
            @test hasmethod(im, (typeof(dssd),))
            rv = im(dssd)
            @test length(rv) == 6
        end
    end
end


text:
## FUNTIONS TO LOAD MODEL IMPLEMENTATION CODE

"""
load(name::String; pkg=nothing, modl=Main, verbosity=1, allow_ambiguous=false)

Load the model implementation code for the model type with specified
`name` into the module `modl`, specifying `pkg` if necesssary, to
resolve duplicate names.

    load(proxy; pkg=nothing, modl=Main, verbosity=1)

In the case that `proxy` is a return value of `traits` (ie, has the
form `(name = ..., package_name = ..., etc)`) this is equivalent to
`load(name=proxy.name, pkg=proxy.package_name)`.

If `allow_ambiguous=true` then multiple models with the same name can be imported, otherwise a warning is raised and the previously available model is kept.

See also [`@load`](@ref)

### Examples

    load("ConstantClassifier")
    load(localmodels()[1])

See also [`@load`](@ref)

"""
function load(proxy::ModelProxy; modl=Main, verbosity=0, allow_ambiguous=false)
    # get name, package and load path:
    name = proxy.name
    pkg = proxy.package_name
    handle = (name=name, pkg=pkg)

    path = INFO_GIVEN_HANDLE[handle][:load_path]
    path_components = split(path, '.')

    # decide what to print
    toprint = verbosity > 0

    # return if model is already loaded unless allow_ambiguous
    localnames = map(p->p.name, localmodels(modl=modl))
    if !allow_ambiguous && name ∈ localnames
        @info "A model type \"$name\" is already loaded. \n"*
        "No new code loaded. "
        return
    end

    verbosity > 1 && @info "Loading into module \"$modl\": "

    # if needed, put MLJModels in the calling module's namespace:
    if path_components[1] == "MLJModels"
        toprint && print("import MLJModels ")
        modl.eval(:(import MLJModels))
        toprint && println('\u2714')
    end

    # load the package (triggering lazy-load of implementation code if
    # this is in MLJModels):
    pkg_ex = Symbol(pkg)
    toprint && print("import $pkg_ex ")
    modl.eval(:(import $pkg_ex))

    toprint && println('\u2714')

    # load the model:
    load_ex = Meta.parse("import $path")
    toprint && print(string(load_ex, " "))
    modl.eval(load_ex)
    toprint && println('\u2714')

    nothing
end

load(name::String; pkg=nothing, kwargs...) =
    load(info(name; pkg=pkg); kwargs...)

"""
    @load name pkg=nothing verbosity=0

Load the model implementation code for the model with specified `name`
into the calling module, provided `pkg` is specified in the case of
duplicate names.

### Examples

    @load DecisionTreeeRegressor
    @load PCA verbosity=1
    @load SVC pkg=LIBSVM

See also [`load`](@ref)

"""
macro load(name_ex, kw_exs...)
    name_ = string(name_ex)

    # parse kwargs:
    warning = "Invalid @load syntax.\n "*
    "Sample usage: @load PCA pkg=\"MultivariateStats\" verbosity=1"
    for ex in kw_exs
        ex.head == :(=) || throw(ArgumentError(warning))
        variable_ex = ex.args[1]
        value_ex = ex.args[2]
        if variable_ex == :pkg
            pkg = string(value_ex)
        elseif variable_ex == :verbosity
            verbosity = value_ex
        else
            throw(ArgumentError(warning))
        end
    end
    (@isdefined pkg) || (pkg = nothing)
    (@isdefined verbosity) || (verbosity = 0)

    # get rid brackets in name_, as in
    # "(MLJModels.Clustering).KMedoids":
    name = filter(name_) do c !(c in ['(',')']) end

    load(name, modl=__module__, pkg=pkg, verbosity=verbosity)

    esc(quote
            try
                $name_ex()
            catch
                try # hack for baremodules that have imported Base.eval:
                    $name_ex = Base.$name_ex
                    $name_ex()
                catch
                    @warn "Code is loaded but no instance returned. "
                    nothing
                end
            end
        end)
end

ext: jl


text:
if !haskey(ENV, "BUILDKITE")
    import Pkg
    Pkg.develop(Pkg.PackageSpec(; path = dirname(@__DIR__)))
end

include("TemperatureProfiles.jl")
include("relations.jl")


text:
module BasisMatricesTests

using BasisMatrices
using Test, LinearAlgebra, SparseArrays, Statistics

tests = ["types.jl", "basis.jl", "util.jl", "spline.jl", "interp.jl",
         "cheb.jl", "lin.jl", "basis_structure.jl", "complete.jl",
         "spline_sparse.jl", "smol.jl"]

if length(ARGS) > 0
    tests = ARGS
end

end_jl(s) = endswith(s, ".jl") ? s : s * ".jl"

for t in tests
    printstyled("* $t\n", color=:green)
    include(end_jl(t))
end

end  # module


text:
s = open("input.txt") do file
    readlines(file);
end

nums = map(x-> parse(Int, x), s);

invalidNum = 248131121;

function FindRange(nums, invalidNum)
    rangeStart = 1;
    rangeEnd = 2;
    curSum = sum(nums[1:2]);
    next = 3;

    while rangeEnd <= length(nums)
        rangeEnd += 1;
        curSum += nums[rangeEnd];

        while curSum > invalidNum
            curSum -= nums[rangeStart];
            rangeStart += 1;
        end

        if curSum == invalidNum
            return nums[rangeStart:rangeEnd]
        end
    end
end

sumRange = FindRange(nums, invalidNum);
println(minimum(sumRange) + maximum(sumRange))


text:
#__precompile__() # this module is safe to precompile
module NeidSolarScripts
#using PyCall

#=
function __init__()
  pyimport("astropy")
end
=#

using Dates, LinearAlgebra

# For reading data from pyrheliometer files
include("pyroheliometer.jl")
export Pyroheliometer

#=
# For computing effects of apparent solar rotation rate
include("solar_rotation.jl")
export SolarRotation
export CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic, CalculateFWHMDifference_SolarRotation_from_loc_Equatorial
export CalculateFWHMDifference_SolarRotation_from_long_lat_alt, CalculateFWHMDifference_SolarRotation_from_obs
#export get_wiyn_loc, get_harpsn_loc

# For computing differential extinction
include("diff_solar_extinction.jl")
export DifferentialExtinction
=#

#include("continuum_AFS.jl")
#export Continuum_AFS
include("continuum_rassine_like.jl")
export Continuum

end


text:
export ELFSegments, ELFSegment, ELFSegment32, ELFSegment64

"""
    ELFSegments

ELF segment header table type, containing information about the number of
segments within the ELF object, the location of the segment headers, etc...
"""
struct ELFSegments{H <: ELFHandle} <: Segments{H}
    handle::H
end
Segments(oh::ELFHandle) = ELFSegments(oh)
handle(segs::ELFSegments) = segs.handle
lastindex(segs::ELFSegments) = header(handle(segs)).e_phnum

"""
    ELFSegment

ELF Segment type, also known as a program header, containing information
about the virtual memory layout of a chunk of the program.
"""
abstract type ELFSegment{H <: ELFHandle} <: Segment{H} end

@io struct ELFSegment32{H <: ELFHandle} <: ELFSegment{H}
    p_type::UInt32
    p_offset::UInt32
    p_vaddr::UInt32
    p_paddr::UInt32
    p_filesz::UInt32
    p_memsz::UInt32
    p_flags::UInt32
    p_align::UInt32
end

@io struct ELFSegment64{H <: ELFHandle} <: ELFSegment{H}
    p_type::UInt32
    p_flags::UInt32
    p_offset::UInt64
    p_vaddr::UInt64
    p_paddr::UInt64
    p_filesz::UInt64
    p_memsz::UInt64
    p_align::UInt64
end

# Construction
function read(oh::ELFHandle, ::Type{ST}) where {ST <: ELFSegment}
    type_func = segment_header_type
    size_func = segment_header_size
    return read_struct(oh, type_func, size_func, "Segment Header")
end

function contents(oh::ELFHandle, seg::ELFSegment)
    # Allocate memory to hold this segment
    data = Array(UInt8, segment_file_size(seg))

    # Seek to where this segment starts
    seek(oh, segment_offset(seg))

    # Read and return it
    read(oh, data)
    return data
end

# Accessors
segment_type(seg::ELFSegment) = seg.p_type
segment_offset(seg::ELFSegment) = seg.p_offset
segment_file_size(seg::ELFSegment) = seg.p_filesz
segment_memory_size(seg::ELFSegment) = seg.p_memsz
segment_physical_address(seg::ELFSegment) = seg.p_paddr
segment_virtual_address(seg::ELFSegment) = seg.p_vaddr
segment_address(seg::ELFSegment) = segment_virtual_address(seg)

"""
    ELFSegmentRef

Provides a reference to a `Segment`, along with a reference to the
`ObjectHandle` this `Segment` comes from.  This should be the primary method by
which users interact with segments inside object files.  The list of available
API operations is given below.  Note that this overlaps heavily with the
`Segment` object API, this is by design as many of the methods are simply
passthroughs to the underlying `Segment` API calls for ease of use.

### Creation:
  - SegmentRef()

### Utility
  - deref()
  - handle()
  - Segments()

### IO-like operations:
  - read()
  - seekstart()
  - seek()
  - eof()

### Format-specific properties:
  - segment_type()
  - segment_offset()
  - segment_file_size()
  - segment_memory_size()
  - segment_physical_address()
  - segment_virtual_address()
  - segment_number()
"""
struct ELFSegmentRef{H<:ELFHandle} <: SegmentRef{H}
    segments::ELFSegments{H}
    segment::ELFSegment{H}
    idx::UInt32
end

function SegmentRef(segments::ELFSegments{H},
                    segment::ELFSegment{H}, idx) where {H <: ELFHandle}
    return ELFSegmentRef(segments, segment, UInt32(idx))
end
deref(seg::ELFSegmentRef) = seg.segment
Segments(seg::ELFSegmentRef) = seg.segments

segment_number(seg::ELFSegmentRef) = seg.idx
@derefmethod segment_physical_address(seg::ELFSegmentRef)
@derefmethod segment_virtual_address(seg::ELFSegmentRef)
@derefmethod segment_type(seg::ELFSegmentRef)

# Printing
"""
    segment_type_string(s::ELFSegment)

Return the given `ELFSegment`'s segment type as a string.
"""
function segment_type_string(s::ELFSegment)
    global P_TYPE

    s_type = segment_type(s)
    if haskey(P_TYPE, s_type)
        return P_TYPE[s_type]
    end
    return string("Unknown Segment Type (0x", string(s_type, base=16), ")")
end
@derefmethod segment_type_string(s::ELFSegmentRef)

function show(io::IO, seg::Union{ELFSegment,ELFSegmentRef})
    print(io, "ELFSegment")

    if !get(io, :compact, false)
        fsz = string(segment_file_size(seg), base=16)
        msz = string(segment_memory_size(seg), base=16)
        paddr = string(segment_physical_address(seg), base=16)
        vaddr = string(segment_virtual_address(seg), base=16)

        println(io)
        println(io, "       Type: $(segment_type_string(seg))")
        println(io, "       Size: (File: 0x$(fsz),   Mem: 0x$(msz))")
        println(io, "     Offset: 0x$(string(segment_offset(seg), base=16))")
        print(io,   "    Address:  (Phy: 0x$(paddr) Virt: 0x$(vaddr))")
    else
        print(io, " $(segment_type_string(seg))")
    end
end


text:
export SimElimination
"""
    SimElimination

 Simulates multiple elimination competitions between a list of players under some performance model.

## Parameters
* `n::Int`: number of n of RR competition to run

## Limitations
* no ties
* no additional factors in competitions
* just win/loss, so scores are trivial
* random ordering of players (no match making)

```
"""
struct SimElimination <: SimulateRule
    n::Int
end
SimElimination(; n::Int=1) = SimElimination(n)

function simulate( r::RatingsList, model::SimElimination, perf_model::ContinuousUnivariateDistribution )
    m = length( r.players )
    L = m-1 # number of matches to play per tournament
    df = DataFrame([String, String, Int, Int, Int], [PlayerA, PlayerB, Outcome, ScoreA, ScoreB], model.n*L)

    # simulate all players playing each other once per round
    c = 1
    winner = ""
    for i=1:model.n
        # random ordering of players
        current_player_list = shuffle( r.players )
        println(" $current_player_list")
        not_end = true
        level = 1
        while not_end
            k = 1
            new_player_list = Array{String,1}(undef, Int(ceil(length(current_player_list)/2)) )
            while k<length(current_player_list)
                # println("  level=$level, k=$k, c=$c, l=$(length(current_player_list))")
                pA = current_player_list[k]
                pB = current_player_list[k+1]
                df[c, PlayerA] = pA
                df[c, PlayerB] = pB
                o = outcome( r.ratings[pA], r.ratings[pB],  perf_model )
                df[c, Outcome] = o
                if o == 1
                    df[c, ScoreA] = 1
                    df[c, ScoreB] = 0
                    new_player_list[ Int(ceil(k/2)) ] = pA
                elseif o == -1
                    df[c, ScoreA] = 0
                    df[c, ScoreB] = 1
                    new_player_list[ Int(ceil(k/2)) ] = pB
                end
                c += 1
                k += 2
            end
            if k==length(current_player_list)
                new_player_list[end] = current_player_list[end]
            end
            # current_player_list = shuffle(new_player_list)
            current_player_list = new_player_list
            if length(current_player_list) == 1
                not_end = false
            end
            level += 1
        end
        winner = current_player_list[1] # winner is last person standing
    end
    return df, winner
end

function outcome( rA::Float64, rB::Float64, perf_model::ContinuousUnivariateDistribution)
    p = cdf( perf_model, rA - rB )
    x = rand()
    if x <= p
        o = 1
    else
        o = -1
    end
    return o
end


text:
using HSA
using HSA.Builtins: get_global_id

code_llvm(get_global_id, (Int32,))

@target hsail function vcopy(a,b)
    i = get_global_id(Int32(0))

	x = b[i+1]
	a[i+1] = x
	return nothing
end

println("------- LLVM -------------------------------")
code_llvm(vcopy, (Ptr{Int64}, Ptr{Int64}))

println("------- SPIR -------------------------------")
code_spir(vcopy, (Ptr{Int64}, Ptr{Int64}))

println("------- HSA --------------------------------")
code_hsa(vcopy, (Ptr{Int}, Ptr{Int}))

println("------- HSAIL ------------------------------")
println(src_hsail(vcopy, (Ptr{Int}, Ptr{Int})))


text:
# This file was generated, do not modify it. # hide
boston3 = coerce(boston, :Chas => OrderedFactor, :Rad => OrderedFactor);


text:
const kDoubleSignificandSize = 53

function FillDigits32FixedLength(n1,requested_len,buffer,len)
    for i = (requested_len-1):-1:0
        buffer[len+i] = 0x30 + n1 % 10
        n1 = div(n1,10)
    end
    return len + requested_len
end

function FillDigits32(n,buffer,len)
    n_len = 0
    while n != 0
        digit = n % 10
        n = div(n,10)
        buffer[len+n_len] = 0x30 + digit
        n_len += 1
    end
    i,j = len, len + n_len - 1
    while i < j
        buffer[i], buffer[j] = buffer[j], buffer[i]
        i += 1
        j -= 1
    end
    return len + n_len
end

function FillDigits64FixedLength(n2,buffer,len)
    kTen7 = 10000000
    part2 = n2 % kTen7
    n2 = div(n2,kTen7)
    part0, part1 = divrem(n2,kTen7)
    len = FillDigits32FixedLength(part0, 3, buffer, len)
    len = FillDigits32FixedLength(part1, 7, buffer, len)
    len = FillDigits32FixedLength(part2, 7, buffer, len)
    return len
end

function FillDigits64(n3,buffer,len)
    kTen7 = 10000000
    part2 = n3 % kTen7
    n3 = div(n3,kTen7)
    part0, part1 = divrem(n3,kTen7)
    if part0 != 0
        len = FillDigits32(part0, buffer, len)
        len = FillDigits32FixedLength(part1, 7, buffer, len)
        len = FillDigits32FixedLength(part2, 7, buffer, len)
    elseif part1 != 0
        len = FillDigits32(part1, buffer, len)
        len = FillDigits32FixedLength(part2, 7, buffer, len)
    else
        len = FillDigits32(part2, buffer, len)
    end
    return len
end

function RoundUp(buffer, len, decimal_point)
    if len == 1
        buffer[1] = 0x31
        decimal_point = 1
        len = 2
        return len, decimal_point
    end
    buffer[len - 1] += 1
    for i = (len-1):-1:2
        buffer[i] != 0x30 + 10 && return len, decimal_point
        buffer[i] = 0x30
        buffer[i - 1] += 1
    end
    if buffer[1] == 0x30 + 10
        buffer[1] = 0x31
        decimal_point += 1
    end
    return len, decimal_point
end

function FillFractionals(fractionals, exponent,
                         fractional_count, buffer,
                         len, decimal_point)
    if -exponent <= 64
        point = -exponent
        for i = 1:fractional_count
            fractionals == 0 && break
            fractionals *= 5
            point -= 1
            digit = fractionals >> point
            buffer[len] = 0x30 + digit
            len += 1
            fractionals -= uint64(digit) << point
        end
        if ((fractionals >> (point - 1)) & 1) == 1
            len, decimal_point = RoundUp(buffer, len, decimal_point)
        end
    else
        fract128 = uint128(fractionals) << 64
        fract128 = shift(fract128,-exponent - 64)
        point = 128
        for i = 1:fractional_count
            fract128 == 0 && break
            fract128 *= 5
            point -= 1
            digit, fract128 = divrem2(fract128,point)
            buffer[len] = 0x30 + digit
            len += 1
        end
        if bitat(fract128,point - 1) == 1
            len, decimal_point = RoundUp(buffer, len, decimal_point)
        end
    end
    return len, decimal_point
end

low(x) = uint64(x)
high(x) = uint64(x >> 64)
bitat(x::Uint128,y) = y >= 64 ? (int32(high(x) >> (y-64)) & 1) : (int32(low(x) >> y) & 1)
function divrem2(x,power)
    h = high(x)
    l = low(x)
    if power >= 64
        result = int32(h >> (power - 64))
        h -= uint64(result) << (power - 64)
        return result, (uint128(h) << 64) + l
    else
        part_low::Uint64 = l >> power
        part_high::Uint64 = h << (64 - power)
        result = int32(part_low + part_high)
        return result, uint128(l - (part_low << power))
    end
end
function shift(x::Uint128,amt)
    if amt == 0
      return x
    elseif amt == -64
        return x << 64
    elseif amt == 64
        return x >> 64
    elseif amt <= 0
        h = high(x); l = low(x)
        h <<= -amt
        h += l >> (64 + amt)
        l <<= -amt
        return (uint128(h) << 64) + l
    else
        h = high(x); l = low(x)
        l >>= amt
        l += h << (64 - amt)
        h >>= amt
        return (uint128(h) << 64) + l
    end
end

function TrimZeros(buffer, len, decimal_point)
    while len > 1 && buffer[len - 1] == 0x30
        len -= 1
    end
    first_non_zero::Int32 = 1
    while first_non_zero < len && buffer[first_non_zero] == 0x30
        first_non_zero += 1
    end
    if first_non_zero != 1
        for i = first_non_zero:(len-1)
            buffer[i - first_non_zero + 1] = buffer[i]
        end
        len -= first_non_zero-1
        decimal_point -= first_non_zero-1
    end
    return len, decimal_point
end

function fastfixedtoa(v,mode,fractional_count,buffer)
    v = float64(v)
    significand::Uint64 = _significand(v)
    exponent = _exponent(v)
    exponent > 20 && return false, 0, 0, buffer
    fractional_count > 20 && return false, 0, 0, buffer
    len = 1
    if exponent + kDoubleSignificandSize > 64
        kFive17 = divisor = 5^17
        divisor_power = 17
        dividend = significand
        if exponent > divisor_power
            dividend <<= exponent - divisor_power
            quotient = div(dividend,divisor)
            remainder = (dividend % divisor) << divisor_power
        else
            divisor <<= divisor_power - exponent
            quotient = div(dividend,divisor)
            remainder = (dividend % divisor) << exponent
        end
        len = FillDigits32(quotient, buffer, len)
        len = FillDigits64FixedLength(remainder, buffer, len)
        decimal_point = len-1
    elseif exponent >= 0
        significand <<= exponent
        len = FillDigits64(significand, buffer, len)
        decimal_point = len-1
    elseif exponent > -kDoubleSignificandSize
        integrals = significand >> -exponent
        fractionals = significand - (integrals << -exponent)
        if integrals > 0xFFFFFFFF
            len = FillDigits64(integrals,buffer,len)
        else
            len = FillDigits32(uint32(integrals),buffer,len)
        end
        decimal_point = len-1
        len, decimal_point = FillFractionals(fractionals,exponent,fractional_count,
                                                    buffer,len, decimal_point)
    elseif exponent < -128
        len = 1
        decimal_point = -fractional_count
    else
        decimal_point = 0
        len, decimal_point = FillFractionals(significand,exponent,fractional_count,
                                                    buffer,len, decimal_point)
    end
    len, decimal_point = TrimZeros(buffer,len,decimal_point)
    buffer[len] = 0
    if (len-1) == 0
        decimal_point = -fractional_count
    end
    return true, len, decimal_point, buffer
end


text:
"""
    metadata(x)

Returns metadata associated with `x`
"""
metadata(x) = _metadata(parent_type(x), x)
metadata(x::AbstractDict) = x
metadata(x::NamedTuple) = x
_metadata(::Type{P}, x::T) where {P,T} = metadata(parent(x))
_metadata(::Type{T}, x::T) where {T} = no_metadata

"""
    metadata(x::AbstractArray; dim)

Returns the metadata associated with dimension `dim` of `x`.
"""
function metadata(x::AbstractArray; dim=nothing)
    if dim === nothing
        return _metadata(parent_type(x), x)
    else
        return _metadata_dim(x, to_dims(x, dim))
    end
end
_metadata_dim(x, dim::Int) = metadata(axes(x, dim))
_metadata_dim(x, dim::StaticInt{D}) where {D} = metadata(axes(x, dim))
_metadata_dim(x::LinearIndices, dim::Int) = metadata(getfield(x.indices, dim))
_metadata_dim(x::LinearIndices, dim::StaticInt{D}) where {D} = metadata(getfield(x.indices, D))
_metadata_dim(x::CartesianIndices, dim::Int) = metadata(getfield(x.indices, dim))
_metadata_dim(x::CartesianIndices, dim::StaticInt{D}) where {D} = metadata(getfield(x.indices, D))

"""
    metadata(x, k)

Returns the value associated with key `k` of `x`'s metadata.
"""
function metadata(x, k)
    out = getmeta(x, k, no_metadata)
    out === no_metadata && throw(KeyError(k))
    return out
end
function metadata(m::Module)
    if isdefined(m, GLOBAL_METADATA)
        return getfield(m, GLOBAL_METADATA)::GlobalMetadata
    else
        return GlobalMetadata(m)
    end
end

"""
    metadata(x::AbstractArray, k; dim)

Returns the value associated with key `k` of `x`'s metadata.
"""
@inline function metadata(x::AbstractArray, k; dim=nothing)
    if dim === nothing
        if has_metadata(x)
            return metadata(metadata(x), k)
        else
            return no_metadata
        end
    else
        return metadata(metadata(x; dim=dim), k)
    end
end

"""
    metadata!(x::AbstractArray, k, val)

Set the value associated with key `k` of `x`'s metadata to `val`.
"""
@inline metadata!(x, k, val) = metadata!(metadata(x), k, val)
metadata!(x::AbstractDict, k, val) = setindex!(x, val, k)
metadata!(m::AbstractDict{String}, k::Symbol, val) = setindex!(m, val, String(k))
metadata!(m::AbstractDict{Symbol}, k::String, val) = setindex!(m, val, Symbol(k))

"""
    metadata!(x::AbstractArray, k, val; dim)

Set the value associated with key `k` of the metadata at dimension `dim` of `x` to `val`.
"""
metadata!(x::AbstractArray, k, val; dim=nothing) = metadata!(metadata(x; dim=dim), k, val)

"""
    metadata_type(::Type{<:AbstractArray}; dim)::Type

Returns the type of the metadata of `x`. If `dim` is specified then returns type of
metadata associated with dimension `dim`.
"""
metadata_type(x::AbstractArray; dim=nothing) = metadata_type(typeof(x); dim=dim)
@inline function metadata_type(::Type{T}; dim=nothing) where {T<:AbstractArray}
    if dim === nothing
        return _metadata_type(parent_type(T), T)
    else
        return _metadata_dim_type(T, to_dims(T, dim))
    end
end
_metadata_type(::Type{T}, ::Type{T}) where {T} = NoMetadata
_metadata_type(::Type{P}, ::Type{T}) where {P,T} = metadata_type(P)
_metadata_dim_type(::Type{T}, dim) where {T} = metadata_type(axes_types(T, dim))

"""
    getmeta(x, key, default)

Return the metadata associated with `key`, or return `default` if `key` is not found.
"""
@inline getmeta(x, k, default) = _getmeta(metadata(x), k, default)
_getmeta(m, k, default) = get(m, k, default)
@inline _getmeta(m::AbstractDict{String}, k::Symbol, default) = get(m, String(k), default)
@inline _getmeta(m::AbstractDict{Symbol}, k::String, default) = get(m, Symbol(k), default)
@inline _getmeta(m::NamedTuple, k::String, default) = get(m, Symbol(k), default)

"""
    getmeta(f::Function, x, key)

Return the metadata associated with `key`, or return `f(x)` if `key` is not found. Note that
this behavior differs from `Base.get(::Function, x, keys)` in that `getmeta` passes `x` to
`f` as an argument (as opposed to `f()`).
"""
@inline function getmeta(f::Union{Function,Type}, x, k)
    m = getmeta(metadata(x), k, no_metadata)
    if m === no_metadata
        return f(x)
    else
        return m
    end
end

"""
    getmeta!(x, key, default)

Return the metadata associated with `key`. If `key` is not found then `default` is returned
and stored at `key`.
"""
@inline getmeta!(x, k, default) = _getmeta!(metadata(x), k, default)
_getmeta!(m, k, default) = get!(m, k, default)
@inline _getmeta!(m::AbstractDict{String}, k::Symbol, default) = get!(m, String(k), default)
@inline _getmeta!(m::AbstractDict{Symbol}, k::String, default) = get!(m, Symbol(k), default)

"""
    getmeta!(f::Function, x, key)

Return the metadata associated with `key`. If `key` is not found then `f(x)` is returned
and stored at `key`. Note that this behavior differs from `Base.get!(::Function, x, keys)` in
that `getmeta!` passes `x` to `f` as an argument (as opposed to `f()`).
"""
@inline function getmeta!(f::Function, x, k)
    m = metadata(x)
    out = get(m, k, no_metadata)
    if out === no_metadata
        out = f(x)
        metadata!(m, k, out)
    end
    return out
end

"""
    metadata_type(::Type{T})

Returns the type of the metadata associated with `x`.
"""
metadata_type(x) = metadata_type(typeof(x))
metadata_type(::Type{T}) where {T} = _metadata_type(parent_type(T), T)
metadata_type(::Type{T}) where {T<:AbstractDict} = T
metadata_type(::Type{T}) where {T<:NamedTuple} = T
metadata_type(::Type{T}) where {T<:Module} = GlobalMetadata
metadata_type(::Type{MetaStruct{P,M}}) where {P,M} = M
@inline function metadata_type(::Type{T}; dim=nothing) where {M,A,T<:MetaArray{<:Any,<:Any,M,A}}
    if dim === nothing
        return M
    else
        return metadata_type(A; dim=dim)
    end
end
@inline function metadata_type(::Type{T}; dim=nothing) where {R,M,T<:MetaUnitRange{<:Any,R,M}}
    if dim === nothing
        return M
    else
        return metadata_type(R; dim=dim)
    end
end

"""
    has_metadata(x)::Bool

Returns `true` if `x` has metadata.
"""
has_metadata(x) = has_metadata(typeof(x))
has_metadata(::Type{T}) where {T} = _has_metadata(metadata_type(T))
_has_metadata(::Type{NoMetadata}) = false
_has_metadata(::Type{T}) where {T} = true

"""
    has_metadata(x::AbstractArray; dim)::Bool

Returns `true` if `x` has metadata associated with dimension `dim`.
"""
has_metadata(x::AbstractArray; dim=nothing) = has_metadata(typeof(x); dim=dim)
function has_metadata(::Type{T}; dim=nothing) where {T<:AbstractArray}
    if dim === nothing
        return _has_metadata(metadata_type(T))
    else
        return _has_metadata(metadata_type(T; dim=dim))
    end
end

"""
    has_metadata(x, k)::Bool

Returns `true` if metadata associated with `x` has the key `k`.
"""
has_metadata(x, k) = haskey(metadata(x), k)

"""
    has_metadata(x::AbstractArray, k; dim)::Bool

Returns `true` if metadata associated with dimension `dim` of `x` has the key `k`.
"""
has_metadata(x::AbstractArray, k; dim=nothing) = haskey(metadata(x; dim=dim), k)

"""
    attach_metadata(x, metadata)

Generic method for attaching metadata to `x`.
"""
attach_metadata(x, m=Dict{Symbol,Any}()) = MetaStruct(x, m)
attach_metadata(x::AbstractArray, m=Dict{Symbol,Any}()) = MetaArray(x, m)
function attach_metadata(x::AbstractRange, m=Dict{Symbol,Any}())
    if known_step(x) === oneunit(eltype(x))
        return MetaUnitRange(x, m)
    else
        return MetaArray(x, m)
    end
end
attach_metadata(x::IO, m=Dict{Symbol,Any}()) = MetaIO(x, m)
attach_metadata(m::METADATA_TYPES) = Base.Fix2(attach_metadata, m)

"""
    share_metadata(src, dst) -> attach_metadata(dst, metadata(src))

Shares the metadata from `src` by attaching it to `dst`.
The returned instance will have properties that are synchronized with `src` (i.e.
modifying one's metadata will effect the other's metadata).

See also: [`copy_metadata`](@ref).
"""
share_metadata(src, dst) = attach_metadata(dst, metadata(src))

"""
    copy_metadata(src, dst) -> attach_metadata(dst, copy(metadata(src)))

Copies the the metadata from `src` and attaches it to `dst`. Note that this method
specifically calls `deepcopy` on the metadata of `src` to ensure that changing the
metadata of `dst` does not affect the metadata of `src`.

See also: [`share_metadata`](@ref).
"""
copy_metadata(src, dst) = attach_metadata(dst, deepcopy(metadata(src)))

"""
    drop_metadata(x)

Returns `x` without metadata attached.
"""
drop_metadata(x) = parent(x)

# This allows dictionaries's keys to be treated like property names
@inline function metadata_keys(x::AbstractArray; dim=nothing)
    if has_metadata(x; dim=dim)
        return metadata_keys(metadata(x; dim=dim))
    else
        return propertynames(x)
    end
end

metadata_keys(x::AbstractDict) = keys(x)
metadata_keys(::NamedTuple{L}) where {L} = L
function metadata_keys(x::X) where {X}
    if has_metadata(X)
        return metadata_keys(metadata(x))
    else
        return propertynames(x)
    end
end

"""
    MetadataPropagation(::Type{T})

Returns type informing how to propagate metadata of type `T`.
See [`DropMetadata`](@ref), [`CopyMetadata`](@ref), [`ShareMetadata`](@ref).
"""
abstract type MetadataPropagation end

"""
    DropMetadata

Informs operations that may propagate metadata to insead drop it.
"""
struct DropMetadata <: MetadataPropagation end

"""
    CopyMetadata

Informs operations that may propagate metadata to attach a copy to any new instance created.
"""
struct CopyMetadata <: MetadataPropagation end

"""
    ShareMetadata

Informs operations that may propagate metadata to attach a the same metadata to
any new instance created.
"""
struct ShareMetadata <: MetadataPropagation end

MetadataPropagation(x) = MetadataPropagation(typeof(x))
MetadataPropagation(::Type{T}) where {T} = MetadataPropagation(metadata_type(T))
MetadataPropagation(::Type{T}) where {T<:AbstractDict} = ShareMetadata()
MetadataPropagation(::Type{T}) where {T<:NamedTuple} = ShareMetadata()
MetadataPropagation(::Type{T}) where {T<:NoMetadata} = DropMetadata()

function propagate_metadata(src, dst)
    return propagate_metadata(MetadataPropagation(src), src, dst)
end
propagate_metadata(::DropMetadata, src, dst) = dst
propagate_metadata(::ShareMetadata, src, dst) = share_metadata(src, dst)
propagate_metadata(::CopyMetadata, src, dst) = copy_metadata(src, dst)

function MetadataPropagation(::Type{T}) where {P,M,T<:MetaStruct{P,M}}
    if P <: Number
        return DropMetadata()
    else
        return ShareMetadata()
    end
end

"""
    metadata_summary([io], x)

Creates summary readout of metadata for `x`.
"""
metadata_summary(x) = metadata_summary(stdout, x)
function metadata_summary(io::IO, x)
    print(io, "$(lpad(Char(0x2022), 3)) metadata:")
    suppress = getmeta(x, :suppress, no_metadata)
    if suppress !== no_metadata
        suppress = metadata(x, :suppress)
        for k in metadata_keys(x)
            if k !== :suppress
                println(io)
                print(io, "     ")
                print(io, "$k")
                print(io, " = ")
                if in(k, suppress)
                    print(io, "<suppressed>")
                else
                    print(io, metadata(x, k))
                end
            end
        end
    else
        for k in metadata_keys(x)
            println(io)
            print(io, "     ")
            print(io, "$k")
            print(io, " = ")
            print(io, metadata(x, k))
        end
    end
end

# this is a currently informal way of changing how showarg displays metadata in
# the argument list. If someone makes a metadata type that's long or complex they
# may want to overload this.
#
# - used within Base.showarg for MetaArray
showarg_metadata(x) = "::$(metadata_type(x))"

function _construct_meta(meta::AbstractDict{Symbol}, kwargs::NamedTuple)
    for (k, v) in pairs(kwargs)
        meta[k] = v
    end
    return meta
end

function _construct_meta(meta, kwargs::NamedTuple)
    if isempty(kwargs)
        return meta
    else
        error("Cannot assign key word arguments to metadata of type $(typeof(meta))")
    end
end

"""
    attach_global_metadata(x, meta, m::Module)

Attach metadata `meta` to the object id of `x` (`objectid(x)`) in global metadata of
module `m`.
"""
function attach_global_metadata(x, meta::MDict, m::Module)
    setindex!(metadata(m), meta, objectid(x))
    if _attach_global_metadata(x, x, m)
        @warn "Cannot create finalizer for $(typeof(x)). Global dictionary must be manually deleted."
    end
    return meta
end

function attach_global_metadata(x, meta, m::Module)
    gm = Dict{Symbol,Any}()
    for (k,v) in pairs(meta)
        gm[k] = v
    end
    return attach_global_metadata(x, gm, m)
end

function _attach_global_metadata(x, xfield, m::Module)
    there_is_no_finalizer = true
    if ismutable(x)
        _assign_global_metadata_finalizer(x, xfield, m)
        there_is_no_finalizer = false
    elseif isstructtype(typeof(x))
        N = fieldcount(typeof(x))
        if N !== 0
            i = 1
            while (there_is_no_finalizer && i <= N)
                there_is_no_finalizer = _assign_global_metadata_finalizer(x, getfield(xfield, i), m)
                i += 1
            end
        end
    end
    return there_is_no_finalizer
end

function _assign_global_metadata_finalizer(x, xfield, m::Module)
    if Base.ismutable(xfield)
        finalizer(xfield) do _
            @async delete_global_metadata!(x, m)
        end
        return false
    else
        return true
    end
end

delete_global_metadata!(x, m::Module) = delete!(metadata(m), objectid(x))

"""
    global_metadata(x, m::Module)
    global_metadata(x, k::Symbol, m::Module)

Retreive metadata associated with the object id of `x` (`objectid(x)`) in stored in the
global metadata of module `m`. If the key `k` is specified only the value associated with
that key is returned.
"""
global_metadata(x, m::Module) = get(metadata(m), objectid(x), no_metadata)
global_metadata(x, k, m::Module) = getindex(global_metadata(x, m), k)

"""
    global_metadata!(x, k, val, m::Module)

Set the value of `x`'s global metadata in module `m` associated with the key `k` to `val`.
"""
global_metadata!(x, k, val, m::Module) = setindex!(global_metadata(x, m), val, k)

"""
    @attach_metadata(x, meta)

Attach metadata `meta` to the object id of `x` (`objectid(x)`) in the current module's
global metadata.

See also: [`GlobalMetadata`](@ref)
"""
macro attach_metadata(x, meta)
    return esc(:(Metadata.attach_global_metadata($x, $meta, @__MODULE__)))
end

"""
    @metadata(x[, k])

Retreive metadata associated with the object id of `x` (`objectid(x)`) in the current
module's global metadata. If the key `k` is specified only the value associated with
that key is returned.
"""
macro metadata(x)
    return esc(:(Metadata.global_metadata($(x), @__MODULE__)))
end

macro metadata(x, k)
    return esc(:(Metadata.global_metadata($(x), $(k), @__MODULE__)))
end

"""
    @metadata!(x, k, val)

Set the value of `x`'s global metadata associated with the key `k` to `val`.
"""
macro metadata!(x, k, val)
    return esc(:(Metadata.global_metadata!($x, $k, $val, @__MODULE__)))
end

"""
    @has_metadata(x)::Bool
    @has_metadata(x, k)::Bool

Does `x` have metadata stored in the curren modules' global metadata? Checks for the
presenece of the key `k` if specified.
"""
macro has_metadata(x)
    return esc(:(Metadata.global_metadata($x, @__MODULE__) !== Metadata.no_metadata))
end
macro has_metadata(x, k)
    return esc(:(haskey(Metadata.global_metadata($x, @__MODULE__), $k)))
end

"""
    @share_metadata(src, dst) -> @attach_metadata(@metadata(src), dst)

Shares the metadata from `src` by attaching it to `dst`. This assumes that metadata
for `src` is stored in a global dictionary (i.e. not part of `src`'s structure) and
attaches it to `dst` through a global reference within the module.

See also: [`@copy_metadata`](@ref), [`share_metadata`](@ref)
"""
macro share_metadata(src, dst)
    return esc(:(Metadata.attach_global_metadata($dst, Metadata.global_metadata($src, @__MODULE__), @__MODULE__)))
end

"""
    @copy_metadata(src, dst) -> attach_metadata(dst, copy(metadata(src)))

Copies the metadata from `src` by attaching it to `dst`. This assumes that metadata
for `src` is stored in a global dictionary (i.e. not part of `src`'s structure) and
attaches a new copy to `dst` through a global reference within the module.

See also: [`@share_metadata`](@ref), [`copy_metadata`](@ref)
"""
macro copy_metadata(src, dst)
    return esc(:(Metadata.attach_global_metadata($dst, deepcopy(Metadata.metadata($src)), @__MODULE__)))
end


text:
using Test
using YaoBlocks.AD
using YaoBlocks
using YaoArrayRegister

@testset "NoParams" begin
    c = chain(put(3, 2 => Ry(0.2)), put(3, 1 => Rx(0.5)))
    np = chain(NoParams(put(3, 2 => Ry(0.2))), put(3, 1 => Rx(0.5)))
    reg = rand_state(3)
    @test content(np[1]) == put(3, 2 => Ry(0.2))
    @test parameters(np) == [0.5]
    @test nparameters(np) == 1
    @test getiparams(np[1]) == ()
    @test niparams(np[1]) == 0
    @test mat(c) ≈ mat(np)
    @test apply!(copy(reg), c) ≈ apply!(copy(reg), np)
    @test np' == chain(put(3, 1 => Rx(-0.5)), NoParams(put(3, 2 => Ry(-0.2))))
    @test dispatch!(np, [0.3]) == chain(NoParams(put(3, 2 => Ry(0.2))), put(3, 1 => Rx(0.3)))
end


text:
#
module PccfMain

export
    Pccf,           # Done
    Adwin,          # In progress
    Adwin3,         # In progress
    BayesChp,       # Done
    BayesianBlocks, # Done
    BasicCounting   # In progress

include("./lib/pccf.jl")

end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


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

const name = "libcgal_julia"
const version = v"0.8.0"

# Collection of sources required to build CGAL
const sources = [
    GitSource("https://github.com/rgcv/libcgal-julia.git",
              "ce3b4312d8ba161d85ff6bd2c4aed788f5523613"),
    # julia binaries
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.1-linux-x86_64.tar.gz",
                  "faa707c8343780a6fe5eaf13490355e8190acf8e2c189b9e7ecbddb0fa2643ad"; unpack_target="julia-x86_64-linux-gnu"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-apple-darwin14.tar.gz",
                  "f2e5359f03314656c06e2a0a28a497f62e78f027dbe7f5155a5710b4914439b1"; unpack_target="julia-x86_64-apple-darwin14"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-w64-mingw32.tar.gz",
                  "c7b2db68156150d0e882e98e39269301d7bf56660f4fc2e38ed2734a7a8d1551"; unpack_target="julia-x86_64-w64-mingw32"),
]

# Dependencies that must be installed before this package can be built
const dependencies = [
    Dependency("CGAL_jll"),
    Dependency("libcxxwrap_julia_jll"),
]

# Bash recipe for building across all platforms
const script = raw"""
## pre-build setup
# exit on error
set -eu

## "find" julia
case $target in
  x86_64-linux-gnu)
    Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/julia-1.3.1
    ;;
  x86_64-apple-darwin14|x86_64-w64-mingw32)
    Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/juliabin
    ;;
esac

## configure build
cmake libcgal-julia*/ -B build \
  `# cmake specific` \
  -DCMAKE_TOOLCHAIN_FILE="$CMAKE_TARGET_TOOLCHAIN" \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_FIND_ROOT_PATH="$prefix" \
  -DCMAKE_INSTALL_PREFIX="$prefix" \
  `# tell libcxxwrap-julia where julia is` \
  -DJulia_PREFIX="$Julia_PREFIX"

## and away we go..
VERBOSE=ON cmake --build build --config Release --target install -- -j$nproc

install_license libcgal-julia*/LICENSE

# HACK: Apparently, this isn't a simple build system anymore..
case $target in
  *mingw32*) mv "$prefix/lib/"*.dll "$prefix/bin" ;;
esac
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
const platforms = [
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:x86_64),
] |> expand_cxxstring_abis

# The products that we will ensure are always built
const products = [
    LibraryProduct("libcgal_julia_exact", :libcgal_julia_exact),
    LibraryProduct("libcgal_julia_inexact", :libcgal_julia_inexact),
]

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

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

"""
    copy_extension_data(data, new_model::AbstractModel, model::AbstractModel)

Return a copy of the extension data `data` of the model `model` to the extension
data of the new model `new_model`. A method should be added for any JuMP
extension storing data in the `ext` field.
"""
function copy_extension_data end

"""
    ReferenceMap

Mapping between variable and constraint reference of a model and its copy. The
reference of the copied model can be obtained by indexing the map with the
reference of the corresponding reference of the original model.
"""
struct ReferenceMap
    model::Model
    index_map::MOIU.IndexMap
end

function Base.getindex(map::ReferenceMap, vref::VariableRef)
    return VariableRef(map.model, map.index_map[index(vref)])
end

function Base.getindex(map::ReferenceMap, cref::ConstraintRef)
    return ConstraintRef(map.model, map.index_map[index(cref)], cref.shape)
end

function Base.getindex(map::ReferenceMap, expr::GenericAffExpr)
    result = zero(expr)
    for (coef, var) in linear_terms(expr)
        add_to_expression!(result, coef, map[var])
    end
    result.constant = expr.constant
    return result
end

function Base.getindex(map::ReferenceMap, expr::GenericQuadExpr)
    aff = map[expr.aff]
    terms = [UnorderedPair(map[key.a], map[key.b]) => val for
     (key, val) in expr.terms]
    return GenericQuadExpr(aff, terms)
end

Base.getindex(map::ReferenceMap, val::AbstractArray) = getindex.(map, val)

Base.broadcastable(reference_map::ReferenceMap) = Ref(reference_map)

# Return a Boolean if the filtering function (1st argument) indicates that the whole value should
# be copied over.
_should_copy_complete_object(_, _) = true
function _should_copy_complete_object(
    filter_constraints::Function,
    value::ConstraintRef,
)
    return filter_constraints(value)
end
function _should_copy_complete_object(
    filter_constraints::Function,
    value::AbstractArray{T},
) where {T<:ConstraintRef}
    return all([filter_constraints(value[i]) for i in eachindex(value)])
end # all(filter_constraints.(value))

"""
    copy_model(model::Model; filter_constraints::Union{Nothing, Function}=nothing)

Return a copy of the model `model` and a [`ReferenceMap`](@ref) that can be used
to obtain the variable and constraint reference of the new model corresponding
to a given `model`'s reference. A [`Base.copy(::AbstractModel)`](@ref) method
has also been implemented, it is similar to `copy_model` but does not return
the reference map.

If the `filter_constraints` argument is given, only the constraints for which
this function returns `true` will be copied. This function is given a
constraint reference as argument.

## Note

Model copy is not supported in `DIRECT` mode, i.e. when a model is constructed
using the [`direct_model`](@ref) constructor instead of the [`Model`](@ref)
constructor. Moreover, independently on whether an optimizer was provided at
model construction, the new model will have no optimizer, i.e., an optimizer
will have to be provided to the new model in the [`optimize!`](@ref) call.

## Examples

In the following example, a model `model` is constructed with a variable `x` and
a constraint `cref`. It is then copied into a model `new_model` with the new
references assigned to `x_new` and `cref_new`.
```julia
model = Model()
@variable(model, x)
@constraint(model, cref, x == 2)

new_model, reference_map = copy_model(model)
x_new = reference_map[x]
cref_new = reference_map[cref]
```
"""
function copy_model(
    model::Model;
    filter_constraints::Union{Nothing,Function} = nothing,
)
    if mode(model) == DIRECT
        error(
            "Cannot copy a model in `DIRECT` mode. Use the `Model` ",
            "constructor instead of the `direct_model` constructor to be ",
            "able to copy the constructed model.",
        )
    end
    caching_mode = backend(model).mode
    new_model = Model(caching_mode = caching_mode)

    # At JuMP's level, filter_constraints should work with JuMP.ConstraintRef,
    # whereas MOI.copy_to's filter_constraints works with MOI.ConstraintIndex.
    function moi_filter_constraints(cref::MOI.ConstraintIndex)
        jump_cref = constraint_ref_with_index(model, cref)
        return filter_constraints(jump_cref)
    end
    filter = filter_constraints !== nothing ? moi_filter_constraints : nothing

    # Copy the MOI backend, note that variable and constraint indices may have
    # changed, the `index_map` gives the map between the indices of
    # `backend(model` and the indices of `backend(new_model)`.
    index_map = MOI.copy_to(
        backend(new_model),
        backend(model),
        copy_names = true,
        filter_constraints = filter,
    )

    new_model.optimize_hook = model.optimize_hook

    # TODO copy NLP data
    if model.nlp_data !== nothing
        error(
            "copy is not supported yet for models with nonlinear constraints",
            " and/or nonlinear objective function",
        )
    end

    reference_map = ReferenceMap(new_model, index_map)

    for (name, value) in object_dictionary(model)
        if _should_copy_complete_object(filter_constraints, value)
            try
                new_model[name] = getindex(reference_map, value)
            catch err
                if err isa MethodError
                    @warn(
                        "Skipping the copy of object `:$(name)` due to " *
                        "unsupported type $(typeof(value)). Please open a " *
                        "GitHub issue at https://github.com/jump-dev/JuMP.jl " *
                        "with this message.",
                    )
                else
                    rethrow(err)
                end
            end
        end
    end

    for (key, data) in model.ext
        new_model.ext[key] = copy_extension_data(data, new_model, model)
    end

    return new_model, reference_map
end

"""
    copy(model::AbstractModel)

Return a copy of the model `model`. It is similar to [`copy_model`](@ref)
except that it does not return the mapping between the references of `model`
and its copy.

## Note

Model copy is not supported in `DIRECT` mode, i.e. when a model is constructed
using the [`direct_model`](@ref) constructor instead of the [`Model`](@ref)
constructor. Moreover, independently on whether an optimizer was provided at
model construction, the new model will have no optimizer, i.e., an optimizer
will have to be provided to the new model in the [`optimize!`](@ref) call.

## Examples

In the following example, a model `model` is constructed with a variable `x` and
a constraint `cref`. It is then copied into a model `new_model` with the new
references assigned to `x_new` and `cref_new`.
```julia
model = Model()
@variable(model, x)
@constraint(model, cref, x == 2)

new_model = copy(model)
x_new = model[:x]
cref_new = model[:cref]
```
"""
function Base.copy(model::AbstractModel)
    new_model, _ = copy_model(model)
    return new_model
end

"""
    copy_conflict(model::Model)

Return a copy of the current conflict for the model `model` and a
[`ReferenceMap`](@ref) that can be used to obtain the variable and constraint
reference of the new model corresponding to a given `model`'s reference.

This is a convenience function that provides a filtering function for
[`copy_model`](@ref).

## Note

Model copy is not supported in `DIRECT` mode, i.e. when a model is constructed
using the [`direct_model`](@ref) constructor instead of the [`Model`](@ref)
constructor. Moreover, independently on whether an optimizer was provided at
model construction, the new model will have no optimizer, i.e., an optimizer
will have to be provided to the new model in the [`optimize!`](@ref) call.

## Examples

In the following example, a model `model` is constructed with a variable `x` and
two constraints `cref` and `cref2`. This model has no solution, as the two
constraints are mutually exclusive. The solver is asked to compute a conflict
with [`compute_conflict!`](@ref). The parts of `model` participating in the
conflict are then copied into a model `new_model`.
```julia
model = Model() # You must use a solver that supports conflict refining/IIS
# computation, like CPLEX or Gurobi
@variable(model, x)
@constraint(model, cref, x >= 2)
@constraint(model, cref2, x <= 1)

compute_conflict!(model)
if MOI.get(model, MOI.ConflictStatus()) != MOI.CONFLICT_FOUND
    error("No conflict could be found for an infeasible model.")
end

new_model, reference_map = copy_conflict(model)
```
"""
function copy_conflict(model::Model)
    filter_constraints =
        (cref) ->
            MOI.get(model, MOI.ConstraintConflictStatus(), cref) !=
            MOI.NOT_IN_CONFLICT
    new_model, reference_map =
        copy_model(model, filter_constraints = filter_constraints)
    return new_model, reference_map
end

# Calling `deepcopy` over a JuMP model is not supported, nor planned to be
# supported, because it would involve making a deep copy of the underlying
# solver (behind a C pointer).
function Base.deepcopy(::Model)
    return error(
        "`JuMP.Model` does not support `deepcopy` as the reference to the underlying solver cannot be deep copied, use `copy` instead.",
    )
end

function MOI.copy_to(dest::MOI.ModelLike, src::Model)
    if src.nlp_data !== nothing
        # Re-set the NLP block in-case things have changed since last
        # solve.
        MOI.set(src, MOI.NLPBlock(), _create_nlp_block_data(src))
    end
    return MOI.copy_to(dest, backend(src))
end

function MOI.copy_to(dest::Model, src::MOI.ModelLike)
    return MOI.copy_to(backend(dest), src)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Read TOML file

"""
    read_home_config(toml_file::String) -> Dict

Read home configuration from a TOML file.
"""
function read_home_config(toml_file::String)::Dict
    h = TOML.parsefile(toml_file)
    # confirm unique ID
    reserved = Set(["load"])
    ids = Set{String}()
    for id in get_all_load_ids(h)
        if id in ids
            error("The ID '$id' is not unique!")
        end
        if id in reserved
            error("The ID '$id' is reserved and should not be used!")
        end
        push!(ids, id)
    end
    return h
end

"""
    get_all_load_ids(h::Dict) -> Vector{String}

Get the ID of all loads in the home `h`.
"""
function get_all_load_ids(h::Dict)::Vector{String}
    return [l["id"] for l in Iterators.flatten((h["CR"], h["AD"], h["SI"], h["SU"]))]
end

"""
    validate_solution(s::AbstractDataFrame, h::AbstractDict)

Given a solution in `s` along with the scenario data, validate whether it is feasible.
If it is feasible, return the total cost of this solution on `s`.
"""
function validate_solution(s::AbstractDataFrame, h::AbstractDict)
    C_SI = C_SU = C_AD = C_CR = 0.0
    # shiftable constraints
    for l in h["SU"]
        id = l["id"]
        if !isapprox(sum(s[!, Symbol("z_", id)]), 1)    # z_ is stored as Float in `s`
            @show sum(s[!, Symbol("z_", id)])
            error("Infeasible solution in SU load $id")
        end
        C_SU += sum(s[!, Symbol("C_", id)])
    end
    for l in h["SI"]
        id = l["id"]
        if !isapprox(sum(s[!, Symbol("z_", id)]), l["L"])
            error("Infeasible solution in SI load $id")
        end
        C_SI += sum(s[!, Symbol("C_", id)])
    end
    # power balance
    consumption = sum(s[!, Symbol("P_", l["id"])] for l in Iterators.Flatten((h["AD"], h["SI"], h["SU"]))) +
                  s.Ps .+ sum(l["P"] for l in h["CR"])
    generation = s.Pb .+ s.PPV
    if !all(consumption .≈ generation)
        @show consumption generation
        error("Power not balanced!")
    end
    for l in h["AD"]
        C_AD += sum(s[!, Symbol("C_", l["id"])])
    end
    for l in h["CR"]
        C_CR += sum(l["P"] .* s.ρ)
    end
    ρ = s.ρ
    ρs = h["UtilityGrid"]["β"] .* ρ
    return C_SI + C_SU + C_AD + C_CR - sum(ρ .* s.PPV) + sum((ρ .- ρs) .* s.Ps)
end

"""
    compute_cost!(df::DataFrame, h::AbstractDict)->Float64

For home `h`, `df` contains price, PV, and the power of each load. Then, the net power exchange is
computed and stored into `df` as `Pb` and `Ps` columns. Finally, the stepwise electricity cost `Cele` is
calculated, and the total cost `C` is also computed, both stored into `df` in place.

The total cost of all steps is returned.
"""
function compute_cost!(df::DataFrame, h::AbstractDict)::Float64
    T = nrow(df)
    c = 0.0
    for l in h["CR"]
        df[!, "P_$(l["id"])"] = fill(l["P"], T)
    end
    ids = [l["id"] for l in Iterators.Flatten((h["CR"], h["AD"], h["SU"], h["SI"]))]
    Pnet = sum(df[!, "P_$id"] for id in ids) .- df.PPV
    df.Pnet = Pnet
    Pb = max.(Pnet, 0.0)
    Ps = -min.(Pnet, 0.0)
    ρb = df.ρ
    ρs = h["UtilityGrid"]["β"] .* ρb
    Cele = @. ρb * Pb - ρs * Ps
    df.Pb = Pb
    df.Ps = Ps
    df.Cele = Cele
    # discomfort cost
    Cdc = zeros(T)
    for l in h["AD"]
        id = l["id"]
        Pmax = l["Pmax"]
        α = l["α"]
        df[!, "C_$(id)_dc"] = @. α * (Pmax - df[!, "P_$id"])^2
        Cdc .+= df[!, "C_$(id)_dc"]
    end
    df.Cdc = Cdc
    df.C = df.Cdc .+ df.Cele
    return sum(df.C)
end

max_stars_repo_name: ShuhuaGao/HybridHEMS

max_issues_repo_name: ShuhuaGao/HybridHEMS

max_forks_repo_name: ShuhuaGao/HybridHEMS


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export common_prefix_distance, generic_levenshtein, hamming_distance, levenshtein_distance, lcs_distance

"""
    common_prefix(a, b)

common_prefix computes the length of the common prefix among
two strings represented as arrays
"""
function common_prefix(a, b)::Int
    len_a::Int = length(a)
    len_b::Int = length(b)
    i::Int = 1
    min_len::Int = min(len_a, len_b)
    @inbounds while i <= min_len && a[i] == b[i]
    	i += 1
    end

    i - 1
end

function common_prefix_distance(a, b)::Float64
    p = min(length(a), length(b))
	1.0 - common_prefix(a, b) / p
end

"""
    generic_levenshtein(a, b, icost::Int, dcost::Int, rcost::Int)::Int

Computes the edit distance between two strings, this is a low level function. Please use dist_lev
"""
function generic_levenshtein(a, b, icost::Int, dcost::Int, rcost::Int)::Int
    if length(a) < length(b)
        a, b = b, a
    end

    alen::Int = length(a)
    blen::Int = length(b)

    alen == 0 && return blen
    blen == 0 && return alen

    C::Vector{Int} = Array(0:blen)

    prevA::Int = 0
    @inbounds for i in 1:alen
        prevA = i
        prevC::Int = C[1]
        j::Int = 1

        while j <= blen
            cost::Int = rcost
            if a[i] == b[j]
               cost = 0
            end
            C[j] = prevA
            j += 1
            prevA = min(C[j]+dcost, prevA+icost, prevC+cost)
            prevC = C[j]
        end

        C[j] = prevA
    end

    prevA
end

"""
    levenshtein_distance(a, b)::Float64

computes edit distance over two strings
"""
function levenshtein_distance(a, b)::Float64
    generic_levenshtein(a, b, 1, 1, 1)
end

"""
    lcs_distance(a, b)::Float64

computes longest common subsequence distance
"""
function lcs_distance(a, b)::Float64
    generic_levenshtein(a, b, 1, 1, 2)
end

# function kerrormatch(a::T1, b::T2, errors::Int)::Bool where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return alen == blen
#     blen == 0 && return true

#     C::Vector{Int} = Vector{Int}(0:blen)

#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
# 	    end

#         C[j] = prevA
#         if prevA <= errors
#             return true
#         end
#     end

#     return false
# end

# function best_match_levenshtein(a::T1, b::T2)::Int where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return blen
#     blen == 0 && return alen

#     C::Vector{Int} = 1:blen |> collect

#     mindist = alen
#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
#         end

#         C[j] = prevA
#         if prevA < mindist
#             mindist = prevA
#         end
#     end

#     return mindist
# end

"""
     hamming_distance(a, b)::Float64

Computes the hamming distance between two sequences of the same length
"""

function hamming_distance(a, b)::Float64
    d::Int = 0

    @inbounds for i = 1:length(a)
        #if a[i] != b[i]
        #    d += 1
        #end
        d += Int(a[i] != b[i])
    end

    d
end


text:
module TestArithmetic

using Test
using DiffCheckNumbers

constant = DiffCheck{Float64, ()}(2.0, NamedTuple{(), NTuple{0, Float64}}(( )))
x = DiffCheck{Float64, (:x, )}(1.0, NamedTuple{(:x, ), NTuple{1, Float64}}((1.0, )))
xy = DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 1.0)))
z = DiffCheck{Float64, (:z, )}(1.0, NamedTuple{(:z, ), NTuple{1, Float64}}((1.0, )))
complex_z = DiffCheck(Complex{Float64}(1.0), :z)
complex_zw = DiffCheck(1.0+im*0, [:z, :w])

### Sum
@test x+constant === DiffCheck{Float64, (:x, )}(3.0, NamedTuple{(:x, ), NTuple{1, Float64}}((1.0, )))
@test x+2.0 === x+constant
@test x+2 === x+constant
@test x+Complex{Float64}(2.0) === x+constant

@test xy+constant === DiffCheck{Float64, (:y, :x)}(3.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 1.0)))
@test xy+2.0 === xy+constant
@test xy+2 === xy+constant
@test xy+Complex{Float64}(2.0) === xy+constant

@test xy+x === DiffCheck{Float64, (:y, :x)}(2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 2.0)))
@test (xy+x)+2 === DiffCheck{Float64, (:y, :x)}(4.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 2.0)))
@test xy+(x+2) === DiffCheck{Float64, (:y, :x)}(4.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 2.0)))

@test complex_z+constant === DiffCheck{Complex{Float64}, (:z, )}(3.0, NamedTuple{(:z, ), NTuple{1, Complex{Float64}}}((1.0, )))
@test complex_z+2.0 === complex_z+constant
@test complex_z+2 === complex_z+constant
@test complex_z+Complex{Float64}(2.0) === complex_z+constant

@test complex_zw+constant === DiffCheck{Complex{Float64}, (:w,:z)}(3.0, NamedTuple{(:w, :z), NTuple{2, Complex{Float64}}}((1.0, 1.0)))
@test complex_zw+2.0 === complex_zw+constant
@test complex_zw+2 === complex_zw+constant
@test complex_zw+Complex{Float64}(2.0) === complex_zw+constant
@test complex_zw+x === DiffCheck{Complex{Float64}, (:w, :z, :x)}(2.0, NamedTuple{(:w, :z, :x), NTuple{3, Complex{Float64}}}((1.0, 1.0, 1.0+0*im)))

### Difference
### Sum
@test x-constant === DiffCheck{Float64, (:x, )}(-1.0, NamedTuple{(:x, ), NTuple{1, Float64}}((1.0, )))
@test x-2.0 === x-constant
@test x-2 === x-constant
@test x-Complex{Float64}(2.0) === x-constant

@test xy-constant === DiffCheck{Float64, (:y, :x)}(-1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 1.0)))
@test xy-2.0 === xy-constant
@test xy-2 === xy-constant
@test xy-Complex{Float64}(2.0) === xy-constant

@test xy-x === DiffCheck{Float64, (:y, :x)}(0.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 0.0)))
@test (xy-x)-2 === DiffCheck{Float64, (:y, :x)}(-2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 0.0)))
@test xy-(x-2) === (xy-x)+2

@test complex_z-constant === DiffCheck{Complex{Float64}, (:z, )}(-1.0, NamedTuple{(:z, ), NTuple{1, Complex{Float64}}}((1.0, )))
@test complex_z-2.0 === complex_z-constant
@test complex_z-2 === complex_z-constant
@test complex_z-Complex{Float64}(2.0) === complex_z-constant

@test complex_zw-constant === DiffCheck{Complex{Float64}, (:w,:z)}(-1.0, NamedTuple{(:w, :z), NTuple{2, Complex{Float64}}}((1.0, 1.0)))
@test complex_zw-2.0 === complex_zw-constant
@test complex_zw-2 === complex_zw-constant
@test complex_zw-Complex{Float64}(2.0) === complex_zw-constant
@test complex_zw-x === DiffCheck{Complex{Float64}, (:w, :z, :x)}(0.0, NamedTuple{(:w, :z, :x), NTuple{3, Complex{Float64}}}((1.0, 1.0, -1.0+0*im)))
@test x-complex_zw === DiffCheck{Complex{Float64}, (:w, :z, :x)}(0.0-0.0*im, NamedTuple{(:w, :z, :x), NTuple{3, Complex{Float64}}}((-1.0-0*im, -1.0-0*im, +1.0+0*im)))

### Multiplication
@test x*constant === DiffCheck{Float64, (:x, )}(2.0, NamedTuple{(:x, ), NTuple{1, Float64}}((2.0, )))
@test x*2.0 === x*constant
@test x*2 === x*constant
@test x*Complex{Float64}(2.0) === x*constant

@test xy*constant === DiffCheck{Float64, (:y, :x)}(2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((2.0, 2.0)))
@test xy*2.0 === xy*constant
@test xy*2 === xy*constant
@test xy*Complex{Float64}(2.0) === xy*constant

@test xy*x === DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 2.0)))
@test (xy*x)*2 === DiffCheck{Float64, (:y, :x)}(2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((2.0, 4.0)))
@test xy*(x*2) === DiffCheck{Float64, (:y, :x)}(2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((2.0, 4.0)))

@test x*pi === DiffCheck{Float64, (:x, )}(pi, NamedTuple{(:x, ), NTuple{1, Float64}}((pi, )))

### Division

### Multiplication
@test x/constant === DiffCheck{Float64, (:x, )}(0.5, NamedTuple{(:x, ), NTuple{1, Float64}}((0.5, )))
@test x/2.0 === x/constant
@test x/2 === x/constant
@test x/Complex{Float64}(2.0) === x/constant

@test constant/x === DiffCheck{Float64, (:x, )}(2.0, NamedTuple{(:x, ), NTuple{1, Float64}}((-2.0, )))
@test 2.0/x === constant/x
@test 2/x === constant/x
@test Complex{Float64}(2.0)/x === constant/x

@test xy/constant === DiffCheck{Float64, (:y, :x)}(0.5, NamedTuple{(:y, :x), NTuple{2, Float64}}((0.5, 0.5)))
@test xy/2.0 === xy/constant
@test xy/2 === xy/constant
@test xy/Complex{Float64}(2.0) === xy/constant

@test constant/xy === DiffCheck{Float64, (:y, :x)}(2.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((-2.0, -2.0)))
@test 2/xy === constant/xy

@test xy/x === DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 0.0)))
@test (complex_z*x)/x === DiffCheck{Complex{Float64}, (:z, :x)}(1.0, NamedTuple{(:z, :x), NTuple{2, Complex{Float64}}}((1.0, 0.0)))
@test (2+im)*complex_z == DiffCheck{Complex{Float64}, (:z,)}(2+im, NamedTuple{(:z, ), NTuple{1, Complex{Float64}}}((2+im, )))

end


text:
# High Level Interface

function sim(;
    ρ=1.0,
    σ_PSF=.13,
    minphotons=50,
    ndatasets=10,
    nframes=1000,
    framerate=50.0,
    pattern::Pattern=Nmer2D(),
    molecule::Molecule=GenericFluor(;q=[0 50; 1e-2 0]),
    camera::Camera=IdealCamera()
)

# Simulation sequence
smld_true=SMLMSim.uniform2D(ρ,pattern,camera.xpixels*camera.pixelsize,camera.ypixels*camera.pixelsize)
smld_model=SMLMSim.kineticmodel(smld_true,molecule,nframes,framerate;ndatasets,minphotons)
smld_noisy=SMLMSim.noise(smld_model,σ_PSF)

return smld_true, smld_model, smld_noisy
end


text:
"""
    subsumes_string(u::String, v::String[, u_indexing])

Check whether stringified variable name `v` describes a sub-range of stringified variable `u`.

This is a very restricted version `subumes(u::VarName, v::VarName)` only really supporting:
- Scalar: `x` subsumes `x[1, 2]`, `x[1, 2]` subsumes `x[1, 2][3]`, etc.

## Note
- To get same matching capabilities as `AbstractPPL.subumes(u::VarName, v::VarName)`
  for strings, one can always do `eval(varname(Meta.parse(u))` to get `VarName` of `u`,
  and similarly to `v`. But this is slow.
"""
function subsumes_string(u::String, v::String, u_indexing=u * "[")
    return u == v || startswith(v, u_indexing)
end

"""
    inargnames(varname::VarName, model::Model)

Statically check whether the variable of name `varname` is an argument of the `model`.

Possibly existing indices of `varname` are neglected.
"""
@generated function inargnames(::VarName{s}, ::Model{_F,argnames}) where {s,argnames,_F}
    return s in argnames
end

"""
    inmissings(varname::VarName, model::Model)

Statically check whether the variable of name `varname` is a statically declared unobserved variable
of the `model`.

Possibly existing indices of `varname` are neglected.
"""
@generated function inmissings(
    ::VarName{s}, ::Model{_F,_a,_T,missings}
) where {s,missings,_F,_a,_T}
    return s in missings
end

# HACK: Type-piracy. Is this really the way to go?
AbstractPPL.getsym(::AbstractVector{<:VarName{sym}}) where {sym} = sym

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(TreeLas, Symbol("##gaplas#19")) && precompile(Tuple{getfield(TreeLas.MGT, Symbol("##gaplas#19")), Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:max_iter, :learn), Tuple{Int64, Float64}}}, typeof(TreeLas.MGT.gaplas), Array{Float64, 2}, Array{Tuple{Int64, Int64}, 1}, Array{Float64, 1}})
    isdefined(TreeLas, Symbol("#gaplas##kw")) && precompile(Tuple{getfield(TreeLas.MGT, Symbol("#gaplas##kw")), NamedTuple{(:max_iter, :learn), Tuple{Int64, Float64}}, typeof(TreeLas.MGT.gaplas), Array{Float64, 2}, Array{Tuple{Int64, Int64}, 1}, Array{Float64, 1}})
    precompile(Tuple{typeof(TreeLas.Dual.dual!), Array{Float64, 1}, Array{Int64, 1}, Array{Int64, 1}, Float64})
    precompile(Tuple{typeof(TreeLas.MGT.duality_check), Array{Float64, 1}, Array{Float64, 1}})
    precompile(Tuple{typeof(TreeLas.MGT.update_tree!), Array{Float64, 1}, Array{Float64, 1}, Array{Int64, 1}, GraphIdx.EdgeGraph, Array{Int64, 1}})
    precompile(Tuple{typeof(TreeLas.Pwl.QueueUnion.merge), Array{TreeLas.Pwl.Event, 1}, Array{TreeLas.Pwl.QueueUnion.Range, 1}, Int64, Int64})
    precompile(Tuple{typeof(TreeLas.Pwl.clip), Array{TreeLas.Pwl.Event, 1}, GraphIdx.Utils.MutRef{TreeLas.Pwl.QueueUnion.Range}, Float64, Float64, Base.Val{false}, Base.Val{true}})
    precompile(Tuple{typeof(TreeLas.Pwl.clip), Array{TreeLas.Pwl.Event, 1}, GraphIdx.Utils.MutRef{TreeLas.Pwl.QueueUnion.Range}, Float64, Float64, Base.Val{true}, Base.Val{true}})
    precompile(Tuple{typeof(TreeLas.TreeDP.tree_dp!), Array{Float64, 2}, Array{Float64, 2}, GraphIdx.Tree.RootedTree, GraphIdx.Vec{Float64}, GraphIdx.Ones{Float64}, TreeLas.TreeDP.TreeDPMem{Float64, Int64}})
    precompile(Tuple{typeof(TreeLas.Utils.primal_objective), Array{Float64, 2}, Array{Float64, 2}, GraphIdx.WeightedGraph, GraphIdx.Const{Float64}})
    precompile(Tuple{typeof(TreeLas.Utils.primal_objective), Array{Float64, 2}, Array{Float64, 2}, GraphIdx.WeightedGraph})
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# module structure:
Base.:*(a::Number, X::AlgebraElement) =
    mul!(similar(X, promote_type(eltype(X), typeof(a))), X, a)
Base.:*(X::AlgebraElement, a::Number) = a * X
Base.:(/)(X::AlgebraElement, a::Number) = inv(a) * X

# TODO: handle this through mul!?
Base.:(//)(X::AlgebraElement, a::Number) = AlgebraElement(coeffs(X) .// a, parent(X))

# ring structure:
Base.:-(X::AlgebraElement) = neg!(similar(X), X)

function _preallocate_output(X::AlgebraElement, Y::AlgebraElement, op)
    T = Base._return_type(op, Tuple{eltype(X),eltype(Y)})
    if coeffs(Y) isa DenseArray
        return similar(Y, T)
    end
    return similar(X, T)
end

Base.:+(X::AlgebraElement, Y::AlgebraElement) = add!(_preallocate_output(X, Y, +), X, Y)
Base.:-(X::AlgebraElement, Y::AlgebraElement) = sub!(_preallocate_output(X, Y, -), X, Y)
Base.:*(X::AlgebraElement, Y::AlgebraElement) = mul!(_preallocate_output(X, Y, *), X, Y)

Base.:^(a::AlgebraElement, p::Integer) = Base.power_by_squaring(a, p)

# mutable API; TODO: replace with MutableArithmetic

function zero!(a::AlgebraElement)
    a.coeffs .= zero(eltype(coeffs(a)))
    return a
end

function neg!(res::AlgebraElement, X::AlgebraElement)
    @assert parent(res) === parent(X)
    res.coeffs .= -coeffs(X)
    return res
end

function add!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
    @assert parent(res) === parent(X) === parent(Y)
    # res = (res === X || res === Y) ? similar(res) : res
    res.coeffs .= coeffs(X) .+ coeffs(Y)
    return res
end

function sub!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
    @assert parent(res) === parent(X) === parent(Y)
    # res = (res === X || res === Y) ? similar(res) : res
    res.coeffs .= coeffs(X) .- coeffs(Y)
    return res
end

function mul!(res::AlgebraElement, X::AlgebraElement, a::Number)
    @assert parent(res) === parent(X)
    # res = (res === X) ? similar(res) : res
    res.coeffs .= a .* coeffs(X)
    return res
end

function mul!(
    res::AbstractVector,
    X::AbstractVector,
    Y::AbstractVector,
    ms::MultiplicativeStructure,
)
    res = (res === X || res === Y) ? zero(res) : (res .= zero(eltype(res)))
    return fmac!(res, X, Y, ms)
end

function mul!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
    @assert parent(res) === parent(X) === parent(Y)
    res = (res === X || res === Y) ? zero(res) : zero!(res)
    return fmac!(res, X, Y)
end

function fmac!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement)
    fmac!(coeffs(res), coeffs(X), coeffs(Y), parent(res).mstructure)
    return res
end

_nzpairs(v::AbstractVector) = pairs(v)
_nzpairs(v::AbstractSparseVector) =
    zip(SparseArrays.nonzeroinds(v), SparseArrays.nonzeros(v))

function fmac!(
    res::AbstractVector,
    X::AbstractVector,
    Y::AbstractVector,
    mstr::MultiplicativeStructure,
)
    @assert res !== X
    @assert res !== Y
    for (j, y) in _nzpairs(Y)
        for (i, x) in _nzpairs(X)
            res[mstr[i, j]] += x * y
        end
    end
    return res
end


text:
using DSGE, SparseArrays, MAT, FileIO

# Script for doing all that MATLAB does
#m = TwoAssetHANK()
#steadystate!(m)
reduceDistribution = true #get_setting(m, :reduceDistribution)
reduceDist_hor     = 100 #get_setting(m, :reduceDist_hor)
reduceV            = 1 #get_setting(m, :reduceV)

nVars = 127504 #get_setting(m, :nVars)
n_v = 60000 #get_setting(m, :n_v)
n_g = 67499 #get_setting(m, :n_g)
n_p = 4 #get_setting(m, :n_p)
n_Z = 1 #get_setting(m, :n_Z)

#Γ0, Γ1, Ψ, Π C = eqcond(m)
test_out = load("data/eqcond_output_matlab.jld2")
Γ0 = test_out["g0"]
Γ1 = test_out["g1"]
Ψ = test_out["psi"]
Π = test_out["pi"]
C = spzeros(Float64, nVars)

# Distribution Reduction
base_mat = my_speye(nVars)
inv_base = my_speye(nVars)
state_red = my_speye(nVars-n_p)
inv_state_red = my_speye(nVars-n_p)

println("=> Reducing Value Function Distribution")
@show reduceDistribution
#error()
@time if reduceDistribution

   state_red, inv_state_red, n_g = stateSpaceReduction(Γ0, Γ1, n_v, n_g, n_p, n_Z, reduceDist_hor)

else
    Γ0, Γ1, C, Π, Ψ, base_mat, inv_base = solve_static_conditions(Γ0, Γ1, Π, Ψ, C)
    # Confirming these do the same thing
    base_mat2, inv_base2, g1, c, pi, psi  = cleanG0Sparse(Γ0, Γ1, C, Π, Ψ)
    @show base_mat2 ≈ base_mat, inv_base ≈ inv_base2, g1 ≈ Γ1, Π ≈ pi, Ψ ≈ psi
end

# Value Function Reduction
from_spline = my_speye(n_g + n_v + n_Z)
to_spline   = my_speye(n_g + n_v + n_Z)
n_splined   = n_v
n_rest      = n_g + n_Z

if reduceV
    #to_spline, from_spline, n_splined, n_rest = reduce
end


text:
# Log Magnitude Approximation (LMA) digital filter

type LMABaseFilter <: Filter
    order::Int
    delay::Vector{Float64}

    function LMABaseFilter(order::Int)
        new(order, zeros(order+1))
    end
end

function filt!(f::LMABaseFilter, x::Float64, coef::Vector{Float64},
               m1::Int, m2::Int)
    @assert length(coef) == f.order+1
    @assert m1 >= 1
    @assert m2 < length(coef)

    d = delay(f)

    for i=length(coef)-1:-1:2
        @inbounds d[i] = d[i-1]
    end
    d[1] = x

    result = 0.0
    for i=m1+1:m2+1
        @inbounds result += coef[i] * d[i-1]
    end

    result
end

# LMACascadeFilter represents a cascade filter which contains LMA base
# filters.
type LMACascadeFilter <: Filter
    filters::Vector{LMABaseFilter}
    padecoef::Vector{Float64}
    delay::Vector{Float64}

    function LMACascadeFilter(order::Int, pade::Int)
        padecoef = Array(Float64, pade+1)

        if pade == 4
            padecoef = padecoef4th
        elseif pade == 5
            padecoef = padecoef5th
        else
            error("LMADF: Order of pade approximation 4 or 5 is only supported.")
        end

        filters = Array(LMABaseFilter, pade+1)
        for i=1:length(filters)
            filters[i] = LMABaseFilter(order)
        end

        new(filters, padecoef, zeros(pade+1))
    end
end

padecoef(f::LMACascadeFilter) = f.padecoef

function filt!(f::LMACascadeFilter, x::Float64, coef::Vector{Float64},
                 m1::Int, m2::Int)
    d = delay(f)
    pade = padecoef(f)
    result, feedback = 0.0, 0.0

    for i=length(pade):-1:2
        @inbounds d[i] = filt!(f.filters[i], d[i-1], coef, m1, m2)
        @inbounds val = d[i] * pade[i]
        if iseven(i)
            feedback += val
        else
            feedback -= val
        end
        result += val
    end

    d[1] = feedback + x
    result += d[1]
    result
end

type LMADF <: MelGeneralizedCepstrumSynthesisFilter
    filters::Vector{LMACascadeFilter}

    function LMADF(order::Int; pade::Int=5)
        filters = Array(LMACascadeFilter, 2)
        filters[1] = LMACascadeFilter(order, pade)
        filters[2] = LMACascadeFilter(order, pade)
        new(filters)
    end
end

first(f::LMADF) = f.filters[1]
last(f::LMADF) = f.filters[2]
allpass_alpha(f::LMADF) = 0.0
glog_gamma(f::LMADF) = 0.0

function filt!(f::LMADF, x::Float64, coef::Vector{Float64})
    m = length(coef)-1
    y = filt!(first(f), x, coef, 1, 1)
    filt!(last(f), y, coef, 2, m)
end

function to_filtcoef(f::LMADF, c::SpectralParamState{LinearCepstrum})
    copy(rawdata(c))
end


text:
normalize_total(p::Profile, library_size::Real=1e6)

"""
Logarithmize the data matrix.
"""
function log1p(p::Profile; base::Real=ℯ)
    log1p.(p.data) ./ log(base)
end


text:
###############################
# Interface for MLDataPattern #
###############################
# Number of objects for data cells
"""
	nobs(data, [obsdim])

Get the number of observations in a `DataCell`.Check out documentation from `MLDataPatterns` for more details.
"""
nobs(::Void, args...) = 0
nobs(::Void, ::LearnBase.ObsDim.Undefined) = 0
nobs(x::T where T<:CellDataVectors, args...) = size(getx!(x),1)
nobs(x::T where T<:CellDataVectors, ::LearnBase.ObsDim.Undefined) = size(getx!(x),1)
nobs(x::T where T<:CellDataMatrices, args...) = size(getx!(x),2)
nobs(x::T where T<:CellDataMatrices, ::LearnBase.ObsDim.Undefined) = size(getx!(x),2)

"""
	getobs(data, [idx], [obsdim])

Get observation from `DataCell`. Check out documentation from `MLDataPatterns` for more details.
"""
getobs(x::T where T<:CellDataVectors, idx) = x[idx]
getobs(x::T where T<:CellDataMatrices, idx) = x[:,idx]
#getobs!{T<:CellData}(buf, x::T, args...) = copy!(buf, datacell(getobs(strip(x), args...)));

"""
	datasubset(data, [idx], [obsdim])

Returns a lazy subset from a `DataCell`. Check out documentation from `MLDataPatterns` for more details.
"""
#datasubset{T<:CellDataVectors}(x::T, idx) = Cell(view(getx!(x),idx))
#datasubset(x::T where T<:CellData, idx; obsdim=LearnBase.ObsDim.Constant{2}()) = DataCell(datasubset(strip(x), idx; obsdim=obsdim))
datasubset(x::T where T<:CellData, idx::Int) = DataCell(datasubset(strip(x),idx:idx))
datasubset(x::T where T<:CellData, idx::Int, ::LearnBase.ObsDim.Undefined) = DataCell(datasubset(strip(x),idx:idx))
datasubset(x::T where T<:CellData, idx, ::LearnBase.ObsDim.Undefined) = DataCell(datasubset(strip(x), idx))

"""
	varsubset(data,[idx])

Similar to `datasubset` however it returns variable subsets
"""
varsubset(x::T where T<:CellDataVec, idx) = DataCell(view(getx!(x),:,idx))				# for unlabeled vector data, ignore index, labels
varsubset(x::T where T<:CellDataVecVec, idx) = DataCell(view(getx!(x),:,idx), view(gety!(x),:))		# for labeled vector data, ignore index
varsubset(x::T where T<:CellDataVecMat, idx) = DataCell(view(getx!(x),:,idx), view(gety!(x),:,:))	# for labeled vector data, ignore idex
varsubset(x::T where T<:CellDataMat, idx) = DataCell(view(getx!(x),idx,:)) 				# For unlabeled matrix data, ignore labels
varsubset(x::T where T<:CellDataMatVec, idx) = DataCell(view(getx!(x),idx,:), view(gety!(x),:))
varsubset(x::T where T<:CellDataMatMat, idx) = DataCell(view(getx!(x),idx,:), view(gety!(x),:,:))
varsubset(x::T where T<:CellData, idx, ::LearnBase.ObsDim.Undefined) = varsubset(x, idx)
varsubset(x::T where T<:AbstractArray, idx) = _variable_(x,idx)

# Data access: idx selects variable
_variable_(x::T where T<:AbstractVector, idx) = view(x, :, idx)
_variable_(x::T where T<:AbstractMatrix, idx) = view(x, idx, :)
_variable_(x::T where T<:CellData, idx) = _variable_(getx!(x), idx)

# Data access when specifying empty static arrays (always return nothing)
_variable_(x::T where T<:AbstractArray, ::StaticArrays.SVector{0}) = nothing
_variable_(::AbstractArray{T,2} where T, ::StaticArrays.SArray{Tuple{0},S,1,0} where S) = nothing

# gettargets and gettarget low-level functions for MLDataPattern integration
LearnBase.gettargets(x::T where T<:CellDataU, idx, obsdim::LearnBase.ObsDimension=LearnBase.ObsDim.Undefined()) = nothing
LearnBase.gettargets(x::T where T<:CellDataL, idx, obsdim::LearnBase.ObsDimension=LearnBase.ObsDim.Undefined()) = gety!(x)[idx]
LearnBase.gettargets(x::T where T<:CellDataLL, idx::Number, obsdim::LearnBase.ObsDimension=LearnBase.ObsDim.Undefined()) = gety!(x)[:,idx:idx]
LearnBase.gettargets(x::T where T<:CellDataLL, idx, obsdim::LearnBase.ObsDimension=LearnBase.ObsDim.Undefined()) = gety!(x)[:,idx]
LearnBase.gettarget(x::T where T<:CellDataU) = nothing
LearnBase.gettarget(f, x::T where T<:CellDataU) = nothing
LearnBase.gettarget(x::T where T<:CellDataL) = gety!(x)[1]
LearnBase.gettarget(f, x::T where T<:CellDataL) = f(gety!(x)[1])
LearnBase.gettarget(x::T where T<:CellDataLL) = gety!(x)
LearnBase.gettarget(f, x::T where T<:CellDataLL) = f(gety!(x))

# Optimization for the targets function
targets(x::T where T<:CellDataU, args...) = nothing
targets(f::Function, x::T where T<:CellDataU, args...) = nothing
targets(f::Function, x::T where T<:CellData, args...) = targets(f, gety!(x), args...)
targets(f::Function, obsdim::LearnBase.ObsDimension=LearnBase.ObsDim.Constant{2}()) =
	FunctionCell((x)->targets(f, x, obsdim),(), "Target generator ($(string(f)))")


text:
using MemoryBuffers
using Test

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


text:
push!(LOAD_PATH,"../src/")
push!(LOAD_PATH,"../docs/")

using Documenter
using LTVSourceReceiverModel

makedocs(;
    modules=[LTVSourceReceiverModel],
    authors="Mamta Dalal, Steven Sandoval, Hasan Al-Shammari",
    sitename="LTVSourceReceiverModel.jl",
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "LTI Omnidirectional" =>
                    ["Basic Forward Model"     => "omniDirectionalLTI/source.md",],
        "Transmitters" => "transmitters/transmitter.md",
        "Receivers" => "receivers/receiver.md",
        "Types/Methods/Functions" => "typesMethodsFunctions.md",
        "Cite" => "cite.md",
    ],
    repo="https://github.com/NMSU-ISA/LTVSourceReceiverModel.jl/blob/{commit}{path}#{line}",
    assets=String[],
)

deploydocs(;
    repo="github.com/NMSU-ISA/LTVSourceReceiverModel",
    devurl = "dev"
)


text:
using StatsBase:sample, Weights

"""
    RandomPolicy(prob::Array{Float64, 2})
    RandomPolicy(prob::Vector{Float64})

The probability of each action is predefined by `prob`.
If `prob` is a vector, then all states share the same `prob`.
"""
struct RandomPolicy{n} <: AbstractPolicy
    prob::Array{Float64, n}
    RandomPolicy(prob::Array{Float64, 2}) = new{2}(prob)
    RandomPolicy(prob::Vector{Float64}) = new{1}(prob)
end

(p::RandomPolicy{2})(s::Int) = sample(Weights(p.prob[s, :], 1.0))
(p::RandomPolicy{2})(s::Int, ::Val{:dist}) = p.prob[s, :]
(p::RandomPolicy{2})(s::Int, a::Int) = p.prob[s, a]

(p::RandomPolicy{1})(s) = sample(Weights(p.prob, 1.0))
(p::RandomPolicy{1})(s, ::Val{:dist}) = p.prob
(p::RandomPolicy{1})(s, a::Int) = p.prob[a]


text:
@testset "mlp" begin
     atype = KnetLayers.arrtype
     x =atype(randn(10,2))
     m = MLP(10,5,2;winit=randn,binit=zeros)
     y = m(x)
     @test true
end


text:
using Plots

w = 3
h = 3
psi = 0
XQ = [-w/2 w/2 w/2 -w/2 -w/2]
YQ = [h/2 h/2 -h/2 -h/2 h/2]
P = [XQ;YQ]
ct = cos(psi)
st = sin(psi)
R = [ct -st;st ct]
P2 = R*P

plot()
function pltCurb(P2,w,h)
    scatter!((P2[1,:]+w,P2[2,:]+h),marker=(:circle,:black,0.,0.),fill=(0,1,:black),leg=true,grid=true,label="")
end

function intersectionPlot()
    pp = plot(0,leg=:false)

    # parameters TODO add to YAML file
    lw = 3
    w = 3
    h = 3
    psi = 0

    XQ = [-w/2 w/2 w/2 -w/2 -w/2]
    YQ = [h/2 h/2 -h/2 -h/2 h/2]
    P = [XQ;YQ]
    ct = cos(psi)
    st = sin(psi)
    R = [ct -st;st ct]
    P2 = R*P

    pltCurb(P2,w+lw/2,h+lw/2)
    pltCurb(P2,w+lw/2,-h-lw/2)
    pltCurb(P2,-w-lw/2,-h-lw/2)
    pltCurb(P2,-w-lw/2,h+lw/2)

    return pp
end

pp = intersectionPlot()
lim_size = 5
xlims!(-lim_size,lim_size)
ylims!(-lim_size,lim_size)


text:
using StructuresKit

#Test moment and shear calculation for a simply-supported beam.

#Ix Iy Ixy J Cw
sectionProperties = [(100.0, 100.0, 0.0, 10.0, 10.0)]

#E  ν
materialProperties = [(29500, 0.30)]

#ax ay
loadLocation = [(0.0, 0.0)]

#kx kϕ
springStiffness = [(0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness
memberDefinitions = [(30*12,3,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
endBoundaryConditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 30*12]

#load  qx   qy
uniformLoad = (0.0, 5.0/1000)

z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, springStiffness, endBoundaryConditions, supports, uniformLoad)

Mxx = InternalForces.moment(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)

Vyy = InternalForces.shear(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)

maxMxx = maximum(Mxx)
maxMxx_AISC = uniformLoad[2]*z[end]^2/8

maxVyy = maximum(Vyy)
maxVyy_AISC = uniformLoad[2]*z[end]/2

MxxError= abs((maxMxx - maxMxx_AISC))/ maxMxx_AISC
VyyError= abs((maxVyy - maxVyy_AISC))/ maxVyy_AISC

#accept 1% error from numerical solution
@test MxxError <= 0.01
@test VyyError <= 0.01


text:
module TailRec
export @tailrec

macro tailrec(func)
    fargs=map(e->isa(e,Expr) ? e.args[1] : e,func.args[1].args)
    fbody=func.args[2]
    fbody=rewrite(fbody,fargs)
    func.args[2]=Expr(:block,:(@label retry),fbody)
    esc(func)
end

function rewrite(expr,args,callflag=false)
    if !isa(expr,Expr)
        expr
    elseif expr.head == :call && expr.args[1] == args[1]
        if callflag
            @warn("Not tail recursive.")
            expr
        else
            newargs=Expr(:tuple)
            newargs.args=args[2:end]
            oldargs=Expr(:tuple)
            oldargs.args=expr.args[2:end]
            Expr(:block, Expr(:(=),newargs,oldargs), :(@goto retry) )
        end
    elseif expr.head == :block
        expr.args[end]=rewrite(expr.args[end],args,expr.head==:call ? true : callflag)
        expr
    else
        expr.args = map(a->rewrite(a,args,expr.head==:call ? true : callflag), expr.args)
        expr
    end
end
end


text:
"""
$(TYPEDEF)

Represents a simple control system with linear dynamics.

# Parameters:

- `ΔT`: the sampling time of the system (se)
- `nx`: the number of states
- `nu`: the number of inputs

# Fields

$(TYPEDFIELDS)
"""
struct LinearSystem{ΔT, nx, nu, TA<:SMatrix{nx, nx}, TB<:SMatrix{nx, nu}} <: ControlSystem{ΔT, nx, nu}
    "The state transition matrix"
    A::TA
    "The control input matrix"
    B::TB
end

LinearSystem{ΔT}(A::TA, B::TB) where {ΔT, nx, nu, TA<:SMatrix{nx, nx}, TB<:SMatrix{nx, nu}} = LinearSystem{ΔT, nx, nu, TA, TB}(A, B)
LinearizationStyle(::LinearSystem) = TrivialLinearization()

dx(ls::LinearSystem, x::SVector, u::SVector, t::AbstractFloat)  = begin @assert !issampled(ls); ls.A*x + ls.B*u end
next_x(ls::LinearSystem, x::SVector, u::SVector) = begin @assert issampled(ls); ls.A*x + ls.B*u end
linearize(ls::LinearSystem, x::SVector, u::SVector, t::AbstractFloat) = ls
function discretize(ls::LinearSystem, vt::Val{ΔT}) where {ΔT}
    @assert !issampled(ls) "Can't discretize a discrete system."
    @assert ΔT > 0 "Discrtization requires finite sampling time ΔT."
    # use the euler approximation to avoid expensive matrix exponential
    return LinearSystem{ΔT}(I+ΔT*ls.A, ΔT*ls.B)
end

# TODO: this should probably know the absolute time (like the `SystemTrajectory`)
struct LTVSystem{ΔT, nx, nu, TD<:AbstractVector{<:LinearSystem{ΔT, nx, nu}}} <: ControlSystem{ΔT, nx, nu}
    "The discrete time series of linear systems."
    dyn::TD
    LTVSystem(dyn::TD) where {ΔT, nx, nu, TD<:AbstractVector{<:LinearSystem{ΔT, nx, nu}}} = begin
        @assert ΔT > 0 "LTVSystem require finite discretization steps."
        new{ΔT, nx, nu, TD}(dyn)
    end
end
LinearizationStyle(::LTVSystem) = TrivialLinearization()
Base.eltype(::Type{<:LTVSystem{ΔT, nx, nu, TD}}) where {ΔT, nx, nu, TD} = eltype(TD)
Base.getindex(ds::LTVSystem, i) = getindex(ds.dyn, i)
Base.setindex!(ds::LTVSystem, v, i) = setindex!(ds.dyn, v, i)
next_x(cs::LTVSystem, xₖ::SVector, uₖ::SVector, k::Int) = next_x(cs.dyn[k], xₖ, uₖ)

# TODO: we could also get away without this type but it seems more clean (?)
struct LTISystem{ΔT,nx,nu,TL<:LinearSystem{ΔT,nx,nu},
                 TXY,TXI}<:ControlSystem{ΔT,nx,nu}
    dyn::TL
    xyids::TXY
    xids::TXI
end
LTISystem(dyn::LinearSystem, xyids) = LTISystem(dyn, xyids, nothing)

LinearizationStyle(::LTISystem) = TrivialLinearization()
xyindex(cs::LTISystem) = cs.xyids
xindex(cs::LTISystem) = cs.xids
Base.eltype(::Type{<:LTISystem{ΔT,nx,nu,TL}}) where {ΔT,nx,nu,TL} = TL
Base.getindex(cs::LTISystem, i) = cs.dyn
next_x(cs::LTISystem, xₖ::SVector, uₖ::SVector, ::Float64) = next_x(cs.dyn, xₖ, uₖ)


text:
struct TrainTest <: Split
    at::Float64

    TrainTest(at::Real = 0.75) = new(at)
end

function load(
    split::TrainTest,
    dataset::Dataset{N, F, T}
) where {N<:Name, F<:Format, T<:Task}

    train = load(dataset, :train; dopostprocess = false)
    if !hassubset(N, :test)
        train, test = data_split(dataset, train, samplesdim(F), split.at)
    else
        test = load(dataset, :test; dopostprocess = false)
    end

    train = postprocess(dataset, train)
    test = postprocess(dataset, test)

    return train, test
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider