text:
# Forward propagation
function forward(model::PyObject, src_coords::Array{Float32}, rcv_coords::Array{Float32},
                 wavelet::Array{Float32}, e::Array{Float32}, space_order::Integer=8, isic::Bool=false)
    # Number of time steps
    nt = size(wavelet, 1)
    r = size(e, 2)

    # Setting forward wavefield
    u = wu.wavefield(model, space_order)

    # Set up PDE expression and rearrange
    pde = ker.wave_kernel(model, u)
    eu, probe_eq = time_probe(e, u; isic=isic)

    # Setup source and receiver
    geom_expr, _, rcv = geom.src_rec(model, u, src_coords=src_coords, nt=nt,
                                     rec_coords=rcv_coords, wavelet=wavelet)

    # Create operator and run
    subs = model.spacing_map
    op = dv.Operator(vcat(pde, geom_expr, probe_eq), subs=subs, name="forwardp$(r)", opt=ut.opt_op(model))

    summary = op()

    # Output
    return rcv.data, eu, summary
end

function backprop(model::PyObject, y::Array{Float32}, rcv_coords::Array{Float32},
                  e::Array{Float32}, space_order::Integer=8)
    # Number of time steps
    nt = size(y, 1)
    r = size(e, 2)

    # Setting adjoint wavefield
    v = wu.wavefield(model, space_order, fw=false)

    # Set up PDE expression and rearrange
    pde = ker.wave_kernel(model, v, fw=false)
    ev, probe_eq = time_probe(e, v; fw=false)

    # Setup source and receiver
    geom_expr, _, _ = geom.src_rec(model, v, src_coords=rcv_coords, nt=nt,
                                   wavelet=y, fw=false)

    # Create operator and run
    subs = model.spacing_map
    op = dv.Operator(vcat(pde, geom_expr, probe_eq), subs=subs, name="adjointp$(r)", opt=ut.opt_op(model))

    # Run operator
    summary = op()

    return ev, summary
end

# Forward propagation
function born(model::PyObject, src_coords::Array{Float32}, rcv_coords::Array{Float32},
              wavelet::Array{Float32}, e::Array{Float32}, space_order::Integer=8, isic::Bool=false)
    # Number of time steps
    nt = size(wavelet, 1)
    r = size(e, 2)

    # Setting forward wavefield
    u = wu.wavefield(model, space_order)
    ul = wu.wavefield(model, space_order, name="l")

    # Set up PDE expression and rearrange
    pde = ker.wave_kernel(model, u)
    pdel = model.dm == 0 ? [] : ker.wave_kernel(model, ul, q=si.lin_src(model, u, isic=isic))
    eu, probe_eq = time_probe(e, u; isic=isic)

    # Setup source and receiver
    geom_expr, _, rcv = geom.src_rec(model, u, src_coords=src_coords, nt=nt,
                            rec_coords=rcv_coords, wavelet=wavelet)
    geom_exprl, _, rcvl = geom.src_rec(model, ul, rec_coords=rcv_coords, nt=nt)

    # Create operator and run
    subs = model.spacing_map
    op = dv.Operator(vcat(pde, pdel, probe_eq, geom_expr, geom_exprl), subs=subs,
                     name="bornp$(r)", opt=ut.opt_op(model))

    summary = op()

    # Output
    return rcv.data, rcvl.data, eu, summary
end

# Forward propagation
function born_with_back(model::PyObject, src_coords::Array{Float32}, rcv_coords::Array{Float32},
                        wavelet::Array{Float32}, e::Array{Float32}, y::Array{Float32, 2},
                        space_order::Integer=8, isic::Bool=false)
    # Number of time steps
    nt = size(wavelet, 1)
    r = size(e, 2)

    # Setting forward wavefield
    u = wu.wavefield(model, space_order)
    ul = wu.wavefield(model, space_order, name="l")
    v = wu.wavefield(model, space_order, name="b")

    # Set up PDE expression and rearrange
    pde = ker.wave_kernel(model, u)
    pdev = ker.wave_kernel(model, v)
    pdel = model.dm == 0 ? [] : ker.wave_kernel(model, ul, q=si.lin_src(model, u, isic=isic))
    eu, ev, probe_eq = time_probe_sim(e, u, v; isic=isic)

    # Setup source and receiver
    geom_expr, _, rcv = geom.src_rec(model, u, src_coords=src_coords, nt=nt,
                                     rec_coords=rcv_coords, wavelet=wavelet)
    geom_exprl, _, rcvl = geom.src_rec(model, ul, rec_coords=rcv_coords, nt=nt)
    geom_exprv, _, _ = geom.src_rec(model, v, src_coords=rcv_coords, nt=nt, wavelet=y)

    # Create operator and run
    subs = model.spacing_map
    op = dv.Operator(vcat(pde, pdel, pdev, probe_eq, geom_expr, geom_exprl, geom_exprv),
                     subs=subs, name="bornbackp$(r)", opt=ut.opt_op(model))

    summary = op()

    # Output
    return rcvl.data, eu, ev, summary
end

function time_probe(e::Array{Float32, 2}, wf::PyObject; fw=true, isic=false, pe=nothing)
    p_e = dv.DefaultDimension(name="p_e", default_value=size(e, 2))
    s = fw ? wf.grid.time_dim.spacing : 1
    # Probing vector
    nt = size(e, 1)
    q = dv.TimeFunction(name="Q", grid=wf.grid, dimensions=(wf.grid.time_dim, p_e),
                        shape=(nt, size(e, 2)), time_order=0, initializer=e)
    # Probed output
    pe = dv.Function(name="$(wf.name)e", grid=wf.grid, dimensions=(wf.grid.dimensions..., p_e),
                     shape=(wf.grid.shape..., size(e, 2)), space_order=wf.space_order)

    if size(e, 2) < 17
        probing = [dv.Inc(pe, s*q*ic(wf, fw, isic)).xreplace(Dict(p_e => i)) for i=1:size(e, 2)]
    else
        probing = [dv.Inc(pe, s*q*ic(wf, fw, isic))]
    end
    return pe, probing
end

function time_probe(e::Array{Float32, 2}, wf::Tuple{PyCall.PyObject, PyCall.PyObject}; fw=true, isic=false, pe=nothing)
    p_e = dv.DefaultDimension(name="p_e", default_value=size(e, 2))
    s = fw ? wf[1].grid.time_dim.spacing : 1
    # Probing vector
    nt = size(e, 1)
    q = dv.TimeFunction(name="Q", grid=wf[1].grid, dimensions=(wf[1].grid.time_dim, p_e),
                        shape=(nt, size(e, 2)), time_order=0, initializer=e)
    # Probed output
    pe = dv.Function(name="$(wf[1].name)e", grid=wf[1].grid, dimensions=(wf[1].grid.dimensions..., p_e),
                     shape=(wf[1].grid.shape..., size(e, 2)), space_order=wf[1].space_order)

    if size(e, 2) < 17
        probing = [dv.Inc(pe, s*q*ic(wf, fw, isic)).xreplace(Dict(p_e => i)) for i=1:size(e, 2)]
    else
        probing = [dv.Inc(pe, s*q*ic(wf, fw, isic))]
    end
    return pe, probing
end

function time_probe_sim(e::Array{Float32, 2}, wfu::PyObject, wfv::PyObject; isic=false)
    p_e = dv.DefaultDimension(name="p_e", default_value=size(e, 2))
    s = wfu.grid.time_dim.spacing
    # Probing vector
    nt = size(e, 1)
    q = dv.TimeFunction(name="Q", grid=wfu.grid, dimensions=(wfu.grid.time_dim, p_e),
                        shape=(nt, size(e, 2)), time_order=0, initializer=e)
    q_r = q._subs(wfu.grid.time_dim, wfu.grid.time_dim.symbolic_max - wfu.grid.time_dim)
    # Probed output
    peu = dv.Function(name="$(wfu.name)e", grid=wfu.grid, dimensions=(wfu.grid.dimensions..., p_e),
                      shape=(wfu.grid.shape..., size(e, 2)), space_order=wfu.space_order)
    pev = dv.Function(name="$(wfv.name)e", grid=wfv.grid, dimensions=(wfv.grid.dimensions..., p_e),
                      shape=(wfv.grid.shape..., size(e, 2)), space_order=wfv.space_order)

    if size(e, 2) < 17
        probing = [dv.Inc(peu, s*q*ic(wfu, true, isic)).xreplace(Dict(p_e => i)) for i=1:size(e, 2)]
        probing = vcat(probing, [dv.Inc(pev, s*q_r*ic(wfv, false, isic)).xreplace(Dict(p_e => i)) for i=1:size(e, 2)])
    else
        probing = [dv.Inc(peu, s*q*ic(wfu, true, isic)), dv.Inc(pev, s*q_r*ic(wfv, false, isic))]
    end
    return peu, pev, probing
end

function ic(wf::PyObject, fw::Bool, isic::Bool)
    return (isic || ~fw) ? wf : wf.dt2
end

function ic(wf::Tuple{PyCall.PyObject, PyCall.PyObject}, fw::Bool, isic::Bool)
    return (isic || ~fw) ? (wf[1] + wf[2]) : (wf[1].dt2 + wf[2].dt2)
end

function combine(eu::PyObject, ev::PyObject, isic::Bool=false)
    g = dv.Function(name="ge", grid=eu.grid, space_order=0)
    eq = isic ? (eu * ev.laplace + si.inner_grad(eu, ev)) : eu * ev
    eq = eq.subs(ev.indices[end], eu.indices[end])
    op = dv.Operator(dv.Inc(g, -eq), subs=eu.grid.spacing_map)
    op()
    g
end

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
include("neural2.jl")
include("imgtransform.jl")

using ArgParse
using Mmap
using Images
using ImageTransformations
using Serialization # TODO: change to better format

names = ["circle", "triangle", "rectangle"]

function parseArgs()::Dict{Symbol, Any}
  parser = ArgParseSettings()

  @add_arg_table! parser begin
    "--train", "-t"
      help = "Training mode. This option takes a path to the targets matrix and the number of data points as arguments."
      metavar = ["TARGETS", "LENGTH"]
      nargs = 2
    "--weights", "-w"
      help = "Set path to the saved weights"
      arg_type = String
      default = "weights"
    "--test", "-T"
      help = "Test the trained neural network"
      metavar = ["TARGETS", "LENGTH"]
      nargs = 2
    "data"
      help = "Eighter image to process or a training data"
      required = true
  end

  parse_args(parser; as_symbols = true)
end

 # TODO: change to better format
function saveNet(net::Net, path::String)
  serialize(path, net)
end

# TODO: change to better format
loadNet(path::String) = deserialize(path)

function loadMmaps(dataPath::String, targetsPath::String, length::Int)
  dataMap = mmap(dataPath, BitMatrix, (900, length))
  targetsMap = mmap(targetsPath, BitMatrix, (3, length))

  return dataMap, targetsMap
end

function trainingMode(net::Neural.Net, dataPath::String, targetsPath::String, length::Int)
  data, targets = loadMmaps(dataPath, targetsPath, length)
  train!(net, data, targets, 150, 100, 1e-2, 10/9)

  saveNet(net, "weights")
end

function testingMode(weightsPath::String, dataPath::String, targetsPath::String, length::Int)
  data, targets = loadMmaps(dataPath, targetsPath, length)
  net = loadNet!(weightsPath)

  println("Testing accuracy is $(Neural.accuracy(net, data, targets) * 100)%")
end

# Preparse the image for the neural network
function loadAndPrepareImg(imgPath::String)
  img = load(imgPath)

  binarized = ImgTransform.imgToBitMatrix(img) # convert image to "negative" (swapped 0 and 1) BitMatrix
  cutout = ImgTransform.cutoutShape(binarized) # locate the shape on image, cut it out
  positive = xor.(cutout, 1) # change back to positive
  resized = imresize(positive, 30, 30) # resize; returns Matrix{Float64}
  resized = BitMatrix(round.(Int, resized))

  return ImgTransform.toVector(resized) # return vector
end

function identify(weightsPath::String, imgPath::String)
  net = loadNet(weightsPath)

  img = loadAndPrepareImg(imgPath)
  probs = forward(net, img)[end]
  id = argmax(probs)

  println("The shape in the image is a $(names[id]) ($(round(probs[id], digits=4) * 100)% confidence).")

  for (l, p) in zip(names, probs)
    print("$l: $p ")
  end
  println()
end

parsed = parseArgs()

if length(parsed[:train]) > 0
  len = parse(Int, parsed[:train][2])
  trainingMode(net, parsed[:data], parsed[:train][1], len)
elseif length(parsed[:test]) > 0
  len = parse(Int, parsed[:test][2])
  testingMode(parsed[:weights], parsed[:data], parsed[:test][1], len)
else
  identify(parsed[:weights], parsed[:data])
end

ext: jl

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition


text:
using Displaz
using StaticArrays
using RoamesGeometry

println("Press ENTER to make a new example")
println("Type any non-empty string and press ENTER to quit")

while(readline(STDIN) == "")
	l1 = Line(rand(SVector{3,Float64}), rand(SVector{3,Float64}))
	l2 = Line(rand(SVector{3,Float64}), rand(SVector{3,Float64}))
	(p1, p2) = closest_points(l1, l2)
	l3 = Line(p1, p2)

	plot3d!(l1, color = [0.6, 0.1, 1.0], label = "Line 1")
	plot3d!(l2, color = [0.1, 0.6, 1.0], label = "Line 2")
	plot3d!(l3, color = [0.8, 0.8, 0.2], label = "Closest points")
end


text:
################################################################################
# parray.jl
# The Persistent Array type, composed using the PTree type.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2019 Noah C. Benson

# ==============================================================================
# #TODO List

# psparse() method for making PArrays similar to SparseArrays.sparse().
# pzeros() similar to zeros() and spzeros()
# prand() similar to rand() and sprand()
# prandn() similar to randn() and sprandn()
# pdiagm() similar to spdiagm()
# blockdiag() method instance.
# permute() method
# broadcast() method
# Many others: overload arithmetic operators? reshape, others.

# ==============================================================================
# PArray definition.

# We expand on some sparse array methods with PArrays (which are implemented as
# efficitn sparse arrays anyway.
import SparseArrays

"""
    PArray{T,N}

The PArray type is a persistent/immutable corrolary to the Array{T,N} type. Like
Array, PArray can store n-dimensional non-ragged arrays. However, unlike Arrays,
PArrays can create duplicates of themselves with finite edits in log-time.

PArrays have a similar interface as Arrays, but instead of the functions
`push!`, `pop!`, and `setindex!`, PArrays use `push`, `pop`, and `setindex`.
PArrays also have efficient implementations of `pushfirst` and `popfirst`.

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"0-element (PArray{Any, ?1}|PVector{Any})"
julia> PArray()
0-element PArray{Any,1}
```

```jldoctest; filter=r"[23]-element (PArray{Int64, ?1}|PVector{Int64}):"
julia> u = PArray{Int,1}([1,2])
2-element PArray{Int64,1}:
 1
 2

julia> push(u, 3)
3-element PArray{Int64,1}:
 1
 2
 3
```

```jldoctest; filter=r"2×3 (PArray{Symbol, ?2}|PMatrix{Symbol}):"
julia> PArray{Symbol,2}(:abc, (2,3))
2×3 PArray{Symbol,2}:
 :abc  :abc  :abc
 :abc  :abc  :abc
```
"""
struct PArray{T,N} <: AbstractPArray{T,N}
    # The initial element of the list; the tree is basically an enormous
    # circular buffer, so it's okay for this to roll around at UInt max.
    _i0::HASH_T
    # The dimensions of the array and the indices into the tree.
    _index::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
    # The data in a long array.
    _tree::PTree{T}
    # The default value, if any (for sparse arrays).
    _default::Union{Nothing, Tuple{T}}
end

# ==============================================================================
# PArray Constructors

PArray{T,N}(default, size::NTuple{N,<:Integer}) where {T,N} = begin
    default = Tuple{T}((default,))
    return PArray{T,N}(0x0, LinearIndices(size), PTree{T}(), default)
end
PArray{T,N}(::UndefInitializer, size::NTuple{N,<:Integer}) where {T,N} =
    return PArray{T,N}(0x0, LinearIndices(size), PTree{T}(), nothing)
PArray{T,N}(default::S, size::Vararg{<:Integer,N}) where {T,N,S} = begin
    return PArray{T,N}(default, NTuple{N,Int}(size))
end
PArray{T,N}(default::S, size::Vector{<:Integer}) where {T,N,S} = begin
    return PArray{T,N}(default, NTuple{N,Int}(size))
end
PArray(default::T, size::NTuple{N,<:Integer}) where {T,N} = PArray{T,N}(default, size)
PArray(default::T, size::Vararg{<:Integer,N}) where {T,N} = begin
    return PArray{T,N}(default, NTuple{N,Int}(size))
end
PArray(default::T, size::Vector{<:Integer}) where {T,N} = begin
    return PArray{T,N}(default, NTuple{N,Int}(size))
end
PArray{T,N}(a::AbstractArray{S,N}) where {T,N,S} = begin
    tree = PTree{T}(a)
    return PArray{T,N}(HASH_T(0x0), _lindex(size(a)...), tree, nothing)
end
PArray{T,N}(p::PArray{T,N}) where {T,N} = p
PArray{T,N}() where {T,N} = PArray{T,N}(undef, (0, [1 for _ in 2:N]...))
PArray(a::AbstractArray{T,N}) where {T,N} = PArray{T,N}(a)
PArray(p::PArray{T,N}) where {T,N} = p
PArray() = PArray{Any,1}()
# Convert function also.
Base.convert(::Type{PArray{T,N}}, x) where {T,N} = PArray{T,N}(x)
Base.convert(::Type{PArray{T,N}}, x::PArray{T,N}) where {T,N} = x

# ==============================================================================
# PArray aliases

"""
    PVector{T}

An alias for `PArray{T,1}`, representing a persistent vector.
"""
const PVector{T} = PArray{T,1} where {T}
PVector(default::T, len::II) where {T,II<:Integer} = PArray{T,1}(default, (len,))
PVector(default::T, len::Tuple{II}) where {T,II<:Integer} =
    PArray{T,1}(default, len)
PVector(a::AbstractArray{T,1}) where {T} = PArray{T,1}(a)
PVector(p::PArray{T,1}) where {T,N} = p
PVector() = PArray{Any,1}()
export PVector

"""
    PMatrix{T}

An alias for `PArray{T,2}`, representing a persistent matrix.
"""
const PMatrix{T} = PArray{T,2} where {T}
PMatrix(args...) = PArray{Any,2}(args...)
PMatrix(default::T, rs::II, cs::JJ) where {T,II<:Integer,JJ<:Integer} =
    PArray{T,2}(default, (rs,cs))
PMatrix(default::T, sz::Tuple{<:Integer,<:Integer}) where {T} =
    PArray{T,2}(default, len)
PMatrix(a::AbstractArray{T,2}) where {T} = PArray{T,2}(a)
PMatrix(p::PArray{T,2}) where {T,N} = p
PMatrix() = PArray{Any,2}()
export PMatrix

# ==============================================================================
# SparseArrays methods.

"""
    nnz(p::PArray)

Yields the number of explicitly set values in the persistent array p, regardless
of the number that are zero. This is different from the sparse-array library
only in that persistent arrays support arbitrary default values instead of
supporting only the value zero. Thus this counts explicit values instead of
non-zero values.

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"4-element (PArray{Int64, ?1}|PVector{Int64}):"
julia> u = PVector{Int}(0, (4,))
4-element PArray{Int64,1}:
 0
 0
 0
 0

julia> v = setindex(u, 2, 3)
4-element PArray{Int64,1}:
 0
 0
 2
 0

julia> nnz(v)
1

julia> nnz(u)
0
```
"""
SparseArrays.nnz(u::PArray{T,N}) where {T,N} = length(u._tree)
_dropzeros(p::PTree{T}, ::Nothing) where {T} = p
_dropzeros(p::PTree{T}, df::Tuple{T}) where {T} = begin
    df = df[1]
    for (k,v) in p
        if v == df
            p = delete(p, k)
        end
    end

    return p
end
"""
    dropzeros(p::PArray)

Drops explicit values of the given array `p` that are equal to the array's
default value. This differs from the SparseArrays implementation of dropzeros()
only in that `PArray`s allow arbitrary default values, while `SparseArray`s
allow only the default value of zero.

Note that under most circumstances, a `PArrray` will not encode explicit zeros,
so this function typically returns the object `p` untouched.

See also: `SparseArrays.nnz`, `SparseArrays.findnz`, [`PArray`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"4-element (PArray{Int64, ?1}|PVector{Int64}):"
julia> u = PVector{Int}([0,1,2,3])
4-element PArray{Int64,1}:
 0
 1
 2
 3

julia> dropzeros(u) === u
true

julia> v = setindex(u, 0, 3)
4-element PArray{Int64,1}:
 0
 1
 0
 3

julia> dropzeros(v) === v
true
```
"""
SparseArrays.dropzeros(p::PArray{T,N}) where {T,N} = begin
    t = _dropzeros(p._tree, p._default)
    (t === p._tree) && return p
    return PArray{T,N}(p._i0, p._index, t, p._default)
end
SparseArrays.dropzeros!(p::PArray{T,N}) where {T,N} = error(
    "dropzeros!: object of type $(typeof(p)) is immutable")
"""
    findnz(p::PArray)

Yields the explicitly set elements of the given persistent array `p`. This
method is identical to the typical SparseArrays implementation of findnz()
except that it respects the arbitrary default-value that persistent arrays are
allowed to have rather than assuming that this value is a zero, as is done in
the `SparseArray`s library.

Note that under most circumstances, a `PArray` will not encode explicit zeros,
so this function typically returns indices and values for all values that aren't
equal to the default value of the array `p` (which is zero by default).

See also: `SparseArrays.nnz`, [`PArray`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"4-element (PArray{Int64, ?1}|PVector{Int64}):"
julia> u = PVector([0,10,20,30])
4-element PArray{Int64,1}:
  0
 10
 20
 30

julia> findnz(u)
([1, 2, 3, 4], [0, 10, 20, 30])
```

```jldoctest; filter=r"4-element (PArray{Float64, ?1}|PVector{Float64}):"
julia> u = setindex(PVector(0.0, 4), 20, 2)
4-element PArray{Float64,1}:
  0.0
 20.0
  0.0
  0.0

julia> findnz(u)
([2], [20.0])
```
"""
SparseArrays.findnz(p::PArray{T,N}) where {T,N} = begin
    sz = size(p._index)
    ndims = length(sz)
    n = SparseArrays.nnz(p)
    iilists = [Vector{Int}(undef, n) for _ in 1:N]
    vals = Vector{T}(undef, n)
    cis = CartesianIndices(size(p))
    for (elno,(k,v)) in enumerate(p._tree)
        ci = cis[Int(k - p._i0) + 1]
        for (iilist,oo) in zip(iilists,Tuple(ci))
            iilist[elno] = oo
        end
        vals[elno] = v
    end
    return (iilists..., vals)
end
"""
    nonzeros(p::PArray)

Yields the explicitly set values of the given persistent array `p`. This method
is identical to the typical `SparseArrays` implementation of `nonzeros()` for
its sparse array classes except that it returns a persistent array of values and
that it respects the arbitrary default-value that persistent arrays are allowed
to have rather than assuming that this value is a zero, as is done in the
`SparseArrays` library.

Note that because `PArray`s don't typically store values equal to their default
value explicitly, this will typically yield a vector of every non-default value
in the array.

See also: `SparseArrays.findnz`, `SparseArrays.nnz`, [`PArray`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"4-element (PArray{Int64, ?1}|PVector{Int64}):"
julia> u = PVector([0,10,20,30])
4-element PArray{Int64,1}:
  0
 10
 20
 30

julia> nonzeros(u)
4-element PArray{Int64,1}:
  0
 10
 20
 30
```

```jldoctest; filter=r"[14]-element (PArray{Float64, ?1}|PVector{Float64}):"
julia> u = setindex(PVector(0.0, 4), 20, 2)
4-element PArray{Float64,1}:
  0.0
 20.0
  0.0
  0.0

julia> nonzeros(u)
1-element PArray{Float64,1}:
 20.0
```
"""
SparseArrays.nonzeros(p::PArray{T,N}) where {T,N} =
    PVector{T}(T[v for (k,v) in p._tree])
# PArrays are always considered sparse.
SparseArrays.issparse(::PArray) = true

# ==============================================================================
# Iterator methods

Base.IteratorSize(::Type{PArray{T,N}}) where {T,N} = Base.HasShape{N}()
Base.IteratorEltype(::Type{PArray{T,N}}) where {T,N} = Base.HasEltype()

Base.length(u::PArray{T,N}) where {T,N} = length(u._index)
Base.size(u::PArray{T,N}) where {T,N} = size(u._index)
Base.eltype(::Type{PArray{T,N}}) where {T,N} = T
Base.eltype(u::PArray{T,N}) where {T,N} = T

Base.iterate(u::PArray{T,N}, k::Int) where {T,N} = begin
    return k > length(u) ? nothing : (u[k], k+1)
end
Base.iterate(u::PArray{T,N}) where {T,N} = iterate(u, 1)

# ==============================================================================
# Indexing methods

IndexStyle(::Type{PArray{T,N}}) where {T,N} = IndexLinear()
_parray_get(u::PTree{T}, ii::HASH_T, ::Nothing) where {T} = begin
    x = get(u, ii, nothing)
    (x === nothing) || return x
    isa(nothing, T) && in(ii => nothing, u) && return x
    #error("PArray has unset values and no default")
    return Array{T}(undef, 1)[1]
end
_parray_get(u::PTree{T}, ii::HASH_T, d::Tuple{T}) where {T} = begin
    return get(u, ii, d[1])
end
Base.getindex(u::PArray{T,N}, k::Vararg{Int,N}) where {T,N} = begin
    if N > 1
        k = u._index[k...]
    else
        k = k[1]
    end
    (k < 1) && throw(BoundsError(u, k))
    (k > length(u)) && throw(BoundsError(u, k))
    return _parray_get(u._tree, HASH_T(k - 1) + u._i0, u._default)
end
Base.setindex!(u::PArray{T,N}, v, k::Int) where {T,N} = error(
    "setindex!: object of type $(typeof(u)) is immutable; see setindex()")
Base.firstindex(u::PArray{T,N}) where {T,N} = 1
Base.lastindex(u::PArray{T,N}) where {T,N} = length(u)

# ==============================================================================
# AbstractArray methods

Base.push!(u::PArray{T,N}, v) where {T,N} = error(
    "push!: object of type $(typeof(u)) is immutable; see push()")
Base.pop!(u::PArray{T,N}, v) where {T,N} = error(
    "pop!: object of type $(typeof(u)) is immutable; see last() and pop()")
#Base.pushfist!(u::PArray{T,N}, v) where {T,N} = error(
#    "pushfirst!: object of type $(typeof(u)) is immutable; see first() pushfirst()")
#Base.popfirst!(u::PArray{T,N}, v) where {T,N} = error(
#    "popfist!: object of type $(typeof(u)) is immutable; see first() and popfist()")
_lindex(u::Vararg{Int,N}) where {N} = LinearIndices{N,NTuple{N,Base.OneTo{Int}}}(
    ((Base.OneTo{Int}.(u))...,))
Base.permutedims(u::PArray{T,N}, dims::NTuple{N,Int}) where {T,N} = begin
    # This is actually pretty easy, code-wise:
    v = PVector{T}(0x0, LinearIndices(()), PTree{T}(), u._default)
    for k in permutedims(u._index, dims)
        v = push(v, u[k])
    end
    sz = size(u)
    idx = _lindex([sz[k] for k in dims]...)
    return PArray{T,N}(v._i0, idx, v._tree, u._default)
end
#Base.broadcast(fn::F, u::PArray{T,N}, args...) where {F<:Function,T,N} = begin
#
#end

# ==============================================================================
# Persistent array methods

_defaultvalue(::Nothing) = undef
_defaultvalue(u::Tuple{T}) where {T} = u[1]
_eqdefault(::Nothing, x) = false
_eqdefault(dflt::Tuple{T}, x::S) where {T,S} = (dflt[1] == x)
defaultvalue(u::PArray{T,N}) where {T,N} = _defaultvalue(u._default)
setindex(u::PArray{T,N}, v::S, ci::CartesianIndex{N}) where {T,N,S} = begin
    return setindex(u, v, u._index[ci])
end
setindex(u::PArray{T,N}, v::S, k::Vararg{Idx,N}) where {T,N,S,Idx<:Integer} = begin
    if N == 1
        k = k[1]
    else
        k = u._index[k...]
    end
    (k < 1) && throw(BoundsError(u,k))
    n = length(u)
    (k > n + 1) && throw(BoundsError(u,k))
    (N == 1) && (k > n) && return push(u, v)
    ii = u._i0 + HASH_T(k - 1)
    t = (_eqdefault(u._default, v) ? delete(u._tree, ii)
                                   : setindex(u._tree, v, ii))
    return t === u._tree ? u : PArray{T,N}(u._i0, u._index, t, u._default)
end
setindex(u::PArray{T,N}, v::S, ii...) where {T,N,S} = begin
    idcs = getindex(u._index, ii...)
    pp = broadcast(Pair{Int,T}, idcs, v)
    if isa(pp, Pair)
        return setindex(u, pp[2], pp[1])
    else
        for (k,v) in pp
            u = setindex(u, v, k)
        end
        return u
    end
end
# Push and pop methods are only defined for vectors
push(u::PVector{T}, x::S) where {T,S} = begin
    n = length(u)
    if _eqdefault(u._default, x)
        tree = u._tree
    else
        ii   = u._i0 + HASH_T(n + 1 - 1)
        tree = setindex(u._tree, x, ii)
    end
    index = _lindex(n+1)
    return PVector{T}(u._i0, index, tree, u._default)
end
pushfirst(u::PVector{T}, x::S) where {T,S} = begin
    n = length(u)
    if _eqdefault(u._default, x)
        tree = u._tree
    else
        tree = setindex(u._tree, x, u._i0 - 0x1)
    end
    index = _lindex(n+1)
    return PVector{T}(u._i0 - 0x1, index, tree, u._default)
end
pop(u::PVector{T}) where {T} = begin
    n = length(u)
    (n == 0) && throw(ArgumentError("PArray must be non-empty"))
    ii = u._i0 + HASH_T(n - 1)
    tree = delete(u._tree, ii)
    return PVector{T}(u._i0, _lindex(n-1), tree, u._default)
end
popfirst(u::PVector{T}) where {T} = begin
    n = length(u)
    (n == 0) && throw(ArgumentError("PArray must be non-empty"))
    tree = delete(u._tree, u._i0)
    return PVector{T}(u._i0 + 0x1, _lindex(n-1), tree, u._default)
end

# pzeros, pones, and other useful utility functions.
"""
    pzeros(dims...)
    pzeros(T, dims...)

Yields a persistent array (`PArray`) of zeros exactly as done by the `zeros()`
function.

See also [`pones`](@ref), [`pfill`](@ref), `zeros`

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"3-element (PArray{Integer, ?1}|PVector{Integer}):"
julia> pzeros(Integer, 3)
3-element PArray{Integer,1}:
 0
 0
 0
```

```jldoctest; filter=r"1×2 (PArray{Bool, ?2}|PMatrix{Bool}):"
julia> pzeros(Bool, (1,2))
1×2 PArray{Bool,2}:
 0  0
```

```jldoctest; filter=r"2×3 (PArray{Float64, ?2}|PMatrix{Float64}):"
julia> pzeros(2, 3)
2×3 PArray{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
```

```jldoctest; filter=r"1×1×1×1 PArray{Float64, ?4}:"
julia> pzeros((1, 1, 1, 1))
1×1×1×1 PArray{Float64,4}:
[:, :, 1, 1] =
 0.0
```
"""
pzeros(::Type{T}, dims::Vararg{Integer}) where {T} = PArray{T,length(dims)}(0, dims...)
pzeros(::Type{T}, dims::Tuple) where {T} = PArray{T,length(dims)}(0, dims...)
pzeros(dims::Tuple) = PArray{Float64,length(dims)}(0.0, dims)
pzeros(dims::Vararg{Integer}) = PArray{Float64,length(dims)}(0.0, dims...)

"""
    pones(dims...)
    pones(T, dims...)

Yields a persistent array (`PArray`) of ones exactly as done by the `ones()`
function.

See also [`pzeros`](@ref), [`pfill`](@ref), `ones`

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"3-element (PArray{Integer, ?1}|PVector{Integer}):"
julia> pones(Integer, 3)
3-element PArray{Integer,1}:
 1
 1
 1
```

```jldoctest; filter=r"1×2 (PArray{Bool, ?2}|PMatrix{Bool}):"
julia> pones(Bool, (1,2))
1×2 PArray{Bool,2}:
 1  1
```

```jldoctest; filter=r"2×3 (PArray{Float64, ?2}|PMatrix{Float64}):"
julia> pones(2, 3)
2×3 PArray{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
```

```jldoctest; filter=r"1×1×1×1 PArray{Float64, ?4}:"
julia> pones((1, 1, 1, 1))
1×1×1×1 PArray{Float64,4}:
[:, :, 1, 1] =
 1.0
```
"""
pones(::Type{T}, dims...) where {T} = PArray{T,length(dims)}(1, dims...)
pones(::Type{T}, dims::Tuple) where {T} = PArray{T,length(dims)}(1, dims...)
pones(dims::Tuple) = PArray{Float64,length(dims)}(1.0, dims)
pones(dims::Vararg{Integer}) = PArray{Float64,length(dims)}(1.0, dims...)

"""
    pfill(val, dims...)

Yields a persistent array (`PArray`) of values exactly as done by the `fill()`
function.

See also [`pones`](@ref), [`pzeros`](@ref), `fill`

# Examples

```@meta
DocTestSetup = quote
    using Air, SparseArrays
end
```

```jldoctest; filter=r"3-element (PArray{Float64, ?1}|PVector{Float64}):"
julia> pfill(NaN, 3)
3-element PArray{Float64,1}:
 NaN
 NaN
 NaN
```

```jldoctest; filter=r"2×3 (PArray{Symbol, ?2}|PMatrix{Symbol}):"
julia> pfill(:abc, 2, 3)
2×3 PArray{Symbol,2}:
 :abc  :abc  :abc
 :abc  :abc  :abc
```
"""
pfill(val::T, dims::Vararg{Integer}) where {T} = PArray{T,length(dims)}(val, dims...)
pfill(val::T, dims::Tuple) where {T} = PArray{T,length(dims)}(val, dims)

export PArray, PVector, pzeros, pones, pfill


text:
import FLOWFarm; const ff = FLOWFarm
using DelimitedFiles

# Defining Turbine Parameters and locations
D = 126.4
locx = 7D
locy = 0
locz = 90
turbine_x = [7D,0]
turbine_y = [0,0]
turbine_z = [0,0]
nturbines = length(turbine_x)

upstream_turbine_id = 2
downstream_turbine_id = 1
hub_height = [90.0,90.0]
rotor_diameter = [D,D]
turbine_ai = [1.0/3.0,1.0/3.0]
turbine_yaw = [0, 0]
dt = rotor_diameter[upstream_turbine_id]
cut_in_speed = [0.0, 0.0]
cut_out_speed = [25.0, 25.0]
rated_speed = [12.0, 12.0]
rated_power = [1.0176371581904552e6 , 1.0176371581904552e6]
#            turbine_yaw = [-40,-20,0,20,40]

# define power parameters
u = 8.0
air_density = 1.225
generator_efficiency = [.944, .944]
rotor_area = .25*pi*D^2
Fd = 0
ai = 1.0/3.0
ambient_ti = .06

turbine_local_ti = [ambient_ti,ambient_ti]

# values defining wake interaction from Gebraad et. al
me = [-.5,.22,1]
ke = .065
MU = [.5,1,5.5]
aU = 5
bU = 1.66

# identify wake models
wakedeficitmodel = ff.Multizone()
wakedeflectionmodel = ff.MultizoneDeflection()
wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
localtimodel = ff.LocalTIModelMaxTI()
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)

# define wind characteristics
winddirections = [0.0]
windspeeds = [u]
windprobabilities = [1.0]
measurementheights = [hub_height[1]]
ambient_tis = [ambient_ti]
shearexponent = 0.15
wind_shear_model = ff.PowerLawWindShear(shearexponent)
wind_resource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)

# Sample wind point at the hub
rotor_sample_points_y = [0]
rotor_sample_points_z = [0]
sorted_turbine_index = sortperm(turbine_x)

eta = .768

# import cp and ct corresponding to wind velocity
powerdata = readdlm("inputfiles/NREL5MWCPCT.txt", ' ', skipstart=1)
vel_points = powerdata[:,1]
cp_points = powerdata[:,2]
ct_points = powerdata[:,3]

# initialize power model
power_model = ff.PowerModelCpPoints(vel_points, cp_points)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# initialize thurst model
ct_model = ff.ThrustModelCtPoints(vel_points, ct_points)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end


text:
type ZeroElementException <: Exception end

#macro debug_stmt(x)
#  x
#end

macro debug_stmt(x)
end

# used to send some size information.
immutable TableSize{N}
  size::NTuple{N,Int}
end
TableSize(arr::AbstractArray) = TableSize(size(arr))

# some statements for compatibility.

iosize_expr = if startswith(string(VERSION), "0.5")
  :(iosize_compat = Base.displaysize)
  #:(iosize_compat = Base.iosize)
else
  :(iosize_compat = Base.tty_size)
end

eval(iosize_expr)


text:
# x is date relative to Jan 1, with 1.00 representing Jan 1 00:00
function xtojandate(x, year)
    basedate = Dates.DateTime(year, 1, 1)
    deltat = Dates.Second(round((x-1)*Planet.day(param_set)))
    date = basedate + deltat
    return date
end

# Earth-Sun distance
function edist(x, year)
    date = xtojandate(x,year)
    _, dist = daily_zenith_angle(date, FT(0), param_set)
    return dist/orbit_semimaj(param_set)
end

years = 1900:2100
days = zeros(length(years))
for (i,year) in enumerate(years)
    f = (x -> edist(x, year))
    res = optimize(f,1.,30)
    days[i] = Optim.minimizer(res)[1]
end

# test mean is about Jan 3.5
@test mean(days) ≈ 3.5 atol=1

# test increasing
@test mean(days[:100]) < mean(days[100:end])


text:
module CESMReader

    using NCDatasets
    using Formatting
    export FileHandler, getData

    include("nanop.jl")

    mutable struct FileHandler
        filename_format :: String
        form   :: Symbol
        LENS_benchmark :: Bool

        function FileHandler(;
            filename_format :: String,
            form        :: Symbol = :YEAR_MONTH,
            LENS_benchmark :: Bool = false,
        )

            if ! (form in [:YEAR, :YEAR_MONTH]  )
                throw(ErrorException("Error: only two forms allowed. Either :YEAR or :YEAR_MONTH. Now we got " * string(form)))
            end

            return new(
                filename_format,
                form,
                LENS_benchmark,
            )
        end

    end

    function getData(
        handler  :: FileHandler,
        varnames :: Union{String, Array},
        year_rng :: Union{Tuple, Array},
        idxes...;
        verbose=true,
    )
        return _getData(handler, varnames, _genYearMonth(year_rng[1], year_rng[2]), collect(idxes)...; verbose=verbose)
    end

    function _genYearMonth(
        beg_y :: Int64,
        end_y :: Int64,
    )
        #=
        beg_t = (beg_y - 1) * 12 + beg_m - 1
        end_t = (end_y - 1) * 12 + end_m - 1

        months = end_t - beg_t + 1

        if months <= 0
            throw(ErrorException("End time should be larger than begin time"))
        end

        if ! ( ( 1 <= beg_m <= 12 ) && ( 1 <= end_m <= 12 ) )
            throw(ErrorException("Invalid month"))
        end
        =#

        time_vec = Array{Any}(undef, 0)

        for y = beg_y:end_y
            for m = 1:12
                push!(time_vec, (y, m))
            end
        end

        return time_vec
    end

    function _getData(
        handler  :: FileHandler,
        varnames :: Union{String, Array},
        time_vec :: Array,
        idxes...;
        verbose=true,
    )

        local only_one_variable = false

        if typeof(varnames) <: AbstractString
            only_one_variable = true
            varnames = [varnames]
        end

        if ! only_one_variable && length(idxes) == 1
            idxes = [ idxes[1] for _ in 1:length(varnames) ]
        end

        if ( only_one_variable && length(idxes) > 1 ) || length(idxes) != length(varnames)
            println(varnames)
            println(idxes)
            throw(ErrorException("Number of idxes should match number of varnames or only 1 (applied to all varnames)."))
        end

        for k = 1:length(idxes)
            if length(idxes[k]) == 0
                throw(ErrorException("Spatial range not specified."))
            end
        end

        local data, new_idxes

        datas = Array{Any}(undef, length(varnames))
        new_idxes = Array{Any}(undef, length(varnames))

        if handler.form == :YEAR_MONTH

            flag_1s = [true for _ in 1:length(varnames)]

            current_t = 0
            for (y, m) in time_vec

                current_t += 1

                filename = format(handler.filename_format, y, m)
                verbose && println("Loading file: ", filename)
                ds = Dataset(filename, "r")
                for (k, varname) in enumerate(varnames)

                    partial_data = replace(ds[varname][idxes[k]..., 1], missing=>NaN)
                    if flag_1s[k]
                        new_size = size(partial_data)
                        new_idxes[k]  = [Colon() for _ in 1:length(new_size)]
                        datas[k] = zeros(Float64, new_size..., length(time_vec))
                        datas[k] .= NaN
                        flag_1s[k] = false
                    end

                    datas[k][new_idxes[k]..., current_t] = partial_data

                end

                close(ds)

            end

        end

        return (only_one_variable) ? datas[1] : datas
    end

    #=

    function _getData(
        handler  :: FileHandler,
        varnames :: Union{String, Array},
        beg_time :: Union{Tuple, Array},
        end_time :: Union{Tuple, Array},
        idxes...;
        verbose=true,
    )

        local only_one_variable = false

        if typeof(varnames) <: AbstractString
            only_one_variable = true
            varnames = [varnames]
        end

        if ! only_one_variable && length(idxes) == 1
            idxes = [ idxes[1] for _ in 1:length(varnames) ]
        end

        if ( only_one_variable && length(idxes) > 1 ) || length(idxes) != length(varnames)
            println(varnames)
            println(idxes)
            throw(ErrorException("Number of idxes should match number of varnames or only 1 (applied to all varnames)."))
        end

        beg_y, beg_m = beg_time
        end_y, end_m = end_time

        beg_t = (beg_y - 1) * 12 + beg_m - 1
        end_t = (end_y - 1) * 12 + end_m - 1

        months = end_t - beg_t + 1

        if months <= 0
            throw(ErrorException("End time should be larger than begin time"))
        end

        if ! ( ( 1 <= beg_m <= 12 ) && ( 1 <= end_m <= 12 ) )
            throw(ErrorException("Invalid month"))
        end

        for k = 1:length(idxes)
            if length(idxes[k]) == 0
                throw(ErrorException("Spatial range not specified."))
            end
        end

        local data, new_idxes

        datas = Array{Any}(undef, length(varnames))
        new_idxes = Array{Any}(undef, length(varnames))

        if handler.form == :YEAR_MONTH

            flag_1s = [true for _ in 1:length(varnames)]

            for y in beg_y:end_y, m in beg_m:end_m

                current_t = (y-1) * 12 + (m-1) - beg_t + 1

                filename = format(handler.filename_format, y, m)
                verbose && println("Loading file: ", filename)
                ds = Dataset(filename, "r")

                for (k, varname) in enumerate(varnames)

                    partial_data = replace(ds[varname][idxes[k]..., 1], missing=>NaN)
                    if flag_1s[k]
                        new_size = size(partial_data)
                        new_idxes[k]  = [Colon() for _ in 1:length(new_size)]
                        datas[k] = zeros(Float64, new_size..., months)
                        datas[k] .= NaN
                        flag_1s[k] = false
                    end

                    datas[k][new_idxes[k]..., current_t] = partial_data

                end

                close(ds)

            end

        elseif handler.form == :YEAR

            flag_1s = [true for _ in 1:length(varnames)]

            for y in beg_y:end_y

                current_t = (y-1) * 12 - beg_t

                filename = format(handler.filename_format, y)
                verbose && println("Loading file: ", filename)
                ds = Dataset(filename, "r")

                for (k, varname) in enumerate(varnames)

                    if y == beg_y
                        rng = (idxes[k]..., beg_m:12)
                    elseif y == end_y
                        rng = (idxes[k]..., 1:end_m)
                    else
                        rng = (idxes[k]..., 1:12)
                    end

                    partial_data = nomissing(ds[varname][rng...])

                    if flag_1s[k]
                        new_size = size(partial_data)[1:end-1]
                        new_idxes[k]  = [Colon() for _ in 1:length(new_size)]
                        datas[k] = zeros(Float64, new_size..., months)
                        datas[k] .= NaN
                        flag_1s[k] = false
                    end

                    if y == beg_y
                        rng = beg_m:12
                    elseif y == end_y
                        rng = 1:end_m
                    else
                        rng = 1:12
                    end

                    offset = (y-beg_y) * 12
                    datas[k][new_idxes[k]..., rng .+ offset] = partial_data

                end

                close(ds)

            end

        end

        return (only_one_variable) ? datas[1] : datas
    end
    =#
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
#=
utilities for complex polynomial interpolation

TODO make interface for complex and real interpolation consistent, use dispatch
=#

"""
$(SIGNATURES)

Compute interpolation data for a complex weighted sum-of-squares conic constraint
on a domain.
"""
function interpolate(
    R::Type{Complex{T}},
    halfdeg::Int,
    n::Int,
    gs::Vector,
    g_halfdegs::Vector{Int};
    sample_factor::Int = 10,
    use_qr::Bool = false,
    ) where {T <: Real}
    # generate interpolation
    # TODO use more numerically-stable basis for columns, and evaluate in a more
    # numerically stable way by multiplying the columns
    L = binomial(n + halfdeg, n)
    U = L^2
    L_basis = [a for t in 0:halfdeg for a in Combinatorics.multiexponents(n, t)]
    mon_pow(z, ex) = prod(z[i]^ex[i] for i in eachindex(ex))
    V_basis = [z -> mon_pow(z, L_basis[k]) * mon_pow(conj(z), L_basis[l]) for
        l in eachindex(L_basis) for k in eachindex(L_basis)]
    @assert length(V_basis) == U

    # sample from domain (inefficient for general domains, only samples from
    # unit box and checks feasibility)
    num_samples = sample_factor * U
    samples = Vector{Vector{Complex{T}}}(undef, num_samples)
    k = 0
    randbox() = 2 * rand(T) - 1
    while k < num_samples
        z = [Complex(randbox(), randbox()) for i in 1:n]
        if all(g -> g(z) > zero(T), gs)
            k += 1
            samples[k] = z
        end
    end

    # select subset of points to maximize abs(det(V)) in heuristic QR-based
    # procedure (analogous to real case)
    V = [b(z) for z in samples, b in V_basis]
    VF = qr(Matrix(transpose(V)), ColumnNorm())
    keep = VF.p[1:U]
    points = samples[keep]
    V = V[keep, :]

    # setup P matrices
    P0 = V[:, 1:L]
    if use_qr
        P0 = Matrix(qr(P0).Q)
    end
    Ps = [P0]
    for i in eachindex(gs)
        gi = gs[i].(points)
        @views P0i = P0[:, 1:binomial(n + halfdeg - g_halfdegs[i], n)]
        Pi = Diagonal(sqrt.(gi)) * P0i
        if use_qr
            Pi = Matrix(qr(Pi).Q)
        end
        push!(Ps, Pi)
    end

    return (points, Ps)
end


text:
function hamming_distance{T}(::Type{T}, a, b)
    @assert length(a) == length(b)
    d = T(0)
    for i in 1:length(a)
        d += ifelse(a[i] == b[i], T(0), T(1))
    end
    return d
end


text:
module SemiDiscretizationMethod

using Reexport
@reexport using LinearAlgebra
@reexport using SparseArrays
@reexport using StaticArrays
@reexport using Arpack
using QuadGK
using Lazy: iterated, take

include("structures_method.jl")
include("structures_input.jl")
include("structures_result.jl")

include("functions_utility.jl")
include("functions_discretization.jl")
include("functions_method.jl")

export SemiDiscretization, NumericSD,
ProportionalMX,
Delay,DelayMX,
Additive,
LDDEProblem,
DiscreteMapping, DiscreteMapping_1step,
fixPointOfMapping, spectralRadiusOfMapping

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiPolicyV1beta1HostPortRange <: SwaggerModel
    max::Any # spec type: Union{ Nothing, Int32 } # spec name: max
    min::Any # spec type: Union{ Nothing, Int32 } # spec name: min

    function IoK8sApiPolicyV1beta1HostPortRange(;max=nothing, min=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1HostPortRange, Symbol("max"), max)
        setfield!(o, Symbol("max"), max)
        validate_property(IoK8sApiPolicyV1beta1HostPortRange, Symbol("min"), min)
        setfield!(o, Symbol("min"), min)
        o
    end
end # type IoK8sApiPolicyV1beta1HostPortRange

const _property_map_IoK8sApiPolicyV1beta1HostPortRange = Dict{Symbol,Symbol}(Symbol("max")=>Symbol("max"), Symbol("min")=>Symbol("min"))
const _property_types_IoK8sApiPolicyV1beta1HostPortRange = Dict{Symbol,String}(Symbol("max")=>"Int32", Symbol("min")=>"Int32")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1HostPortRange }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1HostPortRange))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1HostPortRange }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1HostPortRange[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1HostPortRange }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1HostPortRange[property_name]

function check_required(o::IoK8sApiPolicyV1beta1HostPortRange)
    (getproperty(o, Symbol("max")) === nothing) && (return false)
    (getproperty(o, Symbol("min")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1HostPortRange }, name::Symbol, val)
end


text:
mutable struct Capture{V,T,R} <: State{V}
    value::V
    time::T
    rate::R
end

Capture(; unit, time, timeunit, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    v = unitfy(zero(_type), U)
    #V = promote_type(V, typeof(v))
    TU = value(timeunit)
    t = unitfy(value(time), TU)
    T = typeof(t)
    RU = rateunittype(U, TU)
    R = valuetype(_type, RU)
    Capture{V,T,R}(v, t, zero(R))
end

constructortags(::Val{:Capture}) = (:unit, :time, :timeunit)

@generated rateunit(s::Capture{V,T,R}) where {V,T,R} = unittype(R)

updatetags!(d, ::Val{:Capture}; _...) = begin
    !haskey(d, :time) && (d[:time] = :(context.clock.time))
    #TODO: automatic inference without explicit `timeunit` tag
    !haskey(d, :timeunit) && (d[:timeunit] = @q $C.timeunittype($(d[:unit]), $C.timeunit(__Context__)))
end

genvartype(v::VarInfo, ::Val{:Capture}; N, U, V, _...) = begin
    TU = gettag(v, :timeunit)
    T = @q $C.valuetype(Float64, $TU)
    RU = @q $C.rateunittype($U, $TU)
    R = @q $C.valuetype($N, $RU)
    @q Capture{$V,$T,$R}
end

gendefault(v::VarInfo, ::Val{:Capture}) = nothing

genupdate(v::VarInfo, ::Val{:Capture}, ::MainStep; kw...) = begin
    @gensym s t t0 d
    @q let $s = $(symstate(v)),
           $t = $C.value($(gettag(v, :time))),
           $t0 = $s.time,
           $d = $s.rate * ($t - $t0)
        $C.store!($s, $d)
    end
end

genupdate(v::VarInfo, ::Val{:Capture}, ::PostStep; kw...) = begin
    w = gettag(v, :when)
    f = isnothing(w) ? genbody(v) : @q $C.value($w) ? $(genbody(v)) : zero($(gettag(v, :_type)))
    @gensym s t r
    @q let $s = $(symstate(v)),
           $t = $C.value($(gettag(v, :time))),
           $r = $C.unitfy($f, $C.rateunit($s))
        $s.time = $t
        $s.rate = $r
    end
end


text:
# GENERAL RULES FOR INTERSECTION AND UNIONS

using Base.Iterators: product

struct Intersection{T <: Surface} <: Surface
    content::Vector{T}

    # function Intersection{T}(content) where T
    #     content = unique(content)
    #     if any(isempty(h) for h in content)
    #         return new([])
    #     else
    #         return new(content)
    #     end
    # end
end

==(i1::Intersection, i2::Intersection) = all(h1 == h2 for (h1, h2) in zip(i1.content, i2.content))

struct Reunion{T <: Surface} <: Surface
    content::Vector{T}

    Reunion{T}(content) where T = new(filter(!isempty, unique(content)))
end

# Generic fonctions

in(p, i::Intersection) = all((p in h) for h in i.content)
in(p, u::Reunion) = any((p in h) for h in u.content)

isempty(i::Intersection) = any(isempty(h) for h in i.content)
isempty(u::Reunion) = all(isempty(h) for h in u.content)

complement(i::Intersection) = foldl(union, (complement(h) for h in i.content))
complement(u::Reunion) = foldl(intersect, (complement(h) for h in u.content))

\(c1::Surface, c2::Surface) = c1 ∩ complement(c2)
\(c1::Surface, u2::Reunion) = foldl(\, u2.content, init=c1)

# Union and intersection of half planes

convert(::Type{Intersection{HalfPlane{T}}}, h::HalfPlane) where T = Intersection{HalfPlane{T}}([h])
convert(::Type{Intersection{HalfPlane{T}}}, i::Intersection{HalfPlane{U}}) where {T, U} = Intersection{HalfPlane{T}}([convert(HalfPlane{T}, h) for h in i.content])

convert(::Type{Reunion{HalfPlane{T}}}, h::HalfPlane) where T = Reunion{HalfPlane{T}}([h])
convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, h::HalfPlane) where T = Reunion{Intersection{HalfPlane{T}}}([Intersection{HalfPlane{T}}([h])])
convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, i::Intersection{HalfPlane{U}}) where {T, U} = Reunion{Intersection{HalfPlane{T}}}([i])
convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, u::Reunion{HalfPlane{U}}) where {T, U} = Reunion{Intersection{HalfPlane{T}}}([Intersection{HalfPlane{T}}([h]) for h in u.content])

convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, u::Reunion{Intersection{HalfPlane{U}}}) where {T, U} = Reunion{Intersection{HalfPlane{T}}}([convert(Intersection{HalfPlane{T}}, i) for i in u.content])

promote_rule(::Type{HalfPlane{T}}, ::Type{Intersection{HalfPlane{U}}}) where {T, U} = Intersection{HalfPlane{promote_type(T, U)}}
promote_rule(::Type{HalfPlane{T}}, ::Type{Reunion{HalfPlane{U}}}) where {T, U} = Reunion{HalfPlane{promote_type(T, U)}}
promote_rule(::Type{HalfPlane{T}}, ::Type{Reunion{Intersection{HalfPlane{U}}}}) where {T, U} = Reunion{Intersection{HalfPlane{promote_type(T, U)}}}
promote_rule(::Type{Intersection{HalfPlane{T}}}, ::Type{Reunion{Intersection{HalfPlane{U}}}}) where {T, U} = Reunion{Intersection{HalfPlane{promote_type(T, U)}}}
promote_rule(::Type{Reunion{HalfPlane{T}}}, ::Type{Reunion{Intersection{HalfPlane{U}}}}) where {T, U} = Reunion{Intersection{HalfPlane{promote_type(T, U)}}}
promote_rule(::Type{Reunion{HalfPlane{T}}}, ::Type{Intersection{HalfPlane{U}}}) where {T, U} = Reunion{Intersection{HalfPlane{promote_type(T, U)}}}

union(s1::Surface, s2::Surface) = union(promote(s1, s2)...)
intersect(s1::Surface, s2::Surface) = intersect(promote(s1, s2)...)

union(h1::HalfPlane{T}, h2::HalfPlane{T}) where T = Reunion{HalfPlane{T}}([h1, h2])
union(u1::Reunion{HalfPlane{T}}, u2::Reunion{HalfPlane{T}}) where T = Reunion{HalfPlane{T}}(vcat(u1.content, u2.content))
union(i1::Intersection{HalfPlane{T}}, i2::Intersection{HalfPlane{T}}) where T = Reunion{Intersection{HalfPlane{T}}}([i1, i2])
union(ui1::Reunion{Intersection{HalfPlane{T}}}, ui2::Reunion{Intersection{HalfPlane{T}}}) where T = Reunion{Intersection{HalfPlane{T}}}(vcat(ui1.content, ui2.content))

intersect(h1::HalfPlane{T}, h2::HalfPlane{T}) where T = Intersection{HalfPlane{T}}([h1, h2])
intersect(i1::Intersection{HalfPlane{T}}, i2::Intersection{HalfPlane{T}}) where T = Intersection{HalfPlane{T}}(vcat(i1.content, i2.content))
intersect(u1::Reunion{HalfPlane{T}}, u2::Reunion{HalfPlane{T}}) where T = intersect(convert(Reunion{Intersection{HalfPlane{T}}}, u1), convert(Reunion{Intersection{HalfPlane{T}}}, u2))
function intersect(ui1::Reunion{Intersection{HalfPlane{T}}}, ui2::Reunion{Intersection{HalfPlane{T}}}) where T
    inters = []
    for (i1, i2) in product(ui1.content, ui2.content)
        push!(inters, Intersection{HalfPlane{T}}(vcat(i1.content, i2.content)))
    end
    return Reunion{Intersection{HalfPlane{T}}}(inters)
end

rotate(u::Reunion{T}, ϕ; kw...) where T = Reunion{T}(map(x -> rotate(x, ϕ; kw...), u.content))
rotate(i::Intersection{T}, ϕ; kw...) where T = Intersection{T}(map(x -> rotate(x, ϕ; kw...), i.content))

translate(u::Reunion{T}, v) where T = Reunion{T}(map(x -> translate(x, v), u.content))
translate(i::Intersection{T}, v) where T = Intersection{T}(map(x -> translate(x, v), i.content))

scale(u::Reunion{T}, λ; center=Point(0.0, 0.0)) where T = Reunion{T}(map(x -> scale(x, λ; center), u.content))
scale(i::Intersection{T}, λ; center=Point(0.0, 0.0)) where T = Intersection{T}(map(x -> scale(x, λ; center), i.content))


text:
using PyPlot, Seismic

dt = 0.002
w = Ricker(dt=dt)
wmin = SeisKolmogoroff(w)

nw = length(w)
t = dt*collect(0:1:nw-1)

figure("Kolmogoroff factorization", figsize=(6, 6))

subplot(211)
plot(t, w)
axis("tight")
xlabel("Time (s)")
title("Ricker wavelet")

subplot(212)
plot(t,wmin)
axis("tight")
xlabel("Time (s)")
title("Minimum Phase equivalent")

tight_layout()


text:
const UIntOrChar = Union{Unsigned, AbstractChar}

struct StaticString{N, T<:Unsigned} <: AbstractString
    data::NTuple{N, T}
    function StaticString{N, T}(t::NTuple{M, <:UIntOrChar}) where {N, T, M}
        N == M || throw(DimensionMismatch(
            "cannot construct StaticString{$N, $T} from input of length $M"))
        new{N, T}(t)
    end
end

StaticString{N, T}(cs::UIntOrChar...) where {N, T} = StaticString{N, T}(cs)
StaticString{N, T}(fs::StaticString{N, T}) where {N, T} = fs
StaticString{N, T}(s) where {N, T} = StaticString{N, T}(String(s)...)

_units_to_type(n::Int) = ifelse(n == 1, UInt8, ifelse(n==2, UInt16, UInt32))

macro static_str(s)
    quote
        StaticString{
            length($(esc(s))),
            _units_to_type(maximum(ncodeunits, $(esc(s))))
        }($(esc(s)))
    end
end

convert(::Type{StaticString{N, T}}, s::AbstractString) where {N, T} = StaticString{N, T}(s)

ncodeunits(s::StaticString{N, T}) where {T,N} = N

sizeof(s::StaticString) = sizeof(s.data)

length(s::StaticString) = ncodeunits(s)
length(S::Type{StaticString{N, T}}) where {N, T} = N

lastindex(s::StaticString{N, T}) where {N, T} = N

function iterate(s::StaticString{N, T}, i::Int = 1) where {N, T}
    i > N && return nothing
    return Char(s[i]), i+1
end

codeunit(::StaticString{N, T}) where {N, T} = T
@propagate_inbounds codeunit(s::StaticString, i::Integer) = s.data[i]
@propagate_inbounds getindex(s::StaticString, i::Integer)::Char = s.data[i]

isvalid(s::StaticString, i::Integer) = checkbounds(Bool, s, i)