text:
using FrameFun.FrameFunInterface, FrameFun.Platforms, FrameFun.ApproximationProblems,
    Test, LinearAlgebra, BasisFunctions, FrameFun.ParameterPaths, FrameFun.WeightedSumPlatforms,
    FrameFun.ExtensionFramePlatforms

ap1 = approximationproblem(platform(Fourier(10)),10)
ap2 = approximationproblem(platform(Fourier(10)^2),(10,10))
dict1 = Fourier(10)
dict2 = Fourier(10)^2
plat1 = (platform(Fourier(10)),10)
plat2 = (platform(Fourier(10)^2),(10,10))

aps = (ap1,ap2,)
alls = (dict1, dict2, aps...,)

@testset "samplingparameter" begin
    # adaptive approximations has no samplingparameter
    for (a,sol) in zip(alls ,(10,(10,10),10,(10,10)))
        @test samplingparameter(a)==sol
    end
    @test samplingparameter(plat1...) == 10
    @test samplingparameter(plat2...) == (10,10)
end

@testset "interpolation_grid" begin
    @test interpolation_grid(plat1...) == interpolation_grid(dict1)
    @test interpolation_grid(plat2...) == interpolation_grid(dict2)

    @test interpolation_grid(ap1) == interpolation_grid(dict1)
    @test interpolation_grid(ap2) == interpolation_grid(dict2)
end

@testset "oversampling_grid" begin
    @test oversampling_grid(dict1) == FourierGrid(10)
    @test oversampling_grid(dict1, 11) == FourierGrid(11)
    @test oversampling_grid(dict2) == FourierGrid(10)^2
    @test oversampling_grid(dict2, (11,11)) == FourierGrid(11)^2
    @test oversampling_grid(plat1...) == FourierGrid(10)
    @test oversampling_grid(ap1) == FourierGrid(10)
    @test oversampling_grid(ap2) == FourierGrid(10)^2
end

@testset "samplingoperator" begin
    op1 = samplingoperator(dict1)
    op2 = samplingoperator(dict2)
    op3 = samplingoperator(plat1...)
    op4 = samplingoperator(plat2...)
    op5 = samplingoperator(ap1)
    op6 = samplingoperator(ap2)

    op1a, op1b = elements(op1)
    @test grid(op1a) == FourierGrid(10)
    @test op1b ≈ IdentityOperator(dict1)

    op2a, op2b = elements(op3)
    @test grid(op2a) == FourierGrid(10)
    @test op2b ≈ IdentityOperator(dict1)

    op3a, op3b = elements(op5)
    @test grid(op3a) == FourierGrid(10)
    @test op3b ≈ IdentityOperator(dict1)
    f2 = (x,y)->exp(x*y)
    @test tensorproduct(op1,op1)*f2≈op2*f2
    @test tensorproduct(op3,op3)*f2≈op4*f2
    @test tensorproduct(op5,op5)*f2≈op6*f2
end

@testset "sampling_grid" begin
    @test sampling_grid(dict1) == FourierGrid(10)
    @test sampling_grid(dict2) == FourierGrid(10)^2
    @test sampling_grid(plat1...) == FourierGrid(10)
    @test sampling_grid(plat2...) == FourierGrid(10)^2
    @test sampling_grid(ap1) == FourierGrid(10)
    @test sampling_grid(ap2) == FourierGrid(10)^2
end

@testset "platform_grid" begin
    @test platform_grid(dict1;samplingstyle=GridStyle(),grid=FourierGrid(1)) == FourierGrid(1)
    @test platform_grid(dict2;samplingstyle=ProductSamplingStyle(GridStyle(),GridStyle()),grid=FourierGrid(1)) == FourierGrid(1)^2
    @test platform_grid(dict2;samplingstyle=GridStyle(),grid=FourierGrid(1)) == FourierGrid(1)
    @test platform_grid(plat1...,grid=FourierGrid(1)) == FourierGrid(1)
    @test platform_grid(plat2...,grid=FourierGrid(1)) == FourierGrid(1)^2
    @test platform_grid(ap1;samplingstyle=GridStyle(),grid=FourierGrid(1)) == FourierGrid(1)
    @test platform_grid(ap2;samplingstyle=GridStyle(),grid=FourierGrid(1)) == FourierGrid(1)
    @test sampling_grid(dict2;samplingstyle=GridStyle(),grid=FourierGrid(1))==FourierGrid(1)
    @test sampling_grid(dict2;samplingstyle=ProductSamplingStyle(GridStyle(),GridStyle()),grid=FourierGrid(1)) == FourierGrid(1)^2
end

@testset "discretemeasure" begin
    @test discretemeasure(dict1) == discretemeasure(sampling_grid(dict1))
    @test discretemeasure(dict2)== discretemeasure(sampling_grid(dict2))
    @test discretemeasure(plat1...)== discretemeasure(sampling_grid(plat1...))
    @test discretemeasure(plat2...)== discretemeasure(sampling_grid(plat2...))
    @test discretemeasure(ap1)== discretemeasure(sampling_grid(ap1))
    @test discretemeasure(ap2)== discretemeasure(sampling_grid(ap2))
end

@testset "measure" begin
    @test measure(dict1) == measure(dict1)
    @test measure(dict2)== measure(dict2)
    @test measure(plat1...)== measure(dict1)
    @test measure(plat2...)== measure(dict2)
    @test measure(ap1)== measure(dict1)
    @test measure(ap2)== measure(dict2)
end

@testset "azdual_dict" begin
    @test operator(azdual_dict(dict1))≈operator(azdual_dict(ap1))≈operator(azdual_dict(plat1...))
    op = azdual_dict(dict2)
    @test op isa TensorProductDict
    @test ≈(operator.(elements(op))...)
    op = azdual_dict(plat2...)
    @test op isa TensorProductDict
    @test ≈(operator.(elements(op))...)
    op = azdual_dict(ap2)
    @test op isa TensorProductDict
    @test ≈(operator.(elements(op))...)
    @test operator(element(azdual_dict(dict2),1))≈operator(azdual_dict(dict1))

    @test operator(azdual_dict(dict1;samplingstyle=GramStyle()))≈
        operator(element(azdual_dict(ap2;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
    @test operator(azdual_dict(dict1;samplingstyle=GramStyle()))≈
        operator(element(azdual_dict(plat2...;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
    @test operator(azdual_dict(dict1;samplingstyle=GramStyle()))≈
        operator(element(azdual_dict(dict2;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
end

@testset "discretization" begin
    op1 = discretization(ap1)
    op2 = discretization(ap2)
    op3 = discretization(plat1...)
    op4 = discretization(plat2...)
    op5 = discretization(dict1)
    op6 = discretization(dict2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator
    @test op6 isa TensorProductOperator

    @test op1 ≈ op3 ≈ op5
    @test op2 ≈ op4 ≈ op6

    f1 = exp
    f2 = (x,y)->exp(x*y)
    op1, b1 = discretization(f1,ap1)
    op2, b2 = discretization(f2,ap2)
    op3, b3 = discretization(f1,plat1...)
    op4, b4 = discretization(f2,plat2...)
    op5, b5 = discretization(f1,dict1)
    op6, b6 = discretization(f2,dict2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator
    @test op6 isa TensorProductOperator

    @test op1 ≈ op3 ≈ op5
    @test op2 ≈ op4 ≈ op6

    @test b1≈b3≈b5
    @test b2≈b4≈b6
end

@testset "dualdiscretization" begin

    op1 = dualdiscretization(ap1)
    op2 = dualdiscretization(ap2)
    op3 = dualdiscretization(plat1...)
    op4 = dualdiscretization(plat2...)
    op5 = dualdiscretization(dict1)
    op6 = dualdiscretization(dict2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator
    @test op6 isa TensorProductOperator

    @test op1 ≈ op3 ≈ op5
    @test op2 ≈ op4 ≈ op6

end

@testset "solver" begin
    op1 = solver(dict1)
    op2 = solver(dict2)
    op3 = solver(plat1...)
    op4 = solver(plat2...)
    op5 = solver(ap1)
    op6 = solver(ap2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator
    @test op6 isa TensorProductOperator

    @test op1 ≈ op3 ≈ op5
    @test op2 ≈ op4 ≈ op6

    @test element(op2,1) ≈ op1 ≈ element(op2,2)

    for STYLE in (InverseStyle, DirectStyle, DualStyle, IterativeStyle,AZStyle,AZSmoothStyle)

        op1 = solver(dict1;solverstyle=STYLE())
        op2 = solver(dict2;solverstyle=ProductSolverStyle(STYLE(),STYLE()))
        op3 = solver(plat1...;solverstyle=STYLE())
        op4 = solver(plat2...;solverstyle=ProductSolverStyle(STYLE(),STYLE()))
        op5 = solver(ap1;solverstyle=STYLE())
        op6 = solver(ap2;solverstyle=ProductSolverStyle(STYLE(),STYLE()))

        @test op2 isa TensorProductOperator
        @test op4 isa TensorProductOperator
        @test op6 isa TensorProductOperator

        @test op1 ≈ op3 ≈ op5
        @test op2 ≈ op4 ≈ op6

        @test element(op2,1) ≈ op1 ≈ element(op2,2)
    end
end

@testset "AZ_A, AZ_Z, AZ_Zt, plungeoperator, firstAZstepoperator, plungerank" begin
    A1 = AZ_A(dict1)
    A2 = AZ_A(dict2)
    @test tensorproduct(A1,A1) ≈ A2

    Z1 = AZ_Z(dict1)
    Z2 = AZ_Z(dict2)
    @test tensorproduct(Z1,Z1) ≈ Z2

    Zt1 = AZ_Zt(dict1)
    Zt2 = AZ_Zt(dict2)
    @test tensorproduct(Zt1,Zt1) ≈ Zt2

    @test Zt2 ≈ Z2'

    P1 = plungeoperator(dict1)
    P2 = plungeoperator(dict2)
    @test 1+norm(tensorproduct(P1,P1) -P2)≈1
    @test I-A2*Zt2 ≈ P2

    M1 = firstAZstepoperator(dict1)
    M2 = firstAZstepoperator(dict2)
    @test 1+norm(tensorproduct(M1,M1) - M2)≈1
    @test 1+norm(A2*P2 - M2)≈1

    @test plungerank(dict1) <= 1
    @test plungerank(dict2) <= 1
end

@testset begin "ParameterPath"
    p = IncrementalCartesianParameterPath{2}()
    P = ExtensionFramePlatform(WeightedSumPlatform(platform(ChebyshevT(10)^2), (x,y)->1.,
                (x,y)->sqrt(x^2+y^2)), .9 .* UnitDisk())
    path = HierarchyPath(p,ProductPath(p,p))
    paramP = parametrizedplatform(P, path)

    paramPparam=45
    Pparam = path[paramPparam]

    @test SamplingStyle(P)==SamplingStyle(paramP)
    @test SolverStyle(P, SamplingStyle(P)) == SolverStyle(paramP, SamplingStyle(P))
    @test sampling_grid(P,Pparam)≈sampling_grid(paramP, paramPparam)
    @test dimensions(dictionary(P,Pparam)) == dimensions(dictionary(paramP,paramPparam))
    @test dimensions(azdual_dict(P,Pparam))==dimensions(azdual_dict(paramP, paramPparam))
    @test AZ_A(P,Pparam)≈AZ_A(paramP, paramPparam)
    @test AZ_Zt(P,Pparam)≈AZ_Zt(paramP, paramPparam)
    @test AZ_A(P,Pparam)≈AZ_A(paramP, paramPparam)
    @test typeof(solver(P,Pparam))==typeof(solver(paramP,paramPparam))
end


text:
#Settings

length_min = [50, 50, 50, 50, 50, 50]
length_max = [50, 50, 50, 50, 50, 50]
categories = [:guessed_difficulty]

ext: jl


text:
my_tests = ["utils.jl",
            "transformations.jl",
            "isomap.jl",
            "lle.jl",
            "hlle.jl",
            "ltsa.jl",
            "lem.jl",
            "diffmaps.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
# plot lidar's sensing data grouped by hour

module GroupByHour
    using Plots, DataFrames, CSV, Statistics
    pyplot()

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

        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # sort
        df_sort = sort(df_new, "hour")

        # grouping
        df_grp = groupby(df_sort, "hour")

        # mean by group
        grps = Array(0:23)
        ovsbs = [mean(grp.lidar) for grp in df_grp]
        plot(grps, ovsbs, label="observation")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/group_by_hour/group_by_hour.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
Base.show(io::IO, sp::DiscreteUnipartiteSpeciesPool) = begin
    stack = CrayonStack()
    Base.print(io, stack, "discrete species pool with ")
    Base.print(io, push!(stack, crayon"green"), numspecies(sp))
    Base.print(io, pop!(stack), " species and ")
    Base.print(io, push!(stack, crayon"blue"), 1)
    Base.print(io, pop!(stack), " partition")
end

Base.show(io::IO, sp::DiscreteKpartiteSpeciesPool) = begin
    stack = CrayonStack()

    Base.print(io, stack, "discrete species pool with ")
    Base.print(io, push!(stack, crayon"green"), numspecies(sp))
    Base.print(io, pop!(stack), " species across ")
    Base.print(io, push!(stack, crayon"blue"), sp.partitions)
    Base.print(io, pop!(stack), " partitions")
end


text:
doppler_comp_simple = calc_doppler_component_simple(lambda,obs);
genpca_simple_out = fit_gen_pca_rv(obs,doppler_comp_simple)
rvs_out = est_rvs_from_pc(obs,genpca_simple_out[1],genpca_simple_out[2][:,1])

true_rvs = [-10.0*cos(2pi*phases[i]) for i in 1:25]  # not actually sure what true rvs are

@test maximum(abs(rvs_out .- true_rvs)) < 0.1

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
time_steps = 1:24

base_dir = string(dirname(dirname(pathof(PowerSimulations))))
DATA_DIR = joinpath(base_dir, "test/test_data")
include(joinpath(DATA_DIR, "data_5bus_pu.jl"))
include(joinpath(DATA_DIR, "data_14bus_pu.jl"))

#Base Systems
nodes = nodes5()
c_sys5 = System(nodes, thermal_generators5(nodes), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5))
        add_forecast!(c_sys5, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
end

nodes = nodes14()
c_sys14 = System(nodes, thermal_generators14(nodes), loads14(nodes), branches14(nodes), nothing, 100.0, nothing, nothing)
for (ix, l) in enumerate(get_components(PowerLoad, c_sys14))
    add_forecast!(c_sys14, l, Deterministic("get_maxactivepower", timeseries_DA14[ix]))
end

PTDF5 = PTDF(c_sys5);
PTDF14 = PTDF(c_sys14);

#System with Renewable Energy
nodes = nodes5()
c_sys5_re = System(nodes, vcat(thermal_generators5(nodes), renewable_generators5(nodes)), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_re))
        add_forecast!(c_sys5_re, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
    for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_re))
        add_forecast!(c_sys5_re, r, Deterministic("get_rating", ren_timeseries_DA[t][ix]))
    end
end

nodes = nodes5()
c_sys5_re_only = System(nodes, renewable_generators5(nodes), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_re_only))
        add_forecast!(c_sys5_re_only, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
    for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_re_only))
        add_forecast!(c_sys5_re_only, r, Deterministic("get_rating", ren_timeseries_DA[t][ix]))
    end
end

nodes = nodes5()
# System with HydroPower Energy
c_sys5_hy = System(nodes, vcat(thermal_generators5(nodes), hydro_generators5(nodes)[1]), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
   for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hy))
       add_forecast!(c_sys5_hy, l, Deterministic("maxactivepower", load_timeseries_DA[t][ix]))
   end
   for (ix, h) in enumerate(get_components(HydroGen, c_sys5_hy))
       add_forecast!(c_sys5_hy, h, Deterministic("rating", hydro_timeseries_DA[t][ix]))
   end
end

nodes = nodes5()
c_sys5_hyd = System(nodes, vcat(thermal_generators5(nodes), hydro_generators5(nodes)[2]), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
   for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hyd))
       add_forecast!(c_sys5_hyd, l, Deterministic("maxactivepower", load_timeseries_DA[t][ix]))
   end
   for (ix, h) in enumerate(get_components(HydroGen, c_sys5_hyd))
       add_forecast!(c_sys5_hyd, h, Deterministic("rating", hydro_timeseries_DA[t][ix]))
   end
end

#System with Storage Device
nodes = nodes5()
c_sys5_bat = System(nodes, thermal_generators5(nodes), loads5(nodes), branches5(nodes), battery5(nodes), 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_bat))
        add_forecast!(c_sys5_bat, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
end

#System with Interruptible Load
nodes = nodes5()
c_sys5_il = System(nodes, thermal_generators5(nodes), vcat(loads5(nodes), interruptible(nodes)), branches5(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_il))
        add_forecast!(c_sys5_il, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
    for (ix, i) in enumerate(get_components(InterruptibleLoad, c_sys5_il))
        add_forecast!(c_sys5_il, i, Deterministic("get_maxactivepower", Iload_timeseries_DA[t][ix]))
    end
end

#Systems with HVDC data in the branches
nodes = nodes5()
c_sys5_dc = System(nodes, vcat(thermal_generators5(nodes), renewable_generators5(nodes)), loads5(nodes), branches5_dc(nodes), nothing, 100.0, nothing, nothing)
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_dc))
        add_forecast!(c_sys5_dc, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
    for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_dc))
        add_forecast!(c_sys5_dc, r, Deterministic("get_rating", ren_timeseries_DA[t][ix]))
    end
end

nodes = nodes14()
c_sys14_dc = System(nodes, thermal_generators14(nodes), loads14(nodes), branches14_dc(nodes), nothing, 100.0, nothing, nothing)
for (ix, l) in enumerate(get_components(PowerLoad, c_sys14_dc))
    add_forecast!(c_sys14_dc, l, Deterministic("get_maxactivepower", timeseries_DA14[ix]))
end

PTDF5_dc = PTDF(c_sys5_dc);
PTDF14_dc = PTDF(c_sys14_dc);

# System to test UC Forms
#Park City and Sundance Have non-binding Ramp Limitst at an Hourly Resolution
# Solitude, Sundance and Brighton have binding time_up constraints.
# Solitude and Brighton have binding time_dn constraints.
# Sundance has non-binding Time Down constraint at an Hourly Resolution
# Alta, Park City and Brighton start at 0.
thermal_generators5_uc_testing(nodes5) = [ThermalStandard("Alta", true, nodes5[1], 0.0, 0.0,
           TechThermal(0.5, PowerSystems.ST, PowerSystems.COAL, (min=0.2, max=0.40),  (min = -0.30, max = 0.30), nothing, nothing),
           ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
           ),
           ThermalStandard("Park City", true, nodes5[1], 0.0, 0.0,
               TechThermal(2.2125, PowerSystems.ST, PowerSystems.COAL, (min=0.65, max=1.70), (min =-1.275, max=1.275), (up=0.02, down=0.02), nothing),
               ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
           ),
           ThermalStandard("Solitude", true, nodes5[3], 2.7, 0.00,
               TechThermal(5.20, PowerSystems.ST, PowerSystems.COAL, (min=1.0, max=5.20), (min =-3.90, max=3.90), (up=0.0012, down=0.0012), (up=5.0, down=3.0)),
               ThreePartCost((0.0, 3000.0), 0.0, 3.0, 1.5)
           ),
           ThermalStandard("Sundance", true, nodes5[4], 0.0, 0.00,
               TechThermal(2.5, PowerSystems.ST, PowerSystems.COAL, (min=1.0, max=2.0), (min =-1.5, max=1.5), (up=0.015, down=0.015), (up=2.0, down=1.0)),
               ThreePartCost((0.0, 4000.0), 0.0, 4.0, 2.0)
           ),
           ThermalStandard("Brighton", true, nodes5[5], 6.0, 0.0,
               TechThermal(7.5, PowerSystems.ST, PowerSystems.COAL, (min=3.0, max=6.0), (min =-4.50, max=4.50), (up=0.0015, down=0.0015), (up=5.0, down=3.0)),
               ThreePartCost((0.0, 1000.0), 0.0, 1.5, 0.75)
           )];
nodes = nodes5()
c_sys5_uc = System(nodes, thermal_generators5_uc_testing(nodes), loads5(nodes), branches5(nodes), nothing, 100.0, nothing, nothing);
for t in 1:2
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_uc))
        add_forecast!(c_sys5_uc, l, Deterministic("get_maxactivepower", load_timeseries_DA[t][ix]))
    end
    for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_uc))
        add_forecast!(c_sys5_uc, r, Deterministic("get_rating", ren_timeseries_DA[t][ix]))
    end
    for (ix, i) in enumerate(get_components(InterruptibleLoad, c_sys5_uc))
        add_forecast!(c_sys5_uc, i, Deterministic("get_maxactivepower", Iload_timeseries_DA[t][ix]))
    end
end

c_sys5_ed = System(nodes, vcat(thermal_generators5_uc_testing(nodes), renewable_generators5(nodes)), vcat(loads5(nodes), interruptible(nodes)), branches5(nodes), nothing, 100.0, nothing, nothing);

for t in 1:2 # loop over days
    for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_ed))
        ta = load_timeseries_DA[t][ix]
        for i in 1:length(ta) # loop over hours
            ini_time = timestamp(ta[i]) #get the hour
            data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
            add_forecast!(c_sys5_ed, l, Deterministic("get_maxactivepower", data))
        end
    end
end
for t in 1:2
    for (ix, l) in enumerate(get_components(RenewableGen, c_sys5_ed))
        ta = load_timeseries_DA[t][ix]
        for i in 1:length(ta) # loop over hours
            ini_time = timestamp(ta[i]) #get the hour
            data = when(load_timeseries_RT[t][ix], hour,hour(ini_time[1])) # get the subset ts for that hour
            add_forecast!(c_sys5_ed, l, Deterministic("get_rating", data))
        end
    end
end
for t in 1:2
    for (ix, l) in enumerate(get_components(InterruptibleLoad, c_sys5_ed))
        ta = load_timeseries_DA[t][ix]
        for i in 1:length(ta) # loop over hours
            ini_time = timestamp(ta[i]) #get the hour
            data = when(load_timeseries_RT[t][ix], hour,hour(ini_time[1])) # get the subset ts for that hour
            add_forecast!(c_sys5_ed, l, Deterministic("get_maxactivepower", data))
        end
    end
end

#reserve_uc = reserve5(get_components(ThermalStandard, c_sys5_uc))
#add_component!(c_sys5_uc, reserve_uc[1])
#add_forecast!(c_sys5_uc, reserve_uc[1], Deterministic("get_requirement", Reserve_ts[1][1]))
#add_forecast!(c_sys5_uc, reserve_uc[1], Deterministic("get_requirement", Reserve_ts[2][1]))

function build_init(gens, data)
    init = Vector{InitialCondition}(undef, length(collect(gens)))
    for (ix,g) in enumerate(gens)
        init[ix] = InitialCondition(g,
                    PSI.UpdateRef{Device}(Symbol("P_$(typeof(g))")),
                    data[ix],TimeStatusChange)
    end
    return init
end


text:
using FastRecurFlux
using Flux
using TrackerFlux
using Random
using Statistics
using Test

@testset "FastRecurFlux.jl+TrackerFlux.jl" begin
    Random.seed!(1234)
    x = randn(Float32, 10, 1, 5)
    y = mean(x, dims = 1)
    for layer in (GRU, LSTM)
        model = Chain(layer(10, 100), layer(100, 1)) |> TrackerFlux.track
        function loss(x, y)
            xs = Flux.unstack(x, 3)
            ys = Flux.unstack(y, 3)
            l = 0f0
            for t in 1:length(xs)
                ŷ = model(xs[t])
                l += Flux.mse(ys[t], ŷ)
            end
            return l / length(xs)
        end
        ps = Flux.params(model)
        data = repeat([(x, y)], 100)
        opt = ADAM(1e-3, (0.9, 0.999))
        cb = () -> Flux.reset!(model)
        Flux.@epochs 10 Flux.train!(loss, ps, data, opt, cb = cb)
        @test loss(x, y) < 0.01
    end
end


text:
module Fuzz

using ..DataFlow
import DataFlow: thread!

export grow

grow(::Type{T}, value) where T <: Vertex= T(value())

function grow(::Type{T}, nodes::Integer, edges::Integer = nodes) where T <: Vertex
  vs = [grow(T, ()->i) for i = 1:nodes]
  for _ = 1:edges
    thread!(rand(vs), rand(vs))
  end
  return rand(vs)
end

# using Atom
# function testcase(test, tries, cap = 10)
#   @progress for nodes = 1:cap
#     for _ = 1:tries
#       g = grow(DVertex, nodes)
#       test(g) || return g
#     end
#   end
# end

end


text:
using Distributed
addprocs(Sys.CPU_THREADS)

import PyPlot; const plt = PyPlot
import DSP: conv
import Random: seed!, randperm
import Statistics: std
import LinearAlgebra: norm, dot
import Arpack: eigs
using Printf

import KirchMig

dz, dx = 15, 15
x = 0:dx:1000; nx = length(x)
z = 0:dz:1000; nz = length(z)

rho = 1000ones(nz, nx)
rho[div(330,dx):end,:] .+= 1000
rho[div(670,dz):end,:] .-= 500

vel = 2000 .+ zero(rho)

imp = 1e-3rho.*vel;

refl = [zeros(nx)'
    (imp[2:end,:] .- imp[1:end-1,:])./(imp[2:end,:] .+ imp[1:end-1,:])]

refl[div(nz, 2)-1, div(nx, 2)] += 0.2
refl[div(nz, 2)+1, div(nx, 2)] += 0.2;

nr = 46

rec_z = zeros(nr)
rec_x = range(x[1], stop=x[end], length=nr)

ns = 10
seed!(12)
src_z = zeros(ns)
src_x = sort(rec_x[randperm(nr)][1:ns])

trav_r = KirchMig.eikonal_const_vel([rec_z rec_x], z, x, vel[1]);
trav_s = KirchMig.eikonal_const_vel([src_z src_x], z, x, vel[1]);

t = 0:0.008:1; nt = length(t)
G = KirchMig.KirchMap(t, trav_r, trav_s)

ricker(t0, f) = @. (1 - 2pi^2 * f^2 * t0^2) * exp(-pi^2 * f^2 * t0^2)
rick_dtt = ricker(t .- t[div(nt,5)], 15);
@views rick_dtt[2:end-1] = (rick_dtt[1:end-2] - 2.0*rick_dtt[2:end-1] + rick_dtt[3:end])/(t[2] - t[1])^2;
rick_dtt /= maximum(abs.(rick_dtt));

W = KirchMig.ConvMap(rick_dtt, nr, ns, nt);

L = W*G;

@time d = L*refl[:];

seed!(1)
n = randn(size(d))
d += n*std(d)/std(n);

@time m_mig = L'd;

xran = x[end]-x[1]
mod_mig = reshape(m_mig, nz, nx)
vmin, vmax = maximum(abs.(mod_mig))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_mig, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_mig[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Migrated image")
fig.tight_layout()

J0(m) = norm(L*m - d)^2

function ∇J0!(storage, m)
    storage[1:end] = 2L'L*m - 2m_mig
end

seed!(123);

nL = size(L, 2)
u = randn(nL)
v = randn(nL)
h = 0.001*maximum(abs.(u))

storage = zero(u)
∇J0!(storage, u)

g1 = (J0(u + h*v) - J0(u - h*v))/2h
g2 = dot(storage, v)
err = 100abs(g1 - g2)/((g1+g2)/2)
println("$(@sprintf("%.2f", err))% error")

import Optim

@time res = Optim.optimize(J0, ∇J0!, zeros(size(L, 2)),
                           Optim.ConjugateGradient(),
                           Optim.Options(iterations = 20, show_trace=true))

pclip=0.5
xran = x[end]-x[1]
mod_lsm_0 = reshape(res.minimizer, nz, nx)
vmin, vmax = pclip*maximum(abs.(mod_lsm_0))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_lsm_0, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_lsm_0[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="LSM image: no reg.")
fig.tight_layout()

@time λ, ϕ = eigs(L'L; nev=1, maxiter=2); λ = λ[1]

function TV(x)
    grad = KirchMig.gradient(x, nz, nx) # Returns a nz × nx × 2 array in 2D
    norm_grad = sqrt.(sum(grad.^2, dims=3))
    e = maximum([100eps(), 0.01maximum(norm_grad)])
    return 0.002*λ*sum(sqrt.(norm_grad.^2 .+ e^2))
end
function ∇TV!(storage, x)
    grad = KirchMig.gradient(x, nz, nx)
    norm_grad = sqrt.(sum(grad.^2, dims=3))
    e = maximum([100eps(), 0.01maximum(norm_grad)])
    grad ./= sqrt.(norm_grad.^2 .+ e^2)

    storage[1:end] = -0.002*λ*KirchMig.divergence(grad[:], nz, nx)[:]
end

seed!(123);

nL = size(L, 2)
u = randn(nL)
v = randn(nL)
h = 0.001*maximum(abs.(u))

storage = zero(u)
∇TV!(storage, u)

g1 = (TV(u + h*v) - TV(u - h*v))/2h
g2 = dot(storage, v)
err = 100abs(g1 - g2)/((g1+g2)/2)
println("$(@sprintf("%.2f", err))% error")

JTV(x) = J0(x) + TV(x)

function ∇JTV!(storage, x)
    ∇J0!(storage, x)
    s = storage[:]
    ∇TV!(storage, x)
    storage[1:end] = s + storage[1:end]
end

# Check gradient
seed!(123);

nL = size(L, 2)
u = randn(nL)
v = randn(nL)
h = 0.001*maximum(abs.(u))

storage = zero(u)
∇JTV!(storage, u)

g1 = (JTV(u + h*v) - JTV(u - h*v))/2h
g2 = dot(storage, v)
err = 100abs(g1 - g2)/((g1+g2)/2)
println("$(@sprintf("%.2f", err))% error")

@time res_tv = Optim.optimize(JTV, ∇JTV!, zeros(size(L, 2)),
                              Optim.ConjugateGradient(),
                              Optim.Options(iterations = 20, show_trace=true))

pclip=0.5
xran = x[end]-x[1]
mod_lsm_tv = reshape(res_tv.minimizer, nz, nx)
vmin, vmax = pclip*maximum(abs.(mod_lsm_tv))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_lsm_tv, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_lsm_tv[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="LSM image: TV reg.")
fig.tight_layout()


text:
# Adapted from the maple scripts in libxc
# Look at scripts/maple2c_new.pl in particular
#         maple/util.jl

using SymPy

r0, r1 = symbols("ρ1 ρ2")
s0, s1, s2 = symbols("σ0 σ1 σ2")

# unpol, 3D
DIMENSIONS = 3
RS_FACTOR    = (3/(4*PI))^(1//3)
X_FACTOR_C   = 3//8*(3/PI)^(1//3)*4^(2//3)
r_ws(n) = RS_FACTOR/n^(1//DIMENSIONS)
dens(r0, r1) = r0
zeta(r0, r1) = 0
xs0(r0, r1, sigma0, sigma2) = sqrt(sigma0/4)/((r0/2)^(1 + 1//DIMENSIONS))
xs1(r0, r1, sigma0, sigma2) = sqrt(sigma0/4)/((r0/2)^(1 + 1//DIMENSIONS))
xt(r0, r1, sigma0, sigma1, sigma2) = sqrt(sigma0)/r0^(1 + 1//DIMENSIONS)

f_zeta(z) = ((1 + z)^(4//3) + (1 - z)^(4//3) - 2)/(2^(4//3) - 2)

# Select case to print code for
case = "lda_x"
case = "lda_c_vwn"

include(case * ".jl")

println("#")
println("#-- " * case)
println("#")

# Energy per unit particle
f_expr = f(r_ws(dens(r0, r1)), zeta(r0, r1), xt(r0, r1, s0, s1, s2), xs0(r0, r1, s0, s2), xs1(r0, r1, s0, s2))

# Energy per unit volume
E_expr = f_expr * dens(r0, r1)

# Derivative with respect to r0
dE_expr = diff(E_expr, r0)

symbol = ["ε", "Vρ"]
f_expr = simplify(f_expr)
E_expr = simplify(E_expr)

sexpr, simplified = sympy.cse([f_expr, dE_expr])

for (var, expr) in sexpr
    println(sympy.julia_code(expr, assign_to=var))
end
println()
for (i, simp) in enumerate(simplified[1])
    println(symbol[i], " = ", sympy.julia_code(simp))
end


text:
module Pads

"""
|  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |
|  |  |  |  T  |  |  |  |
|  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |
"""

end # Pads


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export EssencialNeighborhood

struct EssencialNeighborhood <: NeighborhoodAlgorithm
    k::Int
end

function EssencialNeighborhood()
    return EssencialNeighborhood(32)
end

function optimize_neighborhood!(algo::EssencialNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where {T}
end

function neighborhood(algo::EssencialNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where {T}
    nbuffer::Vector{Int32} = Vector{Int}(undef, 0)
    knn = search(index, dist::Function, item, KnnResult(algo.k))
    visible = Set{Int32}()

    @inbounds for p in knn
        in(p.objID, visible) && continue
        for neighbor in index.links[p.objID]
            push!(visible, neighbor)
        end

        push!(nbuffer, p.objID)
    end

    return knn, nbuffer
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    shift_constant(set::MOI.AbstractScalarSet, offset)

Returns a new scalar set `new_set` such that `func`-in-`set` is equivalent to
`func + offset`-in-`new_set`.

Only define this function if it makes sense to!

Use [`supports_shift_constant`](@ref) to check if the set supports shifting:
```Julia
if supports_shift_constant(typeof(old_set))
    new_set = shift_constant(old_set, offset)
    f.constant = 0
    add_constraint(model, f, new_set)
else
    add_constraint(model, f, old_set)
end
```

See also [`supports_shift_constant`](@ref).

## Examples

The call `shift_constant(MOI.Interval(-2, 3), 1)` is equal to
`MOI.Interval(-1, 4)`.
"""
function shift_constant end

"""
    supports_shift_constant(::Type{S}) where {S<:MOI.AbstractSet}

Return `true` if [`shift_constant`](@ref) is defined for set `S`.

See also [`shift_constant`](@ref).
"""
supports_shift_constant(::Type{S}) where {S<:MOI.AbstractSet} = false

function shift_constant(set::MOI.LessThan, offset)
    return MOI.LessThan(MOI.constant(set) + offset)
end
supports_shift_constant(::Type{<:MOI.LessThan}) = true

function shift_constant(set::MOI.GreaterThan, offset)
    return MOI.GreaterThan(MOI.constant(set) + offset)
end
supports_shift_constant(::Type{<:MOI.GreaterThan}) = true

function shift_constant(set::MOI.EqualTo, offset)
    return MOI.EqualTo(MOI.constant(set) + offset)
end
supports_shift_constant(::Type{<:MOI.EqualTo}) = true

function shift_constant(set::MOI.Interval, offset)
    return MOI.Interval(set.lower + offset, set.upper + offset)
end
supports_shift_constant(::Type{<:MOI.Interval}) = true

"""
    ScalarLinearSet{T}

The union of scalar-valued linear sets with element type `T`.

This is used in the vectorize and scalarize bridges.

See also: [`VectorLinearSet`](@ref).
"""
const ScalarLinearSet{T} =
    Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}}

"""
    VectorLinearSet

The union of vector-valued linear cones.

This is used in the vectorize and scalarize bridges.

See also: [`ScalarLinearSet`](@ref).
"""
const VectorLinearSet = Union{MOI.Zeros,MOI.Nonnegatives,MOI.Nonpositives}

"""
    vector_set_type(::Type{S}) where {S}

A utility function to map scalar sets `S` to their vector equivalents.

This is used in the vectorize and scalarize bridges.

See also: [`scalar_set_type`](@ref).
"""
vector_set_type(::Type{<:MOI.EqualTo}) = MOI.Zeros
vector_set_type(::Type{<:MOI.LessThan}) = MOI.Nonpositives
vector_set_type(::Type{<:MOI.GreaterThan}) = MOI.Nonnegatives

"""
    scalar_set_type(::Type{S}, ::Type{T}) where {S,T}

A utility function to map vector sets `S` to their scalar equivalents with
element type `T`.

This is used in the vectorize and scalarize bridges.

See also: [`vector_set_type`](@ref).
"""
scalar_set_type(::Type{<:MOI.Zeros}, T::Type) = MOI.EqualTo{T}
scalar_set_type(::Type{<:MOI.Nonpositives}, T::Type) = MOI.LessThan{T}
scalar_set_type(::Type{<:MOI.Nonnegatives}, T::Type) = MOI.GreaterThan{T}

"""
    is_diagonal_vectorized_index(index::Base.Integer)

Return whether `index` is the index of a diagonal element in a
[`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set.
"""
function is_diagonal_vectorized_index(index::Base.Integer)
    # See https://en.wikipedia.org/wiki/Triangular_number#Triangular_roots_and_tests_for_triangular_numbers
    perfect_square = 1 + 8index
    return isqrt(perfect_square)^2 == perfect_square
end

"""
    side_dimension_for_vectorized_dimension(n::Integer)

Return the dimension `d` such that
`MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))` is `n`.
"""
function side_dimension_for_vectorized_dimension(n::Base.Integer)
    return div(isqrt(1 + 8n), 2)
end

"""
    trimap(row::Integer, column::Integer)

Convert between the row and column indices of a matrix, to the linear index of
the corresponding element in the triangular representation.

This is most useful when mapping between `ConeSquare` and `ConeTriangle` sets,
e.g., as part of an [`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set.
"""
function trimap(row::Integer, column::Integer)
    if row < column
        return trimap(column, row)
    end
    return div((row - 1) * row, 2) + column
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


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

mutable struct TrafficAnalyticsConfigurationProperties <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    workspaceId::Any # spec type: Union{ Nothing, String } # spec name: workspaceId
    workspaceRegion::Any # spec type: Union{ Nothing, String } # spec name: workspaceRegion
    workspaceResourceId::Any # spec type: Union{ Nothing, String } # spec name: workspaceResourceId
    trafficAnalyticsInterval::Any # spec type: Union{ Nothing, Int32 } # spec name: trafficAnalyticsInterval

    function TrafficAnalyticsConfigurationProperties(;enabled=nothing, workspaceId=nothing, workspaceRegion=nothing, workspaceResourceId=nothing, trafficAnalyticsInterval=nothing)
        o = new()
        validate_property(TrafficAnalyticsConfigurationProperties, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(TrafficAnalyticsConfigurationProperties, Symbol("workspaceId"), workspaceId)
        setfield!(o, Symbol("workspaceId"), workspaceId)
        validate_property(TrafficAnalyticsConfigurationProperties, Symbol("workspaceRegion"), workspaceRegion)
        setfield!(o, Symbol("workspaceRegion"), workspaceRegion)
        validate_property(TrafficAnalyticsConfigurationProperties, Symbol("workspaceResourceId"), workspaceResourceId)
        setfield!(o, Symbol("workspaceResourceId"), workspaceResourceId)
        validate_property(TrafficAnalyticsConfigurationProperties, Symbol("trafficAnalyticsInterval"), trafficAnalyticsInterval)
        setfield!(o, Symbol("trafficAnalyticsInterval"), trafficAnalyticsInterval)
        o
    end
end # type TrafficAnalyticsConfigurationProperties

const _property_map_TrafficAnalyticsConfigurationProperties = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("workspaceId")=>Symbol("workspaceId"), Symbol("workspaceRegion")=>Symbol("workspaceRegion"), Symbol("workspaceResourceId")=>Symbol("workspaceResourceId"), Symbol("trafficAnalyticsInterval")=>Symbol("trafficAnalyticsInterval"))
const _property_types_TrafficAnalyticsConfigurationProperties = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("workspaceId")=>"String", Symbol("workspaceRegion")=>"String", Symbol("workspaceResourceId")=>"String", Symbol("trafficAnalyticsInterval")=>"Int32")
Base.propertynames(::Type{ TrafficAnalyticsConfigurationProperties }) = collect(keys(_property_map_TrafficAnalyticsConfigurationProperties))
Swagger.property_type(::Type{ TrafficAnalyticsConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TrafficAnalyticsConfigurationProperties[name]))}
Swagger.field_name(::Type{ TrafficAnalyticsConfigurationProperties }, property_name::Symbol) =  _property_map_TrafficAnalyticsConfigurationProperties[property_name]

function check_required(o::TrafficAnalyticsConfigurationProperties)
    true
end

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


text:
# INITIAL GUESS FOR STOCHASTIC APPROXIMATION ALGORITHM BASED ON A GP MODEL

# returns initial guess given line penalty parameters
function getInitialGuess(gamma_line::Float64,delta_line::Float64,useGurobi::Bool=true)

# set generator bounds violation penalty parameters
	gamma_gen::Float64 = 20.0
	delta_gen::Float64 = 1.0

# generate scenarios
	numScenarios::Int64 = 500
	omega_d, omega_w = generateScenarios(numScenarios)

# compute sum of demand deviations for each scenario
	sumMeanDemands::Float64 = sum(meanDemand[j] for j=1:numLoads)
	sumDemands = (sum(omega_d,1))'
	sumDemandDeviations = sumDemands - sumMeanDemands*ones(Float64,numScenarios)

# compute the maximum wind power outputs for the scenarios
	pmax_scen = pmax*ones(Float64,1,numScenarios)
	pmax_scen[numRegGen+1:end,:] = omega_w

	if(useGurobi)
		ig_mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0,BarConvTol=1E-06,OptimalityTol=1E-06))
	else
		ig_mod = Model(solver=IpoptSolver(tol=1E-06,print_level=0))
	end

# first-stage variables
	@variable(ig_mod, pmin[i] <= p0_ig[i=1:numGen] <= pmax_FS[i])
	@variable(ig_mod, 0 <= rp_ig[i=1:numGen] <= pmax_FS[i] - pmin[i])
	@variable(ig_mod, 0 <= rm_ig[i=1:numGen] <= pmax_FS[i] - pmin[i])
	@variable(ig_mod, 0 <= alpha_ig[i=1:numGen] <= 1)

# recourse variables
	@variable(ig_mod, p_ig[1:numGen,1:numScenarios])
	@variable(ig_mod, theta_ig[1:numBuses,1:numScenarios])

# auxiliary variables
	@variable(ig_mod, posResViol_ig[1:numGen,1:numScenarios] >= 0)
	@variable(ig_mod, negResViol_ig[1:numGen,1:numScenarios] >= 0)
	@variable(ig_mod, lineViol_ig[1:numLines,1:numScenarios] >= 0)
	@variable(ig_mod, genBoundViol_ig[1:numGen,1:numScenarios] >= 0)

	@objective(ig_mod, Min, sum(genCost[i]*p0_ig[i] + posResCost[i]*rp_ig[i] + negResCost[i]*rm_ig[i] for i=1:numGen) + (1.0/numScenarios)*sum(gamma_res*posResCost[i]*posResViol_ig[i,scen] + gamma_res*negResCost[i]*negResViol_ig[i,scen] for scen=1:numScenarios for i = 1:numGen) + (1.0/numScenarios)*sum(gamma_line*(lineViol_ig[k,scen])^2 for k = 1:numLines for scen = 1:numScenarios) + (1.0/numScenarios)*sum(gamma_gen*(genBoundViol_ig[i,scen])^2 for i = 1:numGen for scen = 1:numScenarios))

# first-stage constraints. same as in the reserve saturation model
	@constraint(ig_mod, sum(p0_ig[i] for i=1:numGen) == sumMeanDemands)
	@constraint(ig_mod, [i=1:numRegGen], p0_ig[i] + rp_ig[i] <= pmax[i])
	@constraint(ig_mod, [i=1:numRegGen], p0_ig[i] - rm_ig[i] >= pmin[i])
	@constraint(ig_mod, [i=1:numGen; genGuarRes[i]], alpha_ig[i] >= min_alpha)
	@constraint(ig_mod, [i=1:numRegGen; !genRes[i]], alpha_ig[i] == 0)
	@constraint(ig_mod, sum(alpha_ig[i] for i=1:numGen) == 1)

	if(caseID == "no_wind_reserves")
		@constraint(ig_mod, [i=1:numWindGen], alpha_ig[i+numRegGen] == min_alpha/10.0)
	end

# constraints defining p: affine recourse policy without worrying about generation limits
	@constraint(ig_mod, [i=1:numGen,scen=1:numScenarios], p_ig[i,scen] == p0_ig[i] + alpha_ig[i]*sumDemandDeviations[scen])

# theta constraints
	@constraint(ig_mod, [i=1:numBuses,scen=1:numScenarios], sum(betaAdjBus[i][k]*(theta_ig[i,scen] - theta_ig[linesAdjBus[i][k],scen]) for k = 1:size(linesAdjBus[i],1)) == sum(p_ig[regGenAtBus[i][k],scen] for k = 1:size(regGenAtBus[i],1)) + sum(p_ig[windGenAtBus[i][k],scen] for k = 1:size(windGenAtBus[i],1)) - sum(omega_d[loadAtBus[i][k],scen] for k = 1:size(loadAtBus[i],1)))
# reference bus
	@constraint(ig_mod, [scen=1:numScenarios], theta_ig[numBuses,scen] == 0)

# line violations
	@constraint(ig_mod, [k=1:numLines,scen=1:numScenarios], lineViol_ig[k,scen] >= beta[k]*(theta_ig[lines[k][1],scen] - theta_ig[lines[k][2],scen]) - delta_line*flow_max[k])
	@constraint(ig_mod, [k=1:numLines,scen=1:numScenarios], lineViol_ig[k,scen] >= -beta[k]*(theta_ig[lines[k][1],scen] - theta_ig[lines[k][2],scen]) - delta_line*flow_max[k])

# auxiliary constraints for violations of positive and negative reserves
	@constraint(ig_mod, [i=1:numGen,scen=1:numScenarios], posResViol_ig[i,scen] >= p_ig[i,scen] - p0_ig[i] - rp_ig[i])
	@constraint(ig_mod, [i=1:numGen,scen=1:numScenarios], negResViol_ig[i,scen] >= p0_ig[i] - p_ig[i,scen] - rm_ig[i])

# generator bounds violation
	@constraint(ig_mod, [i=1:numGen,scen=1:numScenarios], genBoundViol_ig[i,scen] >= p_ig[i,scen] - (pmin[i] + delta_gen*(pmax_scen[i,scen]-pmin[i])))
	@constraint(ig_mod, [i=1:numGen,scen=1:numScenarios], genBoundViol_ig[i,scen] >= (pmin[i] + (1-delta_gen)*(pmax_scen[i,scen]-pmin[i])) - p_ig[i,scen])

	solve_status = solve(ig_mod)

	if(solve_status != :Optimal && solve_status != :Suboptimal)
		if(useGurobi)
			println("GP based IG model not solved to (sub)optimality using Gurobi!!!")
		else
			println("GP based IG model not solved to (sub)optimality using Ipopt!!!")
		end
		println("solve_status: ",solve_status)
	end

	p0_initguess = getvalue(p0_ig)
	rp_initguess = getvalue(rp_ig)
	rm_initguess = getvalue(rm_ig)
	alpha_initguess = getvalue(alpha_ig)

	return p0_initguess, rp_initguess, rm_initguess, alpha_initguess, solve_status
end

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
using Oceananigans.Fields: ComputedField
using Statistics

import Oceananigans.Fields: AveragedField

"""
    AveragedField(op::AbstractOperation; dims, data=nothing, operand_data=nothing,
                  recompute_safely=false, recompute_operand_safely=false)

Forms a `ComputedField` to store the result of computing `op`, and returns an
`AveragedField` whose operand is the new `ComputedField`, representing an average over
`dims`.

If `data` is not provided, memory is allocated to store the result of the average.
See `AveragedField(field::AbstractField)`.

The keyword argument `operand_data` can be used to specify memory or scratch space
for the new `ComputedField` data.
"""
function AveragedField(op::AbstractOperation; dims, data=nothing, operand_data=nothing,
                       recompute_safely=false, recompute_operand_safely=false)

    computed = ComputedField(op, data=operand_data, recompute_safely=recompute_operand_safely)

    return AveragedField(computed, dims=dims, data=data, recompute_safely=recompute_safely)
end

"""
    mean(op::AbstractOperation; kwargs...)

Returns an Oceananigans.AveragedField representing the an average over `op`eration.
See `Oceananigans.AbstractField`.
"""
Statistics.mean(op::AbstractOperation; kwargs...) = AveragedField(op; kwargs...)


text:
module PhaseTypeDistributions

using Origin: @origin
using SparseArrays: SparseMatrixCSC, nnz, sparse
using SparseMatrix: SparseCSR, SparseCSC, SparseCOO, spdiag, BlockCOO
using NMarkov: rightbound, poipmf!, poipmf, unif, unifstep!, @scal, @axpy, @dot, convunifstep!, itime
using Deformula: deint

export CF1, GPH
export phllf
export phfit!, phfit
export estep!, mstep!, Estep
export phcdf, phpdf, phccdf
export WeightedSample, PointSample
export initializePH

include("_ph.jl")
include("_phunif.jl")

include("_dist.jl")

include("_data.jl")
include("_emstep.jl")
include("_cf1init.jl")
include("_phfit.jl")

include("_phllf.jl")

end # module


text:
# module TestMutating

using Test, Dates, Temporal

@testset "Mutating" begin
    @testset "General" begin
        A = TS(rand(N,K))
        A[:] = 0
        @test all(A.values .== 0)
    end
    @testset "Single-Element Mutations" begin
        A = TS(rand(N,K))
        @test A.index[end] == today()
        A[1,1] = 1.0
        @test A.values[1,1] == 1.0
        A[1, :B] = 2.0
        s = string(today())
        @test A.values[1,2] == 2.0
        A[s, 3] = 3.0
        @test A.values[end,3] == 3.0
        A[s, :D] = 4.0
        @test A.values[end,end] == 4.0
        t = A.index[1]
        A[t, 1] = 1
        @test A.values[1,1] == 1
        A[t, :A] = 2
        @test A.values[1,1] == 2
    end
    @testset "Column Mutations" begin
        A = TS(rand(N,K))
        A[:,1] = 1.0
        @test all(A.values[:,1] .== 1.0)
        A[:B] = 2.0
        @test all(A.values[:,2] .== 2.0)
        A[:,[false,false,true,false]] = 3.0
        @test all(A.values[:,3] .== 3.0)
        A[:] = 4.0
        @test all(A.values .== 4.0)
        A[:,:] = 5.0
        @test all(A.values .== 5.0)
        A[:,1:2] = 0.0
        @test all(A.values[:,1:2] .== 0.0)
        A[[:A,:B]] = -1.0
        @test all(A.values[:,1:2] .== -1.0)
        A[:,[:A,:B]] = -2.0
        @test all(A.values[:,1:2] .== -2.0)
        idx = falses(K)
        idx[1] = true
        X[:,idx] = -3.0
        @test all(X.values[:,1] .== -3.0)
    end
    @testset "Row Mutations" begin
        A = TS(rand(N,K))
        A[1] = 1.0
        @test all(A.values[1,:] .== 1.0)
        A[2,:] = 2.0
        @test all(A.values[2,:] .== 2.0)
        idx = falses(N)
        idx[3] = true
        A[idx] = 3.0
        @test all(A.values[3,:] .== 3.0)
        idx[3] = false
        idx[4] = true
        A[idx,:] = 4.0
        @test all(A.values[4,:] .== 4.0)
        A[today()] = 0.0
        @test all(A.values[end,:] .== 0.0)
        yr = string(year(today()))
        A[yr] = 3.14
        @test all(A[yr].values .== 3.14)
        A[yr,:] = -3.14
        @test all(A[yr].values .== -3.14)
        A[1:2] = 100.0
        @test all(A.values[1:2,:] .== 100.0)
        A[1:2] = -100.0
        @test all(A.values[1:2,:] .== -100.0)
        t = A.index[1]
        A[t] = 0.0
        @test all(A.values[1,:] .== 0.0)
        A[t,:] = 1.0
        @test all(A.values[1,:] .== 1.0)
        t = A.index[1:2]
        A[t] = 3.0
        @test all(A.values[1:2,:] .== 3.0)
    end
end


text:
# lnprob evaluation for Gaussian

nchild = 3
addprocs(nchild)

@everywhere using parallel
@everywhere using visibilities
@everywhere using image
@everywhere using gauss_model
@everywhere using gridding
@everywhere using constants

# Dimensions of the problem
@everywhere const nx = 128
@everywhere const ny = 128
@everywhere const ra = fftspace(6., nx) # [arcsec]
@everywhere const dec = fftspace(6., ny) # [arcsec]

# convert ra and dec in [arcsec] to radians, and then take the sin to convert to ll, mm
@everywhere ll = sin(ra * arcsec) # direction cosines
@everywhere mm = sin(dec * arcsec)

@everywhere function initfunc(key)
    println("Loading dataset $key")
    return DataVis("data/V4046Sgr_fake.hdf5", key)
end

@everywhere function f(dv::DataVis, key::Int, p::Vector{Float64})
    # p is naturally in arcsec

    # Compute the image
    img = imageGauss(ll, mm, p, key)
    skim = SkyImage(img, ra, dec, dv.lam)

    # Apply the gridding correction function before doing the FFT
    corrfun!(skim, 1.0) #alpha = 1.0

    # FFT the appropriate image channel
    vis_fft = transform(skim)

    # Interpolate the `vis_fft` at the same locations as the DataSet
    mvis = ModelVis(dv, vis_fft)

    # Calculate chi^2 between these two
    return lnprob(dv, mvis)
end

keylist = Int[1, 2, 3]
pipes = initialize(nchild, keylist, initfunc, f)

function fprob(p::Vector{Float64})
    distribute!(pipes, p)
    return gather!(pipes)
end

println(fprob([1.2, 1.0, 1.2, 1.0]))

# Necessary wrapper for NLopt requires gradient as an argument (even if it's not used)
function fgrad(p::Vector, grad::Vector)
    val = fprob(p)
    println(p, " : ", val)
    return val
end
#
# function fp(p::Vector)
#     val = f(p)
#     println(p, " : ", val)
#     return val
# end

# # Now try optimizing the function using NLopt
using NLopt

starting_param = [1.3, 1.2, 0.7, 0.7]

nparam = length(starting_param)
opt = Opt(:LN_COBYLA, nparam)

max_objective!(opt, fgrad)
xtol_rel!(opt,1e-4)

(optf,optx,ret) = optimize(opt, starting_param)
println(optf, " ", optx, " ", ret)

# using LittleMC
#
# using Distributions
# using PDMats
#
# mc = MC(fp, 500, [1.2, 1.0, 1.0, 1.0], PDiagMat([0.03^2, 0.03^2, 0.01^2, 0.01^2]))
#
# start(mc)
#
#
# println(mean(mc.samples, 2))
# println(std(mc.samples, 2))
#
# runstats(mc)
#
# write(mc, "mc.hdf5")

quit!(pipes)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
module DuffyQuadraturesTests

using Test
using Gridap

D = 3
order=2
quad = DuffyQuadrature{D}(order)
@test isa(coordinates(quad),Vector{Point{D,Float64}})
@test isa(weights(quad),Vector{Float64})
@test sum(weights(quad)) ≈ 1/6

D = 2
order=2
quad = DuffyQuadrature{D}(order)
x = coordinates(quad)
r = Point{D,Float64}[
  (0.155051, 0.178559), (0.644949, 0.0750311),
  (0.155051, 0.66639), (0.644949, 0.28002)]
@test isapprox(reinterpret(Float64,x),reinterpret(Float64,r),rtol=1e-3)
@test sum(weights(quad)) ≈ 1/2

D = 1
order=2
quad = DuffyQuadrature{D}(order)
@test sum(weights(quad)) ≈ 1

end # module


text:
#=using Pkg
import Pkg; Pkg.add("MathOptInterface")
import Pkg; Pkg.add("Statistics")
import Pkg; Pkg.add("QuantEcon")
import Pkg; Pkg.add("Interpolations")
import Pkg; Pkg.add("NLsolve")
import Pkg; Pkg.add("Optim")
import Pkg; Pkg.add("Random")
import Pkg; Pkg.add("IterTools")
import Pkg; Pkg.add("JuMP")
import Pkg; Pkg.add("Ipopt")
import Pkg; Pkg.add("BenchmarkTools")
import Pkg; Pkg.add("GaussianProcesses")
=#

using
    MathOptInterface,
    BenchmarkTools,
    LinearAlgebra,
    Statistics,
    QuantEcon,
    Interpolations,
    NLsolve,
    Optim,
    Random,
    IterTools,
    JuMP,
    Ipopt,
    GaussianProcesses,
    TimerOutputs,
    Random;

δ = 0.1
β = 0.86
ϕ = 0.6
ξ = 0.5
μ = 0.5
ψ = 1/ϕ
b = 1/(ϕ^ϕ * (1-ϕ)^(1-ϕ))

numSec = 2
numLin = 10
gridMin = 2
gridMax = 30
gridSiz = numLin^numSec
u(x,y) = (log(x) + log(y))
wInit(x) = (log(x[1]) + log(x[2]))/(1 - β)
function gradwInit(g, x, y)
    g[1] = 1/x/(1-β)
    g[2] = 1/y/(1-β)
end
function hesswInit(h, x, y)
    h[1,1] = -1/x^2/(1-β)
    h[1,2] = 0
    h[2,1] = 0
    h[2,2] = -1/y^2/(1-β)
end
yTargCbd = Vector{Float64}[]
sTargCbd = Vector{Float64}[]
mTargCbd = Matrix{Float64}[]
xTargCbd = Matrix{Float64}[]
for i in 1:gridSiz
    push!(yTargCbd, zeros(numSec))
    push!(sTargCbd, zeros(numSec))
    push!(mTargCbd, zeros(numSec,numSec))
    push!(xTargCbd, zeros(numSec, numSec))
end
wTargCbd = Float64[]
#randvec = 20 .+ 30*rand(rng,numLin)
#for p in product(randvec, randvec)
#    push!(grid, collect(p))
#end
grid = Vector{Float64}[]
for p in product(LinRange(gridMin, gridMax, numLin),LinRange(gridMin, gridMax, numLin))
    push!(grid, collect(p))
end
rng = MersenneTwister(1234)
#for timing
tOCbd = TimerOutput()
#wval = w.(grid)
#wval = w.(grid)

#the following function generates target values (of the value function) for the ML process
##note that (unlike Sargent and Stachursky) the grid is written in terms of kapital
#this simplified the problem substantially as otherwise intermediate variables became dynamic
##also note that this accepts the smooth function
function TargCbd(w, grid, β ; return_policy= false)
    @timeit tOCbd string("Combined-for-", length(grid),"-pts") begin
    modTrial = Model(Ipopt.Optimizer)
    set_silent(modTrial)
    @variable(modTrial,  c[1:numSec, 1:length(grid)] >= 0.0001)
    @variable(modTrial, y[1:numSec, 1:length(grid)] >= 0.0001)
    @variable(modTrial, kp[1:numSec, 1:length(grid)]>= 0.0001)
    @variable(modTrial, x[1:numSec, 1:numSec, 1:length(grid)]>=0.0001)
    @variable(modTrial, xComb[1:numSec, 1:length(grid)]>=0.0001)
    @variable(modTrial, m[1:numSec, 1:numSec, 1:length(grid)]>=0.0001)
    @variable(modTrial, mComb[1:numSec, 1:length(grid)]>=0.0001)
    @variable(modTrial, u[1:length(grid)])
    @variable(modTrial, w[1:length(grid)])
    k = grid
    for n in 1:length(grid)
        for i in 1:numSec
            #output
            @NLconstraint(modTrial, y[i,n] == b * k[n][i]^ϕ * mComb[i,n]^(1-ϕ))
            #future capital
            @constraint(modTrial, kp[i,n] == xComb[i,n] + (1-δ) * k[n][i])
            #supply equals demand (eq'm constraint)
            @constraint(modTrial, 0 == y[i,n] - c[i,n] - sum(m[i,:,n]) - sum(x[i,:,n]))
        end
        @NLconstraint(modTrial, xComb[1,n] == x[1,1,n]^ξ * x[2,1,n]^(1-ξ))
        @NLconstraint(modTrial, xComb[2,n] == x[1,2,n]^ξ * x[2,2,n]^(1-ξ))
        @NLconstraint(modTrial, mComb[1,n] == m[1,1,n]^μ * m[2,1,n]^(1-μ))
        @NLconstraint(modTrial, mComb[2,n] == m[1,2,n]^μ * m[2,2,n]^(1-μ))
        @NLconstraint(modTrial, u[n] == log(c[1,n]) + log(c[2,n]))
        @NLconstraint(modTrial, w[n] == (log(kp[1,n]) + log(kp[2,n]))/(1-β))
#=             register(modTrial, :wInit, 2, wInit, gradwInit)
            @NLconstraint(modTrial, w == wInit(kp[1],kp[2]))
            register(modTrial, :u, 2, u; autodiff = true) =#
    end
    @NLobjective(modTrial, Max, sum(u[i] + β * w[i] for i in 1:length(grid)))
    JuMP.optimize!(modTrial)
    #push!(wTargCbd, value.(w))
    if return_policy
        for i in 1:length(grid)
        sTargCbd[i] = value.(c)[1:2, i]
        yTargCbd[i] = value.(y)[1:2, i]
        mTargCbd[i] = value.(m)[1:2,1:2, i]
        xTargCbd[i] = value.(x)[1:2,1:2, i]
        end
    end
                   end #of timeit
    #global intK[n] = value.(f)
    #global prevK = intK
    return value.(w)
end

function evalBellCbd()
    wTargCbd = TargCbd(wInit, grid, β; return_policy=false)
    show(tOCbd)
    println("\n ---------")
    return wTargCbd
end
#=XTargCbd = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(XTargCbd, sum(xTargCbd[j][:,i] for i in 1:numSec))
       end
MTargCbd = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(MTargCbd, sum(mTargCbd[j][:,i] for i in 1:numSec))
       end
fTargCbd = sTargCbd + MTargCbd + XTargCbd

MPCf = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(MPCf, sTargCbd[j]./fTargCbd[j])
       end
MPC = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(MPC, sTargCbd[j]./yTargCbd[j])
       end
MPM = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(MPM, MTargCbd[j]./yTargCbd[j])
       end
MPX = Vector{Float64}[]
       for j = 1:numLin^numSec
           push!(MPX, XTargCbd[j]./yTargCbd[j])
       end
# a function for turning vectors of vectors into matrices
#function vvm(x)
#           dim1 = length(x)
#           dim2 = length(x[1])
#           matrix = zeros(Float64, dim1, dim2)
#           @inbounds @fastmath for i in 1:dim1, for j in 1:dim2
#                   matrix[i, j] = x[i][j]
#               end
#           end
#           return matrix
#end
#alternatively
gridM = reduce(hcat, grid)

mz = MeanZero()
ml = MeanLin([2., 3.])
mp = MeanPoly([2. 3.; 1. 2.])
kern = SE(0.0, 0.0)
gp = GP(gridM, wTargCbd, mz, kern)
#gp = GP(gridM, wTargCbd, MeanLin([2., 3.]), kern)
GaussianProcesses.optimize!(gp)
predict_f(gp, gridM)[1]
fstar(x) = predict_f(gp,hcat(x))[1]

plot!(contour(gp), heatmap(gp); fmt=:png)
#wTargp(x, y) = predict_f(gp,hcat([x, y]))[1]=#

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
# coding: utf-8

# In[1]:

#fractal is one of the interesting topics in geometry
#it is usually described by a recursive function
#voila,here we are!
using Plots

# In[2]:

#create rectangle shape
rectangle(top_left,top_right,bottom_left,
        bottom_right)=Shape([(top_left[1],top_left[2]),
        (top_right[1],top_right[2]),
        (bottom_right[1],bottom_right[2]),
        (bottom_left[1],bottom_left[2])])

# In[3]:

#compute euclidean distance
function euclidean_distance(point1,point2)
    return √((point1[1]-point2[1])^2+(point1[2]-point2[2])^2)
end

# In[4]:

#simple solution to get coefficients of the equation
function get_line_params(x1,y1,x2,y2)

    slope=(y1-y2)/(x1-x2)
    intercept=y1-slope*x1

    return slope,intercept

end

# In[5]:

#standard solution to quadratic equation
function solve_quadratic_equation(A,B,C)
    x1=(-B+√(B^2-4*A*C))/(2*A)
    x2=(-B-√(B^2-4*A*C))/(2*A)
    return [x1,x2]
end

# In[6]:

#analytic geometry to compute target datapoints
function get_datapoint(pivot,measure,length,direction="inner")

    #for undefined slope
    if pivot[1]==measure[1]
        y1=pivot[2]+length
        y2=pivot[2]-length
        x1=pivot[1]
        x2=pivot[1]

    #for general cases
    else

        #get line equation
        slope,intercept=get_line_params(pivot[1],pivot[2],
                                   measure[1],measure[2],)

        #solve quadratic equation
        A=1
        B=-2*pivot[1]
        C=pivot[1]^2-length^2/(slope^2+1)
        x1,x2=solve_quadratic_equation(A,B,C)

        #get y from line equation
        y1=slope*x1+intercept
        y2=slope*x2+intercept

    end

    if direction=="inner"

        #take the one between pivot and measure points
        if euclidean_distance((x1,y1),measure)<euclidean_distance((x2,y2),measure)
            datapoint=(x1,y1)
        else
            datapoint=(x2,y2)
        end

    else

        #take the one farther away from measure points
        if euclidean_distance((x1,y1),measure)>euclidean_distance((x2,y2),measure)
            datapoint=(x1,y1)
        else
            datapoint=(x2,y2)
        end

    end

    return datapoint

end

# In[7]:

#recursively plot symmetric pythagorean tree at 45 degree
# https//larryriddle.agnesscott.org/ifs/pythagorean/pythTree.htm
function pythagorean_tree(top_left,top_right,bottom_left,
        bottom_right,current_angle,line_len,n)

    #plot square
    plot!(rectangle(top_left,top_right,bottom_left,bottom_right))

    if n==0
        return
    else

        #find mid point
        #midpoint has to satisfy two conditions
        #it has to be on the same line as bottom_left and bottom_right
        #assume this line follows y=kx+b
        #the midpoint is (x,kx+b)
        #bottom_left is (α,kα+b),bottom_right is (δ,kδ+b)
        #the euclidean distance between midpoint and bottom_left should be
        #half of the euclidean distance between bottom_left and bottom_right
        #(x-α)**2+(kx+b-kα-b)**2=((α-δ)**2+(kα+b-kδ-b)**2)/4
        #apart from x,everything else in the equation is constant
        #this forms a simple quadratic solution to get two roots
        #one root would be between bottom_left and bottom_right which yields midpoint
        #and the other would be farther away from bottom_right
        #this function solves the equation via (-B+(B**2-4*A*C)**0.5)/(2*A)
        #alternatively,you can use scipy.optimize.root
        #the caveat is it does not offer both roots
        #a wrong initial guess could take you to the wrong root
        bottom_mid=get_datapoint(bottom_left,bottom_right,line_len/2)
        top_mid=get_datapoint(top_left,top_right,line_len/2)

        #compute the top point of a triangle
        #the computation is similar to midpoint
        #the euclidean distance between triangle_top and top_mid should be
        #half of the distance between top_mid and bottom_mid
        triangle_top=get_datapoint(top_mid,bottom_mid,
                                   line_len/2,"outer")

        #get top left for right square
        #the computation is similar to midpoint
        #the euclidean distance between triangle_top and rightsq_topleft
        #should be the same as the distance between triangle_top and top_left
        rightsq_topleft=get_datapoint(triangle_top,top_left,
                                      line_len/(√(2)),"outer")

        #get midpoint of the diagonal between rightsq_topleft and top_right
        #the computation is similar to midpoint
        #the euclidean distance between rightsq_diag_mid and rightsq_topleft
        #should be half of the distance between rightsq_topleft and top_right
        rightsq_diag_mid=get_datapoint(top_right,rightsq_topleft,line_len/2)
        rightsq_topright=get_datapoint(rightsq_diag_mid,triangle_top,
                                       line_len/2,"outer")

        #get top left and right for left square similar to right square
        leftsq_topleft=get_datapoint(triangle_top,rightsq_topright,
                                     line_len,"outer")
        leftsq_topright=get_datapoint(triangle_top,top_right,
                                      line_len/(√(2)),"outer")

        #recursive do the same for left square
        pythagorean_tree(leftsq_topleft,leftsq_topright,
                         top_left,triangle_top,current_angle+45,
                         line_len/(√(2)),n-1)

        #recursive do the same for right square
        pythagorean_tree(rightsq_topleft,rightsq_topright,
                         triangle_top,top_right,current_angle-45,
                         line_len/(√(2)),n-1)

    end

end

# In[8]:

#initialize
top_left=(0,0)
top_right=(1,0)
bottom_left=(0,-1)
bottom_right=(1,-1)
n=4
current_angle=0
line_len=euclidean_distance(top_left,top_right);

# In[9]:

#viz
gr(size=(250,200))
fig=plot(legend=false,grid=false,axis=false,ticks=false,
    )

pythagorean_tree(top_left,top_right,bottom_left,
                     bottom_right,current_angle,line_len,n)
fig

max_stars_repo_path: pythagorean tree.jl

max_issues_repo_path: pythagorean tree.jl

max_forks_repo_path: pythagorean tree.jl


text:
module StableArrays

using Base
using LinearAlgebra

export StableArray, stabilize!, stabilize, unstabilize

include("StableNumber.jl")

"""
    StableArray(base::AbstractArray, exponent)
"""
struct StableArray{E, A <: AbstractArray}
    base::A
    exponent::E
end

base(x::StableArray) = x.base
exponent(x::StableArray) = x.exponent

"""
    stabilize!(x::AbstractArray)

Construct a [`StableArray`](@ref)
"""
function stabilize!(x::AbstractArray)
    factor = maximum(abs, x)
    if factor == 0
        StableArray(x, factor)
    else
        x ./= factor
        StableArray(x, log(factor))
    end
end

function stabilize(x::AbstractArray)
    factor = maximum(abs, x)
    if factor == 0
        StableArray(x, factor)
    else
        x /= factor
        StableArray(x, log(factor))
    end
end

function unstabilize(x::StableArray)
    base(x) * exp(exponent(x))
end

function unstabilize(x::StableNumber)
    base(x) * exp(exponent(x))
end

const StableMatrix{E, A} = StableArray{E, A} where A<:AbstractMatrix

function Base.:(*)(x1::StableArray{E}, x2::StableArray{E}) where E
    x = x1.base * x2.base
    factor = maximum(abs, x)
    if factor == 0
        return StableArray(x, factor)
    end
    x ./= factor
    xexponent = x1.exponent + x2.exponent + log(convert(E, factor))
    StableArray(x, xexponent)
end

function Base.:(*)(x1::StableArray, x2)
    x1 * stabilize(x2)
end

function Base.:(*)(x1, x2::StableArray)
    stabilize(x1) * x2
end

function Base.:(^)(x::StableArray, n::Integer)
    prod((x for _ in 1:n))
end

function Base.maximum(x::StableArray)
    StableNumber(maximum(base(x)), exponent(x))
end

function Base.adjoint(x::StableArray)
    StableArray(x.base', x.exponent)
end

function Base.size(x::StableArray, args...)
    size(x.base, args...)
end

function LinearAlgebra.tr(x::StableArray)
    StableNumber(tr(base(x)), exponent(x))
end

end


text:
function statsuri(cli::Client, stats_type::String)
    return "http://$(cli.host):$(cli.port)/$(cli.version)/stats/$(stats_type)"
end

"""
    stats(cli, stats_type) -> Dict

Fetches any stats from the cluster (`stats_type` may be "store", "self" or "leader")
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#statistics)
for more details.
"""
function stats(cli::Client, stats_type::String)
    return request(Requests.get, statsuri(cli, stats_type), Dict())
end


text:
module NMRquant

using Dates
using LinearAlgebra: norm
using Mustache
using NMR
using RecipesBase
using Pkg
using Printf
using Sockets
using Statistics
using TimeZones

function write_config(path)
    open(joinpath(path, "config"), "w") do f
        write(f, match(r"inet ([\d\.]+)", read(`ip address show enp0s3`, String)).captures[1])
    end
end

function sorted_kv(dict,args...)
    k = sort(collect(keys(dict)), args...)
    ((x, dict[x]) for x in k)
end

function dark_signals(s::NMR.Spectrum, dark_areas)
    darks = NMR.extract.(s, [area[1] for area in dark_areas])
    for (d,name) in zip(darks, [area[2] for area in dark_areas])
        d[d.default_proc].title = "Name = $name"
    end
    darks
end

function remove_darks!(s::NMR.Spectrum, darks::Array{NMR.Spectrum,1})
    d = NMR.decompose(NMR.lsq_analyze(s, darks)) # find dark signals in s
    s[:] = d[3] # residue (non-dark part)
end

function remove_darks!(lib)
    darks = [l for l in lib if !haskey(sampleinfo(l), "Protons")]
    for l in lib
        if haskey(sampleinfo(l), "Protons")
            remove_darks!(l, darks)
        end
    end
end

"""
    analyze(s::NMR.Spectrum, lib::Array{NMR.Spectrum,1}, stan::Spectrum; out = stdout)

Perform analysis on the NMR spectrum `s` against the library group `lib`, using the internal standard spectrum `stan` for quantitation.

Returns the tuple (`decomposition results`, `quantitation results`) of type (DecompositionResult, Dict{Int, Array{Float64}}).
`Quantitation results` is a dictionary with library reference numbers as keys, and an array of % quantitations as values.
"""
function analyze(s::NMR.Spectrum, lib::Array{NMR.Spectrum,1}, stan::Spectrum; out = stdout, tol = NMR.TOL)
    callback(refs) = println(out, "Found $(join([name(lib[r]) for r in refs], ", ")).")
    d = NMR.lsq_analyze(s, lib; callback=callback, tol=tol)
    quants = quantitate(s, lib, d, stan)
    (d, quants)
end

"""
    quantitate(s::NMR.Spectrum, lib::Vector{NMR.Spectrum}, d::NMR.DecompositionResult, stan::Spectrum, minquant = quant_threshold(s))

Return a Dictionary of quantitation results based on DecompositionResult `d`, reported as percent by weight based internal standard spectrum `stan`,
library spectra `lib` and unknown spectrum `s`. An identification is only reported if the result is above the quantitation threshold `minquant`.
"""
function quantitate(s::NMR.Spectrum, lib::Vector{NMR.Spectrum}, d::NMR.DecompositionResult, stan::Spectrum, minquant = quant_threshold(s))
    results = Dict{Int,Array{Float64}}()
    sampinfo = sampleinfo(s)
    stanint = NMR.integrate(s, stan[stan.default_proc].intrng[1])
    stanconc = blank_conc(s) / MW(stan)
    stanprotons = protons(stan)[1]
    sampleweight = parse(Float64, sampinfo["Sample weight (mg)"])
    volume = parse(Float64, sampinfo["Volume (mL)"])
    for (refno, ref, coeff) in zip(d.refnums, lib[d.refnums], d.coefficients)
        if !isdark(ref)
            pr = protons(ref)
            mw = MW(ref)
            q = NMR.integrate(ref)./pr.*(coeff/stanint*stanprotons*mw*stanconc*volume/sampleweight*100.0)
            if mean(q) > max(minquant, quant_threshold(ref)) # only report more than meanquant %
                results[refno] = q
            end
        end
    end
    results
end

"""
    serve(port, home, libs)

Establish a listening port specified by `port` for receiving NMR spectra data located in a directory defined by `home`.

Valid NMR spectra will be analyzed against a library group `lib`, which must be a Dictionary with library names as keys,
and ([Library Spectrum], [Dark Spectrum], IS Spectrum) as values. A report will be generated at the location of the NMR spectra.
"""
function serve(port, home, libs)
    @async begin
        server = listen(ip"0.0.0.0", port)
        while true
            sock = accept(server)
            @info "Connection opened"
            l = strip(readline(sock))
            @info "Received $l"
            if strip(l) == "bye"
                @info "Shutting down server"
                close(sock)
                close(server)
                return 0
            end
            try
                s = NMR.Spectrum(joinpath(home,l),1)
                sampleinfovalid(s, libs)
                lib, darks, is = getlib(s, libs)

                #3.4
                if lowercase(get(sampleinfo(lib[1]), "Reference Deconvolution", "false")) ∈ ["true", "yes"]
                    tm1 = get(sampleinfo(lib[1]), "tm1", NMR.TM1)
                    tm2 = get(sampleinfo(lib[1]), "tm2", NMR.TM2)
                    s = refdecon(s, lib[3];tm1=tm1, tm2=tm2)
                    s[1].procno = 99
                    s[1].title *= "\r\ntm1 = $tm1\r\ntm2 = $tm2"
                    dump(joinpath(home,l,"pdata","99"),s[1])
                end

                d,q = analyze(s, lib, is; out = sock)
                outfile = report(s, lib, d, q, darks, is, joinpath(home,l))
                println(sock, "Done")
            catch err
                @error err
                @error stacktrace(catch_backtrace())
                println(sock, err)
            finally
                close(sock)
                @info "Connection closed, listening"
            end
        end
    end
end

include("plots.jl")
include("report.jl")
include("util.jl")

export serve, analyze, report, library_report, quantitate, sampleinfo, library

end


text:
export Export
export sourceMap, targetMap, distributor, isLocallyComplete
export permuteToLIDs, permuteFromLIDs, exportLIDs, remoteLIDs, remotePIDs, numSameIDs

"""
Communication plan for data rekg.ribution from a (possibly) multiple-owned to a uniquely owned distribution
"""
struct Export{GID <: Integer, PID <:Integer, LID <: Integer}
    exportData::ImportExportData{GID, PID, LID}
end

#TODO document

## Constructors ##

function Export(source::BlockMap{GID, PID, LID}, target::BlockMap{GID, PID, LID}, remotePIDs::Union{AbstractArray{PID}, Nothing}=nothing; plist...) where {GID <: Integer, PID <: Integer, LID <: Integer}
    Export(source, target,
        Dict(Tuple{Symbol, Any}[pair for pair in plist]))
end

function Export(source::BlockMap{GID, PID, LID}, target::BlockMap{GID, PID, LID},
        plist::Dict{Symbol}) where {GID <: Integer, PID <: Integer, LID <: Integer}
    Export(source, target, nothing, plist)
end

function Export(source::BlockMap{GID, PID, LID}, target::BlockMap{GID, PID, LID},
        remotePIDs::Union{AbstractArray{PID}, Nothing}, plist::Dict{Symbol}) where {
            GID <: Integer, PID <: Integer, LID <: Integer}

    expor = Export(ImportExportData(source, target))

    exportGIDs = setupSamePermuteExport(expor)

    if distributedGlobal(source)
        setupRemote(expor, exportGIDs)
    end

    expor
end

## internal construction methods ##

function setupSamePermuteExport(expor::Export{GID, PID, LID})::AbstractArray{GID} where {GID <: Integer, PID <: Integer, LID <: Integer}

    data = expor.exportData

    source = sourceMap(expor)
    target = targetMap(expor)

    sourceGIDs = myGlobalElements(source)
    targetGIDs = myGlobalElements(target)

    numSrcGIDs = length(sourceGIDs)
    numTgtGIDs = length(targetGIDs)
    numGIDs = min(numSrcGIDs, numTgtGIDs)

    numSameGIDs = 1
    while numSameGIDs <= numGIDs && sourceGIDs[numSameGIDs] == targetGIDs[numSameGIDs]
        numSameGIDs += 1
    end
    numSameGIDs -= 1
    numSameIDs(data, numSameGIDs)

    exportGIDs = Array{GID, 1}(undef, 0)
    permuteToLIDs = JuliaPetra.permuteToLIDs(data)
    permuteFromLIDs = JuliaPetra.permuteFromLIDs(data)
    exportLIDs = JuliaPetra.exportLIDs(data)

    for srcLID = (numSameGIDs+1):numSrcGIDs
        curSrcGID = sourceGIDs[srcLID]
        tgtLID = lid(target, curSrcGID)
        if tgtLID != 0
            push!(permuteToLIDs, tgtLID)
            push!(permuteFromLIDs, srcLID)
        else
            push!(exportGIDs, curSrcGID)
            push!(exportLIDs, srcLID)
        end
    end

    if length(exportLIDs) != 0 && !distributedGlobal(source)
        isLocallyComplete(data, false)
        warn("Source has export LIDs but source not distributed globally.  " *
            "Exporting to a submap of the target map.")
    end

    if distributedGlobal(source)
        #resize!(JuliaPetra.exportPIDs(data), length(exportGIDs))

        (exportPIDs, exportLIDs) = remoteIDList(target, exportGIDs)
        JuliaPetra.exportPIDs(data, exportPIDs)
        missingGIDs = 0
        for i = 1:length(exportPIDs)
            if exportPIDs[i] == 0
                missingGIDs += 1
            end
        end

        if missingGIDs != 0
            warn("The source Map has GIDs not found in the target Map")

            isLocallyComplete(data, false)
            numInvalidExports = missingGIDs
            totalNumExports = length(exportPIDs)

            if numInvalidExports == totalNumExports
                # all exports invalid, can delete all exports
                resize!(exportGIDs, 0)
                resize!(exportLIDs, 0)
                resize!(exportPIDs, 0)
            else
                #some exports are valid, need to keep the valid exports
                numValidExports = 1
                for e = 1:totalNumExports
                    if exportPIDs[e] != 0
                        exportGIDs[numValidExports] = exportGIDs[e]
                        exportLIDs[numValidExports] = exportLIDs[e]
                        exportPIDs[numValidExports] = exportPIDs[e]
                        numValidExports += 1
                    end
                end
                numValidExports -= 1

                resize!(exportGIDs, numValidExports)
                resize!(exportLIDs, numValidExports)
                resize!(exportPIDs, numValidExports)
            end
        end
    end
    exportGIDs
end

function setupRemote(expor::Export{GID, PID, LID}, exportGIDs::AbstractArray{GID, 1}) where {GID <: Integer, PID <: Integer, LID <: Integer}

    data = expor.exportData

    target = targetMap(data)

    exportPIDs = JuliaPetra.exportPIDs(data)

    order = sortperm(exportPIDs)
    permute!(exportPIDs, order)
    permute!(exportLIDs(data), order)
    permute!(exportGIDs, order)

    numRemoteIDs = createFromSends(distributor(data), exportPIDs)

    remoteGIDs = resolve(distributor(data), exportGIDs)

    remoteLIDs = JuliaPetra.remoteLIDs(data)

    resize!(remoteLIDs, numRemoteIDs)

    for i in 1:length(remoteGIDs)
        remoteLIDs[i] = lid(target, remoteGIDs[i])
    end

end
## Getters ##

"""
Get the source map for the given ImportExportData
"""
function sourceMap(impor::Export{GID, PID, LID})::BlockMap{GID, PID, LID} where GID <: Integer where PID <:Integer where LID <: Integer
    impor.exportData.source
end

"""
Get the target map for the given ImportExportData
"""
function targetMap(impor::Export{GID, PID, LID})::BlockMap{GID, PID, LID} where GID <: Integer where PID <:Integer where LID <: Integer
    impor.exportData.target
end

"""
List of elements in the target map that are permuted.
"""
function permuteToLIDs(impor::Export{GID, PID, LID})::Array{LID} where GID <: Integer where PID <:Integer where LID <: Integer
    impor.exportData.permuteToLIDs
end

"""
List of elements in the source map that are permuted.
"""
function permuteFromLIDs(impor::Export{GID, PID, LID})::Array{LID} where GID <: Integer where PID <:Integer where LID <: Integer
    impor.exportData.permuteFromLIDs
end

"""
List of elements in the target map that are coming from other processors
"""
function remoteLIDs(impor::Export{GID, PID, LID})::Array{LID} where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.remoteLIDs
end

"""
List of elements that will be sent to other processors
"""
function exportLIDs(impor::Export{GID, PID, LID})::Array{LID} where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.exportLIDs
end

"""
List of processors to which elements will be sent `exportLID[i]` will be sent to processor `exportPIDs[i]`
"""
function exportPIDs(impor::Export{GID, PID, LID})::Array{PID} where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.exportPIDs
end

"""
Returns the number of elements that are identical between the source and target maps, up to the first different ID
"""
function numSameIDs(impor::Export{GID, PID, LID})::LID where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.numSameIDs
end

"""
Returns the distributor being used
"""
function distributor(impor::Export{GID, PID, LID})::Distributor{GID, PID, LID} where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.distributor
end

"""
Returns whether the import or export is locally complete
"""
function isLocallyComplete(impor::Export{GID, PID, LID})::Bool where GID <: Integer where PID <: Integer where LID <: Integer
    impor.exportData.isLocallyComplete
end


text:
using HDF5
using Compat.Test

@testset "extendtest" begin

fn = tempname()

fid = h5open(fn, "w")
g = g_create(fid, "shoe")
d = d_create(g, "foo", datatype(Float64), ((10, 20), (100, 200)), "chunk", (1, 1))
#println("d is size current $(map(int,HDF5.get_dims(d)[1])) max $(map(int,HDF5.get_dims(d)[2]))")
dims, max_dims = HDF5.get_dims(d)
@test dims == (UInt64(10), UInt64(20))
@test max_dims == (UInt64(100), UInt64(200))
set_dims!(d, (100, 150))
dims, max_dims = HDF5.get_dims(d)
@test dims == (UInt64(100), UInt64(150))
@test max_dims == (UInt64(100), UInt64(200))
d[1, 1:5] = [1.1231, 1.313, 5.123, 2.231, 4.1231]
set_dims!(d, (1, 5))
@test size(d) == (1, 5)
#println("d is size current $(map(int,HDF5.get_dims(d)[1])) max $(map(int,HDF5.get_dims(d)[2]))")
b = d_create(fid, "b", Int, ((1000,), (-1,)), "chunk", (100,)) #-1 is equivalent to typemax(Hsize) as far as I can tell
#println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")
b[1:200] = ones(200)
dims, max_dims = HDF5.get_dims(b)
@test dims == (UInt64(1000),)
@test max_dims == (HDF5.H5S_UNLIMITED,)
set_dims!(b, (10000,))
dims, max_dims = HDF5.get_dims(b)
@test dims == (UInt64(10000),)
@test max_dims == (HDF5.H5S_UNLIMITED,)
#println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")
# b[:] = [1:10000] # gave error no method endof(HDF5Dataset{PlainHDF5File},),
# so I defined endof(dset::HDF5Dataset) = length(dset), and exported endof
# but that didn't fix the error, despite the endof function working
# d[1] produces error ERROR: Wrong number of indices supplied, should datasets support linear indexing?
b[1:10000] = [1:10000;]
#println(b[1:100])

close(fid)

fid = h5open(fn, "r")
d_again = fid["shoe/foo"]
dims, max_dims = HDF5.get_dims(d_again)
@test dims == (UInt64(1), UInt64(5))
@test max_dims == (UInt64(100), UInt64(200))
@test (sum(d_again[1, 1:5]) - sum([1.1231, 1.313, 5.123, 2.231, 4.1231])) == 0
#println("d is size current $(map(int,HDF5.get_dims(re_d)[1])) max $(map(int,HDF5.get_dims(re_d)[2]))")
@test fid["b"][1:10000] == [1:10000;]
b_again = fid["b"]
dims, max_dims = HDF5.get_dims(b_again)
@test dims == (UInt64(10000),)
@test max_dims == (HDF5.H5S_UNLIMITED,)
#println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")

close(fid)
rm(fn)

end # testset extend_test