text:
"""
**Create default parameters**

This function creates model parameters, based on a food web
matrix A. See documentation for more information. Specifically, the default values of the keyword parameters are:

| Parameter         | Default Value | Meaning                                                                                |
|:------------------|:--------------|:--------------------------------------------------------------------------------------------|
| K                 | 1.0           | carrying capacity of producers                                                              |
| Z                 | 1.0           | consumer-resource body mass ratio                                                           |
| bodymass          | [0.0]         | default is to calculates masses from Z unless a vector of bodymass is passed                |
| dry_mass_293      | [0.0]         | dry masses at 20 degrees celsius                                                            |
| TSR_type          | :no_response  | temperature size rule: set the function for masses dependence to temperature, see docs      |
| vertebrates       | [false]       | metabolic status of species (invertebrates or ectotherm vertebrates)                        |
| r                 | 1.0           | growth rate of producers                                                                    |
| c                 | 0             | quantifies the predator interference                                                        |
| h                 | 1             | Hill coefficient                                                                            |
| e_carnivore       | 0.85          | assimilation efficiency of carnivores                                                       |
| e_herbivore       | 0.45          | assimilation efficiency of herbivores                                                       |
| y_invertebrate    | 8             | maximum consumption rate of invertebrate predators relative to their metabolic rate         |
| y_vertebrate      | 4             | maximum consumption rate of vertebrate predators relative to their metabolic rate           |
| Γ                 | 0.5           | half-saturation density                                                                     |
| α                 | 1.0           | interspecific competition relatively to intraspecific competition                           |
| scale_bodymass    | true          | whether to normalize body masses by the mass of the smallest producer                       |
| scale_growth      | false         | whether to normalize growth rates by the growth rate of the smallest producer               |
| scale_metabolism  | false         | whether to normalize metabolic rates by the growth rate of the smallest producer            |
| scale_maxcons     | false         | whether to normalize max. consumption rates by metabolic rates                              |
| productivity      | :species      | type of productivity regulation                                                             |
| dc                | x -> x .* 0.0 | density dependent mortality function for consumers                                          |
| dp                | x -> x .* 0.0 | density dependent mortality function for producers                                          |
| rewire_method     | :none         | method for rewiring the foodweb following extinction events                                 |
| adbm_trigger      | :extinction   | (ADBM) trigger for ADBM rewiring (on extinctions or periodic with :interval)                |
| adbm_interval     | 100           | (ADBM) Δt for periodic rewiring                                                             |
| e                 | 1             | (ADBM) Scaling constant for the net energy gain                                             |
| a_adbm            | 0.0189        | (ADBM) Scaling constant for the attack rate                                                 |
| ai                | -0.491        | (ADBM) Consumer specific scaling exponent for the attack rate                               |
| aj                | -0.465        | (ADBM) Resource specific scaling exponent for the attack rate                               |
| b                 | 0.401         | (ADBM) Scaling constant for handling time                                                   |
| h_adbm            | 1.0           | (ADBM) Scaling constant for handling time                                                   |
| hi                | 1.0           | (ADBM) Consumer specific scaling exponent for handling time                                 |
| hj                | 1.0           | (ADBM) Resource specific scaling constant for handling time                                 |
| n                 | 1.0           | (ADBM) Scaling constant for the resource density                                            |
| ni                | 0.75          | (ADBM) Species-specific scaling exponent for the resource density                           |
| Hmethod           | :ratio        | (ADBM) Method used to calculate the handling time                                           |
| Nmethod           | :original     | (ADBM) Method used to calculate the resource density                                        |
| cost              | 0.0           | (Gilljam) Rewiring cost (a consumer decrease in efficiency when exploiting novel resource)  |
| specialistPrefMag | 0.9           | (Gilljam) Strength of the consumer preference for 1 prey if `preferenceMethod = :specialist`|
| preferenceMethod  | :generalist   | (Gilljam) Scenarios with respect to prey preferences of consumers                           |
| D                 | 0.25          | (nutrient intake) global turnover rate                                                      |
| ν                 | [1.0, 0.5]    | (nutrient intake) conversion rates for nutrients                                            |
| K1                | [0.15]        | (nutient intake) species half saturation densities for nutrient 1                           |
| K2                | [0.15]        | (nutient intake) species half saturation densities for nutrient 2                           |
| T                 | 273.15        | (temperature dependence) temperature                                                        |
| growthrate        | NoEffectTemperature(:r) | Function used to calculate growth rate |
| metabolicrate     | NoEffectTemperature(:x) | Function used to calculate metabolic rate |
| handlingtime      | NoEffectTemperature(:handlintime) | Function used to calculate handling time (and max. consumption rate) |
| attackrate        | NoEffectTemperature(:attackrate) | Function used to calculate attack rate (and Γ) |

All of these values are passed as optional keyword arguments to the function.

~~~ julia
A = [0 1 1; 0 0 0; 0 0 0]
p = model_parameters(A, Z=100.0, productivity=:system)
~~~

The `productivity` keyword can be either `:species` (each species has an
independant carrying capacity equal to `K`), `:system` (the carrying capacity is
K divided by the number of primary producers), or `:competitive` (the species
compete with themselves at rate 1.0, and with one another at rate α).

It is possible for the user to specify a vector of species body-mass, called
`bodymass` -- please do pay attention to the fact that the default behavior is
to assume that primary producers have a bodymass equal to unity, since all
biological rates are expressed relatively. We do not perform any check on
whether or not the user-supplied body-mass vector is correct (mostly because
there is no way of defining correctness for vectors where body-mass of producers
are not equal to unity). It is possible to turn off the scaling (`scale_bodymass
= false`), in that case we strongly recommend setting `scale_growth = true` (
and possibly `scale_metabolism = true` and `scale_maxcons = true`).

The keyword `vertebrates` is an array of `true` or `false`
for every species in the matrix. If the vector is of size 1 (e.g. [false]),
all species will take this metabolic status. By default, all species are invertebrates.

A rewiring method can pe passed to specified if the foodweb should be rewired
following extinctions events, and the method that should be used to perform the
rewiring. This `rewire_method` keyword can be eighter `:none` (no rewiring),
`:ADBM` (allometric diet breadth model as described in Petchey
et al., 2008), `:DS` (rewiring mechanism used by Gilljam et al., 2015, based on diet
similarity) or `:DO` (rewiring mechanism used by Staniczenko et al, 2010, based
on diet overlap).

If `rewire_method`is `:ADBM` or `:DS`, additional keywords can be passed.
See the online documentation and the original references for more details.

Temperature dependence for metabolic rates and body masses can be implemented,
please check the documentation for more information.

"""
function model_parameters(A;
        K::Float64=1.0,
        Z::Float64=1.0,
        c::Float64=0.0,
        h::Number=1.0,
        e_carnivore::Float64=0.85,
        e_herbivore::Float64=0.45,
        α::Float64=1.0,
        productivity::Symbol=:species,
        bodymass::Array{Float64, 1}=[0.0],
        scale_bodymass::Bool=true,
        vertebrates::Array{Bool, 1}=[false],
        dc::Function= (x -> x .* 0.0),
        dp::Function= (x -> x .* 0.0),
        rewire_method::Symbol = :none,
        adbm_trigger::Symbol = :extinction,
        adbm_interval::Int64 = 100,
        e::Float64 = 1.0,
        a_adbm::Float64 = 0.0189,
        ai::Float64 = -0.491,
        aj::Float64 = -0.465,
        b::Float64 = 0.401,
        h_adbm::Float64 = 1.0,
        hi::Float64 = 1.0,
        hj::Float64 = 1.0,
        n::Float64 = 1.0,
        ni::Float64= -0.75,
        Hmethod::Symbol = :ratio,
        Nmethod::Symbol = :original,
        cost::Float64 = 0.0,
        specialistPrefMag::Float64 = 0.9,
        preferenceMethod::Symbol = :generalist,
        D::Float64 = 0.25,
        supply::Array{Float64, 1} = [4.0],
        υ::Array{Float64, 1} = [1.0, 0.5],
        K1::Array{Float64, 1} = [0.15],
        K2::Array{Float64, 1} = [0.15],
        T::Float64 = 273.15,
        handlingtime::Function = NoEffectTemperature(:handlingtime),
        attackrate::Function = NoEffectTemperature(:attackrate),
        metabolicrate::Function = NoEffectTemperature(:metabolism),
        growthrate::Function = NoEffectTemperature(:growth),
        scale_growth::Bool = false,
        scale_metabolism::Bool = false,
        scale_maxcons::Bool = false,
        dry_mass_293::Array{Float64, 1}=[0.0],
        TSR_type::Symbol = :no_response)

  check_food_web(A)

  # Step 1 -- create a dictionnary to store the parameters
  parameters = Dict{Symbol,Any}(
  :K              => K,
  :Z              => Z,
  :c              => c,
  :e_carnivore    => e_carnivore,
  :e_herbivore    => e_herbivore,
  :h              => h,
  :vertebrates    => falses(size(A)[1]),
  :A              => A,
  :α              => α,
  :TSR_type       => TSR_type,
  :dry_mass_293   => dry_mass_293,
  :T              => T
  )
  check_initial_parameters(parameters)

  # Step 2 -- vertebrates ?
  if length(vertebrates) > 1
    if length(vertebrates) == size(A, 1)
      parameters[:vertebrates] = vertebrates
    else
      error("when calling `model_parameters` with an array of values for `vertebrates`, there must be as many elements as rows/columns in the matrix")
    end
  end

  # Step 3 -- body mass and dry mass at 293 K
  parameters[:bodymass] = bodymass
  if length(parameters[:bodymass]) > 1
    if length(parameters[:bodymass]) != size(A, 1)
      error("when calling `model_parameters` with an array of values for `bodymass`, there must be as many elements as rows/columns in the matrix")
    end
  end

  parameters[:dry_mass_293] = dry_mass_293
  if length(parameters[:dry_mass_293]) > 1
    parameters[:dry_mass_293] = dry_mass_293
    if length(parameters[:dry_mass_293]) != size(A, 1)
      error("when calling `model_parameters` with an array of values for `dry_mass_293`, there must be as many elements as rows/columns in the matrix")
    end
  end

  # Step 4 -- TSR type
  if TSR_type ∈ [:no_response, :mean_aquatic, :mean_terrestrial, :maximum, :reverse, :no_response]
    parameters[:TSR_type] = TSR_type
  else
    error("Invalid value for TSR_type -- must be :no_response, :mean_aquatic, :mean_terrestrial, :maximum, :reverse, :no_response")
  end

  # Step 5 -- Identify producers
  is_producer = vec(sum(A, dims = 2) .== 0)
  parameters[:is_producer] = is_producer
  producers_richness = sum(is_producer)

  # Step 6 -- productivity type
  if productivity ∈ [:species, :system, :competitive, :nutrients]
    parameters[:productivity] = productivity
  else
    error("Invalid value for productivity -- must be :system, :species, :competitive or :nutrients")
  end

  # step 7 -- productivity parameters for the NP model
  if parameters[:productivity] == :nutrients
    parameters[:D] = D
    parameters[:supply] = supply
    if length(parameters[:supply]) > 1
      if length(parameters[:supply]) != 2
        error("when calling `model_parameters` with an array of values for `S` (nutrient supply), there must be as many elements as nutrients (2)")
      end
    else
      parameters[:supply] = repeat(supply, 2)
    end
    parameters[:υ] = υ
    if length(parameters[:υ]) != 2
        error("when calling `model_parameters` with an array of values for `υ` (conversion rates), there must be as many elements as nutrients (2)")
    end

    parameters[:K1] = K1
    parameters[:K2] = K2
    if length(parameters[:K1]) > 1
      if length(parameters[:K1]) != size(A, 1)
        error("when calling `model_parameters` with an array of values for `K1` (species half-saturation densities for nutrient 1), there must be as many elements as species")
      end
    else
      parameters[:K1] = is_producer .* K1
    end
    if length(parameters[:K2]) > 1
      if length(parameters[:K2]) != size(A, 1)
        error("when calling `model_parameters` with an array of values for `K2` (species half-saturation densities for nutrient 2), there must be as many elements as species")
      end
    else
      parameters[:K2] = is_producer .* repeat(K2, size(A, 1))
    end

  end

  # Step 8 -- rewire method

 if rewire_method ∈ [:DO, :DS, :stan, :none, :ADBM, :Gilljam, :ADBM_interval]
    parameters[:rewire_method] = rewire_method
    parameters[:extinctions] = []
    parameters[:extinctionstime] = []
    parameters[:tmpA] = []
 else
    error("Invalid method for rewiring -- must be :DO (alternatively :stan), :ADBM, :DS (alternatively :Gilljam) or :none")
 end

 if rewire_method == :ADBM
     if adbm_trigger ∈ [:extinction, :interval]
       parameters[:adbm_trigger] = adbm_trigger
       parameters[:adbm_trigger] == :interval ? parameters[:adbm_interval] = adbm_interval : nothing
       adbm_parameters(parameters, e, a_adbm, ai, aj, b, h_adbm, hi, hj, n, ni, Hmethod, Nmethod)
     else
       error("Invalid trigger for ADBM trigger, must be either :interval or :extinction")
     end
 elseif rewire_method ∈ [:Gilljam, :DS]
     gilljam_parameters(parameters, cost, specialistPrefMag, preferenceMethod)
 end
 check_rewiring_parameters(parameters, parameters[:rewire_method])

  # Setup some objects
  S = size(A)[1]
  parameters[:S] = S
  F = zeros(Float64, size(A))
  efficiency = zeros(Float64, size(A))
  M = zeros(Float64, S)
  x = zeros(Float64, S) # metabolic rate
  r = zeros(Float64, S) # producers growth rate
  attack_r = zeros(Float64, S) # attack rates
  handling_t = zeros(Float64, S) # handling times
  Γ = zeros(Float64, S) #B0, half saturation density
  TR = trophic_rank(A)
  parameters[:trophic_rank] = TR
  is_herbivore = falses(S)

  # Step 9 -- Identify herbivores (Herbivores consume producers)
  get_herbivores(parameters)

  # Step 10 -- Measure generality and extract the vector of 1/n
  getW_preference(parameters)

  # Step 11 -- Get the body mass
  temperature_size_rule(parameters)

  # Step 12 -- Growth rate
  m_producer = minimum(parameters[:bodymass][is_producer])
  parameters[:m_producer] = m_producer
  body_size_relative = parameters[:bodymass] ./ parameters[:m_producer]
  m = scale_bodymass ? body_size_relative : bodymass
  r = growthrate(m, T, parameters)
  rspp = r[sortperm(parameters[:bodymass])[1]]
  r_scaled = r ./ rspp
  parameters[:r] = scale_growth ? r_scaled : r

  # Step 13 -- Metabolic rate
  x = metabolicrate(m, T, parameters)
  x_scaled = x ./ rspp
  parameters[:x] = scale_metabolism ? x_scaled : x

  # Step 14 -- Handling time
  handling_t = handlingtime(m, T, parameters)
  parameters[:ht] = handling_t

  # Step 16 -- Maximum relative consumption rate
  y = 1 ./ handling_t
  parameters[:y] = scale_maxcons == true ? y ./ x : y

  # Step 15 -- Attack rate
  attack_r = attackrate(m, T, parameters)

  # Step 17 -- Half-saturation constant
  Γ = 1 ./ (attack_r .* handling_t)
  Γ[isnan.(Γ)] .= 0.0
  parameters[:Γ] = Γ

  # Step 18 -- Efficiency matrix
  get_efficiency(parameters)

  parameters[:Γh] = parameters[:Γ] .^ parameters[:h]
  parameters[:np] = sum(parameters[:is_producer])
  parameters[:ar] = attack_r

  # Step  19 -- Density dependent mortality
  parameters[:dc] = dc
  parameters[:dp] = dp

  check_parameters(parameters)

  return parameters
end

function get_herbivores(parameters::Dict{Symbol,Any})
  #used internally by model_parameters and update_parameters
  S = size(parameters[:A], 1)
  is_herbivore = falses(S)
  # Identify herbivores
  # Herbivores consume producers
  for consumer in 1:S
    if !parameters[:is_producer][consumer]
      for resource in 1:S
        if parameters[:is_producer][resource]
          if parameters[:A][consumer, resource] == 1
            is_herbivore[consumer] = true
          end
        end
      end
    end
  end
  parameters[:is_herbivore] = is_herbivore
  #return(p)
end

function get_efficiency(parameters::Dict{Symbol,Any})
  #used internally by model_parameters and update_parameters
  S = size(parameters[:A], 1)
  # Efficiency matrix
  efficiency = zeros(Float64,(S, S))
  for consumer in 1:S
    for resource in 1:S
      if parameters[:A][consumer, resource] == 1
        if parameters[:is_producer][resource]
          efficiency[consumer, resource] = parameters[:e_herbivore]
        else
          efficiency[consumer, resource] = parameters[:e_carnivore]
        end
      end
    end
  end
  parameters[:efficiency] = efficiency
  #return(p)
end

function getW_preference(parameters::Dict{Symbol,Any})
  #used internally by model_parameters and update_parameters
  S = size(parameters[:A],1)
  generality = float(vec(sum(parameters[:A], dims = 2)))

  if (parameters[:rewire_method] ∈ [:none, :ADBM, :stan, :DO])

    w = zeros(Float64,(S))
    for i in eachindex(generality)
      if generality[i] > 0.0
        w[i] = 1.0 / generality[i]
      end
    end
    w = w .*parameters[:A]

  else

    if parameters[:preferenceMethod] == :generalist
      w = zeros(Float64,(S))
      for i in eachindex(generality)
        if generality[i] > 0.0
          w[i] = 1.0 / generality[i]
        end
      end
      w = w .*parameters[:A]

    elseif parameters[:preferenceMethod] == :specialist
      w = zeros(Float64,(S,S))
      for i in eachindex(generality)
        if generality[i] > 0.0
          for j = 1:S
            if parameters[:A][i,j] == 1
              if parameters[:specialistPref][i] == j
                if generality[i] > 1
                  w[i,j] = parameters[:specialistPrefMag]
                else
                  w[i,j] = 1.0
                end
              else
                if generality[i] > 1
                  w[i,j] = (1 - parameters[:specialistPrefMag])/(generality[i]-1)
                else
                  w[i,j] = 1.0
                end
              end
            end
          end
        end
      end

    end # ifelse parameters[:preferenceMethod]

  end #ifelse parameters[:rewire_method]

  parameters[:w] = w
  #return(p)
end

function get_specialist_preferences(pr, A)
  specials = zeros(Int64,size(A,1))
  if pr[:preferenceMethod] == :specialist
    for pred = 1:size(A,2)
        prey = findall(x -> x!= 0, A[pred,:])
        if length(prey) > 0
          specials[pred] = sample(prey,1)[1]
        end
    end
  end
  return(specials)
end

function preference_parameters(cost, specialistPrefMag, A, preferenceMethod)
  #Work out the jaccard similarity
  S = size(A,1)
  similarity = zeros(Float64,(S,S))
  for i = 1:S, j = 1:S
    if i == j
    similarity[i,j] = 0
    else
      X = findall(x -> x != 0, A[i,:])
      Y = findall(x -> x != 0, A[j,:])
      if length(X) == 0 && length(Y) == 0
        similarity[i,j] = 0
      else
        similarity[i,j] = size(intersect(X,Y),1) / size(union(X,Y),1)
      end
    end
  end

  similarity_indexes = Vector{Vector{Int}}(undef, S)
  #convert to indexes
  for i = 1:S
    similarity_indexes[i] = sortperm(similarity[i,:])
  end

  preference_parameters = Dict{Symbol,Any}(:similarity => similarity_indexes,
                                          :cost       => cost,
                                          :specialistPrefMag => specialistPrefMag,
                                          #:extinctions => Array{Int,1}(),
                                          :costMat => ones(Float64,(S,S)),
                                          :preferenceMethod => preferenceMethod)
  return(preference_parameters)
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
### Sameq

@mkapprule SameQ  nodefault => true
@doap SameQ(args...) = sameq(args...)
## Mma does the following
@doap SameQ(x) = true
@doap SameQ() = true

### UnSameq

@mkapprule UnsameQ nargs => 2
@doap UnsameQ(x,y) = ! sameq(x,y)

sameq(x,y) = x === y

sameq(x::BigInt,y::BigInt) = (x == y)
sameq(x::BigFloat,y::BigFloat) = (x == y)
sameq(x::String,y::String) = (x == y)
sameq(x) = true
sameq() = true

function sameq(args...)
    for i in 1:length(args)-1
        sameq(args[i],args[i+1]) || return false
    end
    true
end

struct Compare
    result
    known::Bool
end

### Equal

## TODO: Implement Equal(a,b,c), etc.
## TODO: Use Union{Nothing, T} here rather than struct
@mkapprule Equal
@doap function Equal(x,y)
    res = sjequal(x,y)
    res.known == false && return mx
    return res.result
end
@doap Equal() = true

function sjequal(x::Symbol,y)
    x == :Undefined && return Compare(:Undefined,true)
    x == y  && return Compare(true,true)
    Compare(false,false)
end

function sjequal(x,y::Symbol)
    y == :Undefined && return Compare(:Undefined,true)
    x == y  && return Compare(true,true)
    Compare(false,false)
end

function sjequal(x::Symbol,y::Symbol)
    (y == :Undefined || x == :Undefined) && return Compare(:Undefined,true)
    x == y  && return Compare(true,true)
    Compare(false,false)
end

sjequal(x::Real,y::Real) = Compare(x == y, true)
sjequal(x::String,y::String) = Compare(x == y, true)

function sjequal(x,y)
    x == y  && return Compare(true,true)
    Compare(false,false)
end

### Unequal

@mkapprule Unequal nargs => 2
@doap function Unequal(x,y)
    res = sjequal(x,y)
    res.known == false && return mx
    res.result == :Undefined && return :Undefined
    !(res.result)
end

### Less

## TODO: generate methods for Undefined using eval. Or find a better, more general solution

## In Mma, Less is nary
@mkapprule Less
@doap function Less(x,y)
    res = sjless(x,y)
    res.known == false && return mx
    res.result
end

## Mma does this
## TODO: use a macro to generate these. or otherwise organize this.
@doap Less() = true
@doap Less(x) = true
#@doap Less(args...) = mx

sjless(x::Real, y::Real) = Compare(x < y, true)
function sjless(x,y)
    Compare(false,false)
end

### LessEqual

## In Mma, Less is nary
@mkapprule LessEqual
@doap function LessEqual(x,y)
    res = sjlessequal(x,y)
    res.known == false && return mx
    res.result
end

sjlessequal(x::Real, y::Real) = Compare(x <= y, true)
function sjlessequal(x,y)
    Compare(false,false)
end

### Greater

@mkapprule Greater
@doap function Greater(x,y)
    res = sjgreater(x,y)
    res.known == false && return mx
    res.result
end

@doap Greater() = true
@doap Greater(x) = true

sjgreater(x::Real, y::Real) = Compare(x > y, true)
function sjgreater(x,y)
    Compare(false,false)
end

### GreaterEqual

@mkapprule GreaterEqual
@doap function GreaterEqual(x,y)
    res = sjgreaterequal(x,y)
    res.known == false && return mx
    res.result
end

sjgreaterequal(x::Real, y::Real) = Compare(x >= y, true)
function sjgreaterequal(x,y)
    Compare(false,false)
end

#### And

function apprules(mx::Mxpr{:And})
    args = margs(mx)
    length(args) == 0 && return true
    nargs = newargs()
    for arg in args
        arg = doeval(arg) # And has attribute HoldAll
        if isa(arg,Bool)
            arg == true && continue
            arg == false && return false
        end
        push!(nargs, arg)
    end
    length(nargs) == 1 && return nargs[1]
    mxpr(:And,nargs)
end

#### Or

function apprules(mx::Mxpr{:Or})
    args = margs(mx)
    length(args) == 0 && return false
    nargs = newargs()
    for arg in args
        arg = doeval(arg) # Or has attribute HoldAll
        if isa(arg,Bool)
            arg == true && return true
            arg == false && continue
        end
        push!(nargs, arg)
    end
    length(nargs) == 1 && return nargs[1]
    mxpr(:Or,nargs)
end

#### Not

@sjdoc Not """
    Not(expr)

return `False` if `expr` is `True`, and `True` if it is `False`.

`Not` reduces some very simple logical expressions and otherwise remains unevaluated. `Not(expr)` may also be entered `! expr`.
"""

@mkapprule Not nargs => 1

@doap Not(ex::Bool) = ex == true ? false : true

const comparison_negations  = Dict(
                               :<   =>  :>=,
                               :>   =>  :<=,
                               :<=  =>  :>,
                               :>=  =>  :<,
                               :(==)  =>  :!=,
                               :!=  =>  :(==)
                               )

function do_Not(mx::Mxpr{:Not},  ex::Mxpr{:Comparison})
    length(ex) == 3 && return mxpr(:Comparison, ex[1], comparison_negations[ex[2]], ex[3])
    return mx
end

for (a,b) in ( (:Less, :GreaterEqual), (:Greater, :LessEqual), (:Equal, :Unequal) )
    @eval begin
        function do_Not(mx::Mxpr{:Not}, ex::Mxpr{$(QuoteNode(a))})
            length(ex) == 2 && return mxpr($(QuoteNode(b)), ex[1], ex[2])
            return mx
        end
        function do_Not(mx::Mxpr{:Not}, ex::Mxpr{$(QuoteNode(b))})
            length(ex) == 2 && return mxpr($(QuoteNode(a)), ex[1], ex[2])
            return mx
        end
    end
end

### Comparison

@sjdoc Comparison """
    Comparison(expr1,c1,expr2,c2,expr3,...)

performs or represents a chain of comparisons. `Comparison` expressions are usually input and
displayed using infix notation.
"""

@sjexamp(Comparison,
         ("Clear(a,b,c)",""),
         ("a == a","true"),
         ("a == b","false"),
         ("a < b <= c","a < b <= c"),
         ("(a=1,b=2,c=2)","2"),
         ("a < b <= c","true"))

@mkapprule Comparison  nodefault => true

## Following may be a stopgap in transition to binary comparisons
@doap function Comparison(x,op,y)
    mxpr(comparison_translation[op],x,y)
end

# Mma does this a == a != b  --->  a == a && a != b,  and  a == a  -->  True
# Note: Mma 10, at least does this: a == a != b  ---> a != b, in disagreement with the above

# We convert expressions that are not already numbers to floating point numbers, if possible.
# But, not for == or ===.
# We then use these approximations for comparison.
# This gives correct results for Pi>0, Sqrt(2) > 0, etc.
function maybe_N(x,cmp)
    (cmp != :(==)) && (cmp != :(===)) && (! isa(x,Number)) && is_Numeric(x) ? doeval(do_N(x)) : x
end

# FIXME: Don't convert all chained comparisons to conjunctions
# But, Mma  does this a < b < c,  ie. does not always return conjunctions.
# This always returns conjunctions if more than one comparison remains
# after removing true comparisons.
function do_Comparison(mx::Mxpr{:Comparison},args...)
    len = length(args)
    nargs = newargs()
    for i in 2:2:len
        a = args[i-1]
        cmp = args[i]
        b = args[i+1]
        an = maybe_N(a,cmp)
        bn = maybe_N(b,cmp)
        res = _do_Comparison(an,cmp,bn)
        if isa(res, Bool)
            res == false && return res
            push!(nargs,res)
        else
            push!(nargs,(a,cmp,b))
        end
    end
    nargs1 = newargs()
    for i in 1:length(nargs)
        a = nargs[i]
        if a != true
            push!(nargs1, mxpr(:Comparison, a...))
        end
    end
    length(nargs1) == 1 && return nargs1[1]
    mxpr(:And, nargs1)
end

# This does:  1 < 2 < b  -->  1 < 2 < b
function old_do_Comparison(mx::Mxpr{:Comparison},args...)
    len = length(args)
    for i in 2:2:len
        a = args[i-1]
        cmp = args[i]
        b = args[i+1]
        res = _do_Comparison(a,cmp,b)
        res == false && return res
        res != true && return mx
    end
    return true
end

function do_Comparison(mx::Mxpr{:Comparison},a::T,comp::SJSym,b::V) where {T<:Number,V<:Number}
    _do_Comparison(a,comp,b)
end

function _do_Comparison(a::T, comp::SJSym, b::V) where {T<:Number, V<:Number}
    if comp == :<    # Test For loop shows this is much faster than evaling Expr
        return a < b
    elseif comp == :>
        return a > b
    elseif comp == :(==)
        return a == b
    elseif comp == :(>=)
        return a >= b
    elseif comp == :(<=)
        return a <= b
    elseif comp == :(!=)
        return a != b
    elseif comp == :(===)
        return a === b
    end
    eval(Expr(:comparison,a,comp,b)) # This will be slow.
end

## FIXME Uh this is just copied from above. This is required to disambiguate
# from the catchall below
function _do_Comparison(a::T, comp::SJSym, b::T) where T<:Number
    if comp == :<    # Test For loop shows this is much faster than evaling Expr
        return a < b
    elseif comp == :>
        return a > b
    elseif comp == :(==)
        return a == b
    elseif comp == :(>=)
        return a >= b
    elseif comp == :(<=)
        return a <= b
    elseif comp == :(!=)
        return a != b
    elseif comp == :(===)
        return a === b
    end
    eval(Expr(:comparison,a,comp,b)) # This will be slow.
end

# This catches some cases
function _do_Comparison(mx::Mxpr{:DirectedInfinity}, comp::SJSym, n::T) where T<: Number
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return nothing
end

function _do_Comparison(n::T, comp::SJSym, mx::Mxpr{:DirectedInfinity}) where T<: Number
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return nothing
end

# FIXME. duplicated code. Maybe Symata needs its own Boolean type, one that is not <: Number
function _do_Comparison(mx::Mxpr{:DirectedInfinity}, comp::SJSym, n::Bool)
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return nothing
end

function _do_Comparison(n::Bool, comp::SJSym, mx::Mxpr{:DirectedInfinity})
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return nothing
end

# a == a  --> True, etc.  for unbound a
#function _do_Comparison{T<:Union{Mxpr,SJSym,AbstractString,DataType}}(a::T,comp::SJSym,b::T)
function _do_Comparison(a::T,comp::SJSym,b::V) where {T<:Union{Mxpr,SJSym,AbstractString,DataType}, V<:Union{Mxpr,SJSym,AbstractString,DataType}}
    if comp == :(==)
        res = a == b
        res && return res
    elseif comp == :(!=)
        res = a == b
        res && return false
    elseif comp == :(>=)  # Julia says  :a <= :b because symbols are ordred lexicographically
        res = a == b      # We don't want this behavior
        res && return res
    elseif comp == :(<=)
        res = a == b
        res && return res
    elseif comp == :(===)
        return a === b
    end
    return nothing
end

# TODO: Try to find why the Unions don't work and condense these methods
function _do_Comparison(a::Mxpr,comp::SJSym,b::Mxpr)
    if comp == :(==)
        res = a == b
        res && return res
    elseif comp == :(!=)
        res = a == b
        res && return false
    elseif comp == :(===)
        return a === b
    end
    return nothing
end

function _do_Comparison(a::Mxpr,comp::SJSym,b::SJSym)
    if comp == :(==)
        res = a == b
        res && return res
    elseif comp == :(!=)
        res = a == b
        res && return false
    elseif comp == :(===)
        return a === b
    end
    return nothing
end

_do_Comparison(a::SJSym, comp::SJSym, b::T) where {T<:SJReal} = nothing
_do_Comparison(a::T, comp::SJSym, b::SJSym) where {T<:Union{Mxpr,AbstractString,DataType}} = nothing

function _do_Comparison(a::T, comp::SJSym, b::Mxpr) where T<:SJReal
    nothing
end

function _do_Comparison(a::Mxpr, comp::SJSym, b::T) where T<:SJReal
    nothing
end

function _do_Comparison(a::T, comp::SJSym, b::Bool) where T<:Number
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return false
end

function _do_Comparison(a::Bool, comp::SJSym, b::Bool)
    comp == :(==) && return a == b
    comp == :(!=) && return a != b
    comp == :(===) && return a == b
    return false
end

function _do_Comparison(a, comp::SJSym, b::Bool)
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return false
end

_do_Comparison(a::Qsym, comp::Symbol, b::Bool) = nothing
_do_Comparison(a::Qsym, comp::SJSym, b::T) where {T<:Number} = nothing

function _do_Comparison(a, comp::SJSym, b::T) where T<:Number
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return false
end

# Note the asymmetry between this and previous method.
# This one, at least, is correct. and catches 2 < b
function _do_Comparison(a::T, comp::SJSym, b::SJSym) where T<:Number
    comp == :(==) && return false
    comp == :(!=) && return true
    comp == :(===) && return false
    return nothing
end

# used this to search for bug
# _do_Comparison{T<:Number, V<:Mxpr}(mx::V, comp::SJSym, n::T) = false

# Fix bug in (a == b) != False in mxp_test.sj, and similar expressions
function  _do_Comparison(mx::V, comp, n::T) where {T<:Bool, V<:Mxpr}
    comp == :(!=) && return true
    return false
end

function  _do_Comparison(mx::V, comp, n::T) where {T<:Number, V<:Mxpr}
    if typeof(comp) != SJSym
        symerror("_do_Comparison: Comparing with $comp, of type ", typeof(comp))
    else
        symerror("_do_Comparison: (assert error) Got symbol $comp, when expecting non-symbol. mx : $mx, n : $n")
    end
end

_do_Comparison(a::Mxpr, comp::SJSym, b::String) = false

# function _do_Comparison(a::Bool, comp::SJSym, b::Bool)
#     comp == :(==) && return a == b
#     comp == :(!=) && return a != b
#     comp == :(===) && return a == b
#     return false  # I guess this is good
# end

# This is meant to be a catchall for any object.
# But, we should use try catch because == may not be defined.
# Currently, this catches qsym
# function _do_Comparison{T}(a::T, comp::SJSym, b::T)
#     comp == :(==) && return a == b ? true : nothing
#     comp == :(!=) && return a != b ? nothing : true
#     comp == :(===) && return a == b
#     return nothing
# end

# function _do_Comparison{T}(a::T, comp::SJSym, b::T)
#     comp == :(==) && return a == b ? true : nothing
#     comp == :(!=) && return a != b ? nothing : true
#     comp == :(===) && return a == b
#     return nothing
# end

# FIXME. We need >=, <= like this in several places
# Break them out into a function
# NB. Mma leaves  a < a unevaluated.
# This is probably good because a may be of a type for which there is no order
function _do_Comparison(a::T, comp::SJSym, b::T) where T<:Qsym
    (comp == :(==) || comp == :(>=) || comp == :(<=))  && return a == b ? true : nothing
    comp == :(!=) && return a == b ?  false : nothing
    comp == :(===) && return a == b
#    (comp == :(<) || comp == :(>)) && return a == b ? false : nothing
    return nothing
end

#_do_Comparison(a::Qsym, comp::Symbol, b::Bool) = nothing

# function _do_Comparison{T}(a::Qsym, comp::SJSym, b::T)
#     return nothing
# end

function _do_Comparison(args...)
    symerror("No comparison for args ", args)
end

## These allow converting values returned by sympy, although we could do it differntly
apprules(mx::Mxpr{:<}) = mxpr(:Comparison,mx[1],:< ,mx[2])

ext: jl


text:
@enum(Vexity, VEX_CONST, VEX_AFFINE, VEX_CONVEX, VEX_CONCAVE, VEX_NONDCP, VEX_UNSET)

function -(v::Vexity)
    if v == VEX_CONVEX
        return VEX_CONCAVE
    elseif v == VEX_CONCAVE
        return VEX_CONCAVE
    end
    return v
end

function +(v::Vexity, w::Vexity)
    if (v == VEX_CONST) && (w == VEX_CONST)            # constant times other vexity cases
        return v
    elseif (v == VEX_CONST) && (w == VEX_NONDCP)
        return w
    elseif (v == VEX_NONDCP) && (w == VEX_CONST)
        return v
    elseif (v == VEX_CONST)
        return w
    elseif (w == VEX_CONST)
        return v
    elseif (v == VEX_AFFINE) && (w == VEX_AFFINE)      # affine times other vexity cases
        return v
    elseif (v == VEX_AFFINE) && (w == VEX_CONVEX)
        return w
    elseif (v == VEX_CONVEX) && (w == VEX_AFFINE)
        return v
    elseif (v == VEX_AFFINE) && (w == VEX_CONCAVE)
        return w
    elseif (v == VEX_CONCAVE) && (w == VEX_AFFINE)     # remaining vexity cases
        return v
    elseif (v == VEX_CONVEX) && (w == VEX_CONVEX)
        return v
    elseif (v == VEX_CONCAVE) && (w == VEX_CONCAVE)
        return v
    elseif (v == VEX_CONCAVE) && (w == VEX_CONVEX)
        return VEX_NONDCP
    end
    return VEX_NONDCP                                  # (v == VEX_CONVEX) && (w == VEX_CONCAVE)
end

curvature(f, x::Interval) = error("curvature(f,x) for $f is not implemented.")

function curvature(d2f::Interval{T}, x::Interval{T}) where T
    if zero(T) <= d2f
        return VEX_CONCAVE
    elseif zero(T) >= d2f
        return VEX_CONVEX
    elseif isatomic(x)
        return VEX_CONST
    end
    return VEX_NONDCP
end

for F in (positive, negative, lower_bnd, upper_bnd, bnd, rad2deg, deg2rad)
    @eval curvature(::typeof($F), x::Interval) = VEX_AFFINE
end

curvature(::typeof(+), x::Interval, y::Interval, args...) = VEX_CONST
curvature(::typeof(+), x::Interval) = VEX_CONST

curvature(::typeof(-), x::Interval, y::Interval) = VEX_AFFINE
curvature(::typeof(-), x::Interval) = VEX_AFFINE

for F in (abs, abs2, exp, exp2, exp10, expm1, leaky_relu, xlogx, cosh)
    @eval curvature(::typeof($F), x::Interval) = VEX_CONVEX
end

for F in (log, log2, log10, log1p, acosh, sqrt)
    @eval curvature(::typeof($F), x::Interval) = VEX_CONCAVE
end

curvature(::typeof(sin), x::Interval) = curvature(-sin(x), x)
curvature(::typeof(cos), x::Interval) = curvature(-cos(x), x)
curvature(::typeof(tan), x::Interval) = curvature(2*tan(x)*sec(x)^2, x)

for F in (asin, acot, sinh, csch, coth, atanh, acsch, xabsx, erfinv, erfc)
    @eval function curvature(::typeof($F), x::Interval)
        if x <= zero(T)
            return VEX_CONCAVE
        elseif x >= zero(T)
            return VEX_CONVEX
        end
        return VEX_NONDCP
    end
end

for F in (acos, atan, cbrt, tanh, erf, erfcinv)
    @eval function curvature(::typeof($F), x::Interval)
        if x <= zero(T)
            return VEX_CONVEX
        elseif x >= zero(T)
            return VEX_CONCAVE
        end
        return VEX_NONDCP
    end
end

function curvature(::typeof(xexpax), x::Interval{T}, a::Interval{T}) where T
    if isatomic(a)
        inflect_pnt = -2/a.lo
        if inflect_pnt < x
            return VEX_CONVEX
        elseif inflect_pnt > x
            return VEX_CONCAVE
        end
    end
    return VEX_NONDCP
end

# TODO: HANDLE ARGMAX/ARGMIN IS CONSTANT CASE -> VEXAFFINE CASE... (need z interval bounds?)
curvature(::typeof(max), x::Interval{T}, args...) = VEX_CONVEX
curvature(::typeof(min), x::Interval{T}, args...) = VEX_CONCAVE

function curvature(::typeof(relu), x::Interval{T}) where T
    if x >= zero(T)
        return VEX_CONVEX
    end
    return VEX_CONST
end

# TODO: sech
# TODO: maxtanh, pentanh, sigmoid, bisigmoid, softsign, softplus, elu

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
stress = zeros(n,9,3)
for i=1:length(elem[:,1])
    uv = [
    -1 -1
    1 -1
    1 1
    -1 1
    0 -1
    1 0
    0 1
    -1 0
    0 0]

    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]

    u = U[[i1 (i1 .+ Int64(0.5*dofs))]]'
    stress[i,1,:] = E*LN(uv[1,:],xs,ys)*u
    stress[i,2,:] = E*LN(uv[2,:],xs,ys)*u
    stress[i,3,:] = E*LN(uv[3,:],xs,ys)*u
    stress[i,4,:] = E*LN(uv[4,:],xs,ys)*u
    stress[i,5,:] = E*LN(uv[5,:],xs,ys)*u
    stress[i,6,:] = E*LN(uv[6,:],xs,ys)*u
    stress[i,7,:] = E*LN(uv[7,:],xs,ys)*u
    stress[i,8,:] = E*LN(uv[8,:],xs,ys)*u
    stress[i,9,:] = E*LN(uv[9,:],xs,ys)*u
 end
println("calculating nodal stress...")
fstress = zeros(3, size(nodes,1), 2)
fstress[:,:,2] .== 1.
for i = 1:n
    a = elem[i,:]' #global node number
    mapa = [findall(x -> x == a[k], elem[i, :])[1] for k = 1:length(a)] #lobal node number

    fstress[:, a, 1] = fstress[:, a, 1][:, 1, :] .+ stress[i, mapa, :]'
    fstress[:, a, 2] = fstress[:, a, 2][:, 1, :] .+ 1
end
println("averaging stress...")
num_indexed = fstress[:,:,2] #node sharing over elements
fstress = fstress[:,:,1] ./ fstress[:,:,2] #divide to obtain average
stress =1
elem = 1
edges = 1
bc = 1
ix = 1
F = 1
print("Averaging complete!")

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using DashBootstrapComponents

progress = html_div([
    dcc_interval(id = "progress-interval", n_intervals = 0, interval = 500),
    dbc_progress(id = "progress"),
]);

callback!(
    app,
    Output("progress", "value"),
    Output("progress", "label"),
    Input("progress-interval", "n_intervals"),
) do n
    # check progress of some background process, in this example we'll just
    # use n_intervals constrained to be in 0-100
    progress = min(n % 110, 100)
    # only add text after 5% progress to ensure text isn't squashed too much
    return progress >= 5 ? (progress, "$progress %") : (progress, "")
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


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

name = "Gingko"
version = v"1.4.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ginkgo-project/ginkgo/archive/refs/tags/v$(version).tar.gz", "6dcadbd3e93f6ec58ef6cda5b980fbf51ea3c7c13e27952ef38804058ac93f08")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ginkgo-*
mkdir build && cd build

cmake .. \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DGINKGO_BUILD_TESTS=OFF \
-DGINKGO_BUILD_EXAMPLES=OFF \
-DGINKGO_BUILD_BENCHMARKS=OFF \
-DGINKGO_BUILD_REFERENCE=OFF \
-DGINKGO_BUILD_DOC=OFF \
-DGINKGO_BUILD_HWLOC=OFF

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
 platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libginkgo_hip", :libgingko_hip),
    LibraryProduct("libginkgo_dpcpp", :libgingko_dpcpp),
    LibraryProduct("libginkgo_omp", :libgingko_omp),
    LibraryProduct("libginkgo_reference", :libgingko_reference),
    LibraryProduct("libginkgo_cuda", :libgingko_cuda),
    LibraryProduct("libginkgo_device", :libgingko_device),
    LibraryProduct("libginkgo", :libgingko)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
#needs at least gcc 5.5 according to README, BB only has 5.2.0
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"6.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("../examples/baker.jl")
using LinearAlgebra
using PyPlot
using PyCall
@pyimport matplotlib.collections as coll
@pyimport matplotlib.colors as cs
@pyimport matplotlib.cm as cm
function plot_y(m)
	s = [0.3, 0.3, 0.3, 0.0]
	n = 400
	x = LinRange(0.,2*pi,n)
	u = [[x[i] + 2*pi*rand(),x[j] + 2*pi*rand()]
		 for i=1:n,j = 1:n][:]
	q = [rand(2) for i=1:n, j=1:n][:]
	w = [zeros(2) for i=1:n, j=1:n][:]
	y = [zeros(2) for i=1:n, j=1:n][:]
	v = [zeros(2) for i=1:n, j=1:n][:]
	X = [zeros(2) for i=1:n, j=1:n][:]

	D2 = [zeros(2,2) for i=1:n, j=1:n][:]
	DXx = [zeros(2) for i=1:n, j=1:n][:]

	n = length(u)
	z = zeros(n)
	a = zeros(n)
	delta = zeros(n)
	ny = zeros(n)
	segments = zeros(n, 2, 2)
	eps = 1.e-1
	eps1 = 1.e-3
	filepath = string("/home/nishac/Research/PhDThesis/papers/",
					  "Decomposition-of-Linear-Response/figs/")
	for i = 1:20
		@show i
		u .= next.(u, Ref(s))
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
	end
	y_pos = 2.0
	x = Array{Float64}[]
	yx = Array{Array{Float64}}[]
	vx = Array{Array{Float64}}[]

	for i = 1:m
		@show i
		D2 .= pushforward_second_order.(u, q, Ref(s))
		y .= pushforward.(u, y, Ref(s))
		y .+= tensordot(D2, v)
		X .= pert.(u, 1) .+ pert.(u, 2) .+ pert.(u, 3)
		DXx .= dpert_x.(u, Ref(s))
		v .= pushforward.(u, v, Ref(s))
		v .+= X
		w .= pushforward.(u, w, Ref(s))
		w .+= tensordot(D2, q)

		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z

		a .= dot.(v, q)
		v .-= a.*q

		w ./= (z.*z)
		w .= w .- dot.(w, q).*q
		y ./= z
		y .+= -a.*w .+ DXx
		delta .= dot.(y, q) .+ dot.(v, w)
		y .-= delta.*q
		u .= next.(u, Ref(s))
		if i >=5
			inds = [abs(uk[2] - y_pos) < eps1 for uk in u]
			push!(x, [ui[1] for ui in u[inds]])
			push!(yx, y[inds])
			push!(vx, v[inds])
		end

	end
	qperp = [[-qi[2], qi[1]] for qi in q]
	ny = dot.(y, qperp)
	x = vcat(x...)
	yx = vcat(yx...)
	vx = vcat(vx...)
	vxy = [vxi[2] for vxi in vx]
	yxy = [yxi[2] for yxi in yx]
	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)
	plt.tight_layout()
	ax.plot(x, vxy,".", label=L"$v_{x_2}(x_2 = 6.0)$", ms=2.0, color="royalblue")
	#ax.plot([(x .- eps)'[1:10:end] .% (2*pi);
	#		(x .+ eps)'[1:10:end] .% (2*pi)],
	#		[(vxy .- eps*yxy)'[1:10:end];
	#		 (vxy .+ eps*yxy)'[1:10:end]],
#			color="blue")

	#ax.set_aspect(1.5)
	ax.grid(true)
	leg = ax.legend(loc="right", bbox_to_anchor=(0.8,0.2),
					fontsize=30, framealpha=0)
	leg.legendHandles[1]._legmarker.set_markersize(20)
	return x, vxy, yx
end

function plot_v(m)
	s = [0.1, 0.1, 0., 0.1]
	n = 400
	x = LinRange(0.,2*pi,n)
	u = [[x[i] + 2*pi*rand(),x[j] + 2*pi*rand()]
		 for i=1:n,j = 1:n][:]
	q = [rand(2) for i=1:n, j=1:n][:]
	v = [zeros(2) for i=1:n, j=1:n][:]
	X = [zeros(2) for i=1:n, j=1:n][:]
	n = length(u)
	z = zeros(n)
	a = zeros(n)
	segments = zeros(n, 2, 2)
	eps = 1.e-1
	eps1 = 1.e-3
	filepath = string("/home/nishac/Research/PhDThesis/papers/",
					  "Decomposition-of-Linear-Response/figs/")
	y_pos = 4.5
	x = Array{Float64}[]
	qx = Array{Array{Float64}}[]
	vx = Array{Array{Float64}}[]
	for i = 1:20
		@show i
		u .= next.(u, Ref(s))
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
	end
	for i = 1:m
		@show i
		X .= pert.(u, 1) .+ pert.(u, 2) .+ pert.(u, 4)
		v .= pushforward.(u, v, Ref(s))
		v .+= X
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
		a .= dot.(v, q)
		v .-= a.*q
		u .= next.(u, Ref(s))
		if i >=5
			inds = [abs(uk[2] - y_pos) < eps1 for uk in u]
			push!(x, [ui[1] for ui in u[inds]])
			push!(qx, q[inds])
			push!(vx, v[inds])
		end

	end
	x = vcat(x...)
	qx = vcat(qx...)
	vx = vcat(vx...)

	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)
	plt.tight_layout()
	ax.plot(x, vx,".",ms=0.5,label=(L"$v_{x_1}$", L"$v_{x_2}$"))
	#ax.plot(x, vx,".", color="darkcyan")
	fig.legend(fontsize=30)
	#ax.axis("scaled")

	return x, vx
end
function tensordot(A, b)
	return [Ai*b[i] for (i, Ai) in enumerate(A)]
end
function pushforward_second_order(u, v1, s)
	d2u = d2step(u, s)
	d2u_v1 = reshape([dot(d2u[:,i], v1) for i=1:4], 2, 2)
	return d2u_v1
end
function create_line_colls(u, q, eps)
	n = length(u)
	lc = zeros(n, 2, 2)
	for i = 1:n
		lc[i,1,:] = u[i] .- eps*q[i]
		lc[i,2,:] = u[i] .+ eps*q[i]
	end
	return lc
end
function pushforward(u, q, s)
	du = dstep(u, s)
	q = du*q
	return q
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ ffbfa7ea-ffa7-11eb-055a-417ea022d758
using DataFrames, CSV, VegaLite

# ╔═╡ 6acb2176-0f38-4529-ab54-4820d9467d90
md"""
# Stack Overflow Developer Survey Results

Each year, Stack Overflow conducts a [survey](https://insights.stackoverflow.com/survey/) of its users on a variety of aspects.  Two of these aspects are gender identity and the programming languages used.  Given this data, we can look at how the gender diversity of various languages compare to each other.
"""

# ╔═╡ 92ca8227-beff-44e6-ad57-6204b88abe45
# 2020 Data from https://insights.stackoverflow.com/survey/
DATA_PATH="survey_results_public.csv"

# ╔═╡ 0f325b73-65ba-49a8-aa6f-b56c8e4ea548
md"""
R had the highest share of respondents identifying as a gender minority with 12.68%, and Julia has the **third most with 9.25%**.

9.25% of Julia users identifying as a gender minority is notably higher than the Julia Developer Survey found (5%), though lower than the Google Analytics estimates of traffic to julialang.org (~20%).
"""

# ╔═╡ 4789d7b1-7dd4-4b6b-8111-c96f010cf213
md"""
## Correlation of user gender by language pairs

To go a layer deeper, we can look at the percentage of users who identify as gender minorities by language pairs to see the correlation between gender minorities within different communities. To do this, we find the users who use both languages in the pair and collect the count and % gender minority.

The plot is a diagonal as the language pairs values are symmetrical. The diagonal of the plot represents the overall statstics for that language (total count, % gender minority).

We can see that Julia users overlap most with Python, R, and SQL. The Julia/R users have a higher share of gender minorities than Julia alone, but is similar to the overall share for R. The Julia/SQL pair has a slightly higher gender minority share than just Julia.  Then the Julia/Python users have a similar proportion gender minorities to Julia as a whole.
"""

# ╔═╡ 8983a52e-7c3a-44ab-9edd-eaa2bcb4da06
function is_gender_minority(response)
	gender_minorities = Set(["Woman", "Non-binary, genderqueer, or gender non-conforming"])
	return !isempty(intersect(gender_minorities, split(response, ";")))
end

# ╔═╡ bd01939a-4e41-4efb-9ce3-5c1bd728ee0b
begin
	results = CSV.File(DATA_PATH; select=["Respondent", "Gender", "LanguageWorkedWith"], threaded=false) |> DataFrame
	results.languages = split.(results.LanguageWorkedWith, ";")
	results.GenderMinority = is_gender_minority.(results.Gender)
	results
end

# ╔═╡ 223c3d89-e310-41e2-a6c3-bb8103067011
overall_pct_gm = sum(results.GenderMinority) / nrow(results) * 100

# ╔═╡ cc82fc2a-ae81-4f5e-8a54-e091f855b837
unique_languages = filter(x -> x != "NA", unique(reduce(vcat, results.languages)))

# ╔═╡ 5d0bd506-bed4-44e7-8872-85a86922ba3e
begin
	summary = []

	for a in unique_languages
		for b in unique_languages
			language_a = a < b ? a : b
			language_b = a < b ? b : a
			filtered = filter(row -> language_a ∈ row.languages && language_b ∈ row.languages, results)
			count = nrow(filtered)
			push!(summary, (;language_a, language_b, count, gm_pct = sum(filtered.GenderMinority) / count))
		end
	end

	summary = unique(filter(row -> row.count > 0, DataFrame(summary)))
end

# ╔═╡ 091f24ea-b915-4053-b1e1-67bbdf1bdcc8
summary |> @vlplot(
	:circle,
	x={"language_a:n", axis={ticks=false, domain=false, labelColor="grey"}},
	y={"language_b:n", axis={ticks=false, domain=false, labelColor="grey"}},
	size={"count:q", scale={type="log"}},
	color={"gm_pct:q",
		scale={type="threshold", domain=[0.05, 0.075, 0.1, 0.125], scheme="turbo", reverse=true}
	},
	config={view={stroke="transparent"}}
)

# ╔═╡ 57d7edf7-7500-48ac-a598-c23234b72c4d
function gm_pct(language, df)
	filtered = filter(row -> language ∈ row.languages, df)
	return sum(filtered.GenderMinority) / nrow(filtered) * 100
end

# ╔═╡ f8daa230-f572-4382-a579-99e510615802
begin
	df = DataFrame()
	df.language = unique_languages
	df.pct_gm =  map(language -> gm_pct(language, results), unique_languages)
	sort!(df, :pct_gm; rev=true)
	df
end

# ╔═╡ 9b27ba42-0440-4c73-9218-067f4c8c35a3
df |> @vlplot(:bar, x={"language:n", sort="-pct_gm", title="Language"}, y={:pct_gm, title="% Gender Minority"})

# ╔═╡ 8561a887-46ec-4500-942c-04a300472534
round2(n) = round(n, digits=2)

# ╔═╡ 78c81d6c-dcc8-4f36-98b4-267ea8be1c10
md"""
The `Gender` and `LanguageWorkedWith` columns both contain `;` delimited lists of gender identities and languages.  To determine the percent for each language identifying as a gender minority, we summarize the gender identities to a `GenderMinority` value, and for each language in the survey, find the sum of gender minorities using that language vs users using that language of any gender.

The overall percent of respondents identifying as a gender minority is $(round2(overall_pct_gm))%
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
VegaLite = "~2.6.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "344f143fa0ec67e47917848795ab19c6a455f32c"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.32.0"

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[FilePaths]]
deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
version = "0.8.3"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "0f5e8d0cb91a6386ba47bd1527b240bd5725fbae"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.10"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.8.19"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSONSchema]]
deps = ["HTTP", "JSON", "ZipFile"]
git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "0.4.5"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[NodeJS]]
deps = ["Pkg"]
git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a"
uuid = "2bd173c7-0d6d-553b-b6af-13a54713934c"
version = "1.3.0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "a3a337914a035b2d59c9cbe7f1a38aaba1265b02"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.6"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[TableTraitsUtils]]
deps = ["DataValues", "IteratorInterfaceExtensions", "Missings", "TableTraits"]
git-tree-sha1 = "8fc12ae66deac83e44454e61b02c37b326493233"
uuid = "382cd787-c1b6-5bf2-a167-d5b971a19bda"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Vega]]
deps = ["DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "JSONSchema", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "Setfield", "TableTraits", "TableTraitsUtils", "URIParser"]
git-tree-sha1 = "43f83d3119a868874d18da6bca0f4b5b6aae53f7"
uuid = "239c3e63-733f-47ad-beb7-a12fde22c578"
version = "2.3.0"

[[VegaLite]]
deps = ["Base64", "DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "TableTraits", "TableTraitsUtils", "URIParser", "Vega"]
git-tree-sha1 = "3e23f28af36da21bfb4acef08b144f92ad205660"
uuid = "112f6efa-9a02-5b7d-90c0-432ed331239a"
version = "2.6.0"

[[ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─6acb2176-0f38-4529-ab54-4820d9467d90
# ╠═ffbfa7ea-ffa7-11eb-055a-417ea022d758
# ╟─92ca8227-beff-44e6-ad57-6204b88abe45
# ╟─78c81d6c-dcc8-4f36-98b4-267ea8be1c10
# ╟─bd01939a-4e41-4efb-9ce3-5c1bd728ee0b
# ╟─223c3d89-e310-41e2-a6c3-bb8103067011
# ╟─0f325b73-65ba-49a8-aa6f-b56c8e4ea548
# ╟─f8daa230-f572-4382-a579-99e510615802
# ╟─9b27ba42-0440-4c73-9218-067f4c8c35a3
# ╟─4789d7b1-7dd4-4b6b-8111-c96f010cf213
# ╟─091f24ea-b915-4053-b1e1-67bbdf1bdcc8
# ╟─5d0bd506-bed4-44e7-8872-85a86922ba3e
# ╟─cc82fc2a-ae81-4f5e-8a54-e091f855b837
# ╟─8983a52e-7c3a-44ab-9edd-eaa2bcb4da06
# ╟─57d7edf7-7500-48ac-a598-c23234b72c4d
# ╟─8561a887-46ec-4500-942c-04a300472534
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
import Observables
import Makie

export Whirl
export update

"""
    Represents a whirl.

fields: points, gauge1, gauge2, configuration, segments, color and observable.
"""
mutable struct Whirl <: Sprite
    points::Array{ComplexPlane,1}
    gauge1::Array{U1,1}
    gauge2::Array{U1,1}
    configuration::Biquaternion
    segments::Int
    color::Observables.Observable{Array{Makie.ColorTypes.RGBA{Float32},2}}
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Whirl(scene,
          points,
          [gauge1, [gauge2, [configuration, [segments, [color, [transparency]]]]]])

Construct a whirl with the given `scene`, `points`, `section`, `gauge1`, `gauge2`,
`configuration`, `segments`, `color` and `transparency`.
"""
function Whirl(scene::Makie.Scene,
               points::Array{ComplexPlane,1},
               gauge1::Array{U1,1},
               gauge2::Array{U1,1};
               configuration::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
               segments::Int = 36,
               color::Makie.RGBAf = Makie.RGBAf(1.0, 0.2705, 0.0, 0.5),
               transparency::Bool = false)
    whirl = constructwhirl(points, gauge1, gauge2, configuration, segments)
    colorarray = Observables.Observable(fill(color, size(whirl)...))
    observable = buildsurface(scene, whirl, colorarray, transparency = transparency)
    Whirl(points, gauge1, gauge2, configuration, segments, colorarray, observable)
end

"""
    update(whirl, points)

Update a Whirl by changing its observable with the given `whirl` and `points`.
"""
function update(whirl::Whirl, points::Array{ComplexPlane,1})
    whirl.points = points
    value = constructwhirl(whirl.points,
                           whirl.gauge1,
                           whirl.gauge2,
                           whirl.configuration,
                           whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, gauge1, gauge2)

Update a Whirl by changing its observable with the given `whirl`, `gauge1` and `gauge2`.
"""
function update(whirl::Whirl, gauge1::Array{U1,1}, gauge2::Array{U1,1})
    whirl.gauge1 = gauge1
    whirl.gauge2 = gauge2
    value = constructwhirl(whirl.points,
                           whirl.gauge1,
                           whirl.gauge2,
                           whirl.configuration,
                           whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, points, gauge1, gauge2)

Update a Whirl by changing its observable with the given `whirl`, `points`, `gauge1` and `gauge2`.
"""
function update(whirl::Whirl, points::Array{ComplexPlane,1}, gauge1::Array{U1,1}, gauge2::Array{U1,1})
    whirl.points = points
    whirl.gauge1 = gauge1
    whirl.gauge2 = gauge2
    value = constructwhirl(whirl.points,
                           whirl.gauge1,
                           whirl.gauge2,
                           whirl.configuration,
                           whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, configuration)

Update a Whirl by changing its observable with the given `whirl` and `configuration`.
"""
function update(whirl::Whirl, configuration::Biquaternion)
    whirl.configuration = configuration
    value = constructwhirl(whirl.points,
                           whirl.gauge1,
                           whirl.gauge2,
                           whirl.configuration,
                           whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, color)

Update a Whirl by changing its observable with the given `whirl` and `color`.
"""
function update(whirl::Whirl, color::Makie.RGBAf)
    whirl.color[] = fill(color, whirl.segments, length(whirl.points))
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module SaferIntegers

export SafeInteger, SafeSigned, SafeUnsigned,
       SafeInt, SafeInt8, SafeInt16, SafeInt32, SafeInt64, SafeInt128,
       SafeUInt, SafeUInt8, SafeUInt16, SafeUInt32, SafeUInt64, SafeUInt128,
       @saferintegers,
       SafeRational

import Base.Checked: checked_neg, checked_abs, checked_add, checked_sub, checked_mul,
                     checked_div, checked_rem, checked_fld, checked_mod, checked_cld,
                     add_with_overflow, sub_with_overflow, mul_with_overflow

import Base: promote_rule,
    string, bitstring, show, repr, parse, tryparse,
    (~), (&), (|), (⊻), (>>>), (>>), (<<),
    (<), (<=), (==), (!=), (>=), (>),
    isequal, isless,
    abs, abs2,
    (+), (-), (*), (/), (\), (^),
    div, fld, fld1, cld, rem, mod, mod1, divrem, fldmod, fldmod1,
    zero, one, sizeof, typemax, typemin, widen,
    signbit, sign, count_ones, count_zeros, ndigits0z,
    leading_zeros, trailing_zeros, leading_ones, trailing_ones,
    copysign, flipsign,
    gcd, lcm, divgcd

using Random

include("type.jl")
include("construct.jl")
include("promote.jl")
include("int_ops.jl")
include("binary_ops.jl")
include("arith_ops.jl")
include("pow.jl")
include("string_io.jl")
include("parse.jl")
include("rand.jl")
include("cover.jl")

include("SaferIntTypes.jl")
using .SaferIntTypes

include("SaferRationalTypes.jl")

end # module SaferIntegers


text:
using Graphs, SimpleWeightedGraphs, StatsBase

"""
    Convert a dictionary keyed by edge tuples to
    integer weights into 3 vectors for input
    to SimpleWeightedDigraph.

"""
function vectors_from_edgedict(edgedict::Dict{Tuple, Integer})
    src = zeros(Int64, length(edgedict))
    dest = zeros(Int64, length(edgedict))
    weights = zeros(Int64, length(edgedict))
    idx = 1
    for tup in edgedict
        src[idx] = tup[1][1]
        dest[idx] = tup[1][2]
        weights[idx] = tup[2][1]
        idx += 1
    end
    return src, dest, weights
end

"""
    Add a node to a mapping dictionary
    and return the incremented index
    for convenience.

"""
function add_to_mapping!(mapping, node, idx)
    if !haskey(mapping, node)
        mapping[node] = idx
        idx = idx + 1
    end
    return idx
end

"""
"""
function remove_selfloops!(w)
    i = 2
    e = length(w)+1
    while i < e
        if w[i-1] == w[i]
            deleteat!(w, i-1)
            e-=1
        else
            i+=1
        end
    end
end

"""
    Compute a DeBruijn graph of order k and it's
    associated first order graph from an ngram file
    with or without frequency values.

"""
function from_ngram(file::String, frequency::Bool, k::Integer)
    # Initialize empty mappings and edgelists
    # fo = first-order
    # ko = kth-order
    fo_map = Dict{String, Integer}()
    ko_map = Dict{Tuple, Integer}()
    fo_edgelist = Dict{Tuple, Integer}()
    ko_edgelist = Dict{Tuple, Integer}()
    fin = open(file)
    fo_idx = 1
    ko_idx = 1
    for line in eachline(fin)
        s = split(line, ',')
        if frequency
            walk = s[1:end-1]
            freq = parse(Float64, s[end])
        else
            walk = s
            freq = 1
        end
        remove_selfloops!(walk)
        for i in range(1, length(walk))
            if i > 1
                u = walk[i-1]
                v = walk[i]
		if u == v
		    println("Self-loop: $u, $v")
		end
                fo_idx = add_to_mapping!(fo_map, u, fo_idx)
                fo_idx = add_to_mapping!(fo_map, v, fo_idx)
                if haskey(fo_edgelist, (fo_map[u],fo_map[v]))
                    fo_edgelist[(fo_map[u],fo_map[v])] += freq
                else
                    fo_edgelist[(fo_map[u],fo_map[v])] = freq
                end
            end
            if i < length(walk)-k+1
                u = Tuple(walk[i:i+k-1])
                v = Tuple(walk[i+1:i+k])
		if u == v
		    println("Self-loop: $u, $v")
		end
                ko_idx = add_to_mapping!(ko_map, u, ko_idx)
                ko_idx = add_to_mapping!(ko_map, v, ko_idx)

                if haskey(ko_edgelist, (ko_map[u],ko_map[v]))
                    ko_edgelist[(ko_map[u],ko_map[v])] += freq
                else
                    ko_edgelist[(ko_map[u],ko_map[v])] = freq
                end
            end
        end
    end
    close(fin)

    # Convert edge dictionary representations into SimpleWeightedDigraphs
    fo_src, fo_dest, fo_weights = vectors_from_edgedict(fo_edgelist)
    forder = SimpleWeightedDiGraph(fo_src, fo_dest, fo_weights)
    ko_src, ko_dest, ko_weights = vectors_from_edgedict(ko_edgelist)
    korder = SimpleWeightedDiGraph(ko_src, ko_dest, ko_weights)
    return forder, fo_map, korder, ko_map
end

"""
Convert edges of ko to walks of length k. Use edge
weights in ko as weights for walks.
"""
function walks_from_edges(ko, ko_map, fo_map)
    A = adjacency_matrix(ko)
    cids = findall(>(0), A)
    edges_as_walks = Vector{Tuple}(undef, length(cids))
    ko_rev_map = Dict(value => key for (key, value) in ko_map)
    weights = Vector{Int64}(undef, ne(ko))
    Threads.@threads for aidx in range(1, length(cids))
        cidx = cids[aidx]
        u, v = Tuple(cidx)
        ko_u = Tuple([fo_map[j] for j in ko_rev_map[u]])
        ko_v = Tuple([fo_map[j] for j in ko_rev_map[v]])
        walk = Tuple{Vararg{Int64}}([ko_u...,ko_v[end]])
        edges_as_walks[aidx] = walk
        weights[aidx] = A[u, v]
    end
    return edges_as_walks, weights
end

"""
Construct a graph from input walks.

NOTE: Expects nodes to be integers.
"""
function from_walks(walks::Vector{Tuple}, k::Int64, rev_fo_map::Dict{Integer, String})
    ko_map = Dict{Tuple, Integer}()
    fo_edgelist = Dict{Tuple, Integer}()
    ko_edgelist = Dict{Tuple, Integer}()
    new_konode_idx = 1
    for w in walks
        for i in range(2, length(w))
            edge = w[i-1]...,w[i]
            if !haskey(fo_edgelist, edge)
                fo_edgelist[edge] = 0
            fo_edgelist[edge] += 1
            end
        end

        for i in range(1, length(w)-k)
            mapped = Tuple{Vararg{String}}([String(rev_fo_map[u]) for u in w[i:i+k]])
            u = mapped[i:i+k-1]
            v = mapped[i+1:i+k]
            if !haskey(ko_map, u)
                ko_map[u] = new_konode_idx
                new_konode_idx += 1
            end
            if !haskey(ko_map, v)
                ko_map[v] = new_konode_idx
                new_konode_idx += 1
            end
            edge = Tuple{Integer, Integer}([ko_map[u],ko_map[v]])
            if !haskey(ko_edgelist, edge)
               ko_edgelist[edge] = 0
            end
            ko_edgelist[edge] += 1
        end
    end
    fo_src, fo_dest, fo_weights = vectors_from_edgedict(fo_edgelist)
    forder = SimpleWeightedDiGraph(fo_src, fo_dest, fo_weights)
    ko_src, ko_dest, ko_weights = vectors_from_edgedict(ko_edgelist)
    korder = SimpleWeightedDiGraph(ko_src, ko_dest, ko_weights)
    return forder, korder, ko_map
end

"""
Return all k-edge walks based on first-order
graph G.
"""
function get_all_walks(G, k)
    adjacency_dict = get_adj_dict(G, adjacency_matrix(G))
    walklist = Vector{Tuple}()
    for node in range(1, nv(G))
        nodewalks = all_walks(adjacency_dict, node, k)
        append!(walklist, nodewalks)
    end
    return walklist
end

"""
Return all k-edge walks based on first-order
graph G and also return all k-order nodes.
"""
function get_all_walks_with_nodes(G, k)
    adjacency_dict = get_adj_dict(G, adjacency_matrix(G))
    walklist = Vector{Tuple}()
    nodes = Set{Tuple}()
    for node in range(1, nv(G))
        nodewalks = all_walks(adjacency_dict, node, k)
        append!(walklist, nodewalks)
        for walk in nodewalks
            u = walk[1:end-1]
            v = walk[2:end]
            push!(nodes, u)
            push!(nodes, v)
        end
    end
    return walklist, nodes
end


text:
function move(::InPlace, z, x, s, d, α)
    @. s = α * d
    @. z = x + s
    z, s
end
function move(::OutOfPlace, z, x, s, d, α)
    s = @. α * d
    z = @. x + s
    z, s
end
include("static.jl")
include("backtracking.jl")
include("hagerzhangline.jl")

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using AbstractAlgebra

@testset "ff isgroebner" begin
    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:degrevlex)

    @test Groebner.isgroebner([R(1)])
    @test Groebner.isgroebner([x])
    @test Groebner.isgroebner([x, y])
    @test Groebner.isgroebner([x, x])
    @test Groebner.isgroebner([x, y, z])

    @test !Groebner.isgroebner([x + y, x])
    @test Groebner.isgroebner([z*x, z*x, R(1)])

    fs = Groebner.change_ordering(Groebner.rootn(3, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs)
    @test Groebner.isgroebner(Groebner.groebner(fs))

    fs = Groebner.change_ordering(Groebner.rootn(5, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs)
    @test Groebner.isgroebner(Groebner.groebner(fs))
    @test Groebner.isgroebner(Groebner.groebner(fs, reduced=false))

    fs = Groebner.change_ordering(Groebner.noonn(2, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs)
    @test Groebner.isgroebner(Groebner.groebner(fs))

    fs = Groebner.change_ordering(Groebner.noonn(3, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs)
    @test Groebner.isgroebner(Groebner.groebner(fs))
    @test Groebner.isgroebner(Groebner.groebner(fs, reduced=false))

    fs = Groebner.change_ordering(Groebner.noonn(6, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs)
    @test Groebner.isgroebner(Groebner.groebner(fs))

end


text:
module test_julia_replace

using Test

@test replace([1, 2, 3], 2=>0, 3=>0, 0=>5) == [1, 0, 0]
@test replace("abble", 'b'=>'p') == "apple"
@test replace("abble", "bb"=>"pp") == "apple"

end # module test_julia_replace


text:
import DataFrames

"""
    report(df::DataFrames.DataFrame; kwargs...)

Describe the variables in a DataFrame.

# Arguments
- `df`: DataFrame.
- `missing_percentage::Bool`: Show missings by percentage (default) or number.
- `levels_percentage::Bool`: Show factor levels by percentage (default) or number.
- `median::Bool`: Show `mean` and `sd` (default) or `median` and `mad`.
- `dispersion::Bool`: Show dispersion (`sd` or `mad`).
- `range::Bool`: Show range.
- `n_strings::Int`: Number of different string elements to show.

!!! note

    **Ideas / help required:**
    - Add more indices (See [#19](https://github.com/neuropsychology/Psycho.jl/issues/19))
    - Deal with SubDataFrames (See [#20](https://github.com/neuropsychology/Psycho.jl/issues/20))

# Examples
```julia
report(simulate_data_correlation([[0.1], [0.2]]))
```
"""
function report(df::DataFrames.DataFrame; kwargs...)
    text = "The data contains $(DataFrames.nrow(df)) observations of the following variables:"

    for (colname, var) in DataFrames.eachcol(df)
        description = report(var; kwargs...)
        text = text * "\n  - " * String(colname) * " (" * description.text * ")"
    end

    return Report(text=text)
end


text:
using LinearAlgebra
include("../ellpla.jl")
include("../expmap.jl")

x = [0.; 0.; 1.];
a = 1.;
b = 1.;
c = 1.05;
R = Matrix{Float64}(I, 3, 3)
point = [0.; 0.; 0.];
normal = [0.; 0.; 1.];

(conpnt, depth, A, B) = ellipsoid_plane_contact(x, a, b, c, R, point, normal)

println("-------------------------------------------------------------")
println("for x = ", x, "; a, b, c = ", a, ", ", b, ", ", c, "; R = ", R)
println("and point = ", point, ", normal = ", normal)
println("conpnt = ", conpnt, ", depth = ", depth, ", A, B = ", A, ", ", B)

a = 0.1
x[3] = 1.049

(conpnt, depth, A, B) = ellipsoid_plane_contact(x, a, b, c, R, point, normal)

println("-------------------------------------------------------------")
println("for x = ", x, "; a, b, c = ", a, ", ", b, ", ", c, "; R = ", R)
println("and point = ", point, ", normal = ", normal)
println("conpnt = ", conpnt, ", depth = ", depth, ", A, B = ", A, ", ", B)

b = 0.1
x[3] = 1.0499

(conpnt, depth, A, B) = ellipsoid_plane_contact(x, a, b, c, R, point, normal)

println("-------------------------------------------------------------")
println("for x = ", x, "; a, b, c = ", a, ", ", b, ", ", c, "; R = ", R)
println("and point = ", point, ", normal = ", normal)
println("conpnt = ", conpnt, ", depth = ", depth, ", A, B = ", A, ", ", B)

a = 0.5
b = 0.7
x[3] = 0.7
R = expmap([1.0 1.0 1.0])

(conpnt, depth, A, B) = ellipsoid_plane_contact(x, a, b, c, R, point, normal)

println("-------------------------------------------------------------")
println("for x = ", x, "; a, b, c = ", a, ", ", b, ", ", c, "; R = ", R)
println("and point = ", point, ", normal = ", normal)
println("conpnt = ", conpnt, ", depth = ", depth, ", A, B = ", A, ", ", B)

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
"""
    project_reach(Rsets, vars, n, options)

Projection of a reachability analysis result in 2D.

### Input

- `Rsets`   -- reachable states representation
- `vars`    -- variables to plot; two-dimensional index vector
- `n`       -- system dimension
- `options` -- options

### Notes

The `vars` argument is required even if the optional argument
`projection_matrix` is passed, because we also determine whether time is used as
a dimension from this variable.
"""
function project_reach(
        Rsets::Vector{<:ReachSet{<:LazySets.LazySet{numeric_type}}},
        vars::Vector{Int64},
        n::Int64,
        options::AbstractOptions)::Vector{<:ReachSet} where {numeric_type<:Real}

    # parse input
    @assert(length(vars) == 2)
    if n == 2
        return project_2d_reach(Rsets, vars, n, options)
    end

    # first projection dimension
    xaxis = vars[1]
    if xaxis == 0
        got_time = true
        xaxis = n+1 # we add a new dimension for time
    else
        got_time = false
        if (xaxis <= 0 || xaxis > n)
            throw(DomainError("value $xaxis for X variable not allowed"))
        end
    end

    # build projection matrix
    projection_matrix = options[:projection_matrix]
    output_function = !options[:project_reachset]
    m = got_time ? n+1 : n
    if projection_matrix == nothing
        # projection to a state variable
        yaxis = vars[2]
        if (yaxis <= 0 || yaxis > n)
            throw(DomainError("value $yaxis for Y variable not allowed"))
        end
        projection_matrix = sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, m)
    else
        @assert(size(projection_matrix) == (1,n))
        # make vector a 2-row matrix
        projection_matrix = sparse(fill(2, n), 1:n, projection_matrix[1, :], 2, m)
        projection_matrix[1, xaxis] = 1.0
    end

    # apply optional transformation to projection matrix
    if options[:transformation_matrix] != nothing
        transformation_matrix = options[:transformation_matrix]
        if got_time
            # add another dimension for time: block matrix [S 0; 0 1]
            transformation_matrix = sparse(cat([1, 2], transformation_matrix, [1]))
        end
        projection_matrix = projection_matrix * transformation_matrix
    end

    N = length(Rsets)

    # allocate output and define overapproximation function
    ε = options[:ε_proj]
    if ε < Inf
        oa = x -> overapproximate(x, HPolygon, ε)
        RsetsProj = Vector{ReachSet{HPolygon{numeric_type}}}(undef, N)
    else
        set_type = options[:set_type_proj]
        oa = x -> overapproximate(x, set_type)
        RsetsProj = Vector{ReachSet{set_type{numeric_type}}}(undef, N)
    end

    if got_time
        @inbounds for i in 1:N
            t0 = Rsets[i].t_start
            t1 = Rsets[i].t_end
            radius = (t1 - t0)/2.0
            RsetsProj[i] = ReachSet(
                oa(projection_matrix *
                    CartesianProduct(Rsets[i].X, BallInf([t0 + radius], radius))),
                Rsets[i].t_start, Rsets[i].t_end)
        end
    else
        @inbounds for i in 1:N
            RsetsProj[i] = ReachSet(
                oa(projection_matrix * Rsets[i].X),
                Rsets[i].t_start, Rsets[i].t_end)
        end
    end

    return RsetsProj
end

"""
    project_reach(Rsets, vars, n, options)

This function projects a sequence of sets into the time variable, or can be
used to take a linear combination of the given variables.

### Input

- `Rsets`   -- reachable states representation
- `vars`    -- variables to plot; two-dimensional index vector
- `n`       -- system dimension
- `options` -- options

### Notes

The input `Rsets` is an array of sets (instead of a `CartesianProductArray`).
This array contains the collection of reach sets in 2D.

It is assumed that the variable given in vars belongs to the block computed
in the sequence of 2D sets `Rsets`.
"""
function project_2d_reach(
        Rsets::Vector{<:ReachSet{<:LazySets.LazySet{numeric_type}}},
        vars::Vector{Int64},
        n::Int64,
        options::AbstractOptions)::Vector{<:ReachSet} where {numeric_type<:Real}

    # first projection dimension
    xaxis = vars[1]
    if xaxis == 0
        got_time = true
        xaxis = 3  # time is associated to dimension 3
    else
        got_time = false
        if (xaxis <= 0 || xaxis > n)
            throw(DomainError())
        end
        # map back to 2d block
        xaxis = iseven(xaxis) ? 2 : 1
    end

    # build projection matrix
    projection_matrix = options[:projection_matrix]
    output_function = !options[:project_reachset]
    if projection_matrix == nothing
        # projection to a state variable
        yaxis = vars[2]
        if (yaxis <= 0 || yaxis > n)
            throw(DomainError())
        end
        # map back to 2d block
        yaxis = iseven(yaxis) ? 2 : 1
        m = got_time ? 3 : 2
        projection_matrix = sparse([1, 2], [xaxis, yaxis], [1.0, 1.0], 2, m)
    elseif !output_function
        error("projection matrix not allowed for this algorithm")
    end

    N = length(Rsets)

    # allocate output and define overapproximation function
    ε = options[:ε_proj]
    if ε < Inf
        oa = x -> overapproximate(x, HPolygon, ε)
        RsetsProj = Vector{ReachSet{HPolygon{numeric_type}}}(undef, N)
    else
        set_type = options[:set_type_proj]
        oa = x -> overapproximate(x, set_type)
        RsetsProj = Vector{ReachSet{set_type{numeric_type}}}(undef, N)
    end

    if output_function
        @inbounds for i in 1:N
            t0 = Rsets[i].t_start
            t1 = Rsets[i].t_end
            radius = (t1 - t0)/2.0
            RsetsProj[i] = ReachSet(
                oa(CartesianProduct(BallInf([t0 + radius], radius), Rsets[i].X)),
                Rsets[i].t_start, Rsets[i].t_end)
        end
    elseif got_time # x variable is 'time'
        @inbounds for i in 1:N
            t0 = Rsets[i].t_start
            t1 = Rsets[i].t_end
            radius = (t1 - t0)/2.0
            RsetsProj[i] = ReachSet(
                oa(projection_matrix *
                    CartesianProduct(Rsets[i].X, BallInf([t0 + radius], radius))),
                Rsets[i].t_start, Rsets[i].t_end)
        end
    else
        @inbounds for i in 1:N
            RsetsProj[i] = ReachSet(oa(projection_matrix * Rsets[i].X),
                Rsets[i].t_start, Rsets[i].t_end)
        end
    end

    return RsetsProj
end

"""
    project(Rsets, options)

Projects a sequence of sets according to the settings defined in the options.

### Input

- `Rsets`   -- solution of a reachability problem
- `options` -- options structure

### Notes

A projection matrix can be given in the options structure, or passed as a
dictionary entry.
"""
function project(Rsets::Vector{<:ReachSet}, options::AbstractOptions)
    plot_vars = copy(options[:plot_vars])
    for i in 1:length(plot_vars)
        if plot_vars[i] != 0
            plot_vars[i] = options[:inout_map][plot_vars[i]]
        end
    end
    reduced_n = sum(x -> x != 0, options[:inout_map])
    output_function = !options[:project_reachset]
    RsetsProj = project_reach(Rsets, plot_vars, reduced_n, options)
end

project(reach_sol::ReachSolution) = project(reach_sol.Xk, reach_sol.options)

project(Rsets::Vector{<:ReachSet}, options::Pair{Symbol,<:Any}...) =
    project(Rsets, Options(Dict{Symbol,Any}(options)))


text:
"""
    struct RepresentationSpace{G<:Sector} <: AbstractRepresentationSpace{G}

Generic implementation of a representation space, i.e. a complex Euclidean space
with a direct sum structure corresponding to different superselection sections of
type `G<:Sector`, e.g. the irreps of a compact or finite group, or the labels of
a unitary fusion category.
"""
struct RepresentationSpace{G<:Sector} <: AbstractRepresentationSpace{G}
    dims::Dict{G,Int}
    dual::Bool
end
RepresentationSpace(dims::Vararg{Pair{G,Int}}; dual::Bool = false) where {G<:Sector} = RepresentationSpace{G}(Dict{G,Int}(dims...), dual)
RepresentationSpace{G}(dims::Vararg{Pair{<:Any,Int}}; dual::Bool = false) where {G<:Sector} = RepresentationSpace{G}(Dict{G,Int}(convert(G,s)=>d for (s,d) in dims), dual)

# Corresponding methods:
function sectors(V::RepresentationSpace)
    s = collect(keys(V.dims))
    V.dual ? s : map!(conj,s,s)
end

dim(V::RepresentationSpace) = sum(values(V.dims))
dim(V::RepresentationSpace{G}, c::G) where {G<:Sector} = get(V.dims, c, 0)

Base.conj(V::RepresentationSpace) = RepresentationSpace(copy(V.dims), !V.dual)

Base.getindex(::ComplexNumbers, G::Type{<:Sector}) = RepresentationSpace{G}
Base.getindex(::ComplexNumbers, d1::Pair{G,Int}, dims::Vararg{Pair{G,Int}}) where {G<:Sector} = RepresentationSpace{G}(d1, dims...)

# Show methods
function Base.show(io::IO, V::RepresentationSpace{G}) where {G<:Sector}
    print(io, "RepresentationSpace{", G, "}(")
    s = sectors(V)
    c = s[1]
    print(io, sprint(showcompact, c), "=>", dim(V,c))
    for k = 2:length(s)
        c = s[k]
        print(io, ", ", sprint(showcompact, c), "=>", dim(V,c))
    end
    print(io, ")")
end

# direct sum of RepresentationSpaces
function ⊕(V1::RepresentationSpace{G}, V2::RepresentationSpace{G}) where {G<:Sector}
    V1.dual == V2.dual || throw(SpaceMismatch("Direct sum of a vector space and its dual do not exist"))
    dims = Dict{G,Int}()
    for c in union(sectors(V1), sectors(V2))
        dims[c] = dim(V1,c) + dim(V2,c)
    end
    return RepresentationSpace(dims, V1.dual)
end

# indices
Base.indices(V::RepresentationSpace) = Base.OneTo(dim(V))
function Base.indices(V::RepresentationSpace{G}, c::G) where {G}
    offset = 0
    if V.dual
        c = conj(c)
    end
    for c′ in sectors(V)
        c′ == c && break
        offset += dim(V, c′)
    end
    return offset + (1:dim(V, c))
end

Spin(j::Union{Int,Rational{Int}}) = RepresentationSpace{SU₂}(j=>1)
fSpin(j::Union{Int,Rational{Int}}) = RepresentationSpace{fSU₂}(j=>1)

# ZNSpace -- more efficient implemenation for small N
#-----------------------------------------------------
"""
    struct ZNSpace{N} <: AbstractRepresentationSpace{ZNIrrep{N}}

Optimized mplementation of a graded `ℤ_N` space, i.e. a complex Euclidean space graded
by the irreps of type `ZNIrrep{N}`
"""
struct ZNSpace{N} <: AbstractRepresentationSpace{ZNIrrep{N}}
    dims::NTuple{N,Int}
    dual::Bool
end

# various constructors for convenience
ZNSpace{N}(dims::NTuple{N,Int}) where {N} = ZNSpace{N}(dims, false)
ZNSpace{N}(dims::Vararg{Int,N}) where {N} = ZNSpace{N}(dims, false)
ZNSpace(dims::NTuple{N,Int}) where {N} = ZNSpace{N}(dims)
ZNSpace(dims::Vararg{Int,N}) where {N} = ZNSpace{N}(dims)

ZNSpace{N}(dims::Vararg{Pair{<:Any,Int}}) where {N} = ZNSpace{N}((ZNIrrep{N}(c)=>d for (c,d) in dims)...)
function ZNSpace{N}(dims::Vararg{Pair{ZNIrrep{N},Int}}) where {N}
    newdims = zeros(Int, N)
    @inbounds for (c,d) in dims
        newdims[c.n+1] = d
    end
    return ZNSpace{N}(ntuple(n->newdims[n], Val(N)))
end

Base.getindex(::ComplexNumbers, ::Type{ZNIrrep{N}}) where {N} = ZNSpace{N}
Base.getindex(::ComplexNumbers, d1::Pair{ZNIrrep{N},Int}, dims::Vararg{Pair{ZNIrrep{N},Int}}) where {N} = ZNSpace{N}(d1, dims...)

# properties
function sectors(V::ZNSpace{N}) where {N}
    s = ntuple(n->ZNIrrep{N}(n-1), Val(N))
    V.dual ? map(conj, s) : s
end

dim(V::ZNSpace) = sum(V.dims)
dim(V::ZNSpace{N}, c::ZNIrrep{N}) where {N} = V.dims[c.n+1]

Base.conj(V::ZNSpace{N}) where {N} = ZNSpace{N}(V.dims, !V.dual)

# indices
Base.indices(V::ZNSpace) = Base.OneTo(dim(V))
function Base.indices(V::ZNSpace{N}, c::ZNIrrep{N}) where {N}
    dims = V.dims
    n = V.dual ? conj(c).n : c.n
    offset = 0
    @inbounds for k = 1:n
        offset += dims[k]
    end
    return offset + (1:dims[n+1])
end

const ParitySpace = ZNSpace{2}
const ℤ₂Space = ZNSpace{2}
const ℤ₃Space = ZNSpace{3}
const ℤ₄Space = ZNSpace{4}

# Show methods
Base.show(io::IO, ::Type{ZNSpace{2}}) = print(io, "ℤ₂Space")
Base.show(io::IO, ::Type{ZNSpace{3}}) = print(io, "ℤ₃Space")
Base.show(io::IO, ::Type{ZNSpace{4}}) = print(io, "ℤ₄Space")

Base.show(io::IO, V::ZNSpace{2}) = print(io, "ℤ₂Space", V.dims, V.dual ? "'" : "")
Base.show(io::IO, V::ZNSpace{3}) = print(io, "ℤ₃Space", V.dims, V.dual ? "'" : "")
Base.show(io::IO, V::ZNSpace{4}) = print(io, "ℤ₄Space", V.dims, V.dual ? "'" : "")
Base.show(io::IO, V::ZNSpace{N}) where {N} = print(io, "ZNSpace{", N, "}", V.dims, V.dual ? "'" : "")