text:
Base.IndexStyle(::Type{T}) where {T<:DataView} = IndexLinear()
Base.size(A::DataView) = (length(A),)
Base.lastindex(A::DataView) = length(A)
getobs(A::DataView) = map(getobs, A)
getobs(A::DataView, i) = getobs(A[i])
getobs(A::DataView{T}) where {T<:Tuple} = map(i->getobs(A,i), 1:length(A))
getobs(A::DataView{T}, i::Integer) where {T<:Tuple} = map(getobs, A[i])

allowcontainer(fun, ::AbstractObsView) = true
allowcontainer(fun, ::DataView) = false

# for proper dispatch to trump the abstract arrays one
for T in (ObsDim.Constant, ObsDim.Last, Tuple)
    @eval function nobs(A::DataView, obsdim::$T)
        @assert obsdim === default_obsdim(A)
        nobs(A)
    end
    @eval function getobs(A::DataView, idx, obsdim::$T)
        @assert obsdim === default_obsdim(A)
        getobs(A, idx)
    end
end

# --------------------------------------------------------------------

"""
    ObsView(data, [obsdim])

Description
============

Create a view of the given `data` in the form of a vector of
individual observations. Any data access is delayed until
`getindex` is called, and even `getindex` returns the result of
[`datasubset`](@ref) which in general avoids data movement until
[`getobs`](@ref) is called.

If used as an iterator, the view will iterate over the dataset
once, effectively denoting an epoch. Each iteration will return a
lazy subset to the current observation.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Methods
========

Aside from the `AbstractArray` interface following additional
methods are provided:

- **`getobs(data::ObsView, indices::AbstractVector)`** :
    Returns a `Vector` of indivdual observations specified by
    `indices`.

- **`nobs(data::ObsView)`** :
    Returns the number of observations in `data` that the
    iterator will go over.

Details
========

For `ObsView` to work on some data structure, the type of the
given variable `data` must implement the data container
interface. See `?DataSubset` for more info.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

```julia
X, Y = MLDataUtils.load_iris()

A = obsview(X)
@assert typeof(A) <: ObsView <: AbstractVector
@assert eltype(A) <: SubArray{Float64,1}
@assert length(A) == 150 # Iris has 150 observations
@assert size(A[1]) == (4,) # Iris has 4 features

for x in obsview(X)
    @assert typeof(x) <: SubArray{Float64,1}
end

# iterate over each individual labeled observation
for (x,y) in obsview((X,Y))
    @assert typeof(x) <: SubArray{Float64,1}
    @assert typeof(y) <: String
end

# same but in random order
for (x,y) in obsview(shuffleobs((X,Y)))
    @assert typeof(x) <: SubArray{Float64,1}
    @assert typeof(y) <: String
end
```

see also
=========

[`eachobs`](@ref), [`BatchView`](@ref), [`shuffleobs`](@ref),
[`getobs`](@ref), [`nobs`](@ref), [`DataSubset`](@ref)
"""
struct ObsView{TElem,TData,O} <: AbstractObsView{TElem,TData}
    data::TData
    obsdim::O
end

function ObsView(data::T, obsdim::O) where {T,O}
    E = typeof(datasubset(data, 1, obsdim))
    ObsView{E,T,O}(data,obsdim)
end

function ObsView(A::T, obsdim) where T<:DataView
    @assert obsdim == A.obsdim
    @warn string("Trying to nest a ", T.name, " into an ObsView, which is not supported. Returning ObsView(parent(_)) instead")
    ObsView(parent(A), obsdim)
end

function ObsView(A::ObsView, obsdim)
    @assert obsdim == A.obsdim
    A
end

ObsView(data; obsdim = default_obsdim(data)) =
    ObsView(data, convert(LearnBase.ObsDimension,obsdim))

nobs(A::ObsView) = nobs(A.data, A.obsdim)
Base.parent(A::ObsView) = A.data
Base.length(A::ObsView) = nobs(A)
Base.getindex(A::ObsView, i::Int) = datasubset(A.data, i, A.obsdim)
Base.getindex(A::ObsView, i::AbstractVector) =
    ObsView(datasubset(A.data, i, A.obsdim), A.obsdim)

# compatibility with nested functions
default_obsdim(A::ObsView) = A.obsdim

@doc (@doc ObsView)
const obsview = ObsView

function Base.showarg(io::IO, A::ObsView, toplevel)
    print(io, "obsview(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, replace(string(A.obsdim), "LearnBase." => ""))
    print(io, ')')
    toplevel && print(io, " with eltype ", nameof(eltype(A))) # simpify
end

# --------------------------------------------------------------------

default_batch_size(source, obsdim) = clamp(div(nobs(source,obsdim), 5), 2, 100)

"""
Helper function to compute sensible and compatible values for the
`size` and `count`
"""
function _compute_batch_settings(source, size::Int = -1, count::Int = -1, obsdim = default_obsdim(source), upto = false)
    num_observations = nobs(source, obsdim)::Int
    @assert num_observations > 0
    if upto && size > num_observations
        size = num_observations
    end
    size  <= num_observations || throw(ArgumentError("Specified batch-size is too large for the given number of observations"))
    count <= num_observations || throw(ArgumentError("Specified batch-count is too large for the given number of observations"))
    if size > 0 && upto
        while num_observations % size != 0 && size > 1
            size = size - 1
        end
    end
    if size <= 0 && count <= 0
        # no batch settings specified, use default size and as many batches as possible
        size = default_batch_size(source, obsdim)::Int
        count = floor(Int, num_observations / size)
    elseif size <= 0
        # use count to determine size. try use all observations
        size = floor(Int, num_observations / count)
    elseif count <= 0
        # use size and as many batches as possible
        count = floor(Int, num_observations / size)
    else
        # use count just for boundscheck
        max_batchcount = floor(Int, num_observations / size)
        count <= max_batchcount || throw(ArgumentError("Specified number of partitions is too large for the specified size"))
        count = max_batchcount
    end

    # check if the settings will result in all data points being used
    unused = num_observations - size*count
    if unused > 0
        @warn "The specified values for size and/or count will result in $unused unused data points" maxlog=1
    end
    size::Int, count::Int
end

"""
Helper function to translate a batch-index into a range of observations.
"""
function _batchrange(batchsize::Int, batchindex::Int)
    startidx = (batchindex - 1) * batchsize + 1
    startidx:(startidx + batchsize - 1)
end

# --------------------------------------------------------------------

"""
    BatchView(data, [size|maxsize], [count], [obsdim])

Description
============

Create a view of the given `data` that represents it as a vector
of batches. Each batch will contain an equal amount of
observations in them. The number of batches and the batch-size
can be specified using (keyword) parameters `count` and `size`.
In the case that the size of the dataset is not dividable by the
specified (or inferred) `size`, the remaining observations will
be ignored with a warning.

Note that any data access is delayed until `getindex` is called,
and even `getindex` returns the result of [`datasubset`](@ref)
which in general avoids data movement until [`getobs`](@ref) is
called.

If used as an iterator, the object will iterate over the dataset
once, effectively denoting an epoch. Each iteration will return a
mini-batch of constant [`nobs`](@ref), which effectively allows
to iterator over [`data`](@ref) one batch at a time.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`size`** : The batch-size of each batch. I.e. the number of
    observations that each batch must contain.

- **`maxsize`** : The maximum batch-size of each batch. I.e. the
    number of observations that each batch should contain. If the
    number of total observation is not divideable by the size it
    will be reduced until it is.

- **`count`** : The number of batches that the view will contain.

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Methods
========

Aside from the `AbstractArray` interface following additional
methods are provided.

- **`getobs(data::BatchView, batchindices)`** :
    Returns a `Vector` of the batches specified by `batchindices`.

- **`nobs(data::BatchView)`** :
    Returns the total number of observations in `data`. Note that
    unless the batch-size is 1, this number will differ from
    `length`.

Details
========

For `BatchView` to work on some data structure, the type of the
given variable `data` must implement the data container
interface. See `?DataSubset` for more info.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

```julia
using MLDataUtils
X, Y = MLDataUtils.load_iris()

A = batchview(X, size = 30)
@assert typeof(A) <: BatchView <: AbstractVector
@assert eltype(A) <: SubArray{Float64,2}
@assert length(A) == 5 # Iris has 150 observations
@assert size(A[1]) == (4,30) # Iris has 4 features

# 5 batches of size 30 observations
for x in batchview(X, size = 30)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert nobs(x) === 30
end

# 7 batches of size 20 observations
# Note that the iris dataset has 150 observations,
# which means that with a batchsize of 20, the last
# 10 observations will be ignored
for (x,y) in batchview((X,Y), size = 20)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
    @assert nobs(x) === nobs(y) === 20
end

# 10 batches of size 15 observations
for (x,y) in batchview((X,Y), maxsize = 20)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
    @assert nobs(x) === nobs(y) === 15
end

# randomly assign observations to one and only one batch.
for (x,y) in batchview(shuffleobs((X,Y)))
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
end

# iterate over the first 2 batches of 15 observation each
for (x,y) in batchview((X,Y), size=15, count=2)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
    @assert size(x) == (4, 15)
    @assert size(y) == (15,)
end
```

see also
=========

[`eachbatch`](@ref), [`ObsView`](@ref), [`shuffleobs`](@ref),
[`getobs`](@ref), [`nobs`](@ref), [`DataSubset`](@ref)
"""
struct BatchView{TElem,TData,O} <: AbstractBatchView{TElem,TData}
    data::TData
    size::Int
    count::Int
    obsdim::O
end

function BatchView(data::T, size::Int, count::Int, obsdim::O = default_obsdim(data), upto::Bool = false) where {T,O}
    nsize, ncount = _compute_batch_settings(data, size, count, obsdim, upto)
    E = typeof(datasubset(data, 1:nsize, obsdim))
    BatchView{E,T,O}(data, nsize, ncount, obsdim)
end

function BatchView(data::T, size::Int, obsdim::O = default_obsdim(data), upto::Bool = false) where {T,O<:Union{Tuple,ObsDimension}}
    BatchView(data, size, -1, obsdim, upto)
end

function BatchView(A::BatchView, size::Int, count::Int, obsdim, upto::Bool = false)
    @assert obsdim == A.obsdim
    BatchView(parent(A), size, count, obsdim, upto)
end

BatchView(data, obsdim::Union{Tuple,ObsDimension}) =
    BatchView(data, -1, -1, obsdim)

function BatchView(data; size = -1, maxsize = -1, count = -1, obsdim = default_obsdim(data))
    maxsize != -1 && size != -1 && throw(ArgumentError("Providing both \"size\" and \"maxsize\" is not supported"))
    if maxsize != -1
        # set upto to true in order to allow a flexible batch size
        BatchView(data, maxsize, count, convert(LearnBase.ObsDimension,obsdim), true)
    else
        # force given batch size
        BatchView(data, size, count, convert(LearnBase.ObsDimension,obsdim))
    end
end

allowcontainer(::typeof(shuffleobs), ::BatchView) = true
allowcontainer(::typeof(splitobs), ::BatchView) = true

"""
    batchsize(data) -> Int

Return the fixed size of each batch in `data`.
"""
batchsize(A::BatchView) = A.size
nobs(A::BatchView) = A.count * A.size
Base.parent(A::BatchView) = A.data
Base.length(A::BatchView) = A.count
Base.getindex(A::BatchView, batchindex::Int) =
    datasubset(A.data, _batchrange(A.size, batchindex), A.obsdim)

function Base.getindex(A::BatchView, batchindices::AbstractVector)
    obsindices = union((_batchrange(A.size, bi) for bi in batchindices)...)::Vector{Int}
    BatchView(datasubset(A.data, obsindices, A.obsdim), A.size, -1, A.obsdim)
end

# compatibility with nested functions
default_obsdim(A::BatchView) = A.obsdim

@doc (@doc BatchView)
const batchview = BatchView

function Base.showarg(io::IO, A::BatchView, toplevel)
    print(io, "batchview(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, A.size, ", ")
    print(io, A.count, ", ")
    print(io, replace(string(A.obsdim), "LearnBase." => ""))
    print(io, ')')
    toplevel && print(io, " with eltype ", nameof(eltype(A))) # simpify
end

# --------------------------------------------------------------------

# if subsetting a DataView, then DataView the subset instead.
for fun in (:DataSubset, :datasubset), O in (ObsDimension, Tuple)
    @eval @generated function ($fun)(A::T, i, obsdim::$O) where T<:AbstractObsView
        quote
            @assert obsdim == A.obsdim
            ($(T.name.name))(($($fun))(parent(A), i, obsdim), obsdim)
        end
    end
    @eval function ($fun)(A::BatchView, i, obsdim::$O)
        @assert obsdim == A.obsdim
        length(i) < A.size && throw(ArgumentError("The chosen batch-size ($(A.size)) is greater than the number of observations ($(length(i)))"))
        BatchView(($fun)(parent(A), i, obsdim), A.size, -1, obsdim)
    end
end


text:
module AutoGibbs

include("tracking.jl")
include("dependencies.jl")
include("plotting.jl")
include("conditionals.jl")
include("auto_conditional.jl")
include("static_conditional.jl")
include("crp.jl")

export trackdependencies

function slicedependencies(model::Model{F}, args...) where {F}
    trace = trackmodel(model, args...)
    strip = strip_model_layers(F, trace)
    slice = strip_dependencies(strip)
    return slice
end

function trackdependencies(model, args...)
    slice = slicedependencies(model, args...)
    return makegraph(slice)
end

end # module

ext: jl


text:
import Distributions.quantile
import BlackBoxOptim.Utils.latin_hypercube_sampling

"""
Get quantiles from `p` according to weights provided by
Latin Hypercube Sampling `samples`.

"""
function get_samples(p, samples)
    out = try
        quantile.(p, samples)
    catch
        [quantile(p, s) for s in samples]
    end
    out
end

"""
Generate an ensemble of size `nens` using Latin Hypercube Sampling of the
list of distributions or collections provided as `args`.

"""
function lhs_ensemble(nens, args...)
    nvars = length(args)
    usamples = latin_hypercube_sampling([0. for _ in 1:nvars], [1. for _ in 1:nvars], nens)
    [get_samples(args[i], usamples[i, :]) for i in 1:nvars]
end

"""
Generate a prior ensemble of discharge, roughness coefficient, channel shape parameter,
and downstream bed elevation from provided distributions or sample data.

"""
function prior_ensemble(x, Qₚ, nₚ, rₚ, zₚ, nens)
    ze = zeros(length(x), nens)
    Qe, ne, re, ze[1, :] = lhs_ensemble(nens, Qₚ, nₚ, rₚ, zₚ)
    Qe, ne, re, ze
end

"""
Generate an ensemble of water height profiles from Gradually-Varied-Flow simulations
and associated profiles of bed elevation.

- `H`: water surface elevation
- `W`: water width
- `S`: water surface slope
- `x`: channel chainage
- `hbf`: bankfull water surface elevation
- `wbf`: bankfull width
- `Qe`: ensemble discharge
- `ne`: ensemble roughness coefficient
- `re`: ensemble channel shape parameter
- `ze`: ensemble bed elevation profiles
- `perturb_slope`: treat water surface slope stochastically (default is false)

"""
function gvf_ensemble!(H, W, S, x, hbf, wbf, Qe, ne, re, ze)
    nens = length(Qe)
    if ndims(S) > 1
        Se = S
    else
        Se = repeat(S', outer=nens)'
    end
    for j in 2:length(x)
        ze[j, :] = ze[j-1, :] .+ Se[j, :] .* (x[j] - x[j-1]);
    end
    ybf = hbf .- ze
    he = zeros(length(x), nens)
    for i in 1:nens
        he[:, i] = gvf(Qe[i], (H[1]-ze[1, i])*re[i]/(re[i] + 1), W[1], Se[:, i], ne[i], x, wbf, ybf[:, i], [re[i] for _ in 1:length(x)])
    end
    he
end


text:
# function related to calculating areas and volumes

@doc """
### calcVolumeContribution!

Returns the contribution from the given faces to the volume of enclosed region.
The Gauss-Divergence theorem is applied to (1/3) ∫∫∫ ∇⋅(x,y,z) dΩ = vol.

**Inputs**

* `sbpface`: an SBP face operator type
* `xsbp`: SBP-face nodes in physical space; [coord, sbp node, face]
* `nrm`: scaled face-normal at the sbpface nodes; [component, sbp node, face]
* `vol`: volume contribution from boundaries other than this set of faces

**Returns**

* `vol`: the updated value for the volume

"""->
function calcVolumeContribution!(sbpface::AbstractFace{Tsbp},
                                 xsbp::AbstractArray{Tmsh,3},
                                 nrm::AbstractArray{Tmsh,3}) where {Tsbp,Tmsh}
  @assert( sbpface.numnodes == size(xsbp,2) == size(nrm,2) )
  @assert( size(xsbp,3) == size(nrm,3) )
  @assert( size(xsbp,1) == size(nrm,1) )
  vol = zero(Tmsh)
  for f = 1:size(nrm,3)
    for i = 1:sbpface.numnodes
      for di = 1:size(nrm,1)# 3
        vol += (1/3)*sbpface.wface[i]*xsbp[di,i,f]*nrm[di,i,f]
      end
    end
  end
  return vol
end

"""
  This function computes the volume contribution of a set of faces, where the
  set of faces is given by which boundary condition they have applied to them.
"""
function calcVolumeContribution!(mesh::AbstractMesh{Tmsh}, eqn::AbstractSolutionData, bndry_nums::Array{Int, 1}) where Tmsh
# compute the volume contribution for the faces associated with a set of
# boundary conditions

  vol = zero(Tmsh)
  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = sview(mesh.bndryfaces, face_range)

    # calculate normal
    nrm = sview(mesh.nrm_bndry, :, :, face_range)
#    nrm = computeNormal(mesh, eqn, bndry_faces, mesh.dxidx_bndry)
    coords = sview(mesh.coords_bndry, :, :, face_range)

    # compute volume contribution
    vol += calcVolumeContribution!(mesh.sbpface, coords, nrm)
  end

  return vol
end

function calcVolumeContribution_rev!(mesh::AbstractMesh{Tmsh}, eqn::AbstractSolutionData, bndry_nums::Array{Int, 1}, vol_bar::Tmsh) where Tmsh

  vol = zero(Tmsh)
  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = sview(mesh.bndryfaces, face_range)
    coords = sview(mesh.coords_bndry, :, :, face_range)
    coords_bar = zeros(coords)  #TODO: make this a field of mesh?

    nrm = sview(mesh.nrm_bndry, :, :, face_range)
    nrm_bar = sview(mesh.nrm_bndry_bar, :, :, face_range)
#    nrm = computeNormal(mesh, eqn, bndry_faces)
#    nrm_bar = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace, length(bndry_faces))

    # compute volume contribution
    calcVolumeContribution_rev!(mesh.sbpface, coords, coords_bar, nrm, nrm_bar, vol_bar)

    # calculate normal
#    computeNormal_rev(mesh, eqn,  bndry_faces, nrm_bar)

  end

  return nothing
end

"""
  Computes the normal vector at each surface node of the given list of faces.
  Straight-sided meshes only!

  Inputs:
    mesh: an AbstractMesh
    sbp: an SBP operator
    bndryfaces: list of boundary faces

  Outputs: nrm: mesh.dim x mesh.numNodesPerFace x length(bndryfaces) array
                containing the normal vectors
"""
function computeNormal(mesh::AbstractMesh{Tmsh}, eqn::AbstractSolutionData,
                       bndryfaces::AbstractArray{Boundary, 1},
                       dxidx::AbstractArray{Tmsh, 4}) where Tmsh

  if typeof(mesh) <: AbstractDGMesh
    @assert mesh.coord_order == 1
  end

  nfaces = length(bndryfaces)
  nrm = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace, nfaces)

  for i=1:nfaces
    bndry_i = bndryfaces[i]
    nrm_xi = sview(mesh.sbpface.normal, :, bndry_i.face)
    for j=1:mesh.numNodesPerFace
      dxidx_j = sview(dxidx, :, :, j, i)
      nrm_xy_j = sview(nrm, :, j, i)
      calcBCNormal(eqn.params, dxidx_j, nrm_xi, nrm_xy_j)
    end
  end

  return nrm
end

"""
  This function uses reverse mode to back propigate perturbations in nrm to
  dxidx_bndry
"""
function computeNormal_rev(mesh::AbstractMesh, eqn::AbstractSolutionData,
                           bndryfaces::AbstractArray{Boundary, 1},
                           nrm_bar::Abstract3DArray{Tmsh},
                           dxidx::AbstractArray{Tmsh, 4},
                           dxidx_bar::AbstractArray{Tmsh, 4}) where Tmsh
  if typeof(mesh) <: AbstractDGMesh
    @assert mesh.coord_order == 1
  end

  nfaces = length(bndryfaces)
  for i=1:nfaces
    bndry_i = bndryfaces[i]
    nrm_xi = sview(mesh.sbpface.normal, :, bndry_i.face)
    for j=1:mesh.numNodesPerFace
      dxidx_j = sview(dxidx, :, :, j, i)
      dxidxbar_j = sview(dxidx_bar, :, :, j, i)
      nrmbar_xy_j = sview(nrm_bar, :, j, i)
      calcBCNormal_revm(eqn.params, dxidx_j, nrm_xi, nrmbar_xy_j, dxidxbar_j)
    end
  end

  return nothing
end

@doc """
### calcVolumeContribution_rev!

This is the reverse differentiated version of calcVolumeContribution!.  See docs
of calcVolumeContribution! for further details of the primal method.  This
function is differentiated with respect to the primal version's `xsbp` and `nrm`
variables.

**Inputs**

* `sbpface`: an SBP face operator type
* `xsbp`: SBP-face nodes in physical space; [coord, sbp node, face]
* `nrm`: scaled face-normal at the sbpface nodes; [component, sbp node, face]
* `vol_bar`: left multiplies d(vol)/d(nrm) and d(vol)/d(xsbp); it is a scalar

**InOuts**

* `xsbp_bar`: result of vector Jacobian product; [component, sbp node, face]
* `nrm_bar`: result of vector Jacobian product; [component, sbp node, face]

"""->
function calcVolumeContribution_rev!(sbpface::AbstractFace{Tsbp},
                                     xsbp::AbstractArray{Tmsh,3},
                                     xsbp_bar::AbstractArray{Tmsh,3},
                                     nrm::AbstractArray{Tmsh,3},
                                     nrm_bar::AbstractArray{Tmsh,3},
                                     vol_bar::Tmsh) where {Tsbp,Tmsh}
  @assert( sbpface.numnodes == size(xsbp,2) == size(xsbp_bar,2) == size(nrm,2)
           == size(nrm_bar,2) )
  @assert( size(xsbp,3) == size(xsbp_bar,3) == size(nrm,3) == size(nrm_bar,3) )
  @assert( size(xsbp,1) == size(xsbp,1) == size(nrm,1) == size(nrm_bar,1) )
  for f = 1:size(nrm,3)
    for i = 1:sbpface.numnodes
      for di = 1:size(nrm,1) # 3
        # vol += (1/3)*sbpface.wface[i]*xsbp[di,i,f]*nrm[di,i,f]
        fac = vol_bar*(1/3)*sbpface.wface[i]
        xsbp_bar[di,i,f] += fac*nrm[di,i,f]
        nrm_bar[di,i,f] += fac*xsbp[di,i,f]
      end
    end
  end
end

@doc """
### calcProjectedAreaContribution!

  Returns the contribution from the given faces to the projected area onto the
  plane of coordiante `di`

  Because this function uses the absolute value of the normal vector, it is not
  differentiable at zero.

**Inputs**

* `sbpface`: an SBP face operator type
* `di`: the index that indicates the projection plane (e.g. `di`=1 --> yz plane)
* `nrm`: scaled face-normal at the sbpface nodes; [component, sbp node, face]
* `projarea`: projected area of faces from other boundaries.

**Returns**

* `projarea`: the updated value for projected area

"""->
function calcProjectedAreaContribution!(sbpface::AbstractFace{Tsbp},
                                        nrm::AbstractArray{Tmsh,3},
                                        di::Int) where {Tsbp,Tmsh}
  @assert( sbpface.numnodes == size(nrm,2) )
  projarea = zero(Tmsh)
  for f = 1:size(nrm,3)
    for i = 1:sbpface.numnodes
      projarea += 0.5*sbpface.wface[i]*abs(nrm[di,i,f])
    end
  end
  return projarea
end

"""
  This function computes the volume contribution of a set of faces, where the
  set of faces is given by which boundary condition they have applied to them.
"""
function calcProjectedAreaContribution!(mesh::AbstractMesh{Tmsh}, eqn::AbstractSolutionData, bndry_nums::Array{Int, 1}, di::Int) where Tmsh

  proj_area = zero(Tmsh)
  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = sview(mesh.bndryfaces, face_range)

    # calculate normal
    nrm = sview(mesh.nrm_bndry, :, :, face_range)
#    nrm = computeNormal(mesh, eqn, bndry_faces)
#    coords = sview(mesh.coords_bndry, :, :, face_range)

    proj_area += calcProjectedAreaContribution!(mesh.sbpface, nrm, di)
  end

  return proj_area
end

@doc """
### calcProjectedAreaContribution_rev!

This is the reverse differentiated version of calcProjectedAreaContribution!.
See docs of calcProjectedAreaContribution! for further details of the primal
method.  This function is differentiated with respect to the primal version's
`nrm` variable.

**Inputs**

* `sbpface`: an SBP face operator type
* `di`: the index that indicates the projection plane (e.g. `di`=1 --> yz plane)
* `nrm`: scaled face-normal at the sbpface nodes; [component, sbp node, face]
* `projarea_bar`: left multiplies d(projarea)/d(nrm); it is a scalar

**InOuts**

* `nrm_bar`: result of vector Jacobian product; [component, sbp node, face]

"""->
function calcProjectedAreaContribution_rev!(sbpface::AbstractFace{Tsbp}, nrm::AbstractArray{Tmsh,3},
    nrm_bar::AbstractArray{Tmsh,3}, di::Int, projarea_bar::Tmsh) where {Tsbp,Tmsh}
  @assert( sbpface.numnodes == size(nrm,2) == size(nrm_bar,2) )
  for f = 1:size(nrm,3)
    for i = 1:sbpface.numnodes
      # projarea += 0.5*sbpface.wface[i]*abs(nrm[di,i,f])
      #  v1 = nrm[di,i,f]
      #  v2 = abs(v1)
      #  v3 = 0.5*sbpface.wface[i]
      #  v4 = v2*v3

      nrm_bar[di,i,f] += (projarea_bar*0.5*sbpface.wface[i]*
                          absvalue_deriv(nrm[di,i,f]))
    end
  end
end

"""
  Back propigate projarea_bar to mesh.dxidx_bndry_bar
"""
function calcProjectedAreaContribution_rev!(mesh::AbstractMesh, eqn::AbstractSolutionData, bndry_nums::Array{Int, 1}, di::Int,  projarea_bar::Tmsh) where Tmsh

  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = sview(mesh.bndryfaces, face_range)

    nrm = sview(mesh.nrm_bndry, :, :, face_range)
    nrm_bar = sview(mesh.nrm_bndry_bar, :, :, face_range)
#    nrm = computeNormal(mesh, eqn, bndry_faces)
#    nrm_bar = zero(nrm)
    calcProjectedAreaContribution_rev!(mesh.sbpface, nrm, nrm_bar, di, projarea_bar)

    # back propigate to dxidx_bndry_bar
#    computeNormal_rev(mesh, eqn, bndry_faces, nrm_bar)
  end

  return nothing
end

function crossProd(x::AbstractArray{T,1}, y::AbstractArray{T,1},
                   cp::AbstractArray{T,1}) where T
  cp[1] = x[2]*y[3] - y[2]*x[3]
  cp[2] = x[3]*y[1] - y[3]*x[1]
  cp[3] = x[1]*y[2] - y[1]*x[2]
end

function crossProd_rev(x::AbstractArray{T,1}, x_bar::AbstractArray{T,1},
                       y::AbstractArray{T,1}, y_bar::AbstractArray{T,1},
                       cp_bar::AbstractArray{T,1}) where T
  #cp[1] = x[2]*y[3] - y[2]*x[3]
  x_bar[2] +=  y[3]*cp_bar[1]
  y_bar[3] +=  x[2]*cp_bar[1]
  y_bar[2] += -x[3]*cp_bar[1]
  x_bar[3] += -y[2]*cp_bar[1]
  #cp[2] = x[3]*y[1] - y[3]*x[1]
  x_bar[3] +=  y[1]*cp_bar[2]
  y_bar[1] +=  x[3]*cp_bar[2]
  y_bar[3] += -x[1]*cp_bar[2]
  x_bar[1] += -y[3]*cp_bar[2]
  #cp[3] = x[1]*y[2] - y[1]*x[2]
  x_bar[1] +=  y[2]*cp_bar[3]
  y_bar[2] +=  x[1]*cp_bar[3]
  y_bar[1] += -x[2]*cp_bar[3]
  x_bar[2] += -y[1]*cp_bar[3]
end


text:
"""
The `Pfam` module, defines functions to measure the protein contact prediction performance of information measure between column pairs from a Pfam MSA.

**Features**

- Read and download Pfam MSAs
- Obtain PDB information from alignment annotations
- Map between sequence/alignment residues/columns and PDB structures
- Measure of AUC (ROC curve) for contact prediction of MI scores

```julia

using MIToS.Pfam
```
"""
module Pfam

  using MIToS.Utils
  using MIToS.MSA
  using MIToS.SIFTS
  using MIToS.PDB
  using MIToS.Information
  using PairwiseListMatrices
  using DataStructures
  using ROCAnalysis

  export  Stockholm, downloadpfam, getseq2pdb,
          msacolumn2pdbresidue, hasresidues,
          msacontacts, msaresidues,
          getcontactmasks, AUC

  include("download.jl")
  include("pdb.jl")

end


text:
# These are entirely functions relating to the REPL capability of Astellarn.

"""
    monkey!(board)

The monkey AI plays random moves. It is a monkey.
"""
function monkey!(board::Board)
    if isdrawbymaterial(board)
        return :DRAW
    end
    if is50moverule(board)
        return :DRAW
    end
    ml = MoveStack(150)
    gen_moves!(ml, board)
    if length(ml) == 0
        if ischeck(board)
            if board.turn == WHITE
                return :BLACK_WINS
            else
                return :WHITE_WINS
            end
        else
            return :DRAW
        end
    end
    apply_move!(board, rand(ml))
    board
end

"""
    engine!(board, depth)

Looks ahead using alpha-beta pruning plus quiescence, alongside a naive evaluation. Selects a move.
"""
function engine!(board::Board, depth::Int)
    if isdrawbymaterial(board)
        return :DRAW
    end
    if is50moverule(board)
        return :DRAW
    end
    thread = Thread()
    copy!(thread.board, board)

    timeman = TimeManagement(100, true, 0, depth, time(), 0, 0, -1, 0, 0)
    thread.timeman = timeman
    thread.stop = false
    thread.ss.nodes = 0
    thread.ss.depth = 0
    thread.ss.seldepth = 0
    thread.ss.tbhits = 0

    ttable = TT_Table()
    eval = find_best_move(thread, ttable, depth)
    move = thread.pv[1][1]
    if move == MOVE_NONE
        if ischeck(board)
            if board.turn == WHITE
                return :BLACK_WINS
            else
                return :WHITE_WINS
            end
        else
            return :DRAW
        end
    end
    apply_move!(board, move)
    board
end

# clear both squares (in case of captures), and set the second square
function user_move!(board::Board, movefrom::String, moveto::String)
    sqr1 = LABEL_TO_SQUARE[movefrom]
    sqr2 = LABEL_TO_SQUARE[moveto]
    if type(board[sqr1]) == PAWN
        if abs(sqr2 - sqr1) > 10
            f = __DOUBLE_PAWN
        elseif abs(sqr2 - sqr1) !== 8
            if board[sqr2] == BLANK
                f = __ENPASS
            else
                f = __NORMAL_MOVE
            end
        else
            f = __NORMAL_MOVE
        end
    else
        f = __NORMAL_MOVE
    end
    ctype = color(board[sqr1])
    if ctype !== board.turn
        @warn "Illegal move, try another"
        return
    end
    move = Move(sqr1, sqr2, f)
    if verify_move(board, move)
        apply_move!(board, Move(sqr1, sqr2, f))
    else
        @warn "Not a legal move"
    end
    return
end
user_move!(board::Board, move::String) = play!(board, move[1:2], move[3:4])

function user_move!(board::Board, movefrom::String, moveto::String, promo::PieceType)
    sqr1 = LABEL_TO_SQUARE[movefrom]
    sqr2 = LABEL_TO_SQUARE[moveto]
    if promo == KNIGHT
        f = __KNIGHT_PROMO
    elseif promo == BISHOP
        f = __BISHOP_PROMO
    elseif promo == ROOK
        f = __ROOK_PROMO
    elseif promo == QUEEN
        f = __QUEEN_PROMO
    end
    ctype = color(board[sqr1])
    if ctype !== board.turn
        @warn "Illegal move, try another"
        return
    end
    move = Move(sqr1, sqr2, f)
    if verify_move(board, move)
        apply_move!(board, Move(sqr1, sqr2, f))
    else
        @warn "Not a legal move"
    end
    return
end

function user_move!(board::Board, castling::Expr)
    if castling.args == [:-,0,0]
        f = __KING_CASTLE
    elseif castling.args == [:-,:(0-0),0]
        f = __QUEEN_CASTLE
    else
        @warn "Input not understood"
    end
    sqr1 = square(kings(board) & friendly(board))
    if f == __KING_CASTLE
        sqr2 = sqr1 - 2
    else
        sqr2 = sqr1 + 2
    end
    ctype = color(board[sqr1])
    if ctype !== board.turn
        @warn "Illegal move, try another"
        return
    end
    move = Move(sqr1, sqr2, f)
    if verify_move(board, move)
        apply_move!(board, Move(sqr1, sqr2, f))
    else
        @warn "Not a legal move"
    end
    return
end

function verify_move(board::Board, move::Move)
    ml = MoveStack(200)
    gen_moves!(ml, board)
    in(move, ml)
end

# edit the global board
macro move(sqr1::Symbol, sqr2::Symbol)
    user_move!(_globalboard, String(sqr1), String(sqr2))
    _globalboard
end

macro move(sqr1::Symbol, sqr2::Symbol, promo::Symbol)
    user_move!(_globalboard, String(sqr1), String(sqr2), eval(promo))
    _globalboard
end

macro move(sqr12::Symbol)
    user_move!(_globalboard, String(sqr12))
    _globalboard
end

macro move(castling::Expr)
    user_move!(_globalboard, castling)
    _globalboard
end

"""
    @board

Create a new blank 'board' obect, accessible as the `global` object `_globalboard`.
"""
macro board()
    global _globalboard = Board()
end

"""
    @newgame

Create a new game from the starting position (accessible if needed as the `global` object `_globalboard`). Moves can be made using the @move macro.
"""
macro newgame()
    global _globalboard = importfen(START_FEN)
end

macro importfen(fen::String)
    global _globalboard = importfen(fen)
end

"""
    @random

Play a random legal move on the global board.
"""
macro random()
    monkey!(_globalboard)
end

"""
    @random

Play the best engine move on the global board.
"""
macro engine()
    engine!(_globalboard, 10)
end

"""
    @engine n

Play the best engine move on the global board, evaluated at a depth 'n::Int'.
"""
macro engine(depth::Int)
    engine!(_globalboard, depth)
end

macro perft(depth::Int)
    perft(_globalboard, depth)
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
# Julia wrapper for header: /usr/include/scip/scip_validation.h
# Automatically generated using Clang.jl wrap_c

function SCIPvalidateSolve(scip, primalreference, dualreference, reftol, quiet, feasible, primalboundcheck, dualboundcheck)
    ccall((:SCIPvalidateSolve, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, UInt32, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, primalreference, dualreference, reftol, quiet, feasible, primalboundcheck, dualboundcheck)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -20.345339284625858 <= q <= 16.479136897808104)

                     add_NL_constraint(m, :(softplus(0.24719801624822457 + 0.29438950671798914*softplus(0.7291821855933653 + 0.701331929085816*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.7293079760660532*softplus(0.05096633991424904 + 0.31501139336769945*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8601808654748551*softplus(0.18255055897839467 + -0.7759686413687503*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.9720568561360068*softplus(0.5199334546165302 + 0.4081309459110076*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + softplus(-0.49402133773229506 + 0.2465550492028581*softplus(0.7291821855933653 + 0.701331929085816*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.17626257051164496*softplus(0.05096633991424904 + 0.31501139336769945*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8688281780619751*softplus(0.18255055897839467 + -0.7759686413687503*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.3904178182819744*softplus(0.5199334546165302 + 0.4081309459110076*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + softplus(-0.6759660684871824 + 0.9311431633625067*softplus(0.7291821855933653 + 0.701331929085816*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.2015839225926066*softplus(0.05096633991424904 + 0.31501139336769945*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.01676614484829253*softplus(0.18255055897839467 + -0.7759686413687503*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.07861065087286123*softplus(0.5199334546165302 + 0.4081309459110076*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + softplus(-0.08101830092609674 + 0.5708385401734484*softplus(0.7291821855933653 + 0.701331929085816*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.9629021826429458*softplus(0.05096633991424904 + 0.31501139336769945*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.4780289946412233*softplus(0.18255055897839467 + -0.7759686413687503*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.25224290978348307*softplus(0.5199334546165302 + 0.4081309459110076*softplus(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*softplus(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*softplus(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*softplus(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Transits.Orbits: trueanom, kepler_solver

# Compute sin_ν, cos_ν without using arctan function directly
function compute_sincos_ν_no_atan(E, ecc; tol=1e-10)
    sin_E, cos_E = sincos(E)
    denom = 1.0 + cos_E

    # Adjust denominator if necessary to avoid dividing by zero
    m = denom > tol
    m_inv = !m
    denom += 1.0 * m_inv

    # Compute sin ν, cos ν
    x = √((1 + ecc) / (1 - ecc)) * sin_E / denom # where ν = 2 arctan(x)
    x² = x * x
    # Apply trig identites:
    #     sincos(arctan x) = (x, 1) ./ √(1 + x²)
    #     sincos(ν) = 2 sin(arctan x)cos(arctan x), 1 - 2 sin²(arctan x)
    denom = 1.0 / (1.0 + x²)
    sin_ν = 2.0 * x * denom
    cos_ν = (1.0 - x²) * denom
    return sin_ν * m , cos_ν * m - 1.0 * m_inv
end

function compute_E_solver(E, ecc)
    M =  E - ecc * sin(E)
    return kepler_solver(M, ecc) # <- E
end

function _compute_vals(E, ecc)
    # Restrict E to the range [-π, π]
    E = rem2pi(Float64(E), RoundNearest)
    # Computed from solver
    E_solver = compute_E_solver(E, ecc)
    sin_ν_solver, cos_ν_solver = compute_sincos_ν_no_atan(E_solver, ecc)
    # Computed directly
    sin_ν, cos_ν = sincos(trueanom(E, ecc))
    return [E_solver, E, sin_ν_solver, sin_ν, cos_ν_solver, cos_ν]
end
# Returns a generator that can be unpacked into:
#     (E_solver,    E_user,
#     sin_ν_solver, sin_ν_user,
#     cos_ν_solver, cos_ν_user)
compute_summary(Es, eccs) = eachrow(mapreduce(_compute_vals, hcat, Es, eccs))

function test_vals(summary)
    (E_solver,     E_user,
     sin_ν_solver, sin_ν_user,
     cos_ν_solver, cos_ν_user) = summary

    @test all(isfinite.(sin_ν_solver))
    @test all(isfinite.(cos_ν_solver))
    @test allclose(E_solver, E_user)
    @test allclose(sin_ν_solver, sin_ν_user)
    @test allclose(cos_ν_solver, cos_ν_user)
end

# Tests from:
# https://github.com/dfm/kepler.py/blob/main/tests/test_kepler.py

@testset "kepler_solver: M, E edge case" begin
    Es = [0.0, 2 * π, -226.2, -170.4]
    eccs = fill(1.0 - 1e-6, length(Es))
    eccs[end] = 0.9939879759519037
    test_vals(compute_summary(Es, eccs))
end

@testset "kepler_solver: Let them eat π" begin
    eccs = range(0.0, 1.0; length=100)[begin:end-1]
    Es = fill(π, length(eccs))
    test_vals(compute_summary(Es, eccs))
end

@testset "kepler_solver: solver" begin
    ecc_range = range(0, 1; length=500)[begin:end-1]
    E_range = range(-300, 300; length=1_001)
    E_ecc_pairs = Iterators.product(E_range, ecc_range)
    Es = reshape(map(x -> x[1], E_ecc_pairs), :, 1)
    eccs = reshape(map(x -> x[2], E_ecc_pairs), :, 1)
    test_vals(compute_summary(Es, eccs))
end


text:
using SafeTestsets
@safetestset "ChaoticEncryption Tests" begin include("test_chaotic_encryption.jl") end


text:
## draw a GtkIDE logo

function logo()

    circle(t, r, c, opt) = layer(x=r[1]*sin.(2π*t) .+ c[1], y=r[2]*cos.(2π*t) .+ c[2], opt...)
    line(t, v, c, opt)   = layer(x=v[1]*t .+ c[1], y=v[2]*t .+ c[2], opt...)

    c0 = Colors.RGBA(0.15, 0.15, 0.15, 0.9)
    c1 = Colors.RGBA(107/255, 171/255, 91/255, 0.9)
    c2 = Colors.RGBA(0.84, 0.4, 0.38, 0.9)
    c3 = Colors.RGBA(0.67, 0.49, 0.75, 0.9)

    xt = range(0, stop=1, length=25)
    xt = @. xt^4 / (0.3^4 + xt^4)

    opt  = (Geom.point, Theme(default_color=c0, highlight_width=0pt, point_size=3pt))
    opt1 = (Geom.point, Theme(default_color=c1, highlight_width=0pt, point_size=3pt))
    opt2 = (Geom.point, Theme(default_color=c2, highlight_width=0pt, point_size=3pt))
    opt3 = (Geom.point, Theme(default_color=c3, highlight_width=0pt, point_size=3pt))

    for maxt in xt

        t = collect(range(0, stop=maxt, length=90))

        spacing = 0.3 + maxt/15
        h = 1.6

        p = plot(
            #G
            circle(t/2  .+ 0.5, (0.6, 0.8), (0, 0.8), opt),
            line(t, (0, 0.6), (0, 0), opt),
            circle(-t/5 .- 0.5, (0.2, 0.2), (0, 0.8), opt),
            #T
            line(t, (0, h*2/3), (spacing, 0), opt),
            line(t, (0.5, 0), (spacing-0.25, h*2/3), opt),
            #K
            line(t, (0, h*2/3), (2*spacing, 0), opt),
            line(t, (1/3, 1/3), (2*spacing, h/4), opt),
            line(t, (1/3, -0.45), (2*spacing+0.05, h/4+0.05), opt),
            # I
            line(t, (0, h), (3*spacing+0.1, 0), opt1),
            # D
            line(t, (0, h), (4*spacing-0.1, 0), opt2),
            circle(t/2, (0.6, 0.8), (4*spacing-0.1, 0.8), opt2),
            # E
            line(t, (0, h), (6*spacing-0.1, 0), opt3),
            line(t, (0.6, 0), (6*spacing-0.1, 0), opt3),
            line(t, (0.5, 0), (6*spacing-0.1, h/2), opt3),
            line(t, (0.6, 0), (6*spacing-0.1, h), opt3),
            #
            Coord.cartesian(fixed=true, xmin=-0.8, xmax=6*spacing+0.8, ymin=-0.2, ymax=1.8),
            Guide.xticks(ticks=nothing), Guide.yticks(ticks=nothing),
            Guide.xlabel(""), Guide.ylabel(""),

        )
        display(p)
        sleep(0.01)

    end
    sleep(0.01)
    display(plot(x=xt, y=sin.(2π*xt), Geom.line))

end


text:
using Test, LinearAlgebra
using DiffEqSensitivity, StochasticDiffEq
using Random

@info "SDE Adjoints"

seed = 100
Random.seed!(seed)

tstart = 0.0
tend = 0.1
dt = 0.005
trange = (tstart, tend)
t = tstart:dt:tend
tarray = collect(t)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

p2 = [1.01,0.87]

# scalar noise
@testset "SDE inplace scalar noise tests" begin
  using DiffEqNoiseProcess

  dtscalar = tend/1e3

  f!(du,u,p,t) = (du .= p[1]*u)
  σ!(du,u,p,t) = (du .= p[2]*u)

  @info "scalar SDE"

  Random.seed!(seed)
  W = WienerProcess(0.0,0.0,0.0)
  u0 = rand(2)

  linear_analytic_strat(u0,p,t,W) = @.(u0*exp(p[1]*t+p[2]*W))

  prob = SDEProblem(SDEFunction(f!,σ!,analytic=linear_analytic_strat),σ!,u0,trange,p2,
    noise=W
    )
  sol = solve(prob,EulerHeun(), dt=dtscalar, save_noise=true)

  @test isapprox(sol.u_analytic,sol.u, atol=1e-4)

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint())

  @show res_sde_u0, res_sde_p

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test isapprox(res_sde_u0, res_sde_u02,  atol=1e-8)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-8)

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test isapprox(res_sde_u0, res_sde_u02,  atol=1e-8)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-8)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=tend/1e2,adaptive=false,sensealg=InterpolatingAdjoint())

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  rtol=1e-4)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  rtol=1e-4)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  rtol=1e-4)

  @show res_sde_u02, res_sde_p2

  function compute_grads(sol, scale=1.0)
    _sol = deepcopy(sol)
    _sol.W.save_everystep = false
    xdis = _sol(tarray)
    helpu1 = [u[1] for u in xdis.u]
    tmp1 = sum((@. xdis.t*helpu1*helpu1))

    Wtmp = [_sol.W(t)[1][1] for t in tarray]
    tmp2 = sum((@. Wtmp*helpu1*helpu1))

    tmp3 = sum((@. helpu1*helpu1))/helpu1[1]

    return [tmp3, scale*tmp3], [tmp1*(1.0+scale^2), tmp2*(1.0+scale^2)]
  end

  true_grads = compute_grads(sol, u0[2]/u0[1])

  @show  true_grads

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  rtol=1e-4)
  @test isapprox(true_grads[2], res_sde_p', atol=1e-4)
  @test isapprox(true_grads[1], res_sde_u0, rtol=1e-4)
  @test isapprox(true_grads[2], res_sde_p2', atol=1e-4)
  @test isapprox(true_grads[1], res_sde_u02, rtol=1e-4)
end

@testset "SDE oop scalar noise tests" begin
  using DiffEqNoiseProcess

  dtscalar = tend/1e3

  f(u,p,t) = p[1]*u
  σ(u,p,t) = p[2]*u

  Random.seed!(seed)
  W = WienerProcess(0.0,0.0,0.0)
  u0 = rand(2)

  linear_analytic_strat(u0,p,t,W) = @.(u0*exp(p[1]*t+p[2]*W))

  prob = SDEProblem(SDEFunction(f,σ,analytic=linear_analytic_strat),σ,u0,trange,p2,
    noise=W
   )
  sol = solve(prob,EulerHeun(), dt=dtscalar, save_noise=true)

  @test isapprox(sol.u_analytic,sol.u, atol=1e-4)

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint())

  @show res_sde_u0, res_sde_p

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=false))

  @test isapprox(res_sde_u0, res_sde_u02,  atol=1e-8)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-8)

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=BacksolveAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test isapprox(res_sde_u0, res_sde_u02,  atol=1e-8)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-8)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=tend/1e2,adaptive=false,sensealg=InterpolatingAdjoint())

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-4)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=false))

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-4)

  @show res_sde_u02, res_sde_p2

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=dtscalar,adaptive=false,sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP()))

  @test isapprox(res_sde_u0, res_sde_u02,  rtol=1e-4)
  @test isapprox(res_sde_p, res_sde_p2,  atol=1e-4)

  @show res_sde_u02, res_sde_p2

  function compute_grads(sol, scale=1.0)
    _sol = deepcopy(sol)
    _sol.W.save_everystep = false
    xdis = _sol(tarray)
    helpu1 = [u[1] for u in xdis.u]
    tmp1 = sum((@. xdis.t*helpu1*helpu1))

    Wtmp = [_sol.W(t)[1][1] for t in tarray]
    tmp2 = sum((@. Wtmp*helpu1*helpu1))

    tmp3 = sum((@. helpu1*helpu1))/helpu1[1]

    return [tmp3, scale*tmp3], [tmp1*(1.0+scale^2), tmp2*(1.0+scale^2)]
  end

  true_grads = compute_grads(sol, u0[2]/u0[1])

  @show  true_grads

  @test isapprox(true_grads[2], res_sde_p', atol=1e-4)
  @test isapprox(true_grads[1], res_sde_u0, rtol=1e-4)
  @test isapprox(true_grads[2], res_sde_p2', atol=1e-4)
  @test isapprox(true_grads[1], res_sde_u02, rtol=1e-4)

end


text:
module MacroTools

using DataStructures, Markdown, Random
export @match, @capture

include("match/match.jl")
include("match/types.jl")
include("match/union.jl")
include("match/macro.jl")

include("utils.jl")
include("structdef.jl")

include("examples/destruct.jl")
include("examples/threading.jl")
include("examples/forward.jl")

const animals = Symbol[]
const animals_file = joinpath(@__DIR__, "..", "animals.txt")

# Load and initialize animals symbols.
_animals = split(read(animals_file, String))
resize!(animals, length(_animals))
animals .= Symbol.(lowercase.(_animals))

function __init__()
    Random.shuffle!(animals)
end

end # module


text:
# HalidePerovskites
# Uses the newly forked off FeynmanKadanoffOsakaHellwarth.jl package
# Codes by Jarvist Moore Frost, 2017
# These codes were developed with Julia 0.5.0, and requires the Optim and Plots packages.

# This file, is based on doing further calculations beyond:
# https://arxiv.org/abs/1704.05404
# Polaron mobility in halide perovskites
# Jarvist Moore Frost
# (Submitted on 18 Apr 2017 [v1])

push!(LOAD_PATH,"../src/") # load module from local directory

const T=300 # Standard temperature we will calculate for

using PolaronMobility

##### load in library routines... #####
# Plot figures with Plots, which defaults to Pyplot backend
#using Plots
#pyplot()
#gr() # GR backend to Plots
#default(grid=false) # No silly dotted grid lines
#default(size=(400,300)) # Small, for two-column output (i.e. PDF/EPS for article)
#default(size=(800,600)) # Good size for small-ish PNGs for slides

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space

# Hold my beer, I'm going to use a Macro.
macro showprint(ex)
    println("Executing: ",ex)
    return :( $ex )
end

polaronmobility(x...) = PolaronMobility.polaronmobility(x..., verbose=true)

#####
# Call simulation

# CsSnX3 X={Cl,Br,I}
# L. Huang, W.Lambrecht - PRB 88, 165203 (2013)
# Dielectric consts, from TABLE VII
# Effective masses from TABLE VI, mh*
const cm1=2.997e10 # cm-1 to Herz
@showprint CsSnCl=polaronmobility(T, 4.80, 29.4,  243cm1, 0.140) # alpha= 1.386311
@showprint CsSnBr=polaronmobility(T, 5.35, 32.4,  183cm1, 0.082) # alpha= 1.094468
@showprint CsSnI= polaronmobility(T, 6.05, 48.2,  152cm1, 0.069) # alpha= 1.020355

# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus
# CsSnBr3.dat:300 0.877382 511.513 356.358 874.499 8.13686 7.31511 12.6976 0.23729 -3.78992 2.15682 0.84436 0.788735 0.185396
# CsSnCl3.dat:300 1.165048 212.034 147.36 272.993 6.50772 5.58522 11.1557 0.357614 -3.5556 2.42611 0.947833 0.181657 0.143969
# CsSnI3.dat:300 0.728755 703.502 487.085 1448.02 9.64022 8.80508 15.4044 0.198691 -4.01438 2.18761 0.859401 0.967371 0.207864

# So for CsSnSi3, a hole-mobility of 487(Hellwarth)-703(Kadanoff). Compares to expt. 400 (transport) -  585 (Hall effect). Ace!

#Ts,Kμs, Hμs, FHIPμs, ks, Ms, As, Bs, Cs, Fs, Taus
#effectivemass=0.12 # the bare-electron band effective-mass.
# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
@showprint MAPIe=polaronmobility(T, 4.5, 24.1, 2.25E12, 0.12)
@showprint MAPIh=polaronmobility(T, 4.5, 24.1, 2.25E12, 0.15)

# CsPbI3: 3.82, 8.27, 2.83 THz (highest freq mode in cubic)
#    ^-- v. quick Kmesh=3x3x3 PBESol Calc, JMF 2017-04-12. Files: jarvist@titanium:~/phonopy-work/2017-03-Scott-Distortions/1000-Dielectric/3x3x3-kmesh/
# Horribly non-converged!
#                         Kmesh=6x6x6; 7.2/12.1
# Kmesh=9x9x9x, Ediff=10^-9;           6.1/12.0, 2.57 THz
@showprint CsPbI=polaronmobility(T, 6.1,6.1+12.0, 2.57E12, 0.12)

function SendnerCrosscheck()
    cm1=2.997e10 # cm-1 to Herz

    # Rob's / Sendner's paper - values extracted form IR measures
    # https://doi.org/10.1039%2Fc6mh00275g
    #Ts,a,MAPI=polaronmobility("Rob-MAPI", 5.0, 33.5, 40*cm1, 0.104)
    #Ts,a,MAPBr=polaronmobility("Rob-MAPBr", 4.7, 32.3, 51*cm1, 0.117)
    #Ts,a,MAPCl=polaronmobility("Rob-MAPCl", 4.0, 29.8, 70*cm1, 0.2)

    # Private communication. It is not well described in the paper but they
    # used one the Hellwarth effective-mode method to reduce the observed IR
    # oscillators down to a single mode. Rob provided these values by email (as
    # below). The Effetive Masses and dielectric constants are from the paper.
    # Combined, this reproduces their mobilities, and the internal w and
    # v parameters (again, email from Rob).
    # For omega we used: MAPbI/Br/Cl = 112.9/149.4/214.0
    @showprint RobMAPI=polaronmobility(10:10:400, 5.0, 33.5, 112.9*cm1, 0.104)
    @showprint RobMAPBr=polaronmobility(10:10:400, 4.7, 32.3,149.4*cm1, 0.117)
    @showprint RobMAPCl=polaronmobility(10:10:400, 4.0, 29.8, 214.0*cm1, 0.2)

    plot(RobMAPI.T,RobMAPI.Hμ,label="(Rob's values) MAPI",markersize=2,marker=:uptriangle,ylim=(0,400))
    plot!(RobMAPBr.T,RobMAPBr.Hμ,label="(Rob's values) MAPBr",markersize=2,marker=:diamond)
    plot!(RobMAPCl.T,RobMAPCl.Hμ,label="(Rob's values) MAPCl",markersize=2,marker=:diamond)
    savefig("Rob-comparison.png")
#    savefig("Rob-comparison.eps")
end

#SendnerCrosscheck()

println("That's me!")


text:
using Mimi

page_years = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
page_year_0 = 2015

function getpageindexfromyear(year)
    i = findfirst(isequal(year), page_years)
    if i == 0
        error("Invalid PAGE year: $year.")
    end
    return i
end

function getperiodlength(year)      # same calculations made for yagg_periodspan in the model
    i = getpageindexfromyear(year)

    if year==page_years[1]
        start_year = page_year_0
    else
        start_year = page_years[i - 1]
    end

    if year == page_years[end]
        last_year = page_years[end]
    else
        last_year = page_years[i + 1]
    end

    return (last_year - start_year) / 2
end

"""Applies undiscounting factor to get the SCC."""
function undiscount_scc(m::Model, year::Int)
    df = m[:EquityWeighting, :df_utilitydiscountfactor]
    consfocus0 = m[:GDP, :cons_percap_consumption_0][1]
    consfocus = m[:GDP, :cons_percap_consumption][:, 1]
    emuc = m[:EquityWeighting, :emuc_utilityconvexity]
    sccii = getpageindexfromyear(year)

    return df[sccii] * ((consfocus[sccii] / consfocus0)^-emuc)
end

@defcomp ExtraEmissions begin
    e_globalCO2emissions = Parameter(index=[time],unit="Mtonne/year")
    pulse_size = Parameter()
    pulse_year = Parameter()
    e_globalCO2emissions_adjusted = Variable(index=[time],unit="Mtonne/year")

    function run_timestep(p, v, d, t)
        if gettime(t) == p.pulse_year
            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t] + p.pulse_size / getperiodlength(p.pulse_year)
        else
            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t]
        end
    end
end

"""
compute_scc(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided Mimi-PAGE model.
If no model is provided, the default model from main_model.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
"""
function compute_scc(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing, pulse_size = 100000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)

    return scc
end

"""
compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided Mimi-PAGE model.
If no model is provided, the default model from main_model.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
"""
function compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing, pulse_size = 100000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
    scc_disaggregated = mm[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] / undiscount_scc(mm.base, year)

    return (scc = scc, scc_disaggregated = scc_disaggregated, mm = mm)
end

"""
get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing)
Returns a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from main_model.get_model() is used as the base model.
Note that the returned MarginalModel has already been run.
"""
function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing, pulse_size = 100000.)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $page_years.") : nothing

    mm = create_marginal_model(m, pulse_size)

    add_comp!(mm.marginal, ExtraEmissions, :extra_emissions; after=:co2emissions)
    connect_param!(mm.marginal, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
    set_param!(mm.marginal, :extra_emissions, :pulse_size, pulse_size)
    set_param!(mm.marginal, :extra_emissions, :pulse_year, year)

    connect_param!(mm.marginal, :CO2Cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted)

    run(mm)

    return mm
end

function compute_scc_mcs(m::Model, samplesize::Int; year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing, pulse_size = 100000.)
    # Setup of location of final results
    scc_results = zeros(samplesize)

    function mc_scc_calculation(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, ignore::Nothing)
        marginal = sim_inst.models[1]
        marg_damages = marginal[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
        scc_results[trialnum] = marg_damages
    end

    # get simulation
    mcs = getsim()

    # Setup models
    eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run

    # Run
    res = run(mcs, mm, samplesize; post_trial_func=mc_scc_calculation)

    scc_results
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module NetworkResilience

#Library packages
using LinearAlgebra
using Statistics
using Distributed
using Distributions
using DifferentialEquations

#Local files
include("netconstructors.jl")
include("cascadedrivers.jl")
include("networkalgorithms.jl")
include("swing.jl")

#Utilities to export
export NetworkRecorder, Edge, Node, Graph, addnode!, addedge!, numedges, numnodes, NetStore
export incidencematrix, getdegrees, sourcesinkvector, smallworldnetwork, sourcesinklocs
export connectedcomponents, adjacencymatrix, weightedlaplacianmatrix
export nodevoltages, edgecurrents, testnetwork, fracture!
export newrecorder, cascaderesults, smallworldcascade, smallworldensemble, smallworldalphas
export smcascadefixedtop, smcascadeclock, sourcessinkclockface, cascadesequence
export singlecascade, smcascadedist, cascadebisection, sequencedstartbisection
export noisysourcesinkvector, partitionadjacency, avcompontentsize, clustersize
export sourcesinkwellmixed, sourcesinkhalfhalf, austriannetwork
export triggerhist, cascadebisectionrandomtrigger
export edgepower, fswing, fsteadystate, mainjac, centralityvstransient, fracture2!
export smcascadeclock2, cascadebisection2, fracture3!, simplecascade, NetStore2
export swingfracture!, swingcascadevalpha, makegammasourcesinkvector
export fracturerand!, cascadebisection3rand, smcascadeclockrand, swingcascadebisection
export CascadeNetworkSwing, criticalcouplingbisection, swingcascadeandtvalpha, svt
export swingcascadeouterbisection, swingcascadebreakdownvsalpha

end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
### lucy.jl ---  Lucy deconvolution
#
# Copyright (C) 2019 Jakub Wronowski.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: richadson, lucy, deconvolution, signal processing
#
# This file is a part of Deconvolution.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file provides the `lucy' function to perform Richardson-Lucy deconvolution.
#
### Code:

export lucy

function lucy(observed::AbstractArray, psf::AbstractArray; iterations::Integer = 1)
    @assert size(observed) == size(psf)
    @assert iterations >= 0

    psf_ft = fft(psf)
    psf_ft_conj = conj.(psf_ft)

    function lucystep(e)
        return e .* real(ifft(fft(observed ./ ifft(fft(e) .* psf_ft)) .* psf_ft_conj))
    end

    estimated = observed
    for i in 1:iterations
        estimated = lucystep(estimated)
    end

    return estimated
end

"""
    lucy(observed, psf[, iterations])

Returns the Richardson-Lucy deconvolution of `observed`, using the point spread
function `psf`.

All arguments must be arrays in the time/space domain and all of the same size,
they will be converted into the frequency domain internally using the `fft`
function.
"""
lucy


text:
using XRayLabTool
using Test

@testset "XRayLabTool.jl" begin
    data = Refrac(["SiO2","H2O"],[1,2,3.5,5,6],[2.2,1.0])
    si = SubRefrac("Si", [20.0], 2.33)
    @test data["SiO2"]["Dispersion"][3] - 3.814915113516378e-5 < eps()
    @test data["SiO2"]["f1"][1] - 29.448492570232553 < eps()
    @test data["SiO2"]["reSLD"][5] - 1.8979332462750346e-5 < eps()
    @test si["Dispersion"][1] - 1.2096870811637957e-6 < eps()
    @test si["f1"][1] - 14.048053002299712 < eps()
    @test si["reSLD"][1] - 1.977791074150516e-5 < eps()
end


text:
module GeometryWrappersTests

using Gridap
using UnstructuredGrids: RefCell, UGrid
using Test
using UnstructuredGrids: generate_dual_connections
using UnstructuredGrids: find_cell_to_faces

# Polytope to UnstructuredGrid

t = TET_AXIS
polytope = Polytope((t,t,t))

grid = Grid(polytope,1)
grid = Grid(polytope,2)
grid = Grid(polytope,3)
#writevtk(grid,"grid")

t = HEX_AXIS
polytope = Polytope((t,t,t))

grid = Grid(polytope,1)
grid = Grid(polytope,2)
grid = Grid(polytope,3)
#writevtk(grid,"grid")

trian = Triangulation(grid)
quad = CellQuadrature(trian,degree=2)

q = coordinates(quad)
phi = CellGeomap(trian)
x = evaluate(phi,q)

#writevtk(x,"x")

# Polytope to RefCell

t = HEX_AXIS
polytope = Polytope((t,t,t))

refcell = RefCell(polytope)

# UnstructuredGrid to UGrid

c2v = [[1,2,4,5],[2,3,5,6]]
v2x = Point{2,Float64}[(0.,0.),(0.,1.),(0.,2.),(1.,0.),(1.,1.),(1.,2.)]
l = length(c2v)
order = 1
t = (HEX_AXIS,HEX_AXIS)
c2t = ConstantCellValue(t,l)
c2o = ConstantCellValue(order,l)
c2v_data, c2v_ptrs = generate_data_and_ptrs(c2v)

grid = UnstructuredGrid(v2x,c2v_data,c2v_ptrs,c2t,c2o)
ugrid = UGrid(grid)

graph = FullGridGraph(grid)
test_full_grid_graph(graph,2)

t = TET_AXIS
polytope = Polytope((t,t,t))

grid = Grid(polytope,2)

# Others

model = CartesianDiscreteModel(partition=(4,4,3))
graph = FullGridGraph(model)
edge_to_vertices = connections(graph,1,0)
vertex_to_edges = connections(graph,0,1)
r = generate_dual_connections(edge_to_vertices)
@test r == vertex_to_edges

grid2 = Grid(model,2)
vertex_to_edges = connections(graph,0,1)
face_to_edges = connections(graph,2,1)
r = find_cell_to_faces(grid2,vertex_to_edges,1)
@test r == face_to_edges

end # module


text:
using Revise
using Plots
using Pkg
# Pkg.activate(".")
using BenchmarkTools
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer.SolarModel
# using DelimitedFiles
# using Interact
# using Interpolations
# using InstrumentOperator
# using ForwardDiff

## albedo = 0

R_0 = [0.03916, 0.04235, 0.04378, 0.04593, 0.04940, 0.05438, 0.05816, 0.06541, 0.07639, 0.08773, 0.09549, 0.11842, 0.13619, 0.17836, 0.22393, 0.27069]

R_90 = [0.03916, 0.03922, 0.03930, 0.03952, 0.04021, 0.04208, 0.04409, 0.04894, 0.05748, 0.06698, 0.07367, 0.09388, 0.10977, 0.14790, 0.18964, 0.23441]

R_180 = [0.03916, 0.03634, 0.03533, 0.03414, 0.03316, 0.03381, 0.03557, 0.04073, 0.05062, 0.06186, 0.06981, 0.09387, 0.11280, 0.15826, 0.20833, 0.26423]

## albedo = 0.25

# R_0 = [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473]

# R_90 = [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473]

# R_180 = [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473]

R_trues = [R_0, R_90, R_180]

R_modeled = zeros(16, 3)

for i in 1:3 #[0, 90, 180]

    # Load parameters from file
    # parameters = vSmartMOM.default_parameters()
    parameters = vSmartMOM.parameters_from_yaml("/home/rjeyaram/RadiativeTransfer/test/helper/6SV1_1.yaml");

    parameters.vaz = repeat([[0, 90, 180][i]], 16)

    model = model_from_parameters(parameters);

    @show model.τ_rayl[1]
    model.τ_rayl[1] .= 0.1

    R_modeled[:,i] = vSmartMOM.rt_run(model, i_band=1)[:,:,1]

end

p1 = plot(parameters.vza, R_modeled[:,1], xlabel = "VZA", title = "AZ=0", legend=false)
p1 = plot!(parameters.vza, R_trues[3])

p2 = plot(parameters.vza, R_modeled[:,2], xlabel = "VZA", title = "AZ=90", legend=false)
p2 = plot!(parameters.vza, R_trues[2])

p3 = plot(parameters.vza, R_modeled[:,3], xlabel = "VZA", title = "AZ=180", legend=false)
p3 = plot!(parameters.vza, R_trues[1])

p4 = plot(parameters.vza, (R_modeled[:,1]-R_trues[3])./R_trues[3], xlabel = "VZA", title = "AZ=0", legend=false)

p5 = plot(parameters.vza, (R_modeled[:,2]-R_trues[2])./R_trues[2], xlabel = "VZA", title = "AZ=90", legend=false)

p6 = plot(parameters.vza, (R_modeled[:,3]-R_trues[1])./R_trues[1], xlabel = "VZA", title = "AZ=180", legend=false)
plot(p1, p2, p3, p4, p5, p6, layout=(2,3), size=(1200,900))

##

# [0.0000, 11.4783, 16.2602, 23.0739, 32.8599, 43.9455, 50.2082, 58.6677, 66.4218, 71.3371, 73.7398, 78.4630, 80.7931, 84.2608, 86.5602, 88.854]

# "/home/rjeyaram/RadiativeTransfer/test/helper/6SV1_1.yaml"

# params_dict = YAML.load_file("RadiativeTransfer/src/vSmartMOM/ModelParameters/DefaultParameters.yaml")

# vSmartMOM.validate_yaml_parameters(params_dict)

# file = "RadiativeTransfer/src/vSmartMOM/ModelParameters/SampleProfile.yaml"
# profile_nc4 = vSmartMOM.read_atmos_profile(parameters.file, parameters.lat, parameters.lon, parameters.timeIndex);
# profile_yaml = vSmartMOM.read_atmos_profile(file);

# Sets all the "specific" parameters
# parameters = vSmartMOM.default_parameters();

# Generates all the derived attributes from above parameters

#model.params.architecture = RadiativeTransfer.Architectures.GPU();

# function run_auto(x)
#     (x)
#     (x[1] isa ForwardDiff.Dual)
#     (x[1])
#     (typeof(parameters.μ))
#     parameters.μ  = [x[1]]
#     parameters.σ  = [x[2]]
#     parameters.nᵣ = [x[3]]
#     parameters.nᵢ = [x[4]]
#     model = model_from_parameters(parameters);
#     # model.architecture = RadiativeTransfer.Architectures.CPU()
#     model.params.architecture = RadiativeTransfer.Architectures.CPU()
#     R = vSmartMOM.rt_run(model, i_band=1)
#     return R[1,1,:]
# end

# x = [1.3, 2.0, 1.3, 0.000001];

# # @time run_auto(x)
# @time dfdx = ForwardDiff.jacobian(run_auto, x);

# a = 10

# g = x -> ForwardDiff.jacobian(run_auto, x);

# ν_grid = collect((1e7/778):0.015:(1e7/755))

# T = 5777
# black_body = planck_spectrum(T, ν_grid)

# solar = readdlm("RadiativeTransfer/src/solar_merged_20160127_600_26316_100.out")

# solar_idx_start = argmin(abs.(solar[:, 1] .- minimum(ν_grid)))
# solar_idx_end   = argmin(abs.(solar[:, 1] .- maximum(ν_grid)))

# solar_new = solar[(solar_idx_start-10):(solar_idx_end+10), :]

# itp = LinearInterpolation(solar_new[:, 1],
#                           solar_new[:, 2])

# solar_interpolated = itp.(black_body[:,1])

# plot(black_body[:,1], (solar_interpolated .* black_body[:,2]))
# plot!(black_body[:,1], black_body[:,2])

# solar_out = solar_interpolated .* black_body[:,2]

# earth_out = solar_out .* R[1,1,:]

# file = "/home/cfranken/oco2_L1bScND_18688a_180105_B8100r_180206190633.h5"

# ils_Δ, ils_in, dispersion = InstrumentOperator.read_ils_table(file, "/home/rjeyaram/RadiativeTransfer/src/vSmartMOM/ils_oco2.json");

# ils_pixel =

# oco2_kernel = VariableKernelInstrument(ils_pixel, ν, ind_out)

# plot!(1e6 ./ νs[argmin(abs.((νs ./ 100) .- solar[end,1])):(end-1)], solar_interpolated * 1e13)

# 1 ./ νs[argmin(abs.((νs ./ 100) .- solar[end,1])):(end-1)]

# a = ones(1500)
# a[380:(end-1)] .= solar_interpolated

# plot(1:1500, a .* Ls)

# [0.03916, 0.04235, 0.04378, 0.04593, 0.04940, 0.05438, 0.05816, 0.06541, 0.07639, 0.08773, 0.09549, 0.11842, 0.13619, 0.17836, 0.22393, 0.27069]


text:
using NLsolve

ϕ = 0.5
y = 1
β = 0.996
θ = 0.199
α = 0.4
s = 0.0225

function f!(F, v)
    ω = v[1]
    c = v[2]
    F[1] = (1-(1-ϕ)*α)*ω - ϕ*(y + θ*c)
    F[2] = θ^α*(1-β*(1-s))*c - β*(y-ω)
end

soln = nlsolve(f!, ones(2), autodiff = :forward)

@show soln.zero


text:
#Here are parameters for scenarios 1, 2, 3, 4, 5 for Project 2
#For convenience they are stored in a struct, NetworkParameters

using Parameters #You need to install the Parameters.jl package: https://github.com/mauro3/Parameters.jl
using LinearAlgebra

#The @with_kw macro comes from the Parameters package
@with_kw struct NetworkParameters
    L::Int
    gamma_scv::Float64 #This is constant for all scenarios at 3.0
    λ::Float64 #This is undefined for the scenarios since it is varied
    η::Float64 #This is assumed constant for all scenarios at 4.0
    μ_vector::Vector{Float64} #service rates
    P::Matrix{Float64} #routing matrix
    Q::Matrix{Float64} #overflow matrix
    p_e::Vector{Float64} #external arrival distribution
    K::Vector{Int} #-1 means infinity
end

############################
scenario1 = NetworkParameters(  L=3,
                                gamma_scv = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(3),
                                P = [0 1.0 0;
                                    0 0 1.0;
                                    0 0 0],
                                Q = zeros(3,3),
                                p_e = [1.0, 0, 0],
                                K = fill(5,3))
@show scenario1

############################
scenario2 = NetworkParameters(  L=3,
                                gamma_scv = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(3),
                                P = [0 1.0 0;
                                    0 0 1.0;
                                    0.5 0 0],
                                Q = zeros(3,3),
                                p_e = [1.0, 0, 0],
                                K = fill(5,3))
@show scenario2

############################
scenario3 = NetworkParameters(  L=3,
                                gamma_scv = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(3),
                                P = [0 1.0 0;
                                    0 0 1.0;
                                    0.5 0 0],
                                Q = [0 0.5 0;
                                     0 0 0.5;
                                     0.5 0 0],
                                p_e = [1.0, 0, 0],
                                K = fill(5,3))
@show scenario3

############################
scenario4 = NetworkParameters(  L=5,
                                gamma_scv = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = collect(5:-1:1),
                                P = [0   0.5 0.5 0   0;
                                     0   0   0   1   0;
                                     0   0   0   0   1;
                                     0.5 0   0   0   0;
                                     0.2 0.2 0.2 0.2 0.2],
                                Q = [0 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0],
                                p_e = [0.2, 0.2, 0, 0, 0.6],
                                K = [-1, -1, 10, 10, 10])
@show scenario4

############################
scenario5 = NetworkParameters(  L=20,
                                gamma_scv = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(20),
                                P = zeros(20,20),
                                Q = diagm(3=>0.8*ones(17), -17=>ones(3)), #Fixed on Oct 17
                                p_e = vcat(1,zeros(19)),
                                K = fill(5,20))
@show scenario5


text:
#=
  Young tableaux and partition in Julia ConstraintSolver.jl

  See
  http://mathworld.wolfram.com/YoungTableau.html
  and
  http://en.wikipedia.org/wiki/Young_tableau
  """
  The partitions of 4 are
   {4}, {3,1}, {2,2}, {2,1,1}, {1,1,1,1}

  And the corresponding standard Young tableaux are:

  1.   1 2 3 4

  2.   1 2 3         1 2 4    1 3 4
       4             3        2

  3.   1 2           1 3
       3 4           2 4

  4    1 2           1 3      1 4
       3             2        2
       4             4        3

  5.   1
       2
       3
       4
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# Note: count_ctr is defined in constraints_utils.jl
#
function young_tableaux(n=5)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> true,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                        ))

    println("Young tableaux and partitions of order $n")

    # the Young tableuax grid
    @variable(model, 1 <= x[1:n,1:n] <= n+1, Int)

    # the partition structure p
    @variable(model, 0 <= p[1:n] <= n, Int)

    # 1..N is used exactly once (N+1 may be used many times)
    for i in 1:n
        # This is how I would like to model it:
        # @constraint(model,(1 == sum(x[:]==i)))
        # @constraint(model, 1 == sum([x[j,k] == i for j in 1:n, k in 1:n] ))
        # @constraint(model, 1 == sum([xflat[j] == i for j in 1:length(xflat)] ))
        count_ctr(model, x,:(==), i, 1)

    end

    # First position
    @constraint(model, x[1,1] == 1)

    # All rows and columns should be ordered
    for i in 1:n
        increasing(model,x[i,:]) # rows
        increasing(model,x[:,i]) # columns
    end

    # Calculate the structure (the partition)
    for i in 1:n
        # This is what I would like to write
        # @constraint(model, p[i] == sum([x[i,j] for j in 1:n]))
        # This works:
        count_ctr(model,x[i,:],:(<=),n, p[i])

    end

    # P should be ordered
    decreasing(model,p)

    # We just work with n numbers
    @constraint(model, n == sum(p))

    println("solve")
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL

        num_sols = MOI.get(model, MOI.ResultCount())
        println("\nnum_sols:$num_sols\n")

        for sol in 1:num_sols
            println("\nsolution #$sol")
            xx = convert.(Integer,JuMP.value.(x,result=sol))
            pp = convert.(Integer,JuMP.value.(p,result=sol))
            for i in 1:n
                for j in 1:n
                    if xx[i,j] <= n
                        print(xx[i,j]," ")
                    else
                        print("_ ")
                    end
                end
                println()
            end
            println("p:$pp")

        end
    end
end

young_tableaux(5)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module MyPkgTests

using Base.Test
using StaticArrays
import CollisionDetection

include("test_core.jl")
include("test_bloated.jl")
include("test_searcheq.jl")

end # module

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection