text:
# see web page for full listing
# https://code.google.com/p/yahoo-finance-managed/wiki/enumQuoteProperty

AfterHoursChangeRealtime = "c8"
AnnualizedGain = "g3"
Ask = "a0"
AskRealtime = "b2"
AskSize = "a5"
AverageDailyVolume = "a2"
Bid = "b0"
BidRealtime = "b3"
BidSize = "b6"
BookValuePerShare = "b4"
Change = "c1"
Change_ChangeInPercent = "c0"
ChangeFromFiftydayMovingAverage = "m7"
ChangeFromTwoHundreddayMovingAverage = "m5"
ChangeFromYearHigh = "k4"
ChangeFromYearLow = "j5"
ChangeInPercent = "p2"
ChangeInPercentRealtime = "k2"
ChangeRealtime = "c6"
Commission = "c3"
Currency = "c4"
DaysHigh = "h0"
DaysLow = "g0"
DaysRange = "m0"
DaysRangeRealtime = "m2"
DaysValueChange = "w1"
DaysValueChangeRealtime = "w4"
DividendPayDate = "r1"
TrailingAnnualDividendYield = "d0"
TrailingAnnualDividendYieldInPercent = "y0"
DilutedEPS = "e0"
EBITDA = "j4"
EPSEstimateCurrentYear = "e7"
EPSEstimateNextQuarter = "e9"
EPSEstimateNextYear = "e8"
ExDividendDate = "q0"
FiftydayMovingAverage = "m3"
SharesFloat = "f6"
HighLimit = "l2"
HoldingsGain = "g4"
HoldingsGainPercent = "g1"
HoldingsGainPercentRealtime = "g5"
HoldingsGainRealtime = "g6"
HoldingsValue = "v1"
HoldingsValueRealtime = "v7"
LastTradeDate = "d1"
LastTradePriceOnly = "l1"
LastTradeRealtimeWithTime = "k1"
LastTradeSize = "k3"
LastTradeTime = "t1"
LastTradeWithTime = "l0"
LowLimit = "l3"
MarketCapitalization = "j1"
MarketCapRealtime = "j3"
MoreInfo = "i0"
Name = "n0"
Notes = "n4"
OneyrTargetPrice = "t8"
Open = "o0"
OrderBookRealtime = "i5"
PEGRatio = "r5"
PERatio = "r0"
PERatioRealtime = "r2"
PercentChangeFromFiftydayMovingAverage = "m8"
PercentChangeFromTwoHundreddayMovingAverage = "m6"
ChangeInPercentFromYearHigh = "k5"
PercentChangeFromYearLow = "j6"
PreviousClose = "p0"
PriceBook = "p6"
PriceEPSEstimateCurrentYear = "r6"
PriceEPSEstimateNextYear = "r7"
PricePaid = "p1"
PriceSales = "p5"
Revenue = "s6"
SharesOwned = "s1"
SharesOutstanding = "j2"
ShortRatio = "s7"
StockExchange = "x0"
Symbol = "s0"
TickerTrend = "t7"
TradeDate = "d2"
TradeLinks = "t6"
TradeLinksAdditional = "f0"
TwoHundreddayMovingAverage = "m4"
Volume = "v0"
YearHigh = "k0"
YearLow = "j0"
YearRange = "k2"


text:
export FFNNet, propagate, train!

"""
`type FFNNet`

Type representing a Neural Network.

### Fields
* `layers` (`Vector{FFNNLayer}`): Vector containing each layer of the network
* `weights` (`Vector{Matrix{Float64}}`): Vector containing the weight matrices between layers
* `inputsize` (`Int`): Input size accepted by this network
"""
type FFNNet
    layers::Vector{FFNNLayer}
    weights::Vector{Matrix{Float64}}
    inputsize::Int
end

#############################
#        CONSTRUCTORS       #
#############################
"""
`FFNNet(sizes::Int...)`

Construct a network given the input size and the sizes of each layer. By default, the hidden layers have an bias unit and the output layer don't. One the other hand, all the layers have `tanh` as activation function by default.

### Arguments
* `sizes` (`Int...`): Integers specifying the sizes for the network. The first is the input size and the rest is the size of each network layer, from the first one up to the size of the output layer.

### Returns
A Neural Network (`FFNNet{N,I}`)
"""
function FFNNet(sizes::Int...)
    if length(sizes) < 3
        error("Network must have at least one hidden layer")
    end

    # Create an Array of Neural Network Layers of the right sizes
    # The first size corresponds to the input size of the network
    layers = Array(FFNNLayer, length(sizes) - 1)
    @inbounds for i in 2:length(sizes) - 1
        layers[i - 1] = FFNNLayer(sizes[i])
    end
    layers[end] = FFNNLayer(sizes[end], bias=false) # Last layer without bias

    return FFNNet(layers, sizes[1])
end

"""
`FFNNet(sizes::Int...)`

Construct a network given its layers and its input size.

### Arguments
* `layers` (`Vector{FFNNLayer}`): A vector with all the layers of the network (in order, with the last one being the output layer).
* `inputsize` (`Int`): Integer specifying the input size of the layer.

### Returns
A Neural Network (`FFNNet`)
"""
function FFNNet(layers::Vector{FFNNLayer}, inputsize::Int)
    # Create a vector of weight matrices
    weights = Array(Matrix{Float64}, length(layers))

    # This weight initialization is sensible (Glorot & Bengio 2010)
    # W ~ U(∓ √6/√(sum of sizes of the 2 layers connected by W))
    ɛ(a,b) = sqrt(6)/(sqrt(a + b))

    # The first weight matrix is from the input (including bias)
    #  to the first layer (excluding bias unit)
    eps = ɛ(size(layers[1]), inputsize)
    weights[1] = rand(size(layers[1]), inputsize + 1)*2*eps - eps

    # Matrices from layer i-1 (including bias) to layer i
    @inbounds for i in 2:length(layers)
        eps = ɛ(size(layers[i]), size(layers[i-1]))
        weights[i] = rand(size(layers[i]), size(layers[i-1]) + 1)
    end

    return FFNNet(layers, weights, inputsize)
end

############################
#      BASIC FUNCTIONS     #
############################
length(net::FFNNet) = length(net.layers)
size(net::FFNNet) = (net.inputsize, length(net.layers))

function show(io::IO, net::FFNNet)
    I, L = size(net)
    print(io, L, " Layers Feedforward Neural Network:\n  Input Size: ", I)
    for (i, l) in enumerate(net.layers)
        print(io, "\n  Layer ", i, ": ", size(l), " neurons")
        if l.bias
            print(io, " + bias unit")
        end
        print(io, ", ", string(l.activation))
    end
end

############################
#      OTHER FUNCTIONS     #
############################
"""
`propagate(net::FFNNet, x::Vector{Float64})`

Propagate an input `x` through the network `net` and return the output

### Arguments
* `net` (`FFNNet`): A neural network that will process the input and give the output
* `x` (`Vector{Float64}`): The input vector. This must be of the same size as the input size of `net`.

### Returns
The output of the network (`Vector{Float64}`). This is simply the activation of the last layer of the network after forwardpropagating the input.
"""
function propagate(net::FFNNet, x::Vector{Float64})
    if length(x) != net.inputsize
        error("Network does not support input size $length(x), only $net.inputsize")
    end

    # Insert bias unit on input and update first layer
    update!(net.layers[1], net.weights[1] * vcat([1.0], x))

    # Update all remaining layers
    @inbounds for i in 2:length(net)
        update!(net.layers[i], net.weights[i] * activate(net.layers[i-1]))
    end

    # Return the activation of the last layer
    return activate(net.layers[end])
end

function constructbatches(perm::Vector{Int}, batchsize::Int)
    batches = Vector{Int}[]
    while length(perm) > 0
        # Size of this batch
        sizebatch = min(length(perm), batchsize)
        push!(batches, splice!(perm, 1:sizebatch))
    end
    return batches
end

function backpropagate(net::FFNNet,
                       output::Vector{Float64},
                       target::Vector{Float64},
                       cost::Function)

    L = length(net)
    # Vector storing one delta vector for each layer
    δ = Array(Vector{Float64}, L)

    # Compute δ for the last layer
    #   δ^L = ∂E/∂(last.neurons)
    last = net.layers[L]  # Last layer
    δ[L] = (der(last.activation)(last))' * der(cost)(output, target)

    # Find δ of previous layers, backwards
    @inbounds for l in (L-1):-1:1
       layer = net.layers[l] # Current layer
       W = net.weights[l+1]  # W^(l+1)

       # δ^l = ϕ'(s^l) ⊙ W^(l+1)' δ^(l+1)
       δ[l] = (der(layer.activation)(layer))' * (W'δ[l+1])

       # If there is a bias unit, remove the first
       # element of δ^l, that corresponds to it
       layer.bias && deleteat!(δ[l], 1)
    end

    return δ
end

"""
`train!(net::FFNNet,
        inputs::Vector{Vector{Float64}},
        outputs::Vector{Vector{Float64}};
        α::Real = 0.5,
        η::Real = 0.1,
        epochs::Int = 1,
        batchsize::Int = 1,
        cost::Function = quaderror)`

Train the Neural Network using the examples provided in `inputs` and `outputs`.

### Arguments
* `net` (`FFNNet`): Feedforward Neural Network to be trained.
* `inputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This have K elements, K being the number of examples in this dataset. Each input vector must have I elements, where I is the input size of `net`.
* `outputs` (`Vector{Vector{Float64}}`): Vector containing input examples (each one is a vector). This must also have K elements, K being the number of examples in this dataset. Each output vector must have O elements, where O is the size of the output layer of `net`.

### Keyword Arguments
* `α` (`Real`, 0.5 by default): Learning Rate.
* `η` (`Real`, 0.1 by default): Momentum Rate.
* `epochs` (`Int`, 1 by default): Number of iterations of the learning algorithm on this dataset.
* `batchsize` (`Int`, 1 by default): Size of the batch used by the algorithm (1 is simply the default stochastic gradient descent).
* `cost` (`Function`, `quaderror` by default): Cost function to be minimized by the learning algorithm.
"""
function train!(net::FFNNet,
                inputs::Vector{Vector{Float64}},
                outputs::Vector{Vector{Float64}};
                α::Real = 0.5,               # Learning rate
                η::Real = 0.1,               # Momentum rate
                epochs::Int = 1,             # Iterations through entire data
                batchsize::Int = 1,          # Number of examples on each iteration
                cost::Function = quaderror)  # Cost function

    L = length(net)

    # Initialize gradient matrices
    grad      = Matrix{Float64}[zeros(i) for i in net.weights]
    last_grad = Matrix{Float64}[zeros(i) for i in net.weights]

    for τ in 1:epochs
        # Select training order randomly
        perm = randperm(length(inputs))

        # Create batches of the correct size
        batches = constructbatches(perm, batchsize)
        for batch in batches
            # Compute the gradient using this entire batch
            for ex in batch
                input_ex   = vcat([1.0], inputs[ex])     # Example's input with bias
                output_ex  = outputs[ex]                 # Example's output
                output_net = propagate(net, inputs[ex]) # Network's output

                # Find the δs using the backpropagation of this example
                δ = backpropagate(net, output_net, output_ex, cost)

                # Find gradients
                grad[1] += δ[1] ⊗ input_ex # First layer     ∇E = δ^1 ⊗ input
                for l in 2:L               # Other layers    ∇E = δ^l ⊗ x^(l-1)
                    grad[l] += δ[l] ⊗ activate(net.layers[l-1])
                end
            end

            @inbounds for l in 1:L
                # Update Weights using Momentum Gradient Descent
                #  W^(l) = W^(l) - α∇E - η∇E_old
                net.weights[l] -= α*grad[l] + η*last_grad[l]
                # Save last gradients
                last_grad[l][:] = grad[l][:]
                # Reset gradient component for the next batch
                grad[l][:] = 0.0
            end
        end
    end
end

ext: jl


text:
using MDCDL
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
# include("parallelFB.jl")
# include("convolutionalOperator.jl")
include("multidimensionals.jl")
include("permutations.jl")
include("orthonormal.jl")
include("cnsolt.jl")
include("rnsolt.jl")

# include("multiscale.jl")

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# van_der_Pol_Poincare_section_v1.jl
# =============================================

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using StaticArrays
using Random
using Statistics

function van_der_Pol_ODE(du,u,p,t)
    d, a, omega = p
    du[1] =  u[2]
    du[2] =  - u[1] - d*(u[1]*u[1]-1) * u[2] + a*sin(omega*t)
end

# parameters and initial values
# -----------------------------
ntrans = 5000.   # transient periods
nattr = 2000    # no. of points on the attractor used for estimating winding no.
y_init = [0.1 0.1]

d = 5.
a = 1.
omega = 1.

p = [d a omega]
period = 2*pi / omega

 # transient to attractor
 # ----------------------
tspan = (0., ntrans*period)
prob = ODEProblem(van_der_Pol_ODE, y_init, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, maxiters = 1e12 )

# attractor
# ---------
tspan = (0., nattr*period)
y_init = sol.u[end]
prob = ODEProblem(van_der_Pol_ODE, y_init, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, maxiters = 1e12,  saveat = period)

u1vec = [u[1] for u in sol.u]
u2vec = [u[2] for u in sol.u]
nts = length(u1vec)

theta_vec = (atan.( u2vec , u1vec  ) .+ pi ) ./ (2*pi)

# open plot
# ---------
# %%
fig = figure(figsize=(0.6*figx,figy))

ax1 = subplot(1,2,1)
plt.scatter(u1vec,u2vec,s=4, c = COLORS[1])

xticks([-2,0,2], ["-2","","2"])
yticks([-6,0,6], ["-6","","6"])
xlabel(L"u_1";labelpad = -20)
ylabel(L"u_2";labelpad = -20)

ax2 = subplot(1,2,2)
plt.scatter(theta_vec[1:end-1],theta_vec[2:end],s=4, c = COLORS[2])
xlabel(L"\theta_n"; labelpad = -20)
ylabel(L"\theta_{n+1}"; labelpad = -20)
xticks([0,0.5,1], ["0", "", "2π"])
yticks([0,0.5,1], ["0", "", "2π"])
yticks([0,0.5,1])
xlim(0,1)
ylim(0,1)
add_identifiers!(fig)
fig.tight_layout(pad=0.33)

wsave(plotsdir("9", "vanderPol_poincare_section"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Comonicon
using AdvancedHMC: AdvancedHMC

println("Environment variables for testing")
println(ENV)

const DIRECTORY_AdvancedHMC = dirname(dirname(pathof(AdvancedHMC)))
const DIRECTORY_Turing_tests = joinpath(DIRECTORY_AdvancedHMC, "test", "turing")
const GROUP = get(ENV, "AHMC_TEST_GROUP", "AdvancedHMC")

if GROUP == "All" || GROUP == "AdvancedHMC"
    using ReTest, CUDA

    include("metric.jl")
    include("hamiltonian.jl")
    include("integrator.jl")
    include("trajectory.jl")
    include("adaptation.jl")
    include("sampler.jl")
    include("sampler-vec.jl")
    include("demo.jl")
    include("models.jl")
    include("abstractmcmc.jl")

    if CUDA.functional()
        include("cuda.jl")
    else
        @warn "Skipping GPU tests because no GPU available."
    end

    @main function runtests(patterns...; dry::Bool=false)
        retest(patterns...; dry=dry, verbose=Inf)
    end
end

if GROUP == "All" || GROUP == "Downstream"
    using Pkg
    try
        # activate separate test environment
        Pkg.activate(DIRECTORY_Turing_tests)
        Pkg.develop(PackageSpec(; path=DIRECTORY_AdvancedHMC))
        Pkg.instantiate()

        # make sure that the new environment is considered `using` and `import` statements
        # (not added automatically on Julia 1.3, see e.g. PR #209)
        if !(joinpath(DIRECTORY_Turing_tests, "Project.toml") in Base.load_path())
            pushfirst!(LOAD_PATH, DIRECTORY_Turing_tests)
        end

        # Avoids conflicting namespaces, e.g. `NUTS` used in Turing.jl's tests
        # refers to `Turing.NUTS` not `AdvancedHMC.NUTS`.
        include(joinpath("turing", "runtests.jl"))
    catch err
        err isa Pkg.Resolve.ResolverError || rethrow()
        # If we can't resolve that means this is incompatible by SemVer and this is fine
        # It means we marked this as a breaking change, so we don't need to worry about
        # Mistakenly introducing a breaking change, as we have intentionally made one
        @info "Not compatible with this release. No problem." exception = err
    end
end


text:
using ExactOptimalTransport

using LinearAlgebra
using Random
using Test
using Distributions

Random.seed!(100)

@testset "utils.jl" begin
    @testset "FiniteDiscreteMeasure" begin
        @testset "Univariate Finite Discrete Measure" begin
            n = 100
            m = 80
            μsupp = rand(n)
            νsupp = rand(m)
            μprobs = normalize!(rand(n), 1)

            μ = ExactOptimalTransport.discretemeasure(μsupp, μprobs)
            ν = ExactOptimalTransport.discretemeasure(νsupp)
            # check if it vectors are indeed probabilities
            @test isprobvec(μ.p)
            @test isprobvec(probs(μ))
            @test ν.p == ones(m) ./ m
            @test probs(ν) == ones(m) ./ m

            # check if it assigns to DiscreteNonParametric when Vector/Matrix is 1D
            @test μ isa DiscreteNonParametric
            @test ν isa DiscreteNonParametric

            # check if support is correctly assinged
            @test sort(μsupp) == μ.support
            @test sort(μsupp) == support(μ)
            @test sort(vec(νsupp)) == ν.support
            @test sort(vec(νsupp)) == support(ν)
        end
        @testset "Multivariate Finite Discrete Measure" begin
            n = 10
            m = 3
            μsupp = [rand(m) for i in 1:n]
            νsupp = [rand(m) for i in 1:n]
            μprobs = normalize!(rand(n), 1)
            μ = ExactOptimalTransport.discretemeasure(μsupp, μprobs)
            ν = ExactOptimalTransport.discretemeasure(νsupp)
            # check if it vectors are indeed probabilities
            @test isprobvec(μ.p)
            @test isprobvec(probs(μ))
            @test ν.p == ones(n) ./ n
            @test probs(ν) == ones(n) ./ n

            # check if support is correctly assinged
            @test μsupp == μ.support
            @test μsupp == support(μ)
            @test νsupp == ν.support
            @test νsupp == support(ν)
        end
    end
end


text:
### Global constant allowing to chose between forward_diff and zygote_diff for hyperparameter optimization ###
const ADBACKEND = Ref(:Zygote)

function setadbackend(ad_backend::Symbol)
    (ad_backend == :ForwardDiff || ad_backend == :Zygote) ||
        error("Wrong AD Backend symbol, options are :ForwardDiff or :Zygote")
    return ADBACKEND[] = ad_backend
end

# opt(::AbstractInducingPoints) = nothing
# opt(Z::OptimIP) = Z.opt
# data(Z::OptimIP) = Z.Z
# data(Z::AbstractInducingPoints) = Z

## Generic fallback when gradient is nothing
update_kernel!(::Any, ::Any, ::Nothing) = nothing
update_Z!(::Any, ::Any, ::Nothing) = nothing

## Generic fallback when optimizer is nothing
update_kernel!(::Nothing, ::Kernel, ::NamedTuple) = nothing
update_kernel!(::Nothing, ::AbstractVector, ::AbstractArray) = nothing
update_Z!(::Nothing, ::AbstractVector, ::AbstractVector) = nothing
update_Z!(::Nothing, ::AbstractVector, ::NamedTuple) = nothing

## Updating prior mean parameters ##
function update!(μ::PriorMean, g::AbstractVector, X::AbstractVector)
    return update!(μ, g, X)
end

## Updating kernel parameters ##

## ForwardDiff.jl approach (with destructure())
function update!(opt, k::Kernel, Δ::AbstractVector)
    ps = params(k)
    i = 1
    for p in ps
        d = length(p)
        Δp = Δ[i:(i + d - 1)]
        Δlogp = Optimise.apply!(opt, p, p .* Δp)
        @. p = exp(log(p) + Δlogp)
        i += d
    end
end

## Zygote.jl approach with named tuple
function update_kernel!(opt, k::Union{Kernel,Transform}, g::NamedTuple)
    foreach(pairs(g)) do (fieldname, grad)
        update_kernel!(opt, getfield(k, fieldname), grad)
    end
end

function update_kernel!(opt, x::AbstractArray, g::AbstractArray)
    Δ = Optimise.apply!(opt, x, x .* g)
    @. x = exp(log(x) + Δ) # Always assume that parameters need to be positive
end

## Updating inducing points
function update_Z!(opt, Z::Union{ColVecs,RowVecs}, Z_grads::NamedTuple)
        Z.X .+= Optimise.apply!(opt, Z.X, Z_grads.X)
end

function update_Z!(opt, Z::AbstractVector, Z_grads::AbstractVector)
    for (z, zgrad) in zip(Z, Z_grads)
        z .+= Optimise.apply!(opt, z, zgrad)
    end
end


text:
using Bridge, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models

T = 2.0
dt = 1/1000

tt = 0.:dt:T
struct PBIntegratedDiffusion <: ContinuousTimeProcess{Float64}
    γ::Float64
end
struct PBIntegratedDiffusionAux <: ContinuousTimeProcess{Float64}
    γ::Float64
end

PorPtilde = Union{PBIntegratedDiffusion, PBIntegratedDiffusionAux}

βu(t, x::Float64, P::PBIntegratedDiffusion) = - (x+sin(x)) + 1/2
βu(t, x::Float64, P::PBIntegratedDiffusionAux) = -x + 1/2
# not really a 'beta'

Bridge.b(t::Float64, x, P::PorPtilde) = Bridge.b!(t, x, copy(x), P)
function Bridge.b!(t, x, out, P::PorPtilde)
    out[1] = x[2]
    out[2] = βu(t, x[2], P)
    out
end

function Bridge.σ!(t, x, dm, out, P::PorPtilde)
    out[1] = 0.0
    out[2] = dm*P.γ
    out
end
function Bridge.a(t, P::PorPtilde)
    [0.0 0.0; 0.0 P.γ^2]
end
Bridge.a(t, x, P::PorPtilde) = Bridge.a(t, P::PorPtilde)

Bridge.constdiff(::PorPtilde) = true

function Bridge.B!(t, arg, out, P::PBIntegratedDiffusionAux)
    B = [0.0 1.0; 0.0 -1.0]
    out .= (B*arg)
    out
end
function BBt!(t, arg, out, P::PBIntegratedDiffusionAux)
    B = [0.0 1.0; 0.0 -1.0]
    out .= (B*arg + arg*B')
    out
end

function Bridge.dP!(t, p, out, P)
    BBt!(t, p, out, P)
    out[2,2] -= P.γ^2
    out
end

# Generate Data
Random.seed!(1)

P = PBIntegratedDiffusion(0.7)
Pt = PBIntegratedDiffusionAux(0.7)

W1 = sample(tt, Wiener())
x0 = [2.0, 1.0]
X1 = solve(EulerMaruyama!(), x0, W1, P)

L = [1. 0.]
Σnoise = fill(0.01, 1, 1)
v = [2.5]

# Solve Backward Recursion

N = length(tt)
νt = [zero(x0) for i in 1:N]
Σt = [zero(Bridge.outer(x0)) for i in 1:N]
ϵ = 0.02

Bridge.partialbridgeode!(Bridge.R3!(), tt, L, Σnoise, v, νt, Σt, Pt, ϵ)

j = length(tt)÷2
B = [0.0 1.0; 0.0 -1.0]
A = Bridge.a(NaN, Pt)
@test norm((νt[j+1] - νt[j])/dt - (Bridge.b(tt[j+1],  νt[j], Pt))) < 0.01
j = length(tt) - 3
@test norm((B*Σt[j+1] + Σt[j+1]*B' - Bridge.a(tt[j+1], Pt)) - Bridge.dP!(tt[j+1], Σt[j+1], copy(Σt[j+1]), Pt)) < 0.01
@test_broken norm((Σt[j+1] - Σt[j])/dt - (B*Σt[j+1] + Σt[j+1]*B' - Bridge.a(tt[j+1], Pt))) < 0.01

Po1 = Bridge.PartialBridge!(tt, P, Pt, L, v, ϵ, Σnoise)

@test Po1.H == inv.(Σt)

x0 = [2.0, 1.0]
Xo1 = deepcopy(X1)

solve!(Bridge.EulerMaruyama!(), Xo1, x0, W1, Po1)

# Likelihood

ll1 = llikelihood(Bridge.LeftRule(), Xo1, Po1)

lls1 = Float64[]

function mcmc1(x0, tt, Po2)
    # MCMC parameter

    iterations = 1000
    subsamples = 0:100:iterations
    ρ = 0.9

    # initalization
    W = sample(tt, Wiener())
    X = solve(EulerMaruyama!(), x0, W, Po1)
    ll = llikelihood(Bridge.LeftRule(), X, Po1)

    acc = 0

    Wo = deepcopy(W)
    Wrho = deepcopy(W)

    Xo = deepcopy(X)

    @assert !(X.yy[1] === Xo.yy[1])

    XX = Any[]
    if 0 in subsamples
        push!(XX, deepcopy(X))
    end
    for iter in 1:iterations
        # Proposal
        sample!(Wrho, Wiener())
        Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*Wrho.yy

        solve!(EulerMaruyama!(), Xo, x0, Wo, Po1)
        llo = llikelihood(Bridge.LeftRule(), Xo, Po1)
        push!(lls1, llo)

        if log(rand()) <= llo - ll

            X, Xo = Xo, X
            W, Wo = Wo, W
            ll = llo
            acc += 1
        end
        if iter in subsamples
            push!(XX, deepcopy(X))
        end
    end

end

@time @testset "MCMC1" begin
    mcmc1(x0, tt, Po2)

end


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using Test
using CausalityToolsBase
using StateSpaceReconstruction
using PerronFrobenius

# Individual constructors (i.e. rectangularinvariantmeasure(pts, ϵ))
include("test_rectangularmeasure_constructor.jl")
#include("test_inducedmeasure.jl")
#include("test_inducedmeasure_constructor.jl")
#include("test_averagemeasure.jl")
#include("test_averagemeasure_constructor.jl")

# Test the estimators (i.e. invariantmeasure(pts, binningscheme))
include("test_invariantmeasure_estimators.jl")

include("test_gridestimator.jl")
include("test_transferoperator_rectangular.jl")

#include("triangulations.jl")
include("test_triangulationestimators.jl")


text:
module Example

using Glimmer, Glimmer.FlexUI

println("Start [$(splitext(basename(@__FILE__))[1])]")

# ---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
# ---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Grid Types Tree")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

# ---------------------------------------------------------------
# load the data
# ---------------------------------------------------------------
types=Vector{Any}()
# this function do not belong in this package and is just included for this example.
# it returns a table of all the types/functions/modules inside the requsted module (Glimmer in this case).
GridUtils.getTypes(Glimmer, types)
table = GridUtils.toJSONTable(types)

# ---------------------------------------------------------------
# Define Variables
# ---------------------------------------------------------------
type_doc = addVariable!(app, Variable(name="type_doc", type="string", value="No Selected Type"))
table_var = addVariable!(app, Variable(name="table_var", type="aggrid", value=""))

# setup the grid data and options - will get friendlier in the near future
# for tree type tables, we need to supply the name of the "id" firld and the "parent" field.
table_var[] = GridUtils.table2agGridTree(table, "_id", "_parentId")

# setting some grid options - because of it's complexity, the grid's data and visualization properties
# are all stored in the variable and not the UI control.
gridOption!(table_var, "rowSelection", "single")
gridOption!(table_var, "showRowIndex", true)
gridOption!(table_var, "treeViewIndentPixels", 30)
gridOption!(table_var, "suppressScrollOnNewData", true)
gridOption!(table_var, "showTreeViewFilter", true)
forEachGridColumn!(table_var, "editable", false)

# override the default column definitions - don't want to show all the fields in the table
clearGridColumns!(table_var)
addGridColumn!(table_var, gridTreeViewCell(header="Type Name", eventName="myTreeFieldClicked", field="name", width=350))
addGridColumn!(table_var, gridCell(header="Description", field="desc"))
addGridColumn!(table_var, gridCell(header="Type", field="type"))
addGridColumn!(table_var, gridCell(header="Has Documentation", field="hasDocs"))

# define the column default values (note the cellStyle to color a line red according to documantation existance)
gridDefaultColDef!(table_var, gridCell(
    editable=false,
    resizable=true,
    cellStyle="__js__ params => params.data.hasDocs==false ? { color: 'red' } : { color: 'green' }",
))

# ---------------------------------------------------------------
# Define Events
# ---------------------------------------------------------------
# see https://www.ag-grid.com/angular-data-grid/grid-events/ for details

# this is an example of how to convert the grid variable back into a table - not used in this example
on(table_var, :valueChanged) do val
    @show "ON valueChanged"
    new_table = Tables.matrix(table_var[]);
    @show new_table
end

on(table_var, :cellClicked) do val
    @show "ON cellClicked"
end

on(table_var, :cellDoubleClicked) do val
    @show "ON cellDoubleClicked"
end

on(table_var, :rowSelected) do val
    @show "ON rowSelected"
end

# when selection changed we update the documentation display for the new selected type
on(table_var, :selectionChanged) do val
    @show "ON selectionChanged"

    rows = val["selectedRows"];
    if (length(rows) == 0)
        type_doc[] = "No Selected Type";
        return
    end
    if (length(rows) > 1)
        @warn "currently can't handle more than a single selected row - will process first one only"
    end

    doc = rows[1]["docs"]
    type_doc[] = (doc!=="") ? doc : "No Selected Type"
end

# example - will use in the future
on(table_var, :myGridButtonEvent) do val
    @show "ON myGridButtonEvent" val
end

# ---------------------------------------------------------------
# Define Controls
# ---------------------------------------------------------------
ui = VContainer(
    Card(
        title="This Example demonstrate some the the AGGrid component capabilities",
        subtitle="The grid show a tree structure representing the types/functions/modules inside the Glimmer package.",
        content=HContainerStartStart(
            AGGrid(
                style="width: 500px; max-width: 500px; height: 500px;",
                variable="table_var",
            ),
            Card(
                title="Documentation",
                # experementing with different styles because i have no clue about html/css
                contentStyle="max-width: 500px; background-color: lightgreen; border: 1px solid black;overflow:auto;",
                # contentStyle="max-width: 500px; overflow:auto;",
                content=HContainer(
                    Markdown(
                        content="\$(type_doc)",
                    ),
                ),
            ),
        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

# ---------------------------------------------------------------
# Run the application
# ---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
@testset "IndOp" begin
    op = IndOp(Paulis.X, 1)
    @test copy(op) == op

    t1 = PauliTerm("XIIYI", complex(2.0))
    t2 = PauliTerm("IZIYZ", complex(3.0))
    ot1 = OpTerm{IndOp}(t1)
    ot2 = OpTerm{IndOp}(t2)
    @test ot1 * ot2 == OpTerm{IndOp}(t1 * t2)

    ps = PauliSum([PauliTerm("IZYX", 1 + 0im),PauliTerm("XZXX", 1 + 0im),PauliTerm("ZXXX", 1 + 0im),])
    ops = OpSum{IndOp}(ps)
    @test map(weight, ps) == map(weight, ops)
    @test ops isa OpSum{<:IndOp}
    @test typeof(ops) == typeof(OpSum{IndOp}(ops))

    ft = FermiTerm("ENE+")
    oft = OpTerm{IndOp}(ft)
    @test count_bodies(ft) == 3
    @test count_bodies(oft) == 3
end


text:
using DitherPunk
using DitherPunk: gradient_image, test_on_gradient
using Images

# # Gallery
# ## On color spaces
# A simple linear gradient works well to reveal the characteristic patterns of
# different dithering algorithms.
srbg, linear = gradient_image(100, 800);
mosaicview(srbg, linear)
# The pixel intensities in the image `srgb` increase linearly from 0 to 1.
# The second image `linear` has been converted from sRGB to a linear representation,
# which more closely matches our human perception of brightness.
#
# The helper function `test_on_gradient` takes a dithering algorithm and runs it on
# both the `srgb` and the `linear` image.
# It then shows a comparison of both inputs and outputs.
#
# Most dithering algorithms in DitherPunk.jl provide an optional parameter `to_linear`,
# which converts the input image to a linear color space before applying the dithering.
# Select what looks best!
#
# ## Threshold dithering
# #### `ConstantThreshold`
test_on_gradient(ConstantThreshold())

# #### `WhiteNoiseThreshold`
test_on_gradient(WhiteNoiseThreshold())

# ## Ordered dithering
# ### Bayer matrices
# #### `bayer_dithering`
test_on_gradient(Bayer())

# The order of the Bayer-matrix can be specified through the parameter `level`,
# which defaults to `1`.
# **Level 2**
test_on_gradient(Bayer(; level=2))
# **Level 3**
test_on_gradient(Bayer(; level=3))
# **Level 4**
test_on_gradient(Bayer(; level=4))

# ### Clustered / halftone dithering
# The following methods have large characteristic patterns and are therefore
# better suited for large images.
#
# #### `ClusteredDots`
test_on_gradient(ClusteredDots())

# #### `CentralWhitePoint`
test_on_gradient(CentralWhitePoint())

# #### `BalancedCenteredPoint`
test_on_gradient(BalancedCenteredPoint())

# #### `Rhombus`
test_on_gradient(Rhombus())

# #### ImageMagick methods
test_on_gradient(IM_checks())
#
test_on_gradient(IM_h4x4a())
#
test_on_gradient(IM_h6x6a())
#
test_on_gradient(IM_h8x8a())
#
test_on_gradient(IM_h4x4o())
#
test_on_gradient(IM_h6x6o())
#
test_on_gradient(IM_h8x8o())
#
test_on_gradient(IM_c5x5())
#
test_on_gradient(IM_c6x6())
#
test_on_gradient(IM_c7x7())

# ## Error diffusion
# #### `SimpleErrorDiffusion`
test_on_gradient(SimpleErrorDiffusion())

# #### `FloydSteinberg`
test_on_gradient(FloydSteinberg())

# #### `JarvisJudice`
test_on_gradient(JarvisJudice())

# #### `Stucki`
test_on_gradient(Stucki())

# #### `Burkes`
test_on_gradient(Burkes())

# #### `Sierra`
test_on_gradient(Sierra())

# #### `TwoRowSierra`
test_on_gradient(TwoRowSierra())

# #### `SierraLite`
test_on_gradient(SierraLite())

# #### `Fan93()`
test_on_gradient(Fan93())

# #### `ShiauFan()`
test_on_gradient(ShiauFan())

# #### `ShiauFan2()`
test_on_gradient(ShiauFan2())

# #### `Atkinson()`
test_on_gradient(Atkinson())


text:
##################################################
## Classical derivative-based, iterative, root-finding algorithms: Newton and Halley
## Historic, we have derivative free versions of similar order

## Newton
"""

    Roots.Newton()

Implements Newton's [method](http://tinyurl.com/b4d7vls): `x_n1 = xn -
f(xn)/f'(xn)`.  This is a quadratically converging method requiring
one derivative. Two function calls per step.
"""
struct Newton <: AbstractUnivariateZeroMethod
end

function update_state(method::Newton, fs, o, options)
    xn = o.xn1
    fxn = o.fxn1
    fpxn = fs(xn,1)

    if isissue(fpxn)
        o.stopped=true
        return
    end

    xn1 = xn - fxn / fpxn
    fxn1 = fs(xn1)
    incfn(o)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1

end

"""
    newton(f, fp, x0; kwargs...)

Implementation of Newton's method: `x_n1 = x_n - f(x_n)/ f'(x_n)`

Arguments:

* `f::Function` -- function to find zero of

* `fp::Function` -- the derivative of `f`.

* `x0::Number` -- initial guess. For Newton's method this may be complex.

With the `FowardDiff` package derivatives may be computed automatically. For example,  defining
`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` to be used for the first derivative.

Keyword arguments are passed to `find_zero` using the `Roots.Newton()` method.

"""
newton(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Newton(); kwargs...)
@deprecate newton(f, x0; kwargs...)  newton(f, fp, x0; kwargs...)

## Halley

"""
    Roots.Halley()

Implements Halley's [method](http://tinyurl.com/yd83eytb),
`x_n1 = xn - (2 f(xn)*f'(xn)) / (2 f'(xn)^2 - f(xn) * f''(xn))`.
This method is cubically converging, but requires more function calls per step (3) than
other methods.
"""
struct Halley <: AbstractUnivariateZeroMethod
end

function update_state(method::Halley, fs, o::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}
    xn = o.xn1
    fxn = o.fxn1
    fpxn = fs(xn,1); incfn(o)
    fppxn = fs(xn,2); incfn(o)

    xn1 = xn - 2fxn*fpxn / (2*fpxn*fpxn - fxn * fppxn)
    fxn1 = fs(xn1); incfn(o)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1
end

"""
    halley(f, fp, fpp, x0; kwargs...)

Implementation of Halley's method. `xn1 = xn - 2f(xn)*f'(xn) / (2*f'(xn)^2 - f(xn) * f''(xn))`

Arguments:

* `f::Function` -- function to find zero of

* `fp::Function` -- derivative of `f`.

* `fpp:Function` -- second derivative of `f`.

* `x0::Number` -- initial guess

With the `FowardDiff` package derivatives may be computed automatically. For example,  defining
`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second
derivatives, respectively.

Keyword arguments are passed to `find_zero` using the `Roots.Halley()` method.

"""
halley(f, fp, fpp, x0; kwargs...) = find_zero((f, fp, fpp), x0, Halley(); kwargs...)
#halley(f,  x0; kwargs...) = find_zero(f, x0, Halley(); kwargs...) # deprecated
#halley(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Halley(); kwargs...) # deprecated
@deprecate halley(f,  x0; kwargs...)    halley(f, fp, fpp, x0; kwargs...)
@deprecate halley(f, fp, x0; kwargs...) halley(f, fp, fpp, x0; kwargs...)


text:
#--------------------------------------------------------------------------------------------------
# Tuning the constant stepsize for SGD and SMD - E-NET VERSION
#--------------------------------------------------------------------------------------------------

using Random
using PyPlot
using LaTeXStrings, Printf

include("solve_enet.jl");

#-------------------------------------
#   Parameters
#-------------------------------------
maxIter = 7000;
r = 5;  # rank
d = 100;
stoch_err = sqrt(10);  # standard deviation of errors in stochastic measurements b
init_radius = 1.0;  # 2-norm measure of relative deviation between Xtrue and Xinit

num_tests = 5;
subgrad_stepsizes = [1e-6, 2e-6, 3e-6, 4e-6, 5e-6];
mirror_stepsizes = [0.3, 0.4, 0.5, 0.6, 0.7];

Random.seed!(321);  # for reproducibility

#-------------------------------------
#   Initialize plots
#-------------------------------------
subgrad_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel(L"Relative distance to solution set $\min_U \; \frac{\| X_k U - X_{true} \|_2^2 }{\| X_{true} \|_2^2}$");
title(@sprintf("SGD for e-net covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)
ylim(1e-4, 1e1)

subgrad_fun_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel("Function Error");
title(@sprintf("SGD for e-net covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)

mirror_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel(L"Relative distance to solution set $\min_U \; \frac{\| X_k U - X_{true} \|_2^2 }{\| X_{true} \|_2^2}$");
title(@sprintf("SMD for e-net covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)
ylim(1e-4, 1e1)

mirror_fun_plot = figure(figsize=[10,6]);
xlabel(L"Iteration $k$");
ylabel("Function Error");
title(@sprintf("SMD for e-net covariance estimation (r=%i, d=%i)", r,d));
xlim(0,maxIter)

#-----------------------------------------------------------------------------------------
#   Run method(s) for each dimension
#-----------------------------------------------------------------------------------------
for i=1:num_tests
    ηsubgrad = subgrad_stepsizes[i];
    ηmirror = mirror_stepsizes[i];

    # Generate true matrix we are searching for
    Xtrue = randn(d,r);

    # Generate initial point
    pert = randn(d,r);
    Xinit = Xtrue + init_radius * (norm(Xtrue, 2) / norm(pert, 2)) * pert;

    # Run the two methods
    (subgrad_err_hist, subgrad_fun_hist) = solve_enet(Xinit, Xtrue, stoch_err,  maxIter,
                                                            fill(ηsubgrad, maxIter), method="subgradient",
                                                            noiseType="sparse")
    (mirror_err_hist, mirror_fun_hist) = solve_enet(Xinit, Xtrue, stoch_err,  maxIter,
                                                                fill(ηmirror, maxIter), method="mirror",
                                                                noiseType="sparse")

    # Add errors to plots
    plt[:figure](subgrad_plot[:number])
    semilogy(subgrad_err_hist, #color=step_colors[i],
                    label=@sprintf("step=%1.2e", ηsubgrad));

    plt[:figure](subgrad_fun_plot[:number])
    semilogy(subgrad_fun_hist, #color=step_colors[i],
                    label=@sprintf("step=%1.2e", ηsubgrad));

    plt[:figure](mirror_plot[:number])
    semilogy(mirror_err_hist, #color=step_colors[i],
                    label=@sprintf("step=%1.2e", ηmirror));

    plt[:figure](mirror_fun_plot[:number])
    semilogy(mirror_fun_hist, #color=step_colors[i],
                    label=@sprintf("step=%1.2e", ηmirror));
end

# Add legends to plots and save final plots
plt[:figure](subgrad_plot[:number])
legend(loc="upper right")
savefig("plots/sgd_enet_stepsize_test.pdf");

plt[:figure](subgrad_fun_plot[:number])
legend(loc="upper right")
savefig("plots/sgd_enet_stepsize_test_fun_vals.pdf");

plt[:figure](mirror_plot[:number])
legend(loc="upper right")
savefig("plots/smd_enet_stepsize_test.pdf");

plt[:figure](mirror_fun_plot[:number])
legend(loc="upper right")
savefig("plots/smd_enet_stepsize_test_fun_vals.pdf");

max_stars_repo_name: kelliemac/stochastic-model-based

max_issues_repo_name: kelliemac/stochastic-model-based

max_forks_repo_name: kelliemac/stochastic-model-based


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using Flux
using DifferentialEquations: Tsit5

import Random
Random.seed!(1234);

t_start = 0.0
t_step = 0.01
t_stop = 1.0 # 5.0
tData = t_start:t_step:t_stop

# generate traing data
myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
parameters = Dict("mass_s0" => 1.3)
realSimData = fmiSimulateCS(myFMU, t_start, t_stop; parameters=parameters, recordValues=["mass.a"], saveat=tData)

# sine(t) as external force
function extForce(t)
    return [sin(t)]
end
accData = fmi2GetSolutionValue(realSimData, "mass.a")

# loss function for training
function losssum()
    solution = problem(extForce, t_step)

    accNet = fmi2GetSolutionValue(solution, 1; isIndex=true)

    Flux.Losses.mse(accNet, accData)
end

# callback function for training
iterCB = 0
lastLoss = 0.0
function callb()
    global iterCB += 1
    global lastLoss

    if iterCB == 1
        lastLoss = losssum()
    end

    if iterCB % 50 == 0
        loss = losssum()
        @info "Loss: $loss"
        @test loss < lastLoss
        lastLoss = loss
    end
end

# NeuralFMU setup
numInputs = length(myFMU.modelDescription.inputValueReferences)
numOutputs = length(myFMU.modelDescription.outputValueReferences)

net = Chain(inputs -> fmiInputDoStepCSOutput(myFMU, t_step, inputs),
            Dense(numOutputs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numOutputs))

problem = CS_NeuralFMU(myFMU, net, (t_start, t_stop); saveat=tData)
@test problem != nothing

# train it ...
p_net = Flux.params(problem)

optim = ADAM()
Flux.train!(losssum, p_net, Iterators.repeated((), 100), optim; cb=callb)

fmiUnload(myFMU)


text:
function bench_mod(a,b,c)
    for m in b
        for i in eachindex(a)
            c[i] = mod(a[i],m)
        end
    end
end
function bench_fastmod(a,b,c)
    for m in b
        modm = Mod(m)
        for i in eachindex(a)
            c[i] = fastmod(a[i],modm)
        end
    end
end
function bench_divtest(a,b,c)
    for m in b
        for i in eachindex(a)
            c[i] = (mod(a[i],m) == 0)
        end
    end
end
function bench_fast_divtest(a,b,c)
    for m in b
        modm = Mod(m)
        for i in eachindex(a)
            c[i] = divtest(a[i],modm)
        end
    end
end
@testset "unsigned modulus performance" begin
    a = convert.(UInt32,collect(100:100_000))
    b = convert.(UInt32,collect(2:1000))
    c1 = similar(a)
    c2 = similar(a)

    bench_fastmod = @benchmark bench_fastmod($a,$b,$c2)
    bench_mod = @benchmark bench_mod($a,$b,$c1)
    display(bench_fastmod)
    display(bench_mod)
    @test minimum(bench_fastmod.times) < minimum(bench_mod.times)/2 #should be at least twice as fast for these sizes
end

@testset "unsigned divtest performance" begin
    a = convert.(UInt32,collect(100:100_000))
    b = convert.(UInt32,collect(2:1000))
    c1 = BitVector(undef,length(a))
    c2 = BitVector(undef,length(a))

    bench_fastmod = @benchmark bench_fast_divtest($a,$b,$c2)
    bench_mod = @benchmark bench_divtest($a,$b,$c1)
    display(bench_fastmod)
    display(bench_mod)
    @test minimum(bench_fastmod.times) < minimum(bench_mod.times) #shockingly slow
end


text:
# interfacing with NNlib.jl

import ..CuArrays: CuVecOrMat, CuVector

import NNlib: conv!, ∇conv_filter!, ∇conv_data!, stride, dilation, flipkernel,
  maxpool!, meanpool!, ∇maxpool!, ∇meanpool!, spatial_dims, padding, kernel_size,
  softmax, softmax!, ∇softmax!, logsoftmax, logsoftmax!, ∇logsoftmax

# Softmax

const CUDNNFloat = Union{Float16,Float32,Float64}

reshape4D(x::AbstractVector) = reshape(x, 1, 1, length(x), 1)
reshape4D(x::AbstractMatrix) = reshape(x, 1, 1, size(x)...)

function softmax!(out::CuVecOrMat{T}, xs::CuVecOrMat{T}) where T<:CUDNNFloat
  cudnnSoftmaxForward(reshape4D(xs), reshape4D(out))
  return out
end

function ∇softmax!(out::CuVecOrMat{T}, Δ::CuVecOrMat{T}, xs::CuVecOrMat{T}) where T<:CUDNNFloat
  cudnnSoftmaxBackward(reshape4D(softmax(xs)), reshape4D(Δ), reshape4D(out))
  return out
end

function logsoftmax!(out::CuVecOrMat{T}, xs::CuVecOrMat{T}) where T<:CUDNNFloat
  cudnnSoftmaxForward(reshape4D(xs), reshape4D(out), algorithm=CUDNN_SOFTMAX_LOG)
  return out
end

function ∇logsoftmax!(out::CuVecOrMat{T}, Δ::CuVecOrMat{T}, xs::CuVecOrMat{T}) where T<:CUDNNFloat
  cudnnSoftmaxBackward(reshape4D(logsoftmax(xs)), reshape4D(Δ), reshape4D(out);
                       algorithm=CUDNN_SOFTMAX_LOG)
  return out
end

∇logsoftmax(Δ::CuVecOrMat{T}, xs::CuVecOrMat{T}) where T<:CUDNNFloat =
  ∇logsoftmax!(similar(xs), Δ, xs)

# Convolution

# Since CUDNN does not support 1D convolution, Conv in Flux will give a CUDNNError if the size is 1-dimensional.
# We have to reshape the CuArray/PoolDims/DenseConvDims to 4D before feeding to CUDNN.
fix1d(x) = x

fix1d(x::CuArray{T, 3}) where T = reshape(x, size(x, 1), 1, size(x, 2), size(x, 3))

fix1d(cdims::DenseConvDims{1,K,C_in,C_out,S,P,D,F}) where {K,C_in,C_out,S,P,D,F} =
  DenseConvDims{2,(K...,1),C_in,C_out,(S...,1),(P...,0,0),(D...,1),F}((cdims.I...,1))

fix1d(pdims::PoolDims{1,K,S,P,D}) where {K,S,P,D,F} =
  PoolDims{2,(K...,1),(S...,1),(P...,0,0),(D...,1)}((pdims.I..., 1), pdims.C_in)

function conv!(y::CuArray{T}, x::CuArray{T}, w::CuArray{T}, cdims::DenseConvDims;
               alpha=1, algo=0) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end
  cudnnConvolutionForward(fix1d(y), fix1d(x), fix1d(w), fix1d(cdims), alpha=alpha, algo=algo)
  return y
end

function ∇conv_filter!(dw::CuArray{T}, x::CuArray{T}, dy::CuArray{T},
                       cdims::DenseConvDims; alpha=1, algo=0) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  cudnnConvolutionBackwardFilter(fix1d(dw), fix1d(x), fix1d(dy), fix1d(cdims), alpha=alpha, algo=algo)
  return dw
end

function ∇conv_data!(dx::CuArray{T}, dy::CuArray{T}, w::CuArray{T},
                     cdims::DenseConvDims; alpha=1, algo=0) where T<:CUDNNFloat
  if version() < v"6"
    all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end

  cudnnConvolutionBackwardData(fix1d(dx), fix1d(w), fix1d(dy), fix1d(cdims), alpha=alpha, algo=algo)
  return dx
end

∇conv_bias!(db::CuArray{T}, dy::CuArray{T}; alpha=1, beta=0) where T<:CUDNNFloat =
  (cudnnConvolutionBackwardBias(fix1d(db), fix1d(dy), alpha=alpha, beta=beta); return db)

maxpool!(y::CuArray{T}, x::CuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingForward(fix1d(y), fix1d(x), fix1d(pdims); mode=0); return y)

∇maxpool!(dx::CuArray{T}, dy::CuArray{T}, y::CuArray{T}, x::CuArray{T},
          pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingBackward(fix1d(dx), fix1d(dy), fix1d(x), fix1d(y), fix1d(pdims), mode=0); return dx)

meanpool!(y::CuArray{T}, x::CuArray{T}, pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingForward(fix1d(y), fix1d(x), fix1d(pdims), mode=1); return y)

∇meanpool!(dx::CuArray{T}, dy::CuArray{T}, y::CuArray{T}, x::CuArray{T},
           pdims::PoolDims) where T<:CUDNNFloat =
  (cudnnPoolingBackward(fix1d(dx), fix1d(dy), fix1d(x), fix1d(y), fix1d(pdims), mode=1); return dx)


text:
module Simulate_PathPlanning

using  Modia3D

const ptp_path = PTP_path(["angle1", "angle2", "angle3"],
                          positions = [0.0 2.0 3.0;
                                       0.5 3.0 4.0;
                                       0.8 1.5 0.3;
                                       0.2 1.5 0.8],
                          startTime = 0.1,
                          v_max = 2*ones(3),
                          a_max = 3*ones(3))

plotPath(ptp_path)

plotPath(ptp_path, names=("angle2", "angle3"), heading="Test of PTP plots",
         tend = 0.9*ptp_path.Tend, figure=2, ntime=99)

# Test getPosition
pos = getPosition(ptp_path, 2, 0.5)
println("... pos_angle2(time=0.5) = ", pos)

end


text:
type class_nl_function
    n::Int64
    value::Function
    gradient::Function
    hessian::Function

    function class_nl_function(num_var::Int64,val::Function,grad::Function,hess::Function)
        return new(num_var,val,grad,hess)
    end
end

type class_nlp <: internal_AbstractNLPEvaluator
  # linear program with non-linear objective
  # min f(x)
  # g_i(x) = 0
  # x >= 0

	_n::Int64 # number of variables
	_m::Int64 # number of constraints
  obj::class_nl_function
  constraints::Array{class_nl_function,1}

	function class_nlp(number_variables::Int64, obj::class_nl_function, constraints::Array{class_nl_function,1})
		return new(number_variables, length(constraints), obj, constraints);
	end
end

################
# METHODS
################

# evaluate objective
function internal_eval_c(nlp::class_nlp, x::Array{Float64,1})
    try
      return nlp.obj.value(x);
    catch e
        println("ERROR internal_eval_c")
        throw(e)
    end
end

# evalutate constraints
function internal_eval_a(nlp::class_nlp, x::Array{Float64,1})
    try
        num_constr = length(nlp.constraints)
        a = zeros(num_constr)
        for i = 1:num_constr
            a[i] = nlp.constraints[i].value(x)
        end

        return a;
    catch e
        println("ERROR internal_eval_a")
        throw(e)
    end
end

# evaluate gradient of constraints
function internal_eval_jac_a(nlp::class_nlp, x::Array{Float64,1}) # J
    try
      J = spzeros(m(nlp), n(nlp))

      for i = 1:m(nlp)
          J[i,:] = nlp.constraints[i].gradient(x)'
      end

      return J
    catch e
        println("ERROR internal_eval_jac_a")
        throw(e)
    end
end

# hessian of lagrangian
function internal_eval_hesslag_prod(nlp::class_nlp, x::Array{Float64,1}, y::Array{Float64,1})
    try
      H = nlp.obj.hessian(x)

      for i = 1:m(nlp)
          H = H - y[i] * nlp.constraints[i].hessian(x)
      end

      return H
    catch e
        println("ERROR internal_eval_hesslag_prod")
        throw(e)
    end
end

# gradient of lagrangian
function internal_eval_gradlag(nlp::class_nlp, x::Array{Float64,1}, y::Array{Float64,1})
    try
        g = internal_eval_gradc(nlp, x)

        for i = 1:m(nlp)
            g = g - y[i] * nlp.constraints[i].gradient(x)
        end
        return g
    catch e
        println("ERROR internal_eval_gradlag")
        throw(e)
    end
end

# gradient of f
function internal_eval_gradc(nlp::class_nlp, x::Array{Float64,1})
    try
        return nlp.obj.gradient(x);
    catch e
        println("ERROR internal_eval_gradc")
        throw(e)
    end
end

# (\nabla g * x - g)
function internal_eval_b(nlp::class_nlp, x::Array{Float64,1})
    try
        return internal_eval_jac_a(nlp, x) * x - internal_eval_a(nlp, x);
    catch e
        println("ERROR internal_eval_b")
        throw(e)
    end
end

function n(nlp::class_nlp)
    return nlp._n;
end

function m(nlp::class_nlp)
    return nlp._m;
end

max_forks_repo_name: ohinder/Everest


text:
@testset ExtendedTestSet "mask_elements" begin
    @testset ExtendedTestSet "mask_elements" begin
        img = [
            1 10 1 1
            -2 -2 -2 -2
            1 1 1 1
            1 1 -2 1
            1 1 1 1
            -2 -2 -2 -2
        ]
        answer = Bool.([
            1 1 1 1
            0 0 0 0
            1 1 1 1
            1 1 1 1
            1 1 1 1
            0 0 0 0
        ])
        @test mask_elements(img, 0, 2) == answer
    end

    @testset ExtendedTestSet "mask_elements" begin
        img = [
            1 -2 -2 -2
            -2 -2 -2 -2
            1 -2 -2 -2
            -2 -2 -2 -2
            1 -2 -2 -2
            -2 -2 -2 -2
        ]
        answer = Bool.([
            1 0 0 0
            0 0 0 0
            1 0 0 0
            0 0 0 0
            1 0 0 0
            0 0 0 0
        ])
        @test mask_elements(img, 0, 2) == answer
    end
end


text:
# only want to automatically determine the type for holding data T
# when it is not given.
# But do want to assert at compile time that the given type has suitable
# typemin and typemax to hold the range expected

# the use case is making an array of BInt{ 1:10, UInt8 }
# rather than the default BInt{ 1:10, Int64 }

# needs to a Number for promotions to work:
# immutable BInt{ R, T } <: Number

function range_within_range( xl, xu, range_l, range_u )
        range_l <= xl <= range_u && range_l <= xu <= range_u
end

function range_within_range( r1, r2 )
        l1,u1 = get_range_bounds( r1 )
        l2,u2 = get_range_bounds( r2 )
        range_within_range( l1, u1, l2, u2 )
end

function range_union( r1, r2 )
        l1,u1 = get_range_bounds( r1 )
        l2,u2 = get_range_bounds( r2 )
        return make_range( min(l1,l2), max(u1,u2) )
end

function default_int_type( range )
         l,u = get_range_bounds( range )
         T = if range_within_range( l,u, typemin( Int64 ), typemax( Int64 ) )
                 Int64
         elseif range_within_range( l,u, typemin( Int128 ), typemax( Int128 ) )
                 Int128
         else
                 BigInt
         end
         T
end

 # make constructors for above that take a range - and does the
 # number of bits calculation

 function is_signed_2sComplement_range( lo, hi )
         if lo != -hi-1
                 return false
         end
         p = nextpow2( hi )
         p_lo = -p
         p_hi = p-1
         return lo == p_lo && hi == p_hi
 end
 function is_unsigned_binary_range( lo, hi )
         if lo != 0
                 return false
         end
         p = nextpow2( hi )
         p_hi = p-1
         return hi == p_hi
 end

 number_of_bits_unsigned(x::Integer) = x<=0 ? 0 : ndigits(x,2)
 number_of_bits_signed(x::Integer)   = x==0 ? 0 : 1 + (x<0 ? number_of_bits_unsigned(-x-1) : number_of_bits_unsigned(x) )

 function signed_2sComplement_range( lo, hi )
         p1 = lo < 0 ? nextpow2( -lo ) : 0
         p2 = nextpow2( hi )
         p = max( p1, p2 )
         return SignedN(p)
 end

 function unsigned_binary_range( lo, hi )
         assert( lo >= 0 )
         p2 = nextpow2( hi )
         p = max( p1, p2 )
         return UnsignedN(p)
 end

immutable UnsignedN
        n::Int
        UnsignedN( n::Int ) = new(n)
end
immutable SignedN
        n::Int
        SignedN( n::Int ) = new(n)
end

# # Note the difference between the following
# Might need to change this if this is confusing...
# julia> SignedN( 50 )
# SignedN(50)
#
# julia> SignedN( 50:50 )
# SignedN(7)

function UnsignedN( r )
        l,u = get_range_bounds( r )
        UnsignedN( max( number_of_bits_unsigned(l), number_of_bits_unsigned(u) ) )
end

function SignedN( r )
        l,u = get_range_bounds( r )
        SignedN( max( number_of_bits_signed(l), number_of_bits_signed(u) ) )
end

immutable BigRange{ Lower, Upper }
end

get_range_bounds( n::UnsignedN ) = big( 0 ), big(1)<<n.n - 1
get_range_bounds( n::SignedN ) = -big(1)<<(n.n - 1), big(1)<<(n.n - 1) -1
get_range_bounds{ L, U }( r::BigRange{ L, U } ) = big( from_type_constant( L ) ), big( from_type_constant( U ) )
get_range_bounds( n::UnitRange ) = big( n.start ), big( n.stop )
get_range_bounds( n::Integer ) = big( n ), big( n )
get_range_bounds( n::DataType ) = big( typemin(n) ), big( typemax(n) )
get_range_bounds( n::Tuple ) = big( from_type_constant( n ) ), big( from_type_constant( n ) )

function get_range_bounds{T}( r, ::Type{T} )
        l,u = get_range_bounds( r )
        return T(l), T(u)
end

@generated function Base.typemin{R,T}( ::Type{ BInt{R,T} } )
        l,u = get_range_bounds( R )
        return T(l)
end
@generated function Base.typemax{R,T}( ::Type{ BInt{R,T} } )
        l,u = get_range_bounds( R )
        return T(u)
end
@generated function Base.typemin{R}( ::Type{ BInt{R} } )
        T = default_int_type(R)
        l,u = get_range_bounds( R )
        return T( l )
end
@generated function Base.typemax{R}( ::Type{ BInt{R} } )
        T = default_int_type(R)
        l,u = get_range_bounds( R )
        return T( u )
end

 # want this to make a UnitRange if the range values allow
 # else make a BigRange
 # either way - it must be a bitstype
 # so UnitRange{ BigInt } won't work
 function make_range( lower::Integer, upper::Integer )
        if     range_within_range( lower, upper, typemin(Int64 ), typemax(Int64 ) )
                return Int64( lower):Int64( upper)
        elseif range_within_range( lower, upper, typemin(Int128), typemax(Int128) )
                return Int128(lower):Int128(upper)
        else
                return BigRange{ to_type_constant(lower), to_type_constant( upper ) }()
        end
 end

# Is this used?
# # a range specified by a single number is just a constant
# # - can only take on a single value
#  function make_range( n::Integer )
#          return to_type_constant(n)
#  end


text:
############################################################
## joAbstractLinearOperator - extra functions

# elements(jo)
elements(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = joHelpers_etc.elements_helper(A)
elements(A::joMatrix{DDT,RDT}) where {DDT,RDT} = A*jo_speye(DDT,A.n)

# hasinverse(jo)
hasinverse(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = !isnull(A.iop)

# issquare(jo)
issquare(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = (A.m == A.n)

# istall(jo)
istall(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = (A.m > A.n)

# iswide(jo)
iswide(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = (A.m < A.n)

# iscomplex(jo)
iscomplex(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = !(DDT<:Real && RDT<:Real)

# islinear(jo)
function islinear(A::joAbstractLinearOperator{DDT,RDT},samples::Integer=3;tol::Float64=0.,verbose::Bool=false) where {DDT,RDT}
    Test=true
    TEST=Array{Bool,1}(undef,0)
    MYTOL=Array{Float64,1}(undef,0)
    DIF=Array{Float64,1}(undef,0)
    RER=Array{Float64,1}(undef,0)
    RTO=Array{Float64,1}(undef,0)
    for s=1:samples
        x= DDT<:Real ? jo_convert(DDT,randn(A.n)) : jo_convert(DDT,complex.(randn(A.n),randn(A.n)))
        y= DDT<:Real ? jo_convert(DDT,randn(A.n)) : jo_convert(DDT,complex.(randn(A.n),randn(A.n)))
        Axy=A*(x+y)
        AxAy=(A*x+A*y)
        dif=norm(Axy-AxAy); push!(DIF,dif)
        rer=abs(dif/norm(Axy)); push!(RER,rer)
        rto=abs(norm(AxAy)/norm(Axy)); push!(RTO,rto)
        mytol=(tol>0 ? tol : sqrt(max(eps(norm(Axy)),eps(norm(AxAy))))); push!(MYTOL,mytol)
        test=(dif < mytol); push!(TEST,test); Test=Test&&test
        result = test ? "PASSED" : "FAILED"
        if verbose println("Linear test [$s] $result with\n tol=    $mytol: \n diff=   $dif \n relerr= $rer \n ratio=  $rto") end
    end
    return Test,TEST,MYTOL,DIF,RER,RTO
end

# isadjoint(jo)
function isadjoint(A::joAbstractLinearOperator{DDT,RDT},samples::Integer=3;tol::Float64=0.,normfactor::Real=1.,userange::Bool=false,verbose::Bool=false) where {DDT,RDT}
    Test=true
    TEST=Array{Bool,1}(undef,0)
    MYTOL=Array{Float64,1}(undef,0)
    DIF=Array{Float64,1}(undef,0)
    RER=Array{Float64,1}(undef,0)
    RTO=Array{Float64,1}(undef,0)
    for s=1:samples
        if userange
            x= RDT<:Real ? jo_convert(RDT,randn(A.m)) : jo_convert(RDT,complex.(randn(A.m),randn(A.m)))
            y= RDT<:Real ? jo_convert(RDT,randn(A.m)) : jo_convert(RDT,complex.(randn(A.m),randn(A.m)))
            x=adjoint(A)*x
        else
            x= DDT<:Real ? jo_convert(DDT,randn(A.n)) : jo_convert(DDT,complex.(randn(A.n),randn(A.n)))
            y= DDT<:Real ? jo_convert(DDT,randn(A.n)) : jo_convert(DDT,complex.(randn(A.n),randn(A.n)))
            y=A*y
        end
        nfr=convert(RDT, normfactor)
        nfd=convert(DDT, normfactor)
        Axy=dot((A*x)/nfr,y)
        xAty=dot(x,(adjoint(A)*y)*nfd)
        dif=abs(xAty-Axy); push!(DIF,dif)
        rer=abs(dif/Axy); push!(RER,rer)
        rto=abs(xAty/Axy); push!(RTO,rto)
        mytol=(tol>0 ? tol : sqrt(max(eps(abs(Axy)),eps(abs(xAty))))); push!(MYTOL,mytol)
        test=(dif < mytol); push!(TEST,test); Test=Test&&test
        result = test ? "PASSED" : "FAILED"
        if verbose println("Adjoint test [$s] $result with\n tol=    $mytol: \n diff=   $dif \n relerr= $rer \n ratio=  $rto") end
    end
    return Test,TEST,MYTOL,DIF,RER,RTO
end


text:
struct BSONOID
    _wrap_::VoidPtr
    _ref_::Any

    function BSONOID()
        buffer = alloc_buf(12)
        ccall((:bson_oid_init, libbson), Cvoid, (BytePtr, VoidPtr), buffer, C_NULL)
        r = Base.unsafe_convert(BytePtr, buffer)
        new(r, buffer)
    end

    function BSONOID(str::AbstractString)
        cstr = string(str)

        isValid = ccall((:bson_oid_is_valid, libbson), Bool, (BytePtr, Csize_t),
                        cstr, length(cstr))
        isValid || error("'" * str * "': not a valid BSONOID string")

        buffer = alloc_buf(12)
        ccall((:bson_oid_init_from_string, libbson), Cvoid, (BytePtr, BytePtr), buffer, cstr)
        r = Base.unsafe_convert(BytePtr, buffer)
        new(r, buffer)
    end

    BSONOID(_ref_::Any) = new(pointer(_ref_), _ref_)
end

export BSONOID

Base.:(==)(lhs::BSONOID, rhs::BSONOID) =
    ccall((:bson_oid_equal, libbson), Bool, (VoidPtr, VoidPtr), lhs._wrap_, rhs._wrap_)

Base.hash(oid::BSONOID, h::UInt) =
    hash(ccall((:bson_oid_hash, libbson), UInt32, (BytePtr,), oid._wrap_), h)

function Base.convert(::Type{AbstractString}, oid::BSONOID)
    cstr = alloc_buf(25)
    ccall((:bson_oid_to_string, libbson), Cvoid, (BytePtr, BytePtr), oid._wrap_, cstr)
    String(unsafe_string(Base.unsafe_convert(BytePtr, cstr)))
end

Base.string(oid::BSONOID) = convert(AbstractString, oid)

Base.show(io::IO, oid::BSONOID) = print(io, "BSONOID($(convert(AbstractString, oid)))")


text:
@testset "new readme" begin

    mat = (1:4)' .+ rand(2,4)
    @cast rows[r][c] := mat[r,c]
    @cast cols[♜][🚣] := rows[🚣][♜]  # unicode 👍
    @reduce sum_r[c] := sum(r) mat[r,c]
    @test sum_r == sum(rows) # true

    B = rand(2*5, 3)
    @cast A[i,j,k] := B[(i,k),j]  i:2
    @test size(A) == (2,3,5)
    @cast A[i,j,k] = B[(i,k),j];

    imgs = [ rand(8,8) for i=1:16 ];

    @cast G[(i,I), (j,J)] := imgs[(I,J)][i,j] J:4
    @cast G[ i\I,   j\J ] = imgs[ I\J ][i,j] # in-place

    G = rand(16,32);

    @reduce H[a, b] := maximum(α,β)  G[α\a, β\b]  α:2,β:2
    @test size(G) == 2 .* size(H)

    @reduce H4[a, b] := maximum(α:4,β:4)  G[α\a, β\b]
    @test size(G) == 4 .* size(H4)

    W = randn(2,3,5,7);
    @cast Z[_,i,_,k] := W[2,k,4,i]  # equivalent to Z[1,i,1,k] on left
    @test size(Z) == (1,7,1,3)

    B = [ i .* ones(2,2) for i=1:4 ]
    @reduce A[μ,ν,J] := prod(i:2) B[(i,J)][μ,ν]
    @test size(A) == (2, 2, 2)

end
@testset "old readme" begin

    B = rand(3,4,5);
    @cast A[(i,j),k] := B[i,j,k]  # new matrix from tensor B

    B = rand(3*5,4);
    A = zeros(3,4,5);
    @cast A[i,j,k] = B[(i,k),j]   # write into an existing tensor A

    B = rand(3,4,5);
    @cast A[(i,j,k)] := B[i,j,k]  # reshaped view A = vec(B)

    B = [rand(3) for i=1:4];
    @cast A[i,j] := B[i][j]       # hcat a vector of vectors
    @test size(A) == (4,3)

    B = [rand(7) for i=1:3, k=1:4];
    A = zeros(3,7,4);
    @cast A[i,j,k] = B[i,k][j]    # write into A

    B = rand(2,3);
    @cast A[i][j] := B[j,i]       # create views A = collect(eachcol(B))

    B = [rand(3) for i=1:4];
    A = @cast [(i,j)] := B[j][i]  # vcat a vector of vectors
    @test size(A) == (12,)

    B = [rand(3,4,5,6) for i=1:7];
    A = @cast [(i,j),l][k,m] := B[i][j,k,l,m]; # glue then slice then reshape
    @test size(A) == (21,5)

    B = rand(2*5, 3);
    @cast A[i,j,k] := B[(i,k),j]  i:2  # could give (i:2, j:3, k:5)
    @test size(A) == (2,3,5)

    @cast A[i,j,k] := B[(i,k),j]  (i:2, j:3, k:5)
    @test size(A) == (2,3,5)

    @pretty @cast A[(i,j)] = B[i,j]

    @pretty @cast A[k][i,j] := B[i,(j,k)]  k:length(C)

    M = rand(3,4)
    @cast S[i][j] := M[i,j]       # S = julienne(M, (*,:)) creates views, S[i] == M[i,:]
    @cast Z[i,j] := S[i][j]       # Z = align(S, (*,:)) makes a copy
    @test size(Z) == (3,4)

    B = [rand(2,3) for k=1:4, l=1:5];
    @cast A[i,j,k,l] := B[k,l][i,j]
    @test size(A) == (2,3,4,5)

    using StaticArrays
    M = rand(1:99, 2,3)

    @cast S[k]{i} := M[i,k]  i:2  # S = reinterpret(SVector{2,Int}, vec(M)) needs the 2
    @cast N[k,i] := S[k]{i}       # such slices can be reinterpreted back again

    M[1,2]=42; N[2,1]==42          # all views of the original matrix
    @test N[2,1]==42

end
@testset "other random" begin

    A = rand(1:99, 3,3);
    X = randn(4)
    @cast B[i,j,k] := A[i,j] * X[k]
    @test B == A .* reshape(X, 1,1,:)

    G = rand(4, 4);
    two = 2
    @cast G2[k, l] := G[k,2] * G[l,2]
    @cast G3[k, l] := G[k,$two] * G[l,$two]
    @test G2 == G3 == G[:,2] .* G[:,2]'

end


text:
module Util

"""
create annotations for any struct
"""
macro namedargs(StructDec)
  sstruct = string(StructDec)
  typename = sstruct[findfirst(r"(?<=struct).*(?=\n)", sstruct)]
  typename = strip(typename)
  parametric = findfirst(r"{.*}", typename)
  if parametric == nothing
    parametric = ""
  else
    parametric = "where $(typename[parametric])"
  end
  typename = typename[findfirst(r"\w+", typename)]

  function fieldNameAndType(expr)
    return split(string(expr), "::")
  end

  args = StructDec.args[3].args[2:2:end]
  # args = filter(a -> !occursin(r"function\s+\w+", string(a)), args)
  args = filter(a -> !occursin(r"new\(+", string(a)), args)
  args = map(a -> fieldNameAndType(a), args)
  fnames = [a[1] for a in args]
  ftypes = [a[2] for a in args]

  callArgs = join(["$(a[1])" for a in args], ", ")
  oldConstructorCall = "$typename($callArgs)"

  fn_args = join(["$(a[1])::$(a[2])" for a in args], ", ")
  namedArgsConstructor = "$typename(; $fn_args) $parametric"

  return quote
    $(esc(StructDec))
    $(esc(Meta.parse("$namedArgsConstructor = $oldConstructorCall")))
  end
end

end # module Util

#= Test. TODO: Put in runtests.jl
import .Util.@namedargs

@namedargs struct Bab{T <: Number}
  x::Int
  y::T
end

@namedargs struct Bob
  x::Int
  y::String
end
beb = Bob(x=1, y="sos")

@namedargs struct Beb{S <: Number, T <: Int}
  x::S
  y::T
end
beb = Beb(x=1.0, y=2)
beb.x += 1 # should not work (immutable)

# Note that the first version is slower here.
@time for i in 1:100000000
  beb = Beb(x=1.0, y=2)
end

@time for i in 1:100000000
  beb = Beb(1.0, 2)
end

# But the timings here are the same
@time [Beb(x=1.0, y=2) for i in 1:10000];
@time [Beb(1.0, 2) for i in 1:10000];

@namedargs mutable struct Bib{S <: Number, T <: Int}
  x::S
  y::T
end
bib = Bib(x=1, y=2.0)
bib.x += 1 # this should work

@namedargs struct Bub{T <: Number}
  x::Int
  y::T
end;

=#


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

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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 65beba5e-776d-40d5-96a6-d0a2f3a9a4ea
#using Pkg; Pkg.activate("/Users/jj/JuliaProjects/LaserLab/")

# ╔═╡ 8a67d331-35ae-44d6-942c-a61db4e7afa4
begin
	using PlutoUI
	using CSV
	using DataFrames
	using Images
	using ImageBinarization
	using Colors
	using Plots
	using Printf
	using Interpolations
	using QuadGK
	using Markdown
	using InteractiveUtils
	using LsqFit
	using Statistics
	using StatsBase
	import Glob
end

# ╔═╡ d1b57350-bd9d-11ec-21f4-3d7a8656a31c
import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W,
    A, N, mol, mmol, V, L, M

# ╔═╡ 53943a33-9cf3-4c9f-add5-278fdc3bc0fe
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ b829dbf6-486e-4bd1-914c-54fe6d389b9c
lfi = ingredients("../src/LaserLab.jl")

# ╔═╡ 7aafa82d-4fe3-4bdc-a14c-ec48a7178260
PlutoUI.TableOfContents(title="G2 Bold73 analysis", indent=true)

# ╔═╡ 6926c7d6-0c1d-4054-9d02-3b12c21ad81c
collect(10:-1:1 )

# ╔═╡ 9da70d55-28b9-4c12-ba6a-c560832a5194
md"""
# Read Data
"""

# ╔═╡ 6502fdee-cd60-4d32-ac77-c42c845b12ba
begin
	proot ="/Users/jj/JuliaProjects/LaserLab/labdata"
	sroot = "/Users/jj/JuliaProjects/LaserLab/data/G2Sl"
	psample = "G2_BOLD_073_A1"
	xrun = ["G20406","G20407", "G2Ba0408"]
	prun = Dict("G20406"=>"G2SL_20220406",
		        "G20407"=>"G2SL_20220407",
		        "G2Ba0408"=>"G2SL-BaCl2_20220408")
	pfiles ="*.dat"
	fltn = collect(2:11)
	sflts = string.(fltn)
	sreps = string.(collect(0:10))

	md""" Select run type : $(@bind whichxr Select(xrun))"""
end

# ╔═╡ 21743b1b-3a48-4eb2-9dcd-3a8262c45dec
begin
	path = joinpath(proot,psample, prun[whichxr])
	fls = joinpath(path,pfiles)
	nfiles = length(fls)
	files = Glob.glob(pfiles, path);
	#names = [split(f,"/")[end] for f in files];

md"
- Path for glob: $path
- files : $fls
- number of files: $nfiles
"
end

# ╔═╡ a2ce79cf-a447-4cab-8843-a2cd691b4004
md"""
# Select Filter and Rep
"""

# ╔═╡ 90bda847-8afe-4968-9a02-50ad7c7cc572
md""" Select Filter number (2-11): $(@bind whichf Select(sflts))"""

# ╔═╡ 1cd87997-3ff3-40fa-98b7-ee9655848a76
md""" Select repetition number (0-10): $(@bind whichr Select(sreps))"""

# ╔═╡ 21f71c85-9e22-4834-a7d7-76e3d1ab681b
#imgm, drkm, sdrkm, simgm = get_images(files, whichf, whichr);

# ╔═╡ 0387a9f9-4e0b-4568-a698-5fe9e66e5c39
imgnt = lfi.LaserLab.get_images(files, whichf, whichr)

# ╔═╡ 5c283457-fc1d-4249-bd4e-029204e5c57b
lfi.LaserLab.plot_images(imgnt, whichf, whichr);

# ╔═╡ bfc98395-7d39-4589-b007-2f8fbbedeaf8
imgmxnt = lfi.LaserLab.signal_around_maximum(imgnt.cimage, imgnt.cdark; nsigma=3);

# ╔═╡ 60b420bd-9bfd-4361-bec9-11a1ccd205a5
begin
	imgmn = imgnt.cimage ./maximum(imgnt.cimage)
	gimg = Gray.(imgmn);
	imgmxn = imgmxnt.img ./maximum(imgmxnt.img)
	gimgx = Gray.(imgmxn);
	mosaicview(gimg, gimgx; nrow = 1)
end

# ╔═╡ bd3beecc-2154-4c70-a8b4-7334b8cc643d
gimgx

# ╔═╡ 1bfdfd1a-9435-420d-be62-0df760b4bc06
size(imgmxnt.img)

# ╔═╡ 20902f12-03c0-4a98-937a-a9ff8433fb1b
stot = sum(imgmxnt.img)

# ╔═╡ b776cef6-90d6-4897-bd6e-b81a3a58b416
md"""
- Filter = $whichf
- repetition = $whichr
- **signal**
- value of max = $(imgmxnt.max)
- position of max: i = $(imgmxnt.imax)  j = $(imgmxnt.jmax)
- Total light in spot (DC subtracted) = $(round(stot, sigdigits=3))
"""

# ╔═╡ 01aeef0d-12e1-4e32-974f-f5405a661d95
md"""
# Reconstruct spectra
"""

# ╔═╡ 8feb58fd-3cb7-4860-bc22-52bcee93d40a
scs = [string(c) for c in collect(5:10)]

# ╔═╡ 18a0121f-b0a5-46e5-8178-9bc800476b65
@bind zrepx CheckBox()

# ╔═╡ f6c0d5d7-bfdd-49fd-97da-786d49e53211
begin
	fdf = lfi.LaserLab.load_df_from_csv(sroot, "G20406.csv", lfi.LaserLab.enG)
	fdf1 = lfi.LaserLab.load_df_from_csv(sroot, "G20407.csv", lfi.LaserLab.enG)
	fdf2 = lfi.LaserLab.load_df_from_csv(sroot, "G2Ba0408.csv", lfi.LaserLab.enG)
	xn1 = names(fdf1)
end

# ╔═╡ 3aa1788a-131a-402f-9493-32d353020219
if whichxr == "G20406"
	plot(collect(1:10), fdf.r1, lw=2, label="rep1 G20406")
elseif whichxr == "G20407"
	plot(collect(1:10), fdf1.r5, lw=2, label=string(xn1[1],whichxr))
	plot!(collect(1:10), fdf1.r6, lw=2, label=string(xn1[2],whichxr))
elseif whichxr == "G2Ba0408"
	plot(collect(1:10), fdf2.r5, lw=2, label=string(xn1[1],whichxr))
	plot!(collect(1:10), fdf2.r6, lw=2, label=string(xn1[2],whichxr))
end

# ╔═╡ 8740f867-109f-47c1-9e90-57ece3b83886
md"""
# Comparisons
"""

# ╔═╡ 2d8b888a-9cf8-4c53-9f6b-d61f2d2b7514
@bind zcomp CheckBox()

# ╔═╡ 41cf396d-01ba-4e09-9ac9-bdf0a713de9e
if zcomp
	plot(collect(1:10), fdf.r1, lw=2, label="rep1 G20406")
	plot!(collect(1:10), fdf1.r5, lw=2, label=string(xn1[1],whichxr))
	plot!(collect(1:10), fdf1.r6, lw=2, label=string(xn1[2],whichxr))
	plot!(collect(1:10), fdf2.r5, lw=2, label=string(xn1[1],whichxr))
	plot!(collect(1:10), fdf2.r6, lw=2, label=string(xn1[2],whichxr))

end

# ╔═╡ 7eaaa7a9-2a5b-4ed6-b3b8-1a03e0fb8081
if zcomp
	pngn = joinpath(sroot, "G2comp.png")
	png(pngn)
	end

# ╔═╡ cec73687-de04-46e6-a73e-a4d1d3beafc5
md"# Functions"

# ╔═╡ 3f08fd90-ac06-4c29-a66c-cdb782f0717b
function sumflt(files, fltn; rep="0", nsigma=3)
	STOT = Vector{Float64}()
	for iflt in sflts
		imgnt = lfi.LaserLab.get_images(files, iflt, rep)
		imgmxnt = lfi.LaserLab.signal_around_maximum(imgnt.cimage, imgnt.cdark; nsigma=nsigma)
		stot = sum(imgmxnt.img)
		push!(STOT,stot)
	end
	STOT
end

# ╔═╡ 2c71142a-34b2-4804-9af7-dc8147ab30a3
if zrepx
	if whichxr == "G20406"
		sts = sumflt(files, fltn; rep="0", nsigma=3)
		df = DataFrame("r1" => sts)
		dfn = joinpath(sroot, "G20406.csv")
		CSV.write(dfn, df)
		plot(collect(1:10), sts, lw=2, label="rep0 G20406")

	elseif whichxr == "G20407"
		df1 = DataFrame()
		for rep in scs
			sts1 = sumflt(files, fltn; rep=rep, nsigma=3)
			lbl = string("r",rep)
			df1[!,lbl] = sts1
		end
		dfn1 = joinpath(sroot, "G20407.csv")
		CSV.write(dfn1, df1)

		xn = names(df1)
		plot(collect(1:10), df1.r5, lw=2, label=string(xn[1],whichxr))
		plot!(collect(1:10), df1.r6, lw=2, label=string(xn[2],whichxr))
	elseif whichxr == "G2Ba0408"
		df2 = DataFrame()
		for rp in scs
			sts2 = sumflt(files, fltn; rep=rp, nsigma=3)
			lbl = string("r",rp)
			df2[!,lbl] = sts2
		end
		dfn2 = joinpath(sroot, "G2Ba0408.csv")
		CSV.write(dfn2, df2)
		xn2 = names(df2)
		plot(collect(1:10), df2.r5, lw=2, label=string(xn2[1],whichxr))
		plot!(collect(1:10), df2.r6, lw=2, label=string(xn2[2],whichxr))
	else
		@info "not implemented"
	end

end

# ╔═╡ 6fdaa0d2-b6af-40ca-b8e2-5cb8176bbadb
function select_image(inames::Vector{String}, fnumber::String, repetition::String,
	                  itype::String="Imag", inum::String="1")

	names = [split(f,"/")[end] for f in inames]
	if repetition == "0"
		fhead = string("Filter_", fnumber, "_ExpoT_1000000ms_")
	else
		fhead = string("Filter_", fnumber, "_rep_",repetition,"_ExpoT_1000000ms_")
	end

	if itype == "Dark"
		ftail = string(itype, "_", inum, ".dat")
	elseif itype == "Imag" && inum=="1"
		ftail = string("Imag_1.dat")
	elseif itype == "Imag" && inum=="c"
		ftail = string("Imag_corrected.dat")
	else
		@info "not implemented"
	end

	ff = string(fhead, ftail)
	indx = findall([occursin(ff, name) for name in names])[1]
	inames[indx]
end

# ╔═╡ edabab28-72bb-43f2-b90b-7b6ba939f2eb
function dftomat(img::DataFrame)
	mimg = Array{Float64}(undef, size(img)...);
	for i in 1:size(img)[1]
		for j in 1:size(img)[2]
			if img[i,j] >= 0.0 && !isnan(img[i,j])
    			mimg[i, j] = img[i,j]
			else
				mimg[i, j] = 0.0
			end
		end
	end
	mimg
end

# ╔═╡ b8af0d87-f9e1-4ea7-9fa0-1c675233a8d3
function findpeaks(img::Matrix{Float64}; xrange::Int64=10)
	maxx, indxt = findmax(img)
	ixmin = max(1, indxt[1] - xrange)
	ixmax = min(size(img)[1], indxt[1] + xrange -1)
	iymin = max(1, indxt[2] - xrange)
	iymax = min(size(img)[2], indxt[2] + xrange -1)

	img2 = Array{Float64}(undef, ixmax-ixmin+1, iymax-iymin+1);

	stot = 0.0
	for (i, ix) in enumerate(ixmin:ixmax)
		for (j, iy) in enumerate(iymin:iymax)
			img2[i,j]=img[ix,iy]
			stot += img[ix,iy]
		end
	end
	img2, stot
end

# ╔═╡ 5dc25eb7-7841-4c04-aa5e-7d63acc5a312
function get_images(files::Vector{String}, whichf::String, whichr::String)
	ffimg = select_image(files, whichf, whichr, "Imag", "1");
	imgdf = DataFrame(CSV.File(ffimg, header=false,delim="\t"));
	imgm  = dftomat(imgdf)

	ffdrk = select_image(files, whichf, whichr, "Dark", "1");
	drkdf = DataFrame(CSV.File(ffdrk, header=false,delim="\t"));
	drkm  = dftomat(drkdf)

	ffdrk2 = select_image(files, whichf, whichr, "Dark", "2");
	drkdf2 = DataFrame(CSV.File(ffdrk2, header=false,delim="\t"));
	sdrk   = drkdf .- drkdf2
	sdrkm  = dftomat(sdrk)

	simgdf = imgdf .- drkdf
	simgm  = dftomat(simgdf)

	imgm, drkm, sdrkm, simgm
end

# ╔═╡ f41cf5c7-d208-40f9-a4ee-5e992f73ca36
function recospec(files::Vector{String}; rep::String="5", xrange::Int64=30)
	STOT = []
	MAX  = []
	IX   = []
	IY   = []
	STOTd = []
	MAXd  = []
	IXd  = []
	IYd   = []
	for flt in string.(collect(2:11))
		imgm, drkm, sdrkm, simgm = get_images(files, flt, rep)
		pimgx, stot = findpeaks(simgm; xrange=xrange)
		maxx, indxt = findmax(pimgx)
		pimgd, stotd = findpeaks(sdrkm; xrange=xrange)
		maxd, indxd = findmax(pimgd)
		push!(STOT, stot)
		push!(MAX, maxx)
		push!(IX, indxt[1])
		push!(IY, indxt[2])
		push!(STOTd, stotd)
		push!(MAXd, maxd)
		push!(IXd, indxd[1])
		push!(IYd, indxd[2])
	end
	Dict("stot" =>STOT, "smax"=>MAX, "six"=>IX,"siy"=>IY,
	     "dtot" =>STOTd, "dmax"=>MAXd, "dix"=>IXd,"diy"=>IYd)
end

# ╔═╡ 3794ff0b-87bf-4b08-9b33-4016280ebbbf
function plot_images2(imgnt)
	imghm = heatmap(imgnt.image, title=string("rep=", whichr, " filt=", whichf, " Img" ), titlefontsize=10)

	drkhm = heatmap(imgnt.dark, title=string("rep=", whichr, " filt=", whichf, " Dark" ), titlefontsize=10)

	simghm = heatmap(imgnt.corr, title=string("rep=", whichr, " filt=", whichf, "Img - Dark" ))

	plot(imghm, drkhm, simghm, layout=(3,1), titlefontsize=10)
end

# ╔═╡ fca14add-ed54-449b-82a3-bffd03f88cdd
function plot_images(imgm, drkm, sdrkm, simgm)
	imghm = heatmap(imgm, fill_z=imgm, title=string("rep=", whichr, " filt=", whichf, " Img" ), titlefontsize=10)

	drkhm = heatmap(drkm, title=string("rep=", whichr, " filt=", whichf, " Dark" ), titlefontsize=10)

	sdrkhm = heatmap(sdrkm, title=string("rep=", whichr, " filt=", whichf, "Dark - dark2" ), titlefontsize=10)

	simghm = heatmap(simgm, title=string("rep=", whichr, " filt=", whichf, "Img - Dark" ))

	plot(imghm, drkhm, sdrkhm, simghm, layout=(2,2), titlefontsize=10)
end

# ╔═╡ 4480a6a1-3299-43b9-b871-0b9f525577e5
function img_stats(img)
	xmean = describe(img).mean
	xmax = describe(img).max
	xmin = describe(img).min
	hxmean = lfi.LaserLab.hist1d(xmean, 50, minimum(xmean), maximum(xmean))
	hxmin = lfi.LaserLab.hist1d(xmin, 50, minimum(xmin), maximum(xmin))
	hxmax = lfi.LaserLab.hist1d(xmax, 50, minimum(xmax), maximum(xmax))
	(xmean = xmean, xmin = xmin, xmax = xmax,
	hxmean = hxmean, hxmin = hxmin, hxmax = hxmax)
end

# ╔═╡ 454cab5c-e715-46bd-a799-0e9e2e9a7f67
function img_title(label)
	string(label,"_Filter_", whichf, "_rep_", whichr, "_type_", whicht, "_number_", whichn)
end

# ╔═╡ 4c1c878c-fcfb-4069-a27a-5b1288896955
with_terminal() do
	for (i, ix) in enumerate(10:12)
		for (j, iy) in enumerate(20:22)
			#println("i = ", i, " j = ",j ," ix = ", ix, " iy = ", iy)
		end
	end
end

# ╔═╡ Cell order:
# ╠═8a67d331-35ae-44d6-942c-a61db4e7afa4
# ╠═65beba5e-776d-40d5-96a6-d0a2f3a9a4ea
# ╠═d1b57350-bd9d-11ec-21f4-3d7a8656a31c
# ╠═53943a33-9cf3-4c9f-add5-278fdc3bc0fe
# ╠═b829dbf6-486e-4bd1-914c-54fe6d389b9c
# ╠═7aafa82d-4fe3-4bdc-a14c-ec48a7178260
# ╠═6926c7d6-0c1d-4054-9d02-3b12c21ad81c
# ╠═9da70d55-28b9-4c12-ba6a-c560832a5194
# ╠═6502fdee-cd60-4d32-ac77-c42c845b12ba
# ╠═21743b1b-3a48-4eb2-9dcd-3a8262c45dec
# ╠═a2ce79cf-a447-4cab-8843-a2cd691b4004
# ╠═90bda847-8afe-4968-9a02-50ad7c7cc572
# ╠═1cd87997-3ff3-40fa-98b7-ee9655848a76
# ╠═21f71c85-9e22-4834-a7d7-76e3d1ab681b
# ╠═0387a9f9-4e0b-4568-a698-5fe9e66e5c39
# ╠═5c283457-fc1d-4249-bd4e-029204e5c57b
# ╠═bfc98395-7d39-4589-b007-2f8fbbedeaf8
# ╠═60b420bd-9bfd-4361-bec9-11a1ccd205a5
# ╠═bd3beecc-2154-4c70-a8b4-7334b8cc643d
# ╠═1bfdfd1a-9435-420d-be62-0df760b4bc06
# ╠═20902f12-03c0-4a98-937a-a9ff8433fb1b
# ╠═b776cef6-90d6-4897-bd6e-b81a3a58b416
# ╠═01aeef0d-12e1-4e32-974f-f5405a661d95
# ╠═8feb58fd-3cb7-4860-bc22-52bcee93d40a
# ╠═18a0121f-b0a5-46e5-8178-9bc800476b65
# ╠═2c71142a-34b2-4804-9af7-dc8147ab30a3
# ╠═f6c0d5d7-bfdd-49fd-97da-786d49e53211
# ╠═3aa1788a-131a-402f-9493-32d353020219
# ╠═8740f867-109f-47c1-9e90-57ece3b83886
# ╠═2d8b888a-9cf8-4c53-9f6b-d61f2d2b7514
# ╠═41cf396d-01ba-4e09-9ac9-bdf0a713de9e
# ╠═7eaaa7a9-2a5b-4ed6-b3b8-1a03e0fb8081
# ╠═cec73687-de04-46e6-a73e-a4d1d3beafc5
# ╠═3f08fd90-ac06-4c29-a66c-cdb782f0717b
# ╠═6fdaa0d2-b6af-40ca-b8e2-5cb8176bbadb
# ╠═edabab28-72bb-43f2-b90b-7b6ba939f2eb
# ╠═b8af0d87-f9e1-4ea7-9fa0-1c675233a8d3
# ╠═f41cf5c7-d208-40f9-a4ee-5e992f73ca36
# ╠═5dc25eb7-7841-4c04-aa5e-7d63acc5a312
# ╠═3794ff0b-87bf-4b08-9b33-4016280ebbbf
# ╠═fca14add-ed54-449b-82a3-bffd03f88cdd
# ╠═4480a6a1-3299-43b9-b871-0b9f525577e5
# ╠═454cab5c-e715-46bd-a799-0e9e2e9a7f67
# ╠═4c1c878c-fcfb-4069-a27a-5b1288896955

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
using Plots

function bracket_minimum(f, x=0, s=1e-2, k=2.0)
    a, ya = x, f(x);
    b, yb = a + s, f(a + s);

    if yb > ya
        a, b = b, a;
        ya, yb = yb, ya;
        s = -s;
    end

    while true
        c, yc = b + s, f(b + s);
        if yc > yb
            return a < c ? (a, c) : (c, a);
        end

        a, ya, b, yb = b, yb, c, yc;
        s *= k;
    end
end

# test

f1 = x -> (1 .+ x).^2;
f2 = x -> (1 + x)^2;
x = -3:0.1:3;
y = f1(x);

a, b = bracket_minimum(f2, -2, 1e-2, 1.01);

plot(x, y)
vline!([a, b])


text:
# developed with Julia 1.1.1
#
# function to initialize parallel computing

function init_parallel(n::Integer = Sys.CPU_THREADS; kw...)
    if n < 1
        error("number of workers must be greater than 0")
    elseif n == 1 && workers() != [1]
        rmprocs(workers())
    elseif n > nworkers()
        p = addprocs(n - (nprocs() == 1 ? 0 : nworkers()); kw...)
    elseif n < nworkers()
        rmprocs(workers()[n + 1:end])
    end
    @everywhere @eval using EMSx
    return workers()
end


text:
using Base.Test
using ComparativeAutograder

return3() = 3
result = runtestcase(return3, TestCase(()))
println(result)
@test result.result == 3
@test result.exception == nothing
@test result.backtrace == nothing
@test result.stdout == ""
@test result.stderr == ""
@test result.time > 0

return3bad() = nonexistantvariable + 8
result = runtestcase(return3bad, TestCase(()))
println(result)
@test result.result == nothing
@test contains(result.exception, "UndefVarError")
@test result.backtrace != nothing

returnsum(x, y) = x + y
result = runtestcase(returnsum, TestCase((3, 4)))
@test result.result == 7


text:
# This example shows how to use RAFF.jl to detect circles drawn from
# the user. This example was inspired in `GtkReactive.jl` drawing
# example.

# Packages needed
#
# add Gtk
# add GtkReactive
# add Graphics
# add Colors

using Gtk, Gtk.ShortNames, GtkReactive, Graphics, Colors
using RAFF

win = Window("Drawing") |> (bx = Box(:v))

set_gtk_property!(bx, :spacing, 10)

cb = GtkReactive.dropdown(["RAFF", "Least Squares"])
push!(bx, cb)

c = canvas(UserUnit, 200, 200)       # create a canvas with user-specified coordinates
push!(bx, c)

const moddraw = Signal([]) # the model points
const newline = Signal([]) # the in-progress line (will be added to list above)

const drawing = Signal(false)  # this will become true if we're actively dragging

choice = map(x->(push!(moddraw, []);value(x)), cb)

function run_raff()

    n, model, = RAFF.model_list["circle"]

    np, = size(value(newline))

    data = Array{Float64, 2}(undef, np, 3)

    l = value(newline)

    for i in 1:np

        data[i,1] = l[i].x
        data[i,2] = l[i].y
        data[i,3] = 0.0

    end

    r = if value(choice) == "RAFF"

        raff(model, data, n, ftrusted=0.7, MAXMS=10)

    else

        raff(model, data, n, ftrusted=1.0, MAXMS=10)

    end
    println(r)
    # Build points for drawing the circle
    p = (α) -> [abs(r.solution[3]) * cos(α) + r.solution[1], abs(r.solution[3]) * sin(α) + r.solution[2]]

    push!(moddraw, [p(i) for i in LinRange(0.0, 2 * π, 100)])

end

# If has changed the selection box, then run RAFF again.
sigc = map(choice) do c

    run_raff()

end

sigstart = map(c.mouse.buttonpress) do btn
    # This is the beginning of the function body, operating on the argument `btn`
    if btn.button == 1 && btn.modifiers == 0 # is it the left button, and no shift/ctrl/alt keys pressed?
        push!(drawing, true)   # activate dragging
        push!(newline, [btn.position])  # initialize the line with the current position
        push!(moddraw, [])
    end
end

const dummybutton = MouseButton{UserUnit}()
# See the Reactive.jl documentation for `filterwhen`
sigextend = map(filterwhen(drawing, dummybutton, c.mouse.motion)) do btn
    # while dragging, extend `newline` with the most recent point
    push!(newline, push!(value(newline), btn.position))
end

sigend = map(c.mouse.buttonrelease) do btn
    if btn.button == 1
        push!(drawing, false)  # deactivate dragging
        run_raff()
    end
end

function drawline(ctx, l, color)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p.x, p.y)
    set_source(ctx, color)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p.x, p.y)
    end
    stroke(ctx)
end

function drawcircle(ctx, l, color)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p[1], p[2])
    set_source(ctx, color)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p[1], p[2])
    end
    stroke(ctx)
end

# Because `draw` isn't a one-line function, we again use do-block syntax:
redraw = draw(c, moddraw, newline) do cnvs, circ, newl  # the function body takes 3 arguments
    fill!(cnvs, colorant"white")   # set the background to white
    set_coordinates(cnvs, BoundingBox(0, 1, 0, 1))  # set coordinates to 0..1 along each axis
    ctx = getgc(cnvs)   # gets the "graphics context" object (see Cairo/Gtk)
    drawcircle(ctx, circ, colorant"blue")  # draw old lines in blue
    drawline(ctx, newl, colorant"red")    # draw new line in red
end

showall(win)

#If we are not in a REPL
if (!isinteractive())

    # Create a condition object
    c = Condition()

    # Get the window
    # win = guidict["gui"]["window"]

    # Notify the condition object when the window closes
    signal_connect(win, :destroy) do widget
        notify(c)
    end

    # Wait for the notification before proceeding ...
    wait(c)
end


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

include("block.jl")
include("inline.jl")

@flavor common [list, indentcode, blockquote, admonition, footnote, hashheader, horizontalrule,
                paragraph,

                linebreak, escapes, inline_code,
                asterisk_bold, underscore_bold, asterisk_italic, underscore_italic, image, footnote_link, link, autolink]

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
if false
M=Model(:[1-R*β*Expect(λ[+1])/λ
           Uh/λ-η
],:[b       = (-2,10.,8)
    η       = (1,0.9,0.1,1)
],:[b       = (-2,10.,b*0.95)
    h       = (0,1,0.7)
    c       = h*η+R*b[-1]-b
    λ       = c^-σc
    Uh      = ϕh*(1-h)^-σh
    B       = ∫(b,0.0)
    H       = ∫(h*η,0.3)
],:[β       = 0.98
    σc      = 2.5
    ϕh      = 2.0
    σh      = 2.0
    R       = 1.0166])

HADSGE.initdisplay(M)
solve(M,disp=50,crit=1e-12)

end

function initdisplay(M::Model)
    W = Window()
    Blink.addstyle!(W,"""
    body {background-color:#2d2d2d;font-size: 10;color: white}
    .tooltip {
        position: relative;display: inline-block;border-bottom: 1px dotted black;}
    .tooltip .tooltiptext {visibility: hidden;width: 120px;background-color: black;color: #fff;text-align:center;padding: 5px 0;border-radius:6px;position:absolute;z-index:1;}
    .tooltip:hover .tooltiptext {visibility: visible;}
    .field-tip {position:relative;cursor:default}
    .field-tip .tip-content {position:absolute;top:-10px;right:9999px;width:200px;margin-right:-220px;padding:10px;color:#fff;background:#333;opacity:0;
        -webkit-transition:opacity 250ms ease-out;
           -moz-transition:opacity 250ms ease-out;
            -ms-transition:opacity 250ms ease-out;
             -o-transition:opacity 250ms ease-out;
                transition:opacity 250ms ease-out;}
    .field-tip .tip-content:before {
        content:' ';position:absolute;top:50%;
        left:-16px;width:0;height:0;margin-top:-8px;border:8px solid transparent;border-right-color:#333;}
    .field-tip:hover .tip-content {right:-20px;opacity:1;}
    """)
    @js W resizeTo(495,400)
    @js W moveTo(866,5)

    # Summary
    io=IOBuffer()
    for v in State(M.variables)
        if isa(v,Stochastic)
            if isa(v,AR)
                desc="[$(v.name) = $(round(1-v.ρ*v.μ,4)) + $(v.ρ)*$(v.name)[-1] + $(v.σ)*ϵ]  $(v.bounds)"
            elseif isa(v,Markov)
                desc="[$(v.name) = T($(v.name)[-1])] $(v.bounds)"
            end
        else
            desc="Endogenous: $(v.bounds)"
        end
        desc*= " $(length(v.x))"

        print(io,"""
            $(v.name)
            $desc
          """)
    end
    states = takebuf_string(io)
    close(io)

    io=IOBuffer()
    for v in Policy(M.variables)
        desc= " $(map(x->round(x,3),extrema(M[v.name,0])))"
        print(io,"""
            $(v.name)
            $desc
          """)
    end
    policies = takebuf_string(io)
    close(io)

    summary = """

            States
            $states

            Policy
            $policies

    """

    #Equations

    io = IOBuffer()
    print(io,"""
""")
    for i in 1:length(M.summary.equations.args)
        print(io,"""

        $i

        $(M.summary.equations.args[i])

                $(signif(median(M.Fval[:,i]),2))  ($(signif(maximum(M.Fval[:,i]),2)))
""")
    end
    print(io,"
")
    eqs = takebuf_string(io)
    close(io)

    # Progress meter

    progress= """

 Solve Step: <text id="solvestep"></text>

Time: <text id="Telapsed"></text>

Iterations: <text id="iiter">0</text>/<text id="Niter">0</text>

    """

    # progress ="""
    #
Progress (<text id="solvestep"></text>)

    #

    #

    #

    # """
    #

    #

    #     tid
    #     1
    #     2
    #     3
    #     4
    #     5
    #     6
    #     7
    #     8
    #

    #

    #
    #     0
    #     0
    #     0
    #     0
    #     0
    #     0
    #     0
    #     0
    #

    #

    # """

    body!(W,"""

Summary

    $eqs
    $summary
    $progress
    """),
    push!(M.summary.displays,W)
end

function displaysolve(M,iiter,Niter,MaxError,SError,t)
    if !active(M.summary.displays[end])
        initdisplay(M)
    else
        # @js M.summary.displays[end] document.getElementById("solvestep").innerHTML="Equations"
        @js M.summary.displays[end] document.getElementById("iiter").innerHTML=$iiter
        @js M.summary.displays[end] document.getElementById("Niter").innerHTML=$Niter
        @js M.summary.displays[end] document.getElementById("Telapsed").innerHTML=$t
        for i = 1:length(M.summary.equations.args)
            s = """$(signif(MaxError[i],2))  ($(signif(SError[i],2)))"""
            targ = "fval$i"
            @js M.summary.displays[end] document.getElementById($targ).innerHTML=$s
        end
    end
end

# function updatedisplay(M::Model,step,i,n,crit,t)
#     if !active(M.summary.displays[end])
#         initdisplay(M)
#     else
#         @js M.summary.displays[end] document.getElementById("solvestep").innerHTML=$step
#         iter = "$i/$n"
#         @js M.summary.displays[end] document.getElementById("iter").innerHTML=$iter
#         time = "$(signif(t,2))s"
#         @js M.summary.displays[end] document.getElementById("time").innerHTML=$time
#         for i = 1:length(M.summary.equations.args)
#             s = """$(signif(median(M.Fval[:,i]),2))  ($(signif(maximum(M.Fval[:,i]),2)))"""
#             targ = "fval$i"
#             @js M.summary.displays[end] document.getElementById($targ).innerHTML=$s
#         end
#     end
# end

#
# for i = 1:8
#     pn="progress$i"
# @js W document.getElementById($pn).innerHTML=$(i*10)
# end
#
# @js W document.getElementById("solvestep").innerHTML="a"
# @js M.summary.displays[end] document.getElementById("fval1").innerHTML

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE