text:
using MathOptInterface
using Xpress
using XpressPSR
using Test

const MOI = MathOptInterface

XpressPSR.initialize()

#Example of simple Thermal Generation Dispatch and the derivatives of the generation by the demand (dg/dd)
mutable struct DispachModel
    optimizer::MOI.AbstractOptimizer # Optimizer
    g::Vector{MOI.VariableIndex}     # Generation Variable Indexes
    Df::MOI.VariableIndex            # Decifit Variable Index
    c_limit_lower                    # Constraint of lower limit for variables
    c_limit_upper                    # Constraint of upper limite for variables
    c_demand                         # Constraint of the generation by the demand
end

#Simple Model of Thermal Generation Dispatch by Constraints with MOI.SingleVariable
function GenerateModel_SingleVariable()
    # Parameters
    d = 45.0                   # Demand
    I = 3                      # Number of generators
    g_sup = [10.0, 20.0, 30.0] # Upper limit of generation for each generator
    c_g = [1.0, 3.0, 5.0]      # Cost of generation for each generator
    c_Df = 10.0                # Cost of deficit

    # Creation of the Optimizer
    optimizer = Xpress.Optimizer(PRESOLVE=0, logfile = "outputXpress_SV.log")
    # Variables
    g = MOI.add_variables(optimizer, I) # Generation for each generator
    Df = MOI.add_variable(optimizer)    # Deficit

    # Constraints
    c_limit_inf = Vector{Any}(undef, I + 1) # Lower limit constraints
    c_limit_sup = Vector{Any}(undef, I)     # Upper limit constraints
    for i in 1:I
        c_limit_inf[i] = MOI.add_constraint(optimizer, MOI.SingleVariable(g[i]), MOI.GreaterThan(0.0))
        c_limit_sup[i] = MOI.add_constraint(optimizer, MOI.SingleVariable(g[i]), MOI.LessThan(g_sup[i]))
    end
    c_limit_inf[I+1] = MOI.add_constraint(optimizer, MOI.SingleVariable(Df), MOI.GreaterThan(0.0))

    c_demand = MOI.add_constraint(optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(I+1),[g;Df]), 0.0),
        MOI.EqualTo(d)
    ) # Constraint of the Demand

    # Objectives
    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([c_g; c_Df], [g; Df]), 0.0)     # Total cost function
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)  # Pass the objective function to the optimizer
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)                                             # Set the objective sense to MIN

    # Solve the optimizer
    MOI.optimize!(optimizer)

    MOI.write_to_file(optimizer, "LP_SV.lp")

    # Return the solved optimizer with interest variables and constraints
    return DispachModel(optimizer, g, Df, c_limit_inf, c_limit_sup, c_demand)
end

#Simple Model of Thermal Generation Dispatch by Constraints with MOI.ScalarAffineFunction
function GenerateModel_ScalarAffineFunction()
    # Parameters
    d = 45.0                   # Demand
    I = 3                      # Number of generators
    g_sup = [10.0, 20.0, 30.0] # Upper limit of generation for each generator
    c_g = [1.0, 3.0, 5.0]      # Cost of generation for each generator
    c_Df = 10.0                # Cost of deficit

    # Creation of the Optimizer
    optimizer = Xpress.Optimizer(PRESOLVE=0, logfile = "outputXpress_SAF.log")
    # Variables
    g = MOI.add_variables(optimizer, I) # Generation for each generator
    Df = MOI.add_variable(optimizer)    # Deficit

    # Constraints
    c_limit_inf = Vector{Any}(undef, I + 1) # Lower limit constraints
    c_limit_sup = Vector{Any}(undef, I)     # Upper limit constraints
    vectAux = zeros(I + 1)                  # Auxiliar Vector
    for i in 1:I
        vectAux[i] = 1.0
        c_limit_inf[i] = MOI.add_constraint(optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vectAux,[g;Df]), 0.0),
            MOI.GreaterThan(0.0)
        )
        c_limit_sup[i] = MOI.add_constraint(optimizer,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vectAux,[g;Df]), 0.0),
            MOI.LessThan(g_sup[i])
        )
        vectAux[i] = 0.0
    end
    vectAux[I+1] = 1.0
    c_limit_inf[I+1] = MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(vectAux,[g;Df]),0.0), MOI.GreaterThan(0.0))
    vectAux[I+1] = 0.0

    c_demand = MOI.add_constraint(optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(I+1),[g;Df]), 0.0),
        MOI.EqualTo(d)
    ) # Constraint of the Demand

    # Objectives
    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([c_g; c_Df], [g; Df]), 0.0)     # Total cost function
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)  # Pass the objective function to the optimizer
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)                                             # Set the objective sense to MIN

    # Solve the optimizer
    MOI.optimize!(optimizer)

    MOI.write_to_file(optimizer, "LP_SAF.lp")

    # Return the solved optimizer with interest variables and constraints
    return DispachModel(optimizer, g, Df, c_limit_inf, c_limit_sup, c_demand)
end

#It uses a solved (optimized) Dispach Model created by the function GenerateModel
function Forward(model::DispachModel, ϵ::Float64 = 1.0)
    #Initialization of parameters and references to simplify the notation
    optimizer = model.optimizer
    vectRef = [model.g; model.Df]
    I = length(model.g)

    #Get the primal solution of the model
    vect =  MOI.get.(optimizer, MOI.VariablePrimal(), vectRef)

    #Pass the perturbation to the Xpress Framework and set the context to Forward
    MOI.set(optimizer, Xpress.ForwardSensitivityInConstraint(), model.c_demand, ϵ)

    #Get the derivative of the model
    dvect = MOI.get.(optimizer, Xpress.ForwardSensitivityOutVariable(), vectRef)

    #Return the values as a vector
    return [vect;dvect]
end

function Backward(model::DispachModel, ϵ::Float64 = 1.0)
    #Initialization of parameters and references to simplify the notation
    optimizer = model.optimizer
    vectRef = [model.g; model.Df]
    I = length(model.g)

    #Get the primal solution of the model
    vect =  MOI.get.(optimizer, MOI.VariablePrimal(), vectRef)

    #Set variables needed for the Xpress Backward Framework
    dvect = zeros(I + 1)

    #Loop for each primal variable
    for i in 1:I + 1
        #Set the perturbation in the Primal Variables and set the context to Backward
        MOI.set(optimizer, Xpress.BackwardSensitivityInVariable(), vectRef[i], ϵ)

        #Get the value of the derivative of the model
        dvect[i] = MOI.get(optimizer, Xpress.BackwardSensitivityOutConstraint(), model.c_demand)

        #Reset perturbation
        MOI.set(optimizer, Xpress.BackwardSensitivityInVariable(), vectRef[i], 0.0)
    end

    #Return the values as a vector
    return [vect;dvect]
end

#Initialization of Parameters
@testset verbose = true "Results" begin
    @testset verbose = true "Single_Variable" begin
        #Generate model by SingleVariable
        model1 = GenerateModel_SingleVariable()
        #Get the results of models with the DiffOpt Forward context
        @testset "Forward" begin
            resultForward = Forward(model1)
            @test resultForward[1:4] == [10.0;20.0;15.0;0.0]
            @test resultForward[5:8] == [0.0;0.0;1.0;0.0]
        end
        #Get the results of models with the DiffOpt Backward context
        @testset "Backward" begin
           resultBackward = Backward(model1)
           @test resultBackward[1:4] == [10.0;20.0;15.0;0.0]
           @test resultBackward[5:8] == [0.0;0.0;1.0;0.0]
        end
    end
    @testset verbose = true "Scalar_Affine_Function" begin
        #Generate model by ScalarAffineFunction
        model2 = GenerateModel_ScalarAffineFunction()
        #Get the results of models with the DiffOpt Forward context
        @testset verbose = true "Forward" begin
            resultForward = Forward(model2)
            @test resultForward[1:4] == [10.0;20.0;15.0;0.0]
            @test resultForward[5:8] == [0.0;0.0;1.0;0.0]
        end
       #Get the results of models with the DiffOpt Backward context
       @testset verbose = true "Backward" begin
           resultBackward = Backward(model2)
           @test resultBackward[1:4] == [10.0;20.0;15.0;0.0]
           @test resultBackward[5:8] == [0.0;0.0;1.0;0.0]
       end
    end
end

;#END


text:
using Logging

@testset "Travis CI deploy configuration" begin; with_logger(NullLogger()) do
    # Regular tag build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "v1.2.3",
            "TRAVIS_EVENT_TYPE" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
        @test Documenter.authentication_method(cfg) === Documenter.SSH
    end
    # Broken tag build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "not-a-version",
            "TRAVIS_EVENT_TYPE" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) === nothing
    end
    # Regular/broken devbranch build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => nothing,
            "TRAVIS_EVENT_TYPE" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Regular pull request build
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "42",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "something",
            "TRAVIS_TAG" => nothing,
            "TRAVIS_EVENT_TYPE" => nothing,
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Missing/broken environment variables
    withenv("TRAVIS_CI" => "true",
            "TRAVIS_PULL_REQUEST" => "false",
            "TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "TRAVIS_BRANCH" => "master",
            "TRAVIS_TAG" => "v1.2.3",
            "TRAVIS_EVENT_TYPE" => nothing,
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.Travis()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=false) === nothing
    end
end end

@testset "GitHub Actions deploy configuration" begin; with_logger(NullLogger()) do
    # Regular tag build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular tag build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Broken tag build
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/not-a-version",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) === nothing
    end
    # Regular devbranch build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/heads/master",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular devbranch build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/heads/master",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "hello-world"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=true) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Regular pull request build with GITHUB_TOKEN
    withenv("GITHUB_EVENT_NAME" => "pull_request",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/pull/42/merge",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.HTTPS
        @test Documenter.authenticated_repo_url(cfg) === "https://github-actions:SGVsbG8sIHdvcmxkLg==@github.com/JuliaDocs/Documenter.jl.git"
    end
    # Regular pull request build with SSH deploy key (SSH key prioritized)
    withenv("GITHUB_EVENT_NAME" => "pull_request",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/pull/42/merge",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => "SGVsbG8sIHdvcmxkLg==",
            "DOCUMENTER_KEY" => "SGVsbG8sIHdvcmxkLg==",
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) == "previews/PR42"
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="not-master", devurl="hello-world", push_preview=false) === nothing
        @test Documenter.authentication_method(cfg) === Documenter.SSH
        @test Documenter.documenter_key(cfg) === "SGVsbG8sIHdvcmxkLg=="
    end
    # Missing environment variables
    withenv("GITHUB_EVENT_NAME" => "push",
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
            "GITHUB_REF" => "refs/tags/v1.2.3",
            "GITHUB_ACTOR" => "github-actions",
            "GITHUB_TOKEN" => nothing,
            "DOCUMENTER_KEY" => nothing,
        ) do
        cfg = Documenter.GitHubActions()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="hello-world", push_preview=true) === nothing
    end
end end

struct CustomConfig <: Documenter.DeployConfig end
Documenter.deploy_folder(::CustomConfig; kwargs...) = "v1.2.3"
struct BrokenConfig <: Documenter.DeployConfig end

@testset "Custom configuration" begin; with_logger(NullLogger()) do
        cfg = CustomConfig()
        @test Documenter.deploy_folder(cfg; repo="github.com/JuliaDocs/Documenter.jl.git",
                  devbranch="master", devurl="dev", push_preview=true) == "v1.2.3"
        cfg = BrokenConfig()
        @test (@test_logs (:warn, r"Documenter\.deploy_folder\(::BrokenConfig; kwargs\.\.\.\) not implemented") Documenter.deploy_folder(cfg)) === nothing
        @test (@test_logs (:warn, r"Documenter could not auto-detect") Documenter.deploy_folder(nothing)) === nothing
end end

@testset "Autodetection of deploy system" begin
    withenv("TRAVIS_REPO_SLUG" => "JuliaDocs/Documenter.jl",
            "GITHUB_REPOSITORY" => nothing,
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg isa Documenter.Travis
    end
    withenv("TRAVIS_REPO_SLUG" => nothing,
            "GITHUB_REPOSITORY" => "JuliaDocs/Documenter.jl",
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg isa Documenter.GitHubActions
    end
    withenv("TRAVIS_REPO_SLUG" => nothing,
            "GITHUB_REPOSITORY" => nothing,
        ) do
        cfg = Documenter.auto_detect_deploy_system()
        @test cfg === nothing
    end
end

ext: jl


text:
#=
Popis
  Funkcia vypocita hodnotu radiality centrality pre kazdy vrchol.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozene ohodnotenie hrany :: String
  normalize : parameter pre aplikovanie normalizacie :: Bool

Navratova hodnota
  radiality : slovnik s mapovanim (kluc=vrchol_grafu => hodnota=radiality) :: Dict{Int64, Float64}

Poznamky
=#

function centrality_radiality(g::GenericGraph; dist_key::AbstractString="weight", normalize::Bool=true)

  dia = graph_diameter(g)
  epi = AttributeEdgePropertyInspector{Float64}(dist_key)
  radiality = [v.index => 0.0 for v=vertices(g)]
  total_dists = Float64[]
  temp_rad = Float64[]
  temp_sum = Float64[]

  if normalize

    for v in vertices(g)

      temp_dists = dijkstra_shortest_paths(g, epi, [v]).dists
      total_dists = filter(x -> x != Inf, temp_dists)

        for s in total_dists
          push!(temp_sum, dia+1-s)
        end

      temp_rad = ((sum(temp_sum) -dia -1)/(num_vertices(g) - 1))/dia
      empty!(temp_sum)
      radiality[v.index] = temp_rad

    end

  elseif !(normalize)

    for v in vertices(g)

      temp_dists = dijkstra_shortest_paths(g, epi, [v]).dists
      total_dists = filter(x -> x != Inf, temp_dists)

        for s in total_dists
          push!(temp_sum, dia+1-s)
        end

      temp_rad = (sum(temp_sum) -dia -1)/(num_vertices(g) - 1)
      empty!(temp_sum)
      radiality[v.index] = temp_rad

    end

  end

  return radiality

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
"""
showSampling
 List detailed information on sampling in time series (e.g monthly occurrences, gaps etc)

Syntax

 showSampling(in::tsParm,days4gap::Float64)

Description

 List detailed information on sampling in time series (e.g monthly occurrences, gaps etc)
 The detailed information is useful for including in a report and find gaps.  It also checks for
  duplicate dates in the series.

 If value vector is provided then any NaNs will be treated as gaps.

 INPUTS

  in::tsParm       	tsParm structure containing at least sDate and Data
                          sDate will be used to define dates
			  A column (iColumn) will be checked for NaNs which are considered missing data
  iColumn[1]            Integer defining column in Data which will be used to assess
                          missing data
  days4gap              Define minimum gaps size
                          (This is not used for calculating % data return)

 OUTPUTS

   parcels,goodDurationDays,iParcels

   parcels              Structure contain information on Parcels of data separated by the gap size
                             as determined by MetoceanTools.findGaps
   durationGoodDays     Total effective duration [days] the data represent
                          determined by the total duration contained within contiguous blocks of
                          data with no gaps greater than days4gap
   iParcels             Array of length in.sDate containing sequential parcel number or 0 if a gap

Revisions

 Revision 1 26 Dec 2016 Jason McConochie

"""
function showSampling(in::tsParm,iColumn::Int64,days4gap::Float64)

  o=STDOUT
  @printf(o,"==========================================================================================\n")
  @printf(o,"Sampling Analysis Output - showSampling\n")
  @printf(o,"  Generated: %s\n",Dates.format(now(),"dd/uuu/yyyy HH:MM:SS"))
  @printf(o,"==========================================================================================\n")

  t=in.sDate
  z=in.Data[:,iColumn]

  #..check for duplicate dates
      #### Testing tsParm for duplicates, repeated data, should be a function
  usDate=unique(in.sDate)
  dd=length(in.sDate)!=length(usDate)
  if dd>0
     @printf(o,".. there are %d duplicate dates\n",dd)
     #  ADD CODE TO LIST THE UNIQUE DATES
  end

  #..get good data
  isGood=!isnan(z)
  tGood=t[isGood]
  zGood=z[isGood]
  if isempty(tGood)
       @printf(o,"..No dates with good data\n")
       @printf(o,"==========================================================================================\n")
       return
  end

  ms2days=1/(1000*60*60*24)
  ms2hours=1/(1000*60*60)
  ms2minutes=1/(1000*60)
  ms2seconds=1/(1000)

  @printf(o,"Time period covered by good data (i.e. not NaN)\n")
  @printf(o,"  start at %20s\n",Dates.format(tGood[1],"dd/uuu/yyyy HH:MM:SS"))
  @printf(o,"  ends  at %20s\n",Dates.format(tGood[end],"dd/uuu/yyyy HH:MM:SS"))
  tdur=(Dates.value(tGood[end])-Dates.value(tGood[1]))*ms2days
  @printf(o,"Total duration is %g days\n",tdur)

  ##### call s=tsSimpleStats(tsParm) and print more information
  @printf(o,"Simple statistics\n")
  @printf(o,"    Min value=%g\n",minimum(zGood))
  @printf(o,"    Mean value=%g\n",mean(zGood))
  @printf(o,"    Max value=%g\n",maximum(zGood))

  #..calculate the number of values in each year and month
  yy=Dates.year(tGood)
  mm=Dates.month(tGood)
  HH=Dates.hour(tGood)
  uyy=unique(yy)

  @printf(o,"Number of samples in each year\n")
  out=""
  for i=1:length(uyy)
     out=*(out,@sprintf("%4d: %-8d ",uyy[i],sum(uyy[i].==yy)))
     if i/5==round(i/5)
         @printf(o,"   %s\n",out)
         out=""
     end
  end
  if length(uyy)/5!=round(length(uyy)/5)
     @printf(o,"   %s\n",out)
  end

  @printf(o,"Number of samples in each month\n")
  out=""
  for i=1:12
     out=*(out,@sprintf("%3d: %-8d ",i,sum(i.==mm)))
     if i/6==round(i/6)
         @printf(o,"   %s\n",out)
         out=""
     end
  end
  @printf(o,"   %s\n",out)

  @printf(o,"Number of samples in each hour of the day\n")
  out=""
  for i=0:23
    out=*(out,@sprintf("%3d: %-8d ",i,sum(i.==HH)))
    if (i+1)/6==round((i+1)/6)
      @printf(o,"   %s\n",out)
      out=""
    end
  end
  @printf(o,"   %s\n",out)

  if length(zGood)==1
    @printf(o,"==========================================================================================\n")
    return
  end

  @printf(o,"Most common sampling intervals of good data\n")
  @printf(o,"%10s %10s %10s %10s %10s\n","days","hours","minutes","seconds","occurrence")
  #.. calculate in seconds to nearest millisecond
  (uCounts,uValues)=MetoceanTools.countUniqueValues(diff(Dates.value(tGood)))
  isCounts=sortperm(uCounts,rev=true)
  sCounts=uCounts[isCounts]
  sValues=uValues[isCounts]
  for i=1:minimum([length(sValues),5])
    @printf(o,"%10g %10g %10g %10g %10d\n",sValues[i]*ms2days,sValues[i]*ms2hours,sValues[i]*ms2minutes,sValues[i]*ms2seconds,sCounts[i])
  end

  #.. percentage of values at the principal sampling interval for each month of each year
  @printf(o,"Percentage occurrence of most common sampling interval: %g hour(s)\n",uValues[1]*ms2hours)
  @printf(o,"         JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC\n")
  for iyy=1:length(uyy)
    oStr=@sprintf("    %4d ",uyy[iyy])
    for imth=1:12
      mask=(yy.==uyy[iyy]) & (mm.==imth)  #this month and year
      nvalues=sum(mask)
      if nvalues>1
          dd=diff(round(Dates.value(tGood[mask])))
          cnt=sum(dd .== uValues[1])
          pct=round(cnt/nvalues*100)
          oStr=*(oStr,@sprintf("%3d ",pct))
      else
          oStr=*(oStr,@sprintf("%3d ",0))
      end
    end
    @printf(o,"%s \n",oStr)
  end

  #.. percentage of values at the second most common sampling interval for each month of each year
  if length(uValues)>1
    @printf(o,"Percentage occurrence of second most common sampling interval: %g hour(s)\n",uValues[2]*ms2hours)
    @printf(o,"         JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC\n")
    for iyy=1:length(uyy)
      oStr=@sprintf("    %4d ",uyy[iyy])
      for im=1:12
        mask=(yy.==uyy[iyy]) & (mm.==im)  #this month and year
        nvalues=sum(mask)
        if nvalues>1
            dd=diff(round(Dates.value(tGood[mask])))
            cnt=sum(dd .== uValues[2])
            pct=round(cnt/nvalues*100)
            oStr=*(oStr,@sprintf("%3d ",pct))
        else
            oStr=*(oStr,@sprintf("%3d ",0))
        end
      end
      @printf(o,"%s \n",oStr)
    end
  end

  #.. now if requested check for gaps in the data
  parcels,goodDurationDays,iParcels=MetoceanTools.findGaps(in,iColumn,days4gap)
  #..loop over all parcels
  if any(parcels.ID.==0) # Found gaps
     dur=0.
     ig=0
     #@printf(o,"Found %d gaps greater than %g day(s)\n",length(p.sDate)-1,in.days4gap)
     for ip=1:length(parcels.ID)
        if parcels.ID[ip]==0
            ig+=1
            @printf(o,"   Gap %d between %s and %s, duration %8g days\n",ig,
            Dates.format(parcels.sDate[ip,1],"dd/uuu/yyyy HH:MM:SS"),
            Dates.format(parcels.sDate[ip,2],"dd/uuu/yyyy HH:MM:SS"),parcels.Duration[ip])
  	  dur+=parcels.Duration[ip]
         end
     end
     pct=100*sum(dur)./tdur
     @printf(o,"There is %g %% gaps greater than %g day(s) in the data\n",pct,days4gap)
  end

  # Now find out how many samples are missing assuming that we should
  # have all values with the most common sampling value
  @printf(o,"==========================================================================================\n")
  selt=Dates.value(tGood[1]):uValues[1]:Dates.value(tGood[end])
  @printf(o,"If the record had all samples taken at the most common sampling interval\n")
  @printf(o,"   then there should be %d samples.\n",length(selt))
  i=intersect(Dates.value(tGood),selt)
  @printf(o,"   There are %d of these samples found.\n",length(i))
  @printf(o,"   Thus the data represents %g %%.\n",100*length(i)/length(selt))

  @printf(o,"==========================================================================================\n")

  return parcels,goodDurationDays,iParcels

end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
## Versión in-place intermedia para evaluación:
#  capitaliza las matrices de variaciones in-place para ahorrar memoria

# export InflationTotalEvalCPI

Base.@kwdef struct InflationTotalEvalCPI <: InflationFunction
    name::String = "Variación interanual IPC"
end

# Función para bases cuyo índice base es un escalar
function (inflfn::InflationTotalEvalCPI)(base::VarCPIBase{T, T}) where {T <: AbstractFloat}
    base_ipc = convert(IndexCPIBase, base)
    ipc = base_ipc.ipc * base.w / base.baseindex
    varinterm!(ipc, ipc, 100)
    ipc
end

# # Variación interanual
# function (inflfn::InflationTotalEvalCPI)(cs::CountryStructure)
#     vm = mapfoldl(inflfn, vcat, cs.base)
#     capitalize!(vm, vm, 100)
#     varinteran!(vm)
#     vm[12:end]
# end

## Versión in-place extrema para evaluación:
#  capitaliza las matrices de variaciones in-place para ahorrar memoria y
#  guarda en vector los resultados

# Base.@kwdef struct TotalExtremeCPI <: InflationFunction
#     name::String = "Variación interanual IPC"
# end

# # # Función para bases cuyo índice base es un escalar
# # function (inflfn::TotalExtremeCPI)(base::VarCPIBase{T, T}) where {T <: AbstractFloat}
# #     base_ipc = convert(IndexCPIBase, base)
# #     ipc = base_ipc.ipc * base.w / base.baseindex
# #     varinterm!(ipc, ipc, 100)
# #     ipc
# # end

# import LinearAlgebra: mul!

# # Variación interanual
# function (inflfn::TotalExtremeCPI)(tray_infl::AbstractVector, cs::CountryStructure)
#     B = length(cs.base) # número de bases

#     # Para cada base, obtener IPC y variación intermensual
#     ib = 0
#     for b in 1:B
#         l = ib + size(cs[b].v, 1) # número de períodos
#         interm = @view tray_infl[ib+1:l] # vector para guardar inflación intermensual

#         base_ipc = convert(IndexCPIBase, cs.base[b])
#         mul!(interm, base_ipc.ipc, base_ipc.w)
#         varinterm!(interm, interm, 100*base_ipc.baseindex)
#         ib = l
#     end

#     # Capitalizar y obtener var interanual
#     capitalize!(tray_infl, tray_infl, 100)
#     varinteran!(tray_infl)
#     nothing
# end


text:
module RoundFast

import Base:(+),(-),(*),(/),sqrt
export square

include("../../module/AdjacentFloat/AdjacentFloat.jl")
using .AdjacentFloat

include("eftArith.jl")
include("eftRound.jl")
include("rounding.jl")

end # module


text:
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @warnpcfail precompile(HyperRectangle{2,Float32}, (Int, Int, Int, Int))
    @warnpcfail precompile(==, (HyperRectangle{2,Float32}, HyperRectangle{2,Float32}))
    @warnpcfail precompile(normal_mesh, (Tesselation{3,Float32,Cylinder{3,Float32},1},))
    @warnpcfail precompile(normal_mesh, (Tesselation{3,Float32,HyperSphere{3,Float32},1},))
    @warnpcfail precompile(normal_mesh, (HyperSphere{3,Float32},))
    @warnpcfail precompile(Tuple{typeof(*),SMatrix{4, 4, Float32, 16},HyperRectangle{3, Float32}})   # time: 0.11091917
    @warnpcfail precompile(Tuple{typeof(coordinates),HyperRectangle{2, Float32},Tuple{Int64, Int64}})   # time: 0.08693867
    @warnpcfail precompile(union, (HyperRectangle{3, Float32}, HyperRectangle{3, Float32}))
    @warnpcfail precompile(Tuple{typeof(decompose),Type{Point{2, Float32}},HyperRectangle{2, Float32}})   # time: 0.026609203
    @warnpcfail precompile(Tuple{Type{HyperRectangle{3, Float32}},HyperRectangle{2, Float32}})   # time: 0.023717888
    @warnpcfail precompile(Tuple{typeof(+),HyperRectangle{3, Float32},Point{3, Float32}})   # time: 0.006633118
    @warnpcfail precompile(Tuple{Type{Rect2{T} where T},Float32,Float32,Float32,Float32})   # time: 0.001636267
    @warnpcfail precompile(Tuple{typeof(*),HyperRectangle{2, Float32},Float32})   # time: 0.001057589

    if Base.VERSION >= v"1.6.0-DEV.1083"
        @warnpcfail precompile(triangle_mesh, (Polygon{2, Float32, Point2f, LineString{2, Float32, Point2f,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}, false}},
                                        Vector{LineString{2, Float32, Point2f, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}, false}}}},))
    else
        @warnpcfail precompile(triangle_mesh, (Polygon{2, Float32, Point2f, LineString{2, Float32, Point2f,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}}},
                                        Vector{LineString{2, Float32, Point2f, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}}}}},))
    end

    @warnpcfail precompile(split_intersections, (Vector{Point2f},))
end


text:
using Random
using LightGraphs
import Base: ==
using JLD

struct TestStruct{T}
    a::Dict{Symbol,Vector}
    b::Set
    c::Array{T,1}
end

function TestStruct()
    TestStruct(Dict{Symbol,Vector}(), Set(), Int[])
end

==(x::TestStruct, y::TestStruct) = (x.a==y.a) && (x.b==y.b) && (x.c==y.c)

@testset "Persistence" begin
    g = StorageGraph()
    dep = ((a=1,),(b=2,),(c=3,))
    add_nodes!(g, foldr(=>, (dep..., (d=rand(2),))))
    d = g[:d, dep...]
    add_nodes!(g, foldr(=>, (dep..., (d=d,), (x=randstring.(1:2),), (y=TestStruct(),))))

    if VERSION ≥ v"1.1"
        savegraph("test.jls", g, SGNativeFormat())
        g_serialize = loadgraph("test.jls", SGNativeFormat())
        @test g == g_serialize
        @test eltype(g_serialize) == eltype(g)
        @test eltype(g_serialize.data) == eltype(g.data)
        @test eltype(g_serialize.index) == eltype(g.index)
        @test eltype(g_serialize.paths) == eltype(g.paths)
        rm("test.jls")
    end

    savegraph("test.jld", g, "g", SGJLDFormat())
    g_jld = loadgraph("test.jld", "g", SGJLDFormat())
    @test g == g_jld
    @test g.data == g_jld.data
    @test eltype(g_jld) == eltype(g)
    @test eltype(g_jld.data) == eltype(g.data)
    @test eltype(g_jld.index) == eltype(g.index)
    @test eltype(g_jld.paths) == eltype(g.paths)
    rm("test.jld")

    savegraph("test.bson", g, :g, SGBSONFormat())
    g_bson = loadgraph("test.bson", :g, SGBSONFormat())
    @test g == g_bson
    @test eltype(g_bson) == eltype(g)
    @test eltype(g_bson.data) == eltype(g.data)
    @test eltype(g_bson.index) == eltype(g.index)
    @test eltype(g_bson.paths) == eltype(g.paths)

    rm("test.bson")
end


text:
import Term: Tree

tree_dict = Dict(
    "nestedasdasdsadasdasdsadsadasdasdasdsadasasdasdassfsdfdsfdsfdsfsdfsdfdsfsdfd" =>
        Dict("n1" => 1, "n2" => 2),
)

tree_dict_1 =
    Dict("nested" => Dict("n1" => 1, "n2" => 2), "nested2" => Dict("n1" => "a", "n2" => 2))

tree_dict_2 = Dict(
    "nested" => Dict("n1" => 1, "n2" => 2),
    "leaf2" => 2,
    "leaf" => 2,
    "leafme" => "v",
    "canopy" => "test",
)

tree_dict_3 = Dict(
    "nested" => Dict(
        "deeper" => Dict("aleaf" => "unbeliefable", "leaflet" => "level 3"),
        "n2" => Int,
        "n3" => 1 + 2,
    ),
    "nested2" => Dict("n1" => "a", "n2" => 2),
)

tree_dict_4 = Dict(
    "nested" => Dict(
        "deeper" => Dict(
            "aleaf" => "unbeliefable",
            "leaflet" => "level 3",
            "sodeep" => Dict("a" => 4),
        ),
        "n2" => Int,
        "n3" => 1 + 2,
        "adict" => Dict("x" => 2),
    ),
    "nested2" => Dict("n1" => "a", "n2" => 2),
)

@testset "\e[34mTree" begin
    # creation
    testtree(Tree(tree_dict), 50, 6)

    testtree(Tree(tree_dict_1), 15, 9)

    testtree(Tree(tree_dict_2), 18, 11)

    testtree(Tree(tree_dict_3), 33, 13)

    testtree(Tree(tree_dict_4), 33, 19)

    # styling
    for guides_type in (:standardtree, :boldtree, :asciitree)
        @test_nothrow Tree(
            tree_dict;
            title = string(guides_type),
            guides_type = guides_type,
        )
    end

    testtree(
        Tree(
            tree_dict;
            title = "my colors",
            title_style = "bold red",
            node_style = "blue underline",
            leaf_style = "green",
            guides_style = "red dim",
        ),
        52,
        6,
    )

    # test with no errors
    @test_nothrow Tree(Float64)
    @test_nothrow Tree(AbstractFloat)
end


text:
using Dates
using ALGCM1D
#include("./src/ALGCM1D.jl")

stime=DateTime(0,7,1)
ALGCM1D.run(stime)
#include("constants.jl")
#include("holtslag_psis.jl")
#include("holtslag_psim.jl")
#include("VelocityScale.jl")

#VelocityScale.(1:100,100.,-0.01,0.01,300,0.06)


text:
# Julia wrapper for header: /Users/dls/jl/openssl-1.0.1e/include/openssl/md5.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 MD5_Init (Ptr{MD5_CTX},) libcrypto
@c Int32 MD5_Update (Ptr{MD5_CTX}, Ptr{None}, Csize_t) libcrypto
@c Int32 MD5_Final (Ptr{Uint8}, Ptr{MD5_CTX}) libcrypto
@c Ptr{Uint8} MD5 (Ptr{Uint8}, Csize_t, Ptr{Uint8}) libcrypto
@c None MD5_Transform (Ptr{MD5_CTX}, Ptr{Uint8}) libcrypto


text:
using ConvertBNG
using Test

const PREC_EN = 0.002       #  2 mm for easting, northing absolute accuracy
const PREC_LL = 0.0000001   #  8 dp for lat, lon absolute accuracy
const PREC_FL = 0.000000001 # 10 dp for floating point comparison

const AIRY_1830_SEMI_MAJOR = 6377563.396
const AIRY_1830_SEMI_MINOR = 6356256.909

@testset "All tests" begin

    @testset "caister water tower" begin
        # From Transformations and OSGM15 User Guide - Annexe A
        # Caister Water Tower example
        # ETRS89 longitude, latitude: (1.716073973, 52.658007833)
        # ETRS89 eastings, northings: (651307.003, 313255.686)
        # Shifts eastings, northings: (102.801, -78.236)
        # OSGB36 eastings, northings: (651409.804, 313177.450)
        @test convert_etrs89_to_osgb36(651307.003, 313255.686) ≈ [651409.804 313177.450] atol=PREC_EN
        @test convert_bng(1.716073973, 52.658007833) ≈ [651409.804 313177.450] atol=PREC_EN
        @test convert_etrs89(1.716073973, 52.658007833) ≈ [651307.003 313255.686] atol=PREC_EN
        @test convert_etrs89_to_ll(651307.003, 313255.686) ≈ [1.716073973 52.658007833] atol=PREC_LL
        @test convert_lonlat(651409.804, 313177.450) ≈ [1.716073973 52.658007833] atol=PREC_LL
        # Check OSTN15 shift retrival
        @test collect(get_ostn_ref(651, 313)) ≈ [102.787, -78.242, 44.236] atol=PREC_FL
        @test collect(ostn15_shifts(651307.003, 313255.686)) ≈ [102.800790098, -78.235984076, 44.228405219] atol=PREC_FL

        # From Transformations and OSGM15 User Guide - Annexe C
        # Caister Water Tower example
        # E,N to lon,lat using Airy 1830 ellipsoid
        @test convert_to_ll(651409.903, 313177.270, AIRY_1830_SEMI_MAJOR, AIRY_1830_SEMI_MINOR) ≈ [1.717921583 52.65757031] atol=PREC_LL
    end

    @testset "min/max bounds tests" begin
        # failed ostn15 array retrieval
        @test get_ostn_ref(999, 999)[1] === NaN
        @test get_ostn_ref(999, 999)[2] === NaN
        @test get_ostn_ref(999, 999)[3] === NaN

        # Test min/max easting, northing, longitude, latitude
        # bad max easting
        @test convert_lonlat(700000.1, 1240000.0)[1] === NaN
        @test convert_lonlat(700000.1, 1240000.0)[2] === NaN
        # bad max northing
        @test convert_lonlat(600000.0, 1250000.1)[1] === NaN
        @test convert_lonlat(600000.0, 1250000.1)[2] === NaN
        # bad lon
        @test convert_bng(181., 51.44533267)[1] === NaN
        @test convert_bng(181., 51.44533267)[2] === NaN
        # bad lat
        @test convert_bng(-0.32824866, -90.01)[1] === NaN
        @test convert_bng(-0.32824866, -90.01)[2] === NaN

        # Test extent bounds of check function
        # below min lon
        max_lon = 1.768960
        min_lon = -6.379880
        @test check(-6.379881, (min_lon, max_lon)) === NaN
        # below min lat
        max_lat = 55.811741
        min_lat = 49.871159
        @test check(49.871158, (min_lat, max_lat)) === NaN
        # above max lon
        max_lon = 1.768960
        min_lon = -6.379880
        @test check(1.768961, (min_lon, max_lon)) === NaN
        # above max lat
        max_lat = 55.811741
        min_lat = 49.871159
        @test check(55.811742, (min_lat, max_lat)) === NaN
    end

    @testset "OS latlon to OSGB36 testset" begin
        # TP01 to TP40
        @test convert_bng(-6.299777520, 49.92226394) ≈ [  91492.146   11318.804] atol=PREC_EN
        @test convert_bng(-5.203046100, 49.96006138) ≈ [ 170370.718   11572.405] atol=PREC_EN
        @test convert_bng(-4.108645636, 50.43885826) ≈ [ 250359.811   62016.569] atol=PREC_EN
        @test convert_bng(-1.297822772, 50.57563665) ≈ [ 449816.371   75335.861] atol=PREC_EN
        @test convert_bng(-1.450514337, 50.93127938) ≈ [ 438710.920  114792.250] atol=PREC_EN
        @test convert_bng(-3.551283492, 51.40078220) ≈ [ 292184.870  168003.465] atol=PREC_EN
        @test convert_bng( 1.444547304, 51.37447026) ≈ [ 639821.835  169565.858] atol=PREC_EN
        @test convert_bng(-2.544076183, 51.42754743) ≈ [ 362269.991  169978.690] atol=PREC_EN
        @test convert_bng(-0.119925572, 51.48936565) ≈ [ 530624.974  178388.464] atol=PREC_EN
        @test convert_bng(-4.308524770, 51.85890896) ≈ [ 241124.584  220332.641] atol=PREC_EN
        @test convert_bng( 0.897243270, 51.89436637) ≈ [ 599445.590  225722.826] atol=PREC_EN
        @test convert_bng(-2.154586144, 52.25529382) ≈ [ 389544.190  261912.153] atol=PREC_EN
        @test convert_bng(-0.912489570, 52.25160951) ≈ [ 474335.969  262047.755] atol=PREC_EN
        @test convert_bng( 0.401535471, 52.75136687) ≈ [ 562180.547  319784.995] atol=PREC_EN
        @test convert_bng(-1.197476559, 52.96219109) ≈ [ 454002.834  340834.943] atol=PREC_EN
        @test convert_bng(-2.640493208, 53.34480280) ≈ [ 357455.843  383290.436] atol=PREC_EN
        @test convert_bng(-4.289180698, 53.41628516) ≈ [ 247958.971  393492.909] atol=PREC_EN
        @test convert_bng(-4.289177929, 53.41630925) ≈ [ 247959.241  393495.583] atol=PREC_EN
        @test convert_bng(-3.040454907, 53.77911026) ≈ [ 331534.564  431920.794] atol=PREC_EN
        @test convert_bng(-1.663791682, 53.80021520) ≈ [ 422242.186  433818.701] atol=PREC_EN
        @test convert_bng(-4.634521682, 54.08666318) ≈ [ 227778.330  468847.388] atol=PREC_EN
        @test convert_bng(-0.077731332, 54.11685144) ≈ [ 525745.670  470703.214] atol=PREC_EN
        @test convert_bng(-4.388491181, 54.32919541) ≈ [ 244780.636  495254.887] atol=PREC_EN
        @test convert_bng(-2.938277411, 54.89542340) ≈ [ 339921.145  556034.761] atol=PREC_EN
        @test convert_bng(-1.616576852, 54.97912274) ≈ [ 424639.355  565012.703] atol=PREC_EN
        @test convert_bng(-4.296490163, 55.85399953) ≈ [ 256340.925  664697.269] atol=PREC_EN
        @test convert_bng(-3.294792193, 55.92478266) ≈ [ 319188.434  670947.534] atol=PREC_EN
        @test convert_bng(-5.828366919, 57.00606696) ≈ [ 167634.202  797067.144] atol=PREC_EN
        @test convert_bng(-2.048560307, 57.13902519) ≈ [ 397160.491  805349.736] atol=PREC_EN
        @test convert_bng(-4.219263986, 57.48625001) ≈ [ 267056.768  846176.972] atol=PREC_EN
        @test convert_bng(-8.578544561, 57.81351838) ≈ [   9587.909  899448.996] atol=PREC_EN
        @test convert_bng(-7.592555606, 58.21262247) ≈ [  71713.132  938516.404] atol=PREC_EN
        @test convert_bng(-6.260914555, 58.51560361) ≈ [ 151968.652  966483.780] atol=PREC_EN
        @test convert_bng(-3.726310221, 58.58120461) ≈ [ 299721.891  967202.992] atol=PREC_EN
        @test convert_bng(-3.214540011, 59.03743871) ≈ [ 330398.323 1017347.016] atol=PREC_EN
        @test convert_bng(-4.417576746, 59.09335035) ≈ [ 261596.778 1025447.602] atol=PREC_EN
        @test convert_bng(-5.827993398, 59.09671617) ≈ [ 180862.461 1029604.114] atol=PREC_EN
        @test convert_bng(-1.625169661, 59.53470794) ≈ [ 421300.525 1072147.239] atol=PREC_EN
        @test convert_bng(-1.274869104, 59.85409914) ≈ [ 440725.073 1107878.448] atol=PREC_EN
        @test convert_bng(-2.073828228, 60.13308092) ≈ [ 395999.668 1138728.951] atol=PREC_EN
    end

    @testset "OS OSGB36 to latlon testset" begin
        # TP01 to TP40
        @test convert_lonlat(  91492.146,   11318.804) ≈ [-6.299777520 49.92226394] atol=PREC_LL
        @test convert_lonlat( 170370.718,   11572.405) ≈ [-5.203046100 49.96006138] atol=PREC_LL
        @test convert_lonlat( 250359.811,   62016.569) ≈ [-4.108645636 50.43885826] atol=PREC_LL
        @test convert_lonlat( 449816.371,   75335.861) ≈ [-1.297822772 50.57563665] atol=PREC_LL
        @test convert_lonlat( 438710.920,  114792.250) ≈ [-1.450514337 50.93127938] atol=PREC_LL
        @test convert_lonlat( 292184.870,  168003.465) ≈ [-3.551283492 51.40078220] atol=PREC_LL
        @test convert_lonlat( 639821.835,  169565.858) ≈ [ 1.444547304 51.37447026] atol=PREC_LL
        @test convert_lonlat( 362269.991,  169978.690) ≈ [-2.544076183 51.42754743] atol=PREC_LL
        @test convert_lonlat( 530624.974,  178388.464) ≈ [-0.119925572 51.48936565] atol=PREC_LL
        @test convert_lonlat( 241124.584,  220332.641) ≈ [-4.308524770 51.85890896] atol=PREC_LL
        @test convert_lonlat( 599445.590,  225722.826) ≈ [ 0.897243270 51.89436637] atol=PREC_LL
        @test convert_lonlat( 389544.190,  261912.153) ≈ [-2.154586144 52.25529382] atol=PREC_LL
        @test convert_lonlat( 474335.969,  262047.755) ≈ [-0.912489570 52.25160951] atol=PREC_LL
        @test convert_lonlat( 562180.547,  319784.995) ≈ [ 0.401535471 52.75136687] atol=PREC_LL
        @test convert_lonlat( 454002.834,  340834.943) ≈ [-1.197476559 52.96219109] atol=PREC_LL
        @test convert_lonlat( 357455.843,  383290.436) ≈ [-2.640493208 53.34480280] atol=PREC_LL
        @test convert_lonlat( 247958.971,  393492.909) ≈ [-4.289180698 53.41628516] atol=PREC_LL
        @test convert_lonlat( 247959.241,  393495.583) ≈ [-4.289177929 53.41630925] atol=PREC_LL
        @test convert_lonlat( 331534.564,  431920.794) ≈ [-3.040454907 53.77911026] atol=PREC_LL
        @test convert_lonlat( 422242.186,  433818.701) ≈ [-1.663791682 53.80021520] atol=PREC_LL
        @test convert_lonlat( 227778.330,  468847.388) ≈ [-4.634521682 54.08666318] atol=PREC_LL
        @test convert_lonlat( 525745.670,  470703.214) ≈ [-0.077731332 54.11685144] atol=PREC_LL
        @test convert_lonlat( 244780.636,  495254.887) ≈ [-4.388491181 54.32919541] atol=PREC_LL
        @test convert_lonlat( 339921.145,  556034.761) ≈ [-2.938277411 54.89542340] atol=PREC_LL
        @test convert_lonlat( 424639.355,  565012.703) ≈ [-1.616576852 54.97912274] atol=PREC_LL
        @test convert_lonlat( 256340.925,  664697.269) ≈ [-4.296490163 55.85399953] atol=PREC_LL
        @test convert_lonlat( 319188.434,  670947.534) ≈ [-3.294792193 55.92478266] atol=PREC_LL
        @test convert_lonlat( 167634.202,  797067.144) ≈ [-5.828366919 57.00606696] atol=PREC_LL
        @test convert_lonlat( 397160.491,  805349.736) ≈ [-2.048560307 57.13902519] atol=PREC_LL
        @test convert_lonlat( 267056.768,  846176.972) ≈ [-4.219263986 57.48625001] atol=PREC_LL
        @test convert_lonlat(   9587.909,  899448.996) ≈ [-8.578544561 57.81351838] atol=PREC_LL
        @test convert_lonlat(  71713.132,  938516.404) ≈ [-7.592555606 58.21262247] atol=PREC_LL
        @test convert_lonlat( 151968.652,  966483.780) ≈ [-6.260914555 58.51560361] atol=PREC_LL
        @test convert_lonlat( 299721.891,  967202.992) ≈ [-3.726310221 58.58120461] atol=PREC_LL
        @test convert_lonlat( 330398.323, 1017347.016) ≈ [-3.214540011 59.03743871] atol=PREC_LL
        @test convert_lonlat( 261596.778, 1025447.602) ≈ [-4.417576746 59.09335035] atol=PREC_LL
        @test convert_lonlat( 180862.461, 1029604.114) ≈ [-5.827993398 59.09671617] atol=PREC_LL
        @test convert_lonlat( 421300.525, 1072147.239) ≈ [-1.625169661 59.53470794] atol=PREC_LL
        @test convert_lonlat( 440725.073, 1107878.448) ≈ [-1.274869104 59.85409914] atol=PREC_LL
        @test convert_lonlat( 395999.668, 1138728.951) ≈ [-2.073828228 60.13308092] atol=PREC_LL
    end

end


text:
#Script for plotting blockade shift results
#26/07/2017

using Plots, JLD, LaTeXStrings
pyplot()
include("functions.jl")
PyPlot.close("all")

atom = "87Rb"
nn = 50
ll = 1
jj = 0.5
mj = 0.5

RRSI, θ, blockadeshiftmeshGHz, C_6val = BlockadeShift(atom,nn,ll,jj,mj)

PyPlot.figure()
heatmap(RRSI*1e6,θ,blockadeshiftmeshGHz*1e3)
plot!(xlabel=L"R\, (μ\mathrm{m})")
plot!(ylabel=L"θ")
plot!(zlabel=L"ΔW\, (\mathrm{MHz})")
plot!(colorbar=false)
plot!(title=string("Calculated Rydberg blockade shift\n for |n,l,j,mⱼ⟩ = |",string(nn),",",string(ll),",",string(jj),",",string(mj),"⟩ state of ",atom))
gui()
# savefig("SaveFiles/BlockadeSurface.pdf")

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
module lattice

export Lattice,LatticeD2Q9,LatticeD2Q4,LatticeD2Q5

abstract type Lattice end

# D2Q9 lattice
struct LatticeD2Q9 <: Lattice
  Δx::AbstractFloat
  Δt::AbstractFloat
  cc::AbstractFloat
  cs::AbstractFloat
  cssq::AbstractFloat
  f::Array{Float64, 3}
  c::Matrix{Int64}
  w::Vector{Float64}
  nw::Int
  nI::Int               # This should not belong to lattice?
  nJ::Int

  function LatticeD2Q9(dx::AbstractFloat, dt::AbstractFloat, nI::Int, nJ::Int,
                       rho::AbstractFloat=0.0)
    # Default lattice speed vectors and associated weights
    cdef = [1 0; 0 1; -1 0; 0 -1; 1 1; -1 1; -1 -1; 1 -1; 0 0]'
    wdef = [1.0/9.0; 1.0/9.0; 1.0/9.0; 1.0/9.0; 1.0/36.0; 1.0/36.0;
            1.0/36.0; 1.0/36.0; 4.0/9.0]
    nw = length(wdef)
    # Lattice functions for computing `c`, `cs`, and `cssq`
    cf(dx, dt)    = dx/dt
    csf(dx, dt)   = cf(dx, dt) / sqrt(3)
    cssqf(dx, dt) = csf(dx, dt) * csf(dx, dt)
    f = zeros(Float64, (nw, nI, nJ))
    for k=1:nw
      f[k,:,:] = fill(rho * wdef[k], (nI, nJ))
    end
    new(dx, dt, cf(dx, dt), csf(dx, dt), cssqf(dx, dt), f, cdef, wdef,
      nw, nI, nJ)
  end

end

# D2Q4 lattice
struct LatticeD2Q4 <: Lattice
  Δx::AbstractFloat     # Spatial interval
  Δt::AbstractFloat     # Time interval
  cc::AbstractFloat     # Lattice speed
  cs::AbstractFloat     # Sound speed
  cssq::AbstractFloat   # Sound speed squared
  f::Array{Float64, 3}  # Distribution
  c::Matrix{Int64}      # Speed vectors
  w::Vector{Float64}    # Weights
  nw::Int               # Number of velocity directions
  nI::Int               # This should not belong to lattice?
  nJ::Int

  function LatticeD2Q4(dx::AbstractFloat, dt::AbstractFloat, nI::Int, nJ::Int,
                       rho::Float64=0.0)
    # Default lattice speed vectors and associated weights
    cdef = [1 0; -1 0; 0 1; 0 -1]'
    wdef = [1.0/4.0; 1.0/4.0; 1.0/4.0; 1.0/4.0]
    nw = size(cdef, 2)

    # Lattice functions for computing `c`, `cs`, and `cssq`
    cf(dx, dt)    = dx/dt
    csf(dx, dt)   = cf(dx, dt) / sqrt(2)
    cssqf(dx, dt) = csf(dx, dt) * csf(dx, dt)

    f = zeros(Float64, (nw, nI, nJ))
    for k=1:nw
      f[k,:,:] = fill(rho * wdef[k], (nI, nJ))
    end
    new(dx, dt, cf(dx, dt), csf(dx, dt), cssqf(dx, dt), f, cdef, wdef,
      nw, nI, nJ)
  end

#  LatticeD2Q4(dx::AbstractFloat, dt::AbstractFloat, f::Array{Float64, 3}) =
#    new(dx, dt, cf(dx, dt), csf(dx, dt), cssqf(dx, dt), f, cdef, wdef, n)
end

# D2Q4 lattice
struct LatticeD2Q5 <: Lattice
  Δx::AbstractFloat     # Spatial interval
  Δt::AbstractFloat     # Time interval
  cc::AbstractFloat     # Lattice speed
  cs::AbstractFloat     # Sound speed
  cssq::AbstractFloat   # Sound speed squared
  f::Array{Float64, 3}  # Distribution
  c::Matrix{Int64}      # Speed vectors
  w::Vector{Float64}    # Weights
  nw::Int               # Number of velocity directions
  nI::Int               # This should not belong to lattice?
  nJ::Int

  function LatticeD2Q5(dx::AbstractFloat, dt::AbstractFloat, nI::Int, nJ::Int,
                       rho::Float64=0.0)
    # Default lattice speed vectors and associated weights
    cdef = [1 0; -1 0; 0 1; 0 -1; 0 0]'
    wdef = [2.0/12.0; 2.0/12.0; 2.0/12.0; 2.0/12.0; 4.0/12.0]
    nw = size(cdef, 2)

    # Lattice functions for computing `c`, `cs`, and `cssq`
    cf(dx, dt)    = dx/dt
    csf(dx, dt)   = cf(dx, dt) / sqrt(3)
    cssqf(dx, dt) = csf(dx, dt) * csf(dx, dt)

    f = zeros(Float64, (nw, nI, nJ))
    for k=1:nw
      f[k,:,:] = fill(rho * wdef[k], (nI, nJ))
    end
    new(dx, dt, cf(dx, dt), csf(dx, dt), cssqf(dx, dt), f, cdef, wdef,
      nw, nI, nJ)
  end

#  LatticeD2Q4(dx::AbstractFloat, dt::AbstractFloat, f::Array{Float64, 3}) =
#    new(dx, dt, cf(dx, dt), csf(dx, dt), cssqf(dx, dt), f, cdef, wdef, n)
end

# Opposite direction lattice vector index
#
# \param lat Lattice
# \param k Index of lattice vector
# \return Index of opposite direction
function opp_lat_vec(lat::LatticeD2Q9, k::Int)
  if k == 1
    return 3
  elseif k == 2
    return 4
  elseif k == 3
    return 1
  elseif k == 4
    return 2
  elseif k == 5
    return 7
  elseif k == 6
    return 8
  elseif k == 7
    return 5
  elseif k == 8
    return 6
  elseif k == 9
    return 9
  else
    error("$k > 9 for a D2Q9 lattice. Only nine vectors (1-9) possible");
  end
end

# Opposite direction lattice vector index
#
# \param lat Lattice
# \param k Index of lattice vector
# \return Index of opposite direction
function opp_lat_vec(lat::LatticeD2Q4, k::Int)
if k == 1
    return 2
elseif k == 2
    return 1
elseif k == 3
    return 4
elseif k == 4
    return 3
else
    error("$k > 4 for a D2Q4 lattice. Only four vectors (1-4) possible")
end

end

#=
# Helper function for initializing lattice
function _fill_lat(lat::Lattice, i_range::UnitRange{Int}, j_range::UnitRange{Int},
                   rho::Real)
  for k=1:lat.n, j=j_range, i=i_range
    lat.f[k, i, j] = rho * lat.w[k];
  end
end
=#

end


text:
importall Bukdu

type WelcomeController <: ApplicationController
end

first(::WelcomeController) = 1
second(::WelcomeController) = 2

Logger.set_level(false)

using Base.Test
@test isempty(Bukdu.RouterRoute.routes)

@test_throws NoRouteError (Endpoint)("/")

Router() do
    get("/", WelcomeController, first)
end

@test !isempty(Bukdu.RouterRoute.routes)

conn = (Router)(get, "/")
@test 200 == conn.status

@test_throws NoRouteError (Endpoint)("/")

type SecondRouter <: ApplicationRouter
end

SecondRouter() do
end

@test isempty(Bukdu.RouterRoute.routes)

SecondRouter() do
    get("/", WelcomeController, second)
end

@test !isempty(Bukdu.RouterRoute.routes)

Endpoint() do
    plug(Router)
    plug(SecondRouter)
end

conn = (Router)(get, "/")
@test 200 == conn.status

conn = (SecondRouter)(get, "/")
@test 200 == conn.status

conn = (Endpoint)("/")
@test 1 == conn.resp_body

type SecondEndpoint <: ApplicationEndpoint
end

SecondEndpoint() do
    plug(SecondRouter)
    plug(Router)
end

conn = (Endpoint)("/")
@test 1 == conn.resp_body

conn = (SecondEndpoint)("/")
@test 2 == conn.resp_body

type NothingRouter <: ApplicationRouter
end

type NothingEndpoint <: ApplicationEndpoint
end

NothingRouter() do
end

NothingEndpoint() do
end

@test_throws NoRouteError (NothingRouter)(get, "/")

@test_throws NoRouteError (NothingEndpoint)("/")

NothingEndpoint() do
    plug(NothingRouter)
end

@test_throws NoRouteError (NothingEndpoint)("/")

conn = (Endpoint)("/")
@test 200 == conn.status


text:
log_temperature(t::Real) = 1 / log(t)

constant_temperature(t::Real) = 1.0

function default_neighbor!(x::Array, x_proposal::Array)
    @assert size(x) == size(x_proposal)
    for i in 1:length(x)
        @inbounds x_proposal[i] = x[i] + randn()
    end
    return
end

macro satrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
            end
            grnorm = NaN
            update!(tr,
                    iteration,
                    f_x,
                    grnorm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable SimulatedAnnealing <: Optimizer
    neighbor!::Function
    temperature::Function
    keep_best::Bool # not used!?
end

SimulatedAnnealing(; neighbor!::Function = default_neighbor!,
                     temperature::Function = log_temperature,
                     keep_best::Bool = true) =
  SimulatedAnnealing(neighbor!, temperature, keep_best)

function optimize{T}(cost::Function,
                     initial_x::Array{T},
                     mo::SimulatedAnnealing,
                     o::OptimizationOptions)
    # Print header if show_trace is set
    print_header(o)

    # Maintain current and proposed state
    x, x_proposal = copy(initial_x), copy(initial_x)

    # Record the number of iterations we perform
    iteration = 0

    # Track calls to function and gradient
    f_calls = 0

    # Count number of parameters
    n = length(x)

    # Store f(x) in f_x
    f_x = cost(x)
    f_calls += 1

    # Store the best state ever visited
    best_x = copy(x)
    best_f_x = f_x

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @satrace

    # We always perform a fixed number of iterations
    while iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Determine the temperature for current iteration
        t = mo.temperature(iteration)

        # Randomly generate a neighbor of our current state
        mo.neighbor!(x, x_proposal)

        # Evaluate the cost function at the proposed state
        f_proposal = cost(x_proposal)
        f_calls += 1

        if f_proposal <= f_x
            # If proposal is superior, we always move to it
            copy!(x, x_proposal)
            f_x = f_proposal

            # If the new state is the best state yet, keep a record of it
            if f_proposal < best_f_x
                best_f_x = f_proposal
                copy!(best_x, x_proposal)
            end
        else
            # If proposal is inferior, we move to it with probability p
            p = exp(-(f_proposal - f_x) / t)
            if rand() <= p
                copy!(x, x_proposal)
                f_x = f_proposal
            end
        end

        @satrace
    end

    return MultivariateOptimizationResults("Simulated Annealing",
                                           initial_x,
                                           best_x,
                                           Float64(best_f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           false,
                                           NaN,
                                           false,
                                           NaN,
                                           false,
                                           NaN,
                                           tr,
                                           f_calls,
                                           0)
end


text:
using UnsteadyFlowSolvers

#alphadef = SinDef(deg2rad(3), deg2rad(10), 1, 0)
alphadef = ConstDef(deg2rad(1))
hdef = ConstDef(0.)
udef = ConstDef(1.)
full_kinem = KinemDef(alphadef, hdef, udef)

pvt = 0.25
geometry = "FlatPlate"
surf = TwoDSurf(geometry, pvt, full_kinem, rho = 0.02)

Re = 10000
nu = udef.amp  * surf.c / Re
curfield = TwoDVFlowField(nu)
#curfield = TwoDFlowField()
dtstar = find_tstep(alphadef)
t_tot = 3 * dtstar
nsteps = Int(round(t_tot/dtstar))+1
startflag = 0
writeflag = 1
writeInterval = t_tot/10.
delvort = delNone()

mat, surf, curfield = lautat(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort, wakerollup=1)

#makeForcePlots2D()
makeVortPlots2D()

return


text:
@doc doc"""
    AbstractGroupOperation

Abstract type for smooth binary operations $∘$ on elements of a Lie group $G$:
```math
∘: G × G → G
```
An operation can be either defined for a specific [`AbstractGroupManifold`](@ref)
or in general, by defining for an operation `Op` the following methods:

    identity!(::AbstractGroupManifold{Op}, y, x)
    identity(::AbstractGroupManifold{Op}, x)
    inv!(::AbstractGroupManifold{Op}, y, x)
    inv(::AbstractGroupManifold{Op}, x)
    compose(::AbstractGroupManifold{Op}, x, y)
    compose!(::AbstractGroupManifold{Op}, z, x, y)

Note that a manifold is connected with an operation by wrapping it with a decorator,
[`AbstractGroupManifold`](@ref). In typical cases the concrete wrapper
[`GroupManifold`](@ref) can be used.
"""
abstract type AbstractGroupOperation end

@doc doc"""
    AbstractGroupManifold{<:AbstractGroupOperation} <: Manifold

Abstract type for a Lie group, a group that is also a smooth manifold with an
[`AbstractGroupOperation`](@ref), a smooth binary operation. `AbstractGroupManifold`s must
implement at least [`inv`](@ref), [`identity`](@ref), [`compose`](@ref), and
[`translate_diff`](@ref).
"""
abstract type AbstractGroupManifold{O<:AbstractGroupOperation} <: Manifold end

"""
    GroupManifold{M<:Manifold,O<:AbstractGroupOperation} <: AbstractGroupManifold{O}

Decorator for a smooth manifold that equips the manifold with a group operation, thus making
it a Lie group. See [`AbstractGroupManifold`](@ref) for more details.

Group manifolds by default forward metric-related operations to the wrapped manifold.

# Constructor

    GroupManifold(manifold, op)
"""
struct GroupManifold{M<:Manifold,O<:AbstractGroupOperation} <: AbstractGroupManifold{O}
    manifold::M
    op::O
end

show(io::IO, G::GroupManifold) = print(io, "GroupManifold($(G.manifold), $(G.op))")

is_decorator_manifold(::GroupManifold) = Val(true)

is_decorator_group(::AbstractGroupManifold) = Val(true)
is_decorator_group(M::Manifold) = is_decorator_group(M, is_decorator_manifold(M))
is_decorator_group(M::Manifold, ::Val{true}) = is_decorator_group(M.manifold)
is_decorator_group(::Manifold, ::Val{false}) = Val(false)

"""
    base_group(M::Manifold) -> AbstractGroupManifold

Undecorate `M` until an `AbstractGroupManifold` is encountered.
Return an error if the [`base_manifold`](@ref) is reached without encountering a group.
"""
function base_group(M::Manifold)
    is_decorator_group(M) === Val(true) && return base_group(M.manifold)
    error("base_group: manifold $(typeof(M)) with base manifold $(typeof(base_manifold(M))) has no base group.")
end
base_group(G::AbstractGroupManifold) = G

# piping syntax for decoration
if VERSION ≥ v"1.3"
    (op::AbstractGroupOperation)(M::Manifold) = GroupManifold(M, op)
    (::Type{T})(M::Manifold) where {T<:AbstractGroupOperation} = GroupManifold(M, T())
end

########################
# GroupManifold forwards
########################

function check_tangent_vector(G::GroupManifold, x, v; kwargs...)
    return check_tangent_vector(G.manifold, x, v; kwargs...)
end
distance(G::GroupManifold, x, y) = distance(G.manifold, x, y)
exp(G::GroupManifold, x, v) = exp(G.manifold, x, v)
exp!(G::GroupManifold, y, x, v) = exp!(G.manifold, y, x, v)
injectivity_radius(G::GroupManifold) = injectivity_radius(G.manifold)
injectivity_radius(G::GroupManifold, x) = injectivity_radius(G.manifold, x)
function injectivity_radius(G::GroupManifold, x, method::AbstractRetractionMethod)
    return injectivity_radius(G.manifold, x, method)
end
inner(G::GroupManifold, x, v, w) = inner(G.manifold, x, v, w)
inverse_retract(G::GroupManifold, x, y) = inverse_retract(G.manifold, x, y)
function inverse_retract(G::GroupManifold, x, y, method::AbstractInverseRetractionMethod)
    return inverse_retract(G.manifold, x, y, method)
end
inverse_retract!(G::GroupManifold, v, x, y) = inverse_retract!(G.manifold, v, x, y)
function inverse_retract!(
    G::GroupManifold,
    v,
    x,
    y,
    method::AbstractInverseRetractionMethod,
)
    return inverse_retract!(G.manifold, v, x, y, method)
end
function inverse_retract!(G::GroupManifold, v, x, y, method::LogarithmicInverseRetraction)
    return inverse_retract!(G.manifold, v, x, y, method)
end
isapprox(G::GroupManifold, x, y; kwargs...) = isapprox(G.manifold, x, y; kwargs...)
isapprox(G::GroupManifold, x, v, w; kwargs...) = isapprox(G.manifold, x, v, w; kwargs...)
log(G::GroupManifold, x, y) = log(G.manifold, x, y)
log!(G::GroupManifold, v, x, y) = log!(G.manifold, v, x, y)
norm(G::GroupManifold, x, v) = norm(G.manifold, x, v)
project_point(G::GroupManifold, x) = project_point(G.manifold, x)
project_point!(G::GroupManifold, y, x) = project_point!(G.manifold, y, x)
project_tangent(G::GroupManifold, x, v) = project_tangent(G.manifold, x, v)
project_tangent!(G::GroupManifold, w, x, v) = project_tangent!(G.manifold, w, x, v)
retract(G::GroupManifold, x, v) = retract(G.manifold, x, v)
function retract(G::GroupManifold, x, v, method::AbstractRetractionMethod)
    return retract(G.manifold, x, v, method)
end
retract!(G::GroupManifold, y, x, v) = retract!(G.manifold, y, x, v)
function retract!(G::GroupManifold, y, x, v, method::AbstractRetractionMethod)
    return retract!(G.manifold, y, x, v, method)
end
function retract!(G::GroupManifold, y, x, v, method::ExponentialRetraction)
    return retract!(G.manifold, y, x, v, method)
end
function vector_transport_along!(G::GroupManifold, vto, x, v, c, args...)
    return vector_transport_along!(G.manifold, vto, x, v, c, args...)
end
function vector_transport_along(G::GroupManifold, x, v, c, args...)
    return vector_transport_along(G.manifold, x, v, c, args...)
end
function vector_transport_direction!(G::GroupManifold, vto, x, v, vdir, args...)
    return vector_transport_direction!(G.manifold, vto, x, v, vdir, args...)
end
function vector_transport_direction(G::GroupManifold, x, v, vdir, args...)
    return vector_transport_direction(G.manifold, x, v, vdir, args...)
end
function vector_transport_to!(G::GroupManifold, vto, x, v, y, args...)
    return vector_transport_to!(G.manifold, vto, x, v, y, args...)
end
function vector_transport_to(G::GroupManifold, x, v, y, args...)
    return vector_transport_to(G.manifold, x, v, y, args...)
end
zero_tangent_vector(G::GroupManifold, x) = zero_tangent_vector(G.manifold, x)
zero_tangent_vector!(G::GroupManifold, y, x) = zero_tangent_vector!(G.manifold, y, x)

###################
# Action directions
###################

"""
    ActionDirection

Direction of action on a manifold, either [`LeftAction`](@ref) or [`RightAction`](@ref).
"""
abstract type ActionDirection end

"""
    LeftAction()

Left action of a group on a manifold.
"""
struct LeftAction <: ActionDirection end

"""
    RightAction()

Right action of a group on a manifold.
"""
struct RightAction <: ActionDirection end

"""
    switch_direction(::ActionDirection)

Returns a [`RightAction`](@ref) when given a [`LeftAction`](@ref) and vice versa.
"""
switch_direction(::ActionDirection)
switch_direction(::LeftAction) = RightAction()
switch_direction(::RightAction) = LeftAction()

##################################
# General Identity element methods
##################################

@doc doc"""
    Identity(G::AbstractGroupManifold)

The group identity element $e ∈ G$.
"""
struct Identity{G<:AbstractGroupManifold}
    group::G
end

Identity(M::Manifold) = Identity(M, is_decorator_manifold(M))
Identity(M::Manifold, ::Val{true}) = Identity(M.manifold)
Identity(M::Manifold, ::Val{false}) = error("Identity not implemented for manifold $(M)")

show(io::IO, e::Identity) = print(io, "Identity($(e.group))")

(e::Identity)(x) = identity(e.group, x)

# To ensure allocate_result_type works
number_eltype(e::Identity) = Bool

copyto!(e::TE, ::TE) where {TE<:Identity} = e
copyto!(x, ::TE) where {TE<:Identity} = identity!(e.group, x, e)
copyto!(x::AbstractArray, e::TE) where {TE<:Identity} = identity!(e.group, x, e)

isapprox(x, e::Identity; kwargs...) = isapprox(e::Identity, x; kwargs...)
isapprox(e::Identity, x; kwargs...) = isapprox(e.group, e, x; kwargs...)
isapprox(e::E, ::E; kwargs...) where {E<:Identity} = true

function check_manifold_point(M::Manifold, x::Identity; kwargs...)
    if is_decorator_group(M) === Val(true)
        return check_manifold_point(base_group(M), x; kwargs...)
    end
    return DomainError(x, "The identity element $(x) does not belong to $(M).")
end
function check_manifold_point(G::GroupManifold, x::Identity; kwargs...)
    x === Identity(G) && return nothing
    return DomainError(x, "The identity element $(x) does not belong to $(G).")
end
function check_manifold_point(G::GroupManifold, x; kwargs...)
    return check_manifold_point(G.manifold, x; kwargs...)
end

##########################
# Group-specific functions
##########################

@doc doc"""
    inv(G::AbstractGroupManifold, x)

Inverse $x^{-1} ∈ G$ of an element $x ∈ G$, such that
$x \circ x^{-1} = x^{-1} \circ x = e ∈ G$.
"""
inv(M::Manifold, x) = inv(M, x, is_decorator_manifold(M))
inv(M::Manifold, x, ::Val{true}) = inv(M.manifold, x)
function inv(M::Manifold, x, ::Val{false})
    return error("inv not implemented on $(typeof(M)) for points $(typeof(x))")
end
function inv(G::AbstractGroupManifold, x)
    y = allocate_result(G, inv, x)
    return inv!(G, y, x)
end

@doc doc"""
    inv!(G::AbstractGroupManifold, y, x)

Inverse $x^{-1} ∈ G$ of an element $x ∈ G$, such that
$x \circ x^{-1} = x^{-1} \circ x = e ∈ G$.
The result is saved to `y`.
"""
inv!(M::Manifold, y, x) = inv!(M, y, x, is_decorator_manifold(M))
inv!(M::Manifold, y, x, ::Val{true}) = inv!(M.manifold, y, x)
function inv!(M::Manifold, y, x, ::Val{false})
    return error("inv! not implemented on $(typeof(M)) for points $(typeof(x))")
end

@doc doc"""
    identity(G::AbstractGroupManifold, x)

Identity element $e ∈ G$, such that for any element $x ∈ G$, $x \circ e = e \circ x = x$.
The returned element is of a similar type to `x`.
"""
identity(M::Manifold, x) = identity(M, x, is_decorator_manifold(M))
identity(M::Manifold, x, ::Val{true}) = identity(M.manifold, x)
function identity(M::Manifold, x, ::Val{false})
    return error("identity not implemented on $(typeof(M)) for points $(typeof(x))")
end
function identity(G::AbstractGroupManifold, x)
    y = allocate_result(G, identity, x)
    return identity!(G, y, x)
end

identity!(M::Manifold, y, x) = identity!(M, y, x, is_decorator_manifold(M))
identity!(M::Manifold, y, x, ::Val{true}) = identity!(M.manifold, y, x)
function identity!(M::Manifold, y, x, ::Val{false})
    return error("identity! not implemented on $(typeof(M)) for points $(typeof(y)) and $(typeof(x))")
end

isapprox(M::Manifold, x, e::Identity; kwargs...) = isapprox(M, e, x; kwargs...)
function isapprox(M::Manifold, e::Identity, x; kwargs...)
    is_decorator_group(M) === Val(true) && return isapprox(base_group(M), e, x; kwargs...)
    error("isapprox not implemented for manifold $(typeof(M)) and points $(typeof(e)) and $(typeof(x))")
end
function isapprox(M::Manifold, e::E, ::E; kwargs...) where {E<:Identity}
    is_decorator_group(M) === Val(true) && return isapprox(base_group(M), e, e; kwargs...)
    error("isapprox not implemented for manifold $(typeof(M)) and points $(typeof(e)) and $(typeof(e))")
end
function isapprox(G::GT, e::Identity{GT}, x; kwargs...) where {GT<:AbstractGroupManifold}
    return isapprox(G, identity(G, x), x; kwargs...)
end
function isapprox(
    ::GT,
    ::E,
    ::E;
    kwargs...,
) where {GT<:AbstractGroupManifold,E<:Identity{GT}}
    return true
end
function isapprox(G::GT, x, e::Identity{GT}; kwargs...) where {GT<:GroupManifold}
    return isapprox(G, e, x; kwargs...)
end
function isapprox(G::GT, e::Identity{GT}, x; kwargs...) where {GT<:GroupManifold}
    return isapprox(G, identity(G, x), x; kwargs...)
end
isapprox(::GT, ::E, ::E; kwargs...) where {GT<:GroupManifold,E<:Identity{GT}} = true

@doc doc"""
    compose(G::AbstractGroupManifold, x, y)

Compose elements $x,y ∈ G$ using the group operation $x \circ y$.
"""
compose(M::Manifold, x, y) = compose(M, x, y, is_decorator_manifold(M))
compose(M::Manifold, x, y, ::Val{true}) = compose(M.manifold, x, y)
function compose(M::Manifold, x, y, ::Val{false})
    return error("compose not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y))")
end
function compose(G::AbstractGroupManifold, x, y)
    z = allocate_result(G, compose, x, y)
    return compose!(G, z, x, y)
end

compose!(M::Manifold, z, x, y) = compose!(M, z, x, y, is_decorator_manifold(M))
compose!(M::Manifold, z, x, y, ::Val{true}) = compose!(M.manifold, z, x, y)
function compose!(M::Manifold, z, x, y, ::Val{false})
    return error("compose! not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y))")
end

_action_order(x, y, conv::LeftAction) = (x, y)
_action_order(x, y, conv::RightAction) = (y, x)

@doc doc"""
    translate(G::AbstractGroupManifold, x, y[, conv::ActionDirection=LeftAction()])

For group elements $x,y ∈ G$, translate $y$ by $x$ with the specified convention, either
left $L_x$ or right $R_x$, defined as
```math
\begin{aligned}
L_x &: y ↦ x \circ y\\
R_x &: y ↦ y \circ x.
\end{aligned}
```
"""
translate(M::Manifold, x, y) = translate(M, x, y, LeftAction())
function translate(M::Manifold, x, y, conv::ActionDirection)
    return translate(M, x, y, conv, is_decorator_manifold(M))
end
function translate(M::Manifold, x, y, conv::ActionDirection, ::Val{true})
    return translate(M.manifold, x, y, conv)
end
function translate(M::Manifold, x, y, conv::ActionDirection, ::Val{false})
    return error("translate not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y)) and direction $(typeof(conv))")
end
function translate(G::AbstractGroupManifold, x, y, conv::ActionDirection)
    return compose(G, _action_order(x, y, conv)...)
end

@doc doc"""
    translate!(G::AbstractGroupManifold, z, x, y[, conv::ActionDirection=LeftAction()])

For group elements $x,y ∈ G$, translate $y$ by $x$ with the specified convention, either
left $L_x$ or right $R_x$, defined as
```math
\begin{aligned}
L_x &: y ↦ x \circ y\\
R_x &: y ↦ y \circ x.
\end{aligned}
```
Result of the operation is saved in `z`.
"""
translate!(M::Manifold, z, x, y) = translate!(M, z, x, y, LeftAction())
function translate!(M::Manifold, z, x, y, conv::ActionDirection)
    return translate!(M, z, x, y, conv, is_decorator_manifold(M))
end
function translate!(M::Manifold, z, x, y, conv::ActionDirection, ::Val{true})
    return translate!(M.manifold, z, x, y, conv)
end
function translate!(M::Manifold, z, x, y, conv::ActionDirection, ::Val{false})
    return error("translate! not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y)) and direction $(typeof(conv))")
end
function translate!(G::AbstractGroupManifold, z, x, y, conv::ActionDirection)
    return compose!(G, z, _action_order(x, y, conv)...)
end

@doc doc"""
    inverse_translate(G::AbstractGroupManifold, x, y, [conv::ActionDirection=Left()])

For group elements $x,y ∈ G$, inverse translate $y$ by $x$ with the specified convention,
either left $L_x^{-1}$ or right $R_x^{-1}$, defined as
```math
\begin{aligned}
L_x^{-1} &: y ↦ x^{-1} \circ y\\
R_x^{-1} &: y ↦ y \circ x^{-1}.
\end{aligned}
```
"""
inverse_translate(M::Manifold, x, y) = inverse_translate(M, x, y, LeftAction())
function inverse_translate(M::Manifold, x, y, conv::ActionDirection)
    return inverse_translate(M, x, y, conv, is_decorator_manifold(M))
end
function inverse_translate(M::Manifold, x, y, conv::ActionDirection, ::Val{true})
    return inverse_translate(M.manifold, x, y, conv)
end
function inverse_translate(M::Manifold, x, y, conv::ActionDirection, ::Val{false})
    return error("inverse_translate not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y)) and direction $(typeof(conv))")
end
function inverse_translate(G::AbstractGroupManifold, x, y, conv::ActionDirection)
    return translate(G, inv(G, x), y, conv)
end

@doc doc"""
    inverse_translate!(G::AbstractGroupManifold, z, x, y, [conv::ActionDirection=Left()])

For group elements $x,y ∈ G$, inverse translate $y$ by $x$ with the specified convention,
either left $L_x^{-1}$ or right $R_x^{-1}$, defined as
```math
\begin{aligned}
L_x^{-1} &: y ↦ x^{-1} \circ y\\
R_x^{-1} &: y ↦ y \circ x^{-1}.
\end{aligned}
```
Result is saved in `z`.
"""
inverse_translate!(M::Manifold, z, x, y) = inverse_translate!(M, z, x, y, LeftAction())
function inverse_translate!(M::Manifold, z, x, y, conv::ActionDirection)
    return inverse_translate!(M, z, x, y, conv, is_decorator_manifold(M))
end
function inverse_translate!(M::Manifold, z, x, y, conv::ActionDirection, ::Val{true})
    return inverse_translate!(M.manifold, z, x, y, conv)
end
function inverse_translate!(M::Manifold, z, x, y, conv::ActionDirection, ::Val{false})
    return error("inverse_translate! not implemented on $(typeof(M)) for elements $(typeof(x)) and $(typeof(y)) and direction $(typeof(conv))")
end
function inverse_translate!(G::AbstractGroupManifold, z, x, y, conv::ActionDirection)
    return translate!(G, z, inv(G, x), y, conv)
end

@doc doc"""
    translate_diff(G::AbstractGroupManifold, x, y, v[, conv::ActionDirection=LeftAction()])

For group elements $x,y ∈ G$ and tangent vector $v ∈ T_y G$, compute the action of the
differential of the translation by $x$ on $v$, written as $(\mathrm{d}τ_x)_y (v)$, with the
specified left or right convention. The differential transports vectors:
```math
\begin{aligned}
(\mathrm{d}L_x)_y (v) &: T_y G → T_{x \circ y} G\\
(\mathrm{d}R_x)_y (v) &: T_y G → T_{y \circ x} G\\
\end{aligned}
```
"""
translate_diff(M::Manifold, x, y, v) = translate_diff(M, x, y, v, LeftAction())
function translate_diff(M::Manifold, x, y, v, conv::ActionDirection)
    return translate_diff(M, x, y, v, conv, is_decorator_manifold(M))
end
function translate_diff(M::Manifold, x, y, v, conv::ActionDirection, ::Val{true})
    return translate_diff(M.manifold, x, y, v, conv)
end
function translate_diff(M::Manifold, x, y, v, conv::ActionDirection, ::Val{false})
    return error("translate_diff not implemented on $(typeof(G)) for elements $(typeof(vout)), $(typeof(x)) and $(typeof(y)), vector $(typeof(v)), and direction $(typeof(conv))")
end
function translate_diff(G::AbstractGroupManifold, x, y, v, conv::ActionDirection)
    xy = translate(G, x, y, conv)
    vout = zero_tangent_vector(G, xy)
    translate_diff!(G, vout, x, y, v, conv)
    return vout
end

function translate_diff!(M::Manifold, vout, x, y, v)
    return translate_diff!(M, vout, x, y, v, LeftAction())
end
function translate_diff!(M::Manifold, vout, x, y, v, conv::ActionDirection)
    return translate_diff!(M, vout, x, y, v, conv, is_decorator_manifold(M))
end
function translate_diff!(M::Manifold, vout, x, y, v, conv::ActionDirection, ::Val{true})
    return translate_diff!(M.manifold, vout, x, y, v, conv)
end
function translate_diff!(M::Manifold, vout, x, y, v, conv::ActionDirection, ::Val{false})
    return error("translate_diff! not implemented on $(typeof(M)) for elements $(typeof(vout)), $(typeof(x)) and $(typeof(y)), vector $(typeof(v)), and direction $(typeof(conv))")
end

@doc doc"""
    inverse_translate_diff(G::AbstractGroupManifold, x, y, v[, conv::ActionDirection=Left()])

For group elements $x,y ∈ G$ and tangent vector $v ∈ T_y G$, compute the inverse of the
action of the differential of the translation by $x$ on $v$, written as
$((\mathrm{d}τ_x)_y)^{-1} (v) = (\mathrm{d}τ_{x^{-1}})_y (v)$, with the specified left or
right convention. The differential transports vectors:
```math
\begin{aligned}
((\mathrm{d}L_x)_y)^{-1} (v) &: T_y G → T_{x^{-1} \circ y} G\\
((\mathrm{d}R_x)_y)^{-1} (v) &: T_y G → T_{y \circ x^{-1}} G\\
\end{aligned}
```
"""
function inverse_translate_diff(M::Manifold, x, y, v)
    return inverse_translate_diff(M, x, y, v, LeftAction())
end
function inverse_translate_diff(M::Manifold, x, y, v, conv::ActionDirection)
    return inverse_translate_diff(M, x, y, v, conv, is_decorator_manifold(M))
end
function inverse_translate_diff(M::Manifold, x, y, v, conv::ActionDirection, ::Val{true})
    return inverse_translate_diff(M.manifold, x, y, v, conv)
end
function inverse_translate_diff(M::Manifold, x, y, v, conv::ActionDirection, ::Val{false})
    return error("inverse_translate_diff not implemented on $(typeof(M)) for elements $(typeof(vout)), $(typeof(x)) and $(typeof(y)), vector $(typeof(v)), and direction $(typeof(conv))")
end
function inverse_translate_diff(G::AbstractGroupManifold, x, y, v, conv::ActionDirection)
    return translate_diff(G, inv(G, x), y, v, conv)
end

function inverse_translate_diff!(M::Manifold, vout, x, y, v)
    return inverse_translate_diff!(M, vout, x, y, v, LeftAction())
end
function inverse_translate_diff!(M::Manifold, vout, x, y, v, conv::ActionDirection)
    return inverse_translate_diff!(M, vout, x, y, v, conv, is_decorator_manifold(M))
end
function inverse_translate_diff!(
    M::Manifold,
    vout,
    x,
    y,
    v,
    conv::ActionDirection,
    ::Val{true},
)
    return inverse_translate_diff!(M.manifold, vout, x, y, v, conv)
end
function inverse_translate_diff!(
    M::Manifold,
    vout,
    x,
    y,
    v,
    conv::ActionDirection,
    ::Val{false},
)
    return error("inverse_translate_diff! not implemented on $(typeof(M)) for elements $(typeof(vout)), $(typeof(x)) and $(typeof(y)), vector $(typeof(v)), and direction $(typeof(conv))")
end
function inverse_translate_diff!(
    G::AbstractGroupManifold,
    vout,
    x,
    y,
    v,
    conv::ActionDirection,
)
    return translate_diff!(G, vout, inv(G, x), y, v, conv)
end

#################################
# Overloads for AdditionOperation
#################################

"""
    AdditionOperation <: AbstractGroupOperation

Group operation that consists of simple addition.
"""
struct AdditionOperation <: AbstractGroupOperation end

const AdditionGroup = AbstractGroupManifold{AdditionOperation}

+(e::Identity{G}) where {G<:AdditionGroup} = e
+(::Identity{G}, x) where {G<:AdditionGroup} = x
+(x, ::Identity{G}) where {G<:AdditionGroup} = x
+(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

-(e::Identity{G}) where {G<:AdditionGroup} = e
-(::Identity{G}, x) where {G<:AdditionGroup} = -x
-(x, ::Identity{G}) where {G<:AdditionGroup} = x
-(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

*(e::Identity{G}, x) where {G<:AdditionGroup} = e
*(x, e::Identity{G}) where {G<:AdditionGroup} = e
*(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e

zero(e::Identity{G}) where {G<:AdditionGroup} = e

identity(::AdditionGroup, x) = zero(x)

identity!(::AdditionGroup, y, x) = fill!(y, 0)

inv(::AdditionGroup, x) = -x

inv!(::AdditionGroup, y, x) = copyto!(y, -x)

compose(::AdditionGroup, x, y) = x + y

function compose!(::GT, z, x, y) where {GT<:AdditionGroup}
    x isa Identity{GT} && return copyto!(z, y)
    y isa Identity{GT} && return copyto!(z, x)
    z .= x .+ y
    return z
end

translate_diff(::AdditionGroup, x, y, v, ::ActionDirection) = v

translate_diff!(::AdditionGroup, vout, x, y, v, ::ActionDirection) = copyto!(vout, v)

inverse_translate_diff(::AdditionGroup, x, y, v, ::ActionDirection) = v

function inverse_translate_diff!(::AdditionGroup, vout, x, y, v, ::ActionDirection)
    return copyto!(vout, v)
end

#######################################
# Overloads for MultiplicationOperation
#######################################

"""
    MultiplicationOperation <: AbstractGroupOperation

Group operation that consists of multiplication.
"""
struct MultiplicationOperation <: AbstractGroupOperation end

const MultiplicationGroup = AbstractGroupManifold{MultiplicationOperation}

*(e::Identity{G}) where {G<:MultiplicationGroup} = e
*(::Identity{G}, x) where {G<:MultiplicationGroup} = x
*(x, ::Identity{G}) where {G<:MultiplicationGroup} = x
*(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e

/(x, ::Identity{G}) where {G<:MultiplicationGroup} = x
/(::Identity{G}, x) where {G<:MultiplicationGroup} = inv(x)
/(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e

\(x, ::Identity{G}) where {G<:MultiplicationGroup} = inv(x)
\(::Identity{G}, x) where {G<:MultiplicationGroup} = x
\(e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}} = e

inv(e::Identity{G}) where {G<:MultiplicationGroup} = e

one(e::Identity{G}) where {G<:MultiplicationGroup} = e

transpose(e::Identity{G}) where {G<:MultiplicationGroup} = e

LinearAlgebra.det(::Identity{<:MultiplicationGroup}) = 1

LinearAlgebra.mul!(y, e::Identity{G}, x) where {G<:MultiplicationGroup} = copyto!(y, x)
LinearAlgebra.mul!(y, x, e::Identity{G}) where {G<:MultiplicationGroup} = copyto!(y, x)
function LinearAlgebra.mul!(y, e::E, ::E) where {G<:MultiplicationGroup,E<:Identity{G}}
    return identity!(e.group, y, e)
end

identity(::MultiplicationGroup, x) = one(x)

function identity!(G::GT, y, x) where {GT<:MultiplicationGroup}
    isa(x, Identity{GT}) || return copyto!(y, one(x))
    error("identity! not implemented on $(typeof(G)) for points $(typeof(y)) and $(typeof(x))")
end
identity!(::MultiplicationGroup, y::AbstractMatrix, x) = copyto!(y, I)

inv(::MultiplicationGroup, x) = inv(x)

inv!(G::MultiplicationGroup, y, x) = copyto!(y, inv(G, x))

compose(::MultiplicationGroup, x, y) = x * y

# TODO: z might alias with x or y, we might be able to optimize it if it doesn't.
compose!(::MultiplicationGroup, z, x, y) = copyto!(z, x * y)

inverse_translate(::MultiplicationGroup, x, y, ::LeftAction) = x \ y
inverse_translate(::MultiplicationGroup, x, y, ::RightAction) = y / x

function inverse_translate!(G::MultiplicationGroup, z, x, y, conv::ActionDirection)
    return copyto!(z, inverse_translate(G, x, y, conv))
end


text:
struct Word{T}
    alphabet::Vector{T}
    letters::Vector{Int}

    function Word(a::AbstractVector{T}, l::AbstractVector{<:Integer}) where T
        all(i->1<=i<=length(a), l) || throw(ArgumentError("Invalid word over alphabet $a: $w"))
        return new{T}(a,l)
    end
end
Base.show(io::IO, w::Word) = join(io, w.alphabet[w.letters], "·")

Base.:(==)(w::Word, v::Word) = w.alphabet == v.alphabet && w.letters == v.letters
Base.hash(w::Word, h::UInt=UInt(0)) = hash(w.alphabet, hash(w.letters, hash(Word, h)))

struct OnLetters <: SymbolicWedderburn.ByPermutations end
function SymbolicWedderburn.action(::OnLetters, p::PermutationGroups.AbstractPerm, w::Word)
    return Word(w.alphabet, [l^p for l in w.letters])
end

@testset "Extending homomorphism" begin
    words = let A = [:a, :b, :c], radius = 4
        w = Word(A, [1,2,3,2,1])

        # (a·b·c·b·a)^(2,3) == a·c·b·c·a
        @test SymbolicWedderburn.action(OnLetters(), perm"(2,3)", w) == Word(A, [1,3,2,3,1])

        words = [Word(A, [1]), Word(A, [2]), Word(A, [3])]
        for r in 2:radius
            append!(
                words,
                [Word(A, collect(w)) for w in Iterators.product(fill(1:3, r)...)]
            )
        end
        words
    end

    G = PermGroup(perm"(1,2,3)",perm"(1,2)") # G acts on words permuting letters
    basis = words
    action = OnLetters()
    tbl = SymbolicWedderburn.CharacterTable(Rational{Int}, G)
    ehom = SymbolicWedderburn.CachedExtensionHomomorphism(G, action, basis, precompute=true)
    @test all(g ∈ keys(ehom.cache) for g in G) # we actually cached

    ψ = SymbolicWedderburn.action_character(ehom, tbl)
    @test SymbolicWedderburn.constituents(ψ) == [40, 22, 18]
    irr = SymbolicWedderburn.irreducible_characters(tbl)
    multips = SymbolicWedderburn.constituents(ψ)
    @test dot(SymbolicWedderburn.degree.(irr), multips) == length(basis)
    simple = isone.(SymbolicWedderburn.degree.(irr))
    @test simple == [false, true, true]

    inv_vec = SymbolicWedderburn.invariant_vectors(tbl, action, SymbolicWedderburn.basis(ehom))
    @test length(inv_vec) == 22
    @test eltype(inv_vec) == SparseVector{Rational{Int}}

    @testset "semisimple decomposition" begin
        let i = 1
            χ, m, s = irr[i], multips[i], simple[i]
            b = SymbolicWedderburn.image_basis(ehom, χ)
            @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 80
        end

        let i = 2
            χ, m, s = irr[i], multips[i], simple[i]
            b = SymbolicWedderburn.image_basis(ehom, χ)
            @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 22
        end

        let i = 3
            χ, m, s = irr[i], multips[i], simple[i]
            b = SymbolicWedderburn.image_basis(ehom, χ)
            @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 18
        end

        @test symmetry_adapted_basis(G, action, basis, semisimple=true) isa
            Vector{<:SymbolicWedderburn.DirectSummand{<:SymbolicWedderburn.Cyclotomic}}
        @test symmetry_adapted_basis(Rational{Int}, G, action, basis, semisimple=true) isa
            Vector{<:SymbolicWedderburn.DirectSummand{Rational{Int}}}
        @test symmetry_adapted_basis(Float64, G, action, basis, semisimple=true) isa
            Vector{<:SymbolicWedderburn.DirectSummand{Float64}}

        sa_basis = symmetry_adapted_basis(G, action, basis, semisimple=true)

        @test [convert(Matrix{Rational{Int}}, b) for b in sa_basis] isa Vector{Matrix{Rational{Int}}}
        @test [convert(Matrix{Float64}, b) for b in sa_basis] isa Vector{Matrix{Float64}}

        @test length(sa_basis) == 3
        @test multiplicity.(sa_basis) == [40, 22, 18]
        @test SymbolicWedderburn.degree.(sa_basis) == [2,1,1]
        @test size.(sa_basis, 1) ==
            multips .* SymbolicWedderburn.degree.(irr) ==
            [80, 22, 18]
        @test sum(first ∘ size, sa_basis) == length(words)
    end

    @testset "simple decomposition" begin
        RG = let G = G
            b = StarAlgebras.Basis{UInt16}(collect(G))
            StarAlgebra(G, b, (length(b), length(b)))
        end

        let χ = irr[1], m = multips[1]
            (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG)
            @test fl == isone(χ(a))

            µ = AlgebraElement(χ, RG)*a
            mpr = SymbolicWedderburn.image_basis(ehom, µ)
            @test mpr isa AbstractMatrix{eltype(µ)}
            @test size(mpr, 1) == m

            µR = AlgebraElement{Rational{Int}}(µ)
            mpr = SymbolicWedderburn.image_basis(ehom, µR)
            @test mpr isa AbstractMatrix{eltype(µR)}
            @test size(mpr, 1) == m

            µFl = AlgebraElement{Float64}(µ)
            mpr = SymbolicWedderburn.image_basis(ehom, µFl)
            @test mpr isa AbstractMatrix{eltype(µFl)}
            @test size(mpr, 1) == m
        end

        let χ = irr[2], m = multips[2]
            (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG)
            @test fl == isone(χ(a))

            µ = AlgebraElement(χ, RG)*a
            mpr = SymbolicWedderburn.image_basis(ehom, µ)
            @test mpr isa AbstractMatrix{eltype(µ)}
            @test size(mpr, 1) == m
        end

        let χ = irr[3], m = multips[3]
            (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG)
            @test fl == isone(χ(a))

            µ = AlgebraElement(χ, RG)*a
            mpr = SymbolicWedderburn.image_basis(ehom, µ)
            @test mpr isa AbstractMatrix{eltype(µ)}
            @test size(mpr, 1) == m
        end

        @test symmetry_adapted_basis(G, action, basis, semisimple=false) isa
            Vector{<:SymbolicWedderburn.DirectSummand{<:SymbolicWedderburn.Cyclotomic}}
        @test symmetry_adapted_basis(Rational{Int}, G, action, basis, semisimple=false) isa
            Vector{<:SymbolicWedderburn.DirectSummand{Rational{Int}}}
        @test symmetry_adapted_basis(Float64, G, action, basis, semisimple=false) isa
            Vector{<:SymbolicWedderburn.DirectSummand{Float64}}

        sa_basis = symmetry_adapted_basis(G, action, basis)

        @test length(sa_basis) == 3
        @test multiplicity.(sa_basis) == [40, 22, 18]
        @test SymbolicWedderburn.degree.(sa_basis) == [2,1,1]
        @test all(issimple, sa_basis)
        @test size.(sa_basis, 1) == multips == [40, 22, 18]
    end
end


text:
function localbasis(porder,dim,elemtype)

pelem,telem,pface,tface,perm = masternodes(porder,dim,elemtype);

npe = size(pelem,1);
npf = size(pface,1);

if dim==1
    nve = 2;
    nvf = 1;
elseif dim==2 && elemtype==0 # tri
    nve = 3;
    nvf = 2;
elseif dim==2 && elemtype==1 # quad
    nve = 4;
    nvf = 2;
elseif dim==3 && elemtype==0 # tet
    nve = 4;
    nvf = 3;
elseif dim==3 && elemtype==1 # hex
    nve = 8;
    nvf = 4;
end

phielem = zeros(Float64, npe, nve);
phiface = zeros(Float64, npf, nvf);

if dim==1
    xi  = pelem[:,1];
    phielem[:,1] = 1.0 .- xi;
    phielem[:,2] = xi;
    phiface = [1.0];
elseif dim==2 && elemtype==0 # tri
    xi  = pelem[:,1];
    eta = pelem[:,2];
    phielem[:,1] = 1.0 .- xi .- eta;
    phielem[:,2] = xi;
    phielem[:,3] = eta;

    xi = pface[:,1];
    phiface[:,1] = 1.0 .- xi;
    phiface[:,2] = xi;
elseif dim==2 && elemtype==1 # quad
    xi  = pelem[:,1];
    eta = pelem[:,2];
    phielem[:,1] = (1.0 .- xi).*(1.0 .- eta);
    phielem[:,2] = xi.*(1.0 .- eta);
    phielem[:,3] = xi.*eta;
    phielem[:,4] = (1.0 .- xi).*eta;

    xi = pface[:,1];
    phiface[:,1] = 1.0 .- xi;
    phiface[:,2] = xi;
elseif dim==3 && elemtype==0 # tet
    xi   = pelem[:,1];
    eta  = pelem[:,2];
    zeta = pelem[:,3];
    phielem[:,1] = 1.0 .- xi .- eta .- zeta;
    phielem[:,2] = xi;
    phielem[:,3] = eta;
    phielem[:,4] = zeta;

    xi = pface[:,1];
    eta = pface[:,2];
    phiface[:,1] = 1.0 .- xi .- eta;
    phiface[:,2] = xi;
    phiface[:,3] = eta;
elseif dim==3 && elemtype==1 # hex
    xi   = pelem[:,1];
    eta  = pelem[:,2];
    zeta = pelem[:,3];
    phielem[:,1] = (1.0 .-xi).*(1.0 .-eta).*(1.0 .-zeta);
    phielem[:,2] = xi.*(1.0 .-eta).*(1.0 .-zeta);
    phielem[:,3] = xi.*eta.*(1.0 .-zeta);
    phielem[:,4] = (1.0 .-xi).*eta.*(1.0 .-zeta);
    phielem[:,5] = (1.0 .-xi).*(1.0 .-eta).*(zeta);
    phielem[:,6] = xi.*(1.0 .-eta).*(zeta);
    phielem[:,7] = xi.*eta.*(zeta);
    phielem[:,8] = (1.0 .-xi).*eta.*(zeta);

    xi = pface[:,1];
    eta = pface[:,2];
    phiface[:,1] = (1.0 .- xi).*(1.0 .- eta);
    phiface[:,2] = xi.*(1.0 .- eta);
    phiface[:,3] = xi.*eta;
    phiface[:,4] = (1.0 .- xi).*eta;
end

return phielem,phiface,pelem,pface,perm

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
abstract type AbstractDevice end

export CPU
struct CPU <: AbstractDevice end

export GPU
struct GPU <: AbstractDevice end

max_stars_repo_name: JuliaVlasov/VlasovSolvers

max_issues_repo_name: JuliaVlasov/VlasovSolvers

max_forks_repo_name: JuliaVlasov/VlasovSolvers


text:
"""
    randindex(a, b) -> NTuple{4, Float64}

Compute the tuple of Rand-related indices between the clusterings `c1` and `c2`.

`a` and `b` can be either [`ClusteringResult`](@ref) instances or
assignments vectors (`AbstractVector{<:Integer}`).

Returns a tuple of indices:
  - Hubert & Arabie Adjusted Rand index
  - Rand index (agreement probability)
  - Mirkin's index (disagreement probability)
  - Hubert's index (``P(\\mathrm{agree}) - P(\\mathrm{disagree})``)

# References
> Lawrence Hubert and Phipps Arabie (1985). *Comparing partitions.*
> Journal of Classification 2 (1): 193–218

> Meila, Marina (2003). *Comparing Clusterings by the Variation of
> Information.* Learning Theory and Kernel Machines: 173–187.
"""
function randindex(a, b)
    c = counts(a, b)

    n = sum(c)
    nis = sum(abs2, sum(c, dims=2))        # sum of squares of sums of rows
    njs = sum(abs2, sum(c, dims=1))        # sum of squares of sums of columns

    t1 = binomial(n, 2)                    # total number of pairs of entities
    t2 = sum(abs2, c)                      # sum over rows & columnns of nij^2
    t3 = .5*(nis+njs)

    # Expected index (for adjustment)
    nc = (n*(n^2+1)-(n+1)*nis-(n+1)*njs+2*(nis*njs)/n)/(2*(n-1))

    A = t1+t2-t3;        # agreements count
    D = -t2+t3;          # disagreements count

    if t1 == nc
        # avoid division by zero; if k=1, define Rand = 0
        ARI = 0
    else
        # adjusted Rand - Hubert & Arabie 1985
        ARI = (A-nc)/(t1-nc)
    end

    RI = A/t1            # Rand 1971      # Probability of agreement
    MI = D/t1            # Mirkin 1970    # p(disagreement)
    HI = (A-D)/t1        # Hubert 1977    # p(agree)-p(disagree)

    return (ARI, RI, MI, HI)
end