text:
using Random, PyPlot, PyCall, DataFrames

"""
    PlotMapProj(df::DataFrame,nn::Integer)

Plot random subset of size nn trajectories using PyPlot & basemap.
"""
function PlotMapProj(df::DataFrame,nn::Integer)

   #PyPlot.figure()

   # Set up Equidistant cylindrical map projection. Use low resolution coastlines.
   ccrs=pyimport("cartopy.crs")
   #ax = plt.axes(projection=ccrs.Orthographic(central_latitude=-90.0,central_longitude=-160.0))
   ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-160.0))
   ax.coastlines(linewidth=0.3)
   ax.stock_img()

   # Draw trajectories
   IDs = randperm(maximum(df.ID))
   COs=["w" "y" "g" "k"]

   #for global ocean case:
   dMax=90.
   println("dMax=$dMax")

   for ii=1:nn
      tmp=df[df.ID .== IDs[ii], :]
      jj=findall(tmp[!,:lon] .< 20)
      tmp[jj,:lon]=tmp[jj,:lon] .+ 360.0;
      if dMax > 0.
         d=abs.(diff(tmp[!,:lon]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
         d=abs.(diff(tmp[!,:lat]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
      end
      #
      CO=COs[mod(ii,4)+1]
      plt.plot(tmp[!,:lon], tmp[!,:lat], color=CO, linewidth=0.2, transform=ccrs.Geodetic() )
   end

   #to display the figure in JUNO (not needed in REPL or Jupyter):
   #gcf()

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
struct PV <: AbstractTech
    tilt
    array_type
    module_type
    losses
    azimuth
    name
    location
    existing_kw
    min_kw
    max_kw
    cost_per_kw
    om_cost_per_kw
    degradation_pct
    macrs_option_years
    macrs_bonus_pct
    macrs_itc_reduction
    total_itc_pct
    total_rebate_per_kw
    kw_per_square_foot
    acres_per_kw
    inv_eff
    dc_ac_ratio
    prod_factor_series_kw

    function PV(;
        tilt::Real,
        array_type::Int=1,
        module_type::Int=0,
        losses::Real=0.14,
        azimuth::Real=180,
        name::String="PV",
        location::String="both",
        existing_kw::Real=0,
        min_kw::Real=0,
        max_kw::Real=1.0e5,
        cost_per_kw::Real=1600.0,
        om_cost_per_kw::Real=16.0,
        degradation_pct::Real=0.005,
        macrs_option_years::Int = 5,
        macrs_bonus_pct::Float64 = 1.0,
        macrs_itc_reduction::Float64 = 0.5,
        total_itc_pct::Float64 = 0.26,
        total_rebate_per_kw::Float64 = 0.0,
        kw_per_square_foot::Float64=0.01,
        acres_per_kw::Float64=6e-3,
        inv_eff::Float64=0.96,
        dc_ac_ratio::Float64=1.2,
        prod_factor_series_kw::Union{Missing, Array{Real,1}} = missing
        )

        # validate inputs
        invalid_args = String[]
        if !(0 <= azimuth < 360)
            push!(invalid_args, "azimuth must satisfy 0 <= azimuth < 360, got $(azimuth)")
        end
        if !(array_type in [0, 1, 2, 3, 4])
            push!(invalid_args, "array_type must be in [0, 1, 2, 3, 4], got $(array_type)")
        end
        if !(module_type in [0, 1, 2])
            push!(invalid_args, "module_type must be in [0, 1, 2], got $(module_type)")
        end
        if !(0.0 <= losses <= 0.99)
            push!(invalid_args, "losses must satisfy 0.0 <= losses <= 0.99, got $(losses)")
        end
        if !(0 <= tilt <= 90)
            push!(invalid_args, "tilt must satisfy 0 <= tilt <= 90, got $(tilt)")
        end
        if !(location in ["roof", "ground", "both"])
            push!(invalid_args, "location must be in [\"roof\", \"ground\", \"both\"], got $(location)")
        end
        if !(0.0 <= degradation_pct <= 1.0)
            push!(invalid_args, "degradation_pct must satisfy 0 <= degradation_pct <= 1, got $(degradation_pct)")
        end
        if !(0.0 <= inv_eff <= 1.0)
            push!(invalid_args, "inv_eff must satisfy 0 <= inv_eff <= 1, got $(inv_eff)")
        end
        if !(0.0 <= dc_ac_ratio <= 2.0)
            push!(invalid_args, "dc_ac_ratio must satisfy 0 <= dc_ac_ratio <= 1, got $(dc_ac_ratio)")
        end
        # TODO validate additional args
        if length(invalid_args) > 0
            error("Invalid argument values: $(invalid_args)")
        end

        new(
            tilt,
            array_type,
            module_type,
            losses,
            azimuth,
            name,
            location,
            existing_kw,
            min_kw,
            max_kw,
            cost_per_kw,
            om_cost_per_kw,
            degradation_pct,
            macrs_option_years,
            macrs_bonus_pct,
            macrs_itc_reduction,
            total_itc_pct,
            total_rebate_per_kw,
            kw_per_square_foot,
            acres_per_kw,
            inv_eff,
            dc_ac_ratio,
            prod_factor_series_kw
        )
    end
end

ext: jl

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
module RDExample

using ForwardDiff

greet() = print("Hello World!")

include("extrafile_RD.jl")

export my_fun, der_my_fun

end # module


text:
#####
##### `sum(x)`
#####

function frule((_, ẋ), ::typeof(sum), x::Tuple)
    return sum(x), sum(ẋ)
end
function frule((_, ẋ), ::typeof(sum), x; dims=:)
    return sum(x; dims=dims), sum(ẋ; dims=dims)
end

function rrule(::typeof(sum), x::AbstractArray; dims=:)
    project = ProjectTo(x)
    y = sum(x; dims=dims)
    function sum_pullback(dy_raw)
        dy = unthunk(dy_raw)
        x_thunk = InplaceableThunk(
            # Protect `dy` from broadcasting, for when `x` is an array of arrays:
            dx -> dx .+= (dims isa Colon ? Ref(dy) : dy),
            @thunk project(_unsum(x, dy, dims))  # `_unsum` handles Ref internally
        )
        return (NoTangent(), x_thunk)
    end
    return y, sum_pullback
end

# This broadcasts `dy` to the shape of `x`, and should preserve e.g. CuArrays, StaticArrays.
# Ideally this would only need `typeof(x)` not `x`, but `similar` only has a suitable method
# when `eltype(x) == eltype(dy)`, which isn't guaranteed.
_unsum(x, dy, dims) = broadcast(last∘tuple, x, dy)
_unsum(x, dy, ::Colon) = broadcast(last∘tuple, x, Ref(dy))

# Allow for second derivatives of `sum`, by writing rules for `_unsum`:

function frule((_, _, dydot, _), ::typeof(_unsum), x, dy, dims)
    return _unsum(x, dy, dims), _unsum(x, dydot, dims)
end

function rrule(::typeof(_unsum), x, dy, dims)
    z = _unsum(x, dy, dims)
    _unsum_pullback(dz) = (NoTangent(), NoTangent(), sum(unthunk(dz); dims=dims), NoTangent())
    return z, _unsum_pullback
end

#####
##### `sum(f, x)`
#####

# Can't map over Adjoint/Transpose Vector
function rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(sum),
    f,
    xs::Union{Adjoint{<:Number,<:AbstractVector},Transpose{<:Number,<:AbstractVector}};
    kwargs...
)
    op = xs isa Adjoint ? adjoint : transpose
    # since summing a vector we don't need to worry about dims which simplifies adjointing
    vector = parent(xs)
    y, vector_sum_pb = rrule(config, sum, f, vector; kwargs...)
    function covector_sum_pb(ȳ)
        s̄um, f̄, v̄ = vector_sum_pb(ȳ)
        return s̄um, f̄, op(v̄)
    end

    return y, covector_sum_pb
end

function rrule(
    config::RuleConfig{>:HasReverseMode}, ::typeof(sum), f, xs::AbstractArray; dims=:
)
    fx_and_pullbacks = map(x->rrule_via_ad(config, f, x), xs)
    y = sum(first, fx_and_pullbacks; dims=dims)

    pullbacks = last.(fx_and_pullbacks)

    project = ProjectTo(xs)

    function sum_pullback(ȳ)
        call(f, x) = f(x)
        # if dims is :, then need only left-handed only broadcast
        broadcast_ȳ = dims isa Colon  ? (ȳ,) : ȳ
        f̄_and_x̄s = call.(pullbacks, broadcast_ȳ)
        # no point thunking as most of work is in f̄_and_x̄s which we need to compute for both
        f̄ = if fieldcount(typeof(f)) === 0 # Then don't need to worry about derivative wrt f
            NoTangent()
        else
            sum(first, f̄_and_x̄s)
        end
        x̄s = map(unthunk ∘ last, f̄_and_x̄s) # project does not support receiving InplaceableThunks
        return NoTangent(), f̄, project(x̄s)
    end
    return y, sum_pullback
end

# https://github.com/JuliaDiff/ChainRules.jl/issues/522
# The rule above assumes `f` is callable. Arrays are not, this came up when summing
# arrays with weights in StatsBase
@opt_out ChainRulesCore.rrule(
    config::RuleConfig{>:HasReverseMode},
    ::typeof(sum),
    x::AbstractArray,
    y::AbstractArray;
    dims=:
)

function frule(
    (_, _, Δx),
    ::typeof(sum),
    ::typeof(abs2),
    x::AbstractArray{T};
    dims=:,
) where {T<:Union{Real,Complex}}
    ẋ = unthunk(Δx)
    y = sum(abs2, x; dims=dims)
    ∂y = if dims isa Colon
        2 * realdot(x, ẋ)
    elseif VERSION ≥ v"1.2" # multi-iterator mapreduce introduced in v1.2
        mapreduce(+, x, ẋ; dims=dims) do xi, dxi
            2 * realdot(xi, dxi)
        end
    else
        2 * sum(realdot.(x, ẋ); dims=dims)
    end
    return y, ∂y
end

function rrule(
    ::typeof(sum),
    ::typeof(abs2),
    x::AbstractArray{T};
    dims=:,
) where {T<:Union{Real,Complex}}
    y = sum(abs2, x; dims=dims)
    function sum_abs2_pullback(ȳ)
        x_thunk = InplaceableThunk(
            dx -> dx .+= 2 .* real.(ȳ) .* x,
            @thunk(2 .* real.(ȳ) .* x),
        )
        return (NoTangent(), NoTangent(), x_thunk)
    end
    return y, sum_abs2_pullback
end

# Fix dispatch for this pidgeon-hole optimization,
# Rules with RuleConfig dispatch with priority over without (regardless of other args).
# and if we don't specify what do do for one that HasReverseMode then it is ambigious
for Config in (RuleConfig, RuleConfig{>:HasReverseMode})
    @eval function rrule(
        ::$Config, ::typeof(sum), ::typeof(abs2), x::AbstractArray{T}; dims=:,
    ) where {T<:Union{Real,Complex}}
        return rrule(sum, abs2, x; dims=dims)
    end
end

#####
##### `prod`
#####

function rrule(::typeof(prod), x::AbstractArray{T}; dims=:) where {T<:CommutativeMulNumber}
    y = prod(x; dims=dims)
    project_x = ProjectTo(x)
    # vald = dims isa Colon ? nothing : dims isa Integer ? Val(Int(dims)) : Val(Tuple(dims))
    function prod_pullback(ȳ)
        dy = unthunk(ȳ)
        x_thunk = InplaceableThunk(
            # In-place versions -- same branching
            dx -> if dims === (:)
                ∇prod!(dx, x, dy, y)
            elseif any(iszero, x)
                vald = dims isa Colon ? nothing : dims isa Integer ? Val(Int(dims)) : Val(Tuple(dims))
                ∇prod_dims!(dx, vald, x, dy, y)
            else
                dx .+= conj.(y ./ x) .* dy
            end,
            # Out-of-place versions
            @thunk project_x(if dims === (:)
                ∇prod(x, dy, y)
            elseif any(iszero, x)  # Then, and only then, will ./x lead to NaN
                vald = dims isa Colon ? nothing : dims isa Integer ? Val(Int(dims)) : Val(Tuple(dims))
                ∇prod_dims(vald, x, dy, y)  # val(Int(dims)) is about 2x faster than Val(Tuple(dims))
            else
                conj.(y ./ x) .* dy
            end)
        )
        return (NoTangent(), x_thunk)
    end
    return y, prod_pullback
end

function ∇prod_dims(vald::Val{dims}, x, dy, y=prod(x; dims=dims)) where {dims}
    T = promote_type(eltype(x), eltype(dy))
    dx = fill!(similar(x, T, axes(x)), zero(T))
    ∇prod_dims!(dx, vald, x, dy, y)
    return dx
end

function ∇prod_dims!(dx, ::Val{dims}, x, dy, y) where {dims}
    iters = ntuple(d -> d in dims ? tuple(:) : axes(x,d), ndims(x))  # Without Val(dims) this is a serious type instability
    @inbounds for ind in Iterators.product(iters...)
        jay = map(i -> i isa Colon ? 1 : i, ind)
        @views ∇prod!(dx[ind...], x[ind...], dy[jay...], y[jay...])
    end
    return dx
end

function ∇prod(x, dy::Number=1, y::Number=prod(x))
    T = promote_type(eltype(x), eltype(dy))
    dx = fill!(similar(x, T, axes(x)), zero(T)) # axes(x) makes MArray on StaticArrays, Array for structured matrices
    ∇prod!(dx, x, dy, y)
    return dx
end

function ∇prod!(dx, x, dy::Number=1, y::Number=prod(x))
    numzero = iszero(y) ? count(iszero, x) : 0
    if numzero == 0  # This can happen while y==0, if there are several small xs
        dx .+= conj.(y ./ x) .* dy
    elseif numzero == 1
        ∇prod_one_zero!(dx, x, dy)
    else
        # numzero > 1, then all first derivatives are zero
    end
    return dx
end

function ∇prod_one_zero!(dx, x, dy::Number=1)  # Assumes exactly one x is zero
    i_zero = 0
    p_rest = one(promote_type(eltype(x), typeof(dy)))
    for i in eachindex(x)
        xi = @inbounds x[i]
        p_rest *= ifelse(iszero(xi), one(xi), conj(xi))
        i_zero = ifelse(iszero(xi), i, i_zero)
    end
    dx[i_zero] += p_rest * dy
    return
end

#####
##### `cumprod`
#####

function rrule(::typeof(cumprod), x::AbstractVector{<:Real}; dims::Integer=1)
    y = cumprod(x; dims=dims)  # does nothing unless dims == 1
    project_x = ProjectTo(x)
    function cumprod_pullback_1(dy_raw)
        dy = unthunk(dy_raw)
        dx_thunk = InplaceableThunk(
            dx -> if dims == 1
                ∇cumprod!(dx, x, dy, y)
            else
                dx .+= dy
            end
            ,
            @thunk project_x(if dims == 1
                ∇cumprod(x, dy, y)
            else
                dy
            end)
            )
        return (NoTangent(), dx_thunk)
    end
    return y, cumprod_pullback_1
end

function rrule(::typeof(cumprod), x::AbstractArray{<:Real}; dims::Integer)
    y = cumprod(x; dims=dims)
    project_x = ProjectTo(x)
    function cumprod_pullback_2(dy_raw)
        dy = unthunk(dy_raw)
        dx_thunk = InplaceableThunk(
            dx -> if dims <= ndims(x)
                vald = Val(Int(dims))
                ∇cumprod_dim!(dx, vald, x, dy, y)
            else
                dx .+= dy
            end
            ,
            @thunk project_x(if dims <= ndims(x)
                vald = Val(Int(dims))
                ∇cumprod_dim(vald, x, dy, y)
            else
                dy
            end)
            )
        return (NoTangent(), dx_thunk)
    end
    return y, cumprod_pullback_2
end

function ∇cumprod_dim(vald::Val{dim}, x::AbstractArray, dy=fill!(zero(x),1), y=cumprod(x; dims=dim)) where {dim}
     T = promote_type(eltype(x), eltype(dy))
     dx = fill!(similar(x, T, axes(x)), zero(T))
     ∇cumprod_dim!(dx, vald, x, dy, y)
     return dx
 end

@inline function ∇cumprod_dim!(dx::AbstractArray, ::Val{dim}, x::AbstractArray, dy, y) where {dim}
    iters = ntuple(k -> k==dim ? Ref(:) : axes(x,k), ndims(x))
    for ind in Iterators.product(iters...)
        @views ∇cumprod!(dx[ind...], x[ind...], dy[ind...], y[ind...])
    end
    return dx
end

function ∇cumprod(x::AbstractVector, dy=one(x), y=cumprod(x))
    T = promote_type(eltype(x), eltype(dy))  # really needs to allow dy * y / x
    dx = fill!(similar(x, T, axes(x)), zero(T))  # axes(x) makes MArray on StaticArrays, Array for structured matrices
    ∇cumprod!(dx, x, dy, y)
    return dx
end

@inline function ∇cumprod!(dx::AbstractVector, x::AbstractVector, dy, y)
    lo, hi = firstindex(x), lastindex(x)
    z = something(findfirst(iszero, x), hi+1)
    acc = zero(eltype(dy))
    @inbounds for k in z-1:-1:lo
        acc += y[k] * dy[k]
        dx[k] += acc / x[k]
    end
    @inbounds if z != hi+1
        yk = z==1 ? one(eltype(y)) : y[z-1]  # will be prod(x[j] for j=1:k if j!=z)
        dx[z] += yk * dy[z]
        for k in (z+1):hi
            yk *= x[k]
            dx[z] += yk * dy[k]
        end
    end
    return dx
end

#####
##### `foldl`
#####

# `foldl` guarantees to execute `f` in order, left to right. So it makes sense even when
# this `f` is stateful, in which case the gradient must be calculated in the reverse order.

# The implementation aims to be efficient for both tuples and arrays, although using accumulate
# to carry intermediate results along creates arrays of tuples which could be avoided; using a
# loop can be a few times faster. Note also that it does not return a gradient for `init`.

function rrule(
        config::RuleConfig{>:HasReverseMode}, ::typeof(foldl), op::G, x::Union{AbstractArray, Tuple};
        init=_InitialValue()
    ) where {G}
    list, start = if init === _InitialValue()
        _drop1(x), first(x)
    else
        # Case with init keyword is simpler to understand first!
        _reshape1(x, :), init  # (vec is for Julia 1.0, accumulate is fussy)
    end
    hobbits = accumulate(list; init=(start, nothing)) do (a,_), b
        # Here `a` is what we would normally cary forward, and `_` ignores
        # the previous iteration's pullback function (needed later),
        # while `b` is the fresh input from `list` as usual.
        c, back = rrule_via_ad(config, op, a, b)  # LHS is just documentation here!
        # We don't really need to store every `c`, last one is `foldl` output.
        # (The name, BTW, is because "there and back again" is the subtitle of Tolkien's book.)
    end
    y = first(last(hobbits))
    axe = axes(x)
    project = ProjectTo(x)
    function unfoldl(dy)
        trio = accumulate(_reverse1(hobbits); init=(0, dy, 0)) do (_, dc, _), (_, back)
            ds, da, db = back(dc)
            # Don't need to store every `da`, need one for the next iteration + maybe last
        end
        dop = sum(first, trio)
        dx = map(last, _reverse1(trio))
        if init === _InitialValue()
            # `hobbits` is one short
            dx = _vcat1(trio[end][2], dx)
        end
        return (NoTangent(), dop, project(_reshape1(dx, axe)))
    end
    return y, unfoldl
end

#####
##### Iterator-or-Tuple functions
#####

# This zoo of underscore functions helps `foldl` & `accumulate` handle both tuples and arrays,
# and also provides some alternatives for versions of Julia where iterators weren't supported.
# Inspired by `Base._reverse`, used in defn of `foldr`.

# To support 2nd derivatives, some may need their own gradient rules. And _drop1 should perhaps
# be replaced by _peel1 like Iterators.peel

if VERSION >= v"1.6"
    _reverse1(x) = Iterators.reverse(x)
    _drop1(x) = Iterators.drop(x, 1)
    _zip2(x, y) = zip(x, y)  # for `accumulate`, below
else
    # Old versions don't support accumulate(::itr), nor multi-dim reverse
    _reverse1(x) = reverse(vec(x))
    _drop1(x) = vec(x)[2:end]
    _zip2(x, y) = collect(zip(x, y))
end
_reverse1(x::Tuple) = reverse(x)
_drop1(x::Tuple) = Base.tail(x)
_zip2(x::Tuple{Vararg{Any,N}}, y::Tuple{Vararg{Any,N}}) where N = ntuple(i -> (x[i],y[i]), N)

struct _InitialValue end  # Old versions don't have `Base._InitialValue`

_vcat1(x, ys::AbstractVector) = vcat(x, ys)
_vcat1(x::AbstractArray, ys::AbstractVector) = vcat([x], ys)
_vcat1(x, ys::Tuple) = (x, ys...)

_reshape1(x::AbstractArray, axe) = reshape(x, axe)
_reshape1(x::Tuple, axe) = x

_no_tuple_tangent(dx::Tangent) = ChainRulesCore.backing(dx)
_no_tuple_tangent(dx) = dx

#####
##### `accumulate`
#####

# Like `foldl` this by definition works in order, so it makes sense to allow stateful `f`.

function rrule(
        config::RuleConfig{>:HasReverseMode}, ::typeof(accumulate), op::G, x::Union{AbstractArray, Tuple};
        init=_InitialValue(), dims=nothing
    ) where {G}
    isnothing(dims) || dims == 1 && x isa Base.AbstractVecOrTuple || throw(
        "accumulate(op, x; dims) is not currently supported by ChainRules, sorry"
        # It's not supported by AD either, so no point calling back, and no regression:
        # gradient(x -> sum(accumulate(/, x, dims=1)), rand(3,4))
        # ERROR: Mutating arrays is not supported
    )
    list, start = if init === _InitialValue()
        _drop1(x), first(x)
    else
        x, init
    end
    hobbits = accumulate(list; init = (start, nothing)) do (a, _), b
        c, back = rrule_via_ad(config, op, a, b)
    end
    y = map(first, hobbits)
    if init === _InitialValue()
        # `hobbits` is one short, and first one doesn't invoke `op`
        y = _vcat1(first(x), y)
    end
    axe = axes(x)
    project = ProjectTo(x)
    function decumulate(dy)
        dy_plain = _no_tuple_tangent(unthunk(dy))
        rev_list = if init === _InitialValue()
            if VERSION >= v"1.6"
                # Here we rely on `zip` to stop early. Begin explicit with _reverse1(_drop1(...))
                # gets "no method matching iterate(::Base.Iterators.Reverse{Base.Iterators.Drop{Array{"
                _zip2(_reverse1(hobbits), _reverse1(dy_plain))
            else
                # However, on 1.0 and some others, zip does not stop early. But since accumulate
                # also doesn't work on iterators, `_drop1` doesn't make one, so this should work:
                _zip2(_reverse1(hobbits), _reverse1(_drop1(dy_plain)))
                # What an awful tangle.
            end
        else
            _zip2(_reverse1(hobbits), _reverse1(dy_plain))
        end
        trio = accumulate(rev_list; init=(0, ZeroTangent(), 0)) do (_, dc, _), ((_, back), dz)
            ds, da, db = back(dc + dz)
            # Don't need to store every 'da', but need for next iteration, and the last one.
        end
        dop = sum(first, trio)
        dx = map(last, _reverse1(trio))
        if init == _InitialValue()
            # `hobbits` is one short, and the first one is weird
            dx = _vcat1(trio[end][2] + dy_plain[1], dx)
        end
        return (NoTangent(), dop, project(_reshape1(dx, axe)))
    end
    return _reshape1(y, axe), decumulate
end


text:
# Computes the fraction of mutations that do not change the output
# If active_only==true, then only mutations that change active nodes are considered
# If active_only==false, then all mutations are considered
# Note that results are partially random since the result of mutations at a specific location can be different.
function mutational_robustness( c::Chromosome, funcs::Vector{Func}; active_only::Bool=false )
  context = construct_context( c.params.numinputs )
  prev_out = execute_chromosome( c, context )
  count_no_change = 0
  count_no_change = 0.0
  num_mut_locs = 0
  for j = 1:num_mutate_locations( c, funcs )
    cc = deepcopy(c)
    (cc,active) = mutate_chromosome!( cc, funcs, j )
    out_c = execute_chromosome( cc, context )
    #println("j: ",j,"  active: ",active)
    #print_chromosome( cc )
    if !active_only || active # consider all mutations if !active_only, only active mutations if active_only
      #out_c = execute_chromosome( cc, context )
      count_no_change += ((out_c == prev_out) ? 1 : 0)
      if active_only && active
        num_mut_locs += 1
      end
    end
    println("j: ",j,"  active: ",active,"  no_change: ",(out_c == prev_out),"  count_no_change: ",count_no_change)
  end
  println("num_mut_locs: ",num_mut_locs,"  num_mutate_locations(): ",num_mutate_locations(c,funcs),"  count_no_change: ",count_no_change)
  #@assert num_mut_locs == num_mutate_locations(c,funcs)
  if !active_only
    num_mut_locs = num_mutate_locations( c, funcs )
  end
  count_no_change/num_mut_locs
end


text:
function MyFunc(x)
    return x * 1.23456
end

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 67d98a4a-4da0-11eb-3ff8-db5e68b54f8d
begin
    import Pkg
    #Pkg.activate("Mage")

    md"""
    ## BASIC MtG GAME SETTINGS

    #### Edit the elements below to modify basic game options.
    """
end

# ╔═╡ c656c348-3852-11eb-2806-0df5baa0d2e1
begin
    psg = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatSimulator"
    gd = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage"

    import Pkg
    Pkg.activate(psg)

    using Colors
    using PlutoUI
    using Serialization
    using PlaymatSimulator
end

# ╔═╡ 52c3d32e-4e2f-11eb-20f7-4d61ff865e75
US = deserialize("$gd/../../PlaymatSimulator/tmp/user_selection.jls")

# ╔═╡ 3ffcf8a4-385b-11eb-25da-fd9bfc189e06
# length in pixels
const SCREEN_WIDTH = Int32(1920)

# ╔═╡ 519b147e-385b-11eb-30ef-9354da78c3b2
# length in pixels
const SCREEN_HEIGHT = Int32(1080)

# ╔═╡ 731afc68-385b-11eb-1084-1f85f3b05ed0
# length in pixels
const SCREEN_BORDER = Int32(10)

# ╔═╡ a913a358-3930-11eb-104b-a1e5b1b9091c
# length in pixels
const DEFAULT_FONT = "Base/fonts/OpenSans-Regular.ttf"

# ╔═╡ 2936eba6-3942-11eb-14ec-05cab84f472d
# max FPS
const MAX_FPS = Int32(60)

# ╔═╡ 3fe53d58-4d9c-11eb-11c3-1f1ea0505f98
begin
    BKG_NAMES = [f => f for f in readdir("$gd/Base/ui/backgrounds")]

    md"""
    `BKG_NAMES:`
    $(@bind BKG_NAME Select(BKG_NAMES,  default="land.jpg"))
    """
end

# ╔═╡ 089d469c-4da1-11eb-2c82-6132328b548c
BKG_PATH = "$gd/Base/ui/backgrounds/$BKG_NAME"

# ╔═╡ 2b99d0d8-4d9d-11eb-3f90-2fd6c50a9fa8
BACKGROUND_IMG = LocalResource(BKG_PATH)

# ╔═╡ 0c14a126-3852-11eb-2561-9dd3f0f435b0
md"""
`BACKGROUND:`
$(BACKGROUND = colorant"black")

*Available colornames can be found here: http://juliagraphics.github.io/Colors.jl/stable/namedcolors/*
"""

# ╔═╡ 1fd3ee0a-4faa-11eb-1ee9-a5b89af3d078
html"""

"""

# ╔═╡ e3105084-4c1e-11eb-297d-e1654170ec4a
GS = Dict(
    :SCREEN_WIDTH => SCREEN_WIDTH,
    :SCREEN_HEIGHT => SCREEN_HEIGHT,
    :SCREEN_BORDER => SCREEN_BORDER,
    :BACKGROUND => BACKGROUND,
    :DEFAULT_FONT => DEFAULT_FONT,
    :MAX_FPS => MAX_FPS,
);

# ╔═╡ Cell order:
# ╟─67d98a4a-4da0-11eb-3ff8-db5e68b54f8d
# ╟─c656c348-3852-11eb-2806-0df5baa0d2e1
# ╟─52c3d32e-4e2f-11eb-20f7-4d61ff865e75
# ╟─3ffcf8a4-385b-11eb-25da-fd9bfc189e06
# ╟─519b147e-385b-11eb-30ef-9354da78c3b2
# ╟─731afc68-385b-11eb-1084-1f85f3b05ed0
# ╟─a913a358-3930-11eb-104b-a1e5b1b9091c
# ╟─2936eba6-3942-11eb-14ec-05cab84f472d
# ╟─3fe53d58-4d9c-11eb-11c3-1f1ea0505f98
# ╟─089d469c-4da1-11eb-2c82-6132328b548c
# ╟─2b99d0d8-4d9d-11eb-3f90-2fd6c50a9fa8
# ╟─0c14a126-3852-11eb-2561-9dd3f0f435b0
# ╟─1fd3ee0a-4faa-11eb-1ee9-a5b89af3d078
# ╟─e3105084-4c1e-11eb-297d-e1654170ec4a


text:
@compat function (meshtype::Type{T}){T <: AbstractMesh}(c::Pyramid)
    T(decompose(vertextype(T), c), decompose(facetype(T), c))
end

"""
Standard way of creating a mesh from a primitive.
Just walk through all attributes of the mesh and try to decompose it.
If there are attributes missing, just hope it will get managed by the mesh constructor.
(E.g. normal calculation, which needs to have vertices and faces present)
"""
@compat function (meshtype::Type{T}){T <: AbstractMesh}(c::GeometryPrimitive, args...)
    attribs = attributes(T)
    newattribs = Dict{Symbol, Any}()
    for (fieldname, typ) in attribs
        if isdecomposable(eltype(typ), c)
            newattribs[fieldname] = decompose(eltype(typ), c, args...)
        end
    end
    T(homogenousmesh(newattribs))
end

@compat function (meshtype::Type{T}){T <: HMesh,HT}(
        c::Union{HyperCube{3,T}, HyperRectangle{3,HT}}

)    xdir = Vec{3, HT}(widths(c)[1],0f0,0f0)
    ydir = Vec{3, HT}(0f0,widths(c)[2],0f0)
    zdir = Vec{3, HT}(0f0,0f0,widths(c)[3])
    quads = [
        Quad(origin(c) + zdir,   xdir, ydir), # Top
        Quad(origin(c),          ydir, xdir), # Bottom
        Quad(origin(c) + xdir,   ydir, zdir), # Right
        Quad(origin(c),          zdir, ydir), # Left
        Quad(origin(c),          xdir, zdir), # Back
        Quad(origin(c) + ydir,   zdir, xdir) #Front
    ]
    merge(map(meshtype, quads))
end


text:
# Define the AbstractGP type and its API.

"""
    abstract type AbstractGP end

Supertype for various Gaussian process (GP) types. A common interface is provided for
interacting with each of these objects. See [1] for an overview of GPs.

[1] - C. E. Rasmussen and C. Williams. "Gaussian processes for machine learning".
MIT Press. 2006.
"""
abstract type AbstractGP end

"""
    mean(f::AbstractGP, x::AbstractVector)

Computes the mean vector of the multivariate Normal `f(x)`.
"""
Statistics.mean(::AbstractGP, ::AbstractVector)

"""
    cov(f::AbstractGP, x::AbstractVector)

Compute the `length(x)` by `length(x)` covariance matrix of the multivariate Normal `f(x)`.
"""
Statistics.cov(::AbstractGP, x::AbstractVector)

"""
    var(f::AbstractGP, x::AbstractVector)

Compute only the diagonal elements of `cov(f(x))`.
"""
Statistics.var(::AbstractGP, ::AbstractVector)

"""
    cov(f::AbstractGP, x::AbstractVector, y::AbstractVector)

Compute the `length(x)` by `length(y)` cross-covariance matrix between `f(x)` and `f(y)`.
"""
Statistics.cov(::AbstractGP, x::AbstractVector, y::AbstractVector)

"""
    mean_and_cov(f::AbstractGP, x::AbstractVector)

Compute both `mean(f(x))` and `cov(f(x))`. Sometimes more efficient than separately
computation, particularly for posteriors.
"""
StatsBase.mean_and_cov(f::AbstractGP, x::AbstractVector) = (mean(f, x), cov(f, x))

"""
    mean_and_var(f::AbstractGP, x::AbstractVector)

Compute both `mean(f(x))` and the diagonal elements of `cov(f(x))`. Sometimes more efficient
than separately computation, particularly for posteriors.
"""
StatsBase.mean_and_var(f::AbstractGP, x::AbstractVector) = (mean(f, x), var(f, x))

for (m, f) in [
    (:Statistics, :mean),
    (:Statistics, :var),
    (:Statistics, :cov),
    (:StatsBase, :mean_and_cov),
    (:StatsBase, :mean_and_var),
]
    @eval function $m.$f(::AbstractGP)
        return error(
            "`",
            $f,
            "(f::AbstractGP)` is not defined (on purpose!).\n",
            "Please provide an `AbstractVector` of locations `x` at which you wish to compute your ",
            $f,
            $((f === :mean_and_cov || f === :mean_and_var) ? " vectors" : " vector"),
            ", and call `",
            $f,
            "(f(x))`\n",
            "For more details please have a look at the AbstractGPs docs.",
        )
    end
end


text:
using PkgBenchmark
import GibbsTypePriors

results = benchmarkpkg("/home/gkonkamking/Work/GibbsTypePriors", retune = true)
PkgBenchmark.export_markdown("benchmark_results.md", results)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
"""
    SeisWindow(in,out;<keyword arguments>)

Window a seis file using header words.

# Arguments
* `in::String`: filename of input
* `out::String`: filename of output

# Keyword arguments
* `key`
* `minval`
* `maxval`

note that windowing along the time axis is achieved by using the key "t".

*Credits: AS, 2015*
"""
function SeisWindow(in::String,out::String;key=[],minval=[],maxval=[])
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    extent = ReadTextHeader(in)
    tmin = extent.o1
    tmax = extent.o1 + extent.d1*(extent.n1-1)
    for ikey=1:length(key)
	if key[ikey] == "t"
	    tmin = minval[ikey]
	    tmax = maxval[ikey]
	    println("tmin= ",tmin," tmax= ",tmax)
	end
    end

#itmin is the sample at which tmin is saved
    itmin = convert(Int32,round((tmin - extent.o1)/extent.d1) + 1)
    if (itmin < 1)
	     itmin = 1
    end
    itmax = convert(Int32,round((tmax - extent.o1)/extent.d1))
    if (itmax) > extent.n1
	     itmax = extent.n1
    end
println("itmin= ",itmin," itmax= ",itmax)

    SeisWindowHeaders(in, out; key=key, minval=minval, maxval=maxval, tmin=tmin,
                      tmax=tmax)
    FetchTraces(in,out;itmin=itmin,itmax=itmax)
    tmp = join(["tmp_SeisWindow_",string(round(Int,rand()*100000))])
    SeisProcessHeaders(out, tmp, [UpdateHeader],
                               [Dict(:itmin=>itmin,:itmax=>itmax)])
    filename_h_tmp = join([DATAPATH tmp "@headers@"])
    filename_h_out = join([DATAPATH out "@headers@"])
    cp(filename_h_tmp,filename_h_out,force=true);
    rm(filename_h_tmp,force=true);

end

function FetchTraces(in::String, out::String; ntrace=500, itmin=round(Int,1),
                     itmax=round(Int,9999999999))
    NX = GetNumTraces(out)
	println("NX in FetchTraces= ",NX)
    itrace = round(Int,1)
    filename_data_in = ParseDataName(in)
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_data_out = join([DATAPATH out "@data@"])
    stream_in  = open(filename_data_in)
    stream_out  = open(filename_data_out,"w")
    extent = ReadTextHeader(in)
    nt = extent.n1

    if itmin < 1
	itmin = round(Int,1)
    end
    if itmax > nt
	itmax = round(Int,nt)
    end
    while itrace <= NX
	if (itrace > 1)
	    stream_out = open(filename_data_out,"a")
	end
	nx = NX - itrace + 1
	ntrace = nx > ntrace ? ntrace : nx
	h = SeisReadHeaders(out,group="some",itrace=itrace,ntrace=ntrace)
	d = zeros(Float32,itmax-itmin+1,ntrace)
	SeekTraces!(d,stream_in,h,itmin,itmax,nt,round(Int,ntrace))

	write(stream_out,d)
	close(stream_out)
	itrace += ntrace
    end
    close(stream_in)
end

function SeekTraces!(d::AbstractArray{T,2}, stream_in::IOStream,
                     h::Array{Header,1},itmin,itmax,nt,ntrace) where T
    d1 = zeros(Float32,nt)
    for ix = 1 : ntrace
	position = 4*nt*(h[ix].tracenum-1)
	seek(stream_in,position)
	d1 = read!(stream_in,Array{Float32}(undef,nt))
	for it = itmin : itmax
	    d[it - itmin + 1,ix] = d1[it]
	end
    end
    nothing
end

function UpdateHeader(h;itrace=1,itmin=1,itmax=9e9)
    ot = (itmin-1)*h[1].d1 + h[1].o1
    nt = itmax - itmin + 1
    for j = 1 : length(h)
	h[j].o1 = ot
	h[j].n1 = nt
	h[j].tracenum = j + itrace - 1
    end
    return h
end


text:
#
# Proof of concept using real data.
# Test with the real data set
#
include("./lib/metadata.jl")
include("./lib/pccf.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")

using SignalGenerator
using BayesChpOrig
using Metadata
using Roc
using Pccf

function concept_proof_real(datafolder, outputfolder)
    println("... concept_proof_real > data folder   = ", datafolder)
    println("... concept_proof_real > output folder = ", outputfolder)

    # work_path = "./data/activity-recognition-sample/"
    signal = readdlm(datafolder * "sample-sig-10.dat", ',')
    x      = signal[:,1]
    y      = signal[:,2]

    n   = length(x)
    mu0 = 84.6
    sd0 = 9.16

    # l_common = 300
    l_common = 100

    maxes, mat = bayes_detector(x, lambda = l_common)
    detections = extract_changes(maxes)

    probs                = scale01(pccf(n, mu0, sd0))
    pccf_maxes, pccf_mat = pccf_bayes_detector(x, probs, l_common + 100, l_common, 0.3, h_common = l_common, max_d = 500)
    pccf_detections      = extract_changes(pccf_maxes)

    changes = readdlm(projectdir * "src_julia/data/activity-recognition-sample/changes-sig-active-10.dat")
    changes = reshape(changes, length(changes))
    if 1==1
        mid_changes = readdlm(projectdir * "src_julia/data/activity-recognition-sample/mid-changes-sig-active-10.dat")
        mid_changes = reshape(mid_changes, length(mid_changes))
        changes = vcat(changes, mid_changes)
        sort!(changes)
    end

    max_delay = 10
    tpr       = tpRate(detections, changes, max_delay)
    tpr_pccf  = tpRate(pccf_detections, changes, max_delay)
    fpr       = faRateScaled(detections     , changes, max_delay, n, tn_price_scale = 0.1 * 0.5)
    fpr_pccf  = faRateScaled(pccf_detections, changes, max_delay, n, tn_price_scale = 0.1 * 0.5)

    @printf( "... concept_proof4 > tp      = %0.3f; fp      = %0.3f\n", tpr, fpr)
    @printf( "... concept_proof4 > tp_pccf = %0.3f; fp_pccf = %0.3f\n", tpr_pccf, fpr_pccf)

    writedlm(outputfolder * "detections4.dat"     , detections)
    writedlm(outputfolder * "pccf_detections4.dat", pccf_detections)
    writedlm(outputfolder * "pccf4.dat"           , probs)
end

concept_proof_real(projectdir * "src_julia/data/activity-recognition-sample/" , outdir * "concept_proof_real4/")

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency