text:
module Carrot
using ConfParser

include("./bot/post_msg.jl")
include("./forex/fetch_data.jl")

# post_msg.jl
export post_forex_prediction, create_request_text
# fetch_data.jl
export getCandles, create_time_series_metrix

@enum Features batchBot forexPredict web

conf = ConfParse("confs/config.ini")
parse_conf!(conf)

const CarrotConfig = Dict(
    "telegram_token"   => retrieve(conf, "telegram", "token"),
    "telegram_chat_id" => retrieve(conf, "telegram", "chat_id"),
    "oanda_token"      => retrieve(conf, "oanda", "token"),
)

function get_feature_mappings()::Dict
    feature_map = Dict()
    for (index, feature) in enumerate(instances(Features))
        feature_map[string(feature)] = index - 1
    end
    return feature_map
end

function launch(command::String)
    @info "command is $command"
    features = get_feature_mappings()
    if haskey(features, command)
        feature = features[command]
        return launch(Val(Features(feature)))
    end
    return "unknown"
end

function launch(::Val{batchBot})
    return "batch bot"
end

function launch(::Val{forexPredict})
    return "forex predict"
end

function launch(::Val{web})
    return "web"
end

end


text:
# general methods defined

# -------------- general method: create tag string from a tag type instance
function new_tag( d::SingleHtmlTag )
    local res::Vector = String[ "<" * d.tag, ]
    # add id, class, inline CSS style
    isa(d.id, Nothing) ? nothing : push!(res, "id=" * add_doublequotations(d.id) )
    isa(d.class, Nothing) ? nothing : push!(res, "class=" * add_doublequotations(d.class) )
    isa(d.style, Nothing) ? nothing : push!(res, "style=" * add_doublequotations(d.style) )
    # add paired attributes
    append!(res, [ key * "=" * add_doublequotations(value) for (key,value) in d.pairedAttributes ])
    # add single attributes
    append!(res, d.singleAttributes)
    # close the tag
    push!(res, "/>")
    # join the tag as a string, each part is separated by a space
    return join(res, " ")::String
end # new tag
# --------
function new_tag( d::PairedHtmlTag )
    local res::Vector = String[ "<" * d.tag, ]
    # add id, class, inline CSS style
    isa(d.id, Nothing) ? nothing : push!(res, "id=" * add_doublequotations(d.id) )
    isa(d.class, Nothing) ? nothing : push!(res, "class=" * add_doublequotations(d.class) )
    isa(d.style, Nothing) ? nothing : push!(res, "style=" * add_doublequotations(d.style) )
    # add paired attributes
    append!(res, [ key * "=" * add_doublequotations(value) for (key,value) in d.pairedAttributes ])
    # add single attributes
    append!(res, d.singleAttributes)
    # close the tag
    push!(res, ">")
    # add contents between paired tags & the right tag
    if length(d.content) == 0
        # if there is no content, just skip it
        nothing
    elseif (length(d.content) == 1) & !(isa(d.content[1], PairedHtmlTag))
        # if there is only one element and it is not a paired tag (only paired tag has content!)
        push!(res, new_tag(d.content[1]))
    else
        # if there is more than one element in the tag`s content, loop and/or recursive stringlize all elements
        for tmptag in d.content
            if isa(tmptag, BlankHtmlTag)
                push!(res, tmptag.content)
            elseif isa(tmptag, SingleHtmlTag)
                push!(res, new_tag(tmptag))
            else
                push!(res, new_tag(tmptag)) # NOTE: recursive calling
            end
        end
    end
    # add </tag>
    push!(res, "</" * d.tag * ">" )
    # join the tag as a string, each part is separated by a space
    return join(res, " ")::String
end # new tag
# --------
new_tag( d::BlankHtmlTag ) = return (d.content)::String

# ------------------------------ convert a vector of any objects to a list of tags
"""
    vecany2vectag( list::Vector{Any} )

convert a vector of string-lize-able elements to a list of tags
"""
function vecany2vectag( list::Vector{Any} )
    local res = AbstractHtmlTag[]
    for x in list
        if isa(x, AbstractHtmlTag)
            push!(res, x)
        else
            try
                push!(res, BlankHtmlTag(string(x)) )
            catch
                @error("unsupported element type: $(typeof(x))")
            end
        end
    end
    return res::Vector{AbstractHtmlTag}
end  # vecany2vectag

# ------------------------------ add element to a paired tag`s content
"""
    add!( suptag::PairedHtmlTag, subtag::AbstractHtmlTag )

add a new sub-tag member to a given paired tag instance.
"""
function add!( suptag::PairedHtmlTag, subtag::AbstractHtmlTag )
    push!(suptag.content, subtag)
end # add
# ------------------------------ reload add! to add String
"""
    add!( suptag::PairedHtmlTag, subtag::String )

add a new sub-tag member (in the form of string) to a given paired tag instance.
"""
function add!( suptag::PairedHtmlTag, subtag::String )
    push!(suptag.content, BlankHtmlTag(subtag) )
end # add
# ------------------------------ reload add! to allow vector input
"""
    add!( suptag::PairedHtmlTag, v_subtag::Vector{T} where T <: AbstractHtmlTag )

add a list of new sub-tag members to a given paired tag instance.
"""
function add!( suptag::PairedHtmlTag, v_subtag::Vector{T} where T <: AbstractHtmlTag )
    for x in v_subtag; push!(suptag.content, x ); end
end # add
"""
    add!( suptag::PairedHtmlTag, v_subtag::Vector{String} )

add a list of new sub-tag members to a given paired tag instance.
"""
function add!( suptag::PairedHtmlTag, v_subtag::Vector{String} )
    for x in v_subtag; push!(suptag.content, BlankHtmlTag(x) ); end
end # add
"""
    add!( suptag::PairedHtmlTag, v_subtag::Vector{Any} )

add a mixed list of new sub-tag members to a given paired tag instance.
"""
function add!( suptag::PairedHtmlTag, v_subtag::Vector{Any} )
    for x in v_subtag; add!(suptag, x ); end
end # add
# -----------------------------------

# ------------------------------ pop! reload for tag types
"""
    pop!( d::PairedHtmlTag, k::Int )

reload of pop!(). deletes the item with key k in a PairedHtmlTag,
and returns the value associated with k.
"""
pop!( d::PairedHtmlTag, k::Int ) = pop!(d.content, k)::AbstractHtmlTag

# ----------------------- StrOrTag -> Tags
StrOrTag2Tag(content::StrOrTag) = return ( isa(content, String) ? BlankHtmlTag(content) : content )::AbstractHtmlTag
# ----------------------- quick "empty" tags
function quickEmptyPairedTag(tagname::String, content::StrOrTag)
    return PairedHtmlTag(tagname, nothing, nothing, nothing, content = AbstractHtmlTag[StrOrTag2Tag(content),] )
end # quickEmptyPairedTag
function quickEmptyPairedTag(tagname::String, content::Vector)
    return PairedHtmlTag(tagname, nothing, nothing, nothing, content = vecany2vectag(content) )
end # quickEmptyPairedTag
function quickEmptySingleTag(tagname::String)
    return SingleHtmlTag(tagname, nothing, nothing, nothing )
end # quickEmptyPairedTag

# -------------- quickly packaging your contents with a paired tag
"""
    add_pairedtags( tag::String, content::String )

a quick and convenient method to package your content(s) with a paired tag (e.g. <div>)
**WITHOUT ANY ATTRIBUTE**. For example:
add_pairedtags("del", "this is a sample.") --> "<del>this is a sample.</del>";
returns a string.
"""
add_pairedtags( tag::String, content::String ) = "<" * tag * ">" * content * "</" * tag * ">"

# -------------- add \"\" to a string
"""
    add_doublequotations( s::String )

add explicit double quotations to the given string,
e.g. add_doublequotations("fuck") --> "\"fuck\""
"""
add_doublequotations( s::String ) = begin
    if s[1] == '"' == s[end]
        return s::String
    else
        return ( "\"" * s * "\"" )::String
    end
end # add_doublequotations

#

ext: jl

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

# # Computational Modeling of Behavioral Data by Prof. Kentaro Katahira
#
# ## Rescorla-Wagner model

using Plots
using Interact
using Random

# +
"""
Nₜ: number of trials
α: learning rate
Pᵣ: probability of getting reward
"""

@manipulate for Nₜ = 0:1:500, α = 0:0.05:1, Pᵣ = 0:0.05:1

    rng = MersenneTwister(1234) #create a seed for random numbers

    𝐕 = zeros(Nₜ) #strengths of association as Nₜ-length vector
    𝐑 = rand(rng, Nₜ) .< Pᵣ # presence of reinforcement (1 or 0) as Nₜ-length vector

    for t in 1: Nₜ-1

        𝐕[t+1] = 𝐕[t] + α *(𝐑[t]-𝐕[t])
    end

    plot(𝐕, label= string("a ", α))
    plot!([(i, Pᵣ) for i in 1:1:Nₜ], label="expected value of r: " * string(Pᵣ))
    xlabel!("number of trials")
    ylabel!("strength of association")
    ylims!((0, 1))
    title!("Rescorla-Wagner model")
end
# -

# ## Q-learning simulation
# ### softmax function

function softmax(β, Δq)
    return 1 / (1+ exp(-β * (Δq)))
end

@manipulate for β in 0:0.05:5
    plot([(Δq, softmax(β, Δq)) for Δq in -4:0.1:4], m=:o, label=string("beta ", β))
    xlabel!("difference in Q")
    ylabel!("probability")
    ylims!((0, 1))
    title!("Softmax Function")
end

# ### interactive plot of Q-learning model

# +
"""
Nₜ: number of trials
α: learning rate
β: inverse temperature
Pᵣ: probability of getting reward in A
"""

@manipulate for Nₜ in 0:5:200, α in 0:0.05:1, β in 0:0.25:5, Pᵣ in 0:0.05:1

    rng = MersenneTwister(1234)

    𝐐 = zeros(Real, (2, Nₜ)) #initial value of Q in 2 by Nₜ matrix
    𝐜 = zeros(Int, Nₜ) #initial choice in each Nₜ trial
    𝐫 = zeros(Nₜ) # 0 (no reward) or 1 (reward) in each Nₜ trial
    Pₐ = zeros(Nₜ) # probability of choosing A in each trial
    P = (Pᵣ, 1-Pᵣ)

    for t in 1:Nₜ-1
        Pₐ = softmax(β, 𝐐[1, t] - 𝐐[2, t])

        if rand(rng) < Pₐ
            𝐜[t] = 1 #choose A
            𝐫[t] = Int(rand(rng) < P[1])
        else
            𝐜[t] = 2 #choose B
            𝐫[t] = Int(rand(rng) < P[2])
        end

        𝐐[𝐜[t], t+1] = 𝐐[𝐜[t], t] + α * (𝐫[t] - 𝐐[𝐜[t], t])
        𝐐[3 - 𝐜[t], t+1] = 𝐐[3 - 𝐜[t], t] # retain value of unpicked choice
    end

    plot(𝐐[1, :], label="Qt(A)", color="orange")
    plot!([(i, P[1]) for i in 1:1:Nₜ], label="expected value of reward for A:" * string(P[1]), color="darkorange")
    plot!(𝐐[2, :], label="Qt(B)", color="skyblue")
    plot!([(i, P[2]) for i in 1:1:Nₜ], label="expected value of reward for B:" * string(P[2]), color="darkblue")
    xlabel!("number of trials")
    ylabel!("Q (value of behavior?)")
    ylims!((0, 1))
    title!("Q-learning model")
end
# -

# ## Parameter Estimation of Q-learing model
#
# ### Preparation

function generate_qlearning_data(Nₜ, α, β, Pᵣ)

    𝐐 = zeros((2, Nₜ)) #initial value of Q in 2 by Nₜ matrix
    𝐜 = zeros(Int, Nₜ) #initial choice in each Nₜ trial
    𝐫 = zeros(Nₜ) # 0 (no reward) or 1 (reward) in each Nₜ trial
    Pₐ = zeros(Nₜ) # probability of choosing A in each trial
    P = (Pᵣ, 1-Pᵣ)

    for t in 1:Nₜ-1
        Pₐ = softmax(β, 𝐐[1, t] - 𝐐[2, t])

        if rand() < Pₐ
            𝐜[t] = 1 #choose A
            𝐫[t] = (rand() < P[1])
        else
            𝐜[t] = 2 #choose B
            𝐫[t] = Int(rand() < P[2])
        end

        𝐐[𝐜[t], t+1] = 𝐐[𝐜[t], t] + α * (𝐫[t] - 𝐐[𝐜[t], t])
        𝐐[3 - 𝐜[t], t+1] = 𝐐[3 - 𝐜[t], t] # retain value of unpicked choice
    end

    return 𝐜, 𝐫
end

"""
init_values: [α, β]
α: learning rate
β: inverse temperature
𝐜: vector of choices in each Nₜ trial in 1(A) or 2(B)
𝐫: 0 (no reward) or 1 (reward) in each Nₜ trial
"""
function func_qlearning(init_values, 𝐜, 𝐫) #needed for passing list as variables for Optim

    Nₜ = length(𝐜)
    Pₐ = zeros(Nₜ) #probabilities of selecting A
    𝐐 = zeros(Real, (2, Nₜ))
    logl = 0 #initial value of log likelihood

    for t in 1:Nₜ - 1
        Pₐ[t] = softmax(init_values[2], 𝐐[1, t] - 𝐐[2, t])
        logl += (𝐜[t] == 1) * log(Pₐ[t]) + (𝐜[t] == 2) * log(1 - Pₐ[t])
        𝐐[𝐜[t], t + 1] = 𝐐[𝐜[t], t] + init_values[1] * (𝐫[t] - 𝐐[𝐜[t], t])
        𝐐[3 - 𝐜[t], t + 1] =  𝐐[3 - 𝐜[t], t]
    end

    return (negll = -logl, 𝐐 = 𝐐, Pₐ = Pₐ);
end

# ## Parameter Estimation
#
# ### optimization with JuMP and Ipopt

import Pkg
Pkg.add("Pkg")
Pkg.add("Ipopt")
Pkg.build("Ipopt")

# +
using JuMP, Ipopt, ForwardDiff

#@manipulate for Nₜ in 0:50:1000, α1 in 0:0.05:1, β1 in 0:0.25:5, Pᵣ in 0:0.05:1

Nₜ=500
α1 = 0.3
β1 = 0.2
Pᵣ = 0.5

𝐜, 𝐫 = generate_qlearning_data(Nₜ, α1, β1, Pᵣ)
func_qlearning_JuMP(α, β) = func_qlearning((α, β), 𝐜, 𝐫).negll #JuMP requires separate arguments, not a list

m = Model(Ipopt.Optimizer)
register(m, :func_qlearning_JuMP, 2, func_qlearning_JuMP, autodiff=true)

@variable(m, 0.0 <= α <= 1.0, start=rand(), base_name = "learning_rate")
@variable(m, 0.0 <= β <= 5.0, start=5*rand(), base_name = "inverse_temperature")

@NLobjective(m, Min, func_qlearning_JuMP(α, β))
optimize!(m)
print(""," α = ", value(α), " β = ", value(β))
#end
# -

# #### optimization with Optim

# +
using Optim

@manipulate for Nₜ in 0:5:200, α in 0:0.05:1, β in 0:0.25:5, Pᵣ in 0:0.05:1
    𝐜, 𝐫 = generate_qlearning_data(Nₜ, α, β, Pᵣ)

    func_qlearning_opt(init_values) = func_qlearning(init_values, 𝐜, 𝐫).negll

    initial_values = rand(2)
    lower = [0.0, 0.0]
    upper = [1.0, 5.0]
    inner_optimizer = GradientDescent()
    results = optimize(func_qlearning_opt, lower, upper, initial_values, Fminbox(inner_optimizer))
    #@show optimize(func_qlearning_opt, init_values, lower, upper, LBFGS())
end
# -

# #### optimization with BlackBoxOptim, which is designed for blackbox functions, so this part is only for demonstration purpose

# +
using BlackBoxOptim

@manipulate for Nₜ in 0:5:200, α in 0:0.05:1, β in 0:0.25:5, Pᵣ in 0:0.05:1
    𝐜, 𝐫 = generate_qlearning_data(Nₜ, α, β, Pᵣ)

    func_qlearning_opt(init_values) = func_qlearning(init_values, 𝐜, 𝐫).negll

    results = bboptimize(func_qlearning_opt; SearchRange = [(0.0, 1.0), (0.0, 5.0)], NumDimensions = 2);
    best_candidate(results)
end
# -

# #### We can also compare performances when using different optimizers.

# +
#this cell takes a lot time to run, so execute it only if you want to

#𝐜, 𝐫 = generate_qlearning_data(100, 0.3, 1.2, 0.5)
#func_qlearning_opt(init_values) = func_qlearning(init_values, 𝐜, 𝐫).negll
#compare_optimizers(func_qlearning_opt; SearchRange = [(0.0, 1.0), (0.0, 5.0)], NumDimensions = 2);
# -

# ## comparison of models
#
# ### win-stay lose-shift (WSLS) model

"""
Nₜ: number of trials
ϵ: error rate
Pᵣ: probability of getting reward in A
"""
function wsls_simulstion(Nₜ, ϵ, Pᵣ, seed=1234)

    rng = MersenneTwister(seed)

    Pₐ = zeros(Nₜ) #probabilities of selecting A
    Pₐ[1] = 0.5 # probability at initial trial is 0.5
    𝐜 = zeros(Int, Nₜ) #initial choice in each Nₜ trial
    𝐫 = zeros(Nₜ) # 0 (no reward) or 1 (reward) in each Nₜ trial

    for t in 1:Nₜ-1

        chooseAB = rand(rng)
        get_reward = rand(rng)

        #select A with reward
        if chooseAB < Pₐ[t] && get_reward <  Pᵣ

            Pₐ[t + 1] = 1 - ϵ
            𝐜[t] = 1
            𝐫[t] = 1

        #select B with no reward
        elseif chooseAB > Pₐ[t] && get_reward >  Pᵣ

            Pₐ[t + 1] = 1 - ϵ
            𝐜[t] = 2
            𝐫[t] = 0

        #select A with no reward
        elseif chooseAB < Pₐ[t] && get_reward >  Pᵣ

            Pₐ[t + 1] = ϵ
            𝐜[t] = 1
            𝐫[t] = 0
        #select B with reward
        elseif chooseAB > Pₐ[t] && get_reward <  Pᵣ

            Pₐ[t + 1] = ϵ
            𝐜[t] = 2
            𝐫[t] = 1

        end

    end

    return (Pₐ = Pₐ, 𝐜 = 𝐜, 𝐫 = 𝐫);
end

# #### plot

# +
@manipulate for Nₜ in 0:5:200, ϵ in 0:0.05:1, Pᵣ in 0:0.05:1, seed in 1:1:1234

    Pₐ = wsls_simulstion(Nₜ, ϵ, Pᵣ, seed).Pₐ

    plot(Pₐ, label="P(a = A)", color="orange")
    ylabel!("P(a = A)")
    ylims!((0, 1))
    title!("WSLS Model")

end
# -

# ### random selection model

function random_choice_simulation(Nₜ, Pₐ, seed=1234)

    rng = MersenneTwister(seed)

    𝐜 = 2 .- Int.(rand(rng, Nₜ) .< Pₐ) #dot notation in Julia signifies elemnet-wise operation

    return (Pₐ = Pₐ, 𝐜 = 𝐜)
end

# ####plot

# +
@manipulate for Nₜ in 0:5:200, Pₐ in 0:0.05:1

    plot([Pₐ for i in range(1, stop=Nₜ)], label="P(a = A)", color="orange")
    ylabel!("P(a = A)")
    ylims!((0, 1))
    title!("Random Choice Model")

end
# -

# ### model comparison
#
# #### preparation

# +
"""
ϵ: error rate
𝐜: vector of choices in each Nₜ trial in 1(A) or 2(B)
𝐫: 0 (no reward) or 1 (reward) in each Nₜ trial

when given ϵ, 𝐜, and 𝐫, returns log likelihood and Pₐ
"""
function func_wsls(ϵ, 𝐜, 𝐫)

    Nₜ = length(𝐜)
    Pₐ = zeros(Nₜ) #probabilities of selecting A
    Pₐ[1] = 0.5
    logl = 0 #initial value of log likelihood

    for t in 1:Nₜ - 1
        logl += (𝐜[t] == 1) * log(Pₐ[t]) + (𝐜[t] == 2) * log(1 - Pₐ[t])

        #select A with reward
        if 𝐜[t] == 1 &&   𝐫[t] == 1

            Pₐ[t + 1] = 1 - ϵ

        #select B with no reward
        elseif  𝐜[t] == 2 &&   𝐫[t] == 0

            Pₐ[t + 1] = 1 - ϵ

        #select A with no reward
        elseif  𝐜[t] == 1 &&   𝐫[t] == 0

            Pₐ[t + 1] = ϵ

        #select B with reward
        elseif 𝐜[t] == 2 &&   𝐫[t] == 1

            Pₐ[t + 1] = ϵ

        end
    end

    return (ll = logl, Pₐ = Pₐ);
end

"""
Pₐ: probability of choosing A
𝐜: vector of choices in each Nₜ trial in 1(A) or 2(B)
𝐫: 0 (no reward) or 1 (reward) in each Nₜ trial

when given Pₐ, 𝐜, and 𝐫, returns log likelihood and Pₐ
"""
function func_random_choice(Pₐ, 𝐜, 𝐫)

    Nₜ = length(𝐜)
    logl = 0

    for t in 1:Nₜ
        logl += (𝐜[t] == 1) * log(Pₐ) + (𝐜[t] == 2) * log(1 - Pₐ)
    end

    return logl

end
# -

# #### parameter estimation with JuMP

# +
using JuMP, Ipopt, ForwardDiff

@manipulate for Nₜ in 0:50:1000, α1 in 0:0.05:1, β1 in 0:0.25:5, Pᵣ in 0:0.05:1

    𝐜, 𝐫 = generate_qlearning_data(Nₜ, α1, β1, Pᵣ)
    func_qlearning_JuMP(α, β) = func_qlearning((α, β), 𝐜, 𝐫).negll #JuMP requires separate arguments, not a list

    m = Model(Ipopt.Optimizer)
    register(m, :func_qlearning_JuMP, 2, func_qlearning_JuMP, autodiff=true)

    @variable(m, 0.0 <= α <= 1.0, start=rand(), base_name = "learning_rate")
    @variable(m, 0.0 <= β <= 5.0, start=5*rand(), base_name = "inverse_temperature")

    @NLobjective(m, Min, func_qlearning_JuMP(α, β))
    optimize!(m)
    print(""," α = ", value(α), " β = ", value(β))


text:
module CGE

using DelimitedFiles
using Statistics
using StatsBase
using LinearAlgebra
using louvain_jll
using Random

#auxilary
export parseargs

#landmarks
export landmarks

#divergence
export wGCL
export wGCL_directed

#clustering
export louvain_clust

# Include package code
include("auxilary.jl")
include("landmarks.jl")
include("divergence.jl")
include("clustering.jl")
end


text:
include("constants.jl")
include("properties.jl")
include("dedekind.jl")  #dedekind cuts and the like.
include("h-layer.jl")
include("rounding.jl")
include("convert.jl")
include("comparison.jl")

#empty constructor
Valid{N,ES}() where {N,ES} = Valid{N,ES}(reinterpret(Vnum{N,ES}, zero(@UInt)), reinterpret(Vnum{N,ES}, zero(@UInt)))

include("macros.jl")
include("addsub.jl")
include("mul.jl")
include("div.jl")

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
using Literate

for t in readdir(".")
    if occursin(r"tut_.*.jl", t)
        println("\nTutorial $t in $(pwd())\n")
        Literate.markdown(t, "."; documenter=false);
        # Literate.notebook(t, "."; execute=false, documenter=false);
    end
end


text:
using MyFunRep
using Documenter

ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"

DocMeta.setdocmeta!(MyFunRep, :DocTestSetup, :(using MyFunRep); recursive=true)

makedocs(;
    modules=[MyFunRep],
    authors="Rethna Pulikkoonattu et al.",
    repo="https://github.com/geekymode/MyFunRep.jl/blob/{commit}{path}#{line}",
    sitename="MyFunRep.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://geekymode.github.io/MyFunRep.jl",
        assets=String[],
        collapselevel=1,
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => "tutorials/eigenvalues.md",
        "Explanations" => "explanations/preconditioning.md",
    ],
)

deploydocs(;
    repo="github.com/geekymode/MyFunRep.jl",
    devbranch="main",
)


text:
export NearestNeighbor

"""
    NearestNeigbor()

Nearest neighbor interpolation.

"""
struct NearestNeighbor <: InterpolationMethod end

struct NearestNeighborInterpolant{T, TT} <: ScatteredInterpolant
    data::Array{T}
    tree::TT
end

function interpolate(nn::NearestNeighbor,
                     points::AbstractArray{<:Real,2},
                     samples::AbstractArray{<:Number,N};
                     metric = Euclidean()) where {N}

    # Build a kd-tree of the points
    # If we get an adjoint, make a copy (or else KDTree will be sad)
    if points isa LinearAlgebra.Adjoint
        points = copy(points)
    end
    if samples isa LinearAlgebra.Adjoint
        samples = copy(samples)
    end
    tree = KDTree(points, metric)

    return NearestNeighborInterpolant(samples, tree)
end

function evaluate(itp::NearestNeighborInterpolant, points::AbstractArray{<:Real,2})

    # Get the indices for each points closest neighbor
    inds, _ = knn(itp.tree, points, 1)

    m = size(points, 2)
    n = size(itp.data, 2)
    values = zeros(eltype(itp.data), m, n)

    # knn returns a vector of vectors, so we need a loop
    for i in 1:m
        values[i, :] = itp.data[inds[i][1], :]
    end

    return values
end


text:
include("../src/HssMatrices.jl")
using .HssMatrices
using LinearAlgebra
using BenchmarkTools
using Random
Random.seed!(123)

### run benchmarks on Cauchy matrix
K(x,y) = (x-y) != 0 ? 1/(x-y) : 10000.
A = [ K(x,y) for x=-1:0.001:1, y=-1:0.001:1];
b = randn(size(A,2), 5);

# test the simple implementation of cluster trees
m, n = size(A)
#lsz = 64;
lsz = 64;
rcl = bisection_cluster(1:m, leafsize=lsz)
ccl = bisection_cluster(1:n, leafsize=lsz)

hssA = compress(A, rcl, ccl);

# time conversion to dense
println("Benchmarking full...")
@btime full(hssA);

# time access
println("Benchmarking getindex...")
ii = randperm(100); jj = randperm(100)
@btime Aij = hssA[ii,jj];

# time compression
println("Benchmarking compression...")
@btime hssA = compress(A, rcl, ccl);

# time compression
println("Benchmarking randomized compression...")
@btime hssA = randcompress_adaptive(A, rcl, ccl);

println("Benchmarking re-compression...")
hssB = copy(hssA)
@btime hssA = recompress!(hssB; atol=1e-3, rtol=1e-3);

println("Benchmarking proper...")
hssB = hssA + hssA
@btime orthonormalize_generators!(hssB)

println("Benchmarking addition...")
@btime hssC = hssA + hssA

# time matvec
println("Benchmarking matvec...")
x = randn(size(A,2), 10);
@btime y = hssA*x

# time ulvfactsolve
println("Benchmarking ulvfactsolve...")
b = randn(size(A,2), 10);
@btime x = ulvfactsolve(hssA, b);

# time hssldivide
println("Benchmarking hssldivide...")
hssX = compress(1.0*Matrix(I, n, n), ccl, ccl, atol=1e-3, rtol=1e-3); # this should probably be one constructor
@btime hssC = ldiv!(copy(hssA), hssX);


text:
# f = (x,n) -> (x^2 - 1) % n
# println(f(3, 5))

# function gcd(a,b)
# 	while b != 0
# 		a, b = b, (a % b)
# 	end
# 	a
# end

# function pollardrho(n)
# 	f = (x,n) -> (x^2 - 1) % n
# 	x, y, d = 2, 2, 1
# 	while d == 1
# 		x = f(x,n)
# 		y = f(f(y,n),n)
# 		d = gcd(abs(x-y), n)
# 	end
# 	d
# end

function isprime(n)
	if n & 1 == 0
		return false
	end
	for i = 3:sqrt(n)
		if n % i == 0
			return false
		end
	end
	return true
end

N = 600851475143
largest_factor = 0

for candidate = 3:convert(Int64, ceil(sqrt(N)))
	if N % candidate == 0
		if isprime(candidate)
			largest_factor = candidate
		elseif isprime(convert(Int64, N / candidate))
			largest_factor = convert(Int64, N / candidate)
		end
	end
end
println("The answer to Euler Problem 3 is $largest_factor")

max_stars_repo_path: Problem 3.jl

max_stars_repo_name: andrew-christianson/Polyglot-Euler

max_issues_repo_path: Problem 3.jl

max_issues_repo_name: andrew-christianson/Polyglot-Euler

max_forks_repo_path: Problem 3.jl

max_forks_repo_name: andrew-christianson/Polyglot-Euler


text:
n_microbes      = 1
n_reserves      = 2
E               = fill(0.5,n_microbes*n_reserves)
E               = reshape(E, n_reserves, n_microbes)
V               = fill(1.0, n_microbes)

Isolate         = BaseGenome(1e6*ones(n_microbes), 10*ones(n_microbes), 1.0*ones(n_microbes), 1.0*ones(1,n_microbes), 1.0*ones(n_microbes))

p               = PMetabolismCN{BaseGenome, Array{Float64,2}, Array{Float64,1}}(Isolate)
p.α             = fill(0.1, n_microbes)
p.α_n           = fill(0.1, n_microbes)
p.y_EM          = ones(n_reserves, n_microbes)
p.y_EX          = ones(n_reserves, n_microbes)
p.y_EX[1,1]     = 0.1
p.y_EX[2,1]     = 0.3
r               = growth!(zeros(size(E,2)), E, V, p)

x, rM_CO2, rEG_rej, rEX_rej = growth_production(r, E, V, p)


text:
"""
Defines the `Documenter.jl` build "pipeline" named [`DocumentPipeline`](@ref).

Each stage of the pipeline performs an action on a [`Documents.Document`](@ref) object.
These actions may involve creating directory structures, expanding templates, running
doctests, etc.
"""
module Builder

import ..Documenter:
    Anchors,
    Selectors,
    Documents,
    Documenter,
    Utilities

using Compat, DocStringExtensions

# Document Pipeline.
# ------------------

"""
The default document processing "pipeline", which consists of the following actions:

- [`SetupBuildDirectory`](@ref)
- [`ExpandTemplates`](@ref)
- [`CrossReferences`](@ref)
- [`CheckDocument`](@ref)
- [`Populate`](@ref)
- [`RenderDocument`](@ref)

"""
abstract type DocumentPipeline <: Selectors.AbstractSelector end

"""
Creates the correct directory layout within the `build` folder and parses markdown files.
"""
abstract type SetupBuildDirectory <: DocumentPipeline end

"""
Executes a sequence of actions on each node of the parsed markdown files in turn.
"""
abstract type ExpandTemplates <: DocumentPipeline end

"""
Finds and sets URLs for each `@ref` link in the document to the correct destinations.
"""
abstract type CrossReferences <: DocumentPipeline end

"""
Checks that all documented objects are included in the document and runs doctests on all
valid Julia code blocks.
"""
abstract type CheckDocument <: DocumentPipeline end

"""
Populates the `ContentsNode`s and `IndexNode`s with links.
"""
abstract type Populate <: DocumentPipeline end

"""
Writes the document tree to the `build` directory.
"""
abstract type RenderDocument <: DocumentPipeline end

Selectors.order(::Type{SetupBuildDirectory})   = 1.0
Selectors.order(::Type{ExpandTemplates})       = 2.0
Selectors.order(::Type{CrossReferences})       = 3.0
Selectors.order(::Type{CheckDocument})         = 4.0
Selectors.order(::Type{Populate})              = 5.0
Selectors.order(::Type{RenderDocument})        = 6.0

Selectors.matcher(::Type{T}, doc::Documents.Document) where {T <: DocumentPipeline} = true

Selectors.strict(::Type{T}) where {T <: DocumentPipeline} = false

function Selectors.runner(::Type{SetupBuildDirectory}, doc::Documents.Document)
    Utilities.log(doc, "setting up build directory.")

    # Frequently used fields.
    build  = doc.user.build
    source = doc.user.source

    # The .user.source directory must exist.
    isdir(source) || error("source directory '$(abspath(source))' is missing.")

    # We create the .user.build directory.
    # If .user.clean is set, we first clean the existing directory.
    doc.user.clean && isdir(build) && rm(build; recursive = true)
    isdir(build) || mkpath(build)

    # We'll walk over all the files in the .user.source directory.
    # The directory structure is copied over to .user.build. All files, with
    # the exception of markdown files (identified by the extension) are copied
    # over as well, since they're assumed to be images, data files etc.
    # Markdown files, however, get added to the document and also stored into
    # `mdpages`, to be used later.
    mdpages = String[]
    for (root, dirs, files) in walkdir(source)
        for dir in dirs
            d = normpath(joinpath(build, relpath(root, source), dir))
            isdir(d) || mkdir(d)
        end
        for file in files
            src = normpath(joinpath(root, file))
            dst = normpath(joinpath(build, relpath(root, source), file))
            if endswith(file, ".md")
                push!(mdpages, Utilities.srcpath(source, root, file))
                Documents.addpage!(doc, src, dst)
            else
                Compat.cp(src, dst; force = true)
            end
        end
    end

    # If the user hasn't specified the page list, then we'll just default to a
    # flat list of all the markdown files we found, sorted by the filesystem
    # path (it will group them by subdirectory, among others).
    userpages = isempty(doc.user.pages) ? sort(mdpages) : doc.user.pages

    # Populating the .navtree and .navlist.
    # We need the for loop because we can't assign to the fields of the immutable
    # doc.internal.
    for navnode in walk_navpages(userpages, nothing, doc)
        push!(doc.internal.navtree, navnode)
    end

    # Finally we populate the .next and .prev fields of the navnodes that point
    # to actual pages.
    local prev::Union{Documents.NavNode, Nothing} = nothing
    for navnode in doc.internal.navlist
        navnode.prev = prev
        if prev !== nothing
            prev.next = navnode
        end
        prev = navnode
    end
end

"""
$(SIGNATURES)

Recursively walks through the [`Documents.Document`](@ref)'s `.user.pages` field,
generating [`Documents.NavNode`](@ref)s and related data structures in the
process.

This implementation is the de facto specification for the `.user.pages` field.
"""
function walk_navpages(visible, title, src, children, parent, doc)
    # parent can also be nothing (for top-level elements)
    parent_visible = (parent === nothing) || parent.visible
    if src !== nothing
        src = normpath(src)
        src in keys(doc.internal.pages) || error("'$src' is not an existing page!")
    end
    nn = Documents.NavNode(src, title, parent)
    (src === nothing) || push!(doc.internal.navlist, nn)
    nn.visible = parent_visible && visible
    nn.children = walk_navpages(children, nn, doc)
    nn
end

function walk_navpages(hps::Tuple, parent, doc)
    @assert length(hps) == 4
    walk_navpages(hps..., parent, doc)
end

walk_navpages(title::String, children::Vector, parent, doc) = walk_navpages(true, title, nothing, children, parent, doc)
walk_navpages(title::String, page, parent, doc) = walk_navpages(true, title, page, [], parent, doc)

walk_navpages(p::Pair, parent, doc) = walk_navpages(p.first, p.second, parent, doc)
walk_navpages(ps::Vector, parent, doc) = [walk_navpages(p, parent, doc)::Documents.NavNode for p in ps]
walk_navpages(src::String, parent, doc) = walk_navpages(true, nothing, src, [], parent, doc)

function Selectors.runner(::Type{ExpandTemplates}, doc::Documents.Document)
    Utilities.log(doc, "expanding markdown templates.")
    Documenter.Expanders.expand(doc)
end

function Selectors.runner(::Type{CrossReferences}, doc::Documents.Document)
    Utilities.log(doc, "building cross-references.")
    Documenter.CrossReferences.crossref(doc)
end

function Selectors.runner(::Type{CheckDocument}, doc::Documents.Document)
    Utilities.log(doc, "running document checks.")
    Documenter.DocChecks.missingdocs(doc)
    Documenter.DocChecks.doctest(doc)
    Documenter.DocChecks.footnotes(doc)
    Documenter.DocChecks.linkcheck(doc)
end

function Selectors.runner(::Type{Populate}, doc::Documents.Document)
    Utilities.log("populating indices.")
    Documents.doctest_replace!(doc)
    Documents.populate!(doc)
end

function Selectors.runner(::Type{RenderDocument}, doc::Documents.Document)
    count = length(doc.internal.errors)
    if doc.user.strict && count > 0
        error("`makedocs` encountered $(count > 1 ? "errors" : "an error"). Terminating build")
    else
        Utilities.log(doc, "rendering document.")
        Documenter.Writers.render(doc)
    end
end

Selectors.runner(::Type{DocumentPipeline}, doc::Documents.Document) = nothing

end


text:
__precompile__()

module Configs
using JSON, YAML, Dates, Memoize
export  getconfig,
        setconfig!,
        hasconfig,
        initconfig

mutable struct Confstate
    init::Bool
    isimmutable::Bool
    immutable::NamedTuple
    configs::Dict{String, Any}
end
const state = Confstate(false, false, NamedTuple(), Dict())
const defaultpath = normpath(joinpath(dirname(Base.active_project()), "configs"))

"""
    initconfig(; <keyword arguments>)
Optionally override the configs environment.

Automatically called on first [`getconfig`](@ref), [`hasconfig`](@ref) or [`setconfig!`](@ref)

Can only be called once.

Alternatively, set the corresponding `ENV` variables `DEPLOYMENT_KEY` or `CONFIGS_DIRECTORY` to preference

## Arguments
- `deployment_key::String`: The `ENV` key that defines the intended deployment (eg. production, staging, etc.)
Default: `DEPLOYMENT`
- `configs_directory::String`: The directory containing config definitions
May be relative to project root or absolute
Default: `configs`

"""
function initconfig(; deployment_key = "DEPLOYMENT", configs_directory = defaultpath)
    state.init && throw(Configserror("configs has already been initialised"))
    state.init = true
    conf = parseconfigs(deployment_key, configs_directory)
    len = length(conf.files)
    newtrees = asyncmap(readconffile, fill(configs_directory, len), conf.files)
    for newtree in newtrees
        override!(state.configs, newtree)
    end
end

function setconfig!(path::String, value::Union{Number, Bool})
    value = configpathtodict(path, value)
    makeconfig!(state, value)
end
function setconfig!(path::String, value::String)
    try
        value = JSON.parse(value)
    catch err
    end
    value = configpathtodict(path, value)
    makeconfig!(state, value)
end
"""
    setconfig!(x, y)
Sets a configuration value `y` for given configuration path `x`

Can not be called after initial calls to [`getconfig`](@ref) or [`hasconfig`](@ref)
"""
function setconfig!(path::String, value::Union{Tuple, Array, NamedTuple, Dict})
    value = json(value) |> JSON.parse
    value = configpathtodict(path, value)
    makeconfig!(state, value)
end

"""
    getconfig([x])
Returns a configuration value for given configuration path x

The returned value may be absolute or a `NamedTuple`, which can be further accessed with the syntax:
`value.x.y`

"""
@memoize Dict{Tuple{String}, Any} function getconfig(path::String = "")
    !state.init && initconfig()
    !state.isimmutable && (
        state.immutable= makeimmutable(state.configs);
        state.isimmutable = true
    )
    path === "" && return state.immutable
    parseconfigpath(state.immutable, path)
end

"""
    hasconfig(x)
Returns a `Bool` indicating if the given path `x` exists.

"""
@memoize Dict{Tuple{String}, Any} function hasconfig(path::String)::Bool
    path === "" && throw(Configserror("a path is required to query a config"))
    try
        getconfig(path)
        true
    catch err
        false
    end
end

include("utils.jl")
end

max_stars_repo_name: citkane/Config

max_issues_repo_name: citkane/Config

max_forks_repo_name: citkane/Config


text:
module Localizers

using Distributions
using LinearAlgebra
using DocStringExtensions

export NoLocalization, Delta, RBF, BernoulliDropout
export LocalizationMethod, Localizer

abstract type LocalizationMethod end

struct NoLocalization <: LocalizationMethod end
struct Delta <: LocalizationMethod end

"""
    RBF{FT <: Real} <: LocalizationMethod

Radial basis function localization method. Covariance
terms are damped as d(i,j)= |i-j|/l increases, following
a Gaussian.

# Fields

$(TYPEDFIELDS)
"""
struct RBF{FT <: Real} <: LocalizationMethod
    "Length scale defining the RBF kernel"
    lengthscale::FT
end

"""
    BernoulliDropout{FT <: Real} <: LocalizationMethod

Localization method that drops cross-covariance terms with
probability 1-p, retaining a Hermitian structure.

# Fields

$(TYPEDFIELDS)
"""
struct BernoulliDropout{FT <: Real} <: LocalizationMethod
    "Probability of keeping a given cross-covariance term"
    prob::FT
end

struct Localizer{LM <: LocalizationMethod, T}
    kernel::Union{AbstractMatrix{T}, UniformScaling{T}}
end

"Uniform kernel constructor"
function Localizer(localization::NoLocalization, p::IT, d::IT, T = Float64) where {IT <: Int}
    kernel = ones(T, p, d)
    return Localizer{NoLocalization, T}(kernel)
end

"Delta kernel localizer constructor"
function Localizer(localization::Delta, p::IT, d::IT, T = Float64) where {IT <: Int}
    kernel = T(1) * Matrix(I, p, d)
    return Localizer{Delta, T}(kernel)
end

"RBF kernel localizer constructor"
function Localizer(localization::RBF, p::IT, d::IT, T = Float64) where {IT <: Int}
    l = localization.lengthscale
    kernel = zeros(T, p, d)
    for i in 1:p
        for j in 1:d
            @inbounds kernel[i, j] = exp(-(i - j) * (i - j) / (2 * l * l))
        end
    end
    return Localizer{RBF, T}(kernel)
end

"Randomized Bernoulli dropout kernel localizer constructor"
function Localizer(localization::BernoulliDropout, p::IT, d::IT, T = Float64) where {IT <: Int}
    kernel = T(1) * Matrix(I, p, d)

    # Transpose
    kernel = p < d ? kernel' : kernel
    for i in 2:size(kernel, 1)
        for j in 1:min(i - 1, size(kernel, 2))
            @inbounds kernel[i, j] = T(rand(Bernoulli(localization.prob)))
            if i <= size(kernel, 2)
                kernel[j, i] = kernel[i, j]
            end
        end
    end
    # Transpose back
    kernel = p < d ? kernel' : kernel

    return Localizer{BernoulliDropout, T}(kernel)
end

end # module


text:
using Distributions, Compat.Test, Random, LinearAlgebra
using Distributions: Product

@testset "Testing Product distributions" begin
let
    rng, D = MersenneTwister(123456), 11

    # Construct independent distributions and `Product` distribution from these.
    μ = randn(rng, D)
    ds = Normal.(μ, 1.0)
    x = rand.(Ref(rng), ds)
    d_product = Product(ds)

    # Check that methods for `Product` are consistent.
    @test length(d_product) == length(ds)
    @test logpdf(d_product, x) ≈ sum(logpdf.(ds, x))
    @test mean(d_product) == mean.(ds)
    @test var(d_product) == var.(ds)
    @test cov(d_product) == Diagonal(var.(ds))
    @test entropy(d_product) == sum(entropy.(ds))
end
end


text:
export group_by_type
export rectangular_lattice
export triangular_lattice
export wedge
export line
export remove_overlaps!
export remove_shape!
export nearest_neighbor_pairs

"""
    group_by_type(particles; ptype)

Returns a group of particles with type(s) `ptype`.
"""
function group_by_type(particles::Array{Particle, 1}; ptype::Union{Symbol, Array{Symbol, 1}})
    if !(ptype isa Array)
        ptype = [ptype]
    end

    pgroup = Array{Particle, 1}()
    for particle in particles
        if particle.ptype in ptype
            push!(pgroup, particle)
        end
    end
    return pgroup
end

"""
    rectangular_lattice(; s_x, s_y, M_x, M_y)

Generates a rectangular lattice with lattice constants `s_x, s_y` in the x, y
directions.  Each cell is duplicated `M_x, M_y` times in the x, y directions.
"""
function rectangular_lattice(; s_x::Float64, s_y::Float64, M_x::Int64, M_y::Int64)
    positions = Array{Array{Float64, 1}, 1}()
    for i = 0 : M_x - 1, j = 0 : M_y - 1
        push!(positions, [i * s_x, j * s_y])
    end
    return positions, M_x * s_x, M_y * s_y
end

"""
    triangular_lattice(; s, M_x, M_y)

Generates a triangular lattice with lattice constant `s`.  Each cell consists of
points `(0, 0)` and `(s / 2, s * sqrt(3) / 2)`.  The cells are duplicated
`M_x, M_y` times in the x, y directions.  Note that there are other
differently oriented triangular lattices.
"""
function triangular_lattice(; s::Float64, M_x::Int64, M_y::Int64)
    cell = [[0.0, 0.0], [s / 2, s * sqrt(3) / 2]]
    positions = Array{Array{Float64, 1}, 1}()
    for i = 0 : M_x - 1, j = 0 : M_y - 1
        push!(positions, cell[1] + [i * s, j * s * sqrt(3)])
        push!(positions, cell[2] + [i * s, j * s * sqrt(3)])
    end
    return positions, M_x * s, M_y * s * sqrt(3)
end

"""
    line(; from, to, spacing)

Generates a line of points from point `from` to point `to` with approximate
spacing `spacing`.  Note that the spacing is automatically adjusted slightly so
that points are equally spaced and include the start and end points.
"""
function line(; from::Array{Float64, 1}, to::Array{Float64, 1}, spacing::Float64)
    edge_length = sqrt((to[1] - from[1])^2 + (to[2] - from[2])^2)
    M = ceil(Int64, edge_length / spacing)

    edge = [from]
    for i = 1 : M
        push!(edge, from + i / M * (to - from))
    end
    return edge
end

"""
    wedge(; angle, edge_length, spacing, orientation, tip)

Generates a V-shaped wedge.  The angle of the wedge is `angle` with each arm a
length `edge_length`.  The points are separated by approximately `spacing` so
that the start and end points are included.  `orientation` is the angle relative
to the x-axis and tip is the x, y location of the tip of the wedge.
"""
function wedge(; angle::Float64, edge_length::Float64, spacing::Float64, orientation::Float64 = 0.0, tip::Array{Float64, 1})
    M = ceil(Int64, edge_length / spacing)
    approx_spacing = edge_length / M

    funnel = [tip]
    for i = 1 : M
        push!(funnel, tip + i * approx_spacing * [cos(angle / 2 + orientation), sin(angle / 2 + orientation)])
        push!(funnel, tip + i * approx_spacing * [cos(-angle / 2 + orientation), sin(-angle / 2 + orientation)])
    end
    return funnel
end

"""
    remove_overlaps!(positions; fixed_positions, minimum_distance, period_x, period_y)

Removes points from `positions` that are within a `minimum_distance` of points
in `fixed_positions`.  This is used to avoid large forces due to overlaps.
"""
function remove_overlaps!(positions::Array{Array{Float64, 1}, 1}; fixed_positions::Array{Array{Float64, 1}, 1}, minimum_distance::Float64, period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    indices_to_remove = Array{Int64, 1}()
    for (n, (x, y)) in enumerate(positions)
        for (xx, yy) in fixed_positions
            Δx = wrap_displacement(x - xx; period = period_x)
            Δy = wrap_displacement(y - yy; period = period_y)

            if Δx^2 + Δy^2 < minimum_distance^2
                push!(indices_to_remove, n)
                break
            end
        end
    end
    deleteat!(positions, indices_to_remove)
end

"""
    remove_shape!(positions; shape, in_or_out)

Remove positions that are either inside (:in) or outside (:out) of a
parametric shape.
"""
function remove_shape!(positions::Array{Array{Float64, 1}, 1}; shape::Array{Array{Float64, 1}, 1}, in_or_out::Symbol)
    indices_to_remove = Array{Int64, 1}()
    N_shape = length(shape)
    for (n, (x, y)) in enumerate(positions)
        θ = 0.0
        for p = 1 : N_shape
            next_p = p % N_shape + 1
            Δx_1, Δy_1 = shape[p][1] - x, shape[p][2] - y
            Δx_2, Δy_2 = shape[next_p][1] - x, shape[next_p][2] - y
            cross = (Δx_1 * Δy_2 - Δx_2 * Δy_1) / sqrt((Δx_1^2 + Δy_1^2) * (Δx_2^2 + Δy_2^2))

            θ += asin(cross)
        end

        winding_number = round(θ / (2 * pi))
        if (winding_number == 0.0 && in_or_out == :out) || (winding_number != 0.0 && in_or_out == :in)
            push!(indices_to_remove, n)
        end
    end
    deleteat!(positions, indices_to_remove)
end

"""
    nearest_neighbor_pairs(particles; maximum_distance)

Returns pairs of particles that are within `maximum_distance` of eachother.
"""
function nearest_neighbor_pairs(particles::Array{Particle, 1}; maximum_distance::Float64 = 0.0, period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    pairs = Array{NTuple{2, Particle}, 1}()

    N = length(particles)
    for i = 1 : N, j = i + 1 : N
        Δx = wrap_displacement(particles[i].x - particles[j].x; period = period_x)
        Δy = wrap_displacement(particles[i].y - particles[j].y; period = period_y)

        if Δx^2 + Δy^2 < maximum_distance^2
            push!(pairs, (particles[i], particles[j]))
        end
    end

    return pairs
end


text:
# ------------------------------------------------------------------------------------------
# # Orbit Diagrams, Surfaces of Section
# In this tutorial we will be using a `DynamicalSystem` instance to visualize aspects of the
# system.
#
# Topics:
# * `orbitdiagram`
# * `poincaresos`
# * `produce_orbitdiagram`
# * Docstrings
#
# # Orbit Diagram
# An "orbit diagram" is simply a plot that shows the long term behavior of a discrete system
# when a parameter is varied.
#
# This is exactly what the function `orbitdiagram` does:
# 1. Evolves the system for a transient amount of time.
# 2. Evolves & saves the output of the system for a chosen variable.
# 3. Changes/increments a parameter of the equations of motion.
# 4. Repeat steps 1-3 for all given parameter values!
#
# This approach is also related with the (much more advanced) [sensitivity
# analysis](http://docs.juliadiffeq.org/latest/analysis/sensitivity.html) from
# DifferentialEquations.jl
#
# ---
#
# Let's make the super-ultra-famous orbit diagram of the logistic map:
#
# $$x_{n+1} = rx_n(1-x_n)$$
# ------------------------------------------------------------------------------------------

using DynamicalSystems, PyPlot

logimap = Systems.logistic() # Systems module contains pre-defined well-known systems

# ------------------------------------------------------------------------------------------
# ---
#
# The call signature of `orbitdiagram` is:
#
# ```julia
# orbitdiagram(discrete_system, i, p_index, pvalues; n, Ttr, ...)
# ```
# * `i` is the index of the variable we want to save.
# * `p_index` is the index of the parameter we want to change.
# * `pvalues` are the values of the parameter that will change.
# * Keywords `Ttr` and `n` denote for how much transient time to evolve the system and how
# many states to save.
#
# ------------------------------------------------------------------------------------------

i = 1 # which variable to save (of course logistic map is 1D !)
n = 50 # how many values to save
Ttr = 5000 # transient iterations
p_index = 1
pvalues = 2:0.001:4  # parameter values
output = orbitdiagram(logimap, i, p_index, pvalues; n = n, Ttr = Ttr)
typeof(output)

# ------------------------------------------------------------------------------------------
# * The output is a vector of vectors. Each inner vector has length `n` and contains the
# values of the variable `i` at the given parameter value.
#
# Let's plot interactively!
# ------------------------------------------------------------------------------------------

function bf(pvalues, n, Ttr)
    logi = Systems.logistic()
    output = orbitdiagram(logi, 1, 1, pvalues; n = n, Ttr = Ttr)

    figure(figsize=(6,4))
    for (j, p) in enumerate(pvalues)
        plot(p .* ones(output[j]), output[j], linestyle = "None", # linestyle = None
        marker = "o", ms = 0.2, color = "black")
    end
    xlabel("\$r\$"); ylabel("\$x\$");
    xlim(pvalues[1], pvalues[end])
    # ylim(0, 1)
    return
end

bf(linspace(2.0, 4.0, 1000), 200, 2000)

bf(linspace(3.5, 3.6, 1000), 200, 2000)

# ------------------------------------------------------------------------------------------
# We can use `BenchmarkTools` to see how much time it takes to generate this orbit diagram
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@btime output = orbitdiagram($logimap, $i, $p_index, $pvalues; n = $n, Ttr = $Ttr);
println("for total points: $(length(pvalues)*(Ttr+n)), out of which $(length(pvalues)*n) are saved")

# ------------------------------------------------------------------------------------------
# ---
#
# * `orbitdiagram` works with *any* discrete system! Check out the [documentation
# page](https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/orbitdiagram/) for
# more!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Poincaré Surface of Section
# This is a technique to reduce a continuous system into a discrete map with 1 fewer
# dimension.
# The wikipedia entry on [Poincaré map](https://en.wikipedia.org/wiki/Poincar%C3%A9_map) has
# a lot of useful info, but the technique itself is very simple:
#
# 1. Define a hyperplane in the phase-space of the system.
# 2. Evolve the continuous system for long times. Each time the trajectory crosses this
# plane, record the state of the system.
# 3. Only crossings with a specific `direction` (either positive or negative) are allowed.
#
# And that's it! The recorded crossings are the Poincaré Surface of Section!
#
# ## Defining a hyperplane
# Let's say that our phase-space is $D$ dimensional. If the state of the system is
# $\mathbf{u} = (u_1, \ldots, u_D)$ then the equation for a hyperplane is
#
# $$
# a_1u_1 + \dots + a_Du_D = \mathbf{a}\cdot\mathbf{u}=b
# $$
# where $\mathbf{a}, b$ are the parameters that define the hyperplane.
#
# ---
#
# Here is the call signature for a function that does this:
#
# ```julia
# poincaresos(continuous_system, plane, tfinal = 100.0; direction = 1, ...)
# ```
# In code, `plane` can be either:
#
# * A `Tuple{Int, <: Number}`, like `(j, r)` : the hyperplane is defined as when the `j`
# variable of the system crosses the value `r`.
# * An `AbstractVector` of length `D+1`. The first `D` elements of the vector correspond to
# $\mathbf{a}$ while the last element is $b$. The hyperplane is defined with its formal
# equation.
#
# ---
#
# As an example, let's see a section of the Lorenz system:
# $$
# \begin{aligned}
# \dot{X} &= \sigma(Y-X) \\
# \dot{Y} &= -XZ + \rho X -Y \\
# \dot{Z} &= XY - \beta Z
# \end{aligned}
# $$
#
# ------------------------------------------------------------------------------------------

lor = Systems.lorenz()

tr = trajectory(lor, 100.0, dt = 0.01)
figure(figsize = (8,6))
plot3D(columns(tr)...);
xlabel("X"); ylabel("Y"); zlabel("Z");

psos = poincaresos(lor, (2, 0.0), 2000.0) # find where 2nd variable crosses 0.0

figure(figsize = (6,4))
plot(psos[:, 1], psos[:, 3], lw=0.0, marker ="o", ms = 1.0, color = "C1");
xlabel("X"); ylabel("Z");

# ------------------------------------------------------------------------------------------
# * We see that the surface of section is some kind of 1-dimensional object.
# * This is expected, because as we will show in the tutorial "Entropies & Dimensions" the
# Lorenz system (at least for the default parameters) lives in an almost 2-dimensional
# attractor.
#
# * This means that when you take a cut through this object, the result should be
# 1-dimensional!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# **Next, let's visualize the Poincaré Surface of Section in 3D**
# ------------------------------------------------------------------------------------------

"""
    meshgrid(x, y) -> X, Y
Create a meshgrid, such as that used in PyPlot's countour plots.
"""
function meshgrid(vx, vy)
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

# ------------------------------------------------------------------------------------------
# We can get the attractor by calculating trajectories:
# ------------------------------------------------------------------------------------------

tr = trajectory(lor, 200.0, dt=0.01)
x, y, z = columns(tr);

# ------------------------------------------------------------------------------------------
# Next, let's make sure we can color points on our 3D plot based on where they are relative
# to a cut we'll make in the plot.
#
# To do this, we'll create a vector, `c`, that will store labels for the colors of all
# points:
# ------------------------------------------------------------------------------------------

c = Vector{String}(length(y))
for i in 1:length(y) # cut points: red
    if -0.1 < y[i] < 0.1
        c[i] = "C3"
    elseif y[i] < 0
        c[i] = "C0" # in front of cut: blue
    else
        c[i] = "C2" # behind cut: green
    end
end

figure(figsize = (8,6))
# First let's plot the attractor
plot3D(x,y,z, color = "black", lw = 0.2, alpha = 0.25)
scatter3D(x, y, z, color = c, s = 5.0)

# And then plot the PSOS plane:
x = [-20, 20]; z = [0, 40];
X, Z = meshgrid(x, z)
Y = zeros(X)

plot_surface(X, Y, Z, alpha = 0.25, color = "C1");
xlabel("X"); ylabel("Y"); zlabel("Z");

# ------------------------------------------------------------------------------------------
# # Producing an orbit diagram
#
# 1. The `orbitdiagram` function does not make sense for continuous systems. In order for
# `orbitdiagram` to have meaning one must have a map.
#
# 2. We can take advantage of the `poincaresos` function, and reduce a continuous system to
# a map.
# 3. Then, we can formally calculate an orbit diagram for a continuous system!
#
# In this example I will use the Shinriki oscillator, which displays a period-doubling route
# to chaos like the logistic map!
# ------------------------------------------------------------------------------------------

shinriki_voltage(V) = 2.295e-5*(exp(3.0038*V) - exp(-3.0038*V))
function shinriki_eom(u, p, t)
    R1 = p[1]

    du1 = (1/0.01)*(
    u[1]*(1/6.9 - 1/R1) - shinriki_voltage(u[1] - u[2]) - (u[1] - u[2])/14.5
    )

    du2 = (1/0.1)*(
    shinriki_voltage(u[1] - u[2]) + (u[1] - u[2])/14.5 - u[3]
    )

    du3 = (1/0.32)*(-u[3]*0.1 + u[2])
    return SVector{3}(du1, du2, du3)
end

# Be sure to give a mutable container for the parameter container!
shi = ContinuousDynamicalSystem(shinriki_eom, [-2, 0, 0.2], [22.0])

# ------------------------------------------------------------------------------------------
# Now we can use `shi`, our dynamical system, to calculate trajectory steps and then plot
# the result.
# ------------------------------------------------------------------------------------------

set_parameter!(shi, 1, 20.0) # parameter from 19 to 22
tr = trajectory(shi, 2000.0)
figure(figsize = (8,6))
plot3D(columns(tr[1000:end, :])..., color = "C3", alpha = 0.5, marker = "o", ms = 0.5);

# ------------------------------------------------------------------------------------------
# To get a feeling for the system, let's look at a couple `poicaresos`
# ------------------------------------------------------------------------------------------

figure(figsize=(8,4))
subplot(1,2,1)
# the function set_parameter! is useful here! (see docstring!)
R1 = 19.5
set_parameter!(shi, 1, R1)

# here I use `direction = -1` (makes more sense for the specific system)
psos = poincaresos(shi, (2, 0.0), 1000.0, Ttr = 100.0, direction = -1)

plot(psos[:, 1], psos[:, 3], lw=0.0, marker ="o", ms = 2.0, color = "C3");
xlabel("\$V_1\$"); ylabel("\$I_3\$")
title("R1 = $R1")

R1 = 21.0
set_parameter!(shi, 1, R1)
subplot(1,2,2)
psos = poincaresos(shi, (2, 0.0), 1000.0, Ttr = 100.0, direction = -1)
plot(psos[:, 1], psos[:, 3], lw=0.0, marker ="o", ms = 2.0, color = "C3");
xlabel("\$V_1\$");
title("R1 = $R1");

# ------------------------------------------------------------------------------------------
# Hm, this is interesting!
# 1. For some parameters the motion is clearly periodic (due to the distinct number of
# points)
# 2. But at other parameters the motion ***seems to be*** on a 2-dimensional manifold.
#
# **Don't be too quick to judge the second as chaotic though! In the next tutorial
# "Quantifying Chaos" I will show you tools to quantify chaotic behavior!**
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ---
#
# * We would like to see the behavior of the system while varying the parameter, like an
# orbit diagram.
#
# * We can do this by performing successive surfaces of section and then recording the
# values of a chosen variable!
#
# This procedure is bundled in the very convenient function `produce_orbitdiagram`:
#
# ```julia
# produce_orbitdiagram(continuous_system, plane, i, p_index, pvalues; kwargs...)
# ```
#
# The function produces successive `poincaresos` for the `plane`, and records the values of
# the `i`-th variable at the section.
#
# Most other arguments are simply propagated to `poincaresos`.
#
# *Lets use this function for the Shinriki oscillator example:*
#
# First we can store our input arguments for `produce_orbitdiagram` in a few input vars.
# ------------------------------------------------------------------------------------------

pvalues = linspace(19,22,201) # which parameter values to use
p_index = 1 # which parameter to change

i = 1                  # record 1st variable
plane = (2, 0.0)       # find PSOS of 2nd variable when crossing zero
tf = 200.0             # argument passed to poincaresos

# ------------------------------------------------------------------------------------------
# Now we're ready to call `produce_orbitdiagram` and plot the result!
#
# (NOTE: This computation takes some seconds, since it makes 201 PSOS!)
# ------------------------------------------------------------------------------------------

output = produce_orbitdiagram(shi, plane, i, p_index, pvalues;
                              # all keywords are passed to poincaresos:
                              tfinal = tf, Ttr = 200.0, direction = -1)

figure(figsize=(6,4))
for (j, p) in enumerate(pvalues)
    plot(p .* ones(output[j]), output[j], ls = "None",
    marker = "o", ms = 0.5, color = "black")
end
xlabel("\$R_1\$"); ylabel("\$V_1\$");

# ------------------------------------------------------------------------------------------
# ---
#
# ## Example with parameter-dependent plane
#
# Let's see one more case, using the Gissinger system:
#
# $$
# \begin{aligned}
# \dot{Q} &= \mu Q - VD \\
# \dot{D} &= -\nu D + VQ \\
# \dot{V} &= \Gamma -V + QD
# \end{aligned}
# $$
#
# ------------------------------------------------------------------------------------------

gis = Systems.gissinger([2.32865, 2.02514, 1.98312]) # initial condition necessary to see structure

set_parameter!(gis, 1, 0.112) # parameter from 0.1 to 0.15
tr = trajectory(gis, 1000.0)
figure(figsize = (8,6))
plot3D(columns(tr)...);

pvalues = linspace(0.1,0.15,201) # which parameter values to use
p_index = 1 # change 2nd parameter, which is μ

i = 2            # record i variable
j = 1            # find PSOS of j variable
offset = -0.75   # offset is important here

tf = 5000.0

output = produce_orbitdiagram(gis, (j, offset), i, p_index, pvalues;
                              # all keywords are passed to poincaresos:
                              tfinal = tf, Ttr = 2000.0)

figure(figsize=(8,6))
for (j, p) in enumerate(pvalues)
    length(output[j]) == 0 && continue
    plot(p .* ones(output[j]), output[j], ls = "None",
    marker = "o", ms = 0.5, color = "black", alpha = 0.5)
end
ylim(1,3)
xlabel("\$\\mu\$"); ylabel("\$D\$");

# ------------------------------------------------------------------------------------------
# ---
#
# 1. The above is quite cool, but for this specific system, the "optimal" hyperplane to use
# depends on the parameter `μ`.
#
# 2. In addition, it is not "optimal" to record a specific variable during the crossing, but
# instead a function of the variables.
#
# 3. For this approach, it is better to get the entire `poincaresos`s for each parameter,
# like I show in the following example:
# ------------------------------------------------------------------------------------------

# Define appropriate hyperplane for gissinger system
const ν = 0.1
const Γ = 0.9 # default parameters of the system

# I want hyperperplane defined by these two points:
Np(μ) = SVector{3}(sqrt(ν + Γ*sqrt(ν/μ)), -sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))
Nm(μ) = SVector{3}(-sqrt(ν + Γ*sqrt(ν/μ)), sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))

# Create hyperplane using normal vector to vector connecting points:
gis_plane(μ) = (d = (Np(μ) - Nm(μ)); [d[2], -d[1], 0, 0])

gis_plane(0.112)

μ = 0.12
set_parameter!(gis, 1, μ)
figure(figsize = (8,6))
psos = poincaresos(gis, gis_plane(μ), 2000.0, Ttr = 2000.0)
plot3D(columns(psos)..., marker = "o", ls = "None", ms = 2.0);

pvalues = linspace(0.1,0.14,201) # which parameter values to use

xs = Vector{Float64}[] # empty vector of vectors

for (i, μ) in enumerate(pvalues)
    set_parameter!(gis, 1, μ)
    N = Np(μ)
    psos = poincaresos(gis, gis_plane(μ), 4000.0, Ttr = 2000.0)
    push!(xs, [norm(N - k) for k in psos])
end

figure(figsize=(8,6))
for (j, p) in enumerate(pvalues)
    length(xs[j]) == 0 && continue
    plot(p .* ones(xs[j]), xs[j], ls = "None",
    marker = "o", ms = 0.2, color = "black", alpha = 0.5)
end
xlabel("\$\\mu\$"); ylabel("\$x\$");

figure(figsize=(8,6))
for (j, p) in enumerate(pvalues)
    length(xs[j]) == 0 && continue
    plot(p .* ones(xs[j]), xs[j], ls = "None",
    marker = "o", ms = 0.2, color = "black", alpha = 0.5)
end
ylim(2.7, 2.9)
xlabel("\$\\mu\$"); ylabel("\$x\$");

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?orbitdiagram

?poincaresos

?produce_orbitdiagram

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/introduction-to-dynamicalsystems.jl/2. Orbit diagrams, Surfaces of Section.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
    Hamerly()

Hamerly algorithm implementation, based on "Hamerly, Greg. (2010). Making k-means Even Faster.
 Proceedings of the 2010 SIAM International Conference on Data Mining. 130-140. 10.1137/1.9781611972801.12."

This algorithm provides much faster convergence than Lloyd algorithm with realtively small increase in
memory footprint. It is especially suitable for low to medium dimensional input data.

It can be used directly in `kmeans` function

```julia
X = rand(30, 100_000)   # 100_000 random points in 30 dimensions

kmeans(Hamerly(), X, 3) # 3 clusters, Hamerly algorithm
```
"""
struct Hamerly <: AbstractKMeansAlg end

function kmeans!(alg::Hamerly, containers, X, k;
                n_threads = Threads.nthreads(),
                k_init = "k-means++", max_iters = 300,
                tol = 1e-6, verbose = false, init = nothing)
    nrow, ncol = size(X)
    centroids = init == nothing ? smart_init(X, k, n_threads, init=k_init).centroids : deepcopy(init)

    @parallelize n_threads ncol chunk_initialize(alg, containers, centroids, X)

    converged = false
    niters = 0
    J_previous = 0.0
    p = containers.p

    # Update centroids & labels with closest members until convergence
    while niters < max_iters
        niters += 1
        update_containers(alg, containers, centroids, n_threads)
        @parallelize n_threads ncol chunk_update_centroids(alg, containers, centroids, X)
        collect_containers(alg, containers, n_threads)

        J = sum(containers.ub)
        move_centers(alg, containers, centroids)

        r1, r2, pr1, pr2 = double_argmax(p)
        @parallelize n_threads ncol chunk_update_bounds(alg, containers, r1, r2, pr1, pr2)

        if verbose
            # Show progress and terminate if J stops decreasing as specified by the tolerance level.
            println("Iteration $niters: Jclust = $J")
        end

        # Check for convergence
        if (niters > 1) & (abs(J - J_previous) < (tol * J))
            converged = true
            break
        end

        J_previous = J
    end

    @parallelize n_threads ncol sum_of_squares(containers, X, containers.labels, centroids)
    totalcost = sum(containers.sum_of_squares)

    # Terminate algorithm with the assumption that K-means has converged
    if verbose & converged
        println("Successfully terminated with convergence.")
    end

    # TODO empty placeholder vectors should be calculated
    # TODO Float64 type definitions is too restrictive, should be relaxed
    # especially during GPU related development
    return KmeansResult(centroids, containers.labels, Float64[], Int[], Float64[], totalcost, niters, converged)
end

function collect_containers(alg::Hamerly, containers, n_threads)
    if n_threads == 1
        @inbounds containers.centroids_new[end] .= containers.centroids_new[1] ./ containers.centroids_cnt[1]'
    else
        @inbounds containers.centroids_new[end] .= containers.centroids_new[1]
        @inbounds containers.centroids_cnt[end] .= containers.centroids_cnt[1]
        @inbounds for i in 2:n_threads
            containers.centroids_new[end] .+= containers.centroids_new[i]
            containers.centroids_cnt[end] .+= containers.centroids_cnt[i]
        end

        @inbounds containers.centroids_new[end] .= containers.centroids_new[end] ./ containers.centroids_cnt[end]'
    end
end

function create_containers(alg::Hamerly, k, nrow, ncol, n_threads)
    lng = n_threads + 1
    centroids_new = Vector{Array{Float64,2}}(undef, lng)
    centroids_cnt = Vector{Vector{Int}}(undef, lng)

    for i = 1:lng
        centroids_new[i] = zeros(nrow, k)
        centroids_cnt[i] = zeros(k)
    end

    # Upper bound to the closest center
    ub = Vector{Float64}(undef, ncol)

    # lower bound to the second closest center
    lb = Vector{Float64}(undef, ncol)

    labels = zeros(Int, ncol)

    # distance that centroid moved
    p = Vector{Float64}(undef, k)

    # distance from the center to the closest other center
    s = Vector{Float64}(undef, k)

    # total_sum_calculation
    sum_of_squares = Vector{Float64}(undef, n_threads)

    return (
        centroids_new = centroids_new,
        centroids_cnt = centroids_cnt,
        labels = labels,
        ub = ub,
        lb = lb,
        p = p,
        s = s,
        sum_of_squares = sum_of_squares
    )
end

"""
    chunk_initialize(alg::Hamerly, containers, centroids, design_matrix, r, idx)

Initial calulation of all bounds and points labeling.
"""
function chunk_initialize(alg::Hamerly, containers, centroids, X, r, idx)
    centroids_cnt = containers.centroids_cnt[idx]
    centroids_new = containers.centroids_new[idx]

    @inbounds for i in r
        label = point_all_centers!(containers, centroids, X, i)
        centroids_cnt[label] += 1
        for j in axes(X, 1)
            centroids_new[j, label] += X[j, i]
        end
    end
end

"""
    update_containers(::Hamerly, containers, centroids, n_threads)

Calculates minimum distances from centers to each other.
"""
function update_containers(::Hamerly, containers, centroids, n_threads)
    s = containers.s
    s .= Inf
    @inbounds for i in axes(centroids, 2)
        for j in i+1:size(centroids, 2)
            d = distance(centroids, centroids, i, j)
            d = 0.25*d
            s[i] = s[i] > d ? d : s[i]
            s[j] = s[j] > d ? d : s[j]
        end
    end
end

"""
    chunk_update_centroids(::Hamerly, containers, centroids, X, r, idx)

Detailed description of this function can be found in the original paper. It iterates through
all points and tries to skip some calculation using known upper and lower bounds of distances
from point to centers. If it fails to skip than it fall back to generic `point_all_centers!` function.
"""
function chunk_update_centroids(alg::Hamerly, containers, centroids, X, r, idx)

    # unpack containers for easier manipulations
    centroids_new = containers.centroids_new[idx]
    centroids_cnt = containers.centroids_cnt[idx]
    labels = containers.labels
    s = containers.s
    lb = containers.lb
    ub = containers.ub

    @inbounds for i in r
        # m ← max(s(a(i))/2, l(i))
        m = max(s[labels[i]], lb[i])
        # first bound test
        if ub[i] > m
            # tighten upper bound
            label = labels[i]
            ub[i] = distance(X, centroids, i, label)
            # second bound test
            if ub[i] > m
                label_new = point_all_centers!(containers, centroids, X, i)
                if label != label_new
                    labels[i] = label_new
                    centroids_cnt[label_new] += 1
                    centroids_cnt[label] -= 1
                    for j in axes(X, 1)
                        centroids_new[j, label_new] += X[j, i]
                        centroids_new[j, label] -= X[j, i]
                    end
                end
            end
        end
    end
end

"""
    point_all_centers!(containers, centroids, X, i)

Calculates new labels and upper and lower bounds for all points.
"""
function point_all_centers!(containers, centroids, X, i)
    ub = containers.ub
    lb = containers.lb
    labels = containers.labels

    min_distance = Inf
    min_distance2 = Inf
    label = 1
    @inbounds for k in axes(centroids, 2)
        dist = distance(X, centroids, i, k)
        if min_distance > dist
            label = k
            min_distance2 = min_distance
            min_distance = dist
        elseif min_distance2 > dist
            min_distance2 = dist
        end
    end

    ub[i] = min_distance
    lb[i] = min_distance2
    labels[i] = label

    return label
end

"""
    move_centers(::Hamerly, containers, centroids)

Calculates new positions of centers and distance they have moved. Results are stored
in `centroids` and `p` respectively.
"""
function move_centers(::Hamerly, containers, centroids)
    centroids_new = containers.centroids_new[end]
    p = containers.p

    @inbounds for i in axes(centroids, 2)
        d = 0.0
        for j in axes(centroids, 1)
            d += (centroids[j, i] - centroids_new[j, i])^2
            centroids[j, i] = centroids_new[j, i]
        end
        p[i] = d
    end
end

"""
    chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, r, idx)

Updates upper and lower bounds of point distance to the centers, with regard to the centers movement.
"""
function chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, r, idx)
    p = containers.p
    ub = containers.ub
    lb = containers.lb
    labels = containers.labels

    # Since bounds are squred distance, `sqrt` is used to make corresponding estimation, unlike
    # the original paper, where usual metric is used.
    #
    # Using notation from original paper, `u` is upper bound and `a` is `labels`, so
    #
    # `u[i] -> u[i] + p[a[i]]`
    #
    # then squared distance is
    #
    # `u[i]^2 -> (u[i] + p[a[i]])^2 = u[i]^2 + 2 p[a[i]] u[i] + p[a[i]]^2`
    #
    # Taking into account that in our noations `p^2 -> p`, `u^2 -> ub` we obtain
    #
    # `ub[i] -> ub[i] + 2 sqrt(p[a[i]] ub[i]) + p[a[i]]`
    #
    # The same applies to the lower bounds.
    @inbounds for i in r
        label = labels[i]
        ub[i] += 2*sqrt(abs(ub[i] * p[label])) + p[label]
        if r1 == label
            lb[i] += pr2 - 2*sqrt(abs(pr2*lb[i]))
        else
            lb[i] += pr1 - 2*sqrt(abs(pr1*lb[i]))
        end
    end
end

"""
    double_argmax(p)

Finds maximum and next after maximum arguments.
"""
function double_argmax(p)
    r1, r2 = 1, 1
    d1 = p[1]
    d2 = -1.0
    for i in 2:length(p)
        if p[i] > d1
            r2 = r1
            r1 = i
            d2 = d1
            d1 = p[i]
        elseif p[i] > d2
            d2 = p[i]
            r2 = i
        end
    end

    r1, r2, d1, d2
end


text:
"""
	earthtide(cmd0::String=""; kwargs...)

Compute grids or time-series of solid Earth tides.

Full option list at [`earthtide`]($(GMTdoc)earthtide.html)

```julia
	G = earthtide();
	imshow(G)
```
"""
function earthtide(cmd0::String=""; kwargs...)

	(cmd0 != "") && return monolitic("earthtide", cmd0)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd = "earthtide " * parse_common_opts(d, "", [:R :I :V_params :r])[1]
	cmd = parse_opt_range(d, cmd, "T")[1]
	if ((opt_S = add_opt(d, "", 'S', [:S :sun_moon])) != "")
		return finish_PS_module(d, cmd * opt_S, "", true, false, false)
	elseif ((opt_L = add_opt(d, "", 'L', [:L :location])) != "")
		return finish_PS_module(d, cmd * opt_L, "", true, false, false)
	end

	cmd = ((opt_C = add_opt(d, "", 'C', [:C :components])) != "") ? cmd * opt_C : cmd * " -Cz"
	opt_G = add_opt(d, "", 'G', [:G :grid :outgrid])
	(length(opt_G) > 3) && (cmd *= opt_G)		# G=true will give " -G", which we'll ignore  (Have to)

	return (dbg_print_cmd(d, cmd) !== nothing) ? cmd : gmt(cmd)
end


text:
include("DME_solver.jl")

#INPUT IS BASED ON EXPERIMENT 1

perm_val  = 1e-15 # [m^2] permeability
poros_val = 0.256     # [-] porosity
ρ_rock    = 2700     # [kg/m^3]

Kval = c -> 2.0
dKval_dc = c -> zeros(size(c))
a = 2.0*2700.0*1000.0*(1.0-poros_val) # m^2/m^3 = 2 m^2/g * 2700 kg/m^3 * 1000 g/kg * (1-phi)
D_oil = 0.5*10^(-9.0) # m^2/s #0.5e-11#
D_water = 2*10^(-9.0) # m^2/s #2e-11#

## Creating DME parameters
ρ_oil = c_o -> 748.0015057043045-0.0016431963703508173.*c_o+2.9364412314092913.*10.0.^(-7.0).*c_o.^2.0-4.107225978803957.*10.0.^(-11.0).*c_o.^3.0
dρ_dc_oil = c_o -> -0.0016431963703508173+2.0.*2.9364412314092913.*10.0.^(-7.0).*c_o-3.0.*4.107225978803957.*10.0.^(-11.0).*c_o.^2.0
dρ_dp_oil = c_o -> zeros(size(c_o))
ρ_water = c_w -> 983.6319389796079.*ones(size(c_w))
dρ_dp_water = c_w -> zeros(size(c_w))
μ_oil = c_o -> 0.0010131130815849523-2.528969514668829.*10.0.^(-5.0).*c_o+2.3873417406201346.*10.0.^(-11.0).*c_o.^2.0
dμ_dc_oil = c_o -> -2.528969514668829.*10.0.^(-5.0)+2*2.3873417406201346.*10.0.^(-11.0).*c_o
μ_water = c_w -> 0.00040706344808758737+1.47269029810352587.*10.0.^(-7.0).*c_w
dμ_dc_water = c_w -> 1.47269029810352587.*10.0.^(-7.0).*ones(size(c_w))

MT = FluidParametersDME(ρ_oil, dρ_dc_oil, dρ_dp_oil, ρ_water, dρ_dp_water, μ_oil, dμ_dc_oil, μ_water, dμ_dc_water, Kval, dKval_dc, a, D_oil, D_water)

#M_DME = 46.07/1000 #kg/mol
#M_W = 18.01528/1000 #kg/mol
#M_avg = (M_DME+M_W)/2
#vol = 0.104*pi*0.025^2.0/4
#u_inj = 0.00011/(3600*24) # m/day to m/s injection velocity
D_core = 0.025
A_core  = π*D_core^2/4
injection_rate = 4.5202702702702705e-10 # 0.00011/(3600*24) # m^3/day
u_inj   = injection_rate/A_core # Darcy velocity
c_inj1 = 0.0 # mol/m^3 DME in injection water at initial water flooding
c_inj2 = 10#0.1*983.6319389796079/M_avg # mol/m^3 DME in injection waterat DME flooding
c_tracer1 = 0.0
c_tracer2 = 1.0

# Creating realitive permeability
#swc  = 0.392349
#kro0 = 0.121566
#krw0 = 0.269128
#n_w  = 1.88088
#n_o  = 1.91705
(sor_max, swc, kro0min, krw0, n_o, n_w) = [0.198781  0.205555  0.249794  0.149609  3.03704  3.22222]
kro0max = 0.95
sor_min  = 0.01 # equivalent of equilibrium mass frac with 0.1 mol frac DME-water

#sor_min  = 0.03
#sor_max  = 0.250158
m_kro = -(kro0max-kro0min)/(sor_max-sor_min)
kro0_c = c -> m_kro*c + kro0max

c_min = 0.0 # min DME mol/m^3 in oil
c_max = Kval(0)*c_inj2 # max DME mol/m^3 fraction in oil
c_max = c_max[1]
m_sor = -(sor_max-sor_min)/(c_max-c_min)
sor_c = c -> m_sor*c+sor_max
dsor_c = c -> m_sor + 0.0*c #irrelevant

K = RelativePermeabilityDME(kro0_c, krw0, n_o, n_w, sor_c, dsor_c, swc)

# =============================================================================
# Preparing geology/mesh
# =============================================================================
dim = 2

srand(123456543)
λ_dp    = 0.8
cor_x   = 0.1
cor_y   = 0.1
cor_z   = 0.1
Lx      = 0.104# [m]
Ly      = 0.025 # [m]
Lz      = 0.025 # [m]
Nx      = 50 # number of grids in the x direction
Ny      = 12 # number of grids in the y direction
Nz      = 10 # number of grids in the z direction

if dim == 1
    m           = createMesh1D(Nx, Lx)
    perm_mat    = permfieldlogrnde(Nx, perm_val, λ_dp, cor_x)
elseif dim == 2
    m           = createMesh2D(Nx, Ny, Lx, Ly)
    perm_mat    = permfieldlogrnde(Nx, Ny, perm_val, λ_dp, cor_x, cor_y)
else
    m           = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
    perm_mat    = permfieldlogrnde(Nx, Ny, Nz, perm_val, λ_dp, cor_x, cor_y, cor_z)
end

k       = createCellVariable(m, perm_mat)
#k       = createCellVariable(m, perm_val)
k_face  = harmonicMean(k)     # permeability on the cell faces
ϕ       = createCellVariable(m, poros_val)

G = Geology(m, perm_val, poros_val, k, ϕ)

# =============================================================================
# Calculating heuristic for k_o and k_w
# =============================================================================
d = sqrt(perm_val) #Characteristic length of capillary in porous medium
c_o = Kval(0)*c_inj2
c_o = c_o[1]
c_w = c_inj2
v_o = u_inj/poros_val
v_w = u_inj/poros_val

k_o = 0.023*D_oil/d*((ρ_oil(c_o)*v_o*d)/μ_oil(c_o))^(0.83)*(μ_oil(c_o)/(ρ_oil(c_o)*D_oil))^(1.0/3.0)
k_w = 0.023*D_water/d*((ρ_water(c_w)[1]*v_w*d)/μ_water(c_w))^(0.83)*(μ_water(c_w)/(ρ_water(c_w)[1]*D_water))^(1.0/3.0)

k_o = k_o*1e-3 #optimal solution
k_w = k_w*1e-3 #optimal solution

# =============================================================================
# Preparing time-stepping parameters
# =============================================================================
n_pv_water1     = 9.0 # number of injected pore volumes
t_final_water1  = n_pv_water1*Lx/(u_inj/poros_val) # [s] final time for water flood
dt_water1       = t_final_water1/n_pv_water1/Nx/5 # [s] time step

n_pv_DME        = 6.0#9.0
t_final_DME     = t_final_water1+n_pv_DME*Lx/(u_inj/poros_val) # [s] final time for DME flood
dt_DME          = t_final_DME/(n_pv_water1 + n_pv_DME)/Nx/5 # [s] time step

n_pv_water2     = 4.0
t_final_water2  = t_final_DME+n_pv_water2*Lx/(u_inj/poros_val) # [s] final time for DME flood
dt_water2       = t_final_water2/(n_pv_water1 + n_pv_DME + n_pv_water2)/Nx/5 # [s] time step

T_water1 = SimulationTime(t_final_water1, dt_water1)
T_DME    = SimulationTime(t_final_DME, dt_DME)
T_water2 = SimulationTime(t_final_water2, dt_water2)

T = [T_water1, T_DME]#, T_water2]

## Boundary Conditions for water flooding
#T0 = 70 + 273.15 # [K]

c_oil0      = 0.0 # initial DME concentration in the oil phase
c_water0    = 0.0 # initial DME concentration in the water phase
p0          = 2000/14.7*1e5   # [Pa]
sw0         = 1.0-0.667 # [vol frac] initial water saturation
c_t0        = 0.0 # initial tracer concentration

BCc_oil_water     = createBC(m) # concentration (DME in oil) boundary condition
BCc_water_water   = createBC(m) # concentration (DME in water) boundary condition
BCp_water         = createBC(m) # pressure boundary condition
BCs_water         = createBC(m) # saturation boundary condition
BCt_water         = createBC(m) # tracer concentration boundary condition

BCp_water.right.a[:] = 0.0
BCp_water.right.b[:] = 1.0 #Dirichlet
BCp_water.right.c[:] = p0

BCp_water.left.a[:]  = k_face.xvalue[1,:]./μ_water(c_inj1) #Neumann
BCp_water.left.b[:]  = 0.0
BCp_water.left.c[:]  = -u_inj

BCs_water.left.a[:]  = 0.0
BCs_water.left.b[:]  = 1.0 #Dirichlet
BCs_water.left.c[:]  = 1.0 #Value S_w=1

BCc_water_water.left.a[:]  = 0.0 #
BCc_water_water.left.b[:]  = 1.0 #Dirichlet
BCc_water_water.left.c[:]  = c_inj1 #injected concentration of DME in water

BCc_oil_water.left.a[:]  = 0.0 #
BCc_oil_water.left.b[:]  = 1.0 #Dirichlet
BCc_oil_water.left.c[:]  = Kval(0)*c_inj1 #injected concentration of DME in water

BCt_water.left.a[:]  = 0.0
BCt_water.left.b[:]  = 1.0 #Dirichlet
BCt_water.left.c[:]  = c_tracer1

BC_water = BoundaryConditionDME(BCc_oil_water, BCc_water_water, BCp_water, BCs_water, BCt_water)

c_oil_init      = createCellVariable(m, c_oil0, BCc_oil_water)
c_water_init    = createCellVariable(m, c_water0, BCc_water_water)
p_init          = createCellVariable(m, p0, BCp_water)
sw_init         = createCellVariable(m, sw0, BCs_water)
c_t_init        = createCellVariable(m, c_t0, BCt_water)

IC = InitialConditionDME(c_oil_init, c_water_init, p_init, sw_init, c_t_init)

# Creating Boundary conditions for DME flooding
BCc_oil_DME     = createBC(m) # concentration (DME in oil) boundary condition
BCc_water_DME   = createBC(m) # concentration (DME in water) boundary condition
BCp_DME         = createBC(m) # pressure boundary condition
BCs_DME         = createBC(m) # saturation boundary condition
BCt_DME         = createBC(m) # tracer concentration boundary condition

BCc_water_DME.left.a[:] = 0.0
BCc_water_DME.left.b[:] = 1.0 #Dirichlet
BCc_water_DME.left.c[:] = c_inj2 #Injected concentration of DME

BCc_oil_DME.left.a[:] = 0.0
BCc_oil_DME.left.b[:] = 1.0 #Dirichlet
BCc_oil_DME.left.c[:] = Kval(0)*c_inj2 #Created from reaction at boundary

BCp_DME.right.a[:] = 0.0
BCp_DME.right.b[:] = 1.0 #Dirichlet
BCp_DME.right.c[:] = p0

BCp_DME.left.a[:]  = k_face.xvalue[1,:]/μ_water(c_inj2) #Neumann
BCp_DME.left.b[:]  = 0.0
BCp_DME.left.c[:]  = -u_inj

BCs_DME.left.a[:]  = 0.0
BCs_DME.left.b[:]  = 1.0 #Dirichlet
BCs_DME.left.c[:]  = 1.0 #Value S_w=1

BCt_DME.left.a[:] = 0.0
BCt_DME.left.b[:] = 1.0 #Dirichlet
BCt_DME.left.c[:] = c_tracer2

BC_DME = BoundaryConditionDME(BCc_oil_DME, BCc_water_DME, BCp_DME, BCs_DME, BCt_DME)

BC = [BC_water, BC_DME]#, BC_water]

# Simulating
println("Simple model")
t_sim, R_sim = DME_solver_no_gravity_no_pc_incompressible(G, K, MT, T, BC, IC, k_o, k_w, true)
PV = t_sim.*((u_inj./poros_val)./Lx)

data = [0, 0.03883495145631066,
        0.18404907975460127, 0.23300970873786409,
        0.3067484662576685, 0.4271844660194175,
        0.42944785276073616, 0.4951456310679613,
        0.8588957055214725, 0.5048543689320388,
        0.9815950920245402, 0.5242718446601942,
        1.288343558282209, 0.5339805825242718,
        1.5950920245398772, 0.5436893203883496,
        2.2085889570552153, 0.5631067961165048,
        2.883435582822086, 0.5825242718446602,
        3.619631901840491, 0.5825242718446602,
        4.294478527607362, 0.5922330097087378,
        4.969325153374233, 0.5922330097087378,
        5.644171779141105, 0.6019417475728156,
        6.380368098159509, 0.6019417475728156,
        7.0552147239263805, 0.6019417475728156,
        7.730061349693252, 0.6019417475728156,
        8.404907975460121, 0.6019417475728156,
        9.079754601226991, 0.6019417475728156,
        9.570552147239262, 0.6019417475728156,
        9.938650306748466, 0.6116504854368932,
        10.061349693251532, 0.6310679611650485,
        10.184049079754601, 0.6893203883495146,
        10.306748466257668, 0.7475728155339806,
        10.429447852760735, 0.7669902912621359,
        10.736196319018404, 0.8155339805825242,
        10.920245398773005, 0.8252427184466019,
        11.16564417177914, 0.8446601941747572,
        11.411042944785276, 0.854368932038835,
        11.717791411042944, 0.8640776699029126,
        12.085889570552146, 0.8737864077669903,
        12.760736196319018, 0.9029126213592233,
        13.43558282208589, 0.9223300970873787,
        14.478527607361961, 0.941747572815534,
        15.0920245398773, 0.9514563106796117,
        15.766871165644169, 0.9514563106796117,
        16.44171779141104, 0.9514563106796117,
        17.116564417177912, 0.9514563106796117,
        17.791411042944784, 0.9514563106796117,
        18.343558282208587, 0.9514563106796117,
        18.71165644171779, 0.9514563106796117,
        19.01840490797546, 0.9514563106796117,
        19.386503067484657, 0.9514563106796117,
        20.06134969325153, 0.9514563106796117,
        20.7361963190184, 0.9514563106796117,
        21.411042944785272, 0.9514563106796117]

N_d = length(data)
PV_exp = zeros(N_d/2)
R_exp = zeros(N_d/2)
count_PV = 1
count_R = 1
for i = 1:N_d
    if mod(i,2) == 1
        PV_exp[count_PV] = data[i]
        count_PV += 1
    elseif mod(i,2) == 0
        R_exp[count_R] = data[i]
        count_R += 1
    end
end

figure()
plot(PV_exp, R_exp.*100)
plot(PV, R_sim*100)
plot([n_pv_water1, n_pv_water1], [0.0, 100.0], color=:black, linestyle="--")
#plot([n_pv_water1 + n_pv_DME, n_pv_water1 + n_pv_DME], [0.0, 100.0], color=:black, linestyle="--")
plot([18.0, 18.0], [0.0, 100.0], color=:black, linestyle="--")
ylim([0, 100])
grid("on")
title("Recovery factor")
xlabel(L"$PVI$")
ylabel(L"$RF$"" [%]")
legend(["Experimental", "Heuristic"])
savefig("2D Results/recovery_factor.png")

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
function wall_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxp)

    ∑_Δx_sqr = zero(Float64)
    ∑_Δy_sqr = zero(Float64)
    ∑_Δx_Δy = zero(Float64)

    fill!(∑_Δx_Δf, zero(Float64))
    fill!(∑_Δy_Δf, zero(Float64))

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx::Float64 = ny
    ty::Float64 = -nx

    for itm in loc_globaldata[idx].xpos_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        ∑_Δx_sqr = ∑_Δx_sqr + dels*dels_weights
        ∑_Δy_sqr = ∑_Δy_sqr + deln*deln_weights
        ∑_Δx_Δy = ∑_Δx_Δy + dels*deln_weights

        @. qtilde_i = loc_globaldata[idx].q - 0.5*(delx * loc_globaldata[idx].dq1 + dely * loc_globaldata[idx].dq2)
        @. qtilde_k = globaldata_itm.q - 0.5*(delx * globaldata_itm.dq1 + dely * globaldata_itm.dq2)

        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), qtilde_k)
        # end

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata[idx], vl_const, phi_i)
            venkat_limiter(qtilde_k, globaldata_itm, vl_const, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq1 + dely*loc_globaldata[idx].dq2)
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq1 + dely*globaldata_itm.dq2)
        end

        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        qtilde_to_primitive(result, qtilde_i, gamma)
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), result)
        # end
        flux_quad_GxII(G_i, nx, ny, result[1], result[2], result[3], result[4])
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        qtilde_to_primitive(result, qtilde_k, gamma)
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        flux_quad_GxII(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            ∑_Δx_Δf[i] += (G_k[i] - G_i[i]) * dels_weights
            ∑_Δy_Δf[i] += (G_k[i] - G_i[i]) * deln_weights
        end
    end

    det = @. ∑_Δx_sqr*∑_Δy_sqr - ∑_Δx_Δy*∑_Δx_Δy
    one_by_det = 1.0 / det

    for iter in 1:4
        Gxp[iter] = (∑_Δx_Δf[iter]*∑_Δy_sqr - ∑_Δy_Δf[iter]*∑_Δx_Δy)*one_by_det
    end
    return nothing
end

function wall_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxn)

    ∑_Δx_sqr = zero(Float64)
    ∑_Δy_sqr = zero(Float64)
    ∑_Δx_Δy = zero(Float64)

    fill!(∑_Δx_Δf, zero(Float64))
    fill!(∑_Δy_Δf, zero(Float64))

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    for itm in loc_globaldata[idx].xneg_conn
        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        ∑_Δx_sqr = ∑_Δx_sqr + dels*dels_weights
        ∑_Δy_sqr = ∑_Δy_sqr + deln*deln_weights

        ∑_Δx_Δy = ∑_Δx_Δy + dels*deln_weights

        @. qtilde_i = loc_globaldata[idx].q - 0.5*(delx * loc_globaldata[idx].dq1 + dely * loc_globaldata[idx].dq2)
        @. qtilde_k = globaldata_itm.q - 0.5*(delx * globaldata_itm.dq1 + dely * globaldata_itm.dq2)

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata[idx], vl_const, phi_i)
            venkat_limiter(qtilde_k, globaldata_itm, vl_const, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq1 + dely*loc_globaldata[idx].dq2)
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq1 + dely*globaldata_itm.dq2)
        end

        qtilde_to_primitive(result, qtilde_i, gamma)
        flux_quad_GxI(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, gamma)
        flux_quad_GxI(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            ∑_Δx_Δf[i] += (G_k[i] - G_i[i]) * dels_weights
            ∑_Δy_Δf[i] += (G_k[i] - G_i[i]) * deln_weights
        end

    end
    det = ∑_Δx_sqr*∑_Δy_sqr - ∑_Δx_Δy*∑_Δx_Δy
    one_by_det = 1.0 / det
    for iter in 1:4
        Gxn[iter] = (∑_Δx_Δf[iter]*∑_Δy_sqr - ∑_Δy_Δf[iter]*∑_Δx_Δy)*one_by_det
    end
    return nothing
end

function wall_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gyn)

    ∑_Δx_sqr = zero(Float64)
    ∑_Δy_sqr = zero(Float64)
    ∑_Δx_Δy = zero(Float64)

    fill!(∑_Δx_Δf, zero(Float64))
    fill!(∑_Δy_Δf, zero(Float64))

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    for itm in loc_globaldata[idx].yneg_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        ∑_Δx_sqr = ∑_Δx_sqr + dels*dels_weights
        ∑_Δy_sqr = ∑_Δy_sqr + deln*deln_weights

        ∑_Δx_Δy = ∑_Δx_Δy + dels*deln_weights

        @. qtilde_i = loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq1 + dely*loc_globaldata[idx].dq2)
        @. qtilde_k = globaldata_itm.q - 0.5*(delx*globaldata_itm.dq1 + dely*globaldata_itm.dq2)

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata[idx], vl_const, phi_i)
            venkat_limiter(qtilde_k, globaldata_itm, vl_const, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq1 + dely*loc_globaldata[idx].dq2)
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq1 + dely*globaldata_itm.dq2)
        end

        qtilde_to_primitive(result, qtilde_i, gamma)
        flux_Gyn(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, gamma)
        flux_Gyn(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            ∑_Δx_Δf[i] += (G_k[i] - G_i[i]) * dels_weights
            ∑_Δy_Δf[i] += (G_k[i] - G_i[i]) * deln_weights
        end

    end

    det = ∑_Δx_sqr*∑_Δy_sqr - ∑_Δx_Δy*∑_Δx_Δy
    one_by_det = 1.0 / det
    for iter in 1:4
        Gyn[iter] = (∑_Δy_Δf[iter]*∑_Δx_sqr - ∑_Δx_Δf[iter]*∑_Δx_Δy)*one_by_det
    end
    return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver