text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function annuity(years::Int, rate_escalation::Float64, rate_discount::Float64)
    """
        this formulation assumes cost growth in first period
        i.e. it is a geometric sum of (1+rate_escalation)^n / (1+rate_discount)^n
        for n = 1, ..., years
    """
    x = (1 + rate_escalation) / (1 + rate_discount)
    if x != 1
        pwf = round(x * (1 - x^years) / (1 - x), digits=5)
    else
        pwf = years
    end
    return pwf
end

function annuity_escalation(analysis_period::Int, rate_escalation::Float64, rate_discount::Float64)
    """
    :param analysis_period: years
    :param rate_escalation: escalation rate
    :param rate_discount: discount rate
    :return: present worth factor with escalation (inflation, or degradation if negative)
    NOTE: assumes escalation/degradation starts in year 2 (unlike the `annuity` function above)
    """
    pwf = 0
    for yr in range(1, stop=analysis_period + 1)
        pwf += (1 + rate_escalation)^(yr - 1) / (1 + rate_discount)^yr
    end
    return pwf
end

function levelization_factor(years::Int, rate_escalation::Float64, rate_discount::Float64,
    rate_degradation::Float64)
    #=
    NOTE: levelization_factor for an electricity producing tech is the ratio of:
    - an annuity with an escalation rate equal to the electricity cost escalation rate, starting year 1,
        and a negative escalation rate (the tech's degradation rate), starting year 2
    - divided by an annuity with an escalation rate equal to the electricity cost escalation rate (pwf_e).
    Both use the offtaker's discount rate.
    levelization_factor is multiplied by each use of dvRatedProduction in reopt.jl
        (except dvRatedProduction[t,ts] == dvSize[t] ∀ ts).
    This way the denominator is cancelled in reopt.jl when accounting for the value of energy produced
    since each value constraint uses pwf_e.

    :param analysis_period: years
    :param rate_escalation: escalation rate
    :param rate_discount: discount rate
    :param rate_degradation: positive degradation rate
    :return: present worth factor with escalation (inflation, or degradation if negative)
    NOTE: assume escalation/degradation starts in year 2
    =#
    num = 0
    for yr in range(1, stop=years)
        num += (1 + rate_escalation)^(yr) / (1 + rate_discount)^yr * (1 - rate_degradation)^(yr - 1)
    end
    den = annuity(years, rate_escalation, rate_discount)

    return num/den
end

function effective_cost(;
    itc_basis::Float64,
    replacement_cost::Float64,
    replacement_year::Int,
    discount_rate::Float64,
    tax_rate::Float64,
    itc::Float64,
    macrs_schedule::Array{Float64,1},
    macrs_bonus_pct::Float64,
    macrs_itc_reduction::Float64,
    rebate_per_kw::Float64=0.0,
    )

    """ effective PV and battery prices with ITC and depreciation
        (i) depreciation tax shields are inherently nominal --> no need to account for inflation
        (ii) ITC and bonus depreciation are taken at end of year 1
        (iii) battery replacement cost: one time capex in user defined year discounted back to t=0 with r_owner
        (iv) Assume that cash incentives reduce ITC basis
        (v) Assume cash incentives are not taxable, (don't affect tax savings from MACRS)
        (vi) Cash incentives should be applied before this function into "itc_basis".
             This includes all rebates and percentage-based incentives besides the ITC
    """

    # itc reduces depreciable_basis
    depr_basis = itc_basis * (1 - macrs_itc_reduction * itc)

    # Bonus depreciation taken from tech cost after itc reduction ($/kW)
    bonus_depreciation = depr_basis * macrs_bonus_pct

    # Assume the ITC and bonus depreciation reduce the depreciable basis ($/kW)
    depr_basis -= bonus_depreciation

    # Calculate replacement cost, discounted to the replacement year accounting for tax deduction
    replacement = replacement_cost * (1-tax_rate) / ((1 + discount_rate)^replacement_year)

    # Compute savings from depreciation and itc in array to capture NPV
    tax_savings_array = [0.0]
    for (idx, macrs_rate) in enumerate(macrs_schedule)
        depreciation_amount = macrs_rate * depr_basis
        if idx == 1
            depreciation_amount += bonus_depreciation
        end
        taxable_income = depreciation_amount
        push!(tax_savings_array, taxable_income * tax_rate)
    end

    # Add the ITC to the tax savings
    tax_savings_array[2] += itc_basis * itc

    # Compute the net present value of the tax savings
    tax_savings = npv(discount_rate, tax_savings_array)

    # Adjust cost curve to account for itc and depreciation savings ($/kW)
    cap_cost_slope = itc_basis - tax_savings + replacement - rebate_per_kw

    # Sanity check
    if cap_cost_slope < 0
        cap_cost_slope = 0
    end

    return round(cap_cost_slope, digits=4)
end

function dictkeys_tosymbols(d::Dict)
    d2 = Dict()
    for (k, v) in d
        # handling some type conversions for API inputs and JSON
        if k in [
            "loads_kw", "critical_loads_kw",
            "monthly_totals_kwh",
            "prod_factor_series",
            "monthly_energy_rates", "monthly_demand_rates",
            "wholesale_rate", "blended_doe_reference_percents",
            "coincident_peak_load_charge_per_kw", "fuel_cost_per_mmbtu"
            ] && !isnothing(v)
            try
                v = convert(Array{Real, 1}, v)
            catch
                @warn "Unable to convert $k to an Array{Real, 1}"
            end
        end
        if k in [
            "blended_doe_reference_names"
        ]
            try
                v = convert(Array{String, 1}, v)
            catch
                @warn "Unable to convert $k to an Array{String, 1}"
            end
        end
        if k in [
            "coincident_peak_load_active_timesteps"
        ]
            try
                v = convert(Vector{Vector{Int64}}, v)
            catch
                @warn "Unable to convert $k to a Vector{Vector{Int64}}"
            end
        end
        d2[Symbol(k)] = v
    end
    return d2
end

function filter_dict_to_match_struct_field_names(d::Dict, s::DataType)
    f = fieldnames(s)
    d2 = Dict()
    for k in f
        if haskey(d, k)
            d2[k] = d[k]
        else
            @warn "dict is missing struct field $k"
        end
    end
    return d2
end

function npv(rate::Float64, cash_flows::Array)
    npv = cash_flows[1]
    for (y, c) in enumerate(cash_flows[2:end])
        npv += c/(1+rate)^y
    end
    return npv
end

"""
    per_hour_value_to_time_series(x::T, time_steps_per_hour::Int) where T <: Real

Convert a per hour value (eg. dollars/kWh) to time series that matches the settings.time_steps_per_hour
"""
function per_hour_value_to_time_series(x::T, time_steps_per_hour::Int, name::String) where T <: Real
    repeat([x / time_steps_per_hour], 8760 * time_steps_per_hour)
end

"""
    per_hour_value_to_time_series(x::AbstractVector{<:Real}, time_steps_per_hour::Int, name::String)

Convert a monthly or time-sensitive per hour value (eg. dollars/kWh) to a time series that matches the
settings.time_steps_per_hour.
"""
function per_hour_value_to_time_series(x::AbstractVector{<:Real}, time_steps_per_hour::Int, name::String)
    if length(x) == 8760 * time_steps_per_hour
        return x
    end
    vals = Real[]
    if length(x) == 12  # assume monthly values
        for mth in 1:12
            append!(vals, repeat(
                [x[mth] / time_steps_per_hour],
                time_steps_per_hour * 24 * daysinmonth(Date("2017-" * string(mth)))
                )
            )
        end
        return vals
    end
    @error "Cannot convert $name to appropriate length time series."
end

"""
    generate_year_profile_hourly(year::Int64, consecutive_periods::AbstractVector{Dict})

This function creates a year-specific hourly (8760) profile with 1.0 value for timesteps which are defined in `consecutive_periods` based on
    relative (non-year specific) datetime metrics. All other values are 0.0. This functions uses the `Dates` package.

- `year` applies the relative calendar-based `consecutive_periods` to the year's calendar and handles leap years by truncating the last day
- `consecutive_periods` is a list of dictionaries where each dict defines a consecutive period of time which gets a value of 1.0
-- keys for each dict must include "month", "start_week_of_month", "start_day_of_week", "start_hour", "duration_hours
- Returns the `year_profile_hourly` which is an 8760 profile with 1.0 for timesteps defined in consecutive_periods, and 0.0 for all other hours.
"""
function generate_year_profile_hourly(year::Int64, consecutive_periods::AbstractVector{Dict})
    # Create datetime series of the year, remove last day of the year if leap year
    if Dates.isleapyear(year)
        end_year_datetime = DateTime(string(year)*"-12-30T23:00:00")
    else
        end_year_datetime = DateTime(string(year+1)*"-12-31T23:00:00")
    end

    dt_hourly = collect(DateTime(string(year)*"-01-01T00:00:00"):Hour(1):end_year_datetime)

    year_profile_hourly = zeros(8760)

    # Note, day = 1 is Monday, not Sunday
    day_of_week_name = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
    for i in 1:length(consecutive_periods)
        start_month = convert(Int,consecutive_periods[i]["month"])
        start_week_of_month = convert(Int,consecutive_periods[i]["start_week_of_month"])
        start_day_of_week = convert(Int,consecutive_periods[i]["start_day_of_week"])  # Monday - Sunday is 1 - 7
        start_hour = convert(Int,consecutive_periods[i]["start_hour"])
        duration_hours = convert(Int,consecutive_periods[i]["duration_hours"])
        error_start_text = "Error in chp.unavailability_period $(i)."

        try
            start_date_of_month_year = Date(Dates.Year(year), Dates.Month(start_month))
            start_date = Dates.firstdayofweek(start_date_of_month_year) + Dates.Week(start_week_of_month - 1) + Dates.Day(start_day_of_week - 1)
            # Throw an error if start_date is in the previous month when start_week_of_month=1 and there is no start_day_of_week in the first week of the month.
            if Dates.month(start_date) != start_month
                @error "For $(error_start_text), there is no day $(start_day_of_week) ($(day_of_week_name[start_day_of_week])) in the first week of month $(start_month) ($(Dates.monthname(start_date))), $(year)"
            end
            start_datetime = Dates.DateTime(start_date) + Dates.Hour(start_hour - 1)
            if Dates.year(start_datetime + Dates.Hour(duration_hours)) > year
                @error "For $(error_start_text), the start day/time and duration_hours exceeds the end of the year. Please specify two separate unavailability periods: one for the beginning of the year and one for up to the end of the year."
            else
                #end_datetime is the last hour that is 1.0 (e.g. that is still unavailable), not the first hour that is 0.0 after the period
                end_datetime = start_datetime + Dates.Hour(duration_hours - 1)
                year_profile_hourly[findfirst(x->x==start_datetime, dt_hourly):findfirst(x->x==end_datetime, dt_hourly)] .= 1.0
            end
        catch e
            println("For $error_start_text, invalid set for month $start_month (1-12), start_week_of_month $start_week_of_month (1-4, possible 5 and 6), $start_day_of_week (1-7), and $start_hour (1-24) for the year $year.")
        end
    end
    return year_profile_hourly
end

function get_ambient_temperature(latitude::Real, longitude::Real; timeframe="hourly")
    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
        "&lat=", latitude , "&lon=", longitude, "&tilt=", latitude,
        "&system_capacity=1", "&azimuth=", 180, "&module_type=", 0,
        "&array_type=", 0, "&losses=", 14,
        "&timeframe=", timeframe, "&dataset=nsrdb"
    )

    try
        @info "Querying PVWatts for ambient temperature... "
        r = HTTP.get(url)
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from PVWatts: $(response["errors"])")
        end
        @info "PVWatts success."
        tamb = collect(get(response["outputs"], "tamb", []))  # Celcius
        if length(tamb) != 8760
            @error "PVWatts did not return a valid temperature. Got $tamb"
        end
        return tamb
    catch e
        @error "Error occurred when calling PVWatts: $e"
    end
end

function get_pvwatts_prodfactor(latitude::Real, longitude::Real; timeframe="hourly")
    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
        "&lat=", latitude , "&lon=", longitude, "&tilt=", latitude,
        "&system_capacity=1", "&azimuth=", 180, "&module_type=", 0,
        "&array_type=", 0, "&losses=", 14,
        "&timeframe=", timeframe, "&dataset=nsrdb"
    )

    try
        @info "Querying PVWatts for production factor of 1 kW system with tilt set to latitude... "
        r = HTTP.get(url)
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from PVWatts: $(response["errors"])")
        end
        @info "PVWatts success."
        watts = collect(get(response["outputs"], "ac", []) / 1000)  # scale to 1 kW system (* 1 kW / 1000 W)
        if length(watts) != 8760
            @error "PVWatts did not return a valid prodfactor. Got $watts"
        end
        return watts
    catch e
        @error "Error occurred when calling PVWatts: $e"
    end
end


text:
import Symbolics.get_variables; export get_variables

# pretty-printing
display(vars::HarmonicVariable) = display(vars.names)
display(vars::Vector{HarmonicVariable}) = display.(getfield.(vars, Symbol("names")))

# Each keyword gives two names for the new variables
coordinate_names = Dict([["Cartesian", ["u", "v"]], ["polar", ["a", "ϕ"]]])

function _coordinate_transform(new_vars, ω, t)
    Dict([["Cartesian",  new_vars[1] * cos(ω*t) + new_vars[2] * sin(ω*t)],
["polar", new_vars[1] * cos(ω*t + new_vars[2])]])
end

"Generates two variables describing the harmonic ω of a natural variable."
function _rotate_variable(nat_var::Num, ω::Num, time::Num, transform::String; new_symbols::Vector{String})

    types = coordinate_names[transform] # is the transformation "Cartesian" or "polar"
    new_vars = [declare_variable(s, time) for s in new_symbols] # this holds the internal symbols

    var_names_strings = [VDP_name * "_{" * var_name(nat_var) * "," * Base.replace(string(ω),"*"=>"") * "}" for VDP_name in types]
    var_names = Dict(zip([Num(s) for s in new_symbols], var_names_strings))
    repl_rule = _coordinate_transform(new_vars, ω, time)[transform] # a dictionary mapping the coordinate type to an expression
    VDP = HarmonicVariable(new_vars, var_names, types, ω, nat_var)

    repl_rule, VDP
end

###
# Functions for variable substutions and manipulation of HarmonicVariable
###

function substitute_all(vars::HarmonicVariable, rules)
    sym, freq = vars.symbols, vars.ω
    HarmonicVariable(substitute_all(sym, rules), vars.names, vars.types, substitute_all(freq, rules), vars.natural_variable)
end

substitute_all(vars::Vector{HarmonicVariable}, rules) = [substitute_all(var, rules) for var in vars]

"Returns the symbols of a `HarmonicVariable`."
get_variables(vars::Vector{Num}) = unique(flatten([Num.(get_variables(x)) for x in vars]))

get_variables(vars::HarmonicVariable) = unique(get_variables(vars.symbols))

ext: jl


text:
abstract type SingleObjAquisition{S}<:EGOAquisition end

include("single_obj_improvement.jl")


text:
# This file was generated, do not modify it. # hide
mach_pipe = machine(pipe, X, y)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using AxisKeys
using PrettyTables

export AbstractCVMethod, AbstractCV

const POINTWISE_LABELS = (:cv_est, :naive_est, :overfit, :ess, :pareto_k)
const CV_DESC = """
# Fields

  - `estimates::KeyedArray`: A KeyedArray with columns `:total, :se_total, :mean, :se_mean`,
    and rows `:cv_est, :naive_est, :overfit`. See `# Extended help` for more.
      - `:cv_est` contains estimates for the out-of-sample prediction error, as
        predicted using the jackknife (LOO-CV).
      - `:naive_est` contains estimates of the in-sample prediction error.
      - `:overfit` is the difference between the previous two estimators, and estimates
        the amount of overfitting. When using the log probability score, it is equal to
        the effective number of parameters -- a model with an overfit of 2 is "about as
        overfit" as a model with 2 independent parameters that have a flat prior.
  - `pointwise::KeyedArray`: A `KeyedArray` of pointwise estimates with 5 columns --
      - `:cv_est` contains the estimated out-of-sample error for this point, as measured
      using leave-one-out cross validation.
      - `:naive_est` contains the in-sample estimate of error for this point.
      - `:overfit` is the difference in the two previous estimates.
      - `:ess` is the effective sample size, which estimates the simulation error caused by
        using Monte Carlo estimates. It does not measure the accuracy of predictions.
      - `:pareto_k` is the estimated value for the parameter `ξ` of the generalized Pareto
        distribution. Values above .7 indicate that PSIS has failed to approximate the true
        distribution.
  - `psis_object::Psis`: A `Psis` object containing the results of Pareto-smoothed
    importance sampling.

# Extended help

The total score depends on the sample size, and summarizes the weight of evidence for or
against a model. Total scores are on an interval scale, meaning that only differences of
scores are meaningful. *It is not possible to interpret a total score by looking at it.*
The total score is not a relative goodness-of-fit statistic (for this, see the average
score).

The overfit is equal to the difference between the in-sample and out-of-sample predictive
accuracy. When using the log probability score, it is equal to the "effective number of
parameters" -- a model with an overfit of 2 is "about as overfit" as a model with 2
free parameters and flat priors.

The average score is the total score, divided by the sample size. It estimates the expected
log score, i.e. the expectation of the log probability density of observing the next point.
The average score is a relative goodness-of-fit statistic which does not depend on sample
size.

Unlike for chi-square goodness of fit tests, models do not have to be nested for model
comparison using cross-validation methods.
"""

##########################
#### CROSS VALIDATION ####
##########################

"""
    AbstractCV

An abstract type used in cross-validation.
"""
abstract type AbstractCV end

"""
    AbstractCVMethod

An abstract type used to dispatch the correct method for cross validation.
"""
abstract type AbstractCVMethod end

##########################
### BAYESIAN BOOTSTRAP ###
##########################

# """
#     BayesCV{
#         F <: Real,
#         AF <: AbstractArray{F},
#         VF <: AbstractVector{F},
#         I <: Integer,
#         VI <: AbstractVector{I},
#     } <: AbstractCV

# A struct containing the results of cross-validation using the Bayesian bootstrap.

# $CV_DESC

# See also: [`bayes_cv`]@ref, [`psis_loo`]@ref, [`psis`]@ref, [`Psis`]@ref
# """
# struct BayesCV{
#     F <: Real,
#     AF <: AbstractArray{F},
#     VF <: AbstractVector{F},
#     I <: Integer,
#     VI <: AbstractVector{I},
# } <: AbstractCV
#     estimates::KeyedArray
#     posteriors::KeyedArray
#     psis_object::Psis{F, AF, VF, I, VI}
# end

# function Base.show(io::IO, ::MIME"text/plain", cv_object::BayesCV)
#     table = cv_object.estimates
#     post_samples = cv_object.psis_object.posterior_sample_size
#     data_size = cv_object.psis_object.data_size
#     _throw_pareto_k_warning(cv_object.resamples(:pareto_k))
#     println(
#         "Results of Bayesian bootstrap CV with $post_samples Monte Carlo samples and " *
#         "$data_size data points.",
#     )
#     return pretty_table(
#         table;
#         compact_printing=false,
#         header=table.statistic,
#         row_names=table.criterion,
#         formatters=ft_printf("%5.2f"),
#         alignment=:r,
#     )
# end


text:
"""
    ParametersEulerGravity(; background_density=0.0,
                             gravitational_constant=1.0,
                             cfl=1.0,
                             resid_tol=1.0e-4,
                             n_iterations_max=10^4,
                             timestep_gravity=timestep_gravity_erk52_3Sstar!)

Set up parameters for the gravitational part of a [`SemidiscretizationEulerGravity`](@ref).
"""
struct ParametersEulerGravity{RealT<:Real, TimestepGravity}
  background_density    ::RealT # aka rho0
  gravitational_constant::RealT # aka G
  cfl                   ::RealT
  resid_tol             ::RealT
  n_iterations_max      ::Int
  timestep_gravity::TimestepGravity
end

function ParametersEulerGravity(; background_density=0.0,
                                  gravitational_constant=1.0,
                                  cfl=1.0,
                                  resid_tol=1.0e-4,
                                  n_iterations_max=10^4,
                                  timestep_gravity=timestep_gravity_erk52_3Sstar!)
  background_density, gravitational_constant, cfl, resid_tol = promote(background_density, gravitational_constant, cfl, resid_tol)
  ParametersEulerGravity(background_density, gravitational_constant, cfl, resid_tol, n_iterations_max, timestep_gravity)
end

function Base.show(io::IO, parameters::ParametersEulerGravity)
  @nospecialize parameters # reduce precompilation time

  print(io, "ParametersEulerGravity(")
  print(io,   "background_density=", parameters.background_density)
  print(io, ", gravitational_constant=", parameters.gravitational_constant)
  print(io, ", cfl=", parameters.cfl)
  print(io, ", n_iterations_max=", parameters.n_iterations_max)
  print(io, ", timestep_gravity=", parameters.timestep_gravity)
  print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", parameters::ParametersEulerGravity)
  @nospecialize parameters # reduce precompilation time

  if get(io, :compact, false)
    show(io, parameters)
  else
    setup = [
             "background density (ρ₀)" => parameters.background_density,
             "gravitational constant (G)" => parameters.gravitational_constant,
             "CFL (gravity)" => parameters.cfl,
             "max. #iterations" => parameters.n_iterations_max,
             "time integrator" => parameters.timestep_gravity,
            ]
    summary_box(io, "ParametersEulerGravity", setup)
  end
end

"""
    SemidiscretizationEulerGravity

A struct containing everything needed to describe a spatial semidiscretization
of a the compressible Euler equations with self-gravity, reformulating the
Poisson equation for the gravitational potential as steady-state problem of
the hyperblic diffusion equations.
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  "A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics"
  [arXiv: 2008.10593](https://arXiv.org/abs/2008.10593)
"""
struct SemidiscretizationEulerGravity{SemiEuler, SemiGravity,
                                      Parameters<:ParametersEulerGravity, Cache} <: AbstractSemidiscretization
  semi_euler::SemiEuler
  semi_gravity::SemiGravity
  parameters::Parameters
  performance_counter::PerformanceCounter
  gravity_counter::PerformanceCounter
  cache::Cache

  function SemidiscretizationEulerGravity{SemiEuler, SemiGravity, Parameters, Cache}(
      semi_euler::SemiEuler, semi_gravity::SemiGravity,
      parameters::Parameters, cache::Cache) where {SemiEuler, SemiGravity,
                                                   Parameters<:ParametersEulerGravity, Cache}
    @assert ndims(semi_euler) == ndims(semi_gravity)
    @assert typeof(semi_euler.mesh) == typeof(semi_gravity.mesh)
    @assert polydeg(semi_euler.solver) == polydeg(semi_gravity.solver)

    performance_counter = PerformanceCounter()
    gravity_counter = PerformanceCounter()

    new(semi_euler, semi_gravity, parameters, performance_counter, gravity_counter, cache)
  end
end

"""
    SemidiscretizationEulerGravity(semi_euler::SemiEuler, semi_gravity::SemiGravity, parameters)

Construct a semidiscretization of the compressible Euler equations with self-gravity.
`parameters` should be given as [`ParametersEulerGravity`](@ref).
"""
function SemidiscretizationEulerGravity(semi_euler::SemiEuler, semi_gravity::SemiGravity, parameters) where
    {Mesh, SemiEuler<:SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations},
           SemiGravity<:SemidiscretizationHyperbolic{Mesh, <:AbstractHyperbolicDiffusionEquations}}

  u_ode = compute_coefficients(zero(real(semi_gravity)), semi_gravity)
  du_ode     = similar(u_ode)
  u_tmp1_ode = similar(u_ode)
  u_tmp2_ode = similar(u_ode)
  cache = (; u_ode, du_ode, u_tmp1_ode, u_tmp2_ode)

  SemidiscretizationEulerGravity{typeof(semi_euler), typeof(semi_gravity), typeof(parameters), typeof(cache)}(
    semi_euler, semi_gravity, parameters, cache)
end

# TODO: AD, add appropriate method for remake

function Base.show(io::IO, semi::SemidiscretizationEulerGravity)
  @nospecialize semi # reduce precompilation time

  print(io, "SemidiscretizationEulerGravity using")
  print(io,       semi.semi_euler)
  print(io, ", ", semi.semi_gravity)
  print(io, ", ", semi.parameters)
  print(io, ", cache(")
  for (idx,key) in enumerate(keys(semi.cache))
    idx > 1 && print(io, " ")
    print(io, key)
  end
  print(io, "))")
end

function Base.show(io::IO, mime::MIME"text/plain", semi::SemidiscretizationEulerGravity)
  @nospecialize semi # reduce precompilation time

  if get(io, :compact, false)
    show(io, semi)
  else
    summary_header(io, "SemidiscretizationEulerGravity")
    summary_line(io, "semidiscretization Euler", semi.semi_euler |> typeof |> nameof)
    show(increment_indent(io), mime, semi.semi_euler)
    summary_line(io, "semidiscretization gravity", semi.semi_gravity |> typeof |> nameof)
    show(increment_indent(io), mime, semi.semi_gravity)
    summary_line(io, "parameters", semi.parameters |> typeof |> nameof)
    show(increment_indent(io), mime, semi.parameters)
    summary_footer(io)
  end
end

# The compressible Euler semidiscretization is considered to be the main semidiscretization.
# The hyperbolic diffusion equations part is only used internally to update the gravitational
# potential during an rhs! evaluation of the flow solver.
@inline function mesh_equations_solver_cache(semi::SemidiscretizationEulerGravity)
  mesh_equations_solver_cache(semi.semi_euler)
end

@inline Base.ndims(semi::SemidiscretizationEulerGravity) = ndims(semi.semi_euler)

@inline Base.real(semi::SemidiscretizationEulerGravity) = real(semi.semi_euler)

# computes the coefficients of the initial condition
@inline function compute_coefficients(t, semi::SemidiscretizationEulerGravity)
  compute_coefficients!(semi.cache.u_ode, t, semi.semi_gravity)
  compute_coefficients(t, semi.semi_euler)
end

# computes the coefficients of the initial condition and stores the Euler part in `u_ode`
@inline function compute_coefficients!(u_ode, t, semi::SemidiscretizationEulerGravity)
  compute_coefficients!(semi.cache.u_ode, t, semi.semi_gravity)
  compute_coefficients!(u_ode, t, semi.semi_euler)
end

@inline function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationEulerGravity, cache_analysis)
  calc_error_norms(func, u, t, analyzer, semi.semi_euler, cache_analysis)
end

function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)
  @unpack semi_euler, semi_gravity, cache = semi

  u_euler   = wrap_array(u_ode , semi_euler)
  du_euler  = wrap_array(du_ode, semi_euler)
  u_gravity = wrap_array(cache.u_ode, semi_gravity)

  time_start = time_ns()

  # standard semidiscretization of the compressible Euler equations
  @timeit_debug timer() "Euler solver" rhs!(du_ode, u_ode, semi_euler, t)

  # compute gravitational potential and forces
  @timeit_debug timer() "gravity solver" update_gravity!(semi, u_ode)

  # add gravitational source source_terms to the Euler part
  if ndims(semi_euler) == 1
    @views @. du_euler[2, .., :] -= u_euler[1, .., :] * u_gravity[2, .., :]
    @views @. du_euler[3, .., :] -= u_euler[2, .., :] * u_gravity[2, .., :]
  elseif ndims(semi_euler) == 2
    @views @. du_euler[2, .., :] -=  u_euler[1, .., :] * u_gravity[2, .., :]
    @views @. du_euler[3, .., :] -=  u_euler[1, .., :] * u_gravity[3, .., :]
    @views @. du_euler[4, .., :] -= (u_euler[2, .., :] * u_gravity[2, .., :] +
                                     u_euler[3, .., :] * u_gravity[3, .., :])
  elseif ndims(semi_euler) == 3
    @views @. du_euler[2, .., :] -=  u_euler[1, .., :] * u_gravity[2, .., :]
    @views @. du_euler[3, .., :] -=  u_euler[1, .., :] * u_gravity[3, .., :]
    @views @. du_euler[4, .., :] -=  u_euler[1, .., :] * u_gravity[4, .., :]
    @views @. du_euler[5, .., :] -= (u_euler[2, .., :] * u_gravity[2, .., :] +
                                     u_euler[3, .., :] * u_gravity[3, .., :] +
                                     u_euler[4, .., :] * u_gravity[4, .., :])
  else
    error("Number of dimensions $(ndims(semi_euler)) not supported.")
  end

  runtime = time_ns() - time_start
  put!(semi.performance_counter, runtime)

  return nothing
end

# TODO: Taal refactor, add some callbacks or so within the gravity update to allow investigating/optimizing it
function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
  @unpack semi_euler, semi_gravity, parameters, gravity_counter, cache = semi

  # Can be changed by AMR
  resize!(cache.du_ode,     length(cache.u_ode))
  resize!(cache.u_tmp1_ode, length(cache.u_ode))
  resize!(cache.u_tmp2_ode, length(cache.u_ode))

  u_euler    = wrap_array(u_ode,        semi_euler)
  u_gravity  = wrap_array(cache.u_ode,  semi_gravity)
  du_gravity = wrap_array(cache.du_ode, semi_gravity)

  # set up main loop
  finalstep = false
  @unpack n_iterations_max, cfl, resid_tol, timestep_gravity = parameters
  iter = 0
  t = zero(real(semi_gravity.solver))

  # iterate gravity solver until convergence or maximum number of iterations are reached
  @unpack equations = semi_gravity
  while !finalstep
    dt = @timeit_debug timer() "calculate dt" cfl * max_dt(u_gravity, t, semi_gravity.mesh,
                                                           have_constant_speed(equations), equations,
                                                           semi_gravity.solver, semi_gravity.cache)

    # evolve solution by one pseudo-time step
    time_start = time_ns()
    timestep_gravity(cache, u_euler, t, dt, parameters, semi_gravity)
    runtime = time_ns() - time_start
    put!(gravity_counter, runtime)

    # update iteration counter
    iter += 1
    t += dt

    # check if we reached the maximum number of iterations
    if n_iterations_max > 0 && iter >= n_iterations_max
      @warn "Max iterations reached: Gravity solver failed to converge!" residual=maximum(abs, @views du_gravity[1, .., :]) t=t dt=dt
      finalstep = true
    end

    # this is an absolute tolerance check
    if maximum(abs, @views du_gravity[1, .., :]) <= resid_tol
      finalstep = true
    end
  end

  return nothing
end

# Integrate gravity solver for 2N-type low-storage schemes
function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
                              a, b, c)
  G    = gravity_parameters.gravitational_constant
  rho0 = gravity_parameters.background_density
  grav_scale = -4.0*pi*G

  @unpack u_ode, du_ode, u_tmp1_ode = cache
  u_tmp1_ode .= zero(eltype(u_tmp1_ode))
  du_gravity = wrap_array(du_ode, semi_gravity)
  for stage in eachindex(c)
    t_stage = t + dt * c[stage]

    # rhs! has the source term for the harmonic problem
    # We don't need a `@timeit_debug timer() "rhs!"` here since that's already
    # included in the `rhs!` call.
    rhs!(du_ode, u_ode, semi_gravity, t_stage)

    # Source term: Jeans instability OR coupling convergence test OR blast wave
    # put in gravity source term proportional to Euler density
    # OBS! subtract off the background density ρ_0 (spatial mean value)
    @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)

    a_stage = a[stage]
    b_stage_dt = b[stage] * dt
    @timeit_debug timer() "Runge-Kutta step" begin
      @threaded for idx in eachindex(u_ode)
        u_tmp1_ode[idx] = du_ode[idx] - u_tmp1_ode[idx] * a_stage
        u_ode[idx] += u_tmp1_ode[idx] * b_stage_dt
      end
    end
  end

  return nothing
end

function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
  # Coefficients for Carpenter's 5-stage 4th-order low-storage Runge-Kutta method
  a = SVector(0.0, 567301805773.0 / 1357537059087.0,2404267990393.0 / 2016746695238.0,
              3550918686646.0 / 2091501179385.0, 1275806237668.0 / 842570457699.0)
  b = SVector(1432997174477.0 / 9575080441755.0, 5161836677717.0 / 13612068292357.0,
              1720146321549.0 / 2090206949498.0, 3134564353537.0 / 4481467310338.0,
              2277821191437.0 / 14882151754819.0)
  c = SVector(0.0, 1432997174477.0 / 9575080441755.0, 2526269341429.0 / 6820363962896.0,
              2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0)

  timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity, a, b, c)
end

# Integrate gravity solver for 3S*-type low-storage schemes
function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
                                  gamma1, gamma2, gamma3, beta, delta, c)
  G    = gravity_parameters.gravitational_constant
  rho0 = gravity_parameters.background_density
  grav_scale = -4 * G * pi

  @unpack u_ode, du_ode, u_tmp1_ode, u_tmp2_ode = cache
  u_tmp1_ode .= zero(eltype(u_tmp1_ode))
  u_tmp2_ode .= u_ode
  du_gravity = wrap_array(du_ode, semi_gravity)
  for stage in eachindex(c)
    t_stage = t + dt * c[stage]

    # rhs! has the source term for the harmonic problem
    # We don't need a `@timeit_debug timer() "rhs!"` here since that's already
    # included in the `rhs!` call.
    rhs!(du_ode, u_ode, semi_gravity, t_stage)

    # Source term: Jeans instability OR coupling convergence test OR blast wave
    # put in gravity source term proportional to Euler density
    # OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed
    @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)

    delta_stage   = delta[stage]
    gamma1_stage  = gamma1[stage]
    gamma2_stage  = gamma2[stage]
    gamma3_stage  = gamma3[stage]
    beta_stage_dt = beta[stage] * dt
    @timeit_debug timer() "Runge-Kutta step" begin
      @threaded for idx in eachindex(u_ode)
        u_tmp1_ode[idx] += delta_stage * u_ode[idx]
        u_ode[idx]       = (gamma1_stage * u_ode[idx] +
                            gamma2_stage * u_tmp1_ode[idx] +
                            gamma3_stage * u_tmp2_ode[idx] +
                            beta_stage_dt * du_ode[idx])
      end
    end
  end

  return nothing
end

function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
  # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
  # and examples/parameters_hypdiff_lax_friedrichs.toml
  # 5 stages, order 1
  gamma1 = SVector(0.0000000000000000E+00, 5.2910412316555866E-01, 2.8433964362349406E-01, -1.4467571130907027E+00, 7.5592215948661057E-02)
  gamma2 = SVector(1.0000000000000000E+00, 2.6366970460864109E-01, 3.7423646095836322E-01,  7.8786901832431289E-01, 3.7754129043053775E-01)
  gamma3 = SVector(0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00,  8.0043329115077388E-01, 1.3550099149374278E-01)
  beta   = SVector(1.9189497208340553E-01, 5.4506406707700059E-02, 1.2103893164085415E-01,  6.8582252490550921E-01, 8.7914657211972225E-01)
  delta  = SVector(1.0000000000000000E+00, 7.8593091509463076E-01, 1.2639038717454840E-01,  1.7726945920209813E-01, 0.0000000000000000E+00)
  c      = SVector(0.0000000000000000E+00, 1.9189497208340553E-01, 1.9580448818599061E-01,  2.4241635859769023E-01, 5.0728347557552977E-01)

  timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
                           gamma1, gamma2, gamma3, beta, delta, c)
end

function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
  # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
  # and examples/parameters_hypdiff_lax_friedrichs.toml
  # 5 stages, order 2
  gamma1 = SVector(0.0000000000000000E+00, 5.2656474556752575E-01,  1.0385212774098265E+00, 3.6859755007388034E-01, -6.3350615190506088E-01)
  gamma2 = SVector(1.0000000000000000E+00, 4.1892580153419307E-01, -2.7595818152587825E-02, 9.1271323651988631E-02,  6.8495995159465062E-01)
  gamma3 = SVector(0.0000000000000000E+00, 0.0000000000000000E+00,  0.0000000000000000E+00, 4.1301005663300466E-01, -5.4537881202277507E-03)
  beta   = SVector(4.5158640252832094E-01, 7.5974836561844006E-01,  3.7561630338850771E-01, 2.9356700007428856E-02,  2.5205285143494666E-01)
  delta  = SVector(1.0000000000000000E+00, 1.3011720142005145E-01,  2.6579275844515687E-01, 9.9687218193685878E-01,  0.0000000000000000E+00)
  c      = SVector(0.0000000000000000E+00, 4.5158640252832094E-01,  1.0221535725056414E+00, 1.4280257701954349E+00,  7.1581334196229851E-01)

  timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
                           gamma1, gamma2, gamma3, beta, delta, c)
end

function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity)
  # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
  # and examples/parameters_hypdiff_lax_friedrichs.toml
  # 5 stages, order 3
  gamma1 = SVector(0.0000000000000000E+00,  6.9362208054011210E-01, 9.1364483229179472E-01,  1.3129305757628569E+00, -1.4615811339132949E+00)
  gamma2 = SVector(1.0000000000000000E+00,  1.3224582239681788E+00, 2.4213162353103135E-01, -3.8532017293685838E-01,  1.5603355704723714E+00)
  gamma3 = SVector(0.0000000000000000E+00,  0.0000000000000000E+00, 0.0000000000000000E+00,  3.8306787039991996E-01, -3.5683121201711010E-01)
  beta   = SVector(8.4476964977404881E-02,  3.0834660698015803E-01, 3.2131664733089232E-01,  2.8783574345390539E-01,  8.2199204703236073E-01)
  delta  = SVector(1.0000000000000000E+00, -7.6832695815481578E-01, 1.2497251501714818E-01,  1.4496404749796306E+00,  0.0000000000000000E+00)
  c      = SVector(0.0000000000000000E+00,  8.4476964977404881E-02, 2.8110631488732202E-01,  5.7093842145029405E-01,  7.2999896418559662E-01)

  timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
                           gamma1, gamma2, gamma3, beta, delta, c)
end

# TODO: Taal decide, where should specific parts like these be?
@inline function save_solution_file(u_ode, t, dt, iter,
                                    semi::SemidiscretizationEulerGravity, solution_callback,
                                    element_variables=Dict{Symbol,Any}())

  u_euler = wrap_array(u_ode, semi.semi_euler)
  filename_euler = save_solution_file(u_euler, t, dt, iter,
                                      mesh_equations_solver_cache(semi.semi_euler)...,
                                      solution_callback, element_variables, system="euler")

  u_gravity = wrap_array(semi.cache.u_ode, semi.semi_gravity)
  filename_gravity = save_solution_file(u_gravity, t, dt, iter,
                                        mesh_equations_solver_cache(semi.semi_gravity)...,
                                        solution_callback, element_variables, system="gravity")

  return filename_euler, filename_gravity
end

@inline function (amr_callback::AMRCallback)(u_ode,
                                             semi::SemidiscretizationEulerGravity,
                                             t, iter; kwargs...)
  passive_args = ((semi.cache.u_ode, mesh_equations_solver_cache(semi.semi_gravity)...),)
  amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)..., t, iter;
               kwargs..., passive_args=passive_args)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
import Pkg.Types: VersionSpec

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "polymake"
upstream_version = v"4.4"
version = VersionNumber(upstream_version.major*100,upstream_version.minor*100,0)

# Collection of sources required to build polymake
sources = [
    ArchiveSource("https://github.com/polymake/polymake/archive/V$(upstream_version.major).$(upstream_version.minor).tar.gz",
                  "6560445a3f6d359c0d1c3a01eb3ec3f3174b3554bf33865c913082030e72c176")
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
# the prepared config files for non-native architectures
# assume a fixed source directory
mv $WORKSPACE/srcdir/{polymake-*,polymake}
cd $WORKSPACE/srcdir/polymake

perl_version=5.30.3

# to be able to generate a similiar dependency tree at runtime
# we prepare a symlink tree for all dependencies
mkdir -p ${prefix}/deps
for dir in FLINT GMP MPFR PPL Perl bliss boost cddlib lrslib normaliz; do
   ln -s .. ${prefix}/deps/${dir}_jll
done

# adjust for hardcoded /workspace dirs
atomic_patch -p1 ../patches/relocatable.patch

# to unbreak ctrl+c in julia
atomic_patch -p1 ../patches/sigint.patch

if [[ $target == *darwin* ]]; then
  # we cannot run configure and instead provide config files
  mkdir -p build/Opt
  mkdir -p build/perlx/$perl_version/apple-darwin14
  cp ../config/config-$target.ninja build/config.ninja
  cp ../config/build-Opt-$target.ninja build/Opt/build.ninja
  cp ../config/targets.ninja build/targets.ninja
  ln -s ../config.ninja build/Opt/config.ninja
  cp ../config/perlx-config-$target.ninja build/perlx/$perl_version/apple-darwin14/config.ninja
  # for a modified pure perl JSON module and to make miniperl find the correct modules
  export PERL5LIB=$prefix/lib/perl5/$perl_version:$prefix/lib/perl5/$perl_version/darwin-2level:$WORKSPACE/srcdir/patches
  atomic_patch -p1 ../patches/polymake-cross.patch
  atomic_patch -p1 ../patches/polymake-cross-build.patch
else
  ./configure CFLAGS="-Wno-error" CC="$CC" CXX="$CXX" \
              PERL=${prefix}/deps/Perl_jll/bin/perl \
              LDFLAGS="$LDFLAGS -L${prefix}/deps/Perl_jll/lib -Wl,-rpath,${prefix}/deps/Perl_jll/lib" \
              --prefix=${prefix} \
              --with-flint=${prefix}/deps/FLINT_jll \
              --with-gmp=${prefix}/deps/GMP_jll \
              --with-mpfr=${prefix}/deps/MPFR_jll \
              --with-ppl=${prefix}/deps/PPL_jll \
              --with-bliss=${prefix}/deps/bliss_jll \
              --with-boost=${prefix}/deps/boost_jll \
              --with-cdd=${prefix}/deps/cddlib_jll \
              --with-lrs=${prefix}/deps/lrslib_jll \
              --with-libnormaliz=${prefix}/deps/normaliz_jll \
              --without-singular \
              --without-native
fi

ninja -v -C build/Opt -j${nproc}

ninja -v -C build/Opt install

# undo patch needed for building
if [[ $target == *darwin* ]]; then
  atomic_patch -R -p1 ../patches/polymake-cross-build.patch
fi
install -m 644 -D support/*.pl $prefix/share/polymake/support/

# replace miniperl
sed -i -e "s/miniperl-for-build/perl/" ${libdir}/polymake/config.ninja ${bindir}/polymake*
# replace binary path with env
sed -i -e "s#$bindir/perl#/usr/bin/env perl#g" ${libdir}/polymake/config.ninja ${bindir}/polymake*
# remove target and sysroot
sed -i -e "s#--sysroot[ =]\S\+##g" ${libdir}/polymake/config.ninja
sed -i -e "s#-target[ =]\S\+##g" ${libdir}/polymake/config.ninja

# copy of build config that has prefix as a variable
sed -e "s#${prefix}#\${prefix}#g" ${libdir}/polymake/config.ninja > ${libdir}/polymake/config-reloc.ninja

# cleanup symlink tree
rm -rf ${prefix}/deps

# copy julia script to generate dependency-tree at load time
cp ../patches/generate_deps_tree.jl $prefix/share/polymake

install_license COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpolymake", :libpolymake)
    LibraryProduct("libpolymake-apps-rt", :libpolymake_apps_rt)
    ExecutableProduct("polymake", :polymake)
    ExecutableProduct("polymake-config", Symbol("polymake_config"))
    FileProduct("share/polymake/generate_deps_tree.jl", :generate_deps_tree)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("MPFR_jll", v"4.0.2"),
    Dependency("FLINT_jll", compat = "~200.700"),
    Dependency("PPL_jll", compat = "~1.2"),
    Dependency("Perl_jll", compat = "=5.30.3"),
    Dependency("bliss_jll", compat = "~0.73"),
    Dependency("boost_jll", compat = "~1.71"),
    Dependency("cddlib_jll", compat = "~0.94"),
    Dependency("lrslib_jll", compat = "~0.1, ~0.2"),
    Dependency("normaliz_jll", compat = "~300.800.900"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: topolarity/Yggdrasil

max_issues_repo_name: topolarity/Yggdrasil

max_forks_repo_name: topolarity/Yggdrasil


text:
struct HorizontalHeader
    ascender::Int16
    descender::Int16
    line_gap::Int16
    advance_width_max::UInt16
    min_left_side_bearing::Int16
    min_right_side_bearing::Int16
    x_max_extent::Int16
    caret_slope_rise::Int16
    caret_slope_run::Int16
    caret_offset::Int16
    metric_data_format::Int16
    nhmetrics::UInt16
end

function Base.parse(io::IO, ::Type{HorizontalHeader})
    major = read(io, UInt16)
    minor = read(io, UInt16)
    @assert VersionNumber(major, minor) == v"1.0"
    hhea = HorizontalHeader(
        [read(io, T) for T in fieldtypes(HorizontalHeader)[1:10]]...,
        (skip(io, 8); read(io, Int16)),
        read(io, UInt16),
    )
    @assert hhea.metric_data_format == 0 "Metric data format should be zero, got $(hhea.metric_data_format)"
    hhea
end

struct VerticalHeader
    ascent::Int16
    descent::Int16
    line_gap::Int16
    advance_height_max::Int16
    min_top_side_bearing::Int16
    min_bottom_side_bearing::Int16
    y_max_extent::Int16
    caret_slope_rise::Int16
    caret_slope_run::Int16
    caret_offset::Int16
    metric_data_format::Int16
    nvmetrics::UInt16
end

function Base.parse(io::IO, ::Type{VerticalHeader})
    skip(io, 4)
    VerticalHeader(
        [read(io, T) for T in fieldtypes(VerticalHeader)[1:10]]...,
        (skip(io, 8); read(io, Int16)),
        read(io, UInt16),
    )
end

struct HorizontalMetric
    advance_width::UInt16
    left_side_bearing::Int16
end

struct HorizontalMetrics
    metrics::Vector{HorizontalMetric}
    left_side_bearings::Vector{Int16}
end

Base.read(io::IO, ::Type{HorizontalMetric}) = HorizontalMetric(read(io, UInt16), read(io, Int16))

function Base.parse(io::IO, ::Type{HorizontalMetrics}, hhea::HorizontalHeader, maxp::MaximumProfile)
    metrics = [read(io, HorizontalMetric) for _ in 1:hhea.nhmetrics]
    left_side_bearings = [read(io, Int16) for _ in 1:(maxp.nglyphs - hhea.nhmetrics)]
    HorizontalMetrics(metrics, left_side_bearings)
end

struct VerticalMetric
    advance_width::UInt16
    top_side_bearing::Int16
end

struct VerticalMetrics
    metrics::Vector{VerticalMetric}
    top_side_bearings::Vector{Int16}
end

Base.read(io::IO, ::Type{VerticalMetric}) = VerticalMetric(read(io, UInt16), read(io, Int16))

function Base.parse(io::IO, ::Type{VerticalMetrics}, vhea::VerticalHeader, maxp::MaximumProfile)
    metrics = [read(io, VerticalMetric) for _ in 1:vhea.nvmetrics]
    top_side_bearings = [read(io, Int16) for _ in 1:(maxp.nglyphs - vhea.nvmetrics)]
    VerticalMetrics(metrics, top_side_bearings)
end


text:
"""
    Structure(cell::Mat3, atoms::Vector{Atom})

The structure on which the [`Calculations`](@ref Calculation) will be performed.

    Structure(cif_file::String)

Creates a [`Structure`](@ref) from the supplied cif file.
"""
mutable struct Structure
    cell  :: Mat3{typeof(1.0Ang)}
    atoms :: Vector{Atom}
end
Structure() = Structure(Mat3(fill(1.0Ang, 3, 3)), Atom[])

function Structure(cif_file::String)
    str = cif2structure(cif_file)
    @info "Structure extracted from $cif_file\n\tcell parameters: \n\t a = $((str.cell[:,1]...,))\n\t b = $((str.cell[:,2]...,))\n\t c = $((str.cell[:,3]...,))\n\tnat = $(length(str.atoms))\n\telements = $(unique(getfield.(str.atoms, :name)))"
    return str
end

StructTypes.StructType(::Type{Structure}) = StructTypes.Struct()

function Base.:(==)(str1::Structure, str2::Structure)
    return str1.cell == str2.cell && str1.atoms == str2.atoms
end

"Uses cif2cell to Meta.parse a cif file, then returns the parsed structure."
function cif2structure(cif_file::String)
    tmpdir = dirname(cif_file)
    tmpfile = joinpath(tmpdir, "tmp.in")
    @assert splitext(cif_file)[2] == ".cif" error("Please specify a valid cif calculation file")
    if Sys.which("cif2cell") === nothing
        run(`$(DFControl.PYTHONPATH) $(DFControl.CIF2CELLPATH) $cif_file --no-reduce -p quantum-espresso -o $tmpfile`)
    else
        run(`cif2cell $cif_file --no-reduce -p quantum-espresso -o $tmpfile`)
    end
    bla, structure = DFC.FileIO.qe_read_calculation(tmpfile)
    rm(tmpfile)
    return structure
end

# TODO extend base.merge
"Takes a vector of structures and merges all the attributes of the atoms."
function mergestructures(structures::Vector{Structure})
    nonvoid = filter(x -> x != nothing, structures)
    out = nonvoid[1]
    default_at = Atom(; name = :Rh, position_cart = zero(Point3{typeof(1.0Ang)}),
                      position_cryst = zero(Point3{Float64}))
    for structure in nonvoid[2:end]
        for at1 in out.atoms, at2 in structure.atoms
            if at1 == at2
                for fname in fieldnames(typeof(at1))
                    if fname in [:name, :element, :position_cart, :position_cryst]
                        continue
                    end
                    field = getfield(at2, fname)
                    if field != getfield(default_at, fname)
                        setfield!(at1, fname, field)
                    end
                end
            end
        end
    end
    return out
end

"""
    update_geometry!(str1::Structure, str2::Structure)

Updates the spatial parameters of the `atoms` and `cell` of the first structure to those found in the second.
"""
function update_geometry!(str1::Structure, str2::Structure)
    str1.cell = copy(str2.cell)
    ats2 = str2.atoms
    for at1 in str1.atoms
        tats2 = filter(y -> y.name == at1.name, ats2)
        id = findmin(map(x -> norm(x.position_cart - at1.position_cart), tats2))[2]
        if id === nothing
            @error "No atom of the species $(at1.name) found in the second structure"
        end
        set_position!(at1, tats2[id].position_cryst, str1.cell)
    end
    return str1
end
Base.length(str::Structure) = length(str.atoms)
Base.iterate(str::Structure, args...) = iterate(str.atoms, args...)

### CELL ###
"""
    set_cell!(structure::Structure, c::Mat3)

Sets the `cell` of the `Structure` to `c`.
`c` will first be converted to the right units to match those
of the current `cell`.
All the absolute coordinates of the atoms will be recalculated
based on the new `cell`.
"""
function set_cell!(structure::Structure, c::Mat3)
    @warn "Converting cell to the units of the current cell of the structure."
    cell_ = uconvert.(Ang, c)
    for a in structure.atoms
        set_position!(a, a.position_cryst, cell_)
    end
    return structure.cell = cell_
end

"""
    scale_cell!(structure::Structure, scalemat::Matrix)

Rescales the cell of the `structure`.
"""
function scale_cell!(structure::Structure, scalemat::Matrix)
    new_cell = Mat3(structure.cell * scalemat)
    return set_cell!(structure, new_cell)
end

"""
    volume(cell::Mat3)
    volume(str::Structure)

Calculates the volume for the unit cell.
"""
volume(cell::Mat3) = det(cell)
volume(str::Structure) = volume(str.cell)
reciprocal(cell::AbstractMatrix) = inv(cell)

"""
    a(str::Structure)

First lattice vector.
"""
a(str::Structure) = str.cell[:, 1]
"""
    b(str::Structure)

Second lattice vector.
"""
b(str::Structure) = str.cell[:, 2]
"""
    c(str::Structure)

Third lattice vector.
"""
c(str::Structure) = str.cell[:, 3]

"""
    create_supercell(structure::Structure, na::Int, nb::Int, nc::Int; make_afm=false)
    create_supercell(structure::Structure, na::UnitRange, nb::UnitRange, nc::UnitRange; make_afm=false)

Takes a structure and creates a supercell from it with: the given amount of additional cells if (`na::Int, nb::Int, nc::Int`) along the a, b, c direction, or amount of cells specified by the ranges i.e. `-1:1, -1:1, -1:1` would create a 3x3x3 supercell.
If `make_afm` is set to `true` all the labels and magnetizations of the magnetic atoms will be reversed in a checkerboard fashion.
"""
function create_supercell(structure::Structure, na::UnitRange, nb::UnitRange, nc::UnitRange;
                          make_afm = false)
    orig_ats  = structure.atoms
    orig_cell = structure.cell
    scale_mat = diagm(0 => length.([na, nb, nc]))

    orig_uats = unique(orig_ats)

    new_cell  = orig_cell * scale_mat
    new_atoms = eltype(orig_ats)[]
    for ia in na, ib in nb, ic in nc
        transl_vec = orig_cell * [ia, ib, ic]
        factor = isodd(ia + ib + ic) ? -1 : 1
        for at in orig_ats
            cart_pos = at.position_cart + transl_vec
            cryst_pos = inv(new_cell) * cart_pos
            if make_afm && norm(at.magnetization) != 0
                new_magnetization = factor * at.magnetization
                push!(new_atoms,
                      Atom(at; projections = deepcopy(at.projections),
                           magnetization = new_magnetization, position_cart = cart_pos,
                           position_cryst = Point3(cryst_pos)))
            else
                push!(new_atoms,
                      Atom(at; projections = deepcopy(at.projections),
                           position_cart = cart_pos, position_cryst = Point3(cryst_pos)))
            end
        end
    end
    out = Structure(Mat3(new_cell), new_atoms)
    sanitize!(out)
    return out
end

function create_supercell(structure::Structure, na::Int, nb::Int, nc::Int; make_afm = false)
    return create_supercell(structure, 0:na, 0:nb, 0:nc; make_afm = make_afm)
end

### Atoms ###
"""
    getindex(structure::Structure, i::Int)
    getindex(structure::Structure, name::Symbol)
    getindex(structure::Structure, el::Element)

Returns the `i`th atom in `structure`, or all atoms with `name` or are of element `el`.
"""
Base.getindex(str::Structure, args...) = str.atoms[args...]
Base.getindex(str::Structure, el::Symbol) = filter(x -> x.name == el, str.atoms)
Base.getindex(str::Structure, el::Element) = filter(x -> x.element == el, str.atoms)
ismagnetic(str::Structure) = any(x -> norm(x.magnetization) > 0, str.atoms)

## Magnetism ##
function iscolin(str::Structure)
    magats = filter(x -> norm(x.magnetization) > 0, str.atoms)
    return !isempty(magats) &&
           all(x -> isapprox(abs(normalize(x.magnetization) ⋅ [0, 0, 1]), 1.0), magats)
end

function isnoncolin(str::Structure)
    return any(x -> norm(x.magnetization) > 0 &&
                   !isapprox(abs(normalize(x.magnetization) ⋅ [0, 0, 1]), 1.0), str.atoms)
end

function sanitize!(str::Structure)
    magnetic_ats = filter(a -> norm(a.magnetization) != 0 || a.dftu.U != 0.0, str.atoms)
    magnetic_elements = map(x -> x.element, magnetic_ats)
    for e in magnetic_elements
        magnetizations = Vec3[]
        dftus          = DFTU[]
        names          = Symbol[]
        for a in filter(x -> x.element == e, magnetic_ats)
            nameid = findfirst(x -> x == a.name, names)
            if nameid === nothing
                push!(names, a.name)
                push!(magnetizations, a.magnetization)
                push!(dftus, a.dftu)
            else
                mag = magnetizations[nameid]
                dftu_ = dftus[nameid]
                if (a.magnetization != mag || a.dftu != dftu_)
                    id = findfirst(x -> dftus[x] == a.dftu &&
                                       magnetizations[x] == a.magnetization,
                                   1:length(dftus))
                    if id === nothing
                        id = 1
                        tname = Symbol(string(e.symbol) * "$id")
                        while tname ∈ names
                            id += 1
                            tname = Symbol(string(e.symbol) * "$id")
                        end
                        push!(names, tname)
                        push!(magnetizations, a.magnetization)
                        push!(dftus, a.dftu)
                        oldname = a.name
                        a.name = tname
                        @info "Renamed atom from $oldname to $(a.name) in order to distinguish different magnetization species."
                    else
                        oldname = a.name
                        a.name = names[id]
                        @info "Renamed atom from $oldname to $(a.name) in order to distinguish different magnetization species."
                    end
                end
            end
        end
    end
end

"""
    polyhedron(at::Atom, atoms::Vector{Atom}, order::Int)
    polyhedron(at::Atom, str::Structure, order::Int)

Returns a polyhedron around the atom, i.e. the `order` closest atoms.
The returned atoms will be ordered according to their distance to the first one.
In the case of a structure rather than a set of atoms, the search will
be performed over all atoms in the structure.
"""
function polyhedron(at::Atom, atoms::Vector{Atom}, order::Int)
    return sort(atoms; by = x -> norm(x.position_cart - at.position_cart))[1:order]
end
function polyhedron(at::Atom, str::Structure, order::Int)
    return polyhedron(at, create_supercell(str, -1:1, -1:1, -1:1).atoms, order)
end

function set_pseudos!(structure::Structure, pseudos)
    for at in structure.atoms
        pseudo = get(pseudos, at.element.symbol, nothing)
        if pseudo === nothing
            @warn "Pseudo for $(at.name) not found."
        else
            at.pseudo = pseudo
        end
    end
end

### SPGLIB ###
const DEFAULT_TOLERANCE = 1e-5

struct SPGStructure
    lattice::Matrix{Cdouble}
    positions::Matrix{Cdouble}
    species_indices::Vector{Cint}
end

function SPGStructure(s::Structure)
    clattice = convert(Matrix{Cdouble}, ustrip.(s.cell'))
    cpositions = convert(Matrix{Cdouble}, hcat(map(x -> x.position_cryst, s.atoms)...))
    uats = unique(s.atoms)
    species_indices = Cint[findfirst(x -> isequal_species(x, at), uats) for at in s.atoms]
    return SPGStructure(clattice, cpositions, species_indices)
end

"""
    symmetry_operators(s::Structure; maxsize=52, tolerance=$DEFAULT_TOLERANCE)

Finds and returns all the rotations and translations that are symmetry operators of the structure.
"""
symmetry_operators(s::Structure; kwargs...) = symmetry_operators(SPGStructure(s); kwargs...)
function symmetry_operators(s::SPGStructure; maxsize = 52, tolerance = DEFAULT_TOLERANCE)
    rotations    = Array{Cint}(undef, 3, 3, maxsize)
    translations = Array{Cdouble}(undef, 3, maxsize)

    num_ops = ccall((:spg_get_symmetry, SPGLIB), Cint,
                    (Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
                     Cint, Cdouble), rotations, translations, maxsize, s.lattice,
                    s.positions, s.species_indices, length(s.species_indices), tolerance)
    return (rotations = [Mat3{Int}(rotations[:, :, i]) for i in 1:num_ops],
            translations = [Vec3(translations[:, i]) for i in 1:num_ops])
end

"""
    international(s::Structure; tolerance=$DEFAULT_TOLERANCE)

Returns the international symbol of the space group of the structure.
"""
international(s::Structure; kwargs...) = international(SPGStructure(s); kwargs...)
function international(s::SPGStructure; tolerance = DEFAULT_TOLERANCE)
    res = zeros(Cchar, 11)

    num = ccall((:spg_get_international, SPGLIB), Cint,
                (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble), res,
                s.lattice, s.positions, s.species_indices, length(s.species_indices),
                tolerance)
    num == 0 && error("Could not determine the international symbol.")

    return num, join(convert(Vector{Char}, res[1:findfirst(iszero, res)-1]))
end

"""
    niggli_reduce(s::Structure; tolerance=$DEFAULT_TOLERANCE)

Returns the niggli reduced `Structure`.
"""
function niggli_reduce(s::Structure; tolerance = DEFAULT_TOLERANCE)
    reduced_spg_structure = niggli_reduce!(SPGStructure(s); tolerance = DEFAULT_TOLERANCE)
    uats = unique(s.atoms)
    c = Mat3{Float64}(reduced_spg_structure.lattice') .* 1Ang
    return Structure(c, map(1:length(reduced_spg_structure.species_indices)) do i
                         at = deepcopy(uats[reduced_spg_structure.species_indices[i]])
                         pos = reduced_spg_structure.positions[:, i]
                         set_position!(at, pos, c)
                         return at
                     end)
end

function niggli_reduce(c::Mat3{Float64}; tolerance = DEFAULT_TOLERANCE)
    ccell = convert(Matrix{Cdouble}, c)
    numops = ccall((:spg_niggli_reduce, SPGLIB), Cint, (Ptr{Cdouble}, Cdouble), ccell,
                   tolerance)
    return Mat3{Float64}(ccell)
end

function niggli_reduce!(s::SPGStructure; tolerance = DEFAULT_TOLERANCE)
    numops = ccall((:spg_niggli_reduce, SPGLIB), Cint, (Ptr{Cdouble}, Cdouble), s.lattice,
                   tolerance)
    numops == 0 && error("Could not determine the niggli reduced cell.")

    return s
end

function find_primitive!(s::SPGStructure; tolerance = DEFAULT_TOLERANCE)
    numats = ccall((:spg_find_primitive, SPGLIB), Cint,
                   (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint, Cdouble), s.lattice,
                   s.positions, s.species_indices, length(s.species_indices), tolerance)
    numats == 0 && error("Could not find the primitive of the supplied structure.")
    return SPGStructure(s.lattice, s.positions[:, 1:numats], s.species_indices[1:numats])
end

function find_primitive(s::Structure; kwargs...)
    uats = unique(s.atoms)
    spg = find_primitive!(SPGStructure(s))
    new_cell = Mat3{Float64}(spg.lattice) * unit(eltype(s.cell))
    newats = eltype(uats)[]
    for i in 1:length(spg.species_indices)
        tat = deepcopy(uats[spg.species_indices[i]])
        set_position!(tat, Point3(spg.positions[:, i]), new_cell)
        push!(newats, tat)
    end

    return Structure(new_cell, newats)
end

"""
    cell_parameters(cell::Mat3)
    cell_parameters(str::Structure)

Parameters `(a, b, c, α, β, γ)`of the calculation cell returned in a `NamedTuple`.
"""
cell_parameters(s::Structure) = cell_parameters(ustrip.(s.cell))
function cell_parameters(cell::Mat3)
    G = transpose(cell) * cell
    a, b, c = sqrt.(diag(G))
    α = acos(0.5(G[2, 3] + G[3, 2]) / (c * b))
    β = acos(0.5(G[3, 1] + G[1, 3]) / (a * c))
    γ = acos(0.5(G[1, 2] + G[2, 1]) / (a * b))
    return a, b, c, α, β, γ
end

function crystal_kind(s::Structure; tolerance = DEFAULT_TOLERANCE)
    n, sym = international(s)
    f = (i, j) -> i <= n <= j
    cs = [:triclinic => (1, 2), :monoclinic => (3, 15), :orthorhombic => (16, 74),
          :tetragonal => (75, 142), :trigonal => (143, 167), :hexagonal => (168, 194),
          :cubic => (195, 230)]
    for (k, v) in cs
        if f(v...)
            return k
        end
    end
    return error("Crystal kind could not be determined")
end

function lattice_kind(s::Structure; tolerance = DEFAULT_TOLERANCE)
    n, sym = international(s)
    kind = crystal_kind(s; tolerance = tolerance)
    if n ∈ (146, 148, 155, 160, 161, 166, 167)
        return :rhombohedral
    elseif kind == :trigonal
        return :hexagonal
    else
        return kind
    end
end

"""
    high_symmetry_kpath(s::Structure, npoints_per_segment::Int; package=QE, kwargs...)

Generates a QE bands calculation compliant high symmetry kpath, to be used with
e.g. `set_kpoints!(bands_calculation, kpoints)`.
"""
function high_symmetry_kpath(s::Structure, npoints_per_segment::Int; package = QE,
                             kwargs...)
    kpoints, kpath = high_symmetry_kpoints(s)
    out_k = NTuple{4,Float64}[]
    if package === QE
        for path in kpath
            for p in path
                push!(out_k, (kpoints[p]..., float(npoints_per_segment)))
            end
        end
    else
        error("high_symmetry_kpath not implemented for package $package.")
    end
    out_k[end] = (out_k[end][1:3]..., 1.0)
    return out_k
end

"""
    high_symmetry_kpoints(s::Structure; tolerance = 1e-5)

Returns `(kpoints, path)` where `kpoints` are the high-symmetry k-points,
and `path` are the sections of the high symmetry path through the first Brillouin Zone.
"""
function high_symmetry_kpoints(s::Structure; tolerance = DEFAULT_TOLERANCE)
    n, sym    = international(s)
    kind      = lattice_kind(s)
    primitive = find_primitive(s)

    primcell = ustrip.(primitive.cell)

    a, b, c, α, β, γ = cell_parameters(s)

    symerror = () -> error("Unexpected value for international symbol: $sym")
    if kind == :cubic
        if occursin("P", sym)
            return cubic()
        elseif occursin("F", sym)
            return fcc()
        elseif occursin("I", sym)
            return bcc()
        else
            symerror()
        end
    elseif kind == :tetragonal
        if occursin("P", sym)
            return tet()
        elseif occursin("I", sym)
            if c < a
                return bctet1(c, a)
            else
                return bctet2(c, a)
            end
        else
            symerror()
        end
    elseif kind == :orthorhombic
        if occursin("P", sym)
            return orc()
        elseif occursin("F", sym)
            if 1 / a^2 > 1 / b^2 + 1 / c^2
                return orcf1(a, b, c)
            elseif 1 / a^2 < 1 / b^2 + 1 / c^2
                return orcf2(a, b, c)
            else
                return orcf3(a, b, c)
            end

        elseif occursin("I", sym)
            return orci(a, b, c)

        elseif occursin("C", sym) || occursin("A", sym)
            return orcc(a, b, c)
        else
            symerror()
        end

    elseif kind == :hexagonal
        return hex()

    elseif kind == :rhombohedral
        if α < π / 2
            return rhl1(α)
        else
            return rhl2(α)
        end

    elseif kind == :monoclinic
        prim_rec = inv(primcell)
        ka, kb, kc, kα, kβ, kγ = cell_parameters(prim_rec)
        if occursin("P", sym)
            return mcl(b, c, α)

        elseif occursin("C", sym)
            if kγ > π / 2
                return mclc1(a, b, c, α)
            elseif kγ == π / 2
                return mclc2(a, b, c, α)
            elseif kγ < π / 2
                if b * cos(α) / c + b^2 * sin(α)^2 / a^2 < 1
                    return mclc3(a, b, c, α)
                elseif b * cos(α) / c + b^2 * sin(α)^2 / a^2 == 1
                    return mclc4(a, b, c, α)
                elseif b * cos(α) / c + b^2 * sin(α)^2 / a^2 > 1
                    return mclc5(a, b, c, α)
                end
            end
        else
            symerror()
        end

    elseif kind == :triclinic
        prim_rec = inv(primcell)
        ka, kb, kc, kα, kβ, kγ = cell_parameters(prim_rec)
        if kα > π / 2 && kβ > π / 2 && kγ > π / 2
            return tria()
        elseif kα < π / 2 && kβ < π / 2 && kγ < π / 2
            return trib()
        elseif kα > π / 2 && kβ > π / 2 && kγ == π / 2
            return tria()
        elseif kα < π / 2 && kβ < π / 2 && kγ == π / 2
            return trib()
        end

    else
        error("Unknown lattice type $kind")
    end
end

function cubic()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :X => Vec3([0.0, 0.5, 0.0]),
                   :R => Vec3([0.5, 0.5, 0.5]), :M => Vec3([0.5, 0.5, 0.0]))
    path = [[:Γ, :X, :M, :Γ, :R, :X], [:M, :R]]
    return (kpoints = kpoints, path = path)
end

function fcc()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]),
                   :K => Vec3([3.0 / 8.0, 3.0 / 8.0, 3.0 / 4.0]),
                   :L => Vec3([0.5, 0.5, 0.5]),
                   :U => Vec3([5.0 / 8.0, 1.0 / 4.0, 5.0 / 8.0]),
                   :W => Vec3([0.5, 1.0 / 4.0, 3.0 / 4.0]), :X => Vec3([0.5, 0.0, 0.5]))

    path = [[:Γ, :X, :W, :K, :Γ, :L, :U, :W, :L, :K], [:U, :X]]
    return (kpoints = kpoints, path = path)
end

function bcc()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :H => Vec3([0.5, -0.5, 0.5]),
                   :P => Vec3([0.25, 0.25, 0.25]), :N => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :H, :N, :Γ, :P, :H], [:P, :N]]
    return (kpoints = kpoints, path = path)
end

function tet()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([0.5, 0.5, 0.5]),
                   :M => Vec3([0.5, 0.5, 0.0]), :R => Vec3([0.0, 0.5, 0.5]),
                   :X => Vec3([0.0, 0.5, 0.0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :X, :M, :Γ, :Z, :R, :A, :Z], [:X, :R], [:M, :A]]

    return (kpoints = kpoints, path = path)
end

function bctet1(c, a)
    eta = (1 + c^2 / a^2) / 4.0
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :M => Vec3([-0.5, 0.5, 0.5]),
                   :N => Vec3([0.0, 0.5, 0.0]), :P => Vec3([0.25, 0.25, 0.25]),
                   :X => Vec3([0.0, 0.0, 0.5]), :Z => Vec3([eta, eta, -eta]),
                   :Z_1 => Vec3([-eta, 1 - eta, eta]))
    path = [[:Γ, :X, :M, :Γ, :Z, :P, :N, :Z_1, :M], [:X, :P]]
    return (kpoints = kpoints, path = path)
end

function bctet2(c, a)
    eta = (1 + a^2 / c^2) / 4.0
    zeta = a^2 / (2 * c^2)
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :N => Vec3([0.0, 0.5, 0.0]),
                   :P => Vec3([0.25, 0.25, 0.25]), :Sigma => Vec3([-eta, eta, eta]),
                   :Sigma_1 => Vec3([eta, 1 - eta, -eta]), :X => Vec3([0.0, 0.0, 0.5]),
                   :Y => Vec3([-zeta, zeta, 0.5]), :Y_1 => Vec3([0.5, 0.5, -zeta]),
                   :Z => Vec3([0.5, 0.5, -0.5]))
    path = [[:Γ, :X, :Y, :Sigma, :Γ, :Z, :Sigma_1, :N, :P, :Y_1, :Z], [:X, :P]]
    return (kpoints = kpoints, path = path)
end

function orc()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :R => Vec3([0.5, 0.5, 0.5]),
                   :S => Vec3([0.5, 0.5, 0.0]), :T => Vec3([0.0, 0.5, 0.5]),
                   :U => Vec3([0.5, 0.0, 0.5]), :X => Vec3([0.5, 0.0, 0.0]),
                   :Y => Vec3([0.0, 0.5, 0.0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :X, :S, :Y, :Γ, :Z, :U, :R, :T, :Z], [:Y, :T], [:U, :X], [:S, :R]]
    return (kpoints = kpoints, path = path)
end

function orcf1(a, b, c)
    zeta = (1 + a^2 / b^2 - a^2 / c^2) / 4
    eta = (1 + a^2 / b^2 + a^2 / c^2) / 4

    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([0.5, 0.5 + zeta, zeta]),
                   :A_1 => Vec3([0.5, 0.5 - zeta, 1 - zeta]), :L => Vec3([0.5, 0.5, 0.5]),
                   :T => Vec3([1, 0.5, 0.5]), :X => Vec3([0.0, eta, eta]),
                   :X_1 => Vec3([1, 1 - eta, 1 - eta]), :Y => Vec3([0.5, 0.0, 0.5]),
                   :Z => Vec3([0.5, 0.5, 0.0]))
    path = [[:Γ, :Y, :T, :Z, :Γ, :X, :A_1, :Y], [:T, :X_1], [:X, :A, :Z], [:L, :Γ]]
    return (kpoints = kpoints, path = path)
end

function orcf2(a, b, c)
    phi = (1 + c^2 / b^2 - c^2 / a^2) / 4
    eta = (1 + a^2 / b^2 - a^2 / c^2) / 4
    delta = (1 + b^2 / a^2 - b^2 / c^2) / 4
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :C => Vec3([0.5, 0.5 - eta, 1 - eta]),
                   :C_1 => Vec3([0.5, 0.5 + eta, eta]),
                   :D => Vec3([0.5 - delta, 0.5, 1 - delta]),
                   :D_1 => Vec3([0.5 + delta, 0.5, delta]), :L => Vec3([0.5, 0.5, 0.5]),
                   :H => Vec3([1 - phi, 0.5 - phi, 0.5]),
                   :H_1 => Vec3([phi, 0.5 + phi, 0.5]), :X => Vec3([0.0, 0.5, 0.5]),
                   :Y => Vec3([0.5, 0.0, 0.5]), :Z => Vec3([0.5, 0.5, 0.0]))
    path = [[:Γ, :Y, :C, :D, :X, :Γ, :Z, :D_1, :H, :C], [:C_1, :Z], [:X, :H_1], [:H, :Y],
            [:L, :Γ]]
    return (kpoints = kpoints, path = path)
end

function orcf3(a, b, c)
    zeta = (1 + a^2 / b^2 - a^2 / c^2) / 4
    eta = (1 + a^2 / b^2 + a^2 / c^2) / 4
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([0.5, 0.5 + zeta, zeta]),
                   :A_1 => Vec3([0.5, 0.5 - zeta, 1 - zeta]), :L => Vec3([0.5, 0.5, 0.5]),
                   :T => Vec3([1, 0.5, 0.5]), :X => Vec3([0.0, eta, eta]),
                   :X_1 => Vec3([1, 1 - eta, 1 - eta]), :Y => Vec3([0.5, 0.0, 0.5]),
                   :Z => Vec3([0.5, 0.5, 0.0]))
    path = [[:Γ, :Y, :T, :Z, :Γ, :X, :A_1, :Y], [:X, :A, :Z], [:L, :Γ]]
    return (kpoints = kpoints, path = path)
end

function orci(a, b, c)
    zeta = (1 + a^2 / c^2) / 4
    eta = (1 + b^2 / c^2) / 4
    delta = (b^2 - a^2) / (4 * c^2)
    mu = (a^2 + b^2) / (4 * c^2)
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :L => Vec3([-mu, mu, 0.5 - delta]),
                   :L_1 => Vec3([mu, -mu, 0.5 + delta]),
                   :L_2 => Vec3([0.5 - delta, 0.5 + delta, -mu]),
                   :R => Vec3([0.0, 0.5, 0.0]), :S => Vec3([0.5, 0.0, 0.0]),
                   :T => Vec3([0.0, 0.0, 0.5]), :W => Vec3([0.25, 0.25, 0.25]),
                   :X => Vec3([-zeta, zeta, zeta]), :X_1 => Vec3([zeta, 1 - zeta, -zeta]),
                   :Y => Vec3([eta, -eta, eta]), :Y_1 => Vec3([1 - eta, eta, -eta]),
                   :Z => Vec3([0.5, 0.5, -0.5]))
    path = [[:Γ, :X, :L, :T, :W, :R, :X_1, :Z, :Γ, :Y, :S, :W], [:L_1, :Y], [:Y_1, :Z]]
    return (kpoints = kpoints, path = path)
end

function orcc(a, b, c)
    zeta = (1 + a^2 / b^2) / 4
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([zeta, zeta, 0.5]),
                   :A_1 => Vec3([-zeta, 1 - zeta, 0.5]), :R => Vec3([0.0, 0.5, 0.5]),
                   :S => Vec3([0.0, 0.5, 0.0]), :T => Vec3([-0.5, 0.5, 0.5]),
                   :X => Vec3([zeta, zeta, 0.0]), :X_1 => Vec3([-zeta, 1 - zeta, 0.0]),
                   :Y => Vec3([-0.5, 0.5, 0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :X, :S, :R, :A, :Z, :Γ, :Y, :X_1, :A_1, :T, :Y], [:Z, :T]]
    return (kpoints = kpoints, path = path)
end

function hex()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([0.0, 0.0, 0.5]),
                   :H => Vec3([1.0 / 3.0, 1.0 / 3.0, 0.5]),
                   :K => Vec3([1.0 / 3.0, 1.0 / 3.0, 0.0]), :L => Vec3([0.5, 0.0, 0.5]),
                   :M => Vec3([0.5, 0.0, 0.0]))
    path = [[:Γ, :M, :K, :Γ, :A, :L, :H, :A], [:L, :M], [:K, :H]]
    return (kpoints = kpoints, path = path)
end

function rhl1(alpha)
    eta = (1 + 4 * cos(alpha)) / (2 + 4 * cos(alpha))
    nu = 3.0 / 4.0 - eta / 2.0
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :B => Vec3([eta, 0.5, 1.0 - eta]),
                   :B_1 => Vec3([1.0 / 2.0, 1.0 - eta, eta - 1.0]),
                   :F => Vec3([0.5, 0.5, 0.0]), :L => Vec3([0.5, 0.0, 0.0]),
                   :L_1 => Vec3([0.0, 0.0, -0.5]), :P => Vec3([eta, nu, nu]),
                   :P_1 => Vec3([1.0 - nu, 1.0 - nu, 1.0 - eta]),
                   :P_2 => Vec3([nu, nu, eta - 1.0]), :Q => Vec3([1.0 - nu, nu, 0.0]),
                   :X => Vec3([nu, 0.0, -nu]), :Z => Vec3([0.5, 0.5, 0.5]))
    path = [[:Γ, :L, :B_1], [:B, :Z, :Γ, :X], [:Q, :F, :P_1, :Z], [:L, :P]]
    return (kpoints = kpoints, path = path)
end

function rhl2(alpha)
    eta = 1 / (2 * tan(alpha / 2.0)^2)
    nu = 3.0 / 4.0 - eta / 2.0
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :F => Vec3([0.5, -0.5, 0.0]),
                   :L => Vec3([0.5, 0.0, 0.0]), :P => Vec3([1 - nu, -nu, 1 - nu]),
                   :P_1 => Vec3([nu, nu - 1.0, nu - 1.0]), :Q => Vec3([eta, eta, eta]),
                   :Q_1 => Vec3([1.0 - eta, -eta, -eta]), :Z => Vec3([0.5, -0.5, 0.5]))
    path = [[:Γ, :P, :Z, :Q, :Γ, :F, :P_1, :Q_1, :L, :Z]]
    return (kpoints = kpoints, path = path)
end

function mcl(b, c, beta)
    eta = (1 - b * cos(beta) / c) / (2 * sin(beta)^2)
    nu = 0.5 - eta * c * cos(beta) / b
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :A => Vec3([0.5, 0.5, 0.0]),
                   :C => Vec3([0.0, 0.5, 0.5]), :D => Vec3([0.5, 0.0, 0.5]),
                   :D_1 => Vec3([0.5, 0.5, -0.5]), :E => Vec3([0.5, 0.5, 0.5]),
                   :H => Vec3([0.0, eta, 1.0 - nu]), :H_1 => Vec3([0.0, 1.0 - eta, nu]),
                   :H_2 => Vec3([0.0, eta, -nu]), :M => Vec3([0.5, eta, 1.0 - nu]),
                   :M_1 => Vec3([0.5, 1 - eta, nu]), :M_2 => Vec3([0.5, 1 - eta, nu]),
                   :X => Vec3([0.0, 0.5, 0.0]), :Y => Vec3([0.0, 0.0, 0.5]),
                   :Y_1 => Vec3([0.0, 0.0, -0.5]), :Z => Vec3([0.5, 0.0, 0.0]))
    path = [[:Γ, :Y, :H, :C, :E, :M_1, :A, :X, :H_1], [:M, :D, :Z], [:Y, :D]]
    return (kpoints = kpoints, path = path)
end

function mclc1(a, b, c, alpha)
    zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha)^2)
    eta = 0.5 + 2 * zeta * c * cos(alpha) / b
    psi = 0.75 - a^2 / (4 * b^2 * sin(alpha)^2)
    phi = psi + (0.75 - psi) * b * cos(alpha) / c
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :N => Vec3([0.5, 0.0, 0.0]),
                   :N_1 => Vec3([0.0, -0.5, 0.0]),
                   :F => Vec3([1 - zeta, 1 - zeta, 1 - eta]),
                   :F_1 => Vec3([zeta, zeta, eta]), :F_2 => Vec3([-zeta, -zeta, 1 - eta]),
                   # :F_3 => Vec3([1 - zeta, -zeta, 1 - eta]),
                   :I => Vec3([phi, 1 - phi, 0.5]), :I_1 => Vec3([1 - phi, phi - 1, 0.5]),
                   :L => Vec3([0.5, 0.5, 0.5]), :M => Vec3([0.5, 0.0, 0.5]),
                   :X => Vec3([1 - psi, psi - 1, 0.0]), :X_1 => Vec3([psi, 1 - psi, 0.0]),
                   :X_2 => Vec3([psi - 1, -psi, 0.0]), :Y => Vec3([0.5, 0.5, 0.0]),
                   :Y_1 => Vec3([-0.5, -0.5, 0.0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :Y, :F, :L, :I], [:I_1, :Z, :F_1], [:Y, :X_1], [:X, :Γ, :N], [:M, :Γ]]
    return (kpoints = kpoints, path = path)
end

function mclc2(a, b, c, alpha)
    zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha)^2)
    eta = 0.5 + 2 * zeta * c * cos(alpha) / b
    psi = 0.75 - a^2 / (4 * b^2 * sin(alpha)^2)
    phi = psi + (0.75 - psi) * b * cos(alpha) / c
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :N => Vec3([0.5, 0.0, 0.0]),
                   :N_1 => Vec3([0.0, -0.5, 0.0]),
                   :F => Vec3([1 - zeta, 1 - zeta, 1 - eta]),
                   :F_1 => Vec3([zeta, zeta, eta]), :F_2 => Vec3([-zeta, -zeta, 1 - eta]),
                   :F_3 => Vec3([1 - zeta, -zeta, 1 - eta]),
                   :I => Vec3([phi, 1 - phi, 0.5]), :I_1 => Vec3([1 - phi, phi - 1, 0.5]),
                   :L => Vec3([0.5, 0.5, 0.5]), :M => Vec3([0.5, 0.0, 0.5]),
                   :X => Vec3([1 - psi, psi - 1, 0.0]), :X_1 => Vec3([psi, 1 - psi, 0.0]),
                   :X_2 => Vec3([psi - 1, -psi, 0.0]), :Y => Vec3([0.5, 0.5, 0.0]),
                   :Y_1 => Vec3([-0.5, -0.5, 0.0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :Y, :F, :L, :I], [:I_1, :Z, :F_1], [:N, :Γ, :M]]
    return (kpoints = kpoints, path = path)
end

function mclc3(a, b, c, alpha)
    mu = (1 + b^2 / a^2) / 4.0
    delta = b * c * cos(alpha) / (2 * a^2)
    zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha)^2)
    eta = 0.5 + 2 * zeta * c * cos(alpha) / b
    phi = 1 + zeta - 2 * mu
    psi = eta - 2 * delta
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :F => Vec3([1 - phi, 1 - phi, 1 - psi]),
                   :F_1 => Vec3([phi, phi - 1, psi]),
                   :F_2 => Vec3([1 - phi, -phi, 1 - psi]), :H => Vec3([zeta, zeta, eta]),
                   :H_1 => Vec3([1 - zeta, -zeta, 1 - eta]),
                   :H_2 => Vec3([-zeta, -zeta, 1 - eta]), :I => Vec3([0.5, -0.5, 0.5]),
                   :M => Vec3([0.5, 0.0, 0.5]), :N => Vec3([0.5, 0.0, 0.0]),
                   :N_1 => Vec3([0.0, -0.5, 0.0]), :X => Vec3([0.5, -0.5, 0.0]),
                   :Y => Vec3([mu, mu, delta]), :Y_1 => Vec3([1 - mu, -mu, -delta]),
                   :Y_2 => Vec3([-mu, -mu, -delta]), :Y_3 => Vec3([mu, mu - 1, delta]),
                   :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :Y, :F, :H, :Z, :I, :F_1], [:H_1, :Y_1, :X, :Γ, :N], [:M, :Γ]]
    return (kpoints = kpoints, path = path)
end

function mclc4(a, b, c, alpha)
    mu = (1 + b^2 / a^2) / 4.0
    delta = b * c * cos(alpha) / (2 * a^2)
    zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha)^2)
    eta = 0.5 + 2 * zeta * c * cos(alpha) / b
    phi = 1 + zeta - 2 * mu
    psi = eta - 2 * delta
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :F => Vec3([1 - phi, 1 - phi, 1 - psi]),
                   :F_1 => Vec3([phi, phi - 1, psi]),
                   :F_2 => Vec3([1 - phi, -phi, 1 - psi]), :H => Vec3([zeta, zeta, eta]),
                   :H_1 => Vec3([1 - zeta, -zeta, 1 - eta]),
                   :H_2 => Vec3([-zeta, -zeta, 1 - eta]), :I => Vec3([0.5, -0.5, 0.5]),
                   :M => Vec3([0.5, 0.0, 0.5]), :N => Vec3([0.5, 0.0, 0.0]),
                   :N_1 => Vec3([0.0, -0.5, 0.0]), :X => Vec3([0.5, -0.5, 0.0]),
                   :Y => Vec3([mu, mu, delta]), :Y_1 => Vec3([1 - mu, -mu, -delta]),
                   :Y_2 => Vec3([-mu, -mu, -delta]), :Y_3 => Vec3([mu, mu - 1, delta]),
                   :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :Y, :F, :H, :Z, :I], [:H_1, :Y_1, :X, :Γ, :N], [:M, :Γ]]
    return (kpoints = kpoints, path = path)
end

function mclc5(a, b, c, alpha)
    zeta = (b^2 / a^2 + (1 - b * cos(alpha) / c) / sin(alpha)^2) / 4
    eta = 0.5 + 2 * zeta * c * cos(alpha) / b
    mu = eta / 2 + b^2 / (4 * a^2) - b * c * cos(alpha) / (2 * a^2)
    nu = 2 * mu - zeta
    rho = 1 - zeta * a^2 / b^2
    omega = (4 * nu - 1 - b^2 * sin(alpha)^2 / a^2) * c / (2 * b * cos(alpha))
    delta = zeta * c * cos(alpha) / b + omega / 2 - 0.25
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :F => Vec3([nu, nu, omega]),
                   :F_1 => Vec3([1 - nu, 1 - nu, 1 - omega]),
                   :F_2 => Vec3([nu, nu - 1, omega]), :H => Vec3([zeta, zeta, eta]),
                   :H_1 => Vec3([1 - zeta, -zeta, 1 - eta]),
                   :H_2 => Vec3([-zeta, -zeta, 1 - eta]), :I => Vec3([rho, 1 - rho, 0.5]),
                   :I_1 => Vec3([1 - rho, rho - 1, 0.5]), :L => Vec3([0.5, 0.5, 0.5]),
                   :M => Vec3([0.5, 0.0, 0.5]), :N => Vec3([0.5, 0.0, 0.0]),
                   :N_1 => Vec3([0.0, -0.5, 0.0]), :X => Vec3([0.5, -0.5, 0.0]),
                   :Y => Vec3([mu, mu, delta]), :Y_1 => Vec3([1 - mu, -mu, -delta]),
                   :Y_2 => Vec3([-mu, -mu, -delta]), :Y_3 => Vec3([mu, mu - 1, delta]),
                   :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:Γ, :Y, :F, :L, :I], [:I_1, :Z, :H, :F_1], [:H_1, :Y_1, :X, :Γ, :N], [:M, :Γ]]
    return (kpoints = kpoints, path = path)
end

function tria()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :L => Vec3([0.5, 0.5, 0.0]),
                   :M => Vec3([0.0, 0.5, 0.5]), :N => Vec3([0.5, 0.0, 0.5]),
                   :R => Vec3([0.5, 0.5, 0.5]), :X => Vec3([0.5, 0.0, 0.0]),
                   :Y => Vec3([0.0, 0.5, 0.0]), :Z => Vec3([0.0, 0.0, 0.5]))
    path = [[:X, :Γ, :Y], [:L, :Γ, :Z], [:N, :Γ, :M], [:R, :Γ]]
    return (kpoints = kpoints, path = path)
end

function trib()
    kpoints = Dict(:Γ => Vec3([0.0, 0.0, 0.0]), :L => Vec3([0.5, -0.5, 0.0]),
                   :M => Vec3([0.0, 0.0, 0.5]), :N => Vec3([-0.5, -0.5, 0.5]),
                   :R => Vec3([0.0, -0.5, 0.5]), :X => Vec3([0.0, -0.5, 0.0]),
                   :Y => Vec3([0.5, 0.0, 0.0]), :Z => Vec3([-0.5, 0.0, 0.5]))
    path = [[:X, :Γ, :Y], [:L, :Γ, :Z], [:N, :Γ, :M], [:R, :Γ]]
    return (kpoints = kpoints, path = path)
end