text:
# ps3000_run_streaming_ns.jl

export ps3000_run_streaming_ns

function ps3000_run_streaming_ns(handle :: Int16,
                                 sample_interval,
                                 time_unit,
                                 max_samples,
                                 auto_stop,
                                 no_of_samples_per_aggeregate,
                                 overview_buffer_size = 0x4000)
  ps_status = ccall((:ps3000_run_streaming_ns, ps3000driver), Int16,
    (Int16, UInt32, Int16, UInt32, Int16, UInt32, UInt32),
    handle, sample_interval, time_unit, max_samples, auto_stop,
    no_of_samples_per_aggeregate, overview_buffer_size)
  if ps_status == 0
    error("ps3000 Error: problem occurred or value out of range")
  end
  return nothing
end

function ps3000_run_streaming_ns(handle :: Int16,
                                 sample_interval,
                                 time_unit,
                                 max_samples,
                                 auto_stop::Bool,
                                 no_of_samples_per_aggeregate,
                                 overview_buffer_size = 0x4000)
  ps3000_run_streaming_ns(handle,
                          sample_interval,
                          time_unit,
                          max_samples,
                          auto_stop?0x0001:0x0000,
                          no_of_samples_per_aggeregate,
                          overview_buffer_size)

end


text:
function keyprinter(window::GLFW.Window, event::KeyEvent)
    if isaction("PRESS", event)
        println("Key: $(event.key), scancode: $(event.scancode), action: $(event.action), mods: $(event.mods)")
    end
end

function main_keymap(window::GLFW.Window, event::KeyEvent)
    if iskey("ctrl+a", "press", event)
        GLFW.SetWindowShouldClose(window, true)
    else
        keyprinter(window, event)
    end
end

ext: jl


text:
#Demo 13: Broadcast test: Slicing a Paraboloid
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
a_yvsx = cons(:a, labels = set(yaxis="Y", xaxis="X"))
dfltwfrmattr = cons(:a,
	line = set(style=:none, width=2),
	glyph = set(shape=:o),
)

#==Input data
===============================================================================#
nsteps = 100 #Number of steps to sample x/y values of paraboloid
pext = 10 #sample paraboloid up to x/y = +/-pext

#Paraboloid coefficients
a = 1
b = 2
c = 1

#Axis ranges based on input data:
axrange = cons(:a, xyaxes=set(xmin=-pext, xmax=pext, ymin=-pext, ymax=pext))

#==Computations
===============================================================================#
paraboloid(x,y) = ((x/a)^2 + (y/b)^2) * c

#Generate parameter sweeps:
#NOTE: x is outermost sweep (because outermost sweep plots on x-axis)
sweeplist = PSweep[
	PSweep("y", collect(-pext:(2pext/nsteps):pext))
	PSweep("x", collect(-pext:(2pext/nsteps):pext))
]

para = DataHR{DataFloat}(sweeplist) #Create empty pattern
plane1 = DataHR{DataFloat}(sweeplist) #Create empty pattern
plane2 = DataHR{DataFloat}(sweeplist) #Create empty pattern
for inds in subscripts(para)
	(y,x) = coordinates(para, inds)
	para.elem[inds...] = paraboloid(x,y)
	plane1.elem[inds...] = 3x+y/4+2
	plane2.elem[inds...] = -5(x+10)+5y+50
end

slice_z3 = xcross(para-3) #Cross @ z=3: (x-values@x-coord)
slice_z60 = xcross(para-60) #Cross @ z=60
slice_p1 = xcross(para-plane1) #Cross @ plane1
slice_p2 = xcross(para-plane2) #Cross @ plane1

#Keep x-coordinate, but translate up by y-coordinate:
for dsym in (:slice_z3, :slice_z60, :slice_p1, :slice_p2)
@eval begin
_y = parameter($dsym, "y") #Really only needs to be computed once
$dsym = $dsym*0+_y #Keep x-coord, but translate up by y-coordinate of $dsym
end
end

#==Generate plot
===============================================================================#
p1 = push!(cons(:plot, a_yvsx, axrange, title="Slice with planes z=3 & z=60", legend=false),
	cons(:wfrm, slice_z3, dfltwfrmattr, label="z=3"),
	cons(:wfrm, slice_z60, dfltwfrmattr, label="z=60"),
)
p2 = push!(cons(:plot, a_yvsx, axrange, title="Slice with planes 1 & 2", legend=false),
	cons(:wfrm, slice_p1, dfltwfrmattr, label="plane 1"),
	cons(:wfrm, slice_p2, dfltwfrmattr, label="plane 1"),
)

pcoll = push!(cons(:plotcoll, title="Sliced Paraboloid"), p1, p2)
	pcoll.theme.colorscheme = EasyPlot.ColorScheme(EasyPlot.colormap("Blues", 12))
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using Test
using CSV
using JLD2
using Dates
using DataFrames

@testset "Test OGS" begin
    println("Opening dataset")
    data = CSV.read("input/summary_filtered.csv")
    prior_dict = Dict{String,ttt.Player}()
    for h_key in Set([(row.handicap, row.width) for row in eachrow(data) ])
        prior_dict[string(h_key)] = ttt.Player(ttt.Gaussian(0.,25.0/3.),0.,1.0/100)
    end
    prior_copy = copy(prior_dict)
    results = [row.black_win == 1 ? [0.,1.] : [1., 0.] for row in eachrow(data) ]
    events = [ r.handicap<2 ? [[string(r.white)],[string(r.black)]] : [[string(r.white)],[string(r.black),string((r.handicap,r.width))]] for r in eachrow(data) ]
    times = Vector{Int64}()

    println(now())
    println("Creating History")
    h = ttt.History(events, results, times , prior_dict, mu=0.0,sigma=10.,beta=1.,gamma=0.2125,iter=16)
    ts_log_evidence = ttt.log_evidence(h)
    println(now())
    println("Converging History")
    ttt.convergence(h)
    ttt_log_evidence = ttt.log_evidence(h)
    println(now())

    print("TS: ", ts_log_evidence, ", TTT:", ttt_log_evidence)
    @test ts_log_evidence < ttt_log_evidence
    @test isapprox(ttt_log_evidence, -195388.7125823722)

    w_mean = [ ttt.posterior(h.batches[r], string(data[r,"white"])).mu for r in 1:size(data)[1]]
    b_mean = [ ttt.posterior(h.batches[r], string(data[r,"black"])).mu  for r in 1:size(data)[1]]
    w_std = [ ttt.posterior(h.batches[r], string(data[r,"white"])).sigma for r in 1:size(data)[1]]
    b_std = [ ttt.posterior(h.batches[r], string(data[r,"black"])).sigma for r in 1:size(data)[1]]

    h_mean = [ data[r,"handicap"] > 1 ? ttt.posterior(h.batches[r] ,string((data[r,"handicap"],data[r,"width"]))).mu : 0.0 for r in 1:size(data)[1]]
    h_std = [ data[r,"handicap"] > 1 ? ttt.posterior(h.batches[r] ,string((data[r,"handicap"],data[r,"width"]))).sigma : 0.0 for r in 1:size(data)[1]]
    evidence = [ h.batches[r].events[1].evidence for r in 1:size(data)[1]]

    df = DataFrame(id = data[:"id"]
                  ,white = data[:"white"]
                  ,black = data[:"black"]
                  ,handicap = data[:"handicap"]
                  ,width = data[:"width"]
                  ,w_mean = w_mean
                  ,b_mean = b_mean
                  ,w_std = w_std
                  ,b_std = b_std
                  ,h_mean = h_mean
                  ,h_std = h_std
                  ,evidence = evidence)

    CSV.write("output/longtest_output.csv", df; header=true)

end


text:
function get_microdata(outfile = "SPFmicrodata.xlsx")
  url = "https://www.philadelphiafed.org/-/media/research-and-data/real-time-center/survey-of-professional-forecasters/historical-data/spfmicrodata.xlsx?la=en"
  urldownload(url, parser = identity, save_raw = outfile)
  return nothing
end


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# Function computes H-representation for obstacles given their vertices
# it is assumed that the vertices are given in CLOCK-WISE, and that the first vertex is repeated at the end of the vertex list
###############

function obstHrep(nOb, vOb, lOb)

	# do simple checks
	if nOb != length(lOb)
		println("ERROR in number of obstacles")
	end

	# these matrices contain the H-rep
	A_all = zeros(sum(vOb)-nOb,2)
	b_all = zeros(sum(vOb)-nOb,1)

	# counter for lazy people
	lazyCounter = 1;

	for i = 1 : nOb	# building H-rep

		A_i = zeros(vOb[i]-1,2)
		b_i = zeros(vOb[i]-1,1)

		# take two subsequent vertices, and compute hyperplane
		for j = 1 : vOb[i]-1

			# extract two vertices
			v1 = lOb[i][j]		# vertex 1
			v2 = lOb[i][j+1]	# vertex 2

			# find hyperplane passing through v1 and v2
			if v1[1] == v2[1]	# perpendicular hyperplane, not captured by general formula
				if v2[2] < v1[2]   # line goes "down"
					A_tmp = [1 0]
					b_tmp = v1[1]
				else
					A_tmp = [-1 0]
					b_tmp = -v1[1]
				end
			elseif v1[2] == v2[2] # horizontal hyperplane, captured by general formula but included for numerical stability
				if v1[1] < v2[1]
					A_tmp = [0 1]
					b_tmp = v1[2]
				else
					A_tmp = [0 -1]
					b_tmp = -v1[2]
				end
			else   # general formula for non-horizontal and non-vertical hyperplanes
				ab = [v1[1] 1 ; v2[1] 1] \ [v1[2] ; v2[2]]
				a = ab[1]
				b = ab[2]

				if v1[1] < v2[1]  # v1 --> v2 (line moves right)
					A_tmp = [-a 1]
					b_tmp = b
				else  # v2 <-- v1 (line moves left)
					A_tmp = [a -1]
					b_tmp = -b

				end
			end
			# store vertices
			A_i[j,:] = A_tmp
			b_i[j] = b_tmp
		end

		# store everything
		A_all[lazyCounter : lazyCounter+vOb[i]-2,:] = A_i
		b_all[lazyCounter : lazyCounter+vOb[i]-2] = b_i

		# update counter
		lazyCounter = lazyCounter + vOb[i]-1
	end

	return A_all, b_all

end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
import REPL
function REPL.display(d::REPL.REPLDisplay, mime::MIME"text/plain", x)
    io = REPL.outstream(d.repl)
    write(io, OUTPUT_PROMPT)
    Base.have_color && write(io, REPL.answer_color(d.repl))
    show(IOContext(io, :limit => true), mime, x)
    println(io)
end


text:
@doc """
Symbolic Cholesky factorization.  The fill-reducing permutation
is calculated externally, for example using Metis.NodeND
"""->
type CholSymb{T<:Integer}    # symbolic Cholesky
    pinv::Vector{T}          # fill reducing perm for Chol
    cp::Vector{T}            # column pointers
    parent::Vector{T}        # elimination tree
end

@doc """
Find the nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k))
"""->
function ereach{Tv,Ti}(A::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}}, k, parent)
    if A.uplo != 'U'
        error("ereach requires the upper triangle of A to be stored, i.e. A.uplo == 'U'")
    end
    k = convert(Ti,k)
    Ad = A.data
    Ai = rowvals(Ad)
    s = IntSet([k])                     # non-zeros in row k of L
    v = IntSet([k])                     # visited nodes
    for i in sub(Ai,nzrange(Ad,k))
        if i > k
            continue                    # use only upper triangle
        end
        while !in(i,v)                  # traverse up etree
            push!(s,i)
            push!(v,i)
            i = parent[i]
        end
    end
    collect(s)
end

# find nonzero pattern of Cholesky L[k,1:k-1] using etree and triu(A[:,k])
# based on cs_ereach p. 43, "Direct Methods for Sparse Linear Systems"
function ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti})
    m,n = size(A); Ap = A.colptr; Ai = A.rowval
    s = Ti[]; sizehint!(s, n)            # to be used as a stack
    visited = falses(n)
    visited[k] = true
    for p in Ap[k]:(Ap[k+1] - 1)
        i = Ai[p]                # A[i,k] is nonzero
        if i > k continue end    # only use upper triangular part of A
        while !visited[i]        # traverse up etree
            push!(s,i)           # L[k,i] is nonzero
            visited[i] = true
            i = parent[i]
        end
    end
    s
end

@doc """
Determine the pattern of the sparse Cholesky factor of C

Returns the elimination tree and a lower triangular sparse matrix filled with ones.
"""->
function cholpattern{Tv,Ti}(C::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}})
    C.uplo == 'U' || error("cholpattern requires the upper triangle of C to be stored")
    n = size(C,2)
    parent = etree(C)
    Cd = C.data
    Ci = rowvals(Cd)
    cp = ones(Ti,n+1)
    rv = Ti[]
    s = sizehint!(IntSet([]),n) # nonzeros on column j of L' (row j of L)
    v = sizehint!(IntSet([]),n) # visited nodes
    for j in 1:n
        empty!(s)
        empty!(v)
        push!(s,j)
        push!(v,j)
        for i in sub(Ci,nzrange(Cd,j))
            i > j && break              # upper triangle only
            while !in(i,v)
                push!(s,i)
                push!(v,i)
                i = parent[i]
            end
        end
        append!(rv,collect(s))
        cp[j+1] += convert(Ti,length(rv))
    end
    parent,LowerTriangular(transpose(SparseMatrixCSC(n,n,cp,rv,ones(cp[n+1]))))
end

## write a function that copies the lower triangle of the symmetric matrix to L

function LLcholfact{Tv,Ti}(C::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}})
    parent,L = cholpattern(C)
    n = size(C,2)
    Ct = ltri!(transpose(C.data)) # for now evaluate the lower triangle
    Cti = rowvals(Ct)
    Ctx = nonzeros(Ct)                 # these values will be modified
    Ld = L.data
    Li = rowvals(Ld)
    Lx = fill!(0.,nonzeros(Ld))         # initialize contents of L to zero
    for k in 1:n
        ## spread values of k'th column of A into k'th column of L
        Akr = nzrange(Ct,k)
        Lkr = nzrange(Ld,k)
        Lki = sub(Li,Lkr)
        Lkx = sub(Lx,Lkr)
        if length(Akr) == length(Lkr)
            copy!(Lkx,sub(Ax,Akr))
        else
            Aki = sub(Cti,Akr)
            Akx = sub(Ctx,Akr)
            for i in eachindex(Aki)
                Lkx[searchsortedfirst(Lki,Aki[i])] = Akx[i]
            end
        end
    end
end

function cholsymb{Tv,Ti}(S::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}})
    _,pinv = Metis.nodeND(S)
    C = symperm(S,pinv)
    parent,post = etree(C,true)
    cp = colcounts(C,parent,post)
    CholSymb(pinv,convert(Vector{Ti},cp),parent)
end

function cholfact{Tv,Ti}(A::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}}, S::CholSymb{Ti})
    n = size(A,2)
    cp = S.cp; pinv = S.pinv; parent = S.parent
## maybe change the natural case (i.e. no permutation) to pinv = Array(Ti, 0)
    C  = (pinv == 1:n) ? A : symperm(A, pinv)
    Cp = C.colptr; Ci = C.rowval; Cx = C.nzval
    Lp = copy(cp); Li = Array(Ti, cp[n+1]); Lx = Array(Tv, cp[n+1])
    c  = copy(cp[1:n])
    s  = Array(Ti, n)                   # Ti workspace
    x  = Array(Tv, n)                   # Tv workspace
    for k in 1:n                        # compute L[k,:] for L*L' = C
        ## Nonzero pattern of L[k,:]
        top = cs_ereach(C, k, parent, s, c) # find pattern of L(k,:)
        x[k] = 0                            # x[1:k] is now zero
        for p in Cp[k]:(Cp[k+1]-1)          # x = full(triu(C(:,k)))
            if (Ci[p] <= k) x[Ci[p]] = Cx[p] end
        end
        d = x[k]                        # d = C[k,k]
        x[k] = zero(Tv)                 # clear x for k+1st iteration
        ## Triangular solve
        while top <= n             # solve L[0:k-1,0:k-1] * x = C[:,k]
            i = s[top]             # s[top..n-1] is pattern of L[k,:]
            lki = x[i]/Lx[Lp[i]]   # L[k,i] = x[i]/L[i,i]
            x[i] = zero(Tv)        # clear x for k+1st iteration
            for p in (Lp[i]+1):(c[i]-1)
                x[Li[p]] -= Lx[p] * lki
            end
            d -= lki * lki ;            # d = d - L[k,i]*L[k,i]
            p = c[i]
            c[i] += 1
            Li[p] = k                   # store L[k,i] in column i
            Lx[p] = lki
        end
        ## Compute L[k,k]
        if (d <= 0) error("Matrix is not positive definite, detected at row $k") end
        p = c[k]
        c[k] += 1
        Li[p] = k                 # store L[k,k] = sqrt(d) in column k
        Lx[p] = sqrt(d)
    end
    SparseMatrixCSC(n, n, Lp, Li, Lx)
end


text:
# code for howtowrite

r0101,c0101=simpleprover("howto/howto01_01.cnf",3,3)
r0102,c0102=simpleprover("howto/howto01_02.cnf",3,3)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module hetLK
struct VIterable
    upp::Array{Int}
    low::Array{Int}
    function VIterable(upp::Array{Int},low=fill(0,size(upp)))
        if length(upp)!=length(low) error("upp and low has to have the same length") end
        new(upp,low)
    end
end

function viterate(upp::Vector{Int},low=fill(0,size(upp)))
    Iterators.product(range.(low,upp,step=1)...)
end

# Vector Iterator
function Base.iterate(iter::VIterable,state=(iter.low,1))
    el, count = state
    if count>prod(iter.upp.-iter.low.+1) return nothing end
    nextel=copy(el)
    nextel[1]+=1
    j=1
    while j<length(nextel)
        if nextel[j]>iter.upp[j]
            nextel[j+1]+=1
            nextel[j]=iter.low[j]
        else
            break
        end
        j+=1
    end
    return (el, (nextel,count+1))
end
Base.length(iter::VIterable)=prod(iter.upp.-iter.low.+1)
Base.eltype(iter::VIterable)=Array{Int}

using Distributions
using LinearAlgebra
using StatsFuns
function LKlatent(i::Array{Int,1},n::Array{Int,1},lambda::Array{Float64,1},Rmat::Array{Float64,2},gamma::Float64,p::Array{Float64,1},memo::Dict{Pair{Vector{Int},Vector{Float64}},Float64},linememo::Dict{Pair{Vector{Int},Vector{Int}},Float64})
    if all(p.==1)
        return LK1(i,n,lambda,Rmat,gamma,memo,linememo)
    else
        ll=0.0
        p=fill(0.0,size(i)).+p
        for lat_i in VIterable(n,i)
            ll+=LK1(lat_i,n,lambda,Rmat,gamma,memo,linememo)+sum(logpdf.(Binomial.(lat_i,p),i))
            return ll
        end
    end
end

function LKlatent(i::SubArray{Int,1},n::SubArray{Int,1},lambda::Array{Float64,1},Rmat::Array{Float64,2},gamma::Float64,p::Array{Float64,1},memo::Dict{Pair{Vector{Int},Vector{Float64}},Float64},linememo::Dict{Pair{Vector{Int},Vector{Int}},Float64})
    if all(p.==1)
        return LK1(i,n,lambda,Rmat,gamma,memo,linememo)
    else
        ll=0.0
        p=fill(0.0,size(i)).+p
        for lat_i in VIterable(n,i)
            ll+=LK1(lat_i,n,lambda,Rmat,gamma,memo,linememo)+sum(logpdf.(Binomial.(lat_i,p),i))
            return ll
        end
    end
end

function LKnn(n::Array{Int,1},lambda::Array{Float64,1},Rmat::Matrix{Float64},nfactor::Array{Float64,1},gamma::Float64,memo::Dict{Pair{Vector{Int},Vector{Float64}},Float64})
	if haskey(memo,n=>nfactor) return memo[n=>nfactor] end
    if all(n.==0) return 0.0 end

	#foi=Rmat*k.*nfactor.+lambda

    deduct = [sum(log.(binomial.(n,k)))-sum((n-k).* (Rmat*k.*nfactor.+lambda) )+LKnn(k,lambda,Rmat,nfactor,gamma,memo) for k in VIterable(n) if k!=n]
    lpr=log1mexp(logsumexp(deduct))
	memo[n=>nfactor]=lpr
    return(lpr)
end

function LK1(i::Array{Int,1},n::Array{Int,1},lambda::Array{Float64,1},Rmat::Array{Float64,2},gamma::Float64,memo::Dict{Pair{Vector{Int},Vector{Float64}},Float64},linememo::Dict{Pair{Vector{Int},Vector{Int}},Float64})
    if haskey(linememo,n=>i) return linememo[n=>i] end
    Neff=(Rmat*n-diag(Rmat).*min.(n,1))/(sum(Rmat)/size(Rmat,1))
    nfactor=1.0./Neff.^gamma
	nfactor[isnan.(nfactor)].=0.0
    foi=Rmat*i.*nfactor.+lambda
    rtn=sum(log.(binomial.(n,i)))-sum((n-i).*foi)+LKnn(i,lambda,Rmat,nfactor,gamma,memo)
    linememo[n=>i]=rtn
    return rtn
end

function ll(i::Array{Int,2},n::Array{Int,2},lambda::Array{Float64,1},Rmat::Array{Float64,2},gamma:: Float64=0.0,p::Array{Float64,1}=[1.0])
	memo=Dict{Pair{Vector{Int},Vector{Float64}},Float64}() #memoize LKnn by n => nfactor
    linememo=Dict{Pair{Vector{Int},Vector{Int}},Float64}()
    if size(i)!=size(n) error("sizes of n and i don't match") end
    if !(size(i,1)==size(lambda,1)==size(Rmat,1)==size(Rmat,2)) error("the number of classes don't match") end
    llik=0.0
    for c in 1:size(n,2)
        llik+=LKlatent(i[:,c],n[:,c],lambda,Rmat,gamma,p,memo,linememo)
    end
    return llik
end
end


text:
#!/usr/bin/env julia
using DocOpt, Pope

doc = """
Matter Simulator (aka LJH Rewriter)
Usage:
  mattersimulator.jl <ljhpath> [<output>]
  mattersimulator.jl [--maxchannels=4] [--timeout=0] [--fastforward=3.2] <ljhpath> [<output>]
  mattersimulator.jl <ljhpath> [<output>]
  mattersimulator.jl -h | --help

Options:
  --maxchannels=<mc>   Maximum number of channels to re-write. [default: 240].
  --timeout=<to>       Maximum time (in seconds) to wait between writing successive pulses. Will not maintain synchronization between channels. [default: 0.01]
  --fastforward=<ff>   Reduce the delays between pulses by this factor. Use with large timeout to replay fast, but with synchronization. [default: 1.0]

Help exposition:
  Where <ljhpath> points to a single ljh file, and <output> points directory in which to re-write the ljh files.
  Matter Simulator will find other channels of the same ljh group, and include up to <mc> of them. Then it will start
  to write to new ljh files in directory <output> (created if neccesary). The new ljh files will eventually be identical to the
  input ljh files. However, the pulse records will be written at appoximatley the same rate as they were origninally, based on the
  timestamps in the original ljh files. Reduce <to> to make the new files be written faster.

  Matter simualator will update the sentinel file (in ~/.daq) when it starts and finishes writing.

  If output is not passed, a random temporary directory will be used.

  If you have succesfully run `Pkg.test("Pope") from within Julia, this should work:
  ./mattersimulator.jl ~/.julia/v0.6/ReferenceMicrocalFiles/ljh/20150707_D_chan13.ljh
"""

arguments = docopt(doc, version=v"0.0.1")
maxchannels = parse(Int,arguments["--maxchannels"])
timeout_s = parse(Float64,arguments["--timeout"])
fastforward = parse(Float64,arguments["--fastforward"])
ljhpath = expanduser(arguments["<ljhpath>"])
output = arguments["<output>"]
if output == nothing
  output = joinpath(tempdir(), randstring(8))
else
  output = expanduser(output)
end

ulimit() = a=parse(Int,readstring(`bash -c "ulimit -n"`))
if ulimit() <= 500
  println("WARNING: open file limit too low, run `ulimit -n 1000` before opening julia")
  exist()
end

Pope.mattersim(ljhpath, output, timeout_s, fastforward, maxchannels)


text:
import YAML
config = YAML.load(open("config.yaml"))

# Points at which to create images
parameters = readdlm("parameters.dat")

parnames = config["emulator"]["parameters"]

# Go through each of these and plot images, saving with the filename `image_001.out` etc,
# where each index corresponds to the line number in the parameters file.

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using HDF5

grd = config["grid"]
grid = Grid(grd)

# Load the parameters file from the configuration setting, and make this an object.
pars = convert_dict(config["parameters"], config["model"])
vel = pars.vel # [km/s]
npix = config["npix"] # number of pixels
species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

write_grid("", grid)

# Iterate over rows
for i in 1:size(parameters,1)
    row = parameters[i,:]

    println("processing ", row)

    # Frequency will always be the first parameter
    nu, incl = row

    # Write out the Frequency to synthesize at
    lam = cc./(nu * 1e9) * 1e4 # [μm]

    write_lambda([lam], "")

    # Update the parameters
    pars.incl = incl

    # Write the RADMC input files
    write_model(pars, "", grid, species)

    # Perform the synthesis and save it
    (sizeau_desired, sizeau_command) = size_au(config["size_arcsec"], pars.dpc, grid) # [AU]

    tic()
    run(`radmc3d image incl $(pars.incl) posang $(pars.PA) npix $npix loadlambda sizeau $sizeau_command`)
    println("Synthesis runtime")
    toc()

    # Copy image.out to a new filename
    out = @sprintf("image_%3.3i.out", i)
    cp("image.out", "images/" * out, remove_destination=true)
    println("Finished ", pars)

end

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
using Documenter, Flux, NNlib

DocMeta.setdocmeta!(Flux, :DocTestSetup, :(using Flux); recursive=true)
makedocs(modules=[Flux, NNlib],
         doctest = VERSION >= v"1.4",
         sitename = "Flux",
         pages = ["Home" => "index.md",
                  "Building Models" =>
                    ["Basics" => "models/basics.md",
                     "Recurrence" => "models/recurrence.md",
                     "Regularisation" => "models/regularisation.md",
                     "Model Reference" => "models/layers.md",
                     "Advanced Model Building" => "models/advanced.md",
                     "NNlib" => "models/nnlib.md"],
                  "Handling Data" =>
                    ["One-Hot Encoding" => "data/onehot.md",
                     "DataLoader" => "data/dataloader.md"],
                  "Training Models" =>
                    ["Optimisers" => "training/optimisers.md",
                     "Training" => "training/training.md"],
                  "GPU Support" => "gpu.md",
                  "Saving & Loading" => "saving.md",
                  "The Julia Ecosystem" => "ecosystem.md",
                  "Utility Functions" => "utilities.md",
                  "Performance Tips" => "performance.md",
                  "Datasets" => "datasets.md",
                  "Community" => "community.md"],
         format = Documenter.HTML(
             analytics = "UA-36890222-9",
             assets = ["assets/flux.css"],
             prettyurls = get(ENV, "CI", nothing) == "true"),
         )

deploydocs(repo = "github.com/FluxML/Flux.jl.git",
           target = "build",
           push_preview = true)


text:
using StanModels, CSV

df = filter(row -> row[:age] >= 18,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "Howell1.csv"), delim=';'))

m4_1s = "
// Inferring a Rate
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 20);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

sm = SampleModel("m4.1s", m4_1s);

m4_1_data = Dict("N" => length(df[!, :height]), "h" => df[!, :height]);

(sample_file, log_file) = stan_sample(sm, data=m4_1_data);

if !(sample_file == nothing)
  chn = read_samples(sm)
  # Update parameter names
  chn = set_names(chn, Dict("mu" => "μ", "sigma" => "σ"))
  describe(chn)
end


text:
# WARNING: These tests interfere with the global identity map state, and as such will break
# any Node instances that currently exist.

_test_block() = TimeDag.Block(DateTime(2020, 1, 1):Day(1):DateTime(2020, 1, 5), 1:5)

"""
For tests that care about inspecting the internal identity map state, call this first to
ensure that the identity map is emptied.
"""
function _reset_global_identity_map()
    GC.gc()
    id_map = TimeDag.global_identity_map()
    id_map.weak_to_ref = Dict()
    return nothing
end

@testset "identity mapping" begin
    _reset_global_identity_map()

    block = _test_block()
    n1 = TimeDag.block_node(block)
    n2 = TimeDag.block_node(block)

    @test n1 == n2
    @test n1 === n2

    id_map = TimeDag.global_identity_map()
    @test length(id_map) == 1

    n3 = TimeDag.lag(n1, 1)
    n4 = TimeDag.lag(n1, 1)
    @test n3 == n4
    @test n3 === n4
    @test length(id_map) == 2
end

@testset "weakref" begin
    function _create_and_discard()
        id_map = TimeDag.global_identity_map()
        @test length(id_map) == 0
        n1 = TimeDag.block_node(_test_block())
        @test length(id_map) == 1
        n2 = TimeDag.lag(n1, 1)
        @test length(id_map) == 2
        n3 = TimeDag.lag(n2, 1)
        @test length(id_map) == 3
        return nothing
    end

    _reset_global_identity_map()
    _create_and_discard()
    GC.gc()
    # Manual cleanup invocation required, as otherwise it only occurs when we create a new
    # node.
    TimeDag._cleanup!(TimeDag.global_identity_map())
    @test isempty(TimeDag.global_identity_map())
    @test isempty(TimeDag.global_identity_map().weak_to_ref)
end

@testset "many weakref" begin
    function _create_and_discard()
        count = 1
        start = DateTime(2020)
        block = TimeDag.Block(start:Day(1):(start + Day(count - 1)), 1:count)
        n1 = TimeDag.block_node(block)
        x = n1
        for _ in 1:100
            x += n1
        end
        return nothing
    end

    _reset_global_identity_map()
    _create_and_discard()
    _create_and_discard()
    GC.gc()
    TimeDag._cleanup!(TimeDag.global_identity_map())
    @test isempty(TimeDag.global_identity_map())
    @test isempty(TimeDag.global_identity_map().weak_to_ref)
end

@testset "block with missing" begin
    _reset_global_identity_map()

    block1 = Block([DateTime(2021) => 1, DateTime(2022) => missing])
    block2 = Block([DateTime(2021) => 1, DateTime(2022) => missing])

    @test isequal(block1, block2)
    @test ismissing(block1 != block2)

    @test block_node(block1) === block_node(block1)
    @test block_node(block2) === block_node(block2)
    @test block_node(block1) === block_node(block2)
end


text:
# Parallel instance of fwi_objective function # Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

"""
    fwi_objective(model, source, dobs; options=Options())

Evaluate the full-waveform-inversion (reduced state) objective function. Returns a tuple with function value and \\
gradient. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======

    function_value, gradient = fwi_objective(model, source, dobs)

"""
function fwi_objective(model::Model, source::judiVector, dObs::judiVector; options=Options())
# fwi_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.

    results = judipmap(j -> fwi_objective(model, source[j], dObs[j], subsample(options, j)), 1:dObs.nsrc)

    # Collect and reduce gradients
    obj, gradient = reduce((x, y) -> x .+ y, results)

    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end

function fwi_objective(model::Array{Model,1}, source::Array{judiVector{T,Array{T,2}},1}, dObs::Array{judiVector{T,Array{T,2}},1}; options=Options()) where T
# fwi_objective function for multiple sources and multiple vintages. The function distributes the sources and the input data amongst the available workers.

    results = judipmap((m, q, d) -> fwi_objective(m, q, d; options=options), model, source, dObs)

    obj = sum([results[i][1] for i = 1:length(results)])
    gradient = [results[i][2] for i = 1:length(results)]

    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end


text:
#=
Examples of arrays.
=#

# Array of strings, get first item, and range of items
fruits = ["orange", "melon", "apple", "lemon"]
println("fruits = ", fruits)
println("fruits[1] = ", fruits[1])
println("fruits[2:3] = ", fruits[2:3])

# Replace value of first item
numbers = [4.2, 8, 10, 3, 4]
numbers[1] = 20
println("numbers = $numbers")

# Copy an array
numbers2 = [4.2, 8, 10, 3, 4]
numbers3 = copy(numbers2)
numbers2[1] = 100
println("numbers2 = $numbers2")
println("numbers3 = $numbers3")

# An array comprehension
z = [i * 2 for i in 1:5]
println(z)

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
include("hydrogen_fem_julia.jl")
using Printf
using .Hydrogen_FEM
#using Plots

function main()
    param, val = Hydrogen_FEM.construct()
    eigenval = Hydrogen_FEM.do_run(param, val)

    @printf "計算が終わりました: 基底状態のエネルギー固有値E = %.14f (Hartree)\n" eigenval[1]
    Hydrogen_FEM.save_result(param, eigenval, val)
    @printf "計算結果を%sと%sに書き込みました\n" param.EIGENFUNC_FILENAME param.EIGENVAL_FILENAME
    #plot(val.node_r_glo, val.phi)
end
@time main()


text:
module Napari

export napari

using PyCall
using Conda

# Napari module
const napari = PyNULL()
# PyQt5.QtWidgets.QApplication
const qapp_obj = PyNULL()

const layers = (
    :image,
    :points,
    :labels,
    :shapes,
    :surface,
    :vectors,
    :tracks
)
# We want to precompile the layer commands, so we declare this statically
# but the above constant should be effectively the same as below
# const layers = Symbol.( (napari.layers.NAMES...,) )

include("layers.jl")
include("layer_name_macros.jl")
include("viewer.jl")
include("install.jl")

function __init__(;
    use_qt = parse(Bool, get( ENV, "NAPARI_JL_QT", "true") ),
    use_conda = parse(Bool, get(ENV, "NAPARI_JL_USE_CONDA", "true") )
)
    try
        # We need to coordinate between PyCall and Napari
        # https://github.com/JuliaPy/PyCall.jl/blob/master/src/gui.jl#L140
        # https://github.com/napari/napari/blob/master/napari/_qt/event_loop.py

        if use_conda
            napari_local = pyimport_conda("napari", "napari", "conda-forge")
        else
            napari_local = pyimport("napari")
        end
        copy!(napari,napari_local)

        @info "napari version" version = napari.__version__
        @info dirname(napari.__file__)
        if use_qt
            pygui_start(:qt5)
            # Needed to kickstart Qt5 for some reason
            QApplication = pyimport("qtpy.QtWidgets").QApplication
            qapp_obj_local = QApplication.instance()
            if qapp_obj_local === nothing
            #if isequal( qapp_obj_local , PyObject(nothing) )
                # @info "Creating QApplication"
                # See gui_qt in napari._qt.event_loop
                Qt = pyimport("qtpy.QtCore").Qt
                QApplication.setAttribute(Qt.AA_EnableHighDpiScaling)
                qapp_obj_local = QApplication( [] )
                qapp_obj_local.setApplicationName("napari")
            end
            copy!(qapp_obj, qapp_obj_local)
        else
            @warn """pygui_start(:qt5) not run. Start the event loop manually.
                     You may need to initailize qtpy.QtWidgets.QApplication.
                     See environmental variable NAPARI_JL_QT"""
        end
    catch err
        @warn """Napari.jl has failed to import qtpy and napari from Python.
                 Please make sure these are installed. See the
                 Napari.install_with_pip() and Napari.install_with_conda()
                 methods of this package.
        """
        @debug err
        rethrow(err)
    end
end

function install_at_exit_handler()
    atexit( () -> qapp_obj )
end

function startup_logo(time = 1)
    # imports
    QPixmap = pyimport("qtpy.QtGui").QPixmap
    Qt = pyimport("qtpy.QtCore").Qt
    QSplashScreen = pyimport("qtpy.QtWidgets").QSplashScreen
    # Almost verbatim from napari._qt.eventloop
    logopath = joinpath( dirname(napari.__file__), "resources", "logo.png" )
    pm = QPixmap(logopath).scaled(
        360, 360, Qt.KeepAspectRatio, Qt.SmoothTransformation
    )
    splash_widget = QSplashScreen(pm)
    splash_widget.show()
    qapp_obj._splash_widget = splash_widget
    @async begin
        sleep(time)
        splash_widget.close()
    end
    splash_widget
end

function astronaut()
    pyimport("skimage.data").astronaut()
end

end


text:
import Statistics.mean

include("assemble_cell.jl")
# function doassemble(cellvalues_u::CellVectorValues{dim},
#                     cellvalues_p::CellScalarValues{dim},
#                     facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
#                     dh::DofHandler, bcn::Maybe(Dict), mp::StructArray, rheology::Rheology) where {dim}
#
#     f = zeros(ndofs(dh))
#     assembler = start_assemble(K, f)
#     nu = getnbasefunctions(cellvalues_u)
#     np = getnbasefunctions(cellvalues_p)
#
#     fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
#     Ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix
#
#     # assemble_up dispatch on rheology type parameters
#     assemble_up!(Ke, fe, assembler, dh::DofHandler,
#                 cellvalues_u, cellvalues_p, facevalues_u, grid, mp, rheology, bcn)
#
#     return K, f
# end;

function doassemble!(model::Model{dim,1,Nothing,Nothing,E,Nothing},nbasefuncs) where {dim,E}
    assembler = start_assemble(model.K, model.RHS)

    # Only one primitive variable here
    n = nbasefuncs[1]
    cv = model.cellvalues_tuple[1]

    # initialize local balance equation terms
    fe = zeros(n) # local force vector
    Ke = zeros(n,n) # local stiffness matrix

    @inbounds for (i,cell) in enumerate(CellIterator(model.dofhandler))
        fill!(Ke, 0)
        fill!(fe, 0)

        @timeit "assemble cell" assemble_cell!(Ke, fe, model, cell, cv, n)

        assemble!(assembler, celldofs(cell), fe, Ke)
    end
end

function doassemble!(model::Model{dim,2,Nothing,TV,TE,Nothing}, nbasefuncs, nodal_vars, nodal_vars_prev) where {dim,TV,TE}
    assembler = start_assemble(model.K, model.RHS)
    n1 = nbasefuncs[1]
    n2 = nbasefuncs[2]
    fe = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector
    Ke = PseudoBlockArray(zeros(n1 + n2, n1 + n2), [n1, n2], [n1, n2]) # local stiffness matrix

    # unpack things from model
    dh = model.dofhandler
    cvu, cvp, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
    mp = model.material_properties

    # cache ɛdev outside the element routine to avoid some unnecessary allocations
    ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:n1]
    u▄, p▄ = 1, 2

    @inbounds for cell in CellIterator(dh)
        fill!(Ke, 0)
        fill!(fe, 0)

        # assemble_up dispatch on rheology type parameters
        @timeit "assemble cell" assemble_cell!(Ke, fe, model, cell, ɛdev, cvu, cvp, n1, n2, u▄, p▄, nodal_vars_prev)
        # Assemble local terms to global ones
        assemble!(assembler, celldofs(cell), fe, Ke)
    end
end

function doassemble_elast!(res, K, model::Model{dim,2,Nothing,TV,TE,TP}, nbasefuncs, nodal_vars, nodal_vars_prev) where {dim,TV<:Viscosity,TE,TP<:ViscousDruckerPrager}

    # Initialize res vector and K matrix
    #assembler = start_assemble(model.K, model.RHS)
    assembler = start_assemble(K, res)

    # initialize local residual vector
    n1 = nbasefuncs[1]
    n2 = nbasefuncs[2]
    re = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector
    Ke = PseudoBlockArray(zeros(n1 + n2, n1 + n2), [n1, n2], [n1, n2]) # local stiffness matrix

    # TODO preallocate element nodal variables ue, De ...

    # unpack things from model
    dh = model.dofhandler
    cvu, cvD, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
    mp = model.material_properties

    # blocks names
    field1▄, field2▄ = 1, 2

    @inbounds for cell in CellIterator(dh)
        eldofs = celldofs(cell)
        field1_e = nodal_vars[eldofs[1:n1]]
        field2_e = nodal_vars[eldofs[n1+1:end]]
        #De_prev = exp.(nodal_vars_prev[eldofs[n1+1:end]])

        fill!(re, 0)
        fill!(Ke, 0)
        # assemble_up dispatch on rheology type parameters
        @timeit "assemble cell" assemble_cell_elast!(re, Ke, model, cell, cvu, cvD, n1, n2, field1▄, field2▄, field1_e, field2_e, De_prev)

        # increment global residual with local ones
        assemble!(assembler, eldofs, re, Ke)
    end
end
# function doassemble_res!(model::Model{dim,2,TD,TV,TE,TP}, nodal_vars, nodal_vars_prev) where {dim,TD,TV,TE,TP}
#
#     # Initialize res vector
#     fill!(model.RHS, 0)
#
#     # initialize local residual vector
#     n1, n2 = getnbasefunctions.(model.cellvalues_tuple)
#     re = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector
#
#     # TODO preallocate element nodal variables ue, De ...
#
#     # unpack things from model
#     dh = model.dofhandler
#     cvu, cvD, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
#     mp = model.material_properties
#
#     # blocks names
#     u▄ = 1
#     D▄ = 2
#
#     @inbounds for cell in CellIterator(dh)
#         eldofs = celldofs(cell)
#         ue = nodal_vars[eldofs[1:n1]]
#         De = nodal_vars[eldofs[n1+1:end]]
#         De_prev = nodal_vars_prev[eldofs[n1+1:end]]
#
#         fill!(re, 0)
#         # assemble_up dispatch on rheology type parameters
#         @timeit "assemble cell" assemble_cell!(re, model, cell, cvu, cvp, n1, n2, u▄, D▄, ue, De, De_prev)
#         # increment global residual with local ones
#         assemble!(model.RHS, eldofs, re)
#     end
# end
function doassemble_elast!(res, K, model::Model{dim,2,TD,TV,TE,TP}, nbasefuncs, nodal_vars, nodal_vars_prev) where {dim,TD<:Damage,TV,TE,TP}

    #test#
    # println("ue cell 1 : ")
    # display(nodal_vars[1:22])
    ######
    # Initialize res vector and K matrix
    #assembler = start_assemble(model.K, model.RHS)
    assembler = start_assemble(K, res)

    # initialize local residual vector
    n1 = nbasefuncs[1]
    n2 = nbasefuncs[2]
    re = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector
    Ke = PseudoBlockArray(zeros(n1 + n2, n1 + n2), [n1, n2], [n1, n2]) # local stiffness matrix

    # TODO preallocate element nodal variables ue, De ...

    # unpack things from model
    dh = model.dofhandler
    cvu, cvD, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
    mp = model.material_properties

    # blocks names
    u▄, D▄ = 1, 2

    @inbounds for cell in CellIterator(dh)
        eldofs = celldofs(cell)
        ue = nodal_vars[eldofs[1:n1]]
        #test#
        # if cell.current_cellid.x == 1
        #     println("ue cell 1 : ")
        #     display(ue)
        # end
        ######
        De = exp.(nodal_vars[eldofs[n1+1:end]])
        De_prev = exp.(nodal_vars_prev[eldofs[n1+1:end]])

        fill!(re, zero(eltype(re)))
        fill!(Ke, zero(eltype(Ke)))
        # assemble_up dispatch on rheology type parameters
        @timeit "assemble cell" assemble_cell_elast!(re, Ke, model, cell, cvu, cvD, n1, n2, u▄, D▄, ue, De, De_prev)

        # increment global residual with local ones
        assemble!(assembler, eldofs, re, Ke)
    end
end

function doassemble_AD!(res, K, model::Model{dim,2,TD,TV,TE,TP}, nbasefuncs, nodal_vars, nodal_vars_prev) where {dim,TD,TV,TE,TP}

    assembler = start_assemble(K, res)

    # initialize local residual vector
    n1 = nbasefuncs[1]
    n2 = nbasefuncs[2]
    re = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector
    Ke = PseudoBlockArray(zeros(n1 + n2, n1 + n2), [n1, n2], [n1, n2]) # local stiffness matrix

    # TODO preallocate element nodal variables ue, De ...

    # unpack things from model
    dh = model.dofhandler
    cvu, cvD, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
    mp = model.material_properties

    # blocks names
    u▄, D▄ = 1, 2

    @inbounds for cell in CellIterator(dh)
        fill!(re, zero(eltype(re)))
        fill!(Ke, zero(eltype(re)))

        eldofs = celldofs(cell)
        nodal_vars_el = nodal_vars[eldofs]
        nodal_vars_el_prev = nodal_vars_prev[eldofs]
        # assemble_up dispatch on rheology type parameters
        @timeit "assemble cell AD" assemble_cell_AD!(re, Ke, model, cell, cvu, cvD, n1, n2, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)

        # increment global residual with local ones
        assemble!(assembler, eldofs, re, Ke)
    end
end

function doassemble_res!(res, model::Model{dim,2,TD,TV,TE,TP}, nbasefuncs, nodal_vars, nodal_vars_prev) where {dim,TD,TV,TE,TP}

    # Initialize res vector
    fill!(res, 0)

    # initialize local residual vector
    n1 = nbasefuncs[1]
    n2 = nbasefuncs[2]
    re = PseudoBlockArray(zeros(n1 + n2), [n1, n2]) # local force vector

    # TODO preallocate element nodal variables ue, De ...

    # unpack things from model
    dh = model.dofhandler
    cvu, cvD, fvu  = model.cellvalues_tuple[1], model.cellvalues_tuple[2], model.facevalues
    mp = model.material_properties

    # blocks names
    u▄, D▄ = 1, 2

    @inbounds for cell in CellIterator(dh)
        eldofs = celldofs(cell)
        nodal_vars_el = nodal_vars[eldofs]
        nodal_vars_el_prev = nodal_vars_prev[eldofs]
        # assemble_up dispatch on rheology type parameters
        @timeit "assemble cell" assemble_res_cell!(re, model, cell, cvu, cvD, n1, n2, u▄, D▄, nodal_vars_el, nodal_vars_el_prev)

        # increment global residual with local ones
        assemble!(res, eldofs, re)
    end
end

# 1 unknown only

function doassemble!(model::Model{dim,1,D,V,E,P},nbasefuncs, u ; noplast = false) where {dim,D,V,E,P}
    assembler = start_assemble(model.K, model.RHS)

    # Only one primitive variable here
    n = nbasefuncs[1]
    cv = model.cellvalues_tuple[1]

    # initialize local balance equation terms
    re = zeros(n) # local force vector
    Ke = zeros(n,n) # local stiffness matrix

    @inbounds for (i,cell) in enumerate(CellIterator(model.dofhandler))
        fill!(Ke, 0)
        fill!(re, 0)

        eldofs = celldofs(cell)
        ue = u[eldofs]

        #@timeit "assemble cell"
        assemble_cell!(Ke, re, model, cell, cv, n, ue, noplast)

        assemble!(assembler, eldofs, re, Ke)
    end
end

### individual assembly of res function and derivative K used for linesearch
function doassemble_res!(model::Model{dim,1,D,V,E,P},nbasefuncs, u ; noplast = false) where {dim,D,V,E,P}
    assembler = start_assemble(model.K, model.RHS)

    # Only one primitive variable here
    n = nbasefuncs[1]
    cv = model.cellvalues_tuple[1]

    # initialize local balance equation terms
    re = zeros(n) # local force vector
    Ke = zeros(n,n) # local stiffness matrix

    @inbounds for (i,cell) in enumerate(CellIterator(model.dofhandler))
        fill!(Ke, 0)
        fill!(re, 0)

        eldofs = celldofs(cell)
        ue = u[eldofs]

        #@timeit "assemble cell"
        assemble_cell!(Ke, re, model, cell, cv, n, ue, noplast)

        assemble!(assembler, eldofs, re, Ke)
    end
end

doassemble_K!(model::Model{dim,1,D,V,E,P},nbasefuncs, u ; noplast = false) where {dim,D,V,E,P} = doassemble_K!(model.K, u, model, nbasefuncs ; noplast = false)
function doassemble_K!(K, u::AbstractVector, cellvalues::CellVectorValues{dim},
                    facevalues::FaceVectorValues{dim}, grid::Grid,
                    dh::DofHandler,bcn::Maybe(Dict), mp, states, t; noplast = false) where {dim}

    assembler = start_assemble(K)
    # Only one primitive variable here
    n = nbasefuncs[1]
    cv = getnbasefunctions(cellvalues[1])
    Ke = zeros(n,n) # local stiffness matrix

    @inbounds for (i,cell) in enumerate(CellIterator(model.dofhandler))
        fill!(Ke, 0)

        eldofs = celldofs(cell)
        ue = u[eldofs]

        #@timeit "assemble cell"
        assemble_cell!(Ke, re, model, cell, cv, n, ue, noplast)

        assemble!(assembler, eldofs, re, Ke)
    end
end
function doassemble_K(K::SparseMatrixCSC, u::AbstractVector, cellvalues::CellVectorValues{dim},
                    facevalues::FaceVectorValues{dim}, grid::Grid,
                    dh::DofHandler,bcn::Maybe(Dict), mp, states, t ; elastic_only = false) where {dim}
    assembler = start_assemble(K)
    nu = getnbasefunctions(cellvalues)
    ke = zeros(nu, nu) # element tangent matrix

    @inbounds for (i, cell) in enumerate(CellIterator(dh))
        state, r = states[i], mp[i]
        fill!(ke, 0)
        eldofs = celldofs(cell)
        ue = u[eldofs]
        assemble_cell_K!(ke, cell, cellvalues, facevalues, grid, r, bcn,
                       ue, state, t, nu, elastic_only)
        assemble!(assembler, eldofs, ke)
    end
    return K
end

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

"""
    apply_Neumann_bc!(fe, cell::CellIterator, grid::Grid, bc_dicts::Dict, facevalues, n_basefuncs)

Apply traction boundary conditions on the element force vector `fe`. This is done by looping over the element faces and checking if it belong to a Neumann faceset in 'bc_dicts.neumann'. The traction applied on each relevant face is then computed as the average of the traction function of `(x,t)` evaluated at the face nodes.
Time dependency of the traction function is not allowed for now.

"""
function apply_Neumann_bc!(fe, model, cell::CellIterator, n_basefuncs ; inc_sign = +)
    t = (model.clock isa Clock) ? model.clock.current_time : 0.0
    model.neumann_bc != nothing && @inbounds for face in 1:nfaces(cell)
        try
            onboundary(cell, face)
        catch e
            println("cellid : ", cell.current_cellid)
            println("faceid : ", face)
        end
        if onboundary(cell, face)
            for (neumann_set, traction_func) in pairs(model.neumann_bc)
                if (cellid(cell), face) ∈ getfaceset(model.grid, neumann_set)
                    reinit!(model.facevalues, cell, face)
                    face_nodes_coords = get_face_coordinates(cell::CellIterator, face::Int)
                    tractions = traction_func.(face_nodes_coords,Ref(t))
                    face_traction = mean(hcat(tractions...),dims=2)
                    for q_point in 1:getnquadpoints(model.facevalues)
                        dΓ = getdetJdV(model.facevalues, q_point)
                        for i in 1:n_basefuncs
                            δu = shape_value(model.facevalues, q_point, i)
                            (inc_sign == +) && (fe[i] += (δu ⋅ face_traction) * dΓ)
                            (inc_sign == -) && (fe[i] -= (δu ⋅ face_traction) * dΓ)
                        end
                    end
                end
            end
        end
    end
end


text:
# Don't modify this, it's just there to make the tests work.
module DataScienceTutorials
end


text:
module Modbus
	module libmodbus
		using CBinding

		let
			incdir = joinpath(dirname(@__DIR__), "deps/usr/include/modbus")
			libdir = joinpath(dirname(@__DIR__), "deps/usr/lib")

			c`-I$(incdir) -fparse-all-comments -L$(libdir) -lmodbus`
		end

		const c"struct timeval" = Cvoid
		const c"int8_t"  = Int8
		const c"int16_t" = Int16
		const c"int32_t" = Int32
		const c"int64_t" = Int64
		const c"uint8_t"  = UInt8
		const c"uint16_t" = UInt16
		const c"uint32_t" = UInt32
		const c"uint64_t" = UInt64

		c"""
			#include <modbus.h>
			#include <modbus-rtu.h>
			#include <modbus-tcp.h>
		"""ji
	end

	using .libmodbus
	using Sockets

	export ModbusDevice, ModbusRef

	abstract type ModbusKind end
	struct TCP <: ModbusKind end

	mutable struct ModbusDevice
		ptr::Ptr{modbus_t}
		kind::ModbusKind

		function ModbusDevice(ptr::Ptr{modbus_t}, kind::ModbusKind)
			mb = new(ptr, kind)
			finalizer(mb) do x
				ptr == C_NULL || close(x)
			end
			return mb
		end
	end

	function ModbusDevice(addr::IPAddr, port::Integer)
		ptr = modbus_new_tcp(string(addr), port)
		systemerror("modbus_new_tcp", ptr == C_NULL)
		systemerror("modbus_connect", modbus_connect(ptr) != 0)
		return ModbusDevice(ptr, TCP())
	end

	function Base.close(mb::ModbusDevice)
		modbus_close(mb.ptr)
		modbus_free(mb.ptr)
		mb.ptr = C_NULL
	end

	struct ModbusRef{T}
		mb::ModbusDevice
		addr::Int
		count::Int
		words::Vector{UInt16}

		ModbusRef{T}(mb::ModbusDevice, addr::Int, count::Int = 1, numwords::Int = 2) where {T} = new{T}(mb, addr, count, zeros(UInt16, numwords*count))
	end

	Base.length(ref::ModbusRef) = ref.count
	Base.size(ref::ModbusRef) = (ref.count,)
	Base.eltype(ref::ModbusRef{T}) where {T} = T
	Base.iterate(ref::ModbusRef, state = 1) = state > length(ref) ? nothing : (ref[state], state+1)

	function Base.read(ref::ModbusRef)
		# TODO: break up reads into pages (limited number of words can be read at a time)
		systemerror("modbus_read_registers", modbus_read_registers(ref.mb.ptr, ref.addr, length(ref.words), ref.words) != length(ref.words))
	end

	function Base.fetch(ref::ModbusRef, ind::Int = 1)
		val = UInt32(ref.words[(ind-1)*2 + 1]) | (UInt32(ref.words[(ind-1)*2 + 2]) << 16)
		val = eltype(ref) <: AbstractFloat ?
			reinterpret(Float32, val) :
			eltype(ref) <: Signed ?
				reinterpret(Int32, val) :
				eltype(ref) <: Bool ?
					(val != 0) :
					val
		return eltype(ref)(val)
	end

	function Base.write(ref::ModbusRef, val, ind::Int = 1)
		val = convert(eltype(ref), val)
		val = eltype(ref) <: AbstractFloat ?
			Float32(val) :
			eltype(ref) <: Signed ?
				Int32(val) :
				eltype(ref) <: Bool ?
					(val ? ~zero(UInt32) : zero(UInt32)) :
					UInt32(val)
		val = reinterpret(UInt32, val)
		ref.words[(ind-1)*2 + 1] = val & 0xffff
		ref.words[(ind-1)*2 + 2] = (val >> 16) & 0xffff
	end

	function Base.flush(ref::ModbusRef)
		# TODO: break up writes into pages (limited number of words can be written at a time)
		systemerror("modbus_write_registers", modbus_write_registers(ref.mb.ptr, ref.addr, length(ref.words), ref.words) != length(ref.words))
	end

	function Base.getindex(ref::ModbusRef, ind::Int = 1)
		read(ref)
		return fetch(ref, ind)
	end

	function Base.setindex!(ref::ModbusRef, val, ind::Int = 1)
		write(ref, val, ind)
		flush(ref)
	end
end


text:
module AlgencanWrapper

using Libdl

algencan = Libdl.dlopen("./Algencan/libalgencan")
algencan_optimize = Libdl.dlsym(algencan, "c_algencan")

function unsafe_store_array!(dest, convert_to, src)
    for (i, v) in enumerate(src)
        unsafe_store!(dest, convert_to(v), i)
    end
end

function make_evalf(fun)
    return function (n::Cint, x::Ptr{Cdouble}, f::Ptr{Cdouble}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, Int(n))

        unsafe_store!(f, Cdouble(fun(j_x)))

        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function make_evalg(fun)
    return function (n::Cint, x::Ptr{Cdouble}, g::Ptr{Cdouble}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, Int(n))

        unsafe_store_array!(g, Cdouble, fun(j_x))

        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function make_evalh(fun)
    return function (n::Cint, x::Ptr{Cdouble},
                     hrow::Ptr{Cint}, hcol::Ptr{Cint},
                     hval::Ptr{Cdouble}, hnnz::Ptr{Cint},
                     lim::Cint, lmem::Ptr{Cuchar}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, Int(n))

        j_hrow, j_hcol, j_hval = fun(j_x)
        j_hnnz = length(j_hval)

        if j_hnnz > Int(lim)
            unsafe_store!(lmem, Cuchar(1))
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        unsafe_store_array!(hrow, Cint, j_hrow)
        unsafe_store_array!(hcol, Cint, j_hcol)
        unsafe_store_array!(hval, Cdouble, j_hval)
        unsafe_store!(hnnz, Cint(j_hnnz))

        unsafe_store!(lmem, Cuchar(0))

        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function make_evalc(fun)
    return function (n::Cint, x::Ptr{Cdouble}, ind::Cint, c::Ptr{Cdouble}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_n = Int(n)
        j_ind = Int(ind)# + 1 # convert 0-indexed to 1-indexed

        if j_ind > j_n
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, j_n)

        j_c = fun(j_ind, j_x)

        unsafe_store!(c, Cdouble(j_c))
        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function make_evaljac(fun)
    return function (n::Cint, x::Ptr{Cdouble}, ind::Cint,
                     jcvar::Ptr{Cint}, jcval::Ptr{Cdouble}, jcnnz::Ptr{Cint},
                     lim::Cint, lmem::Ptr{Cuchar}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_n = Int(n)
        j_ind = Int(ind)# + 1 # convert 0-indexed to 1-indexed

        if j_ind > j_n
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, Int(n))

        j_jcvar, j_jcval = fun(j_ind, j_x)
        j_jcnnz = length(j_jcval)

        if j_jcnnz > Int(lim)
            unsafe_store!(lmem, Cuchar(1))
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        unsafe_store_array!(jcvar, Cint, map(x -> x - 1, j_jcvar)) # convert back
        unsafe_store_array!(jcval, Cdouble, j_jcval)
        unsafe_store!(jcnnz, Cint(j_jcnnz))

        unsafe_store!(lmem, Cuchar(0))
        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function make_evalhc(fun)
    return function (n::Cint, x::Ptr{Cdouble}, ind::Cint,
                     hcrow::Ptr{Cint}, hccol::Ptr{Cint},
                     hcval::Ptr{Cdouble}, hcnnz::Ptr{Cint},
                     lim::Cint, lmem::Ptr{Cuchar}, flag::Ptr{Cint})
        if fun == nothing
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_n = Int(n)
        j_ind = Int(ind)# + 1 # convert 0-indexed to 1-indexed

        if j_ind > j_n
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        j_x = unsafe_wrap(Array, x, Int(n))

        j_hcrow, j_hccol, j_hcval = fun(j_ind, j_x)
        j_hcnnz = length(j_hcval)

        if j_hcnnz > Int(lim)
            unsafe_store!(lmem, Cuchar(1))
            unsafe_store!(flag, Cint(-1))
            return nothing
        end

        unsafe_store_array!(hcrow, Cint, map(x -> x - 1, j_hcrow)) # convert back
        unsafe_store_array!(hccol, Cint, map(x -> x - 1, j_hccol)) # convert back
        unsafe_store_array!(hcval, Cdouble, j_hcval)
        unsafe_store!(hcnnz, Cint(j_hcnnz))

        unsafe_store!(lmem, Cuchar(0))
        unsafe_store!(flag, Cint(0))
        return nothing
    end
end

function myevalfc(n::Cint, x::Ptr{Cdouble}, f::Ptr{Cdouble}, m::Cint, c::Ptr{Cdouble}, flag::Ptr{Cint})
    unsafe_store!(flag, Cint(-1))
    return nothing
end

function myevalgjac(n::Cint, x::Ptr{Cdouble}, g::Ptr{Cdouble}, m::Cint, jcfun::Ptr{Cint}, jcvar::Ptr{Cint},
                    jcval::Ptr{Cdouble}, jcnnz::Ptr{Cint}, lim::Cint, lmem::Ptr{Cuchar}, flag::Ptr{Cint})
    unsafe_store!(flag, Cint(-1))
    return nothing
end

function myevalgjacp(n::Cint, x::Ptr{Cdouble}, g::Ptr{Cdouble}, m::Cint, p::Ptr{Cdouble}, q::Ptr{Cdouble},
                     work::Cchar, gotj::Ptr{Cuchar}, flag::Ptr{Cint})
    unsafe_store!(flag, Cint(-1))
    return nothing
end

function myevalhl(n::Cint, x::Ptr{Cdouble}, m::Cint, lambda::Ptr{Cdouble}, scalef::Cdouble,
                  scalec::Ptr{Cdouble}, hlrow::Ptr{Cint}, hlcol::Ptr{Cint}, hlval::Ptr{Cdouble},
                  hlnnz::Ptr{Cint}, lim::Cint, lmem::Ptr{Cuchar}, flag::Ptr{Cint})
    unsafe_store!(flag, Cint(-1))
    return nothing
end

function myevalhlp(n::Cint, x::Ptr{Cdouble}, m::Cint, lambda::Ptr{Cdouble}, scalef::Cdouble,
                   scalec::Ptr{Cdouble}, p::Ptr{Cdouble}, hp::Ptr{Cdouble},
                   goth::Ptr{Cuchar}, flag::Ptr{Cint})
    unsafe_store!(flag, Cint(-1))
    return nothing
end

function optimizeVAZI(;
                  n::Int=1,
                  m::Int=0,
                  f::Function,

                  g::Union{Function, Nothing}=nothing,
                  h::Union{Function, Nothing}=nothing,

                  c::Union{Function, Nothing}=nothing,
                  jac::Union{Function, Nothing}=nothing,
                  hc::Union{Function, Nothing}=nothing,
                  equatn::Union{Array{Int}, Nothing}=nothing,
                  linear::Union{Array{Int}, Nothing}=nothing,

                  x::Union{Array{Float64}, Nothing}=nothing,
                  l::Union{Array{Float64}, Nothing}=nothing,
                  u::Union{Array{Float64}, Nothing}=nothing,
                  lambda::Union{Array{Float64}, Nothing}=nothing,

                  jcnnzmax = 1e8,
                  hnnzmax = 1e8,

                  epsfeas = 1e-8,
                  epsopt = 1e-8,
                  efstain = sqrt(epsfeas),
                  eostain = epsopt^1.5,
                  efacc = sqrt(epsfeas),
                  eoacc = sqrt(epsopt),

                  checkder = 0,

                  outputfnm = "algencan.out",
                  specfnm = "",
                  nvparam = 3,
                  vparam = ["ITERATIONS-OUTPUT-DETAIL 10", "LARGEST-PENALTY-PARAMETER-ALLOWED 1.0e8", "OUTER-ITERATIONS-LIMIT 20"]
                 )
    println("FOI")
    return [], []
end

function optimize(;
                  n::Int=1,
                  m::Int=0,
                  f::Function,

                  g::Union{Function, Nothing}=nothing,
                  h::Union{Function, Nothing}=nothing,

                  c::Union{Function, Nothing}=nothing,
                  jac::Union{Function, Nothing}=nothing,
                  hc::Union{Function, Nothing}=nothing,
                  equatn::Union{Array{Int}, Nothing}=nothing,
                  linear::Union{Array{Int}, Nothing}=nothing,

                  x::Union{Array{Float64}, Nothing}=nothing,
                  l::Union{Array{Float64}, Nothing}=nothing,
                  u::Union{Array{Float64}, Nothing}=nothing,
                  lambda::Union{Array{Float64}, Nothing}=nothing,

                  jcnnzmax = 1e8,
                  hnnzmax = 1e8,

                  epsfeas = 1e-8,
                  epsopt = 1e-8,
                  efstain = sqrt(epsfeas),
                  eostain = epsopt^1.5,
                  efacc = sqrt(epsfeas),
                  eoacc = sqrt(epsopt),

                  checkder = 0,

                  outputfnm = "algencan.out",
                  specfnm = "",
                  nvparam = 3,
                  vparam = ["ITERATIONS-OUTPUT-DETAIL 10", "LARGEST-PENALTY-PARAMETER-ALLOWED 1.0e8", "OUTER-ITERATIONS-LIMIT 20"]
                 )

    myevalf = make_evalf(f)
    myevalg = make_evalg(g)
    myevalh = make_evalh(h)
    myevalc = make_evalc(c)
    myevaljac = make_evaljac(jac)
    myevalhc = make_evalhc(hc)

    myevalf_c = @cfunction($myevalf, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}))
    myevalg_c = @cfunction($myevalg, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}))
    myevalh_c = @cfunction($myevalh, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                                             Ptr{Cint}, Cint, Ptr{Cuchar}, Ptr{Cint}))
    myevalc_c = @cfunction($myevalc, Cvoid, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cint}))
    myevaljac_c = @cfunction($myevaljac, Cvoid, (Cint, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cdouble},
                                                 Ptr{Cint}, Cint, Ptr{Cuchar}, Ptr{Cint}))
    myevalhc_c = @cfunction($myevalhc, Cvoid, (Cint, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Cint, Ptr{Cuchar}, Ptr{Cint}))
    myevalfc_c = @cfunction($myevalfc, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cint}))
    myevalgjac_c = @cfunction($myevalgjac, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint},
                                                   Ptr{Cint},
                                                   Ptr{Cdouble}, Ptr{Cint}, Cint, Ptr{Cuchar}, Ptr{Cint}))
    myevalgjacp_c = @cfunction($myevalgjacp, Cvoid, (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble},
                                                     Ptr{Cdouble}, Cchar, Ptr{Cuchar}, Ptr{Cint}))
    myevalhl_c = @cfunction($myevalhl, Cvoid, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cdouble, Ptr{Cdouble},
                                               Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Cint,
                                               Ptr{Cuchar}, Ptr{Cint}))
    myevalhlp_c = @cfunction($myevalhlp, Cvoid, (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cdouble,
                                                 Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                                                 Ptr{Cuchar}, Ptr{Cint}))

    if (x == nothing)
        x = zeros(n)
    end

    if (l == nothing)
        l = zeros(n)
    end

    if (u == nothing)
        u = ones(n)
    end

    if equatn == nothing
        equatn = zeros(Cuchar, m)
    else
        equatn = map(Cuchar, equatn)
    end

    if linear == nothing
        linear = zeros(Cuchar, m)
    else
        linear = map(Cuchar, linear)
    end

    if lambda == nothing
        lambda = zeros(m)
    end

    coded = map(Cuchar, zeros(11)) # same, for being effectively coded by the user
    if f != nothing
        coded[1] = Cuchar(1)
    end
    if g != nothing
        coded[2] = Cuchar(1)
    end
    if h != nothing
        coded[3] = Cuchar(1)
    end
    if c != nothing
        coded[4] = Cuchar(1)
    end
    if jac != nothing
        coded[5] = Cuchar(1)
    end
    if hc != nothing
        coded[6] = Cuchar(1)
    end

    fval = cnorm = snorm = nlpsupn = inform = 0.0

    ccall(algencan_optimize, Cvoid, (
                            Ptr{Nothing},    # *myevalf,
                            Ptr{Nothing},    # *myevalg,
                            Ptr{Nothing},    # *myevalh,
                            Ptr{Nothing},    # *myevalc,
                            Ptr{Nothing},    # *myevaljac,
                            Ptr{Nothing},    # *myevalhc,
                            Ptr{Nothing},    # *myevalfc,
                            Ptr{Nothing},    # *myevalgjac,
                            Ptr{Nothing},    # *myevalgjacp,
                            Ptr{Nothing},    # *myevalhl,
                            Ptr{Nothing},    # *myevalhlp,
                            Cint,            # jcnnzmax,
                            Cint,            # hnnzmax,
                            Ref{Cdouble},    # *epsfeas,
                            Ref{Cdouble},    # *epsopt,
                            Ref{Cdouble},    # *efstain,
                            Ref{Cdouble},    # *eostain,
                            Ref{Cdouble},    # *efacc,
                            Ref{Cdouble},    # *eoacc,
                            Cstring,         # *outputfnm,
                            Cstring,         # *specfnm,
                            Cint,            # nvparam,
                            Ptr{Ptr{Cuchar}},# **vparam,
                            Cint,            # n,
                            Ref{Cdouble},    # *x,
                            Ref{Cdouble},    # *l,
                            Ref{Cdouble},    # *u,
                            Cint,            # m,
                            Ref{Cdouble},    # *lambda,
                            Ref{Cuchar},     # *equatn,
                            Ref{Cuchar},     # *linear,
                            Ref{Cuchar},     # *coded,
                            Cuchar,          # checkder,
                            Ref{Cdouble},    # *f,
                            Ref{Cdouble},    # *cnorm,
                            Ref{Cdouble},    # *snorm,
                            Ref{Cdouble},    # *nlpsupn,
                            Ref{Cint}        # *inform
                           ),
          myevalf_c, myevalg_c, myevalh_c, myevalc_c, myevaljac_c, myevalhc_c, myevalfc_c,
          myevalgjac_c, myevalgjacp_c, myevalhl_c, myevalhlp_c, jcnnzmax, hnnzmax,
          epsfeas, epsopt, efstain, eostain, efacc, eoacc, outputfnm, specfnm,
          nvparam, vparam, n, x, l, u, m, lambda,
          equatn, linear, coded, checkder, fval, cnorm, snorm,
          nlpsupn, inform
         )

    # OUTPUTS ARE NOT WORKING CORRECTLY (JULIA'S PROBABLY PASSING BY VALUE)

    return x, lambda
end

function unload()
    Libdl.dlclose(algencan)
end

end

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
function testupdateparameters(model,p,v)
    updateparameters!(model,p,v)
    return true
end

@test testupdateparameters(twoEQNK,"rho",0.95)
@test testupdateparameters(twoEQNK,["rho","phi"],[0.9,(0.8,0.6)])


text:
# The following variables need to be defined before `include`-ing this file:
# `testset_name_suffix`
# `n_values`
# `k_values`
# `m_values`

@time @testset "Matrix Multiply Float32 $(testset_name_suffix)" begin
    T = Float32
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Float64 $(testset_name_suffix)" begin
    T = Float64
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Int32 $(testset_name_suffix)" begin
    T = Int32
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Int64 $(testset_name_suffix)" begin
    T = Int64
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end


text:
"""
Pre-fab acset schemas+semagrams schemas
"""
module Examples

export TheoryReactionNet, ReactionNet, ReactionNetSema,
  TheoryDirectedPortGraph, DirectedPortGraph, DirectedPortGraphSema,
  TheoryColoredDPG, ColoredDPG, ColoredDPGSema,
  TheoryCircuitGraph, CircuitGraph, CircuitGraphSema,
  TheoryDDS, DDS, DDSSema,
  TheoryUWD, UWD, UWDSema

using Catlab.Present, Catlab.CSetDataStructures
using JSExpr
using ..Schema
using ..Boxes

@present TheoryReactionNet(FreeSchema) begin
  (T,S,I,O)::Ob
  it::Hom(I,T)
  is::Hom(I,S)
  ot::Hom(O,T)
  os::Hom(O,S)
  N::Data
  L::Data
  rate::Attr(T,N)
  concentration::Attr(S,N)
  species_label::Attr(S,L)
  transitions_label::Attr(T,L)
end

const ReactionNet = ACSetType(TheoryReactionNet)

@semagramschema ReactionNetSema(TheoryReactionNet) begin
  @box S Circle label=:species_label
  @box T Square label=:transitions_label
  @wire I(is,it)
  @wire O(ot,os)
  @data N Numeric
end

@present TheoryDirectedPortGraph(FreeSchema) begin
  Box::Ob
  IPort::Ob
  OPort::Ob
  Wire::Ob
  ibox::Hom(IPort,Box)
  obox::Hom(OPort,Box)
  src::Hom(Wire,OPort)
  tgt::Hom(Wire,IPort)

  String::Data
  label::Attr(Box,String)
end

const DirectedPortGraph = ACSetType(TheoryDirectedPortGraph)

@semagramschema DirectedPortGraphSema(TheoryDirectedPortGraph) begin
  @box Box Square label=:label
  @port IPort(ibox) style="Input"
  @port OPort(obox) style="Output"
  @wire Wire(src,tgt)
  @data String Stringlike
end

@present TheoryColoredDPG <: TheoryDirectedPortGraph begin
  Color::Data
  iportcolor::Attr(IPort,Color)
  oportcolor::Attr(OPort,Color)
  wirecolor::Attr(Wire,Color)
end

const ColoredDPG = ACSetType(TheoryColoredDPG)

color(attr) = js"weights => { return { stroke: weights[$attr] }; }"

@semagramschema ColoredDPGSema(TheoryColoredDPG) begin
  @box Box Square label=:label
  @port IPort(ibox) style="Input" style_fn=color(:iportcolor)
  @port OPort(obox) style="Output" style_fn=color(:oportcolor)
  @wire Wire(src,tgt) style_fn=color(:wirecolor)
  @data String Stringlike
  @data Color Stringlike
end

@present TheoryCircuitGraph(FreeSchema) begin
  Box::Ob
  Port::Ob
  Wire::Ob
  box::Hom(Port,Box)
  src::Hom(Wire,Port)
  tgt::Hom(Wire,Port)

  Resistance::Data
  R::Attr(Wire, Resistance)
  Voltage::Data
  V::Attr(Port, Voltage)
end

const CircuitGraph = ACSetType(TheoryCircuitGraph)

@semagramschema CircuitGraphSema(TheoryCircuitGraph) begin
  @box Box Circle
  @port Port(box) style="Circular"
  @wire Wire(src, tgt)
  @data Resistance Stringlike
  @data Voltage Stringlike
end

@present TheoryDDS(FreeSchema) begin
  X::Ob
  next::Hom(X,X)
end

const DDS = ACSetType(TheoryDDS)

@semagramschema DDSSema(TheoryDDS) begin
  @box X Square
end

@present TheoryUWD(FreeSchema) begin
  Box::Ob
  Port::Ob
  Junction::Ob
  OuterPort::Ob
  box::Hom(Port,Box)
  junction::Hom(Port,Junction)
  outer_junction::Hom(OuterPort, Junction)
end

const UWD = ACSetType(TheoryUWD)

@semagramschema UWDSema(TheoryUWD) begin
  @box Box Circle
  @port Port(box) style="Circular"
  @box Junction TinyCircle
  @box OuterPort SmallCircle
end

end


text:
const NCURSES_VERSION_MAJOR = 5
const NCURSES_VERSION_MINOR = 9
const NCURSES_VERSION_PATCH = 20150214

const NCURSES_VERSION = "$(NCURSES_VERSION_MAJOR).$(NCURSES_VERSION_MINOR)"

const NCURSES_MOUSE_VERSION = 2
const NCURSES_OPAQUE = 1
const NCURSES_REENTRANT = 1
const NCURSES_INTEROP_FUNCS = 1
const NCURSES_TPARM_VARARGS = 1
const NCURSES_EVENT_VERSION	= 1
const NCURSES_EXT_COLORS = 20150214
const CCHARW_MAX	= 5
const TRUE	= 1
const FALSE	= 0
const NCURSES_WIDECHAR	= 1
const NCURSES_ATTR_SHIFT = 8

typealias NCURSES_COLOR_T Cshort
typealias NCURSES_PAIRS_T Cshort
typealias NCURSES_SIZE_T Cshort
typealias NCURSES_TPARM_ARG Clong
# typealias NCURSES_BOOL Cuchar
typealias NCURSES_BOOL Bool # should be compiler's option
typealias chtype Cuint
typealias	attr_t chtype
typealias mmask_t Culong

type cchar_t
  attr::Ptr{attr_t}
  # attr::attr_t
  chars::Array{Cwchar_t, CCHARW_MAX}
  ext_color::Cint
end

typealias NCURSES_CH_T cchar_t # cchar_t

# function NCURSES_CAST(type, value)
#   type(value)
# end
#

COLOR_PAIR( n ) = NCURSES_BITS(n, 0)

function NCURSES_BITS( m, shf)
    m << (shf + NCURSES_ATTR_SHIFT)
end

const A_NORMAL     = uint32(0)
const A_ATTRIBUTES = ~uint32(0)
const A_CHARTEXT   = (NCURSES_BITS(uint32(1),0) - uint32(1))
const A_COLOR      = NCURSES_BITS((uint32(1) << 8) - uint32(1),0)
const A_STANDOUT   = NCURSES_BITS(uint32(1),8)
const A_UNDERLINE  = NCURSES_BITS(uint32(1),9)
const A_REVERSE    = NCURSES_BITS(uint32(1),10)
const A_BLINK      = NCURSES_BITS(uint32(1),11)
const A_DIM        = NCURSES_BITS(uint32(1),12)
const A_BOLD       = NCURSES_BITS(uint32(1),13)
const A_ALTCHARSET = NCURSES_BITS(uint32(1),14)
const A_INVIS      = NCURSES_BITS(uint32(1),15)
const A_PROTECT    = NCURSES_BITS(uint32(1),16)
const A_HORIZONTAL = NCURSES_BITS(uint32(1),17)
const A_LEFT       = NCURSES_BITS(uint32(1),18)
const A_LOW        = NCURSES_BITS(uint32(1),19)
const A_RIGHT      = NCURSES_BITS(uint32(1),20)
const A_TOP        = NCURSES_BITS(uint32(1),21)
const A_VERTICAL   = NCURSES_BITS(uint32(1),22)
const A_ITALIC     = NCURSES_BITS(uint32(1),23)

const WA_ATTRIBUTES	= A_ATTRIBUTES
const WA_NORMAL	= A_NORMAL
const WA_STANDOUT	= A_STANDOUT
const WA_UNDERLINE	= A_UNDERLINE
const WA_REVERSE	= A_REVERSE
const WA_BLINK	=A_BLINK
const WA_DIM	=	A_DIM
const WA_BOLD	=	A_BOLD
const WA_ALTCHARSET=	A_ALTCHARSET
const WA_INVIS	=A_INVIS
const WA_PROTECT	=A_PROTECT
const WA_HORIZONTAL=	A_HORIZONTAL
const WA_LEFT=	A_LEFT
const WA_LOW	=	A_LOW
const WA_RIGHT=	A_RIGHT
const WA_TOP	=	A_TOP
const WA_VERTICAL=	A_VERTICAL
const WA_ITALIC=	A_ITALIC

const COLOR_BLACK   = 0
const COLOR_RED     = 1
const COLOR_GREEN   = 2
const COLOR_YELLOW  = 3
const COLOR_BLUE    = 4
const COLOR_MAGENTA = 5
const COLOR_CYAN    = 6
const COLOR_WHITE   = 7

# not really used. See getmouse() hack
NCURSES_MOUSE_MASK(b,m) = m<<((b-1)*6)

const NCURSES_BUTTON_PRESSED = uint32(2)
const BUTTON1_PRESSED = NCURSES_MOUSE_MASK( 1, NCURSES_BUTTON_PRESSED )
const BUTTON2_PRESSED = NCURSES_MOUSE_MASK( 2, NCURSES_BUTTON_PRESSED )
const BUTTON3_PRESSED = NCURSES_MOUSE_MASK( 3, NCURSES_BUTTON_PRESSED )
const BUTTON4_PRESSED = NCURSES_MOUSE_MASK( 4, NCURSES_BUTTON_PRESSED )
const REPORT_MOUSE_POSITION = NCURSES_MOUSE_MASK( 5, 8 )

const ERR = -1
const OK  = 0

const _SUBWIN = 0x01
const _ENDLINE   = 0x02
const _FULLWIN   = 0x04
const _SCROLLWIN = 0x08
const _ISPAD	   = 0x10
const _HASMOVED  = 0x20
const _WRAPPED   = 0x40

const _NOCHANGE = -1
const _NEWINDEX = -1

type ldat

end

type pdat
  _pad_y::NCURSES_SIZE_T
  _pad_x::NCURSES_SIZE_T
  _pad_top::NCURSES_SIZE_T
  _pad_left::NCURSES_SIZE_T
  _pad_bottom::NCURSES_SIZE_T
  _pad_right::NCURSES_SIZE_T
end

immutable Window

  _cury::NCURSES_SIZE_T
  _curx::NCURSES_SIZE_T

  _maxy::NCURSES_SIZE_T
  _maxx::NCURSES_SIZE_T
  _begy::NCURSES_SIZE_T
  _begx::NCURSES_SIZE_T

  _flags::Cshort
  _attrs::attr_t
  _bkgd::chtype

  _notimeout::NCURSES_BOOL
  _clear::NCURSES_BOOL
  _leaveok::NCURSES_BOOL
  _scroll::NCURSES_BOOL
  _idlok::NCURSES_BOOL
  _idcok::NCURSES_BOOL
  _immed::NCURSES_BOOL
  _sync::NCURSES_BOOL
  _use_keypad::NCURSES_BOOL
  _delay::Cint

  _line::Ptr{ldat}

  _regtop::NCURSES_SIZE_T
  _regbottom::NCURSES_SIZE_T

  _parx::Cint
  _pary::Cint
  _parent::Ptr{Void}

  _pad::pdat
  _yoffset::NCURSES_SIZE_T
  _bkgrnd::cchar_t
  _color::Cint

end

typealias WINDOW Window

# type SCREEN
#   # screen
# end

export NCURSES_VERSION_MAJOR, NCURSES_VERSION_MINOR, NCURSES_VERSION_PATCH, NCURSES_VERSION,NCURSES_MOUSE_VERSION, NCURSES_OPAQUE, NCURSES_REENTRANT, NCURSES_INTEROP_FUNCS, NCURSES_TPARM_VARARGS, NCURSES_EXT_COLORS, CCHARW_MAX, TRUE, FALSE, NCURSES_WIDECHAR, NCURSES_ATTR_SHIFT, NCURSES_COLOR_T, NCURSES_PAIRS_T, NCURSES_SIZE_T, NCURSES_TPARM_ARG, NCURSES_BOOL, chtype, attr_t, mmask_t, cchar_t, NCURSES_CH_T, COLOR_PAIR, A_NORMAL, A_ATTRIBUTES, A_CHARTEXT, A_COLOR, A_STANDOUT, A_UNDERLINE, A_REVERSE, A_BLINK, A_DIM, A_BOLD, A_ALTCHARSET, A_INVIS, A_PROTECT, A_HORIZONTAL, A_LEFT, A_LOW, A_RIGHT, A_TOP, A_VERTICAL, A_ITALIC, WA_ATTRIBUTES, WA_NORMAL, WA_STANDOUT, WA_UNDERLINE, WA_REVERSE, WA_BLINK, WA_DIM, WA_BOLD, WA_ALTCHARSET, WA_INVIS, WA_PROTECT, WA_HORIZONTAL, WA_LEFT, WA_LOW, WA_RIGHT, WA_TOP, WA_VERTICAL, WA_ITALIC, COLOR_BLACK, COLOR_RED, COLOR_GREEN, COLOR_YELLOW, COLOR_BLUE, COLOR_MAGENTA, COLOR_CYAN, COLOR_WHITE, NCURSES_MOUSE_MASK, NCURSES_BUTTON_PRESSED, BUTTON1_PRESSED, BUTTON2_PRESSED, BUTTON3_PRESSED, BUTTON4_PRESSED, REPORT_MOUSE_POSITION, ERR, OK, _SUBWIN, _ENDLINE, _FULLWIN, _SCROLLWIN, _ISPAD, _HASMOVED, _WRAPPED, _NOCHANGE, _NEWINDEX, ldat, pdat, _win_st, WINDOW, _NC_EVENT_TIMEOUT_MSEC, _NC_EVENT_FILE, _NC_EVENT_FILE_READABLE, _NC_EVENT_FILE_WRITABLE, _NC_EVENT_FILE_EXCEPTION, _nc_event, _nc_eventlist, wgetch_events, wgetch_events


text:
left_margin=50px

using CSV
using DelimitedFiles
using DataFrames
using HDF5
using Plots
using Dates
using StatsBase
using LaTeXStrings
using LinearAlgebra
using Polynomials
using Optim
using Distributions
using LsqFit
using DifferentialEquations
#using StatsPlots
using Measures
Pkg.update()

# in utc:
#=Tfid=h5open("../data/other/CE_data_collected_by_johannes/input_actual_timezone_utc.h5","r")
Tdgr=Tfid["df"]
#println(names(dgr))
Tdata=read(Tdgr["block0_values"])
Tsces=read(Tdgr["block0_items"])
Tsce1=read(Tdgr["axis0"])
Tsce2=read(Tdgr["axis1"])
Tsce2=unix2datetime.(Tsce2[:]/1000000000)=#

#in timezone berlin:
fid=h5open("../data/other/CE_data_collected_by_johannes/input_actual_timezone_berlin.h5","r")
dgr=fid["df"]
#println(names(dgr))
data=read(dgr["block0_values"])
sces=read(dgr["block0_items"])
sce1=read(dgr["axis0"])
sce2=read(dgr["axis1"])
sce2=unix2datetime.(sce2[:]/1000000000)
#sce2[35065+0:(35065+364*24)]
sce2[26305+0:(26305+364*24+23)]
startd=26305

# 1) use goods2 from model_atPeak/parameter_check_mle

# 2)now filter the intervals that ar both in goods2 as well as complete in Johannes' dataset
cort=zeros(length(goods2))
nanent=zeros(length(goods2))
j=1
nj=1
for i in goods2
    if maximum(isnan.(data[:,startd+i]))==false #iff no Nan entry in this column
        cort[j]=i
        global j+=1
    else
        nanent[nj]=i
        global nj+=1
    end
end
cor=convert.(Int64,cort[1:j-1]) #alle, die rausfliegen
nanen=convert.(Int64,nanent[1:nj-1]) # alle intervalle die gehen

# calculate the correlations
spearman=zeros(49,5)
for jo in 1:49
    for mle in 1:5
        spearman[jo,mle]=corspearman(data[jo,cor.+startd],mleres[cor,2+mle])
    end
end
spearman

# make the heatmap
scespl=["\$\\textrm{load}\$",
    "\$\\textrm{gen \\ biomas}\$",
    "\$\\textrm{gen \\ lignite}\$",
    "\$\\textrm{gen \\ coal-gas}\$",
    "\$\\textrm{gen \\ gas}\$",
    "\$\\textrm{gen \\ hardcoal}\$",
    "\$\\textrm{gen \\ oil}\$",
    "\$\\textrm{gen \\ geothermal}\$",
    "\$\\textrm{gen \\ pumped \\ hydro}\$",
    "\$\\textrm{gen \\ run \\ off \\ hydro}\$",
    "\$\\textrm{gen \\ reservoir \\ hydro}\$",
    "\$\\textrm{gen \\ nuclear}\$",
    "\$\\textrm{gen \\ other \\ renewables}\$",
    "\$\\textrm{gen \\ solar}\$",
    "\$\\textrm{gen \\ waste}\$",
    "\$\\textrm{gen \\ wind \\ off}\$",
    "\$\\textrm{gen \\ wind \\ on}\$",
    "\$\\textrm{gen \\ other}\$",
    "\$\\textrm{gen \\ total}\$",
    "\$\\textrm{gen \\ synchr.}\$",
    "\$\\textrm{load \\ ramp}\$",
    "\$\\textrm{total \\ gen \\ ramp}\$",
    "\$\\textrm{biomass \\ ramp}\$",
    "\$\\textrm{lignite \\ ramp}\$",
    "\$\\textrm{coal \\ gas \\ ramp}\$",
    "\$\\textrm{gas \\ ramp}\$",
    "\$\\textrm{hard \\ coal \\ ramp}\$",
    "\$\\textrm{oil \\ ramp}\$",
    "\$\\textrm{geothermal \\ ramp}\$",
    "\$\\textrm{pumped \\ hydro \\ ramp}\$",
    "\$\\textrm{run \\ off \\ hydro \\ ramp}\$",
    "\$\\textrm{reservoir \\ hydro \\ ramp}\$",
    "\$\\textrm{nuclear \\ ramp}\$",
    "\$\\textrm{other \\ renewables \\ ramp}\$",
    "\$\\textrm{solar \\ ramp}\$",
    "\$\\textrm{waste \\ ramp}\$",
    "\$\\textrm{wind \\ off \\ ramp}\$",
    "\$\\textrm{wind \\ on \\ ramp}\$",
    "\$\\textrm{other \\ ramp}\$",
    "\$\\textrm{fe \\ wind \\ on}\$",
    "\$\\textrm{fe \\ wind \\ off}\$",
    "\$\\textrm{fe \\ solar}\$",
    "\$\\textrm{fe \\ total \\ gen}\$",
    "\$\\textrm{fe \\ load}\$",
    "\$\\textrm{fe \\ load \\ ramp}\$",
    "\$\\textrm{fe \\ total \\ gen \\ ramp}\$",
    "\$\\textrm{fe \\ wind \\ off \\ ramp}\$",
    "\$\\textrm{fe \\ wind \\ on \\ ramp}\$",
    "\$\\textrm{fe \\ solar \\ ramp}\$"]

Plots.heatmap(["\$\\gamma_1\$","\$\\gamma_2\$","\$q\$","\$r\$", "\$\\textrm{D}\$"],
    scespl,
    spearman,
    ticks=:all,
    left_margin=45pt,
    c= :diverging_bwr_40_95_c42_n256,
    size=(800,1200)
)
png("/Users/raphaelbiertz/Documents/masterarbeit/coding/codes/cor_analysis/spearman3")

writedlm("C:/Users/Raphael Biertz/Desktop/spearman.csv",spearman)
writedlm("C:/Users/Raphael Biertz/Desktop/spear.csv",sces)


text:
module DiskData

using Lazy

const AVector = AbstractVector

include("iter/iter.jl")
include("iter/split.jl")
include("vectors/chunked.jl")
include("vectors/cached.jl")
include("vectors/disk.jl")
include("vectors/bigvector.jl")

end # module


text:
using Zygote, NNlib, Test
using Zygote: gradient
using NNlib: conv
import Random

function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

gradcheck(f, xs...) =
  all(isapprox.(ngradient(f, xs...),
                gradient(f, xs...), rtol = 1e-5, atol = 1e-5))

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

Random.seed!(0)

@test gradtest((x, W, b) -> σ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> σ.(W*x .+ b), (5,3), (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), (5,3), (2,5), 2)

@test gradtest((w, x) -> w'*x, randn(10, 2), randn(10))
@test gradtest((w, x) -> transpose(w)*x, randn(5,5), randn(5,5))

@test gradtest(x -> sum(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x), (3,4,5))

@test gradtest(x -> softmax(x).*(1:3), 3)
@test gradtest(x -> softmax(x).*(1:3), (3,5))
@test gradtest(x -> logsoftmax(x).*(1:3), 3)
@test gradtest(x -> logsoftmax(x).*(1:3), (3,5))

@test gradtest(conv, rand(10, 3, 2), randn(Float64,2, 3, 2))
@test gradtest(conv, rand(10, 10, 3, 2), randn(Float64,2, 2, 3, 2))
@test gradtest(conv, rand(10, 10, 10, 3, 2), randn(Float64,2, 2, 2, 3, 2))

@test gradtest(x -> maxpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> maxpool(x, (2,2,2)), rand(10, 10, 10, 3, 2))

@test gradtest(x -> meanpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> meanpool(x, (2,2,2)), rand(5, 5, 5, 3, 2))

@test gradtest(x -> permutedims(x, [3,1,2]), rand(4,5,6))

@test gradtest(x -> repeat(x; inner=2), rand(5))
@test gradtest(x -> repeat(x; inner=2, outer=3), rand(5))
@test gradtest(x -> repeat(x; inner=(2,2,1), outer=(1,1,3)), rand(5,4,3))

@test gradtest(kron, rand(5), rand(3))
@test gradtest(kron, rand(5), rand(3), rand(8))
@test gradtest(kron, rand(5,1), rand(3,1))
@test gradtest(kron, rand(5,1), rand(3,1), rand(8,1))
@test gradtest(kron, rand(5,2), rand(3,2), rand(8,2))

@testset "maximum" begin
  @test gradtest(maximum, rand(2, 3))

  @test gradtest(x -> maximum(x, dims=1), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=2), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> maximum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "minimum" begin
  @test gradtest(minimum, rand(2, 3))

  @test gradtest(x -> minimum(x, dims=1), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=2), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> minimum(x, dims=[1, 2]), rand(2, 3, 4))
end


text:
using Documenter
using RATiLQR

format =
Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
)

makedocs(
    format = format,
    sitename = "RATiLQR.jl",
    authors = "Haruki Nishimura",
    source = "source",
    build = "build",
    clean = true,
    doctest = true,
    modules = [RATiLQR],
    pages = ["Home" => "index.md",
             "getting-started.md",
             "optimal-control.md",
             "solvers.md",
             "references.md"]
)

deploydocs(
    repo = "github.com/StanfordMSL/RATiLQR.jl.git"
)


text:
using eFEM, Test

mesh1 = squareMesh([-2,2,-1,1],12,1)
OperatorType = :Poisson2D

dNodes = Dirichlet(:left,:right,:bottom,:top)
dBCf = Dirichlet((x,y) -> x^4*y^4)
ff   = Forcing((x,y) -> -12*x^2*y^2*(x^2+y^2))

Nodes = [dNodes]; bcfun = [dBCf,ff]
prob = Problem(mesh1,Nodes,bcfun,OperatorType)
sol = solve(prob,mesh1)

L1   = DomainNorm(mesh1.xy,mesh1.cm,sol.u;normID="1")
L2   = DomainNorm(mesh1.xy,mesh1.cm,sol.u;normID="2")
Linf = DomainNorm(mesh1.xy,mesh1.cm,sol.u;normID="Inf")
hCalc(mesh1)

mesh2 = squareMesh([-2,2,-1,1],12,2)
OperatorType = :Poisson2D

dNodes = Dirichlet(:left,:right,:bottom,:top)
dBCf = Dirichlet((x,y) -> x^4*y^4)
ff   = Forcing((x,y) -> -12*x^2*y^2*(x^2+y^2))

Nodes = [dNodes]; bcfun = [dBCf,ff]
prob = Problem(mesh2,Nodes,bcfun,OperatorType)
sol = solve(prob,mesh2)

L1   = DomainNorm(mesh2.xy,mesh2.cm,sol.u;normID="1")
L2   = DomainNorm(mesh2.xy,mesh2.cm,sol.u;normID="2")
Linf = DomainNorm(mesh2.xy,mesh2.cm,sol.u;normID="Inf")
hCalc(mesh2)


text:
####################################
# Rational Krylov utility functions
####################################

export lejabagby, scgendivdiffs, ratnewtoncoeffs, ratnewtoncoeffsm

"""
Generate Leja-Bagby points (a,b) on (A,B), with scaling factors β such that
the uniform norm on the control set C is 1. Greedy search for a minimum is
performed on B. If keepA is true then the points in the output a will be exactly
those of A, otherwise the points in a are also chosen via greedy search on A. If
forceInf is a positive integer, the first forceInf poles in b will be infinity.
"""
function lejabagby(A::AbstractVector{CT}, B::AbstractVector{T}, C::AbstractVector{CT}, m::Int, keepA::Bool=false, forceInf::Int=0) where {T<:Real, CT<:Complex{T}}
    if minimum(abs.(B)) < 1e-9
        @warn "There is at least one pole candidate in B being nearby zero. Consider shifting your problem for stability."
    end

    a = [A[1]]
    b = [forceInf > 0 ? T(Inf) : B[1]]
    β = [T(1.)]

    o = one(eltype(A))
    sA = fill(o, size(A))
    sB = fill(o, size(B))
    sC = fill(o, size(C))

    @inbounds for j = 1:m-1
        binv = 1 / b[j]
        βinv = 1 / β[j]
        for k = 1:length(A); sA[k] = sA[k] * βinv * (A[k] - a[j]) / (1 - A[k] * binv); end
        for k = 1:length(B); sB[k] = sB[k] * βinv * (B[k] - a[j]) / (1 - B[k] * binv); end
        for k = 1:length(C); sC[k] = sC[k] * βinv * (C[k] - a[j]) / (1 - C[k] * binv); end

        push!(a, A[keepA ? j+1 : argmax([isnan(x) ? -Inf : abs.(x) for x in sA])])
        push!(b, forceInf > j ? Inf : B[argmin([isnan(x) ? Inf : abs.(x) for x in sB])])
        push!(β, maximum(abs.(sC)))

        # treat single point case
        if β[j+1] < eps()
            β[j+1] = 1
        end
    end

    return a, b, β
end

"""
Compute scalar generalized divided differences.

# Arguments
- `σ`: Discretization of target set.
- `ξ`: Discretization of singularity set.
- `β`: Scaling factors.
"""
function scgendivdiffs(σ::AbstractVector{CT}, ξ, β, maxdgr, isfunm, pff) where CT<:Complex{<:Real}
    sgdd = zeros(CT, length(pff), maxdgr+2)
    for ii = 1:length(pff)
        if isfunm
            sgdd[ii,:] = ratnewtoncoeffsm(pff[ii], σ, ξ, β)
        else
            sgdd[ii,:] = map(m -> m[1], ratnewtoncoeffs(pff[ii], σ, ξ, β))
        end
    end
    return sgdd
end

"""
Compute rational divided differences for the function fun (can be matrix
valued), using differencing. The σ need to be distinct. For scalar
functions or non-distinct σ it may be better to use `ratnewtoncoeffsm`.
"""
function ratnewtoncoeffs(fun, σ::AbstractVector{CT}, ξ::AbstractVector{T}, β::AbstractVector{T}) where {T<:Real, CT<:Complex{T}}
    m = length(σ)
    D = Vector{Matrix{CT}}(undef, m)

    # compute divided differences D0,D1,...,Dm
    as_matrix(x::Number) = (M = Matrix{eltype(x)}(undef,1,1); M[1] = x; M)
    D[1] = fun(as_matrix(σ[1])) * β[1]

    for j = 2:m
        # evaluate current linearizaion at σ[j]
        Qj = zeros(T, size(D[1]))
        for k = 1:j-1
            Qj += D[k] * evalrat(σ[1:k-1], ξ[1:k-1], β[1:k], [σ[j]])[1]
        end

        # get divided difference from recursion (could be done via Horner)
        D[j] = (fun(as_matrix(σ[j])) .- Qj) / evalrat(σ[1:j-1], ξ[1:j-1], β[1:j], [σ[j]])[1]
    end

    return D
end

"""
Compute rational divided differences for the scalar function fun, using matrix
functions. `fm` has to be a function object representing a matrix function.
"""
function ratnewtoncoeffsm(fm, σ::AbstractVector{CT}, ξ::AbstractVector{T}, β::AbstractVector{T}) where {T<:Real, CT<:Complex{T}}
    m = length(σ) - 1

    σ = σ[:]
    ξ = ξ[:]
    β = β[:]

    # build Hessenberg matrices
    K = Bidiagonal(ones(m+1), β[2:m+1]./ξ[1:m], :L)
    H = Bidiagonal(σ[1:m+1], CT.(β[2:m+1]), :L)

    # column balancing
    P = Diagonal(1 ./ maximum(abs.(K), dims = 1)[:])
    K *= P
    H *= P
    H=Matrix(H) # this matrix is sparse, it is needed to convert before matrix-fun evaluation
    D = fm(H/CT.(K))[:,1] * β[1]
    D = copy(transpose(D))

    return D
end

"Evaluate nodal rational function at the points z."
function evalrat(σ::AbstractVector{CT}, ξ::AbstractVector{T}, β::AbstractVector{T}, z::AbstractVector{CT}) where {T<:Real, CT<:Complex{T}}
    r = ones(CT, size(z)) / β[1]
    for j = 1:length(σ)
        r .*= (z .- σ[j]) ./ (1 .- z/ξ[j]) / β[j+1]
    end
    return r
end