text:
using FiniteDifferences: grad, jacobian, _jvp, _j′vp, jvp, j′vp, to_vec

# Dummy type where length(x::DummyType) ≠ length(first(to_vec(x)))
struct DummyType{TX<:Matrix}
    X::TX
end

function FiniteDifferences.to_vec(x::DummyType)
    x_vec, back = to_vec(x.X)
    return x_vec, x_vec -> DummyType(back(x_vec))
end

Base.:(==)(x::DummyType, y::DummyType) = x.X == y.X
Base.length(x::DummyType) = size(x.X, 1)

@testset "grad" begin

    @testset "grad(::$T)" for T in (Float64, ComplexF64)
        rng, fdm = MersenneTwister(123456), central_fdm(5, 1)
        x = randn(rng, T, 2)
        xc = copy(x)
        @test grad(fdm, x->sin(x[1]) + cos(x[2]), x)[1] ≈ [cos(x[1]), -sin(x[2])]
        @test xc == x
    end

    function check_jac_and_jvp_and_j′vp(fdm, f, ȳ, x, ẋ, J_exact)
        xc = copy(x)
        @test jacobian(fdm, f, x; len=length(ȳ))[1] ≈ J_exact
        @test jacobian(fdm, f, x)[1] == jacobian(fdm, f, x; len=length(ȳ))[1]
        @test _jvp(fdm, f, x, ẋ) ≈ J_exact * ẋ
        @test _j′vp(fdm, f, ȳ, x) ≈ transpose(J_exact) * ȳ
        @test xc == x
    end

    @testset "jacobian / _jvp / _j′vp (::$T)" for T in (Float64, ComplexF64)
        rng, P, Q, fdm = MersenneTwister(123456), 3, 2, central_fdm(5, 1)
        ȳ, A, x, ẋ = randn(rng, T, P), randn(rng, T, P, Q), randn(rng, T, Q), randn(rng, T, Q)
        Ac = copy(A)

        check_jac_and_jvp_and_j′vp(fdm, x->A * x, ȳ, x, ẋ, A)
        @test Ac == A
        check_jac_and_jvp_and_j′vp(fdm, x->sin.(A * x), ȳ, x, ẋ, cos.(A * x) .* A)
        @test Ac == A
    end

    @testset "multi vars jacobian/grad" begin
        rng, fdm = MersenneTwister(123456), central_fdm(5, 1)

        f1(x, y) = x * y + x
        f2(x, y) = sum(x * y + x)
        f3(x::Tuple) = sum(x[1]) + x[2]
        f4(d::Dict) = sum(d[:x]) + d[:y]

        @testset "jacobian" begin
            @testset "check multiple matrices" begin
                x, y = rand(rng, 3, 3), rand(rng, 3, 3)
                jac_xs = jacobian(fdm, f1, x, y)
                @test jac_xs[1] ≈ jacobian(fdm, x->f1(x, y), x)[1]
                @test jac_xs[2] ≈ jacobian(fdm, y->f1(x, y), y)[1]
            end

            @testset "check mixed scalar and matrices" begin
                x, y = rand(3, 3), 2
                jac_xs = jacobian(fdm, f1, x, y)
                @test jac_xs[1] ≈ jacobian(fdm, x->f1(x, y), x)[1]
                @test jac_xs[2] ≈ jacobian(fdm, y->f1(x, y), y)[1]
            end
        end

        @testset "grad" begin
            @testset "check multiple matrices" begin
                x, y = rand(rng, 3, 3), rand(rng, 3, 3)
                dxs = grad(fdm, f2, x, y)
                @test dxs[1] ≈ grad(fdm, x->f2(x, y), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f2(x, y), y)[1]
            end

            @testset "check mixed scalar & matrices" begin
                x, y = rand(rng, 3, 3), 2
                dxs = grad(fdm, f2, x, y)
                @test dxs[1] ≈ grad(fdm, x->f2(x, y), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f2(x, y), y)[1]
            end

            @testset "check tuple" begin
                x, y = rand(rng, 3, 3), 2
                dxs = grad(fdm, f3, (x, y))[1]
                @test dxs[1] ≈ grad(fdm, x->f3((x, y)), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f3((x, y)), y)[1]
            end

            @testset "check dict" begin
                x, y = rand(rng, 3, 3), 2
                d = Dict(:x=>x, :y=>y)
                dxs = grad(fdm, f4, d)[1]
                @test dxs[:x] ≈ grad(fdm, x->f3((x, y)), x)[1]
                @test dxs[:y] ≈ grad(fdm, y->f3((x, y)), y)[1]
            end
        end
    end

    function test_to_vec(x)
        x_vec, back = to_vec(x)
        @test x_vec isa Vector
        @test x == back(x_vec)
        return nothing
    end

    @testset "to_vec(::$T)" for T in (Float64, ComplexF64)
        if T == Float64
            test_to_vec(1.0)
            test_to_vec(1)
        else
            test_to_vec(.7 + .8im)
            test_to_vec(1 + 2im)
        end
        test_to_vec(randn(T, 3))
        test_to_vec(randn(T, 5, 11))
        test_to_vec(randn(T, 13, 17, 19))
        test_to_vec(randn(T, 13, 0, 19))
        test_to_vec([1.0, randn(T, 2), randn(T, 1), 2.0])
        test_to_vec([randn(T, 5, 4, 3), (5, 4, 3), 2.0])
        test_to_vec(reshape([1.0, randn(T, 5, 4, 3), randn(T, 4, 3), 2.0], 2, 2))
        test_to_vec(UpperTriangular(randn(T, 13, 13)))
        test_to_vec(Symmetric(randn(T, 11, 11)))
        test_to_vec(Diagonal(randn(T, 7)))
        test_to_vec(DummyType(randn(T, 2, 9)))

        @testset "$Op" for Op in (Adjoint, Transpose)
            test_to_vec(Op(randn(T, 4, 4)))
            test_to_vec(Op(randn(T, 6)))
            test_to_vec(Op(randn(T, 2, 5)))
        end

        @testset "Tuples" begin
            test_to_vec((5, 4))
            test_to_vec((5, randn(T, 5)))
            test_to_vec((randn(T, 4), randn(T, 4, 3, 2), 1))
            test_to_vec((5, randn(T, 4, 3, 2), UpperTriangular(randn(T, 4, 4)), 2.5))
            test_to_vec(((6, 5), 3, randn(T, 3, 2, 0, 1)))
            test_to_vec((DummyType(randn(T, 2, 7)), DummyType(randn(T, 3, 9))))
            test_to_vec((DummyType(randn(T, 3, 2)), randn(T, 11, 8)))
        end
        @testset "Dictionary" begin
            if T == Float64
                test_to_vec(Dict(:a=>5, :b=>randn(10, 11), :c=>(5, 4, 3)))
            else
                test_to_vec(Dict(:a=>3 + 2im, :b=>randn(T, 10, 11), :c=>(5+im, 2-im, 1+im)))
            end
        end
    end

    @testset "jvp(::$T)" for T in (Float64, ComplexF64)
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, T, N), randn(rng, T, M)
        ẋ, ẏ = randn(rng, T, N), randn(rng, T, M)
        xy, ẋẏ = vcat(x, y), vcat(ẋ, ẏ)
        ż_manual = _jvp(fdm, (xy)->sum(sin, xy), xy, ẋẏ)[1]
        ż_auto = jvp(fdm, x->sum(sin, x[1]) + sum(sin, x[2]), ((x, y), (ẋ, ẏ)))
        ż_multi = jvp(fdm, (x, y)->sum(sin, x) + sum(sin, y), (x, ẋ), (y, ẏ))
        @test ż_manual ≈ ż_auto
        @test ż_manual ≈ ż_multi
    end

    @testset "j′vp(::$T)" for T in (Float64, ComplexF64)
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, T, N), randn(rng, T, M)
        z̄ = randn(rng, T, N + M)
        xy = vcat(x, y)
        x̄ȳ_manual = j′vp(fdm, xy->sin.(xy), z̄, xy)[1]
        x̄ȳ_auto = j′vp(fdm, x->sin.(vcat(x[1], x[2])), z̄, (x, y))[1]
        x̄ȳ_multi = j′vp(fdm, (x, y)->sin.(vcat(x, y)), z̄, x, y)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_auto...)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_multi...)
    end
end


text:
using PyPlot
using RollingFunctions

function plot_acceptance(updts::Vector{<:ParamUpdate{MetropolisHastingsUpdt}})
    fig, ax = plt.subplots(2, length(updts), figsize=(3*length(updts), 5))
    for (i,updt) in enumerate(updts)
        ax[1,i].plot(updt.readjust_history, color="steelblue", linewidth=1.0)
        ax[2,i].plot(rollmean(updt.accpt_history, 100), color="steelblue", linewidth=1.0)
        coord = indices(updt.updt_coord)[1]
        print(coord)
        ax[1,i].set_title(string("Parameter ", string(coord)))
    end
    ax[1,1].set_ylabel("dispersion")
    ax[2,1].set_ylabel("acceptance rate")
    plt.tight_layout()
    ax
end

function plot_acceptance(updts::Vector{<:Imputation}, adaptations=nothing, adpt_intv=100)
    max_num_segm = maximum([length(updt.ρs) for updt in updts])
    fig, ax = plt.subplots(2*length(updts), max_num_segm, figsize=(3*max_num_segm, 5*length(updts)))
    if adaptations !== nothing
        adaptations = cumsum(adaptations)
    end
    for (i,updt) in enumerate(updts)
        for j in 1:length(updt.ρs_history)
            ρ = [ρs[j] for ρs in updt.ρs_history]
            ah = [a[j] for a in updt.accpt_history]
            ax[2*(i-1)+1,j].plot(ρ, color="steelblue", linewidth=1.0)
            ax[2*(i-1)+2,j].plot(rollmean(ah, adpt_intv), color="steelblue", linewidth=1.0)
            ylimA = ax[2*(i-1)+1,j].get_ylim()
            ylimB = ax[2*(i-1)+2,j].get_ylim()
            if adaptations !== nothing
                for adpt in adaptations
                    ax[2*(i-1)+1,j].plot([adpt/adpt_intv, adpt/adpt_intv], ylimA, color="red", linestyle="dashed", linewidth=1.0)
                    ax[2*(i-1)+2,j].plot([adpt, adpt], ylimB, color="red", linestyle="dashed", linewidth=1.0)
                end
            end
        end
        ax[2*(i-1)+1,1].set_ylabel("dispersion")
        ax[2*(i-1)+2,1].set_ylabel("acceptance rate")
    end
    plt.tight_layout()
    ax
end

function plot_chains(ws::MCMCWorkspace,indices=nothing;truth=nothing,figsize=(15,10),
                     ylims=nothing)
    θs = ws.θ_chain
    num_steps = length(θs)
    if indices === nothing
        indices = 1:length(θs[1])
    end
    n_chains = length(indices)

    fig, ax = plt.subplots(n_chains, 1, figsize=figsize, sharex=true)

    for (i,ind) in enumerate(indices)
        ax[i].plot([θ[ind] for θ in θs],color="steelblue", linewidth=0.5)
    end
    if truth !== nothing
        for (i,ind) in enumerate(indices)
            ax[i].plot([1,length(θs)],[truth[ind],truth[ind]],color="black",
                       linestyle="--")
        end
    end
    if ylims !== nothing
        for (i,axis) in enumerate(ax)
            if ylims[i] !== nothing
                axis.set_ylim(ylims[i])
            end
        end
    end
    plt.tight_layout()
    ax
end

function plot_paths(ws::Workspace, ws_mcmc::MCMCWorkspace, schedule, coords=nothing; transf=nothing,
                    figsize=(12,8), alpha=0.5, obs=nothing,
                    path_indices=1:length(ws.paths), ylims=nothing, θ_fixed=nothing)
    yy, tt = ws.paths, ws.time
    if coords === nothing
        coords = 1:length(yy[1][1])
    end
    if transf === nothing
        transf = [(x,θ)->x for _ in coords]
    end
    n_coords = length(coords)
    if θ_fixed !== nothing
        θs = fill(θ_fixed, length(yy))
    else
        θs = θs_for_transform(ws_mcmc, schedule)
    end

    fig, ax = plt.subplots(n_coords, 1, figsize=figsize, sharex=true)
    M = min(length(yy), length(path_indices))
    cmap = plt.get_cmap("BuPu")(range(0, 1, length=M)) #"cividis"
    for (i,ind) in enumerate(coords)
        cc = 1
        for j in 1:length(yy)
            if j in path_indices
                path = [transf[i](y, θs[j])[ind] for y in yy[j]]
                ax[i].plot(tt, path, color=cmap[cc,1:end], alpha=alpha, linewidth=0.5)
                cc += 1
            end
        end
    end

    _i = 1
    if obs !== nothing
        for (i,ind) in enumerate(coords)
            if ind in obs.indices
                ax[i].plot(obs.times, [o[_i] for o in obs.vals], mfc="orange",
                           mec="orange", marker="o", linestyle="none",
                           markersize=4)
                _i += 1
            end
        end
    end
    if ylims !== nothing
        for (i,axis) in enumerate(ax)
            if ylims[i] !== nothing
                axis.set_ylim(ylims[i])
            end
        end
    end
    plt.tight_layout()
    ax
end

function θs_for_transform(ws::Any) # old implementation
    θs = ws.θ_chain.θ_chain
    if ws.θ_chain.counter < 2
        θs = [params(ws.P[1].Target) for i in 1:length(ws.paths)]
    else
        warm_up, save_iter = ws.action_tracker.warm_up, ws.action_tracker.save_iter
        num_updts = ws.accpt_tracker.updt_len
        N = Int64((length(θs)-1)/num_updts) + warm_up
        θs = θs[[num_updts*(i-warm_up) for i in 1:N
                 if (i % save_iter == 0) && i > warm_up ]]
    end
    θs
end

function θs_for_transform(ws::MCMCWorkspace, schedule)
    θs = ws.θ_chain
    warm_up, save_iter = schedule.actions.warm_up, schedule.actions.save
    num_updts = length(schedule.updt_idx[1])-1
    N = Int64((length(θs)-1)/num_updts) + warm_up
    θs = θs[[num_updts*(i-warm_up) for i in 1:N
            if (i % save_iter == 0) && i > warm_up ]]
    θs
end

ext: jl


text:
# Container for a DataFrame row
immutable DataFrameRow{T <: AbstractDataFrame}
    df::T
    row::Int
end

function Base.getindex(r::DataFrameRow, idx::AbstractArray)
    return DataFrameRow(r.df[idx], r.row)
end

function Base.getindex(r::DataFrameRow, idx::Any)
    return r.df[r.row, idx]
end

function Base.setindex!(r::DataFrameRow, value::Any, idx::Any)
    return setindex!(r.df, value, r.row, idx)
end

Base.names(r::DataFrameRow) = names(r.df)
_names(r::DataFrameRow) = _names(r.df)

Base.sub(r::DataFrameRow, c) = DataFrameRow(r.df[[c]], r.row)

index(r::DataFrameRow) = index(r.df)

Base.length(r::DataFrameRow) = size(r.df, 2)

Base.endof(r::DataFrameRow) = size(r.df, 2)

Base.collect(r::DataFrameRow) = @compat Tuple{Symbol, Any}[x for x in r]

Base.start(r::DataFrameRow) = 1

Base.next(r::DataFrameRow, s) = ((_names(r)[s], r[s]), s + 1)

Base.done(r::DataFrameRow, s) = s > length(r)

Base.convert(::Type{Array}, r::DataFrameRow) = convert(Array, r.df[r.row,:])

# hash of DataFrame rows based on its values
# so that duplicate rows would have the same hash
function Base.hash(r::DataFrameRow, h::UInt)
    for col in columns(r.df)
        if isna(col, r.row)
            h = hash(false, h)
        else
            h = hash(true, hash(col[r.row], h))
        end
    end
    return h
end

# compare two elements in the array
_isequalelms(a::Array, i::Int, j::Int) = isequal(a[i], a[j])

# compare the two elements in the data array
function _isequalelms(a::DataArray, i::Int, j::Int)
    if isna(a, i)
        return isna(a, j)
    else
        return !isna(a, j) && isequal(a.data[i], a.data[j])
    end
end

# compare two elements in the pooled array
# NOTE assume there are no duplicated elements in the pool
_isequalelms(a::PooledDataArray, i::Int, j::Int) = isequal(a.refs[i], a.refs[j])

# comparison of DataFrame rows
# only the rows of the same DataFrame could be compared
# rows are equal if they have the same values (while the row indices could differ)
function Base.isequal(r1::DataFrameRow, r2::DataFrameRow)
    if r1.df !== r2.df
        throw(ArgumentError("Comparing rows from different frames not supported"))
    end
    if r1.row == r2.row
        return true
    end
    for col in columns(r1.df)
        if !_isequalelms(col, r1.row, r2.row)
            return false
        end
    end
    return true
end


text:
__precompile__()

module PiecewiseLinearOpt

import JuMP

export PWLFunction, UnivariatePWLFunction, BivariatePWLFunction, piecewiselinear

include("types.jl")
include("jump.jl")

end # module


text:
module PiecewiseSplicedArrays
using PiecewiseIncreasingRanges,Compat
import Base: tail, @propagate_inbounds

export PiecewiseSplicedArray

struct PiecewiseSplicedArray{TA,N,AA<:AbstractArray,R<:AbstractUnitRange} <: AbstractArray{TA,N}
    parent::AA
    axesranges::NTuple{N,PiecewiseUnitRange{Int,R}}
    function PiecewiseSplicedArray{T,N,AA,R}(a,axesranges) where {T,N,AA<:AbstractArray{T,N},R}
    	checknumberofinds(a,axesranges)
    	new{T,N,AA,R}(a,axesranges)
    end
end

function PiecewiseSplicedArray(a::AbstractVector{TA},
	axesranges::PiecewiseUnitRange{Int,R}) where {TA,R<:AbstractUnitRange}

	PiecewiseSplicedArray{TA,1,typeof(a),R}(a,(axesranges,))
end

function PiecewiseSplicedArray(a::AbstractVector{TA},
	axesranges::Vararg{<:AbstractUnitRange}) where {TA}
	ax = PiecewiseUnitRange([axesranges...])
	PiecewiseSplicedArray{TA,1,typeof(a),R}(a,(ax,))
end

function PiecewiseSplicedArray(a::AbstractArray{TA,N},
	axesranges::Vararg{PiecewiseUnitRange{Int,R},N}) where {TA,N,R<:AbstractUnitRange}
	PiecewiseSplicedArray{TA,N,typeof(a),R}(a,axesranges)
end

function PiecewiseSplicedArray(a::AbstractArray{TA,N},
	axesranges::NTuple{N,PiecewiseUnitRange{Int,R}}) where {TA,N,R<:AbstractUnitRange}
	axesranges = Tuple(PiecewiseUnitRange.(axesranges))
	PiecewiseSplicedArray{TA,N,typeof(a),R}(a,axesranges)
end

function PiecewiseSplicedArray(a::AbstractArray{TA,N},
	axesranges::Tuple{Vararg{Vector{R},N}}) where {TA,N,R<:AbstractUnitRange}
	ax = Tuple([PiecewiseUnitRange(axesranges[i]) for i in eachindex(axesranges)])
	PiecewiseSplicedArray{TA,N,typeof(a),R}(a,ax)
end

struct ArraySizeMismatch{T1,T2,T3} <: Exception
	dim::T1
	s :: T2
	ninds :: T3
end

struct DimensionMismatch{T1,T2} <: Exception
	ndims :: T1
	naxes :: T2
end

Base.showerror(io::IO, e::ArraySizeMismatch) = println("size along dimension $(e.dim) is $(e.s) but you have provided $(e.ninds) indices")
Base.showerror(io::IO, e::DimensionMismatch) = println("array has $(e.ndims) dimension but you have provided indices along $(e.naxes) axes")

function checknumberofdims(a,axesranges)
	length(axesranges) == ndims(a) || throw(DimensionMismatch(ndims(a),length(axesranges)))
end

function checknumberofinds(a,axesranges)
	s = size(a)
	for (dim,inds) in enumerate(axesranges)
		(s[dim] == length(inds)) || throw(ArraySizeMismatch(dim,s[dim],length(inds)))
	end
end

Base.parent(A::PiecewiseSplicedArray) = A.parent

Base.size(a::PiecewiseSplicedArray) = size(parent(a))
Base.size(a::PiecewiseSplicedArray,d) = size(parent(a),d)

Base.axes(A::PiecewiseSplicedArray) = A.axesranges
Base.axes(A::PiecewiseSplicedArray{T,N},d) where {T,N} = 1 <= d <= N ? A.axesranges[d] : 1:1

Base.IndexStyle(::Type{OA}) where {OA<:PiecewiseSplicedArray} = IndexStyle(parenttype(OA))
parenttype(::Type{PiecewiseSplicedArray{T,N,AA,R}}) where {T,N,AA,R} = AA
parenttype(A::PiecewiseSplicedArray) = parenttype(typeof(A))

@inline _parentinds(::Tuple{},::Tuple{}) = ()
@inline function _parentinds(ax::NTuple{N,<:PiecewiseUnitRange},inds::NTuple{N,Int}) where {N}
	(searchsortedfirst(ax[1],inds[1]),_parentinds(tail(ax),tail(inds))...)
end
@inline function parentinds(a::PiecewiseSplicedArray{<:Number,N},inds::NTuple{N,Int}) where {N}
	_parentinds(a.axesranges,inds)
end

@inline @propagate_inbounds function Base.getindex(A::PiecewiseSplicedArray{TA,N},
	I::Vararg{Int,N}) where {TA,N}

	@boundscheck checkbounds(A, I...)
	@inbounds ret = parent(A)[parentinds(A,I)...]
	ret
end

@inline @propagate_inbounds function Base.getindex(A::PiecewiseSplicedArray, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds ret = parent(A)[i]
    ret
end

@inline @propagate_inbounds function Base.setindex!(A::PiecewiseSplicedArray{TA,N},val,I::Vararg{Int,N}) where {TA,N}
	@boundscheck checkbounds(A, I...)
	@inbounds parent(A)[parentinds(A,I)...] = val
	val
end

@inline @propagate_inbounds function Base.setindex!(A::PiecewiseSplicedArray{TA,N},val,i::Int) where {TA,N}
	@boundscheck checkbounds(A, i)
	@inbounds parent(A)[i] = val
	val
end

Base.print_array(io::IO,A::PiecewiseSplicedArray) = Base.print_array(io,parent(A))

end


text:
const VT = Vector
const VVT = Vector{VT}

millionnth(coll::VT) = all_perms(coll::VT, 1_000_000)[end]

function all_perms(coll::VT)::VVT
  length(coll) == 0 && return [[]]
  length(coll) == 1 && return [coll]
  length(coll) == 2 && return [coll, coll |> reverse]

  # length of coll ≥ 3
  perms = VVT([])

  for d ∈ coll
    coll_d = filter(x -> x != d, coll)
    _perms = prepend_all!(d, all_perms(coll_d))  # rec. call
    push!(perms, _perms...)
  end

  return perms
end

function all_perms(coll::VT, limit)::VVT
  length(coll) == 0 && return [[]]
  length(coll) == 1 && return [coll]
  length(coll) == 2 && return [coll, coll |> reverse]

  limit == 0 && return [[]]

  # length of coll ≥ 3
  perms = VVT([])

  for d ∈ coll
    coll_d = filter(x -> x != d, coll)
    _perms = prepend_all!(d, all_perms(coll_d, limit - 1))  # rec. call
    push!(perms, _perms...)
  end

  return perms
end

function prepend_all!(item, perms::VVT)::VVT
  for perm ∈ perms
    pushfirst!(perm, item)
  end

  perms
end

prepend_all!(item::T, perms::Vector{Vector{T}}) where T <: Integer = prepend_all!(item, VVT(perms))

prepend_all!(item::T, perms::Vector{Vector{T}}) where T <: Any = prepend_all!(T(item), VVT(perms))

prepend_all!(item::T, perms::Vector{Vector{Any}}) where T <: Any = prepend_all!(item, VVT(perms))

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
immutable CanonicalCorrelation{T} <: RegressionModel
    X::Matrix{T}        # (centered) input X
    Y::Matrix{T}        # (centered) input Y
    coefX::Matrix{T}    # canonical coefficients of X
    coefY::Matrix{T}    # canonical coefficients of Y
    r::Vector{T}        # canonical correlations
end

coef(cc::CanonicalCorrelation) = (cc.coefX, cc.coefY)
Base.cor(cc::CanonicalCorrelation) = cc.r
scores(cc::CanonicalCorrelation) = (cc.X*cc.coefX, cc.Y*cc.coefY)

function canoncor(X::AbstractMatrix, Y::AbstractMatrix; center::Bool=true)
    if center
        X = X .- mean(X, 1)
        Y = Y .- mean(Y, 1)
    end
    qr1 = qrfact(X)
    qr2 = qrfact(Y)
    usv = svdfact(full(qr1[:Q])'*full(qr2[:Q]))
    A = qr1[:R]\usv[:U]
    B = qr2[:R]\usv[:Vt]'

    s = sqrt(size(X, 1) - 1)
    scale!(A, s)
    scale!(B, s)

    CanonicalCorrelation(X, Y, A, B, usv[:S])
end


text:
module tension_compression_examples

using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed2DStrain
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.MatDeforNeohookeanADModule: MatDeforNeohookeanAD
using FinEtoolsDeforNonlinear.MatDeforStVKModule: MatDeforStVK
using FinEtoolsDeforNonlinear.MatDeforStVKADModule: MatDeforStVKAD
using FinEtoolsDeforNonlinear.MatDeforI1RivlinADModule: MatDeforI1RivlinAD
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots

function neohookeanad_q4()
    mr = DeforModelRed2DStrain
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookeanAD(mr, E, nu)
    L= 6/2*phun("mm");
    W = 2/2*phun("mm");
    H = 2/2*phun("mm");
    umag = 1.5*phun("mm");# Magnitude of the displacement
    nincr = 48
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 15
    tolerance = W / 1000

    fens, fes = Q4block(L, W, 2, 1)

    l1  = selectnode(fens; box = [0,0,-Inf,Inf], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    l2  = selectnode(fens; box = [-Inf,Inf,0,0], inflate  =  tolerance)
    e2 = FDataDict("node_list"=>l2, "component"=>2)

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf], inflate  =  tolerance)
    table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
    move(x, lambda) = table(lambda);
    e4 = FDataDict("node_list"=>movel1, "component"=>1, "displacement"=>move)

    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(2, 2), H), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1],  "essential_bcs"=>[e1, e2, e4])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
        # print("-\n")
        # @show modeldata["un1"]
        # vtkexportmesh("junk$(iter).vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
        # print("-\n")
    end
    Ux = FFlt[]; Rx = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
        push!(Ux, mean(modeldata["un1"].values[movel1,1]));
        push!(Rx, sum(modeldata["reactions"].values[movel1,1]));
    end

    modeldata = nonlinearstatics(modeldata);
    # @show Ux, Rx
    @show (minimum(Ux), minimum(Rx), maximum(Ux), maximum(Rx))

    pl = lineplot((L .+ Ux) ./ L, Rx, canvas = DotCanvas)
    display(pl)
end # function neohookeanad_q4

function allrun()
    println("#####################################################")
    println("# neohookeanad_q4 ")
    neohookeanad_q4()
    return true
end # function allrun

end # module tension_compression_examples


text:
f = open("hightemp.txt")
n = int(ARGS[1])
filelength = countlines(f)
f = seekstart(f)
sepa = int(trunc(filelength / n))
mod = filelength % n

for i in 1:n
  out = open(string(i)*".txt", "w")
  lines = sepa
  if i <= mod
    lines = sepa + 1
  end

  for j in 1:lines
    print(out,readline(f))
  end
  close(out)
end


text:
module MaxHelpingHandRainbowSpinnerColorControllerDisabler

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/RainbowSpinnerColorControllerDisabler" RainbowSpinnerColorControllerDisabler(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Rainbow Spinner Color Controller (Disable)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        RainbowSpinnerColorControllerDisabler
    )
)

function Ahorn.selection(entity::RainbowSpinnerColorControllerDisabler)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RainbowSpinnerColorControllerDisabler, room::Maple.Room) = Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/rainbowSpinnerColorControllerDisable", 0, 0)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function calc_perimeter_ngon(a,b)
if length(a) == length(b) && length(a) >= 2
    L = 0;
    L = L + abs(a[1]+im*b[1]);
    for i = 2:length(a)
        L = L + abs((a[i]-a[i-1])+im*(b[i]-b[i-1]));
    end
    L = L + abs(a[length(a)]+im*b[length(a)]);
    return L
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using FredData
using Test

include("test_without_key.jl")

# Normal usage - API key must be present in ENV
if FredData.has_fred_key()
    include("test_with_key.jl")
end


text:
@testset "rotational damper jacobian" begin
    joint_types = [
        :Fixed,
        :Prismatic,
        :Planar,
        :FixedOrientation,
        :Revolute,
        :Cylindrical,
        :PlanarAxis,
        :FreeRevolute,
        :Orbital,
        :PrismaticOrbital,
        :PlanarOrbital,
        :FreeOrbital,
        :Spherical,
        :CylindricalFree,
        :PlanarFree
        ]
    for joint_type in joint_types
        mech = DojoEnvironments.get_snake(gravity=0.0, num_bodies=2, damper=0.3, joint_type=joint_type)
        DojoEnvironments.initialize_snake!(mech)
        function ctrl!(m,k)
            Dojo.set_input!(m, 0.01*ones(Dojo.minimal_dimension(m)))
        end
        storage = Dojo.simulate!(mech, 1.0, ctrl!)
        # Dojo.visualize(mech, storage, vis=vis)

        rot0 = mech.joints[2].rotational
        timestep0 = mech.timestep
        pbody0 = mech.bodies[1]
        cbody0 = mech.bodies[2]
        xa0, va0, qa0, ωa0 = Dojo.current_configuration_velocity(pbody0.state)
        xb0, vb0, qb0, ωb0 = Dojo.current_configuration_velocity(cbody0.state)

        # Configuration
        J0 = Dojo.damper_jacobian_configuration(:parent, :parent, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            xq -> timestep0 * Dojo.damper_force(:parent, rot0,
                Dojo.Quaternion(xq[4:7]...,true), ωa0, qb0, ωb0, timestep0; rotate=true, unitary=false),
            [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_configuration(:parent, :child, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            xq -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, ωa0,
                Dojo.Quaternion(xq[4:7]...,true), ωb0, timestep0; rotate=true, unitary=false),
            [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_configuration(:child, :parent, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            xq -> timestep0 * Dojo.damper_force(:child, rot0,
                Dojo.Quaternion(xq[4:7]...,true), ωa0, qb0, ωb0, timestep0; rotate=true, unitary=false),
            [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_configuration(:child, :child, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            xq -> timestep0 * Dojo.damper_force(:child, rot0, qa0, ωa0,
                Dojo.Quaternion(xq[4:7]...,true), ωb0, timestep0; rotate=true, unitary=false),
            [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        # Velocity
        J0 = Dojo.damper_jacobian_velocity(:parent, :parent, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            vω -> timestep0 * Dojo.damper_force(:parent, rot0, qa0,
                vω[Dojo.SUnitRange(4,6)], qb0, ωb0, timestep0; rotate=true, unitary=false),
            [va0; ωa0])
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_velocity(:parent, :child, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            vω -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, ωa0, qb0,
                vω[Dojo.SUnitRange(4,6)], timestep0; rotate=true, unitary=false),
            [vb0; ωb0])
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_velocity(:child, :parent, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            vω -> timestep0 * Dojo.damper_force(:child, rot0, qa0,
                vω[Dojo.SUnitRange(4,6)], qb0, ωb0, timestep0; rotate=true, unitary=false),
            [va0; ωa0])
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8

        J0 = Dojo.damper_jacobian_velocity(:child, :child, rot0, pbody0, cbody0, timestep0)
        J1 = ForwardDiff.jacobian(
            vω -> timestep0 * Dojo.damper_force(:child, rot0, qa0, ωa0, qb0,
                vω[Dojo.SUnitRange(4,6)], timestep0; rotate=true, unitary=false),
            [vb0; ωb0])
        norm(J0 - J1, Inf)
        @test norm(J0 - J1, Inf) < 1.0e-8
    end
end


text:
using Dates, Distributed

"""
    runbatch(sim, writer, expqueue)

Run the simulation in parallel, using the writer to collect data.

Calls the simulation function `sim` in parallel as `f(expqueue[idx], idx, writer)` for each entry
`expqueue[idx]` in `expqueue`. Individual runs of the simulation are automatically distributed over
all available workers by the master process. The master process uses the `writer` to collect and
store the results.

# Arguments
* `sim::Function`: the simulation function.
* `writer::Writer`: the Writer to use to store the data.
* `expqueue::Vector`: a Vector that the arguments for each run of the simulation.
"""
function runbatch(sim::Function, writer::Writer, expqueue::Vector)
    sink = open_sink(writer)
    try
        next_expidx = 0
        getnext_expidx() = (next_expidx += 1)

        overview = zeros(Int, nworkers())

        function updateoverview(pid, expidx)
            overview[ifelse(nprocs() == 1, 1, pid - 1)] = expidx
        end

        @sync begin
            @async begin
                while true
                    try
                        flushedbuffer = take!(writer.remotequeue)
                        insert!(sink, writer, flushedbuffer)
                    catch err
                        if isa(err, InvalidStateException) && err.state == :closed
                            break # Queue is empty and was closed because workers are done.
                        else
                            rethrow(err)
                        end
                    end
                end
            end
            @sync begin # Prevents race conditions when a worker's async task throws an error.
                println(now(), ": Queued ", length(expqueue), " experiments on ", nworkers(), " cores")
                for pid in workers()
                    @async begin
                        last_expidx = length(expqueue)
                        cur_expidx = getnext_expidx()
                        while cur_expidx <= last_expidx
                            updateoverview(pid, cur_expidx)
                            print("$(now()): Computing experiments $overview...\n")
                            remotecall_fetch(sim, pid, expqueue[cur_expidx], cur_expidx, writer)
                            cur_expidx = getnext_expidx()
                        end
                        updateoverview(pid, 0)
                        print("$(now()): Computing experiments $overview...\n")
                    end
                end
            end
            close(writer.remotequeue)
        end
        println(now(), ": All experiments have been completed")
    finally
        close_sink(sink, writer)
    end
    return get_return(writer)
end


text:
using JuMP
using SCS
using LinearAlgebra

q = 4 # dimension of estimate space
p = 8 # number of experimental vectors
nmax = 3 # upper bound on lambda
n = 12

V = randn(q, p)

eye = Matrix{Float64}(I, q, q);

aOpt = Model(with_optimizer(SCS.Optimizer, verbose = 0))
@variable(aOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(aOpt, u[1:q], lower_bound = 0)

@constraint(aOpt, sum(np) <= n)
for i = 1:q
    @SDconstraint(aOpt, [V * diagm(0 => np ./ n) * V' eye[:, i]; eye[i, :]' u[i]] >= 0)
end

@objective(aOpt, Min, sum(u))

optimize!(aOpt)

@show objective_value(aOpt);
@show value.(np);

eOpt = Model(with_optimizer(SCS.Optimizer, verbose = 0))
@variable(eOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(eOpt, t)

@SDconstraint(eOpt, V * diagm(0 => np ./ n) * V' - (t .* eye) >= 0)
@constraint(eOpt, sum(np) <= n)

@objective(eOpt, Max, t)

optimize!(eOpt)

@show objective_value(eOpt);
@show value.(np);

dOpt = Model(with_optimizer(SCS.Optimizer, verbose = 0))
@variable(dOpt, np[1:p], lower_bound = 0, upper_bound = nmax)
@variable(dOpt, t)
@objective(dOpt, Max, t)
@constraint(dOpt, sum(np) <= n)
E = V * diagm(0 => np ./ n) * V'
@constraint(dOpt, [t, 1, (E[i, j] for i in 1:q for j in 1:i)...] in MOI.LogDetConeTriangle(q))

optimize!(dOpt)

@show objective_value(dOpt);
@show value.(np);


text:
mutable struct parameters
      H::Float64
     Ly::Float64
     f₀::Float64
      g::Float64
     Lj::Float64
     Uj::Float64
     Ny::Int64
     dk::Float64
     kₘ::Float64
     Fr::Float64
     Ro::Float64
end

function parameters(;
    H  = 500,
    Ly = 2e5,
    f₀ = 1e-4,
    g  = 9.81,
    Lj = 1e4,
    Uj = 0.1,
    Ny = 100,
    dk = 5e-6,
    kₘ = 2e-4,
    Fr = 0.0,
    Ro = 0.0
    )
    return parameters(H, Ly, f₀ , g, Lj, Uj, Ny, dk, kₘ, Fr, Ro)
end


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

mutable struct VirtualNetworkGatewayConnectionType <: SwaggerModel

    function VirtualNetworkGatewayConnectionType(;)
        o = new()
        o
    end
end # type VirtualNetworkGatewayConnectionType

const _property_map_VirtualNetworkGatewayConnectionType = Dict{Symbol,Symbol}()
const _property_types_VirtualNetworkGatewayConnectionType = Dict{Symbol,String}()
Base.propertynames(::Type{ VirtualNetworkGatewayConnectionType }) = collect(keys(_property_map_VirtualNetworkGatewayConnectionType))
Swagger.property_type(::Type{ VirtualNetworkGatewayConnectionType }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkGatewayConnectionType[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayConnectionType }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayConnectionType[property_name]

function check_required(o::VirtualNetworkGatewayConnectionType)
    true
end

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


text:
# inspired by Jiahao Chen's post:
# https://groups.google.com/d/msg/julia-dev/iBBSQt1FGHE/KaXNq24Q7l8J

using Traits

@traitdef Group{X,Y} begin  # a trait defining a group
    @constraints begin
        X==Y
    end
    *(X,Y)   -> X
    *(X,Y,X) -> X
    one(X)   -> X
    inv(X)   -> X
end

@show istrait(Group{Int,Int}) # false
@show istrait(Group{Float64, Float64}) #true

#false, should be true:
@show istrait(Group{Array{Float64,2}, Array{Float64,2}}, verbose=true)
# Because inv can error when matrix is not square.
# However, it does not complain about `one`?

@traitfn f{X; Group{X,X}}(x::X, y::X) = inv(x*y*y)

try
    f(5,6)
catch e
    println(e) # ERROR: TraitException("No matching trait found for function f")
end

@show f(5.,6.)

###
# Version 2: use only one parameter, maybe better?
###
@traitdef Group2{X} begin  # a trait defining a group
    *(X,X)   -> X
    *(X,X,X) -> X
    one(X)   -> X
    inv(X)   -> X
end

@show istrait(Group2{Int}) # false
@show istrait(Group2{Float64}) #true

#false, should be true:
@show istrait(Group2{Array{Float64,2}}, verbose=true)
# Because inv can error when matrix is not square.
# However, it does not complain about `one`?

@traitfn f2{X; Group2{X}}(x::X, y::X) = inv(x*y*y)

try
    f2(5,6)
catch e
    println(e) # ERROR: TraitException("No matching trait found for function f")
end

@show f2(5.,6.)


text:
# ir_dump.jl

using MIRT: ir_dump

using Test: @test

x = (int=1, real=2., string="String", shorttuple=(1,2), longtuple=(1,2,3,4), )
@test ir_dump(x, io=isinteractive() ? stdout : IOBuffer()) isa Nothing


text:
@testset "Non-commutative Monomial" begin
    Mod.@ncpolyvar x
    X = constantmonomial(typeof(x))
    @test nvariables(X) == 0
    @test isempty(X.vars)
    @test isempty(X.z)
    X = constantmonomial(x)
    @test nvariables(X) == 1
    @test variables(X)[1] == x
    @test collect(exponents(X)) == [0]
    X = monomial(x)
    @test nvariables(X) == 1
    @test variables(X)[1] == x
    @test collect(exponents(X)) == [1]
    X = x^0
    Y = X * x^2
    @test nvariables(Y) == 1
    @test variables(Y)[1] == x
    @test collect(exponents(Y)) == [2]
    Y = x^2 * X
    @test nvariables(Y) == 1
    @test variables(Y)[1] == x
    @test collect(exponents(Y)) == [2]
end
@testset "Non-commutative MonomialVector" begin
    Mod.@ncpolyvar x y
    X = emptymonovec(typeof(x))
    @test iszero(nvariables(X))
    @test isempty(variables(X))
    X = emptymonovec(x)
    @test nvariables(X) == 1
    @test variables(X)[1] == x
    X = monovec([y, x])
    @test X[1] == x
    @test X[2] == y
    X = monomials([x, y], 2)
    @test variables(X)[1] == x
    @test variables(X)[2] == y
    @test variables(X)[3] == x
    @test collect(exponents(X[1])) == [2, 0, 0]
    @test collect(exponents(X[2])) == [1, 1, 0]
    @test collect(exponents(X[3])) == [0, 2, 0]
    @test collect(exponents(X[4])) == [0, 1, 1]
    X0 = [x^3, x^2*y, x*y^2, x*y*x, y^3, y^2*x, y*x^2, y*x*y]
    X1 = monomials([x, y], 3)
    X2 = monovec(X0)
    Z = [[3,0,0,0],[2,1,0,0],[1,2,0,0],[1,1,1,0],[0,3,0,0],[0,2,1,0],[0,1,2,0],[0,1,1,1]]
    @test X1 isa AbstractVector{<:AbstractMonomial}
    @test X2 isa AbstractVector{<:AbstractMonomial}
    @test length(X1) == length(X2) == length(Z)
    @test issorted(X1, rev=true)
    @test issorted(X2, rev=true)
    @test X1 == X0
    for i in 1:length(Z)
        @test variables(X1[i])[1] == variables(X2[i])[1] == x
        @test variables(X1[i])[2] == variables(X2[i])[2] == y
        @test variables(X1[i])[3] == variables(X2[i])[3] == x
        @test variables(X1[i])[4] == variables(X2[i])[4] == y
        @test collect(exponents(X1[i])) == collect(exponents(X2[i])) == Z[i]
    end
end


text:
using Mocking

copy_backup_1_patch = @patch function copy_backup(a...; b...)
    return Dict{String, Any}("Backup" => Dict{String, Any}("ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:backup/backup-0a3364eded1014b28", "CreationTime" => 1.617954808068e9, "Tags" => Any[Dict{String, Any}("Value" => "MyBackup", "Key" => "Name")], "BackupId" => "backup-0a3364eded1014b28", "Type" => "USER_INITIATED", "KmsKeyId" => "arn:aws:fsx:us-east-1:012345678912:key/d1234e22-543a-12b7-a98f-e12c2b54001a", "Lifecycle" => "COPYING", "OwnerId" => "123456789012", "FileSystem" => Dict{String, Any}("StorageCapacity" => 2400, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:file-system/fs-0f5179e395f597e66", "LustreConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 0, "PerUnitStorageThroughput" => 50, "DeploymentType" => "PERSISTENT_1", "WeeklyMaintenanceStartTime" => "1:05:00"), "StorageType" => "SSD", "FileSystemId" => "fs-0498eed5fe91001ec", "FileSystemType" => "LUSTRE")))
end

create_backup_1_patch = @patch function create_backup(a...; b...)
    return Dict{String, Any}("Backup" => Dict{String, Any}("ProgressPercent" => 0, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:backup/backup-03e3c82e0183b7b6b", "CreationTime" => "1481841524.0", "Tags" => Any[Dict{String, Any}("Value" => "MyBackup", "Key" => "Name")], "BackupId" => "backup-03e3c82e0183b7b6b", "Type" => "USER_INITIATED", "Lifecycle" => "CREATING", "FileSystem" => Dict{String, Any}("StorageCapacity" => 300, "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 30, "DailyAutomaticBackupStartTime" => "05:00", "WeeklyMaintenanceStartTime" => "1:05:00", "ActiveDirectoryId" => "d-1234abcd12"), "FileSystemId" => "fs-0498eed5fe91001ec", "OwnerId" => "012345678912")))
end

create_file_system_1_patch = @patch function create_file_system(a...; b...)
    return Dict{String, Any}("FileSystem" => Dict{String, Any}("Lifecycle" => "CREATING", "OwnerId" => "012345678912", "StorageCapacity" => 3200, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:file-system/fs-0123456789abcdef0", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "VpcId" => "vpc-ab1234cd", "CreationTime" => "1481841524.0", "SubnetIds" => Any["subnet-1234abcd"], "FileSystemId" => "fs-0123456789abcdef0", "DNSName" => "fs-0123456789abcdef0.fsx.com", "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 30, "ThroughputCapacity" => 32, "DailyAutomaticBackupStartTime" => "05:00", "Aliases" => Any[Dict{String, Any}("Lifecycle" => "CREATING", "Name" => "accounting.corp.example.com")], "WeeklyMaintenanceStartTime" => "1:05:00", "ActiveDirectoryId" => "d-1234abcd12"), "StorageType" => "HDD", "KmsKeyId" => "arn:aws:kms:us-east-1:012345678912:key/1111abcd-2222-3333-4444-55556666eeff"))
end

create_file_system_from_backup_1_patch = @patch function create_file_system_from_backup(a...; b...)
    return Dict{String, Any}("FileSystem" => Dict{String, Any}("Lifecycle" => "CREATING", "OwnerId" => "012345678912", "StorageCapacity" => 300, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:file-system/fs-0498eed5fe91001ec", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "VpcId" => "vpc-ab1234cd", "CreationTime" => "1481841524.0", "SubnetIds" => Any["subnet-1234abcd"], "FileSystemId" => "fs-0498eed5fe91001ec", "DNSName" => "fs-0498eed5fe91001ec.fsx.com", "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 30, "ThroughputCapacity" => 8, "DailyAutomaticBackupStartTime" => "05:00", "WeeklyMaintenanceStartTime" => "1:05:00", "ActiveDirectoryId" => "d-1234abcd12"), "KmsKeyId" => "arn:aws:kms:us-east-1:012345678912:key/0ff3ea8d-130e-4133-877f-93908b6fdbd6"))
end

delete_backup_1_patch = @patch function delete_backup(a...; b...)
    return Dict{String, Any}("BackupId" => "backup-03e3c82e0183b7b6b", "Lifecycle" => "DELETED")
end

delete_file_system_1_patch = @patch function delete_file_system(a...; b...)
    return Dict{String, Any}("FileSystemId" => "fs-0498eed5fe91001ec", "Lifecycle" => "DELETING")
end

describe_backups_1_patch = @patch function describe_backups(a...; b...)
    return Dict{String, Any}("Backups" => Any[Dict{String, Any}("ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:backup/backup-03e3c82e0183b7b6b", "CreationTime" => "1481841524.0", "Tags" => Any[Dict{String, Any}("Value" => "MyBackup", "Key" => "Name")], "BackupId" => "backup-03e3c82e0183b7b6b", "Type" => "USER_INITIATED", "Lifecycle" => "AVAILABLE", "FileSystem" => Dict{String, Any}("StorageCapacity" => 300, "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 30, "DailyAutomaticBackupStartTime" => "05:00", "WeeklyMaintenanceStartTime" => "1:05:00", "ActiveDirectoryId" => "d-1234abcd12"), "FileSystemId" => "fs-0498eed5fe91001ec", "OwnerId" => "012345678912"))])
end

describe_file_systems_1_patch = @patch function describe_file_systems(a...; b...)
    return Dict{String, Any}("FileSystems" => Any[Dict{String, Any}("Lifecycle" => "AVAILABLE", "OwnerId" => "012345678912", "StorageCapacity" => 300, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:file-system/fs-0498eed5fe91001ec", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "VpcId" => "vpc-ab1234cd", "CreationTime" => "1481841524.0", "SubnetIds" => Any["subnet-1234abcd"], "FileSystemId" => "fs-0498eed5fe91001ec", "NetworkInterfaceIds" => Any["eni-abcd1234"], "DNSName" => "fs-0498eed5fe91001ec.fsx.com", "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 30, "ThroughputCapacity" => 8, "DailyAutomaticBackupStartTime" => "05:00", "WeeklyMaintenanceStartTime" => "1:05:00", "ActiveDirectoryId" => "d-1234abcd12"), "KmsKeyId" => "arn:aws:kms:us-east-1:012345678912:key/0ff3ea8d-130e-4133-877f-93908b6fdbd6")])
end

list_tags_for_resource_1_patch = @patch function list_tags_for_resource(a...; b...)
    return Dict{String, Any}("Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")])
end

update_file_system_1_patch = @patch function update_file_system(a...; b...)
    return Dict{String, Any}("FileSystem" => Dict{String, Any}("Lifecycle" => "AVAILABLE", "OwnerId" => "012345678912", "StorageCapacity" => 300, "ResourceARN" => "arn:aws:fsx:us-east-1:012345678912:file-system/fs-0498eed5fe91001ec", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "VpcId" => "vpc-ab1234cd", "CreationTime" => "1481841524.0", "SubnetIds" => Any["subnet-1234abcd"], "FileSystemId" => "fs-0498eed5fe91001ec", "DNSName" => "fs-0498eed5fe91001ec.fsx.com", "WindowsConfiguration" => Dict{String, Any}("AutomaticBackupRetentionDays" => 10, "ThroughputCapacity" => 8, "DailyAutomaticBackupStartTime" => "06:00", "WeeklyMaintenanceStartTime" => "3:06:00"), "KmsKeyId" => "arn:aws:kms:us-east-1:012345678912:key/0ff3ea8d-130e-4133-877f-93908b6fdbd6"))
end


text:
"""
  makeblocksize(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}) where {T <: AbstractFloat, N}

Returns Tuple{Int, Int} of the size ob the block in the output matrix

"""

function makeblocksize(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}) where {T <: AbstractFloat, N}
  a = bm.bls
  b = bm.bls
  if !bm.sqr
    if i[1] == bm.bln
      a = bm.dats - (bm.bln-1)*bm.bls
    end
    if i[2] == bm.bln
      b = bm.dats - (bm.bln-1)*bm.bls
    end
  end
  a,b
end

"""
  computeblock(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}, dims::Tuple) where {T <: AbstractFloat, N}

Returns Matrix{T}, the single block of the output higher correlation matrix

"""

function computeblock(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}, dims::Tuple) where {T <: AbstractFloat, N}
  x = bm.bln^(N-2)
  R = zeros(T, makeblocksize(bm, i))
  for j in 1:x
    @inbounds k = Tuple(CartesianIndices(dims)[j])
    for k1 in 1:bm.bln
      @inbounds M1 = unfold(getblock(bm, (i[1],k1, k...)),1)
      @inbounds M2 = unfold(getblock(bm, (k1,i[2],k...)),2)
      @inbounds R += M1*transpose(M2)
    end
  end
  return R
end

"""
  computeblock_p(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}, dims::Tuple) where {T <: AbstractFloat, N}

Returns Matrix{T}, the single block of the output higher correlation matrix
Parallel implementation
"""

function computeblock_p(bm::SymmetricTensor{T, N}, i::Tuple{Int, Int}, dims::Tuple) where {T <: AbstractFloat, N}
  x = bm.bln^(N-2)
  @inbounds R = SharedArray(zeros(T, (makeblocksize(bm, i)..., x)))
  @sync @distributed for j in 1:x
    @inbounds k = Tuple(CartesianIndices(dims)[j])
    for k1 in 1:bm.bln
      @inbounds M1 = unfold(getblock(bm, (i[1],k1, k...)),1)
      @inbounds M2 = unfold(getblock(bm, (k1,i[2],k...)),2)
      @inbounds R[:,:,j] += M1*transpose(M2)
    end
  end
  @inbounds R = Array(R)
  return mapreduce(i -> R[:,:,i], +, 1:size(R,3))
end
"""
  cum2mat(bm::SymmetricTensor{T, N}) where {T <: AbstractFloat, N}

Returns higher order correlation matrix in the form of bm::SymmetricTensor{T, 2}

```jldoctest

julia> Random.seed!(42);

julia> t = rand(SymmetricTensor{Float64, 3}, 4);

julia> cum2mat(t)
SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[7.69432 4.9757; 4.9757 5.72935] [6.09424 4.92375; 5.05157 3.17723]; nothing [7.33094 4.93128; 4.93128 4.7921]], 2, 2, 4, true)

```
"""

function cum2mat(bm::SymmetricTensor{T, N}) where {T <: AbstractFloat, N}
    ret = arraynarrays(T, bm.bln, bm.bln)
    ds = (fill(bm.bln, N-2)...,)
    if nworkers() == 1
      for i in pyramidindices(2, bm.bln)
        @inbounds ret[i...] = computeblock(bm, i, ds)
      end
    else
      for i in pyramidindices(2, bm.bln)
        @inbounds ret[i...] = computeblock_p(bm, i, ds)
      end
    end
    SymmetricTensor(ret; testdatstruct = false)
end


text:
######################################################################
# filters.jl
# Handles the different filters
######################################################################

abstract type AbstractFilter end

"""
`update!(f::AbstractFilter, p::Pose, o::Real)`

`update!(f::AbstractFilter, x::Vehicle, o::Real)`

Updates the belief of filter `f` given a vehicle `x` and observation `o`.
The observation `o` is a real number in [0,360).
If your filter requires this number to be binned into an integer, it must take care of that internally.
"""
function update!(f::AbstractFilter, x::Vehicle, o)
    update!(f, get_pose(x), o)
end

function update!(f::AbstractFilter, p::Pose, o)
    error(typeof(f), " does not yet implement update!(f,x,o).")
end

# multi vehicle update
function update!(f::AbstractFilter,vx::Vector{Vehicle}, vo::Vector{Float64})
    num_vehicles = length(vx)
    for i = 1:num_vehicles
        update!(f, vx[i], vo[i])
    end
end

"""
`centroid(f::AbstractFilter)`

Returns the centroid of the filter's belief. If the belief is a Gaussian, this will return the mean.
"""
function centroid(f::AbstractFilter)
    error(typeof(f), " does not yet implement centroid(f).")
end

"""
`covariance(f::AbstractFilter)`

Returns the covariance of the filter's belief. If the belief is a Gaussian, this will return Sigma.
"""
function covariance(f::AbstractFilter)
    error(typeof(f), " does not yet implement centroid(f).")
end

"""
`entropy(f::AbstractFilter)`

Returns the entropy of the filter's belief.
"""
function entropy(f::AbstractFilter)
    error(typeof(f), " does not yet implement entropy(f).")
end

"""
`reset!(f::AbstractFilter)`

Resets the belief of filter `f`.
"""
function reset!(f::AbstractFilter)
    error(typeof(f), " does not yet implement reset!(f).")
end

# Include the different filters
include("df.jl")
include("kalman.jl")
include("ekf.jl")
include("eif.jl")
include("ukf.jl")
include("pf.jl")


text:
function get_cover(bam::Bam, index::BamIndex)

end


text:
# VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ImagePlot

using FunctionalDataUtils
using Compose, Colors, Compat
import Cairo, FixedPointNumbers
import Compose.compose
import Images

export plot

# M x N  or  M x N x 3 image of float in 0..1
function compose{T<:Real,N}(img::Array{T,N}, c::Context; usealpha = true)
    sm, sn = size(img)

    surface = Cairo.CairoARGBSurface(zeros(UInt32, sn, sm)) # width and height
    draw(PNG(surface), c)
    overlay = reinterpret(UInt8, surface.data, (4, sm, sn)) / 255

    outimg = similar(img, sm, sn, 3)
    imgind = (ndims(img) == 3 && size(img,3) == 3) ? [1,2,3] : [1,1,1]
    for m = 1:sm, n = 1:sn
        alpha = overlay[4,m,n]
        alpha = usealpha ? alpha : convert(typeof(alpha), alpha > 0.5)
        beta = one(alpha) - alpha
        for i = 1:3
            outimg[m,n,i] = beta*img[m,n,imgind[i]] + alpha*overlay[4-i,m,n]
        end
    end

    outimg
end

function plot{T<:Number}(img, points::Array{T,2}, labels::AbstractArray = [];
    compact = false,
    radius = 3px,
    fillcolor = colorant"yellow",
    strokecolor = nothing,
    colormap = [],
    kargs...)
    isempty(points) && return img
    c = context()
    function plotpoints(points, fillcolor)
        newc = context()
        for i = 1:size(points,2)
            newc = compose(newc,  circle(points[1,i],points[2,i], radius))
        end
        if strokecolor != nothing
            newc = compose(newc, stroke(color(strokecolor)))
        end
        newc = compose(newc, fill(fillcolor))
        c = compose(c, newc)
    end
    if isempty(labels)
        plotpoints(points, fillcolor)
    else
        U = sort(unique(labels))
        if all(x->isa(x,Real), U)
            if isempty(colormap)
                colors = @p jetcolormap len(U) | unstack | map x->RGB(x...)
            else
                colors = colormap
            end
        else
            colors = U
        end
        for i in 1:length(U)
            plotpoints(points[:,vec(labels .== U[i])], colors[i])
        end
    end
    width, height = size(img)
    c = compose(context(units=UnitBox(0, 0, width, height)), c)
    r = compose(img, c; kargs...)
    if compact
        mi = round(Int,minimum(points,2))
        ma = round(Int,maximum(points,2))
        r = r[clamp(mi[1]-20:ma[1]+20,1,size(r,1)), clamp(mi[2]-20:ma[2]+20,1,size(r,2)), :]
    end
    r
end

end # module


text:
using VLConstraintBasedModelGenerationUtilities

# build an array of compounds -
formula_array = [
    "C6H12O6"        ;      # 1 glucose
    "C6H13O9P"       ;      # 2 g6p
    "C10H16N5O13P3"  ;      # 3 atp
    "C10H15N5O10P2"  ;      # 4 adp
];

# test -
compound_dict = parse_molecular_formula_string(formula_array[4]) |> check

# compute the atom matrix -
atom_matrix = build_atom_matrix(formula_array) |> check


text:
module ParticleFilter

using RayCast.Bresenham: cast_heading
using Serialization
using RoboLib.Util: rangebearing2point, binarize, img2grid, grid2img, binarize, euclidean
using mushr_pf.SensorModels: LaserScanModel, DiscBeamPDF
using mushr_pf.SensorModels
using mushr_pf.MotionModels: AckerParams, AckerData, step!, BasicStochasticPredictModel, predict!, reset!
using mushr_pf.MotionModels
using RoboLib.Geom: Pose2D, Scale2D, T2D, R2D
using CoordinateTransformations, Rotations
using StaticArrays

include("particlefilter.jl")

struct OccMap{M<:AbstractMatrix{Bool}, TWM, TMW, P, O, S}
    map::M
    T_WM::TWM
    T_MW::TMW
    pos_WM::P
    ori_WM::O
    scale_WM::S
end
function OccMap(map::M, pos_WM::P, ori_WM::O, scale_WM::S) where {M<:AbstractMatrix{Bool},P<:T2D,O<:R2D,S<:Real}
    T_WM = Translation(pos_WM) ∘ LinearMap(ori_WM) ∘ LinearMap(Scale2D(scale_WM))
    T_MW = inv(T_WM)
    OccMap{M, typeof(T_WM), typeof(T_MW), P, O, S}(map, T_WM, T_MW, pos_WM, ori_WM, scale_WM)
end
OccMap(testfn, map, pos_WM, ori_WM, scale_WM) = OccMap(binarize(map, testfn), pos_WM, ori_WM, scale_WM)

# Dangerous indeed!
@inline (m::OccMap)(x::Integer, y::Integer) = @inbounds m.map[x,y]

function setup_predict_model(c::Dict{Symbol, Any})
    mc = c[:predictconf]

    if mc[:type] == :ackermann
        data = AckerData{c[:dtype]}()
        params = AckerParams{c[:dtype]}(mc[:car_length])
    else
        error("\"$mc[:type]\" not a known model")
    end

    if mc[:stochastic]
        process_dist = MvNormal([0,0,0], [mc[:sigma_x], mc[:sigma_y], mc[:sigma_theta]])
        ctrl_dist = MvNormal([0,0], [mc[:vel_sigma], mc[:delta_sigma]])
        init_dist = MvNormal([0,0,0], [mc[:init_sigma_x], mc[:init_sigma_y], mc[:init_sigma_theta]])
        model = BasicStochasticPredictModel(data, params, ctrl_dist, process_dist, init_dist)
    else
        error("not yet implemented")
    end
    model
end

function setup_reweight_model(om::OccMap, c::Dict{Symbol, Any})
    sc = c[:reweightconf]

    if sc[:type] == :beammodel
        #TODO make configurable
        raymethodW = RangeMethodW(om, sc[:max_range_meters])
        beammodel = DiscBeamPDF(om.scale_WM, sc[:max_range_meters]; a_short_lin=sc[:a_short_lin], a_short_exp=sc[:a_short_exp], a_max=sc[:a_max], a_rand=sc[:a_rand], a_hit=sc[:a_hit], sigma_hit=sc[:sigma_hit], lambda_short=sc[:lambda_short])
        reweightmodel = LaserScanModel{c[:dtype]}(beammodel, raymethodW, sc[:inv_squash_factor], sc[:nrays], sc[:method])
    else
        error("not yet implemented")
    end
    reweightmodel
end

struct RangeMethodW{D, M<:OccMap, MR<:Real}
    occmap::M
    mr_M::MR
    function RangeMethodW{D}(map::M, mr_W::MR) where {D,M,MR}
        new{D,M,MR}(map, mr_W/map.scale_WM)
    end
    RangeMethodW(map, mr_W) = RangeMethodW{Int}(map, mr_W)
end

function (rm::RangeMethodW{D})(T_WP::Pose2D{T}, angle) where {D,T}
    T_MP = rm.occmap.T_MW ∘ T_WP
    x_M, y_M, theta_M = T_MP.statev
    x_G, y_G, theta_G = img2grid(x_M, y_M, theta_M)
    xhit_G, yhit_G = cast_heading(D, x_G, y_G, theta_G+angle, rm.mr_M, rm.occmap)
    T(rm.occmap.scale_WM * euclidean(x_G, y_G, xhit_G, yhit_G))::T
end

function setup_resampler(c)
    rs = c[:resampler]
    if rs === :naive
        return resample_naive!
    elseif rs === :lowvar
        error("fix me") #TODO
        return resample_lowvar!
    else
        error("Resampler \"$rs\" unknown!")
    end
end

function setup_pf(occmap::OccMap, c::Dict{Symbol, Any})
    predictmodel = setup_predict_model(c)
    reweightmodel = setup_reweight_model(occmap, c)
    resampler = setup_resampler(c)
    particles = BufferedWeightedParticleBelief{Pose2D{c[:dtype]}}(c[:max_particles])

    StatefulParticleFilter(particles, predictmodel, reweightmodel, resampler, c[:rng])
end

end # module


text:
using Literate

freadme=joinpath(@__DIR__, "page1.jl")
#output_folder=joinpath(@__DIR__, "..","..","..","notebooks","modeling")
output_folder=joinpath(@__DIR__, ".")

#Literate.notebook(freadme, output_folder)
Literate.markdown(freadme, output_folder, documenter=true)

#freadme=joinpath(@__DIR__, "page2.jl")
#Literate.markdown(freadme, output_folder, documenter=true)
#Literate.notebook(freadme, output_folder)


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

using Markdown
using InteractiveUtils

# ╔═╡ 0c4edc09-133b-4cb8-a3e1-1e51cf9141f9
using Pkg, DrWatson

# ╔═╡ 8830efde-f966-4c2f-8370-6e6c45d86098
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers
	using DataFrames, Plots, Interpolations
end

# ╔═╡ cae6336b-054b-4c8e-ae88-93628ec5a6ae
md" ## j-05-01.jl"

# ╔═╡ d365b053-99cc-4775-9763-de6d7a8a9b19
begin
	x = [1.0, 3.0, 6.0, 5.0]
	y = [1.0, 5.0, 10.0, 9.0]
	xi = [2.0, 4.5]
	xint = 1:0.1:6
	indxs = sortperm(x)
end

# ╔═╡ deb9dac5-dfda-458d-9bea-a2659bf67555
md" #####

# ╔═╡ 03a53e2f-95f8-416e-888b-a50cca36fa13
begin
	interp_linear = LinearInterpolation(x[indxs], y[indxs])
	interp_linear.(xi)
end

# ╔═╡ 51fcc406-c130-49d3-9ba9-40437db90e05
begin
	x1 = 1:1:6
	y1 = [1.0, 3.0, 5.0, 7.0, 9.0, 10.0]
	interp = CubicSplineInterpolation(x1, y1)
	interp.(xi)
end

# ╔═╡ 5fded190-bcd2-461f-a7b8-72a4cce4751d
(dfin1, dfxi) = lagrangianpolynomial(length(x), x, y, xi)

# ╔═╡ 0b7ce582-0b49-4d27-b129-4e6e24d9eb15
(dfin2, dfxint) = lagrangianpolynomial(length(x), x, y, collect(xint))

# ╔═╡ 8546534a-55be-4447-9e1a-a722d60624fa
begin
	plot(dfxint.xi, dfxint.yi, line=(:path, 1), label="interpolated curve",
		leg=:topleft)
	scatter!(dfin1.x, dfin1.y, marker=(:circle, 4), label="input points",
		color=:blue)
	scatter!(dfxi.xi, dfxi.yi, marker=(:star, 8), color=:red,
		label="interpolated points")
end

# ╔═╡ 8069f1b3-8a61-436e-91a4-c0ecf0db8685
md" ## End of j-05-01.jl"

# ╔═╡ Cell order:
# ╟─cae6336b-054b-4c8e-ae88-93628ec5a6ae
# ╠═0c4edc09-133b-4cb8-a3e1-1e51cf9141f9
# ╠═8830efde-f966-4c2f-8370-6e6c45d86098
# ╠═d365b053-99cc-4775-9763-de6d7a8a9b19
# ╠═deb9dac5-dfda-458d-9bea-a2659bf67555
# ╠═03a53e2f-95f8-416e-888b-a50cca36fa13
# ╠═51fcc406-c130-49d3-9ba9-40437db90e05
# ╠═5fded190-bcd2-461f-a7b8-72a4cce4751d
# ╠═0b7ce582-0b49-4d27-b129-4e6e24d9eb15
# ╠═8546534a-55be-4447-9e1a-a722d60624fa
# ╟─8069f1b3-8a61-436e-91a4-c0ecf0db8685


text:
stress_scale = 1.0e+5
strain_scale = 1

# tid = parse(Int64, ARGS[1])
force_scale = 5.0
tid = 200
# if Sys.MACHINE=="x86_64-pc-linux-gnu"
#    global tid = parse(Int64, ARGS[1])
#    global force_scale = parse(Float64, ARGS[2])
# end
printstyled("tid=$tid\n", color=:green)

testtype = "NeuralNetwork2D"
nntype = "piecewise"
include("nnutil.jl")

# H0 = [1.04167e6  2.08333e5  0.0
#       2.08333e5  1.04167e6  0.0
#       0.0        0.0        4.16667e5]/stress_scale

H0 = [1.0406424793819175e6 209077.08366547766         0.0
      209077.08366547766   1.0411467691352057e6       0.0
      0.0                  0.0                   419057.32049008965]/stress_scale

#s = ae_to_code("Data/$(nntype)/NNLearn_$(idx)_$(H_function)_ite50.mat", nntype)
#s = ae_to_code("Data/$(nntype)/NNPreLSfit_$(idx)_$(H_function)_2.mat", nntype)
s = ae_to_code("Data/$(nntype)/NN_Train_$(idx)_from_3_$(H_function)_ite50.mat", nntype)

eval(Meta.parse(s))
# density 4.5*(1 - 0.25) + 3.2*0.25
#fiber_fraction = 0.25
#todo
prop = Dict("name"=> testtype, "rho"=> 4.5, "nn"=>post_nn)

T = 200.0
NT = 200

# nx_f, ny_f = 12, 4
# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5

porder = 2
Lx, Ly = 1.0, 0.5
nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny,porder, Lx, Ly; force_scale=force_scale)

ndofs=2
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

ts = LinRange(0, T, NT+1)
adaptive_solver_args = Dict("Newmark_rho"=> 0.0,
                          "Newton_maxiter"=>10,
                          "Newton_Abs_Err"=>1e-4,
                          "Newton_Rel_Err"=>1e-6,
                          "damped_Newton_eta" => 1.0)

globdat, domain, ts = AdaptiveSolver("NewmarkSolver", globdat, domain, T, NT, adaptive_solver_args)

# Δt = T/NT

# ρ_oo = 0.0
# αm = (2*ρ_oo - 1)/(ρ_oo + 1)
# αf = ρ_oo/(ρ_oo + 1)
# for i = 1:NT
#     @info i, "/" , NT
#     solver = NewmarkSolver(Δt, globdat, domain, αm, αf, 1e-4, 1e-6, 100)
#     # error()
#     # close("all")
#     # visσ(domain,-1.5e9, 4.5e9)
#     # savefig("Debug/$i.png")
#     # error()
#     # if i==100
#     #     close("all")
#     #     visσ(domain)
#     #     # visσ(domain,-1.5e9, 4.5e9)
#     #     savefig("Debug/test$(tid)i=75.png")
#     #     # break
#     # end
# end

# plot
close("all")
visσ(domain)
savefig("Debug/test$tid.png")

close("all")
visσ(domain)
axis("equal")
savefig("Debug/order$porder/test_stress$(tid)_$force_scale.png")

close("all")
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(ts, ux)
savefig("Debug/order$porder/test_ux$(tid)_$force_scale.png")

close("all")
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(ts, uy)
savefig("Debug/order$porder/test_uy$(tid)_$force_scale.png")


text:
#
# Test script
#
# author: Atsushi Sakai
#

module PkgTests

using Test

include("../src/trailer_hybrid_a_star.jl")

function test()
    println("Test Start !!!")

    gx = 0.0  # [m]
    gy = 0.0  # [m]
    gyaw0 = deg2rad(90.0)
    gyaw1 = deg2rad(90.0)

    ox = Float64[]
    oy = Float64[]

    for i in -25:25
        push!(ox, Float64(i))
        push!(oy, 15.0)
    end
    for i in -25:-4
        push!(ox, Float64(i))
        push!(oy, 4.0)
    end
    for i in -15:4
        push!(ox, -4.0)
        push!(oy, Float64(i))
    end
    for i in -15:4
        push!(ox, 4.0)
        push!(oy, Float64(i))
    end
    for i in 4:25
        push!(ox, Float64(i))
        push!(oy, 4.0)
    end
    for i in -4:4
        push!(ox, Float64(i))
        push!(oy, -15.0)
    end

    oox = ox[:]
    ooy = oy[:]

    sx = -10.0  # [m]
    sy = 6.0  # [m]
    syaw0 = deg2rad(00.0)
    syaw1 = deg2rad(00.0)

    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = 14.0  # [m]
    sy = 10.0  # [m]
    syaw0 = deg2rad(00.0)
    syaw1 = deg2rad(00.0)

    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = -14.0  # [m]
    sy = 12.0  # [m]
    syaw0 = deg2rad(00.0)
    syaw1 = deg2rad(00.0)
    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = -20.0  # [m]
    sy = 6.0  # [m]
    syaw0 = deg2rad(00.0)
    syaw1 = deg2rad(00.0)
    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = -14.0  # [m]
    sy = 12.0  # [m]
    syaw0 = deg2rad(00.0)
    syaw1 = deg2rad(00.0)
    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = -20.0  # [m]
    sy = 6.0  # [m]
    syaw0 = deg2rad(180.0)
    syaw1 = deg2rad(180.0)
    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION)

    @test length(path.x)>=1

    sx = -20.0  # [m]
    sy = 12.0  # [m]
    syaw0 = deg2rad(180.0)
    syaw1 = deg2rad(180.0)

    @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
                                                               trailer_hybrid_a_star.XY_GRID_RESOLUTION,
                                                               trailer_hybrid_a_star.YAW_GRID_RESOLUTION
                                                               )

    @test length(path.x)>=1

    println("Test Done !!!")

end

@time test()

end

max_stars_repo_name: yinflight/HybridAStarTrailer

max_issues_repo_name: e-schumann/HybridAStarTrailer

max_forks_repo_name: e-schumann/HybridAStarTrailer


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

include("model_AzureAsyncOperationResult.jl")
include("model_Error.jl")
include("model_ErrorDetails.jl")
include("model_Resource.jl")
include("model_SubResource.jl")


text:
#using Pkg
#pkg"activate .."
#push!(LOAD_PATH,"../src/")

using Documenter
using HawkesProcesses

const PAGES = Any[
    "Home" => "index.md",
    "Simulating" => "simulating.md",
    "Inference" => "inference.md",
    "Model Checking" => "modelchecking.md"
]

makedocs(
    sitename = "HawkesProcesses",
    authors = "Dean Markwick",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [HawkesProcesses],
    pages = PAGES
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/dm13450/HawkesProcesses.jl.git"
)


text:
module Ratios

import Base: convert, promote_rule, *, /, +, -, ^, ==, decompose

using Requires

export SimpleRatio

struct SimpleRatio{T<:Integer} <: Real
    num::T
    den::T
end

SimpleRatio(num::Integer, den::Integer) = SimpleRatio(promote(num, den)...)

convert(::Type{BigFloat}, r::SimpleRatio{S}) where {S} = BigFloat(r.num)/r.den
function convert(::Type{T}, r::SimpleRatio{S}) where {T<:AbstractFloat,S}
    P = promote_type(T,S)
    convert(T, convert(P, r.num)/convert(P, r.den))
end
SimpleRatio{T}(i::Integer) where {T<:Integer} = SimpleRatio{T}(convert(T, i), oneunit(T))
SimpleRatio{T}(r::Rational{S}) where {T<:Integer, S<:Integer} = SimpleRatio(convert(T, r.num), convert(T, r.den))
Rational{T}(r::SimpleRatio{S}) where {T<:Integer, S<:Integer} = convert(T, r.num) // convert(T, r.den)

*(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.num, x.den*y.den)
*(x::SimpleRatio, y::Bool) = SimpleRatio(x.num*y, x.den)
*(x::SimpleRatio, y::Integer) = SimpleRatio(x.num*y, x.den)
*(x::Bool, y::SimpleRatio) = SimpleRatio(x*y.num, y.den)
*(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.num, y.den)
/(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den, x.den*y.num)
/(x::SimpleRatio, y::Integer) = SimpleRatio(x.num, x.den*y)
/(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den, y.num)
+(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den + y.num, y.den)
-(x::Integer, y::SimpleRatio) = SimpleRatio(x*y.den - y.num, y.den)
+(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den + x.den*y.num, x.den*y.den)
-(x::SimpleRatio, y::SimpleRatio) = SimpleRatio(x.num*y.den - x.den*y.num, x.den*y.den)
^(x::SimpleRatio, y::Integer) = SimpleRatio(x.num^y, x.den^y)

-(x::SimpleRatio{T}) where {T<:Signed} = SimpleRatio(-x.num, x.den)
-(x::SimpleRatio{T}) where {T<:Unsigned} = throw(VERSION < v"0.7.0-DEV.1269" ? OverflowError() : OverflowError("cannot negate unsigned number"))

promote_rule(::Type{SimpleRatio{T}}, ::Type{S}) where {T<:Integer,S<:Integer} = SimpleRatio{promote_type(T,S)}
promote_rule(::Type{SimpleRatio{T}}, ::Type{SimpleRatio{S}}) where {T<:Integer,S<:Integer} = SimpleRatio{promote_type(T,S)}
promote_rule(::Type{SimpleRatio{T}}, ::Type{S}) where {T<:Integer,S<:AbstractFloat} = promote_type(T,S)
promote_rule(::Type{SimpleRatio{T}}, ::Type{Rational{S}}) where {T<:Integer,S<:Integer} = Rational{promote_type(T,S)}

==(x::SimpleRatio, y::SimpleRatio) = x.num*y.den == x.den*y.num

==(x::SimpleRatio, y::Integer) = x.num == x.den*y
==(x::Integer, y::SimpleRatio) = x*y.den == y.num

function ==(x::AbstractFloat, q::SimpleRatio)
    if isfinite(x)
        (count_ones(q.den) == 1) & (x*q.den == q.num)
    else
        x == q.num/q.den
    end
end

==(q::SimpleRatio, x::AbstractFloat) = x == q

decompose(x::SimpleRatio) = x.num, 0, x.den

function __init__()
    @require FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" begin
        using .FixedPointNumbers: FixedPoint, Fixed, Normed, rawone
        rawone_noerr(::Type{Fixed{T,f}}) where {T,f} = widen(oneunit(T)) << f
        rawone_noerr(::Type{N}) where N<:Normed = rawone(N)
        rawone_noerr(x::FixedPoint) = rawone_noerr(typeof(x))
        Base.promote_rule(::Type{SimpleRatio{S}}, ::Type{<:FixedPoint{T}}) where {S<:Integer,T<:Integer} = SimpleRatio{promote_type(S, T)}
        SimpleRatio{S}(x::FixedPoint) where S<:Integer = SimpleRatio{S}(reinterpret(x), rawone_noerr(x))
        SimpleRatio(x::FixedPoint) = SimpleRatio(reinterpret(x), rawone_noerr(x))
        Base.convert(::Type{S}, x::FixedPoint) where S<:SimpleRatio = S(x)
    end
end

end


text:
# simple example of a harmonic power flow execution
using HarmonicPowerFlow
using StatsPlots

harmonics = [h for h in 1:2:101]
# settings
settings_c = init_settings(true, harmonics)
settings_uc = init_settings(false, harmonics)

# import nodes and lines from csv files
grid = init_power_grid(
    import_nodes_from_csv("examples\\simple_coupled_hpf\\net2_buses.csv"),
    import_lines_from_csv("examples\\simple_coupled_hpf\\net2_lines.csv"),
    settings_uc)

# run harmonic power flow
u, err_h_final, n_iter_h = hpf(grid, settings_uc)

THD = HarmonicPowerFlow.THD(u)

# show barplot of harmonics at bus 4
HarmonicPowerFlow.barplot_THD(u, 4, h_max= 21)


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
using SimpleSDMLayers

coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

## Writing bug as noticed
# Read, write, re-read
ref = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
geotiff("test.tif", ref)
test = geotiff(SimpleSDMPredictor, "test.tif")

isequal(size(ref), size(test)) # not equal, 331 x 571 vs 331 x 570

isequal(ref.left, test.left) # equal
isequal(ref.bottom, test.bottom) # equal
isequal(ref.top, test.top) # equal
isequal(ref.right, test.right) # not equal

ref.right # -50.0
test.right # -50.16666666665

# Check the raster file itself
using ArchGDAL
d = ArchGDAL.read("test.tif") # correct raster size, 331 x 571
gt = ArchGDAL.getgeotransform(d)
gt[1] + gt[2]*ArchGDAL.width(d) # correct right bound, ≈ 50.0

# wcpath = joinpath(ENV["SDMLAYERS_PATH"], "WorldClim", "BioClim", "10", "wc2.1_10m_bio_1.tif")
# geotiff(SimpleSDMPredictor, wcpath; coords...)

## Tests in QGIS
# Loaded the world scale temperature data in QGIS
# Cropped it to my coordinates extent
# Dimensions are 330 x 570, not 331 x 571 as I've had for a while
d = ArchGDAL.read("../../rasters/wc1_cropped.tif")
gt = ArchGDAL.getgeotransform(d)
gt[1] + gt[2]*ArchGDAL.width(d)

nlat = (-50.0 - -145.0)/0.1666666666666666574
nlon = (75.0 - 20.0)/0.1666666666666666574

test2 = geotiff(SimpleSDMPredictor, "../../rasters/wc1_cropped.tif")
geotiff("test2.tif", test2)
geotiff(SimpleSDMPredictor, "test2.tif")

# No data value attempt
d = ArchGDAL.read("../../rasters/wc1_cropped.tif")
d1 = ArchGDAL.getband(d, 1)
d1_nodata = ArchGDAL.getnodatavalue(d1)
geotiff("test2.tif", test2; nodata=Float32(-3.4e+38))
geotiff("test2.tif", test2; nodata=Float32(d1_nodata))
geotiff(SimpleSDMPredictor, "test2.tif")
rm("test2.tif")


text:
import Base.convert
export MPSQuantumRegister, convert, two_qubit_gate_to_mpo, mpo_to_two_qubit_gate
export print_info, compress!, enaglemment_entropy, execute!

"Implementation type MPS quantum register"
struct MPSQuantumRegister{T} <: QuantumRegister
    N::Integer
    state::Array{Array{T, 3}, 1}
    s_values::Array{Union{Nothing, Array{<:AbstractFloat, 1}}, 1}
    function MPSQuantumRegister{T}(N::Integer, conf::String) where T <: Number
        @assert length(conf) == N
        new(N, [init_mps_tensor(c, T) for c in conf], fill(nothing, N-1))
    end
end

MPSQuantumRegister(N::Integer, conf::String) = MPSQuantumRegister{ComplexF64}(N, conf)
"""
    convert(::Type{FullStateQuantumRegister}, x::MPSQuantumRegister)

Convert an MPS quantum register to a full state quantum register
"""
function convert(::Type{FullStateQuantumRegister{T}}, x::MPSQuantumRegister{U}) where T <: Number where U <: Number
    state = ones(U, 1, 1)
    for i in x.N:-1:1
        @tensor state[a, b, c] := state[a, j] * x.state[i][b, c, j]
        dims = size(state)
        state = reshape(state, (dims[1] * dims[2], dims[3]))
    end
    dims = size(state)
    state = reshape(state, (dims[1] * dims[2]))
    return FullStateQuantumRegister{T}(x.N, state)
end

"""
    init_mps_tensor(qubit_value::Char, ::Type{T})

Initialise an array of MPS tensors from a string configuration
"""
function init_mps_tensor(qubit_value::Char, ::Type{T}) where T <: Number
    if qubit_value == '0'
        state = [1., 0.]
    else
        state = [0, 1]
    end
    return reshape(convert(Array{T}, state), (2, 1, 1)) # idx, up, down
end

function to_str(qreg::MPSQuantumRegister{T}) where T <: Number
    to_str(convert(FullStateQuantumRegister{T}, qreg))
end

"""
    apply_1qubit!(qreg::MPSQuantumRegister, gate, i)

Apply the 1 qubit gate to qubit i
"""
function apply_1qubit!(qreg::MPSQuantumRegister, gate, i)
    @tensor qreg.state[i][a, b, c] := qreg.state[i][m, b, c] * gate[a, m]
end

"""
    two_qubit_gate_to_mpo(gate::Array{T, 2}, n::Integer) where T <: Number

Convert a gate into an mpo
"""
function two_qubit_gate_to_mpo(gate::Array{T, 2}, threshold::AbstractFloat=1e-15) where T <: Number
    gate = reshape(gate, (2, 2, 2, 2))  #  r2, r1, l2, l1
    gate = permutedims(gate, (2, 4, 1, 3)) # r1, l1, r2, l2
    gate = reshape(gate, (4, 4)) # (r1, l1), (r2, l2)
    F = svd(gate)
    mpo = Array{Array{T, 4}, 1}()
    chi = sum(F.S .> threshold)
    push!(mpo, reshape(F.U[:,1:chi] * Diagonal(F.S[1:chi]), (2, 2, 1, chi))) # r1, l1, Up, Down
    push!(mpo, permutedims(reshape(F.Vt[1:chi,:], (chi, 2, 2, 1)), (2, 3, 1, 4))) # r2, l2, Up, Down
    mpo
end

"""
    mpo_to_two_qubit_gate(mpo::Array{Array{T, 4}, 1}) where T <: Number

Convert a gate into an mpo
"""
function mpo_to_two_qubit_gate(mpo::Array{Array{T, 4}, 1}) where T <: Number
    a = reshape_tensor(mpo[1], (1, 2, (3, 4))) # r1, l1, down
    b = reshape_tensor(mpo[2], (1, 2, (3, 4))) # r2, l2, up
    @tensor c[a, b, c, d] := a[a, b, x] * b[c, d, x] # r1, l1, r2, l2
    reshape(permutedims(c, (3, 1, 4, 2)), (4, 4))
end

"""
    apply_2qubit!(qreg::MPSQuantumRegister, gate, i, j)

Apply the 2 qubit gate to qubits i and j, assumes |i - j| = 1
"""
function apply_2qubit!(qreg::MPSQuantumRegister, gate, i, j)
    # @assert abs(i -j) == 1
    if j < i
        apply_2qubit!(qreg, swap_2qubits(gate), j, i)
    else
        # reshape gate to mpo
        gate_mpo = two_qubit_gate_to_mpo(gate)
        # remove upward facing dim 1 index to simplify contraction
        gate_mpo_1 = reshape_tensor(gate_mpo[1], [[1, 3], 2, 4])

        # remove downward facing dim 1 index to simplify contraction
        gate_mpo_2 = reshape_tensor(gate_mpo[2], [[1, 4], 2, 3])

        # tensor gate mpo with state mpos
        @tensor state_1[qubit_i, up1, down1, down2] :=  qreg.state[i][qubit_i_c, up1, down1] * gate_mpo_1[qubit_i, qubit_i_c, down2]
        @tensor state_2[qubit_j, up1, up2, down1] :=  qreg.state[j][qubit_j_c, up1, down1] * gate_mpo_2[qubit_j, qubit_j_c, up2]

        # reshape
        qreg.state[i] = reshape_tensor(state_1, [1, 2, [3, 4]])
        qreg.state[j] = reshape_tensor(state_2, [1, [2, 3], 4])

        if abs(i -j) > 1
            chi = size(gate_mpo_1)[3]
            filler_mpo = tensor_identity(chi, chi, 2)
            for k in i+1:j-1
                @tensor tmp[qubit_k, up1, up2, down1, down2] :=
                            qreg.state[k][qubit_k_c, up1, down1] *
                            filler_mpo[up2, down2, qubit_k_c, qubit_k]
                qreg.state[k] = reshape_tensor(tmp, [1, [2, 3], [4, 5]])
            end
        end
    end
    return
end

"""
    function print_info(qreg::MPSQuantumRegister)

Print the dimension of virtual bond dimensions of the given MPS
"""
function print_info(qreg::MPSQuantumRegister)
    for i = 1:qreg.N
        dims = size(qreg.state[i])[2:3]
        println("$(i) has virtual dims $(dims)")
    end
end

"""
    compress!(qreg::MPSQuantumRegister)

Compress the provided MPS quantum register
"""
function compress!(qreg::MPSQuantumRegister, threshold::AbstractFloat=1e-15)
    for i in 1:qreg.N-1
        compress_bond!(qreg, i, threshold)
    end
    for i in qreg.N-1:-1:1
        compress_bond!(qreg, i, threshold)
    end
end

function compress_bond!(qreg::MPSQuantumRegister, index::Integer, threshold::AbstractFloat=1e-15)
    A = qreg.state[index]
    A_dim = size(A)
    B = qreg.state[index+1]
    B_dim = size(B)

    # Contract virtual bond connecting A and B
    @tensor C[idx1, up, idx2, down] := A[idx1, up, c] * B[idx2, c, down]

    # Reshape to a matrix and decompose
    dims = size(C)
    C = reshape(C, (dims[1]*dims[2], dims[3]*dims[4]))
    F = svd(C)

    # Take cutoff and reshape matrices back
    s = F.S
    chi = sum(s .> threshold)
    s = s[1:chi]
    qreg.s_values[index] = s./sqrt(sum(s.^2))
    A = F.U[:,1:chi] * Diagonal(sqrt.(s))
    qreg.state[index] = reshape(A, (2, A_dim[2], chi))
    B = Diagonal(sqrt.(s)) * F.Vt[1:chi,:]
    qreg.state[index+1] = permutedims(reshape(B, (chi, 2, B_dim[3])), (2, 1, 3))
end

function entropy(s::Array{<:AbstractFloat, 1})
    -sum(s.^2 .* log.(s.^2))
end

function enaglemment_entropy(qreg::MPSQuantumRegister)
    if qreg.s_values[1] == nothing
        return
    else
        return [entropy(qreg.s_values[x]) for x in 1:qreg.N-1]
    end
end

"""
    execute!(qc::QuantumCircuit, state::MPSQuantumRegister)

Specialised execute for MPS Quantum Register
"""
function execute!(qc::Gates.QuantumCircuit, state::MPSQuantumRegister,
                  compress_freq::Integer=0)
    for (i, op) in enumerate(qc.ops)
        if op.n == 1
            apply_1qubit!(state, op.gate, op.qubits[1])
        elseif op.n == 2
            apply_2qubit!(state, op.gate, op.qubits[1], op.qubits[2])
        end
        if compress_freq > 0 && i % compress_freq == 0
            compress!(state)
        end
    end
end


text:
abstract Time_Operator
export Time_Operator
include("Smoother.jl")
include("Selector.jl")

function TProcess{N}(sst :: Array{Float64,N}, TS :: Time_Operator,
                      Dim :: Tuple{Int64,Vararg{Int64}})
  error("Still No function for pure Time_Operator");
end

function TProcess{N}(sst :: Array{Float64,N}, TS :: Time_Operator, time :: Array{Int64,1},
                      Dim :: Tuple{Int64,Vararg{Int64}} = TS.DefaultDim )
  error("Still No function for pure Time_Operator");
end

include("Simple_Functions.jl")

export Time_Operator


text:
using FancyDiagnostics
using Base.Test


text:
"""
    EGPnormal(σ, ξ, κ)

*EGPnormal* corresponds to the extended GP model developed by Gamet and Jalbert (2020), with the truncated normal distribution of mean 1 and variance 1/κ^2.

It is a three parameters family: κ controls the shape of the lower tail, σ is a scale parameter, and ξ controls the rate of upper tail decay.

```julia
EGPnormal(σ, ξ, κ)   # EGP of Gamet and Jalbert (2020) with scale parameter σ, rate of upper tail decay ξ and shape of the lower tail κ.

params(d)           # Get the parameters, i.e. (σ, ξ, κ)
```

Reference :

* Unpublished for now
"""
struct EGPnormal{T<:Real} <: ContinuousUnivariateDistribution
    σ::T    # scale parameter
    ξ::T    # rate of upper tail decay
    κ::T    # shape of the lower tail

    function EGPnormal{T}(σ::T, ξ::T, κ::T) where {T <: Real}
        new{T}(σ, ξ, κ)
    end
end

function EGPnormal(σ::T, ξ::T, κ::T; check_args=true) where {T <: Real}
    check_args && @check_args(EGPnormal, σ > zero(σ) && κ >= zero(κ))
    return EGPnormal{T}(σ, ξ, κ)
end

EGPnormal(σ::Real, ξ::Real, κ::Real) = EGPnormal(promote(σ, ξ, κ)...)
EGPnormal(σ::Integer, ξ::Integer, κ::Integer) = EGPnormal(float(σ), float(ξ), float(κ))

minimum(d::EGPnormal) = 0.0
maximum(d::EGPnormal) = Inf * (d.ξ >= 0) - ( d.σ / d.ξ ) * (d.ξ < 0)
insupport(d::EGPnormal, x::Real) = minimum(d) <= x <= maximum(d)

#### Parameters

scale(d::EGPnormal) = d.σ
decay(d::EGPnormal) = d.ξ    # noms à vérifier...
shape(d::EGPnormal) = d.κ

params(d::EGPnormal) = (d.σ, d.ξ, d.κ)
partype(::EGPnormal{T}) where {T} = T

#### Evaluation

function pdf(d::EGPnormal{T}, x::Real) where T<:Real

    (σ, ξ, κ) = params(d)
    η = 1

    if κ == 0

        return Distributions.pdf.(GeneralizedPareto(0, σ, ξ), x)

    else

        return κ/(Distributions.cdf(Normal(), κ*(1-η))-Distributions.cdf(Normal(), -κ*η))*Distributions.pdf(Normal(), κ*(Distributions.cdf(GeneralizedPareto(0, σ, ξ),x)-η))*Distributions.pdf(GeneralizedPareto(0, σ, ξ),x)

    end
end

logpdf(d::EGPnormal, x::Real) = log(pdf(d, x))

function cdf(d::EGPnormal{T}, x::Real) where T<:Real

    (σ, ξ, κ) = params(d)

    return 2/Distributions.erf(κ/sqrt(2)) * (Distributions.cdf(Normal(), κ*(Distributions.cdf(GeneralizedPareto(σ, ξ), x)-1)) - 1 + Distributions.cdf(Normal(), κ))

end

logcdf(d::EGPnormal, x::Real) = log(cdf(d, x))

function quantile(d::EGPnormal{T}, p::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    (σ, ξ, κ) = params(d)

    return Distributions.quantile(GeneralizedPareto(σ,ξ), 1 + 1/κ * Distributions.quantile(Normal(), 0.5*(1 - (1-p)*Distributions.erf(κ/sqrt(2)))))

end

#### Sampling

function rand(rng::AbstractRNG, d::EGPnormal)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    if u == 1

        return maximum(d)

    else

        return quantile(d, u)

    end
end

sampler(d::EGPnormal) = d


text:
import IntervalArithmetic

"""
    IntervalAbs(lo::Real, hi::Real)

Interval abstraction for real-valued numbers.
"""
struct IntervalAbs{T <: Real}
    interval::IntervalArithmetic.Interval{T}
end

IntervalAbs{T}(x::Integer) where {T <: Integer} =
    IntervalAbs{T}(T(x), T(x)+1)
IntervalAbs{T}(x::Union{AbstractFloat,Integer}) where {T <:AbstractFloat} =
    IntervalAbs{T}(T(x), T(x))
IntervalAbs{T}(a::T, b::T) where {T <: Real} =
    IntervalAbs{T}(IntervalArithmetic.Interval{T}(a, b))

Base.copy(x::IntervalAbs) = x

Base.convert(I::Type{IntervalAbs}, x::Real) = I(x)
Base.convert(I::Type{IntervalAbs{T}}, x::Integer) where {T <: Real} = I(x)

Base.show(io::IO, x::IntervalAbs) = Base.show(io, x.interval)

lub(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval ∪ b.interval)
glb(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval ∩ b.interval)

widen(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval ∪ b.interval)

Base.:(==)(a::IntervalAbs, b::IntervalAbs) = a.interval == b.interval

Base.:+(a::IntervalAbs) = a
Base.:-(a::IntervalAbs) = IntervalAbs(-a)

Base.:+(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval + b.interval)
Base.:-(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval - b.interval)
Base.:*(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval * b.interval)
Base.:/(a::IntervalAbs, b::IntervalAbs) = IntervalAbs(a.interval / b.interval)

Base.:+(a::IntervalAbs, b::Real) = IntervalAbs(a.interval + b)
Base.:-(a::IntervalAbs, b::Real) = IntervalAbs(a.interval - b)
Base.:*(a::IntervalAbs, b::Real) = IntervalAbs(a.interval * b)
Base.:/(a::IntervalAbs, b::Real) = IntervalAbs(a.interval / b)

Base.:+(a::Real, b::IntervalAbs) = IntervalAbs(a + b.interval)
Base.:-(a::Real, b::IntervalAbs) = IntervalAbs(a - b.interval)
Base.:*(a::Real, b::IntervalAbs) = IntervalAbs(a * b.interval)
Base.:/(a::Real, b::IntervalAbs) = IntervalAbs(a / b.interval)

equiv(a::IntervalAbs, b::IntervalAbs) =
    !isdisjoint(a.interval, b.interval)
nequiv(a::IntervalAbs, b::IntervalAbs) =
    isdisjoint(a.interval, b.interval)
Base.:<(a::IntervalAbs, b::IntervalAbs) =
    !IntervalArithmetic.precedes(b.interval, a.interval)
Base.:<=(a::IntervalAbs, b::IntervalAbs) =
    !IntervalArithmetic.strictprecedes(b.interval, a.interval)

equiv(a::IntervalAbs, b::Real) = b ⊆ a.interval
nequiv(a::IntervalAbs, b::Real) = !(a.interval.lo == a.interval.hi == b)
Base.:<(a::IntervalAbs, b::Real) = a.interval.lo < b
Base.:<=(a::IntervalAbs, b::Real) = a.interval.lo <= b

equiv(a::Real, b::IntervalAbs) = a ⊆ b.interval
nequiv(a::Real, b::IntervalAbs) = !(a == b.interval.lo == b.interval.hi)
Base.:<(a::Real, b::IntervalAbs) = a < b.interval.hi
Base.:<=(a::Real, b::IntervalAbs) = a <= b.interval.hi


text:
struct DecimalType <: LogicalType
    type::String
    logicalType::String
    size::Int64
    precision::Int64
    scale::Int64
end

DecimalType(type::String, logicalType::String, size::Union{Nothing, Int64}, precision::Int64, scale::Union{Nothing, Int64}) = DecimalType(type, logicalType, something(size, 16), precision, something(scale, 0))
DecimalType(precision::Integer, scale::Union{Nothing, Integer}) = DecimalType("fixed", "decimal", 16, precision, something(scale, 0))

const FIXEDDEC = FixedType("decimal", 16)

struct Decimal{S, P}
    value::Int128
end

juliatype(x::DecimalType) = Decimal{x.scale, x.precision}
schematype(::Type{Decimal{S, P}}) where {S, P} = DecimalType(P, S)

function writevalue(B::Binary, D::DecimalType, x::Decimal, buf, pos, len, opts)
    @assert D.type == "fixed"
    return writevalue(B, FIXEDDEC, _cast(NTuple{16, UInt8}, x.value), buf, pos, len, opts)
end

nbytes(::DecimalType, x::Decimal) = 16

inttype(sz) = sz == 1 ? Int8 : sz == 2 ? Int16 : sz == 4 ? Int32 : sz == 8 ? Int64 : Int128

function readvalue(B::Binary, D::DecimalType, ::Type{Decimal{S, P}}, buf, pos, len, opts) where {S, P}
    if D.type == "fixed"
        x, pos = readvalue(B, FixedType("", D.size), NTuple{D.size, UInt8}, buf, pos, len, opts)
        return Decimal{S, P}(Int128(_cast(inttype(D.size), x))), pos
    elseif D.type == "bytes"
        b, pos = readvalue(B, bytes, Vector{UInt8}, buf, pos, len, opts)
        resize!(b, 16)
        return Decimal{S, P}(reinterpret(Int128, b)[1]), pos
    else
        error("invalid Decimal type = $(D.type)")
    end
end

function skipvalue(B::Binary, D::DecimalType, ::Type{Decimal{S, P}}, buf, pos, len, opts) where {S, P}
    if D.type == "fixed"
        return skipvalue(B, FixedType("", D.size), NTuple{D.size, UInt8}, buf, pos, len, opts)
    elseif D.type == "bytes"
        return skipvalue(B, bytes, Vector{UInt8}, buf, pos, len, opts)
    else
        error("invalid Decimal type = $(D.type)")
    end
end

struct UUIDType <: LogicalType
    type::String
    logicalType::String
end
UUIDType() = UUIDType("string", "uuid")

juliatype(::UUIDType) = UUID
schematype(::Type{UUID}) = UUIDType()

writevalue(B::Binary, ::UUIDType, x::UUID, buf, pos, len, opts) =
    writevalue(B, string, Base.string(x), buf, pos, len, opts)

nbytes(::UUIDType, ::UUID) = nbytes(36) + 36

function readvalue(B::Binary, ::UUIDType, ::Type{UUID}, buf, pos, len, opts)
    x, pos = readvalue(B, string, String, buf, pos, len, opts)
    return UUID(x), pos
end

skipvalue(B::Binary, ::UUIDType, ::Type{UUID}, buf, pos, len, opts) =
    skipvalue(B, string, String, buf, pos, len, opts)

struct DateType <: LogicalType
    type::String
    logicalType::String
end
DateType() = DateType("int", "date")

juliatype(::DateType) = Date
schematype(::Type{Date}) = DateType()

const UNIX_EPOCH_DATE = Dates.value(Dates.Date(1970))
unix(x::Date) = Int32(Dates.value(x) - UNIX_EPOCH_DATE)
date(x::Int32) = Dates.Date(Dates.UTD(Int64(x + UNIX_EPOCH_DATE)))

writevalue(B::Binary, ::DateType, x::Date, buf, pos, len, opts) =
    writevalue(B, int, unix(x), buf, pos, len, opts)

nbytes(::DateType, x::Date) = nbytes(unix(x))

function readvalue(B::Binary, ::DateType, ::Type{Date}, buf, pos, len, opts)
    x, pos = readvalue(B, int, Int32, buf, pos, len, opts)
    return date(x), pos
end

skipvalue(B::Binary, ::DateType, ::Type{Date}, buf, pos, len, opts) =
    skipvalue(B, int, Int32, buf, pos, len, opts)

struct TimeMillisType <: LogicalType
    type::String
    logicalType::String
end
TimeMillisType() = TimeMillisType("int", "time-millis")

juliatype(::TimeMillisType) = Time

writevalue(B::Binary, ::TimeMillisType, x::Time, buf, pos, len, opts) =
    writevalue(B, long, div(Dates.value(x), 1000_000), buf, pos, len, opts)

nbytes(::TimeMillisType, x::Time) = nbytes(div(Dates.value(x), 1000_000))

function readvalue(B::Binary, ::TimeMillisType, ::Type{Time}, buf, pos, len, opts)
    x, pos = readvalue(B, int, Int32, buf, pos, len, opts)
    return Time(Nanosecond(Millisecond(x))), pos
end

skipvalue(B::Binary, ::TimeMillisType, ::Type{Time}, buf, pos, len, opts) =
    skipvalue(B, int, Int32, buf, pos, len, opts)

struct TimeMicrosType <: LogicalType
    type::String
    logicalType::String
end
TimeMicrosType() = TimeMicrosType("long", "time-micros")

juliatype(::TimeMicrosType) = Time
schematype(::Type{Time}) = TimeMicrosType()

writevalue(B::Binary, ::TimeMicrosType, x::Time, buf, pos, len, opts) =
    writevalue(B, long, div(Dates.value(x), 1000), buf, pos, len, opts)

nbytes(::TimeMicrosType, x::Time) = nbytes(div(Dates.value(x), 1000))

function readvalue(B::Binary, ::TimeMicrosType, ::Type{Time}, buf, pos, len, opts)
    x, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return Time(Nanosecond(Microsecond(x))), pos
end

skipvalue(B::Binary, ::TimeMicrosType, ::Type{Time}, buf, pos, len, opts) =
    skipvalue(B, long, Int64, buf, pos, len, opts)

struct TimestampMillisType <: LogicalType
    type::String
    logicalType::String
end
TimestampMillisType() = TimestampMillisType("long", "timestamp-millis")

juliatype(::TimestampMillisType) = DateTime

const UNIX_EPOCH_DATETIME = Dates.value(Dates.DateTime(1970))
unix(x::DateTime) = Int64(Dates.value(x) - UNIX_EPOCH_DATETIME)
datetime(x::Int64) = Dates.DateTime(Dates.UTM(Int64(x + UNIX_EPOCH_DATETIME)))

writevalue(B::Binary, ::TimestampMillisType, x::DateTime, buf, pos, len, opts) =
    writevalue(B, long, unix(x), buf, pos, len, opts)

nbytes(::TimestampMillisType, x::DateTime) = nbytes(unix(x))

function readvalue(B::Binary, ::TimestampMillisType, ::Type{DateTime}, buf, pos, len, opts)
    x, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return datetime(x), pos
end

skipvalue(B::Binary, ::TimestampMillisType, ::Type{DateTime}, buf, pos, len, opts) =
    skipvalue(B, long, Int64, buf, pos, len, opts)

struct TimestampMicrosType <: LogicalType
    type::String
    logicalType::String
end
TimestampMicrosType() = TimestampMicrosType("long", "timestamp-micros")

juliatype(::TimestampMicrosType) = DateTime

writevalue(B::Binary, ::TimestampMicrosType, x::DateTime, buf, pos, len, opts) =
    writevalue(B, long, unix(x) * 1000, buf, pos, len, opts)

nbytes(::TimestampMicrosType, x::DateTime) = nbytes(unix(x) * 1000)

function readvalue(B::Binary, ::TimestampMicrosType, ::Type{DateTime}, buf, pos, len, opts)
    x, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return datetime(div(x, 1000)), pos
end

skipvalue(B::Binary, ::TimestampMicrosType, ::Type{DateTime}, buf, pos, len, opts) =
    skipvalue(B, long, Int64, buf, pos, len, opts)

struct LocalTimestampMillisType <: LogicalType
    type::String
    logicalType::String
end
LocalTimestampMillisType() = LocalTimestampMillisType("long", "local-timestamp-millis")

juliatype(::LocalTimestampMillisType) = DateTime
schematype(::Type{DateTime}) = LocalTimestampMillisType()

writevalue(B::Binary, ::LocalTimestampMillisType, x::DateTime, buf, pos, len, opts) =
    writevalue(B, long, unix(x), buf, pos, len, opts)

nbytes(::LocalTimestampMillisType, x::DateTime) = nbytes(unix(x))

function readvalue(B::Binary, ::LocalTimestampMillisType, ::Type{DateTime}, buf, pos, len, opts)
    x, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return datetime(x), pos
end

skipvalue(B::Binary, ::LocalTimestampMillisType, ::Type{DateTime}, buf, pos, len, opts) =
    skipvalue(B, long, Int64, buf, pos, len, opts)

struct LocalTimestampMicrosType <: LogicalType
    type::String
    logicalType::String
end
LocalTimestampMicrosType() = LocalTimestampMicrosType("long", "local-timestamp-micros")

juliatype(::LocalTimestampMicrosType) = DateTime

writevalue(B::Binary, ::LocalTimestampMicrosType, x::DateTime, buf, pos, len, opts) =
    writevalue(B, long, unix(x) * 1000, buf, pos, len, opts)

nbytes(::LocalTimestampMicrosType, x::DateTime) = nbytes(unix(x) * 1000)

function readvalue(B::Binary, ::LocalTimestampMicrosType, ::Type{DateTime}, buf, pos, len, opts)
    x, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return datetime(div(x, 1000)), pos
end

skipvalue(B::Binary, ::LocalTimestampMicrosType, ::Type{DateTime}, buf, pos, len, opts) =
    skipvalue(B, long, Int64, buf, pos, len, opts)

struct DurationType <: LogicalType
    type::String
    logicalType::String
    size::Int64
end
DurationType() = DurationType("fixed", "duration", 12)

struct Duration
    months::Int32
    days::Int32
    millis::Int32
end

juliatype(::DurationType) = Duration
schematype(::Type{Duration}) = DurationType()
const FIXEDDUR = FixedType("duration", 12)

function writevalue(B::Binary, D::DurationType, x::Duration, buf, pos, len, opts)
    @assert D.type == "fixed"
    return writevalue(B, FIXEDDUR, _cast(NTuple{12, UInt8}, x), buf, pos, len, opts)
end

nbytes(::DurationType, x::Duration) = 12

function readvalue(B::Binary, ::DurationType, ::Type{Duration}, buf, pos, len, opts)
    x, pos = readvalue(B, FIXEDDUR, NTuple{12, UInt8}, buf, pos, len, opts)
    months, days, millis = _cast(NTuple{3, Int32}, x)
    return Duration(months, days, millis), pos
end

skipvalue(B::Binary, ::DurationType, ::Type{Duration}, buf, pos, len, opts) =
    skipvalue(B, FIXEDDUR, NTuple{12, UInt8}, buf, pos, len, opts)


text:
# main.jl
# Code to generate various practice flashcards
# for Julia's || and && operators.

# Change this if you want to generate more or fewer flashcards. There are 127 possibilities in total.
NUM_FLASHCARDS = 30

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

operators = (:(||), :(&&))
bools = (true, false)

# It's good practice to preallocate the arrays.

A = Array{Bool}(undef, NUM_FLASHCARDS, 4)
B = Array{Symbol}(undef, NUM_FLASHCARDS, 3)

# Now we fill them with randomly-selected members from the
# tuples above.

Random.rand!(A, bools)
Random.rand!(B, operators)

# Finally, let's make the flashcards themselves. S will have
# the question slide, and E will have the answer slide.

S = Array{String}(undef, NUM_FLASHCARDS)
E = Array{Bool}(undef, NUM_FLASHCARDS)

for i in 1:size(S)[1]
    S[i] = "$(A[i,1]) $(B[i,1]) $(A[i,2]) $(B[i,2]) $(A[i,3]) $(B[i,3]) $(A[i,4])"
end

@. E = eval(Meta.parse(S))

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

open("..\\Julia-Short-Circuit-Flashcard-Generators\\generated_flashcards.txt", "w") do io
    for i in 1:size(S)[1]
        print(io, "$(S[i])\n$(E[i])\n\n")
    end
end

println("Flashcards generated.~")
# Good work, my friends.

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


text:
using Base.Profile
n = 100
function f(x)
    A = reshape(x[1:n^2], n, n)
    B = reshape(x[n^2 + 1:2n^2], n, n)
    c = x[2n^2+1:end]
    trace(log(A * B .+ c))
end

@show N = 2n^2 + n
x = rand(N)

include("prototype6.jl")

∇f = grad(f, N)

stor = similar(x)
∇f(x, stor)

Profile.clear_malloc_data()
∇f(x, stor)