text:
export BasicDQNLearner

using StatsBase: mean
using Zygote

"""
    BasicDQNLearner(;kwargs...)

See paper: [Playing Atari with Deep Reinforcement Learning](https://www.cs.toronto.edu/~vmnih/docs/dqn.pdf)

This is the very basic implementation of DQN. Compared to the traditional Q learning, the only difference is that,
in the updating step it uses a batch of transitions sampled from an experience buffer instead of current transition.
And the `approximator` is usually a [`NeuralNetworkQ`](@ref).

You can start from this implementation to understand how everything is organized and how to write your own customized algorithm.

# Keywords

- `approximator`::[`AbstractQApproximator`](@ref): used to get Q-values of a state.
- `loss_fun`: the loss function to use. TODO: provide a default [`huber_loss`](@ref)?
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
"""
Base.@kwdef mutable struct BasicDQNLearner{Tq<:AbstractQApproximator,Tf} <: AbstractLearner
    approximator::Tq
    loss_fun::Tf
    γ::Float32 = 0.99f0
    batch_size::Int = 32
    update_horizon::Int = 1
    min_replay_history::Int = 32
end

function update!(learner::BasicDQNLearner{<:NeuralNetworkQ}, batch)
    Q, γ, loss_fun, update_horizon, batch_size = learner.approximator,
        learner.γ,
        learner.loss_fun,
        learner.update_horizon,
        learner.batch_size
    states, rewards, terminals, next_states = map(x->to_device(Q, x), (batch.states, batch.rewards, batch.terminals, batch.next_states))
    actions = CartesianIndex.(batch.actions, 1:batch_size)

    gs = gradient(Q) do
        q = batch_estimate(Q, states)[actions]
        q′ = vec(maximum(batch_estimate(Q, next_states); dims = 1))
        G = rewards .+ γ^update_horizon .* (1 .- terminals) .* q′

        batch_losses = loss_fun(G, q)
        mean(batch_losses)
    end

    update!(Q, gs)
end

function extract_transitions(buffer::CircularTurnBuffer{RTSA}, learner::BasicDQNLearner)
    if length(buffer) > learner.min_replay_history
        inds = sample_indices(buffer, learner.batch_size, learner.update_horizon, nothing)
        extract_SARTS(buffer, inds, learner.γ, learner.update_horizon, nothing)
    else
        nothing
    end
end


text:
module TSML
using Reexport

# reexport common functions to Main
@reexport using CSV
@reexport using Dates
@reexport using DataFrames
@reexport using Random
@reexport using Statistics

include("types.jl")
using .TSMLTypes

include("utils.jl")
using .Utils

include("basefilters.jl")
using .BaseFilters

include("valdatefilters.jl")
using .ValDateFilters

include("baseline.jl")
using .BaselineAlgos

include("mlbase.jl")
using .MLBaseWrapper

include("decisiontree.jl")
using .DecisionTreeLearners

include("statifier.jl")
using .Statifiers

include("monotonicer.jl")
using .Monotonicers

include("cliwrapper.jl")
using .CLIWrappers

include("tsclassifier.jl")
using .TSClassifiers

include("outliernicer.jl")
using .Outliernicers

include("plotter.jl")
using .Plotters

include("timescaledb.jl")
using .TimescaleDBs

include("demo.jl")
using .TSMLDemo

include("argparse.jl")
using .ArgumentParsers

include("ensemble.jl")
using .EnsembleMethods

end # module

ext: jl


text:
using GenericInstruments
GI = GenericInstruments
####################################
resmgr = GI.viOpenDefaultRM()
# Instantiate obj
awg1 = GI.AWG.INSTR(:Agilent33500B, "USB0::0x0957::0x2C07::MY52803073::INSTR")
GI.AWG.connect!(resmgr, awg1)
# IDN
GI.AWG.get_idn(awg1)
#or
GI.AWG.query(awg1,"*IDN?")

# Reset
GI.AWG.reset_instr(awg1)
# Output load INF
GI.AWG.set_load(awg1,ch="1",load="INF")
# Wfm
GI.AWG.set_wfm(awg1,ch="1",func="square")
# Amplitude
GI.AWG.set_amplit(awg1,ch="1",unit="vpp",volt=1)
# Offset
GI.AWG.set_offs(awg1,ch="1",volt=0.4)
# High level
GI.AWG.set_hilev(awg1,ch="1",volt=3)
# Low level
GI.AWG.set_lolev(awg1,ch="1",volt=1)
# Duty cycle
GI.AWG.set_duty(awg1,ch="1",func="square",duty=22)
# Arbitrary waveform
GI.AWG.clear_arbs(awg1,ch="1")
wfm = Float32.([0,0, 0,1, 1,1,1,-1,-1,-1,-1,-1,1,1,1, 0, 0,0,0,0]) # normed to -1:1
GI.AWG.send_arb(awg1,ch="1",wfm=wfm,dt=1e-7,arb="mywfm")
GI.AWG.set_wfm(awg1,ch="1",func="ARB")
GI.AWG.set_arb_wfm(awg1,ch="1",arb="mywfm")
GI.AWG.set_burst_stat(awg1,ch="1",st="on")
GI.AWG.set_trig_sour(awg1,ch="1",sour="BUS")
GI.AWG.set_outp(awg1,ch="1",st="on")
GI.AWG.send_soft_trig(awg1)

GI.AWG.reset_instr(awg1)

GI.AWG.disconnect!(awg1)
# pending...


text:
"""
FunctionTransform
```
    f(x) = abs.(x)
    tr = FunctionTransform(f)
```
Take a function `f` as an argument which is going to act on each vector individually.
Make sure that `f` is supposed to act on a vector by eventually using broadcasting
For example `f(x)=sin(x)` -> `f(x)=sin.(x)`
"""
struct FunctionTransform{F} <: Transform
    f::F
end

transform(t::FunctionTransform,X::T,obsdim::Int=defaultobs) where {T} = mapslices(t.f,X,dims=obsdim)

params(t::FunctionTransform) = t.f


text:
using ITensors
using ITensors.ITensorNetworkMaps

function transfer_matrix(ψ::MPS, n::Int)
  ψ′ = prime(linkinds, ψ)
  lⁿ = linkinds(ψ, n)
  lⁿ⁻¹ = linkinds(ψ, n - 1)
  input_inds_n = (lⁿ..., dag(lⁿ')...)
  output_inds_n = (lⁿ⁻¹..., dag(lⁿ⁻¹')...)
  return ITensorNetworkMap(
    [ψ[n], dag(ψ′[n])]; input_inds=input_inds_n, output_inds=output_inds_n
  )
end

function transfer_matrix(Hψ::Tuple{MPO,MPS}, n::Int)
  H, ψ = Hψ
  ψ′ = ψ'
  lⁿ = linkinds(ψ, n)
  lⁿ⁻¹ = linkinds(ψ, n - 1)
  hⁿ = linkinds(H, n)
  hⁿ⁻¹ = linkinds(H, n - 1)
  input_inds_n = (lⁿ..., hⁿ..., dag(lⁿ')...)
  output_inds_n = (lⁿ⁻¹..., hⁿ⁻¹..., dag(lⁿ⁻¹')...)
  return ITensorNetworkMap(
    [ψ[n], H[n], dag(ψ′[n])]; input_inds=input_inds_n, output_inds=output_inds_n
  )
end

function transfer_matrix(ψ, r::UnitRange)
  return prod([transfer_matrix(ψ, n) for n in r])
end

function cache(T::Vector{<:ITensorNetworkMap})
  v = Dict()
  v[0 => 1] = ITensor(1)
  for n in 1:N
    v[n => n + 1] = v[n - 1 => n] * T[n]
  end
  v[N + 1 => N] = ITensor(1)
  for n in reverse(1:N)
    v[n => n - 1] = T[n] * v[n + 1 => n]
  end
  return v
end

function hessian(H::MPO, cache, n::Int)
  return hessian(H, cache, n:n)
end

function hessian(H::MPO, cache, n⃗::UnitRange)
  n1 = first(n⃗)
  n2 = last(n⃗)
  Hᴸ = vᴴ[n1 - 1 => n1]
  Hᴿ = vᴴ[n2 + 1 => n2]
  Hⁿ⃗ = H[n⃗]
  return ITensorNetworkMap([Hᴸ, Hⁿ⃗..., Hᴿ])
end


text:
using LinearAlgebra
using CartesianGrids
using RigidBodyTools

Δx = 0.02
xlim = (-5.98,5.98)
ylim = (-5.98,5.98)
g = PhysicalGrid(xlim,ylim,Δx)

w = Nodes(Dual,size(g))
q = Edges(Dual,w)
oc = similar(w)
oc .= 1

radius = 1.0
body = Circle(radius,1.5*Δx)
X = VectorData(collect(body))
ds = areas(body)

f = VectorData(X)
ϕ = ScalarData(f)
ϕ .= 1

os = similar(ϕ)
os .= 1

@testset "Inner Products" begin
   w .= 1
   u = similar(w)
   u .= 1
   @test sqrt(dot(w,u,g)) == norm(w,g)

   u = Nodes(Dual,5,5)
   @test_throws MethodError dot(w,u,g)
   @test_throws AssertionError dot(u,u,g)

   @test norm(w,g) == sqrt((xlim[2]-xlim[1])*(ylim[2]-ylim[1]))

   @test isapprox(dot(os,os,ds),2π,atol=1e-3)

   nrm = normals(body)
   @test isapprox(integrate(pointwise_dot(nrm,nrm),areas(body)),2π,atol=1e-3)

end

bl = BodyList()
push!(bl,deepcopy(body))
push!(bl,deepcopy(body))

@testset "Operations on body lists" begin
  Xl = points(bl)

  @test dot(Xl.u,Xl.u,areas(bl),bl,1) == dot(Xl.u,Xl.u,areas(bl),bl,2) == dot(X.u,X.u,ds)

  @test dot(Xl,Xl,areas(bl),bl,1) == dot(Xl,Xl,areas(bl),bl,2) == dot(X,X,ds)

  @test norm(Xl.u,areas(bl),bl,1) == norm(Xl.u,areas(bl),bl,2) == norm(X.u,ds)

  @test norm(Xl,areas(bl),bl,1) == norm(Xl,areas(bl),bl,2) == norm(X,ds)

  nrm = normals(bl)
  @test isapprox(integrate(pointwise_dot(nrm,nrm),areas(bl),bl,1),2π,atol=1e-3)
  @test isapprox(integrate(pointwise_dot(nrm,nrm),areas(bl),bl,2),2π,atol=1e-3)

  fl = ScalarData(Xl)
  gl = ScalarData(Xl)
  gl .= rand(length(gl))

  copyto!(fl,gl,bl,2)

  @test iszero(view(fl,bl,1))
  @test view(fl,bl,2) == view(gl,bl,2)

  # test the version that copies a vector of just the right size
  fl = ScalarData(Xl)
  copyto!(fl,view(gl,bl,2),bl,2)

  @test iszero(view(fl,bl,1))
  @test view(fl,bl,2) == view(gl,bl,2)

  @test_throws AssertionError copyto!(fl,[1.0,2.0,3.0],bl,2)

  fl = VectorData(Xl)
  gl = VectorData(Xl)
  gl .= rand(length(gl))

  copyto!(fl,gl,bl,1)
  @test iszero(view(fl.u,bl,2)) && iszero(view(fl.v,bl,2))
  @test view(fl,bl,1) == view(gl,bl,1)

end

@testset "Regularization and Interpolation" begin
  reg = Regularize(X,Δx,I0=origin(g),weights=dlengthmid(body),ddftype=CartesianGrids.Yang3)

  Rc = RegularizationMatrix(reg,ϕ,oc)
  Ec = InterpolationMatrix(reg,oc,ϕ)

  @test isapprox(abs(dot(oc,Rc*os,g) - dot(os,os,ds)),0,atol=1e-14)

  u = similar(w)
  u .= randn(size(u))
  ϕ .= randn(size(ϕ))
  @test isapprox(dot(u,Rc*ϕ,g),dot(Ec*u,ϕ,ds),atol=1e-14)

  Rf = RegularizationMatrix(reg,f,q);
  Ef = InterpolationMatrix(reg,q,f);

  f.u .= 1;
  q.u .= 1;
  @test isapprox(abs(dot(q.u,(Rf*f).u,g)),sum(ds),atol=1e-14)
  @test isapprox(abs(dot(q.u,(Rf*f).u,g)),2π,atol=1e-3)

end


text:
module NewtonMethod2

using LinearAlgebra, Statistics, ForwardDiff

function newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)
    # setup the algorithm
    x_old = x0
    normdiff = Inf
    iter = 0
    while normdiff > tol && iter <= maxiter
        x_new = x_old-f(x_old)/f_prime(x_old) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter <= maxiter
        return (x_old, normdiff, iter)
    else
        return (nothing, normdiff, iter) # hope it is correct!
    end
end

export newtonroot

#########

function newtonroot2(f; x0, tol=1E-7, maxiter=1000)
    D(f) = x -> ForwardDiff.derivative(f, x)
    f_prime = D(f)
    # setup the algorithm
    x_old = x0
    normdiff = Inf
    iter = 0
    while normdiff > tol && iter <= maxiter
        x_new = x_old-f(x_old)/f_prime(x_old) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter <= maxiter
        return (x_old, normdiff, iter)
    else
        return (nothing, normdiff, iter)
    end
end

export(newtonroot2)

end


text:
# Through a little deduction, you should now be able to determine the remaining
# digits. Consider again the first example above:
#
# acedgfb cdfbe gcdfa fbcad dab cefabd cdfgeb eafb cagedb ab |
# cdfeb fcadb cdfeb cdbaf
#
# After some careful analysis, the mapping between signal wires and segments
# only make sense in the following configuration:
#
#  dddd
# e    a
# e    a
#  ffff
# g    b
# g    b
#  cccc
#
# So, the unique signal patterns would correspond to the following digits:
#
# acedgfb: 8
# cdfbe: 5
# gcdfa: 2
# fbcad: 3
# dab: 7
# cefabd: 9
# cdfgeb: 6
# eafb: 4
# cagedb: 0
# ab: 1
#
# Then, the four digits of the output value can be decoded:
#
# cdfeb: 5
# fcadb: 3
# cdfeb: 5
# cdbaf: 3
#
# Therefore, the output value for this entry is 5353.
#
# Following this same process for each entry in the second, larger example
# above, the output value of each entry can be determined:
#
# fdgacbe cefdb cefbgd gcbe: 8394
# fcgedb cgb dgebacf gc: 9781
# cg cg fdcagb cbg: 1197
# efabcd cedba gadfec cb: 9361
# gecf egdcabf bgf bfgea: 4873
# gebdcfa ecba ca fadegcb: 8418
# cefg dcbef fcge gbcadfe: 4548
# ed bcgafe cdgba cbgef: 1625
# gbdfcae bgc cg cgb: 8717
# fgae cfgab fg bagce: 4315
#
# Adding all of the output values in this larger example produces 61229.
#
# For each entry, determine all of the wire/segment connections and decode the
# four-digit output values. What do you get if you add up all of the output
# values?

digits = Array{String}(undef, 10)
cnt = sum(eachline("input.txt")) do l
  (input, output) = split.(split(l, " | "))
  inputlens = length.(input)

  # find the unique digits - idxs is a vector of indexes that would sort the
  # inputlens. Since all 10 digits must be in the input, we know exactly where
  # to find 1, 4, 7, and 8 because inputlens[idxs] will always be
  # [2, 3, 4, 5, 5, 5, 6, 6, 6, 7]
  idxs = sortperm(inputlens)
  digits[2] = input[idxs[1]]
  digits[5] = input[idxs[3]]
  digits[8] = input[idxs[2]]
  digits[9] = input[idxs[10]]

  fivesegments = input[idxs[4:6]]
  sixsegments = input[idxs[7:9]]

  # of the 5-segment digits, only 3 shares all segments in common with 7
  threeidx = findfirst(x -> issubset(digits[8], x), fivesegments)
  digits[4] = fivesegments[threeidx]
  deleteat!(fivesegments, threeidx)

  # of the 6-segment digits, 6 is the only one that does NOT share all segments
  # in common with 1
  sixidx = findfirst(x -> !issubset(digits[2], x), sixsegments)
  digits[7] = sixsegments[sixidx]
  deleteat!(sixsegments, sixidx)

  # of the remaining 5-segment digits, only 5 will share all segments with 6
  fiveidx = findfirst(x -> issubset(x, digits[7]), fivesegments)
  digits[6] = fivesegments[fiveidx]
  deleteat!(fivesegments, fiveidx)

  # the remaining 5-segment number must be 2
  digits[3] = fivesegments[1]

  # of the remaining 6-segment numbers, only 9 shares all segments with 4
  nineidx = findfirst(x -> issubset(digits[5], x), sixsegments)
  digits[10] = sixsegments[nineidx]
  deleteat!(sixsegments, nineidx)

  # the remaining 6-segment number is 0
  digits[1] = sixsegments[1]

  # now decode output
  mapfoldl((a,x) -> a*10+x, output) do x
    # digit is the index of the entry in digits, minus 1
    findfirst(y -> issetequal(x, y), digits) - 1
  end
end

println("Result: ", cnt)


text:
"""
    wdist(d; h = 2, cri = 4, squared = false)
Compute weights from distances, using a decreasing exponential function.
* `d` : A vector of distances.
* `h` : A scaling positive scalar defining the shape of the function.
* `cri` : A positive scalar defining outliers in the distances vector.
* `squared`: If true, distances are replaced by the squared distances;
    the weight function is then a Gaussian (RBF) kernel function.

Weights are computed by exp(-d / (h * MAD(d))), or are set to 0 for
distances > Median(d) + cri * MAD(d).

The weights decrease with increasing distances. Lower is h, sharper is the decreasing function.
Weights are set to 0 for outliers (extreme distances).
"""
function wdist(d; h = 2, cri = 4, squared = false)
    w = copy(d)
    wdist!(w; h = h, cri = cri, squared = squared)
    w
end

function wdist!(d; h = 2, cri = 4, squared = false)
    # d, out : (n,)
    squared ? d = d.^2 : nothing
    zmed =  Statistics.median(d)
    zmad = Jchemo.mad(d)
    cutoff = zmed + cri * zmad
    d .= map(x -> ifelse(x <= cutoff, exp(-x / (h * zmad)), zero(eltype(d))), d)
    d .= d / maximum(d)
    d[isnan.(d)] .= 1
    return
end


text:
"""
	RegularNetwork <: AbstractNetwork

A regular (ring) network.

# Fields
- N              :: Integer -> Number of nodes
- k              :: Integer -> Node connectivity
"""
mutable struct RegularNetwork <: AbstractNetwork
	N :: Integer
	k :: Integer

	_adjMat :: Union{AbstractMatrix, Nothing}
	_props :: NetworkProperties
end

"""
	RegularNetwork(N::Integer, k::Integer; directed::Bool=false)

Create a regular network with `N` nodes and node connectivity `k`.
"""
function RegularNetwork(N::Integer, k::Integer; directed::Bool=false)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))

	(k % 2 == 1) && (k -= 1)
	(k < 0) && (k = 0)
	(k >= N-1) && (k = N-1)

	RegularNetwork(N, k, nothing, NetworkProperties(N, directed))
end

function show(network::RegularNetwork)
	print(isdirected(network) ? "Directed " : "Undirected ")
	println("Regular network")
	println("- N = $(network.N)")
	println("- k = $(network.k)")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 973169b2-5cb5-11eb-2c27-8115e5bfec1f
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
end

# ╔═╡ 662d102a-5cb6-11eb-1f92-1b7f561b048a
md"""
## The Baird Count Environment
"""

# ╔═╡ 11ba9900-5cb5-11eb-0d89-0f51f890b091
begin
	const DASH_SOLID = (:dashed, :solid)

	Base.@kwdef mutable struct BairdCounterEnv <: AbstractEnv
		current::Int = rand(1:7)
	end

	RLBase.state_space(env::BairdCounterEnv) = Base.OneTo(7)
	RLBase.action_space(env::BairdCounterEnv) = Base.OneTo(length(DASH_SOLID))

	function (env::BairdCounterEnv)(a)
		if DASH_SOLID[a] == :dashed
			env.current = rand(1:6)
		else
			env.current = 7
		end
		nothing
	end

	RLBase.reward(env::BairdCounterEnv) = 0.
	RLBase.is_terminated(env::BairdCounterEnv) = false
	RLBase.state(env::BairdCounterEnv) = env.current
	RLBase.reset!(env::BairdCounterEnv) = env.current = rand(1:6)
end

# ╔═╡ 915ff55a-5cb6-11eb-09bc-2f489a33a737
md"""
## Off Policy
"""

# ╔═╡ b5711c62-5cb6-11eb-2b70-91051ff959d3
# Base.@kwdef struct OffPolicy{P,B} <: AbstractPolicy
#     π_target::P
#     π_behavior::B
# end

# ╔═╡ ba6d1efc-5cb6-11eb-28e5-f73f80f44ad7
# (π::OffPolicy)(env) = π.π_behavior(env)

# ╔═╡ 447be6a0-5cbf-11eb-2daf-9b1b54e403f1
begin

	# const VectorWSARTTrajectory = Trajectory{<:NamedTuple{(:weight, SART...)}}

	# function VectorWSARTTrajectory(;weight=Float64, state=Int, action=Int, reward=Float32, terminal=Bool)
	# 	VectorTrajectory(;weight=Float64, state=state, action=action, reward=reward, terminal=terminal)
	# end

	# function RLBase.update!(
	# 	p::OffPolicy,
	# 	t::VectorWSARTTrajectory,
	# 	e::AbstractEnv,
	# 	s::AbstractStage
	# )
	# 	update!(p.π_target, t, e, s)
	# end

	# function RLBase.update!(
	# 	t::VectorWSARTTrajectory,
	# 	p::OffPolicy,
	# 	env::AbstractEnv,
	# 	s::PreActStage,
	# 	a
	# )
	# 	push!(t[:state], state(env))
	# 	push!(t[:action], a)

	# 	w = prob(p.π_target, s, a) / prob(p.π_behavior, s, a)
	# 	push!(t[:weight], w)
	# end

	# function RLBase.update!(
	# 	t::VectorWSARTTrajectory,
	# 	p::OffPolicy{<:QBasedPolicy{<:TDLearner}},
	# 	env::AbstractEnv,
	# 	s::PreEpisodeStage,
	# )
	# 	empty!(t)
	# end

	# function RLBase.update!(
	# 	t::VectorWSARTTrajectory,
	# 	p::OffPolicy{<:QBasedPolicy{<:TDLearner}},
	# 	env::AbstractEnv,
	# 	s::PostEpisodeStage,
	# )
	# 	action = rand(action_space(env))

	# 	push!(trajectory[:state], state(env))
	# 	push!(trajectory[:action], action)
	# 	push!(t[:weight], 1.0)
	# end
end

# ╔═╡ f01b15c2-5cb6-11eb-1d41-21f632c01a2f
md"""
## Figure 11.2
"""

# ╔═╡ 45c5fdfe-5cb8-11eb-115e-93fcb138261d
world = BairdCounterEnv()

# ╔═╡ c9f52da8-5cb5-11eb-35eb-ad54e81df746
begin
	Base.@kwdef struct RecordWeights <: AbstractHook
		weights::Vector{Vector{Float64}}=[]
	end

	(h::RecordWeights)(::PostActStage, agent, env) = push!(
		h.weights,
		agent.policy.π_target.learner.approximator.weights |> deepcopy
	)
end

# ╔═╡ 30e3d956-5cb8-11eb-05be-ab2d4a182b4e
NW = 8

# ╔═╡ 6396bcb2-5cb8-11eb-2c4a-dd32adf35da6
INIT_WEIGHT = ones(8)

# ╔═╡ 7f294128-5cb8-11eb-3b3f-8b74678e407f
INIT_WEIGHT[7] = 10

# ╔═╡ 5ed1b616-5cb6-11eb-3678-f9795f222f8a
STATE_MAPPING = zeros(NW, length(state_space(world)))

# ╔═╡ 9df6e892-5cb8-11eb-28d3-098d7aff3db0
begin
	for i in 1:6
		STATE_MAPPING[i, i] = 2
		STATE_MAPPING[8, i] = 1
	end
	STATE_MAPPING[7, 7] = 1
	STATE_MAPPING[8, 7] = 2
end

# ╔═╡ bc09d846-5cb8-11eb-0638-7d7a435e5fc4
STATE_MAPPING

# ╔═╡ bff74902-5cb8-11eb-2736-638fcd9a1b5b
π_b = x -> rand() < 6/7 ? 1 : 2

# ╔═╡ 181f9792-5cb9-11eb-2575-6995ea64416d
π_t = VBasedPolicy(
	learner=TDLearner(
		approximator=RLZoo.LinearApproximator(INIT_WEIGHT, Descent(0.01)),
		γ=0.99,
		n=0,
		method=:SRS
	),
	mapping = (env, V) -> 2
)

# ╔═╡ 7ea3658e-5cb9-11eb-297c-1561db9958eb
prob_b = [6/7, 1/7]

# ╔═╡ 9ba05ebc-5cb9-11eb-1a41-cbf495eedbe5
prob_t = [0., 1.]

# ╔═╡ aa855324-5cb9-11eb-1d44-8b3584e93983
md"""
Well, I must admit it is a little tricky here.
"""

# ╔═╡ a1b4b438-5cb9-11eb-0017-4594a2929891
RLBase.prob(::typeof(π_b), s, a::Integer) = prob_b[a]

# ╔═╡ bc514e16-5cb9-11eb-0452-87f09f1bca1a
RLBase.prob(::typeof(π_t), s, a::Integer) = prob_t[a]

# ╔═╡ c165ae40-5cb9-11eb-197c-e10a30939e4e
agent = Agent(
    policy=OffPolicy(
        π_target=π_t,
        π_behavior=π_b
    ),
    trajectory=VectorWSARTTrajectory(state=Any)
)

# ╔═╡ 76cc56d0-5cba-11eb-34b2-771d4a8cb6c6
new_env = StateTransformedEnv(
    BairdCounterEnv(),
    ;state_mapping=s -> @view STATE_MAPPING[:, s]
)

# ╔═╡ a7a82a90-5cba-11eb-0743-2364ec991a79
hook = RecordWeights()

# ╔═╡ b0f5f690-5cba-11eb-1480-b300a5d8aa7e
run(agent, new_env, StopAfterStep(1000),hook)

# ╔═╡ c0f640ea-5cba-11eb-13da-27969e97d064
begin
	p = plot(legend=:topleft)
	for i in 1:length(INIT_WEIGHT)
		plot!(p, [w[i] for w in hook.weights])
	end
	p
end

# ╔═╡ Cell order:
# ╠═973169b2-5cb5-11eb-2c27-8115e5bfec1f
# ╠═662d102a-5cb6-11eb-1f92-1b7f561b048a
# ╠═11ba9900-5cb5-11eb-0d89-0f51f890b091
# ╟─915ff55a-5cb6-11eb-09bc-2f489a33a737
# ╠═b5711c62-5cb6-11eb-2b70-91051ff959d3
# ╠═ba6d1efc-5cb6-11eb-28e5-f73f80f44ad7
# ╠═447be6a0-5cbf-11eb-2daf-9b1b54e403f1
# ╟─f01b15c2-5cb6-11eb-1d41-21f632c01a2f
# ╠═45c5fdfe-5cb8-11eb-115e-93fcb138261d
# ╠═c9f52da8-5cb5-11eb-35eb-ad54e81df746
# ╠═30e3d956-5cb8-11eb-05be-ab2d4a182b4e
# ╠═6396bcb2-5cb8-11eb-2c4a-dd32adf35da6
# ╠═7f294128-5cb8-11eb-3b3f-8b74678e407f
# ╠═5ed1b616-5cb6-11eb-3678-f9795f222f8a
# ╠═9df6e892-5cb8-11eb-28d3-098d7aff3db0
# ╠═bc09d846-5cb8-11eb-0638-7d7a435e5fc4
# ╠═bff74902-5cb8-11eb-2736-638fcd9a1b5b
# ╠═181f9792-5cb9-11eb-2575-6995ea64416d
# ╠═7ea3658e-5cb9-11eb-297c-1561db9958eb
# ╠═9ba05ebc-5cb9-11eb-1a41-cbf495eedbe5
# ╟─aa855324-5cb9-11eb-1d44-8b3584e93983
# ╠═a1b4b438-5cb9-11eb-0017-4594a2929891
# ╠═bc514e16-5cb9-11eb-0452-87f09f1bca1a
# ╠═c165ae40-5cb9-11eb-197c-e10a30939e4e
# ╠═76cc56d0-5cba-11eb-34b2-771d4a8cb6c6
# ╠═a7a82a90-5cba-11eb-0743-2364ec991a79
# ╠═b0f5f690-5cba-11eb-1480-b300a5d8aa7e
# ╠═c0f640ea-5cba-11eb-13da-27969e97d064


text:
#=
Find toplevel items (bind / call)

Downstreams: modules.jl, outline.jl, goto.jl
=#

abstract type ToplevelItem end

struct ToplevelBinding <: ToplevelItem
    expr::EXPR
    bind::Binding
    lines::UnitRange{Int}
end

struct ToplevelCall <: ToplevelItem
    expr::EXPR
    lines::UnitRange{Int}
    verbatim::String
end

struct ToplevelMacroCall <: ToplevelItem
    expr::EXPR
    lines::UnitRange{Int}
    verbatim::String
end

struct ToplevelModuleUsage <: ToplevelItem
    expr::EXPR
    lines::UnitRange{Int}
end

"""
    toplevelitems(text::String; kwargs...)::Vector{ToplevelItem}
    toplevelitems(text::String, expr::EXPR; kwargs...)::Vector{ToplevelItem}

Finds and returns toplevel "item"s (call and binding) in `text`.

keyword arguments:
- `mod::Union{Nothing, String}`: if not `nothing` don't return items within modules
    other than `mod`, otherwise enter into every module.
- `inmod::Bool`: if `true`, don't include toplevel items until it enters into `mod`.
"""
toplevelitems(text::String; kwargs...)::Vector{ToplevelItem} =
    toplevelitems(text, CSTParser.parse(text, true); kwargs...)
function toplevelitems(text::String, expr::EXPR; kwargs...)::Vector{ToplevelItem}
    traverse_expr!(expr)
    return _toplevelitems(text, expr; kwargs...)
end
toplevelitems(text::String, expr::Nothing; kwargs...)::Vector{ToplevelItem} = ToplevelItem[]

function _toplevelitems(
    text::String, expr::EXPR,
    items::Vector{ToplevelItem} = ToplevelItem[], line::Int = 1, pos::Int = 1;
    mod::Union{Nothing, String} = nothing, inmod::Bool = false,
)
    # add items if `mod` isn't specified or in a target modle
    if mod === nothing || inmod
        # binding
        bind = bindingof(expr)
        if bind !== nothing
            lines = line:line+countlines(expr, text, pos, false)
            push!(items, ToplevelBinding(expr, bind, lines))
        end

        lines = line:line+countlines(expr, text, pos, false)

        # destructure multiple returns
        if ismultiplereturn(expr)
            for arg in expr
                (bind = bindingof(arg)) === nothing && continue
                push!(items, ToplevelBinding(arg, bind, lines))
            end
        elseif expr.parent === nothing || !CSTParser.is_assignment(expr.parent)
            # toplevel call
            if iscallexpr(expr)
                push!(items, ToplevelCall(expr, lines, str_value_verbatim(expr, text, pos)))
            # toplevel macro call
            elseif ismacrocall(expr)
                push!(items, ToplevelMacroCall(expr, lines, str_value_verbatim(expr, text, pos)))
            # module usages
            elseif ismoduleusage(expr)
                push!(items, ToplevelModuleUsage(expr, lines))
            end
        end
    end

    if shouldenter(expr, mod)
        # look for more toplevel items in expr:
        if CSTParser.defines_module(expr) && shouldentermodule(expr, mod)
            inmod = true
        end
        for arg in expr
            _toplevelitems(text, arg, items, line, pos; mod = mod, inmod = inmod)
            line += countlines(arg, text, pos)
            pos += arg.fullspan
        end
    end
    return items
end

function shouldenter(expr::EXPR, mod::Union{Nothing, String})
    typof(expr) !== CSTParser.Call &&
    !(hasscope(expr) &&
    !(
        typof(expr) === CSTParser.FileH ||
        (CSTParser.defines_module(expr) && shouldentermodule(expr, mod)) ||
        isdoc(expr)
    ))
end

shouldentermodule(expr::EXPR, mod::Nothing) = true
shouldentermodule(expr::EXPR, mod::String) =
    (bind = bindingof(expr)) === nothing ? false : bind.name == mod


text:
"""
When sample is called with the `save_state` kwarg set to `true`, the chain can be used to plot the tempering swaps that occurred during sampling
"""
function plot_swaps(chain)
    plot(chain.info.samplerstate.Δ_index_history)
end

function plot_ladders(chain)
    plot(chain.info.samplerstate.Δ_history)
end


text:
using ITensors

function main()
  Ny = 6
  Nx = 12

  N = Nx*Ny

  sites = siteinds("S=1/2",N)

  lattice = square_lattice(Nx,Ny,yperiodic=false)
  #lattice = triangular_lattice(Nx,Ny,yperiodic=false)

  ampo = AutoMPO()
  for b in lattice
    add!(ampo,0.5,"S+",b.s1,"S-",b.s2)
    add!(ampo,0.5,"S-",b.s1,"S+",b.s2)
    add!(ampo,    "Sz",b.s1,"Sz",b.s2)
  end
  H = MPO(ampo, sites)

  state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  psi0 = productMPS(sites, state)

  sweeps = Sweeps(10)
  maxdim!(sweeps,10,20,100,100,200,400,800)
  cutoff!(sweeps,1E-8)
  @show sweeps

  energy,psi = dmrg(H,psi0,sweeps)

end
main()


text:
include("../src/LarSurf.jl")

import SparseArrays.spzeros
import SparseArrays.dropzeros!
using Plasm
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

segmentation = zeros(Int8, 5, 6, 7)
segmentation[2:5,2:5,2:6] .= 1
# Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))

filteredFV, Flin, (V, tmodel) = LarSurf.get_surface_grid(segmentation)
(VV, EV, FV, CV) = tmodel
Plasm.View((V,[VV, EV, filteredFV]))


text:
using PortHamiltonian
using PyPlot
using PyCall
@pyimport matplotlib.patches as patch

cfig = figure()
ax = cfig[:add_subplot](1,1,1)

ax[:set_aspect]("equal")
Nint = 5
# weak form
Nelem = 20
p = finelem(Nelem,Nint,0,1)

for i = 1:Nelem
	if (i%2 == 1)
	c = patch.Rectangle(
		(-0.5+(i-1)*(Nint*2+1), -0.5+(i-1)*(Nint*2+1)), Nint*2+1, Nint*2+1,
		alpha=0.3,fc="blue")
	else
	c = patch.Rectangle(
		(-0.5+(i-1)*(Nint*2+1), -0.5+(i-1)*(Nint*2+1)), Nint*2+1, Nint*2+1,
		alpha=0.3,fc="red")
	end
	ax[:add_artist](c)
end
spy(p.J)

pM = PortHamiltonian.new_spectral_element_phs(Nelem,Nint+1,0,1)

figure();
spy(pM.J)
freqp = frequencies(p)
freqpM = frequencies(pM)
[freqp[1:10] freqp[1:10]-[1:2:20]*pi/2 freqpM[1:10] freqpM[1:10]-[1:2:20]*pi/2]

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
import MechanicalSketch
import MechanicalSketch: empty_figure, O, WI, HE, finish
import MechanicalSketch: PALETTE
import MechanicalSketch: circle, scale_to_pt, scale_pt_to_unit
import MechanicalSketch: @import_expand, @layer, Point
import MechanicalSketch: ∙, °
import MechanicalSketch: placeimage, readsvg, readpng, color_with_lumin
import MechanicalSketch: Luxor, Length, set_scale_sketch, text, SVGimage
import Luxor:            Cairo.CairoSurfaceBase, translate, scale

let

if !@isdefined m²
    @import_expand(~m, s, °)
end
include("test_functions_35.jl")
empty_figure(filename = joinpath(@__DIR__, "test_35.png"),
    backgroundcolor = color_with_lumin(PALETTE[8], 30));

totheight = 10.0m
set_scale_sketch(totheight, HE)
totwidth = totheight * WI / HE
Δx = totwidth / 4
Δy = totheight / 4
centers = reshape([Point(x, y) for x in Δx * [-1, 0, 1], y in Δy * [1, 0, -1]], :, 1)
θs = [180, 225, 270]°

for (i, cent) in enumerate(centers)
    circle(cent; d = 1m)
    circle.(cent; r = 2m)
    text(string(i), cent + (1, 1)m, )
end
svfnam = joinpath("..", "resource", "dn69ygSyBh.svg")
isfile(svfnam)
sv = readsvg(svfnam)
place_image_35(centers[1], sv, centered = true)
place_image_35(centers[2], sv, centered = false)
place_image_35(centers[3], sv, centered = true, height = 1m)
place_image_35(centers[4], sv, centered = false, height = 1m)
place_image_35(centers[5], sv, centered = true, width = 1m)
place_image_35(centers[6], sv, centered = false, width = 1m)
place_image_35(centers[7] + 2m .* (cos(θs[1]), sin(θs[1])), sv, centered = true, width = 1m)
place_image_35(centers[8] + 2m .* (cos(θs[2]), sin(θs[2])), sv, centered = true, width = 1m)
place_image_35(centers[9] + 2m .* (cos(θs[2]), sin(θs[2])), sv, centered = false, width = 1m)

end # let

finish()


text:
include("auto_encoder_task.jl")

# ::: NETWORK
networks = Dict(
    "full_sigmoid" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, diminish, Δdiminish),
        DenseLayer(500, 300),
        ActivationLayer(300, diminish, Δdiminish),
        DenseLayer(300, 150),
        ActivationLayer(150, diminish, Δdiminish),
        DenseLayer(150, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
        ],[
        DenseLayer(LATENT_SIZE, 150),
        ActivationLayer(150, diminish, Δdiminish),
        DenseLayer(150, 300),
        ActivationLayer(300, diminish, Δdiminish),
        DenseLayer(300, 500),
        ActivationLayer(500, diminish, Δdiminish),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
   "partial_sigmoid1" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 300),
        ActivationLayer(300, diminish, Δdiminish),
        DenseLayer(300, 150),
        ActivationLayer(150, diminish, Δdiminish),
        DenseLayer(150, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
        ],[
        DenseLayer(LATENT_SIZE, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, 300),
        ActivationLayer(300, diminish, Δdiminish),
        DenseLayer(300, 500),
        ActivationLayer(500, diminish, Δdiminish),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
    "partial_sigmoid2" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 150),
        ActivationLayer(150, diminish, Δdiminish),
        DenseLayer(150, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
       ],[
        DenseLayer(LATENT_SIZE, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 500),
        ActivationLayer(500, diminish, Δdiminish),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
    "partial_sigmoid3" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
       ],[
        DenseLayer(LATENT_SIZE, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]))

# ::: TRAINING
runtime_plotting = false
plot_intervall = 20
begin
    full_sigmoid_tm = [NaN for _ in 1:TRAIN_EPISODES];
    partial_sigmoid_tm1 = [NaN for _ in 1:TRAIN_EPISODES];
    partial_sigmoid_tm2 = [NaN for _ in 1:TRAIN_EPISODES];
    partial_sigmoid_tm3 = [NaN for _ in 1:TRAIN_EPISODES];

    println("---"); @time for e in 1:TRAIN_EPISODES
        t = time()

        z1 = forward!(networks["full_sigmoid"][1], train_data_x[e])
        z2 = forward!(networks["partial_sigmoid1"][1], train_data_x[e])
        z3 = forward!(networks["partial_sigmoid2"][1], train_data_x[e])
        z4 = forward!(networks["partial_sigmoid3"][1], train_data_x[e])

        y1 = forward!(networks["full_sigmoid"][2], z1)
        y2 = forward!(networks["partial_sigmoid1"][2], z2)
        y3 = forward!(networks["partial_sigmoid2"][2], z3)
        y4 = forward!(networks["partial_sigmoid3"][2], z4)

        dloss1 = -Δmse(train_data_x[e], y1)
        dloss2 = -Δmse(train_data_x[e], y2)
        dloss3 = -Δmse(train_data_x[e], y3)
        dloss4 = -Δmse(train_data_x[e], y4)

        dz1 = backward!(reverse(networks["full_sigmoid"][2]), dloss1, opt)
        _ = backward!(reverse(networks["full_sigmoid"][1]), dz1, opt)
        dz2 = backward!(reverse(networks["partial_sigmoid1"][2]), dloss2, opt)
        _ = backward!(reverse(networks["partial_sigmoid1"][1]), dz2, opt)
        dz3 = backward!(reverse(networks["partial_sigmoid2"][2]), dloss3, opt)
        _ = backward!(reverse(networks["partial_sigmoid2"][1]), dz3, opt)
        dz4 = backward!(reverse(networks["partial_sigmoid3"][2]), dloss4, opt)
        _ = backward!(reverse(networks["partial_sigmoid3"][1]), dz4, opt)

        full_sigmoid_tm[e] = mse(train_data_x[e], y1)
        partial_sigmoid_tm1[e] = mse(train_data_x[e], y2)
        partial_sigmoid_tm2[e] = mse(train_data_x[e], y3)
        partial_sigmoid_tm3[e] = mse(train_data_x[e], y4)

        # if runtime_plotting && e % plot_intervall == 0
        #     p = plot(FC_train_metric)
        #     display(p)
        # end

        # println("episode = $e : loss = $loss : dt = $(time() - t)")
    end
end

begin
    gr = 30
    plot([sum(full_sigmoid_tm[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC (full diminish): MSE (lr=$LEARNING_RATE, batch size=$BATCH_SIZE)")
    plot!([sum(partial_sigmoid_tm1[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC (Relu [1] + diminish [3]): MSE (lr=$LEARNING_RATE, batch size=$BATCH_SIZE)")
    plot!([sum(partial_sigmoid_tm2[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC (Relu [2] + diminish [2]): MSE (lr=$LEARNING_RATE, batch size=$BATCH_SIZE)")
    plot!([sum(partial_sigmoid_tm3[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC (Relu [3] + diminish [1]): MSE (lr=$LEARNING_RATE, batch size=$BATCH_SIZE)")
end

generate_images(10, "/Users/BBM2H16AHO/Documents/programming/DNN-Framework/auto_encoder/images/relu_sigmoid_test")


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

using Markdown
using InteractiveUtils

# ╔═╡ 1c7e914c-03c9-11eb-11cf-0d9346305a53
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using NLopt
end

# ╔═╡ 1866ec44-03c9-11eb-0a21-416587e4e0e0
# Solve the problem with NLopt
let

	"""
	Function to be minimized (with its gradient)
	"""
	function myfun(x::Vector, grad::Vector)
		if length(grad) > 0
			grad[1] = -1
			grad[2] = -1
		end
		return -x[1] - x[2]
	end

	"""
	Constraints in form c_i(x)  <= 0 (with gradient)
	"""
	function c1(x::Vector, grad::Vector)
		if length(grad) > 0
			grad[1] = 2*x[1]
			grad[2] = -1
		end
		x[1]^2 - x[2]
	end

	function c2(x::Vector, grad::Vector)
		if length(grad) > 0
			grad[1] = 2*x[1]
			grad[2] = 2*x[2]
		end
		x[1]^2 + x[2]^2 - 1
	end

	# optimisation setup
	opt = Opt(:LD_SLSQP, 2)
	opt.lower_bounds = [-Inf, 0.] # only x_2 > 0
	opt.xtol_rel = 1e-4
	opt.min_objective = myfun
	inequality_constraint!(opt, c1, 1e-8)
	inequality_constraint!(opt, c2, 1e-8)

	# solving it.
	(minf,minx,ret) = NLopt.optimize(opt, [0, 0])
	numevals = opt.numevals # the number of function evaluations
	println("got $minf at $minx after $numevals iterations (returned $ret)")
end

# ╔═╡ 2121e5be-03c9-11eb-2b8b-a10cb146216d
ret

# ╔═╡ Cell order:
# ╠═1c7e914c-03c9-11eb-11cf-0d9346305a53
# ╠═1866ec44-03c9-11eb-0a21-416587e4e0e0
# ╠═2121e5be-03c9-11eb-2b8b-a10cb146216d

max_stars_repo_path: Exercises/PS04 - Optimisation techniques II.jl


text:
"""
    A = JotOpBlend(T, nsamples, shottimes[, nsamples_blended])

where `A` is a shot mixing operation, so that if `d=Am`, then `m::Array{T,2}` is a
common receiver gather with each trace corresponding to a different shot.  The
excitation time of each are the entries of `shottimes::Array{Int,1}` in units of
time samples.  The size of `m` is `(nsamples,length(shottimes))`.  The size of
`d` is either `nsamples_blended` (if provided) or `maximum(shot_times)+nsamples-1`.

# Examples

## blend two receiver gather traces:
```julia
A = JotOpBlend(Float64, 128, [1,64])
m = rand(domain(A))
d = A*m # receiver gather with blended shots
```
"""
function JopBlend(T, nsamples, shottimes, nsamples_blended=0)
    if nsamples_blended == 0
        nsamples_blended = maximum(shottimes) + nsamples - 1
    end

    dom = JetSpace(T, nsamples, length(shottimes))
    rng = JetSpace(T, nsamples_blended)

    JopLn(df! = JopBlend_df!, df′! = JopBlend_df′!, dom = dom, rng = rng, s = (st=shottimes,))
end
export JopBlend

function JopBlend_df!(d::AbstractArray{T,1}, m::AbstractArray{T,2}; st, kwargs...) where {T}
    d .= 0
    nsamples, nshots = size(m)
    nsamples_blended = length(d)
    for ishot = 1:nshots, isample = 1:nsamples
        j = st[ishot] + isample - 1
        if j ∈ 1:nsamples_blended
            d[j] += m[isample,ishot]
        end
    end
    d
end

function JopBlend_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,1}; st, kwargs...) where {T}
    m .= 0
    nsamples, nshots = size(m)
    nsamples_blended = length(d)
    for ishot = 1:nshots, isample = 1:nsamples
        j = st[ishot] + isample - 1
        if j ∈ 1:nsamples_blended
            m[isample,ishot] = d[j]
        end
    end
    m
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Glib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Glib")
JLLWrappers.@generate_main_file("Glib", UUID("7746bdde-850d-59dc-9ae8-88ece973131d"))
end  # module Glib_jll


text:
module HelperFunctions

using XLSX

export benchmarkToExcel

# this array of strings is used to navigate in the excel file using XLSX
const excelXLocation = ["D","E","F","G","H","I","J","K","L","M","N","O","P","Q","R","S","T","U","V","W","X","Y","Z","AA","AB","AC","AD","AE","AF","AG","AH","AI","AJ","AK","AL","AM","AN","AO","AP","AQ","AR","AS","AT","AU","AV","AW","AX","AY","AZ"];
#file to handle benchmarkin process for arbitrary model

# helper functions to transfer memory from bytes to better readable form as in the @benachmark output
# e.g. "11.49 GiB" = helperMemoryEst(12345678910)
function helperMemoryEst(memory::Int)

	strDim = ["B", "KiB", "MeB", "GiB", "TeB", "PeB", "ExB", "ZeB"]
	countDim = 1

	while (memory >= 1024 && countDim < 9)
        memory = memory/1024
		countDim = countDim + 1
	end

	#strMemory = string(memory)
	#len = length(strMemory)
	#strMemory = string(strMemory[1:min(len,5)], " ",strDim[countDim])

	return memory, strDim[countDim]
end

# helper functions to transfer time from nanoseconds to better readable form as in the benchmark output
# e.g. x = helper
function helperTimes(time::Float64)

	strDim = ["ns", "µs", "ms", "s", "m", "h"]
	countDim = 1

	while (time >= 1000 && countDim < 4)
        time = time/1000
		countDim = countDim + 1
	end

	while (time > 60 && countDim >= 4 && countDim < 7)
        time = time/60
		countDim = countDim + 1
	end

	#strTime = string(time)
	#len = length(strTime)
	#strTime = string(strTime[1:min(len,5)], " ",strDim[countDim])

	return time, strDim[countDim]
end

#function that writes the benchmark into an excel file at a specific location
function benchmarkToExcel(excelFile::String, sheetNbr::Int, benchmark::Any, location::String, objValue::Float64)
	# sheet			excel sheet
	# benchmark		benchmark data from @benchmark
	# location 		desired location in excel file eg "A2"

	#transfer benchmark data to array
	numberSamples = length(benchmark.times)
	#times in benchmark.times are ordered ascending
	times = benchmark.times;
	medianTime = times[Int(round(length(times)/2))]
	medTime, dimTime = helperTimes(medianTime)
	memory, dimMemory = helperMemoryEst(benchmark.memory)

	#the zeros are added cause otherwise benchmarks are saved as an Array{Float64,1} and XLSX save this format only as a row vector but we need it as a column vector
	#is there a better workaround? idea is that the zeros get overwritten in every step when this function gets repeated with a shifted location so they are not to distrubing...
	benchmarkvalues = [[memory dimMemory]; [medTime dimTime]; [ "" ""]; [numberSamples ""]; [ ((objValue == 0.0) ? "" : objValue ) ""]]
	XLSX.openxlsx(excelFile, mode="rw") do xf

		#in the excel xf we write on the first sheet
		sheet = xf[sheetNbr]
		#save the array to the excel sheet
		sheet[location] = benchmarkvalues

	end
end

end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
# QCP example
#    minimize x^2
#
#    s.t.  x >= 0.1
#          x ∈ {0, 1}
#
#    solution: (0.1) objv = 0.01

@testset "MIQCP" begin
    @testset "MIQCP 01" begin
        m = CPLEX.CplexMathProgModel(CPX_PARAM_SCRIND=0)
        CPLEX.loadproblem!(m, Matrix{Float64}(undef, 0,1), [0.1], [Inf], [0],
                           Float64[], Float64[], :Min)
        @test CPLEX.get_prob_type(m.inner) == :LP
        MathProgBase.setquadobj!(m, reshape([2],(1,1)))
        @test CPLEX.get_prob_type(m.inner) == :QP
        MathProgBase.optimize!(m)
        @test isapprox(CPLEX.getobjval(m), 0.5*0.1*2*0.1)
    end

    @testset "MIQCP 02" begin
        m = CPLEX.CplexMathProgModel(CPX_PARAM_SCRIND=0)
        CPLEX.loadproblem!(m, Matrix{Float64}(undef, 0,1), [0.1], [Inf], [0], Float64[], Float64[], :Min)
        @test CPLEX.get_prob_type(m.inner) == :LP
        CPLEX.setvartype!(m, [:Bin])
        @test CPLEX.get_prob_type(m.inner) == :MILP
        MathProgBase.setquadobj!(m, reshape([2],(1,1)))
        @test CPLEX.get_prob_type(m.inner) == :MIQP
        MathProgBase.optimize!(m)
        @test isapprox(CPLEX.getobjval(m), 0.5*1*2*1)
    end

    @testset "MIQCP 03" begin
        m = CPLEX.CplexMathProgModel(CPX_PARAM_SCRIND=0)
        CPLEX.loadproblem!(m, Matrix{Float64}(undef, 0,1), [0.1], [Inf], [0], Float64[], Float64[], :Min)
        @test CPLEX.get_prob_type(m.inner) == :LP
        MathProgBase.setquadobj!(m, reshape([2],(1,1)))
        @test CPLEX.get_prob_type(m.inner) == :QP
        CPLEX.setvartype!(m, [:Bin])
        @test CPLEX.get_prob_type(m.inner) == :MIQP
        MathProgBase.optimize!(m)
        @test isapprox(CPLEX.getobjval(m), 0.5*1*2*1)
    end
end


text:
function findperiod(cf::PartialQuotients)
  return 0
end

"""
sqrt(n) = [a0; a1 a2 ... a2 a1 2*a0 a1 a2...] with period "a1 a2 ... a2 a1 2*a0".
"""
function findperiod(cf::PartialQuotients{QuadraticIrrational})
  frac = cf.partialquotients
  len = length(frac)
  (len < 2) && return 0
  (frac[end] != 2 * frac[1]) && return 0
  stop = false
  period_len = (len % 2 == 0) ? 1 : 2
  i = 1
  while !stop
    if 1 + i < len - i
      if frac[1 + i] != frac[len - i]
        return 0
      end
      period_len += 2
    else
      stop = true
    end
    i += 1
  end
  return period_len
end


text:
# Julia wrapper for header: /usr/include/postgresql/libpq-fe.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

module libpq_interface
           # types
    export PGconn,
           PGresult,
           Oid,
           # enums
           ConnStatusType,
           CONNECTION_OK,
           CONNECTION_BAD,
           CONNECTION_STARTED,
           CONNECTION_MADE,
           CONNECTION_AWAITING_RESPONSE,
           CONNECTION_AUTH_OK,
           CONNECTION_SETENV,
           CONNECTION_SSL_STARTUP,
           CONNECTION_NEEDED,
           ExecStatusType,
           PGRES_EMPTY_QUERY,
           PGRES_COMMAND_OK,
           PGRES_TUPLES_OK,
           PGRES_COPY_OUT,
           PGRES_COPY_IN,
           PGRES_BAD_RESPONSE,
           PGRES_NONFATAL_ERROR,
           PGRES_FATAL_ERROR,
           PGRES_COPY_BOTH,
           PGRES_SINGLE_TUPLE,
           PGTransactionStatusType,
           PQTRANS_IDLE,
           PQTRANS_ACTIVE,
           PQTRANS_INTRANS,
           PQTRANS_INERROR,
           PQTRANS_UNKNOWN,
           PGF_TEXT,
           PGF_BINARY,
           # connection functions
           PQconnectdb,
           PQconnectdbParams,
           PQsetdbLogin,
           PQfinish,
           PQreset,
           # connection status functions
           PQdb,
           PQuser,
           PQpass,
           PQhost,
           PQport,
           PQsocket,
           PQstatus,
           PQtransactionStatus,
           PQprotocolVersion,
           PQerrorMessage,
           # command execution functions
           PQexec,
           PQexecParams,
           PQprepare,
           PQexecPrepared,
           PQdescribePrepared,
           PQgetvalue,
           PQgetisnull,
           PQnparams,
           PQparamtype,
           PQntuples,
           PQnfields,
           PQftype,
           PQfname,
           PQgetlength,
           PQgetResult,
           PQresStatus,
           PQresultStatus,
           PQresultErrorMessage,
           PQclear,
           PQputCopyData,
           PQputCopyEnd,
           # misc
           PQescapeLiteral,
           PQescapeIdentifier,
           PQfreemem

    include("libpq_common.jl")

    @c Cint _IO_getc (Ptr{_IO_FILE},) libpq
    @c Cint _IO_putc (Cint, Ptr{_IO_FILE}) libpq
    @c Cint _IO_feof (Ptr{_IO_FILE},) libpq
    @c Cint _IO_ferror (Ptr{_IO_FILE},) libpq
    @c Cint _IO_peekc_locked (Ptr{_IO_FILE},) libpq
    @c Void _IO_flockfile (Ptr{_IO_FILE},) libpq
    @c Void _IO_funlockfile (Ptr{_IO_FILE},) libpq
    @c Cint _IO_ftrylockfile (Ptr{_IO_FILE},) libpq
    @c Cint _IO_vfscanf (Ptr{_IO_FILE}, Ptr{UInt8}, Cint, Ptr{Cint}) libpq
    @c Cint _IO_vfprintf (Ptr{_IO_FILE}, Ptr{UInt8}, Cint) libpq
    # @c __ssize_t _IO_padn (Ptr{_IO_FILE}, Cint, __ssize_t) libpq
    @c Cint _IO_sgetn (Ptr{_IO_FILE}, Ptr{Void}, Cint) libpq
    @c __off64_t _IO_seekoff (Ptr{_IO_FILE}, __off64_t, Cint, Cint) libpq
    @c __off64_t _IO_seekpos (Ptr{_IO_FILE}, __off64_t, Cint) libpq
    @c Void _IO_free_backup_area (Ptr{_IO_FILE},) libpq
    @c Cint remove (Ptr{UInt8},) libpq
    @c Cint rename (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint renameat (Cint, Ptr{UInt8}, Cint, Ptr{UInt8}) libpq
    @c Ptr{FILE} tmpfile () libpq
    @c Ptr{UInt8} tmpnam (Ptr{UInt8},) libpq
    @c Ptr{UInt8} tmpnam_r (Ptr{UInt8},) libpq
    @c Ptr{UInt8} tempnam (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint fclose (Ptr{FILE},) libpq
    @c Cint fflush (Ptr{FILE},) libpq
    @c Cint fflush_unlocked (Ptr{FILE},) libpq
    @c Ptr{FILE} fopen (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Ptr{FILE} freopen (Ptr{UInt8}, Ptr{UInt8}, Ptr{FILE}) libpq
    @c Ptr{FILE} fdopen (Cint, Ptr{UInt8}) libpq
    @c Ptr{FILE} fmemopen (Ptr{Void}, Cint, Ptr{UInt8}) libpq
    @c Ptr{FILE} open_memstream (Ptr{Ptr{UInt8}}, Ptr{Cint}) libpq
    @c Void setbuf (Ptr{FILE}, Ptr{UInt8}) libpq
    @c Cint setvbuf (Ptr{FILE}, Ptr{UInt8}, Cint, Cint) libpq
    @c Void setbuffer (Ptr{FILE}, Ptr{UInt8}, Cint) libpq
    @c Void setlinebuf (Ptr{FILE},) libpq
    @c Cint fprintf (Ptr{FILE}, Ptr{UInt8}) libpq
    @c Cint printf (Ptr{UInt8},) libpq
    @c Cint sprintf (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint vdprintf (Cint, Ptr{UInt8}, Cint) libpq
    @c Cint dprintf (Cint, Ptr{UInt8}) libpq
    @c Cint fscanf (Ptr{FILE}, Ptr{UInt8}) libpq
    @c Cint scanf (Ptr{UInt8},) libpq
    @c Cint sscanf (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint fscanf (Ptr{FILE}, Ptr{UInt8}) libpq
    @c Cint scanf (Ptr{UInt8},) libpq
    @c Cint sscanf (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint fgetc (Ptr{FILE},) libpq
    @c Cint getc (Ptr{FILE},) libpq
    @c Cint getchar () libpq
    @c Cint getc_unlocked (Ptr{FILE},) libpq
    @c Cint getchar_unlocked () libpq
    @c Cint fgetc_unlocked (Ptr{FILE},) libpq
    @c Cint fputc (Cint, Ptr{FILE}) libpq
    @c Cint putc (Cint, Ptr{FILE}) libpq
    @c Cint putchar (Cint,) libpq
    @c Cint fputc_unlocked (Cint, Ptr{FILE}) libpq
    @c Cint putc_unlocked (Cint, Ptr{FILE}) libpq
    @c Cint putchar_unlocked (Cint,) libpq
    @c Cint getw (Ptr{FILE},) libpq
    @c Cint putw (Cint, Ptr{FILE}) libpq
    @c Ptr{UInt8} fgets (Ptr{UInt8}, Cint, Ptr{FILE}) libpq
    @c Ptr{UInt8} gets (Ptr{UInt8},) libpq
    # @c __ssize_t getdelim (Ptr{Ptr{UInt8}}, Ptr{Cint}, Cint, Ptr{FILE}) libpq
    # @c __ssize_t getline (Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{FILE}) libpq
    @c Cint fputs (Ptr{UInt8}, Ptr{FILE}) libpq
    @c Cint puts (Ptr{UInt8},) libpq
    @c Cint ungetc (Cint, Ptr{FILE}) libpq
    @c Cint fread (Ptr{Void}, Cint, Cint, Ptr{FILE}) libpq
    @c Cint fwrite (Ptr{Void}, Cint, Cint, Ptr{FILE}) libpq
    @c Cint fread_unlocked (Ptr{Void}, Cint, Cint, Ptr{FILE}) libpq
    @c Cint fwrite_unlocked (Ptr{Void}, Cint, Cint, Ptr{FILE}) libpq
    @c Cint fseek (Ptr{FILE}, Clong, Cint) libpq
    @c Clong ftell (Ptr{FILE},) libpq
    @c Void rewind (Ptr{FILE},) libpq
    # @c Cint fseeko (Ptr{FILE}, __off_t, Cint) libpq
    # @c __off_t ftello (Ptr{FILE},) libpq
    # @c Cint fgetpos (Ptr{FILE}, Ptr{fpos_t}) libpq
    # @c Cint fsetpos (Ptr{FILE}, Ptr{fpos_t}) libpq
    @c Void clearerr (Ptr{FILE},) libpq
    @c Cint feof (Ptr{FILE},) libpq
    @c Cint ferror (Ptr{FILE},) libpq
    @c Void clearerr_unlocked (Ptr{FILE},) libpq
    @c Cint feof_unlocked (Ptr{FILE},) libpq
    @c Cint ferror_unlocked (Ptr{FILE},) libpq
    @c Void perror (Ptr{UInt8},) libpq
    @c Cint fileno (Ptr{FILE},) libpq
    @c Cint fileno_unlocked (Ptr{FILE},) libpq
    @c Ptr{FILE} popen (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint pclose (Ptr{FILE},) libpq
    @c Ptr{UInt8} ctermid (Ptr{UInt8},) libpq
    @c Void flockfile (Ptr{FILE},) libpq
    @c Cint ftrylockfile (Ptr{FILE},) libpq
    @c Void funlockfile (Ptr{FILE},) libpq
    @c Ptr{PGconn} PQconnectStart (Ptr{UInt8},) libpq
    @c Ptr{PGconn} PQconnectStartParams (Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}}, Cint) libpq
    @c PostgresPollingStatusType PQconnectPoll (Ptr{PGconn},) libpq
    # @c Ptr{PGconn} PQconnectdb (Ptr{UInt8},) libpq
    @c Ptr{PGconn} PQconnectdb (Ptr{UInt8},) libpq
    @c Ptr{PGconn} PQconnectdbParams (Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}}, Cint) libpq
    @c Ptr{PGconn} PQsetdbLogin (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Void PQfinish (Ptr{PGconn},) libpq
    @c Ptr{PQconninfoOption} PQconndefaults () libpq
    @c Ptr{PQconninfoOption} PQconninfoParse (Ptr{UInt8}, Ptr{Ptr{UInt8}}) libpq
    @c Void PQconninfoFree (Ptr{PQconninfoOption},) libpq
    @c Cint PQresetStart (Ptr{PGconn},) libpq
    @c PostgresPollingStatusType PQresetPoll (Ptr{PGconn},) libpq
    @c Void PQreset (Ptr{PGconn},) libpq
    @c Ptr{PGcancel} PQgetCancel (Ptr{PGconn},) libpq
    @c Void PQfreeCancel (Ptr{PGcancel},) libpq
    @c Cint PQcancel (Ptr{PGcancel}, Ptr{UInt8}, Cint) libpq
    @c Cint PQrequestCancel (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQdb (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQuser (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQpass (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQhost (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQport (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQtty (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQoptions (Ptr{PGconn},) libpq
    @c ConnStatusType PQstatus (Ptr{PGconn},) libpq
    @c PGTransactionStatusType PQtransactionStatus (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQparameterStatus (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQprotocolVersion (Ptr{PGconn},) libpq
    @c Cint PQserverVersion (Ptr{PGconn},) libpq
    @c Ptr{UInt8} PQerrorMessage (Ptr{PGconn},) libpq
    @c Cint PQsocket (Ptr{PGconn},) libpq
    @c Cint PQbackendPID (Ptr{PGconn},) libpq
    @c Cint PQconnectionNeedsPassword (Ptr{PGconn},) libpq
    @c Cint PQconnectionUsedPassword (Ptr{PGconn},) libpq
    @c Cint PQclientEncoding (Ptr{PGconn},) libpq
    @c Cint PQsetClientEncoding (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Ptr{Void} PQgetssl (Ptr{PGconn},) libpq
    @c Void PQinitSSL (Cint,) libpq
    @c Void PQinitOpenSSL (Cint, Cint) libpq
    @c PGVerbosity PQsetErrorVerbosity (Ptr{PGconn}, PGVerbosity) libpq
    @c Void PQtrace (Ptr{PGconn}, Ptr{FILE}) libpq
    @c Void PQuntrace (Ptr{PGconn},) libpq
    @c PQnoticeReceiver PQsetNoticeReceiver (Ptr{PGconn}, PQnoticeReceiver, Ptr{Void}) libpq
    @c PQnoticeProcessor PQsetNoticeProcessor (Ptr{PGconn}, PQnoticeProcessor, Ptr{Void}) libpq
    @c pgthreadlock_t PQregisterThreadLock (pgthreadlock_t,) libpq
    @c Ptr{PGresult} PQexec (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Ptr{PGresult} PQexecParams (Ptr{PGconn}, Ptr{UInt8}, Cint, Ptr{Oid}, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}, Cint) libpq
    @c Ptr{PGresult} PQprepare (Ptr{PGconn}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{Oid}) libpq
    @c Ptr{PGresult} PQexecPrepared (Ptr{PGconn}, Ptr{UInt8}, Cint, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}, Cint) libpq
    @c Cint PQsendQuery (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQsendQueryParams (Ptr{PGconn}, Ptr{UInt8}, Cint, Ptr{Oid}, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}, Cint) libpq
    @c Cint PQsendPrepare (Ptr{PGconn}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{Oid}) libpq
    @c Cint PQsendQueryPrepared (Ptr{PGconn}, Ptr{UInt8}, Cint, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}, Cint) libpq
    @c Cint PQsetSingleRowMode (Ptr{PGconn},) libpq
    @c Ptr{PGresult} PQgetResult (Ptr{PGconn},) libpq
    @c Cint PQisBusy (Ptr{PGconn},) libpq
    @c Cint PQconsumeInput (Ptr{PGconn},) libpq
    @c Ptr{PGnotify} PQnotifies (Ptr{PGconn},) libpq
    @c Cint PQputCopyData (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Cint PQputCopyEnd (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQgetCopyData (Ptr{PGconn}, Ptr{Ptr{UInt8}}, Cint) libpq
    @c Cint PQgetline (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Cint PQputline (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQgetlineAsync (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Cint PQputnbytes (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Cint PQendcopy (Ptr{PGconn},) libpq
    @c Cint PQsetnonblocking (Ptr{PGconn}, Cint) libpq
    @c Cint PQisnonblocking (Ptr{PGconn},) libpq
    @c Cint PQisthreadsafe () libpq
    @c PGPing PQping (Ptr{UInt8},) libpq
    @c PGPing PQpingParams (Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}}, Cint) libpq
    @c Cint PQflush (Ptr{PGconn},) libpq
    @c Ptr{PGresult} PQfn (Ptr{PGconn}, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{PQArgBlock}, Cint) libpq
    @c ExecStatusType PQresultStatus (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQresStatus (ExecStatusType,) libpq
    @c Ptr{UInt8} PQresultErrorMessage (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQresultErrorField (Ptr{PGresult}, Cint) libpq
    @c Cint PQntuples (Ptr{PGresult},) libpq
    @c Cint PQnfields (Ptr{PGresult},) libpq
    @c Cint PQbinaryTuples (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQfname (Ptr{PGresult}, Cint) libpq
    @c Cint PQfnumber (Ptr{PGresult}, Ptr{UInt8}) libpq
    @c Oid PQftable (Ptr{PGresult}, Cint) libpq
    @c Cint PQftablecol (Ptr{PGresult}, Cint) libpq
    @c Cint PQfformat (Ptr{PGresult}, Cint) libpq
    @c Oid PQftype (Ptr{PGresult}, Cint) libpq
    @c Cint PQfsize (Ptr{PGresult}, Cint) libpq
    @c Cint PQfmod (Ptr{PGresult}, Cint) libpq
    @c Ptr{UInt8} PQcmdStatus (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQoidStatus (Ptr{PGresult},) libpq
    @c Oid PQoidValue (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQcmdTuples (Ptr{PGresult},) libpq
    @c Ptr{UInt8} PQgetvalue (Ptr{PGresult}, Cint, Cint) libpq
    @c Cint PQgetlength (Ptr{PGresult}, Cint, Cint) libpq
    @c Cint PQgetisnull (Ptr{PGresult}, Cint, Cint) libpq
    @c Cint PQnparams (Ptr{PGresult},) libpq
    @c Oid PQparamtype (Ptr{PGresult}, Cint) libpq
    @c Ptr{PGresult} PQdescribePrepared (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Ptr{PGresult} PQdescribePortal (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQsendDescribePrepared (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Cint PQsendDescribePortal (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Void PQclear (Ptr{PGresult},) libpq
    @c Void PQfreemem (Ptr{Void},) libpq
    @c Ptr{PGresult} PQmakeEmptyPGresult (Ptr{PGconn}, ExecStatusType) libpq
    @c Ptr{PGresult} PQcopyResult (Ptr{PGresult}, Cint) libpq
    @c Cint PQsetResultAttrs (Ptr{PGresult}, Cint, Ptr{PGresAttDesc}) libpq
    @c Ptr{Void} PQresultAlloc (Ptr{PGresult}, Cint) libpq
    @c Cint PQsetvalue (Ptr{PGresult}, Cint, Cint, Ptr{UInt8}, Cint) libpq
    @c Cint PQescapeStringConn (Ptr{PGconn}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{Cint}) libpq
    @c Ptr{UInt8} PQescapeLiteral (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Ptr{UInt8} PQescapeIdentifier (Ptr{PGconn}, Ptr{UInt8}, Cint) libpq
    @c Ptr{Cuchar} PQescapeByteaConn (Ptr{PGconn}, Ptr{Cuchar}, Cint, Ptr{Cint}) libpq
    @c Ptr{Cuchar} PQunescapeBytea (Ptr{Cuchar}, Ptr{Cint}) libpq
    @c Cint PQescapeString (Ptr{UInt8}, Ptr{UInt8}, Cint) libpq
    @c Ptr{Cuchar} PQescapeBytea (Ptr{Cuchar}, Cint, Ptr{Cint}) libpq
    @c Void PQprint (Ptr{FILE}, Ptr{PGresult}, Ptr{PQprintOpt}) libpq
    @c Void PQdisplayTuples (Ptr{PGresult}, Ptr{FILE}, Cint, Ptr{UInt8}, Cint, Cint) libpq
    @c Void PQprintTuples (Ptr{PGresult}, Ptr{FILE}, Cint, Cint, Cint) libpq
    @c Cint lo_open (Ptr{PGconn}, Oid, Cint) libpq
    @c Cint lo_close (Ptr{PGconn}, Cint) libpq
    @c Cint lo_read (Ptr{PGconn}, Cint, Ptr{UInt8}, Cint) libpq
    @c Cint lo_write (Ptr{PGconn}, Cint, Ptr{UInt8}, Cint) libpq
    @c Cint lo_lseek (Ptr{PGconn}, Cint, Cint, Cint) libpq
    @c Oid lo_creat (Ptr{PGconn}, Cint) libpq
    @c Oid lo_create (Ptr{PGconn}, Oid) libpq
    @c Cint lo_tell (Ptr{PGconn}, Cint) libpq
    @c Cint lo_truncate (Ptr{PGconn}, Cint, Cint) libpq
    @c Cint lo_unlink (Ptr{PGconn}, Oid) libpq
    @c Oid lo_import (Ptr{PGconn}, Ptr{UInt8}) libpq
    @c Oid lo_import_with_oid (Ptr{PGconn}, Ptr{UInt8}, Oid) libpq
    @c Cint lo_export (Ptr{PGconn}, Oid, Ptr{UInt8}) libpq
    @c Cint PQlibVersion () libpq
    @c Cint PQmblen (Ptr{UInt8}, Cint) libpq
    @c Cint PQdsplen (Ptr{UInt8}, Cint) libpq
    @c Cint PQenv2encoding () libpq
    @c Ptr{UInt8} PQencryptPassword (Ptr{UInt8}, Ptr{UInt8}) libpq
    @c Cint pg_char_to_encoding (Ptr{UInt8},) libpq
    @c Ptr{UInt8} pg_encoding_to_char (Cint,) libpq
    @c Cint pg_valid_server_encoding_id (Cint,) libpq
end

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
module DotEnv

"""
`DotEnv.parse` accepts a String or an IOBuffer (Any value that
 can be converted into String), and it will return a Dict with
 the parsed keys and values.
"""
function parse( src )
    res = Dict{String,String}()
    for line in split(String(src), '\n')
        m = match(r"^\s*([\w.-]+)\s*=\s*(.*)?\s*$", line)
        if m != nothing
            key = m.captures[1]
            value = string(m.captures[2])

            if (length(value) > 0 && value[1] === '"' && value[end] === '"')
                value = replace(value, r"\\n"m => "\n")
            end

            value = replace(value, r"(^['\u0022]|['\u0022]$)" => "")

            value = strip(value)

            push!(res, Pair(key, value) )
        end
    end
    res
end

"""
`config` reads your .env file, parse the content, stores it to `ENV`,
and finally return a Dict with the content.
"""
function config( path=".env" )
    if (isfile(path))
        parsed = parse(String(read(path)))

        for (k, v) in parsed
            if( !haskey( ENV, k ) )
                ENV[k] = v
            end
        end

        return parsed
    else
        return nothing
    end
end

config( ;path=".env" ) = config(path)

load(opts...) = config(opts...)

end


text:
data_type = Modification
PostgresORM.get_orm(x::Modification) = return(ModificationORM)
gettablename() = "modification"
get_table_name() = gettablename()
get_schema_name() = "public"
# The format of the mapping is: `property name = "column name"`

const columns_selection_and_mapping = Dict(:id => "id",
                                     :entity_type => "entity_type",
                                     :entity_id => "entity_id",
                                     :attrname => "attrname",
                                     :oldvalue => "oldvalue",
                                     :newvalue => "newvalue",
                                     :appuser_id => "user_id",
                                     :action_id => "action_id",
                                     :action_type => "action_type",
                                     :creation_time => "creation_time")

const id_property = :id

# A dictionnary of mapping between fields symbols and overriding types
#   Left hanside is the field symbol ; right hand side is the type override
types_override = Dict()

const track_changes = false
#
# function create_modification(new_object::Modification)
#     dbconn = opendbconn()
#     result = create_entity!(new_object,
#                            dbconn)
#     closedbconn(dbconn)
#     return result
# end
#
# function update_modification(new_object::Modification)
#     dbconn = opendbconn()
#     result = update_entity!(new_object,
#                            dbconn)
#     closedbconn(dbconn)
#     return result
# end
#
# function retrieve_modification(filter_object::Union{Modification,Missing};
#                           include_users = false)
#     dbconn = opendbconn()
#
#     result = retrieve_entity(filter_object,
#                              data_type,
#                              table_name,
#                              columns_selection_and_mapping,
#                              true, # retrieve_complex_props
#                              dbconn)
#
#
#     closedbconn(dbconn)
#     return result
# end
#
# function retrieve_modification()
#     retrieve_modification(missing)
# end


text:
using Printf

function Base.show(io::IO, model::GebModel)
    println(io, "-"^10)

    for y in 1:model.size
        for x in 1:model.size
            print(io, (model.grid[x, y] === nothing) ? " " : "#")
        end
        println(io)
    end

    println(io)

    for organism in model.organisms
        Base.show(io, organism)
    end
end

function Base.show(io::IO, organism::GebOrganism)
    @printf(io, "(%2i, %2i) \t %3i \t %s \t %i\n", organism.coordinates[1], organism.coordinates[2], organism.direction, organism.genotype, length(organism.rules))
end

function Base.show(io::IO, network::Network)
    apply_to_all(network) do node
        Base.show(io, node)
    end
end

function Base.show(io::IO, node::Node)
    if node.type == InputNode()
        println(io, "INP\t", node.string, " ", [objectid(node)], " ", [objectid(x) for x in node.in_inhibitory], [objectid(x) for x in node.in_excitatory],
                [objectid(x) for x in node.out_inhibitory], [objectid(x) for x in node.out_excitatory], " ", node.inhibitory_output, " ", node.excitatory_output, " ", node.io_value, " ", node.reachable_from_input)
    elseif node.type == OutputNode()
        println(io, "OUT\t", node.string, " ", [objectid(node)], " ", [objectid(x) for x in node.in_inhibitory], [objectid(x) for x in node.in_excitatory],
                [objectid(x) for x in node.out_inhibitory], [objectid(x) for x in node.out_excitatory], " ", node.inhibitory_output, " ", node.excitatory_output, " ", node.io_value, " ", node.reachable_from_input)
    else
        println(io, "HID\t", node.string, " ", [objectid(node)], " ", [objectid(x) for x in node.in_inhibitory], [objectid(x) for x in node.in_excitatory],
                [objectid(x) for x in node.out_inhibitory], [objectid(x) for x in node.out_excitatory], " ", node.inhibitory_output, " ", node.excitatory_output, " ", node.io_value, " ", node.reachable_from_input)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using WeightedStats
using Base.Test

include("accuracy.jl")

include("timer.jl")

include("performance.jl")


text:
using SafeTestsets

@safetestset "Rosenbrock" begin include("rosenbrock.jl") end
@safetestset "AD Tests" begin include("ADtests.jl") end


text:
#==============================================================================
    function to return all neighbors in arbitrary dimensions, from an index
        getneighbors is 2D version
        getnbs is N dimensional version

    tested:     YES
    author:     discourse.julialang.org/u/Crown421
 =============================================================================#

using   BenchmarkTools

function getneighbours(M, cartIdx, radius = 1)
    return [M[i,j] for i in cartIdx[1]-radius:cartIdx[1]+radius, j in cartIdx[2]-radius:cartIdx[2]+radius
           if !((i == cartIdx[1]) & (j == cartIdx[2]))&& (1<= i <= size(M, 1)) && (1<= j <= size(M, 2))]
end

function getnbs(M, startIdx::CartesianIndex; radius = 1)
    idx = CartesianIndices(ntuple(i->(-radius:radius), ndims(M))) .+ startIdx
    idx = idx[:]
    deleteat!(idx, ceil(Int, length(idx)/2))
    filter!(id -> checkbounds(Bool, M, id) , idx)
    return M[idx]
end

M = LinearIndices((300,300))

# M = Int.(t)   # t not defined. not sure intention here, but works without it

testIdx = CartesianIndex((30,30))

nabes = getneighbours(M, testIdx)

display(@benchmark getneighbours(M, testIdx))
display(@benchmark getnbs(M, testIdx))

max_stars_repo_name: mahiki/julia-is-great

max_issues_repo_name: mahiki/julia-is-great

max_forks_repo_name: mahiki/julia-is-great


text:
time_start = time()

using PowerModels
using Ipopt
using JuMP

pkg_load_time = time() - time_start

time_start = time()

file_name = "../data/pglib_opf_case5_pjm.m"
#file_name = "../data/pglib_opf_case118_ieee.m"

data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

data_load_time = time() - time_start

time_start = time()

model = Model(Ipopt.Optimizer)
#set_optimizer_attribute(model, "print_level", 0)

@variable(model, va[i in keys(ref[:bus])])
@variable(model, -ref[:bus][i]["vmax"] <= vr[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)
@variable(model, -ref[:bus][i]["vmax"] <= vi[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=0.0)

@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])

@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])

@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]))

for (i,bus) in ref[:ref_buses]
    @constraint(model, vi[i] == 0)
end

for (i,bus) in ref[:bus]
    bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
    bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]

    @constraint(model, bus["vmin"]^2 <= vr[i]^2 + vi[i]^2)
    @constraint(model, bus["vmax"]^2 >= vr[i]^2 + vi[i]^2)

    @constraint(model,
        sum(p[a] for a in ref[:bus_arcs][i]) ==
        sum(pg[g] for g in ref[:bus_gens][i]) -
        sum(load["pd"] for load in bus_loads) -
        sum(shunt["gs"] for shunt in bus_shunts)*(vr[i]^2 + vi[i]^2)
    )

    @constraint(model,
        sum(q[a] for a in ref[:bus_arcs][i]) ==
        sum(qg[g] for g in ref[:bus_gens][i]) -
        sum(load["qd"] for load in bus_loads) +
        sum(shunt["bs"] for shunt in bus_shunts)*(vr[i]^2 + vi[i]^2)
    )

end

# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
    f_idx = (i, branch["f_bus"], branch["t_bus"])
    t_idx = (i, branch["t_bus"], branch["f_bus"])

    p_fr = p[f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
    q_fr = q[f_idx]                     # q_fr is a reference to the optimization variable q[f_idx]
    p_to = p[t_idx]                     # p_to is a reference to the optimization variable p[t_idx]
    q_to = q[t_idx]                     # q_to is a reference to the optimization variable q[t_idx]

    vr_fr = vr[branch["f_bus"]]         # vm_fr is a reference to the optimization variable vm on the from side of the branch
    vr_to = vr[branch["t_bus"]]         # vm_to is a reference to the optimization variable vm on the to side of the branch
    vi_fr = vi[branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
    vi_to = vi[branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch

    g, b = PowerModels.calc_branch_y(branch)
    tr, ti = PowerModels.calc_branch_t(branch)
    tm = tr^2 + ti^2
    g_fr = branch["g_fr"]
    b_fr = branch["b_fr"]
    g_to = branch["g_to"]
    b_to = branch["b_to"]

    # From side of the branch flow
    @constraint(model, p_fr ==  (g+g_fr)/tm*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) )
    @constraint(model, q_fr == -(b+b_fr)/tm*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) )

    # To side of the branch flow
    @constraint(model, p_to ==  (g+g_to)*(vr_to^2 + vi_to^2) + (-g*tr-b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr+g*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) )
    @constraint(model, q_to == -(b+b_to)*(vr_to^2 + vi_to^2) - (-b*tr+g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr-b*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) )

    # Voltage angle difference limit
    @constraint(model, (vi_fr*vr_to - vr_fr*vi_to) <= tan(branch["angmax"])*(vr_fr*vr_to + vi_fr*vi_to))
    @constraint(model, (vi_fr*vr_to - vr_fr*vi_to) >= tan(branch["angmin"])*(vr_fr*vr_to + vi_fr*vi_to))

    # Apparent power limit, from side and to side
    @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
    @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
end

model_build_time = time() - time_start

#println("model building done")
#println(model)

time_start = time()

optimize!(model)
cost = objective_value(model)

solve_time = time() - time_start

println("")
println("\033[1mSummary\033[0m")
println("   case......: $(file_name)")
println("   cost......: $(round(Int, cost))")
println("   pkg time..: $(pkg_load_time)")
println("   data time.: $(data_load_time)")
println("   build time: $(model_build_time)")
println("   solve time: $(solve_time)")
println("")

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
# compare times of direct calculation of c's compared to Crandall's recursive calculation
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot
using SpecialFunctions
import Base: parse
using Printf

Q = Polylogarithms.Q
c = Polylogarithms.c
max_terms = 10
max_n = 5

s = 5.0001;
r = s - n;
μ = -0.6931471805599453 + 3.141592653589793im #  = log(z) = log(-0.5)
ℒ = log(complex(-μ))  # '\u2112'

result_direct = Array{Any}(undef, max_n, max_terms)
result_crandall = Array{Any}(undef, max_n, max_terms)

for k=0:max_n-1
    for j=1:max_terms
       if j<=3
            # b_direct = @benchmarkable Polylogarithms.c($k, $j, $ℒ)
            # result_crandall[k+1, n_terms] = run(b_direct)
            result_crandall[k+1, n_terms] =  minimum(@benchmark Polylogarithms.c($k, $j, $ℒ))
        end
        # b_crandall = @benchmarkable Polylogarithms.c_crandall($k, $j, $ℒ)
        # result_crandall[k+1, n_terms] = run(b_crandall)
        result_crandall[k+1, n_terms] = minimum(@benchmark Polylogarithms.c_crandall($k, $j, $ℒ))
   end
end


text:
using EquilibriumMeasures, StaticArrays
using Plots

"""
This script finds 9 solutions for the equilibrium measure of which 3 are feasible.
"""

# Potential
V = x -> (x-3)*(x-2)*(1+x)*(2+x)*(3+x)*(2x-1)/20

# Initial guess
ic = SVector(-1,1)

# First solution
a = equilibriummeasure(V; a=ic, returnendpoint=true)[2]

# Keep first solution in an array
solns = [a]

for no_sols = 1:2

    # Compute next solution with the same initial guess but deflating
    # all other known solutions
    b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.5, returnendpoint=true)[2]

    # Append array of known solutions with new solution
    push!(solns, b)

end

# Try different damping
b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.3, returnendpoint=true)[2]
push!(solns, b)

# Try different initial guesses
ic = SVector(2,3)
b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.5, returnendpoint=true)[2]
push!(solns, b)

ic = SVector(1,3)
b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.5, returnendpoint=true)[2]
push!(solns, b)

ic = SVector(-1,4)
b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.5, returnendpoint=true)[2]
push!(solns, b)

b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.4, returnendpoint=true)[2]
push!(solns, b)

b = equilibriummeasure(V; a=ic, knownsolutions=solns, dampening=0.5, returnendpoint=true)[2]
push!(solns, b)

# Plot solutions
p = plot(_equilibriummeasure(V, solns[1]...)[2])
if length(solns) > 1
    for i = 2:length(solns)
        p = plot!(_equilibriummeasure(V, solns[i]...)[2], legend=:top)
    end
end

display(p)
# savefig(p, "em.pdf")

max_stars_repo_name: ApproxFun/EquilbriumMeasures

max_issues_repo_name: ApproxFun/EquilbriumMeasures

max_forks_repo_name: ApproxFun/EquilbriumMeasures


text:
#Sets
#Investemnt Options
set_opt_thermalgenerators = 1:5
set_opt_winds = 1:2
#Existing Units
set_thermalgenerators = 1:2
set_winds = 1:3
set_demands = 1:2
set_nodes = 1:3
set_nodes_ref=1
set_nodes_noref=2:3
#Operating conditions-scenarios-clsuters
set_scenarios= 1:3
set_times= 1:3
#Penalty
P=10000
V=10000
#Assign the values of demands according to their locations in the network
#3 operating conditions below
#two time zones below#Maybe no key means and just montecarlo simulation? I do need to select, let's say 876 (10%)
max_demand=150
R=0.15
#max_demand=150
#max_wind=10
p_D_rowdata=[   0.701	    0.7194
                0.3358	  0.3367
                0.6465	  0.6496]
p_D=p_D_rowdata'

p_D_MW=p_D*max_demand
D= maximum(sum(p_D_MW, dims=1))

#Scenarios
a=1/length(set_scenarios)
γ=zeros(length(set_scenarios))
for s in set_scenarios
γ[s]=a
end

Τ=[100
   100
   100]

#2 wind zones
wind_max=10
wind_rowdata= [0.2027	   0.2027	    0.1803
               0.1429	   0.1429	    0.1803
               0.5865	   0.5865	    0.1803]
wind=wind_rowdata'

wind_opt_rowdata= [0.2027	   0.1803
                  0.1429	   0.1803
                  0.5865	   0.1803]
wind_opt=wind_opt_rowdata'

#ϑ
#Change to the weights from k-means

Ns_H=[5000
      3000
      760]

network=[ 1  2  100  500;
          1  3  100  500;
          2  3  100  500;
          2  1  100  500;
          3  1  100  500;
          3  2  100  500]

start_node=network[:,1]
end_node=network[:,2]
F_max=network[:,3]
B=network[:,4]

n_node=3
n_link=6

links = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links[i] = (start_node[i], end_node[i])
end

links_rev = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links_rev[i] = (end_node[i], start_node[i])
end

F_max_dict=Dict()
B_dict=Dict()

for i=1:n_link
  F_max_dict[(start_node[i], end_node[i])] = F_max[i]
  B_dict[(start_node[i], end_node[i])] = B[i]
end

MapG=[(1,1), (2,2)]
MapG_opt=[(1,2),(2,3),(3,1),(4,2),(5,3)]
#Define this according to zones
MapD=[(1,2),(2,3)]
MapW=[(1,1),(2,2),(3,3)]
MapW_opt=[(1,1),(2,3)]

#Main difference with GAMS is that I calculate the variable cost beforehand
#NREL 2020 Data

#Nuclear and Coal in this case
tech_thermal=[        1500      9.20       7185655         10000            1              80               118988       0                 10.461        0.66        60
                      1000      21.76      4149988         10000            1              100              39695        0.0955             8.638        2.01        75]

tech_thermal_opt=[   1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
                     1000      35.84        6830748        10000            1              0               58242        0.0096            12.507        2.01        75
                     500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
                     500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
                     500       26.85        2778138        10000            1              0               26994        0.0053             7.525        2.81        55]

tech_wind=[   300         0.0               1604886        10000            1             10                43205        0                    0          0.00        30
              300         0.0               1604886        10000            1             10                43205        0                    0          0.00        30
              300         0.0               1604886        10000            1             10                43205        0                    0          0.00        30]

tech_wind_opt=[300         0.0              1604886        10000             1             0                43205        0                    0          0.00        30
               300         0.0              1604886        10000             1             0                43205        0                    0          0.00        30]

 #Indices for set_tech
capacity_per_unit=1
varcost=2
invcost=3
maxBuilds=4
ownership=5
capacity_existingunits=6
fixedcost=7
EmissionsRate=8
HeatRate=9
fuelprice=10
life=11

WACC=0.08

CRF_thermal_opt=zeros(5)

for i in 1:5
CRF_thermal_opt[i]=WACC/(1 - (1 / ( (1 + WACC)^tech_thermal_opt[i,life])))
end

CRF_wind_opt=zeros(2)

for i in 1:2
CRF_wind_opt[i]=WACC/(1 - (1 / ( (1 + WACC)^tech_wind_opt[i,life])))
end

#Compute Var Cost
#2020 Data

#=

tech=[1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
      1000      35.84        6830748        10000            1              0               58242        0.0096            12.507        2.01        75
      500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
      500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
      500       26.85        2778138        10000            1              0               26994        0.0053             7.525        2.81        55
      300       0.0          1604886        10000            1              0               43205        0                    0          0.00        30
      300       0.0          1599902        10000            1              0               18760        0                    0          0.00        30
      1500      9.20         7185655        10000            1             10               118988       0                 10.461        0.66        60
      1000      21.76        4149988        10000            1              8               39695        0.0955             8.638        2.01        75]

TABLE Tech(q,GEN_PARAMS)
                    capacity   var_cost      invcost       maxBuilds      ownership   existingunits        fixedcost    EmissionsRate     HeatRate      FPPRICE     LIFE
*                    [MW]      [$/MWh]       [$/MW]   [Number of Units]                                   [$/MW-yr]    [ton/MMBTU]        [MMBTU/MWh] [$/MMBTU]   [years]
Coal_CCS_30_1        1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
Coal_CCS_90_1        1000      35.84        6830748       10000            1              0               58242        0.0096            12.507        2.01        75
Gas_CT_1             500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
Gas_CC_1             500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
Gas_CC_CCS_1         500       26.85       2778138        10000            1              0               26994        0.0053             7.525        2.81        55
Wind_onshore_1       300         0.0        1604886        10000            1              0               43205        0                    0          0.00        30
Solar_Utility_PV_1   300         0.0        1599902        10000            1              0               18760        0                    0          0.00        30
Nuclear_1            1500      9.20       7185655        10000            1             10               118988       0                 10.461        0.66        60
Coal_new_1           1000      21.76       4149988        10000            1              8                39695       0.0955             8.638        2.01        75
;

=#


text:
using Lehmann

println("Generate DLR ...")
d = DLRGrid(Euv = 1.0, β = 1000.0, rtol = 1e-14, isFermi = false) # Initialize DLR object
# β is the inverse temperature, Euv is the ultraviolt energy cutoff of the Green's function. By definition, Λ = Euv*β
tau_k = d.τ  # DLR imaginary time points

println("Prepare the Green's function sample ...")
G_k = Sample.SemiCircle(d, :τ, tau_k)  # Evaluate known G at tau_k
G_x = tau2dlr(d, G_k) # DLR coeffs from G_k

println("Interpolate imaginary-time Green's function ...")
tau_i = collect(LinRange(0, d.β, 40)) # Equidistant tau grid
G_i = dlr2tau(d, G_x, tau_i) # Evaluate DLR at tau_i

G_s = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k

diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
println("Maximum difference: ", diff(G_i, G_s))
# for i in 1:length(G_i)
#     println("$i   $(G_i[i])         $(G_s[i])")
# end


text:
function implicit_keys()::Bool
    CODES = Dict("KEY" => 1)
    return "KEY" in keys(CODES)
end

function explicit_keys()::Bool
    CODES = Dict("KEY" => 1)
    return "KEY" in keys(CODES)
end

function dict_values()::Bool
    CODES = Dict("KEY" => 1)
    return 1 in values(CODES)
end

function return_dict_index_str(key::String)::Int64
    CODES = Dict("KEY" => 1)
    return CODES[key]
end

function return_dict_index_int(key::Int64)::String
    CODES = Dict(1 => "one")
    return CODES[key]
end

function main()
    @assert(implicit_keys())
    @assert(explicit_keys())
    @assert(dict_values())
    @assert(return_dict_index_str("KEY") == 1)
    @assert(return_dict_index_int(1) == "one")
    println(join(["OK"], " "))
end

main()


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

for T in Base.uniontypes(MPI.MPIDatatype)
    # test vector input
    A = ArrayType{T}([rank + 1])
    C = MPI.Allgather(A, comm)
    @test C isa ArrayType{T,1}
    @test C == ArrayType{T,1}(1:size)

    # test scalar input
    A = convert(T, rank + 1)
    C = MPI.Allgather(A, comm)
    @test C isa Vector{T}
    @test C == Vector{T}(1:size)

    val = T(rank + 1)

    # Test passing output buffer with set size
    A = ArrayType(T[val])

    C = ArrayType{T}(undef, size)
    MPI.Allgather!(A, C, comm) # implied size
    @test C == ArrayType{T}(1:size)

    C = ArrayType{T}(undef, size)
    MPI.Allgather!(A, UBuffer(C,1), comm)
    @test C == ArrayType{T}(1:size)

    # Test assertion error
    C = ArrayType{T}(undef, size-1)
    @test_throws AssertionError MPI.Allgather!(A, C, comm)

    # Test explicit IN_PLACE
    C = ArrayType{T}([i == rank ? i : size + 1 for i = 0:size-1])
    MPI.Allgather!(MPI.IN_PLACE, UBuffer(C, 1), comm)
    @test C isa ArrayType{T,1}
    @test C == ArrayType{T}(0:size-1)

    # Test IN_PLACE
    C = ArrayType{T}([i == rank ? i : size + 1 for i = 0:size-1])
    MPI.Allgather!(UBuffer(C, 1), comm)
    @test C isa ArrayType{T,1}
    @test C == ArrayType{T}(0:size-1)
end

MPI.Finalize()
@test MPI.Finalized()


text:
"""
    Model(; mesh, ρ, D, T₂, κ, γ)

Finite element discretized biological model with initial spin densities `ρ`, diffusion
tensors `D`, T₂-relaxation times `T₂`, wall permeabilities `κ`, and gyromacnetic ratio `γ`.
The vectors `ρ`, `D`, `T₂`, are of length `ncompartment`, while `κ` is of length `nboundary`.
"""
Base.@kwdef struct Model{T}
    mesh::FEMesh{T}
    ρ::Vector{Complex{T}}
    D::Vector{SMatrix{3,3,T,9}}
    T₂::Vector{T}
    κ::Vector{T}
    γ::T
end


text:
#--------------------------------------------------------

function sample(RNG::Random.AbstractRNG,w::Vector)
    w /= sum(w)
    p = rand(RNG)
    runningsum = 0
    for i=1:length(w)
        runningsum += w[i]
        if p < runningsum
            return i
        end
    end
end

sample(w::Vector) = sample(Random.globalRNG(),w)

###---- FUNCTIONS FOR GENERATING CONDITIONED RANDOM WALK ----------

lgamma(x) = first(SpecialFunctions.logabsgamma(x))
lbinom(x) = first(SpecialFunctions.logabsbinom(x))
lbinom(n,k) = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1)
p0(a,c,n::BigInt) = mod(n+c-a,2) ≠ 0 ? 0.0 : exp(lbinom(n,(n+c-a)/2) + n * log(1/2))
p0float(a,c,n) = exp(lbinom(n,(n+c-a)/2) + n * log(1/2))
p0precise(a,c,n) = Float64(binomial(big(n),big(div(n+c-a,2))) * big(1/2)^n)
function p0(a,c,n::Int64)
    if mod(n+c-a,2) ≠ 0 || abs(a-c) > n
        return 0.0
    elseif n < 15_000
        return p0float(a,c,n)
    else
        # Use binomial approxmation
        k = div(n+c-a,2)
        return sqrt(2/(pi*n))*exp(-2*(k-n/2)^2/n)
    end
end
p(a,c,n) = p0(a,c,n) - p0(a,-c,n)
q(a,b,n) = p(b,3,n)*p(a,1,n) - p(b,1,n)*p(a,3,n)

function normalize(v)
    return v / sum(v)
end

function transweights(a,b,n)
    return normalize([q(i,j,n) for i=a-1:2:a+1, j=b-1:2:b+1])
end

function takestep(RNG::Random.AbstractRNG, a, b, n)
    return [(i,j) for j=b-1:2:b+1 for i=a-1:2:a+1][
                            sample(RNG,transweights(a,b,n)[:])]
end

takestep(a,b,n) = takestep(Random.GLOBAL_RNG,a,b,n)

function contour_to_tree(v::Array{Int64,1})
    n = div(length(v)-1,2)
    nbs = [Int64[] for _=1:n+3]
    stack = Int64[n+1]
    heads = Int64[]
    num_ones = 0
    bluetree = Dict{Int64,Int64}()
    for s in diff(v)
        if s == 1
            push!(stack,num_ones+1)
            num_ones += 1
        else
            e = pop!(stack)
            push!(heads,e)
            bluetree[e] = stack[end]
            push!(nbs[e],stack[end])
            push!(nbs[stack[end]],e)
        end
    end
    return PlanarMap(map(reverse,nbs)), bluetree, heads
end

function pair_dyck_paths(RNG::Random.AbstractRNG,
                         n::Integer;
                         verbose=false)
    X = Int64[1]
    Y = Int64[3]
    power_of_ten = 10^floor(Integer,log10(n/2))
    for k=1:n
        (a,b) = takestep(RNG,X[end],Y[end],(n-k))
        push!(X,a)
        push!(Y,b)
    end
    X = X .- 1
    Y = Y .- 3
    return X, Y
end

pair_dyck_paths(n::Integer;kwargs...) =
    pair_dyck_paths(Random.GLOBAL_RNG,n)
#----------------------------------------------------

#--- FUNCTIONS FOR GENERATING WOODS -----------------

import Base: keys, getindex
keys(R::RootedTree) = keys(R.ancestors)
getindex(R::RootedTree,k) = getindex(R.ancestors,k)

function add_red_tree(M::PlanarMap,
                      bluetree::Dict{Int64,Int64},
                      Y::Array{Int64,1},
                      heads::Array{Int64,1})
    nbs = deepcopy(M.nbs)
    headDict = Dict(map(reverse,enumerate(heads)))
    n = (length(Y)-1) ÷ 2
    tails = Int64[]
    num_ones = 0
    redtree = Dict{Int64,Int64}()
    tails = cumsum(diff(Y).==1)[findall(diff(Y) .== -1)]
    tails = [t == n ? n+2 : t+1 for t in tails]
    unmatchedHeadInds = Set(1:length(heads))
    for tail in tails
        matchInd = tail == n+2 ? n : findfirst(x -> x ≥ tail, heads) - 1
        while matchInd ∉ unmatchedHeadInds
            matchInd -= 1
        end
        pop!(unmatchedHeadInds,matchInd)
        head = heads[matchInd]
        redtree[head] = tail
        nbs[head] = insertccw(nbs[head],bluetree[head],tail)
        if tail == n+2
            if length(nbs[tail]) == 0
                nbs[tail] = NeighborCycle([head])
            else
                nbs[tail] = insertcw(nbs[tail],nbs[tail][1],head)
            end
        else
            cwneighbor = bluetree[tail]
            while true
                cwneighbor = cw(M,tail,cwneighbor)
                if (cwneighbor == n+1 ||
                    bluetree[cwneighbor] == tail ||
                        bluetree[tail] == cwneighbor)
                    break
                end
            end
            nbs[tail] = insertccw(nbs[tail],cwneighbor,head)
        end
    end
    PlanarMap(nbs), redtree
end

function find_split(F::Face,
                    bluetree::Dict{Int64,Int64},
                    redtree::Dict{Int64,Int64})
    # finds which vertex in each face has two incident faces outgoing red and blue
    f = [F.elements; F.elements[1:1]]
    for (i,v) in enumerate(f)
        prv = f[i == 1 ? length(f)-1 : i - 1]
        nxt = f[i == length(f)-1 ? 1 : i + 1]
        if v in keys(redtree) && redtree[v] == prv && bluetree[v] == nxt
            return v, prv, nxt
        end
    end
    error("Split not found for face "*string(f))
    (0,0,0)
end

function rotate(a::Array{Int64,1},n::Integer)
    vcat(a[n:end],a[1:n-1])
end

import Base.getindex, Base.findfirst
getindex(Z::Base.Iterators.Zip,k::Integer) = (Z.a[k],Z.b[k])
findfirst(C::Face,t::Tuple) = findfirst(isequal(t),pairs(C))

function add_green_tree(M::PlanarMap,
                        bluetree::Dict{Int64,Int64},
                        redtree::Dict{Int64,Int64})
    n = maximum(keys(bluetree))
    nbs = deepcopy(M.nbs)
    greentree = Dict{Int64,Int64}()
    found = false # whether outer face has been found yet
    FL = faces(M;outeredge=(n+1,1))
    # connect edges on left side of the planar map to
    # the green root:
    F = outerface(FL)
    leftedges = F[1:findfirst(isequal(n+2), F.elements)-1]
    for (i,w) in enumerate(leftedges[2:end])
        greentree[w] = n+3
        nbs[w] = insertcw(nbs[w],leftedges[i],n+3)
        if length(nbs[n+3]) == 0
            nbs[n+3] = NeighborCycle([w])
        else
            nbs[n+3] = insertccw(nbs[n+3],nbs[n+3][end],w)
        end
    end
    # triangulate each interior face with green edges
    for F in interiorfaces(FL)
        if length(F) > 3
            v, prv, nxt = find_split(F,bluetree,redtree)
            # idx stores index in the neighbor list of v
            # where incoming green edges should be inserted
            fc = rotate(F,findfirst(isequal(nxt), F))[2:end-1]
            for (i,w) in enumerate(fc[1:end-1])
                greentree[w] = v
                nbs[v] = insertcw(nbs[v],prv,w)
                nbs[w] = insertccw(nbs[w],fc[i+1],v)
            end
        end
    end
    return PlanarMap(nbs), greentree
end

function UWT(RNG::Random.AbstractRNG,n::Integer)
    X,Y = pair_dyck_paths(RNG,2n)
    M, bluetree, heads = contour_to_tree(X)
    M, redtree = add_red_tree(M,bluetree,Y,heads)
    M, greentree = add_green_tree(M,bluetree,redtree)
    add_edge!(M,n+1,n+2)
    add_edge!(M,n+2,n+3)
    add_edge!(M,n+3,n+1)
    return WoodedTriangulation(M,map(RootedTree,
                    (bluetree,redtree,greentree),
                    (n+1,n+2,n+3))...)
end

UWT(n::Integer) = UWT(Random.GLOBAL_RNG,n)

function descendants(d::Dict{Int64,Int64})
    n = maximum(keys(d))
    desc = zeros(Int64,n)
    for j = 1:n
        k = j
        while true
            k = d[k]
            if k < n+1
                desc[k] += 1
            else
                break
            end
        end
    end
    return desc
end

function flowline(v::Int64,d::Dict{Int64,Int64})
    n = maximum(keys(d))
    branch = Int64[v]
    while branch[end] < n+1
        push!(branch,d[branch[end]])
    end
    return branch
end

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


text:
include("src/instance/instance.jl")
include("src/master/master.jl")
include("src/web/web.jl")

# using Instance
using Web
using Master

Backend.App.start_app()

Life.init_life()

life_iteration = Life.construct_life_iteration(10, 180, 300)
autoclose = false

while length(Life.query) > 0 ||
        length(filter((call) -> call.active, Life.calls)) > 0 ||
        !autoclose # With autoclose disabled the program will run even after
        # finishing all tasks.
    life_iteration()
    sleep(300.000)
end

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
module contact_sim

using JuMP
using Ipopt
using Plots
using WGLMakie
using DifferentialEquations
using Interpolations
using DataStructures: CircularBuffer
# using LCPsolve

include("types.jl")
include("simulator.jl")
include("visualize.jl")
include("lcp_solve.jl")
include("foot_collision.jl")
include("controller.jl")
include("dynamics.jl")
include("contacts.jl")

export solve_lcp,
       form_lcp,
       step,
       render_block_trajectory,
       render_ball_trajectory,
       render_walking_trajectory,
       render_kick_trajectory,
       simulate,
       Object,
       Physics,
       Walker,
       collision_cb,
       controller,
       simulate_walking,
       single_stance!,
       footstrike,
       get_ball_contact_jacobian,
       get_block_contact_jacobian

end