text:
##########################
# functions to benchmark #
##########################

function perf_sumelt(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    return s
end

# bounds-checking is deliberately on,
# doesn't use `for a in r`
function perf_sumelt_boundscheck(A)
    s = zero(eltype(A))
    for i = 1:length(A)
        s += A[i]
    end
    return s
end

function perf_sumeach(A)
    s = zero(eltype(A))
    for I in eachindex(A)
        val = Base.unsafe_getindex(A, I)
        s += val
    end
    return s
end

function perf_sumlinear(A)
    s = zero(eltype(A))
    for I in 1:length(A)
        val = Base.unsafe_getindex(A, I)
        s += val
    end
    return s
end
function perf_sumcartesian(A)
    s = zero(eltype(A))
    for I in CartesianIndices(size(A))
        val = Base.unsafe_getindex(A, I)
        s += val
    end
    return s
end

function perf_sumcolon(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    c = Colon()
    @simd for i = 1:ncols
        val = Base.unsafe_getindex(A, c, i)
        s += first(val)
    end
    return s
end

function perf_sumrange(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = 1:nrows
    @simd for i = 1:ncols
        val = Base.unsafe_getindex(A, r, i)
        s += first(val)
    end
    return s
end

function perf_sumlogical(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = falses(nrows)
    r[1:4:end] .= true
    @simd for i = 1:ncols
        val = Base.unsafe_getindex(A, r, i)
        s += first(val)
    end
    return s
end

function perf_sumvector(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = rand(1:nrows, 5)
    @simd for i = 1:ncols
        val = Base.unsafe_getindex(A, r, i)
        s += first(val)
    end
    return s
end

## Views
function perf_sumeach_view(A)
    s = zero(eltype(A))
    @inbounds @simd for I in eachindex(A)
        val = view(A, I)
        s += val[]
    end
    return s
end

function perf_sumlinear_view(A)
    s = zero(eltype(A))
    @inbounds @simd for I in 1:length(A)
        val = view(A, I)
        s += val[]
    end
    return s
end
function perf_sumcartesian_view(A)
    s = zero(eltype(A))
    @inbounds @simd for I in CartesianIndices(size(A))
        val = view(A, I)
        s += val[]
    end
    return s
end

function perf_sumcolon_view(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    c = Colon()
    @inbounds for i = 1:ncols
        val = view(A, c, i)
        s += first(val)
    end
    return s
end

function perf_sumrange_view(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = 1:nrows
    @inbounds for i = 1:ncols
        val = view(A, r, i)
        s += first(val)
    end
    return s
end

function perf_sumlogical_view(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = falses(nrows)
    r[1:4:end] .= true
    @inbounds for i = 1:ncols
        val = view(A, r, i)
        s += first(val)
    end
    return s
end

function perf_sumvector_view(A)
    s = zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = rand(1:nrows, 5)
    @inbounds for i = 1:ncols
        val = view(A, r, i)
        s += first(val)
    end
    return s
end

function perf_sub2ind(sz, irange, jrange, krange)
    linear = LinearIndices(sz)
    s = 0
    for k in krange, j in jrange, i in irange
        ind = linear[i, j, k]
        s += ind
    end
    s
end

function perf_ind2sub(sz, lrange)
    cart = CartesianIndices(sz)
    si = sj = sk = 0
    for l in lrange
        i, j, k = cart[l].I  # TODO: change to Tuple(cart[l])
        si += i
        sj += j
        sk += k
    end
    si, sj, sk
end

function setup_mapr_access(A)
    z = zero(eltype(A))
    zz = sum(z -> z * z, [z]) # z = z*z, with any promotion from sum
    n = minimum(size(A))
    B = Vector{typeof(zz)}(undef, n)
    B, zz, n
end
if VERSION < v"0.7.0-beta-81"
    function perf_mapr_access(A, B, zz, n) #20517
        @inbounds for j in 1:n
            B[j] = mapreduce(k -> A[j,k]*A[k,j], +, zz, 1:j)
        end
        B
    end
else
    function perf_mapr_access(A, B, zz, n) #20517
        @inbounds for j in 1:n
            B[j] = mapreduce(k -> A[j,k]*A[k,j], +, 1:j; init=zz)
        end
        B
    end
end

##########################
# supporting definitions #
##########################

abstract type MyArray{T,N} <: AbstractArray{T,N} end

struct ArrayLS{T,N} <: MyArray{T,N}  # LinearSlow
    data::Array{T,N}
end

struct ArrayLSLS{T,N} <: MyArray{T,N}  # LinearSlow with LinearSlow similar
    data::Array{T,N}
end

struct ArrayLF{T,N} <: MyArray{T,N}  # LinearFast
    data::Array{T,N}
end

struct ArrayStrides{T,N} <: MyArray{T,N}
    data::Array{T,N}
    strides::NTuple{N,Int}
end

ArrayStrides(A::Array) = ArrayStrides(A, strides(A))

struct ArrayStrides1{T} <: MyArray{T,2}
    data::Matrix{T}
    stride1::Int
end

ArrayStrides1(A::Array) = ArrayStrides1(A, size(A,1))

Base.@propagate_inbounds Base.setindex!(A::ArrayLSLS, v, I::Int...) = A.data[I...] = v

Base.@propagate_inbounds Base.unsafe_setindex!(A::ArrayLSLS, v, I::Int...) = Base.unsafe_setindex!(A.data, v, I...)

Base.first(A::ArrayLSLS) = first(A.data)

Base.size(A::MyArray) = size(A.data)
# To ensure that ArrayLS and ArrayLSLS really are LinearSlow even
# after inlining, let's make the size differ from the parent array
# (ref https://github.com/JuliaLang/julia/pull/17355#issuecomment-231748251)
@inline Base.size(A::ArrayLS{T,1}) where {T}   = (sz = size(A.data); (sz[1]-1,))
@inline Base.size(A::ArrayLSLS{T,1}) where {T} = (sz = size(A.data); (sz[1]-1,))
@inline Base.size(A::ArrayLS{T,2}) where {T}   = (sz = size(A.data); (sz[1]-1,sz[2]-1))
@inline Base.size(A::ArrayLSLS{T,2}) where {T} = (sz = size(A.data); (sz[1]-1,sz[2]-1))
@inline Base.size(A::ArrayLS)   = map(n->n-1, size(A.data))
@inline Base.size(A::ArrayLSLS) = map(n->n-1, size(A.data))

@inline Base.similar(A::ArrayLSLS, ::Type{T}, dims::Tuple{Int}) where {T}     = ArrayLSLS(similar(A.data, T, (dims[1]+1,)))
@inline Base.similar(A::ArrayLSLS, ::Type{T}, dims::Tuple{Int,Int}) where {T} = ArrayLSLS(similar(A.data, T, (dims[1]+1,dims[2]+1)))
@inline Base.similar(A::ArrayLSLS, ::Type{T}, dims::Tuple{Vararg{Int}}) where {T} = ArrayLSLS(similar(A.data, T, map(n->n+1, dims)))

Base.@propagate_inbounds Base.getindex(A::ArrayLF, i::Int) = getindex(A.data, i)
Base.@propagate_inbounds Base.getindex(A::ArrayLS{T,2}, i::Int, j::Int) where {T} = getindex(A.data, i, j)
Base.@propagate_inbounds Base.getindex(A::ArrayLS{T,3}, i::Int, j::Int, k::Int) where {T} = getindex(A.data, i, j, k)
Base.@propagate_inbounds Base.getindex(A::ArrayLSLS{T,2}, i::Int, j::Int) where {T} = getindex(A.data, i, j)
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayLF, indx::Int) = Base.unsafe_getindex(A.data, indx)
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayLS{T,2}, i::Int, j::Int) where {T} = Base.unsafe_getindex(A.data, i, j)
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayLS{T,3}, i::Int, j::Int, k::Int) where {T} = Base.unsafe_getindex(A.data, i, j, k)
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayLSLS{T,2}, i::Int, j::Int) where {T} = Base.unsafe_getindex(A.data, i, j)

Base.@propagate_inbounds Base.getindex(A::ArrayStrides{T,2}, i::Real, j::Real) where {T} = getindex(A.data, 1+A.strides[1]*(i-1)+A.strides[2]*(j-1))
Base.@propagate_inbounds Base.getindex(A::ArrayStrides1, i::Real, j::Real) = getindex(A.data, i + A.stride1*(j-1))
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayStrides{T,2}, i::Real, j::Real) where {T} = Base.unsafe_getindex(A.data, 1+A.strides[1]*(i-1)+A.strides[2]*(j-1))
Base.@propagate_inbounds Base.unsafe_getindex(A::ArrayStrides1, i::Real, j::Real) = Base.unsafe_getindex(A.data, i + A.stride1*(j-1))

Base.IndexStyle(::Type{<:ArrayLF}) = IndexLinear()

if !applicable(Base.unsafe_getindex, [1 2], 1:1, 2)
    Base.@propagate_inbounds Base.unsafe_getindex(A::Array, I...) = @inbounds return A[I...]
    Base.@propagate_inbounds Base.unsafe_getindex(A::MyArray, I...) = @inbounds return A[I...]
    Base.@propagate_inbounds Base.unsafe_getindex(A::SubArray, I...) = @inbounds return A[I...]
    Base.@propagate_inbounds Base.unsafe_getindex(A::BitArray, I1::BitArray, I2::Int) = Base.unsafe_getindex(A, Base.to_index(I1), I2)
end

function makearrays(::Type{T}, r::Integer, c::Integer) where T
    A = samerand(T, r, c)
    B = similar(A, r+1, c+1)
    B[1:r, 1:c] = A
    AS = ArrayLS(B)
    ASS = ArrayLSLS(B)
    AF = ArrayLF(A)
    Astrd = ArrayStrides(A)
    Astrd1 = ArrayStrides1(A)
    B = samerand(T, r+1, c+2)
    # And views thereof
    Asub = view(B, 1:r, 2:c+1)
    Asub2 = view(A, :, :)
    Asub3 = view(AS, :, :)
    C = samerand(T, 4, r, c)
    Asub4 = view(C, 1, :, :)
    Asub5 = view(ArrayLS(C), 1, :, :)
    Asub6 = view(reshape(view(C, :, :, :), Val(2)), :, 2:c+1)
    Asub7 = view(reshape(view(ArrayLS(C), :, :, :), Val(2)), :, 2:c+1)

    return (A, AF, AS, ASS, Asub, Asub2, Asub3, Asub4, Asub5, Asub6, Asub7)
end


text: square(x) = exit(0)

ext: jl


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(realpath $(dirname $0))))/examples/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "FaceDetection.jl")) # ../src/FaceDetection.jl

using .FaceDetection
using Test: @testset, @test
using Suppressor: @suppress
# using BenchmarkTools: @btime

@time @testset "FaceDetection.jl" begin
	# constants and variables
	main_data_path = joinpath(@__DIR__, "images")
	pos_training_path = joinpath(main_data_path, "pos")
	neg_training_path = joinpath(main_data_path, "neg")
	pos_testing_path = joinpath(main_data_path, "pos_testing")
	neg_testing_path = joinpath(main_data_path, "neg_testing")
	a, b, c, d, e, f = tuple(rand(Int), rand(Int)), tuple(rand(Int), rand(Int)), rand(Int), rand(Int), rand((0, 1)), rand((0, 1))
    arr = rand(Int, 100, 100)
	int_img = load_image(rand(vcat(filtered_ls.([pos_training_path, neg_training_path, pos_testing_path, neg_testing_path])...)), scale = true, scale_to = (24, 24))
	feature_2v = HaarLikeObject(feature_types.two_vertical, (1, 1), 10, 10, 100000, 1)
	feature_2h = HaarLikeObject(feature_types.two_horizontal, (1, 1), 10, 10, 100000, 1)
	feature_3h = HaarLikeObject(feature_types.three_horizontal, (1, 1), 10, 10, 100000, 1)
	feature_3v = HaarLikeObject(feature_types.two_vertical, (1, 1), 10, 10, 100000, 1)
	feature_4 = HaarLikeObject(feature_types.four, (1, 1), 10, 10, 100000, 1)
	left_area = sum_region(int_img, (1, 1), (24, 12))
	right_area = sum_region(int_img, (1, 12), (24, 24))
	left_area_3h = sum_region(int_img, (1, 1), (8, 24))
	middle_area_3h = sum_region(int_img, (8, 1), (16, 24))
	right_area_3h = sum_region(int_img, (16, 1), (24, 24))
	left_area_3v = sum_region(int_img, (1, 1), (24, 8))
	middle_area_3v = sum_region(int_img, (1, 8), (24, 16))
	right_area_3v = sum_region(int_img, (1, 16), (24, 24))
	top_left_area = sum_region(int_img, (1, 1), (12, 12))
	top_right_area = sum_region(int_img, (12, 1), (24, 12))
	bottom_left_area = sum_region(int_img, (1, 12), (12, 24))
	bottom_right_area = sum_region(int_img, (12, 12), (24, 24))
	expected_2v = feature_2v.threshold * feature_2v.polarity > left_area - right_area ? 1 : 0
	expected_2v_fail = feature_2h.threshold * -1 > left_area - right_area ? 1 : 0
	expected_2h = feature_2h.threshold * feature_2h.polarity > left_area - right_area ? 1 : 0
	expected_3h = feature_3h.threshold * feature_3h.polarity > left_area_3h - middle_area_3h + right_area_3h ? 1 : 0
	expected_3v = feature_3v.threshold * feature_3v.polarity > left_area_3v - middle_area_3v + right_area_3v ? 1 : 0
	expected_4 = feature_4.threshold * feature_4.polarity > top_left_area - top_right_area - bottom_left_area + bottom_right_area ? 1 : 0
	classifiers = []
	features = []
	p, n = 0, 0
	random_img = load_image(rand(vcat(filtered_ls.([pos_training_path, neg_training_path, pos_testing_path, neg_testing_path])...)))

    # IntegralImage.jl
    @test isequal(to_integral_image([17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; 11 18 25 2 9]), [17 41 42 50 65; 40 69 77 99 130; 44 79 100 142 195; 54 101 141 204 260; 65 130 195 260 325])
    @test isequal(sum_region(to_integral_image([1 7 4 2 9; 7 2 3 8 2; 1 8 7 9 1; 3 2 3 1 5; 2 9 5 6 6]), (4,4), (5,5)), 18)
    @test typeof(sum_region(to_integral_image([1 7 4 2 9.9; 7 2 3 8 2; 1 8 7 9 1; 3 2 3 1 5; 2 9 5 6 6]), (4,4), (5,5))) <: AbstractFloat
    @test sum_region(to_integral_image([1 7 4 2 9.9; 7 2 3 8 2; 1 8 7 9 1; 3 2 3 1 5; 2 9 5 6 6]), (4,4), (5,5)) isa AbstractFloat

    # HaarLikeFeature.jl
    @test HaarLikeObject(a, b, c, d, e, f) isa HaarLikeObject
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).feature_type isa Tuple{Integer, Integer}
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).position isa Tuple{Integer, Integer}
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).top_left isa Tuple{Integer, Integer}
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).bottom_right isa Tuple{Integer, Integer}
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).width isa Integer
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).height isa Integer
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).threshold ∈ [0, 1]
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).polarity ∈ [0, 1]
    @test HaarLikeObject((1,3), (1,3), 10, 8, 0, 1).weight ∈ [0, 1]
    @test get_vote(HaarLikeObject(a, b, c, d, e, f), arr) ∈ [-1, 1]
	@test get_vote(feature_2v, int_img) == expected_2v
	@test get_vote(feature_2v, int_img) != expected_2v_fail
	@test get_vote(feature_2h, int_img) == expected_2h
	@test get_vote(feature_3h, int_img) == expected_3h
	@test get_vote(feature_3v, int_img) == expected_3v
	@test get_vote(feature_4, int_img) == expected_4

    # AdaBoost.jl
	@suppress begin
	    classifiers = learn(pos_training_path, neg_training_path, 10, 8, 10, 8, 10)
		features = FaceDetection.create_features(19, 19, 8, 10, 8, 10)
	end
	@test length(features) == 4520

    # Utils.jl
	@test determine_feature_size(pos_training_path, neg_training_path) == (10, 10, 8, 8, (19, 19))
	@test get_faceness(classifiers[rand(1:length(classifiers))], random_img) isa Real
	@suppress begin
		num_faces = length(filtered_ls(pos_testing_path))
		num_non_faces = length(filtered_ls(neg_testing_path))
		p = sum(ensemble_vote_all(pos_testing_path, classifiers)) / num_faces
		n = (num_non_faces - sum(ensemble_vote_all(neg_testing_path, classifiers))) / num_non_faces
	end
	@test isapprox(p, 0.496, atol=1e-1) # these tests implicitly test the whole algorithm
	@test isapprox(n, 0.536, atol=1e-1)
end # end tests


text:
function learn_fitness(cont_f::Function, map_f::Function, meta_f::Function;
                       neps::Int64=6)
    steps = Array{Int64}(undef, neps)
    max_step = -Inf
    for i in 1:neps
        e = Episode(Grid(map_f()); reward=food_reward, meta=meta_f())
        run!(e, cont_f)
        if terminate(e)
            if i == 1
                return Float64(-neps)
            else
                return Float64(i - neps - steps[i-1] / e.meta["max_step"])
            end
        end
        steps[i] = e.meta["step"]
        max_step = e.meta["max_step"]
    end
    (steps[1] - steps[neps]) / max_step
end


text:
using POMDPs
using POMCPOWRAVE
using ProfileView
using POMDPModels

#=
using Gallium
breakpoint(Pkg.dir("POMCPOWRAVE", "src", "solver.jl"), 40)
=#

solver = POMCPOWRAVESolver(tree_queries=50_000,
                     eps=0.01,
                     c=10.0,
                     enable_action_pw=false,
                     alpha_observation=1/8,
                     rng=MersenneTwister(2))

problem = LightDark1D()
policy = solve(solver, problem)
ib = initial_state_distribution(problem)
a = action(policy, ib)

@time a = action(policy, ib)

Profile.clear()
@profile a = action(policy, ib)
ProfileView.view()


text:
# using Printf

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKdegenerate{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT
    cache::IntegratorCacheVPRK{DT,D,S}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersVPRKdegenerate(equ::ET, tab::TableauVPRK{TT}, Δt::TT,
                                  cache::IntegratorCacheVPRK{DT,D,S}) where {DT, TT, ET <: IODE{DT,TT}, D, S}

    @assert D == equ.d
    @assert S == tab.s

    q = zeros(DT,D)
    p = zeros(DT,D)

    ParametersVPRKdegenerate{DT, TT, ET, D, S}(equ, tab, Δt, cache, 0, q, p)
end

"Variational partitioned Runge-Kutta integrator."
mutable struct IntegratorVPRKdegenerate{DT, TT,
                SPT <: ParametersVPRK{DT,TT},
                PPT <: ParametersVPRKdegenerate{DT,TT},
                SST <: NonlinearSolver{DT},
                STP <: NonlinearSolver{DT},
                IT  <: InitialGuessPODE{DT,TT}, D, S} <: AbstractIntegratorVPRK{DT,TT}

    sparams::SPT
    pparams::PPT

    solver::SST
    projector::STP
    iguess::IT

    cache::IntegratorCacheVPRK{DT,D,S}
end

function IntegratorVPRKdegenerate(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create solver params
    sparams = ParametersVPRK(equation, tableau, Δt)

    # create projector params
    pparams = ParametersVPRKdegenerate(equation, tableau, Δt, cache)

    # create nonlinear solver
    solver = create_nonlinear_solver(DT, D*S, sparams)

    # create projector
    projector = create_nonlinear_solver(DT, D, pparams)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheVPRK{DT,D,S}(true)

    # create integrator
    IntegratorVPRKdegenerate{DT, TT, typeof(sparams), typeof(pparams), typeof(solver), typeof(projector), typeof(iguess), D, S}(
                sparams, pparams, solver, projector, iguess, cache, q, p)
end

"Compute solution of degenerate symplectic partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKdegenerate{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    # create temporary vectors
    q̅ = zeros(ST,D)
    p̅ = zeros(ST,D)

    quote
        @assert length(x) == length(b)

        compute_solution!(x, $q̅, $p̅, params)

        # compute b = - [q̅-q-BV]
        for k in 1:div(D,2)
            b[k] = params.q[k] - $q̅[k]
            for i in 1:S
                b[k] += params.Δt * params.tab.q.b[i] * params.cache.V[k,i]
            end
        end

        # compute b = - [p̅-p-BF]
        for k in 1:div(D,2)
            b[div(D,2)+k] = params.p[k] - $p̅[k]
            for i in 1:S
                b[div(D,2)+k] += params.Δt * params.tab.p.b[i] * params.cache.F[k,i]
            end
        end
    end
end

function compute_solution!(
                x::Vector{ST}, q̅::Vector{ST}, p̅::Vector{ST},
                params::ParametersVPRKdegenerate{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    @assert length(x) == length(q̅) == length(p̅)

    # copy x to q
    q̅ .= x

    # compute p̅ = ϑ(q)
    params.equ.ϑ(params.t + params.Δt, q̅, p̅)

end

function initial_guess!(int::IntegratorVPRKdegenerate, sol::AtomicSolutionPODE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).q.c[i])
        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = int.cache.ṽ[k]
        end
    end
end

function initial_guess_projection!(int::IntegratorVPRKdegenerate, sol::AtomicSolutionPODE)
    for k in eachdim(int)
        int.projector.x[k] = sol.q[k]
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKdegenerate{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.sparams, sol)
    update_params!(int.pparams, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.P, int.cache.F, int.sparams)

    # compute unprojected solution
    update_solution!(int, sol)

    # set initial guess for projection
    initial_guess_projection!(int, sol)

    # call solution solver
    solve!(int.projector)

    # print solver status
    print_solver_status(int.projector.status, int.projector.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.projector.status, int.projector.params)

    # compute projection vector fields
    compute_projection!(int.projector.x, sol.q, sol.p, int.pparams)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
"""
    diffusion_time(f)

Diffusion time of time profile `f`.
"""
function diffusion_time end

diffusion_time(f::PGSE) = f.Δ - f.δ / 3
diffusion_time(f::DoublePGSE) = 2(f.Δ - f.δ / 3)
diffusion_time(f::CosOGSE) = 1 / 8 * f.δ / f.nperiod
diffusion_time(f::SinOGSE) = 3 / 8 * f.δ / f.nperiod


text:
import Statistics

function dm_vals(e, closest, max_f_vals, min_f_vals)
	results =  Array{Float64, 2}(undef, size(closest))
	diff_vals = abs.(e .- closest)./(max_f_vals .- min_f_vals .+ eps(Float64))
	for i = 1:length(e)
		results[:,i] = Statistics.mean(diff_vals[:, 1:end .!= i], dims=2)
	end

	return results
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
using SymPy
using LinearAlgebra
λ = μ = 1.0

function gen_code(u)
    s = replace(replace(replace(sympy.julia_code(simplify(u)), ".*"=>"*"), ".^"=>"^"), "./"=>"/")
    "  (x,y)->begin;$s;end,\n"
end

function gen_code(u, v, g, h)
    s1 = gen_code(u)
    s2 = gen_code(v)
    s3 = gen_code(g)
    s4 = gen_code(h)
    "(\n" * s1 * s2 * s3 * s4 * ")"
end

x, y = @vars x y
# u = x*(1-x)*y*(1-y)
# v = x*(1-x)*y*(1-y)

# u = x^2 * (1-x) * y^2 * (1-y^2)
# v = u

u = x^2 * (1-x) * y^2 * (1-y^2)
v = x*(1-x)*y*(1-y)

ux = diff(u, x)
uy = diff(u, y)
vx = diff(v, x)
vy = diff(v, y)

ε = [
    ux (uy+vx)/2
    (uy+vx)/2 vy
]

σ = 2μ * ε + λ*(ux + vy) * diagm(0=>ones(2))

g = diff(σ[1,1], x) + diff(σ[1,2], y)
h = diff(σ[2,1], x) + diff(σ[2,2], y)

println(gen_code(u, v, g, h))


text:
## fit logistic regression

using DataFrames, GLM, Plots
temp = [53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81]
failure = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0]

data = DataFrame(temp = temp, failure = failure)

logit_fit = glm(@formula(failure ~ temp), data, Binomial(), LogitLink())
#plot(temp, predict(logit_fit), legend = false, xlabel = "Temperature", ylab = "Probability")
#scatter!(temp, predict(logit_fit))

## metropolis-hastings
using Distributions
γ = 0.57721
function ll(α::Float64, β::Float64)
    a = exp.(α .+ β*temp)
    return prod( (a ./ (1 .+ a) ).^failure .* (1 ./ (1 .+ a)).^(1 .- failure) )
end
function mh_logit(T::Int, α_hat::Float64, β_hat::Float64, σ_hat::Float64)
    φ = Normal(β_hat, σ_hat)
    π = Exponential(exp(α_hat+γ))
    Α = ones(T)
    Β = ones(T)
    for t = 1:T-1
        α = log(rand(π))
        β = rand(φ)
        r = ( ll(α, β) / ll(Α[t], Β[t]) ) * ( pdf(φ, Β[t]) / pdf(φ, β) )
        if rand() < r
            Α[t+1] = α
            Β[t+1] = β
        else
            Α[t+1] = Α[t]
            Β[t+1] = Β[t]
        end
    end
    return Α, Β
end

# trace plot

Α, Β = mh_logit(10000, 15.04, -0.233, 0.108)

p1 = plot(Α, legend = false, xlab = "Intercept")
hline!([15.04])

p2 = plot(Β, legend = false, xlab = "Slope")
hline!([-0.233])

plot(p1, p2, layout = (1,2))

# mean trace plot

Αmean = cumsum(Α) ./ collect(1:length(Α))
Βmean = cumsum(Β) ./ collect(1:length(Β))

p1 = plot(Αmean, legend = false, xlab = "Intercept")
hline!([15.04])

p2 = plot(Βmean, legend = false, xlab = "Slope")
hline!([-0.233])

plot(p1, p2, layout = (1,2))

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
"""
    RMSEA(sem_fit::SemFit)

Return the RMSEA.
"""
function RMSEA end

RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) =
    RMSEA(df(sem_fit), χ²(sem_fit), n_obs(sem_fit))

RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}) =
    sqrt(length(sem_fit.model.sems))*RMSEA(df(sem_fit), χ²(sem_fit), n_obs(sem_fit))

function RMSEA(df, chi2, n_obs)
    rmsea = (chi2 - df) / (n_obs*df)
    rmsea > 0 ? nothing : rmsea = 0
    return sqrt(rmsea)
end


text:
#derivative logic

#ForwardDiff compiles one method per function, so separating the
#differentiation logic from the property logic allows the differentials
#to be compiled only once

"""
    ∂f∂T(model,V,T,z=SA[1.0])

returns `f` and `∂f/∂T` at constant total volume and composition, where f is the total helmholtz energy, given by `eos(model,V,T,z)`

"""
function ∂f∂T(model,V,T,z=SA[1.0])
    return ForwardDiff.derivative(∂T -> eos(model,V,∂T,z),T)
end

"""
    ∂f∂V(model,V,T,z=SA[1.0])

returns `f` and `∂f/∂V` at constant temperature and composition, where f is the total helmholtz energy, given by `eos(model,V,T,z)`, and V is the total volume

"""
function ∂f∂V(model,V,T,z)
    return ForwardDiff.derivative(∂V -> eos(model,∂V,T,z),V)
end

#returns a tuple of the form ([∂f∂V,∂f∂T],f),using the least amount of computation
"""
    ∂f(model,V,T,z)

returns zeroth order (value) and first order derivative information of the total helmholtz energy (given by `eos(model,V,T,z)`).
the result is given in two values:

```julia
grad_f,fval = ∂2f(model,V,T,z)
```

where:

```julia
fval   = f(V,T) = eos(model,V,T,z)

grad_f = [ ∂f/∂V; ∂f/∂T]

```

Where `V` is the total volume, `T` is the temperature and `f` is the total helmholtz energy.
"""
function ∂f(model,V,T,z)
    f(w) = eos(model,first(w),last(w),z)
    V,T = promote(V,T)
    VT_vec = SVector(V,T)
    ∂result = DiffResults.GradientResult(VT_vec)
    res_∂f =  ForwardDiff.gradient!(∂result, f,VT_vec)
    _f =  DiffResults.value(res_∂f)
    _∂f = DiffResults.gradient(res_∂f)
    return (_∂f,_f)
end

#returns p and ∂p∂V at constant T
#it doesnt do a pass over temperature, so its
#faster that d2f when only requiring d2fdV2

"""
    p∂p∂V(model,V,T,z=SA[1.0])

returns `p` and `∂p/∂V` at constant temperature, where p is the pressure = `pressure(model,V,T,z)` and `V` is the total Volume.

"""
function p∂p∂V(model,V,T,z=SA[1.0])
    V_vec =   SVector(V)
    f(∂V) = pressure(model,only(∂V),T,z)
    ∂result = DiffResults.GradientResult(V_vec)
    res_∂f =  ForwardDiff.gradient!(∂result, f,V_vec)
    _p =  DiffResults.value(res_∂f)
    _∂p∂V = only(DiffResults.gradient(res_∂f))
    return _p,_∂p∂V
end

"""
    ∂2f(model,V,T,z)

returns zeroth order (value), first order and second order derivative information of the total helmholtz energy (given by `eos(model,V,T,z)`).
the result is given in three values:

```
hess_f,grad_f,fval = ∂2f(model,V,T,z)
```

where:
```
fval   = f(V,T) = eos(model,V,T,z)

grad_f = [ ∂f/∂V; ∂f/∂T]

hess_f = [ ∂²f/∂V²; ∂²f/∂V∂T
          ∂²f/∂V∂T; ∂²f/∂V²]
 ```

Where `V` is the total volume, `T` is the temperature and `f` is the total helmholtz energy.
"""
function ∂2f(model,V,T,z)
    f(w) = eos(model,first(w),last(w),z)
    V,T = promote(V,T)
    VT_vec =   SVector(V,T)
    ∂result = DiffResults.HessianResult(VT_vec)
    res_∂f =  ForwardDiff.hessian!(∂result, f,VT_vec)
    _f =  DiffResults.value(res_∂f)
    _∂f = DiffResults.gradient(res_∂f)
    _∂2f = DiffResults.hessian(res_∂f)
    return (_∂2f,_∂f,_f)
end

"""
    ∂2p(model,V,T,z)

returns zeroth order (value), first order and second order derivative information of the pressure.
the result is given in three values:

```
hess_p,grad_p,pval = ∂2p(model,V,T,z)
```

where:
```
pval   = p(V,T) = pressure(model,V,T,z)

grad_p = [ ∂p/∂V; ∂p/∂T]

hess_p = [ ∂²p/∂V²; ∂²p/∂V∂T
          ∂²p/∂V∂T; ∂²p/∂V²]
 ```

Where `V` is the total volume, `T` is the temperature and `p` is the pressure.
"""
function ∂2p(model,V,T,z)
    f(w) = pressure(model,first(w),last(w),z)
    V,T = promote(V,T)
    VT_vec =   SVector(V,T)
    ∂result = DiffResults.HessianResult(VT_vec)
    res_∂f =  ForwardDiff.hessian!(∂result, f,VT_vec)
    _f =  DiffResults.value(res_∂f)
    _∂f = DiffResults.gradient(res_∂f)
    _∂2f = DiffResults.hessian(res_∂f)
    return (_∂2f,_∂f,_f)
end

"""
    f_hess(model,V,T,z)

returns the second order volume (`V`) and temperature (`T`) derivatives of the total helmholtz energy (given by `eos(model,V,T,z)`). the result is given in a 2x2 `SMatrix`, in the form:

```
[ ∂²f/∂V²  ∂²f/∂V∂T
 ∂²f/∂V∂T  ∂²f/∂T²]
 ```

use this instead of the ∂2f if you only need second order information. ∂2f also gives zeroth and first order derivative information, but due to a bug in the used AD, it allocates more than necessary.
"""
function f_hess(model,V,T,z)
    f(w) = eos(model,first(w),last(w),z)
    V,T = promote(V,T)
    VT_vec = SVector(V,T)
    return ForwardDiff.hessian(f,VT_vec)
end

"""
    ∂²³f(model,V,T,z=SA[1.0])

returns `∂²A/∂V²` and `∂³A/∂V³`, in a single ForwardDiff pass. used mainly in `crit_pure` objective function

"""
function ∂²³f(model,V,T,z=SA[1.0])
    V_vec =   SVector(V)
    f(∂A∂V) = ForwardDiff.derivative(∂²A∂V² -> pressure(model,only(∂²A∂V²),T,z),only(∂A∂V))
    ∂result = DiffResults.GradientResult(V_vec)
    res_∂f =  ForwardDiff.gradient!(∂result,f,V_vec)
    _∂²A∂V² =  DiffResults.value(res_∂f)
    _∂³A∂V³ = only(DiffResults.gradient(res_∂f))
    return _∂²A∂V², _∂³A∂V³
end


text:
################################################################################
## Linear Regression
##   y ~ N(b0 + b1 * x, s2)
##   b0, b1 ~ N(0, 1000)
##   s2 ~ invgamma(0.001, 0.001)
################################################################################

using Mamba

## Data
data = Dict(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)

## Log-transformed Posterior(b0, b1, log(s2)) + Constant
logf = function(x::DenseVector)
   b0 = x[1]
   b1 = x[2]
   logs2 = x[3]
   r = data[:y] - b0 - b1 * data[:x]
   (-0.5 * length(data[:y]) - 0.001) * logs2 -
     (0.5 * dot(r, r) + 0.001) / exp(logs2) -
     0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
end

## MCMC Simulation with Slice Sampling
## With multivariate (1) and univariate (2) updating
n = 5000
sim1 = Chains(n, 3, names = ["b0", "b1", "s2"])
sim2 = Chains(n, 3, names = ["b0", "b1", "s2"])
width = [1.0, 1.0, 2.0]
theta1 = SliceUnivariate([0.0, 0.0, 0.0], width, logf)
theta2 = SliceMultivariate([0.0, 0.0, 0.0], width, logf)
for i in 1:n
  sample!(theta1)
  sample!(theta2)
  sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])]
  sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])]
end
describe(sim1)
describe(sim2)


text:
include("../../number_theo/src/number_module.jl")

module quasirandom
using LinearAlgebra
using ..number

function halton(i; base = 2)
    # von der Corput sequence
    result = 0.0
    f = 1.0
    while i > 0
        f = f / base
        result += f * mod(i, base)
        i = floor(Int, i / base)
    end
    return result
end

function halton_filling_seq(n_points; base = 2)
    return [halton(i, base = base) for i = 1:n_points]
end

function halton_filling_seq(n_points, n_dim)
    upto = max(n_dim * (log(n_dim) + log(log(n_dim))), 6) |> ceil |> Int
    bases = number.sieve_of_eratosthenes(upto)
    mat = zeros(n_points, n_dim)
    for dd = 1:n_dim
        mat[:,dd] = halton_filling_seq(n_points; base = bases[dd])
    end
    return mat
end

function pairwise_distances(X)
    # columns are the different dimensions
    # rows are the different points

    (N, n_dim)  = size(X)
    Nn          = div(N * (N-1), 2)
    n           = zeros(Nn)
    kk = 0
    @inbounds for ii = 1:(N-1)
        @inbounds for jj = (ii+1):N
            kk += 1
            v1 = view(X, ii, 1:n_dim)
            v2 = view(X, jj, 1:n_dim)
            n[kk] = norm(v1 - v2)
        end
    end
    return n
end

function morris_mitchell_criterion(X, q = 1.0)
    dists = pairwise_distances(X)
    s = sum( dists.^(-q) )^(1/q)
    return s
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text: maxmultiple(divisor, bound) = (bound÷divisor)*divisor

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
    sx = fitslopex(fx; kmin=1)

Return spectral slope (minus 1) of spatial spectrum `fx`.
"""
	function fitslopex(fx; kmin=1)
        kmin = Int(kmin);
        xr = collect(1:length(fx));
        x = log.(xr[kmin:end])
        y = log.(fx[kmin:end])
        sx, a = [ x[:] ones(length(x))] \ y
	sx = abs(sx)
	sx = sx-1
	return sx
	end


text:
@doc raw"""
    PrimalDualProblem {mT <: Manifold, nT <: Manifold} <: PrimalDualProblem} <: Problem

Describes a Problem for the linearized Chambolle-Pock algorithm.

# Fields

* `M`, `N` – two manifolds $\mathcal M$, $\mathcal N$
* `cost` $F + G(Λ(⋅))$ to evaluate interims cost function values
* `forward_oprator` the operator for the forward operation in the algorthm, either $Λ$ (exact) or $DΛ$ (linearized).
* `linearized_adjoint_operator` The adjoint differential $(DΛ)^* \colon \mathcal N \to T\mathcal M$
* `prox_F` the proximal map belonging to $f$
* `prox_G_dual` the proximal map belonging to $g_n^*$
* `Λ` – (`fordward_operator`) for the linearized variant, this has to be set to the exact forward operator.
  This operator is required in several variants of the linearized algorithm.
  Since the exact variant is the default, `Λ` is by default set to `forward_operator`.

# Constructor

    LinearizedPrimalDualProblem(M, N, cost, prox_F, prox_G_dual, forward_operator, adjoint_linearized_operator,Λ=forward_operator)

"""
mutable struct PrimalDualProblem{mT<:Manifold,nT<:Manifold} <: Problem
    M::mT
    N::nT
    cost::Function
    prox_F::Function
    prox_G_dual::Function
    forward_operator::Function
    adjoint_linearized_operator::Function
    Λ::Union{Function,Missing}
end
function PrimalDualProblem(
    M::mT,
    N::nT,
    cost,
    prox_F,
    prox_G_dual,
    forward_operator,
    adjoint_linearized_operator,
    Λ=forward_operator,
) where {mT<:Manifold,nT<:Manifold}
    return PrimalDualProblem{mT,nT}(
        M, N, cost, prox_F, prox_G_dual, forward_operator, adjoint_linearized_operator, Λ
    )
end

@doc raw"""
    PrimalDualOptions

A general type for all primal dual based options to be used within primal dual
based algorithms
"""
abstract type PrimalDualOptions <: Options end

@doc raw"""
    ChambollePockOptions <: PrimalDualOptions

stores all options and variables within a linearized or exact Chambolle Pock.
The following list provides the order for the constructor, where the previous iterates are
initialized automatically and values with a default may be left out.

* `m` - base point on $ \mathcal M $
* `n` - base point on $ \mathcal N $
* `x` - an initial point on $x^{(0)} \in \mathcal M$ (and its previous iterate)
* `ξ` - an initial tangent vector $\xi^{(0)}\in T^*\mathcal N$ (and its previous iterate)
* `xbar` - the relaxed iterate used in the next dual update step (when using `:primal` relaxation)
* `ξbar` - the relaxed iterate used in the next primal update step (when using `:dual` relaxation)
* `Θ` – factor to damp the helping $\tilde x$
* `primal_stepsize` – (`1/sqrt(8)`) proximal parameter of the primal prox
* `dual_stepsize` – (`1/sqrt(8)`) proximnal parameter of the dual prox
* `acceleration` – (`0.`) acceleration factor due to Chambolle & Pock
* `relaxation` – (`1.`) relaxation in the primal relaxation step (to compute `xbar`)
* `relax` – (`_primal`) which variable to relax (`:primal` or `:dual`)
* `stop` - a [`StoppingCriterion`](@ref)
* `type` – (`exact`) whether to perform an `:exact` or `:linearized` Chambolle-Pock
* `update_primal_base` (`(p,o,i) -> o.m`) function to update the primal base
* `update_dual_base` (`(p,o,i) -> o.n`) function to update the dual base
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use
* `inverse_retraction_method` - (`LogarithmicInverseRetraction()`) an inverse retraction to use.
* `vector_transport_method` - (`ParallelTransport()`) a vector transport to use

where for the last two the functions a [`Problem`](@ref)` p`,
[`Options`](@ref)` o` and the current iterate `i` are the arguments.
If you activate these to be different from the default identity, you have to provide
`p.Λ` for the algorithm to work (which might be `missing` in the linearized case).

# Constructor
    ChambollePockOptions(m::P, n::Q, x::P, ξ::T, primal_stepsize::Float64, dual_stepsize::Float64;
        acceleration::Float64 = 0.0,
        relaxation::Float64 = 1.0,
        relax::Symbol = :primal,
        stopping_criterion::StoppingCriterion = StopAfterIteration(300),
        variant::Symbol = :exact,
        update_primal_base::Union{Function,Missing} = missing,
        update_dual_base::Union{Function,Missing} = missing,
        retraction_method = ExponentialRetraction(),
        inverse_retraction_method = LogarithmicInverseRetraction(),
        vector_transport_method = ParallelTransport(),
    )
"""
mutable struct ChambollePockOptions{
    P,
    Q,
    T,
    RM<:AbstractRetractionMethod,
    IRM<:AbstractInverseRetractionMethod,
    VTM<:AbstractVectorTransportMethod,
} <: PrimalDualOptions
    m::P
    n::Q
    x::P
    xbar::P
    ξ::T
    ξbar::T
    primal_stepsize::Float64
    dual_stepsize::Float64
    acceleration::Float64
    relaxation::Float64
    relax::Symbol
    stop::StoppingCriterion
    variant::Symbol
    update_primal_base::Union{Function,Missing}
    update_dual_base::Union{Function,Missing}
    retraction_method::RM
    inverse_retraction_method::IRM
    vector_transport_method::VTM

    function ChambollePockOptions(
        m::P,
        n::Q,
        x::P,
        ξ::T,
        primal_stepsize::Float64=1 / sqrt(8),
        dual_stepsize::Float64=1 / sqrt(8);
        acceleration::Float64=0.0,
        relaxation::Float64=1.0,
        relax::Symbol=:primal,
        stopping_criterion::StoppingCriterion=StopAfterIteration(300),
        variant::Symbol=:exact,
        update_primal_base::Union{Function,Missing}=missing,
        update_dual_base::Union{Function,Missing}=missing,
        retraction_method::RM=ExponentialRetraction(),
        inverse_retraction_method::IRM=LogarithmicInverseRetraction(),
        vector_transport_method::VTM=ParallelTransport(),
    ) where {
        P,
        Q,
        T,
        RM<:AbstractRetractionMethod,
        IRM<:AbstractInverseRetractionMethod,
        VTM<:AbstractVectorTransportMethod,
    }
        return new{P,Q,T,RM,IRM,VTM}(
            m,
            n,
            x,
            deepcopy(x),
            ξ,
            deepcopy(ξ),
            primal_stepsize,
            dual_stepsize,
            acceleration,
            relaxation,
            relax,
            stopping_criterion,
            variant,
            update_primal_base,
            update_dual_base,
            retraction_method,
            inverse_retraction_method,
            vector_transport_method,
        )
    end
end
get_solver_result(o::ChambollePockOptions) = o.x
@doc raw"""
    primal_residual(p, o, x_old, ξ_old, n_old)

Compute the primal residual at current iterate $k$ given the necessary values $x_{k-1},
ξ_{k-1}, and $n_{k-1}$ from the previous iterate.
```math
\Bigl\lVert
\frac{1}{σ}\operatorname{retr}^{-1}_{x_{k}}x_{k-1} -
V_{x_k\gets m_k}\bigl(DΛ^*(m_k)\bigl[V_{n_k\gets n_{k-1}}ξ_{k-1} - ξ_k \bigr]
\Bigr\rVert
```
where $V_{\cdot\gets\cdot}$ is the vector transport used in the [`ChambollePockOptions`](@ref)
"""
function primal_residual(p::PrimalDualProblem, o::ChambollePockOptions, x_old, ξ_old, n_old)
    return norm(
        p.M,
        o.x,
        1 / o.primal_stepsize *
        inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method) -
        vector_transport_to(
            p.M,
            o.m,
            p.adjoint_linearized_operator(
                o.m,
                vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) -
                o.ξ,
            ),
            o.x,
            o.vector_transport_method,
        ),
    )
end
@doc raw"""
    dual_residual(p, o, x_old, ξ_old, n_old)

Compute the dual residual at current iterate $k$ given the necessary values $x_{k-1},
ξ_{k-1}, and $n_{k-1}$ from the previous iterate. The formula is slightly different depending
on the `o.variant` used:

For the `:lineaized` it reads
```math
\Bigl\lVert
\frac{1}{τ}\bigl(
V_{n_{k}\gets n_{k-1}}(ξ_{k-1})
- ξ_k
\bigr)
-
DΛ(m_k)\bigl[
V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1}
\bigr]
\Bigr\rVert
```

and for the `:exact` variant

```math
\Bigl\lVert
\frac{1}{τ} V_{n_{k}\gets n_{k-1}}(ξ_{k-1})
-
\operatorname{retr}^{-1}_{n_{k}}\bigl(
Λ(\operatorname{retr}_{m_{k}}(V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1}))
\bigr)
\Bigr\rVert
```

where in both cases $V_{\cdot\gets\cdot}$ is the vector transport used in the [`ChambollePockOptions`](@ref).
"""
function dual_residual(p::PrimalDualProblem, o::ChambollePockOptions, x_old, ξ_old, n_old)
    if o.variant === :linearized
        return norm(
            p.N,
            o.n,
            1 / o.dual_stepsize *
            (vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) - o.ξ) -
            p.forward_operator(
                o.m,
                vector_transport_to(
                    p.M,
                    o.x,
                    inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method),
                    o.m,
                    o.vector_transport_method,
                ),
            ),
        )
    elseif o.variant === :exact
        return norm(
            p.N,
            o.n,
            1 / o.dual_stepsize *
            (vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) - o.n) -
            inverse_retract(
                p.N,
                o.n,
                p.Λ(
                    retract(
                        p.M,
                        o.m,
                        vector_transport_to(
                            p.M,
                            o.x,
                            inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method),
                            o.m,
                            ParallelTransport(),
                        ),
                        o.retraction_method,
                    ),
                ),
                o.inverse_retraction_method,
            ),
        )
    else
        throw(
            DomainError(
                o.variant,
                "Unknown Chambolle–Pock variant, allowed are `:exact` or `:linearized`.",
            ),
        )
    end
end
#
# Special Debuggers
#
@doc raw"""
    DebugDualResidual <: DebugAction

A Debug action to print the dual residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugDualResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugDualResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "Dual Residual: ", a)
    end
    function DebugDualResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,T,Q}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "Dual Residual: ", a)
    end
end
function (d::DebugDualResidual)(p::PrimalDualProblem, o::ChambollePockOptions, i::Int)
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(d.io, d.prefix * string(dual_residual(p, o, xOld, ξOld, nOld)))
    end
    return d.storage(p, o, i)
end

@doc raw"""
    DebugPrimalResidual <: DebugAction

A Debug action to print the primal residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugPrimalResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugPrimalResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "Primal Residual: ", a)
    end
    function DebugPrimalResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,T,Q}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "Primal Residual: ", a)
    end
end
function (d::DebugPrimalResidual)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(d.io, d.prefix * string(primal_residual(p, o, xOld, ξOld, nOld)))
    end
    return d.storage(p, o, i)
end
@doc raw"""
    DebugPrimalDualResidual <: DebugAction

A Debug action to print the primaldual residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugPrimalDualResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugPrimalDualResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "PD Residual: ", a)
    end
    function DebugPrimalDualResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,Q,T}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "PD Residual: ", a)
    end
end
function (d::DebugPrimalDualResidual)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(
            d.io,
            d.prefix * string(
                (
                    primal_residual(p, o, xOld, ξOld, nOld) +
                    dual_residual(p, o, xOld, ξOld, nOld)
                ) / manifold_dimension(p.M),
            ),
        )
    end
    return d.storage(p, o, i)
end

#
# Debugs
#
"""
    DebugPrimalChange(opts...)

Print the change of the primal variable by using [`DebugChange`](@ref),
see their constructors for detail.
"""
DebugPrimalChange(opts...) = DebugChange(opts[1], "Primal Change: ", opts[2:end]...)

"""
    DebugPrimalIterate(opts...)

Print the change of the primal variable by using [`DebugIterate`](@ref),
see their constructors for detail.
"""
DebugPrimalIterate(opts...) = DebugIterate(opts...)

"""
    DebugDualIterate(e)

Print the dual variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.ξ`.
"""
DebugDualIterate(opts...) = DebugEntry(:ξ, "ξ:", opts...)

"""
    DebugDualChange(opts...)

Print the change of the dual variable, similar to [`DebugChange`](@ref),
see their constructors for detail, but with a different calculation of the change,
since the dual variable lives in (possibly different) tangent spaces.
"""
mutable struct DebugDualChange <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugDualChange(
        a::StoreOptionsAction=StoreOptionsAction((:ξ, :n)), io::IO=stdout
    )
        return new(io, "Dual Change: ", a)
    end
    function DebugDualChange(
        values::Tuple{T,P},
        a::StoreOptionsAction=StoreOptionsAction((:ξ, :n)),
        io::IO=stdout,
    ) where {P,T}
        update_storage!(a, Dict{Symbol,Any}(k => v for (k, v) in zip((:ξ, :n), values)))
        return new(io, "Dual Change: ", a)
    end
end
function (d::DebugDualChange)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:ξ, :n])) && i > 0 # all values stored
        ξOld, nOld = get_storage.(Ref(d.storage), [:ξ, :n]) #fetch
        print(
            d.io,
            d.prefix * string(
                norm(
                    p.N,
                    o.n,
                    vector_transport_to(p.N, nOld, ξOld, o.n, o.vector_transport_method) -
                    o.ξ,
                ),
            ),
        )
    end
    return d.storage(p, o, i)
end

"""
    DebugDualBaseIterate(io::IO=stdout)

Print the dual base variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.n`.
"""
DebugDualBaseIterate(io::IO=stdout) = DebugEntry(:n, "n:", io)
"""
    DebugDualChange(a=StoreOptionsAction((:ξ)),io::IO=stdout)

Print the change of the dual base variable by using [`DebugEntryChange`](@ref),
see their constructors for detail, on `o.n`.
"""
function DebugDualBaseChange(a::StoreOptionsAction=StoreOptionsAction((:n)), io::IO=stdout)
    return DebugEntryChange(
        :n, (p, o, x, y) -> distance(p.N, x, y), a, "Dual Base Change:", io
    )
end

"""
    DebugPrimalBaseIterate(io::IO=stdout)

Print the primal base variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.m`.
"""
DebugPrimalBaseIterate(io::IO=stdout) = DebugEntry(:m, "m:", io)

"""
    DebugPrimalBaseChange(a::StoreOptionsAction=StoreOptionsAction((:m)),io::IO=stdout)

Print the change of the primal base variable by using [`DebugEntryChange`](@ref),
see their constructors for detail, on `o.n`.
"""
function DebugPrimalBaseChange(
    a::StoreOptionsAction=StoreOptionsAction((:m)), io::IO=stdout
)
    return DebugEntryChange(
        :m, (p, o, x, y) -> distance(p.M, x, y), a, "Primal Base Change:", io
    )
end

#
# Records
#

# Primals are just the entries
"""
    RecordPrimalChange(a)

Create an [`RecordAction`](@ref) that records the primal value change,
i.e. [`RecordChange`](@ref), since we just redord the change of `o.x`.
"""
RecordPrimalChange() = RecordChange()

"""
    RecordDualBaseIterate(x)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordIterate`](@ref), i.e. `o.x`.
"""
RecordPrimalIterate(x) = RecordIterate(x)

"""
    RecordDualIterate(ξ)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordEntry`](@ref) of `o.ξ`, so .
"""
RecordDualIterate(ξ) = RecordEntry(ξ, :ξ)

"""
    RecordDualChange()

Create the action either with a given (shared) Storage, which can be set to the
`values` Tuple, if that is provided).
"""
RecordDualChange() = RecordEntryChange(:ξ, (p, o, x, y) -> distance(p.N, x, y))

"""
    RecordDualBaseIterate(n)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordEntry`](@ref) of `o.n`.
"""
RecordDualBaseIterate(n) = RecordEntry(n, :n)

"""
    RecordDualBaseChange(e)

Create an [`RecordAction`](@ref) that records the dual base point change,
i.e. [`RecordEntryChange`](@ref) of `o.n` with distance to the last value to store a value.
"""
RecordDualBaseChange() = RecordEntryChange(:n, (p, o, x, y) -> distance(p.N, x, y))

"""
    RecordPrimalBaseIterate(x)

Create an [`RecordAction`](@ref) that records the primal base point,
i.e. [`RecordEntry`](@ref) of `o.m`.
"""
RecordPrimalBaseIterate(m) = RecordEntry(m, :m)
"""
    RecordPrimalBaseChange()

Create an [`RecordAction`](@ref) that records the primal base point change,
i.e. [`RecordEntryChange`](@ref) of `o.m` with distance to the last value to store a value.
"""
RecordPrimalBaseChange() = RecordEntryChange(:m, (p, o, x, y) -> distance(p.M, x, y))


text:
@testset "flag" begin
    @testset "basic" begin
        @system SFlag(Controller) begin
            a => true ~ flag
            b => false ~ flag
        end
        s = instance(SFlag)
        @test s.a' == true && s.b' == false
    end

    @testset "flag vs track" begin
        @system SFlagVsTrack(Controller) begin
            a => 1 ~ accumulate
            f1(a) => a >= 1 ~ flag
            f2(a) => a >= 1 ~ track::Bool
            x1(f1) => (f1 ? 1 : 0) ~ track
            x2(f2) => (f2 ? 1 : 0) ~ track
        end
        s = instance(SFlagVsTrack)
        @test s.a' == 0
        @test s.f1' == false && s.f2' == false
        @test s.x1' == 0 && s.x2' == 0
        update!(s)
        @test s.a' == 1
        @test s.f1' == true && s.f2' == true
        @test s.x1' == 1 && s.x2' == 1
        update!(s)
        @test s.a' == 2
        @test s.f1' == true && s.f2' == true
        @test s.x1' == 1 && s.x2' == 1
    end

    @testset "flag parameter" begin
        @system SFlagParameter(Controller) begin
            a => true ~ flag(parameter)
            b => false ~ flag(parameter)
        end
        s1 = instance(SFlagParameter)
        @test s1.a' == true && s1.b' == false
        c = :SFlagParameter => (a = false, b = true)
        s2 = instance(SFlagParameter; config=c)
        @test s2.a' == false && s2.b' == true
    end
end


text:
# utilities for x ↔ scale factor ↔ redshift
a2z(a::T) where T = one(T)/a - one(T)
z2a(z::T) where T = one(T)/(z + one(T))
a2x(a) = log(a)
x2a(x) = exp(x)
z2x(z) = a2x(z2a(z))
x2z(x) = a2z(x2a(x))

# utility function for constructing an interpolator
spline(f, x_grid) = scale(interpolate(f, BSpline(Cubic(Line(OnGrid())))), x_grid)
spline_∂ₓ(f, x_grid) = spline([Interpolations.gradient(f, x)[1] for x in x_grid], x_grid)
spline_∂ₓ²(f, x_grid) = spline([Interpolations.hessian(f, x)[1] for x in x_grid], x_grid)

δ_kron(i, j) = (i == j) ? 1 : 0

#utilities for mapping between comoving momenta and unit interval
to_ui(lq,lqmi,lqma) = -1 + (1- (-1)) / (lqma-lqmi) * (lq-lqmi)
from_ui(x,lqmi,lqma) = lqmi + (lqma- lqmi) / (1- (-1)) * (x- (-1))
dxdq(q,logqmin,logqmax) = (1+to_ui(1+logqmin,logqmin,logqmax))/(q*log(10))
xq2q(x,logqmin,logqmax) = 10.0 ^ from_ui(x,logqmin,logqmax)

# utilities for applying an FFTLog transformation ===

struct FFTLogPlan{T, OT, AA<:AbstractArray{Complex{T},1},
                  AAR<:AbstractArray{T,1}, PT<:Plan, IPT<:Plan}
    L::T
    N::Int
    μ::OT
    q::T
    r₀::T
    k₀r₀::T
    uₘ::AA
    r::AAR
    k::AAR
    fftplan!::PT
    ifftplan!::IPT
end

function plan_fftlog(r::AA, μ, q, k₀r₀=1.0;
                     kropt=true) where {T, AA<:AbstractArray{T,1}}
    logrmin = log(first(r))
    logrmax = log(last(r))
    r₀ = exp((logrmin + logrmax)/2)
    @assert logrmin < logrmax
    N = length(r)
    L = logrmax - logrmin
    dlnr = L / (N - 1)
    if kropt
        k₀r₀ = k₀r₀_low_ringing(N, μ, q, L, k₀r₀)
    end
    k₀ = k₀r₀ / r₀
    Nhalf = N ÷ 2
    n = range(-Nhalf,Nhalf,length=N)
    k = reverse(k₀ .* exp.(n .* L / N))

    m = fftfreq(N, N)  # get indicies that go from [-N/2] to [N/2]
    uₘ_coeff = similar(r, Complex{T})
    for i in eachindex(m)
        uₘ_coeff[i] = uₘ(m[i], μ, q, dlnr, k₀r₀, N)
    end
    uₘ_coeff[N÷2+1] = real(uₘ_coeff[N÷2+1])  # eq 19
    fftplan! = plan_fft!(uₘ_coeff)
    ifftplan! = plan_ifft!(uₘ_coeff)
    return FFTLogPlan(L, N, μ, q, r₀, k₀r₀, uₘ_coeff, r, k, fftplan!, ifftplan!)
end

U_μ(μ, x) = exp(x * log(2.) - lgamma(0.5 * (μ + 1 - x)) + lgamma(0.5 * (μ + 1 + x)))
uₘ(m, μ, q, dlnr, k₀r₀, N) = (k₀r₀)^(-2π * im * m / (dlnr * N)) * U_μ(μ, q + 2π * im * m / (dlnr * N) )

function k₀r₀_low_ringing(N, μ, q, L, k₀r₀=1.0)
    # from pyfftlog
    dlnr = L / (N-1)
    xp = (μ + 1 + q) / 2
    xm = (μ + 1 - q) / 2
    y = π * im / 2 / dlnr
    zp = lgamma(xp + y)
    zm = lgamma(xm + y)
    arg = log(2 / k₀r₀) / dlnr + imag(zp + zm) / π
    return k₀r₀ * exp((arg - round(arg))* dlnr)
end

function mul!(Y, pl::FFTLogPlan, A)
    Y .= A
    Y .*= (pl.r).^(-pl.q)
    pl.fftplan! * Y
    Y .*= pl.uₘ
    pl.ifftplan! * Y
    Y .*= (pl.r).^(pl.q)
end

function ldiv!(Y, pl::FFTLogPlan, A)
    Y .= A
    Y .*= (pl.r).^(-pl.q)
    pl.fftplan! * Y
    Y ./= pl.uₘ
    pl.ifftplan! * Y
    Y .*= (pl.r).^(pl.q)
end


text:
using JuPOT

name = ARGS[1]
age = ARGS[2]
cats = ARGS[3]

message = "Hello after still working " * name * " you are " * age * " and own " * cats * " cats!!"

outfile = "output.txt"
f = open(outfile, "w")
println(f, message)
close(f)


text:
module GMG
export V_cycle!
function jacobi_relaxation!(level::Int,nx::Int,ny::Int,u::Array{Float64,2},f::Array{Float64,2},iters::Int,pre::Bool)

    dx=1.0/nx
    dy=1.0/ny
    Ax=1/dx^2;Ay=1/dy^2
    Ap=1.0/(2.0*(Ax+Ay))

    #Dirichlet BC
    u[  1,:]=-u[    2,:]
    u[end,:]=-u[end-1,:]
    u[  :,1]=-u[    :,2]
    u[:,end]=-u[:,end-1]

  #if it is a pre-sweep, u is fully zero (on the finest grid depends on BC, always true on coarse grids)
  # we can save some calculation, if doing only one iteration, which is typically the case.
    if ((pre) && (level>1))
        u[2:nx+1,2:ny+1] = -Ap*f[2:nx+1,2:ny+1]

        #Dirichlet BC
        u[  1,:]=-u[    2,:]
        u[end,:]=-u[end-1,:]
        u[  :,1]=-u[    :,2]
        u[:,end]=-u[:,end-1]
        iters=iters-1
    end

    for it = 1:iters
        u[2:nx+1,2:ny+1] = Ap*(Ax*(u[3:nx+2,2:ny+1] + u[1:nx,2:ny+1])
                             + Ay*(u[2:nx+1,3:ny+2] + u[2:nx+1,1:ny])
                             - f[2:nx+1,2:ny+1])
        #Dirichlet BC
        u[  1,:]=-u[    2,:]
        u[end,:]=-u[end-1,:]
        u[  :,1]=-u[    :,2]
        u[:,end]=-u[:,end-1]
    end
    #  if(not pre):
    #    return u,None

    res=zeros(Float64,nx+2,ny+2)
    res[2:nx+1,2:ny+1]=f[2:nx+1,2:ny+1]-(( Ax*(u[3:nx+2,2:ny+1]+u[1:nx,2:ny+1])
                                         + Ay*(u[2:nx+1,3:ny+2]+u[2:nx+1,1:ny])
                                         - 2.0*(Ax+Ay)*u[2:nx+1,2:ny+1]))
    return u,res
end

function restrict(nx,ny,v)
  """
  restrict 'v' to the coarser grid
  """
  v_c=zeros(Float64,nx+2,ny+2)

  v_c[2:nx+1,2:ny+1]=0.25*(v[2:2:2*nx,2:2:2*ny]+v[2:2:2*nx,3:2:2*ny+1]+v[3:2:2*nx+1,2:2:2*ny]+v[3:2:2*nx+1,3:2:2*ny+1])

  return v_c
end

function prolong(nx,ny,v)
  """
  interpolate 'v' to the fine grid
  """
  v_f=zeros(Float64,2*nx+2,2*ny+2)

  a=0.5625; b=0.1875; c= 0.0625

  v_f[2:2:2*nx  ,2:2:2*ny  ] = a*v[2:nx+1,2:ny+1]+b*(v[1:nx  ,2:ny+1]+v[2:nx+1,1:ny]  )+c*v[1:nx  ,1:ny  ]
  v_f[3:2:2*nx+1,2:2:2*ny  ] = a*v[2:nx+1,2:ny+1]+b*(v[3:nx+2,2:ny+1]+v[2:nx+1,1:ny]  )+c*v[3:nx+2,1:ny  ]
  v_f[2:2:2*nx  ,3:2:2*ny+1] = a*v[2:nx+1,2:ny+1]+b*(v[1:nx  ,2:ny+1]+v[2:nx+1,3:ny+2])+c*v[1:nx  ,3:ny+2]
  v_f[3:2:2*nx+1,3:2:2*ny+1] = a*v[2:nx+1,2:ny+1]+b*(v[3:nx+2,2:ny+1]+v[2:nx+1,3:ny+2])+c*v[3:nx+2,3:ny+2]

  return v_f
end

function V_cycle!(nx::Int,ny::Int,num_levels::Int,u::Array{Float64,2},f::Array{Float64,2},level::Int)
  """
  V cycle
  """
  if (level==num_levels)#bottom solve
    u,res=jacobi_relaxation!(level,nx,ny,u,f,50,true)
    return u,res
  end
  #Step 1: Relax Au=f on this grid
  u,res=jacobi_relaxation!(level,nx,ny,u,f,2,true)

  #Step 2: Restrict residual to coarse grid
  res_c=restrict(nx÷2,ny÷2,res)

  #Step 3:Solve A e_c=res_c on the coarse grid. (Recursively)
  e_c=zero(res_c)
  e_c,res_c=V_cycle!(nx÷2,ny÷2,num_levels,e_c,res_c,level+1)

  #Step 4: Interpolate(prolong) e_c to fine grid and add to u
  u=u+prolong(nx÷2,ny÷2,e_c)

  #Step 5: Relax Au=f on this grid
  u,res=jacobi_relaxation!(level,nx,ny,u,f,1,false)
  return u,res
end

end

max_stars_repo_name: AbhilashReddyM/Multigrid-in-Julia

max_issues_repo_name: AbhilashReddyM/Multigrid-in-Julia

max_forks_repo_name: AbhilashReddyM/Multigrid-in-Julia


text:
import Base: show

# TODO: real expression printing
Base.show(io::IO, x::Variable) = print(io, "%", x.id)

function show(io::IO, b::Branch)
  if b == unreachable
    print(io, "unreachable")
  elseif isreturn(b)
    print(io, "return $(repr(b.args[1]))")
  else
    print(io, "br $(b.block)")
    if !isempty(b.args)
      print(io, " (")
      join(io, repr.(b.args), ", ")
      print(io, ")")
    end
    b.condition != nothing && print(io, " unless $(b.condition)")
  end
end

function show(io::IO, b::Block)
  bb = basicblock(b)
  print(io, b.id, ":")
  if !isempty(bb.args)
    print(io, " (")
    for i = 1:length(bb.args)
      print(io, bb.args[i])
      bb.argtypes[i] != Any && print(io, " :: ", bb.argtypes[i])
      i != length(bb.args) && print(io, ", ")
    end
    print(io, ")")
  end
  println(io)
  for (x, st) in b
    print(io, "  ")
    x == nothing || print(io, string("%", x.id), " = ")
    st.expr == nothing ? print(io, "nothing") :
      print(io, st.expr)
    st.type == Any || print(io, " :: ", st.type)
    println(io)
  end
  for br in bb.branches
    println(io, "  ", br)
  end
end

show(io::IO, ir::IR) = foreach(b -> show(io, b), blocks(ir))


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastLoop
using GeometricExamples

const Δt     = 5.
const ntime  = 10000
const nloop  = 200
const nsave  = 10
const nplot  = 10
const run_id = "poincare_1st_fast_firk_dt5"

include("../guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

pinv = guiding_center_4d_ode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_1st.jl")


text:
function main()
    S = readline()
    N = length(S)

    cnt = 0
    for i = 1:N, j = i:N
        space_left = max(0, i - 2)
        space_right = max(0, N - j - 1)
        cnt += parse(Int, S[i:j]) * 2^(space_left + space_right)
    end

    println(cnt)
end

main()


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.

function timestamp(quasigradient::AbstractQuasiGradient, ::AbstractQuasiGradientExecution)
    return quasigradient.data.iterations
end

function current_decision(quasigradient::AbstractQuasiGradient, ::AbstractQuasiGradientExecution)
    return quasigradient.x
end

function start_workers!(::AbstractQuasiGradient, ::AbstractQuasiGradientExecution)
    return nothing
end

function close_workers!(::AbstractQuasiGradient, ::AbstractQuasiGradientExecution)
    return nothing
end

function subobjectives(quasigradient::AbstractQuasiGradient, execution::AbstractQuasiGradientExecution)
    return execution.subobjectives
end

function set_subobjectives(quasigradient::AbstractQuasiGradient, Qs::AbstractVector, execution::AbstractQuasiGradientExecution)
    execution.subobjectives .= Qs
    return nothing
end

function solve_master!(quasigradient::AbstractQuasiGradient, ::AbstractQuasiGradientExecution)
    try
        MOI.optimize!(quasigradient.master)
    catch err
        status = MOI.get(quasigradient.master, MOI.TerminationStatus())
        # Master problem could not be solved for some reason.
        @unpack Q,θ = quasigradient.data
        gap = abs(θ-Q)/(abs(Q)+1e-10)
        # Always print this warning
        @warn "Master problem could not be solved, solver returned status $status. The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
        rethrow(err)
    end
    return MOI.get(quasigradient.master, MOI.TerminationStatus())
end

function iterate!(quasigradient::AbstractQuasiGradient, execution::AbstractQuasiGradientExecution)
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(quasigradient)
    if Q == Inf
        # Early termination log
        log!(quasigradient; status = MOI.INFEASIBLE)
        return MOI.INFEASIBLE
    end
    if Q == -Inf
        # Early termination log
        log!(quasigradient; status = MOI.DUAL_INFEASIBLE)
        return MOI.DUAL_INFEASIBLE
    end
    sense = MOI.get(quasigradient.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    if coeff*Q <= coeff*quasigradient.data.Q
        quasigradient.data.Q = Q
        quasigradient.ξ .= quasigradient.x
    end
    # Determine stepsize
    γ = step(quasigradient,
             quasigradient.data.iterations,
             Q,
             quasigradient.x,
             quasigradient.gradient)
    # Proximal subgradient update
    prox!(quasigradient,
          quasigradient.x,
          quasigradient.gradient,
          γ)
    # Log progress
    log!(quasigradient)
    # Check optimality
    if terminate(quasigradient,
                 quasigradient.data.iterations,
                 Q,
                 quasigradient.x,
                 quasigradient.gradient)
        # Optimal, final log
        log!(quasigradient; optimal = true)
        return MOI.OPTIMAL
    end
    # Calculate time spent so far and check perform time limit check
    time_spent = quasigradient.progress.tlast - quasigradient.progress.tinit
    if time_spent >= quasigradient.parameters.time_limit
        log!(quasigradient; status = MOI.TIME_LIMIT)
        return MOI.TIME_LIMIT
    end
    # Dont return a status as procedure should continue
    return nothing
end