text:
function Base.real{DD}(f::Fun{Laurent{DD}})
    n=ncoefficients(f)
    cfs=f.coefficients

    ret=Array(Float64,iseven(n)?n+1:n)
    ret[1]=real(cfs[1])

    for k=2:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]=imag(cfs[k])
        ret[k+1]=real(cfs[k])
    end
    for k=3:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]+=real(cfs[k])
        ret[k-1]-=imag(cfs[k])
    end

    Fun(ret,Fourier(domain(f)))
end

function Base.imag{DD}(f::Fun{Laurent{DD}})
    n=ncoefficients(f)
    cfs=f.coefficients

    ret=Array(Float64,iseven(n)?n+1:n)
    ret[1]=imag(cfs[1])

    for k=2:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]=-real(cfs[k])
        ret[k+1]=imag(cfs[k])
    end
    for k=3:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]+=imag(cfs[k])
        ret[k-1]+=real(cfs[k])
    end

    Fun(ret,Fourier(domain(f)))
end


text:
# Analogous to elixir_euler_source_terms.jl, but uses a locally refined cubed sphere mesh
# of the size of the Earth's atmosphere (using an atmospheric height of 30km).
# The initial condition and source terms have also been rescaled to planetary size.

using OrdinaryDiffEq
using Trixi
using LinearAlgebra

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations3D(gamma)

function initial_condition_convergence_test_sphere(x, t, equations::CompressibleEulerEquations3D)
  x_scaled = x / 6.371229e6
  t_scaled = t / 6.371229e6

  return initial_condition_convergence_test(x_scaled, t_scaled, equations)
end

@inline function source_terms_convergence_test_sphere(u, x, t, equations::CompressibleEulerEquations3D)
  x_scaled = x / 6.371229e6
  t_scaled = t / 6.371229e6

  return source_terms_convergence_test(u, x_scaled, t_scaled, equations) / 6.371229e6
end

function indicator_test(u::AbstractArray{<:Any,5},
                        mesh, equations, dg::DGSEM, cache;
                        kwargs...)
  alpha = zeros(Int, nelements(dg, cache))

  for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, k, element)
      lambda, phi, r = cart_to_sphere(x)
      if 0.22 < lambda < 3.3 && 0.45 < phi < 1.3
        alpha[element] = 1
      end
    end
  end

  return alpha
end

function cart_to_sphere(x)
  r = norm(x)
  lambda = atan(x[2], x[1])
  if lambda < 0
    lambda += 2 * pi
  end
  phi = asin(x[3] / r)

  return lambda, phi, r
end

function Trixi.get_element_variables!(element_variables, indicator::typeof(indicator_test), ::AMRCallback)
  return nothing
end

initial_condition = initial_condition_convergence_test_sphere

boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict(
  :inside  => boundary_condition,
  :outside => boundary_condition
)

surface_flux = flux_hll
# Note that a free stream is not preserved if N < 2 * N_geo, where N is the
# polydeg of the solver and N_geo is the polydeg of the mesh.
# However, the FSP error is negligible in this example.
solver = DGSEM(polydeg=4, surface_flux=surface_flux)

# For performance reasons, only one face of the cubed sphere can be used:

# One face of the cubed sphere
# mapping(xi, eta, zeta) = Trixi.cubed_sphere_mapping(xi, eta, zeta, 6.371229e6, 30000.0, 1)

# trees_per_dimension = (8, 8, 4)
# mesh = P4estMesh(trees_per_dimension, polydeg=4,
#                  mapping=mapping,
#                  initial_refinement_level=0,
#                  periodicity=false)

trees_per_cube_face = (6, 2)
mesh = Trixi.P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0,
                                  polydeg=4, initial_refinement_level=0)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_convergence_test_sphere,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1e5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_controller = ControllerThreeLevel(semi, indicator_test,
                                      base_level=0,
                                      max_level=1, max_threshold=0.6)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=0, # Only initial refinement
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_solution,
                        amr_callback)

###############################################################################
# run the simulation

# Use a Runge-Kutta method with automatic (error based) time step size control
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary

ext: jl


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

mutable struct ContainerServiceServicePrincipalProfile <: SwaggerModel
    clientId::Any # spec type: Union{ Nothing, String } # spec name: clientId
    secret::Any # spec type: Union{ Nothing, String } # spec name: secret

    function ContainerServiceServicePrincipalProfile(;clientId=nothing, secret=nothing)
        o = new()
        validate_property(ContainerServiceServicePrincipalProfile, Symbol("clientId"), clientId)
        setfield!(o, Symbol("clientId"), clientId)
        validate_property(ContainerServiceServicePrincipalProfile, Symbol("secret"), secret)
        setfield!(o, Symbol("secret"), secret)
        o
    end
end # type ContainerServiceServicePrincipalProfile

const _property_map_ContainerServiceServicePrincipalProfile = Dict{Symbol,Symbol}(Symbol("clientId")=>Symbol("clientId"), Symbol("secret")=>Symbol("secret"))
const _property_types_ContainerServiceServicePrincipalProfile = Dict{Symbol,String}(Symbol("clientId")=>"String", Symbol("secret")=>"String")
Base.propertynames(::Type{ ContainerServiceServicePrincipalProfile }) = collect(keys(_property_map_ContainerServiceServicePrincipalProfile))
Swagger.property_type(::Type{ ContainerServiceServicePrincipalProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ContainerServiceServicePrincipalProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceServicePrincipalProfile }, property_name::Symbol) =  _property_map_ContainerServiceServicePrincipalProfile[property_name]

function check_required(o::ContainerServiceServicePrincipalProfile)
    (getproperty(o, Symbol("clientId")) === nothing) && (return false)
    (getproperty(o, Symbol("secret")) === nothing) && (return false)
    true
end

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


text:
@testset "test ac polar pf" begin
    @testset "3-bus case" begin
        result = run_ac_pf("../test/data/case3.m", ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
        @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.10000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92617; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.90000; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["qf"], -0.403045; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["qt"],  0.0647562; atol = 1e-5)
    end
    @testset "5-bus transformer swap case" begin
        result = run_pf("../test/data/case5.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "5-bus asymmetric case" begin
        result = run_pf("../test/data/case5_asym.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "5-bus case with hvdc line" begin
        result = run_ac_pf("../test/data/case5_dc.m", ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["3"]["pg"], 3.336866; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.0635; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 1.0808; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 1.1; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.0641; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["va"], -0; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["5"]["vm"], 1.0530; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.089; atol = 1e-5)

    end
    @testset "6-bus case" begin
        result = run_pf("../test/data/case6.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.10000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.10000; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["va"], 0.00000; atol = 1e-3)
    end
    @testset "24-bus rts case" begin
        result = run_pf("../test/data/case24.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
end

@testset "test ac rect pf" begin
    @testset "5-bus asymmetric case" begin
        result = run_pf("../test/data/case5_asym.m", ACRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    #=
    # numerical issues with ipopt, likely div. by zero issue in jacobian
    @testset "5-bus case with hvdc line" begin
        result = run_pf("../test/data/case5_dc.m", ACRPowerModel, ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["3"]["pg"], 3.336866; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.0635; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 1.0808; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 1.1; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.0641; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["va"], -0; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["5"]["vm"], 1.0530; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.089; atol = 1e-5)
    end
    =#
end

@testset "test ac tan pf" begin
    @testset "5-bus asymmetric case" begin
        result = run_pf("../test/data/case5_asym.m", ACTPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "5-bus case with hvdc line" begin
        result = run_pf("../test/data/case5_dc.m", ACTPowerModel, ipopt_solver, setting = Dict("output" => Dict("branch_flows" => true)))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["3"]["pg"], 3.336866; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.0635; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 1.0808; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 1.1; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.0641; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["va"], -0; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["5"]["vm"], 1.0530; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-5)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.089; atol = 1e-5)
    end
end

@testset "test dc pf" begin
    @testset "3-bus case" begin
        result = run_dc_pf("../test/data/case3.m", ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.54994; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["va"],  0.00000; atol = 1e-5)
        @test isapprox(result["solution"]["bus"]["2"]["va"],  0.09147654582; atol = 1e-5)
        @test isapprox(result["solution"]["bus"]["3"]["va"], -0.28291891895; atol = 1e-5)
    end
    @testset "5-bus asymmetric case" begin
        result = run_dc_pf("../test/data/case5_asym.m", ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "6-bus case" begin
        result = run_dc_pf("../test/data/case6.m", ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 1e-5)
        @test isapprox(result["solution"]["bus"]["4"]["va"], 0.00000; atol = 1e-5)
    end
    @testset "24-bus rts case" begin
        result = run_pf("../test/data/case24.m", DCPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
end

@testset "test soc pf" begin
    @testset "3-bus case" begin
        result = run_pf("../test/data/case3.m", SOCWRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test result["solution"]["gen"]["1"]["pg"] >= 1.480

        @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
        @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92616; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.89999; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-4)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-4)
    end
    @testset "5-bus asymmetric case" begin
        result = run_pf("../test/data/case5_asym.m", SOCWRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "6-bus case" begin
        result = run_pf("../test/data/case6.m", SOCWRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.09999; atol = 1e-3)
    end
    @testset "24-bus rts case" begin
        result = run_pf("../test/data/case24.m", SOCWRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
end

@testset "test sdp pf" begin
    @testset "3-bus case" begin
        result = run_pf("../test/data/case3.m", SDPWRMPowerModel, scs_solver)

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test result["solution"]["gen"]["1"]["pg"] >= 1.480

        @test isapprox(result["solution"]["gen"]["2"]["pg"], 1.600063; atol = 1e-3)
        @test isapprox(result["solution"]["gen"]["3"]["pg"], 0.0; atol = 1e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 0.92616; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.89999; atol = 1e-3)

        @test isapprox(result["solution"]["dcline"]["1"]["pf"],  0.10; atol = 1e-4)
        @test isapprox(result["solution"]["dcline"]["1"]["pt"], -0.10; atol = 1e-4)
    end
    @testset "5-bus asymmetric case" begin
        result = run_pf("../test/data/case5_asym.m", SDPWRMPowerModel, scs_solver)

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
    @testset "6-bus case" begin
        result = run_pf("../test/data/case6.m", SDPWRMPowerModel, scs_solver)

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.09999; atol = 1e-3)
        @test isapprox(result["solution"]["bus"]["4"]["vm"], 1.09999; atol = 1e-3)
    end
    @testset "24-bus rts case" begin
        result = run_pf("../test/data/case24.m", SDPWRMPowerModel, scs_solver)

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 0; atol = 1e-2)
    end
end


text:
addprocs(7);

@everywhere import WAV
@everywhere import Speechbox

@everywhere using WAV
@everywhere using Speechbox

fpath = "/home/vidhya/work/data/SWB2_part/seg_wav_files";
wav_flist = filter(x->endswith(x,".wav"),readdir(fpath));

@everywhere function ex_mfcc(fname)
  ext = ".mfc"
  fpath = "/home/vidhya/work/data/SWB2_part/seg_wav_files";
  opath = "/home/vidhya/work/workspace/swb_mfcc";
  bname,~ = splitext(fname);
  ifile = joinpath(fpath,fname);
  ofile = joinpath(opath,bname*ext);
  #println(ifile);

  x,fs = wavread(ifile);
  sx = normalise_energy(x,fs);
  mf = melfcc(sx,fs);
  write_binary(ofile,mf);
end

println(nprocs())

function run_ex_mfcc(wav_flist)
  @sync @parallel for i=1:10000
    ex_mfcc(wav_flist[i]);
  end
  return nothing
end

#@time pmap(ex_mfcc,wav_flist[1:1000])
@time run_ex_mfcc(wav_flist);

rmprocs(workers());
println(nprocs())


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG
export AbstractRNGSeq, ZRNGSeq, next!, nextn!, skipTo

@inline rand(r::AbstractRNG{UInt64}, ::Type{Float64}) = (Float64(Random.rand(r, UInt64)  >> 12 )+0.5) / 0x1p52
#@inline rand(r::AbstractRNG{UInt64}, ::Type{Float64}) = ldexp(Float64(Random.rand(r, UInt64)  >> 11 )+0.5, ~52)

@inline rand(r::AbstractRNG{UInt32}, ::Type{Float64}) = (Float64(Random.rand(r, UInt32) )+0.5) /0x1p32


text:
@testset "igmrf.jl" begin
    @testset "Base.show(io, obj)" begin
        # does not throw
        G = GMRF.GridStructure((1,2), [[2], [1]], [[1],[2]], spzeros(1, 2), spzeros(1, 2))
        igmrf = iGMRF(G, 1, 1.0)
        io = IOBuffer()
        @test_logs Base.show(io, igmrf)

    end

    @testset "iGMRF(m₁, m₂, order, κ)" begin
        # order not 1 or 2 throws
        @test_throws AssertionError iGMRF(1, 1, 0, 1.0)

        # Simple grid of order 1
        igmrf = iGMRF(1, 1, 1, 1.0)

        @test igmrf.G.gridSize == (1, 1)
        @test igmrf.rankDeficiency == 1
        @test igmrf.κ ≈ 1.0
        # TODO : Test that W̄ was derived from W

        # Simple grid of order 2
        igmrf = iGMRF(1, 1, 2, 1.0)

        @test igmrf.G.gridSize == (1, 1)
        @test igmrf.rankDeficiency == 3
        @test igmrf.κ ≈ 1.0
        # TODO : Test that W̄ was derived from W

    end

    @testset "fo_nbs(m₁, m₂)" begin
        # Grid 1 x 1
        nbs, W = GMRF.fo_nbs(1, 1)

        @test nbs == [[]]
        # TODO : Test W

        # Grid 2 x 2
        nbs, W = GMRF.fo_nbs(2, 2)

        @test nbs == [[2, 3], [1, 4], [1, 4], [2, 3]]
        # TODO : Test W

    end

    @testset "so_nbs(m₁, m₂)" begin
        # Grid 1 x 1
        nbs, W = GMRF.so_nbs(1, 1)

        @test nbs == [[]]
        # TODO : Test W

        # Grid 2 x 2
        nbs, W = GMRF.so_nbs(2, 2)

        # TODO : Test nbs
        # TODO : Test W

    end

    @testset "fo_condindsubsets(m₁, m₂)" begin
        # Grid 1 x 1
        cond = GMRF.fo_condindsubsets(1, 1)

        # TODO : Test cond

        # Grid 2 x 2
        cond = GMRF.fo_condindsubsets(2, 2)

        # TODO : Test cond

    end

    @testset "so_condindsubsets(m₁, m₂)" begin
        # Grid 1 x 1
        cond = GMRF.so_condindsubsets(1, 1)

        # TODO : Test cond

        # Grid 2 x 2
        cond = GMRF.so_condindsubsets(2, 2)

        # TODO : Test cond

    end

    @testset "rand(F)" begin
        # rankdeficiency != 1 or 3 throws
        G = GMRF.GridStructure((1,2), [[2], [1]], [[1],[2]], spzeros(1, 2), spzeros(1, 2))
        igmrf = iGMRF(G, 2, 1.0)
        @test_throws AssertionError rand(igmrf)

        # returns plausible data (rankdeficiency == 1)
        # TODO: Test if all insupport ?

        # returns plausible data (rankdeficiency == 3)
        # TODO: Test if all insupport ?

    end

    @testset "logpdf(F, y)" begin
        # TODO : Test with known values

    end

    @testset "fullconditionals(F, y)" begin
        # TODO : Test with known values

    end

    @testset "fullcondlogpdf(F, y)" begin
        # TODO : Test with known values

    end

    @testset "getconditional(F, B, x)" begin
        # TODO : Test with known values

    end

end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    Optimizer(; <keyword arguments>)

Return a quasi-gradient optimizer.
...
# Arguments
- `master_optimizer::AbstractOptimizer`: MathOptInterface solver capable of solving linear (and possibly quadratic) programs.
- `subproblem_optimizer::AbstractOptimizer`: Optionally specify a different solver for the subproblems.
- `execution::Execution = Serial`: Specify how algorithm should be executed (Serial, Synchronous, Asynchronous). Distributed variants requires worker cores.
- `subproblems::AbstractSubProblemState = Unaltered()`: Specify if a smoothing procedure should be applied.
- `prox::AbstractProx = Polyhedron()`: Specify proximal step.
- `step::AbstractStep = Constant()`: Specify step-size
- `termination::AbstractTermination = AfterMaximumIterations()`: Specify termination criterion

- <keyword arguments>: Algorithm specific parameters, See `?LShaped` for list of possible arguments and default values.
...
"""
mutable struct Optimizer <: AbstractStructuredOptimizer
    master_optimizer
    subproblem_optimizer
    master_params::Dict{MOI.AbstractOptimizerAttribute, Any}
    sub_params::Dict{MOI.AbstractOptimizerAttribute, Any}
    execution::AbstractExecution
    subproblems::AbstractSubProblemState
    prox::AbstractProx
    step::AbstractStepSize
    termination::AbstractTermination
    parameters::QuasiGradientParameters{Float64}

    status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    dual_status::MOI.ResultStatusCode
    raw_status::String
    solve_time::Float64

    quasigradient::Union{AbstractQuasiGradient, Nothing}

    function Optimizer(; master_optimizer = nothing,
                       execution::AbstractExecution = nworkers() == 1 ? Serial() : Synchronous(),
                       subproblems::AbstractSubProblemState = Unaltered(),
                       prox::AbstractProx = Polyhedron(),
                       step::AbstractStepSize = Constant(),
                       terminate::AbstractTermination = AfterMaximumIterations(),
                       subproblem_optimizer = nothing, kw...)
        return new(master_optimizer,
                   subproblem_optimizer,
                   Dict{MOI.AbstractOptimizerAttribute, Any}(),
                   Dict{MOI.AbstractOptimizerAttribute, Any}(),
                   execution,
                   subproblems,
                   prox,
                   step,
                   terminate,
                   QuasiGradientParameters{Float64}(; kw...),
                   MOI.OPTIMIZE_NOT_CALLED,
                   MOI.NO_SOLUTION,
                   MOI.NO_SOLUTION,
                   "Quasi-gradient optimizer has not been run.",
                   NaN,
                   nothing)
    end
end

# Interface #
# ========================== #
function supports_structure(optimizer::Optimizer, ::StageDecompositionStructure)
    return true
end

function default_structure(::UnspecifiedInstantiation, optimizer::Optimizer)
    if optimizer.execution isa Serial && nworkers() == 1
        return StageDecomposition()
    else
        return DistributedStageDecomposition()
    end
end

function check_loadable(optimizer::Optimizer, ::StageDecompositionStructure)
    if optimizer.prox isa Polyhedron && optimizer.master_optimizer === nothing
        msg = "Polyhedron proximal operator requires setting `MasterOptimizer` attribute."
        throw(UnloadableStructure{Optimizer, StageDecompositionStructure}(msg))
    end
    return nothing
end

function ensure_compatible_execution!(optimizer::Optimizer, ::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}})
    if !(optimizer.execution isa Serial)
        @warn "Distributed execution policies are not compatible with a single-core stage-decomposition structure. Switching to `Serial` execution by default."
        MOI.set(optimizer, Execution(), Serial())
    end
    return nothing
end

function ensure_compatible_execution!(optimizer::Optimizer, ::StageDecompositionStructure{2, 1, <:Tuple{DistributedScenarioProblems}})
    if optimizer.execution isa Serial
        @warn "Serial execution not compatible with distributed stage-decomposition structure. Switching to `Synchronous` execution by default."
        MOI.set(optimizer, Execution(), Synchronous())
    end
    return nothing
end

function load_structure!(optimizer::Optimizer, structure::StageDecompositionStructure, x₀::AbstractVector)
    # Sanity check
    check_loadable(optimizer, structure)
    # Restore structure if optimization has been run before
    restore_structure!(optimizer)
    # Ensure that execution policy is compatible
    ensure_compatible_execution!(optimizer, structure)
    # Set master optimizer
    set_master_optimizer!(structure, optimizer.master_optimizer)
    # Set subproblem optimizers
    set_subproblem_optimizer!(structure, optimizer.subproblem_optimizer)
    # Set any given prox/sub optimizer attributes
    for (attr, value) in optimizer.master_params
        MOI.set(backend(structure.first_stage), attr, value)
    end
    for (attr, value) in optimizer.sub_params
        MOI.set(scenarioproblems(structure), attr, value)
    end
    # Create new L-shaped algorithm
    optimizer.quasigradient = QuasiGradientAlgorithm(structure,
                                                     x₀,
                                                     optimizer.execution,
                                                     optimizer.subproblems,
                                                     optimizer.prox,
                                                     optimizer.step,
                                                     optimizer.termination;
                                                     type2dict(optimizer.parameters)...)
    return nothing
end

function restore_structure!(optimizer::Optimizer)
    if optimizer.quasigradient !== nothing
        restore_master!(optimizer.quasigradient)
        restore_subproblems!(optimizer.quasigradient)
    end
    return nothing
end

function MOI.optimize!(optimizer::Optimizer)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    start_time = time()
    optimizer.status = optimizer.quasigradient()
    if optimizer.status == MOI.OPTIMAL
        optimizer.primal_status = MOI.FEASIBLE_POINT
        optimizer.dual_status = MOI.FEASIBLE_POINT
        optimizer.raw_status = "Quasi-gradient procedure converged to optimal solution."
    end
    optimizer.solve_time = time() - start_time
    return nothing
end

function num_iterations(optimizer::Optimizer)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return num_iterations(optimizer.quasigradient)
end

function optimizer_name(optimizer::Optimizer)
    return "$(str(optimizer.execution))Quasi-gradient"
end

# MOI #
# ========================== #
function MOI.get(optimizer::Optimizer, ::MOI.Silent)
    return !MOI.get(optimizer, MOI.RawParameter("log"))
end

function MOI.set(optimizer::Optimizer, attr::MOI.Silent, flag::Bool)
    MOI.set(optimizer, MOI.RawParameter("log"), !flag)
    optimizer.master_params[attr] = flag
    optimizer.sub_params[attr] = flag
    return nothing
end

function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(QuasiGradientParameters))
        error("Unrecognized parameter name: $(name).")
    end
    return getfield(optimizer.parameters, name)
end

function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(QuasiGradientParameters))
        error("Unrecognized parameter name: $(name).")
    end
    setfield!(optimizer.parameters, name, value)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
    limit = MOI.get(optimizer, MOI.RawParameter("time_limit"))
    return isinf(limit) ? nothing : limit
end

function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, limit::Union{Real, Nothing})
    limit = limit === nothing ? Inf : limit
    MOI.set(optimizer, MOI.RawParameter("time_limit"), limit)
    return
end

function MOI.get(optimizer::Optimizer, ::RelativeTolerance)
    return MOI.get(optimizer, MOI.RawParameter("τ"))
end

function MOI.set(optimizer::Optimizer, ::RelativeTolerance, limit::Real)
    MOI.set(optimizer, MOI.RawParameter("τ"), limit)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::MasterOptimizer)
    if optimizer.master_optimizer === nothing
        return nothing
    end
    return MOI.OptimizerWithAttributes(optimizer.master_optimizer, collect(optimizer.master_params))
end

function MOI.set(optimizer::Optimizer, ::MasterOptimizer, optimizer_constructor)
    optimizer.master_optimizer = optimizer_constructor
    return nothing
end

function MOI.get(optimizer::Optimizer, ::MasterOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute)
    if !haskey(optimizer.master_params, attr)
        error("Master optimizer attribute $(attr) has not been set.")
    end
    return optimizer.master_params[attr]
end

function MOI.set(optimizer::Optimizer, ::MasterOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute, value)
    optimizer.master_params[attr] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, param::RawMasterOptimizerParameter)
    moi_param = MOI.RawParameter(param.name)
    if !haskey(optimizer.master_params, moi_param)
        error("Master optimizer attribute $(param.name) has not been set.")
    end
    return optimizer.master_params[moi_param]
end

function MOI.set(optimizer::Optimizer, param::RawMasterOptimizerParameter, value)
    moi_param = MOI.RawParameter(param.name)
    optimizer.master_params[moi_param] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, ::SubProblemOptimizer)
    if optimizer.subproblem_optimizer === nothing
        return MOI.get(optimizer, MasterOptimizer())
    end
    return MOI.OptimizerWithAttributes(optimizer.subproblem_optimizer, collect(optimizer.sub_params))
end

function MOI.set(optimizer::Optimizer, ::SubProblemOptimizer, optimizer_constructor)
    optimizer.subproblem_optimizer = optimizer_constructor
    return nothing
end

function MOI.get(optimizer::Optimizer, ::SubProblemOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute)
    if !haskey(optimizer.sub_params, attr)
        error("Subproblem optimizer attribute $(attr) has not been set.")
    end
    return optimizer.sub_params[attr]
end

function MOI.set(optimizer::Optimizer, ::SubProblemOptimizerAttribute, attr::MOI.AbstractOptimizerAttribute, value)
    optimizer.sub_params[attr] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, param::RawSubProblemOptimizerParameter)
    moi_param = MOI.RawParameter(param.name)
    if !haskey(optimizer.sub_params, moi_param)
        error("Subproblem optimizer attribute $(param.name) has not been set.")
    end
    return optimizer.sub_params[moi_param]
end

function MOI.set(optimizer::Optimizer, param::RawSubProblemOptimizerParameter, value)
    moi_param = MOI.RawParameter(param.name)
    optimizer.sub_params[moi_param] = value
    return nothing
end

function MOI.get(optimizer::Optimizer, ::Execution)
    return optimizer.execution
end

function MOI.set(optimizer::Optimizer, ::Execution, execution::AbstractExecution)
    optimizer.execution = execution
    return nothing
end

function MOI.get(optimizer::Optimizer, ::SubProblems)
    return optimizer.subproblems
end

function MOI.set(optimizer::Optimizer, ::SubProblems, subproblems::AbstractSubProblemState)
    optimizer.subproblems = subproblems
    return nothing
end

function MOI.get(optimizer::Optimizer, ::Prox)
    return optimizer.prox
end

function MOI.set(optimizer::Optimizer, ::Prox, prox::AbstractProx)
    optimizer.prox = prox
    return nothing
end

function MOI.get(optimizer::Optimizer, ::StepSize)
    return optimizer.step
end

function MOI.set(optimizer::Optimizer, ::StepSize, step::AbstractStepSize)
    optimizer.step = step
    return nothing
end

function MOI.get(optimizer::Optimizer, ::Termination)
    return optimizer.termination
end

function MOI.set(optimizer::Optimizer, ::Termination, termination::AbstractTermination)
    optimizer.termination = termination
    return nothing
end

function MOI.get(optimizer::Optimizer, param::ExecutionParameter)
    return MOI.get(optimizer.execution, param)
end

function MOI.set(optimizer::Optimizer, param::ExecutionParameter, value)
    MOI.set(optimizer.execution, param, value)
    return nothing
end

function MOI.get(optimizer::Optimizer, param::ProxParameter)
    return MOI.get(optimizer.prox, param)
end

function MOI.set(optimizer::Optimizer, param::ProxParameter, value)
    MOI.set(optimizer.prox, param, value)
    return nothing
end

function MOI.get(optimizer::Optimizer, param::StepParameter)
    return MOI.get(optimizer.step, param)
end

function MOI.set(optimizer::Optimizer, param::StepParameter, value)
    MOI.set(optimizer.step, param, value)
    return nothing
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    return optimizer.status
end

function MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus)
    return optimizer.primal_status
end

function MOI.get(optimizer::Optimizer, ::MOI.DualStatus)
    return optimizer.dual_status
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return optimizer.raw_status
end

function MOI.get(optimizer::Optimizer, attr::MOI.ListOfVariableIndices)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    list = MOI.get(optimizer.quasigradient.structure, attr)
    # Remove any auxilliary variables from prox
    filter_variables!(optimizer.quasigradient.prox, list)
    return list
end

function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, index::MOI.VariableIndex)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return decision(optimizer.quasigradient, index)
end

function MOI.get(optimizer::Optimizer, attr::MOI.ListOfConstraintIndices)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    list = MOI.get(optimizer.quasigradient.structure, attr)
    # Remove any prox constraints
    filter_constraints!(optimizer.quasigradient.prox, list)
    return list
end

function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue)
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return objective_value(optimizer.quasigradient)
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
    return optimizer.solve_time
end

function MOI.get(optimizer::Optimizer, attr::Union{MOI.AbstractOptimizerAttribute, MOI.AbstractModelAttribute})
    # Fallback to first-stage optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.get(optimizer.quasigradient.structure, attr)
end

function MOI.get(optimizer::Optimizer, attr::MOI.AbstractConstraintAttribute, ci::MOI.ConstraintIndex)
    # Fallback to first-stage optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.get(optimizer.quasigradient.structure, attr, ci)
end

function MOI.set(optimizer::Optimizer, attr::Union{MOI.AbstractOptimizerAttribute, MOI.AbstractModelAttribute}, value)
    # Fallback to first-stage optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, value)
end

function MOI.set(optimizer::Optimizer, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex, value)
    # Fallback to first-stage optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, index, value)
end

function MOI.set(optimizer::Optimizer, attr::MOI.AbstractConstraintAttribute, ci::MOI.ConstraintIndex, value)
    # Fallback to first-stage optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, ci, value)
end

function MOI.get(optimizer::Optimizer, attr::ScenarioDependentModelAttribute)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.get(optimizer.quasigradient.structure, attr)
end

function MOI.get(optimizer::Optimizer, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.get(optimizer.quasigradient.structure, attr, index)
end

function MOI.get(optimizer::Optimizer, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.get(optimizer.quasigradient.structure, attr, ci)
end

function MOI.set(optimizer::Optimizer, attr::ScenarioDependentModelAttribute, value)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, value)
end

function MOI.set(optimizer::Optimizer, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex, value)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, index, value)
end

function MOI.set(optimizer::Optimizer, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex, value)
    # Fallback to subproblem optimizer through structure
    if optimizer.quasigradient === nothing
        throw(UnloadedStructure{Optimizer}())
    end
    return MOI.set(optimizer.quasigradient.structure, attr, ci, value)
end

function MOI.is_empty(optimizer::Optimizer)
    return optimizer.quasigradient === nothing
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.RawParameter) = true
MOI.supports(::Optimizer, ::AbstractStructuredOptimizerAttribute) = true
MOI.supports(::Optimizer, ::RawInstanceOptimizerParameter) = true
MOI.supports(::Optimizer, ::AbstractQuasiGradientAttribute) = true

# High-level attribute setting #
# ========================== #
"""
    get_prox_attribute(stochasticprogram::StochasticProgram, name::String)

Return the value associated with the prox-specific attribute named `name` in `stochasticprogram`.

See also: [`set_prox_attribute`](@ref), [`set_prox_attributes`](@ref).
"""
function get_prox_attribute(stochasticprogram::StochasticProgram, name::String)
    return return MOI.get(optimizer(stochasticprogram), RawProxParameter(name))
end
"""
    set_prox_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)

Sets the prox-specific attribute identified by `name` to `value`.

"""
function set_prox_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticprogram, RawProxParameter(String(name)), value)
end
"""
    set_prox_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)

Given a list of `attribute => value` pairs or a collection of keyword arguments, calls
`set_prox_attribute(stochasticprogram, attribute, value)` for each pair.

"""
function set_prox_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)
    for (name, value) in pairs
        set_prox_attribute(stochasticprogram, name, value)
    end
end
function set_prox_attributes(stochasticprogram::StochasticProgram; kw...)
    for (name, value) in kw
        set_prox_attribute(stochasticprogram, name, value)
    end
end
"""
    get_step_attribute(stochasticprogram::StochasticProgram, name::String)

Return the value associated with the step-specific attribute named `name` in `stochasticprogram`.

See also: [`set_step_attribute`](@ref), [`set_step_attributes`](@ref).
"""
function get_step_attribute(stochasticprogram::StochasticProgram, name::String)
    return MOI.get(optimizer(stochasticprogram), RawStepParameter(name))
end
"""
    set_step_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)

Sets the step-specific attribute identified by `name` to `value`.

"""
function set_step_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticprogram, RawStepParameter(String(name)), value)
end
"""
    set_step_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)

Given a list of `attribute => value` pairs or a collection of keyword arguments, calls
`set_step_attribute(stochasticprogram, attribute, value)` for each pair.

"""
function set_step_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)
    for (name, value) in pairs
        set_step_attribute(stochasticprogram, name, value)
    end
end
function set_step_attributes(stochasticprogram::StochasticProgram; kw...)
    for (name, value) in kw
        set_step_attribute(stochasticprogram, name, value)
    end
end
"""
    get_termination_attribute(stochasticprogram::StochasticProgram, name::String)

Return the value associated with the termination-specific attribute named `name` in `stochasticprogram`.

See also: [`set_termination_attribute`](@ref), [`set_termination_attributes`](@ref).
"""
function get_termination_attribute(stochasticprogram::StochasticProgram, name::String)
    return MOI.get(optimizer(stochasticprogram), RawTerminationParameter(name))
end
"""
    set_termination_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)

Sets the termination-specific attribute identified by `name` to `value`.

"""
function set_termination_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticprogram, RawTerminationParameter(String(name)), value)
end
"""
    set_termination_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)

Given a list of `attribute => value` pairs or a collection of keyword arguments, calls
`set_termination_attribute(stochasticprogram, attribute, value)` for each pair.

"""
function set_termination_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)
    for (name, value) in pairs
        set_termination_attribute(stochasticprogram, name, value)
    end
end
function set_termination_attributes(stochasticprogram::StochasticProgram; kw...)
    for (name, value) in kw
        set_termination_attribute(stochasticprogram, name, value)
    end
end


text:
using StatPlots, RDatasets
iris = dataset("datasets", "iris")
@df iris scatter(:SepalLength, :SepalWidth, group=:Species,
    m=(0.5, [:+ :h :star7], 4), bg=RGB(1.0,1.0,1.0))


text:
using LinearAlgebra

function reverse_gauss_1(A,b)
    x = similar(b)
    N = size(A, 1)
    for k in 0:N-1
        x[N-k] = (b[N-k] - sum(@view(A[N-k,N-k+1:N]) .* @view(x[N-k+1:N])))/A[N-k,N-k]
    end
    return x
end

function reverse_gauss(Ab)
#Ab - транспонированная расширенная матрица
    N = size(Ab, 2)
    x = Vector{eltype(Ab)}(undef,N)

    for k in 0:N-1
        x[N-k] = (Ab[N-k, N] - dot(@view(Ab[N-k+1:N, N-k]), @view(x[N-k+1:N])))/Ab[N-k,N-k]
    end
    return x
end

function upper_triangular(N::Integer)
    A = randn(N, N)
    _, A = lu(A)
    return A
end

function upper_triangular_bad(N::Integer)
    A = zeros(N,N)
    for i in 1:N
        A[i,i+1:N] = randn(N-i)
        A[i,i] = 1.0
    end
    return A
end

######
N=10_000
A = upper_triangular(N)
b = ones(N)
Abt=[A b]

println("1:")
@time reverse_gauss_1(A,b);

println("2:")
@time reverse_gauss(Abt);
println("...................")
########

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

function add_1(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A) == size(B)
    C = similar(A)
    for i in 1:size(A,1), j in 1:size(A,2) #
        C[i,j] = A[i,j] + B[i,j]
    end
    #= это эквивалентно:
    for i in 1:size(A,1)
        for j in 1:size(A,2)
            C[i,j] = A[i,j] + B[i,j]
        end
    end
    =#
    return C
end

function add_2(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A) == size(B)
    C = similar(A)
    for  j in 1:size(A,2), i in 1:size(A,1) # - это более педпочтиельный вариант (преимущество проявляется на больших матрицах)
        C[i,j] = A[i,j] + B[i,j]
    end
    return C
end

using LinearAlgebra

function mul(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A,2) == size(B, 1)
    C = Matrix{T}(undef, size(A,1), size(B,2))
    At=A'
    for j in 1:size(C,2), i in 1:size(C,1)
        p=0.0
        for k in 1:size(A,2) # size(A,2) == size(B, 1)
            p += At[k,i]*B[k,j]
        end
        C[i,j] = p
    end
    return C
end

function mul_1(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A,2) == size(B, 1)
    C = Matrix{T}(undef, size(A,1), size(B,2))
    for i in 1:size(C,1), j in 1:size(C,2)
        @views C[i,j] = dot(A[i,:],B[:,j])
    end
    #= это эквивалентно:
    for i in 1:size(C,1), j in 1:size(C,2)
        C[i,j] = T(0)
        for k in 1:size(A,2) # size(A,2) == size(B, 1)
            C[i,j] += A[i,k]*conj(B[k,j])
        end
    end
    =#
    return C
end

function mul_2(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A,2) == size(B, 1)
    C = Matrix{T}(undef, size(A,1), size(B,2))
    for j in 1:size(C,2), i in 1:size(C,1)
        @views C[i,j] = dot(A[i,:], B[:,j])
    end
    return C
end

function mul_3(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T<:Number
    @assert size(A,2) == size(B, 1)
    C = Matrix{T}(undef, size(A,1), size(B,2))
    At = A'
    for j in 1:size(C,2), i in 1:size(C,1)
        @views C[i,j] = dot(At[:,i], B[:,j])
    end
    return C
end
#----------------------------------------------------
print("N = ")
N = parse(Int, readline())
#=
При N=10_000 выигрыш примерно на порядок
при N=1000 он всего раза в два
=#
A = randn(N,N)
B = randn(N,N)
println("Варианты реализации сложения матриц:")
@time add_1(A,B);
@time add_2(A,B);
println("Варианты реализации перемножения матриц:")
@time mul(A,B);
#=
@time mul_1(A,B);
@time mul_2(A,B);
@time mul_3(A,B);
=#
println()

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using SafeTestsets
@safetestset "Dummy testing" begin include("dummy_tests.jl") end


text:
abstract type BeliefNode{B,A,O} end
abstract type AbstractActNode end

# Note: links to parents were taken out because they hadn't been used in anything we've done so far
# Note: probably don't need the labels, but they don't seem like they would really kill performance

type ActNode{A, O, BNodeType <: BeliefNode} <: AbstractActNode
    label::A # for keeping track of which action this corresponds to
    N::Int
    V::Float64
    children::Dict{O, BNodeType} # maps observations to ObsNodes
end

type ObsNode{Belief,A,O} <: BeliefNode{Belief,A,O}
    label::O
    N::Int # for dpw, this is the number of times we have transitioned from parent to this from the parent
    B::Belief # belief/state distribution
    children::Dict{A,ActNode{A,O,ObsNode{Belief,A,O}}}
end

type RootNode{RootBelief,A,ANodeType} <: BeliefNode{RootBelief}
    N::Int
    B::RootBelief # belief/state distribution
    children::Dict{A,ANodeType} # ActNode not parameterized here to make initialize_belief more flexible
end
RootNode{RootBelief}(b::RootBelief) = RootNode{RootBelief,Any,ActNode}(0, b, Dict{Any,ActNode}())

"""
    init_V(initializer, problem::POMDPs.POMDP, h::BeliefNode, action)

Return the initial value (V) associated with a new action node when it is created. This can be used in concert with `init_N` to incorporate prior experience into the solver.
"""
function init_V end
init_V(n::Number, problem::POMDPs.POMDP, h::BeliefNode, action) = convert(Float64, n)
init_V(f::Function, problem::POMDPs.POMDP, h::BeliefNode, action) = f(problem, h, action)

"""
    init_N(initializer, problem::POMDPs.POMDP, h::BeliefNode, action)

Return the initial number of queries (N) associated with a new action node when it is created.
"""
function init_N end
init_N(n::Number, problem::POMDPs.POMDP, h::BeliefNode, action) = convert(Int, n)
init_N(f::Function, problem::POMDPs.POMDP, h::BeliefNode, action) = f(problem, h, action)


text:
# :PyFloat_FromDouble => (Cdouble,) => PyPtr,
# :PyFloat_AsDouble => (PyPtr,) => Cdouble,

"""
    pyfloat(x=0.0)

Convert `x` to a Python `float`.
"""
pyfloat(x::Real=0.0) = pynew(errcheck(C.PyFloat_FromDouble(x)))
pyfloat(x) = @autopy x pynew(errcheck(C.PyNumber_Float(getptr(x_))))
export pyfloat

pyisfloat(x) = pytypecheck(x, pybuiltins.float)

pyfloat_asdouble(x) = errcheck_ambig(@autopy x C.PyFloat_AsDouble(getptr(x_)))

function pyconvert_rule_float(::Type{T}, x::Py) where {T<:Number}
    val = pyfloat_asdouble(x)
    if T in (Float16, Float32, Float64, BigFloat)
        pyconvert_return(T(val))
    else
        pyconvert_tryconvert(T, val)
    end
end

# NaN is sometimes used to represent missing data of other types
# so we allow converting it to Nothing or Missing
function pyconvert_rule_float(::Type{Nothing}, x::Py)
    if isnan(pyfloat_asdouble(x))
        pyconvert_return(nothing)
    else
        pyconvert_unconverted()
    end
end

function pyconvert_rule_float(::Type{Missing}, x::Py)
    if isnan(pyfloat_asdouble(x))
        pyconvert_return(missing)
    else
        pyconvert_unconverted()
    end
end

pyconvert_rule_fast(::Type{Float64}, x::Py) =
    if pyisfloat(x)
        pyconvert_return(pyfloat_asdouble(x))
    else
        pyconvert_unconverted()
    end


text:
module CytofResearch
#= Note:
Julia uses Gamma(shape, scale) and InverseGamma(shape, scale).
Note that InverseGamma(shape, scale) IS my preferred parameterization.
It has a mean of scale / (shape - 1) for shape > 1.
=#

export MCMC

include("Util/Util.jl")
include("MCMC/MCMC.jl")
include("Model/Model.jl")
include("VB/VB.jl")

end # module

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
@doc raw"""
    GeneralizedBurresWassertseinMetric{T<:AbstractMatrix} <: AbstractMetric

The generalized Bures Wasserstein metric for symmetric positive definite matrices, see[^HanMishraJawanpuriaGao2021].

This metric internally stores the symmetric positive definite matrix ``M`` to generalise the metric,
where the name also follows the mentioned preprint.

[^HanMishraJawanpuriaGao2021]:
    > Han, A., Mishra, B., Jawanpuria, P., Gao, J.:
    > _Generalized Bures-Wasserstein geometry for positive definite matrices_.
    > arXiv: [2110.10464](https://arxiv.org/abs/2110.10464).
"""
struct GeneralizedBuresWassersteinMetric{T<:AbstractMatrix} <: RiemannianMetric
    M::T
    GeneralizedBuresWassersteinMetric(MM::TT) where {TT<:AbstractMatrix} = new{TT}(MM)
end

@doc raw"""
    change_representer(M::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, E::EuclideanMetric, p, X)

Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function on the tangent
space at `p` with respect to the [`EuclideanMetric`](@ref) `g_E`,
this is turned into the representer with respect to the (default) metric,
the [`GeneralizedBuresWassersteinMetric`](@ref) on the [`SymmetricPositiveDefinite`](@ref) `M`.

To be precise we are looking for ``Z∈T_p\mathcal P(n)`` such that for all ``Y∈T_p\mathcal P(n)``
it holds

```math
⟨X,Y⟩ = \operatorname{tr}(XY) = ⟨Z,Y⟩_{\mathrm{BW}}
```
for all ``Y`` and hence we get ``Z = 2pXM + 2MXp``.
"""
change_representer(
    ::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric},
    ::EuclideanMetric,
    p,
    X,
)

function change_representer!(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric},
    Y,
    ::EuclideanMetric,
    p,
    X,
)
    Y .= 2 .* (p * X * M.metric.M + M.metric.M * X * p)
    return Y
end

@doc raw"""
    distance(::MatricManifold{SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q)

Compute the distance with respect to the [`BuresWassersteinMetric`](@ref) on [`SymmetricPositiveDefinite`](@ref) matrices, i.e.

```math
d(p,q) = \operatorname{tr}(M^{-1}p) + \operatorname{tr}(M^{-1}q)
       - 2\operatorname{tr}\bigl( (p^{\frac{1}{2}}M^{-1}qM^{-1}p^{\frac{1}{2}} \bigr)^{\frac{1}{2}},
```
"""
function distance(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric},
    p,
    q,
)
    luM = lu(M.metric.M)
    luMp = luM \ p
    luMq = luM \ q
    return sqrt(tr(luMp) + tr(luMq) - 2 * tr(sqrt(luMq * luMp)))
end

@doc raw"""
    exp(::MatricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X)

Compute the exponential map on [`SymmetricPositiveDefinite`](@ref) with respect to
the [`GeneralizedBuresWassersteinMetric`](@ref) given by

```math
    \exp_p(X) = p+X+\mathcal ML_{p,M}(X)pML_{p,M}(X)
```

where ``q=L_{M,p}(X)`` denotes the generalized Lyapunov operator, i.e. it solves ``pqM + Mqp = X``.
"""
exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X)

function exp!(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric},
    q,
    p,
    X,
)
    m = M.metric.M
    Y = lyapc(p, m, -X) #lyap solves qpM + Mpq - X =0
    q .= p .+ X .+ m * Y * p * Y * m
    return q
end

@doc raw"""
    inner(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X, Y)

Compute the inner product [`SymmetricPositiveDefinite`](@ref) with respect to
the [`GeneralizedBuresWassersteinMetric`](@ref) given by

```math
    ⟨X,Y⟩ = \frac{1}{2}\operatorname{tr}(L_{p,M}(X)Y)
```

where ``q=L_{M,p}(X)`` denotes the generalized Lyapunov operator, i.e. it solves ``pqM + Mqp = X``.
"""
function inner(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric},
    p,
    X,
    Y,
)
    return dot(lyapc(p, M.metric.M, -X), Y) / 2
end

@doc raw"""
    log(::MatricManifold{SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q)

Compute the logarithmic map on [`SymmetricPositiveDefinite`](@ref) with respect to
the [`BuresWassersteinMetric`](@ref) given by

```math
    \log_p(q) = M(M^{-1}pM^{-1}q)^{\frac{1}{2}} + (qM^{-1}pM^{-1})^{\frac{1}{2}}M - 2 p.
```
"""
log(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q)

function log!(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric},
    X,
    p,
    q,
)
    m = M.metric.M
    lum = lu(m)
    lum_p_lum = lum \ p / lum
    X .= real.(Symmetric(m * sqrt(lum_p_lum * q) + sqrt(q * lum_p_lum) * m) - 2 * p)
    return X
end


text:
module TOML

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    # All TOML errors inherit from the abstract Error
    abstract Error <: Exception

    # Errors that arise in parsing input
    type ParserError <: Error
        message::String
    end
    Base.showerror(io::IO, e::ParserError) = print(io,
        join(["TOML.ParserError: ", e.message]))

    # Errors that arise in Table structure
    type TableError <: Error
        message::String
    end
    Base.showerror(io::IO, e::TableError) = print(io,
        join(["TOML.TableError: ", e.message]))

    # Errors that arise in handling types within TOML
    type TypeConversionError <: Error
        message::String
    end
    Base.showerror(io::IO, e::TypeConversionError) = print(io,
        join(["TOML.TypeConversionError: ", e.message]))

    # Errors that arise in dealing with Values
    type ValueError <: Error
        message::String
    end
    Base.showerror(io::IO, e::ValueError) = print(io,
        join(["TOML.ValueError: ", e.message]))

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    # The parsing routines, especially for Value, return the value in the
    # desired format.  However, numbers can be either integers, floats, or
    # both; thus we use this struct for returning numbers.
    type NumberValue
        integer_value::Int64
        float_value::Float64
        valid_integer::Bool
        valid_float::Bool
    end # type NumberValue

    function NumberValue(; integer_value::Int64 = 0,
        float_value::Float64 = 0.0, valid_integer::Bool = false,
        valid_float::Bool = false)
        return(NumberValue(integer_value, float_value,
            valid_integer, valid_float))
    end

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    # Check if the character is a valid letter (a-zA-Z).
    function is_letter(c::Char)
        if c == 'a' || c == 'A'
            return(true)
        elseif c == 'b' || c == 'B'
            return(true)
        elseif c == 'c' || c == 'C'
            return(true)
        elseif c == 'd' || c == 'D'
            return(true)
        elseif c == 'e' || c == 'E'
            return(true)
        elseif c == 'f' || c == 'F'
            return(true)
        elseif c == 'g' || c == 'G'
            return(true)
        elseif c == 'h' || c == 'H'
            return(true)
        elseif c == 'i' || c == 'I'
            return(true)
        elseif c == 'j' || c == 'J'
            return(true)
        elseif c == 'k' || c == 'K'
            return(true)
        elseif c == 'l' || c == 'L'
            return(true)
        elseif c == 'm' || c == 'M'
            return(true)
        elseif c == 'n' || c == 'N'
            return(true)
        elseif c == 'o' || c == 'O'
            return(true)
        elseif c == 'p' || c == 'P'
            return(true)
        elseif c == 'q' || c == 'Q'
            return(true)
        elseif c == 'r' || c == 'R'
            return(true)
        elseif c == 's' || c == 'S'
            return(true)
        elseif c == 't' || c == 'T'
            return(true)
        elseif c == 'u' || c == 'U'
            return(true)
        elseif c == 'v' || c == 'V'
            return(true)
        elseif c == 'w' || c == 'W'
            return(true)
        elseif c == 'x' || c == 'X'
            return(true)
        elseif c == 'y' || c == 'Y'
            return(true)
        elseif c == 'z' || c == 'Z'
            return(true)
        else
            return(false)
        end
    end # function is_letter

    # =========================================================================

    # Check if the character is a valid digit (0-9).
    function is_digit(c::Char)
        if c == '0'
            return(true)
        elseif c == '1'
            return(true)
        elseif c == '2'
            return(true)
        elseif c == '3'
            return(true)
        elseif c == '4'
            return(true)
        elseif c == '5'
            return(true)
        elseif c == '6'
            return(true)
        elseif c == '7'
            return(true)
        elseif c == '8'
            return(true)
        elseif c == '9'
            return(true)
        else
            return(false)
        end
    end # function is_digit

    # =========================================================================

    # Convert valid digit characters to the correct number
    function to_digit(c::Char)
        if c == '0'
            return(0::Int64)
        elseif c == '1'
            return(1::Int64)
        elseif c == '2'
            return(2::Int64)
        elseif c == '3'
            return(3::Int64)
        elseif c == '4'
            return(4::Int64)
        elseif c == '5'
            return(5::Int64)
        elseif c == '6'
            return(6::Int64)
        elseif c == '7'
            return(7::Int64)
        elseif c == '8'
            return(8::Int64)
        elseif c == '9'
            return(9::Int64)
        else
            throw(ParserError(join(["Character \"", string(c),
                "\" is not a digit."])))
        end
    end # function to_digit

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    # skip ahead to end of whitespace -- return number of whitespace characters
    function consume_whitespace(s::String)
        idx = 1
        while idx <= length(s) && (s[idx] == ' ' || s[idx] == '\t')
            idx += 1
        end
        return(s[idx:end])
    end # function consume_whitespace

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    # Attempt to parse a string into a string Value (which follows a particular
    # format).  Return the string and the trailing characters.  Throw a
    # ParserError if parsing fails.
    function parse_string_value(s::String)
        idx = 1
        if s[idx] != '"'
            throw(ParserError("Unable to parse as a string."))
        end
        idx += 1
        temp_str::UTF8String = ""
        while idx <= length(s)
            if s[idx] == '\\'
                idx += 1
                if s[idx] == '"'
                    temp_str = join([temp_str, '"'])
                    idx += 1
                elseif s[idx] == '\\'
                    temp_str = join([temp_str, '\\'])
                    idx += 1
                elseif s[idx] == 'b'
                    temp_str = join([temp_str, "\b"])
                    idx += 1
                elseif s[idx] == 't'
                    temp_str = join([temp_str, "\t"])
                    idx += 1
                elseif s[idx] == 'n'
                    temp_str = join([temp_str, "\n"])
                    idx += 1
                elseif s[idx] == 'f'
                    temp_str = join([temp_str, "\f"])
                    idx += 1
                elseif s[idx] == 'r'
                    temp_str = join([temp_str, "\r"])
                    idx += 1
                else
                    throw(ParserError(join(["Unknown escale character \"\\",
                        s[idx], "\"."])))
                end
            elseif s[idx] == '"'
                break
            else
                temp_str = join([temp_str, s[idx]])
                idx += 1
            end
        end
        if s[idx] != '"'
            throw(ParserError("Unable to parse as a string."))
        end
        idx += 1
        return(temp_str, s[idx:end])
    end # function parse_string_value

    # =========================================================================

    # Attempt to parse a string as a number.  Return the NumberValue struct or
    # raise a ParserError.
    function parse_number_value(s::String)
        temp_num::NumberValue = NumberValue()
        temp_num.valid_integer = false
        temp_num.integer_value = 0
        temp_num.valid_float = false
        temp_num.float_value = 0.0
        idx = 1
        # sign
        sign::Int64
        if s[idx] == '-'
            sign = -1
            idx += 1
        elseif s[idx] == '+'
            sign = 1
            idx += 1
        elseif s[idx] == '.' || is_digit(s[idx])
            sign = 1
        else
            throw(ParserError("Unable to parse as a number."))
        end
        # integer part
        ipart::Int64 = 0
        while idx <= length(s) && is_digit(s[idx])
            next_ipart = 10 * ipart + to_digit(s[idx])
            idx += 1
            if next_ipart <= ipart
                throw(ParserError("Unable to parse as a number."))
            end
            ipart = next_ipart
        end
        # decimal
        dpart::Float64 = 0.0
        shift::Float64 = 0.1
        if idx <= length(s) && s[idx] == '.'
            idx += 1
            while idx <= length(s) && is_digit(s[idx])
                dpart += shift * to_digit(s[idx])
                shift *= 0.1
                idx += 1
            end
        end
        # exponent (scientific notation)
        e_sign::Int64
        exponent::Int64 = 0
        if idx <= length(s) && (s[idx] == 'e' || s[idx] == 'E')
            idx += 1
            if s[idx] == '-'
                e_sign = -1
                idx += 1
            elseif s[idx] == '+'
                e_sign = 1
                idx += 1
            elseif is_digit(s[idx])
                e_sign = 1
            else
                throw(ParserError("Invalid exponent in number."))
            end
            while idx <= length(s) && is_digit(s[idx])
                exponent = 10 * exponent + to_digit(s[idx])
                idx += 1
            end
            exponent *= e_sign
        end
        # construct the number
        as_integer::Int64 = 0
        as_float::Float64 = 0.0
        if dpart == 0 && exponent == 0
            # This is really an integer, and may also be a float
            as_integer = sign * ipart
            as_float = 0.0
            try
                as_float = convert(Float64, as_integer)
                back_to_integer::Int64 = convert(Int64, as_float)
                if back_to_integer == as_integer
                    temp_num.float_value = as_float
                    temp_num.valid_float = true
                end
            catch e
                isa(e, InexactError) || rethrow()
            end
            temp_num.integer_value = as_integer
            temp_num.valid_integer = true
        else
            # This is really a float, and may also be an integer
            as_float = sign * (convert(Float64, ipart) + dpart) * 10.0^exponent
            as_integer = 0
            try
                as_integer = convert(Int64, as_float)
                back_to_float::Float64 = convert(Float64, as_integer)
                if back_to_float == as_float
                    temp_num.integer_value = as_integer
                    temp_num.valid_integer = true
                end
            catch e
                isa(e, InexactError) || rethrow()
            end
            temp_num.float_value = as_float
            temp_num.valid_float = true
        end
        return(temp_num, s[idx:end])
    end # function parse_number_value

    # =========================================================================

    # Attempt to parse a string as a boolean.  Return the boolean or raise a
    # ParserError.
    function parse_boolean_value(s::String)
        temp_bool::Bool = false
        if length(s) >= 4 && s[1:4] == "true"
            idx = 5
            temp_bool = true
        elseif length(s) >= 5 && s[1:5] == "false"
            idx = 6
            temp_bool = false
        else
            throw(ParserError("Unable to parse as a boolean."))
        end
        return(temp_bool, s[idx:end])
    end # function parse_boolean_value

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    type Value
        # The value in different formats
        value_as_string::UTF8String
        value_as_integer::Int64
        value_as_float::Float64
        value_as_boolean::Bool
        # Is the value available in different formats?
        is_conformable_to_string::Bool
        is_conformable_to_integer::Bool
        is_conformable_to_float::Bool
        is_conformable_to_boolean::Bool
    end # type Value

    # =========================================================================

    function Value()
        return(Value("", 0, 0.0, false, false, false, false, false))
    end # function Value()

    # =========================================================================

    function clear(v::Value)
        # not conformable to anything
        v.is_conformable_to_string = false
        v.is_conformable_to_integer = false
        v.is_conformable_to_float = false
        v.is_conformable_to_boolean = false
        # wipe the value(s)
        v.value_as_string = ""
        v.value_as_integer = 0
        v.value_as_float = 0.0
        v.value_as_boolean = false
    end # function clear

    # =========================================================================

    # Set the Value based on parsing the String
    function analyze(v::Value, s::String)
        # Clear the current internal values and flags
        clear(v)

        # Trim leading whitespace
        s = consume_whitespace(s)

        # Ensure there is something (non-comment) left in the string
        if length(s) == 0 || s[1] == '#'
            throw(ParserError("Empty value."))
        end

        # Choose which type to parse
        if s[1] == '"'
            # This is either a String or nothing
            temp_string::UTF8String, s = parse_string_value(s)
            # Save it
            v.value_as_string = temp_string
            v.is_conformable_to_string = true
        elseif s[1] == 't' || s[1] == 'f'
            # This is either a Boolean or nothing
            temp_boolean::Bool, s = parse_boolean_value(s)
            # Save it
            v.value_as_boolean = temp_boolean
            is_conformable_to_boolean = true
        elseif s[1] == '-' || s[1] == '+' || s[1] == '.' || is_digit(s[1])
            # This is either an Integer, a Float, both, or nothing
            temp_number::NumberValue, s = parse_number_value(s)
            # Save it
            v.value_as_integer = temp_number.integer_value
            v.value_as_float = temp_number.float_value
            v.is_conformable_to_integer = temp_number.valid_integer
            v.is_conformable_to_float = temp_number.valid_float
        else
            # This is nothing
            throw(ParserError(join(["Unable to parse \"", s,
                "\" to a value."])))
        end

        # Return the remainder of the string
        return s
    end # function analyze

    # =========================================================================

    # Set the Value from a String
    function set(v::Value, s::UTF8String)
        clear(v)
        v.value_as_string = s
        v.is_conformable_to_string = true
    end # function set(Value, UTF8String)

    # =========================================================================

    # Set the Value from an Integer
    function set(v::Value, i::Int64)
        clear(v)
        v.value_as_integer = i
        v.is_conformable_to_integer = true
        as_float::Float64 = 0.0
        try
            as_float = convert(Float64, i)
            back_to_integer::Int64 = convert(Int64, as_float)
            if back_to_integer == i
                v.value_as_float = as_float
                v.is_conformable_to_float = true
            end
        catch e
            isa(e, InexactError) || rethrow()
        end
    end # function set(Value, Int64)

    # =========================================================================

    # Set the Value from a Float
    function set(v::Value, f::Float64)
        clear(v)
        v.value_as_flaot = f
        v.is_conformable_to_float = true
        as_integer::Int64 = 0
        try
            as_integer = convert(Int64, f)
            back_to_float::Float64 = convert(Float64, as_integer)
            if back_to_float == f
                v.value_as_integer = as_integer
                v.is_conformable_to_integer = true
            end
        catch e
            isa(e, InexactError) || rethrow()
        end
    end # function set(Value, Float)

    # =========================================================================

    # Set the Value from a Boolean
    function set(v::Value, b::Bool)
        clear(v)
        v.value_as_boolean = b
        v.is_conformable_to_boolean = true
    end # function set(Value, Boolean)

    # =========================================================================

    # Return the Value as a String
    function as_string(v::Value)
        if v.is_conformable_to_string
            return(v.value_as_string)
        else
            throw(TypeConversionError(
                "Value cannot be converted to a string."))
        end
    end # function as_string

    # =========================================================================

    # Return the Value as an Integer
    function as_integer(v::Value)
        if v.is_conformable_to_integer
            return(v.value_as_integer)
        else
            throw(TypeConversionError(
                "Value cannot be converted to an integer."))
        end
    end # function as_integer

    # =========================================================================

    # Return the Value as a Float
    function as_float(v::Value)
        if v.is_conformable_to_float
            return(v.value_as_float)
        else
            throw(TypeConversionError("Value cannot be converted to a float."))
        end
    end # function as_float

    # =========================================================================

    # Return the Value as a Boolean
    function as_boolean(v::Value)
        if v.is_conformable_to_boolean
            return(v.value_as_boolean)
        else
            throw(TypeConversionError(
                "Value cannot be converted to a boolean."))
        end
    end # function as_boolean

    # =========================================================================

    # Convert the Value to a string as if writing a new TOML file
    function serialize(v::Value)
        if v.is_conformable_to_boolean
            # Write as a Boolean
            if v.value_as_boolean
                return("true")
            else
                return("false")
            end
        elseif v.is_conformable_to_integer
            # Write as an Integer (anything conformable to both Integer and
            # Float will appear as an Integer because of the ordering here)
            return(string(v.value_as_integer))
        elseif v.is_conformable_to_float
            # Write as a Float
            return(string(v.value_as_float))
        elseif v.is_conformable_to_string
            idx = 1
            output = "\""
            while idx <= length(v.value_as_string)
                if v.value_as_string[idx] == '"'
                    output = join([output, "\\\""])
                elseif v.value_as_string[idx] == '\\'
                    output = join([output, "\\\\"])
                elseif v.value_as_string[idx] == '\b'
                    output = join([output, "\\\b"])
                elseif v.value_as_string[idx] == '\t'
                    output = join([output, "\\\t"])
                elseif v.value_as_string[idx] == '\n'
                    output = join([output, "\\\n"])
                elseif v.value_as_string[idx] == '\f'
                    output = join([output, "\\\f"])
                elseif v.value_as_string[idx] == '\r'
                    output = join([output, "\\\r"])
                else
                    output = join([output, v.value_as_string[idx]])
                end
                idx += 1
            end
            output = join([output, "\""])
            return(output)
        else
            # Not actually a valid Value
            throw(ValueError("Value cannot be serialized."))
        end
    end # function serialize

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    type ValueArray
        array::Array{Value,1}
        is_conformable_to_string::Bool
        is_conformable_to_integer::Bool
        is_conformable_to_float::Bool
        is_conformable_to_boolean::Bool
    end # type ValueArray

    # =========================================================================

    function ValueArray()
        return(ValueArray([], false, false, false, false))
    end # function ValueArray

    # =========================================================================

    function size(va::ValueArray)
        return(length(va.array))
    end # function size

    # -------------------------------------------------------------------------
    # =========================================================================
    # -------------------------------------------------------------------------

    function clear(va::ValueArray)
        while length(va.array) > 0
            deleteat!(va.array, 1)
        va.is_conformable_to_string = false
        va.is_conformable_to_integer = false
        va.is_conformable_to_float = false
        va.is_conformable_to_boolean = false
    end # function clear

    # =========================================================================

    function add(va::ValueArray, v::Value)
        if length(va) == 0
            push!(va, v)
            va.is_conformable_to_string = v.is_conformable_to_string
            va.is_conformable_to_integer = v.is_conformable_to_integer
            va.is_conformable_to_float = v.is_conformable_to_float
            va.is_conformable_to_boolean = v.is_conformable_to_boolean
        else
            if va.is_conformable_to_string && v.is_conformable_to_string
                push!(va, v)
            elseif va.is_conformable_to_integer && v.is_conformable_to_integer
                push!(va, v)
            elseif va.is_conformable_to_float && v.is_conformable_to_float
                push!(va, v)
            elseif va.is_conformable_to_boolean && v.is_conformable_to_boolean
                push!(va, v)
            else
                throw(ValueError(
                    "Value with invalid type cannot be added to ValueArray."))
            end
        end
    end # function add(ValueArray, Value)

    # =========================================================================

    function remove(va::ValueArray, index::Integer)
        # The deleteat! function already handles bounds checking for me
        deleteat!(va.array, index)
    end # function remove

    # =========================================================================

    function as_string(va::ValueArray)
        arr::Array{UTF8String,1} = []
        for index = 1:length(va)
            push!(arr, va[index])
        end
        return(arr)
    end # function as_string(ValueArray)

    # =========================================================================

end # module TOML

max_stars_repo_name: bkkrueger/config

max_issues_repo_name: bkkrueger/TOML

max_forks_repo_name: bkkrueger/TOML


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   DTCFILE.txt
#
#   This file stores the temperature variation `ΔTc` caused by the `Dst`.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

################################################################################
#                       Private Structures and Variables
################################################################################

"""
    _DTCFILE_Structure

Structure to store the interpolations of the data in `DTCFILE.TXT` file.

# Fields

* `DstΔTc`: Temperature variation due to Dst [K].

"""
struct _DTCFILE_Structure
    DstΔTc::_space_indices_itp_linear{Float64}
end

# Remote file: DTCFILE.TXT
_dtcfile = @RemoteFile(
    "http://sol.spacenvironment.net/jb2008/indices/DTCFILE.TXT",
    file="DTCFILE.TXT",
    dir="$path_si_root",
    updates=:daily
)

# Optional variable that will store the `DTCFILE.TXT` data.
@OptionalData(_dtcfile_data, _DTCFILE_Structure,
              "Run `init_space_indices()` with `:dtcfile` in `enabled_files` array to initialize required data.")

################################################################################
#                               Public Functions
################################################################################

#                                   Getters
# ==============================================================================

# Create a function `get_S` for each `S` field in `_DTCFILE_Structure`.
for field in fieldnames(_DTCFILE_Structure)
    field_name = string(field)
    func_name  = "get_" * field_name
    func       = Symbol("get_" * string(field))
    qfield     = Meta.quot(field)

    @eval begin
        export $func
"""
    $($func_name)(DT::Number)

Get the value of the index `$($field_name)` at Julian Day `DT`.

This function requires the initialization of the variable `_dtcfile_data`.
Otherwise, an exception will be raised. To initialize it, run
`init_space_indices()`.

"""
        ($func)(DT::Number) = getfield(get(_dtcfile_data), $qfield)(DT)
    end
end

################################################################################
#                              Private Functions
################################################################################

"""
    _init_dctfile(;force_download = false, local_path = nothing)

Initialize the data in the file `DTCFILE.TXT` by creating `_dtcfile_data`. The
initialization process is composed of:

1. Download the file, if it is necessary;
2. Parse the file;
3. Create the interpolations and the structures.

If the keyword `force_download` is `true`, then the file will always be
downloaded.

The user can also specify a location for the file using the keyword
`local_path`. If it is `nothing`, which is the default, then the file will be
downloaded.

"""
function _init_dtcfile(;force_download = false, local_path = nothing)
    # Update the remote files if no path is given.
    if local_path == nothing
        download(_dtcfile; force = force_download)
        local_path = path(_dtcfile)
    end

    push!(_dtcfile_data,   _parse_dtcfile(local_path))

    nothing
end

"""
    _parse_dtcfile(path::AbstractString)

Parse the `DTCFILE.TXT` file in `path` and return an instance of the structure
`_DTCFILE_Structure` with the initialized interpolations.

The format of the file `DTCFILE.TXT` must be:

    DTC YYYY DOY DTC_0h DTC_1h DTC_2h ... DTC_22h DTC_23h

in which `DOY` is the day of the year and `DTC_Xh` is the `ΔTc` at hour `X`.

"""
function _parse_dtcfile(path::AbstractString)
    # Allocate the raw data.
    DT     = DateTime[]
    DstΔTc = Float64[]

    open(path) do file
        line_num = 0

        for ln in eachline(file)
            line_num += 1

            # Ignore comments.
            (ln[1] == '#') && continue

            tokens   = split(ln)
            num_toks = length(tokens)

            # Ignore blank lines.
            (num_toks == 0) && continue

            # Ignore and print warning about lines with bad format.
            if (num_toks != 27) || (tokens[1] != "DTC")
                @warn "Line $line_num of file DTCFILE.TXT has invalid format!"
                continue
            end

            # Compute the Julian day at midnight of the day.
            year = parse(Int64, tokens[2])
            doy  = parse(Int64, tokens[3])

            DT_0h = DateTime(year,) + Day(doy-1)

            # Parse the data.
            for k = 1:24
                DstΔTc_h = parse(Float64, tokens[k+3])
                push!(DT,  DT_0h + Hour(k-1))
                push!(DstΔTc, DstΔTc_h)
            end

        end
    end

    # Create the interpolations for each parameter.
    knots      = (datetime2julian.(DT),)
    itp_DstΔTc = interpolate(knots, DstΔTc, Gridded(Linear()))

    _DTCFILE_Structure(itp_DstΔTc)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    GitRevWalker(repo::GitRepo)

A `GitRevWalker` *walks* through the *revisions* (i.e. commits) of
a git repository `repo`. It is a collection of the commits
in the repository, and supports iteration and calls to [`map`](@ref LibGit2.map)
and [`count`](@ref LibGit2.count) (for instance, `count` could be used to determine
what percentage of commits in a repository were made by a certain
author).

```julia
cnt = LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
    count((oid,repo)->(oid == commit_oid1), walker, oid=commit_oid1, by=LibGit2.Consts.SORT_TIME)
end
```
Here, `count` finds the number of commits along the walk with a certain `GitHash`.
Since the `GitHash` is unique to a commit, `cnt` will be `1`.
"""
function GitRevWalker(repo::GitRepo)
    ensure_initialized()
    w_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_revwalk_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}), w_ptr, repo.ptr)
    return GitRevWalker(repo, w_ptr[])
end

function Base.iterate(w::GitRevWalker, state=nothing)
    ensure_initialized()
    id_ptr = Ref(GitHash())
    err = ccall((:git_revwalk_next, :libgit2), Cint,
                (Ptr{GitHash}, Ptr{Cvoid}), id_ptr, w.ptr)
    if err == Cint(Error.GIT_OK)
        return (id_ptr[], nothing)
    elseif err == Cint(Error.ITEROVER)
        return nothing
    else
        throw(GitError(err))
    end
end

Base.IteratorSize(::Type{GitRevWalker}) = Base.SizeUnknown()

"""
    LibGit2.push_head!(w::GitRevWalker)

Push the HEAD commit and its ancestors onto the [`GitRevWalker`](@ref)
`w`. This ensures that HEAD and all its ancestor commits will be encountered
during the walk.
"""
function push_head!(w::GitRevWalker)
    ensure_initialized()
    @check ccall((:git_revwalk_push_head, :libgit2), Cint, (Ptr{Cvoid},), w.ptr)
    return w
end

"""
    LibGit2.push!(w::GitRevWalker, cid::GitHash)

Start the [`GitRevWalker`](@ref) `walker` at commit `cid`. This function can be used
to apply a function to all commits since a certain year, by passing the first commit
of that year as `cid` and then passing the resulting `w` to [`map`](@ref LibGit2.map).
"""
function push!(w::GitRevWalker, cid::GitHash)
    ensure_initialized()
    @check ccall((:git_revwalk_push, :libgit2), Cint, (Ptr{Cvoid}, Ptr{GitHash}), w.ptr, Ref(cid))
    return w
end

function push!(w::GitRevWalker, range::AbstractString)
    ensure_initialized()
    @check ccall((:git_revwalk_push_range, :libgit2), Cint, (Ptr{Cvoid}, Ptr{UInt8}), w.ptr, range)
    return w
end

function Base.sort!(w::GitRevWalker; by::Cint = Consts.SORT_NONE, rev::Bool=false)
    ensure_initialized()
    rev && (by |= Consts.SORT_REVERSE)
    ccall((:git_revwalk_sorting, :libgit2), Cvoid, (Ptr{Cvoid}, Cint), w.ptr, by)
    return w
end

repository(w::GitRevWalker) = w.owner

"""
    LibGit2.map(f::Function, walker::GitRevWalker; oid::GitHash=GitHash(), range::AbstractString="", by::Cint=Consts.SORT_NONE, rev::Bool=false)

Using the [`GitRevWalker`](@ref) `walker` to "walk" over every commit in the repository's history,
apply `f` to each commit in the walk. The keyword arguments are:
    * `oid`: The [`GitHash`](@ref) of the commit to begin the walk from. The default is to use
      [`push_head!`](@ref) and therefore the HEAD commit and all its ancestors.
    * `range`: A range of `GitHash`s in the format `oid1..oid2`. `f` will be
      applied to all commits between the two.
    * `by`: The sorting method. The default is not to sort. Other options are to sort by
      topology (`LibGit2.Consts.SORT_TOPOLOGICAL`), to sort forwards in time
      (`LibGit2.Consts.SORT_TIME`, most ancient first) or to sort backwards in time
      (`LibGit2.Consts.SORT_REVERSE`, most recent first).
    * `rev`: Whether to reverse the sorted order (for instance, if topological sorting is used).

# Examples
```julia
oids = LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
    LibGit2.map((oid, repo)->string(oid), walker, by=LibGit2.Consts.SORT_TIME)
end
```
Here, `map` visits each commit using the `GitRevWalker` and finds its `GitHash`.
"""
function map(f::Function, walker::GitRevWalker;
                  oid::GitHash=GitHash(),
                  range::AbstractString="",
                  by::Cint = Consts.SORT_NONE,
                  rev::Bool=false,
                  count::Int=0)
    res = []
    sort!(walker, by=by, rev=rev)
    if !iszero(oid)
        push!(walker, oid)
    elseif !isempty(range)
        push!(walker, range)
    else
        push_head!(walker)
    end

    repo = repository(walker)
    for val in (count == 0 ? walker : Iterators.take(walker, count))
        Base.push!(res, f(val, repo))
    end
    return res
end

"""
    LibGit2.count(f::Function, walker::GitRevWalker; oid::GitHash=GitHash(), by::Cint=Consts.SORT_NONE, rev::Bool=false)

Using the [`GitRevWalker`](@ref) `walker` to "walk" over every commit in the repository's history,
find the number of commits which return `true` when `f` is applied to them. The keyword arguments
are:
    * `oid`: The [`GitHash`](@ref) of the commit to begin the walk from. The default is to use
      [`push_head!`](@ref) and therefore the HEAD commit and all its ancestors.
    * `by`: The sorting method. The default is not to sort. Other options are to sort by
      topology (`LibGit2.Consts.SORT_TOPOLOGICAL`), to sort forwards in time
      (`LibGit2.Consts.SORT_TIME`, most ancient first) or to sort backwards in time
      (`LibGit2.Consts.SORT_REVERSE`, most recent first).
    * `rev`: Whether to reverse the sorted order (for instance, if topological sorting is used).

# Examples
```julia
cnt = LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
    count((oid, repo)->(oid == commit_oid1), walker, oid=commit_oid1, by=LibGit2.Consts.SORT_TIME)
end
```
`count` finds the number of commits along the walk with a certain `GitHash` `commit_oid1`, starting
the walk from that commit and moving forwards in time from it. Since the `GitHash` is unique to
a commit, `cnt` will be `1`.
"""
function count(f::Function, walker::GitRevWalker;
                  oid::GitHash=GitHash(),
                  by::Cint = Consts.SORT_NONE,
                  rev::Bool=false)
    c = 0
    sort!(walker, by=by, rev=rev)
    if !iszero(oid)
        push!(walker, oid)
    else
        push_head!(walker)
    end

    repo = repository(walker)
    for val in walker
        c += f(val, repo) == true
    end
    return c
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# ---
# title: 1563. Stone Game V
# id: problem1563
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/stone-game-v/description/>
# hidden: true
# ---
#
# There are several stones  **arranged in a row** , and each stone has an
# associated value which is an integer given in the array `stoneValue`.
#
# In each round of the game, Alice divides the row into **two non-empty rows**
# (i.e. left row and right row), then Bob calculates the value of each row which
# is the sum of the values of all the stones in this row. Bob throws away the
# row which has the maximum value, and Alice's score increases by the value of
# the remaining row. If the value of the two rows are equal, Bob lets Alice
# decide which row will be thrown away. The next round starts with the remaining
# row.
#
# The game ends when there is only **one stone remaining**. Alice's is initially
# **zero**.
#
# Return _the maximum score that Alice can obtain_.
#
#
#
# **Example 1:**
#
#
#
#     Input: stoneValue = [6,2,3,4,5,5]
#     Output: 18
#     Explanation: In the first round, Alice divides the row to [6,2,3], [4,5,5]. The left row has the value 11 and the right row has value 14. Bob throws away the right row and Alice's score is now 11.
#     In the second round Alice divides the row to [6], [2,3]. This time Bob throws away the left row and Alice's score becomes 16 (11 + 5).
#     The last round Alice has only one choice to divide the row which is [2], [3]. Bob throws away the right row and Alice's score is now 18 (16 + 2). The game ends because only one stone is remaining in the row.
#
#
# **Example 2:**
#
#
#
#     Input: stoneValue = [7,7,7,7,7,7,7]
#     Output: 28
#
#
# **Example 3:**
#
#
#
#     Input: stoneValue = [4]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= stoneValue.length <= 500`
#   * `1 <= stoneValue[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end