text:
export load_wt_janelia, parse_wt_janelia, add_wt_a

function load_wt_janelia(gui,vid_name,data_name)

    meas=matread(string(gui.folder_path,data_name))

    gui.vid_path = string(gui.folder_path,vid_name)

    gui.tif=load(gui.vid_path);

    parse_wt_janelia(gui,meas)

    nothing
end

function parse_wt_janelia(gui,xx)
   gui.whiskers=[Whisker() for i=1:size(gui.tif,3), j=1:2]

    for i=1:length(xx["whiskers"]["x"])

        wid = xx["measurements"]["label"][i]+1
        t = xx["whiskers"]["time"][i]+1

        if (wid<size(gui.whiskers,2)+1)&(wid>0)

            pos_x = xx["whiskers"]["x"][i]
            pos_y = 480 - xx["whiskers"]["y"][i]

            gui.whiskers[t,wid].pos = [Point2f0(pos_x[j],pos_y[j]) for j=1:length(pos_x)]

            gui.whiskers[t,wid].ang = xx["measurements"]["angle"][i]
        end
    end

    nothing
end

function add_wt_a(gui,wnum,cnum)

    gui.y_data=zeros(Float32,size(gui.tif,3),3)

    for i=1:size(gui.whiskers,1)
       gui.y_data[i,cnum] = gui.whiskers[i,wnum].ang
    end

    #Find first nonzero value, calculate mean of after values and set first values equal to it
    f_ind = findfirst(gui.y_data[:,cnum] .!= 0.0)

    if f_ind>1
        mymean = mean(gui.y_data[f_ind:end,cnum])
        gui.y_data[1:(f_ind-1),cnum] = mymean
    end

    #Interpolate between missing slider_values
    data_inds=find(gui.y_data[:,cnum].!=0.0);
    mydata=gui.y_data[data_inds,cnum];

    itp=interpolate((data_inds,),mydata,Gridded(Linear()))

    for i=1:size(gui.y_data,1)
        gui.y_data[i,cnum]=itp[i]
        gui.y_data[i,cnum] = gui.y_data[i,cnum] * -1
    end

    gui.y_data[:,cnum] = gui.y_data[:,cnum] - mean(gui.y_data[:,cnum])

    nothing
end

function add_wt_callbacks(gui)

end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots
using StatsPlots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 5E1
"""
(1) Normal regression
(2) <- Multiple regression without intercept
(3) Normal regression of the difference
(4) <- "Regular" multiple regression
"""

# %% 5E2
"""

"""

# %% 5E3
"""

"""

# %% 5E4
"""
(1)
μ(A = 1) = α + β_A
μ(B = 1) = α + β_B
μ(C = 1) = α
μ(D = 1) = α + β_D

(2)
μ(A = 1) = α + β_A
μ(B = 1) = α + β_B
μ(C = 1) = α + β_C
μ(D = 1) = α + β_D

(3)
μ(A = 1) = α
μ(B = 1) = α + β_B
μ(C = 1) = α + β_C
μ(D = 1) = α + β_D

(4)
μ(A = 1) = α_A
μ(B = 1) = α_B
μ(C = 1) = α_C
μ(D = 1) = α_D

(5)
μ(A = 1) = α_A * (1 - 0 - 0 - 0)
μ(B = 1) = α_B
μ(C = 1) = α_C
μ(D = 1) = α_D

(2) is the odd one out because it has 5 parameters instead of the 4 of all the other models. I'm not sure if there going to be differences in uncertainty between (1)/(3) and (4)/(5).
"""

ext: jl


text:
using SymPy
import PyCall
using Test

## Test of importing an external module using `pyimport_conda` and `import_from`.
## from issue #262

stats = PyCall.pyimport_conda("sympy.stats", "sympy")

# this brings in all functions into the session using introspection, specialized on a
# symbolic first argument
import_from(stats)

@testset "Test stats module" begin

    @vars mu z x
    @vars sigma positive=true
    X = stats.Normal("x", mu, sigma)

    # some operations:
    @test variance(X) == sigma^2
    @test E(X) == mu
    P(Gt(X,3)) ## or P(X ≫ 3)

    Z = stats.Normal("Z", 0, 1)  # needs `stats.` as "Z" is not symbolic,
    # could also do this:
    @vars Z
    Z = Normal(Z, 0, 1)
    @test E(Z) == 0
    @test variance(Z) == 1

    simplify(P(Gt(Z, 1)))

    @test 0.68 <= 1 - 2*N(P(Gt(Z, 1))) <= 0.69

    @vars Y
    @vars a b c
    Y = DiscreteUniform(Y, (a,b,c))
    density(Y).dict
    @test density(Y).dict[a] == 1//3

    # And density
    density(X)  #NormalDistribution(mu, sigma)
    density(X).pdf(z)  # need the `pdf` call here

end


text:
# This file is a part of SimilaritySearch.jl

export allknn

"""
    allknn(g::AbstractSearchContext, k::Integer; parallel=false, pools=getpools(g)) -> knns, dists

Computes all the k nearest neighbors (all vs all) using the index `g`. It removes self references.

Parameters:

- `g`: the index
- `k`: the number of neighbors to retrieve
- `parallel`: If true, the construction will use all threads available threads
- `pools`: A pools object, dependent of `g`

Returns:

- `knns` a (k, n) matrix of identifiers; the i-th column corresponds to the i-th object in the dataset.
    Zeros can happen to the end of each column meaning that the retrieval was less than the desired `k`
- `dists` a (k, n) matrix of distances; the i-th column corresponds to the i-th object in the dataset.
    Zero values in `knns` should be ignored in `dists`

"""
function allknn(g::AbstractSearchContext, k::Integer; parallel=false, pools=getpools(g))
    n = length(g)
    knns = zeros(Int32, k, n)
    dists = Matrix{Float32}(undef, k, n)
    allknn(g, knns, dists; parallel, pools)
end

function _allknn_loop(g::SearchGraph, i, knns, dists, pools)
    k = size(knns, 1) + 1
    res = getknnresult(k, pools)
    vstate = getvstate(length(g), pools)
    # the loop helps to overcome when the current nn is in a small clique (smaller the the desired k)
    for h in g.links[i] # hints
        visited(vstate, h) && continue
        @inbounds search(g.search_algo, g, g[i], res, h, pools; vstate)
        length(res) == k && break
    end

    # again for the same issue
    if length(res) < k
        for j in 1:2k
            h = rand(1:length(g))
            visited(vstate, h) && continue
            @inbounds search(g.search_algo, g, g[i], res, h, pools; vstate)
            length(res) == k && break
        end
    end

    _allknn_inner_loop(res, i, knns, dists)
end

function _allknn_loop(g, i, knns, dists, pools)
    k = size(knns, 1) + 1
    res = getknnresult(k, pools)
    @inbounds search(g, g[i], res)
    _allknn_inner_loop(res, i, knns, dists)
end

@inline function _allknn_inner_loop(res, i,  knns, dists)
    j = 0
    @inbounds for (id, dist) in res
        i == id && continue
        j += 1
        knns[j, i] = id
        dists[j, i] = dist
    end
end

"""
allknn(g, knns, dists; parallel=false, pools=getpools(g)) -> knns, dists

Computes all the k nearest neighbors (all vs all) using the index `g`. It removes self references.

Arguments:

- `g`: the index
- `knns`: an uninitialized integer matrix of (k, n) size for storing the `k` nearest neighbors of the `n` elements
- `dists`: an uninitialized floating point matrix of (k, n) size for storing the `k` nearest distances of the `n` elements
- `parallel`: If true, the construction will use all threads available threads
- `pools`: A pools object, dependent of `g`

Results:

- `knns` and `dists` are returned. Note that the index can retrieve less than `k` objects, and these are represented as zeros at the end of each column (can happen)
"""
function allknn(g::AbstractSearchContext, knns::AbstractMatrix{Int32}, dists::AbstractMatrix{Float32}; parallel=false, pools=getpools(g))
    n = length(g)
    @assert n > 0

    if parallel
        Threads.@threads for i in 1:n
            _allknn_loop(g, i, knns, dists, pools)
        end
    else
        for i in 1:n
            _allknn_loop(g, i, knns, dists, pools)
        end
    end

    knns, dists
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
# SPD-License-Identifier: MIT

function load_filelist(uttids, feadir, ext)
    filelist= []
    open(uttids, "r") do f
        for line in eachline(f)
            push!(filelist, joinpath(feadir, line*ext))
        end
    end
    filelist
end

function makemodel(T, modelconf, μ, Σ)
    C = modelconf["ncomponents"]
    pstrength = modelconf["pstrength"]
    σ = modelconf["noise_init"]
    D = length(μ)
    μ₀ = Array{T}(μ)
    W₀ = Array{T}(Hermitian(inv(Σ)))
    diagΣ₀ = Array{T}(diag(Σ))

    @debug "creating model with $C components (covtype = $(modelconf["covtype"]))"
    components = []
    for c in 1:C
        if modelconf["covtype"] == "full"
            push!(components, Normal(T, D; μ₀, W₀, pstrength, σ))
        elseif modelconf["covtype"] == "diagonal"
            push!(components, NormalDiag(T, D; μ₀, diagΣ₀, pstrength, σ))
        else
            error("unknown cov type $(modelconf["covtype"])")
        end
    end
    Mixture(T, components = Tuple(components))
end

function savemodel(dir, model, modelconf, inferenceconf, step = -1)
    mtype = typeof(model)
    data = Dict(
        "modelconf" => modelconf,
        "inferenceconf" => inferenceconf,
        "mtype" => mtype,
        "mstate" => todict(model)
    )
    modelname = step > 0 ? "$step.jld2" : "final.jld2"

    save(joinpath(dir, modelname), data)
end

function loadmodel(path)
    dict = load(path)
    mtype = dict["mtype"]
    mstate = dict["mstate"]
    model = fromdict(mtype, mstate)
    modelconf = dict["modelconf"]
    model, modelconf
end

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
function extract(input)
    f = open(input, "r")
    #CHECK FILE TYPE
    occursin("ooTextFile", readline(f)) ? print("Successuflly opened ooTextFile") : print("File type NOT: \"ooTextFile")#line#1
    #CHECK OBJECT CLASS
    if !(readline(f) == "Object class = \"TextGrid\"")#line#2
     throw(UndefVarError(:NotTextGridFile))
    end
    readline(f)#SKIP BLANK LINE (line#3)
    #READ START TIME
    xmin = parse(Float64,readline(f)[8:end])#line#4
    #READ END TIME
    xmax = parse(Float64,readline(f)[8:end])#line#5
    #READ TIERS INFO
    tiers = occursin("tiers? <exists>", readline(f))#line#6
    tierSize = parse(Int64,readline(f)[8:end])#line#7
    #CREATE VECTOR VARIABLES
    class = Array{String}(undef, 1, tierSize)
    name = Array{String}(undef, 1, tierSize)
    xmin = Array{Float64}(undef, 1, tierSize)
    xmax = Array{Float64}(undef, 1, tierSize)
    intervalSize = Array{Int64}(undef, 1, tierSize)
    interval = Vector{Vector}(undef, tierSize)
    readline(f)#SKIP "item[]" (line#8)
    for i in 1:tierSize
        readline(f)#SKIP "item[i]" (line#9....)
        class[i] = readline(f)[18:end-2]#
        #println("class: " * readline(f)[18:end-2])
        #name[i] = readline(f)[17:end-2]#
        println("name: " * readline(f)[17:end-2])
        #xmin[i] = parse(Float64,readline(f)[15:end])#
        println("xmin: " * readline(f)[15:end])
        #xmax[i] = parse(Float64,readline(f)[15:end])#
        println("xmax: " * readline(f)[15:end])
        #intervalSize[i] = parse(Int64,readline(f)[26:end])#
        if occursin("IntervalTier", class[i])
            intervalSize[i] = parse(Int64,readline(f)[26:end])#println("intervalSize: " * readline(f)[26:end])
        else
            intervalSize[i] = parse(Int64,readline(f)[23:end])
        end
        #println("intervalSize: " * intervalSize[i])
        intervalLines = Vector{Any}(undef,intervalSize[i])
        println("tier#",i)
        for j in 1:intervalSize[i]
            #println(readline(f),readline(f),readline(f),readline(f),"\n")
            readline(f)#SKIP "intervals [j]:" (line#15)
            #println([readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
            if occursin("IntervalTier", class[i])
                intervalLines[j] = [parse(Float64,readline(f)[19:end]),
                                    parse(Float64,readline(f)[19:end]),
                                    readline(f)[21:end-2]]
            else
                intervalLines[j] = [parse(Float64,readline(f)[21:end]),
                                    readline(f)[21:end-2]]
            end
            #append!(interval,[readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
            if(j == intervalSize[i])
                interval[i] = copy(intervalLines)
            end
        end
    end
    interval end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using SLiDE

UNITS = "billions of us dollars (USD)"

function SLiDE.sum_over(df::DataFrame, col::Array{Symbol,1}; values_only = true, keepkeys = false)

    inp_keys = df[:,find_oftype(df, Not(AbstractFloat))]
    val_cols = find_oftype(df, AbstractFloat)
    by_cols = setdiff(propertynames(df), [col; val_cols])

    gd = groupby(df, by_cols);
    df = combine(gd, val_cols .=> sum)

    # df = by(df, by_cols, Pair.(val_cols, sum))
    df = edit_with(df, Rename.(setdiff(propertynames(df), by_cols), val_cols))

    keepkeys && (df = leftjoin(inp_keys, df, on = by_cols))
    return values_only ? df[:,val_cols[1]] : df
end

# ******************************************************************************************
#   READ BLUENOTE DATA -- For benchmarking!
# ******************************************************************************************
BLUE_DIR = joinpath("data", "windc_output", "2a_build_national_cgeparm_raw")
bluenote_lst = [x for x in readdir(joinpath(SLIDE_DIR, BLUE_DIR)) if occursin(".csv", x)]
bluenote = Dict(Symbol(k[1:end-4]) => sort(edit_with(
    read_file(joinpath(BLUE_DIR, k)), Rename(:Val, :value))) for k in bluenote_lst)

# Add supply/use info for checking.
BLUE_DIR_IN = joinpath("data", "windc_output", "1b_stream_windc_base")
[bluenote[k] = sort(edit_with(read_file(joinpath(BLUE_DIR_IN, string(k, "_units.csv"))), [
    Rename.([:Dim1,:Dim2,:Dim3,:Dim4,:Val], [:yr,:i,:j,:units,:value]);
    Replace.([:i,:j], "upper", "lower")])) for k in [:supply, :use]]

[bluenote[k][!,:value] .= round.(bluenote[k][:,:value]*1E-3, digits=3) # convert millions -> billions USD
    for k in [:supply,:use]]

# ******************************************************************************************
#   READ SETS AND SLiDE SUPPLY/USE DATA.
# ******************************************************************************************
y = read_file(joinpath("dev","buildstream","setlist.yml"));
set = Dict(Symbol(k) => sort(read_file(joinpath(y["SetPath"]..., ensurearray(v)...)))[:,1]
    for (k,v) in y["SetInput"])

# Read supply/use data.
DATA_DIR = joinpath("data", "input")
io_lst = convert_type.(Symbol, ["supply", "use"])
io = Dict(k => read_file(joinpath(DATA_DIR, string(k, ".csv"))) for k in io_lst)

# Perform some minor edits that will apply to all parameters in order to maintain
# consistency with the data in build_national_cgeparm_raw.gdx for easier benchmarking.

for k in keys(io)
    global io[k] = edit_with(io[k], Drop.([:value, :units], [0., "all"], "=="))
    global io[k] = io[k] |> @filter(_.yr in set[:yr]) |> DataFrame

    # global io[k][!,:value] .= round.(io[k][:,:value]*1E-3, digits=3)
    # global io[k][!,:units] .= UNITS
end

io[:supply], io[:use] = fill_zero(io[:supply], io[:use]; permute_keys = true)

# ******************************************************************************************
#   PARTITION DATA INTO PARAMETERS.
# ******************************************************************************************

# Read from use data.
#!  id0(yr,i(ir_use),j(jc_use)) = use(yr,ir_use,jc_use);
#!  fd0(yr,i(ir_use),fd(jc_use)) = use(yr,ir_use,jc_use);
#!  va0(yr,va(ir_use),j(jc_use)) = use(yr,ir_use,jc_use);
#!  ts0(yr,ts(ir_use),j(jc_use)) = use(yr,ir_use,jc_use);
#!  x0(yr,i(ir_use)) = use(yr,ir_use,"exports");
io[:id0] = io[:use] |> @filter(.&(_.i in set[:i],  _.j in set[:j] ))  |> DataFrame
io[:fd0] = io[:use] |> @filter(.&(_.i in set[:i],  _.j in set[:fd]))  |> DataFrame
io[:va0] = io[:use] |> @filter(.&(_.i in set[:va], _.j in set[:j] ))  |> DataFrame
io[:ts0] = io[:use] |> @filter(.&(_.i in set[:ts], _.j in set[:j] ))  |> DataFrame
io[:x0]  = io[:use] |> @filter(.&(_.i in set[:i],  _.j == "exports")) |> DataFrame

io[:va0] = edit_with(io[:va0], Rename(:i, :va));
io[:fd0] = edit_with(io[:fd0], Rename(:j, :fd));

# Read from supply data.
#!  ys0(yr,j(jc_supply),i(ir_supply)) = supply(yr,ir_supply,jc_supply);
#!  m0(yr,i(ir_supply)) = supply(yr,ir_supply,"imports");
#!  mrg0(yr,i(ir_supply)) = supply(yr,ir_supply,"margins");
#!  trn0(yr,i(ir_supply)) = supply(yr,ir_supply,"trncost");
#!  cif0(yr,i(ir_supply)) = supply(yr,ir_supply,"ciffob");
#!  duty0(yr,i(ir_supply)) = supply(yr,ir_supply,"duties");
#!  tax0(yr,i(ir_supply)) = supply(yr,ir_supply,"tax");
#!  sbd0(yr,i(ir_supply)) = - supply(yr,ir_supply,"subsidies");
io[:ys0]   = io[:supply] |> @filter(.&(_.i in set[:i], _.j in set[:j]))   |> DataFrame
io[:m0]    = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "imports")) |> DataFrame
io[:mrg0]  = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "margins")) |> DataFrame
io[:trn0]  = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "trncost")) |> DataFrame
io[:cif0]  = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "ciffob"))  |> DataFrame
io[:duty0] = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "duties"))  |> DataFrame
io[:tax0]  = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "tax"))     |> DataFrame
io[:sbd0]  = io[:supply] |> @filter(.&(_.i in set[:i], _.j == "subsidies")) |> DataFrame

# Treat negative inputs as outputs.
#!  ys0(yr,j,i) = ys0(yr,j,i) - min(0,id0(yr,i,j));
#!  id0(yr,i,j) = max(0,id0(yr,i,j));
#!  ts0(yr,'subsidies',j) = - ts0(yr,'subsidies',j);
#!  sbd0(yr,i(ir_supply)) = - supply(yr,ir_supply,"subsidies");
io[:ys0][!,:value] = io[:ys0][:,:value] - min.(0, io[:id0][:,:value])
io[:id0][!,:value] = max.(0, io[:id0][:,:value])
io[:ts0][io[:ts0][:,:i] .== "subsidies", :value] *= -1
io[:sbd0][!,:value] *= -1

# Adjust transport margins for transport sectors according to CIF/FOB adjustments.
# Insurance imports are specified as net of adjustments.
#!  trn0(yr,i)$(cif0(yr,i) AND NOT SAMEAS(i,'ins')) = trn0(yr,i) + cif0(yr,i);
#!  m0(yr,i)$(SAMEAS(i,'ins')) = m0(yr,i) + cif0(yr,i);
#!  cif0(yr,i) = 0;
i_ins = io[:cif0][:,:i] .== "ins"
io[:trn0][.!i_ins, :value] .= io[:trn0][.!i_ins,:value] + io[:cif0][.!i_ins,:value]
io[:m0  ][  i_ins, :value] .= io[:m0][i_ins,:value] + io[:cif0][i_ins,:value]
io[:cif0][      !, :value] .= 0.0

# Aggregate supply and gross output
#!  s0(yr,j) = sum(i,ys0(yr,i,j));
#!  y0(yr,i) = sum(j, ys0(yr,j,i));
io[:s0] = sum_over(io[:ys0], :i; values_only = false)   # aggregate supply
io[:y0] = sum_over(io[:ys0], :j; values_only = false)   # gross output

# ******************************************************************************************

# Initialize empty DataFrames where values will be calculated.
io[:a0]  = fill_zero((yr = set[:yr], i = set[:i]))
io[:tm0] = fill_zero((yr = set[:yr], i = set[:i]))
io[:ta0] = fill_zero((yr = set[:yr], i = set[:i]))
io[:ms0] = fill_zero((yr = set[:yr], i = set[:i], m = set[:m]))
io[:md0] = fill_zero((yr = set[:yr], m = set[:m], i = set[:i]))

# Balance of payments deficit
#!  bopdef(yr) = 0;
io[:bopdef] = fill_zero((yr = set[:yr], ))

# Margin supply
#!  ms0(yr,i,"trd") = max(-mrg0(yr,i),0);
#!  ms0(yr,i,'trn') = max(-trn0(yr,i),0);
io[:ms0][io[:ms0][:,:m] .== "trd", :value] .= max.(-io[:mrg0][:,:value], 0)
io[:ms0][io[:ms0][:,:m] .== "trn", :value] .= max.(-io[:trn0][:,:value], 0)

# Margin demand
#!  md0(yr,"trd",i) = max(mrg0(yr,i),0);
#!  md0(yr,'trn',i) = max(trn0(yr,i),0);
io[:md0][io[:md0][:,:m] .== "trd", :value] .= max.(io[:mrg0][:,:value], 0)
io[:md0][io[:md0][:,:m] .== "trn", :value] .= max.(io[:trn0][:,:value], 0)

# Household supply
# Move household supply of recycled goods into the domestic output market
# from which some may be exported. Net out margin supply from output.
#!  fs0(yr,i) = -min(0, fd0(yr,i,'pce'));
#!  y0(yr,i) = sum(j,ys0(yr,j,i)) + fs0(yr,i) - sum(m,ms0(yr,i,m));
io[:fs0] = io[:fd0] |> @filter(_.fd == "pce") |> DataFrame
io[:fs0][!,:value] .= - min.(io[:fs0][:,:value], 0)
io[:y0][!,:value]  .= sum_over(io[:ys0], :j) + io[:fs0][:,:value] - sum_over(io[:ms0], :m)

# Armington supply
#!  a0(yr,i) = sum(fd, fd0(yr,i,fd)) + sum(j, id0(yr,i,j));
io[:a0][!,:value] .= sum_over(io[:fd0], :fd) + sum_over(io[:id0], :j)

# Remove commodity taxes and subsidies on the goods which are produced solely
# for supplying retail sales margin:
#!  y0(yr,imrg) = 0;
#!  a0(yr,imrg) = 0;
#!  tax0(yr,imrg) = 0;
#!  sbd0(yr,imrg) = 0;
#!  x0(yr,imrg) = 0;
#!  m0(yr,imrg) = 0;
#!  md0(yr,m,imrg) = 0;
#!  duty0(yr,imrg) = 0;
# Here's how to do this: https://discourse.julialang.org/t/dataframes-obtaining-the-subset-of-rows-by-a-set-of-values/15923/10
[io[k][findall(in(set[:imrg]), io[k][:,:i]), :value] .= 0.0
    for k in [:y0, :a0, :tax0, :sbd0, :x0, :m0, :md0, :duty0]]

# Tax net subsidy rate on intermediate demand.
#!  tm0(yr,i)$duty0(yr,i) = duty0(yr,i)/m0(yr,i);
i_div = io[:m0][:,:value] .!= 0.0
io[:tm0][i_div, :value] .=  io[:duty0][i_div,:value] ./ io[:m0][i_div,:value]

# Import tariff
#!  ta0(yr,i)$(tax0(yr,i)-sbd0(yr,i)) = (tax0(yr,i) - sbd0(yr,i))/a0(yr,i);
i_div = io[:a0][:,:value] .!= 0.0
io[:ta0][i_div, :value] .= (io[:tax0][i_div,:value] - io[:sbd0][i_div,:value]) ./ io[:a0][i_div,:value]

# include("partitionbea_check.jl")
# println("Done.")

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using JuMP, Gurobi, Plots, LaTeXStrings, Distributions, LinearAlgebra
using BPPrivInfo

lb = 0.0 + eps(0.0)
ub = 1.0 - eps(1.0)
dist = MixtureModel([Uniform(lb, 1/2 + 1/100),
                     truncated(Normal(ub, 0.1), 1/2 - 1/100, ub)], [35/100, 65/100])

typedist = discretise(dist, 1000)
prob = probs(typedist)
belief = support(typedist)

function solve_discretised(dist::DiscreteNonParametric)
    # Create problem data
    prob = probs(dist)
    belief = support(dist)
    n = length(belief)

    # Initalise model
    BPPI = Model(Gurobi.Optimizer)

    # Create variables
    @variable(BPPI, 0 ≤ πG[1:n] ≤ 1)
    @variable(BPPI, 0 ≤ πB[1:n] ≤ 1)

    # Create objective
    @expression(BPPI, exppayoff, dot(prob, belief .* πG) + dot(prob, (1 .- belief) .* πB))
    @objective(BPPI, Max, exppayoff)

    # Create constraints
    @constraint(BPPI, obG, belief .* πG .≥ (1 .- belief) .* πB)
    @constraint(BPPI, obB, (1 .- belief) .* (1 .- πB) .≥ belief .* (1 .- πG))
    @constraint(BPPI,
                report[i = 1:n, j = 1:n; i ≠ j],
                belief[i] * πG[i] + (1 - belief[i]) * (1 - πB[i]) ≥ belief[i] * πG[j] + (1 - belief[i]) * (1 - πB[j]))

    # Solve model
    optimize!(BPPI)

    termstatus = termination_status(BPPI)
    termstatus == MOI.OPTIMAL || @warn("Termination status = $termstatus")

    # Generate plots
    plot(belief, value.(πG), label = L"\pi_G", legend = :outertopright, show = true)
    plot!(belief, value.(πB), label = L"\pi_B")
    # plot(belief, -dual.(UpperBoundRef.(πG)), label = L"\pi_G \le 1", legend = :outertopright, show = true)

    return BPPI, πG, πB, report
end

function plot_exp_payoff(dist::DiscreteNonParametric)
    x = support(dist)
    plot!(x, t -> expected_payoff(t, dist), label = L"V", legend = :outertopright)
end

model, πG, πB, report = solve_discretised(typedist)

get_report_duals(report, i, n) = dual.(report[i,j] for j in 1:n if j ≠ i)

function plot_report_duals(dist, report)
    belief = support(dist)
    n = length(belief)
    plt = plot(belief[2:end], get_report_duals(report, 1, n), legend = false)
    for j in 2:n
        plot!(plt, belief[1:end .≠ j], get_report_duals(report, j, n))
    end
    display(plt)
end

max_stars_repo_name: carlocab/BPPrivInfo

max_issues_repo_name: carlocab/BPPrivInfo

max_forks_repo_name: carlocab/BPPrivInfo


text:
function fold(r)
# [rw] = fold(r)

# Fold left wing of vector in "FFT buffer format"
# onto right wing
# J.O. Smith, 1982-2002

   sz = size(r);
   n = maximum(sz)
   if prod(sz) != n
     error("fold.m: input must be a vector");
   end

   if n < 3; return r;
   elseif mod(n,2)==1
       nt = Int((n+1)/2);
       rw = [ r[1], r[2:nt] + conj(r[n:-1:nt+1]), 0*ones(n-nt) ];
   else
       nt = Int(n/2);
       rf = cat(1, r[2:nt] , 0);
       rf = rf + conj(r[n:-1:nt+1]);
       rw = cat(1, r[1] , rf , 0*ones(n-nt-1) );
   end

   rw = reshape(rw,sz)
   return rw
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayFirewallDisabledRuleGroup <: SwaggerModel
    ruleGroupName::Any # spec type: Union{ Nothing, String } # spec name: ruleGroupName
    rules::Any # spec type: Union{ Nothing, Vector{Int32} } # spec name: rules

    function ApplicationGatewayFirewallDisabledRuleGroup(;ruleGroupName=nothing, rules=nothing)
        o = new()
        validate_property(ApplicationGatewayFirewallDisabledRuleGroup, Symbol("ruleGroupName"), ruleGroupName)
        setfield!(o, Symbol("ruleGroupName"), ruleGroupName)
        validate_property(ApplicationGatewayFirewallDisabledRuleGroup, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type ApplicationGatewayFirewallDisabledRuleGroup

const _property_map_ApplicationGatewayFirewallDisabledRuleGroup = Dict{Symbol,Symbol}(Symbol("ruleGroupName")=>Symbol("ruleGroupName"), Symbol("rules")=>Symbol("rules"))
const _property_types_ApplicationGatewayFirewallDisabledRuleGroup = Dict{Symbol,String}(Symbol("ruleGroupName")=>"String", Symbol("rules")=>"Vector{Int32}")
Base.propertynames(::Type{ ApplicationGatewayFirewallDisabledRuleGroup }) = collect(keys(_property_map_ApplicationGatewayFirewallDisabledRuleGroup))
Swagger.property_type(::Type{ ApplicationGatewayFirewallDisabledRuleGroup }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayFirewallDisabledRuleGroup[name]))}
Swagger.field_name(::Type{ ApplicationGatewayFirewallDisabledRuleGroup }, property_name::Symbol) =  _property_map_ApplicationGatewayFirewallDisabledRuleGroup[property_name]

function check_required(o::ApplicationGatewayFirewallDisabledRuleGroup)
    (getproperty(o, Symbol("ruleGroupName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ApplicationGatewayFirewallDisabledRuleGroup }, name::Symbol, val)
end


text:
#1
natural1(x, y) = strip(x) < strip(y)

#2
natural2(x, y) = replace(x, r"\s+" => " ") < replace(y, r"\s+" => " ")

#3
natural3(x, y) = lowercase(x) < lowercase(y)

#4
splitbynum(x) = split(x, r"(?<=\D)(?=\d)|(?<=\d)(?=\D)")
numstringtonum(arr) = [(n = tryparse(Float32, e)) != nothing ? n : e for e in arr]
function natural4(x, y)
    xarr = numstringtonum(splitbynum(x))
    yarr = numstringtonum(splitbynum(y))
    for i in 1:min(length(xarr), length(yarr))
        if typeof(xarr[i]) != typeof(yarr[i])
            a = string(xarr[i]); b = string(yarr[i])
        else
             a = xarr[i]; b = yarr[i]
        end
        if a == b
            continue
        else
            return a < b
        end
    end
    return length(xarr) < length(yarr)
end

#5
deart(x) = replace(x, r"^The\s+|^An\s+|^A\s+"i => "")
natural5(x, y) = deart(x) < deart(y)

#6
const accentdict = Dict(
'À'=> 'A', 'Á'=> 'A', 'Â'=> 'A', 'Ã'=> 'A', 'Ä'=> 'A',
'Å'=> 'A', 'Ç'=> 'C', 'È'=> 'E', 'É'=> 'E',
'Ê'=> 'E', 'Ë'=> 'E', 'Ì'=> 'I', 'Í'=> 'I', 'Î'=> 'I',
'Ï'=> 'I', 'Ñ'=> 'N', 'Ò'=> 'O', 'Ó'=> 'O',
'Ô'=> 'O', 'Õ'=> 'O', 'Ö'=> 'O', 'Ù'=> 'U', 'Ú'=> 'U',
'Û'=> 'U', 'Ü'=> 'U', 'Ý'=> 'Y', 'à'=> 'a', 'á'=> 'a',
'â'=> 'a', 'ã'=> 'a', 'ä'=> 'a', 'å'=> 'a', 'è'=> 'e',
'é'=> 'e', 'ê'=> 'e', 'ë'=> 'e', 'ì'=> 'i', 'í'=> 'i',
'î'=> 'i', 'ï'=> 'i', 'ð'=> 'd', 'ñ'=> 'n', 'ò'=> 'o',
'ó'=> 'o', 'ô'=> 'o', 'õ'=> 'o', 'ö'=> 'o', 'ù'=> 'u',
'ú'=> 'u', 'û'=> 'u', 'ü'=> 'u', 'ý'=> 'y', 'ÿ'=> 'y')
function tr(str, dict=accentdict)
    for (i, ch) in enumerate(str)
        if haskey(dict, ch)
            arr = split(str, "")
            arr[i] = string(dict[ch])
            str = join(arr)
        end
    end
    str
end

natural6(x, y) = tr(x) < tr(y)

#7
const ligaturedict = Dict('œ' => "oe", 'Œ' => "OE", 'æ' => "ae", 'Æ' => "AE", 'Ĳ' => "IJ")
natural7(x, y) = tr(x, ligaturedict) < tr(y, ligaturedict)

#8
const altsdict = Dict('ß' => "ss", 'ſ' => 's', 'ʒ' => 's')
natural8(x, y) = tr(x, altsdict) < tr(y, altsdict)

preprocessors = [natural1, natural2, natural2, natural3, natural4, natural5, natural6, natural7, natural8]

const testarrays = Vector{Vector{String}}([
["ignore leading spaces: 2-2", " ignore leading spaces: 2-1", "  ignore leading spaces: 2+0", "   ignore leading spaces: 2+1"],
["ignore m.a.s spaces: 2-2", "ignore m.a.s  spaces: 2-1", "ignore m.a.s   spaces: 2+0", "ignore m.a.s    spaces: 2+1"],
["Equiv. spaces: 3-3", "Equiv.\rspaces: 3-2", "Equiv.\x0cspaces: 3-1", "Equiv.\x0bspaces: 3+0", "Equiv.\nspaces: 3+1", "Equiv.\tspaces: 3+2"],
["cASE INDEPENENT: 3-2", "caSE INDEPENENT: 3-1", "casE INDEPENENT: 3+0", "case INDEPENENT: 3+1"],
["foo100bar99baz0.txt", "foo100bar10baz0.txt", "foo1000bar99baz10.txt", "foo1000bar99baz9.txt"],
["The Wind in the Willows", "The 40th step more", "The 39 steps", "Wanda"],
["Equiv. ý accents: 2-2", "Equiv. Ý accents: 2-1", "Equiv. y accents: 2+0", "Equiv. Y accents: 2+1"],
["Ĳ ligatured ij", "no ligature"],
["Start with an ʒ: 2-2", "Start with an ſ: 2-1", "Start with an ß: 2+0", "Start with an s: 2+1"]])

for (i, ltfunction) in enumerate(preprocessors)
    println("Testing sorting mod number $i. Sorted is: $(sort(testarrays[i], lt=ltfunction)).")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract NoDecision <: Exception
Base.show(io::IO, nd::NoDecision) = print(io, """
    POMCP failed to choose an action because the following exception was thrown:
    $nd

    To specify an action for this case, use the default_action solver parameter.
    """)

immutable AllSamplesTerminal <: NoDecision
    belief
end
Base.show(io::IO, ast::AllSamplesTerminal) = print(io, """
    POMCP failed to choose an action because all states sampled from the belief were terminal.

    To see the belief, catch this exception as ex and see ex.belief.

    To specify an action for this case, use the default_action solver parameter.
    """)

immutable ExceptionRethrow end

default_action(::ExceptionRethrow, belief, ex) = rethrow(ex)
default_action(f::Function, belief, ex) = f(belief, ex)
default_action(p::POMDPs.Policy, belief, ex) = action(p, belief)
default_action(a, belief, ex) = a


text:
module CDMrdata
    using CodecXz
    using CodecBzip2
    using RData
    export load_data, list_datasets
    include("dataloaders.jl")
end


text:
module Stardates

import Dates
import Printf
using TimeZones

include("constants.jl")
include("startofyear.jl")
include("sd.jl")

export Stardate, stardate, defaulttimezone, mediawiki, mstardate,
    nyse_closing_stardate

end # module


text:
function RenderObject(
        data::Dict{Symbol}, program, pre,
        bbs=Observable(Rect3f(Vec3f(0), Vec3f(1))),
        main=nothing
    )
    RenderObject(convert(Dict{Symbol,Any}, data), program, pre, bbs, main)
end

function Base.show(io::IO, obj::RenderObject)
    println(io, "RenderObject with ID: ", obj.id)
end

Base.getindex(obj::RenderObject, symbol::Symbol) = obj.uniforms[symbol]
Base.setindex!(obj::RenderObject, value, symbol::Symbol) = obj.uniforms[symbol] = value

Base.getindex(obj::RenderObject, symbol::Symbol, x::Function) = getindex(obj, Val(symbol), x)
Base.getindex(obj::RenderObject, ::Val{:prerender}, x::Function) = obj.prerenderfunctions[x]
Base.getindex(obj::RenderObject, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x]

Base.setindex!(obj::RenderObject, value, symbol::Symbol, x::Function)     = setindex!(obj, value, Val(symbol), x)
Base.setindex!(obj::RenderObject, value, ::Val{:prerender}, x::Function)  = obj.prerenderfunctions[x] = value
Base.setindex!(obj::RenderObject, value, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x] = value

"""
Represents standard sets of function applied before rendering
"""
struct StandardPrerender
    transparency::Observable{Bool}
    overdraw::Observable{Bool}
end

function (sp::StandardPrerender)()
    if sp.overdraw[]
        # Disable depth testing if overdrawing
        glDisable(GL_DEPTH_TEST)
    else
        glEnable(GL_DEPTH_TEST)
        glDepthFunc(GL_LEQUAL)
    end

    # Disable cullface for now, untill all rendering code is corrected!
    glDisable(GL_CULL_FACE)
    # glCullFace(GL_BACK)

    if sp.transparency[]
        # disable depth buffer writing
        glDepthMask(GL_FALSE)

        # Blending
        glEnable(GL_BLEND)
        glBlendEquation(GL_FUNC_ADD)

        # buffer 0 contains weight * color.rgba, should do sum
        # destination <- 1 * source + 1 * destination
        glBlendFunci(0, GL_ONE, GL_ONE)

        # buffer 1 is objectid, do nothing
        glDisablei(1, GL_BLEND)

        # buffer 2 is color.a, should do product
        # destination <- 0 * source + (1 - source) * destination
        glBlendFunci(2, GL_ZERO, GL_ONE_MINUS_SRC_COLOR)

    else
        glDepthMask(GL_TRUE)
        enabletransparency()
    end
end

struct StandardPostrender
    vao::GLVertexArray
    primitive::GLenum
end

function (sp::StandardPostrender)()
    render(sp.vao, sp.primitive)
end

struct StandardPostrenderInstanced{T}
    main::T
    vao::GLVertexArray
    primitive::GLenum
end

function (sp::StandardPostrenderInstanced)()
    renderinstanced(sp.vao, to_value(sp.main), sp.primitive)
end

struct EmptyPrerender end

(sp::EmptyPrerender)() = nothing

export EmptyPrerender
export prerendertype

function instanced_renderobject(data, program, bb=Observable(Rect3f(Vec3f(0), Vec3f(1))), primitive::GLenum=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(convert(Dict{Symbol,Any}, data), program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrenderInstanced(main, robj.vertexarray, primitive)
    robj
end

function std_renderobject(data, program, bb=Observable(Rect3f(Vec3f(0), Vec3f(1))), primitive=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(convert(Dict{Symbol,Any}, data), program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrender(robj.vertexarray, primitive)
    robj
end

prerendertype(::Type{RenderObject{Pre}}) where {Pre} = Pre
prerendertype(::RenderObject{Pre}) where {Pre} = Pre


text:
decode_id(c::IO) =
  let id = decode_int(c)
    if id == -1
      error("Backend Error: $(decode_String(c))")
    else
      id
    end
  end

encode_BIMLevel = encode_int
decode_BIMLevel = decode_int_or_error
encode_FloorFamily = encode_int
decode_FloorFamily = decode_int_or_error

function request_operation(conn::IO, name)
  write(conn, Int32(0))
  encode_String(conn, name)
  op = read(conn, Int32)
  if op == -1
    error(name * " is not available")
  else
    op
  end
end

interrupt_processing(conn::IO) = write(conn, Int32(-1))
#=
create_op(name::ASCIIString, argtypes::Array{DataType}, rettype::DataType) = (
  op = request_operation(name);
  (args...) -> (conn = current_connection();
                write(conn, Int32(op));
                @assert length(args) == length(argtypes) "Incorrect number of args";
                for actual_arg in [argtype(arg) for (arg, argtype) in zip(args, argtypes)]
                  write(conn, actual_arg)
                end;
                read(conn, rettype)))

macro def_op(name, argtypes, rettype)
  :($(esc(name)) = create_op($(string(name)), $argtypes, $rettype))
end

@def_op(Circle, [XYZ, XYZ, Float64], Int32)
@def_op(Sphere, [XYZ, Float64], Int32)

circle = (op = request_operation("Circle");
          conn = current_connection();
          (c, n, r)-> (send_data(conn, Int8(op));
                       send_data(conn, c);
                       send_data(conn, n);
                       send_data(conn, Float64(r);
                       read(conn, Int32))))
=#

function parse_c_signature(sig)
  m = match(r"^ *(public|) *(\w+) *([\[\]]*) +((?:\w|:|<|>)+) *\( *(.*) *\)", sig)
  ret = Symbol(m.captures[2])
  array_level = count(c -> c=='[', something(m.captures[3], ""))
  name = Symbol(m.captures[4])
  params = split(m.captures[5], r" *, *", keepempty=false)
  function parse_c_decl(decl)
    m = match(r"^ *((?:\w|:|<|>)+) *([\[\]]*) *(\w+)$", decl)
    (Symbol(m.captures[1]), count(c -> c=='[', something(m.captures[2], "")), Symbol(m.captures[3]))
  end
  (name, [parse_c_decl(decl) for decl in params], (ret, array_level))
end

#=
const julia_type_for_c_type = Dict(
  :byte => :Int8,
  :double => :Float64,
  :float => :Float64,
  :int => :Int,
  :bool => :Bool,
  :Point3d => :XYZ,
  :Point2d => :XYZ,
  :Vector3d => :VXYZ,
  :string => :ASCIIString,
  :ObjectId => :Int32,
  :Entity => :Int32,
  :BIMLevel => :Int32,
  :FloorFamily => :Int32,
    ) #Either Int32 or Int64, depending on the architecture

julia_type(ctype, is_array) = is_array ? :(Vector{$(julia_type(ctype, false))}) : julia_type_for_c_type[ctype]
=#
export show_rpc, step_rpc
const show_rpc = Parameter(false)
const step_rpc = Parameter(false)
function initiate_rpc_call(conn, opcode, name)
    if step_rpc()
        print(stderr, "About to call $(name) [press ENTER]")
        readline()
    end
    if show_rpc()
        print(stderr, name)
    end
end
function complete_rpc_call(conn, opcode, result)
    if show_rpc()
        println(stderr, result == nothing ? "-> nothing" : "-> $(result)")
    end
    result
end

#=
This does not play well with reset_backend because the opcode are in global
scope and completely separated from the backend and, thus, the operation
still tries to use the same opcode even after a reset. Must think about
storing the generated functions on the backend itself, so that we can drop
it and restart.

Another option would be to have the table of generated methods separated
from the channel on the backend side, but I'm afraid that will prevent multiple
clients (not sure they are useful, though).
=#

function my_symbol(prefix, name)
  Symbol(prefix * replace(string(name), ":" => "_"))
end

function translate(sym)
  Symbol(replace(string(sym), r"[:<>]" => "_", ))
end

function rpc(prefix, str)
    name, params, ret = parse_c_signature(str)
    func_name = my_symbol(prefix, name)
    #Expr(:(=), Expr(:call, esc(name), [p[1] for p in params]...), Expr(:call, :+, params[1][1], 2))
    esc(quote
        $func_name =
          let opcode = -1
              buf = IOBuffer()
              (conn, $([:($(p[3])) for p in params]...)) -> begin
                opcode = Int32(request_operation(conn, $(string(name))))
                global $func_name = (conn, $([:($(p[3])) for p in params]...)) -> begin
                  initiate_rpc_call(conn, opcode, $(string(name)))
                  take!(buf) # Reset the buffer just in case there was an encoding error on a previous call
                  write(buf, opcode)
                  $([:($(Symbol("encode_", translate(p[1]), "_array"^p[2]))(buf, $(p[3]))) for p in params]...)
                  write(conn, take!(buf))
                  complete_rpc_call(conn, opcode, $(Symbol("decode_", translate(ret[1]), "_array"^ret[2]))(conn))
                end
                $func_name(conn, $([:($(p[3])) for p in params]...))
          end
        end
        #export $func_name
      end)
end


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_csv.jl - https://github.com/strickek/ReadWriteDlm2.jl

"""

    readcsv2(source, T::Type=Any; opts...)

Equivalent to `readdlm2()` with delimiter `','` and `decimal='.'`. Default Type
`Any` activates parsing of `Bool`, `Int`, `Float64`, `Complex`, `Rational`,
`Missing`, `DateTime`, `Date` and `Time`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> B = Any[1 complex(1.5,2.7);1.0 1//3];

julia> writecsv2("test.csv", B)

julia> readcsv2("test.csv")
2×2 Array{Any,2}:
 1    1.5+2.7im
 1.0    1//3
```
"""
readcsv2(input; opts...) =
    readdlm2auto(input, ',', Nothing, '\n', false; decimal='.', opts...)

readcsv2(input, T::Type; opts...) =
    readdlm2auto(input, ',', T, '\n', false; decimal='.', opts...)

"""

    writecsv2(f, A; opts...)

Equivalent to `writedlm2()` with fixed delimiter `','` and `decimal='.'`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> B = Any[1 complex(1.5,2.7);1.0 1//3];

julia> writecsv2("test.csv", B)

julia> read("test.csv", String)
"1,1.5+2.7im\\n1.0,1//3\\n"
```
"""
writecsv2(f, a; opts...) = writedlm2auto(f, a, ','; decimal='.', opts...)


text:
export getSensMatVec

function getSensMatVec(v::Vector,m::Vector,pFor::EikonalInvParam)

    # extract pointers
    Q     = pFor.Sources
    P     = pFor.Receivers

    nsrc = size(Q,2)

    # allocate space for matvec product
    Jv   = zeros(size(P,2),nsrc)

	if length(v)!=length(m)
		warn("wrong size of vector v in EikonalInv::getSensMatVec.jl");
	end

    # derivative of mass matrix
	Tlin = zeros(size(v));

	if nsrc != length(pFor.eikonalParams)
		error("EikonalInv: Calculating derivatives without having fields... Run getData() to pFor");
	end

	pEiks = pFor.eikonalParams;
    for i=1:nsrc
		pEik_i = pEiks[i];
		if pFor.useFilesForFields
			file = matopen(getFieldsFileName());
			pEik_i.T1 		= read(file,string("T1_",i));
			pEik_i.ordering = read(file,string("ordering_",i));
			pEik_i.OP 		= read(file,string("OP_",i));
			close(file);
		end
		getSensMatVecEik(v,Tlin,pEik_i);
		if pFor.useFilesForFields
			FactoredEikonalFastMarching.clear!(pEik_i);
		end
		Jv[:,i] = P'*Tlin;
    end
    return vec(Jv)
end


text:
const function2latex = Dict{Symbol, String}(
    ## Greek alphabet
    :alpha => "\\alpha",
    :beta => "\\beta",
    :gamma => "\\gamma",
    :delta => "\\delta",
    :epsilon => "\\epsilon",
    :zeta => "\\zeta",
    :eta => "\\eta",
    :theta => "\\theta",
    :iota => "\\iota",
    :kappa => "\\kappa",
    :lambda => "\\lambda",
    :mu => "\\mu",
    :nu => "\\nu",
    :xi => "\\xi",
    :pi => "\\pi",
    :rho => "\\rho",
    :sigma => "\\sigma",
    :tau => "\\tau",
    :upsilon => "\\upsilon",
    :phi => "\\phi",
    :chi => "\\chi",
    :psi => "\\psi",
    :omega => "\\omega",
    :Gamma => "\\Gamma",
    :Delta => "\\Delta",
    :Theta => "\\Theta",
    :Lambda => "\\Lambda",
    :Xi => "\\Xi",
    :Pi => "\\Pi",
    :Sigma => "\\Sigma",
    :Upsilon => "\\Upsilon",
    :Phi => "\\Phi",
    :Psi => "\\Psi",
    :Omega => "\\Omega",
    ## Trinogometry
    :sin => "\\sin",
    :cos => "\\cos",
    :tan => "\\tan",
    :cot => "\\cot",
    :sec => "\\sec",
    :csc => "\\csc",
    :sinh => "\\sinh",
    :cosh => "\\cosh",
    :tanh => "\\tanh",
    :coth => "\\coth",
    :asin => "\\arcsin",
    :acos => "\\arccos",
    :atan => "\\arctan",
    :atan2 => "\\arctan",
    :asinh => "\\mathrm{arcsinh}",
    :sinc => "\\mathrm{sinc}",
    :acosh => "\\mathrm{arccosh}",
    :cosc => "\\mathrm{cosc}",
    :atanh => "\\mathrm{arctanh}",
    :acot => "\\mathrm{arccot}",
    :acoth => "\\mathrm{arccoth}",
    :asec => "\\mathrm{arcsec}",
    :sech => "\\mathrm{sech}",
    :asech => "\\mathrm{arcsech}",
    :acsc => "\\mathrm{arccsc}",
    :csch => "\\mathrm{csch}",
    :acsch => "\\mathrm{arccsch}",
    ## Misc
    :log => "\\log",
    :log10 => "\\log_{10}",
    :log2 => "\\log_{2}",
    :gamma => "\\Gamma", # The Gamma function
    )

const trigonometric_functions = [
    :sin,
    :cos,
    :tan,
    :cot,
    :sec,
    :csc,
    :sinh,
    :cosh,
    :tanh,
    :coth,
    :asin,
    :acos,
    :atan,
    :atan2,
    :asinh,
    :sinc,
    :acosh,
    :cosc,
    :atanh,
    :acot,
    :acoth,
    :asec,
    :sech,
    :asech,
    :acsc,
    :csch,
    :acsch,
    ]


text:
module ImageTransformations

using ImageCore
using CoordinateTransformations
using Rotations
using StaticArrays
using Interpolations, AxisAlgorithms
using OffsetArrays
using ColorVectorSpace

import Base: eltype, size, length
using Base: tail, Indices
using Base.Cartesian
using .ColorTypes: AbstractGray, TransparentGray, TransparentRGB

# these two symbols previously live in ImageTransformations
import ImageBase: restrict, restrict!

export

    restrict,
    imresize,
    center,
    warp,
    WarpedView,
    warpedview,
    InvWarpedView,
    invwarpedview,
    imrotate,
    CenterPoint

include("centerpoint.jl")
include("autorange.jl")
include("interpolations.jl")
include("warp.jl")
include("warpedview.jl")
include("invwarpedview.jl")
include("resizing.jl")
include("compat.jl")
include("deprecated.jl")

# TODO: move to warp.jl
@inline _getindex(A, v::StaticVector) = A[Tuple(v)...]
@inline _getindex(A::AbstractInterpolation, v::StaticVector) = A(Tuple(v)...)
@inline _getindex(A, v) = A[v...]
@inline _getindex(A::AbstractInterpolation, v) = A(v...)

center(img::AbstractArray{T,N}) where {T,N} = SVector{N}(map(_center, axes(img)))
_center(ind::AbstractUnitRange) = (first(ind)+last(ind))/2

end # module


text:
module Types

export EMImage, Segmentation, AffinityMap, ParamDict, SegMST, SegmentPairs, SegmentPairAffinities
export Timg, Tseg, Tsgm, Tec, Tecs

# type of raw image
const EMImage = Array{UInt8,3}

# type of segmentation
const Segmentation = Array{UInt32,3}

# type of affinity map
const AffinityMap = Array{Float32,4}

# type of parameter dictionary
const ParamDict = Dict{Symbol, Dict{Symbol, Any}}

const SegmentPairs = Array{UInt32,2}

const SegmentPairAffinities = Vector{Float32}

mutable struct SegMST
    segmentation::Segmentation
    segmentPairs::SegmentPairs
    segmentPairAffinities::SegmentPairAffinities
end

# defined for backward compatibility
const Timg  =  EMImage
const Tseg  =  Segmentation
const Taff  =  AffinityMap
const Tsgm  =  SegMST

end # end of module


text:
module DataUtils
#
# THIS IS NOT USED
# FIXME messing with Generics and Show() in particular
#
import DataFrames: DataFrame
import Statistics: mean, median, std, quantile
import Parameters: @with_kw

using  ScottishTaxBenefitModel
using  .Definitions

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function findallmax(A)
    max = A[firstindex(A)]
    max_num = []
    for i in (firstindex(A):lastindex(A))
        if A[i] > max
            max = A[i]
            max_num = []
        end
        if A[i] == max
            append!(max_num, i)
        end
    end
    return max_num
end

print(findallmax([3,4,9,2,7,0,1,9,2,3,9,4,5,6]))


text:
using Base: gc_alloc_count
using Test

using AssemblyView
using DZLinearAlgebra
using DZOptimization
using MultiFloats
using RungeKuttaToolKit

function rmk(args...)::Nothing
    print("\33[2K\r")
    print(args...)
    flush(stdout)
end

function say(args...)::Nothing
    print("\33[2K\r")
    println(args...)
    flush(stdout)
end

################################################################################

const KNOWN_SAFE_FUNCTIONS = Regex[]

add_safe_function!(func::String) = push!(KNOWN_SAFE_FUNCTIONS,
    Regex("\"?(julia_)?" * func * "(_[0-9]+)?\"?"))

add_safe_function!("o_memset")
add_safe_function!("gemv!")
add_safe_function!("generic_matvecmul!")

is_safe(func::String)::Bool = any(
    match(safe_func, func) != nothing for safe_func in KNOWN_SAFE_FUNCTIONS)

replace_underscores(expr::Expr, t::Symbol) = Expr(
    replace_underscores(expr.head, t),
    [replace_underscores(arg, t) for arg in expr.args]...)

replace_underscores(sym::Symbol, t::Symbol) = (sym == :_) ? t : sym

replace_underscores(n::Int, t::Symbol) = n

replace_underscores(qn::QuoteNode, t::Symbol) =
    QuoteNode(replace_underscores(qn.value, t))

function test_asm_calls(@nospecialize(func), type_template::Expr)
    for T in (:Float32, :Float64, :Float64x2, :Float64x4, :Float64x8)
        rmk("    ", rpad(func, 40, ' '), ": ", T)
        types = eval(replace_underscores(type_template, T))
        view_asm(IOBuffer(), func, types)
        for call in asm_offsets(func, types)
            if !is_safe(call)
                say("\nCOMPILATION TEST FAILURE:")
                say("Unsafe function call <", call,
                    "> in function <", func, ">.\n")
                exit()
            end
        end
    end
    say("    ", rpad(func, 40, ' '), ": PASS")
end

macro test_block(name::String, block::Expr)
    Expr(:block, :(say("Running ", $name, " tests...")),
        block, :(say("All ", $name, " tests passed.\n")))
end

say()

################################################################################

@test_block "DZMisc compilation" begin
    test_asm_calls(norm, :(Vector{_},))
    test_asm_calls(norm2, :(Vector{_},))
    test_asm_calls(normalize!, :(Vector{_},))
    test_asm_calls(dot, :(Vector{_}, Vector{_}))
    test_asm_calls(identity_matrix!, :(Matrix{_},))
end

@test_block "RKTK2 compilation" begin
    test_asm_calls(populate_explicit!, :(Matrix{_}, Vector{_}, Vector{_}, Int))
    test_asm_calls(populate_explicit!, :(Vector{_}, Matrix{_}, Vector{_}, Int))
    add_safe_function!("populate_explicit!")
    test_asm_calls(RungeKuttaToolKit.compute_butcher_weights!,
        :(Matrix{_}, Matrix{_}, Vector{Pair{Int,Int}}))
    test_asm_calls(RungeKuttaToolKit.backprop_butcher_weights!,
        :(Matrix{_}, Matrix{_}, Vector{_}, Matrix{_}, Vector{_},
            Vector{Int}, Vector{Tuple{Int,Int}}, Vector{Pair{Int,Int}}))
    add_safe_function!("compute_butcher_weights!")
    add_safe_function!("backprop_butcher_weights!")
    test_asm_calls(evaluate_residual2,
        :(Matrix{_}, Vector{_}, RKOCBackpropEvaluator{_}))
    test_asm_calls(evaluate_gradient!,
        :(Matrix{_}, Vector{_}, Matrix{_}, Vector{_}, RKOCBackpropEvaluator{_}))
    add_safe_function!("evaluate_residual2")
    add_safe_function!("evaluate_gradient!")
end

@test_block "DZOptimization compilation" begin
    test_asm_calls(DZOptimization.update_inverse_hessian!,
        :(Matrix{_}, _, Vector{_}, Vector{_}, Vector{_}))
    test_asm_calls(DZOptimization.quadratic_line_search,
        :(DZOptimization.StepObjectiveFunctor{
            RKOCExplicitBackpropObjectiveFunctor{_}, _, 1}, _, _))
    add_safe_function!("update_inverse_hessian!")
    add_safe_function!("quadratic_line_search")
    test_asm_calls(step!,
        :(BFGSOptimizer{RKOCExplicitBackpropObjectiveFunctor{_},
            RKOCExplicitBackpropGradientFunctor{_}, _},))
    add_safe_function!("step!")
end

################################################################################

function do_f64x8_ops(a::Vector{Float64x8}, b::Vector{Float64x8}, n::Int)
    @simd ivdep for i = 1 : n
        @inbounds a[i] = i
    end
    @simd ivdep for i = 1 : n
        @inbounds b[i] = sqrt(a[i])
    end
    @simd ivdep for i = 1 : n
        @inbounds a[i] = b[i] * b[i] + b[i] * b[i]
    end
    @inbounds [Float64(abs(a[i] - 2.0 * i)) for i = 1 : n]
end

function test_f64x8_ops(n::Int)
    a = Vector{Float64x8}(undef, n)
    b = Vector{Float64x8}(undef, n)
    do_f64x8_ops(a, b, 10)
    result, total_time, _, gc_time, gc_diff = @timed do_f64x8_ops(a, b, n)
    for i = 1 : length(result)
        @test result[i] < 1.0e-126 * i
    end
    @test total_time < 0.5
    @test gc_time == 0.0
    @test gc_alloc_count(gc_diff) < 100
end

@test_block "Float64x8 performance" begin
    test_f64x8_ops(100000)
end

################################################################################

@inline function vdpol!(f, y)
    @inbounds f[1] = y[2]
    @inbounds f[2] = 1000.0 * (1.0 - y[1] * y[1]) * y[2] - y[1]
end

function run_rksolver(::Type{T}, method, n::Int) where {T <: Real}
    y0 = T[2, 0]
    h = T(2000) / T(10000000)
    sol = RungeKuttaToolKit.RKSolver{T}(method(T), 2)
    for _ = 1 : n
        RungeKuttaToolKit.runge_kutta_step!(vdpol!, y0, h, sol)
    end
    y0
end

function time_rksolver(::Type{T}, method, n::Int) where {T <: Real}
    run_rksolver(T, method, 10)
    @time run_rksolver(T, method, n)
end

function test_rksolver(::Type{T}, method, n::Int) where {T <: Real}
    run_rksolver(T, method, 10)
    _, total_time, _, gc_time, gc_diff = @timed run_rksolver(T, method, n)
    @test total_time < 0.5
    @test gc_time == 0.0
    @test gc_alloc_count(gc_diff) < 100
end

using RungeKuttaToolKit.ExampleMethods
@test_block "RKSolver performance" begin
    test_rksolver(Float32, rk4_table,   100000)
    test_rksolver(Float64, rk4_table,   100000)
    test_rksolver(Float32, rkck5_table, 100000)
    test_rksolver(Float64, rkck5_table, 100000)
    test_rksolver(Float32, rkf8_table,  100000)
    test_rksolver(Float64, rkf8_table,  100000)
    test_rksolver(Float64x2, rk4_table, 10000)
    test_rksolver(Float64x3, rk4_table, 10000)
    test_rksolver(Float64x4, rk4_table, 10000)
    test_rksolver(Float64x5, rk4_table, 10000)
    test_rksolver(Float64x6, rk4_table, 10000)
    test_rksolver(Float64x7, rk4_table, 10000)
    test_rksolver(Float64x8, rk4_table, 10000)
end

################################################################################

function test_backprop_evaluator()
    x = Float64x4.(rand(BigFloat, 136))
    A = zeros(Float64x4, 16, 16)
    b = zeros(Float64x4, 16)
    k = 0
    for i = 2 : 16
        for j = 1 : i-1
            A[i,j] = x[k += 1]
        end
    end
    for i = 1 : 16
        b[i] = x[k += 1]
    end

    ev1 = RKOCEvaluator{Float64x4}(10, 16)
    ev2 = RKOCBackpropEvaluator{Float64x4}(10, 16)
    residual = Vector{Float64x4}(undef, 1205)
    jacobian = Matrix{Float64x4}(undef, 1205, 136)
    evaluate_residual!(residual, x, ev1)
    evaluate_jacobian!(jacobian, x, ev1)

    gA = zeros(Float64x4, 16, 16)
    gb = zeros(Float64x4, 16)
    evaluate_gradient!(gA, gb, A, b, ev2)
    result, total_time, _, gc_time, gc_diff =
        @timed evaluate_gradient!(gA, gb, A, b, ev2)

    @test total_time < 0.05
    @test gc_time == 0.0
    @test gc_alloc_count(gc_diff) == 0
    @test Float64(result) == Float64(sum(residual.^2))

    g1 = jacobian' * residual
    g1 += g1
    g2 = vcat([gA[i,j] for i = 2 : 16 for j = 1 : i-1], gb)
    for (x1, x2) in zip(g1, g2)
        @test Float64(x1) == Float64(x2)
    end
end

@test_block "Backprop evaluator" begin
    test_backprop_evaluator()
end


text:
make() =
  withenv("LINES"=>11, "COLUMNS"=>72) do
    runtests(joinpath("doc/src", ENV["DOCFILE"]), mod=Main)
    Base.active_repl.options.iocontext[:displaysize] =
      (convert(Integer, trunc(displaysize(Base.stdout)[1] * 2/3)), 72)
    nothing
  end

atreplinit() do repl
    try
        @eval using BenchmarkTools: @btime
        @eval using NarrativeTest
        @eval using Revise
        @eval using DataKnots
        @eval using DataKnots4Postgres
        @eval make()
        @async Revise.wait_steal_repl_backend()
    catch
        @warn "startup error?"
    end
end


text:
"""
    AbstractModelTrace

Common base class for various trace or "variable info" types.
"""
abstract type AbstractModelTrace end


text:
mutable struct VariableDomain
    lb::Float64;
    ub::Float64;
end

"""
    get_lb(domain)

Returns the lower bound in `domain`

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.get_lb(domain)
1.0
```
"""
get_lb(d::VariableDomain)::Float64 = d.lb;

"""
    set_lb!(domain, val)

Sets the lower bound in `domain` to `val`

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.set_lb!(domain, 3.0)
julia> OREnvironment.get_lb(domain)
3.0
```
"""
set_lb!(d::VariableDomain, lb::Float64) = d.lb = lb;

"""
    get_ub(domain)

Returns the upper bound in `domain`

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.get_ub(domain)
5.0
```
"""
get_ub(d::VariableDomain)::Float64 = d.ub;

"""
    set_ub!(domain, val)

Sets the upper bound in `domain` to `val`

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.set_ub!(domain, 9.0)
julia> OREnvironment.get_lb(domain)
9.0
```
"""
set_ub!(d::VariableDomain, ub::Float64) = d.ub = ub;

"""
    get_middle_point(domain)

Computes the middle point corresponding to the domain.

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.get_middle_point(domain)
3.0
```
"""
get_middle_point(d::VariableDomain)::Float64 = (d.ub + d.lb)/2

"""
    is_value_within_the_domain(domain, val)

Checks if `val` is between the lower and upper bound of `domain`.

# Example
```jldoctest
julia> domain = OREnvironment.VariableDomain(1.0,5.0);
julia> OREnvironment.is_value_within_the_domain(domain, 3.0)
true
julia> OREnvironment.is_value_within_the_domain(domain, 5)
true
julia> OREnvironment.is_value_within_the_domain(domain, 6)
false
```
"""
@inline function is_value_within_the_domain(d::VariableDomain, value::T)::Bool where {T<:Real}
    return d.lb <= value <= d.ub;
end

"""
    is_solution_within_bounds(arrayDomains, solution)

Checks if all the values of `solution` are within the domains. Each variable will have 1 domain.
"""
function is_solution_within_bounds(variablesDomain::Array{VariableDomain,1},
                                   solution::Solution)::Bool
    local nVariables::Int = length(variablesDomain);
    for i in 1:nVariables
        local value = get_solution(solution, i);
        local feasible::Bool = is_value_within_the_domain(variablesDomain[i], value);
        if feasible == false return false end
    end
    return true;
end

mutable struct DefaultProblem <: Problem
    costs::Array{Float64,1};
    constraints::Array{<:Constraint,1};
    variablesConstraints::Array{Array{Int,1},1};
    objSense::Symbol;
    variablesDomain::Array{VariableDomain,1};
end

"""
    constructProblem(c, vconstr, objSense, vdomains)

Constructs a `DefaultProblem` struct by providing the vector of cost coefficients `c`, the vector of constraint of the problem `vconstr`, the objective sense, and the vector of domains of the variables.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> OREnvironment.constructProblem(cost, constraints, :max, domain)
```
"""
function constructProblem(costs::Array{Float64,1},
                          constraints::Array{<:Constraint,1},
                          objSense::Symbol,
                          domain::Array{VariableDomain,1})
    variablesConstraints = get_relationship_variables_constraints(constraints, length(costs));
    return DefaultProblem(costs, constraints, variablesConstraints, objSense, domain);
end

"""
    get_lb_variable(p, var)

Returns the lower bound of variable `var`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_lb_variable(p, 1);
0.0
```
"""
get_lb_variable(p::Problem, variable::Int)::Float64 = get_lb(p.variablesDomain[variable]);

"""
    set_lb_variable!(p, var, val)

Set the lower bound of variable `var` to value `val`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.set_lb_variable!(p, 1, 1.0);
julia> OREnvironment.get_lb_variable(p, 1);
1.0
```
"""
@inline function set_lb_variable!(p::Problem, variable::Int, lb::Float64)
    set_lb!(p.variablesDomain[variable], lb);
end

"""
    get_ub_variable(p, var)

Returns the upper bound of variable `var`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_ub_variable(p, 1);
1.0
```
"""
get_ub_variable(p::Problem, variable::Int)::Float64 = get_ub(p.variablesDomain[variable]);

"""
    set_ub_variable!(p, var, val)

Set the upper bound of variable `var` to value `val`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.set_ub_variable!(p, 1, 3.0);
julia> OREnvironment.get_ub_variable(p, 1);
3.0
```
"""
@inline function set_ub_variable!(p::Problem, variable::Int, ub::Float64)
    set_ub!(p.variablesDomain[variable], ub);
end

"""
    is_value_within_the_domain(p, var, val)

Checks if for variable `var`, value `val` is within its domain.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.is_value_within_the_domain(p, 1, 3.0)
false
julia> OREnvironment.is_value_within_the_domain(p, 1, 0.5)
true
```
"""
@inline function is_value_within_the_domain(p::Problem,
                                            variable::Int,
                                            value::T)::Bool where {T<:Real}
    return is_value_within_the_domain(p.variablesDomain[variable], value)
end

"""
    get_middle_point_variable_domain(p, var)

Computes the middle point corresponding to the domain of variable `var`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_middle_point_variable_domain(p, 1);
0.5
```
"""
@inline function get_middle_point_variable_domain(p::Problem, variable::Int)::Float64
    return get_middle_point(p.variablesDomain[variable])
end

"""
    get_number_of_variables(p)

Returns the number of variables in problem `p`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_number_of_variables(p);
6
```
"""
get_number_of_variables(p::Problem)::Int = length(p.costs)

"""
    get_number_of_constraints(p)

Returns the number of constraints in problem `p`. This does not include variables bounds.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_number_of_constraints(p);
2
```
"""
get_number_of_constraints(p::Problem)::Int = length(p.constraints)

"""
    get_cost(p, var)

Returns the cost coefficient of variable `var`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_cost(p, 1);
1.0
```
"""
get_cost(p::Problem, variable::Int)::Float64 = @inbounds p.costs[variable];

"""
    set_cost!(p, var, coef)

Sets `coef` as cost coefficient of variable `var`.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.set_cost!(p, 1, 1.5);
julia> OREnvironment.get_cost(p, 1);
1.5
```
"""
set_cost!(p::Problem, variable::Int, value::Float64) = @inbounds p.costs[variable] = value;

"""
    get_obj_sense(p)

Returns the objetive sense of the problem.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_obj_sense(p);
:max
```
"""
get_obj_sense(p::Problem)::Symbol = p.objSense;

"""
    set_obj_sense!(p, objSense)

Sets `objSense` as the objetive sense of the problem.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.set_obj_sense!(p, :min);
julia> OREnvironment.get_obj_sense(p);
:min
```
"""
set_obj_sense!(p::Problem, value::Symbol) = p.objSense = value;

"""
    get_constraint(p, i)

Returns the i-th constraint of the problem.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_constraint(p, 1)
```
"""
get_constraint(p::Problem, idxConstraint::Int) = @inbounds p.constraints[idxConstraint];

"""
    get_constraints(p)

Returns the constraints of the problem.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_constraints(p)
```
"""
get_constraints(p::Problem) = p.constraints;

"""
    get_constraints_of_variable(p, var)

Returns the indexes of the constraints where variable `var` appears in the problem.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.get_constraints_of_variable(p, 1)
[1, 2] # variable 1 appears in both constraints
```
"""
get_constraints_of_variable(p::Problem, variable::Int)::Array{Int,1} = @inbounds p.variablesConstraints[variable];

"""
    add_constraint!(p, c)

Adds constraint `c` to problem `p`. It is also possible to provide a vector of constraints instead of a single constraint.

**Note:** The vector indicating which variables appears in which constraints is also updated here.

# Example
```jldoctest
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> c = OREnvironment.constructConstraint(3.9, :equal, [5, 6], [1.4, 11.0]);
julia> OREnvironment.add_constraint!(p, c);
```
"""
function add_constraint!(p::Problem, c::Constraint)
    push!(p.constraints, c);
    local idx::Int = length(p.constraints);
    add_constraint_index_to_variables!(c, idx, p.variablesConstraints);
end

function add_constraint!(p::Problem, constraints::Array{<:Constraint,1})
    for c in constraints
        add_constraint!(p,c);
    end
end

"""
    is_feasible(p, s)

Checks if the current status of solution `s` is feasible for problem `p`, regarding constraints and variables bounds.

**Note:** This function DOES NOT compute the current consumption of solution `s`. It just checks the last values saved in memory and compare them with the right-hand side values of the constraints.

# Example
```jldoctest
julia> typeVariables = Int; numVariables = 6;  numConstraints = 2;
julia> status = OREnvironment.constructStatus(numConstraints);
julia> args = (typeVariables, numVariables, status);
julia> s = OREnvironment.constructSolution(:FixedLengthArray, args)
julia> cost = collect(1.0:6.0);
julia> variables1 = [1, 3, 4, 6];
julia> variables2 = [1, 3, 5, 6];
julia> coefs1 = [2.3, 3.2, 3.1, 12.34];
julia> coefs2 = coefs1 .+ 1.0;
julia> constraint1 = OREnvironment.constructConstraint(15.0, :lessOrEq, variables1, coefs1);
julia> constraint2 = OREnvironment.constructConstraint(9.0, :lessOrEq, variables2, coefs2);
julia> constraints = [constraint1, constraint2];
julia> domain = [OREnvironment.VariableDomain(0.0,1.0) for i in 1:6];
julia> p = OREnvironment.constructProblem(cost, constraints, :max, domain);
julia> OREnvironment.is_feasible(p,s)
true
julia> OREnvironment.add_solution!(s, 1, 2) # value 2 violates domain
julia> OREnvironment.is_feasible(p,s)
false
```
"""
function is_feasible(p::Problem, s::Solution)::Bool
    local lhsFeasible::Bool = is_current_consumption_feasible(s, p.constraints);
    if lhsFeasible == false
        return false;
    end
    local boundFeasible::Bool = is_solution_within_bounds(p.variablesDomain, s);
    return boundFeasible;
end

"""
    gap(LB, UB, objSense)

Returns the gap between the lower bound `LB` and the upper bound `UB`. The computation changes depending on the objective sense of the problem.

# Example
```jldoctest
julia> LB = 7.0;
julia> UB = 9.0;
julia> OREnvironment.gap(LB, UB, :min) == (UB-LB)/LB
true
julia> OREnvironment.gap(LB, UB, :max) == (UB-LB)/UB
true
```
"""
gap(lowerBound::Float64, upperBound::Float64, objSense::Symbol)::Float64 = return gap(lowerBound, upperBound, Val(objSense));
gap(lowerBound::Float64, upperBound::Float64, ::Val{:min})::Float64 = return (upperBound-lowerBound)/lowerBound;
gap(lowerBound::Float64, upperBound::Float64, ::Val{:max})::Float64 = return (upperBound-lowerBound)/upperBound;


text:
# function simple_multiply(x, y)
#    product = 0
#    for i in 1:x
#        product += y
#    end
#    product
# end

# function simpledivide(dividend, divisor)
#     quotient = 0
#     remainder = dividend

#     while remainder - divisor >= 0
#         remainder -= divisor
#         quotient += 1
#     end
#     quotient, remainder
# end

addend = 0
accumulator = 0
counter = 0
shifter = 0

function clear()
    global addend = 0
    global accumulator = 0
    global counter = 0
    global shifter = 0
    nothing
end

function setaddend!(x)
   global addend = x
end

function leftshift!()
   global shifter += 1
end

function rightshift!()
   global shifter -= 1
end

function crank!(n)
    for _ in 1:n
        global accumulator += addend * 10^shifter
        global counter += 10^shifter
    end

    accumulator, counter
end

function uncrank!(n)
    for _ in 1:n
        global accumulator -= addend * 10^shifter
        global counter -= 10^shifter
    end

    accumulator, counter
end

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using HillsEigenvalues
using Test

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

max_stars_repo_name: Mikejmnez/HillsEigenvalues

max_issues_repo_name: Mikejmnez/HillsEigenvalues

max_forks_repo_name: Mikejmnez/HillsEigenvalues


text:
"maybe of type `T` or nothing"
const Maybe{T} = Union{Nothing, T}

"""
    PartialDefault{F}

Type for non-constant default value, it depends on the value of another field that has default value.
"""
struct PartialDefault{F}
    lambda::F
    vars::Vector{Symbol}
    expr::Expr
end

(f::PartialDefault)(x) = f.lambda(x)

function Base.show(io::IO, x::PartialDefault)
    print(io, JLFunction(;head=:->, args=x.vars, body=x.expr))
end

"""
    Reflect

Placeholder type for reflected type string.

# Type Alias

if the corresponding type has a [`type_alias`](@ref) defined,
serialization and parsing will use the [`type_alias`](@ref) instead
of the type name, this only works on concrete types since the alias
cannot contain any type var information.

# Example

the following option struct

```julia
@option struct MyOption
    type::Reflect
    name::String = "Sam"
end
```

would be equivalent to

```toml
type = "MyOption"
name = "Sam"
```

this is useful for defining list of different types etc.
"""
struct Reflect end

"""
    create(::Type{T}; kwargs...) where T

Create an instance of option type `T` from `kwargs`. Similar
to the default keyword argument constructor, but one can use this to create
custom keyword argument constructor with extra custom keywords.
"""
function create(::Type{T}; kwargs...) where T
    error("$T is not an option type")
end


text:
#
# Factors-DataFrames conversion
#
# Bridge between the two

"""
Convert a Factor to a DataFrame
"""
function Base.convert(::Type{DataFrame}, ϕ::Factor)
    df = DataFrames.DataFrame(pattern(ϕ))
    DataFrames.rename!(df, names(df), names(ϕ))
    df[:potential] = ϕ.potential[:]

    return df
end


text:
# Run the benchmarks for Julia (run.jl, and boid/run.jl) and Python (run.py and boid/run.py). Copy the results below and plot

using VegaLite
using DataFrames

## Agents.jl benchmark results

# forest fire
jlresults = [0.0639249, 0.142405601, 0.257741999, 0.371782101]

# boid
jlboid = [0.101076999, 0.270401899, 0.5469489, 0.8079563]

# Mesa benchmark results

#forest fire
pyresults = [0.8553307999998196, 2.0069307999999637, 3.3087123000000247, 4.781681599999956]

# boid
pyboid = [0.8770560000000387, 2.4145003999999517, 4.189664500000049, 6.677871200000027]

dd = DataFrame(runtime=vcat(pyresults./jlresults, pyboid./jlboid),
              model = vcat(fill("Forest fire", 4), fill("Boid flocking", 4)))

p = @vlplot(
  width=300, height=90,
  data = dd,
  mark = :circle,
  y = {"model:n", axis={title=""}},
  x = {:runtime, axis={title="Mesa/Agents run time"}},
  color = {"model:n", legend=false}
)
save("benchmark_mesa.svg", p)


text:
#=

  Magic squares in Julia ConstraintSolver.jl

  This models solves n=3:7 easily but it has problem from n=8

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function magic_square(n=4,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            # "branch_strategy" => :IMPS, # default
                                                            "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 10, # default, 10
                                                            "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            "time_limit"=> 6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    s = round(Int,n*(n^2 + 1) / 2)
    ss = s
    println("n:$n n*n:$(n^2) s:$s")
    @variable(model, 1 <= x[1:n,1:n] <= n^2, Int)
    # @variable(model, ss <= s <= ss, Int) # It works if s is a variable

    @constraint(model, x[:] in CS.AllDifferent())

    # rows and cols
    for i in 1:n
        # Rows
        # @constraint(model, s == sum(x[i,:]))  # This works
        # The following yield INFEASIBLE
        # I reported this in https://github.com/Wikunia/ConstraintSolver.jl/issues/226
        # This was fixed in ConstraintSolver.jl v0.6.0.
        @constraint(model, sum(x[i,:]) == s)

        # Columns
        # @constraint(model, sum(x[:,i]) == s)
        @constraint(model, s == sum(x[:,i]))
    end

    # diagonals
    # @constraint(model, s == sum([x[i,i] for i in 1:n]))
    @constraint(model, sum([x[i,i] for i in 1:n]) == s)

    # @constraint(model, s == sum([x[i,n-i+1] for i in 1:n]))
    @constraint(model, sum([x[i,n-i+1] for i in 1:n]) == s)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                println("s:$s")
                display(xx)

            end
        end
    else
        println("status:$status")
    end

    return status
end

print_solutions=true
all_solutions=false
for n in 3:10  # 8.. probably time out
    println("\nn:$n")
    @time status = magic_square(n,print_solutions,all_solutions)
    if status != MOI.OPTIMAL
        break
    end
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function chain_feedforward(chain::Chain, all_t::Array{Float64, 1})
  const L = chain.x_coords[end]
  const rho = 2.3
  const g = 9.81
  const p0 = rho*g*L #mass of the chain
  const M = 6.1511 #mass of the car
  const h = 0.5 #how far the cart will drive

  #initialise outputs
  a0_target = 0.0 #acceleration at the top
  F_ff = Array{Float64}(length(all_t))
  xL_target = Array{Float64}(length(all_t)) #position of the bottom end of the chain (flacher Ausgang)
  x0_target = Array{Float64}(length(all_t)) #position of the car/top of the chain
  v0_target = Array{Float64}(length(all_t)) #velocity of the car/top of the chain
  phi0_target = Array{Float64}(length(all_t)) #chain angle at top of the chain
  phi0_dot_target = Array{Float64}(length(all_t)) #angular velocity of the chain at top of the chain

  poly1_coeffs, poly2_coeffs  = xL_poly_coeffs(h)

  #calculate (?flacher Ausgang?)
  for i = 1:length(all_t)
    xL_target[i] = y_w(all_t[i], poly1_coeffs, poly2_coeffs)
    x0_target[i] = x0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    v0_target[i] = v0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    phi0_target[i] = phi0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    phi0_dot_target[i] = dphi0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)

    F_ff[i] = M*a0_w(all_t[i], poly1_coeffs, poly2_coeffs, L) - p0*sin(phi0_target[i])
  end

  return F_ff, xL_target, x0_target, v0_target, phi0_target, phi0_dot_target
end

function chain_feedforward(all_t::Array{Float64, 1})
  const L = 1.178
  const rho = 2.3
  const g = 9.81
  const p0 = rho*g*L #mass of the chain
  const M = 6.1511 #mass of the car
  const h = 0.5 #how far the cart will drive

  #initialise outputs
  a0_target = 0.0 #acceleration at the top
  F_ff = Array{Float64}(length(all_t))
  xL_target = Array{Float64}(length(all_t)) #position of the bottom end of the chain (flacher Ausgang)
  x0_target = Array{Float64}(length(all_t)) #position of the car/top of the chain
  v0_target = Array{Float64}(length(all_t)) #velocity of the car/top of the chain
  phi0_target = Array{Float64}(length(all_t)) #chain angle at top of the chain
  phi0_dot_target = Array{Float64}(length(all_t)) #angular velocity of the chain at top of the chain

  poly1_coeffs, poly2_coeffs  = xL_poly_coeffs(h)

  #calculate (?flacher Ausgang?)
  for i = 1:length(all_t)
    xL_target[i] = y_w(all_t[i], poly1_coeffs, poly2_coeffs)
    x0_target[i] = x0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    v0_target[i] = v0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    phi0_target[i] = phi0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)
    phi0_dot_target[i] = dphi0_w(all_t[i], poly1_coeffs, poly2_coeffs, L)

    F_ff[i] = M*a0_w(all_t[i], poly1_coeffs, poly2_coeffs, L) - p0*sin(phi0_target[i])
  end

  return F_ff, xL_target, x0_target, v0_target, phi0_target, phi0_dot_target
end

#calculate polynomial coefficients for the (?flacher Ausgang?)
function xL_poly_coeffs(h)
  const t1 = 0.7
  const t2 = t1+0.3465
  const t3 = 3.0
  const t4 = t3+0.3465

  const c4 = -2*(t1+t2)
  const c3 = (t1^2+t2^2+4*t1*t2)
  const c2 = -2*t1*t2*(t1+t2)
  const c1 = (t1*t2)^2
  poly1(t) = 1/5*t^5 + 1/4*t^4*c4 + 1/3*t^3*c3 + 1/2*t^2*c2 + t*c1 #integrated polynomial

  const k4 = -2*(t3+t4)
  const k3 = (t3^2+t4^2+4*t3*t4)
  const k2 = -2*t3*t4*(t3+t4)
  const k1 = (t3*t4)^2
  poly2(t) = 1/5*t^5 + 1/4*t^4*k4 + 1/3*t^3*k3 + 1/2*t^2*k2 + t*k1

  Out1 = Array{Float64}(6)
  Out1[1] = -poly1(t1)/(poly1(t2)-poly1(t1))
  Out1[2] = c1/(poly1(t2)-poly1(t1))
  Out1[3] = 1/2*c2/(poly1(t2)-poly1(t1))
  Out1[4] = 1/3*c3/(poly1(t2)-poly1(t1))
  Out1[5] = 1/4*c4/(poly1(t2)-poly1(t1))
  Out1[6] = 1/5/(poly1(t2)-poly1(t1))

  Out2 = Array{Float64}(6)
  Out2[1] = 1.0+poly2(t3)/(poly2(t4)-poly2(t3))
  Out2[2] = -k1/(poly2(t4)-poly2(t3))
  Out2[3] = -1/2*k2/(poly2(t4)-poly2(t3))
  Out2[4] = -1/3*k3/(poly2(t4)-poly2(t3))
  Out2[5] = -1/4*k4/(poly2(t4)-poly2(t3))
  Out2[6] = -1/5/(poly2(t4)-poly2(t3))

  return Out1*h, Out2*h
end

#Solltrajektorie
function y_w(t, poly1_coeffs, poly2_coeffs)
  const t1 = 0.7
  const t2 = t1+0.3465
  const t3 = 3.0
  const t4 = t3+0.3465
  const t_period = 5.0

  poly1(t) = poly1_coeffs[1] + poly1_coeffs[2]*t + poly1_coeffs[3]*t^2+ poly1_coeffs[4]*t^3+ poly1_coeffs[5]*t^4 + poly1_coeffs[6]*t^5
  poly2(t) = poly2_coeffs[1] + poly2_coeffs[2]*t + poly2_coeffs[3]*t^2+ poly2_coeffs[4]*t^3+ poly2_coeffs[5]*t^4 + poly2_coeffs[6]*t^5

  t = mod(t, t_period)
  if t<=t1
    Out = poly1(t1)
  elseif t1<t<=t2
    Out = poly1(t)
  elseif t2<t<=t3
    Out = poly1(t2)
  elseif t3<t<=t4
    Out = poly2(t)
  else
    Out = poly2(t4)
  end
  return Out::Float64
end

function dydt_w(t, poly1_coeffs, poly2_coeffs)
  const t1 = 0.7
  const t2 = t1+0.3465
  const t3 = 3.0
  const t4 = t3+0.3465
  const t_period = 5.0

  poly1(t) = poly1_coeffs[2] + 2*poly1_coeffs[3]*t+ 3*poly1_coeffs[4]*t^2+ 4*poly1_coeffs[5]*t^3 + 5*poly1_coeffs[6]*t^4
  poly2(t) = poly2_coeffs[2] + 2*poly2_coeffs[3]*t+ 3*poly2_coeffs[4]*t^2+ 4*poly2_coeffs[5]*t^3 + 5*poly2_coeffs[6]*t^4

  t = mod(t, t_period)
  if t<=t1
    Out = 0.0
  elseif t1<t<=t2
    Out = poly1(t)
  elseif t2<t<=t3
    Out = 0.0
  elseif t3<t<=t4
    Out = poly2(t)
  else
    Out = 0.0
  end
  return Out::Float64
end

function ddydt2_w(t, poly1_coeffs, poly2_coeffs)
  const t1 = 0.7
  const t2 = t1+0.3465
  const t3 = 3.0
  const t4 = t3+0.3465
  const t_period = 5.0

  poly1(t) = 2*poly1_coeffs[3] + 6*poly1_coeffs[4]*t+ 12*poly1_coeffs[5]*t^2 + 20*poly1_coeffs[6]*t^3
  poly2(t) = 2*poly2_coeffs[3] + 6*poly2_coeffs[4]*t+ 12*poly2_coeffs[5]*t^2 + 20*poly2_coeffs[6]*t^3

  t = mod(t, t_period)
  if t<=t1
    Out = 0.0
  elseif t1<t<=t2
    Out = poly1(t)
  elseif t2<t<=t3
    Out = 0.0
  elseif t3<t<=t4
    Out = poly2(t)
  else
    Out = 0.0
  end
  return Out::Float64
end

function x0_w(t::Number, poly1_coeffs::Array{Float64, 1}, poly2_coeffs::Array{Float64, 1}, L::Number)
  const g = 9.81 #N/kg acceleration through gravity
  #return 1/pi*quadgk(x->y_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, abstol = 1e-8)[1] #integral
  return 1/pi*integral_simpson(x->y_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, 300) #integral
end

function v0_w(t::Number, poly1_coeffs::Array{Float64, 1}, poly2_coeffs::Array{Float64, 1}, L::Number)
  const g = 9.81 #N/kg acceleration through gravity
  #return 1/pi*quadgk(x->dydt_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, abstol = 1e-8)[1] #integral
  return 1/pi*integral_simpson(x->dydt_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, 300) #integral
end

function a0_w(t::Number, poly1_coeffs::Array{Float64, 1}, poly2_coeffs::Array{Float64, 1}, L::Number)
  const g = 9.81 #N/kg acceleration through gravity
  #return 1/pi*quadgk(x->ddydt2_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, abstol = 1e-8)[1] #integral
  return 1/pi*integral_simpson(x->ddydt2_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs), -pi/2, pi/2, 300) #integral
end

function phi0_w(t::Number, poly1_coeffs::Array{Float64, 1}, poly2_coeffs::Array{Float64, 1}, L::Number)
  const g = 9.81 #N/kg acceleration through gravity
  #return 1/pi*quadgk(x->(-dydt_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs)*sin(x)/(sqrt(L*g))), -pi/2, pi/2, abstol = 1e-8)[1] #integral
  return 1/pi*integral_simpson(x->(-dydt_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs)*sin(x)/(sqrt(L*g))), -pi/2, pi/2, 300) #integral
end

function dphi0_w(t::Number, poly1_coeffs::Array{Float64, 1}, poly2_coeffs::Array{Float64, 1}, L::Number)
  const g = 9.81 #N/kg acceleration through gravity
  #return 1/pi*quadgk(x->(-ddydt2_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs)*sin(x)/(sqrt(L*g))), -pi/2, pi/2, abstol = 1e-8)[1] #integral
  return 1/pi*integral_simpson(x->(-ddydt2_w(t+2*sqrt(L/g)*sin(x), poly1_coeffs, poly2_coeffs)*sin(x)/(sqrt(L*g))), -pi/2, pi/2, 300) #integral
end

#source: https://rosettacode.org/wiki/Numerical_integration
function integral_simpson(f::Function, xmin::Number, xmax::Number, N::Int)
  h = (xmax-xmin)/N
  sum1 = f(xmin+h/2)
  sum2 = 0.0
  for i = 1:(N-1)
    sum1 += f(xmin+h*(i+0.5))
    sum2 += f(xmin+h*i)
  end
  return h/6*(f(xmin)+f(xmax)+4*sum1+2*sum2)
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA