text:
##################################################
#### Check if all possible solutions are found  ##
##################################################

@testset "all solutions and set" begin
    lhs(x) = sum(x) >= 1
    rhs(x) = sum(x) <= 2

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {sum(x) >= 1 && sum(x) <= 2})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test (lhs(xval) && rhs(xval))
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        if lhs([i,j]) && rhs([i,j])
            @test (i,j) in results
        end
    end
end

@testset "all solutions or set" begin
    lhs(x) = sum(x) >= 1
    rhs(x) = sum(x) <= 2

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {sum(x) >= 1 || sum(x) <= 2})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test (lhs(xval) || rhs(xval))
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        if lhs([i,j]) || rhs([i,j])
            @test (i,j) in results
        end
    end
end

@testset "all solutions xor set" begin
    lhs(x) = sum(x) >= 1
    rhs(x) = sum(x) <= 2

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {(sum(x) >= 1) ⊻ (sum(x) <= 2)})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test (lhs(xval) ⊻ rhs(xval))
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        if lhs([i,j]) ⊻ rhs([i,j])
            @test (i,j) in results
        end
    end
end

@testset "all solutions xnor set" begin
    lhs(x) = sum(x) >= 1
    rhs(x) = sum(x) <= 2

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {!((sum(x) >= 1) ⊻ (sum(x) <= 2))})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test !((lhs(xval) ⊻ rhs(xval)))
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        if !(lhs([i,j]) ⊻ rhs([i,j]))
            @test (i,j) in results
        end
    end
end

@testset "all solutions xnor set 2" begin
    lhs(x) = !(sum(x) > 3)
    rhs(x) = sum(x) >= 2

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {!((!(sum(x) > 3)) ⊻ (sum(x) >= 2))})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test !(lhs(xval) ⊻ rhs(xval))
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        if !(lhs([i,j]) ⊻ rhs([i,j]))
            @test (i,j) in results
        end
    end
end

@testset "all solutions variable as constraint" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, a, Bin)
    @variable(m, b, Bin)
    @variable(m, c, Bin)
    @constraint(m, a || b || !c)
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Bool,Bool,Bool}}()
    for i in 1:nresults
        a_val, b_val, c_val = convert.(Bool, round.(JuMP.value.([a,b,c]; result=i)))
        @test a_val || b_val || !c_val
        push!(results, (a_val, b_val, c_val))
    end
    for i in [true, false], j in [true, false], k in [true, false]
        if i || j || !k
            @test (i,j, k) in results
        end
    end
end

#=
@testset "all solutions combined" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_optimal_solutions"=>true))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, b, Bin)
    @constraint(m, b >= 1)
    @constraint(m, b := {(!((!(sum(x) > 3)) ⊻ (sum(x) >= 2))) && (x[1]+2x[2] <= 5)})
    optimize!(m)

    com = CS.get_inner_model(m)
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xval = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test !((!(sum(xval) > 3)) ⊻ (sum(xval) >= 2)) && (xval[1]+2xval[2] <= 5)
        push!(results, (xval[1], xval[2]))
    end
    for i in 0:5, j in 0:5
        x = [i,j]
        if !((!(sum(x) > 3)) ⊻ (sum(x) >= 2)) && (x[1]+2x[2] <= 5)
            @test (i,j) in results
        end
    end
end
=#


text:
# prompt.jl

using MIRT: prompt
using Test: @test

_tmp = prompt(:state) # save current state
prompt(:draw)
@test prompt(:state) === :draw
prompt()
@test prompt(_tmp) isa Symbol # return to original state

ext: jl


text:
using Test: @test
import Trixi

# Use a macro to avoid world age issues when defining new initial conditions etc.
# inside an elixir.
"""
    @test_trixi_include(elixir; l2=nothing, linf=nothing,
                                atol=10*eps(), rtol=0.001,
                                parameters...)

Test Trixi by calling `trixi_include(elixir; parameters...)`.
By default, only the absence of error output is checked.
If `l2` or `linf` are specified, in addition the resulting L2/Linf errors
are compared approximately against these reference values, using `atol, rtol`
as absolute/relative tolerance.
"""
macro test_trixi_include(elixir, args...)

  local l2   = get_kwarg(args, :l2, nothing)
  local linf = get_kwarg(args, :linf, nothing)
  local atol = get_kwarg(args, :atol, 500*eps())
  local rtol = get_kwarg(args, :rtol, sqrt(eps()))
  local kwargs = Pair{Symbol, Any}[]
  for arg in args
    if arg.head == :(=) && !(arg.args[1] in (:l2, :linf, :atol, :rtol))
      push!(kwargs, Pair(arg.args...))
    end
  end

  quote
    Trixi.mpi_isroot() && println("═"^100)
    Trixi.mpi_isroot() && println($elixir)

    # if `maxiters` is set in tests, it is usually set to a small numer to
    # run only a few steps - ignore possible warnings coming from that
    if any(==(:maxiters) ∘ first, $kwargs)
      additional_ignore_content = [
        r"┌ Warning: Interrupted\. Larger maxiters is needed\.\n└ @ SciMLBase .+\n",
        r"┌ Warning: Interrupted\. Larger maxiters is needed\.\n└ @ Trixi .+\n"]
    else
      additional_ignore_content = []
    end

    # evaluate examples in the scope of the module they're called from
    @test_nowarn_mod trixi_include(@__MODULE__, $elixir; $kwargs...) additional_ignore_content

    # if present, compare l2 and linf errors against reference values
    if !isnothing($l2) || !isnothing($linf)
      l2_measured, linf_measured = analysis_callback(sol)

      if Trixi.mpi_isroot() && !isnothing($l2)
        @test length($l2) == length(l2_measured)
        for (l2_expected, l2_actual) in zip($l2, l2_measured)
          @test isapprox(l2_expected, l2_actual, atol=$atol, rtol=$rtol)
        end
      end

      if Trixi.mpi_isroot() && !isnothing($linf)
        @test length($linf) == length(linf_measured)
        for (linf_expected, linf_actual) in zip($linf, linf_measured)
          @test isapprox(linf_expected, linf_actual, atol=$atol, rtol=$rtol)
        end
      end
    end

    Trixi.mpi_isroot() && println("═"^100)
    Trixi.mpi_isroot() && println("\n\n")
  end
end

# Get the first value assigned to `keyword` in `args` and return `default_value`
# if there are no assignments to `keyword` in `args`.
function get_kwarg(args, keyword, default_value)
  val = default_value
  for arg in args
    if arg.head == :(=) && arg.args[1] == keyword
      val = arg.args[2]
      break
    end
  end
  return val
end

# Modified version of `@test_nowarn` that prints the content of `stderr` when
# it is not empty. This is useful for debugging failing tests.
macro test_nowarn_debug(expr)
  quote
    let fname = tempname()
      try
        ret = open(fname, "w") do f
          redirect_stderr(f) do
            $(esc(expr))
          end
        end
        stderr_content = read(fname, String)
        if !isempty(stderr_content)
          println("Content of `stderr`:\n", stderr_content)
        end
        @test isempty(stderr_content)
        ret
      finally
        rm(fname, force=true)
      end
    end
  end
end

# Modified version of `@test_nowarn` that prints the content of `stderr` when
# it is not empty and ignnores module replacements.
macro test_nowarn_mod(expr, additional_ignore_content=String[])
  quote
    let fname = tempname()
      try
        ret = open(fname, "w") do f
          redirect_stderr(f) do
            $(esc(expr))
          end
        end
        stderr_content = read(fname, String)
        if !isempty(stderr_content)
          println("Content of `stderr`:\n", stderr_content)
        end

        # Patterns matching the following ones will be ignored. Additional patterns
        # passed as arguments can also be regular expressions, so we just use the
        # type `Any` for `ignore_content`.
        ignore_content = Any[
          # TODO: Upstream (PlotUtils). This should be removed again once the
          #       deprecated stuff is fixed upstream.
          "WARNING: importing deprecated binding Colors.RGB1 into PlotUtils.\n",
          "WARNING: importing deprecated binding Colors.RGB4 into PlotUtils.\n",
        ]
        append!(ignore_content, additional_ignore_content)
        for pattern in ignore_content
          stderr_content = replace(stderr_content, pattern => "")
        end

        # We also ignore simple module redefinitions for convenience. Thus, we
        # check whether every line of `stderr_content` is of the form of a
        # module replacement warning.
        @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content)
        ret
      finally
        rm(fname, force=true)
      end
    end
  end
end

"""
    @timed_testset "name of the testset" #= code to test #=

Similar to `@testset`, but prints the name of the testset and its runtime
after execution.
"""
macro timed_testset(name, expr)
  @assert name isa String
  quote
    local time_start = time_ns()
    @testset $name $expr
    local time_stop = time_ns()
    if Trixi.mpi_isroot()
      flush(stdout)
      @info("Testset " * $name * " finished in "
            * string(1.0e-9 * (time_stop - time_start)) * " seconds.\n")
      flush(stdout)
    end
  end
end

"""
    @trixi_testset "name of the testset" #= code to test #=

Similar to `@testset`, but wraps the code inside a temporary module to avoid
namespace pollution. It also `include`s this file again to provide the
definition of `@test_trixi_include`. Moreover, it records the execution time
of the testset similarly to [`timed_testset`](@ref).
"""
macro trixi_testset(name, expr)
  @assert name isa String
  # TODO: `@eval` is evil
  # We would like to use
  #   mod = gensym(name)
  #   ...
  #   module $mod
  # to create new module names for every test set. However, this is not
  # compatible with the dirty hack using `@eval` to get the mapping when
  # loading structured, curvilinear meshes. Thus, we need to use a plain
  # module name here.
  quote
    local time_start = time_ns()
    @eval module TrixiTestModule
      using Test
      using Trixi
      include(@__FILE__)
      # We define `EXAMPLES_DIR` in (nearly) all test modules and use it to
      # get the path to the elixirs to be tested. However, that's not required
      # and we want to fail gracefully if it's not defined.
      try
        import ..EXAMPLES_DIR
      catch
        nothing
      end
      @testset $name $expr
    end
    local time_stop = time_ns()
    if Trixi.mpi_isroot()
      flush(stdout)
      @info("Testset " * $name * " finished in "
            * string(1.0e-9 * (time_stop - time_start)) * " seconds.\n")
    end
    nothing
  end
end


text:
# This file contains types and functions for Choropleth Maps

# Abstract Map Classificator
abstract type AbstractMapClassificator end

abstract type AbstractGraduatedMapClassificator <: AbstractMapClassificator end

abstract type AbstractStatisticalMapClassificator  <: AbstractGraduatedMapClassificator end

abstract type AbstractUniqueMapClassificator <: AbstractMapClassificator end

abstract type AbstractCoLocationMapClassificator <: AbstractMapClassificator end

# Default color palettes
# This functions returns the default color palette and a boolean indicating if colors are reversed
defaultpalette(::AbstractGraduatedMapClassificator) = (:YlOrBr, false)

defaultpalette(::AbstractStatisticalMapClassificator) = (:RdBu, true)

defaultpalette(::AbstractUniqueMapClassificator) = (:Paired, false)

defaultpalette(mcr::AbstractCoLocationMapClassificator) = defaultpalette(mcr.mcr)

# Abstract Map Classification
abstract type AbstractMapClassification end

"""
    mapclassify (mcr::AbstractMapClassificator, x::Vector)
Classify observations in variable `x` in classess using the criterion specified in `mcr`.
"""
function mapclassify end

"""
    maplabels (mcAbstractMapClassification)
Get the labels of the classification `mc`.
# Optional Arguments
- `digits=2`: number of decimal digits.
- `sep=", "`: separator between class lower and upper bounds.
- `counts=true`: include the total number of observations on each class.
"""
function maplabels end

"""
    counts(mc::AbstractMapClassification)
Get the vector of classess sizes.

`counts(mc)[k]` is the number of observations assigned to the ``k``-class.
"""
counts(mc::AbstractMapClassification)::Vector{Int} = return mc.ngroup;

"""
    bounds(mc::AbstractMapClassification)
Get the lower and upper bounds of the classes.

`bounds(mc)[1][k]` is the lower bound of the ``k``-class.

`bounds(mc)[2][k]` is the upper bound of the ``k``-class.
"""
function bounds end

"""
    levels(mc::AbstractMapClassification)
Get the levels of the classes.

`levels(mc)[k]` is the level of the ``k``-class.
"""
function level end

"""
    assignments(mc::AbstractMapClassification)
Get the vector of classes indices for each observation.

`assignments(mc)[i]` is the index of the class to which the ``i``-th observation is assigned.
"""
assignments(mc::AbstractMapClassification)::Vector{Int} = return mc.group;

length(mc::AbstractMapClassification)::Int = return mc.k;


text:
using MAT, HDF5, JLD, DelimitedFiles, Statistics, Dates, TimeZones, Pkg.TOML

# If options[:datafolder] is set, then use that path. If not, use path listed in HOMEDIR/.GlobalEnergyGIS_config if available.
# Otherwise, fall back to Supergrid/inputdata where data for Europe8 is supplied.
function getdatafolder(options)
    if options[:datafolder] == ""
        configfile = joinpath(homedir(), ".GlobalEnergyGIS_config")
        if isfile(configfile)
            return joinpath(TOML.parsefile(configfile)["datafolder"], "output")
        else
            return abspath(dirname(@__FILE__), "..", "inputdata")
        end
    else
        return options[:datafolder]
    end
end

# eurasia21 regions for reference
# :NOR,:FRA,:GER,:UK,:MED,:BAL,:SPA,:CEN,:BUK,:TCC,:KZK,:CAS,:RU_C,:RU_SW,:RU_VL,:CH_N,:CH_NE,:CH_E,:CH_SC,:CH_SW,:CH_NW

function makesets(hourinfo, options)
    @unpack regionset, disableregions = options
    inputdata = getdatafolder(options)
    distancevars = matread(joinpath(inputdata, "distances_$regionset.mat"))
    dataregions = Symbol.(vec(distancevars["regionlist"]))
    regionlist = setdiff(dataregions, disableregions)
    makesets(regionlist, dataregions, hourinfo, options)
end

function makesets(REGION, dataregions, hourinfo, options)
    @unpack datayear, regionset = options
    techdata = Dict(
        :name => [:pv,  :pvroof, :csp,     :wind, :offwind, :hydro,    :coal,    :gasGT,   :gasCCGT, :bioGT,   :bioCCGT, :nuclear, :battery],
        :type => [:vre, :vre,    :storage, :vre,  :vre,     :storage,  :thermal, :thermal, :thermal, :thermal, :thermal, :thermal, :storage],
        :fuel => [:_,   :_,      :_,       :_,    :_,       :_,        :coal,    :gas,     :gas,     :biogas,  :biogas,  :uranium, :_]
    )

    inputdata = getdatafolder(options)
    hydrovars = matread(joinpath(inputdata, "GISdata_hydro_$regionset.mat"))
    _, ncostclasses, nreservoirclasses = size(hydrovars["potentialcapac"])
    nvreclasses = 5

    numtechs = length(techdata[:name])
    reservoirs = collect('a':'z')[1:nreservoirclasses]
    vreclass = [Symbol("$letter$number") for letter in ["a", "b"] for number = 1:nvreclasses]
    hydroclass = [:x0;  [Symbol("$letter$number") for letter in reservoirs for number = 1:ncostclasses]]
    noclass = [:_]
    techtype = Dict(techdata[:name][i] => techdata[:type][i] for i=1:numtechs)
    techfuel = Dict(techdata[:name][i] => techdata[:fuel][i] for i=1:numtechs)

    TECH = techdata[:name]
    FUEL = [:_, :coal, :gas, :biogas, :uranium]
    CLASS = Dict(k =>
        k == :hydro ? hydroclass :
        k == :csp || techtype[k] == :vre ? vreclass :
        noclass for k in TECH)
    CLASS[:transmission] = noclass
    STORAGECLASS = Dict(k => [:_] for k in TECH)
    STORAGECLASS[:hydro] = [:x0;  Symbol.(reservoirs)]
    STORAGECLASS[:csp] = vreclass

    reservoirclass = Dict(r => [Symbol("$r$number") for number = 1:ncostclasses] for r in Symbol.(reservoirs))
    reservoirclass[:x0] = [:x0]
    reservoirclass[:_] = [:_]
    for vc in vreclass
        reservoirclass[vc] = [vc]
    end

    HOUR = 1:(24*Dates.daysinyear(datayear)÷hourinfo.hours)

    return Sets(REGION, FUEL, TECH, CLASS, STORAGECLASS, HOUR, techtype, techfuel, reservoirclass, dataregions)
end

# resample hour dimension of array a (indicated by hourdim) using hourindexes in hourinfo structure,
# then reduce hours further by hours
function reducehours(a, hourdim, hourinfo)
    hours = hourinfo.hours
    aa = copy(selectdim(a, hourdim, hourinfo.hourindexes))
    out = copy(selectdim(aa, hourdim, 1:hours:size(aa,hourdim)))    # sample every nth hour
    if true     # true: averaging   false: sampling
        for i = 2:hours
            out += copy(selectdim(aa, hourdim, i:hours:size(aa,hourdim)))
        end
        out = out / hours
    end
    return out
end

# reduce regions from 10 (in Lina's input data) to 8 (in model)
# :MED = :IT + :GR,    :BAL (new) = :BAL (old) + :POL
function ten2eight(a)
    # REG10 = [:NOR, :IT, :FRA, :GER, :UK, :GR, :BAL, :POL, :SPA, :CEN]
    # REGION = [:NOR, :FRA, :GER, :UK, :MED, :BAL, :SPA, :CEN]
    out = a[[1,3,4,5,6,7,9,10],:]
    out[5,:] = a[2,:] + a[6,:]
    out[6,:] = a[7,:] + a[8,:]
    return out
end

CRF(r,T) = r / (1 - 1/(1+r)^T)

function makeparameters(sets, options, hourinfo)
    @unpack REGION, FUEL, TECH, CLASS, HOUR, dataregions = sets
    @unpack discountrate, datayear, regionset, solarwindarea, islandindexes,
            inputdatasuffix, sspscenario, sspyear = options

    hoursperyear = 24 * Dates.daysinyear(datayear)
    hoursperperiod = Int(hourinfo.hoursperperiod)

    initialstoragelevel = 0.7       # make this tech dependent later
    minflow_existinghydro = 0.4
    cspsolarmultiple = 3.0          # peak capacity of collectors divided by turbine power
    cspthermalstoragehours = 12

    numregions = length(REGION)
    nhours = length(HOUR)
    nhydro = length(CLASS[:hydro])
    nclasses = length(CLASS[:pv])

    activeregions = [r in REGION for r in dataregions]
    demand = AxisArray(zeros(numregions, nhours), REGION, HOUR)     # GW

    inputdata = getdatafolder(options)

    # read synthetic demand data (in UTC)
    gisdemand = JLD.load(joinpath(inputdata,
        "SyntheticDemand_$(regionset)_$sspscenario-$(sspyear)_$datayear.jld"), "demand")
    for i = 1:numregions
        demand[i,:] = reducehours(gisdemand[:,i], 1, hourinfo) / 1000       # GW
    end

    hydrovars = matread(joinpath(inputdata, "GISdata_hydro_$regionset.mat"))
    hydrocapacity = AxisArray(zeros(numregions,nhydro), REGION, CLASS[:hydro])
    hydroeleccost = AxisArray(zeros(numregions,nhydro), REGION, CLASS[:hydro])
    monthlyinflow = AxisArray(zeros(numregions,nhydro,12), REGION, CLASS[:hydro], 1:12)
    cfhydroinflow = AxisArray(zeros(numregions,nhydro,nhours), REGION, CLASS[:hydro], HOUR)
    dischargetime = AxisArray(zeros(numregions,3,1+nhydro+nclasses), REGION, [:hydro,:battery,:csp], [CLASS[:hydro]; CLASS[:csp]; :_])

    hydrocapacity[:,:x0] = typeof(hydrovars["existingcapac"]) == Float64 ? [hydrovars["existingcapac"]] : hydrovars["existingcapac"][activeregions]
    hydrocapacity[:,2:end] = reshape(hydrovars["potentialcapac"][activeregions,:,:], numregions, nhydro-1)
    hydrocapacity[isnan.(hydrocapacity)] = zeros(sum(isnan.(hydrocapacity)))

    # eleccost = capcost * crf / (CF * 8760)  =>   eleccost2/eleccost1 = crf2/crf1
    # 1$ = 0.9€ (average 2015-2017)
    hydroeleccost[:,2:end] = reshape(hydrovars["potentialmeancost"][activeregions,:,:], numregions, nhydro-1)   # $/kWh with 10% discount rate
    hydroeleccost[:,:] = hydroeleccost[:,:] * CRF(discountrate,40)/CRF(0.1,40) * 0.9 * 1000     # €/MWh    (0.9 €/$)
    hydroeleccost[isnan.(hydroeleccost)] = fill(999, sum(isnan.(hydroeleccost)))

    monthlyinflow[:,:x0,:] = hydrovars["existinginflowcf"][activeregions,:]
    monthlyinflow[:,2:end,:] = reshape(hydrovars["potentialinflowcf"][activeregions,:,:,:], numregions, nhydro-1, 12)
    monthlyinflow[isnan.(monthlyinflow)] = zeros(sum(isnan.(monthlyinflow)))

    dischargetime[:,:hydro,:x0] .= 168*6    # assume average discharge time 6 weeks by default
    # if we have data on national storage capacity then use it
    hydrostoragecapacity = Dict(:NOR => 121.43, :FRA => 3.59, :MED => 9.2, :SPA => 16.6, :CEN => 7.4)   # TWh
    for (reg, hydroenergy) in hydrostoragecapacity
        if reg in REGION
            dischargetime[reg,:hydro,:x0] = hydroenergy/hydrocapacity[reg,:x0] * 1000
        end
    end
    # these countries have mostly run-of-river hydro so discharge times are lower
    dischargetime[intersect([:GER,:UK,:BAL],REGION),:hydro,:x0] .= 300
    # use the GIS data for discharge time for future hydro
    dischargetime[:,:hydro,2:end-1-nclasses] = reshape(hydrovars["potentialmeandischargetime"][activeregions,:,:], numregions, nhydro-1)

    dischargetime[:,:battery,:_] .= 1
    dischargetime[:,:csp,:] .= cspthermalstoragehours
    dischargetime[isnan.(dischargetime)] = fill(10000, sum(isnan.(dischargetime)))
    dischargetime[dischargetime .> 10000] = fill(10000, sum(dischargetime .> 10000))

    # monthly to hourly hydro inflow
    dayspermonth = Dates.daysinmonth.(Date.(datayear, 1:12))
    lasthour = 24 ÷ hoursperperiod * cumsum(dayspermonth)
    firsthour = [1; 1 .+ lasthour[1:end-1]]
    for m = 1:12
        for i = firsthour[m]:lasthour[m]
            cfhydroinflow[:,:,i] = monthlyinflow[:,:,m]
        end
    end
    cfhydroinflow[cfhydroinflow .< 0.01] = zeros(sum(cfhydroinflow .< 0.01))

    # read regional distances from Matlab file
    distancevars = matread(joinpath(inputdata, "distances_$regionset.mat"))

    if typeof(distancevars["distances"]) == Float64
        distances = fill(distancevars["distances"], 1, 1)
        connected = fill(distancevars["connected"], 1, 1)
        connectedoffshore = fill(distancevars["connectedoffshore"], 1, 1)
    else
        distances = distancevars["distances"][activeregions,activeregions]
        connected = distancevars["connected"][activeregions,activeregions]
        connectedoffshore = distancevars["connectedoffshore"][activeregions,activeregions]
    end

    # add connection from south-central China to Germany
    indexGER = findfirst(REGION .== :GER)
    indexCHSC = findfirst(REGION .== :CH_SC)
    if indexGER != nothing && indexCHSC != nothing
        connected[indexGER,indexCHSC] = connected[indexCHSC,indexGER] = true
    end

    # from Bogdanov & Breyer (2016) "North-East Asian Super Grid..."
    transmissioncostdata = connected .* (180 .+ 0.612*distances) .+ connectedoffshore .* (180 .+ 0.992*distances)
    transmissionfixedcostdata = connected .* (1.8 .+ 0.0075*distances) .+ connectedoffshore .* (1.8 .+ 0.0010*distances)
    transmissioninvestcost = AxisArray(transmissioncostdata, REGION, REGION)        # €/kW
    transmissionfixedcost = AxisArray(transmissionfixedcostdata, REGION, REGION)        # €/kW
    transmissionlossdata = (connected .| connectedoffshore) .* (0.014 .+ 0.016*distances/1000)
    transmissionlosses = AxisArray(transmissionlossdata, REGION, REGION)
    smalltransmissionpenalty = 0.1      # €/MWh elec

    islands = zeros(Bool, length(dataregions), length(dataregions))
    for ndx in islandindexes
        islands[ndx,ndx] .= true
    end
    transmissionislands = AxisArray(islands[activeregions,activeregions], REGION, REGION)

    # Efficiencies for storage technologies are round trip efficiencies.
    # CSP costs are adjusted for solar field size and storage capacity further below.
    techtable = [
        #               investcost  variablecost    fixedcost   lifetime    efficiency  rampingrate
        #               €/kW        €/MWh elec      €/kW/year   years                   share of capacity per hour
        :gasGT          500         1               10          30          0.4         1
        :gasCCGT            800         1               16          30          0.6         0.3
        :coal           1600        2               48          30          0.45        0.15
        :bioGT          500         1               10          30          0.4         1
        :bioCCGT        800         1               16          30          0.6         0.3
        :nuclear        5000        3               150         50          0.4         0.05
        :wind           1200        0               43          25          1           1
        :offwind        2300        0               86          25          1           1
        :transmission   NaN         0               NaN         50          NaN         1
        :battery        150         0.1             1.5         10          0.9         1   # 1h discharge time, 150 €/kW = 150 €/kWh
        :pv             600         0               16          25          1           1
        :pvroof         900         0               20          25          1           1
        :csp            6000        0               35          30          1           1   # for solar multiple=3, storage=12 hours
        :hydro          10          0               0           80          1           1   # small artificial investcost so it doesn't overinvest in free capacity
    ]
    techs = techtable[:,1]
    techdata = Float64.(techtable[:,2:end])
    baseinvestcost = AxisArray(techdata[:,1], techs)    # €/kW
    variablecost = AxisArray(techdata[:,2], techs)      # €/MWh elec
    fixedcost = AxisArray(techdata[:,3], techs)         # €/kW/year
    lifetime = AxisArray(techdata[:,4], techs)          # years
    efficiency = AxisArray(techdata[:,5], techs)
    rampingrate = AxisArray(techdata[:,6], techs)

    # fuel cost references:
    # 1. https://www.gov.uk/government/statistical-data-sets/prices-of-fuels-purchased-by-major-power-producers
    # 2. https://oilprice.com/Energy/Natural-Gas/European-Natural-Gas-Prices-Are-Set-To-Rise-Further.html
    # 3. https://www.bloomberg.com/news/articles/2018-09-03/coal-nears-100-in-europe-as-china-s-power-demand-draws-in-fuel
    # 4. https://www.frisch.uio.no/ressurser/LIBEMOD/pdf/phase_out_28mai2015_golombek_powerpoint.pdf
    # 5. http://www.world-nuclear.org/information-library/economic-aspects/economics-of-nuclear-power.aspx
    # 6. https://www.gasforclimate2050.eu/files/files/Ecofys_Gas_for_Climate_Report_Study_March18.pdf  (page 19, section 3.4, figure 7)
    # Sepulveda/Jenkins:  gas & biogas 21 €/MWh, uranium 3 €/MWh
    # 1: coal 11-12 €/MWh, gas 19-21 €/MWh,  2: gas 20-30 €/MWh,  3: coal 70-90 €/ton = (25 MJ/kg) = 2.8-3.6 €/GJ = 10-13 €/MWh
    # 4: gas 45*.5 = 22 €/MWh, coal 22.3*.4 = 9 €/MWh, bio 26.4*.4 = 11 €/MWh, nuclear 6.7*.35 = 2.3 €/MWh
    # 5: nuclear fuel pellets 0.39 USD/kWh = 3.2 €/MWh    (see also Stuff/Nuclear fuel costs.xlsx)
    # 6: biomethane anaerobic digestion 96 €/MWh (2015), 60 €/MWh (2050), thermal gasification 37 €/MWh (2050)

    fuelcost = AxisArray(Float64[0, 11, 22, 37, 3.2], [:_, :coal, :gas, :biogas, :uranium])     # €/MWh fuel

    crf = AxisArray(discountrate ./ (1 .- 1 ./(1+discountrate).^lifetime), techs)

    emissionsCO2 = AxisArray(zeros(length(FUEL)), FUEL)
    emissionsCO2[[:coal,:gas]] = [0.330, 0.202]     # kgCO2/kWh fuel (or ton/MWh or kton/GWh)

    # do something with B classes (and pvrooftop) later
    windvars = matread(joinpath(inputdata, "GISdata_wind$(datayear)_$regionset$inputdatasuffix.mat"))
    solarvars = matread(joinpath(inputdata, "GISdata_solar$(datayear)_$regionset$inputdatasuffix.mat"))

    allclasses = union(sets.CLASS[:pv], sets.CLASS[:hydro], [:_])
    cf = AxisArray(ones(numregions,length(TECH),length(allclasses),nhours), REGION, TECH, allclasses, HOUR)
    classlimits = AxisArray(zeros(numregions,5,length(CLASS[:pv])), REGION, [:wind, :offwind, :pv, :pvroof, :csp], CLASS[:pv])
    solarcombinedarea = AxisArray(zeros(numregions,length(CLASS[:pv]),length(CLASS[:csp])), REGION, CLASS[:pv], CLASS[:csp])

    cf[:,:wind,1:5,:] = permutedims(reducehours(windvars["CFtime_windonshoreA"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:offwind,1:5,:] = permutedims(reducehours(windvars["CFtime_windoffshore"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:pv,1:5,:] = permutedims(reducehours(solarvars["CFtime_pvplantA"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:pvroof,1:5,:] = permutedims(reducehours(solarvars["CFtime_pvrooftop"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:csp,1:5,:] = permutedims(reducehours(solarvars["CFtime_cspplantA"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:wind,6:10,:] = permutedims(reducehours(windvars["CFtime_windonshoreB"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:pv,6:10,:] = permutedims(reducehours(solarvars["CFtime_pvplantB"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[:,:csp,6:10,:] = permutedims(reducehours(solarvars["CFtime_cspplantB"][:,activeregions,:], 1, hourinfo), [2,3,1])
    cf[isnan.(cf)] = zeros(sum(isnan.(cf)))
    cf[cf .< 0.01] = zeros(sum(cf .< 0.01))     # set small values to 0 for better numerical stability

    classlimits[:,:wind,1:5] = windvars["capacity_onshoreA"][activeregions,:] * solarwindarea
    classlimits[:,:offwind,1:5] = windvars["capacity_offshore"][activeregions,:] * (1 + 0.5*(solarwindarea - 1))    # only half because default offshore area is 33%
    classlimits[:,:pv,1:5] = solarvars["capacity_pvplantA"][activeregions,:] * solarwindarea
    classlimits[:,:pvroof,1:5] = solarvars["capacity_pvrooftop"][activeregions,:] * solarwindarea
    classlimits[:,:csp,1:5] = solarvars["capacity_cspplantA"][activeregions,:] * solarwindarea
    classlimits[:,:wind,6:10] = windvars["capacity_onshoreB"][activeregions,:] * solarwindarea
    classlimits[:,:pv,6:10] = solarvars["capacity_pvplantB"][activeregions,:] * solarwindarea
    classlimits[:,:csp,6:10] = solarvars["capacity_cspplantB"][activeregions,:] * solarwindarea

    solarcombinedarea[:,1:5,1:5] = solarvars["solar_overlap_areaA"][activeregions,:,:] * solarwindarea
    solarcombinedarea[:,6:10,6:10] = solarvars["solar_overlap_areaB"][activeregions,:,:] * solarwindarea
    pv_density = solarvars["pv_density"]
    csp_density = solarvars["csp_density"]

    investcost = AxisArray(zeros(length(techs),length(allclasses)), techs, allclasses)  # €/kW
    for k in techs, c in CLASS[k]
        investcost[k,c] = baseinvestcost[k]
    end
    for k in [:wind,:pv,:csp]
        investcost[k,6:10] .+= 200
    end

    # CSP solar tower cost analysis (see CSP cost analysis.xlsx):
    # Costs vary with the solar multiple (collector field size) and thermal storage capacity (in hours).
    # Assume a solar tower plant with solar multiple = 3 and storage = 12 hours costs 9200 USD/kW (2010).
    # Then assume solar field costs are 35% of total costs and storage about 10% of total costs (roughly inline with IRENA 2012, fig 4.4)
    # for this plant, and that costs for plants with other parameters vary linearly with solar multiple and storage size.
    # Finally convert to euro (1.15 USD/EUR) and assume a cost reduction of 25% to 2050. Resulting cost for a 3/12 solar tower: 6000 €/kW.
    investcost[:csp,:] = investcost[:csp,:] * (0.55 + 0.35*cspsolarmultiple/3 + 0.10*cspthermalstoragehours/12)

    cf[:,:csp,:,:] = cf[:,:csp,:,:] * cspsolarmultiple                          # OK if this surpasses 100%
    classlimits[:,:csp,:] = classlimits[:,:csp,:] / cspsolarmultiple            # GIS data calculated as peak solar power
    # The Capacity variable of the model should now be correct (turbine capacity, with a larger solar collector field)

    # # check demand/solar synchronization
    # plotly()
    # for r = 1:numregions
    #   tt1 = 8760÷2÷hoursperperiod     # test winter, spring, summer (÷12, ÷3, ÷2)
    #   tt = tt1:tt1+2*24÷hoursperperiod
    #   qq = [demand[r,tt] maximum(cf[r,:pv,1:5,tt],dims=1)']
    #   display(plot(qq./maximum(qq,dims=1), size=(1850,950)))
    # end

    return Params(cf, transmissionlosses, demand, hydrocapacity, cfhydroinflow, classlimits, transmissionislands,
        efficiency, rampingrate, dischargetime, initialstoragelevel, minflow_existinghydro, emissionsCO2, fuelcost,
        variablecost, smalltransmissionpenalty, investcost, crf, fixedcost, transmissioninvestcost, transmissionfixedcost,
        hydroeleccost, solarcombinedarea, pv_density, csp_density, cspsolarmultiple)
end

# Run fix_timezone_error() if an error like this is produced (the build step should take care if this for most people):
# (I'm not even sure if this helps...)
#
# julia> r, annualelec, capac, tcapac, chart = runmodel(regionset=:europe8, carboncap=0.1, hours=3);
#
# Reading input data...
#   1.257888 seconds (1.77 M allocations: 92.844 MiB, 2.44% gc time)
# ERROR: UnhandledTimeError: TimeZone Europe/Oslo does not handle dates on or after 2038-03-28T01:00:00 UTC
fix_timezone_error() = TimeZones.TZData.compile(max_year=2200)

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5
@testset "ic0" begin
    @testset for elty in [Float32,Float64]
        @testset "csr" begin
            A    = rand(elty, m ,m)
            A    = A + transpose(A)
            A   += m * eye(elty ,m)
            d_A  = Symmetric(CudaSparseMatrixCSR(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'S', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'S', d_A, info, 'O')
            h_A  = to_host(d_B)
            Ac   = sparse(full(cholfact(A)))
            h_A  = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A    = rand(elty, m ,m)
            A    = A + transpose(A)
            A   += m * eye(elty ,m)
            d_A  = Symmetric(CudaSparseMatrixCSC(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'S', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'S', d_A, info, 'O')
            h_A  = to_host(d_B)
            Ac   = sparse(full(cholfact(A)))
            h_A  = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
    @testset for elty in [Complex64,Complex128]
        @testset "csr" begin
            A    = rand(elty,m,m)
            A    = A + ctranspose(A)
            A   += m * eye(elty,m)
            d_A  = Hermitian(CudaSparseMatrixCSR(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'H', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'H', d_A, info, 'O')
            h_A  = to_host(d_B)
            Ac   = sparse(full(cholfact(A)))
            h_A  = ctranspose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A    = rand(elty,m,m)
            A    = A + ctranspose(A)
            A   += m * eye(elty,m)
            d_A  = Hermitian(CudaSparseMatrixCSC(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'H', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'H', d_A, info, 'O')
            h_A  = to_host(d_B)
            Ac   = sparse(full(cholfact(A)))
            h_A  = ctranspose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end
@testset "ic2" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "csr" begin
            A   = rand(elty, m, m)
            A  += ctranspose(A)
            A  += m * eye(elty, m)
            d_A = CudaSparseMatrixCSR(sparse(tril(A)))
            d_B = CUSPARSE.ic02(d_A, 'O')
            h_A = to_host(d_B)
            Ac  = sparse(full(cholfact(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            A   = rand(elty,m,n)
            d_A = CudaSparseMatrixCSR(sparse(tril(A)))
            @test_throws DimensionMismatch CUSPARSE.ic02(d_A, 'O')
        end
        @testset "csc" begin
            A   = rand(elty, m, m)
            A  += ctranspose(A)
            A  += m * eye(elty, m)
            d_A = CudaSparseMatrixCSC(sparse(tril(A)))
            d_B = CUSPARSE.ic02(d_A, 'O')
            h_A = to_host(d_B)
            Ac  = sparse(full(cholfact(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            A   = rand(elty,m,n)
            d_A = CudaSparseMatrixCSC(sparse(tril(A)))
            @test_throws DimensionMismatch CUSPARSE.ic02(d_A, 'O')
        end
    end
end


text:
using Test
using MultivariateBases
using DynamicPolynomials
@polyvar x y

@testset "Linear" begin
    basis = MonomialBasis([x, y])
    @test polynomialtype(basis, Int) == polynomialtype(x, Int)
    @test polynomial(i -> i^2, basis) == x + 4y
    @test coefficients(x + 4y, MonomialBasis) == [1, 4]
end
@testset "Affine" begin
    # It will be sorted and 1 will be moved at the end
    basis = MonomialBasis([1, x, y])
    @test polynomialtype(basis, Int) == polynomialtype(x, Int)
    @test polynomial(i -> i^2, basis) == x + 4y + 9
    @test coefficients(9 + x + 4y, MonomialBasis) == [1, 4, 9]
end
@testset "Quadratic" begin
    basis = MonomialBasis([x^2, x*y, y^2])
    @test polynomialtype(basis, Int) == polynomialtype(x, Int)
    @test polynomial(i -> i^2, basis) == x^2 + 4x*y + 9y^2
    @test coefficients(x^2 + 4x*y + 9y^2, MonomialBasis) == [1, 4, 9]
end
@testset "API degree = $degree" for degree in 0:3
    api_test(MonomialBasis, degree)
end

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
using Test, LinearAlgebra, Polyhedra, XPORTA

@testset "test/unit/LocalPolytope.jl" begin

using BellScenario

@testset "running all tests for module LocalPolytope.jl" begin

    @testset "LocalPolytope.vrep()" begin
        @testset "no arguments" begin
            scenario = LocalSignaling(3,3,2)
            local_poly = vrep(scenario)

            @test local_poly isa XPORTA.Polyhedron
            @test npoints(local_poly) == 21

            vertices = map(p -> convert.(Int64,p), points(local_poly))
            @test LocalPolytope.vertices(scenario) == vertices
        end

        @testset "passing arguments" begin
            scenario = LocalSignaling(3,3,2)
            local_poly = vrep(scenario, rank_d_only=true)

            @test local_poly isa XPORTA.Polyhedron
            @test npoints(local_poly) == 18

            vertices = map(p -> convert.(Int64,p), points(local_poly))
            @test LocalPolytope.vertices(scenario, rank_d_only=true) == vertices
        end
    end

    println("running LocalPolytope.jl unit tests:")
    for test in readdir("./test/unit/LocalPolytope/")
    # run only julia files in test directory
    if occursin(r"^.*\.jl$", test)
        println("./unit/LocalPolytope/$test")
        @time include("./LocalPolytope/$test")
    end
    end

end

end


text:
### contractile component
# using Parameters
# CC parameter struct
@with_kw struct CCParameters
    @deftype Float64
    Tmax
    T0
    ωmax
    ωc
    Amin
    ω1
    ωr
    θopt
    R
end
# create from dictionary
CCParameters(d::Dict{Symbol,Float64}) = CCParameters(Tmax = d[:Tmax], T0 = d[:T0], ωmax = d[:wmax], ωc = d[:wc], Amin = d[:Amin], ω1 = d[:w1], ωr = d[:wr], θopt = d[:theta_opt], R = d[:R])

# CC struct
mutable struct CC
    α::Float64  # activation level
    τ::Float64  # torque
    θ::Float64  # CC angle
    ω::Float64  # CC angular velocity
    τ_θ::Float64 # normalised torque angle value for current CC angle
    τ_ω::Float64 # normalised torque angular velocity value for current CC ang vel
    cc_p::CCParameters # 9 CC parameters
    α_p::ActivationProfile # activation parameters

    function CC(θ, ω, cc_p::CCParameters, α_p::ActivationProfile)
        α = α_p.a0
        τ_θ = torque_angle(θ, cc_p.θopt, cc_p.R)
        τ_ω = torque_angvel_norm(-ω, cc_p)
        τ = α * τ_θ * τ_ω * cc_p.T0

        return new(α, τ, θ, ω, τ_θ, τ_ω, cc_p, α_p)
    end
end

# calculate torque from CC angle
torque_angle(θ, θopt, R) = exp((-(θ - θopt)^2) / (2.0 * R^2))
torque_angle(cc::CC) = torque_angle(cc.θ, cc.cc_p.θopt, cc.cc_p.R)
torque_angle!(cc::CC) = cc.τ_θ = torque_angle(cc)
torque_angle(θ, p::CCParameters) = torque_angle(θ, p.θopt, p.R)

# caluate torque given CC angular velocity
function torque_angvel(ω, Tmax, T0, ωmax, ωc, Amin, ω1, ωr)
    # constants
    k = 4.3
    Amax = 1.0 # for differential activation

    # intermediate values
    Tc = T0 * ωc / ωmax
    C = Tc * (ωmax + ωc)
    ωe = ((Tmax - T0) / (k * T0)) * (ωmax * ωc / (ωmax + ωc))
    E = -(Tmax - T0) * ωe

    # differential activation
    α = Amin + (Amax - Amin) / (1 + exp(-(ω - ω1) / ωr))

    # torque zero when greater than ωmax
    ω > ωmax && return 0.0

    # otherwise calculate torque
    if ω > 0.0
        T = C / (ωc + ω) - Tc   # concentric torque
    else
        T = E / (ωe - ω) + Tmax
    end

    return T * α
end

# reverse sign for τ-ω function if passing CC, leave as is if passing ω
torque_angvel(cc::CC) = torque_angvel(-cc.ω, cc.cc_p.Tmax, cc.cc_p.T0, cc.cc_p.ωmax, cc.cc_p.ωc, cc.cc_p.Amin, cc.cc_p.ω1, cc.cc_p.ωr)
torque_angvel(ω, p::CCParameters) = torque_angvel(ω, p.Tmax, p.T0, p.ωmax, p.ωc, p.Amin, p.ω1, p.ωr)
torque_angvel_norm(cc::CC) = torque_angvel(-cc.ω, cc.cc_p.Tmax / cc.cc_p.T0, 1.0, cc.cc_p.ωmax, cc.cc_p.ωc, cc.cc_p.Amin, cc.cc_p.ω1, cc.cc_p.ωr)
torque_angvel_norm(ω, p::CCParameters) = torque_angvel(ω, p.Tmax / p.T0, 1.0, p.ωmax, p.ωc, p.Amin, p.ω1, p.ωr)
torque_angvel!(cc::CC) = cc.τ_ω = torque_angvel_norm(cc)

# calculate CC angular velocity given NORMALISED torque - NO DIFFERENTIAL ACTIVATION
function angvel_torque(τ_ω, Tmax, T0, ωmax, ωc)
    # constants
    k = 4.3

    # intermediate values
    Tc = T0 * ωc / ωmax
    C = Tc * (ωmax + ωc)
    ωe = ((Tmax - T0) / (k * T0)) * (ωmax * ωc / (ωmax + ωc))
    E = -(Tmax - T0) * ωe

    # error check
    τ_ω < 0 && error("normalised torque less than zero: $τ_ω")

    # find angular velocity
    if τ_ω ≤ T0
        ω = C / (τ_ω + Tc) - ωc   # concentric, +ve angvel
    elseif τ_ω > T0 && τ_ω < Tmax
        ω = ωe - E / (τ_ω - Tmax) # eccentric, -ve angvel
        ω < -ωmax && (ω = -ωmax)  # limit max eccentric angvel
    else
        ω = -ωmax   # τ_ω > Tmax, set to ωmax to prevent it blowing up
    end

    return ω
end

# reverse sign of ω if passing CC, keep as if passing τ_ω
angvel_torque(cc::CC) = -angvel_torque(cc.τ_ω, cc.cc_p.Tmax, cc.cc_p.T0, cc.cc_p.ωmax, cc.cc_p.ωc)
angvel_torque(τ_ω, p::CCParameters) = angvel_torque(τ_ω, p.Tmax, p.To, p.ωmax, p.ωc)
angvel_torque!(cc::CC) = cc.ω = angvel_torque(cc) # not negated as called function negates
angvel_torque_norm(cc::CC) = -angvel_torque(cc.τ_ω, cc.cc_p.Tmax / cc.cc_p.T0, 1.0, cc.cc_p.ωmax, cc.cc_p.ωc)
angvel_torque_norm(τ_ω, p::CCParameters) = angvel_torque(τ_ω, p.Tmax / p.T0, 1.0, p.ωmax, p.ωc)
angvel_torque_norm!(cc::CC) = cc.ω = angvel_torque_norm(cc) # not negated as called function negates

# calculate CC activation
activation(cc::CC, t) = activation(t, cc.α_p)
activation!(cc::CC, t) = cc.α = activation(cc, t)

# calculate contractile component torque from its current state (angle and angvel) and activation
torque(cc::CC, t) = activation(t, cc.α_p) * torque_angle(cc) * torque_angvel(cc)
torque!(cc::CC, t) = cc.τ = torque(cc, t)
# ω sign not reversed as a specific value being called
torque(t, θ, ω, cc_p::CCParameters, α_p::ActivationProfile) = activation(t, α_p) * torque_angle(θ, cc_p) * torque_angvel(ω, cc_p)
torque(α, θ, ω, p::CCParameters) = α * torque_angle(θ, p) * torque_angvel(ω, p)
torque_norm_tet(θ, ω, p::CCParameters) = torque_angle(θ, p) * torque_angvel_norm(ω, p)

# extend plots
function Plots.plot(p::CCParameters)
    p1 = plot(ω -> torque_angvel_norm(ω, p), -p.ωmax:0.01:p.ωmax, xaxis = "angular velocity", yaxis = "normalised torque", legend = :none)
    p2 = plot(θ -> torque_angle(θ, p), (p.θopt-3p.θopt):0.01:(p.θopt+3p.θopt), xaxis = "angle", yaxis = "normalised torque", legend = :none)
    p3 = plot(τ -> angvel_torque_norm(τ, p), 0.0:0.01:2.0, xaxis = "normalised torque", yaxis = "CC angular velocity", legend = :none)
    plot(p1, p2, p3, layout = (1, 3))
end

Plots.surface(p::CCParameters) = surface((p.θopt-3p.θopt):0.01:(p.θopt+3p.θopt), -p.ωmax:0.01:p.ωmax, (θ, ω) -> torque_norm_tet(θ, ω, p), xaxis = "angle", yaxis = "angular velocity", zaxis = "normalised torque")


text:
# This example came from the R package ggm.

using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir)

df = DataFrame!(CSV.File(scm_path("..", "data", "marks.csv")));

#=
R_dag = "
           mechanics vectors algebra statistics analysis
mechanics          0       0       0          0        0
vectors            1       0       0          0        0
algebra            1       1       0          1        1
statistics         0       0       0          0        0
analysis           0       0       0          1        0
";
=#

d = from_ggm("DAG(
    mechanics ~ vectors+algebra,
    vectors ~ algebra,
    statistics ~ algebra+analysis,
    analysis ~ algebra)"
)
display(d)

# DAG accepts either an OrderedDict, an adjacency_matrix or a ggm/dagitty string.
# Below d_string holds a ggm DAG definition.

d_string = "DAG(
    mechanics ~ vectors+algebra,
    vectors ~ algebra,
    statistics ~ algebra+analysis,
    analysis ~ algebra)"

dag = DAG("marks", d_string, df=df);
show(dag)

fname = ProjDir * "/marks.dot"
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

display(dag.s); println()

bs = basis_set(dag)
display(bs); println()

t = shipley_test(dag)
display("shipley_test = $t"); println()

pt = pcor_test(dag, [:analysis, :statistics, :mechanics], 1, 88)
display("pcor_test = $pt"); println()

f = [:statistics]; s = [:mechanics]; sel = vcat(f, s)

e = d_separation(dag, f, s; c=:algebra)
println("d_separation($(dag.name), $f, $s; c=:algebra) = $e")

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e")

print("d_separation($(dag.name), [:statistics], [:mechanics]; c=[:vectors]) = ")
println(d_separation(dag, [:statistics], [:mechanics]; c=[:vectors]))

print("d_separation($(dag.name), [:statistics], [:mechanics]; c=[:analysis, :vectors]) = ")
println(d_separation(dag, [:statistics], [:mechanics]; c=[:analysis, :vectors]))

print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics]; c=[:algebra]) = ")
println(d_separation(dag, [:statistics, :analysis], [:mechanics]; c=[:algebra]))

print("d_separation($(dag.name), [:statistics], [:mechanics, :vectors]; c=[:algebra]) = ")
println(d_separation(dag, [:statistics], [:mechanics, :vectors]; c=[:algebra]))

print("d_separation($(dag.name), [:statistics], [:mechanics, :analysis]; c=[:algebra]) = ")
println(d_separation(dag, [:statistics], [:mechanics, :analysis]; c=[:algebra]))

print("d_separation($(dag.name), [:analysis], [:vectors]) = ")
println(d_separation(dag, [:analysis], [:vectors]))

print("d_separation($(dag.name), [:analysis], [:vectors]; c=[:algebra]) = ")
println(d_separation(dag, [:analysis], [:vectors]; c=[:algebra]))

print("d_separation($(dag.name), [:vectors], [:statistics]; c=[:algebra]) = ")
println(d_separation(dag, [:analysis], [:vectors]; c=[:algebra]))

print("d_separation($(dag.name), [:statistics], [:algebra]; c=[:analysis]) = ")
println(d_separation(dag, [:statistics], [:algebra]; c=[:analysis]))

print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors]) = ")
println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors]))

print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors]; c=[:algebra]) = ")
println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors]; c=[:algebra]))

adjustmentsets = adjustment_sets(dag, :statistics, :mechanics)
println("\nAdjustment sets:")
adjustmentsets |> display

adjustmentsets = adjustment_sets(dag, :mechanics, :statistics)
println("\nAdjustment sets:")
adjustmentsets |> display

#end


text:
const _vars_dict = Dict(:vars => Symbol[])

"""
    AffineExpression{T}

Data structure to represent affine expressions, such as ``x+2y+z+4``.

### Examples

```jldoctest
julia> @linvars x y z
3-element Vector{AffineExpression{Int64}}:
 x
 y
 z

julia> p1 = x + 2y + z + 4
x+2y+z+4
```
"""
struct AffineExpression{T}
    coeffs::Vector{T}
end

function AffineExpression(x::T) where {T<:Number}
    c = zeros(promote_type(T, Int), length(_vars_dict[:vars]) + 1)
    c[end] = x
    return AffineExpression(c)
end

## VARIABLES CONSTRUCTION

_get_vars(x::Symbol) = [x]

function _get_vars(x...)

    if length(x) == 1
        s = x[1].args[1]
        start = x[1].args[2].args[2]
        stop = x[1].args[2].args[3]
        return [Symbol(s, i) for i in start:stop]
    else
        return collect(x)
    end
end

"""
    @linvars(x...)

Macro to construct the variables used to represent [`AffineExpression`](@ref).

### Examples

```jldoctest
julia> @linvars x
1-element Vector{AffineExpression{Int64}}:
 x

julia> @linvars x y z
3-element Vector{AffineExpression{Int64}}:
 x
 y
 z

julia> @linvars x[1:4]
4-element Vector{AffineExpression{Int64}}:
 x1
 x2
 x3
 x4
```
"""
macro linvars(x...)
    vars = _get_vars(x...)
    _vars_dict[:vars] = vars

    ex = quote end
    vars_ex = Expr(:vect)
    for (i, s) in enumerate(vars)
        c = zeros(Int, length(vars) + 1)
        c[i] = 1
        push!(ex.args, :($(esc(s)) = AffineExpression($c)))
        push!(vars_ex.args, s)
    end
    push!(ex.args, esc(vars_ex))
    return ex
end

(ae::AffineExpression)(p::Vector{<:Number}) = dot(ae.coeffs[1:end-1], p) + ae.coeffs[end]

## printing
function _tostring(ae::AffineExpression)
    iszero(ae.coeffs) && return "0"
    str = ""
    @inbounds for (i, x) in enumerate(_vars_dict[:vars])
        c = ae.coeffs[i]
        iszero(c) && continue
        sgn = c > 0 ? "+" : "-"
        c = abs(c) == 1 ? "" : "$(abs(c))"
        str *= sgn * c * "$(x)"
    end

    c = last(ae.coeffs)
    if !iszero(c)
        sgn = c > 0 ? "+" : ""
        str *= sgn * "$(c)"
    end
    str = (str[1] == '+' ? str[2:end] : str)
    return str
end

show(io::IO, ae::AffineExpression) = print(io, _tostring(ae))

#########################
# BASIC FUNCTIONS       #
#########################

function zero(::AffineExpression{T}) where {T}
    return AffineExpression(zeros(T, length(_vars_dict[:vars]) + 1))
end

function zero(::Type{AffineExpression{T}}) where {T}
    return AffineExpression(zeros(T, length(_vars_dict[:vars]) + 1))
end

one(::Type{AffineExpression{T}}) where {T} = zero(AffineExpression{T}) + 1
one(::AffineExpression{T}) where {T} = zero(AffineExpression{T}) + 1

# coefficients(ae::AffineExpression) = ae.coeffs[1:end-1]
# coefficient(ae::AffineExpression, i) = ae.coeffs[i]
# offset(ae::AffineExpression) = ae.coeffs[end]

#########################
# ARITHMETIC OPERATIONS #
#########################

for op in (:+, :-)
    @eval function $op(ae1::AffineExpression, ae2::AffineExpression)
        return AffineExpression($op(ae1.coeffs, ae2.coeffs))
    end

    @eval function $op(ae::AffineExpression{T}, n::S) where {T<:Number, S<:Number}
        TS = promote_type(T, S)
        c = similar(ae.coeffs, TS)
        c .= ae.coeffs
        c[end] = $op(c[end], n)
        return AffineExpression(c)
    end
end

+(ae::AffineExpression) = ae
-(ae::AffineExpression) = AffineExpression(-ae.coeffs)

function Base.:(*)(ae1::AffineExpression, n::Number)
    return AffineExpression(ae1.coeffs * n)
end

function Base.:(/)(ae1::AffineExpression, n::Number)
    return AffineExpression(ae1.coeffs / n)
end

for op in (:+, :*)
    @eval $op(n::Number, ae::AffineExpression) = $op(ae, n)
end

-(n::Number, ae::AffineExpression) = n + (-ae)

==(ae1::AffineExpression, ae2::AffineExpression) = ae1.coeffs == ae2.coeffs
==(ae1::AffineExpression, n::Number) = iszero(ae1.coeffs[1:end-1]) && ae1.coeffs[end] == n

function *(ae::AffineExpression, A::AbstractArray{T, N}) where {T<:Number, N}
    return map(x -> x * ae, A)
end

function *(A::AbstractArray{T, N}, ae::AffineExpression) where {T<:Number, N}
    return map(x -> x * ae, A)
end

## Convetion and promotion

function promote_rule(::Type{AffineExpression{T}}, ::Type{S}) where {T<:Number, S<:Number}
    AffineExpression{promote_type(T, S)}
end

function promote_rule(::Type{AffineExpression{T}}, ::Type{AffineExpression{S}}) where {T<:Number, S<:Number}
    AffineExpression{promote_type(T, S)}
end

convert(::Type{AffineExpression}, x::Number) = AffineExpression(x)
convert(::Type{AffineExpression{T}}, x::Number) where T = AffineExpression(convert(T, x))
convert(::Type{AffineExpression{T}}, ae::AffineExpression) where {T<:Number} = AffineExpression{T}(ae.coeffs)

max_stars_repo_name: lucaferranti/IntervalLinearAlgebra

max_issues_repo_name: lucaferranti/IntervalLinearAlgebra

max_forks_repo_name: lucaferranti/IntervalLinearAlgebra


text:
module TestObjectiveContainer

using Test

import VecMathOptInterface

const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_ObjectiveSense()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    @test MOI.get(o, MOI.ListOfModelAttributesSet()) == []
    @test MOI.supports(o, MOI.ObjectiveSense())
    for val in (MOI.MIN_SENSE, MOI.MAX_SENSE, MOI.FEASIBILITY_SENSE)
        MOI.set(o, MOI.ObjectiveSense(), val)
        @test MOI.get(o, MOI.ObjectiveSense()) == val
    end
    return
end

function test_FEASIBILITY_SENSE_clears_objective()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    MOI.set(o, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    @test MOI.get(o, MOI.ObjectiveFunction{MOI.VariableIndex}()) ≈ x
    MOI.set(o, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test !MOI.is_empty(o)
    @test o.is_function_set == false
    return
end

function _test_basic_objective(F, T)
    o = MOI.Utilities.ObjectiveContainer{T}()
    @test MOI.is_empty(o)
    @test MOI.supports(o, MOI.ObjectiveFunction{F}())
    x = MOI.VariableIndex(1234)
    f = convert(F, x)
    MOI.set(o, MOI.ObjectiveFunction{F}(), f)
    @test MOI.get(o, MOI.ObjectiveFunctionType()) == F
    @test !MOI.is_empty(o)
    @test MOI.get(o, MOI.ListOfModelAttributesSet()) ==
          Any[MOI.ObjectiveFunction{F}()]
    MOI.empty!(o)
    @test MOI.is_empty(o)
    return
end

function test_basic_objective()
    _test_basic_objective(MOI.VariableIndex, Float32)
    _test_basic_objective(MOI.ScalarAffineFunction{Float32}, Float32)
    _test_basic_objective(MOI.ScalarQuadraticFunction{Float32}, Float32)
    return
end

function test_delete_VariableIndex()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    MOI.set(o, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.delete(o, x)
    @test MOI.is_empty(o)
    return
end

function test_delete_ScalarAffineFunction()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    f = convert(MOI.ScalarAffineFunction{Float16}, x)
    MOI.set(o, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.delete(o, x)
    @test MOI.get(o, MOI.ObjectiveFunctionType()) == typeof(f)
    @test !(MOI.get(o, MOI.ObjectiveFunction{typeof(f)}()) ≈ f)
    return
end

function test_delete_ScalarQuadraticFunction()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    f = convert(MOI.ScalarQuadraticFunction{Float16}, x)
    MOI.set(o, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.delete(o, x)
    @test MOI.get(o, MOI.ObjectiveFunctionType()) == typeof(f)
    @test !(MOI.get(o, MOI.ObjectiveFunction{typeof(f)}()) ≈ f)
    return
end

function test_delete_VariableIndex_plural()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    MOI.set(o, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.delete(o, [x])
    @test MOI.is_empty(o)
    return
end

function test_delete_ScalarAffineFunction_plural()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    f = convert(MOI.ScalarAffineFunction{Float16}, x)
    MOI.set(o, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.delete(o, [x])
    @test MOI.get(o, MOI.ObjectiveFunctionType()) == typeof(f)
    @test !(MOI.get(o, MOI.ObjectiveFunction{typeof(f)}()) ≈ f)
    return
end

function test_delete_ScalarQuadraticFunction_plural()
    o = MOI.Utilities.ObjectiveContainer{Float16}()
    x = MOI.VariableIndex(1234)
    f = convert(MOI.ScalarQuadraticFunction{Float16}, x)
    MOI.set(o, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.delete(o, [x])
    @test MOI.get(o, MOI.ObjectiveFunctionType()) == typeof(f)
    @test !(MOI.get(o, MOI.ObjectiveFunction{typeof(f)}()) ≈ f)
    return
end

end  # module

TestObjectiveContainer.runtests()


text:
extract_val(::Type{Val{X}}) where X = X
function rewrite_backend_specific(m::CLMethod, f::F, types, expr) where F
    li = m.parent
    if f == (ifelse) && length(types) == 3 && types[1] == Bool && all(x-> x <: cli.Types || is_fixedsize_array(x), types[2:3])
        _types = (types[2], types[3], types[1])
        sm = LazyMethod(li, cli.cl_select, _types)
        expr.args = [sm, expr.args[3], expr.args[4], expr.args[2]]
        return true, expr
    elseif f == abs && length(types) == 1 && all(t-> t <: cli.Floats, types)
        expr.args[1] = LazyMethod(li, fabs, types)
        return true, expr
    elseif f == mod && length(types) == 2 && all(t-> t <: cli.Ints, types)
        expr.args[1] = LazyMethod(li, %, types)
        return true, expr
    elseif f == rem && length(types) == 2 && all(t-> t <: cli.Ints, types)
        expr.args[1] = LazyMethod(li, %, types)
        return true, expr
    elseif f == muladd && length(types) == 3 && all(t-> t <: cli.Floats, types)

        expr.args[1] = LazyMethod(li, fma, types)
        return true, expr
    elseif Symbol(f) == :synchronize_threads
        expr.args[1] = LazyMethod(li, cli.barrier, (Cuint,))
        expr.args[2] = :CLK_LOCAL_MEM_FENCE
        return true, expr
    elseif Symbol(f) == Symbol("GPUArrays.LocalMemory")
        if length(expr.args) != 5
            error("LocalMemory needs to be used as: LocalMemory(state, T, N), with T & N being constants / static function parameters")
        end
        T = Sugar.rewrite_ast.(m, expr.args[3])
        N = expr.args[4].args[1]
        if !(isa(T, Type)) #&& (isa(N, Val{<:Integer}) || isa(N, Val{NTuple{N, Integer}} where N)))
             error("LocalMemory needs to be used as: LocalMemory(T, N), with T being a type and Val{N} where N <: Union{Integer, Tuple}. Found: $T $N")
        end
        name = gensym("localmem")
        ret = Expr(:local_memory, T, name)
        ret.typ = cli.LocalPointer{T}
        inline = Sugar.InlineNode([Expr(:local_memory_init, T, prod(extract_val(N)), name)], ret)
        return true, inline
    else
        return false, expr
    end
end


text:
function _numargs(f)
    typ = Tuple{Any, Val{:analytic}, Vararg}
    typ2 = Tuple{Any, Type{Val{:analytic}}, Vararg} # This one is required for overloaded types
    typ3 = Tuple{Any, Val{:jac}, Vararg}
    typ4 = Tuple{Any, Type{Val{:jac}}, Vararg} # This one is required for overloaded types
    typ5 = Tuple{Any, Val{:tgrad}, Vararg}
    typ6 = Tuple{Any, Type{Val{:tgrad}}, Vararg} # This one is required for overloaded types
    numparam = maximum([(m.sig<:typ || m.sig<:typ2 || m.sig<:typ3 || m.sig<:typ4 || m.sig<:typ5 || m.sig<:typ6) ? 0 : num_types_in_tuple(m.sig) for m in methods(f)])
    return (numparam-1) #-1 in v0.5 since it adds f as the first parameter
end

#Get the number of parameters of a Tuple type, i.e. the number of fields.

function num_types_in_tuple(sig)
  length(sig.parameters)
end

function num_types_in_tuple(sig::UnionAll)
  length(Base.unwrap_unionall(sig).parameters)
end

function num_vecjac!(
    du,
    f,
    x,
    v,
    cache1 = similar(v),
    cache2 = similar(v);
    compute_f0 = true,
)
    if _numargs(f) != 2
        du .= num_jacvec(f, x, v)
        return du
    end
    compute_f0 && (f(cache1, x))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    vv = reshape(v, size(x))
    for i = 1:length(x)
        x[i] += ϵ
        f(cache2, x)
        x[i] -= ϵ
        du[i] = (((cache2 .- cache1) ./ ϵ)'*vv)[1]
    end
    return du
end

function num_vecjac(f, x, v, f0 = nothing)
    vv = reshape(v, axes(x))
    f0 === nothing ? _f0 = f(x) : _f0 = f0
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    du = similar(x)
    for i = 1:length(x)
        x[i] += ϵ
        f0 = f(x)
        x[i] -= ϵ
        du[i] = (((f0 .- _f0) ./ ϵ)'*vv)[1]
    end
    return vec(du)
end


text:
###############################################################################
#
# Query MODIS ungridded data folder
#
###############################################################################
"""
    query_RAW(dt::AbstractUngriddedData,
              year::Int,
              days::Int) where {FT<:AbstractFloat}
    query_RAW(dt::MOD15A2Hv006LAI,
              year::Int) where {FT<:AbstractFloat}

Prepare parameters (file name and etc) to work on, given
- `dt` [`AbstractUngriddedData`](@ref) type ungridded data type
- `year` AD year using format XXXX
- `days` How many days per cycle

Note that `query_RAW(dt, year)` function returns the follwings
- Folder in on Fluo server that contains files to process
- File name pattern, such as "MOD15A2H.AYYYYDOY.*.hdf"
- File name prefix, such as "MOD15A2H.A"
- Band name to process, such as "Lai_500m"
- Cache folder location to store generated CSV files
- Scaling and offset array [a,b] as in y=ax+b
- Mask of minimum and maximum values
"""
function query_RAW(
            dt::AbstractMODIS500m{FT},
            year::Int,
            days::Int
) where {FT<:AbstractFloat}
    # Querry the file information using data type
    _folder, _naming, _pref, _band, _cache, _mm = query_RAW(dt, year);

    # set time information for the gridding process
    date_start = DateTime("$(year)-01-01");
    date_stop  = DateTime("$(year)-12-31");
    date_dday  = Day(days);
    date_list  = date_start:date_dday:date_stop;

    # generate file list to process
    @info "Processing data for year $(year)";
    @info "Please wait while fetching file names...";
    params = [];
    for layer in eachindex(date_list)
        date_temp = date_list[layer];
        for _date in date_temp:Day(1):date_temp+Day(days-1)
            file_dict = ["YYYY" => lpad(Dates.year(_date)     , 4, "0"),
                         "DOY"  => lpad(Dates.dayofyear(_date), 3, "0")];
            file_temp = reduce(replace, file_dict, init=_naming);
            for _path in glob(file_temp, _folder)
                push!(params, [dt, layer, _path, _pref, _band, _cache, _mm]);
            end
        end
    end

    return params
end

function query_RAW(dt::MOD09A1v006NIRv{FT}, year::Int) where {FT<:AbstractFloat}
    _folder = joinpath("$(MODIS_HOME)/MOD09A1.006/original", string(year));
    _naming = "MOD09A1.AYYYYDOY.*.hdf";
    _prefix = "MOD09A1.A";
    _cache  = joinpath("$(MODIS_HOME)/MOD09A1.006/cache", string(year));
    _band   = ["sur_refl_b01", "sur_refl_b02"];

    return _folder, _naming, _prefix, _band, _cache, FT[-0.01,1.7]
end

function query_RAW(dt::MOD15A2Hv006LAI{FT}, year::Int) where {FT<:AbstractFloat}
    _folder = joinpath("$(MODIS_HOME)/MOD15A2H.006/original", string(year));
    _naming = "MOD15A2H.AYYYYDOY.*.hdf";
    _prefix = "MOD15A2H.A";
    _cache  = joinpath("$(MODIS_HOME)/MOD15A2H.006/cache", string(year));
    _band   = "Lai_500m";

    return _folder, _naming, _prefix, _band, _cache, FT[0,15]
end


text:
s = open("input.txt") do file
    readlines(file)
end

function getGoodCount(s)
    goodCount = 0;
    for line in s
        m = match(r"(\d+)-(\d+) (.): (.*)", line);

        pos1 = parse(Int32, m.captures[1]);
        pos2 = parse(Int32, m.captures[2]);
        letter = Char(m.captures[3][1]);

        password = m.captures[4];

        if (password[pos1] == letter) ⊻ (password[pos2] == letter)
            goodCount += 1;
        end
    end

    return goodCount;
end

println(getGoodCount(s));


text:
export isintegral

################################################################################
#
#  Copy
#
################################################################################

Base.copy(d::nf_elem) = deepcopy(d)

################################################################################
#
#  Basis matrix
#
################################################################################

function basis_matrix(A::Vector{nf_elem}, ::Type{FakeFmpqMat})
  @assert length(A) > 0
  n = length(A)
  d = degree(parent(A[1]))

  M = zero_matrix(FlintZZ, n, d)

  deno = one(FlintZZ)
  dummy = one(FlintZZ)

  for i in 1:n
    deno = lcm(deno, denominator(A[i]))
  end

  for i in 1:n
    elem_to_mat_row!(M, i, dummy, A[i])
    temp_den = divexact(deno, dummy)
    for j in 1:d
      M[i, j] = mul!(M[i, j], M[i, j], temp_den)
    end
  end
  return FakeFmpqMat(M, deno)
end

function basis_matrix(A::Vector{nf_elem})
  @assert length(A) > 0
  n = length(A)
  d = degree(parent(A[1]))

  M = zero_matrix(FlintQQ, n, d)

  for i in 1:n
    for j in 1:d
      M[i, j] = coeff(A[i], j - 1)
    end
  end
  return M
end

################################################################################
#
#  Characteristic polynomial
#
################################################################################

function charpoly(Qx::FmpqPolyRing, a::nf_elem)
  f = charpoly(Qx, representation_matrix(a))
  return f
end

function charpoly(a::nf_elem)
  f = charpoly(parent(parent(a).pol), a)
  return f
end

function charpoly(a::nf_elem, ::FlintRationalField)
  return charpoly(a)
end

function charpoly(Zx::FmpzPolyRing, a::nf_elem)
  f = charpoly(a)
  if !isone(denominator(f))
    throw(error("Element is not integral"))
  end
  return Zx(f)
end

function charpoly(a::nf_elem, Z::FlintIntegerRing)
  return charpoly(PolynomialRing(Z, cached = false)[1], a)
end

################################################################################
#
#  Minimal polynomial
#
################################################################################

@doc Markdown.doc"""
    minpoly(a::nf_elem) -> fmpq_poly

The minimal polynomial of $a$.
"""
function minpoly(Qx::FmpqPolyRing, a::nf_elem)
  f = minpoly(Qx, representation_matrix(a))
  return f
end

function minpoly(a::nf_elem)
  f = minpoly(parent(parent(a).pol), a)
  return f
end

function minpoly(a::nf_elem, ::FlintRationalField)
  return minpoly(a)
end

function minpoly(a::nf_elem, ZZ::FlintIntegerRing)
  return minpoly(PolynomialRing(ZZ, cached = false)[1], a)
end

function minpoly(Zx::FmpzPolyRing, a::nf_elem)
  f = minpoly(a)
  if !isone(denominator(f))
    throw(error("Element is not integral"))
  end
  return Zx(f)
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function sub!(a::nf_elem, b::nf_elem, c::nf_elem)
   ccall((:nf_elem_sub, libantic), Nothing,
         (Ref{nf_elem}, Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}),
         a, b, c, a.parent)
end

function set_den!(a::nf_elem, d::fmpz)
  ccall((:nf_elem_set_den, libantic), Nothing,
        (Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
        a, d, parent(a))
end

function set_num_coeff!(a::nf_elem, i::Int, b::fmpz)
  ccall((:_nf_elem_set_coeff_num_fmpz, libantic), Nothing,
        (Ref{nf_elem}, Int, Ref{fmpz}, Ref{AnticNumberField}),
        a, i, b, parent(a))
end

function divexact!(z::nf_elem, x::nf_elem, y::fmpz)
  ccall((:nf_elem_scalar_div_fmpz, libantic), Nothing,
        (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
        z, x, y, parent(x))
  return z
end

function gen!(r::nf_elem)
   a = parent(r)
   ccall((:nf_elem_gen, libantic), Nothing,
         (Ref{nf_elem}, Ref{AnticNumberField}), r, a)
   return r
end

function one!(r::nf_elem)
   a = parent(r)
   ccall((:nf_elem_one, libantic), Nothing,
         (Ref{nf_elem}, Ref{AnticNumberField}), r, a)
   return r
end

function one(r::nf_elem)
   a = parent(r)
   return one(a)
end

function zero(r::nf_elem)
   return zero(parent(r))
end

################################################################################
#
#  Norm div
#
################################################################################

@doc Markdown.doc"""
    norm_div(a::nf_elem, d::fmpz, nb::Int) -> fmpq

Computes `divexact(norm(a), d)` provided the result has at most `nb` bits.

Typically, `a` is an element of some ideal with norm `d`.
"""
function norm_div(a::nf_elem, d::fmpz, nb::Int)
   z = fmpq()
   # TODO:
   #CF the resultant code has trouble with denominators,
   #   this "solves" the problem, but it should probably be
   #   adressed in c
   @assert nb > 0
   n = degree(parent(a))
   if n > 20
     p = p_start
     pp = fmpz(1)
     no = 1
     while nbits(pp) < nb
       p = next_prime(p)
       R = GF(Int(p), cached = false)
       Rt, t = PolynomialRing(R, cached = false)
       np = R(divexact(resultant(Rt(parent(a).pol), Rt(a)), R(d)))
       if isone(pp)
         no = lift(np)
         pp = fmpz(p)
       else
         no = crt(no, pp, lift(np), fmpz(p))
         pp *= fmpz(p)
       end
     end
     no = mod_sym(no, pp)
     return no//1
   end
   de = denominator(a)
   ccall((:nf_elem_norm_div, libantic), Nothing,
         (Ref{fmpq}, Ref{nf_elem}, Ref{AnticNumberField}, Ref{fmpz}, UInt),
         z, (a*de), a.parent, (d*de^n), UInt(nb))
   return z
end

################################################################################
#
#  Is norm divisible by
#
################################################################################

# TODO: Use fits(Int, n) and then split into fmpz_mod/nmod case
@doc Markdown.doc"""
    isnorm_divisible(a::nf_elem, n::fmpz) -> Bool

Checks if the norm of $a$ is divisible by $n$, assuming that the norm of $a$ is
an integer.
"""
function isnorm_divisible(a::nf_elem, n::fmpz)
  K = parent(a)
  if !iscoprime(denominator(K.pol), n)
    na = norm(a)
    @assert isone(denominator(na))
    return divides(numerator(na), n)[1]
  end
  s, t = ppio(denominator(a), n)
  if !isone(s)
    m = n*s^degree(K)
  else
    m = n
  end
  if fits(Int, m)
    R1 = ResidueRing(FlintZZ, Int(m), cached = false)
    R1x = PolynomialRing(R1, "x", cached = false)[1]
    el = resultant_ideal(R1x(numerator(a)), R1x(K.pol))
    return iszero(el)
  end
  R = ResidueRing(FlintZZ, m, cached = false)
  Rx = PolynomialRing(R, "x", cached = false)[1]
  el = resultant_ideal(Rx(numerator(a)), Rx(K.pol))
  return iszero(el)
end

#In this version, n is supposed to be a prime power
function isnorm_divisible_pp(a::nf_elem, n::fmpz)
  K = parent(a)
  if !iscoprime(denominator(K.pol), n)
    na = norm(a)
    @assert isone(denominator(na))
    return divides(numerator(na), n)[1]
  end
  s, t = ppio(denominator(a), n)
  if !isone(s)
    m = n*s^degree(K)
  else
    m = n
  end
  if fits(Int, m)
    R1 = ResidueRing(FlintZZ, Int(m), cached = false)
    R1x = PolynomialRing(R1, "x", cached = false)[1]
    el = resultant_ideal_pp(R1x(numerator(a)), R1x(K.pol))
    return iszero(el)
  end
  R = ResidueRing(FlintZZ, m, cached = false)
  Rx = PolynomialRing(R, "x", cached = false)[1]
  el = resultant_ideal_pp(Rx(numerator(a)), Rx(K.pol))
  return iszero(el)
end

################################################################################
#
#  Numerator
#
################################################################################

@doc Markdown.doc"""
    numerator(a::nf_elem) -> nf_elem

For an element $a\in K = Q[t]/f$ write $a$ as $b/d$ with
$b\in Z[t]$, $\deg(a) = \deg(b)$ and $d>0$ minimal in $Z$.
This function returns $b$.
"""
function numerator(a::nf_elem)
   _one = one(FlintZZ)
   z = deepcopy(a)
   ccall((:nf_elem_set_den, libantic), Nothing,
         (Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
         z, _one, a.parent)
   return z
end

################################################################################
#
#  Chinese remaindering
#
################################################################################

function induce_crt(a::nf_elem, p::fmpz, b::nf_elem, q::fmpz, signed::Bool = false)
  c = parent(a)()
  pi = invmod(p, q)
  mul!(pi, pi, p)
  pq = p*q
  if signed
    pq2 = div(pq, 2)
  else
    pq2 = fmpz(0)
  end
  return induce_inner_crt(a, b, pi, pq, pq2), pq
end

function inner_crt(a::fmpz, b::fmpz, up::fmpz, pq::fmpz, pq2::fmpz = fmpz(0))
  #1 = gcd(p, q) = up + vq
  # then u = modinv(p, q)
  # vq = 1-up. i is up here
  #crt: x = a (p), x = b(q) => x = avq + bup = a(1-up) + bup
  #                              = (b-a)up + a
  if !iszero(pq2)
    r = mod(((b-a)*up + a), pq)
    if r > pq2
      return r-pq
    else
      return r
    end
  else
    return mod(((b-a)*up + a), pq)
  end
end

function induce_inner_crt(a::nf_elem, b::nf_elem, pi::fmpz, pq::fmpz, pq2::fmpz = fmpz(0))
  c = parent(a)()
  ca = fmpz()
  cb = fmpz()
  for i=0:degree(parent(a))-1
    Nemo.num_coeff!(ca, a, i)
    Nemo.num_coeff!(cb, b, i)
    Hecke._num_setcoeff!(c, i, inner_crt(ca, cb, pi, pq, pq2))
  end
  return c
end

################################################################################
#
#  Norm
#
################################################################################

@doc Markdown.doc"""
    norm(f::PolyElem{nf_elem}) -> fmpq_poly

>The norm of $f$, that is, the product of all conjugates of $f$ taken
>coefficientwise.
"""
function norm(f::PolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(f)
  f, i = deflate(f)
  if degree(f) == 1 && ismonic(f)
    N = charpoly(-constant_coefficient(f))
  elseif degree(f) > 10 # TODO: find a good cross-over,
                         # do this using CRT modular?
    P = polynomial_to_power_sums(f, degree(f)*degree(K))
    PQ = fmpq[tr(x) for x in P]
    N = power_sums_to_polynomial(PQ)
  else
    Qx = PolynomialRing(FlintQQ, "x", cached = false)[1]
    Qxy = PolynomialRing(Qx, "y", cached = false)[1]

    T = change_base_ring(Qx, K.pol, parent = Qxy)
    h = nf_poly_to_xy(f, Qxy, Qx)
    N = resultant(T, h)
  end
  return inflate(N, i)
end

################################################################################
#
#  Factorization
#
################################################################################

@doc Markdown.doc"""
    factor(f::fmpz_poly, K::NumberField) -> Fac{Generic.Poly{nf_elem}}
    factor(f::fmpq_poly, K::NumberField) -> Fac{Generic.Poly{nf_elem}}

The factorisation of $f$ over $K$.
"""
function factor(f::fmpq_poly, K::AnticNumberField)
  f1 = change_base_ring(K, f)
  return factor(f1)
end

function factor(f::fmpz_poly, K::AnticNumberField)
  f1 = change_base_ring(K, f)
  return factor(f1)
end

function nice(f::PolyElem{nf_elem})
  if degree(f) < 10
    return "$f"
  end
  if ismonic(f)
    return "$(gen(parent(f))^degree(f)) + ... "
  else
    return "$(leading_coefficient(f))$(gen(parent(f))^degree(f)) + ... "
  end
end

@doc Markdown.doc"""
    factor(f::PolyElem{nf_elem}) -> Fac{Generic.Poly{nf_elem}}

The factorisation of $f$.
"""
function factor(f::PolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(f)

  iszero(f) && error("poly is zero")

  if degree(f) == 0
    r = Fac{typeof(f)}()
    r.fac = Dict{typeof(f), Int}()
    r.unit = Kx(leading_coefficient(f))
    return r
  end
  sqf = factor_squarefree(f)
  fac = Dict{typeof(f), Int}()
  for (k, v) in sqf
    if degree(k) == 1
      fac[k] = v
      continue
    end
    el = k
    if iszero(coeff(k, 0))
      el = shift_right(el, 1)
      fac[gen(Kx)] = v
    end
    @vprint :PolyFactor 1 "Factoring $(nice(el))\n"
    lf = _factor(el)
    for g in lf
      fac[g] = v
    end
  end
  r = Fac{typeof(f)}()
  r.fac = fac
  #The unit is just the leading coefficient of f
  r.unit = Kx(leading_coefficient(f))
  return r
end

  #assumes that f is a squarefree polynomial
function _factor(f::PolyElem{nf_elem})

  K = base_ring(f)
  f = f*(1//leading_coefficient(f))

  if degree(f) < degree(K)
    lf = factor_trager(f)::Vector{typeof(f)}
  else
    lf = factor_new(f)::Vector{typeof(f)}
  end
return lf
end

function factor_trager(f::PolyElem{nf_elem})
  k = 0
  g = f
  @vprint :PolyFactor 1 "Using Trager's method\n"
  p = p_start
  F = GF(p)

  Kx = parent(f)
  K = base_ring(Kx)

  Zx = Hecke.Globals.Zx
  @vtime :PolyFactor Np = norm_mod(g, p, Zx)
  while isconstant(Np) || !issquarefree(map_coefficients(F, Np))
    k = k + 1
    g = compose(f, gen(Kx) - k*gen(K))
    @vtime :PolyFactor 2 Np = norm_mod(g, p, Zx)
  end

  @vprint :PolyFactor 2 "need to shift by $k, now the norm\n"
  if any(x -> denominator(x) > 1, coefficients(g)) ||
     !isdefining_polynomial_nice(K)
    @vtime :PolyFactor 2 N = Hecke.Globals.Qx(norm(g))
  else
    @vtime :PolyFactor 2 N = norm_mod(g, Zx)
    @hassert :PolyFactor 1 N == Zx(norm(g))
  end

  while isconstant(N) || !issquarefree(N)
    error("should not happen")
    k = k + 1
    g = compose(f, gen(Kx) - k*gen(K))
    @vtime :PolyFactor 2 N = norm_mod(g)
  end
  @vtime :PolyFactor 2 fac = factor(N)

  res = typeof(f)[]

  for i in keys(fac.fac)
    t = change_base_ring(K, i, parent = Kx)
    t = compose(t, gen(Kx) + k*gen(K))
    @vtime :PolyFactor 2 t = gcd(f, t)
    push!(res, t)
  end

  return res
end

function isirreducible(f::PolyElem{nf_elem})
  isresult_right, result = isirreducible_easy(f)
  if isresult_right
    return result
  end
  fac = _factor(f)
  return length(fac) == 1
end

function isirreducible_easy(f::PolyElem{nf_elem})
  if degree(f) == 1
    return true, true
  end
  if iszero(coeff(f, 0))
    return true, false
  end
  if !issquarefree(f)
    return true, false
  end

  s = Set{Int}()
  i = 1
  for p = PrimesSet(degree(f)+1, -1)
    try
      d = _degset(f, p)
      if length(s) == 0
        s = d
      else
        s = Base.intersect(s, d)
      end
      if length(s) == 1
        return true, true
      end
      i += 1
      if i > 2*degree(f)
        break
      end
    catch e
      if !isa(e, BadPrime)
        rethrow(e)
      end
    end
  end
  return false, false
end

function _ds(fa)
  @assert all(x->x == 1, values(fa.fac))
  T = Int[degree(x) for x = keys(fa.fac)]
  M = MSet(T)
  return Set(sum(s) for s = subsets(M) if length(s) > 0)
end

function _degset(f::fmpz_poly, p::Int)
  F = GF(p, cached = false)
  Ft, t = PolynomialRing(F, cached = false)
  @assert issquarefree(Ft(f))
  g = Ft(f)
  if !issquarefree(g)
    throw(BadPrime(p))
  end
  fa = factor(g)
  return _ds(fa)
end

function (R::GFPPolyRing)(f::fq_nmod_poly)
  g = R()
  for i=0:degree(f)
    setcoeff!(g, i, coeff(coeff(f, i), 0))
  end
  return g
end

function _degset(f::PolyElem{nf_elem}, p::Int, normal::Bool = false)
  K = base_ring(f)

  me = modular_init(K, p, deg_limit = 1)
  #to be competitive, we need to have Fp, not Fq of degree 1
  if isempty(me)
    return Set(1:degree(f))
  end
  fp = modular_proj(f, me)
  R = GF(p, cached = false)
  Rt = PolynomialRing(R, cached = false)[1]
  if !issquarefree(fp[1])
    throw(BadPrime(p))
  end

  s = _ds(factor(Rt(fp[1])))
  if normal
    return s
  end
  for i=2:length(fp)
    if !issquarefree(fp[i])
      throw(BadPrime(p))
    end
    s = Base.intersect(s, _ds(factor(Rt(fp[i]))))
    if length(s) == 1
      return s
    end
  end
  return s
end

################################################################################
#
#  Roots
#
################################################################################

@doc Markdown.doc"""
    roots(f::fmpz_poly, K::AnticNumberField) -> Vector{nf_elem}

Computes all roots in $K$ of a polynomial $f$. It is assumed that $f$ is non-zero,
squarefree and monic.
"""
function roots(f::fmpz_poly, K::AnticNumberField; kw...)
  f1 = change_base_ring(K, f)
  return roots(f1; kw...)
end

@doc Markdown.doc"""
    roots(f::fmpq_poly, K::AnticNumberField) -> Vector{nf_elem}

Computes all roots in $K$ of a polynomial $f$. It is assumed that $f$ is non-zero,
squarefree and monic.
"""
function roots(f::fmpq_poly, K::AnticNumberField; kw...)
  f1 = change_base_ring(K, f)
  return roots(f1; kw...)
end

function elem_in_nf(a::nf_elem)
  return a
end

@doc Markdown.doc"""
    roots(f::Generic.Poly{nf_elem}; max_roots = degree(f),
                                    ispure = false,
                                    isnormal = false)       -> Vector{nf_elem}

Computes the roots of a polynomial $f$. It is assumed that $f$ is non-zero,
squarefree and monic.

- `max_roots` controls the maximal number of roots the function returns.
- `ispure` indicates whether $f$ is of the form $x^d + c$, where $d$ is the
  degree and $c$ the constant coefficient of $f$.
- `isnormal` indicates that the field contains no or all the roots of $f$.
"""
function roots(f::Generic.Poly{nf_elem}; max_roots::Int = degree(f),
                                         ispure::Bool = false,
                                         isnormal::Bool = false)

  iszero(f) && error("Polynomial must be non-zero")

  if max_roots == 0
    return nf_elem[]
  end

  if max_roots <= 1 && iszero(coeff(f, 0))
    return nf_elem[zero(base_ring(f))]
  end

  if degree(f) == 0
    return nf_elem[]
  end

  if degree(f) == 1
    return nf_elem[-coeff(f, 0)//coeff(f, 1)]
  end

  return _roots_hensel(f, max_roots = max_roots, ispure = ispure, isnormal = isnormal)
end

@doc Markdown.doc"""
    hasroot(f::PolyElem{nf_elem}) -> Bool, nf_elem

Tests if $f$ has a root and return it.
"""
function hasroot(f::PolyElem{nf_elem})
  rt = roots(f, max_roots = 1)
  if length(rt) == 0
    return false, zero(base_ring(f))
  else
    return true, rt[1]
  end
end

@doc Markdown.doc"""
    hasroot(f::fmpz_poly, K::AnticNumberField) -> Bool, nf_elem
    hasroot(f::fmpq_poly, K::AnticNumberField) -> Bool, nf_elem

Tests if $f$ has a root in $K$, and return it.
"""
function hasroot(f::fmpz_poly, K::AnticNumberField)
  f1 = change_base_ring(K, f)
  return hasroot(f1)
end

function hasroot(f::fmpq_poly, K::AnticNumberField)
  f1 = change_base_ring(K, f)
  return hasroot(f1)
end

################################################################################
#
#  Roots
#
################################################################################

@doc Markdown.doc"""
    ispower(a::nf_elem, n::Int; with_roots_unity::Bool = false) -> Bool, nf_elem

Determines whether $a$ has an $n$-th root. If this is the case,
the root is returned.

If the field $K$ is known to contain the $n$-th roots of unity,
one can set `with_roots_unity` to `true`.
"""
function ispower(a::nf_elem, n::Int; with_roots_unity::Bool = false, isintegral::Bool = false, trager = false)
#  @req isdefining_polynomial_nice(parent(a)) "Defining polynomial must be integral and monic"
  @assert n > 0
  if n == 1
    return true, a
  end
  if iszero(a)
    return true, a
  end

  if trager
    return ispower_trager(a, n)
  end

  K = parent(a)
  if isintegral
    d = fmpz(1)
  else
    if ismaximal_order_known(K)
      OK = maximal_order(K)
      d = denominator(a, OK)
    else
      d = denominator(a)
    end
  end
  Ky, y = PolynomialRing(K, "y", cached = false)

  if n == 2 || with_roots_unity
    rt = roots(y^n - a*d^n, max_roots = 1, ispure = true, isnormal = true)
  else
    rt = roots(y^n - a*d^n, max_roots = 1, ispure = true)
  end

  if length(rt) > 0
    return true, rt[1]//d
  else
    return false, zero(a)
  end
end

function ispower_trager(a::nf_elem, n::Int)
  # This is done using Trager factorization, but we can do some short cuts
  # The norm will be the minpoly_a(x^n), which will always be squarefree.
  K = parent(a)
  @vprint :PolyFactor 1 "Computing the minpoly\n"
  @vtime :PolyFactor 1 f = minpoly(a)
  b = K(1)
  c = a*b
  if degree(f) < degree(K)
    i = 0
    while true
      @vprint :PolyFactor 1 "Need to shift it\n"
      b = (gen(K)+i)
      c = a*b^n
      f = minpoly(c)
      if degree(f) == degree(K)
        break
      end
      i += 1
    end
  end
  Qx = parent(f)
  x = gen(Qx)
  N = inflate(f, n)
  @vprint :PolyFactor 1 "Factoring the minpoly\n"
  @vtime :PolyFactor 1 fac = factor(N)
  Kt, t = PolynomialRing(K, "a", cached = false)
  for (p, _) in fac
    if degree(p) == degree(f)
      @vprint :PolyFactor 1 "Computing final gcd\n"
      t = gcd(change_base_ring(K, p, parent = Kt), t^n - c)
      @assert degree(t) == 1
      return true, -divexact(coeff(t, 0), coeff(t, 1))//b
    end
  end

  return false, a
end

function _height(a::nf_elem)
  h = fmpz(1)
  for i in 1:degree(parent(a))
    h = max(h, height(coeff(a, i - 1)))
  end
  return h
end

issquare(a::nf_elem) = ispower(a, 2)

issquare_with_sqrt(a::NumFieldElem) = issquare(a)

sqrt(a::nf_elem) = root(a, 2)

function root(a::nf_elem, n::Int)
  fl, rt = ispower(a, n)
  if fl
    return rt
  end

  error("$a has no $n-th root")
end

function roots(a::nf_elem, n::Int)
  @assert n > 0
  if n == 1 || iszero(a)
    return nf_elem[a]
  end

  d = denominator(a)
  Ky, y = PolynomialRing(parent(a), "y", cached = false)
  rt = roots(y^n - a*d^n, ispure = true)

  return nf_elem[x//d for x = rt]
end

function root(a::NfOrdElem, n::Int)
  fl, rt = ispower(a.elem_in_nf, n)
  if fl
    O = parent(a)
    if denominator(rt, O) == 1
      return O(rt)
    end
  end

  error("$a has no $n-th root")
end

################################################################################
#
#  Inversion via lifting
#
################################################################################

function inv_lift_recon(a::nf_elem)  # not competitive....reconstruction is too slow
  p = next_prime(2^60)
  K = parent(a)
  me = modular_init(K, p)
  ap = Hecke.modular_proj(a, me)
  bp = Hecke.modular_lift([inv(x) for x = ap], me)
  pp = fmpz(p)

  fl, b = Hecke.rational_reconstruction(bp, pp)
  t = K()
  while !fl
#    @assert mod_sym(a*bp - 1, pp) == 0
    mul!(pp, pp, pp)
    mul!(t, a, bp)
    rem!(bp, pp)
    sub!(t, 2, t)
    mul!(bp, bp, t)
    rem!(bp, pp)
#    @assert mod_sym(a*bp - 1, pp) == 0
    fl, b = rational_reconstruction(bp, pp)
    if fl
      if b*a == 1
        return b
      end
      fl = false
    end
  end
  return b
end

function inv_lift(a::nf_elem)  # better, but not enough
  p = next_prime(2^60)
  K = parent(a)
  me = modular_init(K, p)
  ap = modular_proj(a, me)
  bp = modular_lift([inv(x) for x = ap], me)
  pp = fmpz(p)
  fl, b = Hecke.rational_reconstruction(bp, pp)
  t = K()
  n = norm(a)
  while !fl
    Hecke.mul!(t, a, bp)
    Hecke.mul!(pp, pp, pp)
    rem!(t, pp)
    Hecke.sub!(t, 2, t)
    Hecke.mul!(bp, bp, t)
    rem!(t, pp)
    mul!(t, bp, n)
    mod_sym!(t, pp)
    if t*a == n
      return t//n
    end
  end
  return b
end

################################################################################
#
#  Mod function
#
################################################################################

function __mod(a::nf_elem, b::fmpz, fl::Bool = true)#, sym::Bool = false) # Not yet
  z = parent(a)()
  ccall((:nf_elem_mod_fmpz_den, libantic), Nothing, (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}, Cint), z, a, b, parent(a), Cint(fl))
  return z
end

function coprime_denominator(a::nf_elem, b::fmpz)
  z = parent(a)()
  ccall((:nf_elem_coprime_den, libantic), Nothing, (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}), z, a, b, parent(a))
  return z
end

function mod_sym!(a::nf_elem, b::fmpz)
  ccall((:nf_elem_smod_fmpz, libantic), Nothing,
        (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
        a, a, b, parent(a))
  return a
end

function mod(b::nf_elem, p::fmpz)
  K = parent(b)
  if isdefining_polynomial_nice(parent(b))
    return coprime_denominator(b, p)
  else
    m = lcm([p, denominator(K.pol), numerator(coeff(K.pol, degree(K.pol)))])
    return coprime_denominator(b, m)
  end
end

mod(x::nf_elem, y::Integer) = mod(x, fmpz(y))

#Assuming that the denominator of a is one, reduces all the coefficients modulo p
# non-symmetric (positive) residue system
function mod!(a::nf_elem, b::fmpz)
  ccall((:nf_elem_mod_fmpz, libantic), Nothing,
        (Ref{nf_elem}, Ref{nf_elem}, Ref{fmpz}, Ref{AnticNumberField}),
        a, a, b, parent(a))
  return a
end

function rem(a::nf_elem, b::fmpz)
  c = deepcopy(a)
  rem!(c, b)
  return c
end

function mod_sym(a::nf_elem, b::fmpz, b2::fmpz)
  return mod_sym(a, b)
  return z
end

function mod_sym(a::nf_elem, b::fmpz)
  c = deepcopy(a)
  mod_sym!(c, b)
  return c
end

################################################################################
#
#  Conversion to nmod_poly and fmpz_mod_poly
#
################################################################################

function nf_elem_to_nmod_poly!(r::nmod_poly, a::nf_elem, useden::Bool = true)
  ccall((:nf_elem_get_nmod_poly_den, libantic), Nothing,
        (Ref{nmod_poly}, Ref{nf_elem}, Ref{AnticNumberField}, Cint),
        r, a, a.parent, Cint(useden))
  return nothing
end

function nf_elem_to_gfp_poly!(r::gfp_poly, a::nf_elem, useden::Bool = true)
  ccall((:nf_elem_get_nmod_poly_den, libantic), Nothing,
        (Ref{gfp_poly}, Ref{nf_elem}, Ref{AnticNumberField}, Cint),
        r, a, a.parent, Cint(useden))
  return nothing
end

function (R::Nemo.NmodPolyRing)(a::nf_elem)
  r = R()
  nf_elem_to_nmod_poly!(r, a)
  return r
end

function (R::Nemo.GFPPolyRing)(a::nf_elem)
  r = R()
  nf_elem_to_gfp_poly!(r, a)
  return r
end

function nf_elem_to_fmpz_mod_poly!(r::fmpz_mod_poly, a::nf_elem, useden::Bool = true)
  ccall((:nf_elem_get_fmpz_mod_poly_den, libantic), Nothing,
        (Ref{fmpz_mod_poly}, Ref{nf_elem}, Ref{AnticNumberField}, Cint, Ref{fmpz_mod_ctx_struct}),
        r, a, a.parent, Cint(useden), r.parent.base_ring.ninv)
  return nothing
end

function nf_elem_to_gfp_fmpz_poly!(r::gfp_fmpz_poly, a::nf_elem, useden::Bool = true)
  ccall((:nf_elem_get_fmpz_mod_poly_den, libantic), Nothing,
        (Ref{gfp_fmpz_poly}, Ref{nf_elem}, Ref{AnticNumberField}, Cint, Ref{fmpz_mod_ctx_struct}),
        r, a, a.parent, Cint(useden), r.parent.base_ring.ninv)
  return nothing
end

function (R::Nemo.FmpzModPolyRing)(a::nf_elem)
  r = R()
  nf_elem_to_fmpz_mod_poly!(r, a)
  return r
end

################################################################################
#
#  Complex conjugate
#
################################################################################

function conjugate_quad(a::nf_elem)
  k = parent(a)
  @assert degree(k) == 2
  #fallback: tr(a) = a + bar(a), so tr(a) - a = bar(a)...
  #in the easy case: tr(gen(k)) = -r
  #if there are dens then it gets messy and I am not sure it is worth while:
  # x + y gen(k) -> x + y bar(k) and bar(k) = tr(k) - gen(k), but
  # tr(k) = -b/a (for the polyomial ax^2 + bx + c), hence:
  # (x+y gen(k)) / d -> (ax - by - ay gen(k))/(ad)
  # and there we might have to do simplification.
  #TODO: on 2nd thought: we might have to simplify in the easy case as well?
  isone(k.pol_den) || return tr(a) - a
  # we have
  # a = x + y gen(k), so bar(a) = x + y bar(k)
  # assume pol(k) is monic: x^2 + rx + t, then
  # bar(gen(k)) = -gen(k) - r
  # since (pq-formula) gen(k) = - r/2 \pm 1/2(sqrt(r^2-4t)
  # so bar(a) = x + y (-bar(k) - r) = (x-ry) - y gen(k)
  b = k()
  q = fmpz()
  GC.@preserve b begin
    a_ptr = reinterpret(Ptr{Int}, pointer_from_objref(a))
    p_ptr = reinterpret(Ptr{Int}, k.pol_coeffs)
    s = sizeof(Ptr{fmpz})

    ccall((:fmpz_mul, libflint), Cvoid, (Ref{fmpz}, Ptr{Int}, Ptr{Int}), q, p_ptr+s, a_ptr +s)
    ccall((:fmpz_sub, libflint), Cvoid, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), reinterpret(Ptr{Int}, pointer_from_objref(b)), a_ptr, q)
    ccall((:fmpz_neg, libflint), Cvoid, (Ptr{fmpz}, Ptr{fmpz}), reinterpret(Ptr{Int}, pointer_from_objref(b))+1*s, a_ptr + s)
    ccall((:fmpz_set, libflint), Cvoid, (Ptr{fmpz}, Ptr{fmpz}), reinterpret(Ptr{Int}, pointer_from_objref(b))+3*s, a_ptr+3*s)
  end
  #TODO:
  # - write in c?
  # - use Ref and Ref(, i) instead of pointers
  # - deal with non-monic fields
  return b
end

function complex_conjugate(a::nf_elem)
  d = degree(parent(a))
  if d == 1
    return a
  end
  if d == 2
    if discriminant(parent(a)) < 0
      return a
    end
    return conjugate_quad(a)
  end
  error("Not implemented yet: element must be in an at most quadratic field")
end


text:
# profile/test conversion utils

using TransformUtils
using Test
using BenchmarkTools, ProfileView

R = SO3(0)
E = Euler(0)

RR = convert(SO3, so3(randn(3)))
E = convert(Euler, RR)

# pvt: 371.464 ms (5000000 allocations: 534.06 MiB)
@btime for i in 1:10^6 convert!(R, E); end
@profiler for i in 1:10^6 convert!(R, E); end

T1 = SE3(zeros(3),R)
T2 = SE3(0)

# pvt: 520.225 ns (20 allocations: 1.39 KiB - norm()
@btime T3 = T1 \ T2

# pvt: 519.044 μs (20000 allocations: 1.36 MiB)
@btime for i in 1:10^3 T3 = T1 \ T2; end
using Profile
Profile.clear()
@profile for i in 1:10^3 T3 = T1 \ T2; end

Juno.profiler()
Juno.profiletree()

@testset "From: quaternion" begin
@test
end

function fast_norm_sonar(u::::Vector{Float64})
  s = 0.0
  @fastmath @inbounds @simd for i in 1:n
      s += u[i]^2
  end
  @fastmath @inbounds return sqrt(s/n)
end

import Base: normalize!

function normalize!(q::Quaternion, tol=1e-6)

  s = q.s^2
  @fastmath @inbounds @simd for i in 1:3
      s += q.v[i]^2
  end
  @fastmath mag = sqrt(s)
  @fastmath q.v ./= mag
  @fastmath q.s /= mag
  nothing
end

@code_warntype normalize!(q)

@code_llvm normalize!(q)

# code_native(normalize!, (Quaternion,))
@code_native normalize!(q)

@btime normalize!(q) # 17ns

@btime for i in 1:1000 normalize!(q); end

@time normalize!(q)

q = convert(Quaternion, so3(randn(3)))
E = convert(Euler, convert(SO3, so3(randn(3))))

@btime convert!(q, E)

#
Profile.clear()
@profiler convert!(q, E)

Profile.clear()
@profiler for i in 1:10^6 convert!(q, E); end

@code_warntype convert!(q, E)

# julia --track-allocation=user testInPlaceConvert.jl