text:
###### activate environment #####
using Pkg
Pkg.activate(pwd())

##### packages #####
using Dates
using DataFrames, CSV
using Distributions

##### defining structs #####

"""
Struct taking different values for summer and winter seasons
"""
mutable struct demand
    summer::AbstractVector
    winter::AbstractVector
end

"""
Struct for data on each appliance
"""
mutable struct appliance
    name::String
    wattage::Float64
    duration::demand
    start::demand
end

"""
Struct for appliance data on each household
"""
mutable struct household
    name::String
    appliances::AbstractVector
end

##### helping functions #####

"""
Draw appliance start from a uniform distribution depending on the length of the day segment
"""
function sample_start(i::Int64)

    # sampling for segments of the day from uniform distribution
    if 0 < i < length(segments_duration[:,1])
        s = rand(1:segments_duration[i,2]) - 1
        return s

    # no sampling if aplliance not in use
    elseif i == 7
        s = 0
        return s

    else
        error("out of range")
    end
end

"""
Draw appliance uptime duration from a truncated normal distribution
"""
function sample_duration(i::Int64)

    # sampling from truncated normal distribution with mean at max duration and truncation at mean
    if i > 0
        d = round(rand(TruncatedNormal(i,1,0,i)),digits=0)
        return d

    # no sampling if appliance not in use
    elseif i == 0
        d = 0
        return d
    else
        error("out of range")
    end
end

"""
Determine type of day (weekday, Friday, Saturday)
"""
function test_day(day::Date)

    # weekday
    if 0 < dayofweek(day) < 5 || dayofweek(day) == 7
        daytype = 1

    # Friday
    elseif dayofweek(day) == 5
        daytype = 2

    # Saturday
    elseif dayofweek(day) == 6
        daytype = 3

    else
        error("out of bounds")
    end
    return daytype
end

"""
Generate an hourly timecode between a starting day s and an ending day e
"""
function generate_timecode(s::Date,e::Date)
    # add hours 00:00 to 23:00
    ts_time_int = s+Time(0):Hour(1):e+Time(23)
    ts_timecode = collect(ts_time_int)
    return ts_timecode
end

##### main functions #####

"""
Generate starting and stopping hour for an appliance a for a specific day d
"""
function app_duration(a::appliance,d::Date)

    # determine start and stop in summer
    if seasons[month(d)] == 1
        start = d + Time(Dates.Hour(segments_start[a.start.summer[test_day(d)],2] + sample_start(a.start.summer[test_day(d)])))
        stop = start + Dates.Hour(sample_duration(a.duration.summer[test_day(d)]))
        return start, stop

    # determine start and stop in winter
    elseif seasons[month(d)] == 2
        start = d + Time(Dates.Hour(segments_start[a.start.summer[test_day(d)],2] + sample_start(a.start.winter[test_day(d)])))
        stop = start + Dates.Hour(sample_duration(a.duration.winter[test_day(d)]))
        return start, stop

    else
        error("out of bounds")
    end
end

"""
Generate load time series for an appliance a between a starting day s and an ending day e with hourly timecode tc
"""
function generate_app_load_ts(a::appliance,s::Date,e::Date,tc::Vector)

    # initialize dataframe with time code
    app_load_ts = DataFrame()
    app_load_ts.Timecode = tc
    app_load_ts[!, :load] .= 0.0

    # add hourly device wattage between start and stop
    for i in collect(s:Day(1):e)
        start_stop = app_duration(a,i)
        for i in 1:length(tc)
            if start_stop[1] <= tc[i] <= start_stop[2]
                app_load_ts.load[i] = app_load_ts.load[i] + a.wattage
            end
        end
    end

    return app_load_ts
end

"""
Generate load time series for a specific household hh based on a starting day s and an ending day e
"""
function generate_hh_load_ts(hh::household,s::Date,e::Date)

    # initialize dataframe with time code
    hh_load_ts = DataFrame()
    hh_load_ts.Timecode = generate_timecode(s,e)

    # simulate every device
    for i in 1:length(hh.appliances)
        if hh.appliances[i] == 1
            hh_load_ts.appl = generate_app_load_ts(appl[i],s,e,hh_load_ts.Timecode).load
            rename!(hh_load_ts,:appl => appl[i].name)
        elseif hh.appliances[i] > 1
            hh_load_ts.appl = generate_app_load_ts(appl[i],s,e,hh_load_ts.Timecode).load
            for n in 2:hh.appliances[i]
                hh_load_ts.appl = hh_load_ts.appl + (generate_app_load_ts(appl[i],s,e,hh_load_ts.Timecode).load * rand(d))
            end
            rename!(hh_load_ts,:appl => appl[i].name)
        end
    end

    # add one column of total demand (sum) to dataframe
    hh_load_ts.total_load = sum(eachcol(hh_load_ts[!,2:end]))
    return hh_load_ts
end

"""
Generate load time series for a swarm grid based on a start day s and an end day e
"""
function generate_sg_load_ts(s::Date,e::Date)

    # initialize dataframe with time code
    sg_load_ts = DataFrame()
    sg_load_ts.Timecode = generate_timecode(s,e)

    # simulate every household
    for i in 1:length(hh)
        @info "simulating household $(hh[i].name)"
        sg_load_ts.hh = generate_hh_load_ts(hh[i],start_ts,end_ts).total_load
        rename!(sg_load_ts,:hh => hh[i].name)
    end

    return sg_load_ts
end

##### loading data #####

# appliance data
appl_dat = CSV.File("_input/data_appliances.csv") |> DataFrame
appl = []

for i in 1:nrow(appl_dat)
    push!(appl,appliance(appl_dat.name[i],appl_dat.wattage[i]/1000,demand([appl_dat.h_sum_wd[i],appl_dat.h_sum_fri[i],appl_dat.h_sum_sat[i]],[appl_dat.h_win_wd[i],appl_dat.h_win_fri[i],appl_dat.h_win_sat[i]]),demand([appl_dat.s_sum_wd[i],appl_dat.s_sum_fri[i],appl_dat.s_sum_sat[i]],[appl_dat.s_win_wd[i],appl_dat.s_win_fri[i],appl_dat.s_win_sat[i]])))
end

# household data
hh_dat = CSV.File("_input/data_households.csv") |> DataFrame
hh = []

for i in 1:nrow(hh_dat)
    push!(hh,household(hh_dat[i,1],Array(hh_dat[i,2:end])))
end

# define segments of the day
"""
# 1 = night (00-06)
# 2 = morning (06-09)
# 3 = late morning (09-12)
# 4 = afternoon (12-15)
# 5 = late afternoon (15-18)
# 6 = evening (18-00)
# 7 = not in use
"""
segments_duration = [1 6;2 3;3 3;4 3;5 3;6 6; 7 0]
segments_start = [1 0;2 6;3 9;4 12;5 15;6 18; 7 0]

# define seasons of the months
"""
1 = summer
2 = winter
"""
seasons = [1,1,1,1,2,2,2,2,2,1,1,1]

# distribution and probability of use for cases where more than one of the same appliance exists in a household
p = 0.3
d = Binomial(1,p)

##### generate smart grid load profile #####

start_ts = Date("2019-06-01") # start day of time series
end_ts = Date("2019-06-07") # end day of time series

sg_load_ts = generate_sg_load_ts(start_ts,end_ts)

CSV.write("_output\\sg_load_ts-$(start_ts)_$(end_ts).csv", sg_load_ts)

##### plotting results #####

using Plots
using TimeSeries

sg_load_ta = TimeArray(sg_load_ts, timestamp = :Timecode) # convert dataframe to timearray
plot(sg_load_ta)
savefig("_output\\sg_load_ts-$(start_ts)_$(end_ts).png")

max_stars_repo_name: weibezahn/SGDemTool

max_issues_repo_name: weibezahn/SGDemTool

max_forks_repo_name: weibezahn/SGDemTool


text:
module IdrisForeignCollections
using CamJulia.IdrisList

FList = IdrisList

const HVECT_NIL = Symbol("Data.HVect.Nil")
const HVECT_CONS = Symbol("Data.HVect.::")

const VECT_NIL = Symbol("Data.Vect.Nil")
const VECT_CONS = Symbol("Data.Vect.::")

const LIST_NIL = Symbol("Prelude.List.Nil")
const LIST_CONS = Symbol("Prelude.List.::")

export from_flist, to_flist
export to_flist!
export map_hvect
export to_fvect, from_fvect
export to_fhvect, from_fhvect

function from_flist(jl_lst :: Vector{T}) where T
    foldr(jl_lst, init=LIST_NIL) do each, prev
         (LIST_CONS, each, prev)
    end
end

function to_flist!(c :: Vector{T}, flist::Symbol) where T
    nothing
end

function to_flist!(c :: Vector{T}, flist::Tuple{Symbol, T, A}) where {T, A}
    while flist !== LIST_NIL
        push!(c, flist[2])
        flist = flist[3]
    end
end

function to_flist(flist :: Symbol)
    []
end

function to_flist(flist::Tuple{Symbol, T, A}) where {T, A}
    c = T[]
    to_flist!(c, flist)
    c
end

function map_hvect(f, a::Symbol)
    a
end

function map_hvect(f, a :: Tuple{Symbol, T, A}) where {T, A}
    (a[1], f(a[2]), map_hvect(f, a[3]))
end

function from_fvect(n :: Int64, vec :: Vector{T}) where T
    prev = VECT_NIL
    for i = n:-1:1
        prev = (VECT_CONS, vec[i], prev)
    end
    prev
end

function to_fvect(n :: Int64, vec :: Symbol)
    []
end

function to_fvect(n :: Int64, vec :: Tuple{Symbol, T, A}) where {T, A}
    v = Vector{T}(UndefInitializer(), n)
    i = 1
    while vec[3] !== VECT_NIL
        v[i] = vec[2]
        vec = vec[3]
        i = i + 1
    end
    @assert vec[3] === VECT_NIL
    v[i] = vec[2]
    v
end

function from_fhvect(n :: Int64, vec :: Tuple) where T
    prev = HVECT_NIL
    for i = n:-1:1
        prev = (HVECT_CONS, vec[i], prev)
    end
    prev
end

function to_fhvect(n :: Int64, vec :: Symbol)
    ()
end

@generated function to_fhvect(n :: Int64, vec :: Tuple{Symbol, T, A}) where {T, A}
    node = vec
    vars = []
    while node <: Tuple
        ps = node.parameters
        push!(vars, gensym("elt" * string(ps[2])))
        node = ps[end]
    end

    tp = gensym("tp")
    expr = foldr(vars, init=Expr(:tuple, vars...)) do cur_var, prev
        :(let $cur_var = $tp[2], $tp=$tp[3]; $prev end)
    end
    :(let $tp = vec; $expr end)
end

end

ext: jl

max_stars_repo_name: thautwarm/idris-cam

max_issues_repo_name: thautwarm/idris-cam

max_forks_repo_name: thautwarm/idris-cam


text:
## test_spectral.jl : test spectral method for GRF generation

@testset "spectral decomposition   " begin

    ## 1d Mat\'ern ##
    cov = CovarianceFunction(1,Matern(0.3,1))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,Spectral(),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## 2d Exponential ##
    cov = CovarianceFunction(2,Exponential(1))
    pts1 = 0:0.05:1
    pts2 = 0:0.05:1
    grf = GaussianRandomField(cov,Spectral(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Exponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## 3d Gaussian ##
    cov = CovarianceFunction(3,Matern(1,2.5))
    pts1 = 0:0.1:1
    pts2 = 0:0.1:1
    pts3 = 0:0.1:1
    grf = GaussianRandomField(cov,Spectral(),pts1,pts2,pts3)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 3
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts) == 3
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test length(grf.pts[3]) == length(pts3)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)
    @test size(sample(grf),3) == length(pts3)

    ## non-SPD covariance matrix ##
    cov = CovarianceFunction(2,SquaredExponential(0.5))
    pts1 = 0:0.05:1
    pts2 = 0:0.05:1
    @suppress grf = GaussianRandomField(cov,Spectral(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,SquaredExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## test non-equidistant structured grid  ##
    cov = CovarianceFunction(2,Matern(1.,2.5,σ=1.,p=2))
    pts1 = sin.(-pi/2:0.2:pi/2)
    pts2 = sin.(-pi/2:0.1:pi/2)
    grf = GaussianRandomField(cov,Spectral(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## test anisotropic
    cov = CovarianceFunction(2,AnisotropicExponential([1000 0; 0 1000]))
    pts1 = range(0,stop = 10,length = 64)
    pts2 = range(0,stop = 10,length = 64)
    grf = GaussianRandomField(cov,Spectral(),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,AnisotropicExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

end


text:
using Ito
using Base.Test

using Calendar
using Ito.Time

#2010 UK LSE Holidays
c=UKLSECalendar(); year=2010
@assert !isWorkingDay(c, ymd(2010, January, 1))
@assert !isWorkingDay(c, ymd(2010, April, 2))
@assert !isWorkingDay(c, ymd(2010, April, 5))
@assert !isWorkingDay(c, ymd(2010, May, 3))
@assert !isWorkingDay(c, ymd(2010, May, 31))
@assert !isWorkingDay(c, ymd(2010, August, 30))
@assert !isWorkingDay(c, ymd(2010, December, 25)) #Saturday
@assert !isWorkingDay(c, ymd(2010, December, 27))
@assert !isWorkingDay(c, ymd(2010, December, 28))

#US Calendars
@assert !isWorkingDay(USSettlementCalendar(), ymd(2012, July, 4))
@assert !isWorkingDay(USNYSECalendar(), ymd(2012, July, 4))
@assert !isWorkingDay(USGovernmentBondCalendar(), ymd(2012, July, 4))
@assert !isWorkingDay(USNERCCalendar(), ymd(2012, July, 4))

#2012 India LSE Calendar

c=IndiaNSECalendar()

@assert !isWorkingDay(c, ymd(2012, January, 1)) #New Year, Sunday
@assert !isWorkingDay(c, ymd(2012, January, 26)) #Republic Day
@assert !isWorkingDay(c, ymd(2012, February, 20)) #Shivratri
@assert !isWorkingDay(c, ymd(2012, March, 8)) #Holi
@assert !isWorkingDay(c, ymd(2012, April, 1)) #Ram Navami, Sunday
@assert !isWorkingDay(c, ymd(2012, April, 5)) #Mahavir Jayanti
@assert !isWorkingDay(c, ymd(2012, April, 6)) #Good Friday
@assert !isWorkingDay(c, ymd(2012, April, 14)) #Ambedkar Jayanti, Saturday
@assert !isWorkingDay(c, ymd(2012, May, 1)) #May Day
@assert !isWorkingDay(c, ymd(2012, August, 15)) #Independence Day
@assert !isWorkingDay(c, ymd(2012, August, 20)) #Ramzan
@assert !isWorkingDay(c, ymd(2012, September, 19)) #Ganesh Chaturthi
@assert !isWorkingDay(c, ymd(2012, October, 2)) #Gandhi Jayanti
@assert !isWorkingDay(c, ymd(2012, October, 24)) #Dussera
@assert !isWorkingDay(c, ymd(2012, October, 27)) #Bakri Id, Saturday
@assert !isWorkingDay(c, ymd(2012, November, 14)) #Diwali
@assert !isWorkingDay(c, ymd(2012, November, 25)) #Moharram, Sunday
@assert !isWorkingDay(c, ymd(2012, November, 28)) #Guru Nanak Jayanti
@assert !isWorkingDay(c, ymd(2012, December, 25)) #Christmas

#DayCount

@test_approx_eq 0.497724380567 yearfraction(ISDAActualActual(), ymd(2003,November,1), ymd(2004, May, 1))
@test_approx_eq 0.50 yearfraction(ISMAActualActual(), ymd(2003,November,1), ymd(2004, May, 1))
@test_approx_eq 0.497267759563 yearfraction(AFBActualActual(), ymd(2003,November,1), ymd(2004, May, 1))

@test_approx_eq 0.410958904110 yearfraction(ISDAActualActual(), ymd(1999,February,1), ymd(1999, July, 1))
@test_approx_eq 0.410958904110 yearfraction(ISMAActualActual(), ymd(1999,February,1), ymd(1999, July, 1), ymd(1998, July, 1), ymd(1999, July, 1))
@test_approx_eq 0.410958904110 yearfraction(AFBActualActual(), ymd(1999,February,1), ymd(1999, July, 1))

@test_approx_eq_eps 1.001377348600 yearfraction(ISDAActualActual(), ymd(1999,July,1), ymd(2000, July, 1)) 1e-10
@test_approx_eq 1 yearfraction(ISMAActualActual(), ymd(1999,July,1), ymd(2000, July, 1), ymd(1999, July, 1), ymd(2000, July, 1))
@test_approx_eq 1 yearfraction(AFBActualActual(), ymd(1999,July,1), ymd(2000, July, 1))

@test_approx_eq 0.915068493151 yearfraction(ISDAActualActual(), ymd(2002,August,15), ymd(2003, July, 15))
@test_approx_eq 0.915760869565 yearfraction(ISMAActualActual(), ymd(2002,August,15), ymd(2003, July, 15), ymd(2003, January, 15), ymd(2003, July, 15))
@test_approx_eq 0.915068493151 yearfraction(AFBActualActual(), ymd(2002,August,15), ymd(2003, July, 15))

@test_approx_eq 0.504004790778 yearfraction(ISDAActualActual(), ymd(2003,July,15), ymd(2004, January, 15))
@test_approx_eq 0.5 yearfraction(ISMAActualActual(), ymd(2003,July,15), ymd(2004, January, 15), ymd(2003, July, 15), ymd(2004, January, 15))
@test_approx_eq 0.504109589041 yearfraction(AFBActualActual(), ymd(2003,July,15), ymd(2004, January, 15))

@test_approx_eq 0.503892506924 yearfraction(ISDAActualActual(), ymd(1999,July,30), ymd(2000, January, 30))
@test_approx_eq 0.5 yearfraction(ISMAActualActual(), ymd(1999,July,30), ymd(2000, January, 30), ymd(1999,July,30), ymd(2000, January, 30))
@test_approx_eq 0.504109589041 yearfraction(AFBActualActual(), ymd(1999,July,30), ymd(2000, January, 30))

@test_approx_eq 0.415300546448 yearfraction(ISDAActualActual(), ymd(2000,January,30), ymd(2000, June, 30))
@test_approx_eq 0.417582417582 yearfraction(ISMAActualActual(), ymd(2000,January,30), ymd(2000, June, 30), ymd(2000,January,30), ymd(2000, July, 30))
@test_approx_eq 0.415300546448 yearfraction(AFBActualActual(), ymd(2000,January,30), ymd(2000, June, 30))


text:
"""
    _tmap(f, xs...)

An internal helper for `tmap` that handles the `basesize == 1` case.
"""
function _tmap(f, xss...)
    # Set the number of threads that BLAS uses to 1 here in case it's using up too many
    nthreads() > 1 && LinearAlgebra.BLAS.set_num_threads(1)
    futures = map(xss...) do x...
        Threads.@spawn f(x...)
    end
    results = map(fetch, futures)

     # Bump up the BLAS threads once we're done
    nthreads() > 1 && LinearAlgebra.BLAS.set_num_threads(typemax(Int32))

    return results
end

"""
    _tmap_with_partition(f, xss...; basesize)

An internal helper for `tmap` that handles the `basesize > 1` case.
Works for `basesize == 1`, but less efficent; since it breaks things up into single
item slices then stiches them back together again.
"""
function _tmap_with_partition(f, xss...; basesize)
    partitioned_xss = Iterators.partition.(xss, basesize)
    partioned_ys = _tmap(partitioned_xss...) do xss...
        map(f, xss...)
    end
    return reduce(vcat, partioned_ys)
end

"""
    tmap(f, xs...; basesize=1)

Multithreaded version of `map`.
`basesize` controls the minimum number of items from `xs` to process per `@spawn`ed task.

!!! tip
    `basesize` should be set high enough that proccessing that
    many items takes about ~1ms. This is to counter the ~50μs overhead
    it takes to dispatch work to a thread. If the function takes >1ms per
    call, then `basesize=1` is recommended.
"""
function tmap(f, xss...; basesize=1)
    if basesize == 1
        _tmap(f, xss...)
    else
        _tmap_with_partition(f, xss...; basesize=basesize)
    end
end

"""
    tmap_with_warmup(f, xs...; basesize=1)

Similar to [`tmap`](@ref), but runs the first call single threaded, before multithreading
the remainnder.
This is useful for dealing with things that benifit from something happening on first run.
Which might be related to caching values, or some compilation troubles.

`basesize` controls the minimum number of items from `xs` to process per `@spawn`ed task.

See [`tmap`](@ref) for more details
"""
function tmap_with_warmup(f, xss...; basesize=1)
    xs = first.(xss)
    xs_tails = Iterators.rest.(xss, 2)

    y = f(xs...)
    ys_tail = tmap(f, xs_tails...; basesize=basesize)
    ys = pushfirst!(ys_tail, y)
    return ys
end


text:
# Support functionality for amend_coverage_from_src!

isevaldef(x) = Base.Meta.isexpr(x, :(=)) && Base.Meta.isexpr(x.args[1], :call) &&
               x.args[1].args[1] == :eval

function_body_lines(ast) = function_body_lines!(Int[], ast, false)
function_body_lines!(flines, arg, infunction) = flines
function function_body_lines!(flines, node::LineNumberNode, infunction)
    line = node.line
    if infunction
        push!(flines, line)
    end
    flines
end
function function_body_lines!(flines, ast::Expr, infunction)
    if ast.head == :line
        line = ast.args[1]
        if infunction
            push!(flines, line)
        end
        return flines
    elseif ast.head == :module
        # Ignore automatically added eval definitions
        args = ast.args[end].args
        if length(args) >= 2 && isevaldef(args[1]) && isevaldef(args[2])
            args = args[3:end]
        end
    else
        args = ast.args
    end
    infunction |= isfuncexpr(ast)
    for arg in args
        flines = function_body_lines!(flines, arg, infunction)
    end
    flines
end


text:
using subpop
using Test

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


text:
#' ---
#' title: Variables, Constraints and Objective
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' While the last tutorial introduced you to basics of of JuMP code, this tutorial will go in depth focusing on how to work
#' with different parts of a JuMP program.

using JuMP
model = Model();

#' ## Variables

#' ### Variable Bounds
#' All of the variables we have created till now have had a bound. We can also create a free variable.

@variable(model, free_x)

#' While creating a variable, instead of using the <= and >= syntax, we can also use the `lower_bound` and `upper_bound` keyword arguments.

@variable(model, keyword_x, lower_bound = 1, upper_bound = 2)

#' We can query whether a variable has a bound using the `has_lower_bound` and `has_upper_bound` functions. The values of the bound can be obtained
#' using the `lower_bound` and `upper_bound` functions.

has_upper_bound(keyword_x)

#+

upper_bound(keyword_x)

#' Note querying the value of a bound that does not exist will result in an error.
#+ tangle = false

lower_bound(free_x)

#' JuMP also allows us to change the bounds on variable. We will learn this in the problem modification tutorial.

#' ### Containers
#' We have already seen how to add a single variable to a model using the `@variable` macro. Let's now look at more ways to add
#' variables to a JuMP model. JuMP provides data structures for adding collections of variables to a model. These data
#' structures are reffered to as Containers and are of three types - `Arrays`, `DenseAxisArrays`, and `SparseAxisArrays`.

#' #### Arrays
#' JuMP arrays are created in a similar syntax to Julia arrays with the addition of specifying that the indices start with 1. If
#' we do not tell JuMP that the indices start at 1, it will create a DenseAxisArray instead.

@variable(model, a[1:2, 1:2])

#' An n-dimensional variable $x \in {R}^n$ having a bound $l \preceq x \preceq u$ ($l, u \in {R}^n$) is added in the following
#' manner.

n = 10
l = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10]
u = [10; 11; 12; 13; 14; 15; 16; 17; 18; 19]

@variable(model, l[i] <= x[i = 1:n] <= u[i])

#' Note that while working with Containers, we can also create variable bounds depending upon the indices

@variable(model, y[i = 1:2, j = 1:2] >= 2i + j)

#' #### DenseAxisArrays
#' DenseAxisArrays are used when the required indices are not one-based integer ranges. The syntax is similar except with an
#' arbitrary vector as an index as opposed to a one-based range.

#' An example where the indices are integers but do not start with one.

@variable(model, z[i = 2:3, j = 1:2:3] >= 0)

#' Another example where the indices are an arbitrary vector.

@variable(model, w[1:5,["red", "blue"]] <= 1)

#' #### SparseAxisArrays
#' SparseAxisArrays are created when the indices do not form a rectangular set. For example, this applies when indices have a
#' dependence upon previous indices (called triangular indexing).

@variable(model, u[i = 1:3, j = i:5])

#' We can also conditionally create variables by adding a comparison check that depends upon the named indices and is separated
#' from the indices by a semi-colon (;).

@variable(model, v[i = 1:9; mod(i, 3) == 0])

#' ### Variable Types

#' The last arguement to the `@variable` macro is usually the variable type. Here we'll look at how to specifiy he variable type.

#' #### Integer Variables
#' Integer optimization variables are constrained to the set $x \in {Z}$
#+ tangle = false

@variable(model, integer_x, Int)

#' or

@variable(model, integer_z, integer = true)

#' #### Binary Variables
#' Binary optimization variables are constrained to the set $x \in \{0, 1\}$.
#+ tangle = false

@variable(model, binary_x, Bin)

#' or

@variable(model, binary_z, binary = true)

#' #### Semidefinite variables
#' JuMP also supports modeling with semidefinite variables. A square symmetric matrix X is positive semidefinite if all eigenvalues
#' are nonnegative.

@variable(model, psd_x[1:2, 1:2], PSD)

#' We can also impose a weaker constraint that the square matrix is only symmetric (instead of positive semidefinite) as follows:

@variable(model, sym_x[1:2, 1:2], Symmetric)

#' ## Constraints

model = Model()
@variable(model, x)
@variable(model, y)
@variable(model, z[1:10]);

#' ### Constraint References
#' While calling the `@constraint` macro, we can also set up a constraint reference. Such a refference is useful for obtaining
#' additional information about the constraint such as its dual.

@constraint(model, con, x <= 4)

#' ### Containers
#' Just as we had containers for variables, JuMP also provides `Arrays`, `DenseAxisArrays`, and `SparseAxisArrays` for storing
#' collections of constraints. Examples for each container type are given below.

#' #### Arrays

@constraint(model, [i = 1:3], i * x <= i + 1)

#' #### DenseAxisArrays

@constraint(model, [i = 1:2, j = 2:3], i * x <= j + 1)

#' #### SparseAxisArrays

@constraint(model, [i = 1:2, j = 1:2; i != j], i * x <= j + 1)

#' ### Constraints in a Loop
#' We can add constraints using regular Julia loops

for i in 1:3
    @constraint(model, 6x + 4y >= 5i)
end

#' or use for each loops inside the `@constraint` macro.

@constraint(model, [i in 1:3], 6x + 4y >= 5i)

#' We can also create constraints such as $\sum _{i = 1}^{10} z_i \leq 1$

@constraint(model, sum(z[i] for i in 1:10) <= 1)

#' ## Objective
#' While the recommended way to set the objective is with the @objective macro, the functions `set_objective_sense` and
#' `set_objective_function` provide an equivalent lower-level interface.

using GLPK

model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
set_objective_sense(model, MOI.MIN_SENSE)
set_objective_function(model, x + y)

optimize!(model)

@show objective_value(model);

#' To query the objective function from a model, we use the `objective_sense`, `objective_function`, and `objective_function_type`
#' functions.

objective_sense(model)

#+

objective_function(model)

#+

objective_function_type(model)

#' # Vectorized Constraints and Objective
#' We can also add constraints and objective to JuMP using vectorized linear algebra. We'll illustrate this by solving an LP in
#' standard form i.e.

#' $$
#' \begin{align*}
#' & \min & c^T x \\
#' & \;\;\text{s.t.} & A x = b \\
#' & & x \succeq 0 \\
#' & & x \in \mathbb{R}^n
#' \end{align*}
#' $$

vector_model = Model(GLPK.Optimizer)

A= [ 1 1 9  5;
     3 5 0  8;
     2 0 6 13]

b = [7; 3; 5]

c = [1; 3; 5; 2]

@variable(vector_model, x[1:4] >= 0)
@constraint(vector_model, A * x .== b)
@objective(vector_model, Min, c' * x)

optimize!(vector_model)

@show objective_value(vector_model);


text:
module TestTrials

using Test
using TreeParzen

space = Dict(:fun => HP.Uniform(:param1, 0.0, 5.0))

function_fail = Dict(
    :fun => (
        HP.Uniform(:param1, 0.0, 5.0),
        HP.Uniform(:param2, 0.0, 1.0)
    )
)
config = Config(0.25, 25, 24, 1001, 1.0)

trials = TreeParzen.API.run(TreeParzen.Trials.Trial[], x -> x[:fun], space, 3, config)

@test isa(trials[1], TreeParzen.Trials.Trial)
@test keys(trials[1].hyperparams) == keys(trials[2].hyperparams)
@test keys(trials[1].vals) == keys(trials[2].vals)

# Catches error when an invalid function returns a Tuple{Float64,Float64} where Float64 is
# expected
@test_throws TypeError TreeParzen.API.run(TreeParzen.Trials.Trial[], x ->[:fun], function_fail, 1, config)

end #module TestTrials
true


text:
using StringDistances, Random
Random.seed!(2)
x = map(Random.randstring, rand(5:25,500_000))
y = map(Random.randstring, rand(5:25,500_000))

function f(t, x, y; min_score = 0.0)
    [compare(x[i], y[i], t; min_score = min_score) for i in 1:length(x)]
end

function g(dist, x, y)
    [dist(x[i], y[i]) for i in 1:length(x)]
end

@time f(Jaro(), x, y);
#0.3s
@time f(Levenshtein(), x, y);
# 0.4s
@time f(Levenshtein(), x, y, min_score = 0.8);
# 0.11
@time f(DamerauLevenshtein(), x, y);
# 0.58s.
@time f(DamerauLevenshtein(), x, y, min_score = 0.8);
# 0.08 (now 0.09)
@time f(RatcliffObershelp(), x, y);
# 1.35s

@time findnearest(x[1], y, Levenshtein());
# 0.02
@time findnearest(x[1], y, DamerauLevenshtein());
# 0.05
@time findnearest(x[1], y, QGram(2));
# 0.75

@time findall(x[1], y, Levenshtein());
# 0.05
@time findall(x[1], y, DamerauLevenshtein());
# 0.05
@time findall(x[1], y, Partial(DamerauLevenshtein()));
# 0.96
@time findall(x[1], y, QGram(2));
# 0.81
@time findall(x[1], y, TokenSort(DamerauLevenshtein()));
# 0.27 (now 0.32)
@time findall(x[1], y, TokenSet(DamerauLevenshtein()));
# 0.55
@time findall(x[1], y, TokenMax(DamerauLevenshtein()));
# 2.25 (now 3.6)

x = map(Random.randstring, rand(5:25,1000))
y = map(Random.randstring, rand(5:25,1000))
@time pairwise(Levenshtein(), x, y)
# 0.25 seconds
@time pairwise(QGram(2), x, y, preprocess = false)
# 2.126829
@time pairwise(QGram(2), x, y, preprocess = true)
# 0.12

#= Rcode
library(stringdist)
x <- sapply(sample(5:25,5 * 1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
y <- sapply(sample(5:25,5 * 1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
system.time(stringdist(x,y,method='lv', nthread = 1))
system.time(stringdist(x,y,method='dl', nthread = 1))
#  0.472
system.time(stringdist(x,y,method='jaccard', nthread = 1))
# 0.739
system.time(stringdist(x,y,method='cosine', nthread = 1))
system.time(stringdist(x,y,method='qgram', nthread = 1))

=#


text:
export rand_supremacy2d

# control-Z entangler.
cz_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>Z) for (ctrl, target) in pairs)

"""
    rand_supremacy2d(nx::Int, ny::Int, depth::Int) -> AbstactBlock

The circuit proposed for realizing quantum supermacy in a near-term device.

References
-------------------------------
* Boixo, Sergio, et al. "Characterizing quantum supremacy in near-term devices." Nature Physics 14.6 (2018): 595-600.

!!! note

    Some restrictions are loosed, please check this circuit carefully.
"""
function rand_supremacy2d(nx::Int, ny::Int, depth::Int)
    nbits = nx*ny
    entanglers = map(pair->cz_entangler(nbits, pair), pair_supremacy(nx, ny))
    gateset = [ConstGate.T, SqrtX, SqrtY]
    c = chain(nbits, repeat(nbits, H, 1:nbits))
    pre_occ = (1:nbits...,)
    pre_gates = Vector{PrimitiveBlock}(fill(H,nbits))
    hastgate = zeros(Bool, nbits)
    for i=1:depth-2
        ent = entanglers[mod1(i,length(entanglers))]
        unit = chain(nbits, [ent])
        occ = occupied_locs(ent)
        for loc in setdiff(pre_occ, occ)
            g1 = hastgate[loc] ? rand(setdiff(gateset, (pre_gates[loc],))) : (hastgate[loc]=true; T)
            push!(unit, put(nbits, loc=>g1))
            pre_gates[loc] = g1
        end
        pre_occ = occ
        push!(c, unit)
    end
    depth!=1 && push!(c, repeat(nbits, H, 1:nbits))
    return c
end

# obtain supremacy pairing patterns
function pair_supremacy(nx::Int, ny::Int; periodic=false)
    Kx = [0.25 0.5; 0.25 -0.5]
    Ky = [0.5 0.25; -0.5 0.25]
    out = Vector[]
    li = LinearIndices((nx, ny))

    for (dx, dy, isxdirection) in [(2,0,true), (0,0,true), (0,3,false), (0,1,false),
                                    (3,0,true), (1,0,true), (0,0,false), (0,2,false)]
        res = Pair{Int, Int}[]
        for i = 1:nx, j=1:ny
            if (periodic || (isxdirection ? i<nx : j<ny)) && all(mod.((isxdirection ? Kx : Ky)*[i-dx-1, j-dy-1], 1) .≈ 0)
                i_ = isxdirection ? mod1(i+1, nx) : i
                j_ = isxdirection ? j : mod1(j+1, ny)
                push!(res, li[i, j] => li[i_, j_])
            end
        end
        push!(out, res)
    end

    return out
end


text:
@testset "Inner products" begin
    @testset "Vectors" begin
        t = range(0,stop=1,length=11)
        B = FEDVR(t, 4)
        B̃ = B[:,1:5]

        uv = ones(size(B,2))
        vv = ones(size(B̃,2))

        @testset "Direct inner products" begin
            u = B*uv
            v = B̃*vv

            @test u'v == 5
            @test (uv'*(B'B̃)*vv)[1] == 5
            @test norm(v) ≈ √5
        end
    end

    @testset "Matrices" begin
        t = range(-1.0,stop=1.0,length=5)
        R = FEDVR(t, 4)[:,2:end-1]

        n = 3
        Φ = rand(ComplexF64, size(R,2), n)
        for i = 1:n
            # Orthogonalize against all previous vectors
            for j = 1:i-1
                c = Φ[:,j]'Φ[:,i]/(Φ[:,j]'Φ[:,j])
                Φ[:,i] -= c*Φ[:,j]
            end
            Φ[:,i] /= norm(R*Φ[:,i])
        end

        @testset "Direct inner products" begin
            Φ̂ = R*Φ
            Φ̂'Φ̂
            @test Φ̂'Φ̂ ≈ I
        end
    end
end


text:
module BackwardKernel

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

using Distributions, SpecialFunctions
import DataStructures: PriorityQueue, peek, enqueue!, dequeue!
import Base: iterate, length, sum, eltype
import Base.Iterators: product

export EventQueue, sum, length, eltype, iterate

mutable struct EventQueue{ℜ<:Real, N}
    T::Int
    X::Array{ℜ,N}
    pq::PriorityQueue{NTuple{N,Int},Float64,Base.Order.ForwardOrdering}
    function EventQueue(X::Array{ℜ,N}) where {ℜ<:Real,N}
        I = size(X)
        T::Int = Int(round(sum(X)))
        pq = PriorityQueue([i => rand(Dirichlet([1.0,X[i...]]))[1]
                           for i ∈ product(map(Iₙ -> 1:Iₙ, I)...) if X[i...]>0])
        return new{ℜ,N}(T,copy(X),pq)
    end
end

Base.sum(L::EventQueue) = sum(L.X)
Base.length(L::EventQueue) = L.T
Base.eltype(L::EventQueue{ℜ,N}) where {ℜ<:Real, N} = NTuple{N,Integer}
Base.size(L::EventQueue) = size(L.X)

function Base.iterate(L::EventQueue, state::Ƶ=1) where {Ƶ<:Integer}
    if L.T < state
        return nothing
    end

    t = peek(L.pq)[2]
    i = dequeue!(L.pq)

    L.X[i...] -= 1
    if L.X[i...] > 0.0
        L.pq[i] = t + (1.0 - t)*rand(Dirichlet([1.0, L.X[i...]]))[1]
    end

    return i, state+1
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Stage-decomposition spconstructs #
# ================================ #
function EWS(stochasticprogram::StochasticProgram,
             structure::StageDecompositionStructure{2,1,Tuple{SP}}) where SP <: DistributedScenarioProblems
    partial_ews = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch(
                w,
                scenarioproblems(structure, 2)[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                subproblem_optimizer(stochasticprogram)) do sp, gen_one, gen_two, one_params, two_params, opt
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0
                    subproblems = map(scenarios(scenarioproblems)) do scenario
                        ws = _WS(gen_one,
                                 gen_two,
                                 one_params,
                                 two_params,
                                 scenario,
                                 DecisionMap(),
                                 DecisionMap(),
                                 opt)
                        return ws
                    end
                    return outcome_mean(subproblems, probability.(scenarios(scenarioproblems)))
                end
        end
    end
    return sum(partial_ews)
end

function statistical_EWS(stochasticprogram::StochasticProgram,
                         structure::StageDecompositionStructure{2,1,Tuple{SP}}) where SP <: DistributedScenarioProblems
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch(
                w,
                scenarioproblems(structure, 2)[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                subproblem_optimizer(stochasticprogram)) do sp, gen_one, gen_two, one_params, two_params, opt
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0, 0.0, 0.0, 0
                    ws_models = map(scenarios(scenarioproblems)) do scenario
                        ws = _WS(gen_one,
                                 gen_two,
                                 one_params,
                                 two_params,
                                 scenario,
                                 DecisionMap(),
                                 DecisionMap(),
                                 opt)
                        return ws
                    end
                    return welford(ws_models, probability.(scenarios(scenarioproblems)))
                end
        end
    end
    𝔼WS, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼WS, sqrt(σ²)
end


text:
import LearningStrategies
const L = LearningStrategies

export GAStrategy

mutable struct GAStrategy{G, P, K,
                          Fs<:SelectionOperator{>:G, P, K},
                          Fc<:CrossoverOperator{>:G, K, P},
                          Fm<:MutationOperator{>:G}} <: L.LearningStrategy
    selection::Fs
    crossover::Fc
    mutation::Fm
    cache::Population{G}

    function GAStrategy(selection::SelectionOperator{U, P, K},
                        crossover::CrossoverOperator{V, K, P},
                        mutation::MutationOperator{W}) where {U, V, W, P, K}
        G = typejoin(U, V, W)
        S, C, M = typeof(selection), typeof(crossover), typeof(mutation)
        new{G, P, K, S, C, M}(selection, crossover, mutation, Population{G}())
    end

    function GAStrategy{G}(selection::SelectionOperator{>:G, P, K},
                           crossover::CrossoverOperator{>:G, K, P},
                           mutation::MutationOperator{>:G}) where {G, P, K}
        S, C, M = typeof(selection), typeof(crossover), typeof(mutation)
        new{G, P, K, S, C, M}(selection, crossover, mutation, Population{G}())
    end
end

function L.setup!(strategy::GAStrategy{G}, model::AbstractPopulationModel{G}) where G
    # setup!(model.fitness, model)
    setup!(strategy.selection, model)
    setup!(strategy.mutation, model)
    setup!(strategy.crossover, model)
    assessfitness!(model)
end

function L.update!(model::AbstractPopulationModel{G}, strategy::GAStrategy{G, P, K},
                   generation, _item) where {G, P, K}
    preparecache!(strategy.cache, length(model.population))

    # TODO: log timing information
    for parents::Family{G, K} in selection(model.population, strategy.selection, generation)
        for child in crossover!(copy.(parents), strategy.crossover, generation)::Family{G, P}
            push!(strategy.cache, mutate!(child, strategy.mutation, generation))
        end
    end

    # swap parents and children -- saves reallocations
    model.population, strategy.cache = strategy.cache, model.population
    assessfitness!(model)

    model
end

function L.finished(verbose_strategy::L.Verbose{<:GAStrategy}, model::AbstractPopulationModel,
                    data, generation)
    done = L.finished(verbose_strategy.strategy, model, data, generation)
    done && @info "Evolved $generation generations, final population size $(length(model.population))"
    done
end


text:
using Documenter, ShapML

makedocs(
    sitename = "ShapML",
    authors = "Nickalus Redell",
    doctest = false,
    format = Documenter.HTML(
        highlights = ["yaml"]
        ),
    pages = [
        "Introduction" => "index.md",
        "Vignettes" => Any["Stochastic vs. TreeSHAP" => "vignettes/consistency.md"],
        "Functions" => "functions/functions.md"
    ]
)

deploydocs(
    repo = "github.com/nredell/ShapML.jl.git"
)


text:
using PCM
using Test


text:
type State
    tokens::Vector{Int}
    headsets::Vector{Vector{Int}}
    step::Int
end

State(tokens, headsets) = State(tokens, headsets, 0)

function Base.next(s::State)::Vector{State}
    s.step == length(tokens) && return State[]
    map(headsets[s.step+1]) do h
        State(s.tokens, s.headsets, h)
    end
end

type Model
    nn
end

function Model()
    T = Float32
    x = Var()
    y = Embedding(T,10000,100)(x)
    y = Linear(T,100,1)(y)
    Model(Graph(y,x))
end

function getscore(m::Model, s::State)::Float64
    1.0
end


text:
function coords!(xsi, x1, y1, x2, y2, orient, x, incx, y, incy)

    Np = length(xsi)
    s = orient?1.0:-1.0
    idx = 1; idy = 1
    for i = 1:Np
        xx = s * xsi[i]
        x[idx] = (1-xx)/2 * x1 + (1+xx)/2 * x2
        y[idy] = (1-xx)/2 * y1 + (1+xx)/2 * y2
        idx += incx
        idy += incy
    end

end

function coords(xsi, x1, y1, x2, y2, orient=true)
    Np = length(xsi)
    x = zeros(Np)
    y = zeros(Np)
    coords!(xsi, x1, y1, x2, y2, orient, x, 1, y, 1)

    return x, y
end

coords(s::Segment, xsi, x1, y1, x2, y2, orient=true) = coords(xsi, x1, y1, x2, y2, orient)
coords!(s::Segment, xsi, x1, y1, x2, y2, orient, x, incx, y, incy) = coords!(xsi, x1, y1, x2, y2, orient, x, incx, y, incy)

function normal(x1, y1, x2, y2)

    dx = x2-x1
    dy = y2-y1

    len = hypot(dx, dy)

    return [dy/len, -dx/len]
end

normal(s::Segment, x, y) = normal(s.x1, s.y1, s.x2, s.y2)

function coords!(crc::Circle, xsi, x1, y1, x2, y2, orient, x, incx, y, incy)
    Np = length(xsi)
    const PI2 = 2*pi
    xc = crc.xc
    yc = crc.yc
    r = crc.rad
    inside = crc.inside

    xr1 = x1 - xc; yr1 = y1-yc
    xr2 = x2 - xc; yr2 = y2-yc

    th1 = atan2(yr1, xr1)
    th2 = atan2(yr2, xr2)

    if inside
        if th2 < th1
            th2 += PI2
        end
    else
        if th2 > th1
            th2 -= PI2
        end
    end

    s = orient?1.0:-1.0

    idx = 1; idy = 1
    for i = 1:Np
        xx = s * xsi[i]
        th = ( (1-xx)*th1 + (1+xx)*th2 ) / 2
        x[idx] = xc + r * cos(th)
        y[idy] = yc + r * sin(th)
        idx += incx
        idy += incy
    end

end

function normal(crc::Circle, x, y)

    theta = atan2(y, x)

    s = crc.inside?-1.0:1.0

    return [s*cos(theta), s*sin(theta)]
end

function coords(crc::Circle, xsi, x1, y1, x2, y2, orient=true)
    Np = length(xsi)
    x = zeros(Np)
    y = zeros(Np)
    coords!(crc::Circle, xsi, x1, y1, x2, y2, orient, x, 1, y, 1)
    return x,y
end

function gordon_hall(x, y, xsi)
    Q = length(xsi)
    Q2 = Q*Q
    vx = zeros(Q,Q)
    vy = zeros(Q,Q)
    ex = zeros(Q,Q)
    ey = zeros(Q,Q)

    x0 = x[1,1]; x1 = x[Q,1]; x2 = x[Q,Q]; x3 = x[1,Q]
    y0 = y[1,1]; y1 = y[Q,1]; y2 = y[Q,Q]; y3 = y[1,Q]

    c = 1
    for j=1:Q, i=1:Q
        xi = xsi[i]
        xj = xsi[j]
	vx[c] = 0.25 * ( (1.0-xi)*(1.0-xj) * x0 + (1.0+xi)*(1.0-xj)*x1 +
		      (1.0+xi)*(1.0+xj)*x2 + (1.0-xi)*(1.0+xj)*x3)
	vy[c] = 0.25 * ( (1.0-xi)*(1.0-xj) * y0 + (1.0+xi)*(1.0-xj)*y1 +
		      (1.0+xi)*(1.0+xj)*y2 + (1.0-xi)*(1.0+xj)*y3);
        c += 1
    end
    c = 1
    for j=1:Q, i=1:Q
        xi = xsi[i]
        xj = xsi[j]
	ex[c] = vx[c] + 0.5*( (1.0-xi)*(x[1,j] - vx[1,j]) +
			     (1.0+xi)*(x[Q,j] - vx[Q,j]) +
			     (1.0-xj)*(x[i,1] - vx[i]) +
			     (1.0+xj)*(x[i,Q] - vx[i,Q]))

	ey[c] = vy[c] + 0.5*( (1.0-xi)*(y[1,j] - vy[1,j]) +
			     (1.0+xi)*(y[Q,j] - vy[Q,j]) +
			     (1.0-xj)*(y[i,1] - vy[i]) +
			     (1.0+xj)*(y[i,Q] - vy[i,Q]))
        c += 1
    end

    for i = 1:Q2
        x[i] = ex[i]
        y[i] = ey[i]
    end
end

function node_coords(elem::Quad, xsi, crvs)

    Q = length(xsi)

    x = zeros(Q,Q)
    y = zeros(Q,Q)
    xv = xverts(elem)
    yv = yverts(elem)
    ic = elem.curve
    if ic[1] > 0
        x[:,1], y[:,1] = coords(crvs[ic[1]], xsi, xv[1], yv[1], xv[2], yv[2], true)
    else
        x[:,1], y[:,1] = coords(xsi, xv[1], yv[1], xv[2], yv[2], true)
    end

    if ic[2] > 0
        x[Q,:], y[Q,:] = coords(crvs[ic[2]], xsi, xv[2], yv[2], xv[3], yv[3], true)
    else
        x[Q,:], y[Q,:] = coords(xsi, xv[2], yv[2], xv[3], yv[3], true)
    end

    if ic[3] > 0
        x[:,Q], y[:,Q] = coords(crvs[ic[3]], xsi, xv[3], yv[3], xv[4], yv[4], false)
    else
        x[:,Q], y[:,Q] = coords(xsi, xv[3], yv[3], xv[4], yv[4], false)
    end

    if ic[4] > 0
        x[1,:], y[1,:] = coords(crvs[ic[4]], xsi, xv[4], yv[4], xv[3], yv[3], false)
    else
        x[1,:], y[1,:] = coords(xsi, xv[4], yv[4], xv[1], yv[1], false)
    end
    gordon_hall(Q, x, y, xsi)

    return x, y
end


text:
##########
##########
# Functions to perform perturbations on the food webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

using DifferentialEquations
using DelimitedFiles

# functions ====================================================================

"""
    biomass_perturb_all(du, u, p, t)

Calling biomass functions, create the actual equation for each species, to be
passed to the simulation. Do this by perturbing all species.
"""
function biomass_perturb_all(du, u, p, t)
    # pick set of random
    # extract all values to be used w/in function from the p value
    species_list, network, base_mortality, z = p
    for i in 1:size(species_list, 1)
        du[i] = growth_rate(i, species_list, u) +
                biomass_gain(species_list, network, u, 2, i) -
                biomass_loss(species_list, network, u, 2, i) -
                mortality_change(i, u, base_mortality, z)
    end
end

"""
    biomass_perturb_random(du, u, p, t)

Calling biomass functions, create the actual equation for each species, to be
passed to the simulation. Do this by randomly picking ~1/3 of the species and
changing the mortality rate. The species to be perturbed are defined before this
function is used, as are the ones NOT to be perturbed.
"""
function biomass_perturb_random!(du, u, p, t)
    # extract all values to be used w/in function from the p value
    # extract all values to be used w/in function from the p value
    species_list, network, base_mortality, amount, perturb_sp, non_perturb = p
    if amount === nothing
        for i in 1:size(species_list, 1)
            du[i] = growth_rate(i, species_list, u) +
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        end
    else
        for i in perturb_sp
            du[i] = growth_rate(i, species_list, u) +
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        end
        for i in non_perturb
            du[i] = growth_rate(i, species_list, u) +
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        end
    end
end

"""
    extinction_cond(out, x, t, integrator)

Condition asking if any of the values in the integrator go below the extinction
threshold.
"""
function extinction_threshold(out,u,t,integrator)
    # loop through all species to make the condition check all of them
    for i in 1:size(species_list, 1)
        out[i] = 1e-05 - u[i]
    end
end

"""
    extinction_affect!(integrator, event_idx)

Affect function to use in callback vector, this will allow us to check
if any of the values of the threshold condition are true, and if so,
set the value to zero in the integrator.
"""
function extinction_affect!(integrator, event_idx)
    # loop again through all species
    for i in 1:size(species_list, 1)
        if event_idx == i
            integrator.u[i] = 0
        end
    end
end

"""
    inc_mort_cond(u,t,integrator)

Condition to find the stop value and when it's hit, increase the parameter value
specified. Do this whenever a value of stop1 is reached.
"""
function inc_mort_cond(u,t,integrator)
  t in stop1
end

"""
    dec_mort_cond(u,t,integrator)

Condition to find the stop value and when it's hit, decrease the parameter value
specified. Do this whenever a value of stop2 is reached.
"""
function dec_mort_cond(u,t,integrator)
  t in stop2
end

"""
    affect_inc!(integrator)

Affect for the condition to work on. Increase the mortality by making the amount
(the fourth element in p) equal to z.
"""
function affect_inc!(integrator, z)
  integrator.p[4] = z
end

"""
    affect_dec!(integrator)

Affect for the condition to work on. Decrease the fourth element in p,
returning it to 1.0.
"""
function affect_dec!(integrator)
  integrator.p[4] = nothing
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
module SoilWater

using Cropbox

@system Pedotransfer begin
    Ψ_wp: tension_wilting_point => 1500 ~ preserve(u"kPa", parameter)
    Ψ_fc: tension_field_capacity => 33 ~ preserve(u"kPa", parameter)
    Ψ_sat: tension_saturation => 0.01 ~ preserve(u"kPa", parameter)

    θ_wp: vwc_wilting_point ~ hold
    θ_fc: vwc_field_capacity ~ hold
    θ_sat: vwc_saturation ~ hold

    K_at(; vwc): hydraulic_conductivity_at ~ hold
    Hm_at: matric_head_at ~ hold
end

@system TabularPedotransfer(Pedotransfer) begin
    vwc2hc => [
        0.005 3.46e-9;
        0.050 4.32e-8;
        0.100 1.3e-7;
        0.150 6.91e-7;
        0.200 4.32e-6;
        0.250 2.59e-5;
        0.300 0.000173;
        0.350 0.001037;
        0.400 0.006912;
        0.450 0.0432;
        1.000 0.0432;
    ] ~ interpolate(u"m/d")
    #hc_to_vwc(vwc_to_hc) ~ interpolate(reverse)

    vwc2ss => [
        0.005 10000;
        0.010 3500;
        0.025 1000;
        0.050 200;
        0.100 40;
        0.150 10;
        0.200 6;
        0.250 3.5;
        0.300 2.2;
        0.350 1.4;
        0.400 0.56;
        0.450 0.001; #HACK: no duplicate 0
        1.000 0;
    ] ~ interpolate(u"m")
    ss2vwc(vwc2ss) ~ interpolate(reverse)

    head(; Ψ(u"kPa")) => (Ψ * u"m" / 9.8041u"kPa") ~ call(u"m")

    θ_wp(ss2vwc, head, Ψ_wp): vwc_wilting_point => ss2vwc(head(Ψ_wp)) ~ preserve # 0.02? 0.06
    θ_fc(ss2vwc, head, Ψ_fc): vwc_field_capacity => ss2vwc(head(Ψ_fc)) ~ preserve # 0.11? 0.26
    θ_sat(ss2vwc, head, Ψ_sat): vwc_saturation => ss2vwc(head(Ψ_sat)) ~ preserve # 0.45

    K_at(vwc2hc; θ): hydraulic_conductivity_at => vwc2hc(θ) ~ call(u"m/d")
    Hm_at(vwc2ss; θ): matric_head_at => vwc_ss(θ) ~ call(u"m")

    # vwc_airdry_water => 0.01 ~ preserve(parameter)
    # vwc_wilting_point => 0.07 ~ preserve(parameter)
    # initial_vwc => 0.4 ~ preserve(parameter)
    # rooting_depth => 0.2 ~ preserve(u"m", parameter)
    # iteration_per_time_step => 100 ~ preserve(parameter)
end

@system Texture begin
    S: sand => 0.29 ~ preserve(parameter)
    C: clay => 0.32 ~ preserve(parameter)
    OM: organic_matter => 1.5 ~ preserve(u"percent", parameter)
end

@system CharacteristicTransfer(Pedotransfer, Texture) begin
    DF: density_factor => 1.0 ~ preserve(parameter)

    # volumetric soil water content at permanent wilting point

    θ_1500t(S, C, OM): vwc_1500_first => begin
        -0.024S + 0.487C + 0.006OM + 0.005S*OM - 0.013C*OM + 0.068S*C + 0.031
    end ~ preserve # theta_1500t (%v)

    θ_1500(θ=θ_1500t): vwc_1500 => begin
        θ + (0.14θ - 0.02)
    end ~ preserve # theta_1500 (%v)

    θ_wp(θ_1500): vwc_wilting_point ~ preserve

    # volumetric soil water content at field capacity

    θ_33t(S, C, OM): vwc_33_first => begin
        -0.251S + 0.195C + 0.011OM + 0.006S*OM - 0.027C*OM + 0.452S*C + 0.299
    end ~ preserve # theta_33t (%v)

    θ_33(θ=θ_33t): vwc_33_normal => begin
        θ + (1.283θ^2 - 0.374θ - 0.015)
    end ~ preserve # theta_33 (%v)

    θ_33_DF(θ_33, θ_s, θ_s_DF): vwc_33_adjusted => begin
        θ_33 - 0.2(θ_s - θ_s_DF)
    end ~ preserve # theta_33_DF (%v)

    θ_fc(θ_33_DF): vwc_field_capacity ~ preserve

    # volumetric soil water content between saturation and field capacity

    θ_s_33t(S, C, OM): vwc_gravitation_first => begin
        0.278S + 0.034C +0.022OM - 0.018S*OM - 0.027C*OM - 0.584S*C + 0.078
    end ~ preserve # theta_s_33t (%v)

    θ_s_33(θ=θ_s_33t): vwc_gravitation_normal => begin
        θ + (0.636θ - 0.107)
    end ~ preserve # theta_s_33 (%v)

    θ_s_33_DF(θ_s_DF, θ_33_DF): vwc_gravitation_adjusted => begin
        θ_s_DF - θ_33_DF
    end ~ preserve # theta_s_33_DF (%v)

    # volumetric soil water content at saturation

    θ_s(θ_33, θ_s_33, S): vwc_saturation_normal => begin
        θ_33 + θ_s_33 - 0.097S + 0.043
    end ~ preserve # theta_s (%v)

    θ_s_DF(θ_s, ρ_DF, ρ_P): vwc_saturation_adjusted => begin
        1 - ρ_DF / ρ_P
    end ~ preserve # theta_s_DF (%v)

    θ_sat(θ_s_DF): vwc_saturation ~ preserve

    # density effects

    ρ_DF(ρ_N, DF): matric_density => begin
        ρ_N * DF
    end ~ preserve(u"g/cm^3") # rho_DF (g cm-3)

    ρ_N(θ_s, ρ_P): normal_density => begin
        (1 - θ_s) * ρ_P
    end ~ preserve(u"g/cm^3") # rho_N (g cm-3)

    ρ_P: particle_density => begin
        2.65
    end ~ preserve(u"g/cm^3") # (g cm-3)

    # hydraulic conductivity (moisture - conductivity)

    # coefficients of moisture-tension, Eq. 11 of Saxton and Rawls 2006
    A(B, θ_33): moisture_tension_curve_coeff_A => begin
        exp(log(33) + B*log(θ_33))
    end ~ preserve

    B(θ_33, θ_1500): moisture_tension_curve_coeff_B => begin
        (log(1500) - log(33)) / (log(θ_33) - log(θ_1500))
    end ~ preserve

    # slope of logarithmic tension-moisture curve
    λ(B): pore_size_distribution => begin
        1 / B
    end ~ preserve

    K_s(θ_s, θ_33, λ): saturated_hydraulic_conductivity => begin
        1930(θ_s - θ_33)^(3-λ)
    end ~ preserve(u"mm/hr") # K_s,i (m day-1)

    K_at(K_s, θ_s, λ; θ): hydraulic_conductivity_at => begin
        #TODO: need bounds check?
        # θ = min(θ, θ_s)
        # (Ψ_at(vwc) < Ψ_ae) && (θ = θ_s)
        K_s * (θ / θ_s)^(3 + 2/λ)
    end ~ call(u"mm/hr") # K_theta,i (m day-1)

    # soil matric suction (moisture - tension)

    Ψ_et(S, C, OM, θ=θ_s_33): tension_air_entry_first => begin
        -21.674S - 27.932C - 81.975θ + 71.121S*θ + 8.294C*θ + 14.05S*C + 27.161
    end ~ preserve(u"kPa") # psi_et (kPa)

    Ψ_e(Ψ=nounit(Ψ_et, u"kPa")): tension_air_entry => begin
        Ψ_e = Ψ + (0.02Ψ^2 - 0.113Ψ - 0.70)
        #TODO: need bounds check?
        # max(Ψ_e, zero(Ψ_e))
    end ~ preserve(u"kPa") # psi_e (kPa)

    Ψ_at(θ_s, θ_33, θ_1500, Ψ_e, A, B; θ): tension_at => begin
        if θ_s <= θ
            Ψ_e
        elseif θ_33 <= θ
            33u"kPa" - (θ - θ_33) * (33u"kPa" - Ψ_e) / (θ_s - θ_33)
        elseif θ_1500 <= θ
            A*θ^-B
        else
            #@show "too low θ = $θ < θ_1500 = $θ_1500"
            A*θ^-B
        end
    end ~ call(u"kPa") # psi_theta (kPa)

    Hm_at(Ψ_at; θ): matric_head_at => begin
        Ψ_at(θ) * u"m" / 9.8041u"kPa"
    end ~ call(u"m") # H_mi (m)
end

#TODO: support convenient way to set up custom Clock
#TODO: support unit reference again?
@system SoilClock(Clock) <: Clock begin
    step => 15u"minute" ~ preserve(u"hr", parameter)
end
@system SoilContext(Context) <: Context begin
    context ~ ::Context(override)
    clock(config) ~ ::SoilClock
end

#TODO: implement LayeredTexture for customization
@system Layer(CharacteristicTransfer) begin
    context ~ ::SoilContext(override)

    i: index ~ ::int(override)
    θ_i: vwc_initial => 0.4 ~ preserve(extern)

    # Soil layer depth and cumulative thickness (2.4.2)
    z: depth ~ preserve(u"m", extern) # z_i (m)
    d_r: rooting_depth ~ track(u"m", override) # d_root (m)

    s: thickness ~ preserve(u"m", extern) # s_i (m)
    ss: cumulative_thickness ~ preserve(u"m", extern) # S_i (m)

    s_r(s, ss, d_r): root_zone_thickness => begin
        z = zero(d_r)
        max(s - max(ss - d_r, z), z)
    end ~ track(u"m") # s_i | s_i - (S_i - d_root) (m)

    𝚯_r(θ, s_r): water_content_root_zone => θ * s_r ~ track(u"m") # Theta_root,i (m) (Eq. 2.95)
    𝚯_r_wp(θ_wp, s_r): water_content_root_zone_wilting_point => θ_wp * s_r ~ track(u"m")
    𝚯_r_fc(θ_fc, s_r): water_content_root_zone_field_capacity => θ_fc * s_r ~ track(u"m")
    𝚯_r_sat(θ_sat, s_r): water_content_root_zone_saturation => θ_sat * s_r ~ track(u"m")

    # Root extraction of water (2.4.5)
    ϕ(z, d_r): water_extraction_ratio => begin
        cj = iszero(d_r) ? 0 : min(1, z / d_r)
        1.8cj - 0.8cj^2
    end ~ track # phi_i

    # Hydraulic conductivity (2.4.6)
    K(K_at, θ): hydraulic_conductivity => K_at(θ) ~ track(u"m/d") # k_i (m day-1)

    # Matric suction head (2.4.7)
    Hm(Hm_at, θ): matric_head => Hm_at(θ) ~ track(u"m") # H_mi (m)

    # Gravity head (2.4.8)
    Hg(z): gravity_head ~ preserve(u"m") # H_gi (m)

    # Total head
    H(Hm, Hg): total_head => Hm + Hg ~ track(u"m") # H_i (m)

    # Water content (2.4.10)
    qi: water_flux_in => 0 ~ track(u"m/d", ref, skip) # q_i (m day-1)
    qo: water_flux_out => 0 ~ track(u"m/d", ref, skip) # q_o (m day-1)
    q̂(qi, qo): water_flux_net => qi - qo ~ track(u"m/d") # q^hat_i (m day-1)
    𝚯(q̂): water_content ~ accumulate(init=𝚯_i, u"m") # Theta_i (m)

    𝚯_i(θ_i, s): water_content_initial => θ_i * s ~ preserve(u"m")
    𝚯_wp(θ_wp, s): water_content_wilting_point => θ_wp * s ~ track(u"m")
    𝚯_fc(θ_fc, s): water_content_field_capacity => θ_fc * s ~ track(u"m")
    𝚯_sat(θ_sat, s): water_content_saturation => θ_sat * s ~ track(u"m")

    # Volumetric water content (-)
    θ(i, 𝚯, 𝚯_wp, 𝚯_sat, s): volumetric_water_content => begin
        #FIXME: remove clamping?
        #HACK: clamping only for vwc
        # Teh uses 0.005 m3/m3 instead of wilting point
        #𝚯 = clamp(𝚯, 𝚯_wp, 𝚯_sat)
        θ = min(𝚯, 𝚯_sat) / s
        θ = max(θ, 0.005)
    end ~ track # Theta_v,i (m3 m-3)
end

@system SurfaceInterface begin
    context ~ ::SoilContext(override)
    l: layer ~ ::Layer(override)

    R: precipitation ~ track(u"m/d", override)
    Ea: evaporation_actual ~ track(u"m/d", override)
    Ta: transpiration_actual ~ track(u"m/d", override)

    Tai(Ta, ϕ=l.ϕ): water_extraction => begin
        Ta * ϕ
    end ~ track(u"m/d")

    q(R, Ea, Tai): flux => begin
        R - Ea - Tai
    end ~ track(u"m/d")

    _q(l, q) => begin
        Cropbox.setvar!(l, :qi, q)
    end ~ ::Nothing
end

@system SoilInterface begin
    context ~ ::SoilContext(override)
    l1: upper_layer ~ ::Layer(override)
    l2: lower_layer ~ ::Layer(override)

    Ta: transpiration_actual ~ track(u"m/d", override)

    K(K1=l1.K, K2=l2.K, s1=l1.s, s2=l2.s): hydraulic_conductivity => begin
        ((K1*s1) + (K2*s2)) / (s1 + s2)
    end ~ track(u"m/d") # k^bar (m day-1)

    # Hydraulic gradient (2.4.9)
    ΔH(H1=l1.H, H2=l2.H): hydraulic_gradient => begin
        H2 - H1
    end ~ track(u"m") # (m)

    Δz(z1=l1.z, z2=l2.z): depth_gradient => begin
        z2 - z1
    end ~ track(u"m") # (m)

    Tai(Ta, ϕ1=l1.ϕ, ϕ2=l2.ϕ): water_extraction => begin
        Ta * (ϕ2 - ϕ1)
    end ~ track(u"m/d")

    q(K, ΔH, Δz, Tai): flux => begin
        K * (ΔH / Δz) - Tai
    end ~ track(u"m/d") # q_i (m day-1)

    _q(l1, l2, q) => begin
        Cropbox.setvar!(l1, :qo, q)
        Cropbox.setvar!(l2, :qi, q)
    end ~ ::Nothing
end

@system BedrockInterface begin
    context ~ ::SoilContext(override)
    l: layer ~ ::Layer(override)

    q(l.K): flux ~ track(u"m/d")

    _q(l, q) => begin
        Cropbox.setvar!(l, :qo, q)
    end ~ ::Nothing
end

@system SoilWeather begin
    s: store ~ provide(index=:timestamp, init=1u"d", parameter)
    R: precipitation ~ drive(from=s, by=:precipitation, u"mm/d")
    T: transpiration ~ drive(from=s, by=:transpiration, u"mm/d")
    E: evaporation ~ drive(from=s, by=:evaporation, u"mm/d")
end

# w = instance(SoilWeather, config=(
#     :Clock => (:step => 24),
#     :SoilWeather => (:filename => "test/PyWaterBal.csv")
# ))

#FIXME: not just SoilClock, but entire Context should be customized for sub-timestep handling
#TODO: implement sub-timestep advance
# 2.4.11 Daily integration
# iterations=100
# Theta_i,t+1 (m day-1) (Eq. 2.105)
@system SoilModule begin
    context ~ ::SoilContext(override)
    w: weather ~ ::SoilWeather(override)

    d_r: rooting_depth ~ track(u"m", override) # d_root (m)

    # Partitioning of soil profile into several layers (2.4.1)
    L(context, d_r): layers => begin
        # Soil layer depth and cumulative thickness (2.4.2)
        n = 5
        s = 0.2u"m" # thickness
        ss = 0u"m" # cumulative_thickness
        θ = 0.4 # vwc_initial
        L = Layer[]
        for i in 1:n
            z = ss + s/2 # depth
            l = Layer(; context, i, θ_i=θ, z, d_r, s, ss)
            push!(L, l)
            ss += s
        end
        L
    end ~ ::Vector{Layer}

    surface_interface(context, layer=L[1], R, Ea, Ta) ~ ::SurfaceInterface

    soil_interfaces(context, L, Ta) => begin
        [SoilInterface(; context, l1=a, l2=b, Ta) for (a, b) in zip(L[1:end-1], L[2:end])]
    end ~ ::Vector{SoilInterface}

    bedrock_interface(context, layer=L[end]) ~ ::BedrockInterface

    interfaces(L, surface_interface, soil_interfaces, bedrock_interface) => begin
        [surface_interface, soil_interfaces..., bedrock_interface]
    end ~ ::Vector{System}(skip)

    # Actual evaporation (2.4.3)
    RD_e(θ=L[1].θ, θ_sat=L[1].θ_sat): evaporation_reduction_factor => begin
        1 / (1 + (3.6073 * (θ / θ_sat))^-9.3172)
    end ~ track # R_D,e
    Ep(w.E): evaporation_potential ~ track(u"m/d")
    Ea(Ep, RD_e): evaporation_actual => Ep * RD_e ~ track(u"m/d") # E_a (m day-1)

    # Actual transpiration (2.4.4)
    θ_r(L, d_r): volumetric_water_content_root_zone => begin
        sum(l.𝚯_r' for l in L) / d_r
    end ~ track # Theta_v,root (m3 m-3)

    θ_r_wp(L, d_r): volumetric_water_content_root_zone_wilting_point => begin
        sum(l.𝚯_r_wp' for l in L) / d_r
    end ~ track # (m3 m-3)

    θ_r_fc(L, d_r): volumetric_water_content_root_zone_field_capacity => begin
        sum(l.𝚯_r_fc' for l in L) / d_r
    end ~ track # (m3 m-3)

    θ_r_sat(L, d_r): volumetric_water_content_root_zone_saturation => begin
        sum(l.𝚯_r_sat' for l in L) / d_r
    end ~ track # (m3 m-3)

    RD_t(θ_r, θ_r_wp, θ_r_sat): transpiration_reduction_factor => begin
        θ_cr = (θ_r_wp + θ_r_sat) / 2
        (θ_r - θ_r_wp) / (θ_cr - θ_r_wp)
        #FIXME: 0 instead of 0.01?
    end ~ track(min=0.01, max=1) # R_D,t
    Tp(w.T): transpiration_potential ~ track(u"m/d")
    Ta(Tp, RD_t): transpiration_actual => Tp * RD_t ~ track(u"m/d") # T_a (m day-1)

    R(w.R): precipitation ~ track(u"m/d")
end

@system SoilController(Controller) begin
    w(context, config): weather ~ ::SoilWeather
    sc(context, config): soil_context ~ ::SoilContext(context)
    d_r: rooting_depth => 0.3 ~ track(u"m")
    s(context=soil_context, w, d_r): soil ~ ::SoilModule
end

end

@testset "soil" begin
    r = simulate(SoilWater.SoilController, stop=80u"d",
        config=(
            :Clock => (:step => 1u"d"),
            :SoilClock => (:step => 15u"minute"),
            :SoilWeather => (:store => "$(@__DIR__)/data/soil/PyWaterBal.csv"),
        ),
        target=(
            :v1 => "s.L[1].θ",
            :v2 => "s.L[2].θ",
            :v3 => "s.L[3].θ",
            :v4 => "s.L[4].θ",
            :v5 => "s.L[5].θ",
        ),
    )
    @test r.time[end] == 80u"d"
    visualize(r, :time, [:v1, :v2, :v3, :v4, :v5], ylim=(0.2, 0.45)) |> display
end


text:
using BinaryBuilder

# Collection of sources required to build Gettext
name = "Gettext"
version = v"0.21.0"

sources = [
    ArchiveSource("https://ftp.gnu.org/pub/gnu/gettext/gettext-$(version.major).$(version.minor).tar.xz",
                  "d20fcbb537e02dcf1383197ba05bd0734ef7bf5db06bdb241eb69b7d16b73192"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gettext-*/

export CFLAGS="-O2"
export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"

if [[ "${target}" == *-mingw* ]]; then
    # Apply patch from https://lists.gnu.org/archive/html/bug-gettext/2020-07/msg00035.html
    #      ../woe32dll/.libs/libgettextsrc_la-c++format.o: In function `__static_initialization_and_destruction_0':
    #      /workspace/srcdir/gettext-0.21/gettext-tools/src/../woe32dll/../src/format.c:67: undefined reference to `__imp_formatstring_ruby'
    atomic_patch -p1 ../patches/0001-build-Fix-build-failure-on-mingw-formatstring_ruby.patch
fi

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-static \
    --enable-relocatable \
    --with-libiconv-prefix=${prefix} \
    am_cv_lib_iconv=yes \
    am_cv_func_iconv=yes
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgettextlib", "libgettextlib-$(version.major)"], :libgettext)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("Libiconv_jll"),
    Dependency("XML2_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#!/usr/bin/env julia

using Luxor, Test, Random

Random.seed!(42)

function image_testing(fname)
    width, height = 4000, 4000
    margin = 500
    Drawing(width, height, fname)
    origin()
    background("grey25")
    setline(5)
    sethue("green")

    # try a nonpng
    @test_throws ErrorException readpng(dirname(dirname(pathof(Luxor))) * "/test/test-image.jl")

    image = readpng(dirname(dirname(pathof(Luxor))) * "/test/stackoverflow.png")
    w = image.width
    h = image.height
    pagetiles = Tiler(width, height, 7, 9)
    tw = pagetiles.tilewidth/2
    for (pos, n) in pagetiles
        circle(pos, tw, :stroke)
        circle(pos, tw, :clip)
        gsave()
        translate(pos)
        scale(25, 25)
        rotate(rand(0.0:pi/8:2pi))
        placeimage(image, O, centered=true)
        grestore()
        clipreset()
    end
    @test finish() == true
end

function placetesting(fname)
    # make svg file
    Drawing(200, 200, "svgimage.svg")
    origin()
    setline(40)
    sethue("rebeccapurple")
    squircle(O, 80, 80, :strokepreserve)
    sethue("white")
    fillpath()
    fontsize(40)
    sethue("black")
    text("SVG", halign=:center, valign=:middle)
    finish()

    # make png file
    Drawing(200, 200, "pngimage.png")
    origin()
    setline(40)
    sethue("red")
    squircle(O, 80, 80, :strokepreserve)
    sethue("white")
    fillpath()
    fontsize(40)
    sethue("black")
    text("PNG", halign=:center, valign=:middle)
    finish()

    Drawing(1200, 1200, fname)
    origin()
    svgimg = readsvg("svgimage.svg")
    pngimg = readpng("pngimage.png")
    for (pt, n) in Tiler(1200, 1200, 4, 4)
        @layer begin
            translate(pt)
            isodd(n) && placeimage(svgimg, O, centered=true)
            iseven(n) && placeimage(pngimg, O, centered=true)
        end
    end
    finish()
end

fname = "test-image.png"
image_testing(fname)

fname = "test-image-place.png"
placetesting(fname)

println("...finished test: output in $(fname)")


text:
abstract type Dist end

"""
    GetDistante(dist::Dist, p1::AbstractArray, p2::AbstractArray)

Returns the distance "dist" between array "p1" and array "p2".

# Examples
    GetDistance(Euclidean(), [2,2], [2,4])
    GetDistance(CityBlock(), [1, 2, 3], [4, 5, 6])
"""
function GetDistance(dist::Dist, p1::AbstractArray, p2::AbstractArray)
    if length(p1) != length(p2)
        throw(DimensionMismatch("Arrays with different dimensions, first array: $(length(p1)) != second array: $(length(p2))."))
    end
    if length(p1) == 0
        throw(DimensionMismatch("Can not calculate distance for empty arrays"))
    end

    return calculate(dist, p1, p2)
end

struct Euclidean <: Dist end
"""
    calculate(::Euclidean, p1::AbstractArray, p2::AbstractArray)

Return the euclidean distance between arrays "p1" and "p2".
"""
calculate(::Euclidean, p1::AbstractArray, p2::AbstractArray) = sqrt(sum((p1 - p2) .^ 2))

struct CityBlock <: Dist end
"""
    calculate(::CityBlock, p1::AbstractArray, p2::AbstractArray)

Return the city block distance between arrays "p1" and "p2".
"""
calculate(::CityBlock, p1::AbstractArray, p2::AbstractArray) = sum(abs.(p1 - p2))

struct TotalVariation <: Dist end
"""
    calculate(::TotalVariation, p1::AbstractArray, p2::AbstractArray)

Return the total variation distance between arrays "p1" and "p2".
"""
calculate(::TotalVariation, p1::AbstractArray, p2::AbstractArray) = sum(abs.(p1 - p2)) / 2

struct Chebyshev <: Dist end
"""
    calculate(::Chebyshev, p1::AbstractArray, p2::AbstractArray)

Return the Chebyshev distance between arrays "p1" and "p2".
"""
calculate(::Chebyshev, p1::AbstractArray, p2::AbstractArray) = maximum(abs.(p1-p2))

struct Jaccard <: Dist end
"""
    calculate(::Jaccard, p1::AbstractArray, p2::AbstractArray)

Return the Jaccard distance between arrays "p1" and "p2".
"""
calculate(::Jaccard, p1::AbstractArray, p2::AbstractArray) = 1 - sum(min(p1, p2)) / sum(max(p1,p2))

struct BrayCurtis <: Dist end
"""
    calculate(::BrayCurtis, p1::AbstractArray, p2::AbstractArray)

Return the Bray Curtis distance between arrays "p1" and "p2".
"""
calculate(::BrayCurtis, p1::AbstractArray, p2::AbstractArray) = sum(abs.(p1 - p2)) / sum(abs.(p1 + p2))

struct CosineDist <: Dist end
"""
    calculate(::CosineDist, p1::AbstractArray, p2::AbstractArray)

Return the cosine distance between arrays "p1" and "p2".
"""
calculate(::CosineDist, p1::AbstractArray, p2::AbstractArray) = 1 - sum(p1 .* p2) / ((sqrt(sum(p1 .^ 2))) * (sqrt(sum(p2 .^ 2))))

struct SpanNormDist <: Dist end
"""
    calculate(::SpanNormDist, p1::AbstractArray, p2::AbstractArray)

Return the span-norm distance between arrays "p1" and "p2".
"""
calculate(::SpanNormDist, p1::AbstractArray, p2::AbstractArray) = maximum(p1 - p2) - minimum(p1 - p2)


text:
function operar()
    println(10+10)
end
function operar(numeroUno::Int, numeroDos::Int)
    println(numeroUno+numeroDos)
end

operar()
operar(100,100)


text:
module Services
include("bucket.jl")
include("clibucket.jl")
include("awscliprovider.jl")
include("gcscliprovider.jl")

include("queue.jl")
include("awsqueue.jl")

include("cache.jl")
include("filesystemcache.jl")

include("datasource.jl")
include("bucketcachedatasource.jl")

include("daemon.jl")
end # module Services

max_forks_repo_name: seung-lab/SimpleTasks


text:
import NonlinearSolvers
const NS = NonlinearSolvers

export calculate_sfc_fluxes_energy
export DryBulkFormulaWithRadiation, DryBulkFormula, LinearRelaxation, DryMonin

"""
Options for calulating surface fluxes
"""
abstract type SurfaceFluxType end

struct LinearRelaxation <: SurfaceFluxType end
struct DryBulkFormula <: SurfaceFluxType end
struct DryBulkFormulaWithRadiation <: SurfaceFluxType end
struct DryMonin <: SurfaceFluxType end

SF.surface_conditions(
    args...;
    wθ_flux_star::Union{Nothing, FT} = nothing,
    universal_func::Union{Nothing, F} = SF.UniversalFunctions.Businger,
    sol_type::NS.SolutionType = NS.CompactSolution(),
    tol::NS.AbstractTolerance = NS.ResidualTolerance{FT}(sqrt(eps(FT))),
    maxiter::Int = 10_000,
) where {FT <: AbstractFloat, APS, F} =
    SF.surface_conditions(args..., wθ_flux_star, universal_func, sol_type, tol, maxiter)

calculate_sfc_fluxes_energy(formulation::LinearRelaxation, p, θ_sfc, θ_1) = p.λ .* (θ_sfc .- θ_1)

"""
    calculate_sfc_fluxes_energy(formulation::DryBulkFormula, parameters, θ_sfc, θ_1, u_1, v_1, ρ_1, t )

- calculates momentum and sensible heat fluxes using a bulk formula with constant drag/transfer coefficients
"""
calculate_sfc_fluxes_energy(formulation::DryBulkFormula, p, θ_sfc, θ_1, u_1, v_1, ρ_1) =
    p.Ch * p.C_p * ρ_1 * sqrt(u_1^2 + v_1^2) * (θ_sfc .- θ_1)

"""
    calculate_sfc_fluxes_energy(formulation::DryBulkFormulaWithRadiation, parameters, θ_sfc, θ_1, u_1, v_1, ρ_1, t )

- calculates momentum and sensible heat fluxes using a bulk formula with constant drag/transfer coefficients
- includes a simple representation of surface radiative fluxes with a temporal dependency
"""
function calculate_sfc_fluxes_energy(formulation::DryBulkFormulaWithRadiation, parameters, θ_sfc, θ_1, u_1, v_1, ρ_1, t)

    R_LW_incoming = 0.0
    p = parameters
    R_SW = (1 - p.α) * p.τ * p.F_sol * (1 .+ sin(t * 2π / p.τ_d))
    R_LW = p.ϵ * (p.σ * θ_sfc .^ 4) - R_LW_incoming
    SH = p.Ch * p.C_p * ρ_1 * sqrt(u_1^2 + v_1^2) * (θ_sfc - θ_1) # convert to theta!!
    #p_sfc = p.pₒ
    #LH   = p.lambda * p.g_w * (q_sat(θ_sfc, p_sfc) - q_a)

    F_tot = -(R_SW - R_LW - SH)#- LH
end

"""
    calculate_sfc_fluxes_energy(formulation::DryMonin, parameters, θ_sfc, θ_1, windspeed_1, ρ_1, t )

- calculates momentum and sensible heat fluxes using drag/transfer coefficients calculated in CliMA's `SurfaceFluxes.jl` module. This uses the Monin Obukhov theory and Nishizawa and Kitamura (2018) method for the finite volume discretization
"""
function calculate_sfc_fluxes(formulation::DryMonin, parameters, θ_sfc, θ_1, uv_1, ρ_1, t)

    domain_atm = Domains.IntervalDomain(0, parameters.atmos_Lz, x3boundary = (:bottom, :top))
    mesh_atm = Meshes.IntervalMesh(domain_atm, nelems = parameters.atmos_Nz)
    center_space_atm = Spaces.CenterFiniteDifferenceSpace(mesh_atm)
    z_centers = Fields.coordinate_field(center_space_atm)
    z_in = parent(z_centers)[1]

    windspeed_1 = LinearAlgebra.norm(uv_1)
    p = parameters

    z_0m = 1e-5 # eventually DryMonin.z_0m? (use cases: where atmos will get it with prescribed land, or vice versa, or coupled)
    z_0c = z_0m # z_0c = function(z0m, ustar...) eventually
    z_0 = [z_0m, z_0c]
    x_in = [windspeed_1, θ_1]

    x_s = [windspeed_1 * 0.0, θ_sfc] # we are assuming θ_sfc = T(-\Delta z).
    θ_basic = deepcopy(θ_1) # for buoyancy calculation

    # Initial guesses for MO parameters, these should be a function of state.
    LMO_init = z_in * 100 # Initial value so that ξ_init<<1
    u_star_init = 0.1 * windspeed_1
    th_star_init = 1.0 * θ_1
    MO_param_guess = [LMO_init, u_star_init, th_star_init] # guess based on current atmos state

    output = SF.surface_conditions(
        CLIMAparam_set,
        MO_param_guess,
        x_in,
        x_s,
        z_0,
        θ_basic,
        z_in,
        SF.FVScheme(),
        maxiter = 500, # often unconverged + coeffs v small, check with FMS (esp Mo_params_guess), update for conditional states
    )

    C_exchange = output.C_exchange
    Cd = C_exchange[1] # Cd
    Ch = C_exchange[2] # Ch

    # upward sensible potential temperature flux
    SH = Ch * ρ_1 * windspeed_1 * (θ_1 - θ_sfc)

    # wind stress
    τ = Geometry.Cartesian3Vector(Cd * windspeed_1) ⊗ uv_1

    fluxes = (τ = τ, SH = SH)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Test
using RadioHead

c = Coax(15.0, 0.8)

@testset "units" begin

    m = Meter(2.0)
    @test toFeet(m).value == 3.2808f0 * 2
    @test toMeters(m).value == 1.0f0 * 2

    m = Foot(3.0)
    @test toMeters(m).value == 0.3048f0 * 3
    @test toFeet(m).value == 1.0f0 * 3

end

@testset "coax" begin

    @test c.centerFreq === 15.0f0
    @test c.velocityFactor === 0.8f0
    @test c.wavelengthUnit === 16.0

end


text:
export
Bounds2, Bounds3,
Bounds2i, Bounds2f, Bounds2d,
Bounds3i, Bounds3f, Bounds3d

struct Bounds{N, T<:Number, U<:Point{N, T}} <:FieldVector{2, T}
    p_min::U
    p_max::U

    Bounds{N,T,U}(p_min, p_max) where {N,T,U} = new{N,T,U}(min(p_min, p_max), max(p_min, p_max))

    Bounds{N,T,U}(p) where {N,T,U} = new{N,T,U}(p, p)

    Bounds{N,T,U}() where {N,T,U} = new{N,T,U}(typeof(U)(typemax(T)), typeof(U)(typemin(T)))
end

Bounds2{T<:Number} = Bounds{2, T, Point2{T}}
Bounds3{T<:Number} = Bounds{3, T, Point3{T}}

Bounds(p_min::Point{N,T}, p_max::Point{N,U}) where {N,T,U} = Bounds{N,promote_type(T,U)}(p_min, p_max)
Bounds2(p_min::Point2{T}, p_max::Point2{U}) where {T,U} = Bounds2{promote_type(T,U)}(p_min, p_max)
Bounds3(p_min::Point3{T}, p_max::Point3{U}) where {T,U} = Bounds3{promote_type(T,U)}(p_min, p_max)

Base.promote_rule(::Type{Bounds{N,T}}, ::Type{Bounds{N,U}}) where {N,T,U} = Bounds{N, promote_type(T,U)}
Base.promote_rule(::Type{Bounds{N,T}}, ::Type{Point{N,U}}) where {N,T,U} = Bounds{N, promote_type(T,U)}

Bounds2i = Bounds2{Int}
Bounds2s = Bounds2{Float32}
Bounds2d = Bounds2{Float64}

Bounds3i = Bounds3{Int}
Bounds3s = Bounds3{Float32}
Bounds3d = Bounds3{Float64}


text: funcs = [ () -> i^2 for i = 1:10 ]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PRS
using Test

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


text:
Nx,Ny,Nz = 10,20,11
G = ones(Bool,Nx+2,Ny+2,Ny+2)
G[1,:,:] = G[end,:,:] = false
G[:,1,:] = G[:,end,:] = false
G[:,:,1] = G[:,:,end] = false

I = AcFDTD.getNodeProperty(G)
@test all(I[2:end-1,2:end-1,2:end-1] .== 0x01)
#walls
@test all(I[1  ,2:end-1,2:end-1].==0x02)
@test all(I[end,2:end-1,2:end-1].==0x03)
@test all(I[2:end-1,  1,2:end-1].==0x04)
@test all(I[2:end-1,end,2:end-1].==0x05)
@test all(I[2:end-1,2:end-1,  1].==0x06)
@test all(I[2:end-1,2:end-1,end].==0x07)
#edges
@test all(I[  1,  1,2:end-1].==0x08)
@test all(I[  1,end,2:end-1].==0x09)
@test all(I[  1,2:end-1,  1].==0x0a)
@test all(I[  1,2:end-1,end].==0x0b)

@test all(I[end,  1,2:end-1].==0x0c)
@test all(I[end,end,2:end-1].==0x0d)
@test all(I[end,2:end-1,  1].==0x0e)
@test all(I[end,2:end-1,end].==0x0f)

@test all(I[2:end-1,   1,  1].==0x10)
@test all(I[2:end-1,   1,end].==0x11)
@test all(I[2:end-1, end,  1].==0x12)
@test all(I[2:end-1, end,end].==0x13)

#corners
@test I[  1,  1,  1] == 0x14
@test I[  1,  1,end] == 0x15
@test I[  1,end,  1] == 0x16
@test I[  1,end,end] == 0x17

@test I[end,  1,  1] == 0x18
@test I[end,  1,end] == 0x19
@test I[end,end,  1] == 0x1a
@test I[end,end,end] == 0x1b

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
export traits, Trait

"Returns the traits of an objec"
function traits end

"Meta data to attach to ω::Ω"
const Tags = NamedTuple

struct Trait{T} end

symtotrait(x::Symbol) = Trait{x}

@generated function traits(k::Tags{K, V}) where {K, V}
  traits_ = map(symtotrait, K)
  AndTrait{Union{traits_...}}()
end

@generated function traits(k::Type{Tags{K, V}}) where {K, V}
  traits_ = map(symtotrait, K)
  AndTrait{Union{traits_...}}()
end


text:
using MagneticReadHead
using Test
@testset "breadcrumbs" begin
    iob = IOBuffer()
    filename = joinpath(@__DIR__, "demo.jl")
    MagneticReadHead.breadcrumbs(iob, filename, 3)
    @test String(take!(iob)) ==
        " function eg1()\n     z = eg2(2)\n➧    eg_last(z)\n end\n \n"
end


text:
module SymSemiseparableMatrices

# Importing Relevant Packages
using LinearAlgebra

# Must be imported here, as it is relevant before "syntax.jl"
import LinearAlgebra: inv!, tr, mul!, logdet, ldiv!, cholesky, det
import Base: inv, size, eltype, adjoint, *, \, size, eltype, Matrix, getindex, getproperty

# Creating Types
struct SymSemiseparableMatrix{T,UU<:AbstractMatrix{T},VV<:AbstractMatrix{T}} <: AbstractMatrix{T}
    n::Int64
    p::Int64
    Ut::UU
    Vt::VV
end
struct SymSemiseparableCholesky{T,UU<:AbstractMatrix{T},WW<:AbstractMatrix{T}} <: AbstractMatrix{T}
    n::Int64
    p::Int64
    Ut::UU
    Wt::WW
end
struct DiaSymSemiseparableMatrix{T,UU<:AbstractMatrix{T},VV<:AbstractMatrix{T},dd<:AbstractVector{T}} <: AbstractMatrix{T}
    n::Int64
    p::Int64
    Ut::UU
    Vt::VV
    d::dd
end
struct DiaSymSemiseparableCholesky{T,UU<:AbstractMatrix{T},WW<:AbstractMatrix{T},dd<:AbstractVector{T}} <: AbstractMatrix{T}
    n::Int64
    p::Int64
    Ut::UU
    Wt::WW
    d::dd
end

# Utility functions
include("utility_functions.jl")

# Matrices
include("SymSemiseparableMatrix.jl")
include("SymSemiseparableCholesky.jl")
include("DiaSymSemiseparableMatrix.jl")
include("DiaSymSemiseparableCholesky.jl")

# Spline Kernel
include("spline_kernel.jl")

# Exporting Relevant
export SymSemiseparableMatrix, SymSemiseparableCholesky
export DiaSymSemiseparableMatrix, DiaSymSemiseparableCholesky
export trinv

end # module


text:
include("JuliaDe3DFitting.jl")

X = [1.2 1 1; 2 2.3 2; 3 3 3.1; 4 4.1 4; 5 5 5]
Y = planeFit(X)
@show X
@show Y

X = [1 0 0; 0 1 0 ; 0 0 1; -1 0 0; 0 -1 0; 0 0 -1]
Y = sphereFit(X, 100)
@show X
@show Y

X = [1 0 0; 0 1 0 ; -1 0 0; 0 -1 0; 1 0 1; 0 1 1; -1 0 1; 0 -1 1; 1 0 2; 0 1 2; -1 0 2; 0 -1 2]
Y = cylinderFit(X)
@show X
@show Y


text:
@testset "53.maximum-subarray.jl" begin
    @test maximum_subarray([-2, 1, -3, 4, -1, 2, 1, -5, 4]) == 6
    @test maximum_subarray([-1]) == -1
    @test maximum_subarray([5, 4, -1, 7, 8]) == 23
    @test maximum_subarray([-2, 1]) == 1
    @test maximum_subarray([
        8, -2, -4, -1, -8, 3, 8, 8, 3, 4, 2, -9, -1, -3, -6, 8, -3, 9
    ]) == 28
end


text:
module FastPolys

using StaticArrays

export fpoly,
       fpoly_d,
       fpoly_ed

push_str!(ex::Vector{Expr}, s::String) = push!(ex, Meta.parse(s))
append_str!(ex::Vector{Expr}, s::Vector{String}) = append!(ex, Meta.parse.(s))

# ACCUMULATOR FUNCTIONS
# ----------------------

# 1st order

_mon_muladd_(I::NTuple{1, Int}) =
   "m += x1_$(I[1])"

_mon_muladd_d_(I::NTuple{1, Int}) =
   ["dm_$(I[1]) += dx1_$(I[1])",]

_mon_muladd_ed_(I::NTuple{1, Int}) =
   ["m += x1_$(I[1])","dm_$(I[1]) += dx1_$(I[1])",]

# 2nd order

_mon_muladd_(I::NTuple{2, Int}) =
   "m = muladd(x1_$(I[1]), x2_$(I[2]), m)"

_mon_muladd_d_(I::NTuple{2, Int}) =
   ["dm_$(I[1]) = muladd(dx1_$(I[1]),  x2_$(I[2]), dm_$(I[1]))",
    "dm_$(I[2]) = muladd( x1_$(I[1]), dx2_$(I[2]), dm_$(I[2]))"]

_mon_muladd_ed_(I::NTuple{2, Int}) =
   ["m = muladd(x1_$(I[1]), x2_$(I[2]), m)",
    "dm_$(I[1]) = muladd(dx1_$(I[1]),  x2_$(I[2]), dm_$(I[1]))",
    "dm_$(I[2]) = muladd( x1_$(I[1]), dx2_$(I[2]), dm_$(I[2]))"]

# 3rd order

_mon_muladd_(I::NTuple{3, Int}) =
   "m = muladd(x1_$(I[1]), x2_$(I[2]) * x3_$(I[3]), m)"

_mon_muladd_d_(I::NTuple{3, Int}) =
   ["dm_$(I[1]) = muladd(dx1_$(I[1]),  x2_$(I[2]) *  x3_$(I[3]), dm_$(I[1]))",
    "dm_$(I[2]) = muladd( x1_$(I[1]), dx2_$(I[2]) *  x3_$(I[3]), dm_$(I[2]))",
    "dm_$(I[3]) = muladd( x1_$(I[1]),  x2_$(I[2]) * dx3_$(I[3]), dm_$(I[3]))"]

 _mon_muladd_ed_(I::NTuple{3, Int}) =
    ["a23 = x2_$(I[2]) * x3_$(I[3])",
     "m = muladd(x1_$(I[1]), a23, m)",
     "dm_$(I[1]) = muladd(dx1_$(I[1]),  a23, dm_$(I[1]))",
     "dm_$(I[2]) = muladd( x1_$(I[1]), dx2_$(I[2]) *  x3_$(I[3]), dm_$(I[2]))",
     "dm_$(I[3]) = muladd( x1_$(I[1]),  x2_$(I[2]) * dx3_$(I[3]), dm_$(I[3]))"]

# 4th order

_mon_muladd_(I::NTuple{4, Int}) =
   "m = muladd(x1_$(I[1]) * x2_$(I[2]), x3_$(I[3]) * x4_$(I[4]), m)"

_mon_muladd_d_(I::NTuple{4, Int}) =
   ["a12 = x1_$(I[1]) * x2_$(I[2])",
    "a34 = x3_$(I[3]) * x4_$(I[4])",
    "dm_$(I[1]) = muladd(dx1_$(I[1]) *  x2_$(I[2]), a34, dm_$(I[1]))",
    "dm_$(I[2]) = muladd( x1_$(I[1]) * dx2_$(I[2]), a34, dm_$(I[2]))",
    "dm_$(I[3]) = muladd( a12, dx3_$(I[3]) *  x4_$(I[4]), dm_$(I[3]))",
    "dm_$(I[4]) = muladd( a12,  x3_$(I[3]) * dx4_$(I[4]), dm_$(I[4]))"]

_mon_muladd_ed_(I::NTuple{4, Int}) =
   ["a12 = x1_$(I[1]) * x2_$(I[2])",
   "a34 = x3_$(I[3]) * x4_$(I[4])",
   "m = muladd(a12, a34, m)",
   "dm_$(I[1]) = muladd(dx1_$(I[1]) *  x2_$(I[2]), a34, dm_$(I[1]))",
   "dm_$(I[2]) = muladd( x1_$(I[1]) * dx2_$(I[2]), a34, dm_$(I[2]))",
   "dm_$(I[3]) = muladd( a12, dx3_$(I[3]) *  x4_$(I[4]), dm_$(I[3]))",
   "dm_$(I[4]) = muladd( a12,  x3_$(I[3]) * dx4_$(I[4]), dm_$(I[4]))"]

# 5th order

_mon_muladd_(I::NTuple{5, Int}) =
   "m = muladd(x1_$(I[1]) * x2_$(I[2]), x3_$(I[3]) * x4_$(I[4]) * x5_$(I[5]), m)"

_mon_muladd_d_(I::NTuple{5, Int}) =
   ["a12 = x1_$(I[1]) * x2_$(I[2])",
    "a123 = a12 * x3_$(I[3])",
    "a45 = x4_$(I[4]) * x5_$(I[5])",
    "a345 = x3_$(I[3]) * a45",
    "dm_$(I[1]) = muladd(dx1_$(I[1]) *  x2_$(I[2]), a345, dm_$(I[1]))",
    "dm_$(I[2]) = muladd( x1_$(I[1]) * dx2_$(I[2]), a345, dm_$(I[2]))",
    "dm_$(I[3]) = muladd( a12, dx3_$(I[3]) *  a45, dm_$(I[3]))",
    "dm_$(I[4]) = muladd( a123,  dx4_$(I[4]) *  x5_$(I[5]), dm_$(I[4]))",
    "dm_$(I[5]) = muladd( a123,   x4_$(I[4]) * dx5_$(I[5]), dm_$(I[5]))" ]

_mon_muladd_ed_(I::NTuple{5, Int}) =
   ["a12 = x1_$(I[1]) * x2_$(I[2])",
   "a123 = a12 * x3_$(I[3])",
   "a45 = x4_$(I[4]) * x5_$(I[5])",
   "a345 = x3_$(I[3]) * a45",
   "m = muladd(a12, a345, m)",
   "dm_$(I[1]) = muladd(dx1_$(I[1]) *  x2_$(I[2]), a345, dm_$(I[1]))",
   "dm_$(I[2]) = muladd( x1_$(I[1]) * dx2_$(I[2]), a345, dm_$(I[2]))",
   "dm_$(I[3]) = muladd( a12, dx3_$(I[3]) *  a45, dm_$(I[3]))",
   "dm_$(I[4]) = muladd( a123,  dx4_$(I[4]) *  x5_$(I[5]), dm_$(I[4]))",
   "dm_$(I[5]) = muladd( a123,   x4_$(I[4]) * dx5_$(I[5]), dm_$(I[5]))" ]

function _check_input_(D, A)
   @assert length(A) == D
   @assert typeof(A[1]) <: NTuple
   @assert eltype(A[1]) <: Integer
   for n = 2:length(A)
      @assert typeof(A[1]) == typeof(A[n])
   end
end

"""
`function fpoly(X, ::Val{A})`

TODO: write documentation
"""
@generated function fpoly(x::NTuple{D, SVector{N, T}},
                          ::Val{A}) where {D, N, T, A}
   # make some assumptions about what the input is
   _check_input_(D, A)

   # generate the code
   # -----------------
   exprs = Expr[]
   # read all the data from the StaticArrays
   #     for very small monomials, this could be the bottleneck, but for
   #     reasonably large ones it should make no difference
   for d = 1:D, n = 1:N
      push_str!(exprs, "x$(d)_$(n) = x[$(d)][$(n)]")
   end
   # also initialise the output
   push_str!(exprs, "m = zero($T)")
   # evaluate the monomial
   for I in zip(A...)
      push_str!(exprs, _mon_muladd_(I))
   end
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, exprs...))
      return m
   end
end

"""
`function fpoly(X, dX, ::Val{A})`

TODO: write documentation
"""
@generated function fpoly_d(
         x::NTuple{D, SVector{N, T}},
         dx::NTuple{D, SVector{N, T}},
         ::Val{A}
      ) where {D, N, T, A}
   # make some assumptions about what the input is
   _check_input_(D, A)

   # generate the code
   # -----------------
   exprs = Expr[]
   # read all the data from the StaticArrays
   #     for very small monomials, this could be the bottleneck, but for
   #     reasonably large ones it should make no difference
   for d = 1:D, n = 1:N
      push_str!(exprs, "x$(d)_$(n) = x[$(d)][$(n)]")
      push_str!(exprs, "dx$(d)_$(n) = dx[$(d)][$(n)]")
   end
   # initialise the output
   push_str!(exprs, "m = zero($T)")
   for n = 1:N
      push_str!(exprs, "dm_$(n) = zero($T)")
   end
   # evaluate the monomial and its derivative
   for I in zip(A...)
      append_str!(exprs, _mon_muladd_d_(I))
   end
   # collect the dm variables into an SVector
   coll = "SVector(" * prod("dm_$n, " for n = 1:N) * ")"
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, exprs...))
      $(Meta.parse(coll))
   end
end

@generated function fpoly_ed(
         x::NTuple{D, SVector{N, T}},
         dx::NTuple{D, SVector{N, T}},
         ::Val{A}
      ) where {D, N, T, A}
   # make some assumptions about what the input is
   _check_input_(D, A)

   # generate the code
   # -----------------
   exprs = Expr[]
   # read all the data from the StaticArrays
   #     for very small monomials, this could be the bottleneck, but for
   #     reasonably large ones it should make no difference
   for d = 1:D, n = 1:N
      push_str!(exprs, "x$(d)_$(n) = x[$(d)][$(n)]")
      push_str!(exprs, "dx$(d)_$(n) = dx[$(d)][$(n)]")
   end
   # initialise the output
   push_str!(exprs, "m = zero($T)")
   for n = 1:N
      push_str!(exprs, "dm_$(n) = zero($T)")
   end
   # evaluate the monomial and its derivative
   for I in zip(A...)
      append_str!(exprs, _mon_muladd_ed_(I))
   end
   # # collect the m and dm variables into an SVector
   # coll = "SVector( m ) ,
   # SVector(" * prod("dm_$n, " for n = 1:N) * ")"
   coll = "SVector(" * prod("dm_$n, " for n = 1:N) * ")"
   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, exprs...))
      return m, $(Meta.parse(coll))
   end
end

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
@recipe function f(m::Walk1D, h::SimHistory)
    xlim := (-1, m.p.t_max)
    ylim := (-m.p.threshx-1, m.p.threshx+1)
    @series begin
        seriestype := :line
        color := :black
        label := "thresh+"
        x = [0.0, m.p.t_max]
        y = [m.p.threshx, m.p.threshx]
        x, y
    end
    @series begin
        seriestype := :line
        color := :black
        label := "thresh-"
        x = [0.0, m.p.t_max]
        y = [-m.p.threshx, -m.p.threshx]
        x, y
    end
    @series begin
        label := "path"
        color := :blue
        x = [s.t for s in state_hist(h)[1:end-1]]
        y = [s.x for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        seriestype := :scatter
        label := "current position"
        t = state_hist(h)[end-1].t
        x = state_hist(h)[end-1].x
        if abs(x) < abs(m.p.threshx)
            color := :blue
        else
            color := :red
        end
        [t], [x]
    end
end