text:
# TODO change operations to avoid nonlinear expressions obtained from operations done on measures

## GeneralVariableRef--GeneralVariableRef extensions for different types
# var1 + var2
function Base.:+(lhs::V, rhs::W)::JuMP.GenericAffExpr where {V <: GeneralVariableRef,
                                                             W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericAffExpr(0.0, Pair{type, Float64}(lhs, 1.0),
                               Pair{type, Float64}(rhs, 1.0))
end

# var1 - var2
function Base.:-(lhs::V, rhs::W)::JuMP.GenericAffExpr where {V <: GeneralVariableRef,
                                                             W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    if lhs == rhs
        return zero(JuMP.GenericAffExpr{Float64, type})
    else
        return JuMP.GenericAffExpr(0.0, Pair{type, Float64}(lhs, 1.0),
                                   Pair{type, Float64}(rhs, -1.0))
    end
end

# var1 * var2
function Base.:*(lhs::V, rhs::W)::JuMP.GenericQuadExpr where {V <: GeneralVariableRef,
                                                              W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(JuMP.GenericAffExpr{Float64, type}(),
                                JuMP.UnorderedPair{type}(lhs, rhs) => 1.0)
end

# GeneralVariableRef--GenericAffExpr
# var + aff
function Base.:+(lhs::V,rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericAffExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    # For the variables to have the proper order in the result,
    # we need to add the lhs first.
    type = _var_type_parser(V, W)
    result = zero(JuMP.GenericAffExpr{C, type})
    result.constant = rhs.constant
    sizehint!(result, length(JuMP.linear_terms(rhs)) + 1)
    JuMP.add_to_expression!(result, one(C), lhs)
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, coef, var)
    end
    return result
end

# var - aff
function Base.:-(lhs::V, rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericAffExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    # For the variables to have the proper order in the result,
    # we need to add the lhs first.
    type = _var_type_parser(V, W)
    result = zero(JuMP.GenericAffExpr{C, type})
    result.constant = -rhs.constant
    sizehint!(result, length(JuMP.linear_terms(rhs)) + 1)
    JuMP.add_to_expression!(result, one(C), lhs)
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, -coef, var)
    end
    return result
end

# var * aff
function Base.:*(lhs::V, rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    if !iszero(rhs.constant)
        result = JuMP.GenericQuadExpr{C, type}(convert(JuMP.GenericAffExpr{C, type},
                                                       lhs * rhs.constant))
    else
        result = zero(JuMP.GenericQuadExpr{C, type})
    end
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, coef, lhs, var)
    end
    return result
end

# GeneralVariableRef--GenericQuadExpr
# var + quad
function Base.:+(v::V, q::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    new_q = convert(JuMP.GenericQuadExpr{C, type}, q)
    return JuMP.GenericQuadExpr(v + new_q.aff, copy(new_q.terms))
 end

 # var - quad
 function Base.:-(v::V, q::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                  V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = convert(JuMP.GenericQuadExpr{C, type}, -q)
    # This makes an unnecessary copy of aff, but it's important for v to appear
    # first.
    result.aff = v + result.aff
    return result
end

## Extend operators for GenericAffExpr--GeneralVariableRef of different types
# aff + var
function Base.:+(lhs::JuMP.GenericAffExpr{C, V}, rhs::W)::JuMP.GenericAffExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.add_to_expression!(JuMP.GenericAffExpr{C, type}(lhs.constant,
                                   lhs.terms), one(C), rhs)
end

# aff - var
function Base.:-(lhs::JuMP.GenericAffExpr{C, V}, rhs::W)::JuMP.GenericAffExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.add_to_expression!(JuMP.GenericAffExpr{C, type}(lhs.constant,
                                   lhs.terms), -one(C), rhs)
end

# aff * var
function Base.:*(lhs::JuMP.GenericAffExpr{C, V}, rhs::W)::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    if !iszero(lhs.constant)
        result = JuMP.GenericQuadExpr{C, type}(convert(JuMP.GenericAffExpr{C,
                                                     type}, lhs.constant * rhs))
    else
        result = zero(JuMP.GenericQuadExpr{C, type})
    end
    for (coef, var) in JuMP.linear_terms(lhs)
        JuMP.add_to_expression!(result, coef, var, rhs)
    end
    return result
end

## Extend operators for AffExpr--AffExpr for different types
# aff + aff
function Base.:+(lhs::JuMP.GenericAffExpr{C, V},
                 rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericAffExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    if length(JuMP.linear_terms(lhs)) > 50 || length(JuMP.linear_terms(rhs)) > 50
        if length(JuMP.linear_terms(lhs)) > 1
            JuMP.operator_warn(JuMP.owner_model(first(JuMP.linear_terms(lhs))[2]))
        end
    end
    type = _var_type_parser(V, W)
    return JuMP.add_to_expression!(copy(convert(JuMP.GenericAffExpr{C, type}, lhs)),
                                   convert(JuMP.GenericAffExpr{C, type}, rhs))
end

# aff - aff
function Base.:-(lhs::JuMP.GenericAffExpr{C, V},
                 rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericAffExpr  where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = copy(convert(JuMP.GenericAffExpr{C, type}, lhs))
    result.constant -= rhs.constant
    sizehint!(result, length(JuMP.linear_terms(lhs)) +
                      length(JuMP.linear_terms(rhs)))
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, -coef, var)
    end
    return result
end

# aff * aff
function Base.:*(lhs::JuMP.GenericAffExpr{C, V},
                 rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = zero(JuMP.GenericQuadExpr{C, type})
    JuMP.add_to_expression!(result, convert(JuMP.GenericAffExpr{C, type}, lhs),
                            convert(JuMP.GenericAffExpr{C, type}, rhs))
    return result
end

## Extend for operators on GenericAffExpr--GenericQuadExpr of different types
# aff + quad
function Base.:+(lhs::JuMP.GenericAffExpr{C, V},
                 rhs::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(lhs + rhs.aff,
                                copy(convert(JuMP.GenericQuadExpr{C, type}, rhs).terms))
end

# aff - quad
function Base.:-(lhs::JuMP.GenericAffExpr{C, V},
                 rhs::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = -convert(JuMP.GenericQuadExpr{C, type}, rhs)
    # This makes an unnecessary copy of aff, but it's important for a to appear
    # first.
    result.aff = lhs + result.aff
    return result
end

## Extend operators for GenericQuadExpr--GeneralVariableRef with different types
# quad + var
function Base.:+(lhs::JuMP.GenericQuadExpr{C, V}, rhs::W)::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(lhs.aff + rhs,
                                copy(convert(JuMP.GenericQuadExpr{C, type}, lhs).terms))
end

# quad - var
function Base.:-(lhs::JuMP.GenericQuadExpr{C, V}, rhs::W)::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(lhs.aff - rhs,
                                copy(convert(JuMP.GenericQuadExpr{C, type}, lhs).terms))
end

# Extend operators for GenericQuadExpr--GenericAffExpr with different types
# quad + aff
function Base.:+(lhs::JuMP.GenericQuadExpr{C, V},
                 rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(lhs.aff + rhs,
                                copy(convert(JuMP.GenericQuadExpr{C, type}, lhs).terms))
end

# quad - aff
function Base.:-(lhs::JuMP.GenericQuadExpr{C, V},
                 rhs::JuMP.GenericAffExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    return JuMP.GenericQuadExpr(lhs.aff - rhs,
                                copy(convert(JuMP.GenericQuadExpr{C, type}, lhs).terms))
end

## Extend operators for GenericQuadExpr--GenericQuadExpr with differnt types
# quad + quad
function Base.:+(lhs::JuMP.GenericQuadExpr{C, V},
                 rhs::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = copy(convert(JuMP.GenericQuadExpr{C, type}, lhs))
    for (coef, var1, var2) in JuMP.quad_terms(rhs)
        JuMP.add_to_expression!(result, coef, var1, var2)
    end
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, coef, var)
    end
    result.aff.constant += rhs.aff.constant
    return result
end

# quad - quad
function Base.:-(lhs::JuMP.GenericQuadExpr{C, V},
                 rhs::JuMP.GenericQuadExpr{C, W})::JuMP.GenericQuadExpr where {C,
                 V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    result = copy(convert(JuMP.GenericQuadExpr{C, type}, lhs))
    for (coef, var1, var2) in JuMP.quad_terms(rhs)
        JuMP.add_to_expression!(result, -coef, var1, var2)
    end
    for (coef, var) in JuMP.linear_terms(rhs)
        JuMP.add_to_expression!(result, -coef, var)
    end
    result.aff.constant -= rhs.aff.constant
    return result
end

ext: jl


text:
module TestExampleDownload

using CSV
using Test
using WorldBankData
using DataFrames

# reset the country and indicator caches from search_wdi.jl tests
# otherwise downloads will fail since the data is not in the cache
WorldBankData.reset_country_cache()
WorldBankData.reset_indicator_cache()

refdf = CSV.read(joinpath(dirname(@__FILE__), "example_data.csv"))

# the data gets frequently updated on the World Bank site use this to update the example_data.csv file
function update_example_data()
    dfnref = wdi(["NY.GNP.PCAP.CD", "AG.LND.ARBL.HA.PC"], ["US", "BR"], 1980, 2008, extra=true, verbose=true)
    CSV.write(joinpath(dirname(@__FILE__), "example_data.csv"), dfnref)
end

# download example case from documentation and compare to csv file
# retry 5 times if no data
function try_download(cntr=5)
    dfweb = ""
    while cntr > 0
        dfweb = wdi(["NY.GNP.PCAP.CD", "AG.LND.ARBL.HA.PC"], ["US", "BR"];
                    startyear = 1980, endyear = 2008, extra=true, verbose=true)
        if size(dfweb)[1] != 0
            break
        end
        println("      retry download")
        cntr = cntr - 1
        if cntr == 0
            error("unable to download data")
        end
    end
    dfweb
end

# data contains NA which breaks == check
function rm_na(df)
    df[map(ismissing, df[!, :NY_GNP_PCAP_CD]), :NY_GNP_PCAP_CD] = -123456
    df[map(ismissing, df[!, :AG_LND_ARBL_HA_PC]), :AG_LND_ARBL_HA_PC] = -123456
end

dfweb = try_download()

@test dfweb[!, :year] == refdf[!, :year]

rm_na(dfweb)
rm_na(refdf)
@test dfweb == refdf

end


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rastername, rasterpath, zipurl, zipname, zippath

@testset "WorldClim Weather" begin
    raster_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Weather", "prec", "wc2.1_2.5m_prec_2001-01.tif")
    @test rasterpath(WorldClim{Weather}, :prec; date=Date(2001, 1)) == raster_file
    @test rastername(WorldClim{Weather}, :prec; date=Date(2001, 1)) == "wc2.1_2.5m_prec_2001-01.tif"

    zip_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Weather", "zips", "wc2.1_2.5m_prec_2010-2018.zip")
    @test zippath(WorldClim{Weather}, :prec; decade=Date(2010)) == zip_file
    @test zipurl(WorldClim{Weather}, :prec; decade=Date(2010)) ==
        URI(scheme="https", host="biogeo.ucdavis.edu", path="/data/worldclim/v2.1/hist/wc2.1_2.5m_prec_2010-2018.zip")
    @test zipname(WorldClim{Weather}, :prec; decade=Date(2010)) ==
        "wc2.1_2.5m_prec_2010-2018.zip"

    # These files are 3GB each. Too big to test in CI.
    if !haskey(ENV, "CI")
        @test getraster(WorldClim{Weather}, :prec; date=Date(2001, 1)) == raster_file
        @test getraster(WorldClim{Weather}, (:prec,); date=Date(2001, 1)) == (prec=raster_file,)
        raster_files = [joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Weather", "prec", "wc2.1_2.5m_prec_2001-$(lpad(string(x), 2, '0')).tif") for x in 1:12]
        @test getraster(WorldClim{Weather}, :prec; date=Date(2001):Month(1):Date(2001, 12)) == raster_files
        @test getraster(WorldClim{Weather}, (:prec,); date=Date(2001):Month(1):Date(2001, 12)) == map(f -> (prec=f,), raster_files)
        @test getraster(WorldClim{Weather}, [:prec]; date=Date(2001):Month(1):Date(2001, 12)) == map(f -> (prec=f,), raster_files)
    end
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
function test_cross_validation()
    println("-- Testing cross validation")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    # 5-fold, top-4 recommendation with MF(data, 2)
    @test cross_validation(5, Recall, 4, MF, data, 2) <= 0.5

    # MostPopular(data)
    @test cross_validation(5, Recall, 4, MostPopular, data) <= 0.5

    # UserKNN(data, 2, true)
    @test cross_validation(5, Recall, 4, UserKNN, data, 2, true) <= 0.5
end

test_cross_validation()


text:
module solve_dist

include("steady_state.jl")
using .steady_state

function logit(x, fixedParam)
    return 1 ./ (1.0 .+ exp.(-fixedParam.lambda .* x))
end
function solveDist(x0, set, fxp, par, SmNC, SmPL, SmYC, SmD, SfNC, SfPL, SfYC, SfD, stages)
    #setting up
    x = copy(x0);
    v = x[length(x)-set.Ny*set.Na+1:length(x)];
    x = deleteat!(x, length(x)-set.Ny*set.Na+1:length(x));
    x = reshape(x,set.Nx*set.Ne,:);
    umNC    = x[:, 1];
    umPL    = x[:,2];
    umYC    = x[:,3];
    umD     = x[:,4];
    ufNC    = x[:,5];
    ufPL    = x[:,6];
    ufYC    = x[:,7];
    ufD     = x[:,8];
    hmNC    = x[:,9:29];
    hmYC    = x[:,30:50];
    hmPL    = x[:,51:71];
    hmD     = x[:,72:92];
    hfNC    = x[:,93:113];
    hfYC    = x[:,114:134];
    hfPL    = x[:,135:155];
    hfD     = x[:,156:176];
    D_f     = stages[4];
    D_m     = stages[8];
    gem = par.d1m .+ par.d2m .* fxp.p;
    gef = par.d1f .+ par.d2f .* fxp.p;
    #HC evoluation Stage
    #Male

    umNC_I  = umNC;
    umPL_I  = umPL;
    umYC_I  = umYC;
    umD_I   = umD;
    hmNC_I = [geth_I(x,y,gem,hmNC) for x in 1:21, y in 1:21];
    hmPL_I  = hmPL;
    hmYC_I = [geth_I(x,y,gem,hmYC) for x in 1:21, y in 1:21];
    hmD_I = [geth_I(x,y,gem,hmD) for x in 1:21, y in 1:21];
    #Female
    ufNC_I  = ufNC;
    ufPL_I  = ufPL;
    ufYC_I  = ufYC;
    ufD_I   = ufD;
    hfNC_I = [geth_I(x,y,gef,hfNC) for x in 1:21, y in 1:21];
    hfYC_I = [geth_I(x,y,gef,hfYC) for x in 1:21, y in 1:21];
    hfPL_I  = hfPL;
    hfD_I = [geth_I(x,y,gef,hfD) for x in 1:21, y in 1:21];

    #(II) Searching and Matching Stage
    #initial distribution of (x,ϵ)
    #initm = vec([x*y for x in fxp.phi0m, y in fxp.eden]);
    initm = vec([0.0259367011071109
    0.0749429899681410
    0.0448519457139780
    0.00904904918382681
    0.00246010759171049
    0.000858637573705890
    0.000351184093983957
    0.111558866891707
    0.322344580669522
    0.192917064551742
    0.0389217452607447
    0.0105814079527510
    0.00369316955142811
    0.00151051205137649
    0.0259367011071109
    0.0749429899681410
    0.0448519457139780
    0.00904904918382681
    0.00246010759171049
    0.000858637573705890
    0.000351184093983957]);
    #initf = vec([x*y for x in fxp.phi0f, y in fxp.eden]);
    initf = vec([0.0507800707384540
    0.0768648059334113
    0.0228879898217145
    0.00513562671455189
    0.00156507309681993
    0.000845274203814434
    0.000445523255319775
    0.218415099471133
    0.330610690171056
    0.0984457583635939
    0.0220893433859680
    0.00673168806483479
    0.00363569106190128
    0.00191628339055626
    0.0507800707384540
    0.0768648059334113
    0.0228879898217145
    0.00513562671455189
    0.00156507309681993
    0.000845274203814434
    0.000445523255319775]);

    umNC_II = [umNC_I[i]*(1 - fxp.gamma -fxp.chi- fxp.lambdau_NC*sum(v.*(SmNC[i,:].>0)))+fxp.phi*D_m*initm[i]+ fxp.deltam_NC*sum(hmNC_I[i,:]) for i in 1:21];
    ufNC_II = [ufNC_I[i]*(1 - fxp.gamma -fxp.chi- fxp.lambdau_NC*sum(v.*(SfNC[i,:].>0)))+fxp.phi*D_f*initf[i]+ fxp.deltaf_NC*sum(hfNC_I[i,:]) for i in 1:21];

    umPL_II = [umPL_I[i]*(1 - fxp.gamma -fxp.etam)+ fxp.chi*(umNC_I[i]+umYC_I[i]) + fxp.deltam_YC*sum(hmPL_I[i,:]) for i in 1:21];
    ufPL_II = [ufPL_I[i]*(1 - fxp.gamma -fxp.etaf)+ fxp.chi*(ufNC_I[i]+ufYC_I[i]) + fxp.deltaf_YC*sum(hfPL_I[i,:]) for i in 1:21];

    umYC_II = [umYC_I[i]*(1 - fxp.gamma -fxp.chi - fxp.lambdau_YC*sum(v.*(SmYC[i,:].>0)))+ fxp.etam*umPL_I[i] + fxp.deltam_YC*sum(hmYC_I[i,:]) for i in 1:21];
    ufYC_II = [ufYC_I[i]*(1 - fxp.gamma -fxp.chi - fxp.lambdau_YC*sum(v.*(SfYC[i,:].>0)))+ fxp.etaf*ufPL_I[i] + fxp.deltaf_YC*sum(hfYC_I[i,:]) for i in 1:21];

    umD_II = [umD_I[i]*(1 - fxp.phi - fxp.lambdau_YC*sum(v.*(SmD[i,:].>0)) )+ fxp.gamma*umNC_I[i] + fxp.gamma*umYC_I[i] + fxp.delta*sum(hmD_I[i,:]) for i in 1:21];
    ufD_II = [ufD_I[i]*(1 - fxp.phi - fxp.lambdau_YC*sum(v.*(SfD[i,:].>0)) )+ fxp.gamma*ufNC_I[i] + fxp.gamma*ufYC_I[i] + fxp.delta*sum(hfD_I[i,:]) for i in 1:21];

    hmNC_II = [get_hgNC_II_ij("m",i,j,umNC_I,SmNC,hmNC_I,fxp,v) for i in 1:21, j in 1:21];

    hmYC_II = [get_hgYC_II_ij("m",i,j,umYC_I,SmYC,hmYC_I,hmPL_I,fxp,v) for i in 1:21, j in 1:21];

    hmPL_II = [hmPL_I[i,j]+fxp.chi*(hmNC_I[i,j] + hmYC_I[i,j])- hmPL_I[i,j]*(fxp.gamma + fxp.deltam_YC + fxp.etam) for i in 1:21, j in 1:21];

    hmD_II = [get_hgD_II_ij(i,j,SmD,umD_I,hmNC_I,hmYC_I,hmPL_I,hmD_I,fxp,v) for i in 1:21, j in 1:21];

    hfNC_II = [get_hgNC_II_ij("f",i,j,ufNC_I,SfNC,hfNC_I,fxp,v) for i in 1:21, j in 1:21];

    hfYC_II = [get_hgYC_II_ij("f",i,j,ufYC_I,SfYC,hfYC_I,hfPL_I,fxp,v) for i in 1:21, j in 1:21];

    hfPL_II = [hfPL_I[i,j]+fxp.chi*(hfNC_I[i,j] + hfYC_I[i,j])- hfPL_I[i,j]*(fxp.gamma + fxp.deltaf_YC + fxp.etaf) for i in 1:21, j in 1:21];

    hfD_II = [get_hgD_II_ij(i,j,SfD,ufD_I,hfNC_I,hfYC_I,hfPL_I,hfD_I,fxp,v) for i in 1:21, j in 1:21];

    # Stage (III) Endogenous Quits
    # after endogenous quits, all measures must equal measures at the beginning
    umNCIII  = vec(umNC_II + sum(hmNC_II.*(SmNC .<=0), dims=2));
    umPLIII  = vec(umPL_II + sum(hmPL_II.*(SmPL .<=0),dims= 2));
    umYCIII  = vec(umYC_II + sum(hmYC_II.*(SmYC .<=0),dims= 2));
    umDIII   = vec(umD_II + sum(hmD_II.*(SmD .<=0), dims=2));
    ufNCIII  = vec(ufNC_II + sum(hfNC_II.*(SfNC .<=0),dims= 2));
    ufPLIII  = vec(ufPL_II + sum(hfPL_II.*(SfPL .<=0),dims= 2));
    ufYCIII  = vec(ufYC_II + sum(hfYC_II.*(SfYC .<=0),dims= 2));
    ufDIII   = vec(ufD_II + sum(hfD_II.*(SfD .<=0),dims= 2));
    hmNCIII  = hmNC_II - hmNC_II.*(SmNC .<=0);
    hmYCIII  = hmYC_II - hmYC_II.*(SmYC .<=0);
    hmPLIII  = hmPL_II - hmPL_II.*(SmPL .<=0);
    hmDIII   = hmD_II - hmD_II.*(SmD .<=0);
    hfNCIII  = hfNC_II - hfNC_II.*(SfNC .<=0);
    hfYCIII  = hfYC_II - hfYC_II.*(SfYC .<=0);
    hfPLIII  = hfPL_II - hfPL_II.*(SfPL .<=0);
    hfDIII   = hfD_II - hfD_II.*(SfD .<=0);

    #calculate new vacancies
    vnew    = fxp.gamma_yα - vec(sum(hmNC,dims=1)) - vec(sum(hmYC,dims=1)) -vec(sum(hmPL,dims=1)) - vec(sum(hmD,dims=1)) - vec(sum(hfNC,dims=1)) - vec(sum(hfYC,dims=1)) -vec(sum(hfPL,dims=1)) - vec(sum(hfD,dims=1));
    result = vcat(umNCIII, umPLIII, umYCIII, umDIII, ufNCIII, ufPLIII, ufYCIII, ufDIII, vec(hmNCIII), vec(hmYCIII), vec(hmPLIII), vec(hmDIII),vec(hfNCIII), vec(hfYCIII), vec(hfPLIII), vec(hfDIII),vnew);

    return result;
end
export solveDist
end


text:
@discord_object struct ActivityParty
    id::String
    size::Vector{Int}
end


text:
"""
$(TYPEDEF)

A `Storage` can be used to store simulation results for a [`Mechanism`](@ref).

Indexing example:
    storage.x[body_id][time_step_number][2] # returns the y-position of the specified body at the specified time step.

# Important attributes
* `x`: Contains the position for each body and each time step.
* `q`: Contains the orientation for each body and each time step.
* `v`: Contains the velocity for each body and each time step.
* `ω`: Contains the angluar velocity for each body and each time step.

# Constructors
    Storage{T}(step_range, nbodies)
    Storage(step_range, nbodies)
    Storage{T}(nsteps, nbodies)
    Storage(nsteps, nbodies)
"""
struct Storage{T,N}
    x::Vector{Vector{SVector{3,T}}}
    q::Vector{Vector{UnitQuaternion{T}}}
    v::Vector{Vector{SVector{3,T}}}
    ω::Vector{Vector{SVector{3,T}}}

    function Storage{T}(steps, nbodies) where T
        x = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        q = [[one(UnitQuaternion{T}) for i = steps] for j = 1:nbodies]
        v = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        ω = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        new{T,length(steps)}(x, q, v, ω)
    end

    Storage(steps, nbodies) = Storage{Float64}(steps, nbodies)
    Storage{T}(nsteps::Integer, nbodies) where T = Storage{T}(Base.OneTo(nsteps), nbodies)
    Storage(nsteps::Integer, nbodies) = Storage{Float64}(nsteps, nbodies)

    function Storage(x::Vector{<:Vector{<:AbstractVector{T}}},q::Vector{Vector{UnitQuaternion{T}}}) where T
        steps = Base.OneTo(length(x[1]))
        nbodies = length(x)

        v = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        ω = [[szeros(T, 3) for i = steps] for j = 1:nbodies]

        new{T,length(steps)}(x, q, v, ω)
    end

    Storage{T}() where T = Storage{T}(Base.OneTo(0),0)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, storage::Storage{T,N}) where {T,N}
    summary(io, storage)
end


text:
#using Compat.Serialization
#using Compat: @info

# Sets up a variety of different Gauss-Legendre data lengths
# for resolving different regimes of oscillatory function
#const GLH_NODES_20, GLH_WEIGHTS_20 = gausslegendre(20)
const GLH_NODES_100, GLH_WEIGHTS_100 = gausslegendre(100)
const GLH_NODES_200, GLH_WEIGHTS_200 = gausslegendre(200)
const GLH_NODES_500, GLH_WEIGHTS_500 = gausslegendre(500)
const GLH_NODES_1000, GLH_WEIGHTS_1000 = gausslegendre(1000)
const GLH_NODES_2000, GLH_WEIGHTS_2000 = gausslegendre(2000)
const GLH_NODES_2500, GLH_WEIGHTS_2500 = gausslegendre(2500)
const GLH_NODES_3500, GLH_WEIGHTS_3500 = gausslegendre(3500)
const GLH_NODES_10000, GLH_WEIGHTS_10000 = gausslegendre(10000)
const GLH_N = [100,200,500,1000,2000,2500,3500,10000]

const LGFH_DIR  = joinpath(pwd(), "cache")

using ProgressMeter

alpha_to_string(α::Float64) = string(100000+α*10000)[2:6]

lgfh_file(α) = joinpath(LGFH_DIR,"lgfhtable_alpha"*alpha_to_string(α)*".dat")

function _get_gl_data(n;i=0)
  # Set the number of Gauss points so that it will resolve the Nyquist frequency
  # on the oscillatory part of the Fourier integral
  while 2*n > GLH_N[i += 1] end
  return getfield(Fields,Symbol("GLH_NODES_",string(GLH_N[i]))), getfield(Fields,Symbol("GLH_WEIGHTS_",string(GLH_N[i])))
end

function load_lgf_helmholtz(N,α)
    if isfile(lgfh_file(α))
        G = deserialize(open(lgfh_file(α), "r"))
        if size(G,1) ≥ N
            return G
        end
    end
    build_lgf_helmholtz(N,α)
end

function build_lgf_helmholtz(N,α)
    @info "Building and caching Helmholtz LGF table for α = "*string(α)

    asymptotic_radius = _get_lgf_switch(α)
    @info "Switch to asymptotic formula at i = "*string(asymptotic_radius)

    lgf_support = _get_lgf_support(α)
    @info "Support radius = "*string(lgf_support)

    g = zeros(ComplexF64,N, N)
    @showprogress for y in 0:N-1, x in 0:y
      ir = _index_radius(x,y)
      if ir < asymptotic_radius
        g[x+1,y+1] = lgf_helmholtz(x,y,α)
      elseif ir < lgf_support
        g[x+1,y+1] = lgf_helmholtz_asymptotic(x,y,α)
      end
    end

    G = Symmetric(g)
    mkpath(LGFH_DIR)
    serialize(open(lgfh_file(α),"w"),G)
    G
end

#quadgauss_helmholtz(f::Function) = dot(GLH_WEIGHTS, f(GLH_NODES))

function lgf_helmholtz(i :: Integer, j :: Integer, α::Real)

  nodes_i, weights_i = _get_gl_data(abs(i))
  nodes_j, weights_j = _get_gl_data(abs(j))

  quadgauss_j(f::Function) = dot(weights_j, f(nodes_j))

  if i ≥ j
      v = ComplexF64(0.0)
      for (idex,y) in enumerate(nodes_i)
          xi1 = 0.5*π*(y+1)
          expixi1 = cos(i*xi1)
          cosxi1 = cos(xi1)
          v += weights_i[idex]*quadgauss_j() do x
              xi2 = 0.5*π*(x .+ 1)
              expixi2 = cos.(j*xi2)
              return expixi2.*expixi1./(im*α .+ 4 .- 2*cos.(xi2) .- 2*cosxi1)
          end
      end
      return v/4
  else
      return lgf_helmholtz(j,i,α)
  end

end

lgf_helmholtz_asymptotic(i :: Integer, j :: Integer, α::Real) =
      conj(im*0.25*hankelh1(0,exp(im*π/4)*sqrt(α)*sqrt(i^2+j^2)))

function _get_lgf_support(α::Number;tol=eps(Float64))
  i = 0
  while abs(lgf_helmholtz(i,0,α)) > tol
     i += 1
  end
  return i
end

function _get_lgf_switch(α::Number;tol=1e-9)
  i = 1
  while abs(lgf_helmholtz(i,0,α)-lgf_helmholtz_asymptotic(i,0,α)) > tol
     i += 1
  end
  return i
end

_index_radius(i::Int,j::Int) = ceil(Int,sqrt(i^2+j^2))


text:
# # [Variational Quantum Eigensolver](@id VQE)

# ## References
# * [A variational eigenvalue solver on a quantum processor](https://arxiv.org/abs/1304.3061)
# * [Variational Quantum Eigensolver with Fewer Qubits](https://arxiv.org/abs/1902.02663)

# ## Define a hamiltonian

# construct a 5-site heisenberg hamiltonian

using Yao, YaoExtensions

N = 5
hami = heisenberg(N)

# The ground state can be obtained by a sparse matrix ground state solver.
# The high performance `mat` function in `Yao.jl` makes computation time lower than `10s`
# to construct a `20` site Heisenberg hamltonian
using KrylovKit: eigsolve

function ed_groundstate(h::AbstractBlock)
    E, V = eigsolve(h |> mat, 1, :SR, ishermitian=true)
    E[1], V[1]
end

ed_groundstate(hami)

# Here we use the `heisenberg` hamiltonian that defined in [`YaoExtensions.jl`](https://github.com/QuantumBFS/YaoExtensions.jl),
# for tutorial purpose, we pasted the code for construction here.
# ```julia
# function heisenberg(nbit::Int; periodic::Bool=true)
#    sx = i->put(nbit, i=>X)
#    sy = i->put(nbit, i=>Y)
#    sz = i->put(nbit, i=>Z)
#    map(1:(periodic ? nbit : nbit-1)) do i
#        j=i%nbit+1
#        sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j)
#     end |> sum
# end
# ```

# ## Define an ansatz
# As an ansatz, we use the canonical circuit for demonstration [`variational_circuit`](@ref YaoExtensions.variational_circuit)
# defined in [`YaoExtensions.jl`](https://github.com/QuantumBFS/YaoExtensions.jl).
c = variational_circuit(N)
dispatch!(c, :random)

# ## Run
# Use the [`Adam`](@ref) optimizer for parameter optimization,
# we provide a poorman's implementation in `QuAlgorithmZoo`
using QuAlgorithmZoo: Adam, update!

optimizer = Adam(lr=0.01)
params = parameters(c)
niter = 100

for i = 1:niter
    ## `expect'` gives the gradient of an observable.
    grad_input, grad_params = expect'(hami, zero_state(N) => c)

    ## feed the gradients into the circuit.
    dispatch!(c, update!(params, grad_params, optimizer))
    println("Step $i, Energy = $(expect(hami, zero_state(N) |> c))")
end


text:
function update_weights(data, e, closest_same, closest_other, weights, weights_mult, max_f_vals, min_f_vals)

	# Penalty term
	penalty = sum(abs.(e .- closest_same)./(max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

	# Reward term
	reward = sum(weights_mult .* abs.(e .- closest_other)./(max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

	# Weights update
	weights = weights .- penalty./(size(data, 1)*size(closest_same, 1) + eps(Float64)) .+
		reward./(size(data, 1)*size(closest_other, 1) + eps(Float64))

	# Return updated weights.
	return vec(weights)
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function annuity(years::Int, rate_escalation::Float64, rate_discount::Float64)
    """
        this formulation assumes cost growth in first period
        i.e. it is a geometric sum of (1+rate_escalation)^n / (1+rate_discount)^n
        for n = 1, ..., years
    """
    x = (1 + rate_escalation) / (1 + rate_discount)
    if x != 1
        pwf = round(x * (1 - x^years) / (1 - x), digits=5)
    else
        pwf = years
    end
    return pwf
end

function levelization_factor(years::Int, rate_escalation::Float64, rate_discount::Float64,
    rate_degradation::Float64)
    #=
    NOTE: levelization_factor for an electricity producing tech is the ratio of:
    - an annuity with an escalation rate equal to the electricity cost escalation rate, starting year 1,
        and a negative escalation rate (the tech's degradation rate), starting year 2
    - divided by an annuity with an escalation rate equal to the electricity cost escalation rate (pwf_e).
    Both use the offtaker's discount rate.
    levelization_factor is multiplied by each use of dvRatedProduction in reopt.jl
        (except dvRatedProduction[t,ts] == dvSize[t] ∀ ts).
    This way the denominator is cancelled in reopt.jl when accounting for the value of energy produced
    since each value constraint uses pwf_e.

    :param analysis_period: years
    :param rate_escalation: escalation rate
    :param rate_discount: discount rate
    :param rate_degradation: positive degradation rate
    :return: present worth factor with escalation (inflation, or degradation if negative)
    NOTE: assume escalation/degradation starts in year 2
    =#
    num = 0
    for yr in range(1, stop=years)
        num += (1 + rate_escalation)^(yr) / (1 + rate_discount)^yr * (1 - rate_degradation)^(yr - 1)
    end
    den = annuity(years, rate_escalation, rate_discount)

    return num/den
end

function effective_cost(;
    itc_basis::Float64,
    replacement_cost::Float64,
    replacement_year::Int,
    discount_rate::Float64,
    tax_rate::Float64,
    itc::Float64,
    macrs_schedule::Array{Float64,1},
    macrs_bonus_pct::Float64,
    macrs_itc_reduction::Float64,
    rebate_per_kw::Float64=0.0,
    )

    """ effective PV and battery prices with ITC and depreciation
        (i) depreciation tax shields are inherently nominal --> no need to account for inflation
        (ii) ITC and bonus depreciation are taken at end of year 1
        (iii) battery replacement cost: one time capex in user defined year discounted back to t=0 with r_owner
        (iv) Assume that cash incentives reduce ITC basis
        (v) Assume cash incentives are not taxable, (don't affect tax savings from MACRS)
        (vi) Cash incentives should be applied before this function into "itc_basis".
             This includes all rebates and percentage-based incentives besides the ITC
    """

    # itc reduces depreciable_basis
    depr_basis = itc_basis * (1 - macrs_itc_reduction * itc)

    # Bonus depreciation taken from tech cost after itc reduction ($/kW)
    bonus_depreciation = depr_basis * macrs_bonus_pct

    # Assume the ITC and bonus depreciation reduce the depreciable basis ($/kW)
    depr_basis -= bonus_depreciation

    # Calculate replacement cost, discounted to the replacement year accounting for tax deduction
    replacement = replacement_cost * (1-tax_rate) / ((1 + discount_rate)^replacement_year)

    # Compute savings from depreciation and itc in array to capture NPV
    tax_savings_array = [0.0]
    for (idx, macrs_rate) in enumerate(macrs_schedule)
        depreciation_amount = macrs_rate * depr_basis
        if idx == 1
            depreciation_amount += bonus_depreciation
        end
        taxable_income = depreciation_amount
        push!(tax_savings_array, taxable_income * tax_rate)
    end

    # Add the ITC to the tax savings
    tax_savings_array[2] += itc_basis * itc

    # Compute the net present value of the tax savings
    tax_savings = npv(discount_rate, tax_savings_array)

    # Adjust cost curve to account for itc and depreciation savings ($/kW)
    cap_cost_slope = itc_basis - tax_savings + replacement - rebate_per_kw

    # Sanity check
    if cap_cost_slope < 0
        cap_cost_slope = 0
    end

    return round(cap_cost_slope, digits=4)
end

function dictkeys_tosymbols(d::Dict)
    d2 = Dict()
    for (k, v) in d
        if k in ["loads_kw", "prod_factor_series_kw"] && !isempty(v)
            try
                v = convert(Array{Real, 1}, v)
            catch
                @warn "Unable to convert $k to an Array{Real, 1}"
            end
        end
        d2[Symbol(k)] = v
    end
    return d2
end

function filter_dict_to_match_struct_field_names(d::Dict, s::DataType)
    f = fieldnames(s)
    d2 = Dict()
    for k in f
        if haskey(d, k)
            d2[k] = d[k]
        else
            @warn "dict is missing struct field $k"
        end
    end
    return d2
end

function npv(rate::Float64, cash_flows::Array)
    npv = cash_flows[1]
    for (y, c) in enumerate(cash_flows[2:end])
        npv += c/(1+rate)^y
    end
    return npv
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
__precompile__()
module IndexedTables

using Compat
using NamedTuples, PooledArrays

import Base:
    show, eltype, length, getindex, setindex!, ndims, map, convert, keys, values,
    ==, broadcast, empty!, copy, similar, sum, merge, merge!, mapslices,
    permutedims, reducedim, serialize, deserialize, sort, sort!

export NDSparse, flush!, aggregate!, aggregate_vec, where, pairs, convertdim, columns, column, rows,
    itable, update!, aggregate, reducedim_vec, dimlabels, collect_columns

const Tup = Union{Tuple,NamedTuple}
const DimName = Union{Int,Symbol}

include("utils.jl")
include("columns.jl")
include("table.jl")
include("ndsparse.jl")
include("collect.jl")

#=
# Poor man's traits

# These support `colnames` and `columns`
const TableTrait = Union{AbstractVector, NextTable, NDSparse}

# These support `colnames`, `columns`,
# `pkeynames`, `permcache`, `cacheperm!`
=#

const Dataset = Union{NextTable, NDSparse}

# no-copy convert
_convert(::Type{NextTable}, x::NextTable) = x
function _convert(::Type{NDSparse}, t::NextTable)
    NDSparse(rows(t, pkeynames(t)), rows(t, excludecols(t, pkeynames(t))),
             copy=false, presorted=true)
end

function _convert(::Type{NextTable}, x::NDSparse)
    convert(NextTable, x.index, x.data;
            perms=x._table.perms,
            presorted=true, copy=false)
end

ndsparse(t::NextTable; kwargs...) = _convert(NDSparse, t; kwargs...)
table(t::NDSparse; kwargs...) = _convert(NextTable, t; kwargs...)

include("sortperm.jl")
include("indexing.jl") # x[y]
include("selection.jl")
include("reduce.jl")
include("flatten.jl")
include("join.jl")
include("reshape.jl")

# TableTraits.jl integration
include("tabletraits.jl")

end # module


text:
make_globals(s) = s
lift_global(s::Symbol) = GlobalRef(@__MODULE__, s)
function lift_global(ex::Expr)
    if ex.head == :.
        return GlobalRef(
                         CodeInfoTools.resolve(lift_global(ex.args[1])),
                         ex.args[2].value)
    end
    return ex
end
function make_globals(ex::Expr)
    if ex.head == :call
        return Expr(ex.head,
                    lift_global(ex.args[1]),
                    map(make_globals, ex.args[2 : end])...)
    elseif ex.head == :(::)
        return Expr(ex.head, ex.args[1],
                    lift_global(ex.args[2]))
    elseif ex.head == :.
        return lift_global(ex)
    else
        return Expr(ex.head, map(make_globals, ex.args)...)
    end
end

make_patterm(s::Number) = PatLiteral(s)
make_patterm(s::Symbol) = PatVar(s)
function make_patterm(ex::Expr)
    if ex.head == :call
        return PatTerm(:call, [PatLiteral(ex.args[1]), map(make_patterm, ex.args[2 : end])...])
    elseif ex.head == :(::)
        @assert(ex.args[1] isa Symbol)
        return PatTypeAssertion(PatVar(ex.args[1]), ex.args[2])
    end
end

create_patterm(e) = e
create_patterm(e::Int) = PatLiteral(e)
function create_patterm(expr::Expr)
    resolved = make_globals(expr)
    resolved = walk(CodeInfoTools.resolve, resolved)
    return make_patterm(resolved)
end

function _methodtheory(ex::Expr)
    @assert(@capture(ex, begin body__ end))
    new = map(body) do sub
        if @capture(sub, lhs_ => rhs_)
            return Expr(:call, :push!, :arr,
                        RewriteRule(create_patterm(lhs), create_patterm(rhs)))
        elseif @capture(sub, lhs_ |> rhs_)
            return Expr(:call, :push!, :arr,
                        DynamicRule(create_patterm(lhs), rhs))
        end
    end
    return quote
        arr = Rule[]
        $(new...)
        arr
    end
end

macro methodtheory(ex)
    new = _methodtheory(ex)
    new = postwalk(rmlines ∘ unblock, new)
    esc(new)
end


text:
#TODO: Document interface to support static size
# Define maybestaticsize, maybestaticlength, and maybestaticfirstindex

struct Static{N} end
Base.@pure Static(N) = Static{N}()

Base.CartesianIndex(I::Tuple{<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)

abstract type AbstractStaticUnitRange <: AbstractUnitRange{Int} end
# rank 2, but 3 unknowns; 5 types can express all different posibilities
# 1: all unknown, UnitRange; 2: all three known (only two must be made explicit):
struct StaticUnitRange{L,U} <: AbstractStaticUnitRange end
Base.@pure StaticUnitRange(L,U) = StaticUnitRange{L,U}() # Do I use this definition anywhere?
# Then each specifying one of the three parameters
struct StaticLowerUnitRange{L} <: AbstractStaticUnitRange
    U::Int
end
struct StaticUpperUnitRange{U} <: AbstractStaticUnitRange
    L::Int
end
struct StaticLengthUnitRange{N} <: AbstractStaticUnitRange
    L::Int
end

@inline Base.first(::StaticUnitRange{L}) where {L} = L
@inline Base.first(::StaticLowerUnitRange{L}) where {L} = L
@inline Base.first(r::StaticUpperUnitRange) = r.L
@inline Base.first(r::StaticLengthUnitRange) = r.L
@inline Base.last(::StaticUnitRange{L,U}) where {L,U} = U
@inline Base.last(r::StaticLowerUnitRange) = r.U
@inline Base.last(::StaticUpperUnitRange{U}) where {U} = U
@inline Base.last(r::StaticLengthUnitRange{N}) where {N} = r.L + N - 1

@inline Base.iterate(x::AbstractStaticUnitRange) = (i = unwrap(first(x)); (i,i))

@inline Base.:(:)(::Static{L}, ::Static{U}) where {L,U} = StaticUnitRange{L,U}()
@inline Base.:(:)(::Static{L}, U::Int) where {L} = StaticLowerUnitRange{L}(U)
@inline Base.:(:)(L::Int, ::Static{U}) where {U} = StaticUpperUnitRange{U}(L)
@inline unwrap(x) = x
@inline unwrap(::Val{N}) where {N} = N
@inline unwrap(::Type{Val{N}}) where {N} = N
@inline unwrap(::Static{N}) where {N} = N
@inline unwrap(::Type{Static{N}}) where {N} = N

@inline StaticLowerUnitRange{L}(::Static{U}) where {L,U} = StaticUnitRange{L,U}()
@inline StaticUpperUnitRange{U}(::Static{L}) where {L,U} = StaticUnitRange{L,U}()
@generated StaticLengthUnitRange{N}(::Static{L}) where {L,N} = StaticUnitRange{L,L+N-1}()

@inline maybestaticsize(A, ::Val{I}) where {I} = size(A, I)
@inline maybestaticsize(A::AbstractArray{<:Any,0}, ::Val{1:2}) = (1, 1)
@inline maybestaticsize(A::AbstractVector, ::Val{1:2}) = (length(A), 1)
@inline maybestaticsize(A::AbstractMatrix, ::Val{1:2}) = size(A)
@inline maybestaticsize(A::AbstractArray, ::Val{1:2}) = (size(A,1),size(A,2))
@inline maybestaticlength(A) = length(A)

@inline maybestaticfirstindex(A::AbstractArray, ::Val{I}) where {I} = firstindex(A, I)
@inline maybestaticfirstindex(A::Array, ::Val{I}) where {I} = Static{1}()

# @inline maybestaticeachindex(A::AbstractArray) = maybestaticrange(eachindex(A))

@inline maybestaticrange(r::Base.OneTo) = StaticLowerUnitRange{1}(last(r))
@inline maybestaticrange(r::UnitRange) = r
@inline maybestaticrange(r::AbstractStaticUnitRange) = r
@inline maybestaticrange(r) = maybestaticfirst(r):maybestaticlast(r)
# @inline maybestaticaxis(A::AbstractArray, ::Val{I}) where {I} = maybestaticfirstindex(A, Val{I}()):maybestaticsize(A, Val{I}())

@inline maybestaticsize(::NTuple{N}, ::Val{1}) where {N} = Static{N}() # should we assert that i == 1?
@inline maybestaticlength(::NTuple{N}) where {N} = Static{N}()
@inline maybestaticsize(::Adjoint{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = Static{1}()
@inline maybestaticsize(::Transpose{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = Static{1}()
@inline maybestaticlength(B::Adjoint) = maybestaticlength(parent(B))
@inline maybestaticlength(B::Transpose) = maybestaticlength(parent(B))
@inline maybestaticsize(B::Adjoint{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
@inline maybestaticsize(B::Adjoint{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())
@inline maybestaticsize(B::Transpose{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
@inline maybestaticsize(B::Transpose{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())

@inline Base.iszero(::Static{0}) = true
@inline Base.iszero(::Static) = false
@generated function Base.divrem(N::Integer, ::Static{L}) where {L}
    if ispow2(L)
        quote
            $(Expr(:meta,:inline))
            d = N >>> $(intlog2(L))
            r = N & $(L-1)
            d, r
        end
    else
        quote
            $(Expr(:meta,:inline))
            divrem(N, $L)
        end
    end
end
@generated Base.divrem(::Static{N}, ::Static{D}) where {N,D} = divrem(N, D)

@inline Base.:+(::Static{N}, i) where {N} = N + i
@inline Base.:+(i, ::Static{N}) where {N} = N + i
@inline Base.:+(::Static{M}, ::Static{N}) where {M,N} = M + N
@inline Base.:*(::Static{N}, i) where {N} = N * i
@inline Base.:*(i, ::Static{N}) where {N} = N * i
@inline Base.:*(::Static{M}, ::Static{N}) where {M,N} = M * N
@inline Base.:-(::Static{N}, i) where {N} = N - i
@inline Base.:-(i, ::Static{N}) where {N} = i - N
@inline Base.:-(::Static{M}, ::Static{N}) where {M,N} = M - N
@inline Base.checked_add(::Static{N}, i) where {N} = Base.checked_add(N, i)
@inline Base.checked_add(i, ::Static{N}) where {N} = Base.checked_add(i, N)
@generated Base.checked_add(::Static{M}, ::Static{N}) where {M,N} = Static{Base.checked_add(M, N)}()
@inline Base.checked_sub(::Static{N}, i) where {N} = Base.checked_sub(N, i)
@inline Base.checked_sub(i, ::Static{N}) where {N} = Base.checked_sub(i, N)
@generated Base.checked_sub(::Static{M}, ::Static{N}) where {M,N} = Static{Base.checked_sub(M, N)}()
@inline Base.:>>(::Static{N}, i) where {N} = N >> i
@inline Base.:>>(i, ::Static{N}) where {N} = i >> N
@inline Base.:>>(::Static{M}, ::Static{N}) where {M,N} = M >> N
@inline Base.:<<(::Static{N}, i) where {N} = N << i
@inline Base.:<<(i, ::Static{N}) where {N} = i << N
@inline Base.:<<(::Static{M}, ::Static{N}) where {M,N} = M << N
@inline Base.:>>>(::Static{N}, i) where {N} = N >>> i
@inline Base.:>>>(i, ::Static{N}) where {N} = i >>> N
@inline Base.:>>>(::Static{M}, ::Static{N}) where {M,N} = M >>> N
@inline Base.:&(::Static{N}, i) where {N} = N & i
@inline Base.:&(i, ::Static{N}) where {N} = N & i
@inline Base.:&(::Static{M}, ::Static{N}) where {M,N} = M & N
@inline Base.:>(::Static{N}, i) where {N} = N > i
@inline Base.:>(i, ::Static{N}) where {N} = i > N
@inline Base.:>(::Static{M}, ::Static{N}) where {M,N} = M > N
@inline Base.:<(::Static{N}, i) where {N} = N < i
@inline Base.:<(i, ::Static{N}) where {N} = i < N
@inline Base.:<(::Static{M}, ::Static{N}) where {M,N} = M < N
@inline Base.:(==)(::Static{M}, i) where {M} = M == i
@inline Base.:(==)(i, ::Static{M}) where {M} = M == i
@inline Base.:(==)(::Static{M}, ::Static{N}) where {M,N} = false
@inline Base.:(==)(::Static{M}, ::Static{M}) where {M} = true

function static_promote(i, j)
    i == j || throw("$i ≠ $j")
    i
end
function static_promote(::Static{M}, i) where {M}
    M == i || throw("$M ≠ $i")
    Static{M}()
end
function static_promote(i, ::Static{M}) where {M}
    M == i || throw("$M ≠ $i")
    Static{M}()
end
static_promote(::Static{M}, ::Static{N}) where {M, N} = throw("$M ≠ $N")
static_promote(::Static{M}, ::Static{M}) where {M} = Static{M}()

@generated staticp1(::Static{N}) where {N} = Static{N+1}()
@inline staticp1(N::Integer) = N + 1
@generated staticm1(::Static{N}) where {N} = Static{N-1}()
@inline staticm1(N::Integer) = N - 1
@inline staticm1(i::Tuple{}) = tuple()
@inline staticm1(i::Tuple{I}) where {I} = @inbounds (i[1] - 1,)
@inline staticm1(i::Tuple{I1,I2}) where {I1,I2} = @inbounds (i[1] - 1, i[2] - 1)
@inline staticm1(i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = @inbounds (i[1] - 1, staticm1(Base.tail(i))...)
@inline Base.ntuple(f::F, ::Static{N}) where {F,N} = ntuple(f, Val{N}())

@inline maybestaticfirst(A) = first(A)
@inline maybestaticfirst(::StaticUnitRange{L}) where {L} = Static{L}()
@inline maybestaticfirst(::StaticLowerUnitRange{L}) where {L} = Static{L}()
@inline maybestaticfirst(::Base.OneTo) = Static{1}()

@inline _maybestaticfirst(A::Tuple{}) = tuple()
@inline _maybestaticfirst(A::Tuple{I}) where {I} = (maybestaticfirst(@inbounds(A[1])),)
@inline _maybestaticfirst(A::Tuple{I1,I2}) where {I1,I2} = @inbounds (maybestaticfirst(A[1]), maybestaticfirst(A[2]))
@inline _maybestaticfirst(A::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = (maybestaticfirst(@inbounds A[1]), _maybestaticfirst(Base.tail(A))...)
@inline maybestaticfirst(A::CartesianIndices) = CartesianVIndex(_maybestaticfirst(A.indices))

@inline maybestaticlast(A) = last(A)
@inline maybestaticlast(::StaticUnitRange{L,U}) where {L,U} = Static{U}()
@inline maybestaticlast(::StaticUpperUnitRange{U}) where {U} = Static{U}()
@inline _maybestaticlast(A::Tuple{}) = tuple()
@inline _maybestaticlast(A::Tuple{I}) where {I} = (maybestaticlast(@inbounds(A[1])),)
@inline _maybestaticlast(A::Tuple{I1,I2}) where {I1,I2} = @inbounds (maybestaticlast(A[1]), maybestaticlast(A[2]))
@inline _maybestaticlast(A::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = (maybestaticlast(@inbounds A[1]), _maybestaticlast(Base.tail(A))...)
@inline maybestaticlast(A::CartesianIndices) = CartesianVIndex(_maybestaticlast(A.indices))


text:
using Lazy
import Lazy: cycle, range, drop, take
using Test

# dummy function to test threading macros on
function add_things(n1, n2, n3)
    100n1 + 10n2 + n3
end

# dummy macro to test threading macros on
macro m_add_things(n1, n2, n3)
    quote
        100 * $(esc(n1)) + 10 * $(esc(n2)) + $(esc(n3))
    end
end

# define structs for @forward macro testing below (PR #112)
struct Foo112 end
struct Bar112 f::Foo112 end

@testset "Lazy" begin

if VERSION >= v"1.0.0"
    @test isempty(detect_ambiguities(Base, Core, Lazy))
end

@testset "Lists" begin
    @test list(1, 2, 3)[2] == 2
    @test prepend(1, list(2,3,4)) == 1:list(2, 3, 4)
    @test seq([1, 2, 3]) == list(1, 2, 3)
    @test seq(1:3) == list(1, 2, 3)
    @test constantly(1)[50] == 1
    testfn() = 1
    @test repeatedly(testfn)[50] == 1
    @test cycle([1, 2, 3])[50] == 2
    @test iterated(x->x^2, 2)[3] == 16
    @test range(1, 5)[3] == 3
    @test range(1, 5)[10] == nothing
    @test range(1, 5)[-1] == 1
    @test list(1, 2, 3) * list(4, 5, 6) == list(1, 2, 3, 4, 5, 6)
    @test first(list(1, 2, 3)) == 1
    @test tail(list(1, 2, 3)) == list(2, 3)
    @test flatten(list(1,2,list(3,4))) == list(1, 2, 3, 4)
    @test list(1,2,list(3,4))[3] == list(3, 4)
    @test list(list(1), list(2))[1] == list(1)
    @test reductions(+, 0, list(1, 2, 3)) == list(1, 3, 6)
    @test [i for i in @lazy[1,2,3]] == [1,2,3]

    l = list(1, 2, 3)
    @test l:7:l == list(list(1, 2, 3), 7, 1, 2, 3)   # ambiguity test
end

@testset "Fibs" begin
    fibs = @lazy 0:1:(fibs + drop(1, fibs));
    @test fibs[20] == 4181
    @test take(5, fibs) == list(0, 1, 1, 2, 3)
end

@testset "Primes" begin
    isprime(n) =
        @>> primes begin
            take_while(x -> x<=sqrt(n))
            map(x -> n % x == 0)
            any; !
        end
    primes = filter(isprime, range(2));
end

@testset "Even squares" begin
    esquares = @>> range() map(x->x^2) filter(iseven);
    @test take(5, esquares) == list(4, 16, 36, 64, 100)
end

@testset "Threading macros" begin
    temp = @> [2 3] sum
    @test temp == 5
    # Reverse from after index 2
    temp = @>> 2 reverse([1, 2, 3, 4, 5])
    @test temp == [1, 5, 4, 3, 2]
    temp = @as x 2 begin
        x^2
        x + 2
    end
    @test temp == 6

    # test that threading macros work with functions
    temp = @> 1 add_things(2,3)
    @test temp == 123

    temp = @>> 3 add_things(1,2)
    @test temp == 123

    temp = @as x 2 add_things(1,x,3)
    @test temp == 123

    # test that threading macros work with macros
    temp = @> 1 @m_add_things(2,3)
    @test temp == 123

    temp = @>> 3 @m_add_things(1,2)
    @test temp == 123

    temp = @as x 2 @m_add_things(1,x,3)
    @test temp == 123

end

@testset "Forward macro" begin
    play(x::Foo112; y) = y                        # uses keyword arg
    play(x::Foo112, z) = z                        # uses regular arg
    play(x::Foo112, z1, z2; y) = y + z1 + z2      # uses both

    @forward Bar112.f play                        # forward `play` function to field `f`

    let f = Foo112(), b = Bar112(f)
        @test play(f, y = 1) === play(b, y = 1)
        @test play(f, 2) === play(b, 2)
        @test play(f, 2, 3, y = 1) === play(b, 2, 3, y = 1)
    end
end

@testset "getindex" begin
    l = Lazy.range(1,10)
    @test l[1] == 1
    @test collect(l[1:5]) == collect(1:5)
end

@testset "Listables" begin
    @test_throws MethodError sin()
end

@static VERSION ≥ v"1.2" && @testset "avoid stackoverflow" begin
    @test (length(takewhile(<(10), Lazy.range(1))); true)
    @test (length(takewhile(<(100000), Lazy.range(1))); true)
end

@testset "any/all" begin
    let xs = list(true, false, false)
        @test any(identity, xs) == true
        @test any(xs) == true
        @test all(identity, xs) == false
        @test all(xs) == false
    end
    let yy = list(1, 0, 1)
        @test any(Bool, yy) == true
        @test all(Bool, yy) == false
    end
    # Base method--ensures no ambiguity with methods here
    @test all([true true; true true], dims=1) == [true true]
end

end


text:
import Pkg
Pkg.add("Conda")
using Conda
Conda.pip_interop(true)
Conda.pip("install", "webio_jupyter_extension")
using PlotlyJS

function batman()
    # reference: https://github.com/alanedelman/18.337_2015/blob/master/Lecture01_0909/The%20Bat%20Curve.ipynb
    σ(x) = @. √(1-x.^2)
    el(x) = @. 3*σ(x/7)
    s(x) = @. 4.2 - 0.5*x - 2.0*σ(0.5*x-0.5)
    b(x) = @. σ(abs(2-x)-1) - x.^2/11 + 0.5x - 3
    c(x) = [1.7, 1.7, 2.6, 0.9]

    p(i, f; kwargs...) = PlotlyJS.scatter(;x=[-i; 0.0; i], y=[f(i); NaN; f(i)],
                                  marker_color="black", showlegend=false,
                                  kwargs...)
    traces = vcat(p(3:0.1:7, el; name="wings 1"),
                  p(4:0.1:7, t->-el(t); name="wings 2"),
                  p(1:0.1:3, s; name="Shoulders"),
                  p(0:0.1:4, b; name="Bottom"),
                  p([0, 0.5, 0.8, 1], c; name="head"))

    PlotlyJS.plot(traces, Layout(title="Batman"))
end
batman()

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
export generate_deffile

# todo: add options (files for compiling, precompiling, and just using the work folder)
# last one is yet to be figured out
function generate_deffile(; excludepkgs = [], commit = "master", script = [])
    ppath = dirname(Base.active_project(false))
    cpath = joinpath(ppath, "container")
    if !isdir(cpath)
        mkdir(cpath)
        gitignore = joinpath(ppath, ".gitignore")
        io = open(gitignore, "a+");
        println(io, (raw"""

        ########################################
        #               Singularity                #
        ########################################

        *.sif
        container/containerhome
        """))
        close(io)

    end

    singjl_path = joinpath(cpath, "Singularity.pack")

    open(singjl_path, "w") do depsjl_file

        println(depsjl_file, strip(raw"""
        Bootstrap: library
        From: crown421/default/juliabase:latest

        %setup
            dir=`pwd`
            git clone \
            "file://$dir" \
            ${SINGULARITY_ROOTFS}/Project
        """))

        println(depsjl_file, strip(raw"""
        %post
            export JULIA_DEPOT_PATH=/user/.julia
            export PATH=/opt/julia/bin:$PATH

            cd Project
        """))

        # these are the variable things, unfortunately singularity build does not take arguments
        println(depsjl_file, ("""
            git checkout $commit -- src/ scripts/ Project.toml Manifest.toml
        """))

        println(depsjl_file, ("""
            julia --project -e 'using Pkg; Pkg.rm.($excludepkgs)'
        """))

        print(depsjl_file, (raw"""
            julia --project -e 'using Pkg; Pkg.instantiate()'

            julia --project -e 'using Pkg; Pkg.precompile()'

            chmod -R a+rX $JULIA_DEPOT_PATH
            chmod -R a+rX /Project/scripts

        %runscript
            if [ -z "$@" ]; then
        """))
        if isempty(script)
        print(depsjl_file, (raw"""
                # if theres none, test julia:
                julia --project=/Project -e 'using Pkg;  Pkg.status()'
        """))
        else
        print(depsjl_file, ("""
            SCRIPT=$script
        """))
        println(depsjl_file, (raw"""
                # set specific script if given
                julia --project=/Project -e "include(\\\"/Project/scripts/$SCRIPT\\\")" > "$SCRIPT-$(date +"%FT%H%M%S").log"
        """))
        end
        println(depsjl_file, (raw"""
            else
                # if theres an argument, then run it! and hope its a julia script :)
                julia --project=/Project -e "include(\\\"/Project/scripts/$@\\\")" > "$@-$(date +"%FT%H%M%S").log"
            fi
        """))
    end

end


text:
#  _|              _|
#      _|_|_|    _|_|_|_|    _|_|      _|_|_|    _|_|    _|  _|_|
#  _|  _|    _|    _|      _|_|_|_|  _|    _|  _|_|_|_|  _|_|
#  _|  _|    _|    _|      _|        _|    _|  _|        _|
#  _|  _|    _|      _|_|    _|_|_|    _|_|_|    _|_|_|  _|
#                                          _|
#                                      _|_|

function jacobi_pols(n, a, b, z)
    if (n+a > length(logfact) || n+b > length(logfact))
        error("Error: j is too high, please check the implementation of jacobi polynomials!")
    end

    # special case when I can not calculate log
    if (z ≈ 1) || (z ≈ -1)
        return sum((s % 2 == 0 ? 1.0 : -1.0) *
                exp(logfact[n+a+1] + logfact[n+b+1]-logfact[n-s+1]-logfact[a+s+1]-logfact[s+1]-logfact[n+b-s+1])*
                ((1-z)/2.0)^s*((1+z)/2.0)^(n-s) for s = 0:n)
    end

    # general case
    ls = log((1.0-z)/2.0);
    lc = log((1.0+z)/2.0);
    res = 0.0;
    for s = 0:n
        logs = logfact[(n+a+1)::Int] + logfact[(n+b+1)::Int]-
                    logfact[(n-s+1)::Int]-logfact[(a+s+1)::Int]-
                    logfact[(s+1)::Int]-logfact[(n+b-s+1)::Int];
        args = s*ls + (n-s)*lc;
        res += (s % 2 == 0 ? 1.0 : -1.0) * exp(logs+args);
    end
    return res;
end

function wignerd_hat(j, m1, m2, z)
        if abs(m1) > j || abs(m2) > j
            return zero(z)
        end
        factor = ((abs(m1-m2)+m1-m2)/2) % 2 == 0 ? one(z) : -one(z);
        am1 = abs(m1); am2 = abs(m2);
        M = (am1 > am2) ? am1 : am2;
        N = (am1 < am2) ? am1 : am2;
        gammas = logfact[Int(j-M+1)]+logfact[Int(j+M+1)]-(logfact[Int(j-N+1)]+logfact[Int(j+N+1)])
        return factor / 2^M * exp(gammas/2)*
               jacobi_pols(Int(j-M), Int(abs(m1-m2)), Int(abs(m1+m2)), z);
end

"""
    wignerd(j, m1, m2, cosβ)

Small wigner d-function for representation `j` with indices `m1`, `m2`, the argument `cosβ` is the cosine of the rotation angle.
The function gives the value of the matrix element
```
    ⟨ j m1 | exp(-i Jy β)| j m2 ⟩.
```
The input values are expected to be __integers__.
For a general case including half-integers see `wignerd_doublearg`.
"""
function wignerd(j, m1, m2, z)
    (z ≈ 1) && return m1 == m2 ? one(z) : zero(z)
    (z ≈ -1) && return m1 == -m2 ? (iseven(j-m2) ? one(z) : -one(z)) : zero(z)
    #
    hat = wignerd_hat(j, m1, m2, z);
    xi = sqrt(1-z)^abs(Int(m1-m2))*sqrt(1+z)^abs(Int(m1+m2));
    return hat*xi;
end

"""
    wignerD(j, m1, m2, α, cosβ, γ)

Wigner D-function for representation `j` with indices `m1`, `m2`, α, β, and γ are the rotation angles
The function gives the value of the matrix element
```
    ⟨ j m1 | exp(-i Jz α) exp(-i Jy β) exp(-i Jz γ) | j m2 ⟩.
```
The input values are expected to be __integers__.
For a general case including half-integers see `wignerD_doublearg`.
"""
wignerD(j, m1, m2, α, cosβ, γ)  = wignerd(j, m1, m2, cosβ) * cis(-m1*α-m2*γ)

#
#  _|                  _|      _|_|              _|              _|
#  _|_|_|      _|_|_|  _|    _|                      _|_|_|    _|_|_|_|    _|_|      _|_|_|    _|_|    _|  _|_|
#  _|    _|  _|    _|  _|  _|_|_|_|  _|_|_|_|_|  _|  _|    _|    _|      _|_|_|_|  _|    _|  _|_|_|_|  _|_|
#  _|    _|  _|    _|  _|    _|                  _|  _|    _|    _|      _|        _|    _|  _|        _|
#  _|    _|    _|_|_|  _|    _|                  _|  _|    _|      _|_|    _|_|_|    _|_|_|    _|_|_|  _|
#                                                                                _|
#                                                                                _|_|

function wignerd_hat_doublearg(two_j, two_m1, two_m2, z)
        # @show j, m1, m2
        if abs(two_m1) > two_j || abs(two_m2) > two_j
            return zero(z)
        end
        factor = (abs(two_m1-two_m2)+two_m1-two_m2) % 8 == 4 ? -one(z) : one(z);
        two_am1 = abs(two_m1); two_am2 = abs(two_m2);
        (two_M, two_N) = (two_am1 > two_am2) ? (two_am1,two_am2) : (two_am2,two_am1);
        #
        j_mnus_M = div(two_j-two_M,2)
        j_plus_M = div(two_j+two_M,2)
        j_mnus_N = div(two_j-two_N,2)
        j_plus_N = div(two_j+two_N,2)
        #
        gammas = logfact[j_mnus_M+1] +
                 logfact[j_plus_M+1] -
                 logfact[j_mnus_N+1] -
                 logfact[j_plus_N+1]
        return factor / 2^(two_M/2) * exp(gammas/2) *
               jacobi_pols(j_mnus_M, div(abs(two_m1-two_m2),2), div(abs(two_m1+two_m2),2), z);
end

"""
    wignerd_doublearg(j, m1, m2, cosβ)

Small wigner d-function for representation `j` with indices `m1`, `m2`, the argument `cosβ` is the cosine of the rotation angle.
The function gives the value of the matrix element
```
    ⟨ j m1 | exp(-i Jy β)| j m2 ⟩.
```
The function requires __doubled value of the momenta__ for the input.
"""
function wignerd_doublearg(two_j, two_m1, two_m2, z)
    (z ≈ 1) && return two_m1 == two_m2 ? one(z) : zero(z)
    (z ≈ -1) && return two_m1 == -two_m2 ? (iseven(div(two_j-two_m2,2)) ? one(z) : -one(z)) : zero(z)
    #
    hat = wignerd_hat_doublearg(two_j, two_m1, two_m2, z);
    xi = (1-z)^(abs(two_m1-two_m2)/4)*(1+z)^(abs(two_m1+two_m2)/4);
    return hat*xi;
end

"""
    wignerD_doublearg(two_j, two_m1, two_m2, α, cosβ, γ)

Wigner D-function for representation `j` with indices `m1`, `m2`, α, β, and γ are the rotation angles
The function gives the value of the matrix element
```
    ⟨ j m1 | exp(-i Jz α) exp(-i Jy β) exp(-i Jz γ) | j m2 ⟩.
```
The function requires __doubled value of the momenta__ for the input.
"""
wignerD_doublearg(two_j, two_m1, two_m2, α, cosβ, γ) =
    wignerd_doublearg(two_j, two_m1, two_m2, cosβ) * cis(-two_m1*α/2-two_m2*γ/2);

"""
    kronecker(i, j) = (i==j) ? 1 : 0
"""
kronecker(i, j) = (i==j) ? 1 : 0


text:
"""
    cylinder()

Return a reference cylinder (`RefMesh`) of 1m radius and 1m length. The cylinder axis is oriented
parallel to the z axis, and the origin (center of the bottom cap of the cylinder) is positionned
at x = 0, y = 0 and z = 0.
"""
function cylinder()
    file = joinpath(dirname(dirname(pathof(PlantGeom))), "data", "cylinder_124faces.ply")

    return RefMesh("cylinder", read_ply(file))
end

"""
    read_ply(file)

Read a ply file into a `Meshes.Meshes.SimpleMesh` format. Code taken from
[Meshes.jl](https://juliageometry.github.io/Meshes.jl/stable/algorithms/smoothing.html)
documentation.

Note: Remove when we have cylinder triangulation in Meshes.jl (and remove the dependency to PlyIO)
"""
function read_ply(file)
    ply = load_ply(file)
    x = ply["vertex"]["x"]
    y = ply["vertex"]["y"]
    z = ply["vertex"]["z"]
    points = Meshes.Point3.(x, y, z)
    connec = [Meshes.connect(Tuple(c .+ 1)) for c in ply["face"]["vertex_indices"]]
    Meshes.SimpleMesh(points, connec)
end


text:
include("divisor.jl")

function compute()

	triangle=0
	for i=1:100000000
		triangle+=i;
		if num_divisors(triangle)>500
			return triangle
		end
	end
end

print(compute())


text:
struct MTData
    src::TextReader        # reader for source language data
    tgt::TextReader        # reader for target language data
    batchsize::Int         # desired batch size
    maxlength::Int         # skip if source sentence above maxlength
    batchmajor::Bool       # batch dims (B,T) if batchmajor=false (default) or (T,B) if true.
    bucketwidth::Int       # batch sentences with length within bucketwidth of each other
    buckets::Vector        # sentences collected in separate arrays called buckets for each length range
    batchmaker::Function   # function that turns a bucket into a batch.
end

function MTData(src::TextReader, tgt::TextReader; batchmaker = arraybatch, batchsize = 128, maxlength = typemax(Int),
                batchmajor = false, bucketwidth = 10, numbuckets = min(128, maxlength ÷ bucketwidth))
    buckets = [ [] for i in 1:numbuckets ] # buckets[i] is an array of sentence pairs with similar length
    MTData(src, tgt, batchsize, maxlength, batchmajor, bucketwidth, buckets, batchmaker)
end

Base.IteratorSize(::Type{MTData}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{MTData}) = Base.HasEltype()
Base.eltype(::Type{MTData}) = NTuple{2}

function Base.iterate(d::MTData, state=nothing)
    if state === nothing
        for b in d.buckets; empty!(b); end
    end
    bucket,ibucket = nothing,nothing
    while true
        if state === nothing
            iter1,iter2 = iterate(d.src),iterate(d.tgt)
        else
            iter1,iter2 = iterate(d.src, state[1]), iterate(d.tgt, state[2])
        end
        if iter1 === nothing || iter2 === nothing
            ibucket = findfirst(x -> !isempty(x), d.buckets)
            bucket = (ibucket === nothing ? nothing : d.buckets[ibucket])
            break
        else
            sent1, state1 = iter1; sent2, state2 = iter2
            state = (state1, state2)
            ## Filter and sort based on src sentence length
            if length(sent1) > d.maxlength || length(sent1) == 0; continue; end
            ibucket = min(1 + (length(sent1)-1) ÷ d.bucketwidth, length(d.buckets))
            bucket = d.buckets[ibucket]
            push!(bucket, (sent1, sent2))
            if length(bucket) === d.batchsize; break; end
        end
    end
    if bucket === nothing; return nothing; end
    batch = d.batchmaker(d, bucket)
    empty!(bucket)
    return batch, state
end

function arraybatch(d::MTData, bucket)
    src,tgt = (x->x[1]).(bucket), (x->x[2]).(bucket)
    srclen,tgtlen = maximum(length.(src)), maximum(length.(tgt))+2
    batchsize = length(bucket)
    srcout = fill(d.src.vocab.eos, batchsize, srclen)
    tgtout = fill(d.tgt.vocab.eos, batchsize, tgtlen)
    for i in 1:batchsize
        srcout[i,end-length(src[i])+1:end] = src[i]
        tgtout[i, 1] = d.tgt.vocab.eos
        tgtout[i, 2:(1+length(tgt[i]))] = tgt[i]
        tgtout[i, 2+length(tgt[i])] = d.tgt.vocab.eos
    end
    if d.batchmajor
        srcout, tgtout = permutedims.((srcout, tgtout))
    end
    return (srcout, tgtout)
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
export RecordOptions
export RecordAction
export RecordGroup, RecordEvery
export RecordChange, RecordCost, RecordIterate, RecordIteration
export RecordEntry, RecordEntryChange
export get_record, has_record
export RecordActionFactory, RecordFactory

#
#
# record Options Decorator
#
#
@doc raw"""
    RecordAction

A `RecordAction` is a small functor to record values.
The usual call is given by `(p,o,i) -> s` that performs the record based on
a [`Problem`](@ref) `p`, [`Options`](@ref) `o` and the current iterate `i`.

By convention `i<=0` is interpreted as "For Initialization only", i.e. only
initialize internal values, but not trigger any record, the same holds for
`i=typemin(Inf)` which is used to indicate `stop`, i.e. that the record is
called from within [`stop_solver!`](@ref) which returns true afterwards.

# Fields (assumed by subtypes to exist)
* `recordedValues` an `Array` of the recorded values.
"""
abstract type RecordAction <: AbstractOptionsAction end

@doc raw"""
    RecordOptions <: Options

append to any [`Options`](@ref) the decorator with record functionality,
Internally a `Dict`ionary is kept that stores a [`RecordAction`](@ref) for
several occasions using a `Symbol` as reference.
The default occasion is `:All` and for example solvers join this field with
`:Start`, `:Step` and `:Stop` at the beginning, every iteration or the
end of the algorithm, respectively

The original options can still be accessed using the [`get_options`](@ref) function.

# Fields
* `options` – the options that are extended by debug information
* `recordDictionary` – a `Dict{Symbol,RecordAction}` to keep track of all
  different recorded values

# Constructors
    RecordOptions(o,dR)

construct record decorated [`Options`](@ref), where `dR` can be

* a [`RecordAction`](@ref), then it is stored within the dictionary at `:All`
* an `Array` of [`RecordAction`](@ref)s, then it is stored as a
  `recordDictionary`(@ref) within the dictionary at `:All`.
* a `Dict{Symbol,RecordAction}`.
"""
mutable struct RecordOptions{O <: Options} <: Options
    options::O
    recordDictionary::Dict{Symbol, <: RecordAction}
    RecordOptions{O}(o::O, dR::Dict{Symbol, <: RecordAction}) where {O <: Options} = new(o,dR)
end
RecordOptions(o::O, dR::D) where {O <: Options, D <: RecordAction} = RecordOptions{O}(o,Dict(:All => dR))
RecordOptions(o::O, dR::Array{ <: RecordAction,1}) where {O <: Options} = RecordOptions{O}(o,Dict(:All => RecordGroup(dR)))
RecordOptions(o::O, dR::Dict{Symbol, <: RecordAction}) where {O <: Options} = RecordOptions{O}(o,dR)
RecordOptions(o::O, format::Array{<:Any,1}) where {O <: Options} = RecordOptions{O}(o, RecordFactory(get_options(o),format))

dispatch_options_decorator(o::RecordOptions) = Val(true)

"""
    has_record(o)

check whether the [`Options`](@ref)` o` are decorated with
[`RecordOptions`](@ref)
"""
has_record(o::RecordOptions) = true
has_record(o::Options) = has_record(o,dispatch_options_decorator(o))
has_record(o::Options, ::Val{true}) = has_record(o.options)
has_record(o::Options, ::Val{false}) = false

# default - stored in the recordedValues field of the RecordAction
@doc raw"""
    get_record(o[,s=:Step])

return the recorded values from within the [`RecordOptions`](@ref) `o` that where
recorded with respect to the `Symbol s` as an `Array`. The default refers to
any recordings during an Iteration represented by the Symbol `:Step`
"""
function get_record(o::RecordOptions,s::Symbol=:Step)
    if haskey(o.recordDictionary,s)
        return get_record(o.recordDictionary[s])
    elseif haskey(o.recordDictionary,:All)
        return get_record(o.recordDictionary[:All])
    else
        error("No record known for key found, since neither :$s nor :All are present.")
    end
end
get_record(o::Options, s::Symbol=:Step) = get_record(o,s,dispatch_options_decorator(o))
get_record(o::Options, s, ::Val{true}) = get_record(o.options,s)
get_record(o::Options, s, ::Val{false}) = error("No Record decoration found")

@doc raw"""
    get_record(r)

return the recorded values stored within a [`RecordAction`](@ref) `r`.
"""
get_record(r::R) where {R <: RecordAction} = r.recordedValues

"""
    record_or_eset!(r,v,i)

either record (`i>0` and not `Inf`) the value `v` within the [`RecordAction`](@ref) `r`
or reset (`i<0`) the internal storage, where `v` has to match the internal
value type of the corresponding Recordaction.
"""
function record_or_eset!(r::R,v,i::Int) where {R <: RecordAction}
    if i > 0
        push!(r.recordedValues,v)
    elseif i < 0 && i > typemin(Int) # reset if negative but not stop indication
        r.recordedValues = Array{typeof(v),1}()
    end
end
"""
    RecordGroup <: RecordAction

group a set of [`RecordAction`](@ref)s into one action, where the internal prints
are removed by default and the resulting strings are concatenated

# Constructor
    RecordGroup(g)

construct a group consisting of an Array of [`RecordAction`](@ref)s `g`,
that are recording `en bloque`; the method does not perform any record itself,
but keeps an array of records. Accessing these yields a `Tuple` of the recorded
values per iteration
"""
mutable struct RecordGroup <: RecordAction
  group::Array{RecordAction,1}
  RecordGroup(g::Array{<:RecordAction,1}) = new(g)
  RecordGroup() = new(Array{RecordAction,1}())
end
function (d::RecordGroup)(p::P,o::O,i::Int) where {P <: Problem, O <: Options}
    for ri in d.group
        ri(p,o,i)
    end
end
get_record(r::RecordGroup) = [zip( get_record.(r.group)...)...]

@doc raw"""
    RecordEvery <: RecordAction

record only every $i$th iteration.
Otherwise (optionally, but activated by default) just update internal tracking
values.

This method does not perform any record itself but relies on it's childrens methods
"""
mutable struct RecordEvery <: RecordAction
    record::RecordAction
    every::Int
    alwaysUpdate::Bool
    RecordEvery(r::RecordAction,every::Int=1,alwaysUpdate::Bool=true) = new(r,every,alwaysUpdate)
end
function (d::RecordEvery)(p::P,o::O,i::Int) where {P <: Problem, O <: Options}
    if i<=0
        d.record(p,o,i)
    elseif (rem(i,d.every)==0)
        d.record(p,o,i)
    elseif d.alwaysUpdate
        d.record(p,o,0)
    end
end
get_record(r::RecordEvery) = get_record(r.record)
#
# Special single ones
#
@doc raw"""
    RecordChange <: RecordAction

debug for the amount of change of the iterate (stored in `o.x` of the [`Options`](@ref))
during the last iteration.

# Additional Fields
* `storage` a [`StoreOptionsAction`](@ref) to store (at least) `o.x` to use this
  as the last value (to compute the change)
"""
mutable struct RecordChange <: RecordAction
    recordedValues::Array{Float64,1}
    storage::StoreOptionsAction
    RecordChange(a::StoreOptionsAction=StoreOptionsAction( (:x,) ) ) = new(Array{Float64,1}(),a)
    function RecordChange(x0, a::StoreOptionsAction=StoreOptionsAction((:x,)))
        update_storage!(a,Dict(:x=>x0))
        return new(Array{Float64,1}(),a)
    end
end
function (r::RecordChange)(p::P,o::O,i::Int) where {P <: Problem, O <: Options}
    record_or_eset!(r,
        has_storage(r.storage, :x) ? distance(p.M,o.x, get_storage(r.storage,:x) ) : 0.0,
        i
    )
    r.storage(p,o,i)
    return r.recordedValues
end

@doc raw"""
    RecordEntry{T} <: RecordAction

record a certain fields entry of type {T} during the iterates

# Fields
* `recordedValues` – the recorded Iterates
* `field` – Symbol the entry can be accessed with within [`Options`](@ref)

"""
mutable struct RecordEntry{T} <: RecordAction
    recordedValues::Array{T,1}
    field::Symbol
    RecordEntry{T}(f::Symbol) where T = new(Array{T,1}(),f)
end
RecordEntry(e::T,f::Symbol) where T = RecordEntry{T}(f)
RecordEntry(d::DataType,f::Symbol) = RecordEntry{d}(f)
(r::RecordEntry{T})(p::Pr,o::O,i::Int) where {T, Pr <: Problem, O <: Options} = record_or_eset!(r, getfield(o, r.field), i)

@doc raw"""
    RecordEntryChange{T} <: RecordAction

record a certain entries change during iterates

# Additional Fields
* `recordedValues` – the recorded Iterates
* `field` – Symbol the field can be accessed with within [`Options`](@ref)
* `distance` – function (p,o,x1,x2) to compute the change/distance between two values of the entry
* `storage` – a [`StoreOptionsAction`](@ref) to store (at least) `getproperty(o, d.field)`
"""
mutable struct RecordEntryChange <: RecordAction
    recordedValues::Array{Float64,1}
    field::Symbol
    distance::Function
    storage::StoreOptionsAction
    RecordEntryChange(
            f::Symbol,
            d::Function,
            a::StoreOptionsAction=StoreOptionsAction( (f,) )
        ) = new(Array{Float64,1}(),f,d,a)
    function RecordEntryChange(v::T where T, f::Symbol, d::Function,
            a::StoreOptionsAction=StoreOptionsAction( (f,) )
        )
        update_storage!(a,Dict(f=>v))
        return new(Array{Float64,1}(),f, d, a)
    end
end
function (r::RecordEntryChange)(p::P,o::O,i::Int) where {P <: Problem, O <: Options}
    record_or_eset!(r,
        has_storage(r.storage, r.field) ? r.distance(p,o, getfield(o, r.field), get_storage(r.storage, r.field) ) : 0.0,
    i)
    r.storage(p,o,i)
end

@doc raw"""
    RecordIterate <: RecordAction

record the iterate

# Constructors
    RecordIterate(x0)

initialize the iterate record array to the type of `x0`, e.g. your initial data.

    RecordIterate(P)

initialize the iterate record array to the data type `T`.
"""
mutable struct RecordIterate{T} <: RecordAction
    recordedValues::Array{T,1}
    RecordIterate{T}() where {T} = new(Array{T,1}())
end
RecordIterate(p::T) where {T} = RecordIterate{T}()
RecordIterate() = throw(ErrorException("The iterate's data type has to be provided, i.e. RecordIterate(x0)."))

function (r::RecordIterate{T})(p::Pr,o::O,i::Int) where {Pr <: Problem, O <: Options,T}
    return record_or_eset!(r, o.x, i)
end

@doc raw"""
    RecordIteration <: RecordAction

record the current iteration
"""
mutable struct RecordIteration <: RecordAction
    recordedValues::Array{Int,1}
    RecordIteration() = new(Array{Int,1}())
end
(r::RecordIteration)(p::P,o::O,i::Int) where {P <: Problem, O <: Options} = record_or_eset!(r, i, i)

@doc raw"""
    RecordCost <: RecordAction

record the current cost function value, see [`get_cost`](@ref).
"""
mutable struct RecordCost <: RecordAction
    recordedValues::Array{Float64,1}
    RecordCost() = new(Array{Float64,1}())
end
(r::RecordCost)(p::P,o::O,i::Int) where {P <: Problem, O <: Options} = record_or_eset!(r, get_cost(p,o.x), i)

@doc raw"""
    RecordFactory(a)

given an array of `Symbol`s and [`RecordAction`](@ref)s and `Ints`

* The symbol `:Cost` creates a [`RecordCost`](@ref)
* The symbol `:iteration` creates a [`RecordIteration`](@ref)
* The symbol `:Change` creates a [`RecordChange`](@ref)
* any other symbol creates a [`RecordEntry`](@ref) of the corresponding field in [`Options`](@ref)
* any [`RecordAction`](@ref) is directly included
* an Integer `k` introduces that record is only performed every `k`th iteration
"""
function RecordFactory(o::O, a::Array{<:Any,1} ) where {O <: Options}
    # filter out every
    group = Array{RecordAction,1}()
    for s in filter(x -> !isa(x,Int), a) # filter ints and stop
        push!(group,RecordActionFactory(o,s) )
    end
    record = RecordGroup(group)
    # filter ints
    e = filter(x -> isa(x,Int),a)
    if length(e) > 0
        record = RecordEvery(record,last(e))
    end
    dictionary = Dict{Symbol,RecordAction}(:All => record)
    return dictionary
end
@doc raw"""
    RecordActionFactory(s)

create a [`RecordAction`](@ref) where

* a [`RecordAction`](@ref) is passed through
* a [`Symbol`] creates [`RecordEntry`](@ref) of that symbol, with the exceptions
  of `:Change`, `:Iterate`, `:Iteration`, and `:Cost`.
"""
RecordActionFactory(o::O,a::A) where {O <: Options, A <: RecordAction} = a
function RecordActionFactory(o::O,s::Symbol) where {O <: Options}
    if (s==:Change)
        return RecordChange()
    elseif (s==:Iteration)
        return RecordIteration()
    elseif (s==:Iterate)
        return RecordIterate(o.x)
    elseif (s==:Cost)
        return RecordCost()
    end
        return RecordEntry(getfield(o,s),s)
end


text:
"""
  Methods for three-dimensional visualization of wind farm layout, wind turbine
  design, and wake models.

  # AUTHORSHIP
    * Author    : Eduardo J Alvarez
    * Email     : Edo.AlvarezR@gmail.com
    * Created   : Jul 2018
    * License   : MIT License
"""
module WFVisual
export generate_windfarm

# ------------ GENERIC MODULES -------------------------------------------------
import JLD
import CSV
import Dierckx
import PyPlot
import PyCall

const plt = PyPlot

# ------------ FLOW MODULES ----------------------------------------------------
# https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
const gt = GeometricTools

# ------------ GLOBAL VARIABLES ------------------------------------------------
const module_path = splitdir(@__FILE__)[1]          # Path to this module
const def_data_path = joinpath(module_path, "../data")  # Path to data files

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["turbine", "farm"]
  include("WFVisual_"*header_name*".jl")
end

end # END OF MODULE

max_issues_repo_name: byuflowlab/wind-farm-utilities

max_forks_repo_name: byuflowlab/wind-farm-utilities


text:
# SplineSpace represents a Spline, right now piecewise constant HeavisideSpace is only implemented case
struct SplineSpace{order,T,R} <: Space{PiecewiseSegment{T},R}
    domain::PiecewiseSegment{T}
end

SplineSpace{m,T}(d::PiecewiseSegment{T}) where {m,T} = SplineSpace{m,T,real(eltype(T))}(d)
SplineSpace{m,T}(d::AbstractVector) where {m,T} = SplineSpace{m}(PiecewiseSegment(sort(d)))

SplineSpace{m}(d::PiecewiseSegment{T}) where {m,T} = SplineSpace{m,T,real(eltype(T))}(d)
SplineSpace{m}(d::AbstractVector) where {m} = SplineSpace{m}(PiecewiseSegment(sort(d)))

const HeavisideSpace{T,R} = SplineSpace{0,T,R}
dimension(h::SplineSpace{λ}) where {λ} = length(h.domain.points)+λ-1

convert(::Type{HeavisideSpace},d::PiecewiseSegment) = HeavisideSpace{eltype(d)}(d)

convert(::Type{HeavisideSpace},d::AbstractVector) =
    HeavisideSpace(PiecewiseSegment(sort(d)))

spacescompatible(a::SplineSpace{λ},b::SplineSpace{λ}) where {λ} = domainscompatible(a,b)

function evaluate(c::AbstractVector{T}, s::HeavisideSpace{<:Real}, x::Real) where T
    p = domain(s).points
    for k=1:length(c)
        if p[k] ≤ x ≤ p[k+1]
            return c[k]
        end
    end
    return zero(T)
end

function evaluate(c::AbstractVector{T}, s::SplineSpace{1,<:Real}, x::Real) where T
    p = domain(f).points
    c = f.coefficients
    for k=1:length(p)-1
        if p[k] ≤ x ≤ p[k+1]
            return (x-p[k])*c[k+1]/(p[k+1]-p[k]) + (p[k+1]-x)*c[k]/(p[k+1]-p[k])
        end
    end
    return zero(T)
end

function points(sp::HeavisideSpace,n)
    x=sp.domain.points
    (x[1:end-1] + diff(x)/2)[1:n]
end

points(sp::SplineSpace{1},n) = sp.domain.points[1:n]

for λ = [0,1]
    @eval begin
        function transform(S::SplineSpace{$λ},vals::AbstractVector,plan...)
            @assert length(vals) ≤ dimension(S)
            vals
        end
        itransform(S::SplineSpace{$λ},cfs::AbstractVector,plan...) = pad(cfs,dimension(S))
    end
end

bandwidths(D::ConcreteDerivative{H}) where {H<:HeavisideSpace} = (0,1)
rangespace(D::ConcreteDerivative{H}) where {H<:HeavisideSpace} = DiracSpace(domain(D).points[2:end-1])

function getindex(D::ConcreteDerivative{H,<:Any,T},k::Integer,j::Integer) where {H<:HeavisideSpace,T}
    if k==j
        -one(T)
    elseif j==k+1
        one(T)
    else
        zero(T)
    end
end

Base.sum(f::Fun{HS}) where {HS<:HeavisideSpace} = dotu(f.coefficients,diff(space(f).domain.points))
function Base.sum(f::Fun{SplineSpace{1,T,R}}) where {T,R}
    vals=pad(f.coefficients,dimension(space(f)))
    dfs=diff(space(f).domain.points)
    ret=vals[1]*dfs[1]/2
    for k=2:length(vals)-1
        ret+=vals[k]*(dfs[k]+dfs[k+1])/2
    end
    ret+=vals[end]*dfs[end]/2
    ret
end

#diffentiate HeavisideSpace
function differentiate(f::Fun{<:HeavisideSpace})
    dp=domain(f).points
    cfs=f.coefficients
    diff=0.0
    for n=1:length(cfs)
        diff=diff+cfs[n]*(DiracDelta(dp[n])-DiracDelta(dp[n+1]))
    end
    return diff
end

#Derivative Operator for HeavisideSpace
function Derivative(H::HeavisideSpace, k::Int)
    @assert k == 1
    ConcreteDerivative(H)
end

differentiate(f::Fun{SplineSpace{1,T,R}}) where {T,R} =
    Fun(HeavisideSpace(space(f).domain),
        diff(pad(f.coefficients,dimension(space(f))))./diff(space(f).domain.points))

integrate(f::Fun{HeavisideSpace{T,R}}) where {T,R} =
    Fun(SplineSpace{1,T,R}(space(f).domain),
        [0;cumsum(f.coefficients).*diff(space(f).domain.points)])


text:
a = Vec4f0(0)
b = Vec2f0(2)
c = Vec4f0(b, 1,2)
d = Vec{4, Int}(b, 1,2)
d = Vec{4, UInt}(b, 1,2)

m  = Mat{4,3,Float32}(9f0)
m2 = Mat{3,3,Float32}(9f0)
m2 = Mat{3,3,Float32}(ntuple(x->0f0, 3), ntuple(x->0f0, 3), ntuple(x->0f0, 3))

gluniform(Int32(1), a)
gluniform(Int32(1), [a,a])
gluniform(Int32(1), m)
gluniform(Int32(1), [m, m])


text:
module PDDL

using Base: @kwdef
using Julog

export Domain, Problem, Action, Event, State
export Term, Compound, Var, Const
export parse_domain, parse_problem, parse_pddl, pddl, pddl_str
export write_domain, write_problem, write_pddl
export load_domain, load_problem, preprocess
export save_domain, save_problem
export get_value
export GoalSpec
export has_term_in_state
export flatten_goal
export get_static_predicates, get_static_functions
export satisfy, evaluate, find_matches
export init_state, goal_state, initialize, transition, simulate
export get_preconditions, get_effect
export effect_diff, precond_diff, effect_dist, update!, update
export available, relevant, execute, execpar, execseq, regress, trigger
export use_available_action_cache!, use_relevant_action_cache!
export clear_available_action_cache!, clear_relevant_action_cache!
export compute_relaxed_actions
export regularize_clauses, compute_hsp_axioms, filter_negative_preconds
export compute_costs_one_step_derivation, compute_costs_one_step_derivation
export get_arbitrary_cost
export length_set, init_facts_costs, get_facts_and_state, convert_dict_any_to_term
export convert_dict_any_to_preconds, cache_global_preconditions

include("requirements.jl")
include("structs.jl")
include("parser.jl")
include("goals.jl")
include("writer.jl")

using .Parser, .Writer

include("utils.jl")
include("core.jl")
include("preprocess.jl")
include("states.jl")
include("effects.jl")
include("actions.jl")
include("events.jl")

end # module


text:
export ImportanceSamplingApproximation

using Random
using StatsBase

"""
    ImportanceSamplingApproximation

This structure stores all information needed to perform an importance sampling procedure and provides
convenient functions to generate samples and weights to approximate expectations

# Fields
- `rng`: random number generator objects
- `nsamples`: number of samples generated by default
"""
struct ImportanceSamplingApproximation{T, R}
    rng        :: R
    nsamples   :: Int
    bsamples   :: Vector{T}
    bweights   :: Vector{T}
    resampling :: Bool
    rsamples   :: Vector{T}
end

function ImportanceSamplingApproximation(rng::R, nsamples::Int; resampling::Bool = true) where { R }
    return ImportanceSamplingApproximation(Float64, rng, nsamples; resampling = resampling)
end

function ImportanceSamplingApproximation(::Type{T}, rng::R, nsamples::Int; resampling::Bool = true) where { T, R }
    return ImportanceSamplingApproximation{T, R}(rng, nsamples, Vector{T}(undef, nsamples), Vector{T}(undef, nsamples), resampling, Vector{T}(undef, nsamples))
end

getsamples(approximation::ImportanceSamplingApproximation, distribution)           = getsamples(approximation, distribution, approximation.nsamples)
getsamples(approximation::ImportanceSamplingApproximation, distribution, nsamples) = rand(approximation.rng, distribution, nsamples)

function approximate_meancov(approximation::ImportanceSamplingApproximation, g::Function, distribution)

    # We use preallocated arrays to sample and compute transformed samples and weightd
    rand!(approximation.rng, distribution, approximation.bsamples)
    map!(g, approximation.bweights, approximation.bsamples)

    if approximation.resampling
        n_eff = 1 / sum(Base.Generator(abs2, approximation.bweights))

        # Here we assume that lengths of bweights ans bsamples vectors are the same
        N = length(approximation.bweights)

        if n_eff < N / 10
            # We use rsamples as a temporary buffer here
            copyto!(approximation.rsamples, 1, approximation.bsamples, 1, N)
            # Here rsamples are equal to bsamples, but during sampling bsamples will be overwritten
            sample!(approximation.rng, approximation.rsamples, Weights(approximation.bweights), approximation.bsamples)
            fill!(approximation.bweights, 1 / N)
        end
    end

    normalization = sum(approximation.bweights)

    if iszero(normalization)
        return mean(distribution), var(distribution)
    end

    map!(Base.Fix2(/, normalization), approximation.bweights, approximation.bweights)

    m = mapreduce(prod, +, zip(approximation.bweights, approximation.bsamples))

    _v = let m = m
        (r) -> r[1] * (r[2] - m) ^ 2
    end

    v = mapreduce(_v, +, zip(approximation.bweights, approximation.bsamples))

    if isnan(m) || isnan(v) || isinf(m) || iszero(v)
        return mean(distribution), var(distribution)
    end

    return m, v

end