text:
@doc raw"""
    energy_balance(psi_container::PSIContainer,
                        initial_conditions::Vector{InitialCondition},
                        efficiency_data::Tuple{Vector{String}, Vector{InOut}},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})

Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple

# Constraints

If t = 1:

``` varenergy[name, 1] == initial_conditions[ix].value + varin[name, 1]*eff_in*fraction_of_hour - varout[name, 1]*fraction_of_hour/eff_out ```

If t > 1:

``` varenergy[name, t] == varenergy[name, t-1] + varin[name, t]*eff_in*fraction_of_hour - varout[name, t]*fraction_of_hour/eff_out ```

# LaTeX

`` x^{energy}_1 == x^{energy}_{init} + frhr \eta^{in} x^{in}_1 - \frac{frhr}{\eta^{out}} x^{out}_1, \text{ for } t = 1 ``

`` x^{energy}_t == x^{energy}_{t-1} + frhr \eta^{in} x^{in}_t - \frac{frhr}{\eta^{out}} x^{out}_t, \forall t \geq 2 ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* initial_conditions::Vector{InitialCondition} : for time zero 'varenergy'
* efficiency_data::Tuple{Vector{String}, Vector{InOut}} :: charging/discharging efficiencies
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_names[1] : varin
- : var_names[2] : varout
- : var_names[3] : varenergy

"""
function energy_balance(
    psi_container::PSIContainer,
    initial_conditions::Vector{InitialCondition},
    efficiency_data::Tuple{Vector{String}, Vector{InOut}},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
)
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / 60
    name_index = efficiency_data[1]

    varin = get_variable(psi_container, var_names[1])
    varout = get_variable(psi_container, var_names[2])
    varenergy = get_variable(psi_container, var_names[3])

    constraint = add_cons_container!(psi_container, cons_name, name_index, time_steps)

    for (ix, name) in enumerate(name_index)
        eff_in = efficiency_data[2][ix].in
        eff_out = efficiency_data[2][ix].out

        constraint[name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            varenergy[name, 1] ==
            initial_conditions[ix].value + varin[name, 1] * eff_in * fraction_of_hour -
            (varout[name, 1]) * fraction_of_hour / eff_out
        )

    end

    for t in time_steps[2:end], (ix, name) in enumerate(name_index)
        eff_in = efficiency_data[2][ix].in
        eff_out = efficiency_data[2][ix].out

        constraint[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            varenergy[name, t] ==
            varenergy[name, t - 1] + varin[name, t] * eff_in * fraction_of_hour -
            (varout[name, t]) * fraction_of_hour / eff_out
        )
    end

    return
end

@doc raw"""
    energy_balance_external_input(psi_container::PSIContainer,
                        initial_conditions::Vector{InitialCondition},
                        efficiency_data::Tuple{Vector{String}, Vector{InOut}},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})
Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple
# Constraints
If t = 1:
``` varenergy[name, 1] == initial_conditions[ix].value + (paraminflow[name, t] - varspill[name, 1] - varout[name, 1])*fraction_of_hour ```
If t > 1:
``` varenergy[name, t] == varenergy[name, t-1] + (paraminflow[name, t] - varspill[name, t] - varout[name, t])*fraction_of_hour ```
# LaTeX
`` x^{energy}_1 == x^{energy}_{init} + frhr  (x^{in}_1 - x^{spillage}_1 -  x^{out}_1), \text{ for } t = 1 ``
`` x^{energy}_t == x^{energy}_{t-1} + frhr (x^{in}_t - x^{spillage}_t - x^{out}_t), \forall t \geq 2 ``
# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* initial_conditions::Vector{InitialCondition} : for time zero 'varenergy'
* inflow_data::Vector{DeviceTimeSeries} :: Inflow energy forecast information
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_names[1] : varspill
- : var_names[2] : varout
- : var_names[3] : varenergy
* param_reference::UpdateRef : UpdateRef to access the inflow parameter
"""
function energy_balance_external_input_param(
    psi_container::PSIContainer,
    initial_conditions::Vector{InitialCondition},
    inflow_data::Vector{DeviceTimeSeries},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
    param_reference::UpdateRef,
)
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / 60
    name_index = (d.name for d in inflow_data)

    varspill = get_variable(psi_container, var_names[1])
    varout = get_variable(psi_container, var_names[2])
    varenergy = get_variable(psi_container, var_names[3])
    container = add_param_container!(psi_container, param_reference, name_index, time_steps)
    paraminflow = get_parameter_array(container)
    multiplier = get_multiplier_array(container)
    constraint = add_cons_container!(psi_container, cons_name, name_index, time_steps)

    for (ix, d) in enumerate(inflow_data)
        multiplier[d.name, 1] = d.multiplier
        paraminflow[d.name, 1] = PJ.add_parameter(psi_container.JuMPmodel, d.timeseries[1])
        exp =
            initial_conditions[ix].value +
            (
                multiplier[d.name, 1] * paraminflow[d.name, 1] - varspill[d.name, 1] -
                varout[d.name, 1]
            ) * fraction_of_hour
        constraint[d.name, 1] =
            JuMP.@constraint(psi_container.JuMPmodel, varenergy[d.name, 1] == exp)

        for t in time_steps[2:end]
            paraminflow[d.name, t] =
                PJ.add_parameter(psi_container.JuMPmodel, d.timeseries[t])
            exp =
                varenergy[d.name, t - 1] +
                (
                    d.multiplier * paraminflow[d.name, 1] - varspill[d.name, t] -
                    varout[d.name, t]
                ) * fraction_of_hour
            constraint[d.name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, varenergy[d.name, t] == exp)
        end
    end
    return
end

@doc raw"""
    energy_balance_external_input(psi_container::PSIContainer,
                        initial_conditions::Vector{InitialCondition},
                        efficiency_data::Tuple{Vector{String}, Vector{InOut}},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})
Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple
# Constraints
If t = 1:
``` varenergy[name, 1] == initial_conditions[ix].value + (paraminflow[name, t] - varspill[name, 1] - varout[name, 1])*fraction_of_hour ```
If t > 1:
``` varenergy[name, t] == varenergy[name, t-1] + (paraminflow[name, t] - varspill[name, t] - varout[name, t])*fraction_of_hour ```
# LaTeX
`` x^{energy}_1 == x^{energy}_{init} + frhr  (x^{in}_1 - x^{spillage}_1 -  x^{out}_1), \text{ for } t = 1 ``
`` x^{energy}_t == x^{energy}_{t-1} + frhr (x^{in}_t - x^{spillage}_t - x^{out}_t), \forall t \geq 2 ``
# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* initial_conditions::Vector{InitialCondition} : for time zero 'varenergy'
* inflow_data::TVector{DeviceTimeSeries} :: Inflow energy forecast information
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_names[1] : varspill
- : var_names[2] : varout
- : var_names[3] : varenergy
"""
function energy_balance_external_input(
    psi_container::PSIContainer,
    initial_conditions::Vector{InitialCondition},
    inflow_data::Vector{DeviceTimeSeries},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
)
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / 60
    name_index = (d.name for d in inflow_data)

    varspill = get_variable(psi_container, var_names[1])
    varout = get_variable(psi_container, var_names[2])
    varenergy = get_variable(psi_container, var_names[3])

    constraint = add_cons_container!(psi_container, cons_name, name_index, time_steps)

    for (ix, d) in enumerate(inflow_data)
        constraint[d.name, 1] = JuMP.@constraint(
            psi_container.JuMPmodel,
            varenergy[d.name, 1] ==
            initial_conditions[ix].value +
            (d.multiplier * d.timeseries[1] - varspill[d.name, 1] - varout[d.name, 1]) *
            fraction_of_hour
        )

        for t in time_steps[2:end]
            constraint[d.name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                varenergy[d.name, t] ==
                varenergy[d.name, t - 1] +
                (d.multiplier * d.timeseries[t] - varspill[d.name, t] - varout[d.name, t]) * fraction_of_hour
            )
        end
    end
    return
end


text:
function generic_tests_ccmatrix(C)
    r = length(C)
    for i = 1:r
        ccm = SymbolicWedderburn.CMMatrix(C, i)
        l = length(C[i])
        @test all([sum(ccm[:,j]) == l for j = 1:r])
    end
    ccm1 = SymbolicWedderburn.CMMatrix(C, 1)
    @test all([isone(ccm1[i,i]) for i = 1:r])
end

@testset "ConjClassMatrix" begin
    @testset "example: Symmetric group (4)" begin
        G = PermutationGroups.SymmetricGroup(4)
        C = conjugacy_classes(G)
        generic_tests_ccmatrix(C)
        @test SymbolicWedderburn.CMMatrix(C, 2)  == [0 1 0 0 0;
                                                     6 0 3 0 2;
                                                     0 4 0 4 0;
                                                     0 0 3 0 4;
                                                     0 1 0 2 0]
        @test SymbolicWedderburn.CMMatrix(C, 3)  == [0 0 1 0 0;
                                                     0 4 0 4 0;
                                                     8 0 4 0 8;
                                                     0 4 0 4 0;
                                                     0 0 3 0 0]
        @test SymbolicWedderburn.CMMatrix(C, 4)  == [0 0 0 1 0;
                                                     0 0 3 0 4;
                                                     0 4 0 4 0;
                                                     6 0 3 0 2;
                                                     0 2 0 1 0]
        @test SymbolicWedderburn.CMMatrix(C, 5)  == [0 0 0 0 1;
                                                     0 1 0 2 0;
                                                     0 0 3 0 0;
                                                     0 2 0 1 0;
                                                     3 0 0 0 2]
    end
    @testset "example: Alternating group (4)" begin
        a = perm"(1,2,3)(4)"
        b = perm"(1,2,4)"
        G = PermGroup([a,b])
        C = [
            Orbit([a,b], one(G)),
            Orbit([a,b], perm"(1,2)(3,4)"),
            Orbit([a,b], perm"(1,2,3)(4)"),
            Orbit([a,b], perm"(1,3,2)(4)"),
            ]

        generic_tests_ccmatrix(C)
        @test SymbolicWedderburn.CMMatrix(C, 1) == [1 0 0 0
                                                    0 1 0 0
                                                    0 0 1 0
                                                    0 0 0 1]
        @test SymbolicWedderburn.CMMatrix(C, 2) == [0 1 0 0
                                                    3 2 0 0
                                                    0 0 3 0
                                                    0 0 0 3]
        @test SymbolicWedderburn.CMMatrix(C, 3) == [0 0 1 0
                                                    0 0 3 0
                                                    0 0 0 4
                                                    4 4 0 0]
        @test SymbolicWedderburn.CMMatrix(C, 4) == [0 0 0 1
                                                    0 0 0 3
                                                    4 4 0 0
                                                    0 0 4 0]
    end

    @testset "example: C₅ ⋊ C₄" begin
        S = [perm"(1,2,4,5,3)", perm"(2,5,3,4)"]
        G = PermGroup(S);

        ccG = [
            Orbit(S, one(G)),
            Orbit(S, perm"(2,3)(4,5)"),
            Orbit(S, perm"(2,4,3,5)"),
            Orbit(S, perm"(2,5,3,4)"),
            Orbit(S, perm"(1,2,4,5,3)"),
            ]
        # the order of cclasses is taken from GAP

        @assert sum(length, ccG) == order(G)

        @test SymbolicWedderburn.CMMatrix(ccG, 1) ==
        [ 1  0  0  0  0
          0  1  0  0  0
          0  0  1  0  0
          0  0  0  1  0
          0  0  0  0  1 ]
        @test SymbolicWedderburn.CMMatrix(ccG, 2) ==
        [ 0  1  0  0  0
          5  0  0  0  5
          0  0  0  5  0
          0  0  5  0  0
          0  4  0  0  0 ]
        @test SymbolicWedderburn.CMMatrix(ccG, 3) ==
        [ 0  0  1  0  0
          0  0  0  5  0
          0  5  0  0  0
          5  0  0  0  5
          0  0  4  0  0 ]
        @test SymbolicWedderburn.CMMatrix(ccG, 4) ==
        [ 0  0  0  1  0
          0  0  5  0  0
          5  0  0  0  5
          0  5  0  0  0
          0  0  0  4  0 ]
        @test SymbolicWedderburn.CMMatrix(ccG, 5) ==
        [ 0  0  0  0  1
          0  4  0  0  0
          0  0  4  0  0
          0  0  0  4  0
          4  0  0  0  3 ]

        generic_tests_ccmatrix(ccG)
        generic_tests_ccmatrix(conjugacy_classes(G)) # might be in different order
    end

    @testset "random subgroups of SymetricGroup(N)" begin
        for i in 2:6
            G = PermutationGroups.SymmetricGroup(i)
            for _ in 1:5
                PG = PermGroup(rand(G, 2))
                generic_tests_ccmatrix(SymbolicWedderburn.conjugacy_classes(PG))
            end
        end
    end
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.02"

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
end

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
md"##### Immutable:"

# ╔═╡ f9d73026-2386-11eb-2bdc-bde5f34e7694
begin
	a1 = 10
	b1 = a1
	b1 = 20
	a1
end

# ╔═╡ f9d759ca-2386-11eb-269c-2540e0b8d985
md"##### No copy:"

# ╔═╡ f9d7db52-2386-11eb-2170-3be42892f232
begin
	a2 = [10]
	b2 = a2
	b2[1] = 20
	a2
end

# ╔═╡ f9df3fa0-2386-11eb-3b6c-fd77875fc054
md"##### Copy:"

# ╔═╡ f9ea5298-2386-11eb-2ba2-3fe09f34fc4c
begin
	a3 = [10]
	b3 = copy(a3)
	b3[1] = 20
	a3
end

# ╔═╡ f9ec6ec8-2386-11eb-2d82-edb2e6c6fa76
md"##### Shallow copy:"

# ╔═╡ f9f2877c-2386-11eb-256d-19b1583fa1b0
begin
	a4 = [[10]]
	b4 = copy(a4)
	b4[1][1] = 20
	a4
end

# ╔═╡ f9f83140-2386-11eb-3c17-7b680d570db8
md"##### Deep copy:"

# ╔═╡ fa03c456-2386-11eb-0265-576e2afde307
begin
	a5 = [[10]]
	b5 = deepcopy(a5)
	b5[1][1] = 20
	a5
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.02"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╟─fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═f9d73026-2386-11eb-2bdc-bde5f34e7694
# ╟─f9d759ca-2386-11eb-269c-2540e0b8d985
# ╠═f9d7db52-2386-11eb-2170-3be42892f232
# ╟─f9df3fa0-2386-11eb-3b6c-fd77875fc054
# ╠═f9ea5298-2386-11eb-2ba2-3fe09f34fc4c
# ╟─f9ec6ec8-2386-11eb-2d82-edb2e6c6fa76
# ╠═f9f2877c-2386-11eb-256d-19b1583fa1b0
# ╟─f9f83140-2386-11eb-3c17-7b680d570db8
# ╠═fa03c456-2386-11eb-0265-576e2afde307
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
# calculate covariance

module CalcCovariance
    using DataFrames, CSV, Statistics

    function main()
        # input
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_700.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        # extract between 12 and 16 o'clock
        df_ext = df_org[(df_org.time.>=120000).&(df_org.time.<160000),:]

        # calculate
        println("Variance of ir: $(var(df_ext.ir))")
        println("Variance of lidar: $(var(df_ext.lidar))")

        diff_ir = df_ext.ir .- mean(df_ext.ir)
        diff_lidar = df_ext.lidar .- mean(df_ext.lidar)

        println("Mean of ir: $(mean(df_ext.ir))")
        println("Mean of lidar: $(mean(df_ext.lidar))")

        a = diff_ir .* diff_lidar
        println("Covariance: $(sum(a)/size(df_ext)[1])")
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module _EllipticCurves

using Base.Test
using Brainstorm
using Nemo

function test_wnf_ZZ()
    @test curve(17, 1) == Brainstorm._Math._EllipticCurves.WNFCurve(17, 1)
    @test_throws ErrorException curve(0, 0)

    c = curve(1, 2)
    p = point(c, 1, 2)
    i = ideal(c)
    @test p == Brainstorm._Math._EllipticCurves.ConcretePoint{typeof(c), Int}(c, 1, 2)
    @test_throws ErrorException point(c, 1, 1)
    @test i == Brainstorm._Math._EllipticCurves.IdealPoint{typeof(c), Int}(c, 0, 0)

    @test -p == Brainstorm._Math._EllipticCurves.ConcretePoint{typeof(c), Int}(c, 1, -2)
    @test -i == i

    @test i + i == i
    @test p + i == p
    @test i + p == p
    @test p + (-p) == i
    @test p - p == i

    @test point(c, 1) == p
end

function test_wnf_QQ()
    c = curve(big(-2//1), big(4//1))
    p = point(c, 3, 5)
    q = point(c, -2, 0)

    @test p-q == point(c, 0, -2)
    @test p+p+p == point(c, -237//121, 845 // 1331)
    @test p+p+p+p+p == point(c, 2312883//1142761, -3507297955//1221611509)
    @test 5p == point(c, 2312883//1142761, -3507297955//1221611509)
    @test q-3p == point(c, 240, 3718)
    @test -20p == point(c,
        872171688955240345797378940145384578112856996417727644408306502486841054959621893457430066791656001//520783120481946829397143140761792686044102902921369189488390484560995418035368116532220330470490000,
        -27483290931268103431471546265260141280423344817266158619907625209686954671299076160289194864753864983185162878307166869927581148168092234359162702751//11884621345605454720092065232176302286055268099954516777276277410691669963302621761108166472206145876157873100626715793555129780028801183525093000000)

    @test point(c, big(3//1)) == p
    @test point(c, big(-2//1)) == q
end

function test_wnf_ZZp()
    R = ResidueRing(ZZ, 5)

    c = curve(R(1), R(1))
    p = point(c, R(2), R(1))

    @test 2p == point(c, R(2), R(4))
    @test 3p == ideal(c)

    @test Brainstorm._Math._EllipticCurves.logpx(p, 2) == (Dict(R(2) => 1), 0, point(c, R(2), R(4)))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 3) == (Dict(R(2) => 1), 3, ideal(c))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 4) == (Dict(R(2) => 1), 3, ideal(c))
end

function test_wnf_polyZZp()
    R = ResidueRing(ZZ, 5)
    S, s = PolynomialRing(R, "s")
    T, t = FiniteField(3 + s^2, "t")

    c = curve(T(1), T(1))
    p = point(c, 2 + t, 2t)

    @test  -p == point(c, 2 +  t,     3t)
    @test p+p == point(c, 3 +  t, 2 + 0t)
    @test  4p == point(c, 3 + 2t, 4 + 4t)
    @test  9p == ideal(c)

    @test Brainstorm._Math._EllipticCurves.logpx(p, 8)  == (Dict((i*p).x => i for i in 1:4), 0, 8p)
    @test Brainstorm._Math._EllipticCurves.logpx(p, 9)  == (Dict((i*p).x => i for i in 1:4), 9, ideal(c))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 10) == (Dict((i*p).x => i for i in 1:4), 9, ideal(c))
end

function test_wnf_group()
    T, t = FiniteField(37, 1, "t")
    c = curve(T(-1), T(3))
    N = 42
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(point(c, T(2), T(3)), N) == 7
    @test order(gen(c, N, 3),         N) == 3
    @test order(gen(c, N, 7),         N) == 7

    T, t = FiniteField(97, 1, "t")
    c = curve(T(2), T(3))
    N = 100
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(gen(c, N, 5), N) == 5

    T, t = FiniteField(29, 1, "t")
    c = curve(T(-1), T(1))
    N = 37
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(gen(c, N, 37), N) == 37
end

function test_gwnf_ZZ()
    @test curve(0, 0, 0, 17, 1) == Brainstorm._Math._EllipticCurves.GWNFCurve(0, 0, 0, 17, 1)
    @test_throws ErrorException curve(0, 0, 0, 0, 0)

    c = curve(0, 0, 0, 1, 2)
    p = point(c, 1, 2)
    i = ideal(c)
    @test p == Brainstorm._Math._EllipticCurves.ConcretePoint{typeof(c), Int}(c, 1, 2)
    @test_throws ErrorException point(c, 1, 1)
    @test i == Brainstorm._Math._EllipticCurves.IdealPoint{typeof(c), Int}(c, 0, 0)

    @test -p == Brainstorm._Math._EllipticCurves.ConcretePoint{typeof(c), Int}(c, 1, -2)
    @test -i == i

    @test i + i == i
    @test p + i == p
    @test i + p == p
    @test p + (-p) == i
    @test p - p == i

    #@test point(c, 1) == p
end

function test_gwnf_QQ()
    c = curve(0, 0, 0, big(-2//1), big(4//1))
    p = point(c, 3, 5)
    q = point(c, -2, 0)

    @test p-q == point(c, 0, -2)
    @test p+p+p == point(c, -237//121, 845 // 1331)
    @test p+p+p+p+p == point(c, 2312883//1142761, -3507297955//1221611509)
    @test 5p == point(c, 2312883//1142761, -3507297955//1221611509)
    @test q-3p == point(c, 240, 3718)
    @test -20p == point(c,
        872171688955240345797378940145384578112856996417727644408306502486841054959621893457430066791656001//520783120481946829397143140761792686044102902921369189488390484560995418035368116532220330470490000,
        -27483290931268103431471546265260141280423344817266158619907625209686954671299076160289194864753864983185162878307166869927581148168092234359162702751//11884621345605454720092065232176302286055268099954516777276277410691669963302621761108166472206145876157873100626715793555129780028801183525093000000)

    #@test point(c, big(3//1)) == p
    #@test point(c, big(-2//1)) == q
end

function test_gwnf_ZZp()
    R = ResidueRing(ZZ, 5)

    c = curve(R(0), R(0), R(0), R(1), R(1))
    p = point(c, R(2), R(1))

    @test 2p == point(c, R(2), R(4))
    @test 3p == ideal(c)

    @test Brainstorm._Math._EllipticCurves.logpx(p, 2) == (Dict(R(2) => 1), 0, point(c, R(2), R(4)))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 3) == (Dict(R(2) => 1), 3, ideal(c))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 4) == (Dict(R(2) => 1), 3, ideal(c))
end

function test_gwnf_polyZZp()
    R = ResidueRing(ZZ, 5)
    S, s = PolynomialRing(R, "s")
    T, t = FiniteField(3 + s^2, "t")

    c = curve(T(0), T(0), T(0), T(1), T(1))
    p = point(c, 2 + t, 2t)

    @test  -p == point(c, 2 +  t,     3t)
    @test p+p == point(c, 3 +  t, 2 + 0t)
    @test  4p == point(c, 3 + 2t, 4 + 4t)
    @test  9p == ideal(c)

    @test Brainstorm._Math._EllipticCurves.logpx(p, 8)  == (Dict((i*p).x => i for i in 1:4), 0, 8p)
    @test Brainstorm._Math._EllipticCurves.logpx(p, 9)  == (Dict((i*p).x => i for i in 1:4), 9, ideal(c))
    @test Brainstorm._Math._EllipticCurves.logpx(p, 10) == (Dict((i*p).x => i for i in 1:4), 9, ideal(c))
end

#=
function test_gwnf_group()
    T, t = FiniteField(37, 1, "t")
    c = curve(T(0), T(0), T(0), T(-1), T(3))
    N = 42
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(point(c, T(2), T(3)), N) == 7
    @test order(gen(c, N, 3),         N) == 3
    @test order(gen(c, N, 7),         N) == 7

    T, t = FiniteField(97, 1, "t")
    c = curve(T(0), T(0), T(0), T(2), T(3))
    N = 100
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(gen(c, N, 5), N) == 5

    T, t = FiniteField(29, 1, "t")
    c = curve(T(0), T(0), T(0), T(-1), T(1))
    N = 37
    @test order(c, _GroupAlgorithm.BabyStepGiantStep) == N
    @test order(gen(c, N, 37), N) == 37
end
=#

function test_all()
    print(rpad("Math.EllipticCurves...", 50, ' '))

    test_wnf_ZZ()
    test_wnf_QQ()
    test_wnf_ZZp()
    test_wnf_polyZZp()
    test_wnf_group()

    test_gwnf_ZZ()
    test_gwnf_QQ()
    test_gwnf_ZZp()
    test_gwnf_polyZZp()
    #test_gwnf_group()

    println("PASS")
end

end


text:
function computePrecessionNutation(ttt, Ax, Ay, As)

    # Define Constants
    ARCSEC2RAD = π / 64800.0
    MUARCSEC2RAD = π / 648000000000

    # Precompute powers of ttt
    ttt2 = ttt*ttt
    ttt3 = ttt2*ttt
    ttt4 = ttt3*ttt

    # Compute Delaunay Arguments
    DelArgs = @SVector[
        rem((485868.249036 + 1717915923.2178*ttt + 31.8792*ttt2 + 0.051635*ttt3 - 0.00024470*ttt4)*ARCSEC2RAD,2*pi),
        rem((1287104.79305 + 129596581.0481*ttt - 0.5532*ttt2 + 0.000136*ttt3 - 0.00001149*ttt4)*ARCSEC2RAD,2*pi),
        rem((335779.526232 + 1739527262.8478*ttt - 12.7512*ttt2 - 0.001037*ttt3 + 0.00000417*ttt4)*ARCSEC2RAD,2*pi),
        rem((1072260.70369 + 16029616001.2090*ttt - 6.3706*ttt2 + 0.006593*ttt3 - 0.00003169*ttt4)*ARCSEC2RAD,2*pi),
        rem((450160.398036 - 6962890.5431*ttt + 7.4722*ttt2 + 0.007702*ttt3 - 0.00005939*ttt4)*ARCSEC2RAD,2*pi),
        rem(4.402608842 + 2608.7903141574*ttt,2*pi),
        rem(3.176146698 + 1021.3285546211*ttt,2*pi),
        rem(1.753470314 + 628.3075849991*ttt,2*pi),
        rem(6.203480913 + 334.0612426700*ttt,2*pi),
        rem(0.599546497 + 52.9690962641*ttt,2*pi),
        rem(0.874016757 + 21.3299104960*ttt,2*pi),
        rem(5.481293872 + 7.4781598567*ttt,2*pi),
        rem(5.311886287 + 3.8133035638*ttt,2*pi),
        rem(0.02438175*ttt + 0.00000538691*ttt2,2*pi)]

    # Compute terms
    XT0 = 0.0; YT0 = 0.0; ST0 = 0.0
    XT1 = 0.0; YT1 = 0.0; ST1 = 0.0
    XT2 = 0.0; YT2 = 0.0; ST2 = 0.0
    XT3 = 0.0; YT3 = 0.0; ST3 = 0.0
    XT4 = 0.0; YT4 = 0.0; ST4 = 0.0
    for i in 1:15
        if i == 1
            for j in 1:1306
                ap   = transpose(view(Ax,j,4:17))*DelArgs
                XT0 += Ax[j,2]*sin(ap) + Ax[j,3]*cos(ap)
            end
        elseif i == 2
            for j in 1:253
                ap   = transpose(view(Ax,1306 + j,4:17))*DelArgs;
                XT1 += Ax[1306 + j,2]*sin(ap) + Ax[1306 + j,3]*cos(ap);
            end
        elseif i == 3
            for j in 1:36
                ap   = transpose(view(Ax,1559 + j,4:17))*DelArgs;
                XT2 += Ax[1559 + j,2]*sin(ap) + Ax[1559 + j,3]*cos(ap);
            end
        elseif i == 4
            for j = 1:4
                ap   = transpose(view(Ax,1595 + j,4:17))*DelArgs;
                XT3 += Ax[1595 + j,2]*sin(ap) + Ax[1595 + j,3]*cos(ap);
            end
        elseif i == 5
            ap   = transpose(view(Ax,1600,4:17))*DelArgs;
            XT4 += Ax[1600,2]*sin(ap) + Ax[1600,3]*cos(ap);
        elseif i == 6
            for j in 1:962
                ap   = transpose(view(Ay,j,4:17))*DelArgs;
                YT0 += Ay[j,2]*sin(ap) + Ay[j,3]*cos(ap);
            end
        elseif i == 7
            for j in 1:277
                ap   = transpose(view(Ay,962 + j,4:17))*DelArgs;
                YT1 += Ay[962 + j,2]*sin(ap) + Ay[962 + j,3]*cos(ap);
            end
        elseif i == 8
            for j in 1:30
                ap   = transpose(view(Ay,1239 + j,4:17))*DelArgs;
                YT2 += Ay[1239 + j,2]*sin(ap) + Ay[1239 + j,3]*cos(ap);
            end
        elseif i == 9
            for j in 1:5
                ap   = transpose(view(Ay,1269 + j,4:17))*DelArgs;
                YT3 += Ay[1269 + j,2]*sin(ap) + Ay[1269 + j,3]*cos(ap);
            end
        elseif i == 10
            ap   = transpose(view(Ay,1275,4:17))*DelArgs;
            YT4 += Ay[1275,2]*sin(ap) + Ay[1275,3]*cos(ap);
        elseif i == 11
            for j in 1:33
                ap   = transpose(view(As,j,4:17))*DelArgs;
                ST0 += As[j,2]*sin(ap) + As[j,3]*cos(ap);
            end
        elseif i == 12
            for j in 1:3
                ap   = transpose(view(As,33 + j,4:17))*DelArgs;
                ST1 += As[33 + j,2]*sin(ap) + As[33 + j,3]*cos(ap);
            end
        elseif i == 13
            for j in 1:25
                ap   = transpose(view(As,36 + j,4:17))*DelArgs;
                ST2 += As[36 + j,2]*sin(ap) + As[36 + j,3]*cos(ap);
            end
        elseif i == 14
            for j in 1:4
                ap  = transpose(view(As,61 + j,4:17))*DelArgs;
                ST3 += As[61 + j,2]*sin(ap) + As[61 + j,3]*cos(ap);
            end
        else
            ap   = transpose(view(As,66,4:17))*DelArgs;
            ST4 += As[66,2]*sin(ap) + As[66,3]*cos(ap);
        end
    end

    # Get pole corrdinate corrections
    JD = ttt*36525 + 2451545.0
    ΔX = getdx(JD)
    ΔY = getdy(JD)

    # Compute X, Y, and S
    X = (-0.016617 + 2004.191898*ttt - 0.4297829*ttt2 - 0.19861834*ttt3 +
         0.000007578*ttt4 + 0.0000059285*ttt4*ttt)*ARCSEC2RAD +
         (XT0 + XT1*ttt + XT2*ttt2 + XT3*ttt3 + XT4*ttt4)*MUARCSEC2RAD + ΔX;

    Y = (-0.006951 - 0.025896*ttt - 22.4072747*ttt2 + 0.00190059*ttt3 +
         0.001112526*ttt4 + 0.0000001358*ttt4*ttt)*ARCSEC2RAD +
         (YT0 + YT1*ttt + YT2*ttt2 + YT3*ttt3 + YT4*ttt4)*MUARCSEC2RAD + ΔY;

    S = -(X*Y)/2 + (0.000094 + 0.00380865*ttt - 0.00012268*ttt2 - 0.072574*ttt3 +
         0.00002798*ttt4 + 0.00001562*ttt4*ttt)*ARCSEC2RAD +
         (ST0 + ST1*ttt + ST2*ttt2 + ST3*ttt3 + ST4*ttt4)*MUARCSEC2RAD;

    # Compute rotation matrix
    d = atan(sqrt((X^2 + Y^2)/(1 - X^2 - Y^2)))
    a = 1/(1 + cos(d))

    m1 = @SMatrix [1-a*X^2  -a*X*Y  X;
                   -a*X*Y  1-a*Y^2  Y;
                   -X  -Y  1-a*(X^2+Y^2)]
    m2 = @SMatrix [ cos(S) sin(S) 0;
                   -sin(S) cos(S) 0;
                    0 0 1]
    PN = SMatrix{3,3}(m1*m2)

    return PN
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
abstract type AbstractVirtualVolume{T} end

in(pt::AbstractCoordinatePoint{T, 3}, avv::AbstractVirtualVolume{T}) where {T <: SSDFloat} = in(pt, avv.geometry)

function construct_virtual_volume(T, pass::Dict, input_units::NamedTuple, transformations::Transformations)
    construct_virtual_volume(T, pass, input_units, Val{Symbol(pass["model"])}, transformations)
end

"""
    struct DeadVolume{T, G} <: AbstractVirtualVolume{T}

Volume inside which the charge drift is set to zero.

## Parametric types
* `T`: Precision type.
* `G`: Type of `geometry`.

## Fields
* `name::String`: Name of the dead volume, relevant for plotting.
* `geometry::G`: Geometry of the dead volume, see [Constructive Solid Geometry (CSG)](@ref).

## Definition in Configuration File
A `DeadVolume` is defined through an entry in the `virtual_drift_volumes`
array of a detector with `model: dead`. It needs a `geometry` and can optionally be given a `name`.

An example definition of dead volumes looks like this:
```yaml
virtual_drift_volume:
  - name: Volume 1
    model: dead
    geometry: # ...
  - name: Volume 2
    model: dead
    geometry: # ...
```
"""
struct DeadVolume{T, G} <: AbstractVirtualVolume{T}
    name::String
    geometry::G
end

function modulate_driftvector(sv::CartesianVector{T}, pt::CartesianPoint{T}, tl::DeadVolume{T})::CartesianVector{T} where {T <: SSDFloat}
    return CartesianVector{T}(0,0,0)
end

function DeadVolume{T}(dict::Dict, input_units::NamedTuple, transformations::Transformations) where T <: SSDFloat
    n = haskey(dict, "name") ? dict["name"] : "external part"
    g = Geometry(T, dict["geometry"], input_units, transformations)
    return DeadVolume{T, typeof(g)}(n, g)
end

function construct_virtual_volume(T, pass::Dict, input_units::NamedTuple, ::Type{Val{:dead}}, transformations::Transformations)
    DeadVolume{T}(pass, input_units, transformations)
end

struct ArbitraryDriftModificationVolume{T} <: AbstractVirtualVolume{T}
    name::String
    id::Int
    geometry::AbstractGeometry{T}
end

function modulate_driftvector(sv::CartesianVector{T}, pt::CartesianPoint{T}, tl::ArbitraryDriftModificationVolume{T})::CartesianVector{T} where {T <: SSDFloat}
    modulate_driftvector(sv, pt, tl, Val{tl.id})
end
function modulate_driftvector(sv::CartesianVector{T}, pt::CartesianPoint{T}, tl::ArbitraryDriftModificationVolume{T}, ::Type{Val{id}})::CartesianVector{T} where {T <: SSDFloat, id}
    error("""
        This function needs to be overwritten by the user. Use `::Type{Val{<id>}}` for the last argument.
        <id> is the corresponding id specified in the configuration file of the detector.
    """)
end

function ArbitraryDriftModificationVolume{T}(dict::Dict, input_units::NamedTuple, transformations::Transformations) where T <: SSDFloat
    n = haskey(dict, "name") ? dict["name"] : "external part"
    g = Geometry(T, dict["geometry"], input_units, transformations)
    id = Int(dict["id"])
    return ArbitraryDriftModificationVolume{T}(n, id, g)
end

function construct_virtual_volume(T, pass::Dict, input_units::NamedTuple, ::Type{Val{:arbitrary}}, transformations::Transformations)
    ArbitraryDriftModificationVolume{T}(pass, input_units, transformations)
end


text:
using BAI
using Test
using LinearAlgebra

tolerance = 1e-14

@test BAI.reduce_index(10, 3) == 9
@test BAI.reduce_index(10, 11) == 10
@test BAI.recover_index(10, 3) == 11
@test BAI.recover_index(10, 10) == 11
@test BAI.recover_index(10, 11) == 10

# testing reduce_omega
omega = [7, 13, 2, 7, 12]
i, j = 3, 1
omega_tilde = [13, 9, 7, 12]
@test BAI.reduce_omega(omega, i, j) == omega_tilde

# testing reduce_lambda
lambda = [1.2, 3.5, 1.123, 7.32, 10.23, 1.5]
j = 4
lambda_tilde = [1.2, 3.5, 1.123, 10.23, 1.5]
@test BAI.reduce_lambda(lambda, j) == lambda_tilde

# testing recover_lambda
lambda_tilde = [2.3, 7.41, 9.2, 10, 13.5]
i, j = 4, 2
epsilon = 0.01
lambda = [2.3, 9.21, 7.41, 9.2, 10, 13.5]
@test norm(BAI.recover_lambda(lambda_tilde, i, j, epsilon) - lambda) < tolerance

# testing reduce_mu
mu = [1.2, 3.4, 5.6, 7.8]
omega = [0.1, 0.4, 0.3, 0.2]
i, j = 2, 3
epsilon = 0.01
mu_tilde = [1.2, 4.337142857142857, 7.8]
@test norm(BAI.reduce_mu(mu, omega, i, j, epsilon) - mu_tilde) < tolerance

# testing reduce_laplacian
laplacian = [2 -1 -1; -1 1.5 -0.5; -1 -0.5 1.5]
i, j = 2, 1
laplacian_tilde = [1.5 -1.5; -1.5 1.5]
@test BAI.reduce_laplacian(laplacian, i, j) == laplacian_tilde

# testing reduce_laplacian_j
laplacian = [2 -1 -1; -1 1.5 -0.5; -1 -0.5 1.5]
i, j = 2, 1
laplacian_tilde_j = [1; -1]
@test BAI.reduce_laplacian_j(laplacian, i, j) == laplacian_tilde_j

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
function setupcalibrateserine!(p_set, m_set, w_set)
    project_name = "L-Serine-700"
    molecule_names = ["L-Serine"; "DSS"]
    w = [20.0/46; 1.0] # BMRB: DSS is 1 % => 46 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)
end

function setupcalibrateentriesNamJan2022!(p_set, m_set, w_set)

    project_name = "L-Glutamine-700"
    molecule_names = ["L-Glutamine"; "DSS"]
    w = [20.0/0.5; 1.0] # BMRB: DSS is 500 uM => 0.5 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Serine-700"
    molecule_names = ["L-Serine"; "DSS"]
    w = [20.0/46; 1.0] # BMRB: DSS is 1 % => 46 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "D-(+)-Glucose-NRC-600"
    molecule_names = ["D-(+)-Glucose";]
    w = [1.0; ]
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "Nam2022_Serine"
    molecule_names = ["D-(+)-Glucose"; "L-Serine";]
    w = [1.0; 1.0] #
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Leucine-500-2mM"
    molecule_names = ["L-Leucine";]
    w = [1.0]
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Isoleucine-700"
    molecule_names = ["L-Isoleucine"; "DSS"]
    w = [20.0/46; 1.0] # BMRB: DSS is 1 % => 46 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

end

function setupcalibrateentries1!(p_set, m_set, w_set)

    # project_name = "L-Serine-700"
    # molecule_names = ["L-Serine"; "DSS"]
    # w = [20.0/46; 1.0] # BMRB: DSS is 1 % => 46 mM
    # push!(p_set, project_name)
    # push!(m_set, molecule_names)
    # push!(w_set, w)

    project_name = "L-Alanine-700"
    molecule_names = ["L-Alanine"; "DSS"]
    w = [20.0/0.5; 1.0] # BMRB: DSS is 500uM => 0.5 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Threonine-700"
    molecule_names = ["L-Threonine"; "DSS"]
    w = [20.0/0.5; 1.0] # BMRB: DSS is 500uM => 0.5 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Tryptophan-700"
    molecule_names = ["L-Tryptophan"; "DSS"]
    w = [20.0/0.5; 1.0] # BMRB: DSS is 500uM => 0.5 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Valine-700"
    molecule_names = ["L-Valine"; "DSS"]
    w = [20.0/0.5; 1.0] # BMRB: DSS is 500uM => 0.5 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

    project_name = "L-Histidine-700"
    molecule_names = ["L-Histidine"; "DSS"]
    w = [20.0/46; 1.0] # BMRB: DSS is 1 % => 46 mM
    push!(p_set, project_name)
    push!(m_set, molecule_names)
    push!(w_set, w)

end


text:
# JIT compilation of Julia code to PTX

const to = Ref{TimerOutput}()

function timings!(new=TimerOutput())
  global to
  to[] = new
  return
end

timings() = (TimerOutputs.print_timer(to[]); println())

include(joinpath("compiler", "common.jl"))
include(joinpath("compiler", "irgen.jl"))
include(joinpath("compiler", "optim.jl"))
include(joinpath("compiler", "validation.jl"))
include(joinpath("compiler", "rtlib.jl"))
include(joinpath("compiler", "mcgen.jl"))
include(joinpath("compiler", "debug.jl"))
include(joinpath("compiler", "driver.jl"))

function __init_compiler__()
    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("--nvptx-fma-level=1")

    timings!()
end


text:
mutable struct DQNPolicy <: Policy
    q::Tensor # Q network
    s::Tensor # placeholder
    env::AbstractEnvironment
    sess
end

function get_action(policy::DQNPolicy, o::Array{Float64})
    # cannot take a batch of observations
    o = reshape(o, (1, size(o)...))
    TensorFlow.set_def_graph(policy.sess.graph)
    q_val = run(policy.sess, policy.q, Dict(policy.s => o))
    ai = indmax(q_val)
    return actions(policy.env)[ai]
end

"""
    method for getting aciton with corrected policy
"""
function normalize_qvals(qvals::Vector{Float64})
    return qvals / sqrt(sum(qvals .^ 2))
end

function get_action(policy::DQNPolicy, o::Array{Float64}, low_fi_state::Vector{Float64})
    if policy.env.nn_policy == false
        q_low_fi = vec(lowfi_values(policy.env, o))

        if policy.env.correction == true
            # cannot take a batch of observations
            o = reshape(o, (1, size(o)...))
            TensorFlow.set_def_graph(policy.sess.graph)
            q_correction = run(policy.sess, policy.q, Dict(policy.s => o) ) # Array{Float32,2}
            # reshape q_correction
            temp_arr = zeros(n_actions(policy.env)) # Array{Float64,1}
            for i in 1 : n_actions(policy.env)
                temp_arr[i] = q_correction[1, i]
            end
            q_correction = temp_arr # non-normalized
            # q_correction = normalize_qvals(temp_arr) # Array{Float64,1}
            # corrected qvals:
            q_decision = policy.env.correction_weight * q_correction + (1 - policy.env.correction_weight) * q_low_fi
            # q_decision = q_correction + q_low_fi

        else # no correction
            q_decision = q_low_fi
        end
    else # if policy.env.nn_policy == true
        o = reshape(o, (1, size(o)...))
        TensorFlow.set_def_graph(policy.sess.graph)
        q_decision = run(policy.sess, policy.q, Dict(policy.s => o) ) # Array{Float32,2}
    end
    return actions(policy.env)[indmax(q_decision)]
end

function get_value(policy::DQNPolicy, o::Array{Float64})
    o = reshape(o, (1, size(o)...))
    TensorFlow.set_def_graph(policy.sess.graph)
    q_val = run(policy.sess, policy.q, Dict(policy.s => o))
    return q_val
end

function get_value(sess, q, o)
    TensorFlow.set_def_graph(policy.sess.graph)
    o = reshape(o, (1, size(o)...))
    q_val = run(sess, q, Dict(s => o) )
    return q_val
end

function get_value_batch(sess, q, o)
    TensorFlow.set_def_graph(policy.sess.graph)
    q_val = run(sess, q, Dict(s => o) )
    return q_val
end

function POMDPs.action(policy::DQNPolicy, s::S) where S
    obs = nothing
    if isa(policy.env.problem, POMDP)
        obs = convert_o(Vector{Float64}, s, policy.env.problem)
    else
        obs = convert_s(Vector{Float64}, s, policy.env.problem)
    end
    return get_action(policy, obs)
end

function POMDPs.value(policy::DQNPolicy, s)
    obs = nothing
    if isa(policy.env.problem, POMDP)
        obs = convert_o(Array{Float64}, s, policy.env.problem)
    else
        obs = convert_s(Array{Float64}, s, policy.env.problem)
    end
    return get_value(policy, obs)
end

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

mutable struct LSTMPolicy <: Policy
    q::Tensor # Q network
    state::LSTMStateTuple
    s::Tensor # placeholder
    state_ph::LSTMStateTuple # hidden state placeholder
    state_val::LSTMStateTuple
    env::AbstractEnvironment
    sess
end

function LSTMPolicy(sess, env, arch, dueling)
    obs_dim = obs_dimensions(env)
    n_outs = n_actions(env)
    init_c = zeros(1, arch.lstm_size)
    init_h = zeros(1, arch.lstm_size)
    # bs x trace_length x dim
    s1 = placeholder(Float32, shape=[-1, 1, obs_dim...])
    q1, q1_init, q1_state = build_recurrent_q_network(s1,
                                                   arch.convs,
                                                   arch.fc_in,
                                                   arch.fc_out,
                                                   arch.lstm_size,
                                                   n_actions(env),
                                                   final_activation = identity,
                                                   scope= Q_SCOPE,
                                                   reuse=true,
                                                   dueling = dueling)
    state_val = LSTMStateTuple(init_c, init_h)
    return LSTMPolicy(q1, q1_state, s1, q1_init, state_val, env, sess)
end

function get_action!(policy::LSTMPolicy, o::Array{Float64}, sess) # update hidden state
    # cannot take a batch of observations
    o = reshape(o, (1, 1, size(o)...))
    feed_dict = Dict(policy.s => o, policy.state_ph => policy.state_val)
    q_val, state_val = run(sess, [policy.q, policy.state], feed_dict)
    policy.state_val = state_val
    ai = indmax(q_val)
    return actions(policy.env)[ai]
end

function POMDPs.action(policy::LSTMPolicy, o::Array{Float64})
    return get_action!(policy, o, policy.sess)
end

function get_value(policy::LSTMPolicy, o::Array{Float64}, sess) # update hidden state
    # cannot take a batch of observations
    o = reshape(o, (1, 1, size(o)...))
    feed_dict = Dict(policy.s => o, policy.state_ph => policy.state_val)
    q_val, state_val = run(sess, [policy.q, policy.state], feed_dict)
    return q_val
end

function get_value!(policy::LSTMPolicy, o::Array{Float64}, sess) # update hidden state
    # cannot take a batch of observations
    o = reshape(o, (1, 1, size(o)...))
    feed_dict = Dict(policy.s => o, policy.state_ph => policy.state_val)
    q_val, state_val = run(sess, [policy.q, policy.state], feed_dict)
    policy.state_val = state_val
    return q_val
end

function reset_hidden_state!(policy::LSTMPolicy)# could use zero_state from tf
    hidden_size = get(get_shape(policy.state_ph.c).dims[end])
    init_c = zeros(1, hidden_size)
    init_h = zeros(1, hidden_size)
    policy.state_val = LSTMStateTuple(init_c, init_h)
end

# XXX be careful!! this change the hidden state of the policy
function POMDPs.value(policy::LSTMPolicy, o::Array{Float64})
    return get_value!(policy, o, policy.sess)
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
module SDPA

using CxxWrap

# This 'using' is required to suppress a warning about SDPA not having Libdl in its
# dependencies (Libdl is used by BinaryProvider), e.g.: bicycle1885/CodecZlib.jl#26.
using Libdl

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("SDPA not properly installed. Please run Pkg.build(\"SDPA\")")
end

@wrapmodule(joinpath(dirname(dirname(pathof(SDPA))), "deps", "usr", "lib",
                     "libsdpawrap"))

function __init__()
    @initcxx
end

include("blockdiag.jl")
include("blockmat.jl")
include("options.jl")
include("MOI_wrapper.jl")
include("MPB_wrapper.jl")

end # module


text:
function event_occurred(obj::EventHandler, arg0::SpacecraftState, arg1::EventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, EventDetector, jboolean), arg0, arg1, arg2)
end


text:
# Full Training Test (needs GPU for reasonable training time)using MuZero
using MuZero
import Zygote
using Grids
include("../../examples/tictactoe/state.jl")
include("../../examples/tictactoe/game.jl")
include("./old.jl")
include("./inference_tests.jl")
using Base.Threads: @spawn, @threads
config = make_tictactoe_config()
config.batch_size = 16
config.num_actors = 3
config.weight_decay = 1e-4
config.lr_init=0.002f0
config.save_interval = 10000
config.num_simulations = 100
config.num_unroll_steps = 1
losses, dynamics_stats = [], []

p_tests = getpolicytests(TicTacToe)
p_tests.period = 1000

policy_perfs = [];
function calc_policy_perf(l)
    push!(policy_perfs, policy_perf(l.test))
end

test_cb = LearnerCallback(
    "test_cb",
    calc_policy_perf,
    period=1000
)

learner = Learner(
    TicTacToe;
    config=config,
    dir="output",
    cb=[test_cb],
    tests=[p_tests]
)

#selfplay!(learner; limit=config.batch_size)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=33)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=33)
# @spawn run_selfplay!(config, storage, replay_buffer; limit=34)
# rand(1:length(replay_buffer))

for i = 1:(config.num_actors)
    @spawn selfplay!(learner)
end
learner.replay_buffer |> length

#cb = [(1000, 1, [calc_policy_perf]),]

#config.training_steps = 3
#train_network!(learner, losses=losses, callbacks=cb)

config.training_steps = 2000
@show @time fit!(learner; steps=3)
@spawn fit!(learner)
#@spawn train_network!(learner, losses=losses, callbacks=cb)
#@spawn train_network!(config, storage, replay_buffer, losses=losses, callbacks=cb)

@show 1
learner.stats.losses |> length
@show 2


text:
# Define a struct subtype of `AbstractNodeParams`
mutable struct VoltageSource <: AbstractNodeParams
    voltage::Float64
end

# The type-symbol should uniquely identify the type of node
NodeGraphEditor.type_symbol(::Type{VoltageSource}) = :voltage_source
# The number of inputs (left-hand side connectors)
NodeGraphEditor.inputs(::Val{:voltage_source}) = 1
# The number of outputs (right-hand side connectors)
NodeGraphEditor.outputs(::Val{:voltage_source}) = 1
# Returns a default voltage source
NodeGraphEditor.create(::Val{:voltage_source}) = VoltageSource(1.0)

# Similarly for a resistance
mutable struct Resistance <: AbstractNodeParams
    resistance::Float64
end

NodeGraphEditor.type_symbol(::Type{Resistance}) = :resistance
NodeGraphEditor.inputs(::Val{:resistance}) = 1
NodeGraphEditor.outputs(::Val{:resistance}) = 1
NodeGraphEditor.create(::Val{:resistance}) = Resistance(1.0)


text:
import LightGraphs as LG
import MetaGraphs

# Distribution on directed forests, where the probability mass assigned to trees with the final node
# as the root is TREE_PROBABILITY_MASS (uniform over all directed forests that satisfy this property)
# and then the remainder of probability mass (1.0 - TREE_PROBABILITY_MASS) is spread uniformly over all
# directed forests that do not satisfy the property.

TREE_PROBABILITY_MASS = (1.0 - 1e-10)

struct TreePreferenceDiForest <: Gen.Distribution{LG.SimpleDiGraph} end
const tree_preference_diforest = TreePreferenceDiForest()
(::TreePreferenceDiForest)(n) = Gen.random(tree_preference_diforest, n)

function Gen.random(::TreePreferenceDiForest, n::Int)
    graph_with_edges(n,[])
end

# The number of directed forests on n nodes is (n+1) ^ (n-1) as given by Cayley's formula.
function Gen.logpdf(::TreePreferenceDiForest, diforest::LG.SimpleDiGraph, n::Int)
    n != nv(diforest) && return -Inf

    num_directed_forests = (n+1) ^ (n-1)
    num_trees_with_root_last_node = (n ^ (n-2)) / n

    if isTree(diforest) && (n ∈ S.rootsOfForest(diforest))
        return log(TREE_PROBABILITY_MASS / num_trees_with_root_last_node)
    else
        return log((1.0 - TREE_PROBABILITY_MASS) / (num_directed_forests - num_trees_with_root_last_node))
    end
end

export tree_preference_diforest


text:
function Center(P,h,n)
    s=zeros(Float64,1,n-1)
    Pc=zeros(Float64,1,n-1)
    s2=0.0

     for i in 1:n
        if i != h
            for j in 1:n-1
                s[j]+=P[i][j]
             end
            end
        end
    for i in 1:n-1
        Pc[i]=s[i]/n
    end
    return Pc

end

function maxx(P,f,n)
    arg_f=zeros(Float64,1,n-1)
    arg_h=zeros(Float64,1,n-1)
    h=1
        for i in 2:n
            for j in 1:n-1
                arg_f[j]=P[i][j]
                arg_h[j]=P[h][j]
            end
            if f(arg_f...)>f(arg_h...)
              h=i
            end
        end
    return h
end

function newmax(P,Pr,f,h,n)
    arg_f=zeros(Float64,1,n-1)
    arg_h=zeros(Float64,1,n-1)
    arg_pr=zeros(Float64,1,n-1)
    if h!=1
      h2=1
    else
       h2=2
    end

         for i in 1:n
            for j in 1:n-1
                arg_f[j]=P[i][j]
                arg_h[j]=P[h2][j]
                arg_pr[j]=Pr[j]
            end
            if n!=h
              if f(arg_f...)>f(arg_h...)
               h2=i
              end
            end
          end
            if f(arg_pr...)>f(arg_h...)
             return true
            else
             return false
            end

end

function min(P,f,n)
    arg_f=zeros(Float64,1,n-1)
    arg_h=zeros(Float64,1,n-1)
    h=1
        for i in 2:n
            for j in 1:n-1
                arg_f[j]=P[i][j]
                arg_h[j]=P[h][j]
            end
           if f(arg_f...)<f(arg_h...)
            h=i
            end
        end
    return h
end

function reduction(P,l,n)
    PL=P[l]

        for i in 1:n
            for j in 1:n-1
            P[i][j]=0.5*(P[i][j]+PL[j])
            end
        end
end

function reflection(Ph,Pc,P,n)
    a=1.0
    Pr=zeros(Float64,1,n-1)
    for i in 1:n-1
    Pr[i]=(1.0+a)*Pc[i]-a*Ph[i]
    end

    return Pr
end

function expansion(Pr,Pc,P,n)
    c=2.0

    Pe=zeros(Float64,1,n-1)
    for i in 1:n-1
      Pe[i]=(1.0+c)*Pr[i]-c*Pc[i]
    end

    return Pe
end

function contraction(Ph,Pc,n)
    b=0.5
     Pcon=zeros(Float64,1,n-1)
    for i in 1:n-1
      Pcon[i]=b*Ph[i]+(1.0-b)*Pc[i]
    end
    return Pcon
end

function Check(P2, P, x,f,n)

    arg_p=zeros(Float64,1,n-1)
    arg_x=zeros(Float64,1,n-1)

        for i in 1:n-1
            arg_p[i]=P2[i]
            arg_x[i]=P[x][i]
        end
       if f(arg_p...)<f(arg_x...)
            return true
       else
            return false
       end
end

function count_e(P,n)

    wynik= P[1][1]
    for i in 2:n
        wynik=wynik-P[i][1]
        end
    wynik=wynik*wynik
    return wynik

end

function create_P(n)
  if n==1
       P=[[1.0],[2.0]]
  else
   P = [rand(n) for i in 1:n+1]
   end
        return P
end

function m_Nelder(f,n)
    n=n+1
step=0
t=0.00000001
e=1.0
P=create_P(n)
while t<e
ZW=0
h=maxx(P,f,n)
l=min(P,f,n)
Pcon=0.0
Pc=Center(P,h,n)
Pr=reflection(P[h],Pc,P,n)
if Check(Pr,P,l,f,n)
            #f(Pr[1],Pr[2])<f(P[l][1],P[l][2])
    Pe=expansion(Pr,Pc,P,n)
    if Check(Pe,P,h,f,n)   #f(Pe[1],Pe[2])<f(P[h][1],P[h][2])
            for i in 1:n-1
            P[h][i]=Pe[i]
            end
    else
            for i in 1:n-1
            P[h][i]=Pr[i]
            end
    end
elseif newmax(P,Pr,f,h,n)
          if Check(Pr,P,h,f,n)   #f(Pr[1],Pr[2])<f(P[h][1],P[h][2])
          Ph=Pr
          Pcon=contraction(Ph,Pc,n)
          else
            Pcon=contraction(P[h],Pc,n)
          end

          if Check(Pcon,P,h,f,n)   #f(Pcon[1],Pcon[2])<f(P[h][1],P[h][2])
            for i in 1:n-1
            P[h][i]=Pcon[i]
            end
          else
            reduction(P,l,n)
            ZW=1
          end
else
     for i in 1:n-1
       P[h][i]=Pr[i]
      end
end
if ZW==1
end
e=count_e(P,n)
step+=1
    if(step>100)
        e=t^2-1
    end
end

 for i in 1:n-1
 print("Zmienna z")
 print(i)
 print(" = ")
 println(P[1][i])
    end
print("Steps = ")
print(step)
end

max_stars_repo_name: Cooltura/Method-Neldera-Meada

max_issues_repo_name: Cooltura/Method-Neldera-Meada

max_forks_repo_name: Cooltura/Method-Neldera-Meada


text:
#!/bin/julia
using XLSX;
using CSV;
using DataFrames;
using Plots;
using Dates;
using Statistics;
using StatsPlots
import PyPlot;
using LinearAlgebra

function free()
    quality =
        "file/xian_guangdian.csv" |>
        CSV.File |>
        DataFrame |>
        data ->
            begin
                combine(nrow, groupby(select(data, :客户等级), :客户等级)) |>
                data -> rename(data, :nrow => "用户数量") |> println
                combine(nrow, groupby(select(data, [:客户等级, :网络类型]), [:客户等级, :网络类型]))
            end |>
            data ->
                rename(data, :nrow => :quantity, :网络类型 => :net_kind, :客户等级 => :user_level)
    data = combine(groupby(quality, :net_kind), [:user_level, :quantity])
    dict = Dict(
        "5星ABD客户" => "star_5ABD",
        "离线" => "out_link",
        "3星AB客户" => "star_3AB",
        "1星D客户" => "star_1D",
        "1星A客户" => "star_1A",
        "VIP商业个人客户" => "vip",
        "3星AD客户" => "start_3AD",
    )
    1:(data|>nrow) .|>
    i -> begin
        data[i, :net_kind] =
            Dict("农网用户" => "village", "城网用户" => "city", " " => "unknown")[data[
                i,
                :net_kind,
            ]]
        data[i, :user_level] = dict[data[i, :user_level]]
    end
    gp = groupby(data, :net_kind)
    gp |>
    keys .|>
    kind -> @df combine(gp[kind], [:user_level, :quantity]) plot(
        :user_level,
        :quantity,
        label = "$kind",
    ) |> fig -> savefig(fig, "images/first_$kind")
end

function draw_plot()
    file_path = "file/xian_beijing_salary.xlsx"
    salarys = DataFrame(XLSX.readdata(file_path, "Sheet1!C3:D14"), :auto) .|> identity
    salary = [salarys.x1, salarys.x2]

    println("mean:$(salary .|> Statistics.mean)")
    println("var:$(salary .|> Statistics.var)")
    println("std:$(salary .|> Statistics.std)")
    println("cov:$(salary |> Statistics.cov)")
    println("cor:$(salary .|> Statistics.cor)")

    violin(["Xi'an"], salarys.x1, label = "Xi'an")
    violin!(["Beijing"], salarys.x2, label = "Beijing") |>
    fig -> savefig(fig, "images/violin")
    boxplot(["Xi'an"], salarys.x1, label = "Xi'an")
    boxplot!(["Beijing"], salarys.x2, label = "Beijing") |>
    fig -> savefig(fig, "images/box")

end

function map_transform()
    file_path = "file/3movie_metadata.csv"
    movie_metadata = file_path |> CSV.File |> DataFrame
    dropmissing!(movie_metadata, :director_name)
    dict = Dict(
        :color => "Color",
        :num_critic_for_reviews => 0,
        :duration => 0,
        :director_facebook_likes => 0,
        :actor_3_facebook_likes => 0,
        :actor_2_name => "",
        :actor_1_facebook_likes => 0,
        :gross => 0,
        :genres => "",
        :actor_1_name => "",
        :movie_title => "",
        :num_voted_users => 0,
        :cast_total_facebook_likes => 0,
        :actor_3_name => "",
        :facenumber_in_poster => 0,
        :plot_keywords => "",
        :movie_imdb_link => "",
        :num_user_for_reviews => 0,
        :language => "",
        :country => "",
        :content_rating => 0,
        :budget => 0,
        :title_year => 0,
        :actor_2_facebook_likes => 0,
        :imdb_score => 0,
        :aspect_ratio => 0,
        :movie_facebook_likes => 0,
    )
    dict |>
    keys .|>
    key -> select!(
        movie_metadata,
        :,
        key => (col -> col .|> each -> if ismissing(each)
            dict[key]
        else
            each
        end) => key,
    )
end

function join_compine()
    ["file/4ReaderInformation.csv", "file/4ReaderRentRecode.csv"] .|>
    CSV.File .|>
    DataFrame |>
    dates -> begin
        println(dates...)
        innerjoin(dates..., on = :num) |>
        file -> begin
            file |> println
            CSV.write("file/join.csv", file)
        end
    end
end

function self_pca()
    mat = "file/5iris.csv" |> CSV.File |> DataFrame
    gp = groupby(mat, :Species)
    names = [:Sepal_length, :Sepal_width, :Petal_length, :Petal_width, :Species]
    @df gp[1] plot(
        :Sepal_length,
        :Sepal_width,
        :Petal_length,
        zcolor = :Petal_width,
        m = (10, 0.2, :blues, Plots.stroke(0)),
        fontfamily = "Yahei",
        xlabel = "Sepal_length",
        ylabel = "Sepal_width",
        zlabel = "Petal_length",
        title = "山鸢尾",
        label = "山鸢尾",
        w = 0,
    )
    @df gp[2] plot(
        :Sepal_length,
        :Sepal_width,
        :Petal_length,
        zcolor = :Petal_width,
        m = (10, 0.2, :blues, Plots.stroke(0)),
        fontfamily = "Yahei",
        xlabel = "Sepal_length",
        ylabel = "Sepal_width",
        zlabel = "Petal_length",
        title = "变色鸢尾",
        label = "变色鸢尾",
        w = 0,
    )
    @df gp[3] plot(
        :Sepal_length,
        :Sepal_width,
        :Petal_length,
        zcolor = :Petal_width,
        m = (10, 0.2, :blues, Plots.stroke(0)),
        fontfamily = "Yahei",
        xlabel = "Sepal_length",
        ylabel = "Sepal_width",
        zlabel = "Petal_length",
        title = "维吉尼亚鸢尾",
        label = "维吉尼亚鸢尾",
        w = 0,
    )
    data = copy(mat)
    data = select!(data, Not([:Species]))
    transform!(
        data,
        :,
        :Sepal_length => (x -> x .- mean(x)) => :Sepal_length,
        :Sepal_width => (x -> x .- mean(x)) =>:Sepal_width,
        :Petal_length => (x -> x .- mean(x)) =>:Petal_length,
        :Petal_width => (x -> x .- mean(x)) =>:Petal_width,
    )
    data = data |> Matrix
    (values, vectors) = data |> Statistics.cov |> eigen
    # sortperm 会将元素从小到大排序
    p = last(sortperm(values), 2) |> x -> vectors[:, x]
    data * p
end

free();
draw_plot();
map_transform();
join_compine();
self_pca();

max_stars_repo_name: lovebaihezi/lab

max_issues_repo_name: lovebaihezi/lab

max_forks_repo_name: lovebaihezi/lab


text:
#compare parallel Dykstra vs PARSDMM
#look at the total number of PARDMM iterations if
#we use ARADMM to solve Dykstra subproblems.

using SetIntersectionProjection
include(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Dykstra_prox_parallel.jl"))
include(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Dykstra_prox_parallel2.jl"))
using MAT
using LinearAlgebra

ENV["MPLBACKEND"]="qt5agg"
using PyPlot

mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#read velocity model
file = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/compass_velocity.mat"))
m=read(file, "Data")
close(file)

m = m[1:341,200:600];
m = permutedims(m,[2,1])

#PARSDMM options:
options          = PARSDMM_options()
options.FL       = Float32
options.feas_tol = 0.001f0
options.evol_rel_tol = 0.0001f0
set_zero_subnormals(true)

#select working precision
if options.FL==Float64
  TF = Float64
  #TI = Int64
elseif options.FL==Float32
  TF = Float32
  #TI = Int32
end

comp_grid = compgrid((TF(25), TF(6)),(size(m,1), size(m,2)))

m  = convert(Vector{TF},vec(m))
m2 = deepcopy(m)
m3 = deepcopy(m)

#constraints (bounds, vertical monotonicity, total-variation)
constraint = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min     = 1500.0
m_max     = 4000.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

(TV,dummy1,dummy2,dummy3)=get_TD_operator(comp_grid,"TV",options.FL)
m_min     = 0.0
m_max     = 0.1*norm(TV*m,1)
set_type  = "l1"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

BLAS.set_num_threads(2) #2 is fine for a small problem
(P_sub,TD_OPS,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OPS,AtA,l,y) = PARSDMM_precompute_distribute(TD_OPS,set_Prop,comp_grid,options)

println("PARSDMM serial (bounds, bounds on D_z and TV):")
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OPS,set_Prop,P_sub,comp_grid,options);

cg_axis=[0 ; cumsum(log_PARSDMM.cg_it)];
cg_axis=cg_axis[1:10:end]

feas_axis=0:10:(10-1)*size(log_PARSDMM.set_feasibility,1)

#define axis limits and colorbar limits
xmax = comp_grid.d[1]*comp_grid.n[1]
zmax = comp_grid.d[2]*comp_grid.n[2]
vmi=1500
vma=4500
figure();imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection PARSDMM (bounds, bounds on D_z, TV))")

# parallel Dykstra: first set up projectors onto complicated sets
#use same stopping conditions as for PARSDMM
maxit_dyk=20
dyk_feas_tol = deepcopy(options.feas_tol)
obj_dyk_tol  = deepcopy(options.obj_tol)

feas_tol_target= deepcopy(options.feas_tol)
obj_tol_target = deepcopy(options.obj_tol)

#use more accurate ARADMM sub-problem solutions
# we can play with these a bit to see if we can decrease total computational cost
# if set too loose, we will not see overall convergence
options.obj_tol      = 0.8*options.obj_tol
options.feas_tol     = 0.8*options.feas_tol
options.evol_rel_tol = 10*eps(TF)#0.2*options.feas_tol

P=Vector{Any}(undef,3)
P[1]=P_sub[1] #projector onto bounds is easy.

# Set up ARADMM to project onto transform-domain bounds (PARSDMM with 1 constraint set is equivalent to ARADMM)
constraint = Vector{SetIntersectionProjection.set_definitions}()

m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

(P_sub2,TD_OPS2,set_Prop2) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OPS2,AtA2,l,y) = PARSDMM_precompute_distribute(TD_OPS2,set_Prop2,comp_grid,options)
P[2] = inp -> PARSDMM(inp,AtA2,TD_OPS2,set_Prop2,P_sub2,comp_grid,options) #projector onto transform-domain bounds

# Set up ARADMM to project onto anisotropic-TV set (PARSDMM with 1 constraint set is equivalent to ARADMM)
constraint = Vector{SetIntersectionProjection.set_definitions}()

m_min     = 0.0
m_max     = 0.1*norm(TV*m,1)
set_type  = "l1"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

(P_sub3,TD_OPS3,set_Prop3) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OPS3,AtA3,l,y) = PARSDMM_precompute_distribute(TD_OPS3,set_Prop3,comp_grid,options)
P[3] = inp -> PARSDMM(inp,AtA3,TD_OPS3,set_Prop3,P_sub3,comp_grid,options) #projector onto transform-domain bounds

closed_form=Vector{Bool}(undef,3)
closed_form[1]=true
closed_form[2]=false
closed_form[3]=false

figure();imshow(permutedims(reshape(m2,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")

@time (x,obj,feasibility_err_dyk,cg_it,ARADMM_it,l1_P,bounds_P,evol_x_PDyk)=Dykstra_prox_parallel(m2,P,P_sub,TD_OPS,closed_form,maxit_dyk,dyk_feas_tol,obj_dyk_tol)

figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection Dykstra (bounds, bounds on D_z, TV)")

fig, ax = subplots()
ax[:semilogy](feas_axis,log_PARSDMM.set_feasibility[:,1],color="b",label="PARSDMM - bounds",linewidth=5);
ax[:semilogy](feas_axis,log_PARSDMM.set_feasibility[:,2],color="r",label="PARSDMM - monotonicity",linewidth=5);
ax[:semilogy](feas_axis,log_PARSDMM.set_feasibility[:,3],color="k",label="PARSDMM - TV",linewidth=5);
ax[:semilogy]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,1],linestyle=":",color="b",label="Parallel Dykstra - bounds",linewidth=5);
ax[:semilogy]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,2],linestyle=":",color="r",label="Parallel Dykstra - monotonicity",linewidth=5);
ax[:semilogy]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,3],linestyle=":",color="k",label="Parallel Dykstra - TV",linewidth=5);
# ax[:loglog](feas_axis,log_PARSDMM.set_feasibility[:,3],color="k",label="PARSDMM - TV");
# ax[:loglog](feas_axis,log_PARSDMM.set_feasibility[:,1],color="b",label="PARSDMM - bounds");
# ax[:loglog](feas_axis,log_PARSDMM.set_feasibility[:,2],color="r",label="PARSDMM - monotonicity");
# ax[:loglog]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,1],linestyle="--",color="b",label="Parallel Dykstra - bounds");
# ax[:loglog]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,2],linestyle="--",color="r",label="Parallel Dykstra - monotonicity");
# ax[:loglog]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,3],linestyle="--",color="k",label="Parallel Dykstra - TV");
#ax[:semilogy]([0 ; cumsum(ARADMM_it)],ones(size(feasibility_err_dyk,1)).*feas_tol_target,label="target");
title("relative set feasibility error", fontsize=16)
xlabel(L"number of sequential $\ell_1$ projections", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_feasibility.eps",bbox_inches="tight",dpi=1200)

fig, ax = subplots()
ax[:semilogy](cg_axis,log_PARSDMM.set_feasibility[1:end-1,1],color="b",label="PARSDMM - bounds",linewidth=5);
ax[:semilogy](cg_axis,log_PARSDMM.set_feasibility[1:end-1,2],color="r",label="PARSDMM - monotonicity",linewidth=5);
ax[:semilogy](cg_axis,log_PARSDMM.set_feasibility[1:end-1,3],color="k",label="PARSDMM - TV",linewidth=5);
ax[:semilogy]([0;cumsum(cg_it)],feasibility_err_dyk[:,1],linestyle=":",color="b",label="Parallel Dykstra - bounds",linewidth=5);
ax[:semilogy]([0;cumsum(cg_it)],feasibility_err_dyk[:,2],linestyle=":",color="r",label="Parallel Dykstra - monotonicity",linewidth=5);
ax[:semilogy]([0;cumsum(cg_it)],feasibility_err_dyk[:,3],linestyle=":",color="k",label="Parallel Dykstra - TV",linewidth=5);
# ax[:loglog](cg_axis,log_PARSDMM.set_feasibility[1:end-1,1],color="b",label="PARSDMM - bounds");
# ax[:loglog](cg_axis,log_PARSDMM.set_feasibility[1:end-1,2],color="r",label="PARSDMM - monotonicity");
# ax[:loglog](cg_axis,log_PARSDMM.set_feasibility[1:end-1,3],color="k",label="PARSDMM - TV");
# ax[:loglog]([0;cumsum(cg_it)],feasibility_err_dyk[:,1],linestyle="--",color="b",label="Parallel Dykstra - bounds");
# ax[:loglog]([0;cumsum(cg_it)],feasibility_err_dyk[:,2],linestyle="--",color="r",label="Parallel Dykstra - monotonicity");
# ax[:loglog]([0;cumsum(cg_it)],feasibility_err_dyk[:,3],linestyle="--",color="k",label="Parallel Dykstra - TV");

#ax[:semilogy]([0 ; cumsum(ARADMM_it)],ones(size(feasibility_err_dyk,1)).*feas_tol_target,label="target");
title("relative set feasibility error", fontsize=16)
xlabel("number of sequential CG iterations", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_feasibility_CG.eps",bbox_inches="tight",dpi=1200)

fig, ax = subplots()
#ax[:semilogy](abs.(diff(log_PARSDMM.obj)./log_PARSDMM.obj[1:end-1]),color="b",label="PARSDMM",linewidth=5);
ax[:semilogy](cumsum(ARADMM_it[1:end]),evol_x_PDyk[1:length(ARADMM_it[1:end]),1],color="r",linestyle=":",label="Parallel Dykstra",linewidth=5);
ax[:semilogy](log_PARSDMM.evol_x[3:end],color="b",label="PARSDMM",linewidth=5);
#ax[:semilogy](ones(sum(ARADMM_it[2:end])).*obj_tol_target,label="target");
#title(L"relative change in $|| \mathbf{m} - \mathbf{x} ||$", fontsize=16)
title(L"relative change in $|| x^k - x^{k-1} ||$", fontsize=16)
xlabel(L"number of sequential $\ell_1$ projections", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_x_evol.eps",bbox_inches="tight",dpi=1200)

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

#now with a different set of constraints:
# transform-domain rank and bounds
constraint = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min        = 1450.0
m_max        = 4000.0
set_type     = "bounds"
TD_OP        = "identity"
app_mode     = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#rank constraint on vertical derivative of the image
m_min        = 0
m_max        = 15
set_type     = "rank"
TD_OP        = "D_z"
app_mode     = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

BLAS.set_num_threads(2) #2 is fine for a small problem
(P_sub,TD_OPS,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OPS,AtA,l,y) = PARSDMM_precompute_distribute(TD_OPS,set_Prop,comp_grid,options)

println("PARSDMM serial (bounds and rank constraint on D_z):")
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OPS,set_Prop,P_sub,comp_grid,options);

cg_axis=[0 ; cumsum(log_PARSDMM.cg_it)];
cg_axis=cg_axis[1:10:end]

feas_axis=0:10:(10-1)*size(log_PARSDMM.set_feasibility,1)

#define axis limits and colorbar limits
xmax = comp_grid.d[1]*comp_grid.n[1]
zmax = comp_grid.d[2]*comp_grid.n[2]
vmi=1500
vma=4500
figure();imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")
figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection PARSDMM (bounds, rank)")

# parallel Dykstra: first set up projectors onto complicated sets
#use same stopping conditions as for PARSDMM
maxit_dyk=12
dyk_feas_tol = deepcopy(options.feas_tol)
obj_dyk_tol  = deepcopy(options.obj_tol)

feas_tol_target = deepcopy(options.feas_tol)
obj_tol_target  = deepcopy(options.obj_tol)

#use more accurate ARADMM sub-problem solutions
# we can play with these a bit to see if we can decrease total computational cost
# if set too loose, we will not see overall convergence
options.obj_tol      = 0.5*options.obj_tol
options.feas_tol     = 0.5*options.feas_tol
options.evol_rel_tol = 10*eps(TF)

P=Vector{Any}(undef,2)
P[1]=P_sub[1] #projector onto bounds is easy.

# Set up ARADMM to project onto set of transform-domain rank (PARSDMM with 1 constraint set is equivalent to ARADMM)
constraint = Vector{SetIntersectionProjection.set_definitions}()

m_min        = 0
m_max        = 15
set_type     = "rank"
TD_OP        = "D_z"
app_mode     = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

(P_sub2,TD_OPS2,set_Prop2) = setup_constraints(constraint,comp_grid,options.FL)
(TD_OPS2,AtA2,l,y) = PARSDMM_precompute_distribute(TD_OPS2,set_Prop2,comp_grid,options)
P[2] = inp -> PARSDMM(inp,AtA2,TD_OPS2,set_Prop2,P_sub2,comp_grid,options) #projector onto set of transform-domain rank

closed_form=Vector{Bool}(undef,2)
closed_form[1]=true
closed_form[2]=false

figure();imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("model to project")

@time (x,obj,feasibility_err_dyk,cg_it,ARADMM_it,svd_P,~,evol_x_PDyk)=Dykstra_prox_parallel(m,P,P_sub,TD_OPS,closed_form,maxit_dyk,dyk_feas_tol,obj_dyk_tol)

figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]); title("Projection Dykstra (bounds, bounds on D_z, TV)")

fig, ax = subplots()
ax[:semilogy](feas_axis,log_PARSDMM.set_feasibility[:,1],color="b",label="PARSDMM - bounds",linewidth=5);
ax[:semilogy](feas_axis,log_PARSDMM.set_feasibility[:,2],color="r",label="PARSDMM - TD rank",linewidth=5);
ax[:semilogy]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,1],linestyle=":",color="b",label="Parallel Dykstra - bounds",linewidth=5);
ax[:semilogy]([0 ; cumsum(ARADMM_it)],feasibility_err_dyk[:,2],linestyle=":",color="r",label="Parallel Dykstra - TD rank",linewidth=5);
#ax[:semilogy]([0 ; cumsum(ARADMM_it)],ones(size(feasibility_err_dyk,1)).*feas_tol_target,label="target");
title("relative set feasibility error", fontsize=16)
xlabel("number of sequential SVDs", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_feasibility2.eps",bbox_inches="tight",dpi=1200)

fig, ax = subplots()
ax[:semilogy](cg_axis,log_PARSDMM.set_feasibility[1:end-1,1],color="b",label="PARSDMM - bounds",linewidth=5);
ax[:semilogy](cg_axis,log_PARSDMM.set_feasibility[1:end-1,2],color="r",label="PARSDMM - TD rank",linewidth=5);
ax[:semilogy]([0;cumsum(cg_it)],feasibility_err_dyk[:,1],linestyle=":",color="b",label="Parallel Dykstra - bounds",linewidth=5);
ax[:semilogy]([0;cumsum(cg_it)],feasibility_err_dyk[:,2],linestyle=":",color="r",label="Parallel Dykstra - TD rank",linewidth=5);
#ax[:semilogy]([0 ; cumsum(ARADMM_it)],ones(size(feasibility_err_dyk,1)).*feas_tol_target,label="target");
title("relative set feasibility error", fontsize=16)
xlabel("number of sequential CG iterations", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_feasibility_CG2.eps",bbox_inches="tight",dpi=1200)

fig, ax = subplots()
ax[:semilogy](log_PARSDMM.evol_x[3:end],color="b",label="PARSDMM",linewidth=5);
ax[:semilogy](cumsum(ARADMM_it[1:end]),evol_x_PDyk[1:end-1,1],color="r",label="Parallel Dykstra",linewidth=5);
#ax[:semilogy](ones(sum(ARADMM_it[2:end])).*obj_tol_target,label="target");
title(L"relative change in $|| x^k - x^{k-1} ||$", fontsize=16)
xlabel("number of sequential SVDs", fontsize=16)
ax[:legend](fontsize=10)
ax[:tick_params]("both",labelsize=12)
savefig("Dykstra_vs_PARSDMM_x_evol2.eps",bbox_inches="tight",dpi=1200)

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
function create_tax_report_bots(usr, Y, M, client)
    files = glob(joinpath(Pinguu.base_dir, "users/" * usr["user_id"] * "/*/closed_orders/$Y/$M.json"))
    if length(files) > 0
        repeat = 1
        while repeat > 0
            @info("Repeat: " * string(repeat))
            repeat = 1
            for file in files
                tmp = JSON.parsefile(file)
                config = JSON.parsefile(joinpath(split(file, "closed_orders")[1], "config.json"))
                for (k,v) in tmp
                    if !(k in ["trade_num", "profit"])
                        for (k1,v1) in v["trades"]
                            filename = joinpath(joinpath(split(file, "closed_orders")[1], "tax/$Y/$M"), k1 * ".json")
                            if !isfile(filename)
                                entry = Dict()
                                entry["buy"] = []
                                entry["sell"] = []
                                skip = false
                                for (k2,v2) in v1["orders"]
                                    sleep(0.5)
                                    if !haskey(v2, "trades")
                                        trades = try client.get_my_trades(symbol=v2["symbol"], orderId=v2["orderId"])
                                            catch
                                            sleep(2)
                                            try client.get_my_trades(symbol=v2["symbol"], orderId=v2["orderId"])
                                            catch
                                                sleep(2)
                                                nothing
                                            end
                                        end
                                    else
                                        trades = v2["trades"]
                                    end
                                    if trades == nothing
                                        skip = true
                                    else
                                        if v2["side"] == "BUY"
                                            append!(entry["buy"], trades)
                                        else
                                            append!(entry["sell"], trades)
                                        end
                                    end
                                end
                                if !skip
                                    entry["profit"] = 0
                                    entry["commission"] = 0
                                    for (side, trades) in entry
                                        if side in ["buy", "sell"]
                                            for trade in trades
                                                commission = parse(Float64, trade["commission"])
                                                if commission != 0
                                                    if trade["commissionAsset"] == config["farm_coin"]
                                                        commission *= parse(Float64, trade["price"])
                                                    elseif trade["commissionAsset"] != config["base_coin"]
                                                        translator = string(trade["commissionAsset"], config["base_coin"])
                                                        from = split(string(unix2datetime(trade["time"]/1000)), ".")[1]
                                                        to = split(string(unix2datetime(trade["time"]/1000) + Minute(1)), ".")[1]

                                                        klines = client.get_historical_klines(translator, client.KLINE_INTERVAL_1MINUTE, from, to)
                                                        tt = Pinguu.klines2table(klines, time_warp=0)

                                                        commission *= tt[1].close
                                                    end
                                                end
                                                buy_sign = side == "buy" ? -1 : 1
                                                entry["profit"] += buy_sign * parse(Float64, trade["quoteQty"])
                                                entry["profit"] -= commission
                                                entry["commission"] += commission
                                            end
                                        end
                                        entry["sell_time"] = string(unix2datetime(entry["sell"][end]["time"]/1000))
                                        !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
                                        open(filename, "w") do f
                                            JSON.print(f, entry, 4)
                                        end
                                    end
                                else
                                    repeat += 1
                                end
                            end
                        end
                    end
                end
            end
            repeat -= 1
        end
    end
    files = glob(joinpath(Pinguu.base_dir, "users/" * usr["user_id"] * "/*/tax/$Y/$M/*.json"))

    df = DataFrame(Handelsende = [], Profit = [], Gebuehren = [], Symbol = [], A = [], B = [], Kaeufe = [], Verkaeufe = [])
    push!(df, ["", "", "", "", "", "", "", "", ])
    profit = 0
    for file in files
        entry = JSON.parsefile(file)
        push!(df, [
                entry["sell_time"],
                entry["profit"],
                entry["commission"],
                entry["sell"][end]["symbol"],
                "",
                "",
                entry["buy"],
                entry["sell"],
                ])
        profit += entry["profit"]
    end
    push!(df, ["", "", "", "", "", "", "", "", ])
    push!(df, ["Gesamtprofit:", round(profit, digits=3), "EUR", "", "", "", "", "", ]);

    filename = joinpath(Pinguu.base_dir, "stats/tax/user/" * usr["user_id"] * "/bots/$Y/$Y-$M.csv")
    !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
    if length(files) > 0
        CSV.write(filename, df, delim=";")
    end
end
export create_tax_report_bots

function create_tax_report_grid(usr, Y, M)
    config_files = glob(joinpath(Pinguu.base_dir, "grid_bots/user/" * usr["user_id"] * "/*/config.json"))
    bot_paths = []
    for file in config_files
        tmp = JSON.parsefile(file)
        if tmp["testing"] == "0"
            push!(bot_paths, dirname(file))
        end
    end
    df = DataFrame(Handelsende = [], Profit = [], Gebuehren = [], Symbol = [], A = [], B = [], Kaeufe = [], Verkaeufe = [])
    push!(df, ["", "", "", "", "", "", "", "", ])
    profit = 0
    for bot_path in bot_paths
        files = glob(joinpath(bot_path, "history/$Y/$M.json"))
        tmp_config = JSON.parsefile(joinpath(bot_path, "config.json"))
        for file in files
            tmp = JSON.parsefile(file)
            for (day, trades) in tmp
                for trade in trades
                    sell_time = typeof(trade["order"]["time"]) == Int64 ? unix2datetime(trade["order"]["time"]/1000) : trade["order"]["time"]
                    entry = Dict(
                        "sell_time" => sell_time,
                        "profit" => trade["profit"],
                        "commission" => trade["fee"],
                        "symbol" => trade["order"]["symbol"],
                        "sell" => trade["order"],
                    )
                    if trade["type"] == "normal"
                        entry["buy"] = trade["buy_order"]
                    else
                        entry["buy"] = tmp_config["init_state"]["order"]
                    end

                    push!(df, [
                        entry["sell_time"],
                        entry["profit"],
                        entry["commission"],
                        entry["symbol"],
                        "",
                        "",
                        entry["buy"],
                        entry["sell"],
                    ])
                    profit += trade["profit"]
                end
            end
        end
    end
    push!(df, ["", "", "", "", "", "", "", "", ])
    push!(df, ["Gesamtprofit:", round(profit, digits=3), "EUR", "", "", "", "", "", ]);
    filename = joinpath(Pinguu.base_dir, "stats/tax/user/" * usr["user_id"] * "/grid_bots/$Y/$Y-$M.csv")
    !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
    if profit != 0
        CSV.write(filename, df, delim=";")
    end
end
export create_tax_report_grid

function tax_summary(usr, Y)
    profit = 0
    trade_num = 0
    months = []

    for M in lpad.(1:12, 2, "0")
        #
        ## Bots
        files = glob(joinpath(Pinguu.base_dir, "users/" * usr["user_id"] * "/*/tax/$Y/$M/*.json"))
        tmp_profit = 0
        for file in files
            entry = JSON.parsefile(file)
            tmp_profit += entry["profit"]
            trade_num += 1
        end

        #
        ## Grid Bots
        config_files = glob(joinpath(Pinguu.base_dir, "grid_bots/user/" * usr["user_id"] * "/*/config.json"))
        bot_paths = []
        for file in config_files
            tmp = JSON.parsefile(file)
            if tmp["testing"] == "0"
                push!(bot_paths, dirname(file))
            end
        end

        for bot_path in bot_paths
            files = glob(joinpath(bot_path, "history/$Y/$M.json"))
            for file in files
                tmp = JSON.parsefile(file)
                for (day, trades) in tmp
                    for trade in trades
                        sell_time = typeof(trade["order"]["time"]) == Int64 ? unix2datetime(trade["order"]["time"]/1000) : trade["order"]["time"]
                        trade_num += 1
                        tmp_profit += trade["profit"]
                    end
                end
            end
        end

        profit += tmp_profit
        if tmp_profit != 0
            push!(months, tmp_profit)
        end
    end
    summary = Dict(
        "profit" => profit,
        "trade_num" => trade_num,
        "avg_month" => mean(months)
    )
    filename = joinpath(Pinguu.base_dir, "stats/tax/user/" * usr["user_id"] * "/$Y-summary.json")
    open(filename, "w") do f
        JSON.print(f, summary, 4)
    end
end
export tax_summary

function create_tax_report(usr, Y, M)
    client = usr["client"]
    Pinguu.create_tax_report_grid(usr, Y, M)
    Pinguu.create_tax_report_bots(usr, Y, M, client)
    Pinguu.tax_summary(usr, Y)
end
export create_tax_report


text:
# # Visual Transformations

# The function `visual` can be used to give data-independent visual information about the plot
# (plotting function or attributes).

using AlgebraOfGraphics, CairoMakie, PalmerPenguins, DataFrames
set_aog_theme!()

penguins = dropmissing(DataFrame(PalmerPenguins.load()))

data(penguins) * visual(Violin) *
    mapping(:species, :bill_depth_mm, color=:sex, dodge=:sex) |> draw

#

plt = data(penguins) * visual(Violin)
plt *= mapping(:species, :bill_depth_mm, color=:sex, side=:sex, dodge=:island)
draw(plt, axis=(limits=((0.5, 3.5), nothing),))

#

data(penguins) * visual(BoxPlot, show_notch=true) *
    mapping(:species, :bill_depth_mm, color=:sex, dodge=:sex) |> draw

#

data(penguins) *
    mapping(:bill_length_mm, :bill_depth_mm, col=:sex) *
    visual(QQPlot, qqline=:fit) |> draw


text:
"""
# Plotting

```julia
plot(ts::TS, cols::Vector{Int} = collect(1:TSx.ncol(ts)))
plot(ts::TS, cols::Vector{T}) where {T<:Union{String, Symbol}}
plot(ts::TS, cols::T) where {T<:Union{Int, String, Symbol}}
```

Plots a TS object with the index on the x-axis and selected `cols` on
the y-axis. By default, plot all the columns. Columns can be selected
using Int indexes, String(s), or Symbol(s).

# Example
```jldoctest; setup = :(using TSx, DataFrames, Dates, Plots, Random, Statistics)
julia> using Random;
julia> random(x) = rand(MersenneTwister(123), x);
julia> dates = Date("2022-01-01"):Month(1):Date("2022-01-01")+Month(11);

julia> df = DataFrame(Index = dates,
        val1 = random(12),
        val2 = random(12),
        val3 = random(12));

julia> ts = TS(df)
julia> show(ts)
(12 x 3) TS with Dates.Date Index

 Index       val1        val2        val3
 Date        Float64     Float64     Float64
────────────────────────────────────────────────
 2022-01-01  -0.319954    0.974594   -0.552977
 2022-02-01  -0.0386735  -0.171675    0.779539
 2022-03-01   1.67678    -1.75251     0.820462
 2022-04-01   1.69702    -0.0130037   1.0507
 2022-05-01   0.992128    0.76957    -1.28008
 2022-06-01  -0.315461   -0.543976   -0.117256
 2022-07-01  -1.18952    -1.12867    -0.0829082
 2022-08-01   0.159595    0.450044   -0.231828
 2022-09-01   0.501436    0.265327   -0.948532
 2022-10-01  -2.10516    -1.11489     0.285194
 2022-11-01  -0.781082   -1.20202    -0.639953
 2022-12-01  -0.169184    1.34879     1.33361

julia> using Plots

julia> # plot(ts)

# plot first 6 rows with selected columns
julia> # plot(ts[1:6], [:val1, :val3]);

# plot columns 1 and 2 on a specified window size
julia> # plot(ts, [1, 2], size=(600, 400));
```
"""
@recipe function f(ts::TS, cols::Vector{Int} = collect(1:TSx.ncol(ts)))
    seriestype := :line
    size --> (1200, 1200)
    xlabel --> :Index
    ylabel --> join(TSx.names(ts)[cols], ", ")
    legend := true
    label := permutedims(TSx.names(ts)[cols])
    (TSx.index(ts), Matrix(ts.coredata[!, cols.+1])) # increment to account for Index
end

@recipe function f(ts::TS, cols::Vector{T}) where {T<:Union{String, Symbol}}
    colindices = [DataFrames.columnindex(ts.coredata, i) for i in cols]
    colindices .-= 1            # decrement to account for Index
    (ts, colindices)
end

@recipe function f(ts::TS, cols::T) where {T<:Union{Int, String, Symbol}}
    (ts, [cols])
end


text:
# f(x,y)=x+y

function Fxpy(x::Float64,y::Float64)
    return ((x+y)*ones(1,1),nothing)
end

function Fxpy(x::Float64,y)
    return (x*ones(size(y))+y,nothing)
end

function Fxpy(x,y::Float64)
    return (y*ones(size(y))+x,nothing)
end

function Fxpy(x,y)
    if size(x)==(1,1)
        return (x[1]*ones(size(y))+y,nothing)
    elseif size(y)==(1,1)
        return (x+y[1]*ones(size(x)),nothing)
    else
        return (x.+y,nothing)
    end
end

function Fxpy_inplace(value,auxvalue,x,y)
    if size(x)==(1,1)
        copy!(value,x[1]*ones(size(y))); axpy!(1.0,y,value)
    elseif size(y)==(1,1)
        copy!(value,y[1]*ones(size(x))); axpy!(1.0,x,value)
    else
        copy!(value,x); axpy!(1.0,y,value)
    end
end

function Dxpy(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if size(x)==(1,1) && derivativeIDX==1
        axpy!(1.0,[sum(grad_c)],grad_n)
    elseif size(y)==(1,1) && derivativeIDX==2
        axpy!(1.0,[sum(grad_c)],grad_n)
    else
        axpy!(1.0,grad_c,grad_n)
    end
end

if PROC=="GPU"

    function Fxpy(x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            tmp=CudaArray(Float64,size(y))
            gfill!(tmp,x)
            axpy!(1.0,y,tmp)
        elseif size(y)==(1,1)
            tmp=CudaArray(Float64,size(x))
            gfill!(tmp,y)
            axpy!(1.0,x,tmp)
        else
            tmp=CudaArray(zeros(size(x))); copy!(tmp,x); axpy!(1.0,y,tmp)
        end
        return (tmp,nothing)
    end

    function Fxpy_inplace(value::CudaArray,auxvalue,x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            gfill!(value,x); axpy!(1.0,y,value)
        elseif size(y)==(1,1)
            gfill!(value,y); axpy!(1.0,x,value)
        else
            copy!(value,x); axpy!(1.0,y,value)
        end
    end

    function Dxpy(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)

        if size(x)==(1,1) && derivativeIDX==1
            axpy!(1.0,sum(grad_c),grad_n)
        elseif size(y)==(1,1) && derivativeIDX==2
            axpy!(1.0,sum(grad_c),grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    end
end

Derivative[Fxpy]=Dxpy
Inplace[Fxpy]=Fxpy_inplace

import Base.+
+(A::ADnode,B::ADnode)=ADnode(Fxpy,[A B])
+(A::Real,B::ADnode)=ADnode(Fxpy,[ADconst(A) B])
+(A::ADnode,B::Real)=ADnode(Fxpy,[A ADconst(B)])

if PROC=="GPU"
    function +(A::CudaArray,B::CudaArray)
#        tmp=CudaArray(zeros(size(A))); copy!(tmp,A); axpy!(1.0,B,tmp);
#        return tmp
        return CUBLAS.geam('N','N',1.0,A,1.0,B)
    end

    function +(A::CudaArray,B::CudaArray,out::CudaArray)
#        fill!(out,0.0); axpy!(1.0,A,out);axpy!(1.0,B,out)
        CUBLAS.geam!('N','N',1.0,A,1.0,B,out)
    end

end

export Fxpy, +

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Test
using BrainFlow
using JSON

@testset "BoardShim constructor" begin
    # for most boards the master_id == the regular id
    params = BrainFlow.BrainFlowInputParams()
    board_shim = BrainFlow.BoardShim(BrainFlow.CYTON_BOARD, params)
    @test board_shim.master_board_id == Integer(BrainFlow.CYTON_BOARD)

    # for some special boards you need to provide the master_id in the other_info as string...
    params = BrainFlow.BrainFlowInputParams(other_info = "1")
    board_shim = BrainFlow.BoardShim(BrainFlow.STREAMING_BOARD, params)
    @test board_shim.master_board_id == 1

    # test the error if the other_info cannot be parsed as Int32
    params = BrainFlow.BrainFlowInputParams(other_info = "")
    @test_throws BrainFlow.BrainFlowError BrainFlow.BoardShim(BrainFlow.STREAMING_BOARD, params)
end

@testset "prepare_session error" begin
    # use a board that is not connected
    params = BrainFlow.BrainFlowInputParams()
    board_shim = BrainFlow.BoardShim(BrainFlow.CYTON_BOARD, params)
    @test_throws BrainFlow.BrainFlowError("Error in prepare_session INVALID_ARGUMENTS_ERROR", 13) BrainFlow.prepare_session(board_shim)
end

@testset "generated functions" begin
    @test BrainFlow.get_eeg_channels isa Function
end

@testset "data filtering" begin
    @test Int32(BrainFlow.CONSTANT) == 1
    @test Int32(BrainFlow.LINEAR) == 2

    @test Int32(BrainFlow.MEAN) == 0
    @test Int32(BrainFlow.MEDIAN) == 1
    @test Int32(BrainFlow.EACH) == 2

    @test Int32(BrainFlow.BUTTERWORTH) == 0
    @test Int32(BrainFlow.CHEBYSHEV_TYPE_1) == 1
    @test Int32(BrainFlow.BESSEL) == 2

    @test Int32(BrainFlow.NO_WINDOW) == 0
    @test Int32(BrainFlow.HANNING) == 1
    @test Int32(BrainFlow.HAMMING) == 2
    @test Int32(BrainFlow.BLACKMAN_HARRIS) == 3
end

@testset "model params" begin
    params = BrainFlowModelParams(BrainFlow.RESTFULNESS, BrainFlow.DEFAULT_CLASSIFIER)
    @test params.metric == BrainFlow.RESTFULNESS
    @test params.classifier == BrainFlow.DEFAULT_CLASSIFIER
end


text:
const libexpat = Libdl.find_library("libexpat")

include("libexpat/common.jl")
include("libexpat/output.jl")

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module Processor
    using ...Utils
    using ...Functions.CommonStaticFunctions.Single

    export None, ArgMaxIdx, SingleFunc, Chain

    @Fun(Func, output::Any, input::Any)

    None() = MemoryWrapper(Sig(:None), sett -> Func(in -> in))

    ArgMaxIdx() = MemoryWrapper(Sig(:ArgMaxIdx),
                  function (sett)
                        T = sett.outputType
                        Func(
                            function (in)
                                maxIdx = 1
                                maxVal::T = -Inf
                                for i in eachindex(in)
                                    if in[i] > maxVal
                                        maxVal = in[i]
                                        maxIdx = i
                                    end
                                end
                                return T(maxIdx)
                            end)
                end)

    SingleFunc(singleFunc::Single.MemoryWrapper) = MemoryWrapper(Sig(:SingleFunc),
        function (sett)
            return Func(singleFunc(sett))
        end)

    "Chain together processors. f(g(...h(x))) = Chain(f, g, ...h)."
    Chain(processors...) = MemoryWrapper(Sig(:Chain),
        function (sett)
            procs = [processors[i](sett) for i in length(processors):-1:1]
            return Func(
                function (input)
                    out = input
                        for proc in procs
                            out = proc(out)
                        end
                        out
                    end)
        end)

end


text:
let mode = :ast
    function ast_mode_parse(s)
        quote dump(Meta.parse($s)) end
    end

    register() do
        modes[mode] = initrepl(
            ast_mode_parse;
            prompt_text="ast> ",
            prompt_color = :light_blue,
            start_key="",
            mode_name=string(mode)
        )
    end
end


text:
using Faker
using Test

@testset "Faker" begin
include("test_default.jl");
include("test_es_mx.jl");
end;


text:
@testset "cnn" begin
    nfs = [8,16,32,32]
    cnn_model = get_model(nfs)
    X = randn(Float32,(28,28,1,12))
    Y = cnn_model(X)

    @test size(Y) == (10,12)
end


text:
using PyCall,Statistics,Dates
@pyimport matplotlib.pyplot as plt

function check_cell(i,j,array)
    total=0
    for di in -1:1
        for dj in -1:1
            if (0<i+di<=n && 0<j+dj<=n && (i+di,j+dj)!=(i,j))
                total+=array[i+di,j+dj]
            end
        end
    end
    if array[i,j]==1 && 1<total<4
        live=1
    elseif array[i,j]==0 && total==3
        live=1
    else
        live=0
    end
    return live
end

n=128*100
lim=0.8
t_end=2

array=rand(Float64, (n,n))

#I'm sure this isn't the correct way todo this
Threads.@threads for i in 1:n
    Threads.@threads for j in 1:n
        array[i,j]=round(array[i,j])
    end
end

fig, ax_lst = plt.subplots()
im = ax_lst.imshow(array)
plt.pause(0.0001)

times=[]

for t in 1:t_end
    new_array=Array{Int, 2}(undef, n, n)
    ts=Dates.value(Dates.now())
    Threads.@threads for i in 1:n
        Threads.@threads for j in 1:n
            #tried just multithreading the outer loop and it made no difference
            #(i.e. not being slowed down by the dispatch like in python)
            #deminishing returns are seen e.g. 0.24s for n=128 and threads=16 but only 0.25s for n=128 and threads=4
            new_array[i,j]=check_cell(i,j,array)
        end
    end
    push!(times,Dates.value(Dates.now())-ts)
    global array=new_array
    im.set_data(array)
    plt.pause(.0001)
    println(t)
end
println(mean(times)/1000,",",std(times)/1000)

max_stars_repo_name: jagoosw/Conways-Game

max_issues_repo_name: jagoosw/Conways-Game

max_forks_repo_name: jagoosw/Conways-Game


text:
################################################################################
#
#  Free part of the unit group
#
################################################################################

@doc Markdown.doc"""
    regulator(x::Array{T, 1}, abs_tol::Int) -> arb

Compute the regulator $r$ of the elements in $x$, such that the radius of $r$
is less then `-2^abs_tol`.
"""
function regulator(x::Array{T, 1}, abs_tol::Int) where T
  K = _base_ring(x[1])
  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  p = 32

  conlog = Vector{Vector{arb}}(undef, r)

  while true
    q = 2
    for i in 1:r
      conlog[i] = conjugates_arb_log(x[i], p)
      for j in 1:r
        q = max(q, bits(conlog[i][j]))
      end
    end

    A = zero_matrix(ArbField(q, false), r, r)

    for i in 1:r
      for j in 1:r
        A[i, j] = conlog[i][j]
      end
    end

    z = abs(det(A))

    if isfinite(z) && radiuslttwopower(z, -abs_tol)
      return z
    end

    p = 2*p
  end
end

function lower_regulator_bound(K::AnticNumberField)
  return ArbField(64, false)("0.054")
end