text:
##################################################
"""
    symbols(name(s), assumptions...)

Calls `sympy.symbols` to produce symbolic variables and symbolic functions. An alternate to the recommended `@syms`, (when applicable)

In sympy `sympy.symbols` and `sympy.Symbol` both allow the construction of symbolic variables and functions. The `Julia` function `symbols` is an alias for `sympy.symbols`.

* Variables are created through `x=symbols("x")`;
* Assumptions on variables by `x=symbols("x", real=true)`;
* Multiple symbols `x1,x2 = symbols("x[1:3]")` can be created. Unlike `@syms`, the number of variables can be specified with a variable through interpolation.
* Symbolic functions can be created py passing `cls=sympy.Function`, `symbols("F", cls=sympy.Function, real=true)`

"""
function symbols(x::AbstractString; kwargs...)
    out = sympy.symbols(x; kwargs...)
end
symbols(x::Symbol; kwargs...) = symbols(string(x); kwargs...)
symbols(xs::T...; kwargs...) where {T <: SymbolicObject} = xs

"""
    @vars x y z

Define symbolic values, possibly with names and assumptions

Examples:
```
@vars x y
@vars a1=>"α₁"
@vars a b real=true
```
!!! Note:
    The `@syms` macro is recommended as it has a more flexible syntax
"""
macro vars(x...)
    q = Expr(:block)
    as = []    # running list of assumptions to be applied
    ss = []    # running list of symbols created
    for s in reverse(x)
        if isa(s, Expr)    # either an assumption or a named variable
            if s.head == :(=)
                s.head = :kw
                push!(as, s)
            elseif s.head == :call && s.args[1] == :(=>)
                push!(ss, s.args[2])
                push!(q.args, Expr(:(=), esc(s.args[2]), Expr(:call, :symbols, s.args[3], map(esc,as)...)))
            end
        elseif isa(s, Symbol)   # raw symbol to be created
            push!(ss, s)
            push!(q.args, Expr(:(=), esc(s), Expr(:call, :symbols, Expr(:quote, s), map(esc,as)...)))
        else
            throw(AssertionError("@vars expected a list of symbols and assumptions"))
        end
    end
    push!(q.args, Expr(:tuple, map(esc,reverse(ss))...)) # return all of the symbols we created
    q
end

"""
    @syms a n::integer x::(real,positive)=>"x₀" y[-1:1] u() v()::real w()::(real,positive) y()[1:3]::real

Construct symbolic variables or functions along with specified assumptions. Similar to `@vars`, `sympy.symbols`, and `sympy.Function`, but the specification of the assumptions is more immediate than those interfaces which follow sympy's constructors.

Allows the specification of assumptions on the variables and functions.

* a type-like annontation, such as `n::integer` is equivalent to `sympy.symbols("n", integer=true)`. Multiple assumptions are combined using parentheses (e.g., `n::(integer,nonnegative)`.

The possible [values](https://docs.sympy.org/latest/modules/core.html#module-sympy.core.assumptions) for assumptions are: "commutative", "complex", "imaginary", "real", "integer", "odd", "even", "prime", "composite", "zero", "nonzero", "rational", "algebraic", "transcendental", "irrational", "finite", "infinite", "negative", "nonnegative", "positive", "nonpositive", "hermitian", "antihermetian".

* a tensor declaration form is provided to define arrays of variables, e.g. `x[-1:1]` or `y[1:4, 2:5]`.

* a symbolic function can be specified using a pair of parentheses after the name, as in `u()`.

* The return type of a function can have assumptions specified, as with a variable. E.g., `h()::complex`. How the symbolic function prints can be set as with a variable, e.g. `h()::complex=>"h̄"`.

* multiple definitions can be separated by commas

* How the symbol prints (the `__str__()` value) can be specified using the syntax `=>"name"`, as in `x=>"xₒ"`

## Examples:

```jldoctest constructors
julia> using SymPy

julia> @syms a b::nonnegative
(a, b)

julia> sqrt(a^2), sqrt(b^2)
(sqrt(a^2), b)
```

```jldoctest constructors
julia> @syms x::prime
(x,)

julia> ask(𝑄.negative(x)), ask(𝑄.integer(x)), ask(𝑄.even(x))  # (false, true, nothing)
(false, true, nothing)
```

```jldoctest constructors
julia> @syms a[0:5], x
(Sym[a₀, a₁, a₂, a₃, a₄, a₅], x)

julia> sum( aᵢ*x^(i) for (i,aᵢ) ∈ zip(0:5, a)) |> print
a₀ + a₁*x + a₂*x^2 + a₃*x^3 + a₄*x^4 + a₅*x^5
```

```jldoctest constructors
julia> @syms x u() v()::nonnegative
(x, u, v)

julia> sqrt(u(x)^2), sqrt(v(x)^2) # sqrt(u(x)^2), Abs(v(x))
(sqrt(u(x)^2), Abs(v(x)))
```

!!! Note:
    Many thanks to `@matthieubulte` for this contribution.
"""
macro syms(xs...)
    # If the user separates declaration with commas, the top-level expression is a tuple
    if length(xs) == 1 && isa(xs[1], Expr) && xs[1].head == :tuple
        _gensyms(xs[1].args...)
    elseif length(xs) > 0
        _gensyms(xs...)
    end
end

function _gensyms(xs...)
    asstokw(a) = Expr(:kw, esc(a), true)

    # Each declaration is parsed and generates a declaration using `symbols`
    symdefs = map(xs) do expr
        decl = parsedecl(expr)
        symname = sym(decl)
        symname, gendecl(decl)
    end
    syms, defs = collect(zip(symdefs...))

    # The macro returns a tuple of Symbols that were declared
    Expr(:block, defs..., :(tuple($(map(esc,syms)...))))
end

## avoid PyObject conversion as possible
Sym(x::T) where {T <: Number} = sympify(x)
Sym(x::Rational{T}) where {T} = Sym(numerator(x))/Sym(denominator(x))
function Sym(x::Complex{Bool})
    !x.re && x.im && return IM
    !x.re && !x.im && return zero(Sym)
    x.re && !x.im && return Sym(1)
    x.re && x.im && return Sym(1) + IM
end
Sym(x::Complex{T}) where {T} = Sym(real(x)) + Sym(imag(x)) * IM
Sym(xs::Symbol...) = Tuple(Sym.((string(x) for x in xs)))
Sym(x::AbstractString) = sympy.symbols(x)
Sym(s::SymbolicObject) = s
Sym(x::Irrational{T}) where {T} = convert(Sym, x)

convert(::Type{Sym}, s::AbstractString) = Sym(s)

sympify(s, args...; kwargs...) = pycall(sympy.sympify::PyCall.PyObject, Sym, s) #sympy.sympify(s, args...; kwargs...)

SymMatrix(s::SymMatrix) = s
SymMatrix(s::Sym) = sympy.ImmutableMatrix([s])
SymMatrix(A::Matrix) = sympy.ImmutableMatrix([A[i,:] for i in 1:size(A)[1]])


text:
import Base: +, -, *, /, %, abs
import Base: &, |, !, Bool, ==, !=
import Base: <, <=, >, >=
import Base: exp, log, sqrt, ^
import Base: sin, cos, tan, asin, acos, atan
import ResumableFunctions: @yield, @resumable

@resumable function sin(z, x_::Out)
  x = asin(z)
  # TODO: this produces duplicates if x == pi/2
  while true
    @yield x
    @yield (pi - x)
    x += 2*pi
  end
end

@resumable function cos(z, x_::Out)
  x = acos(z)
  while true
    @yield x
    @yield (2*pi - x)
    x += 2*pi
  end
end

@resumable function tan(z, x_::Out)
  x = atan(z)
  while true
    @yield x
    x += pi
  end
end

@resumable function ^(z, x_::Out, y)
  @yield exp(log(z) / y)
end

@resumable function ^(z, x, y_::Out)
  @yield log(z) / log(x)
end

@resumable function sqrt(z, x_::Out)
  @yield (z*z)
end

@resumable function exp(z, x_::Out)
  @yield log(z)
end

@resumable function log(z, x_::Out)
  @yield exp(z)
end

@resumable function abs(z, x_::Out)
  @yield z
  @yield (-z)
end

@resumable function +(z, x_::Out)
  # +x == z ==> x = z
  @yield z
end

@resumable function +(z, x, y_::Out)
  # x + y == z  ==>  y = z - x
  @yield (z - x)
end

@resumable function +(z, x_::Out, y)
  # x + y == z  ==>  x = z - y
  @yield (z - y)
end

@resumable function +(z::T, x_::Out, y_::Out) where {T <: Integer}
  x = typemin(T)
  while true
    for y in +(z, x, Out())
      @yield (x, y)
    end
    if x == typemax(T)
      break
    end
    x += 1
  end
end

@resumable function -(z, x_::Out)
  # unary -
  @yield (-z)
end

@resumable function -(z, x, y_::Out)
  # x - y == z  ==>  y = x - z
  @yield (x - z)
end

@resumable function -(z, x_::Out, y)
  # x - y == z  ==>  x = z + y
  @yield (z + y)
end

@resumable function -(z::T, x_::Out, y_::Out) where {T <: Integer}
  x = typemin(T)
  while true
    for y in -(z, x, Out())
      @yield (x, y)
    end
    if x == typemax(T)
      break
    end
    x += 1
  end
end

# For multiplication, distinguish between integers and floats.
# For integers, we only produce a result if the result is an integer.
@resumable function *(z::T, x::T, y_::Out) where {T <: Integer}
  # x * y == z  ==>  y = z / x
  if mod(z, x) == 0
    @yield div(z, x)
  end
end

@resumable function *(z, x, y_::Out)
  # x * y == z  ==>  y = z / x
  @yield z / x
end

@resumable function *(z::T, x_::Out, y::T) where {T <: Integer}
  # x * y == z  ==>  x = z / y
  if mod(z, y) == 0
    @yield div(z, y)
  end
end

@resumable function *(z, x_::Out, y)
  # x * y == z  ==>  x = z / y
  @yield z / y
end

@resumable function *(z::T, x_::Out, y_::Out) where {T <: Integer}
  x = typemin(T)
  while true
    for y in *(z, x, Out())
      @yield (x, y)
    end
    if x == typemax(T)
      break
    end
    x += 1
  end
end

@resumable function *(z, x_::Out, y_::Out)
  error("* cannot be materialized")
end

@resumable function /(z, x, y_::Out)
  # x / y == z  ==> y = x / z
  @yield x / z
end

@resumable function /(z, x_::Out, y)
  # x / y == z  ==> x = y * z
  @yield y * z
end

@resumable function /(z, x_::Out, y_::Out)
  error("/ cannot be materialized")
end

@resumable function %(z::T, x::T, y_::Out) where {T <: Integer}
  # x % y == z ==> x == ay + z for all a
  a = typemin(T)
  while true
    if a != 0
      @yield div(x - z, a)
    end
    if a == typemax(T)
      break
    end
    a += 1
  end
end

@resumable function %(z::T, x_::Out, y::T) where {T <: Integer}
  # x % y == z ==> x == ay + z for all a
  a = typemin(T)
  while true
    @yield (a*y) + z
    if a == typemax(T)
      break
    end
    a += 1
  end
end

@resumable function %(z, x_::Out, y_::Out)
  error("% cannot be materialized")
end

@resumable function <(z::Bool, x_::Out, y::T) where {T <: Integer}
  for x in >(z, y, Out())
    @yield x
  end
end

@resumable function <=(z::Bool, x_::Out, y::T) where {T <: Integer}
  for x in >=(z, y, Out())
    @yield x
  end
end

@resumable function >(z::Bool, x_::Out, y::T) where {T <: Integer}
  for x in <(z, y, Out())
    @yield x
  end
end

@resumable function >=(z::Bool, x_::Out, y::T) where {T <: Integer}
  for x in <=(z, y, Out())
    @yield x
  end
end

@resumable function <(z::Bool, x::T, y_::Out) where {T <: Integer}
  if z
    y = x+1
    while x < y
      @yield y
      y += 1
    end
  else
    y = x
    while !(x < y)
      @yield y
      y -= 1
    end
  end
end

@resumable function <=(z::Bool, x::T, y_::Out) where {T <: Integer}
  if z
    y = x
    while x <= y
      @yield y
      y += 1
    end
  else
    y = x-1
    while !(x <= y)
      @yield y
      y -= 1
    end
  end
end

@resumable function >(z::Bool, x::T, y_::Out) where {T <: Integer}
  if z
    y = x-1
    while x > y
      @yield y
      y -= 1
    end
  else
    y = x
    while !(x > y)
      @yield y
      y += 1
    end
  end
end

@resumable function >=(z::Bool, x::T, y_::Out) where {T <: Integer}
  if z
    y = x
    while x >= y
      @yield y
      y -= 1
    end
  else
    y = x+1
    while !(x >= y)
      @yield y
      y += 1
    end
  end
end

@resumable function ==(z::Bool, x, y_::Out)
  if z
    @yield x
  else
    error("== cannot be materialized")
  end
end

@resumable function ==(z::Bool, x_::Out, y)
  if z
    @yield y
  else
    error("== cannot be materialized")
  end
end

@resumable function ==(z::Bool, x_::Out, y_::Out)
  error("== cannot be materialized")
end

@resumable function !=(z::Bool, x, y_::Out)
  if !z
    @yield x
  else
    error("!= cannot be materialized")
  end
end

@resumable function !=(z::Bool, x_::Out, y)
  if !z
    @yield y
  else
    error("!= cannot be materialized")
  end
end

@resumable function !=(z::Bool, x_::Out, y_::Out)
  error("!= cannot be materialized")
end

@resumable function Bool(z::Bool, x_::Out)
  @yield z
end

@resumable function !(z::Bool, x_::Out)
  @yield (!z)
end

@resumable function and(z::Bool, x::Bool, y_::Out)
  if z
    if x
      @yield true
    end
  else
    if x
      @yield false
    else
      @yield false
      @yield true
    end
  end
end

@resumable function and(z::Bool, x_::Out, y::Bool)
  for x in and(z, y, Out())
    @yield x
  end
end

@resumable function and(z::Bool, x_::Out, y_::Out)
  if z
    @yield (true, true)
  else
    @yield (false, false)
    @yield (false, true)
    @yield (true, false)
  end
end

@resumable function |(z::Bool, x::Bool, y_::Out)
  if z
    if x
      @yield false
      @yield true
    else
      @yield true
    end
  else
    if !x
      @yield false
    end
  end
end

@resumable function |(z::Bool, x_::Out, y::Bool)
  for x in |(z, y, Out())
    @yield x
  end
end

@resumable function |(z::Bool, x_::Out, y_::Out)
  if !z
    @yield (false, false)
  else
    @yield (false, true)
    @yield (true, false)
    @yield (true, true)
  end
end

ext: jl


text:
@testset "src/io/geojson.jl" begin
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module API

using UUIDs
using Printf
import Random
import Dates
import LibGit2

import ..depots, ..depots1, ..logdir, ..devdir
import ..Operations, ..Display, ..GitTools, ..Pkg, ..UPDATED_REGISTRY_THIS_SESSION
using ..Types, ..TOML
using Pkg.Types: VersionTypes

preview_info() = printstyled("───── Preview mode ─────\n"; color=Base.info_color(), bold=true)

include("generate.jl")

function check_package_name(x::AbstractString, mode=nothing)
    if !(occursin(Pkg.REPLMode.name_re, x))
        message = "$x is not a valid packagename."
        if mode !== nothing && any(occursin.(['\\','/'], x)) # maybe a url or a path
            message *= "\nThe argument appears to be a URL or path, perhaps you meant " *
                "`Pkg.$mode(PackageSpec(url=\"...\"))` or `Pkg.$mode(PackageSpec(path=\"...\"))`."
        end
        pkgerror(message)
    end
    return PackageSpec(x)
end

develop(pkg::Union{AbstractString, PackageSpec}; kwargs...) = develop([pkg]; kwargs...)
develop(pkgs::Vector{<:AbstractString}; kwargs...) =
    develop([check_package_name(pkg, :develop) for pkg in pkgs]; kwargs...)
develop(pkgs::Vector{PackageSpec}; kwargs...)      = develop(Context(), pkgs; kwargs...)
function develop(ctx::Context, pkgs::Vector{PackageSpec};
                 shared::Bool=true, strict::Bool=false, kwargs...)
    pkgs = deepcopy(pkgs) # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)

    for pkg in pkgs
        # if julia is passed as a package the solver gets tricked
        pkg.name != "julia" || pkgerror("Trying to develop julia as a package")
        pkg.repo.rev === nothing || pkgerror("git revision can not be given to `develop`")
        pkg.name !== nothing || pkg.uuid !== nothing || pkg.repo.url !== nothing ||
            pkgerror("A package must be specified by `name`, `uuid`, `url`, or `path`.")
        pkg.version == VersionSpec() ||
            pkgerror("Can not specify version when tracking a repo.")
    end

    ctx.preview && preview_info()

    new_git = handle_repos_develop!(ctx, pkgs, shared)

    any(pkg -> Types.collides_with_project(ctx.env, pkg), pkgs) &&
        pkgerror("Cannot `develop` package with the same name or uuid as the project")

    Operations.develop(ctx, pkgs, new_git; strict=strict)
    ctx.preview && preview_info()
    return
end

add(pkg::Union{AbstractString, PackageSpec}; kwargs...) = add([pkg]; kwargs...)
add(pkgs::Vector{<:AbstractString}; kwargs...) =
    add([check_package_name(pkg, :add) for pkg in pkgs]; kwargs...)
add(pkgs::Vector{PackageSpec}; kwargs...)      = add(Context(), pkgs; kwargs...)
function add(ctx::Context, pkgs::Vector{PackageSpec}; strict::Bool=false, kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)

    for pkg in pkgs
        # if julia is passed as a package the solver gets tricked; this catches the error early on
        pkg.name == "julia" && pkgerror("Trying to add julia as a package")
        pkg.name !== nothing || pkg.uuid !== nothing || pkg.repo.url !== nothing ||
            pkgerror("A package must be specified by `name`, `uuid`, `url`, or `path`.")
        if (pkg.repo.url !== nothing || pkg.repo.rev !== nothing)
            pkg.version == VersionSpec() ||
                pkgerror("Can not specify version when tracking a repo.")
        end
    end

    ctx.preview && preview_info()
    Types.update_registries(ctx)

    repo_pkgs = [pkg for pkg in pkgs if (pkg.repo.url !== nothing || pkg.repo.rev !== nothing)]
    new_git = handle_repos_add!(ctx, repo_pkgs)
    # repo + unpinned -> name, uuid, repo.rev, repo.url, tree_hash
    # repo + pinned -> name, uuid, tree_hash

    project_deps_resolve!(ctx.env, pkgs)
    registry_resolve!(ctx.env, pkgs)
    stdlib_resolve!(ctx, pkgs)
    ensure_resolved(ctx.env, pkgs, registry=true)

    any(pkg -> Types.collides_with_project(ctx.env, pkg), pkgs) &&
        pkgerror("Cannot add package with the same name or uuid as the project")

    Operations.add(ctx, pkgs, new_git; strict=strict)
    ctx.preview && preview_info()
    return
end

rm(pkg::Union{AbstractString, PackageSpec}; kwargs...) = rm([pkg]; kwargs...)
rm(pkgs::Vector{<:AbstractString}; kwargs...)          = rm([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
rm(pkgs::Vector{PackageSpec}; kwargs...)               = rm(Context(), pkgs; kwargs...)

function rm(ctx::Context, pkgs::Vector{PackageSpec}; mode=PKGMODE_PROJECT, kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    foreach(pkg -> pkg.mode = mode, pkgs)

    for pkg in pkgs
        pkg.name !== nothing || pkg.uuid !== nothing ||
            pkgerror("Must specify package by either `name` or `uuid`.")
        if !(pkg.version == VersionSpec() && pkg.pinned == false &&
             pkg.tree_hash === nothing && pkg.repo.url === nothing &&
             pkg.repo.rev === nothing && pkg.path === nothing)
            pkgerror("Package may only be specified by either `name` or `uuid`")
        end
    end

    Context!(ctx; kwargs...)
    ctx.preview && preview_info()

    project_deps_resolve!(ctx.env, pkgs)
    manifest_resolve!(ctx.env, pkgs)
    ensure_resolved(ctx.env, pkgs)

    Operations.rm(ctx, pkgs)
    ctx.preview && preview_info()
    return
end

up(ctx::Context; kwargs...)                            = up(ctx, PackageSpec[]; kwargs...)
up(; kwargs...)                                        = up(PackageSpec[]; kwargs...)
up(pkg::Union{AbstractString, PackageSpec}; kwargs...) = up([pkg]; kwargs...)
up(pkgs::Vector{<:AbstractString}; kwargs...)          = up([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
up(pkgs::Vector{PackageSpec}; kwargs...)               = up(Context(), pkgs; kwargs...)

function up(ctx::Context, pkgs::Vector{PackageSpec};
            level::UpgradeLevel=UPLEVEL_MAJOR, mode::PackageMode=PKGMODE_PROJECT,
            update_registry::Bool=true, kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    foreach(pkg -> pkg.mode = mode, pkgs)

    Context!(ctx; kwargs...)
    ctx.preview && preview_info()
    if update_registry
        Types.clone_default_registries()
        Types.update_registries(ctx; force=true)
    end
    if isempty(pkgs)
        if mode == PKGMODE_PROJECT
            for (name::String, uuid::UUID) in ctx.env.project.deps
                push!(pkgs, PackageSpec(name=name, uuid=uuid))
            end
        elseif mode == PKGMODE_MANIFEST
            for (uuid, entry) in ctx.env.manifest
                push!(pkgs, PackageSpec(name=entry.name, uuid=uuid))
            end
        end
    else
        project_deps_resolve!(ctx.env, pkgs)
        manifest_resolve!(ctx.env, pkgs)
        ensure_resolved(ctx.env, pkgs)
    end
    Operations.up(ctx, pkgs, level)
    ctx.preview && preview_info()
    return
end

resolve(ctx::Context=Context()) =
    up(ctx, level=UPLEVEL_FIXED, mode=PKGMODE_MANIFEST, update_registry=false)

pin(pkg::Union{AbstractString, PackageSpec}; kwargs...) = pin([pkg]; kwargs...)
pin(pkgs::Vector{<:AbstractString}; kwargs...)          = pin([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
pin(pkgs::Vector{PackageSpec}; kwargs...)               = pin(Context(), pkgs; kwargs...)

function pin(ctx::Context, pkgs::Vector{PackageSpec}; kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)
    ctx.preview && preview_info()

    for pkg in pkgs
        pkg.name !== nothing || pkg.uuid !== nothing ||
            pkgerror("Must specify package by either `name` or `uuid`.")
        pkg.repo.url === nothing || pkgerror("Can not specify `repo` url")
        pkg.repo.rev === nothing || pkgerror("Can not specify `repo` rev")
    end

    foreach(pkg -> pkg.mode = PKGMODE_PROJECT, pkgs)
    project_deps_resolve!(ctx.env, pkgs)
    ensure_resolved(ctx.env, pkgs)
    Operations.pin(ctx, pkgs)
    return
end

free(pkg::Union{AbstractString, PackageSpec}; kwargs...) = free([pkg]; kwargs...)
free(pkgs::Vector{<:AbstractString}; kwargs...)          = free([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
free(pkgs::Vector{PackageSpec}; kwargs...)       = free(Context(), pkgs; kwargs...)

function free(ctx::Context, pkgs::Vector{PackageSpec}; kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)
    ctx.preview && preview_info()

    for pkg in pkgs
        pkg.name !== nothing || pkg.uuid !== nothing ||
            pkgerror("Must specify package by either `name` or `uuid`.")
        if !(pkg.version == VersionSpec() && pkg.pinned == false &&
             pkg.tree_hash === nothing && pkg.repo.url === nothing &&
             pkg.repo.rev === nothing && pkg.path === nothing)
            pkgerror("Package may only be specified by either `name` or `uuid`")
        end
    end

    foreach(pkg -> pkg.mode = PKGMODE_MANIFEST, pkgs)
    manifest_resolve!(ctx.env, pkgs)
    ensure_resolved(ctx.env, pkgs)

    find_registered!(ctx.env, UUID[pkg.uuid for pkg in pkgs])
    Operations.free(ctx, pkgs)
    return
end

test(;kwargs...)                                         = test(PackageSpec[]; kwargs...)
test(pkg::Union{AbstractString, PackageSpec}; kwargs...) = test([pkg]; kwargs...)
test(pkgs::Vector{<:AbstractString}; kwargs...)          = test([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
test(pkgs::Vector{PackageSpec}; kwargs...)               = test(Context(), pkgs; kwargs...)
function test(ctx::Context, pkgs::Vector{PackageSpec};
              coverage=false, test_fn=nothing,
              julia_args::Cmd=``,
              test_args::Cmd=``, kwargs...)
    pkgs = deepcopy(pkgs) # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)
    ctx.preview && preview_info()
    if isempty(pkgs)
        ctx.env.pkg === nothing && pkgerror("trying to test unnamed project") #TODO Allow this?
        push!(pkgs, ctx.env.pkg)
    else
        project_resolve!(ctx.env, pkgs)
        project_deps_resolve!(ctx.env, pkgs)
        manifest_resolve!(ctx.env, pkgs)
        ensure_resolved(ctx.env, pkgs)
    end
    Operations.test(ctx, pkgs; coverage=coverage, test_fn=test_fn, julia_args=julia_args, test_args=test_args)
    return
end

installed() = __installed(PKGMODE_PROJECT)
function __installed(mode::PackageMode=PKGMODE_MANIFEST)
    diffs = Display.status(Context(), PackageSpec[], mode=mode, use_as_api=true)
    version_status = Dict{String, Union{VersionNumber,Nothing}}()
    diffs == nothing && return version_status
    for entry in diffs
        version_status[entry.name] = entry.new.ver
    end
    return version_status
end

function gc(ctx::Context=Context(); kwargs...)
    Context!(ctx; kwargs...)
    ctx.preview && preview_info()
    env = ctx.env

    # If the manifest was not used
    usage_file = joinpath(logdir(), "manifest_usage.toml")

    # Collect only the manifest that is least recently used
    manifest_date = Dict{String, Dates.DateTime}()
    for (manifest_file, infos) in TOML.parse(String(read(usage_file)))
        for info in infos
            date = info["time"]
            manifest_date[manifest_file] = haskey(manifest_date, date) ? max(manifest_date[date], date) : date
        end
    end

    # Find all reachable packages through manifests recently used
    new_usage = Dict{String, Any}()
    paths_to_keep = String[]
    printpkgstyle(ctx, :Active, "manifests:")
    for (manifestfile, date) in manifest_date
        !isfile(manifestfile) && continue
        println("        $(Types.pathrepr(manifestfile))")
        manifest = try
            read_manifest(manifestfile)
        catch e
            @warn "Reading manifest file at $manifestfile failed with error" exception = e
            nothing
        end
        manifest == nothing && continue
        new_usage[manifestfile] = [Dict("time" => date)]
        for (uuid, entry) in manifest
            if entry.tree_hash !== nothing
                push!(paths_to_keep,
                      Operations.find_installed(entry.name, uuid, entry.tree_hash))
            end
        end
    end

    # Collect the paths to delete (everything that is not reachable)
    paths_to_delete = String[]
    for depot in depots()
        packagedir = abspath(depot, "packages")
        if isdir(packagedir)
            for name in readdir(packagedir)
                if isdir(joinpath(packagedir, name))
                    for slug in readdir(joinpath(packagedir, name))
                        versiondir = joinpath(packagedir, name, slug)
                        if !(versiondir in paths_to_keep) && isdir(versiondir)
                            push!(paths_to_delete, versiondir)
                        end
                    end
                end
            end
        end
    end

    pretty_byte_str = (size) -> begin
        bytes, mb = Base.prettyprint_getunits(size, length(Base._mem_units), Int64(1024))
        return @sprintf("%.3f %s", bytes, Base._mem_units[mb])
    end

    # Delete paths for noreachable package versions and compute size saved
    function recursive_dir_size(path)
        size = 0
        for (root, dirs, files) in walkdir(path)
            for file in files
                size += lstat(joinpath(root, file)).size
            end
        end
        return size
    end

    sz = 0
    for path in paths_to_delete
        sz_pkg = recursive_dir_size(path)
        if !ctx.preview
            try
                Base.rm(path; recursive=true)
            catch
                @warn "Failed to delete $path"
            end
        end
        printpkgstyle(ctx, :Deleted, Types.pathrepr(path) * " (" * pretty_byte_str(sz_pkg) * ")")
        sz += sz_pkg
    end

    # Delete package paths that are now empty
    for depot in depots()
        packagedir = abspath(depot, "packages")
        if isdir(packagedir)
            for name in readdir(packagedir)
                name_path = joinpath(packagedir, name)
                if isdir(name_path)
                    if isempty(readdir(name_path))
                        !ctx.preview && Base.rm(name_path)
                    end
                end
            end
        end
    end

    # Write the new condensed usage file
    if !ctx.preview
        open(usage_file, "w") do io
            TOML.print(io, new_usage, sorted=true)
        end
    end
    ndel = length(paths_to_delete)
    byte_save_str = ndel == 0 ? "" : (" (" * pretty_byte_str(sz) * ")")
    printpkgstyle(ctx, :Deleted, "$(ndel) package installation$(ndel == 1 ? "" : "s")$byte_save_str")

    ctx.preview && preview_info()
    return
end

build(pkgs...; kwargs...) = build([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
build(pkg::Array{Union{}, 1}; kwargs...) = build(PackageSpec[]; kwargs...)
build(pkg::PackageSpec; kwargs...) = build([pkg]; kwargs...)
build(pkgs::Vector{PackageSpec}; kwargs...) = build(Context(), pkgs; kwargs...)
function build(ctx::Context, pkgs::Vector{PackageSpec}; verbose=false, kwargs...)
    pkgs = deepcopy(pkgs)  # deepcopy for avoid mutating PackageSpec members
    Context!(ctx; kwargs...)

    ctx.preview && preview_info()
    if isempty(pkgs)
        if ctx.env.pkg !== nothing
            push!(pkgs, ctx.env.pkg)
        else
            for (uuid, entry) in ctx.env.manifest
                push!(pkgs, PackageSpec(entry.name, uuid))
            end
        end
    end
    project_resolve!(ctx.env, pkgs)
    foreach(pkg -> pkg.mode = PKGMODE_MANIFEST, pkgs)
    manifest_resolve!(ctx.env, pkgs)
    ensure_resolved(ctx.env, pkgs)
    Operations.build(ctx, pkgs, verbose)
end

#####################################
# Backwards compatibility with Pkg2 #
#####################################
function clone(url::String, name::String = "")
    @warn "Pkg.clone is only kept for legacy CI script reasons, please use `add`" maxlog=1
    ctx = Context()
    if !isempty(name)
        ctx.old_pkg2_clone_name = name
    end
    develop(ctx, [Pkg.REPLMode.parse_package_identifier(url; add_or_develop=true)])
end

function dir(pkg::String, paths::AbstractString...)
    @warn "`Pkg.dir(pkgname, paths...)` is deprecated; instead, do `import $pkg; joinpath(dirname(pathof($pkg)), \"..\", paths...)`." maxlog=1
    pkgid = Base.identify_package(pkg)
    pkgid === nothing && return nothing
    path = Base.locate_package(pkgid)
    path === nothing && return nothing
    return abspath(path, "..", "..", paths...)
end

precompile() = precompile(Context())
function precompile(ctx::Context)
    printpkgstyle(ctx, :Precompiling, "project...")

    pkgids = [Base.PkgId(uuid, name) for (name, uuid) in ctx.env.project.deps if !(uuid in keys(ctx.stdlibs))]
    if ctx.env.pkg !== nothing && isfile( joinpath( dirname(ctx.env.project_file), "src", ctx.env.pkg.name * ".jl") )
        push!(pkgids, Base.PkgId(ctx.env.pkg.uuid, ctx.env.pkg.name))
    end

    # TODO: since we are a complete list, but not topologically sorted, handling of recursion will be completely at random
    for pkg in pkgids
        paths = Base.find_all_in_cache_path(pkg)
        sourcepath = Base.locate_package(pkg)
        sourcepath == nothing && continue
        # Heuristic for when precompilation is disabled
        occursin(r"\b__precompile__\(\s*false\s*\)", read(sourcepath, String)) && continue
        stale = true
        for path_to_try in paths::Vector{String}
            staledeps = Base.stale_cachefile(sourcepath, path_to_try)
            staledeps === true && continue
            # TODO: else, this returns a list of packages that may be loaded to make this valid (the topological list)
            stale = false
            break
        end
        if stale
            printpkgstyle(ctx, :Precompiling, pkg.name)
            try
                Base.compilecache(pkg, sourcepath)
            catch
            end
        end
    end
    nothing
end

instantiate(; kwargs...) = instantiate(Context(); kwargs...)
function instantiate(ctx::Context; manifest::Union{Bool, Nothing}=nothing,
                     update_registry::Bool=true, kwargs...)
    Context!(ctx; kwargs...)
    if (!isfile(ctx.env.manifest_file) && manifest == nothing) || manifest == false
        up(ctx; update_registry=update_registry)
        return
    end
    if !isfile(ctx.env.manifest_file) && manifest == true
        pkgerror("manifest at $(ctx.env.manifest_file) does not exist")
    end
    Operations.prune_manifest(ctx.env)
    for (name, uuid) in ctx.env.project.deps
        get(ctx.env.manifest, uuid, nothing) === nothing || continue
        pkgerror("`$name` is a direct dependency, but does not appear in the manifest.",
                 " If you intend `$name` to be a direct dependency, run `Pkg.resolve()` to populate the manifest.",
                 " Otherwise, remove `$name` with `Pkg.rm(\"$name\")`.",
                 " Finally, run `Pkg.instantiate()` again.")
    end
    Operations.is_instanitated(ctx) && return
    Types.update_registries(ctx)
    pkgs = PackageSpec[]
    Operations.load_all_deps!(ctx, pkgs)
    Operations.check_registered(ctx, pkgs)
    new_git = UUID[]
    for pkg in pkgs
        pkg.repo.url !== nothing || continue
        sourcepath = Operations.source_path(pkg)
        isdir(sourcepath) && continue
        # download repo at tree hash
        push!(new_git, pkg.uuid)
        clonepath = Types.clone_path!(pkg.repo.url)
        tmp_source = Types.repo_checkout(clonepath, string(pkg.tree_hash))
        mkpath(sourcepath)
        mv(tmp_source, sourcepath; force=true)
    end
    new_apply = Operations.download_source(ctx, pkgs)
    Operations.build_versions(ctx, union(new_apply, new_git))
end

@deprecate status(mode::PackageMode) status(mode=mode)

status(; kwargs...) = status(PackageSpec[]; kwargs...)
status(pkg::Union{AbstractString,PackageSpec}; kwargs...) = status([pkg]; kwargs...)
status(pkgs::Vector{<:AbstractString}; kwargs...) =
    status([check_package_name(pkg) for pkg in pkgs]; kwargs...)
status(pkgs::Vector{PackageSpec}; kwargs...) = status(Context(), pkgs; kwargs...)
function status(ctx::Context, pkgs::Vector{PackageSpec}; diff::Bool=false, mode=PKGMODE_PROJECT)
    project_resolve!(ctx.env, pkgs)
    project_deps_resolve!(ctx.env, pkgs)
    if mode === PKGMODE_MANIFEST
        foreach(pkg -> pkg.mode = PKGMODE_MANIFEST, pkgs)
    end
    manifest_resolve!(ctx.env, pkgs)
    ensure_resolved(ctx.env, pkgs)
    Pkg.Display.status(ctx, pkgs, diff=diff, mode=mode)
    return nothing
end

function activate()
    Base.ACTIVE_PROJECT[] = nothing
    p = Base.active_project()
    p === nothing || printpkgstyle(Context(), :Activating, "environment at $(pathrepr(p))")
    return nothing
end
function _activate_dep(dep_name::AbstractString)
    Base.active_project() === nothing && return
    env = nothing
    try
        env = EnvCache()
    catch err
        err isa PkgError || rethrow()
        return
    end
    uuid = get(env.project.deps, dep_name, nothing)
    if uuid !== nothing
        entry = manifest_info(env, uuid)
        if entry.path !== nothing
            return joinpath(dirname(env.project_file), entry.path)
        end
    end
end
function activate(path::AbstractString; shared::Bool=false)
    if !shared
        # `pkg> activate path`/`Pkg.activate(path)` does the following
        # 1. if path exists, activate that
        # 2. if path exists in deps, and the dep is deved, activate that path (`devpath` above)
        # 3. activate the non-existing directory (e.g. as in `pkg> activate .` for initing a new env)
        if Types.isdir_windows_workaround(path)
            fullpath = abspath(path)
        else
            fullpath = _activate_dep(path)
            if fullpath === nothing
                fullpath = abspath(path)
            end
        end
    else
        # initialize `fullpath` in case of empty `Pkg.depots()`
        fullpath = ""
        # loop over all depots to check if the shared environment already exists
        for depot in Pkg.depots()
            fullpath = joinpath(Pkg.envdir(depot), path)
            isdir(fullpath) && break
        end
        # this disallows names such as "Foo/bar", ".", "..", etc
        if basename(abspath(fullpath)) != path
            pkgerror("not a valid name for a shared environment: $(path)")
        end
        # unless the shared environment already exists, place it in the first depots
        if !isdir(fullpath)
            fullpath = joinpath(Pkg.envdir(Pkg.depots1()), path)
        end
    end
    Base.ACTIVE_PROJECT[] = Base.load_path_expand(fullpath)
    p = Base.active_project()
    if p !== nothing
        n = ispath(p) ? "" : "new "
        printpkgstyle(Context(), :Activating, "$(n)environment at $(pathrepr(p))")
    end
    return nothing
end

function setprotocol!(;
    domain::AbstractString="github.com",
    protocol::Union{Nothing, AbstractString}=nothing
)
    GitTools.setprotocol!(domain=domain, protocol=protocol)
    return nothing
end

@deprecate setprotocol!(proto::Union{Nothing, AbstractString}) setprotocol!(protocol = proto) false

# API constructor
function Package(;name::Union{Nothing,AbstractString} = nothing,
                 uuid::Union{Nothing,String,UUID} = nothing,
                 version::Union{VersionNumber, String, VersionSpec, Nothing} = nothing,
                 url = nothing, rev = nothing, path=nothing, mode::PackageMode = PKGMODE_PROJECT)
    path !== nothing && url !== nothing &&
        pkgerror("cannot specify both a path and url")
    url !== nothing && version !== nothing &&
        pkgerror("`version` can not be given with `url`, use `rev` instead")
    repo = Types.GitRepo(rev = rev, url = url !== nothing ? url : path)
    version = version === nothing ? VersionSpec() : VersionSpec(version)
    uuid isa String && (uuid = UUID(uuid))
    PackageSpec(;name=name, uuid=uuid, version=version, mode=mode, path=nothing,
                special_action=PKGSPEC_NOTHING, repo=repo, tree_hash=nothing)
end
Package(name::AbstractString) = PackageSpec(name)
Package(name::AbstractString, uuid::UUID) = PackageSpec(name, uuid)
Package(name::AbstractString, uuid::UUID, version::VersionTypes) = PackageSpec(name, uuid, version)

end # module


text:
using OPCSPs
using Base.Test

include("apply_measurement.jl")
include("simple_op.jl")
include("simple_opcsp.jl")


text:
struct IdentityBasis <: AbstractFuncApprox
    num_inputs
    state_ranges
    IdentityBasis(num_inputs) = new(num_inputs, nothing)
    IdentityBasis(state_ranges::Array{T,2}) where {T} = new(size(state_ranges)[1], state_ranges)
end

function get_num_outputs(f::IdentityBasis)
    return f.num_inputs
end

function get_output_range(f::IdentityBasis)
    if f.state_ranges == nothing
        ranges = zeros((f.num_inputs, 2))
        ranges[:, 1] .= -Inf
        ranges[:, 2] .= Inf
        return ranges
    else
        return f.state_ranges
    end
end

function call(f::IdentityBasis, x::T)::T where T
    return x
end

function call!(out, f::IdentityBasis, x)
    out .= x
end

function clone(f::IdentityBasis)::IdentityBasis
    return IdentityBasis(f.num_inputs)
end

export IdentityBasis

include("fourierbasis.jl")
include("randomrbf.jl")
include("concat.jl")
export FourierBasis, RandomRBFBasis, ConcatBasis

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
function simurecwf(
    model::Union{elmod2d,acmod2d},
    source::Array{MTsource},
    rec::receiver,
    RecDataPath::String,
    Slices::Union{UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64},
    WFDataPath::String,
    OptCpnt::Array{String},
    )

    nT = model.medium.nT

    #============  write both Recordings and Wavefields ===========#
    fid_wf = open(RecDataPath, "a+")
    fid_record = open(WFDataPath, "a+")
    # wIB(fid_record,rec,model.medium.nT,OptCpnt)
    # wIB(fid_wf,model,Slices,OptCpnt)
    for it = 1 : nT
        for isn in 1:length(sou)
            addmt!(model.wf, sou[isn], it)
        end
          run!(model)
          wDB(fid_record,fid_wf,model.wf,rec.BDnloc,
            model.nwf,OptCpnt,it,Slices)
    end
    close(fid_wf)
end


text:
using MatpowerCases, Base.Test

println("Running tests...")

include("misc.jl")
include("brch.jl")
include("bus.jl")
include("gen.jl")
include("gencost.jl")
include("case.jl")


text:
# some common functions
# Query operations are located in the derived.jl file.

abstract AbstractIndexableBitVector <: AbstractVector{Bool}
typealias AbstractBitVector Union{BitVector,AbstractIndexableBitVector}

size(b::AbstractIndexableBitVector) = (length(b),)


text:
using HorizonSideRobots

function find_passage(r::Robot)
    side = Ost
    while isborder(r, Nord) == true
        putmarker!(r)
        while ismarker(r) == true
            if isborder(r, side) == true
                side = HorizonSide(mod(Int(side) + 2, 4))
            end
            move!(r, side)
        end
        side = HorizonSide(mod(Int(side)+2, 4))
    end
end


text:
"""
    find_first()

# Examples

```jldoctest

r = 10:-1:1

find_first(<(5), r)

```
"""
@propagate_inbounds find_first(f, x) = find_first(f, x, order(x))

# ==, isequal
function find_first(f::Fix2{<:Union{typeof(==),typeof(isequal)}}, r, ::Ordering)
    if isempty(r)
        return nothing
    else
        idx = unsafe_findvalue(f.x, r)
        @boundscheck if (firstindex(r) > idx || idx > lastindex(r)) || @inbounds(!f(r[idx]))
            return nothing
        end
        return idx
    end
end

# <, isless
function find_first(f::Fix2{<:Union{typeof(<),typeof(isless)}}, r, ::ForwardOrdering)
    @boundscheck if first(r) < f.x
        return firstindex(r)
    end
    return nothing
end
find_first(f::Fix2{<:Union{typeof(<),typeof(isless)}}, r, ::UnorderedOrdering) = nothing
function find_first(f::Fix2{<:Union{typeof(<),typeof(isless)}}, r, ::ReverseOrdering)
    idx = unsafe_findvalue(f.x, r)
    @boundscheck if firstindex(r) > idx
        return 1
    end
    @boundscheck if lastindex(r) < idx
        return nothing
    end
    return f(@inbounds(r[idx])) ? idx : ((idx != lastindex(r)) ? idx + 1 : nothing)
end

# <=
function find_first(f::Fix2{typeof(<=)}, r, ::ForwardOrdering)
    @boundscheck if first(r) <= f.x
        return firstindex(r)
    end
    return nothing
end
find_first(f::Fix2{typeof(<=)}, r, ::UnorderedOrdering) = nothing
function find_first(f::Fix2{typeof(<=)}, r, ::ReverseOrdering)
    idx = unsafe_findvalue(f.x, r)
    @boundscheck if lastindex(r) < idx
        return nothing
    end
    @boundscheck if firstindex(r) >= idx
        return 1
    end
    return @inbounds(f(r[idx])) ? idx : return idx - 1
end

# >
function find_first(f::Fix2{typeof(>)}, r, ::ForwardOrdering)
    idx = unsafe_findvalue(f.x, r)
    @boundscheck if lastindex(r) < idx
        return nothing
    end
    @boundscheck if firstindex(r) > idx
        return 1
    end
    return f(@inbounds(r[idx])) ? idx : (idx != lastindex(r) ? idx + 1 : nothing)
end
function find_first(f::Fix2{typeof(>)}, r, ::ReverseOrdering)
    @boundscheck if first(r) > f.x
        return firstindex(r)
    end
    return nothing
end
find_first(f::Fix2{typeof(>)}, r, ::UnorderedOrdering) = nothing

# >=
function find_first(f::Fix2{typeof(>=)}, r, ::ForwardOrdering)
    idx = unsafe_findvalue(f.x, r)
    @boundscheck if lastindex(r) < idx
        return nothing
    end
    @boundscheck if firstindex(r) > idx
        return 1
    end
    return f(@inbounds(r[idx])) ? idx : idx + 1
end
find_first(f::Fix2{typeof(>=)}, r, ::ReverseOrdering) = first(r) >= f.x ? firstindex(r) : nothing
find_first(f::Fix2{typeof(>=)}, r, ::UnorderedOrdering) = nothing


text:
const CuDense{ElT,VecT}                 = Dense{ElT,VecT} where {VecT<:CuVector}
const CuDenseTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:CuDense}

Dense{T, SA}(x::Dense{T, SB}) where {T<:Number, SA<:CuArray, SB<:Array} = Dense{T, SA}(CuArray(x))
Dense{T, SA}(x::Dense{T, SB}) where {T<:Number, SA<:Array, SB<:CuArray} = Dense{T, SA}(collect(x.data))
Dense{T, S}(size::Integer) where {T, S<:CuArray{<:T}} = Dense{T, S}(CUDA.zeros(T, size))
function Dense{T, S}(x::T, size::Integer) where {T, S<:CuArray{<:T}}
    arr = CuArray{T}(undef, size)
    fill!(arr, x)
    Dense{T, S}(arr)
end
cpu(x::CuDense{T}) where {T<:Number} = Dense(collect(x.data))
Base.complex(::Type{Dense{ElT, VT}}) where {ElT, VT<:CuArray} = Dense{complex(ElT),CuVector{complex(ElT)}}

similartype(::Type{<:CuArray{<:Any,N}}, eltype::Type) where {N} = CuArray{eltype,N}
NDTensors.similar(::Type{<:CuArray{T}}, dims) where {T} = CuArray{T,length(dims)}(undef, dims)

CuArray(x::CuDense{ElT}) where {ElT} = CuVector{ElT}(data(x))
CuArray{ElT, N}(x::CuDenseTensor{ElT, N}) where {ElT, N} = CuArray{ElT, N}(reshape(data(store(x)), dims(inds(x))...))
CuArray(x::CuDenseTensor{ElT, N}) where {ElT, N} = CuArray{ElT, N}(x)

*(D::Dense{T, AT},x::S) where {T,AT<:CuArray,S<:Number} = Dense(x .* data(D))

Base.:(==)(::Type{<:CuDense{ElT1,CVec1}}, ::Type{<:CuDense{ElT2,CVec2}}) where {ElT1,ElT2,CVec1,CVec2} = (ElT1 == ElT2)
Base.getindex(D::CuDense{<:Number})          = collect(data(D))[]
Base.getindex(D::CuDenseTensor{<:Number, 0}) = store(D)[]
LinearAlgebra.norm(T::CuDenseTensor) = norm(data(store(T)))

# This is for type promotion for Scalar*Dense
function Base.promote_rule(::Type{<:Dense{ElT1,CuVector{ElT1}}},
                           ::Type{ElT2}) where {ElT1,
                                                ElT2<:Number}
  ElR  = promote_type(ElT1,ElT2)
  VecR = CuVector{ElR}
  return Dense{ElR,VecR}
end

function Base.permutedims(T::CuDenseTensor{<:Number,N},
                          perm::NTuple{N,Int}) where {N}
  Tp = NDTensors.similar(T,ITensors.NDTensors.permute(inds(T),perm))
  #Tp = permute(T,perm; always_copy=true)
  permute!(Tp, T)
  return Tp
end

function Base.permutedims!(R::CuDenseTensor{<:Number,N},
                           T::CuDenseTensor{<:Number,N},
                           perm::NTuple{N,Int}) where {N}
    return permutedims!!(R, T, perm)
end

function permutedims!!(B::Tensor{ElT,N,StoreT,IndsB},
                       A::Tensor{ElT,N,StoreT,IndsA},
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->permute!(r,t)) where {N,ElT,IndsB,IndsA,StoreT<:CuDense{ElT}}
  Ais = inds(A)
  Bis = ITensors.NDTensors.permute(inds(A), perm)
  B = f(B, A)
  return B
end

function Base.similar(::Type{<:CuDenseTensor{ElT}},
                      inds::Dims) where {ElT}
    storage_arr = CuVector{ElT}(undef,dim(inds))
    return Tensor(Dense(storage_arr),inds)
end

function Base.similar(::Type{<:CuDenseTensor{ElT}},
                      inds) where {ElT}
    storage_arr = CuVector{ElT}(undef,dim(inds))
    return Tensor(Dense(storage_arr),inds)
end

function outer!(R::CuDenseTensor,
                T1::CuDenseTensor,
                T2::CuDenseTensor)
  R_dat = vec(array(T1))*transpose(vec(array(T2)))
  copyto!(data(store(R)), vec(R_dat))
  inds_outer = unioninds(inds(T1),inds(T2))
  return R
end

function contract!!(R::CuDenseTensor{<:Number,NR},
                    labelsR::NTuple{NR},
                    T1::CuDenseTensor{<:Number,N1},
                    labelsT1::NTuple{N1},
                    T2::CuDenseTensor{<:Number,N2},
                    labelsT2::NTuple{N2}) where {NR,N1,N2}
  if N1==0
    # TODO: replace with an add! function?
    # What about doing `R .= T1[] .* PermutedDimsArray(T2,perm)`?
    perm = getperm(labelsR,labelsT2)
    newT2 = Tensor(Dense(data(store(T1)).*data(store(T2))), inds(T2))
    permute!(R,newT2)
  elseif N2==0
    perm = getperm(labelsR,labelsT1)
    newT1 = Tensor(Dense(data(store(T2)).*data(store(T1))), inds(T1))
    permute!(R,newT1)
  elseif N1+N2==NR
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R,T1,T2)
    inds_outer = unioninds(inds(T1),inds(T2))
    R = Tensor(store(R), inds_outer)
  else
    R = _contract!!(R,labelsR,T1,labelsT1,T2,labelsT2)
  end
  return R
end

function permutedims!!(B::CuDenseTensor{ElT,0},
                       A::CuDenseTensor{ElT,0},
                       perm::NTuple{0,Int},
                       f=(r,t)->permute!(r,t)) where {ElT<:Number}
    Cs = f(B, A)
    return Tensor(Dense(vec(Cs)), IndexSet{0}())
end

function permutedims!!(B::CuDenseTensor{ElT,N},
                       A::CuDenseTensor{ElT,0},
                       perm::NTuple{N,Int},
                       f=(r,t)->permute!(r,t)) where {N, ElT<:Number}
    Cis = ITensors.NDTensors.permute(inds(B), perm)
    Cs = f(B, A)
    return Tensor(Dense(vec(Cs)), Cis)
end

function _contract_scalar!(R::CuDenseTensor{ElR}, labelsR,
                           T1::Number, labelsT1,
                           T2::Number, labelsT2,
                           α = one(ElR), β = zero(ElR)) where {ElR}
  if iszero(β)
    copyto!(data(R), [α * T1 * T2])
  elseif iszero(α)
    copyto!(data(R), β.*data(R))
  else
    copyto!(data(R), [α * T1 * T2] .+ β.*data(R))
  end
  return R
end

function _contract_scalar!(R::CuDenseTensor{ElR,NR}, labelsR,
                           T₁::CuDenseTensor, labelsT₁,
                           T₂::CuDenseTensor, labelsT₂,
                           α = one(ElR), β=zero(ElR)) where {ElR,NR}
    if nnz(T₁) == nnz(T₂) == 1
        new_R = Tensor(Dense(data(store(T₁)).*data(store(T₂))), inds(R))
        copyto!(store(R), store(new_R))
    elseif nnz(T₁) == 1
        props = ContractionProperties(labelsT₁, labelsT₂, labelsR)
        compute_contraction_properties!(props, T₁, T₂, R)
        R = _contract!(R, T₁, T₂, props, α, β)
        #perm = getperm(labelsR,labelsT₂)
        #newT2 = Tensor(Dense(data(store(T₁)).*data(store(T₂))), inds(T₂))
        #permute!(R,newT2)
    elseif nnz(T₂) == 1
        props = ContractionProperties(labelsT₁, labelsT₂, labelsR)
        compute_contraction_properties!(props, T₁, T₂, R)
        R = _contract!(R, T₁, T₂, props, α, β)
        #perm = getperm(labelsR,labelsT₁)
        #newT1 = Tensor(Dense(data(store(T₁)).*data(store(T₂))), inds(T₁))
        #permute!(R,newT1)
    else
        error("In _contract_scalar!, one tensor must be a scalar")
    end
    return R
end

function _gemm_contract!(CT::DenseTensor{El,NC},
                         AT::DenseTensor{El,NA},
                         BT::DenseTensor{El,NB},
                         props::ContractionProperties,
                         α::Number=one(El),
                         β::Number=zero(El)) where {El,NC,NA,NB}
    # TODO: directly use Tensor instead of Array
    C = array(CT)
    A = array(AT)
    B = array(BT)

    tA = 'N'
    if props.permuteA
        pA = NTuple{NA,Int}(props.PA)
        Ap = permutedims(A, pA)
        AM = reshape(Ap, props.dmid, props.dleft)
        tA = 'T'
    else
        #A doesn't have to be permuted
        if Atrans(props)
            AM = reshape(A,props.dmid,props.dleft)
            tA = 'T'
        else
            AM = reshape(A,props.dleft,props.dmid)
        end
    end

    tB = 'N'
    if props.permuteB
        pB = NTuple{NB,Int}(props.PB)
        Bp = permutedims(B, pB)
        BM = reshape(Bp, props.dmid, props.dright)
    else
        if Btrans(props)
            BM = reshape(B,props.dright,props.dmid)
            tB = 'T'
        else
            BM = reshape(B,props.dmid,props.dright)
        end
    end

    #TODO: this logic may be wrong
    if props.permuteC
        #Need to copy here since we will be permuting
        #into C later
        CM = reshape(copy(C), props.dleft, props.dright)
    else
        if Ctrans(props)
            CM = reshape(C,props.dleft,props.dright)
            (AM,BM) = (BM,AM)
            if tA==tB
                tA = tB = (tA == 'T' ? 'N' : 'T')
            end
        else
            CM = reshape(C, props.dleft, props.dright)
        end
    end

    CM = CUBLAS.gemm!(tA,tB,El(α),AM,BM,El(β),CM)

    if props.permuteC
        pC = NTuple{NC,Int}(props.PC)
        Cr = reshape(CM,props.newCrange...)
        @strided C .= permutedims(Cr, pC)
    end
    return C
end

function _contract!(CT::CuDenseTensor{El1,NC},
                    AT::CuDenseTensor{El2,NA},
                    BT::CuDenseTensor{El3,NB},
                    props::ContractionProperties,
                    α::Number=one(Ela),β::Number=zero(El3)) where {El1,El2,El3,NC,NA,NB}
  if ndims(CT) > 12 || ndims(BT) > 12 || ndims(AT) > 12
    return _gemm_contract!(CT, AT, BT, props, α, β)
  end
  Ainds = inds(AT)
  Adims = dims(Ainds)
  Binds = inds(BT)
  Bdims = dims(Binds)
  Cinds = inds(CT)
  Cdims = dims(Cinds)
  Adata = reshape(data(store(AT)),Adims...)
  Bdata = reshape(data(store(BT)),Bdims...)
  Cdata = reshape(data(store(CT)),Cdims...)
  contracted = commoninds(Ainds, Binds)
  A_only = uniqueinds(Ainds, Binds)
  B_only = uniqueinds(Binds, Ainds)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(contracted)
      push!(ind_dict, i)
  end
  if length(A_only) > 0
      for (idx, i) in enumerate(A_only)
          push!(ind_dict, i)
      end
  end
  if length(B_only) > 0
      for (idx, i) in enumerate(B_only)
          push!(ind_dict, i)
      end
  end
  ctainds = zeros(Int, length(Ainds))
  ctbinds = zeros(Int, length(Binds))
  ctcinds = zeros(Int, length(Cinds))
  for (ii, ia) in enumerate(Ainds)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Binds)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end
  for (ii, ic) in enumerate(Cinds)
      ctcinds[ii] = findfirst(x->x==ic, ind_dict)
  end
  id_op    = CUDA.CUTENSOR.CUTENSOR.CUTENSOR_OP_IDENTITY
  dict_key = ""
  for cc in zip(ctcinds, Cdims)
      dict_key *= string(cc[1]) * "," * string(cc[2]) * ","
  end
  for aa in zip(ctainds, Adims)
      dict_key *= string(aa[1]) * "," * string(aa[2]) * ","
  end
  for bb in zip(ctbinds, Bdims)
      dict_key *= string(bb[1]) * "," * string(bb[2]) * ","
  end
  if haskey(ENV, "CUTENSOR_AUTOTUNE") && tryparse(Int, ENV["CUTENSOR_AUTOTUNE"]) == 1
      if haskey(ContractionPlans, dict_key)
          dict_val = ContractionPlans[dict_key]
          algo  = dict_val
          #plan  = dict_val[2]
          Cdata = CUDA.CUTENSOR.contraction!(α, Adata, Vector{Char}(ctainds), id_op, Bdata, Vector{Char}(ctbinds), id_op, β, Cdata, Vector{Char}(ctcinds), id_op, id_op; algo=algo)
      else
          # loop through all algos
          # pick the fastest one
          # store that plan!
          best_time = 1e6
          best_plan = nothing
          best_algo = nothing
          max_algos = Ref{Int32}(C_NULL)
          CUDA.CUTENSOR.cutensorContractionMaxAlgos(max_algos)
          # fix once the other options are documented
          #algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(max_algos[] - 1))
          algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(-1))
          for algo in reverse(algos)
              try
                  Cdata, this_time, bytes, gctime, memallocs = @timed CUDA.CUTENSOR.contraction!(α, Adata, Vector{Char}(ctainds), id_op, Bdata, Vector{Char}(ctbinds), id_op, β, Cdata, Vector{Char}(ctcinds), id_op, id_op; algo=CUDA.CUTENSOR.cutensorAlgo_t(algo))
                  if this_time < best_time
                      best_time = this_time
                      #best_plan = this_plan
                      best_algo = CUDA.CUTENSOR.cutensorAlgo_t(algo)
                  end
              catch err
                  @warn "Algorithm $algo not supported"
              end
          end
          ContractionPlans[dict_key] = best_algo
      end
  else
      Cdata = CUDA.CUTENSOR.contraction!(α, Adata, Vector{Char}(ctainds), id_op, Bdata, Vector{Char}(ctbinds), id_op, β, Cdata, Vector{Char}(ctcinds), id_op, id_op)
  end
  return parent(Cdata)
end

function Base.:+(B::CuDenseTensor, A::CuDenseTensor)
  opC  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opA  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(inds(A))
      push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  CUTENSOR.elementwiseBinary!(one(eltype(Adata)), reshapeAdata, ctainds, opA, one(eltype(Bdata)), reshapeBdata, ctbinds, opC, C, ctcinds, opAC)
  copyto!(data(store(B)), vec(C))
  return B
end

function Base.:+(B::CuDense, Bis::IndexSet, A::CuDense, Ais::IndexSet)
  opA  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opC  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
      push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  Cis = Bis
  C = CUTENSOR.elementwiseBinary!(1, reshapeAdata, ctainds, opA, 1, reshapeBdata, ctbinds, opC, C, ctcinds, opAC)
  copyto!(data(B), vec(C))
  return C
end

function Base.:-(B::CuDenseTensor, A::CuDenseTensor)
  opC  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opA  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(inds(A))
      push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis))
  CUTENSOR.elementwiseBinary!(one(eltype(Adata)), reshapeAdata, ctainds, opA, -one(eltype(Bdata)), reshapeBdata, ctbinds, opC, C, ctcinds, opAC)
  copyto!(data(store(B)), vec(C))
  return B
end

function Base.:-(A::CuDense, Ais::IndexSet, B::CuDense, Bis::IndexSet)
  opA  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opC  = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
      push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  Cis = Bis
  C = CUTENSOR.elementwiseBinary!(one(eltype(Adata)), reshapeAdata, ctainds, opA, -one(eltype(Bdata)), reshapeBdata, ctbinds, opC, C, ctcinds, opAC)
  copyto!(data(B), vec(C))
  return C
end

function Base.permute!(B::CuDenseTensor, A::CuDenseTensor)
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
      push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  if ndims(A) < 12 # use CUTENSOR
      ctainds = zeros(Int, length(Ais))
      ctbinds = zeros(Int, length(Bis))
      for (ii, ia) in enumerate(Ais)
          ctainds[ii] = findfirst(x->x==ia, ind_dict)
      end
      for (ii, ib) in enumerate(Bis)
          ctbinds[ii] = findfirst(x->x==ib, ind_dict)
      end
      CUDA.CUTENSOR.permutation!(one(eltype(Adata)), reshapeAdata, Vector{Char}(ctainds), reshapeBdata, Vector{Char}(ctbinds))
  else # use GPUArrays
      perm = Int[]
      for aix in Ais
        b_pos = findfirst(bix->bix==aix, Bis)
        push!(perm, b_pos)
      end
      @assert isperm(perm)
      permutedims!(reshapeBdata, reshapeAdata, invperm(perm))
  end
  return vec(reshapeBdata)
end

function Base.permute!(B::CuDense, Bis::IndexSet, A::CuDense, Ais::IndexSet)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
      push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata,dims(Bis)...)
  reshapeAdata = reshape(Adata,dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x->x==ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x->x==ib, ind_dict)
  end

  CUDA.CUTENSOR.permutation!(one(eltype(Adata)), reshapeAdata, Vector{Char}(ctainds), reshapeBdata, Vector{Char}(ctbinds))
  return vec(reshapeBdata)
end

Base.:/(A::CuDenseTensor, x::Number) = A * inv(x)


text:
using NetCDF
using DataFrames
using CSV
using IcanProj
using Statistics
using StatsBase
using PyPlot

function forcings_mean(path_simulations, target_variable, df_links)

  file = joinpath(path_simulations, "forcings_st/$(target_variable)_1km.nc")

  nspace = length(ncread(file, "dim_space"))
  ntime = length(ncread(file, "dim_time"))

  data, tmp = zeros(nspace), zeros(nspace)

  for i in 1:ntime

    ncread!(file, target_variable, tmp, start=[1, i], count=[-1, 1])

    data += tmp / ntime

  end

  df_tmp = DataFrame(ind_senorge = Int.(ncread(file, "id")), data = data)

  df_links = join(df_links, df_tmp, on = :ind_senorge, kind = :left)

  return df_links

end

function compute_spread(df_res)

  df_std = DataFrame(resolution = String[], y = Float64[], lowererror = Float64[], uppererror = Float64[])

  for resolution in [:ind_5km, :ind_10km, :ind_25km, :ind_50km]

    df_aggregated = aggregate(df_res, resolution, std)

    data_std = df_aggregated.data_std

    ikeep = .!isnan.(data_std)

    data_std = data_std[ikeep]

    y = percentile(data_std, 50)

    lowererror = percentile(data_std, 50) .- percentile(data_std, 10)

    uppererror = percentile(data_std, 90) .- percentile(data_std, 50)

    push!(df_std, (String(resolution), y, lowererror, uppererror))

  end

  return df_std

end

function compute_variability(path_simulations, target_variable)

  df_links = CSV.read(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

  df_res = forcings_mean(path_simulations, target_variable, df_links)

  df_std = compute_spread(df_res)

  return df_std

end

# Settings

path_simulations = "/data04/jmg/fsm_simulations/netcdf"

path_results = joinpath(dirname(pathof(IcanProj)), "..", "plots", "forcings")

# Plot results

fig, axes = plt[:subplots](nrows = 3, ncols = 2)

fig[:set_size_inches](8, 7)

# Air temperature

target_variable = "tair"

df_std = compute_variability(path_simulations, target_variable)

axes[1][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[1][:xaxis][:set_ticklabels]([])

axes[1][:set_ylabel]("Air temperature (C)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Relative humidity

target_variable = "rhum"

df_std = compute_variability(path_simulations, target_variable)

axes[4][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[4][:xaxis][:set_ticklabels]([])

axes[4][:yaxis][:tick_right]()

axes[4][:yaxis][:set_label_position]("right")

axes[4][:set_ylabel]("Relative humidity (%)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Wind speed

target_variable = "wind"

df_std = compute_variability(path_simulations, target_variable)

axes[2][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[2][:xaxis][:set_ticklabels]([])

axes[2][:set_ylabel]("Wind speed (m/s)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Snowfall

target_variable = "snowf"

df_std = compute_variability(path_simulations, target_variable)

axes[5][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[5][:xaxis][:set_ticklabels]([])

axes[5][:yaxis][:tick_right]()

axes[5][:yaxis][:set_label_position]("right")

axes[5][:set_ylabel]("Snowfall (mm/3h)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Incoming shortwave radiation

target_variable = "iswr"

df_std = compute_variability(path_simulations, target_variable)

axes[3][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[3][:set_xticks](collect(1:4))

axes[3][:set_xticklabels](["5km", "10km", "25km", "50km"])

axes[3][:set_ylabel]("Incoming shortwave\nradiation (W/m2)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Incoming longwave radiation

target_variable = "ilwr"

df_std = compute_variability(path_simulations, target_variable)

axes[6][:errorbar](collect(1:4), df_std.y, yerr=permutedims([df_std.lowererror df_std.uppererror]), fmt="o")

axes[6][:set_xticks](collect(1:4))

axes[6][:set_xticklabels](["5km", "10km", "25km", "50km"])

axes[6][:yaxis][:tick_right]()

axes[6][:yaxis][:set_label_position]("right")

axes[6][:set_ylabel]("Incoming longwave\nradiation (W/m2)")

df_std |> CSV.write(joinpath(path_results, "$(target_variable)_std.csv"))

# Save figure

fig[:savefig](joinpath(path_results, "forcings_variability.png"), dpi=600)


text:
@testset "LazyNode" begin
    ss = ["Hello world.", "Make peace.", "Make food.", "Eat penam."]

    @test LazyNode{:Sentence}(ss).data == ss
    @test LazyNode(:Sentence, ss).data == ss
    @test LazyNode{:Sentence}(ss).data == ss
    @test nobs(LazyNode{:Sentence}(ss)) == 4
    @test (LazyNode{:Sentence}(ss))[[1,3]].data == ss[[1,3]]
    @test (LazyNode{:Sentence}(ss))[2].data == [ss[2]]
    @test catobs((LazyNode{:Sentence}(ss)[i] for i in [1:2,3:4])...).data == ss
    @test catobs((LazyNode{:Sentence}(ss)[i] for i in [1,2,3,4])...).data == ss
    @test reduce(catobs, [LazyNode{:Sentence}(ss)[i] for i in [1:2,3:4]]).data == ss
    @test reduce(catobs, [LazyNode{:Sentence}(ss)[i] for i in [1,2,3,4]]).data == ss

    ds = LazyNode{:Sentence}(ss)
    m = Mill.reflectinmodel(ds)
    @test m(ds).data ≈ m.m(Mill.unpack2mill(ds)).data
end


text:
#=
This is the place, where I put functions, which are so annoying in OpenGL, that I felt the need to wrap them and make them more "Julian"
Its also to do some more complex error handling, not handled by the debug callback
=#

function glGetShaderiv(shaderID::GLuint, variable::GLenum)
    result = Ref{GLint}(-1)
    glGetShaderiv(shaderID, variable, result)
    result[]
end
function glShaderSource(shaderID::GLuint, shadercode::Vector{UInt8})
    shader_code_ptrs = Ptr{UInt8}[pointer(shadercode)]
    len              = Ref{GLint}(length(shadercode))
    glShaderSource(shaderID, 1, shader_code_ptrs, len)
end
glShaderSource(shaderID::GLuint, shadercode::String) = glShaderSource(shaderID, Vector{UInt8}(shadercode))
function glGetAttachedShaders(program::GLuint)
    shader_count   = glGetProgramiv(program, GL_ATTACHED_SHADERS)
    length_written = GLsizei[0]
    shaders        = zeros(GLuint, shader_count)

    glGetAttachedShaders(program, shader_count, length_written, shaders)
    shaders[1:first(length_written)]
end

get_attribute_location(program::GLuint, name) = get_uniform_location(program, ascii(name))
get_attribute_location(program::GLuint, name::Symbol) = get_attribute_location(program, string(name))
function get_attribute_location(program::GLuint, name::String)
    const location::GLint = glGetAttribLocation(program, name)
    if location == -1
        warn(
            "Named attribute (:$(name)) is not an active attribute in the specified program object or\n
            the name starts with the reserved prefix gl_\n"
        )
    elseif location == GL_INVALID_OPERATION
        error(
            "program is not a value generated by OpenGL or\n
            program is not a program object or\n
            program has not been successfully linked"
        )
    end
    location
end

get_uniform_location(program::GLuint, name::Symbol) = get_uniform_location(program, String(name))
function get_uniform_location(program::GLuint, name::String)
    location = glGetUniformLocation(program, name)
    if location == -1
        error(
            """Named uniform (:$(name)) is not an active attribute in the specified program object or
            the name starts with the reserved prefix gl_"""
        )
    elseif location == GL_INVALID_OPERATION
        error("""program is not a value generated by OpenGL or
            program is not a program object or
            program has not been successfully linked"""
        )
    end
    location
end

function glGetActiveUniform(programID::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    uniformSize    = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(programID, GL_ACTIVE_UNIFORM_MAX_LENGTH)
    name           = Vector{GLchar}(maxcharsize)

    glGetActiveUniform(programID, index, maxcharsize, actualLength, uniformSize, typ, name)

    actualLength[1] <= 0 &&  error("No active uniform at given index. Index: ", index)

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]", "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], uniformSize[1])
end
function glGetActiveAttrib(programID::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    attributeSize  = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(programID, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH)
    name           = Vector{GLchar}(maxcharsize)

    glGetActiveAttrib(programID, index, maxcharsize, actualLength, attributeSize, typ, name)

    actualLength[1] <= 0 && error("No active uniform at given index. Index: ", index)

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]", "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], attributeSize[1])
end
function glGetProgramiv(programID::GLuint, variable::GLenum)
    result = Ref{GLint}(-1)
    glGetProgramiv(programID, variable, result)
    result[]
end
function glGetIntegerv(variable::GLenum)
    result = Ref{GLint}(-1)
    glGetIntegerv(UInt32(variable), result)
    result[]
end

function glGenBuffers(n=1)
    result = GLuint[0]
    glGenBuffers(1, result)
    id = result[]
    if id <= 0
        error("glGenBuffers returned invalid id. OpenGL Context active?")
    end
    id
end
function glGenVertexArrays()
    result = GLuint[0]
    glGenVertexArrays(1, result)
    id = result[1]
    if id <=0
        error("glGenVertexArrays returned invalid id. OpenGL Context active?")
    end
    id
end
function glGenTextures()
    result = GLuint[0]
    glGenTextures(1, result)
    id = result[1]
    if id <= 0
        error("glGenTextures returned invalid id. OpenGL Context active?")
    end
    id
end
function glGenFramebuffers()
    result = GLuint[0]
    glGenFramebuffers(1, result)
    id = result[1]
    if id <= 0
        error("glGenFramebuffers returned invalid id. OpenGL Context active?")
    end
    id
end

function glDeleteTextures(id::GLuint)
  arr = [id]
  glDeleteTextures(1, arr)
end
function glDeleteVertexArrays(id::GLuint)
  arr = [id]
  glDeleteVertexArrays(1, arr)
end
function glDeleteBuffers(id::GLuint)
  arr = [id]
  glDeleteBuffers(1, arr)
end

function glGetTexLevelParameteriv(target::GLenum, level, name::GLenum)
  result = GLint[0]
  glGetTexLevelParameteriv(target, level, name, result)
  result[1]
end

glViewport(x::SimpleRectangle) = glViewport(x.x, x.y, x.w, x.h)
glScissor(x::SimpleRectangle) = glScissor(x.x, x.y, x.w, x.h)

function glGenRenderbuffers(format::GLenum, attachment::GLenum, dimensions)
    renderbuffer = GLuint[0]
    glGenRenderbuffers(1, renderbuffer)
    glBindRenderbuffer(GL_RENDERBUFFER, renderbuffer[1])
    glRenderbufferStorage(GL_RENDERBUFFER, format, dimensions...)
    glFramebufferRenderbuffer(GL_FRAMEBUFFER, attachment, GL_RENDERBUFFER, renderbuffer[1])
    renderbuffer[1]
end

function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, h::Integer, d::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    glTexImage3D(GL_PROXY_TEXTURE_3D, level, internalFormat, w, h, d, border, format, datatype, C_NULL)
    for l in  0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_WIDTH)
        if result == 0
            error("glTexImage 3D: width too large. Width: ", w)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l,GL_TEXTURE_HEIGHT)
        if result == 0
            error("glTexImage 3D: height too large. height: ", h)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_DEPTH)
        if result == 0
            error("glTexImage 3D: depth too large. Depth: ", d)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            error("glTexImage 3D: internal format not valid. format: ", GLENUM(internalFormat).name)
        end
    end
    glTexImage3D(ttype, level, internalFormat, w, h, d, border, format, datatype, data)
end
function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, h::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    maxsize = glGetIntegerv(GL_MAX_TEXTURE_SIZE)
    glTexImage2D(GL_PROXY_TEXTURE_2D, level, internalFormat, w, h, border, format, datatype, C_NULL)
    for l in 0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_WIDTH)
        if result == 0
            error("glTexImage 2D: width too large. Width: ", w)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_HEIGHT)
        if result == 0
            error("glTexImage 2D: height too large. height: ", h)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            error("glTexImage 2D: internal format not valid. format: ", GLENUM(internalFormat).name)
        end
    end
    glTexImage2D(ttype, level, internalFormat, w, h, border, format, datatype, data)
end
function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    glTexImage1D(GL_PROXY_TEXTURE_1D, level, internalFormat, w, border, format, datatype, C_NULL)
    for l in 0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_1D, l, GL_TEXTURE_WIDTH)
        if result == 0
            error("glTexImage 1D: width too large. Width: ", w)
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_1D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            error("glTexImage 1D: internal format not valid. format: ", GLENUM(internalFormat).name)
        end
    end
    glTexImage1D(ttype, level, internalFormat, w, border, format, datatype, data)
end


text:
using ModelingToolkit, Optimization, OptimizationOptimJL

@variables t x(t) v(t) u(t)
@parameters p[1:2]
D = Differential(t)

loss = (4 - x)^2 + 2v^2 + u^2
eqs = [D(x) ~ v - p[2] * x
       D(v) ~ p[1] * u^3 + v]

@named sys = ControlSystem(loss, eqs, t, [x, v], [u], p)
dt = 0.1
tspan = (0.0, 1.0)
sys = runge_kutta_discretize(sys, dt, tspan)

u0 = rand(length(states(sys))) # guess for the state values
prob = OptimizationProblem(sys, u0, [0.1, 0.1], grad = true)
sol = solve(prob, BFGS())

# issue #819
@testset "Combined system name collisions" begin
    eqs_short = [D(x) ~ -p[2] * x
                 D(v) ~ p[1] * u^3]
    sys1 = ControlSystem(loss, eqs_short, t, [x, v], [u], p, name = :sys1)
    sys2 = ControlSystem(loss, eqs_short, t, [x, v], [u], p, name = :sys1)
    @test_throws ArgumentError ControlSystem(loss, [sys2.v ~ sys1.v], t, [], [], [],
                                             systems = [sys1, sys2], name = :foo)
end


text:
if BACKEND == "GPU"
    using CUDArt
    using CUBLAS
end
function init(material=train_data, num_buf_factor=20, prob_raise=0.75)
    global sampleFactory
    global genBufFac
    global init_model

    info(@sprintf("Using %s-%d", BACKEND, DEVICE))
    idx2wf = jldopen(init_model) do f
         read(f, "idx2wf")
    end

    for i = 1:length(idx2wf)
        @inbounds idx2wf[i] = idx2wf[i]^prob_raise
    end
    # sumx = sum(idx2wf)
    sumx = sum(idx2wf[i] for i = length(idx2wf):-1:1)
    for i = 1:length(idx2wf)
        @inbounds idx2wf[i] /= sumx
    end
    # info(length(idx2wf))

    info("Using $material as training data")
    set!(sampleFactory, SampleFactory(material, idx2wf))
    set!(genBufFac, convert(UInt32,num_buf_factor))
    if BACKEND == "GPU"
        # CUDArt.device(DEVICE)  # already set in config_global.jl
        CUDArt.init(DEVICE)
    end
end

function close_device()
    if BACKEND == "GPU"
        CUBLAS.destroy()
        CUDArt.close(DEVICE)
        CUDArt.device_reset(DEVICE)
    end
end


text:
# coding: utf-8

# In[1]:

#this has nothing to do with recursion algorithm
#it happened to be in the recursion chapter in my book
#so i kept it under recursion
#its about dynamic programming
#its kinda tricky to understand
#if you are familiar with convex optimization or lagrangian
#its better to use em instead of this

#knapsack problem is to maximize the value
#while having a weight constraint
#each value has a different weight
#the knapsack has a maximum capacity which is the constraint

#to solve the problem,we have to use recursive thinking
#lets create a list from 1 to the maximum capacity
#for each capacity in the list
#we try to reach the optimal allocation of weight at the given capacity
#say we have c as the maximum capacity
#we remove the last item i
#we get weight capacity-weight[i]
#we wanna make sure our allocation at capacity-weight[i] is still the optimal
#by optimal,we mean for the same weight we can achieve the highest value
#if capacity-weight[i] is not the optimal
#we can find another combo with the same weight but higher value
#we add the item i back into the knapsack then
#the new total value we get would be larger than the previous so-called optimal
#it will contradict the definition of optimal
#hence,for each capacity,we keep removing items
#until we reach base case 0,and it always stays the optimal at the given capacity

#to get the optimal status
#we shall do a traversal on all items
#we create a matrix with (number of items) * (maximum capacity)
#for each capacity level,we try to add a new item
#if adding new item causes the overall weight larger than the current capacity level
#the knapsack reverts to the previous status without item i which is matrix[i-1][j]
#if adding new item doesnt cause the overall weight bigger than the current capacity level
#we try to see whether adding item i would be the new optimal case
#so we compare the previous status with the status after adding item i
#the status after adding item i shall be matrix[i-1][j-weight[i-1]]+value[i-1]
#we use j-weight[i-1] cuz adding item i would reduce the capacity we have
#we have to use the current constraint level j to minus item i weight

# In[2]:

function knapsack(value,weight,capacity)

    #in this section,we create a nested list with size of (number of items+1)*(capacity+1)
    matrix=[[0 for _ in 1:(capacity+1)] for _ in 1:(length(value)+1)]

    #now we begin our traversal on all elements in matrix
    #i starts from 2 cuz we would be using i-1 to imply item i
    for i in 2:(length(value)+1)

        for j in 2:(capacity+1)

            #this is the part to check if adding item i-1 would exceed the current capacity j
            #if it does,we go back to the previous status
            #if not,we shall find out whether adding item i-1 would be the new optimal
            if weight[i-1]>j

                matrix[i][j]=matrix[i-1][j]

            else

                #julia index starts from 1
                #which is a pain in the ass
                #when current capacity==the new item s weight
                #it creates an issue
                if j==weight[i-1]

                    ind=1

                else

                    ind=j-weight[i-1]

                end

                #we use max funcion to see if adding item i-1 would be the new optimal
                matrix[i][j]=max(matrix[i-1][j],
                matrix[i-1][ind]+value[i-1])

            end

        end

    end

    return matrix[length(value)+1][capacity+1]

end

# In[3]:

println(knapsack([0,50,60,60,120],[0,10,15,20,40],50))


text:
#These are the commands to define the Fem class, and assemble the Matrix in Julia
#full documentation of the API from FEniCS can be found in the link below
#http://fenics.readthedocs.io/projects/UFL/en/latest/api-doc/Expression.html
#Tests for these can be found in the test_jfem.jl file.

using FEniCS

@fenicsclass FunctionSpace

#Use Functionspace for scalar fields
FunctionSpace(mesh::Mesh, family::StringOrSymbol, degree::Int) = FunctionSpace(fenics.FunctionSpace(mesh.pyobject, family, degree))
#Use VectorFunctionSpace for vector fields
VectorFunctionSpace(mesh::Mesh, family::StringOrSymbol,degree::Int) = FunctionSpace(fenics.VectorFunctionSpace(mesh.pyobject, family, degree))

export FunctionSpace, VectorFunctionSpace

@fenicsclass Expression
Argument(V,number,part::StringOrSymbol = nothing) = Expression(fenics.Argument(V.pyobject, number, part=part))
TrialFunction(V::FunctionSpace) = Expression(fenics.TrialFunction(V.pyobject))

TestFunction(V::FunctionSpace) = Expression(fenics.TestFunction(V.pyobject))

function TrialFunctions(V::FunctionSpace)
     vec=fenics.TrialFunctions(V.pyobject)
     expr_vec = [Expression(elem) for elem in vec]
     return expr_vec
end

function TestFunctions(V::FunctionSpace)
     vec=fenics.TestFunctions(V.pyobject)
     expr_vec = [Expression(elem) for elem in vec]
     return expr_vec
end

#Below are attributes for the argument *class*

export Argument,TrialFunction,TrialFunctions,TestFunction,TestFunctions

@fenicsclass Constant
Constant(x::Union{Real,Tuple}) = Expression(fenics.Constant(x, name="Constant($x)"))
export Constant

@fenicsclass FeFunction
function FeFunction(V::FunctionSpace; name::String="")
    if name == ""
        return FeFunction(fenics.Function(V.pyobject))
    else
        return FeFunction(fenics.Function(V.pyobject, name=name))
    end
end
assign(solution1::FeFunction,solution2 )=fenicspycall(solution1,:assign,solution2.pyobject)

function assign(solution::FeFunction, data::AbstractArray)
    solution.pyobject.vector().set_local(data)
end

geometric_dimension(expr::Union{FeFunction,Expression}) = fenicspycall(expr, :geometric_dimension)
export geometric_dimension

function split(fun::FeFunction)
     vec = fenics.split(fun.pyobject)
     expr_vec = [FeFunction(spl) for spl in vec]
     return expr_vec
end

function py_split(fun::FeFunction)
     vec = fun.pyobject.split()
     expr_vec = [FeFunction(spl) for spl in vec]
     return expr_vec
end

export FeFunction,assign,split,py_split

function Expression(cppcode; kw...)
    Expression(fenics.Expression(cppcode; kw...))
end
Identity(dim::Int) = Expression(fenics.Identity(dim))
inner(u::Union{Expression,FeFunction}, v::Union{Expression,FeFunction}) = Expression(fenics.inner(u.pyobject, v.pyobject))
outer(u::Union{Expression,FeFunction}, v::Union{Expression,FeFunction}) = Expression(fenics.outer(u.pyobject, v.pyobject))
dot(u::Union{Expression,FeFunction}, v::Union{Expression,FeFunction}) = Expression(fenics.dot(u.pyobject, v.pyobject))
grad(u::Union{Expression,FeFunction}) = Expression(fenics.grad(u.pyobject))
∇(u::Union{Expression,FeFunction}) = Expression(fenics.grad(u.pyobject))
nabla_grad(u::Union{Expression,FeFunction}) = Expression(ufl.nabla_grad(u.pyobject))
nabla_div(u::Union{Expression,FeFunction}) = Expression(ufl.nabla_div(u.pyobject))
div(u::Union{Expression,FeFunction}) = Expression(fenics.div(u.pyobject))
cross(u::Union{Expression,FeFunction}, v::Union{Expression,FeFunction}) = Expression(fenics.cross(u.pyobject, v.pyobject))
tr(u::Union{Expression,FeFunction}) = Expression(fenics.tr(u.pyobject))
sqrt(u::Union{Expression,FeFunction}) = Expression(fenics.sqrt(u.pyobject))
sym(u::Union{Expression,FeFunction}) = Expression(fenics.sym(u.pyobject))
len(U::Union{Expression,FeFunction}) = length(U.pyobject)

interpolate(solution1::FeFunction,solution2::Expression) = FeFunction(fenicspycall(solution1,:interpolate,solution2.pyobject))

Expression(x::FEniCS.Expression) = convert(Expression,x)
export Expression,Identity,inner,grad, nabla_grad, nabla_div,div, outer,dot,cross, tr, sqrt, sym, len, interpolate
export ∇

#Below are attributes for the Expression and FeFunction types
#Computes values at vertex of a mesh for a given Expression / FeFunction, and returns an array
compute_vertex_values(expr::Expression, mesh::Mesh) = fenicspycall(expr, :compute_vertex_values, mesh.pyobject)
compute_vertex_values(expr::FeFunction, mesh::Mesh) = fenicspycall(expr, :compute_vertex_values, mesh.pyobject)

export compute_vertex_values

@fenicsclass Measure
directional_derivative(solution1::FeFunction,direction)=FeFunction(fenicspycall(solution1,:dx,direction))

# some of these constant are initialized in __init__
export dx,ds,dS,dP,directional_derivative

#https://github.com/FEniCS/Expression/blob/master/Expression/measure.py
@fenicsclass Form

*(expr::Union{Expression,FeFunction}, measure::Measure) = Expression(measure.pyobject.__rmul__(expr.pyobject) )

*(expr::Union{Expression,FeFunction}, expr2::Union{Expression,FeFunction}) = Expression(expr.pyobject.__mul__(expr2.pyobject) )
*(expr::Real, expr2::Union{Expression,FeFunction}) = Expression(expr2.pyobject.__mul__(expr) )
*(expr::Union{Expression,FeFunction}, expr2::Real) = Expression(expr.pyobject.__mul__(expr2) )

+(expr::Union{Expression,FeFunction}, expr2::Real) = Expression(expr.pyobject.__add__(expr2) )
+(expr::Real, expr2::Union{Expression,FeFunction}) = Expression(expr2.pyobject.__add__(expr) )
+(expr::Union{Expression,FeFunction}, expr2::Union{Expression,FeFunction}) = Expression(expr.pyobject.__add__(expr2.pyobject) )

-(expr::Union{Expression,FeFunction}, expr2::Real) = Expression(expr.pyobject.__sub__(expr2) )
-(expr::Real, expr2::Union{Expression,FeFunction}) = -1*(Expression(expr2.pyobject.__sub__(expr) ))
-(expr::Union{Expression,FeFunction}, expr2::Union{Expression,FeFunction}) = Expression(expr.pyobject.__sub__(expr2.pyobject) )

/(expr::Union{Expression,FeFunction}, expr2::Real) = Expression(expr.pyobject.__div__(expr2) )
/(expr::Union{Expression,FeFunction}, expr2::Union{Expression,FeFunction}) = Expression(expr.pyobject.__div__(expr2.pyobject) )

function /(expr::Real,expr2::Union{Expression,FeFunction})
    x = expr2*expr2
    y = x/expr
    z = expr2/y
    return Expression(z)
end

function Transpose(object::Expression)
    x = object.pyobject.T
    y = Expression(x)
    return y
end
export Transpose

"""
rhs(equation::Expression)
Given a combined bilinear and linear form,
    extract the right hand side (negated linear form part).
    Example::
           a = u*v*dx + f*v*dx
           L = rhs(a) -> -f*v*dx
"""
rhs(equation::Expression)=Expression(fenics.rhs(equation.pyobject))
"""
lhs(equation::Expression)
    Given a combined bilinear and linear form,
    extract the left hand side (bilinear form part).
    Example::
        a = u*v*dx + f*v*dx
        a = lhs(a) -> u*v*dx
"""
lhs(equation::Expression)=Expression(fenics.lhs(equation.pyobject))

export lhs, rhs
#this assembles the matrix from a fenics form
@fenicsclass Matrix

assemble(assembly_item::Union{Form,Expression};tensor=nothing, form_compiler_parameters=nothing, add_values=false, finalize_tensor=true, keep_diagonal=false, backend=nothing) = Matrix(fenics.assemble(assembly_item.pyobject,
tensor=tensor,form_compiler_parameters=form_compiler_parameters,add_values=add_values,finalize_tensor=finalize_tensor,keep_diagonal=keep_diagonal,backend=backend))
export assemble

#I have changed this to Function+Form

#https://fenicsproject.org/olddocs/dolfin/1.6.0/python/programmers-reference/cpp/fem/DirichletBC.html
@fenicsclass sub_domain
@fenicsclass BoundaryCondition
"""
DirichletBC works in a slightly altered way to the original FEniCS.
Instead of defining the function with the required return command,
we define the return command directly (
ie instead of function(x)
                 return x
we simple write "x"
"""
DirichletBC(V::FunctionSpace,g,sub_domain)=BoundaryCondition(fenics.DirichletBC(V.pyobject,g.pyobject,sub_domain))#look this up with example also removed type from g(Should be expression)
DirichletBC(V::FunctionSpace,g::Number,sub_domain)=BoundaryCondition(fenics.DirichletBC(V.pyobject,g,sub_domain))#look this up with example also removed type from g(Should be expression)
DirichletBC(V::FunctionSpace,g::Tuple,sub_domain)=BoundaryCondition(fenics.DirichletBC(V.pyobject,g,sub_domain))#look this up with example also removed type from g(Should be expression)

export DirichletBC
#https://fenicsproject.org/olddocs/dolfin/2016.2.0/python/programmers-reference/compilemodules/subdomains/CompiledSubDomain.html
CompiledSubDomain(cppcode::String)  = sub_domain(fenics.CompiledSubDomain(cppcode))
export CompiledSubDomain

apply(bcs::BoundaryCondition, matrix::Matrix) = fenicspycall(bcs, :apply, matrix.pyobject)
apply(bcs, matrix::Matrix) = fenicspycall(BoundaryCondition(bcs), :apply, matrix.pyobject)

export apply

function assemble_system_julia(a::Expression,L::Expression)
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject)
    A  = A_fenics[:array]()
    b = b_fenics[:array]()
    return A,b
end

function assemble_system_julia(a::Expression,L::Expression,bc)
    bcs_py = [bcs.pyobject for bcs in bc]
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject,bcs_py)
    A  = A_fenics[:array]()
    b = b_fenics[:array]()
    return A,b
end

function assemble_system_julia(a::Expression,L::Expression,bc::BoundaryCondition)
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject,bc.pyobject)
    A  = A_fenics[:array]()
    b = b_fenics[:array]()
    return A,b
end

##add assemble_system to FEniCS objects for solving

function assemble_system(a::Expression,L::Expression)
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject)
    A  = Matrix(A_fenics)
    b = Matrix(b_fenics)
    return A,b
end

function assemble_system(a::Expression,L::Expression,bc)
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject,bc)
    A  = Matrix(A_fenics)
    b = Matrix(b_fenics)
    return A,b
end

function assemble_system(a::Expression,L::Expression,bc::BoundaryCondition)
    A_fenics,b_fenics = fenics.assemble_system(a.pyobject,L.pyobject,bc)
    A  = Matrix(A_fenics)
    b = Matrix(b_fenics)
    return A,b
end

export assemble_system, assemble_system_julia

"""
 For a full list of supported arguments, and their usage
please refer to http://matplotlib.org/api/pyplot_api.html
not all kwargs have been imported. Should you require any that are not imported
open as issue, and I will attempt to add them.
Deprecate this in a future version
"""
Plot(in_plot::Union{Mesh,FunctionSpace,FeFunction};alpha=1,animated=false,antialiased=true,color="grey"
,dash_capstyle="butt",dash_joinstyle="miter",dashes="",drawstyle="default",fillstyle="full",label="s",linestyle="solid",linewidth=1
,marker="",markeredgecolor="grey",markeredgewidth="",markerfacecolor="grey"
,markerfacecoloralt="grey",markersize=1,markevery="none",visible=true,title="") =fenics.common.plotting.plot(in_plot.pyobject,
alpha=alpha,animated=animated,antialiased=antialiased,color=color,dash_capstyle=dash_capstyle,dash_joinstyle=dash_joinstyle
,dashes=dashes,drawstyle=drawstyle,fillstyle=fillstyle,label=label,linestyle=linestyle,linewidth=linewidth,marker=marker,markeredgecolor=markeredgecolor
,markeredgewidth=markeredgewidth,markerfacecolor=markerfacecolor,markerfacecoloralt=markerfacecoloralt,markersize=markersize,markevery=markevery
,visible=visible,title=title)#the first is the keyword argument, the second is the value
#export Plot

@fenicsclass FiniteElement
"""
This is the FiniteElement class. Attributes/function inputs can be found below,
or via the FEniCS documentation.
This is only a basic implementation. Many of the expected methods are still missing
Arguments
|          family (string/symbol)
|             The finite element family
|          cell
|             The geometric cell
|          degree (int)
|             The polynomial degree (optional)
|          form_degree (int)
|             The form degree (FEEC notation, used when field is
|             viewed as k-form)
|          quad_scheme
|             The quadrature scheme (optional)
|          variant
|             Hint for the local basis function variant (optional)
"""
FiniteElement(family::StringOrSymbol,cell=nothing,degree=nothing,form_degree=nothing,quad_scheme=nothing,variant=nothing) =
FiniteElement(fenics.FiniteElement(family=family, cell=cell, degree=degree, form_degree=form_degree, quad_scheme=quad_scheme,variant=variant))
export FiniteElement

"""
Methods for the FiniteElement class
mapping(self)
 |
 |  reconstruct(self, family=None, cell=None, degree=None)
 |      Construct a new FiniteElement object with some properties
 |      replaced with new values.
 |
 |  shortstr(self)
 |      Format as string for pretty printing.
 |
 |  sobolev_space(self)
 |      Return the underlying Sobolev space.
 |
 |  variant(self)
 |
"""

# some of these constant are initialized in __init__
export hexahedron, tetrahedron, quadrilateral, triangle

family(finiteelement::FiniteElement) = fenicspycall(finiteelement, :family)

cell(finiteelement::FiniteElement) = fenicspycall(finiteelement, :cell)

degree(finiteelement::FiniteElement) = fenicspycall(finiteelement, :degree)

#form_degree(finiteelement::FiniteElement) = fenicspycall(finiteelement, :form_degree)

#quad_scheme(finiteelement::FiniteElement) = fenicspycall(finiteelement, :quad_scheme)

variant(finiteelement::FiniteElement) = fenicspycall(finiteelement, :variant)

reconstruct(finiteelement::FiniteElement, ;family=nothing,cell=nothing,degree=nothing) =
FiniteElement(fenicspycall(finiteelement, :reconstruct, family,cell,degree))

sobolev_space(finiteelement::FiniteElement) = fenicspycall(finiteelement, :sobolev_space)
export family, cell, degree, variant, reconstruct, sobolev_space

FacetNormal(mesh::Mesh)  = Expression(fenics.FacetNormal(mesh.pyobject))

export FacetNormal

@fenicsclass MixedElement

function MixedElement(vec::Array{FEniCS.FiniteElementImpl,1})
    pyvec = [elem.pyobject for elem in vec]
    me = MixedElement(fenics.MixedElement(pyvec))
end

export MixedElement

FunctionSpace(mesh::Mesh, element::Union{FiniteElement,MixedElement}) = FunctionSpace(fenics.FunctionSpace(mesh.pyobject, element.pyobject))

export FunctionSpace


text:
using StochasticNeighborhoodEmbeddings

d = 16 # relatively high-dimensional
n = 64
X1 = .1randn(d, n) .+ ones(d)
X2 = .1randn(d, n) .- ones(d)
X = hcat(X1, X2)
u = 30 # perplexity: higher value means distribution P[:, i] has larger entropy -> more uneven
Y = tsne(X, u, max_iter = 64, early_iter = 0, verbose = true)

using Plots
@views scatter(Y[1,:], Y[2,:])
gui()


text:
# Two step entry decision

make_test_entry_two_step(J, nc) =
    EntryTwoStepSwitches{Float64}(nTypes = J, nColleges = nc);

function init_entry_decision(objId :: ObjectId,
    switches :: EntryTwoStepSwitches{F1},
    st :: SymbolTable) where F1

    pEntryPref = init_entry_prefscale(switches, st);

    collPrefScale = switches.collPrefScale;
    pCollPref = Param(:collPrefScale,
        LatexLH.description(st, :collPresScale), latex(st, :collPrefScale),
        collPrefScale, collPrefScale, 0.1, 3.0,
        switches.calEntryPrefScale);

    pvec = ParamVector(objId, [pEntryPref]);
    return EntryTwoStep(objId,  pvec, ModelParams.value(pEntryPref),
        collPrefScale, switches)
end

"""
    $(SIGNATURES)

Entry probability for a student who is admitted to colleges in `admitV`.
Value of work argument is not used in this case, but provided so that interface is consistent with one step entry decision.

OUT:
- Entry prob by [type, college]. Conditional on entry
- Expected value by type
"""
function entry_probs(e :: EntryTwoStep{F1}, vWork_jV :: Vector{F1},
    vCollege_jcM :: Matrix{F1}, admitV;
    prefShocks :: Bool = true) where F1 <: Real

    J, nc = size(vCollege_jcM);
    prob_jxM = zeros(F1, J, nc);

    if isempty(admitV)
        eVal_jV = fill(F1(-1e8), J);
    else
        if prefShocks
            probM, eVal_jV = EconLH.extreme_value_decision(
                vCollege_jcM[:, admitV], entry_pref_scale(e); demeaned = true);
            prob_jxM[:, admitV] .= probM;
        else
            eVal_jV, icV = max_choices(fill(F1(-1e8), J),
                vCollege_jcM, admitV);
            prob_jxM[:, icV] .= one(F1);
        end
    end
    return prob_jxM, eVal_jV
end

# The same for one individual
function entry_probs(e :: EntryTwoStep{F1},
    vWork :: F1, vCollege_cV :: Vector{F1}, admitV;
    profShocks :: Bool = true) where F1 <: Real

    nc = length(vCollege_cV);
    prob_cV = zeros(F1, nc);

    @assert !isempty(admitV)  "Two step entry cannot have empty admission set"
    if prefShocks
        # eVal is a one element vector and probV is a Matrix
        probV, eVal = EconLH.extreme_value_decision_one(
            vCollege_cV[admitV], entry_pref_scale(e); demeaned = true);
        prob_cV[admitV] .= probV;
    else
        eVal, ic = max_choice(F1(-1e8), vCollege_cV, admitV);
        (ic > 0)  &&  (prob_cV[ic] = one(F1));
    end

    @assert check_float(eVal)
    return prob_cV, eVal
end

function entry_decisions(
    entryS :: EntryTwoStep{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    vWork_jV :: AbstractVector{F1}, vCollege_jcM :: AbstractMatrix{F1},
    endowPctV :: AbstractVector{F1}, rank_jV)  where {I1, F1}

    error("Not implemented for two step entry yet")
end

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