text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates, Pipe

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_VLAIO = CSV.read(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", DataFrame)
transform!(df_VLAIO, ["PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE", "PARTNER_ISO", "FLOW"] .=> ByRow(string), renamecols=false)
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> convert(Union{Missing, Float64}, x)), renamecols=false)

# remove missings
df_VLAIO.PRODUCT_NC = lpad.(string.(df_VLAIO.PRODUCT_NC), 8, '0') # needs to be done again
subset!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-16027346/21667553 ~26% observations

# MAD adjustment
outlier_cutoff = 3
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
gdf = groupby(df_VLAIO, "PRODUCT_NC")
df_VLAIO = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
subset!(df_VLAIO, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-13110297/16027346 ~18% observations (~40% of total)
df_VLAIO = df_VLAIO[:, Not([:MAD, :UNIT_PRICE])]

# add EU/WORLD again
EU27 = ["Roemenië", "Griekenland", "Oostenrijk", "Polen", "Duitsland", "Spanje", "Hongarije", "Slovakije", "Italië", "Nederland",
       "Frankrijk", "Letland", "Kroatië", "Cyprus", "Malta", "Litouwen", "Slovenië", "Estland", "Portugal", "Finland", "Tsjechië",
       "Luxemburg", "Zweden", "Denemarken", "Bulgarije", "Ierland", "België"]

# slightly modified functions
function append_WORLD(df::DataFrame)

    # sum over PARTNER
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_WORLD.TRADE_TYPE .= "total"
    df_WORLD.PARTNER_ISO .= "WORLD"

    df = vcat(df, df_WORLD)

    return df
end

function append_EU(df::DataFrame, EU::Vector{String})

    # subset and sum over EU ctrys
    df_EU = subset(df, :PARTNER_ISO => ByRow(x -> x in EU))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df_EU, cols_grouping)
    df_EU = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_EU.TRADE_TYPE .= "intra"
    df_EU.PARTNER_ISO .= "EU"

    df = vcat(df, df_EU)

    return df
end

df_VLAIO = append_WORLD(df_VLAIO)
df_VLAIO = append_EU(df_VLAIO, EU27)

# double check if no zeros/missing introduced by WORLD/EU
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false)
subset!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-15228494/15228494 ~0% observations

cols_subset = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC", "VALUE_IN_EUROS", "QUANTITY_IN_KG"]
df_VLAIO = df_VLAIO[:, cols_subset]
sort!(df_VLAIO, :PERIOD)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 1
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# # prepare data
# df = transform(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
# subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-18211714/24661479 ~26% observations

# # remove outliers based on medians abs dev.
# outlier_cutoff = 3
# transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
# gdf = groupby(df, "PRODUCT_NC")
# df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
# subset!(df, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-14861374/18211714 ~18% observations (~40% of total)

# produce data for plotting
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
df = data_fig1(df_VLAIO, ["Vlaanderen"], partners)

# -----------------------

# formatting
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> (x -> HP(x, λ)) .=> [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP])

# create unweighted unit prices
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
transform!(df, [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP] => ByRow((v,q) -> v/q) => :UNIT_PRICE_HP)

# create 3MMA as observations are very volatile
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :PRICE_INDEX => (x -> movingaverage(x,3)) => :PRICE_INDEX_3MMA)

# plotting
sort!(df, :DATE)

for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS/1e9,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros (billion)", title="Flemish "* flow*": values")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/" * "fig1_" * flow * "_values" * ".png") # export image dropbox

    # price index (weighted unit prices)
    # Notes:
    #   - unit prices are extremely volatile, 3 MMA would be useful
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted)", title="Flemish "*flow*": price index (weighted)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/" * "fig1_" * flow * "_prices_weighted" * ".png")

    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX_3MMA,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted, 3MMA)", title="Flemish "*flow*": price index (weighted, 3MMA)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/" * "fig1_" * flow * "_prices_weighted_3MMA" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :UNIT_PRICE,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (unweighted)", title="Flemish "*flow*": price index (unweighted)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/" * "fig1_" * flow * "_prices_unweighted" * ".png")

end

# plotting of HP
for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS_HP/1e9,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros (billion)", title="Flemish "*flow*": values (HP, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/HP/" * "fig1_" * flow * "_values" * "_HP" * ".png") # export image dropbox

    # price index (weighted unit prices)
    # Notes:
    #   - unit prices are extremely volatile, 3 MMA would be useful
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX_HP,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted)", title="Flemish "*flow*": price index (weighted, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/HP/" * "fig1_" * flow * "_prices_weighted" * "_HP" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :UNIT_PRICE_HP,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (unweighted)", title="Flemish "*flow*": price index (unweighted, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig1/HP/" * "fig1_" * flow * "_prices_unweighted" * "_HP" * ".png")

end

# VLA vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :UNIT_PRICE_HP,
        group=:FLOW, lw=2, legend=:topleft, ylabel="euros", title="Price index (unweighted, λ=$λ): \n Flanders vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/fig1/HP/" * "fig1_" * "VLA_GB" * "_prices_unweighted" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 4
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# rolling standard deviation
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX, :UNIT_PRICE, :VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP, :UNIT_PRICE_HP]
                .=> (x -> rolling_std(x, 6)) .=>
                [:STD_VALUE, :STD_QUANTITY, :STD_PRICE, :STD_UNIT, :STD_VALUE_HP, :STD_QUANTITY_HP, :STD_PRICE_HP, :STD_UNIT_HP])

# plotting HP
for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :STD_UNIT_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros", title="Flemish "* flow*": std unit prices (unweighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig4/HP/" * "fig4_" * flow * "_std" * "_prices_unweighted" * "_HP" * ".png") # export image dropbox

end

# VLA vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :STD_UNIT_HP,
        group=:FLOW, lw=2, legend=:topleft, ylabel="euros", title="6 months STD price index (unweighted, λ=$λ): \n Flanders vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/fig4/HP/" * "fig4_" * "VLA_GB" * "_STD_prices_unweighted" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 2
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# # prepare data
# df = transform(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
# subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-18211714/24661479 ~26% observations

# # remove outliers based on medians abs dev.
# outlier_cutoff = 3
# transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
# gdf = groupby(df, "PRODUCT_NC")
# df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
# subset!(df, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-14861374/18211714 ~18% observations (~40% of total)

# produce data for plotting
df = data_fig2(df_VLAIO, ["Vlaanderen"], ["Verenigd Koninkrijk"], 2)

# ------------------

# formatting
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, x -> nrow(x) < 3 ? DataFrame() : x) # remove groups with only two observations otherwise we cannot apply HP filter
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> (x -> HP(x, λ)) .=> [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP])

# create unweighted unit prices
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
transform!(df, [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP] => ByRow((v,q) -> v/q) => :UNIT_PRICE_HP)

# plotting
# subsetting products for figures
sort!(df) # for some reason needed to have correct DATE axis
products = 1:8:89

for flow in ["imports", "exports"]
    for i in products

    product_range = lpad.(i:i+6,2,'0')

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="euros", title="Flemish "*flow*": values")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/values/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_values" * ".png") # export image

    # price index (weighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :PRICE_INDEX,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (weighted)", title="Flemish "*flow*": unit prices (weighted)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/prices_weighted/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_prices_weighted" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :UNIT_PRICE,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (unweighted)", title="Flemish "*flow*": unit prices (unweighted)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/prices_unweighted/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_prices_unweighted" * ".png")

    end
end

for flow in ["imports", "exports"]
    for i in products

    product_range = lpad.(i:i+6,2,'0')

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS_HP/1e6,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="euros (million)", title="Flemish "*flow*": values (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/HP/values/" * flow * "/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_values" * "_HP" * ".png") # export image

    # price index (weighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :PRICE_INDEX_HP,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (weighted)", title="Flemish "*flow*": unit prices (weighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/HP/prices_weighted/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_prices_weighted" * "_HP" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :UNIT_PRICE_HP,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (unweighted)", title="Flemish "*flow*": unit prices (unweighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig2/HP/prices_unweighted/" * "fig2_" * flow * "_" * product_range[1] * "_" * product_range[end] * "_prices_unweighted" * "_HP" * ".png")

    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 3
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# # prepare data
# df = transform(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
# subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-18211714/24661479 ~26% observations

# # remove outliers based on medians abs dev.
# outlier_cutoff = 3
# transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
# gdf = groupby(df, "PRODUCT_NC")
# df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
# subset!(df, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-14861374/18211714 ~18% observations (~40% of total)

# produce data for plotting
partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
df = data_fig3(df_VLAIO, ["Vlaanderen"], partners)

# function to create data for figure 3
function data_fig3(df::DataFrame, declarants::Vector{String}, partners::Vector{String})

    # clean df
    df = subset(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL"), :DECLARANT_ISO => ByRow(x -> x in declarants),
                    :PARTNER_ISO => ByRow(x -> x in [partners; "WORLD"])) # take out TOTAL

    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

    df_join = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")),
                    subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD")), on=[:DECLARANT_ISO, :FLOW, :PERIOD], makeunique=true)

    transform!(df_join, [:VALUE_IN_EUROS, :VALUE_IN_EUROS_1] => ByRow((x,s) -> x/s*100) => :VALUE_SHARE)
    transform!(df_join, [:QUANTITY_IN_KG, :QUANTITY_IN_KG_1] => ByRow((x,s) -> x/s*100) => :QUANTITY_SHARE)

    cols_name = ["VALUE_PARTNER", "QUANTITY_PARTNER", "VALUE_WORLD", "QUANTITY_WORLD"]
    rename!(df_join, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :VALUE_IN_EUROS_1, :QUANTITY_IN_KG_1] .=> cols_name)

    return df_join
end

df = data_fig3(df_VLAIO, ["Vlaanderen"], partners)

# ------------------

# formatting
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
transform!(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> (x -> HP(x, λ)) .=> [:VALUE_SHARE_HP, :QUANTITY_SHARE_HP])

# plotting
sort!(df, :DATE)

for flow in ["imports", "exports"]

    # value share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :VALUE_SHARE,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": value share")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig3/" * "fig3_" * flow * "_value_share" * ".png") # export image

    # quantity share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :QUANTITY_SHARE,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="kg", title="Flemish "*flow*": quantity share")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig3/" * "fig3_" * flow * "_quantity_share" * ".png") # export image

end

# plotting HP
for flow in ["imports", "exports"]

    # value share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :VALUE_SHARE_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "*flow*": value share (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig3/HP/" * "fig3_" * flow * "_value_share" * "_HP" * ".png") # export image

    # quantity share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :QUANTITY_SHARE_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="kg", title="Flemish "*flow*": quantity share (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig3/HP/" * "fig3_" * flow * "_quantity_share" * "_HP" * ".png") # export image

end

# VLA vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :VALUE_SHARE_HP,
        group=:FLOW, lw=2, legend=:bottomleft, ylabel="percentages", title="Total trade share (values, λ=$λ): \n Flanders vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/fig3/HP/" * "fig3_" * "VLA_GB" * "_value_share" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Table 1
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
df = tab1(df_VLAIO, ["Vlaanderen"], partners, 2)

# aggregate to yearly data
transform!(df, :PERIOD => ByRow(x -> string(x)[1:4]) => :YEAR)
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

df_join = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")),
                    subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD")), on=[:DECLARANT_ISO, :FLOW, :PRODUCT_NC_digits, :YEAR], makeunique=true)

transform!(df_join, [:VALUE_IN_EUROS, :VALUE_IN_EUROS_1] => ByRow((x,s) -> x/s) => :SHARE_VALUE)
transform!(df_join, [:QUANTITY_IN_KG, :QUANTITY_IN_KG_1] => ByRow((x,s) -> x/s) => :SHARE_QUANTITY)

cols_name = ["VALUE_PARTNER", "QUANTITY_PARTNER", "VALUE_WORLD", "QUANTITY_WORLD"]
rename!(df_join, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :VALUE_IN_EUROS_1, :QUANTITY_IN_KG_1] .=> cols_name)
cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "VALUE_PARTNER", "QUANTITY_PARTNER", "VALUE_WORLD", "QUANTITY_WORLD", "SHARE_VALUE", "SHARE_QUANTITY"]
df = df_join[:, cols_name]

# ----------------

df_tab1 = subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"))
cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "SHARE_VALUE"]
df_tab1 = df_tab1[:, cols_name]

# # only keep top 10
# df_tab1 = @pipe df_tab1 |>
#     groupby(_, cols_grouping) |>
#     combine(_) do sdf
#         sorted = sort(sdf, order(:SHARE_VALUE, rev=true))
#         first(sorted, 10)
#     end

sort!(df_tab1)
transform!(df_tab1, :SHARE_VALUE => ByRow(x -> round(x*100, digits=2)), renamecols=false)
df_tab1_wide = unstack(df_tab1, :YEAR, :SHARE_VALUE)

XLSX.writetable(dir_dropbox * "results/images/VLAIO/tab/" * "table1_within_prod_importance" * ".xlsx", df_tab1_wide, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Table 2
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

partners = ["Verenigd Koninkrijk", "EU", "Duitsland", "Nederland", "Frankrijk"]
df = tab1(df_VLAIO, ["Vlaanderen"], partners, 2)

# aggregate to yearly data
transform!(df, :PERIOD => ByRow(x -> string(x)[1:4]) => :YEAR)
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum, renamecols=false)

# create total world exports/imports per year
df_WORLD = subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD"))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR"]
gdf = groupby(df_WORLD, cols_grouping)
df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> sum .=> [:VALUE_WORLD, :QUANTITY_WORLD])

df_join = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")), df_WORLD[:,Not(:PARTNER_ISO)], on=[:DECLARANT_ISO, :FLOW, :YEAR])

transform!(df_join, [:VALUE_IN_EUROS, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)
transform!(df_join, [:QUANTITY_IN_KG, :QUANTITY_WORLD] => ByRow((x,s) -> x/s) => :SHARE_QUANTITY)

# ----------------

df_tab2 = subset(df_join, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"))
cols_name = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits", "YEAR", "SHARE_VALUE"]
df_tab2 = df_tab2[:, cols_name]

# # only keep top 10
# df_tab2 = @pipe df_tab2 |>
#     groupby(_, cols_grouping) |>
#     combine(_) do sdf
#         sorted = sort(sdf, order(:SHARE_VALUE, rev=true))
#         first(sorted, 10)
#     end

sort!(df_tab2)
transform!(df_tab2, :SHARE_VALUE => ByRow(x -> round(x*100, digits=2)), renamecols=false)

df_tab2_wide = unstack(df_tab2, :YEAR, :SHARE_VALUE)

XLSX.writetable(dir_dropbox * "results/images/VLAIO/tab/" * "table2_overall_prod_importance" * ".xlsx", df_tab2_wide, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 5
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# function to create data for figure 5
function fig5(df::DataFrame, declarants::Vector{String}, partners::Vector{String})

    # clean df
    subset!(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :PARTNER_ISO => ByRow(x -> x in partners))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out TOTAL

    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, nrow => :COUNT)

    return df
end

partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk"]
df = fig5(df_VLAIO, ["Vlaanderen"], partners)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

for flow in ["imports", "exports"]

    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :COUNT,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="number of products", title="Flemish "* flow*": number of products")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/fig5/" * "fig5_" * flow * "_product_count" * ".png") # export image dropbox

end


text:
@doc raw"""
`ODE`: Ordinary Differential Equation

Defines an initial value problem
```math
\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,
```
with vector field ``v``, initial condition ``q_{0}`` and the solution
``q`` taking values in ``\mathbb{R}^{d}``.

### Fields

* `d`: dimension of dynamical variable ``q`` and the vector field ``v``
* `n`: number of initial conditions
* `v`: function computing the vector field
* `h`: function computing the Hamiltonian
* `t₀`: initial time
* `q₀`: initial condition

The function `v` providing the vector field must have the interface
```julia
    function v(t, q, v)
        v[1] = ...
        v[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, and
`v` is the vector which holds the result of evaluating the vector field ``v``
on `t` and `q`.

"""
struct ODE{dType <: Number, tType <: Number, vType <: Function,
           hType <: Union{Function,Nothing}, pType <: Union{Tuple,Nothing}, N} <: AbstractEquationODE{dType, tType}

    d::Int
    n::Int
    v::vType
    h::hType
    t₀::tType
    q₀::Array{dType,N}
    parameters::pType
    periodicity::Vector{dType}

    function ODE(DT::DataType, N::Int, d::Int, n::Int, v::vType, t₀::tType, q₀::AbstractArray{dType};
                 h::hType=nothing, parameters=nothing, periodicity=zeros(DT,d)) where {
                        dType <: Number, tType <: Number, vType <: Function, hType <: Union{Function,Nothing}}

        @assert d == size(q₀,1)
        @assert n == size(q₀,2)
        @assert ndims(q₀) == N ∈ (1,2)

        new{DT, tType, vType, hType, typeof(parameters), N}(d, n, v, h, t₀,
                convert(Array{DT}, q₀), parameters, periodicity)
    end
end

function ODE(v, t₀, q₀::AbstractArray{DT}; kwargs...) where {DT}
    ODE(DT, ndims(q₀), size(q₀,1), size(q₀,2), v, t₀, q₀; kwargs...)
end

function ODE(v, q₀; kwargs...)
    ODE(v, zero(eltype(q₀)), q₀; kwargs...)
end

Base.hash(ode::ODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.v, hash(ode.t₀,
        hash(ode.q₀, hash(ode.periodicity, hash(ode.parameters, h)))))))

Base.:(==)(ode1::ODE, ode2::ODE) = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.v == ode2.v
                             && ode1.h == ode2.h
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.parameters == ode2.parameters
                             && ode1.periodicity == ode2.periodicity)

function Base.similar(ode::ODE, q₀; kwargs...)
    similar(ode, ode.t₀, q₀; kwargs...)
end

function Base.similar(ode::ODE, t₀::TT, q₀::AbstractArray{DT};
                      h=ode.h, parameters=ode.parameters, periodicity=ode.periodicity) where {DT  <: Number, TT <: Number}
    @assert ode.d == size(q₀,1)
    ODE(ode.v, t₀, q₀; h=h, parameters=parameters, periodicity=periodicity)
end

@inline Base.ndims(ode::ODE) = ode.d

@inline periodicity(equation::ODE) = equation.periodicity

ext: jl


text:
using ITensors
using ITensorInfiniteMPS
using ITensorInfiniteMPS.HDF5
using Test

@testset "HDF5 Read and Write" begin
  function space_shifted(::Model"ising", q̃sz)
    return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
  end

  N = 2
  model = Model"ising"()
  space_ = fill(space_shifted(model, 0), N)
  s = infsiteinds("S=1/2", N; space=space_)
  initstate(n) = "↑"
  ψ = InfMPS(s, initstate)

  @testset "InfiniteCanonicalMPS" begin
    fo = h5open("data.h5", "w")
    write(fo, "ψ", ψ)
    close(fo)

    fi = h5open("data.h5", "r")
    ψr = read(fi, "ψ", InfiniteCanonicalMPS)
    close(fi)

    @test ψ.AL.data == ψr.AL.data
    @test ψ.C.data == ψr.C.data
    @test ψ.AR.data == ψr.AR.data
  end

  # Clean up the test hdf5 file
  rm("data.h5"; force=true)
end


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

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "GAP_lib"
version = v"400.1100.0"
upstream_version = v"4.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/gap-system/gap.git", "7a2dfbec79d3d20a4ac4d5f174eaf35aa34387d7"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(upstream_version)/gap-$(upstream_version)-core.tar.bz2",
#                  "6637f66409bc91af21eaa38368153270b71b13b55b75cc1550ed867c629901d1"),
    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(upstream_version)/packages-required-v$(upstream_version).tar.gz",
                  "29ab6e2752f39d22e3f0a19e5bcbfec1710993b4cfd52337a0ca10fda6a76537";
                  unpack_target="pkg"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

mv ../pkg .
find pkg -name '._*' -exec rm \{\} \; # unwanted files

# run autogen.sh if compiling from it source and/or if configure was patched
./autogen.sh

# compile a native version of GAP so we can use it to generate the manual
# (the manual is only in FULL gap release tarballs, not in the -core tarball
# nor in git snapshots)
mkdir native-build
cd native-build
#apk add gmp-dev zlib-dev
../configure --build=${MACHTYPE} --host=${MACHTYPE} CC=${CC_BUILD} CXX=${CXX_BUILD} --with-zlib=${prefix} --with-gmp=${prefix}
make -j${nproc}
make html   # build the manual (only HTML and txt; for PDF we'd need LaTeX)
cd ..

# remove the native build, it has done its job
rm -rf native-build

# the license
install_license LICENSE

# "install" most of the files
rm -rf autom4te.cache dev etc extern hpcgap/extern pkg
mkdir -p ${prefix}/share/gap/
mv * ${prefix}/share/gap/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = Product[
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Zlib_jll"),
    BuildDependency("GMP_jll"),
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
This example was adapted from Daniel Shiffman's "Snake Game Redux" video,
from 'The Coding Train' Youtube channel.

https://www.youtube.com/watch?v=OMoVcohRgZA&t=1658s
"""
module SnakeGame

using Alexya

@init "Snake Game" 800 600

include("snake.jl")
include("utils.jl")
include("draw.jl")
include("controls.jl")

start()

end


text:
using DelaySSAdocs
using Documenter

DocMeta.setdocmeta!(DelaySSAdocs, :DocTestSetup, :(using DelaySSAdocs); recursive=true)

makedocs(;
    modules=[DelaySSAdocs],
    authors="Xiaoming Fu",
    repo="https://github.com/palmtree2013/DelaySSAdocs.jl/blob/{commit}{path}#{line}",
    sitename="DelaySSAdocs.jl",
    format=Documenter.HTML(;
        mathengine=Documenter.Writers.HTMLWriter.MathJax2(),
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://palmtree2013.github.io/DelaySSAdocs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => [
            "tutorials/tutorials.md",
            "tutorials/bursty.md",
            "tutorials/delay_degradation.md",
            "tutorials/heterogeneous_delay.md",
            "tutorials/delay_oscillator.md",
            "tutorials/stochastic_delay.md",
        ],
        "Algorithm" => [
            "algorithms/notations.md",
            "algorithms/delayrejection.md",
            "algorithms/delaydirect.md",
            "algorithms/delaymnrm.md",
        ],
        "Theory" => "theory.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/palmtree2013/DelaySSAdocs.jl",
    devbranch="main",
)


text:
function compute()

	a::BigInt=sum(NUMBERS)
	ans=string(a)[1:10]
	return ans
end

NUMBERS = BigInt[
	37107287533902102798797998220837590246510135740250,
	46376937677490009712648124896970078050417018260538,
	74324986199524741059474233309513058123726617309629,
	91942213363574161572522430563301811072406154908250,
	23067588207539346171171980310421047513778063246676,
	89261670696623633820136378418383684178734361726757,
	28112879812849979408065481931592621691275889832738,
	44274228917432520321923589422876796487670272189318,
	47451445736001306439091167216856844588711603153276,
	70386486105843025439939619828917593665686757934951,
	62176457141856560629502157223196586755079324193331,
	64906352462741904929101432445813822663347944758178,
	92575867718337217661963751590579239728245598838407,
	58203565325359399008402633568948830189458628227828,
	80181199384826282014278194139940567587151170094390,
	35398664372827112653829987240784473053190104293586,
	86515506006295864861532075273371959191420517255829,
	71693888707715466499115593487603532921714970056938,
	54370070576826684624621495650076471787294438377604,
	53282654108756828443191190634694037855217779295145,
	36123272525000296071075082563815656710885258350721,
	45876576172410976447339110607218265236877223636045,
	17423706905851860660448207621209813287860733969412,
	81142660418086830619328460811191061556940512689692,
	51934325451728388641918047049293215058642563049483,
	62467221648435076201727918039944693004732956340691,
	15732444386908125794514089057706229429197107928209,
	55037687525678773091862540744969844508330393682126,
	18336384825330154686196124348767681297534375946515,
	80386287592878490201521685554828717201219257766954,
	78182833757993103614740356856449095527097864797581,
	16726320100436897842553539920931837441497806860984,
	48403098129077791799088218795327364475675590848030,
	87086987551392711854517078544161852424320693150332,
	59959406895756536782107074926966537676326235447210,
	69793950679652694742597709739166693763042633987085,
	41052684708299085211399427365734116182760315001271,
	65378607361501080857009149939512557028198746004375,
	35829035317434717326932123578154982629742552737307,
	94953759765105305946966067683156574377167401875275,
	88902802571733229619176668713819931811048770190271,
	25267680276078003013678680992525463401061632866526,
	36270218540497705585629946580636237993140746255962,
	24074486908231174977792365466257246923322810917141,
	91430288197103288597806669760892938638285025333403,
	34413065578016127815921815005561868836468420090470,
	23053081172816430487623791969842487255036638784583,
	11487696932154902810424020138335124462181441773470,
	63783299490636259666498587618221225225512486764533,
	67720186971698544312419572409913959008952310058822,
	95548255300263520781532296796249481641953868218774,
	76085327132285723110424803456124867697064507995236,
	37774242535411291684276865538926205024910326572967,
	23701913275725675285653248258265463092207058596522,
	29798860272258331913126375147341994889534765745501,
	18495701454879288984856827726077713721403798879715,
	38298203783031473527721580348144513491373226651381,
	34829543829199918180278916522431027392251122869539,
	40957953066405232632538044100059654939159879593635,
	29746152185502371307642255121183693803580388584903,
	41698116222072977186158236678424689157993532961922,
	62467957194401269043877107275048102390895523597457,
	23189706772547915061505504953922979530901129967519,
	86188088225875314529584099251203829009407770775672,
	11306739708304724483816533873502340845647058077308,
	82959174767140363198008187129011875491310547126581,
	97623331044818386269515456334926366572897563400500,
	42846280183517070527831839425882145521227251250327,
	55121603546981200581762165212827652751691296897789,
	32238195734329339946437501907836945765883352399886,
	75506164965184775180738168837861091527357929701337,
	62177842752192623401942399639168044983993173312731,
	32924185707147349566916674687634660915035914677504,
	99518671430235219628894890102423325116913619626622,
	73267460800591547471830798392868535206946944540724,
	76841822524674417161514036427982273348055556214818,
	97142617910342598647204516893989422179826088076852,
	87783646182799346313767754307809363333018982642090,
	10848802521674670883215120185883543223812876952786,
	71329612474782464538636993009049310363619763878039,
	62184073572399794223406235393808339651327408011116,
	66627891981488087797941876876144230030984490851411,
	60661826293682836764744779239180335110989069790714,
	85786944089552990653640447425576083659976645795096,
	66024396409905389607120198219976047599490197230297,
	64913982680032973156037120041377903785566085089252,
	16730939319872750275468906903707539413042652315011,
	94809377245048795150954100921645863754710598436791,
	78639167021187492431995700641917969777599028300699,
	15368713711936614952811305876380278410754449733078,
	40789923115535562561142322423255033685442488917353,
	44889911501440648020369068063960672322193204149535,
	41503128880339536053299340368006977710650566631954,
	81234880673210146739058568557934581403627822703280,
	82616570773948327592232845941706525094512325230608,
	22918802058777319719839450180888072429661980811197,
	77158542502016545090413245809786882778948721859617,
	72107838435069186155435662884062257473692284509516,
	20849603980134001723930671666823555245252804609722,
	53503534226472524250874054075591789781264330331690,
]

println(compute())


text:
using BrainFlow

# enable all possible logs from all three libs
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)
BrainFlow.enable_dev_logger(BrainFlow.DATA_HANDLER)
BrainFlow.enable_dev_logger(BrainFlow.ML_MODULE)

params1 = BrainFlowInputParams()
params1.other_info = "board1"
board_shim1 = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params1)
params2 = BrainFlowInputParams()
params2.other_info = "board2"
board_shim2 = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params2)

BrainFlow.prepare_session(board_shim1)
BrainFlow.prepare_session(board_shim2)
BrainFlow.start_stream(board_shim1)
sleep(2)
data = BrainFlow.get_current_board_data(10, board_shim1)
println(BrainFlow.calc_stddev(data[2, :]))

BrainFlow.release_all_sessions()

# calc concentration
model_params = BrainFlowModelParams(metric = "concentration", classifier = "KNN")
BrainFlow.prepare(model_params)
BrainFlow.release_all()

max_issues_repo_name: HDPNC/brainflow

max_forks_repo_name: HDPNC/brainflow


text:
@doc raw"""
    inverse_retract(M, p, q method::NLSolveInverseRetraction; kwargs...)

Approximate the inverse of the retraction specified by `method.retraction` from `p` with
respect to `q` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M` using NLsolve. This inverse retraction is
not guaranteed to succeed and probably will not unless `q` is close to `p` and the initial
guess `X0` is close.

If the solver fails to converge, an [`OutOfInjectivityRadiusError`](@ref) is raised.
See [`NLSolveInverseRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.NLSolveInverseRetraction) for configurable parameters.
"""
inverse_retract(::AbstractManifold, p, q, ::NLSolveInverseRetraction; kwargs...)

function inverse_retract_nlsolve!(
    M::AbstractManifold,
    X,
    p,
    q,
    m::NLSolveInverseRetraction;
    kwargs...,
)
    X0 = m.X0 === nothing ? zero_vector(M, p) : m.X0
    res = _inverse_retract_nlsolve(M, p, q, m; kwargs...)
    return copyto!(X, res.zero)
end

function _inverse_retract_nlsolve(
    M::AbstractManifold,
    p,
    q,
    m::NLSolveInverseRetraction;
    kwargs...,
)
    X0 = m.X0 === nothing ? zero_vector(M, p) : m.X0
    function f!(F, X)
        m.project_tangent && project!(M, X, p, X)
        retract!(M, F, p, project(M, p, X), m.retraction; kwargs...)
        m.project_point && project!(M, q, q)
        F .-= q
        return F
    end
    isdefined(Manifolds, :NLsolve) ||
        @warn "To use NLSolveInverseRetraction, NLsolve must be loaded using `using NLsolve`."
    res = NLsolve.nlsolve(f!, X0; m.nlsolve_kwargs...)
    if !res.f_converged
        @debug res
        throw(OutOfInjectivityRadiusError())
    end
    return res
end


text:
# Convenience functions that used to be in Julia, moved to LegacyStrings.jl
# Copyright 2017 Scott P. Jones and other contributors to the Julia language
# Licensed under MIT License, see LICENSE.md

# Renamed from ascii to to_ascii, to prevent issues with Base.ascii
to_ascii(str) = convert(ASCIIStr, str)
to_ascii(pnt::Ptr{UInt8}) = convert(ASCIIStr, unsafe_string(p))
to_ascii(pnt::Ptr{UInt8}, len::Integer) = convert(ASCIIStr, unsafe_string(p, len))

"""
    utf8(s)

Create a UTF-8 string from a byte array, array of `UInt8`, or any other string type. (Data
must be valid UTF-8. Conversions of byte arrays check for a byte-order marker in the first
two bytes, and do not include it in the resulting string.)
"""
utf8(x) = convert(UTF8Str, x)

"""
    utf8(::Union{Ptr{UInt8}, Ptr{Int8}} [, length])

Create a string from the address of a NUL-terminated UTF-8 string. A copy is made; the
pointer can be safely freed. If `length` is specified, the string does not have to be
NUL-terminated.
"""
utf8(p::Ptr{UInt8}) = utf8(unsafe_string(p))
utf8(p::Ptr{UInt8}, len::Integer) = utf8(unsafe_string(p, len))

"""
    utf16(s)

Create a UTF-16 string from a byte array, array of `UInt16`, or any other string type. (Data
must be valid UTF-16. Conversions of byte arrays check for a byte-order marker in the first
two bytes, and do not include it in the resulting string.)
"""
utf16(x) = convert(UTF16Str, x)

"""
    utf16(::Union{Ptr{UInt16}, Ptr{Int16}} [, length])

Create a string from the address of a NUL-terminated UTF-16 string. A copy is made; the
pointer can be safely freed. If `length` is specified, the string does not have to be
NUL-terminated.
"""
function utf16(pnt::Ptr{UInt16}, len::Unsigned)
    buf, out = _allocate(UInt16, len)
    _memcpy(out, pnt, len)
    Str(UTF16CSE, buf)
end
function utf16(pnt::Ptr{UInt16})
    len = 0
    while (ch = unsafe_load(pnt, len += 1)) != 0
        check_valid(ch, len)
    end
    utf16(pnt, len%UInt)
end
utf16(pnt, len::Signed) = len < 0 ? error("negative length") : utf16(pnt, len%Unsigned)

"""
    utf32(s)

Create a UTF-32 string from a byte array, array of `Char` or `UInt32`, or any other string
type. (Conversions of byte arrays check for a byte-order marker in the first four bytes, and
do not include it in the resulting string.)
"""
utf32(x) = convert(UTF32Str, x)

# These are broken by #24999

#utf32(p::Union{Ptr{Char}, Ptr{Int32}}, len::Integer) = utf32(reinterpret(Ptr{UInt32}, p), len)
#utf32(p::Union{Ptr{Char}, Ptr{Int32}}) = utf32(reinterpret(Ptr{UInt32}, p))

"""
    utf32(::Ptr{UInt32}, [, length])

Create a string from the address of a NUL-terminated UTF-32 string. A copy is made; the
pointer can be safely freed. If `length` is specified, the string does not have to be
NUL-terminated.
"""
function utf32(pnt::Ptr{UInt32}, len::Unsigned)
    buf, out = _allocate(UInt32, len)
    _memcpy(out, pnt, len)
    Str(UTF32CSE, buf)
end
function utf32(pnt::Ptr{UInt32})
    len = 0
    while (ch = unsafe_load(pnt, len += 1)) != 0
        check_valid(ch, len)
    end
    utf32(pnt, len%UInt)
end
utf32(pnt, len::Signed) = len < 0 ? error("negative length") : utf32(pnt, len%Unsigned)


text:
# see: A Gretton, e.g. http://www.gatsby.ucl.ac.uk/~gretton/coursefiles/lecture4_introToRKHS.pdf

export
  mmd!,  # KED
  mmd

"""
    $SIGNATURES

Normal kernel used for Hilbert space embeddings.
"""
ker(M::MB.AbstractManifold, p, q, sigma::Real=0.001) = exp( -sigma*(distance(M, p, q)^2) )

# overwrite non-symmetric with alternate implementations
# ker(M::MB.AbstractManifold, p, q, sigma::Real=0.001) = exp( -sigma*(distance(M, p, q)^2) )

"""
    $SIGNATURES

MMD disparity (i.e. 'distance') measure based on Kernel Hilbert Embeddings between two beliefs.

Notes:
- This is the in-place version (well more in-place than mmd)

DevNotes:
- TODO dont assume equally weighted particles
- TODO profile SIMD vs SLEEF
- TODO optimize memory
- TODO make multithreaded

See also: [`mmd`](@ref), [`ker`](@ref)
"""
function mmd!(MF::MB.AbstractManifold,
              val::AbstractVector{<:Real},
              a::AbstractVector,
              b::AbstractVector,
              N::Integer=length(a), M::Integer=length(b);
              bw::AbstractVector{<:Real}=[0.001;] )
  #
  # TODO allow unequal data too
  _N = 1.0/N
  _M = 1.0/M
  _val1 = 0.0
  _val2 = 0.0
  _val3 = 0.0
  @inbounds @fastmath for i in 1:N
    @simd for j in 1:M
      _val1 -= ker(MF, a[i], b[j], bw[1])
    end
  end
  _val1 *= 2.0*_N*_M
  @inbounds @fastmath for i in 1:N
    @simd for j in 1:N
      _val2 += ker(MF, a[i], a[j], bw[1])
    end
  end
  _val2 *= (_N*_N)
  @inbounds @fastmath for i in 1:M
    @simd for j in 1:M
      _val3 += ker(MF, b[i], b[j], bw[1])
    end
  end
  _val3 *= (_M*_M)

  # accumulate all terms
  val[1] = _val1 + _val2 + _val3
  return val
end

"""
    $SIGNATURES

MMD disparity (i.e. 'distance') measure based on Kernel Hilbert Embeddings between two beliefs.

Notes:
- This is a wrapper to the in-place `mmd!` function.

Related

mmd!, ker
"""
function mmd( MF::MB.AbstractManifold,
              a::AbstractVector,
              b::AbstractVector,
              N::Int=length(a), M::Int=length(b);
              bw::AbstractVector{<:Real}=[0.001;])
  #
  val = [0.0;]
  mmd!( MF, val,
        a,b,
        N, M; bw=bw )
  #
  return val[1]
end

function mmd(a::ManifoldKernelDensity{M}, b::ManifoldKernelDensity{M}; bw::Vector{<:Real}=[0.001;]) where M <: MB.AbstractManifold
  # @assert a.manifold == b.manifold "Manifolds not the same $(a.manifold), $(b.manifold)"
  aPts = getPoints(a)
  bPts = getPoints(b)
  mmd(a.manifold, aPts, bPts, bw=bw)
end

function isapprox(a::ManifoldKernelDensity, b::ManifoldKernelDensity; atol::Real=0.1)
  mmd(a,b) < atol
end


text:
module Random

import GreyDecision.GreyNumbers: GreyNumber

function Base.rand(::Type{GreyNumber})::GreyNumber{Float64}
    return GreyNumber(rand(), rand())
end

function Base.rand(::Type{GreyNumber{T}})::GreyNumber{T} where {T <: Real}
    return GreyNumber(rand(T), rand(T))
end

function Base.rand(::Type{GreyNumber{T}}, n::Int)::Array{GreyNumber{T}, 1} where {T <: Real}
    return map(x -> rand(GreyNumber{T}), 1:n)
end

function Base.rand(::Type{GreyNumber}, n::Int)::Array{GreyNumber{Float64}, 1}
    return rand(GreyNumber{Float64}, n)
end

function Base.rand(::Type{GreyNumber{T}}, n::Int, m::Int)::Array{GreyNumber{T}, 2} where {T <: Real}
    mat = Array{GreyNumber{T}, 2}(undef, n, m)
    for i in 1:n
        for j in 1:m
            mat[i, j] = rand(GreyNumber{T})
        end
    end
    return mat
end

function Base.rand(::Type{GreyNumber}, n::Int, m::Int)::Array{GreyNumber{Float64}, 2}
    return rand(GreyNumber{Float64}, n, m)
end

end # End of module Random2


text:
#=

  The BUGS Book

  Example 3.3.3 Trihalomethanes in tap water

  Cf ~/jags/bugs_book_3_3_3.R
  """
  model {
    for (i in 1:n) {
      y[i] ~ dnorm(mu, inv.sigma.squared)
    }
    mu ~ dnorm(gamma, inv.omega.squared)
    inv.omega.squared <- n0/sigma.squared
    inv.sigma.squared <- 1/sigma.squared
    y.pred ~ dnorm(mu, inv.sigma.squared)
    P.crit <- step(y.pred - y.crit)
  }

  Data: list(n=2, y=c(128, 132), gamma=120, n0=0.25,sigma.squared=25, y.crit=145)

  Output:
            Mean      SD  Naive SE Time-series SE
P.crit 3.213e-03 0.05659 0.0002001      0.0002034
mu     1.289e+02 3.33096 0.0117767      0.0117332
y.pred 1.289e+02 5.95716 0.0210618      0.0209726
  """

  Summary Statistics
  parameters       mean       std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol    Float64   Float64    Float64   Float64      Float64   Float64       Float64

          mu   128.9489    3.3703     0.0337    0.0502    4816.4635    0.9999     2519.0709
      y_pred   128.9964    6.0306     0.0603    0.0779    7496.1334    1.0000     3920.5719
      p_crit     0.0046    0.0677     0.0007    0.0006   10096.6867    0.9999     5280.6939

  Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5%
      Symbol    Float64    Float64    Float64    Float64    Float64

          mu   122.2672   126.7344   128.9249   131.2583   135.4088
      y_pred   117.2568   124.8223   128.9742   133.1282   140.7995
      p_crit     0.0000     0.0000     0.0000     0.0000     0.0000

  Cf ~/webppl/bugs_book_3_3_3.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_3_3_3(y=[128.0,132.0])
    n = 2
    gamma = 120
    n0 = 0.25
    sigma_squared = 25
    y_crit = 145.0

    omega_squared = sigma_squared/n0
    mu ~ Normal(gamma, sqrt(omega_squared))
    #  random Real inv_omega_squared ~ n0/sigma_squared;
    y ~ filldist(Normal(mu, sqrt(sigma_squared)),2)

    #  random Real inv_sigma_squared ~ 1/sigma_squared;
    y_pred ~ Normal(mu, sqrt(sigma_squared));
    p_crit ~ Dirac(y_pred > y_crit)

 end

model = bugs_book_3_3_3()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(10_000), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,5), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using BenchmarkTools

lines = readlines("3/input.txt");

function find_trees(input, right, down)
    grid = map(==('#'), hcat(collect.(input)...))';
    nrows, row_length = size(grid);

    # fold (reduce) function (result, element)
    # init res = (pos = 1, sum = 0)
    function fld_fnc((pos, sum), row_idx)
        # return our next pos and current sum
        return (mod1(pos + right, row_length), sum + grid[row_idx, pos])
    end

    # fold over each row, skip by down amount
    return foldl(fld_fnc, 1:down:nrows, init=(1, 0))[2]
end

function part1(input)
    return find_trees(input, 3, 1);
end

function part2(input)
    slopes = [
        (1, 1)
        (3, 1)
        (5, 1)
        (7, 1)
        (1, 2)
    ]
    return prod(map(x -> find_trees(input, x[1], x[2]), slopes))
end

println("Part 1: $(part1(lines))");
println("Part 2: $(part2(lines))");

@benchmark part2(lines)


text:
module AnalysisWrapper

using JavaCall

const BivariateFunction = @jimport org.hipparchus.analysis.BivariateFunction
const Class = @jimport java.lang.Class
const Field = @jimport org.hipparchus.Field
const FieldBivariateFunction = @jimport org.hipparchus.analysis.FieldBivariateFunction
const FieldUnivariateFunction = @jimport org.hipparchus.analysis.FieldUnivariateFunction
const FieldUnivariateMatrixFunction = @jimport org.hipparchus.analysis.FieldUnivariateMatrixFunction
const FieldUnivariateVectorFunction = @jimport org.hipparchus.analysis.FieldUnivariateVectorFunction
const FunctionUtils = @jimport org.hipparchus.analysis.FunctionUtils
const JString = @jimport java.lang.String
const MultivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.MultivariateDifferentiableFunction
const MultivariateFunction = @jimport org.hipparchus.analysis.MultivariateFunction
const MultivariateMatrixFunction = @jimport org.hipparchus.analysis.MultivariateMatrixFunction
const MultivariateVectorFunction = @jimport org.hipparchus.analysis.MultivariateVectorFunction
const Object = @jimport java.lang.Object
const ParametricUnivariateFunction = @jimport org.hipparchus.analysis.ParametricUnivariateFunction
const RealFieldBivariateFunction = @jimport org.hipparchus.analysis.RealFieldBivariateFunction
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RealFieldUnivariateFunction = @jimport org.hipparchus.analysis.RealFieldUnivariateFunction
const RealFieldUnivariateMatrixFunction = @jimport org.hipparchus.analysis.RealFieldUnivariateMatrixFunction
const RealFieldUnivariateVectorFunction = @jimport org.hipparchus.analysis.RealFieldUnivariateVectorFunction
const TrivariateFunction = @jimport org.hipparchus.analysis.TrivariateFunction
const UnivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction
const UnivariateMatrixFunction = @jimport org.hipparchus.analysis.UnivariateMatrixFunction
const UnivariateVectorFunction = @jimport org.hipparchus.analysis.UnivariateVectorFunction

export BivariateFunction
export FieldBivariateFunction
export FieldUnivariateFunction
export FieldUnivariateMatrixFunction
export FieldUnivariateVectorFunction
export FunctionUtils
export MultivariateFunction
export MultivariateMatrixFunction
export MultivariateVectorFunction
export ParametricUnivariateFunction
export RealFieldBivariateFunction
export RealFieldUnivariateFunction
export RealFieldUnivariateMatrixFunction
export RealFieldUnivariateVectorFunction
export TrivariateFunction
export UnivariateFunction
export UnivariateMatrixFunction
export UnivariateVectorFunction
export add
export collector
export combine
export compose
export derivative
export equals
export fix1st_argument
export fix2nd_argument
export get_class
export gradient
export hash_code
export multiply
export notify
export notify_all
export sample
export to_differentiable
export to_real_field_bivariate_function
export to_real_field_univariate_function
export to_real_field_univariate_matrix_function
export to_real_field_univariate_vector_function
export to_string
export value
export wait

include("DifferentiationWrapper/DifferentiationWrapper.jl")
include("FunctionWrapper/FunctionWrapper.jl")
include("IntegrationWrapper/IntegrationWrapper.jl")
include("InterpolationWrapper/InterpolationWrapper.jl")
include("PolynomialsWrapper/PolynomialsWrapper.jl")
include("SolversWrapper/SolversWrapper.jl")
include("bivariate_function.jl")
include("field_bivariate_function.jl")
include("field_univariate_function.jl")
include("field_univariate_matrix_function.jl")
include("field_univariate_vector_function.jl")
include("function_utils.jl")
include("multivariate_function.jl")
include("multivariate_matrix_function.jl")
include("multivariate_vector_function.jl")
include("parametric_univariate_function.jl")
include("real_field_bivariate_function.jl")
include("real_field_univariate_function.jl")
include("real_field_univariate_matrix_function.jl")
include("real_field_univariate_vector_function.jl")
include("trivariate_function.jl")
include("univariate_function.jl")
include("univariate_matrix_function.jl")
include("univariate_vector_function.jl")

end


text:
module LightBayes

using Distributions
using ArrayViews

import Base: length, mean, var, rand, rand!
import Base.LinAlg: axpy!
import Distributions: suffstats

export
    LikelihoodModel,
    IsoGaussModel,

    withparams,
    logpar,
    posterior,
    posterior_mode,
    posterior_mode!,
    posterior_rand,
    posterior_rand!

## source files

include("common.jl")
include("gauss.jl")

end # module


text:
include("../src/prrqr2.jl")
using BenchmarkTools
using LinearAlgebra
using LowRankApprox
using Random

Random.seed!(0)

tol = 1e-3

n = 1000; k = 10; m = 1000
s = zeros(k)
s[1:k] = @. exp(-(1:k));
U = randn(n,k)
V = randn(m,k)
A = U*diagm(s)*V';

#println("Benchmarking rank-revealing QR...")
#@btime _ = prrqr!(randn(n,100), tol; reltol=true)

println("Testing custom implementation...")
X = copy(A)
Q, R = _compress_block!(X; tol, reltol=false)
X = copy(A)
@time Q, R = _compress_block!(X; tol, reltol=false)
println("error in the truncated QR decomposition: ", norm(A - Q*R)/norm(A))

println("Testing LowRankApprox.jl implementation...")
# LowRankApprox.jl implementation
X = copy(A)
F = pqrfact(X; atol = tol)
X = copy(A)
@time F = pqrfact(X; atol = tol)
println("error in the LowRankApprox.jl PQR decomposition: ", norm(A - F.Q*F.R)/norm(A))

#_, σ, _ = svd(A);
#plot(1:n, σ, yaxis=:log, label="singular values")


text:
# SCIP enums
typealias SCIP_BIVAR_CONVEXITY Int8
const SCIP_BIVAR_ALLCONVEX = int8(0)
const SCIP_BIVAR_1CONVEX_INDEFINITE = int8(1)
const SCIP_BIVAR_CONVEX_CONCAVE = int8(2)
const SCIP_BIVAR_UNKNOWN = int8(3)

typealias SCIP_LinearConsType Int8
const SCIP_LINEARCONSTYPE_INVALIDCONS = int8(-1)
const SCIP_LINEARCONSTYPE_LINEAR = int8(0)
const SCIP_LINEARCONSTYPE_LOGICOR = int8(1)
const SCIP_LINEARCONSTYPE_KNAPSACK = int8(2)
const SCIP_LINEARCONSTYPE_SETPPC = int8(3)

typealias SCIP_SetppcType Int8
const SCIP_SETPPCTYPE_PARTITIONING = int8(0)
const SCIP_SETPPCTYPE_PACKING = int8(1)
const SCIP_SETPPCTYPE_COVERING = int8(2)

typealias SCIP_ClockType Int8
const SCIP_CLOCKTYPE_DEFAULT = int8(0)
const SCIP_CLOCKTYPE_CPU = int8(1)
const SCIP_CLOCKTYPE_WALL = int8(2)

typealias SCIP_DispStatus Int8
const SCIP_DISPSTATUS_OFF = int8(0)
const SCIP_DISPSTATUS_AUTO = int8(1)
const SCIP_DISPSTATUS_ON = int8(2)

typealias SCIP_ExprOp Int8
const SCIP_EXPR_VARIDX = int8(1)
const SCIP_EXPR_CONST = int8(2)
const SCIP_EXPR_PARAM = int8(3)
const SCIP_EXPR_PLUS = int8(8)
const SCIP_EXPR_MINUS = int8(9)
const SCIP_EXPR_MUL = int8(10)
const SCIP_EXPR_DIV = int8(11)
const SCIP_EXPR_SQUARE = int8(12)
const SCIP_EXPR_SQRT = int8(13)
const SCIP_EXPR_REALPOWER = int8(14)
const SCIP_EXPR_INTPOWER = int8(15)
const SCIP_EXPR_SIGNPOWER = int8(16)
const SCIP_EXPR_EXP = int8(17)
const SCIP_EXPR_LOG = int8(18)
const SCIP_EXPR_SIN = int8(19)
const SCIP_EXPR_COS = int8(20)
const SCIP_EXPR_TAN = int8(21)
const SCIP_EXPR_MIN = int8(24)
const SCIP_EXPR_MAX = int8(25)
const SCIP_EXPR_ABS = int8(26)
const SCIP_EXPR_SIGN = int8(27)
const SCIP_EXPR_SUM = int8(64)
const SCIP_EXPR_PRODUCT = int8(65)
const SCIP_EXPR_LINEAR = int8(66)
const SCIP_EXPR_QUADRATIC = int8(67)
const SCIP_EXPR_POLYNOMIAL = int8(68)
const SCIP_EXPR_LAST = int8(69)

typealias SCIP_ExprCurv Int8
const SCIP_EXPRCURV_UNKNOWN = int8(0)
const SCIP_EXPRCURV_CONVEX = int8(1)
const SCIP_EXPRCURV_CONCAVE = int8(2)
const SCIP_EXPRCURV_LINEAR = int8(SCIP_EXPRCURV_CONVEX | SCIP_EXPRCURV_CONCAVE)

typealias SCIP_BranchDir Int8
const SCIP_BRANCHDIR_DOWNWARDS = int8(0)
const SCIP_BRANCHDIR_UPWARDS = int8(1)
const SCIP_BRANCHDIR_FIXED = int8(2)
const SCIP_BRANCHDIR_AUTO = int8(3)

typealias SCIP_LPSolStat Int8
const SCIP_LPSOLSTAT_NOTSOLVED = int8(0)
const SCIP_LPSOLSTAT_OPTIMAL = int8(1)
const SCIP_LPSOLSTAT_INFEASIBLE = int8(2)
const SCIP_LPSOLSTAT_UNBOUNDEDRAY = int8(3)
const SCIP_LPSOLSTAT_OBJLIMIT = int8(4)
const SCIP_LPSOLSTAT_ITERLIMIT = int8(5)
const SCIP_LPSOLSTAT_TIMELIMIT = int8(6)
const SCIP_LPSOLSTAT_ERROR = int8(7)

typealias SCIP_BoundType Int8
const SCIP_BOUNDTYPE_LOWER = int8(0)
const SCIP_BOUNDTYPE_UPPER = int8(1)

typealias SCIP_SideType Int8
const SCIP_SIDETYPE_LEFT = int8(0)
const SCIP_SIDETYPE_RIGHT = int8(1)

typealias SCIP_RowOriginType Int8
const SCIP_ROWORIGINTYPE_UNSPEC = int8(0)
const SCIP_ROWORIGINTYPE_CONS = int8(1)
const SCIP_ROWORIGINTYPE_SEPA = int8(2)

typealias SCIP_LPAlgo Int8
const SCIP_LPALGO_PRIMALSIMPLEX = int8(0)
const SCIP_LPALGO_DUALSIMPLEX = int8(1)
const SCIP_LPALGO_BARRIER = int8(2)
const SCIP_LPALGO_BARRIERCROSSOVER = int8(3)

typealias SCIP_VerbLevel Int8
const SCIP_VERBLEVEL_NONE = int8(0)
const SCIP_VERBLEVEL_DIALOG = int8(1)
const SCIP_VERBLEVEL_MINIMAL = int8(2)
const SCIP_VERBLEVEL_NORMAL = int8(3)
const SCIP_VERBLEVEL_HIGH = int8(4)
const SCIP_VERBLEVEL_FULL = int8(5)

typealias SCIP_NlpParam Int8
const SCIP_NLPPAR_FROMSCRATCH = int8(0)
const SCIP_NLPPAR_VERBLEVEL = int8(1)
const SCIP_NLPPAR_FEASTOL = int8(2)
const SCIP_NLPPAR_RELOBJTOL = int8(3)
const SCIP_NLPPAR_LOBJLIM = int8(4)
const SCIP_NLPPAR_INFINITY = int8(5)
const SCIP_NLPPAR_ITLIM = int8(6)
const SCIP_NLPPAR_TILIM = int8(7)
const SCIP_NLPPAR_OPTFILE = int8(8)
const SCIP_NLPPAR_FASTFAIL = int8(9)

typealias SCIP_NlpSolStat Int8
const SCIP_NLPSOLSTAT_GLOBOPT = int8(0)
const SCIP_NLPSOLSTAT_LOCOPT = int8(1)
const SCIP_NLPSOLSTAT_FEASIBLE = int8(2)
const SCIP_NLPSOLSTAT_LOCINFEASIBLE = int8(3)
const SCIP_NLPSOLSTAT_GLOBINFEASIBLE = int8(4)
const SCIP_NLPSOLSTAT_UNBOUNDED = int8(5)
const SCIP_NLPSOLSTAT_UNKNOWN = int8(6)

typealias SCIP_NlpTermStat Int8
const SCIP_NLPTERMSTAT_OKAY = int8(0)
const SCIP_NLPTERMSTAT_TILIM = int8(1)
const SCIP_NLPTERMSTAT_ITLIM = int8(2)
const SCIP_NLPTERMSTAT_LOBJLIM = int8(3)
const SCIP_NLPTERMSTAT_UOBJLIM = int8(4)
const SCIP_NLPTERMSTAT_NUMERR = int8(5)
const SCIP_NLPTERMSTAT_EVALERR = int8(6)
const SCIP_NLPTERMSTAT_MEMERR = int8(7)
const SCIP_NLPTERMSTAT_LICERR = int8(8)
const SCIP_NLPTERMSTAT_OTHER = int8(9)

typealias SCIP_ParamType Int8
const SCIP_PARAMTYPE_BOOL = int8(0)
const SCIP_PARAMTYPE_INT = int8(1)
const SCIP_PARAMTYPE_LONGINT = int8(2)
const SCIP_PARAMTYPE_REAL = int8(3)
const SCIP_PARAMTYPE_CHAR = int8(4)
const SCIP_PARAMTYPE_STRING = int8(5)

typealias SCIP_ParamSetting Int8
const SCIP_PARAMSETTING_DEFAULT = int8(0)
const SCIP_PARAMSETTING_AGGRESSIVE = int8(1)
const SCIP_PARAMSETTING_FAST = int8(2)
const SCIP_PARAMSETTING_OFF = int8(3)

typealias SCIP_ParamEmphasis Int8
const SCIP_PARAMEMPHASIS_DEFAULT = int8(0)
const SCIP_PARAMEMPHASIS_CPSOLVER = int8(1)
const SCIP_PARAMEMPHASIS_EASYCIP = int8(2)
const SCIP_PARAMEMPHASIS_FEASIBILITY = int8(3)
const SCIP_PARAMEMPHASIS_HARDLP = int8(4)
const SCIP_PARAMEMPHASIS_OPTIMALITY = int8(5)
const SCIP_PARAMEMPHASIS_COUNTER = int8(6)

typealias SCIP_Objsense Int8
const SCIP_OBJSENSE_MAXIMIZE = int8(-1)
const SCIP_OBJSENSE_MINIMIZE = int8(+1)

typealias SCIP_Result Int8
const SCIP_DIDNOTRUN = int8(1)
const SCIP_DELAYED = int8(2)
const SCIP_DIDNOTFIND = int8(3)
const SCIP_FEASIBLE = int8(4)
const SCIP_INFEASIBLE = int8(5)
const SCIP_UNBOUNDED = int8(6)
const SCIP_CUTOFF = int8(7)
const SCIP_SEPARATED = int8(8)
const SCIP_NEWROUND = int8(9)
const SCIP_REDUCEDDOM = int8(10)
const SCIP_CONSADDED = int8(11)
const SCIP_CONSCHANGED = int8(12)
const SCIP_BRANCHED = int8(13)
const SCIP_SOLVELP = int8(14)
const SCIP_FOUNDSOL = int8(15)
const SCIP_SUSPENDED = int8(16)
const SCIP_SUCCESS = int8(17)

typealias SCIP_Retcode Int8
const SCIP_OKAY = int8(+1)
const SCIP_ERROR = int8(0)
const SCIP_NOMEMORY = int8(-1)
const SCIP_READERROR = int8(-2)
const SCIP_WRITEERROR = int8(-3)
const SCIP_NOFILE = int8(-4)
const SCIP_FILECREATEERROR = int8(-5)
const SCIP_LPERROR = int8(-6)
const SCIP_NOPROBLEM = int8(-7)
const SCIP_INVALIDCALL = int8(-8)
const SCIP_INVALIDDATA = int8(-9)
const SCIP_INVALIDRESULT = int8(-10)
const SCIP_PLUGINNOTFOUND = int8(-11)
const SCIP_PARAMETERUNKNOWN = int8(-12)
const SCIP_PARAMETERWRONGTYPE = int8(-13)
const SCIP_PARAMETERWRONGVAL = int8(-14)
const SCIP_KEYALREADYEXISTING = int8(-15)
const SCIP_MAXDEPTHLEVEL = int8(-16)
const SCIP_BRANCHERROR = int8(-17)

typealias SCIP_Efficiacychoice Int8
const SCIP_EFFICIACYCHOICE_LP = int8(0)
const SCIP_EFFICIACYCHOICE_RELAX = int8(1)
const SCIP_EFFICIACYCHOICE_NLP = int8(2)

typealias SCIP_Stage Int8
const SCIP_STAGE_INIT = int8(0)
const SCIP_STAGE_PROBLEM = int8(1)
const SCIP_STAGE_TRANSFORMING = int8(2)
const SCIP_STAGE_TRANSFORMED = int8(3)
const SCIP_STAGE_INITPRESOLVE = int8(4)
const SCIP_STAGE_PRESOLVING = int8(5)
const SCIP_STAGE_EXITPRESOLVE = int8(6)
const SCIP_STAGE_PRESOLVED = int8(7)
const SCIP_STAGE_INITSOLVE = int8(8)
const SCIP_STAGE_SOLVING = int8(9)
const SCIP_STAGE_SOLVED = int8(10)
const SCIP_STAGE_EXITSOLVE = int8(11)
const SCIP_STAGE_FREETRANS = int8(12)
const SCIP_STAGE_FREE = int8(13)

typealias SCIP_Setting Int8
const SCIP_UNDEFINED = int8(0)
const SCIP_DISABLED = int8(1)
const SCIP_AUTO = int8(2)
const SCIP_ENABLED = int8(3)

typealias SCIP_SolOrigin Int8
const SCIP_SOLORIGIN_ORIGINAL = int8(0)
const SCIP_SOLORIGIN_ZERO = int8(1)
const SCIP_SOLORIGIN_LPSOL = int8(2)
const SCIP_SOLORIGIN_NLPSOL = int8(3)
const SCIP_SOLORIGIN_RELAXSOL = int8(4)
const SCIP_SOLORIGIN_PSEUDOSOL = int8(5)
const SCIP_SOLORIGIN_UNKNOWN = int8(6)

typealias SCIP_Status Int8
const SCIP_STATUS_UNKNOWN = int8(0)
const SCIP_STATUS_USERINTERRUPT = int8(1)
const SCIP_STATUS_NODELIMIT = int8(2)
const SCIP_STATUS_TOTALNODELIMIT = int8(3)
const SCIP_STATUS_STALLNODELIMIT = int8(4)
const SCIP_STATUS_TIMELIMIT = int8(5)
const SCIP_STATUS_MEMLIMIT = int8(6)
const SCIP_STATUS_GAPLIMIT = int8(7)
const SCIP_STATUS_SOLLIMIT = int8(8)
const SCIP_STATUS_BESTSOLLIMIT = int8(9)
const SCIP_STATUS_OPTIMAL = int8(10)
const SCIP_STATUS_INFEASIBLE = int8(11)
const SCIP_STATUS_UNBOUNDED = int8(12)
const SCIP_STATUS_INFORUNBD = int8(13)

typealias SCIP_NodeType Int8
const SCIP_NODETYPE_FOCUSNODE = int8(0)
const SCIP_NODETYPE_PROBINGNODE = int8(1)
const SCIP_NODETYPE_SIBLING = int8(2)
const SCIP_NODETYPE_CHILD = int8(3)
const SCIP_NODETYPE_LEAF = int8(4)
const SCIP_NODETYPE_DEADEND = int8(5)
const SCIP_NODETYPE_JUNCTION = int8(6)
const SCIP_NODETYPE_PSEUDOFORK = int8(7)
const SCIP_NODETYPE_FORK = int8(8)
const SCIP_NODETYPE_SUBROOT = int8(9)
const SCIP_NODETYPE_REFOCUSNODE = int8(10)

typealias SCIP_Varstatus Int8
const SCIP_VARSTATUS_ORIGINAL = int8(0)
const SCIP_VARSTATUS_LOOSE = int8(1)
const SCIP_VARSTATUS_COLUMN = int8(2)
const SCIP_VARSTATUS_FIXED = int8(3)
const SCIP_VARSTATUS_AGGREGATED = int8(4)
const SCIP_VARSTATUS_MULTAGGR = int8(5)
const SCIP_VARSTATUS_NEGATED = int8(6)

typealias SCIP_Vartype Int8
const SCIP_VARTYPE_BINARY = int8(0)
const SCIP_VARTYPE_INTEGER = int8(1)
const SCIP_VARTYPE_IMPLINT = int8(2)
const SCIP_VARTYPE_CONTINUOUS = int8(3)

typealias SCIP_DomchgType Int8
const SCIP_DOMCHGTYPE_DYNAMIC = int8(0)
const SCIP_DOMCHGTYPE_BOTH = int8(1)
const SCIP_DOMCHGTYPE_BOUND = int8(2)

typealias SCIP_BoundchgType Int8
const SCIP_BOUNDCHGTYPE_BRANCHING = int8(0)
const SCIP_BOUNDCHGTYPE_CONSINFER = int8(1)
const SCIP_BOUNDCHGTYPE_PROPINFER = int8(2)

typealias SCIP_VBCColor Int8
const SCIP_VBCCOLOR_UNSOLVED = int8(3)
const SCIP_VBCCOLOR_SOLVED = int8(2)
const SCIP_VBCCOLOR_CUTOFF = int8(4)
const SCIP_VBCCOLOR_CONFLICT = int8(15)
const SCIP_VBCCOLOR_MARKREPROP = int8(11)
const SCIP_VBCCOLOR_REPROP = int8(12)
const SCIP_VBCCOLOR_SOLUTION = int8(14)
const SCIP_VBCCOLOR_NONE = int8(-1)


text:
@testset "avoid overflow during add/subtract" begin
    x = 2*10^11
    @test FastRational(1, 2x) + FastRational(1, 3x) == FastRational(5, 6x)
    @test FastRational(1, 2x) - FastRational(1, 3x) == FastRational(1, 6x)
end

@testset "constructors1" begin
    @test FastRational(2,4).num == 1
    @test FastRational(2,4).den == 2
    @test FastRational((2,4)).num == 2
    @test FastRational((2,4)).den == 4
end

@testset "convert from Rational" begin
    @test FastRational(1//2) == FastRational(1, 2)
    @test FastRational(2, 4) == FastRational(1, 2)
end

@testset "conversion type" begin
    @test typeof(convert(FastRational{Int64,FastRationals.MayReduce}, FastRational(1,2))) ==
    FastRational{Int64,FastRationals.MayReduce}
end

@testset "inverse matrix" begin
    B = FastRational.([(1,2) (1,3); (1,3) (1,4)])
    @test inv(B) == FastRational.([(18,1) (-24,1); (-24,1) (36,1)])
end


text:
using SimilaritySearch, LinearAlgebra

function create_database()
    X = rand(Float32, 4, 100_000)
    for c in eachcol(X) normalize!(c) end
    MatrixDatabase(X)
end

function main()
    db = create_database()
    H = SearchGraph(; db, dist=NormalizedCosineDistance(), verbose=false)
    index!(H; parallel_block=128)
    @time hI, hD = allknn(H, 32; parallel=true)
end

main()

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
#Kacper Kulczak 279079
#Program do zadania 2.12 na pracownie z przedmiotu Analiza Numeryczna M

using Plotly;
include("interps.jl")

#interpolacja dwoma wielomianami Newtona rybki
ryba_x1 =[3.,13.,26.,38.,49.,65.,77.,93.,106.,123.,131.,138.,144.,149.,155.]
ryba_y1 = [81.,72.,65.,60.,57.,55.,56.,58.,63.,71.,75.,81.,86.,92.,99.]
ryba_x2 =[3.,13.,26.,38.,49.,65.,77.,93.,106.,123.,131.,138.,144.,149.,155.]
ryba_y2 =[81.,90.,97.,102.,105.,107.,106.,104.,99.,91.,87.,81.,76.,70.,63.]
Newton1 = interpNewton(ryba_x1,ryba_y1)
Newton2 = interpNewton(ryba_x2,ryba_y2)
xss = linspace(ryba_x1[1], ryba_x1[end], 1000)
trace1 = scatter(;x = xss, y = map(Newton1,xss) , mode = "lines" ,name = "wielomian 1")
trace2 = scatter(;x = xss, y = map(Newton2,xss) , mode = "lines" ,name = "wielomian 2")
trace3 = scatter(;x=ryba_x1, y=ryba_y1, mode = "markers", name = "wezły interpolacyjne 1")
trace4 = scatter(;x=ryba_x2, y=ryba_y2, mode = "markers", name = "węzy interpolacyjne 2")
layout = Layout(;title="Interpolacja wielomianami interpolacyjnymi Newtona",
    xaxis=attr(title="", ),
    yaxis=attr(title="",))
plot([trace1,trace2,trace3,trace4],layout)

#interpolacja dwoma funckjami sklejanymi rybki
Spline1 = interpspline(ryba_x1,ryba_y1)
Spline2 = interpspline(ryba_x2,ryba_y2)
xss = linspace(ryba_x1[1], ryba_x1[end], 1000)
trace1 = scatter(;x = xss, y = map(Spline1,xss) , mode = "lines", name="funkcja 1")
trace2 = scatter(;x = xss, y = map(Spline2,xss) , mode = "lines", name="funckja 2")
trace3 = scatter(;x=ryba_x1, y=ryba_y1, mode = "markers", name="węzły interpolacyjne")
trace4 = scatter(;x=ryba_x2, y=ryba_y2, mode = "markers", name ="węzły interpolacyjne")
layout = Layout(;title="Interpolacja naturalnymi funkcjami sklejanymi",
    xaxis=attr(title="", ),
    yaxis=attr(title="",))
plot([trace1,trace2,trace3,trace4],layout)

#rybka za pomocą interpolowanej krzywej parametrycznej
nryba_y1 = myreverse(ryba_y1)
pop!(nryba_y1)
nryba_x1 = myreverse(ryba_x1)
pop!(nryba_x1)
rybap_x = append!(nryba_x1 , ryba_x2 )
rybap_y = append!(nryba_y1, ryba_y2)
a = 1
b = length(rybap_x)
t = linspace(a, b, b)
Splinex = interpspline(t,rybap_x)
Spliney = interpspline(t,rybap_y)
xss = linspace(a, b, 1000)
trace1 = scatter(;x = map(Splinex,xss), y = map(Spliney,xss), mode = "lines",name= "funkcja")
trace3 = scatter(;x=rybap_x, y=rybap_y, mode="markers",name = "węzły interpolacji")
layout = Layout(;title = "Interpolacja krzywej parametrycznej naturalną funckją sklejaną")
plot([trace1,trace3], layout)
plot([trace1,trace3], layout)

#mapa Polski z punktow ze list na pracowni
tajem_xs, tajem_ys = ([3.7,3.2,2.7,2.1,1.7,1.1,0.7,0.4,0.4,0.5,0.3,0.6,0.6,0.7,0.6,0.8,0.8,0.6,0.8,0.9,1.1,1.4,1.8,1.7,1.9,2.2,2.1,2.7,2.6,3.3,3.5,3.7,3.9,4.2,4.3,4.5,4.7,5.0,5.5,5.9,6.2,6.4,6.3,6.5,6.8,7.2,7.1,7.2,6.8,6.7,6.8,6.4,6.2,6.9,6.8,6.6,6.5,6.1,5.5,5.0,4.6,4.1,3.7,3.4,3.2,3.7],[6.4,6.7,6.5,6.4,6.0,5.9,5.7,5.7,5.4,5.0,4.6,4.3,4.0,3.7,3.2,2.9,2.6,2.4,2.3,2.4,2.2,2.1,2.0,1.8,1.4,1.5,1.8,1.6,1.4,1.3,0.9,0.6,0.8,0.7,0.4,0.5,0.7,0.6,0.8,0.6,0.4,0.3,0.7,1.2,1.7,2.0,2.2,2.4,2.8,3.2,3.6,3.9,4.2,4.5,5.1,5.6,6.0,6.2,6.1,6.2,6.2,6.3,6.0,6.1,6.5,6.4])
tempxs, tempys  =interpparametric(tajem_xs, tajem_ys)
tajem_trace0 = scatter(; x=tempxs ,  y=tempys , mode = "lines", name = "Interpolowana linia")
trace1= scatter(;x = tajem_xs, y = tajem_ys, mode = "markers", name = "Dane punkty")
plot([tajem_trace0,trace1])

#spirala
a=0
b=10*pi
aa,bb,cc =make_param_curve(t->t*cos(t),t->t*sin(t),linspace(a,b,40) ,linspace(a,b,4000), "")

#Cykloida wydluzona
a=-3*pi
b=3*pi
aa,bb,cc =make_param_curve(t->1/3*t - sin(t),t->1/3 -cos(t),linspace(a,b,40) ,linspace(a,b,4000), "Cykloida wydłużona")

#Krzywa z petlami
a=-2*pi
b=2*pi
aa,bb,cc =make_param_curve(t->t+sin(2*t),t->t+sin(3*t),linspace(a,b,40) ,linspace(a,b,4000), "Krzywa z pętlami")

#Globus
a=0
b=9*pi
aa,bb,cc =make_param_curve(t->t*cos(t),t->sin(t),linspace(a,b,40) ,linspace(a,b,4000), "Globus")

#Inna krzywa
a=0
<<<<<<< HEAD
b=15*pi
aa,bb,cc =make_param_curve(t->t*cos(t),t->t,linspace(a,b,40) ,linspace(a,b,4000), "")
=======
b=10*3.14
make_param_curve(t->t*cos(t),t->t*sin(t),linspace(a,b,50) ,linspace(a,b,4000), "")

>>>>>>> 31480e236ad47585b5bccec604fa1446572e7ba5


text:
import Base.Sort, Base.Sort.Algorithm
import Base:sort!
struct NOOPSortAlg <: Base.Sort.Algorithm end
const NOOPSort = NOOPSortAlg()

sort!(x, ::Integer, ::Integer, ::ShortestPaths.NOOPSortAlg, ::Base.Sort.Ordering) = x

"""
    struct BFS <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [Breadth-First Search algorithm](https://en.m.wikipedia.org/wiki/Breadth-first_search).

An optional sorting algorithm may be specified (default = no sorting).
Sorting helps maintain cache locality and will improve performance on
very large graphs; for normal use, sorting will incur a performance
penalty.

`BFS` is the default algorithm used when a source is specified
but no distance matrix is specified.

### Implementation Notes
`BFS` supports the following shortest-path functionality:
- (optional) multiple sources
- all destinations
"""
struct BFS{T<:Base.Sort.Algorithm} <: ShortestPathAlgorithm
    sort_alg::T
end

BFS() = BFS(NOOPSort)

struct BFSResult{U<:Integer} <: ShortestPathResult
    parents::Vector{U}
    dists::Vector{U}
end

function shortest_paths(
    g::AbstractGraph{U},
    ss::AbstractVector{U},
    alg::BFS,
    ) where U<:Integer

    n = nv(g)
    dists = fill(typemax(U), n)
    parents = zeros(U, n)
    visited = falses(n)
    n_level = one(U)
    cur_level = Vector{U}()
    sizehint!(cur_level, n)
    next_level = Vector{U}()
    sizehint!(next_level, n)
    @inbounds for s in ss
        dists[s] = zero(U)
        visited[s] = true
        push!(cur_level, s)
    end
    while !isempty(cur_level)
        @inbounds for v in cur_level
            @inbounds @simd for i in outneighbors(g, v)
                if !visited[i]
                    push!(next_level, i)
                    dists[i] = n_level
                    parents[i] = v
                    visited[i] = true
                end
            end
        end
        n_level += one(U)
        empty!(cur_level)
        cur_level, next_level = next_level, cur_level
        sort!(cur_level, alg=alg.sort_alg)
    end
    return BFSResult(parents, dists)
end

shortest_paths(g::AbstractGraph{U}, ss::AbstractVector{<:Integer}, alg::BFS) where {U<:Integer} = shortest_paths(g, U.(ss), alg)
shortest_paths(g::AbstractGraph{U}, s::Integer, alg::BFS) where {U<:Integer} = shortest_paths(g, Vector{U}([s]), alg)


text:
using Test, YaoBase, YaoBlocks, BitBasis, YaoArrayRegister, Random
using StatsBase: mean

@testset "measure ghz" begin
    # GHZ state
    st = normalize!(ArrayReg(bit"0000") + ArrayReg(bit"1111"))
    Random.seed!(1234)

    # measure it at 1, 2
    # should collapse to 0000 or 1111 since entangled
    g = Measure(4; locs=(1, 2))
    st |> g

    @test g.results[1] == 0 ? st.state[end] == 0 : st.state[1] == 0
end

@testset "collapseto" begin
    Random.seed!(1234)

    st = rand_state(5; nbatch=3)
    g = Measure(5; locs=(1, 2), collapseto=bit"11")
    st |> g
    for k in 1:32
        if !(st.state[k] ≈ 0.0)
            @test all(bit(k-1; len=5)[1:2] .== 1)
        end
    end
end

@testset "error handling" begin
    @test_throws ErrorException Measure(5; locs=(1, 2), collapseto=bit"11", remove=true)
    @test_throws ErrorException mat(Measure(5; locs=(1, 2), collapseto=bit"11"))
end

@testset "op-measures" begin
    Random.seed!(1234)

    reg = rand_state(6, nbatch=10)
    op = repeat(3, X)

    # measure!
    reg2 = reg |> copy
    res = measure!(op, reg2, 2:4)
    res2 = measure!(op, reg2, 2:4)
    @test size(res) == (10,)
    @test res2 == res

    # measure_collapseto!
    reg2 = reg |> copy
    res = measure_collapseto!(op, reg2, 2:4)
    reg2 |> repeat(6, H, 2:4)
    res2 = measure_collapseto!(op, reg2, 2:4)
    @test size(res) == (10,) == size(res2)
    @test all(res2 .== 1)

    # measure_remove!
    reg2 = reg |> copy
    res = measure_remove!(op, reg2, 2:4)
    reg2 |> repeat(3, H, 2:3)
    @test size(res) == (10,)
    @test nqubits(reg2) == 3

    reg = repeat(ArrayReg(ComplexF64[1,-1]/sqrt(2.0)), 10)
    @test measure!(X, reg) |> mean ≈ -1
    reg = repeat(ArrayReg(ComplexF64[1.0,0]), 1000)
    @test abs(measure!(X, reg) |> mean) < 0.1
end


text:
using RCall, JLD

grid_k = [collect(1:25:10000); 10000]

Vnk_numerical_accuracy = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Vnk_10000.jld", "Vnk_numerical_accuracy") |> x -> Float64.(x)
Pkn_numerical_accuracy = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Pkn_10000.jld", "Pkn_numerical_accuracy") |> x -> Float64.(x)
Cnk_numerical_accuracy = load("test/graphical_tests/saves_for_graphical_tests/accuracy_Cnk_10000.jld", "Cnk_numerical_accuracy") |> x -> Float64.(x)

R"library(tidyverse)"
R"p = tibble(Pkn = $Pkn_numerical_accuracy, Cnk = $Cnk_numerical_accuracy, Vnk = $Vnk_numerical_accuracy, k = $(collect(grid_k)) ) %>%
    gather(Var, Accuracy, -k) %>%
    mutate(Var = factor(Var, levels = c('Vnk', 'Cnk', 'Pkn'))) %>%
    ggplot(aes(x = k, y = Accuracy, colour = Var, group = Var)) +
    theme_bw() +
    geom_point() +
    geom_line() +
    geom_hline(yintercept = 64, linetype = 'dotted', colour = 'red') +
    scale_colour_discrete(name='')"

R"pdf('test/graphical_tests/figures_graphical_tests/accuracy_PknCnkVnk_10000.pdf', height = 5, width = 5)
    plot(p)
    dev.off()"

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
# tensortoolbox.jl
#
# Main file for module TensorToolbox, a Julia package for working with
# with tensors and tensor operations

# module TensorToolbox
#
# # Exports
# #---------
# # Types:
# export VectorSpace, Field, ElementarySpace, InnerProductSpace, EuclideanSpace
# export ComplexSpace, CartesianSpace, GeneralSpace, GradedSpace, ZNSpace
# export Parity, ZNCharge, U1Charge
# export CompositeSpace, ProductSpace, InvariantSpace
# export IndexSpace, TensorSpace, AbstractTensor, Tensor, InvariantTensor
# export TruncationScheme
# export AbstractTensorMap
# export SpaceError, IndexError
#
# # general vector space methods
# export space, dim, dual, fieldtype, sectortype #, cnumber, iscnumber, directsum, fuse
# export ⊗, ⊕, ℂ, ℝ # some unicode
#
# # vector spaces with symmetries
# export sectors, invariant
#
# # tensor characteristics
# export spacetype, tensortype, numind, order
# # tensor constructors
# export tensor, tensorcat
# # index manipulations
# export insertind, deleteind, fuseind, splitind
# # tensor operations
# export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct
# export tensorcopy!, tensoradd!, tensortrace!, tensorcontract!, tensorproduct!
# # tensor factorizations
# export leftorth, rightorth
# export leftorth!, rightorth!, svd!
#
# # truncation schemes
# export notrunc, truncerr, truncdim, truncspace
#
# # tensor maps
# export domain, codomain, hermitian, posdef, tensormap
#
# # tensor networks
# export AbstractTensorNetwork, TensorNetwork
# export network, optimizecontract

using Base: tuple_type_head, tuple_type_tail, tuple_type_cons, tail, front, setindex,
            Iterators.product, ImmutableDict
include("auxiliary/auxiliary.jl")

# Exception types:
#------------------
abstract type TensorException <: Exception end

# Exception type for all errors related to sector mismatch
struct SectorMismatch{S<:Union{Void,String}} <: TensorException
    message::S
end
SectorMismatch()=SectorMismatch{Void}(nothing)
Base.show(io::IO, ::SectorMismatch{Void}) = print(io, "SectorMismatch()")

# Exception type for all errors related to vector space mismatch
struct SpaceMismatch{S<:Union{Void,String}} <: TensorException
    message::S
end
SpaceMismatch()=SpaceMismatch{Void}(nothing)
Base.show(io::IO, ::SpaceMismatch{Void}) = print(io, "SpaceMismatch()")

# Exception type for all errors related to invalid tensor index specification.
struct IndexError{S<:Union{Void,String}} <: TensorException
    message::S
end
IndexError()=IndexError{Void}(nothing)
Base.show(io::IO, ::IndexError{Void}) = print(io, "IndexError()")

# Tensor product operator
#-------------------------
⊗(a, b, c, d...)=⊗(a, ⊗(b, c, d...))

# Definitions and methods for superselection sectors (quantum numbers)
#----------------------------------------------------------------------
include("sectors/sectors.jl")

# Definitions and methods for vector spaces
#-------------------------------------------
include("spaces/vectorspaces.jl")

# Constructing and manipulating fusion trees and iterators thereof
#------------------------------------------------------------------
include("fusiontrees/fusiontrees.jl")

# # Definitions and methods for tensors
# #-------------------------------------
# import TensorOperations
# intentionally shadow original TensorOperation methods for StridedArray objects

# define truncation schemes for tensors
# include("tensors/truncation.jl")

# general definitions
#include("tensors/abstracttensor.jl")

# specific implementation
# include("tensors/tensor.jl")

#end


text:
function compute_recall(neighbors, ground_truth)
    """
    neighbors:(n, topk)
    ground_truth:(n, gt)
    """
    total = 0
    n = size(neighbors)[1]
    for i in 1:n
        gt_row = ground_truth[i,:]
        row = neighbors[i,:]
        total += length(intersect(gt_row, row))
    end

    return total / length(ground_truth)

end

function  recall_atN_(neighbors_matrix, ground_truth)
    """
    Args:
        neighbors_matrix:(nq,topK=512)
        ground_truth: 1d or 2d
    """
    topk = size(neighbors_matrix)[2]
    if length(size(ground_truth)) == 1
        ng = 1
    end

    if length(size(ground_truth)) == 2
        ng = size(ground_truth)[2]
    end

    N = [1,2,4,8,10,16,20,32,64,100,128,256,512]

    recall_list=[]
    N_list = []
    for tn in N
        if ng <= tn <= topk
            neighbors_matrix_topN = neighbors_matrix[:,1:tn]

            recall = compute_recall(neighbors_matrix_topN, ground_truth)
            append!(recall_list, recall)
            append!(N_list,tn)

            println("recall ",ng,"@",tn," = ",recall)
        end
    end
    return recall_list,N_list
end

function recall_atN(neighbors_matrix,ground_truth)
    """
    Args:
        neighbors_matrix:(nq,topK=512)
        ground_truth:(nq,) or (nq,>=10)
    """
    if length(size(ground_truth)) == 2
        ground_truth_1 = ground_truth[:,1]
        ground_truth_10 = ground_truth[:,1:10]
        r1, N1 = recall_atN_(neighbors_matrix, ground_truth_1)
        print("\n")
        r10, N10 = recall_atN_(neighbors_matrix, ground_truth_10)
        print("\n")
        println("N=",N1)
        println("recall1@N:",r1)
        println("N=",N10)
        println("recall10@N:",r10)
    end
    if length(size(ground_truth)) == 1
        r1, N1 = recall_atN_(neighbors_matrix, ground_truth)
        print("\n")
        println("N=",N1)
        println("recall1@N:",r1)
    end
end


text:
abstract type ABasic end

struct Node
    data::String
    class::Type{<:ABasic}
    args::NTuple{N,Node} where N  # array does not work since Node([]]) == Node([]]) returns false
end
Node(data::String, class::Type{<:ABasic}) = Node(data, class, ())
data(self::Node) = self.data
class(self::Node) = self.class
args(self::Node) = self.args

function Base.show(io::IO, self::Node)
    _show(io, self, class(self))
end
_show(io::IO, self::Node, ::Type{<:Any}) = print(io, "$(Meta.parse("($(args(self)[1])) $(data(self)) ($(args(self)[2]))"))")

abstract type ANode <: ABasic end
abstract type AAtom <: ANode end
abstract type AOperator <: ANode end

abstract type AUnknown <: AAtom end
abstract type AVariable <: AUnknown end
struct Variable{T} <: AVariable where T
end
_show(io::IO, self::Node, ::Type{<:AAtom}) = print(io, data(self))

# Only used for truth table purposes
abstract type AWild <: AUnknown end
struct Wild{T} <: AWild where T
end

abstract type ACall <: AOperator end
function (f::Node)(x::Node)
    return Node("()", ACall, (f, x))
end
function (f::Node)(Vararg...)
    return Node("()", f, Vararg)
end
_show(io::IO, self::Node, ::Type{<:ACall}) = print(io, "$(Meta.parse("$(args(self)[1])($(args(self)[2]))"))")

abstract type ASingleton <: AAtom end
abstract type ASingletonStatement <: ASingleton end
abstract type ASingletonSet <: ASingleton end

abstract type ASingletonExpression <: ASingleton end

wild_symbol(name::String, set::Type{<:ASingletonSet}) = Node(name, Wild{set})
wild_symbol(name::String, set::Node) = Node(name, Wild{class(set)})

symbol(name::String, set::Type{<:ASingletonSet}) = Node(name, Variable{set})
symbol(name::String, set::Node) = Node(name, Variable{class(set)})

# External types
abstract type ANonBasic <: ABasic end
convert(::Type{Node}, x::Number) = Node(String(x), ANonBasic)

abstract type ALogicOperator <: AOperator end
abstract type AAnd <: ALogicOperator end
struct And <: AAnd
end
∧(x::Node, y::Node) = Node("∧", AAnd, (x, y))

abstract type AForAll <: ALogicOperator end
ForAll(variable::Node, set::Node, statement::Node) = Node("A", AForAll, (variable, set, statement))

abstract type AOr <: ALogicOperator end
struct Or <: AOr
end
∨(x::Node, y::Node) = Node("∨", AOr, (x, y))

abstract type AExists <: ALogicOperator end
Exists(variable::Node, set::Node, statement::Node) = Node("E", AOr, (variable, set, statement))

abstract type ANot <: ALogicOperator end
struct Not <: ANot
end
~(x::Node) = Node("~", ANot, (x,))
_show(io::IO, self::Node, ::Type{<:ANot}) = print(io, "$(Meta.parse("$(data(self))($(args(self)[1]))"))")

abstract type AImplies <: ALogicOperator end
struct Implies <: AImplies
end
⟹(x::Node, y::Node) = Node("⟹", AImplies, (x, y))

abstract type AEquivalent <: ALogicOperator end
struct Equivalent <: AEquivalent
end
≡(x::Node, y::Node) = Node("≡", AEquivalent, (x, y))

abstract type ASubset <: ALogicOperator end
struct Subset <: ASubset
end
⊂(x::Node, y::Node) = Node("⊂", ASubset, (x, y))

abstract type AElementOf <: ALogicOperator end
struct ElementOf <: AElementOf
end
∈(x::Node, y::Node) = Node("∈", AElementOf, (x, y))

abstract type ASetOperator <: AOperator end
abstract type ASetUnion <: ASetOperator end
struct SetUnion <: ASetUnion
end
∪(x::Node, y::Node) = Node("∪", ASetUnion, (x, y))

abstract type ASetIntersection <: ASetOperator end
struct SetIntersection <: ASetIntersection
end
∩(x::Node, y::Node) = Node("∩", ASetIntersection, (x, y))

abstract type ASetDifference <: ASetOperator end
struct SetDifference <: ASetDifference
end

abstract type AConditionSet <: ASetOperator end
struct ConditionSet <: AConditionSet
end
ConditionSet(x::Node, condition::Node) = Node("|", AConditionSet, (x, condition))  # {x | condition}
_show(io::IO, self::Node, ::Type{<:AConditionSet}) = print(io, "$(Meta.parse("ConditionSet($(args(self)[1]), $(args(self)[2]))"))")

abstract type AExpressionOperator <: AOperator end