text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test
using Compat

# Wrap the functions defined in C++
wrap_modules("build/cxxwrap/libparametric")

import ParametricTypes.TemplateType, ParametricTypes.NonTypeParam

p1 = TemplateType{ParametricTypes.P1, ParametricTypes.P2}()
p2 = TemplateType{ParametricTypes.P2, ParametricTypes.P1}()

println("Dumping object p1:")
dump(p1)

@testset "Parametric tests" begin
    @test ParametricTypes.get_first(p1) == 1
    @test ParametricTypes.get_second(p2) == 1
    @test typeof(ParametricTypes.get_first(p1)) == Int32
    @test typeof(ParametricTypes.get_second(p2)) == Int32

    @test ParametricTypes.get_first(p2) == 10.
    @test ParametricTypes.get_second(p1) == 10.
    @test typeof(ParametricTypes.get_first(p2)) == Float64
    @test typeof(ParametricTypes.get_second(p1)) == Float64

    @test ParametricTypes.TemplateDefaultType{ParametricTypes.P1}() != nothing

    nontype1 = ParametricTypes.NonTypeParam{Int32, Int32(1)}()
    @test ParametricTypes.get_nontype(nontype1) == 1

    nontype2 = ParametricTypes.NonTypeParam{UInt32, UInt32(2)}()
    @test ParametricTypes.get_nontype(nontype2) == UInt32(2)

    nontype3 = ParametricTypes.NonTypeParam{Int32, Int32(1)}(3)
    @test ParametricTypes.get_nontype(nontype3) == 3

    nontype4 = ParametricTypes.NonTypeParam{Int64, Int64(64)}()
    @test ParametricTypes.get_nontype(nontype4) == Int64(64)

    concr = ParametricTypes.ConcreteTemplate{Float64}()
    @test isa(concr, ParametricTypes.AbstractTemplate{Float64})
    @test isa(concr, ParametricTypes.AbstractTemplate)
    @test isa(concr, ParametricTypes.ConcreteTemplate)
    abst = ParametricTypes.to_base(concr)
    @test isa(abst, ParametricTypes.AbstractTemplate{Float64})
    @test isa(abst, ParametricTypes.AbstractTemplate)

    f3 = ParametricTypes.Foo3{Int32, Bool, Float32}()
    @test length(methods(ParametricTypes.foo3_method)) == 6
    f2 = ParametricTypes.Foo2{Float64}()
    @test length(methods(ParametricTypes.foo2_method)) == 2

    @test length(methods(ParametricTypes.foo3_free_method)) == 6
    ParametricTypes.foo3_free_method(f3)

    @test supertype(ParametricTypes.Foo3{Float64,ParametricTypes.P1,Float32}) == ParametricTypes.AbstractTemplate{Float64}

    darr = [1.0, 2.0, 3.0]
    carr = Complex{Float32}[1+2im, 3+4im]
    vec1 = ParametricTypes.CppVector{Float64}(pointer(darr), 3)
    vec2 = ParametricTypes.CppVector2{Float64, Float32}()
    vec3 = ParametricTypes.CppVector{Complex{Float32}}(pointer(carr), 2)
    @test isa(vec1, AbstractVector{Float64})
    @test isa(vec2, AbstractVector{Float64})
    @test isa(vec3, AbstractVector{Complex{Float32}})
    @test ParametricTypes.get(vec1,0) == 1.0
    @test ParametricTypes.get(vec1,1) == 2.0
    @test ParametricTypes.get(vec1,2) == 3.0
    @test ParametricTypes.get(vec3,0) == 1+2im
    @test ParametricTypes.get(vec3,1) == 3+4im
end


text:
function inv(a::RelSeriesElem{<:Nemo.FieldElem})
#function inv(a::RelSeriesElem{nf_elem})
  @assert valuation(a)==0
  # x -> x*(2-xa) is the lifting recursion
  x = parent(a)(inv(coeff(a, 0)))
  set_prec!(x, 1)
  p = precision(a)
  la = [p]
  while la[end]>1
    push!(la, div(la[end]+1, 2))
  end

  two = parent(a)(base_ring(a)(2))
  set_prec!(two, p)

  n = length(la)-1
  y = parent(a)()
  while n>0
    set_prec!(x, la[n])
    set_prec!(y, la[n])
#    y = mul!(y, a, x)
#    y = two-y #sub! is missing...
#    x = mul!(x, x, y)
    x = x*(two-x*a)
    n -=1
  end
  return x
end

function log(a::RelSeriesElem{<:Nemo.FieldElem})
  @assert valuation(a)==0
  return integral(derivative(a)*inv(a))
end

function exp(a::RelSeriesElem{<:Nemo.FieldElem})
  @assert valuation(a) >0
  R = base_ring(parent(a))
  x = parent(a)([R(1)], 1, 2, 0)
  p = precision(a)
  la = [p]
  while la[end]>1
    push!(la, div(la[end]+1, 2))
  end

  one = parent(a)([R(1)], 1, 2, 0)

  n = length(la)-1
  # x -> x*(1-log(a)+a) is the recursion
  while n>0
    set_prec!(x, la[n])
    set_prec!(one, la[n])
    x = x*(one - log(x) + a) # TODO: can be optimized...
    n -=1
  end
  return x
end

doc"""
    derivative(f::RelSeriesElem{T}) -> RelSeriesElem
> Return the derivative of the power series $f$.
"""
function derivative(f::RelSeriesElem{T}) where T
  g = parent(f)()
  set_prec!(g, precision(f)-1)
  v = valuation(f)
  if v==0
    for i=1:Nemo.pol_length(f)
      setcoeff!(g, i-1, (i+v)*Nemo.polcoeff(f, i))
    end
    Nemo.set_val!(g, 0)
  else
    for i=0:Nemo.pol_length(f)
      setcoeff!(g, i, (i+v)*Nemo.polcoeff(f, i))
    end
    Nemo.set_val!(g, v-1)
  end
  Nemo.renormalize!(g)
  return g
end

doc"""
    integral(f::RelSeriesElem{T}) -> RelSeriesElem
> Return the integral of the power series $f$.
"""
function Nemo.integral(f::RelSeriesElem{T}) where T
  g = parent(f)()
  set_prec!(g, precision(f)+1)
  v = valuation(f)
  for i=0:Nemo.pol_length(f)
    setcoeff!(g, i, divexact(Nemo.polcoeff(f, i), i+v+1))
  end
  Nemo.set_val!(g, v+1)
  Nemo.renormalize!(g)
  return g
end

ext: jl


text:
using LinearAlgebra # for the norm and dot function
using ForwardDiff # for automatic differentiation

#=

## Gradient decent with backtracking line search

Write a function `gradient_decent(f, x0, α, β)` which implements the gradient decent with backtracking line search.

Below you find some pseudo-code and a few template for the functions you could create.
(You can change name and arguments of the functions in your implementation of course.)

The pseudo algorithms are:

Gradient decent:
given a starting point x_0

x = x_0
repeat until stopping criterium
	Compute decent direction Δx = ∇f(x).
 	Line search: Choose a step size t > 0.
    Update: x = x + t * Δx

 	Stopping criterium: ‖ ∇f(x) ‖² < ε
end

Backtracking line search:
given x, f, Δx, α, β

t = 1
repeat until f(x + t Δx) < f(x) + α t dot( ∇f(x), Δx )
 	Udate: t = β t
end

Source: [https://web.stanford.edu/class/ee364a/lectures/unconstrained.pdf](https://web.stanford.edu/class/ee364a/lectures/unconstrained.pdf)

=#

function linesearch(f, Δx, x, α, β; t_min = 1e-5)
	# write your code here

	return NaN
end

function gradientdecent(f, x0, α, β, ε; max_steps = 1e4)
	@assert β < 1
	@assert α < 0.5
	# write your code here

	return NaN
end

# test case

begin
	f(x) = exp(x[1]+3x[2]-0.1) + exp(x[1]-3x[2]-0.1) + exp(-x[1]-0.1)
	α = 0.25
	β = 0.5
	ε = 1e-8
	x0 = [0.0, 0.0]

	gradientdecent(f, x0, α, β, ε)
end

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
using Base.Threads

#=
"""
	 index(a::Int64,b::Int64)
Summary
======
Triangular indexing determination.

Arguments
======
a = row index

b = column index
"""
=#
@inline function triangular_index(a::Int,b::Int)
  return (muladd(a,a,-a) >> 1) + b
end

@inline function triangular_index(a::Int)
  return muladd(a,a,-a) >> 1
end

@inline function decompose(input::Int)
  #return ceil(Int,(-1.0+√(1+8*input))/2.0)
  return Base.fptosi(
    Int,
    Base.ceil_llvm(
      0.5*(
        Base.Math.sqrt_llvm(
          Base.sitofp(Float64, muladd(8,input,1))
        ) - 1.0
      )
    )
  )
end

function compute_enuc(mol::Molecule)
  E_nuc = 0.0
  for iatom in 1:length(mol), jatom in 1:(iatom-1)
    ix = mol[iatom].atom_center[1]
    jx = mol[jatom].atom_center[1]

    iy = mol[iatom].atom_center[2]
    jy = mol[jatom].atom_center[2]

    iz = mol[iatom].atom_center[3]
    jz = mol[jatom].atom_center[3]

    distance = √((jx-ix)^2 + (jy-iy)^2 + (jz-iz)^2)

    E_nuc += mol[iatom].atom_id*mol[jatom].atom_id/distance
  end

  return E_nuc
end

function compute_overlap(S::Matrix{Float64}, basis::Basis,
  jeri_oei_engine)

  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    S_block_JERI = zeros(Float64,(abas*bbas,))
    JERI.compute_overlap_block(jeri_oei_engine, S_block_JERI, ash, bsh,
      length(S_block_JERI))
    axial_normalization_factor(S_block_JERI, basis[ash], basis[bsh])

    idx = 1
    for ibas in 0:abas-1, jbas in 0:bbas-1
      iorb = apos + ibas
      jorb = bpos + jbas

      S[max(iorb,jorb),min(iorb,jorb)] = S_block_JERI[idx]

      idx += 1
    end
  end

  for iorb in 1:basis.norb, jorb in 1:iorb
    if iorb != jorb
      S[min(iorb,jorb),max(iorb,jorb)] = S[max(iorb,jorb),min(iorb,jorb)]
    end
  end
end

function compute_ke(T::Matrix{Float64}, basis::Basis,
  jeri_oei_engine)

  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    T_block_JERI = zeros(Float64,(abas*bbas,))
    JERI.compute_kinetic_block(jeri_oei_engine, T_block_JERI, ash, bsh,
      length(T_block_JERI))
    axial_normalization_factor(T_block_JERI, basis[ash],
      basis[bsh])

    idx = 1
    for ibas in 0:abas-1, jbas in 0:bbas-1
      iorb = apos + ibas
      jorb = bpos + jbas

      T[max(iorb,jorb),min(iorb,jorb)] = T_block_JERI[idx]

      idx += 1
    end
  end

  for iorb in 1:basis.norb, jorb in 1:iorb
    if iorb != jorb
      T[min(iorb,jorb),max(iorb,jorb)] = T[max(iorb,jorb),min(iorb,jorb)]
    end
  end
end

function compute_nah(V::Matrix{Float64}, mol::Molecule,
  basis::Basis, jeri_oei_engine)

  #== define ncenter ==#
  #=
  ncenter::Int64 = length(mol)

  Z = Vector{Float64}([])
  x = Vector{Float64}([])
  y = Vector{Float64}([])
  z = Vector{Float64}([])

  for atom in mol
    push!(Z, convert(Float64,atom.atom_id))
    push!(x, atom.atom_center[1])
    push!(y, atom.atom_center[2])
    push!(z, atom.atom_center[3])
  end
  =#
  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    V_block_JERI = zeros(Float64,(abas*bbas,))
    JERI.compute_nuc_attr_block(jeri_oei_engine, V_block_JERI, ash, bsh,
      length(V_block_JERI))
    axial_normalization_factor(V_block_JERI, basis[ash],
      basis[bsh])

    idx = 1
    for ibas in 0:abas-1, jbas in 0:bbas-1
      iorb = apos + ibas
      jorb = bpos + jbas

      V[max(iorb,jorb),min(iorb,jorb)] = V_block_JERI[idx]

      idx += 1
    end
  end

  for iorb in 1:basis.norb, jorb in 1:iorb
    if iorb != jorb
      V[min(iorb,jorb),max(iorb,jorb)] = V[max(iorb,jorb),min(iorb,jorb)]
    end
  end
end

function sad_guess(mol::Molecule, basis::Basis)
  basis_symbol = basis.model

  sad_guess = zeros(Float64, (basis.norb, basis.norb))
  h5open(joinpath(@__DIR__, "../../../records/sadgss.h5"),"r") do sadgss
    anchor = 1
    for atom in mol
      atom_symbol = atom.symbol

      sadgss_buf = read(sadgss["$atom_symbol/$basis_symbol"])
      #println("$anchor, $atom")
      #display(sadgss_buf); println()

      sqrt_nbas_guess = trunc(Int,sqrt(length(sadgss_buf)))

      sadgss_idx = 1
      for i in anchor:(anchor+sqrt_nbas_guess-1)
        for j in anchor:(anchor+sqrt_nbas_guess-1)
          sad_guess[i,j] = sadgss_buf[sadgss_idx]
          sadgss_idx += 1
        end
      end
      anchor += sqrt_nbas_guess
    end
  end

  #display(sad_guess)
  return sad_guess
end

function compute_schwarz_bounds(schwarz_bounds::Matrix{Float64},
  basis::Basis, nsh::Int64)

  max_am = max_ang_mom(basis)
  eri_quartet_batch = Vector{Float64}(undef,eri_quartet_batch_size(max_am))
  jeri_schwarz_engine = JERI.TEIEngine(basis.basis_cxx, basis.shpdata_cxx)

  for ash in 1:nsh, bsh in 1:ash
    #fill!(eri_quartet_batch, 0.0)
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas
    abshp = triangular_index(ash, bsh)

    JERI.compute_eri_block(jeri_schwarz_engine, eri_quartet_batch,
      ash, bsh, ash, bsh, abshp, abshp, abas*bbas, abas*bbas)

    #= axial normalization =#
    axial_normalization_factor(eri_quartet_batch, basis[ash], basis[bsh],
      basis[ash], basis[bsh], abas, bbas, abas, bbas)

    #== compute schwarz bound ==#
    max_index = BLAS.iamax(abas*bbas*abas*bbas,
      eri_quartet_batch, 1)
    schwarz_bounds[ash, bsh] = sqrt(abs(eri_quartet_batch[max_index]))
  end

  for ash in 1:nsh, bsh in 1:ash
    if ash != bsh
      schwarz_bounds[min(ash,bsh),max(ash,bsh)] =
        schwarz_bounds[max(ash,bsh),min(ash,bsh)]
    end
  end
end

#=
"""
		get_oei_matrix(oei::Array{Float64,2})
Summary
======
Extract one-electron integrals from data file object. Kinetic energy integrals,
overlap integrals, and nuclear attraction integrals can all be extracted.

Arguments
======
oei = array of one-electron integrals to extract
"""
=#
function read_in_oei(oei::Vector{T}, nbf::Int) where T
	nbf2 = (nbf*(nbf+1)) >> 1

	oei_matrix = Matrix{Float64}(undef,(nbf,nbf))
	for ibf in 1:nbf2
    i = decompose(ibf)
    j = ibf - triangular_index(i)

		oei_matrix[i,j] = float(oei[ibf])
		oei_matrix[j,i] = oei_matrix[i,j]
	end

	return oei_matrix
end

function DIIS(F::Matrix{Float64}, e_array::Vector{Matrix{Float64}},
  F_array::Vector{Matrix{Float64}}, B_dim::Int64)

  B = Matrix{Float64}(undef,B_dim+1,B_dim+1)
  for i in 1:B_dim, j in 1:B_dim
    B[i,j] = LinearAlgebra.BLAS.dot(length(e_array[i]), e_array[i], 1,
      e_array[j], 1)

	  B[i,B_dim+1] = -1.0
	  B[B_dim+1,i] = -1.0
	  B[B_dim+1,B_dim+1] = 0.0
  end
  #DIIS_coeff::Vector{Float64} = [ fill(0.0,B_dim)..., -1.0 ]
  DIIS_coeff::Vector{Float64} = vcat(zeros(B_dim), [-1.0])

  #DIIS_coeff[:], B[:,:], ipiv = LinearAlgebra.LAPACK.gesv!(B, DIIS_coeff)
  DIIS_coeff[:], B[:,:], ipiv = LinearAlgebra.LAPACK.sysv!('U', B, DIIS_coeff)

  #fill!(F, zero(Float64))
  LinearAlgebra.BLAS.scal!(length(F), 0.0, F, 1)
  for index in 1:B_dim
    #F .+= DIIS_coeff[index] .* F_array[index]
    F .+= DIIS_coeff[index] .* F_array[index]
  end
end

function axial_normalization_factor(oei, ash, bsh)
  ama = ash.am
  amb = bsh.am

  na = ash.nbas
  nb = bsh.nbas

  ab = 0
  for asize::Int64 in 0:(na-1), bsize::Int64 in 0:(nb-1)
    ab += 1

    anorm = axial_norm_fact[asize+1,ama]
    bnorm = axial_norm_fact[bsize+1,amb]

    abnorm = anorm*bnorm
    oei[ab] *= abnorm
  end
end

@inline function axial_normalization_factor(eri_quartet_batch::Vector{Float64},
  μsh::JCModules.Shell, νsh::JCModules.Shell,
  λsh::JCModules.Shell, σsh::JCModules.Shell,
  nμ::Int, nν::Int, nλ::Int, nσ::Int)

  #eri_quartet_batch_simint = similar(eri_quartet_batch)
  #println(ish, ",", jsh, ",", ksh, ",", lsh)
  amμ = μsh.am
  amν = νsh.am
  amλ = λsh.am
  amσ = σsh.am

  if amμ < 3 && amν < 3 && amλ < 3 && amσ < 3
    return
  else
    μνλσ = 0
    for μsize::Int64 in 0:(nμ-1), νsize::Int64 in 0:(nν-1)
      μνλσ = nσ*nλ*νsize + nσ*nλ*nν*μsize

      μnorm = axial_norm_fact[μsize+1,amμ]
      νnorm = axial_norm_fact[νsize+1,amν]

      μνnorm = μnorm*νnorm

      for λsize::Int64 in 0:(nλ-1), σsize::Int64 in 0:(nσ-1)
        μνλσ += 1

        λnorm = axial_norm_fact[λsize+1,amλ]
        σnorm = axial_norm_fact[σsize+1,amσ]

        λσnorm = λnorm*σnorm

        eri_quartet_batch[μνλσ] *= μνnorm*λσnorm
      end
    end
  end
end

function eri_quartet_batch_size(max_am)
  return am_to_nbas_cart(max_am)^4
end


text:
using DataFrames
using CSV
using DelimitedFiles
using Statistics
using Glob
using LinearAlgebra
include("../src/StochasticBlockModel.jl")

all_datasets = glob("instances/S1/in/*.in")
output_path = "experiments/experiment-S1-heuristics.txt"

heuristic_methods = ["ls1","ls2","exact"]

time_limit = 20.0
N_trials = 50

df = DataFrame(instance_path = String[],
               dataset_name = String[],
               n = Int[],
               m = Int[],
               q = Int[],
               avg_degree = Float64[],
               method = String[],
               LB = Float64[],
               UB = Float64[],
               status = String[],
               solvetime = Float64[],
               trial = Int[],
               seed = Int[],
               strong_assort = Bool[],
               weak_assort = Bool[],
               z = Matrix{Int}[],
               w = Matrix{Float64}[]
)

# Write column names of DataFrame to output file
df_columns = string(join(string.(names(df)), "\t"),"\n")
open(output_path, "a+") do io
   write(io, df_columns)
end

for filepath in all_datasets
    dataset_path = split(filepath, "/")[end]
    println(filepath)

    # Load Dataset
    dataset = StochasticBlockModel.Dataset(filepath)

    for method in heuristic_methods
        println("method: $method")
        for trial in 1:N_trials
            println("Trial: $trial")
            # Sample a random seed
            seed = abs(rand(Int))
            println("seed: $seed")

            # Estimate SBM
            sbm, z, opt_results = StochasticBlockModel.expectationMaximization(dataset,
                                                                               method,
                                                                               time_limit=time_limit,
                                                                               seed=seed,
                                                                               z=nothing)

            # Average degree
            avg_degree = mean(dataset.k)

            # Strong assortativity
            strong_assort = (minimum(diag(sbm.w)) >= maximum([sbm.w[i,j]
                             for i=1:dataset.n_communities, j=1:dataset.n_communities if j>i]))

            # Weak assortativity
            weak_assort = all([sbm.w[i,i] >= maximum([sbm.w[i,j] for j=1:dataset.n_communities if i!=j])
                               for i=1:dataset.n_communities])

            # Append to data frame
            row = (filepath, dataset_path, dataset.n, dataset.m, dataset.n_communities, avg_degree, method,
                   opt_results.LB, opt_results.UB, string(opt_results.status), opt_results.solvetime,
                   trial, seed, strong_assort, weak_assort, z, sbm.w)
            push!(df, row)

            # Write new line to output file
            CSV.write(output_path, df[end:end,:], delim="\t", append=true)

            # Save prob matrix w and assignments z
            res_folder = joinpath(join(split(output_path,"/")[1:end-1], "/"), string("n",dataset.n))
            if ~isdir(res_folder)
                mkdir(res_folder)
            end

            # Save w
            w_savepath = replace(joinpath(res_folder, dataset_path),".in"=>"-t$trial-$method.w")
            writedlm(w_savepath, sbm.w, '\t')

            # Save z
            z_savepath = replace(joinpath(res_folder, dataset_path),".in"=>"-t$trial-$method.z")
            writedlm(z_savepath, z, '\t')
        end
    end
end

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Test

# =========================================================================================
# Test tess_pipeline(Image List File).
@testset "tess_pipeline(Image List File)" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local result   = tess_pipeline_text(pipeline)

    local pix1, io = mktemp(;cleanup=false); close(io)
    local pix2, io = mktemp(;cleanup=false); close(io)
    local pix3, Io = mktemp(;cleanup=false); close(io)
    local output, io = mktemp(;cleanup=false);

    pix_write(pix1, pix_with([Line("One 1")]), IFF_TIFF); println(io, pix1);
    pix_write(pix2, pix_with([Line("Two 2")]), IFF_TIFF); println(io, pix2);
    pix_write(pix3, pix_with([Line("Three 3")]), IFF_TIFF); println(io, pix3);

    close(io)

    @test tess_run_pipeline(pipeline, output) == true

    @test result[] == "One 1\n\fTwo 2\n\fThree 3\n"

    rm(pix1)
    rm(pix2)
    rm(pix3)
    rm(output)
end

# =========================================================================================
# Test tess_pipeline(TIFF File).
@testset "tess_pipeline(TIFF File)" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local result   = tess_pipeline_text(pipeline)
    local pix, io = mktemp(;cleanup=false); close(io)

    pix_write_tiff(pix, pix_with([Line("One 1")]); append = true)
    pix_write_tiff(pix, pix_with([Line("Two 2")]); append = true)
    pix_write_tiff(pix, pix_with([Line("Three 3")]); append = true)

    @test tess_run_pipeline(pipeline, pix) == true

    @test result[] == "One 1\n\fTwo 2\n\fThree 3\n"

    rm(pix)
end

# =========================================================================================
# Test tess_run_pipeline(IImage List File) /w no renderers.
@testset "tess_run_pipeline(Image List File) /w no renderers" begin
    local inst       = TessInst("eng", datadir)
    local pipeline   = TessPipeline(inst)
    local err        = "No renderers attached to the pipeline."
    local pix1, io = mktemp(;cleanup=false); close(io)
    local pix2, io = mktemp(;cleanup=false); close(io)
    local pix3, Io = mktemp(;cleanup=false); close(io)
    local output, io = mktemp(;cleanup=false);

    pix_write(pix1, pix_with(), IFF_TIFF); println(io, pix1);
    pix_write(pix2, pix_with(), IFF_TIFF); println(io, pix2);
    pix_write(pix3, pix_with(), IFF_TIFF); println(io, pix3);

    close(io)

    @test (@test_logs (:error, err) tess_run_pipeline(pipeline, output)) == false
    rm(pix1)
    rm(pix2)
    rm(pix3)
    rm(output)
end

# =========================================================================================
# Test tess_pipeline(Image List File) delete instance before run.
@testset "tess_pipeline(Image List File) delete instance before run" begin
    local inst     = TessInst("eng", datadir)
    local pipeline = TessPipeline(inst)
    local err      = "Instance has been freed."

    @test tess_pipeline_text(pipeline) != nothing

    local pix1, io = mktemp(;cleanup=false); close(io)
    local pix2, io = mktemp(;cleanup=false); close(io)
    local pix3, Io = mktemp(;cleanup=false); close(io)
    local output, io = mktemp(;cleanup=false);

    pix_write(pix1, pix_with(), IFF_TIFF); println(io, pix1);
    pix_write(pix2, pix_with(), IFF_TIFF); println(io, pix2);
    pix_write(pix3, pix_with(), IFF_TIFF); println(io, pix3);

    close(io)

    tess_delete!(inst)

    @test (@test_logs (:error, err) tess_run_pipeline(pipeline, output)) == false

    rm(pix1)
    rm(pix2)
    rm(pix3)
    rm(output)
end


text:
"""
    _check(csv::String) -> Nothing

Check the structural integrity of the CSV file (see the exceptions list).

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Throws
- [`IncorrectHeader`](@ref): $(INCORRECT_HEADER_EXCEPTION)
- [`LastElementIsNotAnInteger`](@ref): $(LAST_ELEMENT_IS_NOT_AN_INTEGER_EXCEPTION)
- [`NotFiveColumns`](@ref): $(NOT_FIVE_COLUMNS_EXCEPTION)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, io = mktemp()
println(io, HEADER)
flush(io)
file = _file * ".csv"
mv(_file, file)

PDFHighlights.Internal.CSV._check(file)

# output

```
"""
function _check(csv::String)::Nothing

    lines = filter(!isempty, readlines(csv))

    # Check if the header is correct
    rstrip(lines[1]) != HEADER && throw(IncorrectHeader(csv))

    for (index, line) in enumerate(lines[2:end])

        # Check if the last element is parsing as an integer

        j = lastindex(line)
        i = findprev(',', line, j)

        if i != j
            try
                parse(Int32, line[(i+1):end])
            catch
                throw(LastElementIsNotAnInteger(csv, index + 1))
            end
        end

        # Check that each row represents values of 5 columns

        inside_quotes = false
        commas = 0

        for (char_index, char) in pairs(line)
            if char == ','
                if !inside_quotes
                    commas += 1
                end
            else
                if char == '"'
                    inside_quotes ? inside_quotes = false : inside_quotes = true
                end
            end
        end

        commas != 4 && throw(NotFiveColumns(csv, index + 1))

    end

    return nothing

end


text:
module ThreadSafeAccumulatorsRuntests

using Test
include("../examples/ThreadSafeAccumulators.jl")
using .ThreadSafeAccumulators

@test collect(Accumulator()) == Any[]
@test collect(Accumulator{Int}()) == Int[]
@test collect(push!(Accumulator{Int}(), 2)) == [2]
@test collect(push!(push!(Accumulator{Int}(), 1,2),3,4)) == 1:4
# constructors / convert
@test collect(Accumulator([])) == collect(Accumulator())
@test collect(Accumulator([1,2,3])) == 1:3
@test collect(Accumulator{Float64}([1])) == Float64[1]

@test eltype(Accumulator{Int}()) == Int
@test eltype(Accumulator([1,2,3])) == Int

end


text:
# * Scalar operators

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        D::QuasiDiagonal,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))
        Ai,Bi = indices(A,2), indices(B,2)
        if Ai == Bi
            Diagonal(Vector{T}(undef, length(Ai)))
        else
            m,n = length(Ai),length(Bi)
            offset = Ai[1]-Bi[1]
            BandedMatrix{T}(undef, (m,n), (-offset,offset))
        end
    end
    dest::Diagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))
        for (i,d) in enumerate(locs(B))
            dest.diag[i] = D.diag[d]
        end
    end
    dest::BandedMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))
        for (i,d) in enumerate(locs(B))
            dest.data[i] = D.diag[d]
        end
    end
end


text:
module AVFormat
  include(joinpath(Pkg.dir("VideoIO"),"src","init.jl"))
  w(f) = joinpath(avformat_dir, f)

  using AVUtil
  using AVCodecs

  include(w("LIBAVFORMAT.jl"))
end


text:
using NamedColors
using Test
using ColorTypes
using ColorTypes: N0f8

include("loading.jl")
include("namelookup.jl")


text:
# Author: Lukasz Bratos

include("blocksys.jl")

using Test
using LinearAlgebra

sizes = [16, 10000, 50000]

@testset "Tests for matrix $size" for size in sizes
    A, n, l = Blocksys.readMatrix("./examples/$(size)/A.txt")
    b = Blocksys.readVector("./examples/$(size)/b.txt")

    @testset "Gauss elimination" begin
        @test isapprox(Blocksys.solveGauss(A, n, l, b), \(A, b))
    end

    @testset "Gauss elimination with choice" begin
        @test isapprox(Blocksys.solveGaussWithChoice(A, n, l, b), \(A, b))
    end

    @testset "Solve LU" begin
        L, U = Blocksys.lu(A, n, l)
        @test isapprox(Blocksys.solveLu(L, U, n, l, b), \(A, b))
    end

    @testset "Solve LU with choice" begin
        L, U, perm = Blocksys.luWithChoice(A, n, l)
        @test isapprox(Blocksys.solveLuWithChoice(L, U, n, l, b, perm), \(A, b))
    end
end


text:
using ITensors
using ITensorInfiniteMPS
using Random

Random.seed!(1234)

#
# Example
#

N = 10
s = siteinds("S=1/2", N; conserve_szparity=true)
χ = 6
@assert iseven(χ)
space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
ψ = InfiniteMPS(ComplexF64, s; space=space)
randn!.(ψ)

ψ = orthogonalize(ψ, :)
@show norm(prod(ψ.AL[1:N]) * ψ.C[N] - ψ.C[0] * prod(ψ.AR[1:N]))


text:
"""
    AbstractInjectLocation

An abstract type that represents all inject locations. The type hierarchy
is as follows:

```
AbstractInjectLocation
   |_ LevelLocation
   |_ KwargsLocation
   |_ MessageLocation
         |_ BeginningMessageLocation
         |_ EndingMessageLocation
         |_ WholeMessageLocation
```
"""
abstract type AbstractInjectLocation end

struct LevelLocation <: AbstractInjectLocation end

abstract type MessageLocation <: AbstractInjectLocation end
struct BeginningMessageLocation <: MessageLocation end
struct EndingMessageLocation <: MessageLocation end
struct WholeMessageLocation <: MessageLocation end

struct KwargsLocation <: AbstractInjectLocation end

"""
    AbstractLogProperty

An abstract type that represents all log properties.

```
AbstractLogProperty
   |_ LevelProperty
   |_ MessageProperty
   |_ KwargsProperty
```
"""
abstract type AbstractLogProperty end
struct LevelProperty <: AbstractLogProperty end
struct MessageProperty <: AbstractLogProperty end
struct KwargsProperty <: AbstractLogProperty end

struct ColorSpec
    color::Union{Symbol,Int}
    bold::Bool
end


text:
include("CreateMeshFCFV.jl")
using LoopVectorization
using CairoMakie
import AbstractPlotting.GeometryBasics
using SparseArrays, LinearAlgebra
import UnicodePlots
import Plots

# # ENV["MPLBACKEND"]="Qt5Agg"
# # ENV["MPLBACKEND"] = "module://gr.matplotlib.backend_gr"
# # Problem #1 was to find a simple way to plot piecewise constant fields on triangles:
# # solution from: https://robertsweeneyblanco.github.io/Programming_for_Mathematical_Applications/Computational_Geometry/Triangulations.html
# # Problem #2 was to get any figures popping out in VScode
# # and: https://github.com/JuliaPy/PyPlot.jl/issues/418
# # import PyCall, PyPlot
# using Base.Threads
# using LoopVectorization
# import TriangleMesh
# using Printf
# # PyCall.pygui(:qt5)
# # PyPlot.pygui(true)

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

# function tplot(mesh, v)

#     # Prepare mesh for visalisation
#     p = fill(Float64[], mesh.nv)
#     for i=1:mesh.nv
#         p[i] = [mesh.xv[i], mesh.yv[i]]
#     end
#     t = fill(Int64[], mesh.nel)
#     for i=1:mesh.nel
#         t[i] = [mesh.e2v[i,1], mesh.e2v[i,2], mesh.e2v[i,3]]
#     end

#     # Plot triangular mesh with nodes `p` and triangles `t`
#     PyPlot.clf()
#     tris =  PyPlot.convert(Array{Int64}, hcat(t...)')
#     display(PyPlot.tripcolor(first.(p), last.(p), tris .- 1, v,
#               cmap="viridis", edgecolors="none", linewidth=0) )
#     PyPlot.axis("equal")
#     PyPlot.ylim([0, 1])
#     PyPlot.xlim([0, 1])
#     PyPlot.title("Low res.")
#     PyPlot.xlabel("x")
#     PyPlot.ylabel("y")
#     PyPlot.colorbar()
#     PyPlot.show()
#     return
# end

@views function PlotMakie( mesh, v )
    p   = [AbstractPlotting.GeometryBasics.Polygon( Point2f0[ (mesh.xv[mesh.e2v[i,j]], mesh.yv[mesh.e2v[i,j]]) for j=1:mesh.nf_el] ) for i in 1:mesh.nel]
    scene = poly(p, color = v, colormap = :jet1, strokewidth = 0, strokecolor = :black, markerstrokewidth=0, markerstrokecolor = (0, 0, 0, 0), aspect_ratio=:equal, clims=[0.7 1.4])
    display(scene)

    # red = 20
    # iel = 100
    # p = Plots.plot()
    # for iel=1:10#mesh.nel
    #         Plots.plot!(  [mesh.xv[mesh.e2v[iel,2]], mesh.xv[mesh.e2v[iel,3]]], [mesh.yv[mesh.e2v[iel,2]], mesh.yv[mesh.e2v[iel,3]]], leg = false, aspect_ratio=:equal   )
    #         Plots.plot!( [mesh.xv[mesh.e2v[iel,3]], mesh.xv[mesh.e2v[iel,1]]], [mesh.yv[mesh.e2v[iel,3]], mesh.yv[mesh.e2v[iel,1]]]  )
    #         Plots.plot!( [mesh.xv[mesh.e2v[iel,1]], mesh.xv[mesh.e2v[iel,2]]], [mesh.yv[mesh.e2v[iel,1]], mesh.yv[mesh.e2v[iel,2]]]  )

    #     Plots.scatter!( [mesh.xf[mesh.e2f[iel,1]]], [mesh.yf[mesh.e2f[iel,1]]]  )
    #     Plots.plot!(    [mesh.xf[mesh.e2f[iel,1]], mesh.xf[mesh.e2f[iel,1]]+mesh.n_x[iel,1]/red], [mesh.yf[mesh.e2f[iel,1]], mesh.yf[mesh.e2f[iel,1]]+mesh.n_y[iel,1]/red ] )
    #     Plots.scatter!( [mesh.xf[mesh.e2f[iel,2]]], [mesh.yf[mesh.e2f[iel,2]]]  )
    #     Plots.plot!(    [mesh.xf[mesh.e2f[iel,2]], mesh.xf[mesh.e2f[1iel,2]]+mesh.n_x[iel,2]/red], [mesh.yf[mesh.e2f[iel,2]], mesh.yf[mesh.e2f[iel,2]]+mesh.n_y[iel,2]/red ] )
    #     Plots.scatter!( [mesh.xf[mesh.e2f[iel,3]]], [mesh.yf[mesh.e2f[iel,3]]]  )
    #     Plots.plot!(    [mesh.xf[mesh.e2f[iel,3]], mesh.xf[mesh.e2f[iel,3]]+mesh.n_x[iel,3]/red], [mesh.yf[mesh.e2f[iel,3]], mesh.yf[mesh.e2f[iel,3]]+mesh.n_y[iel,3]/red ] )
    # end
    # display(p)
    return
end

function Tanalytic2!( mesh, T, Tdir, Tneu, se, a, b, c, d, alp, bet )
    # Evaluate T analytic on cell faces
    @avx for in=1:mesh.nf
        x        = mesh.xf[in]
        y        = mesh.yf[in]
        Tdir[in] = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
    end
    # Evaluate T analytic on barycentres
    @avx for iel=1:mesh.nel
        x       = mesh.xc[iel]
        y       = mesh.yc[iel]
        T       = exp(alp*sin(a*x + c*y) + bet*cos(b*x + d*y))
        T[iel]  = T
        se[iel] = T*(-a*alp*cos(a*x + c*y) + b*bet*sin(b*x + d*y))*(a*alp*cos(a*x + c*y) - b*bet*sin(b*x + d*y)) + T*(a^2*alp*sin(a*x + c*y) + b^2*bet*cos(b*x + d*y)) + T*(-alp*c*cos(a*x + c*y) + bet*d*sin(b*x + d*y))*(alp*c*cos(a*x + c*y) - bet*d*sin(b*x + d*y)) + T*(alp*c^2*sin(a*x + c*y) + bet*d^2*cos(b*x + d*y))
    end
    return
end

function StabParam(tau, dA, Vol)
    taui = tau*dA
    return taui
end

function ComputeFCFV(mesh, se, Tdir, tau)
    ae = zeros(mesh.nel)
    be = zeros(mesh.nel)
    ze = zeros(mesh.nel,2)

    # Assemble FCFV elements
    @avx for iel=1:mesh.nel

        be[iel] = be[iel]   + mesh.vole[iel]*se[iel]

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel])                              # Stabilisation parameter for the face

            # Assemble
            ze[iel,1] += (bc==1) * dAi*ni_x*Tdir[nodei]  # Dirichlet
            ze[iel,2] += (bc==1) * dAi*ni_y*Tdir[nodei]  # Dirichlet
            be[iel]   += (bc==1) * dAi*taui*Tdir[nodei]  # Dirichlet
            ae[iel]   +=           dAi*taui

        end
    end
    return ae, be, ze
end

function ComputeElementValues(mesh, uh, ae, be, ze, Tdir, tau)
    ue          = zeros(mesh.nel);
    qx          = zeros(mesh.nel);
    qy          = zeros(mesh.nel);

    @avx for iel=1:mesh.nel

        ue[iel]  =  be[iel]/ae[iel]
        qx[iel]  = -1.0/mesh.vole[iel]*ze[iel,1]
        qy[iel]  = -1.0/mesh.vole[iel]*ze[iel,2]

        for ifac=1:mesh.nf_el

            # Face
            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel])      # Stabilisation parameter for the face

            # Assemble
            ue[iel] += (bc!=1) *  dAi*taui*uh[mesh.e2f[iel, ifac]]/ae[iel]
            qx[iel] -= (bc!=1) *  1.0/mesh.vole[iel]*dAi*ni_x*uh[mesh.e2f[iel, ifac]]
            qy[iel] -= (bc!=1) *  1.0/mesh.vole[iel]*dAi*ni_y*uh[mesh.e2f[iel, ifac]]
         end
    end
    return ue, qx, qy
end

@views function main()

    # Create sides of mesh
    xmin, xmax = 0, 1
    ymin, ymax = 0, 1
    nx, ny     = 2, 2
    # mesh_type  = "Quadrangles"
    mesh_type  = "UnstructTriangles"

    if mesh_type=="Quadrangles"
        tau = 100
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax )
    elseif mesh_type=="UnstructTriangles"
        tau = 100
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax )
    end

    println("Number of elements: ", mesh.nel)

    # Source term and BCs etc...
    Tanal  = zeros(mesh.nel)
    se     = zeros(mesh.nel)
    Tdir   = zeros(mesh.nf)
    Tneu   = zeros(mesh.nf)
    alp = 0.1; bet = 0.3; a = 5.1; b = 4.3; c = -6.2; d = 3.4;
    @time Tanalytic2!(mesh , Tanal, Tdir, Tneu, se, a, b, c, d, alp, bet)

    # Compute some mesh vectors
    ae, be, ze = ComputeFCFV(mesh, se, Tdir, tau)

    # Assemble stiffness matrix
    rows = Int64[]
    cols = Int64[]
    vals = Float64[]
    f    = zeros(mesh.nf);

    for iel=1:mesh.nel

        Ke = zeros(mesh.nf_el,mesh.nf_el);
        fe = zeros(mesh.nf_el,1);

        for ifac=1:mesh.nf_el

            nodei  = mesh.e2f[iel,ifac]
            bci    = mesh.bc[nodei]

            if bci != 1

                dAi  = mesh.dA[iel,ifac]
                ni_x = mesh.n_x[iel,ifac]
                ni_y = mesh.n_y[iel,ifac]
                taui = StabParam(tau,dAi,mesh.vole[iel])

                for jfac=1:mesh.nf_el

                    nodej  = mesh.e2f[iel,jfac]
                    bcj    = mesh.bc[nodej]

                    if bcj!= 1

                        dAj  = mesh.dA[iel,jfac]
                        nj_x = mesh.n_x[iel,jfac]
                        nj_y = mesh.n_y[iel,jfac]
                        tauj = StabParam(tau,dAj,mesh.vole[iel])

                        # Delta
                        del = 0.0
                        if ifac==jfac; del = 1.0; end

                        # Element matrix
                        nitnj         = ni_x*nj_x + ni_y*nj_y;
                        Ke[ifac,jfac] = dAi * (1.0/ae[iel] * dAj * taui*tauj - 1.0/mesh.vole[iel]*dAj*nitnj - taui*del);
                    end

                end

                # RHS vector
                Xi = 0.0;
                if bci == 2; Xi = 1.0; end # indicates Neumann dof
                ti = Tneu[nodei]
                nitze     = ni_x*ze[iel,1] + ni_y*ze[iel,2]
                fe[ifac]  = dAi * (1.0/mesh.vole[iel]*nitze - ti*Xi - 1.0/ae[iel]*be[iel]*taui)

            end
        end
        # println(Ke)

        for ifac=1:mesh.nf_el
            nodei  = mesh.e2f[iel,ifac]
            bci    = mesh.bc[nodei]

            if bci != 1

                for jfac=1:mesh.nf_el

                    nodej  = mesh.e2f[iel,jfac]
                    bcj    = mesh.bc[nodej]

                    if bcj != 1
                        push!(rows, mesh.e2f[ iel,ifac])
                        push!(cols, mesh.e2f[ iel,jfac])
                        push!(vals,      -Ke[ifac,jfac])
                    end
                end
                f[mesh.e2f[ iel,ifac]] -= fe[ifac]
            else
                push!(rows, mesh.e2f[ iel,ifac])
                push!(cols, mesh.e2f[ iel,ifac])
                push!(vals,                 1.0)
                f[mesh.e2f[ iel,ifac]] = Tdir[nodei]
            end
        end

    end
    K = sparse(rows, cols, vals, mesh.nf, mesh.nf)
    droptol!(K, 1e-6)
    uh   = K\f

    println(uh)

    # Reconstruct element values
    Te, qx, qy = ComputeElementValues(mesh, uh, ae, be, ze, Tdir, tau)

    # display(UnicodePlots.spy(K))
    # println(uh[:])

    # for iel=1:mesh.nel
    #     println(ue[iel])
    # end

    # for i=1:length(uh)
    #     println(uh[i])
    # end

    # println(mesh.e2f)

    # import MAT
    # file   = MAT.matopen("/Users/imac/ownCloud/FCFV/Mat100_ue.mat")
    # ue_mat = MAT.read(file, "ue")

    # Visualise
    @time PlotMakie( mesh, Te )
    println(maximum(Te))
    println(maximum(Tanal))
    # tplot(mesh, ue)
    # function qplot(x, y, v)
    # xc=LinRange(xmin,xmax,nx)
    # yc=LinRange(ymin,ymax,ny)
    # display( heatmap(xc,yc,reshape(ue,nx,ny), clim=[0.7 1.5]) )
        # clf()
        # display( pcolor(xc, yc, reshape(ue,nx,ny)) )
        # colorbar()
        # xlabel("x")
        # ylabel("y")
        # show()
    # end

end

main()


text:
module VCFTools

using CodecZlib
using Distributions
using HypothesisTests
using ProgressMeter
using DelimitedFiles
using Dates
using SpecialFunctions
using VariantCallFormat
import VariantCallFormat.findgenokey

export conformgt_by_id, conformgt_by_pos,
    gtstats, geno_ismissing,
    nrecords, nsamples, openvcf,
    sampleID,
    grm,
    # convert functions
    save_snpinfo,
    convert_gt, copy_gt!,
    convert_ht, copy_ht!,
    convert_ds, copy_ds!,
    copy_gt_trans!, copy_ht_trans!, copy_ds_trans!,
    # filter functions
    filter_genotype,
    filter, filter_header,
    filter_chr, filter_range,
    mask_gt, find_duplicate_marker,
    # aim selection
    aim_select,
    # write routines
    write_vcf

# package code goes here
include("gtstats.jl")
include("conformgt.jl")
include("convert.jl")
include("filter.jl")
include("aim_select.jl")
include("grm.jl")
include("write.jl")

# test data directory
datadir(parts...) = joinpath(@__DIR__, "..", "test", parts...)

end # module


text:
#!/usr/bin/env julia
using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3, Twist
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport gazebo_msgs.msg: ModelState, ModelStates
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties, GetLinkState, SetLinkState
@rosimport obstacle_detector.msg: Obstacles, CircleObstacle
@rosimport mavs_control.msg: Control

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.srv
using gazebo_msgs.msg
using obstacle_detector.msg
using mavs_control.msg

using NLOptControl
using VehicleModels
using DataFrames
using MAVs
using KNITRO
using PyCall
@pyimport tf.transformations as tf

# TODO
# 1) detect crash
# 3) get the Gazebo sim time
# 5) pass these parameters to a service or something -> needs to happen on a seperate node!
# 6) update OCP based off of current gazebo model position
# 7) get the current position of Gazebo model => need all of the states!!
# 8) get time from the Gazebo world
# 9) populate obstacle field in Gazebo based off of c.o

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/6/2017, Last Modified: 8/2/2017 \n
--------------------------------------------------------------------------------------\n
"""
function setObstacleData(params)

  if RobotOS.has_param("obstacle_radius")
    Q = params[2][7];           # number of obstacles the algorithm can handle

    if isnan(RobotOS.get_param("obstacle_radius"))
      for i in 1:Q
        setvalue(params[2][1][i],1.0);
        setvalue(params[2][2][i],1.0);
        setvalue(params[2][3][i],-100.0 - 3*i);
        setvalue(params[2][4][i],-100.0);
        setvalue(params[2][5][i],0.0);
        setvalue(params[2][6][i],0.0);
      end
    else
      r=deepcopy(RobotOS.get_param("obstacle_radius"))
      x=deepcopy(RobotOS.get_param("obstacle_x"))
      y=deepcopy(RobotOS.get_param("obstacle_y"))
      vx=deepcopy(RobotOS.get_param("obstacle_vx"))
      vy=deepcopy(RobotOS.get_param("obstacle_vy"))

      L = length(r)               # number of obstacles detected
      N = Q - L;
      if N < 0
        warn(" \n The number of obstacles detected exceeds the number of obstacles the algorithm was designed for! \n
                  Consider increasing the number of obstacles the algorithm can handle \n!")
      end

      for i in 1:Q
        if i <= L          # add obstacle
          setvalue(params[2][1][i],r[i]);
          setvalue(params[2][2][i],r[i]);
          setvalue(params[2][3][i],x[i]);
          setvalue(params[2][4][i],y[i]);
          setvalue(params[2][5][i],vx[i]);
          setvalue(params[2][6][i],vy[i]);
        else              # set obstacle off field
          setvalue(params[2][1][i],1.0);
          setvalue(params[2][2][i],1.0);
          setvalue(params[2][3][i],-100.0 - 3*i);
          setvalue(params[2][4][i],-100.0);
          setvalue(params[2][5][i],0.0);
          setvalue(params[2][6][i],0.0);
        end
      end
    end
  end

  return nothing
end

function main()
  init_node("rosjl_obstacles")
  #linkName = "base_footprint"  # TODO make this a parameter

  modelName = RobotOS.get_param("robotName")

  ###############################
  # set up services and messages

  # pause simulation
  const pause_physics = ServiceProxy("/gazebo/pause_physics",Empty)
  println("Waiting for '/gazebo/pause_physics' service...")
  wait_for_service("/gazebo/pause_physics")

  # unpause simulation
  const unpause_physics = ServiceProxy("/gazebo/unpause_physics",Empty)
  println("Waiting for '/gazebo/unpause_physics' service...")
  wait_for_service("/gazebo/unpause_physics")

  # Set up service to get Gazebo model state
  const get_state = ServiceProxy("/gazebo/get_model_state", GetModelState)
  println("Waiting for 'gazebo/get_model_state' service...")
  wait_for_service("gazebo/get_model_state")

  # Set up service to set Gazebo model state
  const set_state = ServiceProxy("/gazebo/set_model_state", SetModelState)
  println("Waiting for 'gazebo/set_model_state' service...")
  wait_for_service("gazebo/set_model_state")

  # message for solution to optimal control problem
  const pub = Publisher{Control}("/mavs/optimal_control", queue_size=10)

  # set up services and messages
  ###############################

  c=defineCase(;(:mode=>:autoGazebo));

  # Get the current position of the Gazebo model
  gs = GetModelStateRequest()
  gs.model_name = modelName
  gs_r = get_state(gs)

  if !gs_r.success
      error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
  end

  # Define position to move robot
  vehPose = gs_r.pose  # use the current position
  vehPose.position.x = c.m.X0[1]
  vehPose.position.y = c.m.X0[2]
  roll = 0; pitch = 0; yaw = c.m.X0[5];
  Q = tf.quaternion_from_euler(roll, pitch, yaw)
  vehPose.orientation.x = Q[1]
  vehPose.orientation.y = Q[2]
  vehPose.orientation.z = Q[3]
  vehPose.orientation.w = Q[4]

  # Define the robot state
  ms = ModelState()
  ms.model_name = modelName
  ms.pose = vehPose

  # Set the state of the Gazebo model
  ss = SetModelStateRequest()
  ss.model_state = ms
  println("Calling 'gazebo/set_model_state' service...")
  ss_r = set_state(ss)

  n=initializeAutonomousControl(c);

  driveStraight!(n)

  loop_rate = Rate(2.0)

  # unpause physics
  up = EmptyRequest()
  unpause_physics(up)

  while !is_shutdown()

    if ((n.r.dfs_plant[end][:x][end]-c.g.x_ref)^2 + (n.r.dfs_plant[end][:y][end]-c.g.y_ref)^2)^0.5 < 2*n.XF_tol[1]
       println("Goal Attained! \n"); n.mpc.goal_reached=true;
       pp = EmptyRequest()
       pause_physics(pp)
       break;
    end

    println("Running model for the: ", n.r.eval_num," time");
    setObstacleData(n.params)
    updateAutoParams!(n,c)  # check to see if goal is in range

    # get the cur6rent position of Gazebo model
    gs = GetModelStateRequest()
    gs.model_name = modelName
    gs_r = get_state(gs)
    if !gs_r.success
        error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
    end
    X0 = zeros(8)
    X0[1] = gs_r.pose.position.x
    X0[2] = gs_r.pose.position.y
    X0[3] = n.mpc.X0p[3]
    X0[4] = n.mpc.X0p[4]
    Q = [gs_r.pose.orientation.x,gs_r.pose.orientation.y,gs_r.pose.orientation.z,gs_r.pose.orientation.w]
    roll, pitch, yaw = tf.euler_from_quaternion(Q)
    X0[5] = yaw
    X0[6] = n.mpc.X0p[6]
    X0[7] = n.mpc.X0p[7]
    X0[8] = n.mpc.X0p[8]
    println(string(X0,"\n"))
    updateX0!(n,X0;(:userUpdate=>true))
    status=autonomousControl!(n)                # rerun optimization
  #  n.mpc.t0_actual=(n.r.eval_num-1)*n.mpc.tex  # external so that it can be updated easily in PathFollowing
    n.mpc.t0_actual = to_sec(get_rostime())
    msg = Control()
    msg.t = n.mpc.t0_actual + n.r.t_st
    msg.x = n.r.X[:,1]
    msg.y = n.r.X[:,2]
    msg.psi = n.r.X[:,5]
    msg.sa = n.r.X[:,6]
    msg.vx = n.r.X[:,7]
    # TODO consider buffering the message here..
    publish(pub,msg)

  #  spinOnce() # call the callbacks
    rossleep(loop_rate) # NOTE might not want to do this ...
  end
end

if !isinteractive()
    main()
end


text:
#' ---
#' title: Working with Data Files
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' In many cases we might need to read data available in an external file rather than type it into Julia ourselves.
#' This tutorial is concerned with reading tabular data into Julia and using it for a JuMP model.
#' We'll be reading data using the DataFrames.jl package and some other packages specific to file types.

#' Note: There are multiple ways to read the same kind of data intto Julia.
#' However, this tutorial only focuses on DataFrames.jl as
#' it provides the ecosystem to work with most of the required file types in a straightforward manner.

#' ### DataFrames.jl

#' The DataFrames package provides a set of tools for working with tabular data.
#' It is available through the Julia package system.
#+ tangle = false

using Pkg
Pkg.add("DataFrames")

#' ### What is a DataFrame?

#' A DataFrame is a data structure like a table or spreadsheet. You can use it for storing and exploring a set of related data values.
#' Think of it as a smarter array for holding tabular data.

#' ## Reading Tabular Data into a DataFrame
#' We will begin by reading data from different file formats into a DataFrame object.
#' The example files that we will be reading are present in the data folder.

#' ### Excel Sheets
#' Excel files can be read using the ExcelFiles.jl package.
#+ tangle = false

Pkg.add("XLSX")

#' To read a Excel file into a DataFrame, we use the following julia code.
#' The first arguement to the `readtable` function is the file to be read and the second arguement is the name of the sheet.

using DataFrames
using XLSX

#+

data_dir = joinpath(@__DIR__, "data")
excel_df = DataFrame(XLSX.readtable(joinpath(data_dir, "SalesData.xlsx"), "SalesOrders")...)

#' ### CSV Files
#' CSV and other delimited text files can be read the CSV.jl package.
#+ tangle = false

Pkg.add("CSV")

#' To read a CSV file into a DataFrame, we use the `CSV.read` function.

using CSV
csv_df = CSV.read(joinpath(data_dir, "StarWars.csv"))

#' ### Other Delimited Files
#' We can also use the CSV.jl package to read any other delimited text file format.
#' By default, CSV.File will try to detect a file's delimiter from the first 10 lines of the file.
#' Candidate delimiters include `','`, `'\t'`, `' '`, `'|'`, `';'`, and `':'`. If it can't auto-detect the delimiter, it will assume `','`.
#' Let's take the example of space separated data.

ss_df = CSV.read(joinpath(data_dir, "Cereal.txt"))

#' We can also specify the delimiter by passing the `delim` arguement.

delim_df = CSV.read(joinpath(data_dir, "Soccer.txt"), delim = "::")

#' Note that by default, are read-only. If we wish to make changes to the data read, we pass the `copycols = true` arguement in the function call.
#+ tangle = false

ss_df = CSV.read(joinpath(data_dir, "Cereal.txt"), copycols = true)

#' ## Working with DataFrames
#' Now that we have read the required data into a DataFrame, let us look at some basic operations we can perform on it.

#' ### Querying Basic Information
#' The `size` function gets us the dimensions of the DataFrame.

size(ss_df)

#' We can also us the `nrow` and `ncol` functions to get the number of rows and columns respectively.

nrow(ss_df), ncol(ss_df)

#' The `describe` function gives basic summary statistics of data in a DataFrame.

describe(ss_df)

#' Names of every column can be obtained by the `names` function.

names(ss_df)

#' Corresponding data types are obtained using the broadcasted `eltype` function.

eltype.(ss_df)

#' ### Accessing the Data
#' Similar to regular arrays, we use numerical indexing to access elements of a DataFrame.

csv_df[1,1]

#' The following are different ways to access a column.

csv_df[!, 1]

#+

csv_df[!, :Name]

#+

csv_df.Name

#+

csv_df[:, 1] # note that this creates a copy

#' The following are different ways to access a row.

csv_df[1:1, :]

#+

csv_df[1, :] # this produces a DataFrameRow

#' We can change the values just as we normally assign values.

#' Assign a range to scalar.

excel_df[1:3, 5] .= 1

#' Vector to equal length vector.

excel_df[4:6, 5] = [4, 5, 6]

#' Subset of the DataFrame to another data frame of matching size.

excel_df[1:2, 6:7] =  DataFrame([-2 -2; -2 -2], [Symbol("Unit Cost"), :Total])

#+

excel_df

#' There are a lot more things which can be done with a DataFrame.
#' See the [docs](https://juliadata.github.io/DataFrames.jl/stable/) for more information.

#' ## A Complete Modelling Example - Passport Problem

#' Let's now apply what we have learnt to solve a real modelling problem.

#' The [Passport Index Dataset](https://github.com/ilyankou/passport-index-dataset)
#' lists travel visa requirements for 199 countries, in .csv format.
#' Our task is to find out the minimum number of passports required to visit all countries.

#' In this dataset, the first column represents a passport (=from) and each remaining column represents a foreign country (=to).
#' The values in each cell are as follows:
#' * 3 = visa-free travel
#' * 2 = eTA is required
#' * 1 = visa can be obtained on arrival
#' * 0 = visa is required
#' * -1 is for all instances where passport and destination are the same

#' Our task is to find out the minimum number of passports needed to visit every country without requiring a visa.
#' Thus, the values we are interested in are -1 and 3. Let us modify the data in the following manner -

passportdata = CSV.read(joinpath(data_dir, "passport-index-matrix.csv"), copycols = true)

for i in 1:nrow(passportdata)
    for j in 2:ncol(passportdata)
        if passportdata[i,j] == -1 || passportdata[i,j] == 3
            passportdata[i,j] = 1
        else
            passportdata[i,j] = 0
        end
    end
end

#' The values in the cells now represent:
#' * 1 = no visa required for travel
#' * 0 = visa required for travel

#' Let us assossciate each passport with a decision variable $pass_{cntr}$ for each country.
#' We want to minize the sum $\sum pass_{cntr}$ over all countries.

#' Since we wish to visit all the countries, for every country,
#' we should own atleast one passport that lets us travel to that country visa free.
#' For one destination, this can be mathematically represented as $\sum_{cntr \in world} passportdata_{cntr,dest} \cdot pass_{cntr} \geq 1$.

#' Thus, we can represent this problem using the following model:

#' $$
#' \begin{align*}
#' \min && \sum_{cntr \in World} pass_{cntr} \\
#' \text{s.t.} && \sum_{cntr \in World} passportdata_{cntr,dest} \cdot pass_{cntr} \geq 1 && \forall dest \in World \\
#' && pass_{cntr} \in \{0,1\} && \forall cntr \in World
#' \end{align*}
#' $$

#' We'll now solve the problem using JuMP.

using JuMP, GLPK

# Finding number of countries
n = ncol(passportdata) - 1 # Subtract 1 for column representing country of passport

model = Model(GLPK.Optimizer)
@variable(model, pass[1:n], Bin)
@constraint(model, [j = 2:n], sum(passportdata[i,j] * pass[i] for i in 1:n) >= 1)
@objective(model, Min, sum(pass))
optimize!(model)

println("Minimum number of passports needed: ", objective_value(model))

#+

countryindex = findall(value.(pass) .== 1 )

print("Countries: ")
for i in countryindex
    print(names(passportdata)[i+1], " ")
end


text:
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# FILE CONTENTS:
#     SUMMARY: Integration schemes
#
#     INTEGRATION METHODS
#         midpoint: midpoint or trapezoidal integration
#         rk4: Runge-Kutta 4
#         rk3: Runge-Kutta 3
#     OTHER METHODS
#         f_augmented!: Create function with augmented state and control input
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# TODO: Change S to Z

"""
$(SIGNATURES)
In place Midpoint integration

Defines methods for both separated and augmented forms. Returns a discrete version
of a continuous dynamics function.

# Arguments
* f!: in place dynamics function, i.e. `f!(x?,x,u)`
* dt: time step

"""
function midpoint(f!::Function, dt::Float64)
    fd!(xdot,x,u,dt=dt) = begin
        f!(xdot,x,u)
        xdot .*= dt/2.
        f!(xdot, x + xdot, u)
        copyto!(xdot,x + xdot*dt)
    end
end

function midpoint(f_aug!::Function)
    fd_aug!(dS, S) = begin
        dt = S[end]^2
        f_aug!(dS, S)
        dS .*= dt/2.
        f_aug!(dS, S + dS)
        copyto!(dS,S + dS*dt)
    end
end

"""
$(SIGNATURES)
In place Runge Kutta 4 integration

Defines methods for both separated and augmented forms. Returns a discrete version
of a continuous dynamics function.

# Arguments
* f!: in place dynamics function, i.e. `f!(x?,x,u)` for separate or `f!(S?,S)` for augmented dynamics
* dt: time step
"""
function rk4(f!::Function, dt::Float64)
    # Runge-Kutta 4
    fd!(xdot,x,u,dt=dt) = begin
        k1 = k2 = k3 = k4 = zero(xdot)
        f!(k1, x, u);         k1 *= dt;
        f!(k2, x + k1/2, u); k2 *= dt;
        f!(k3, x + k2/2, u); k3 *= dt;
        f!(k4, x + k3, u);    k4 *= dt;
        copyto!(xdot, x + (k1 + 2*k2 + 2*k3 + k4)/6)
    end
end

function rk4(f_aug!::Function)
    # Runge-Kutta 4
    fd!(dS,S::Array) = begin
        dt = S[end]^2
        k1 = k2 = k3 = k4 = zero(S)
        f_aug!(k1,S);         k1 *= dt;
        f_aug!(k2,S + k1/2); k2 *= dt;
        f_aug!(k3,S + k2/2); k3 *= dt;
        f_aug!(k4,S + k3);    k4 *= dt;
        copyto!(dS, S + (k1 + 2*k2 + 2*k3 + k4)/6)
    end
end

"""
$(SIGNATURES)
In place Runge Kutta 3 integration

Defines methods for both separated and augmented forms. Returns a discrete version
of a continuous dynamics function.

# Arguments
* f!: in place dynamics function, i.e. `f!(x?,x,u)` for separate or `f!(S?,S)` for augmented dynamics
* dt: time step
"""
function rk3(f!::Function, dt::Float64)
        # Runge-Kutta 3 (zero order hold)
    fd!(xdot,x,u,dt=dt) = begin
        k1 = k2 = k3 = zero(x)
        f!(k1, x, u);               k1 *= dt;
        f!(k2, x + k1/2, u);       k2 *= dt;
        f!(k3, x - k1 + 2*k2, u);  k3 *= dt;
        copyto!(xdot, x + (k1 + 4*k2 + k3)/6)
    end
end

function rk3(f_aug!::Function)
    # Runge-Kutta 3 augmented (zero order hold)
    fd!(dS,S::Array) = begin
        dt = S[end]^2
        k1 = k2 = k3 = zero(S)
        f_aug!(k1,S);              k1 *= dt;
        f_aug!(k2,S + k1/2);      k2 *= dt;
        f_aug!(k3,S - k1 + 2*k2); k3 *= dt;
        copyto!(dS, S + (k1 + 4*k2 + k3)/6)
    end
end

"""
$(SIGNATURES)
Converts a separated dynamics function into an augmented dynamics function
"""
function f_augmented!(f!::Function, n::Int, m::Int)
    f_aug!(dS::AbstractArray, S::Array) = f!(dS, S[1:n], S[n+1:n+m])
end

function f_augmented(f::Function, n::Int, m::Int)
    f_aug(S::Array) = f(S[1:n], S[n+1:n+m])
end

function ZeroOrderHoldInterpolation(t,X)
    itr = interpolate(X,BSpline(Constant()))
    dt = t[2] - t[1]
    function zoh(t2)
        i2 = t2./dt .+ 1
        itr(floor.(i2))
    end
end

function MidpointInterpolation(t,X)
    Xm = [(X[i] + X[i+1])/2 for i = 1:length(X)-1]
    push!(Xm,Xm[end])
    ZeroOrderHoldInterpolation(t,Xm)
end


text:
using PackageCompiler

## Instruções para gerar sysimage com pacotes úteis pré-compilados, de forma a reduzir enormemente a penalidade
## de latência na carga de pacotes e compilação JIT de funções na 1a utilização - ver
## https://www.youtube.com/watch?v=d7avhSuK2NA para mais infos

PackageCompiler.create_sysimage(
    [:DataFrames, :CSV, :SQLite, :Tables, :StatsBase,
     :Queryverse, :Gadfly, :Cairo, :Colors, :Statistics];
    precompile_statements_file="precompile.jl"
    #, replace_default=false    # descomente se quiser substituir a sysimage padrão do Julia -
                                # ela pode ser restaurada com PackageCompiler::restore_default_sysimage, se necessário
    )

max_stars_repo_name: llrt/eda-qsa-rfb-jl

max_issues_repo_name: llrt/eda-qsa-rfb-jl

max_forks_repo_name: llrt/eda-qsa-rfb-jl


text:
using ArrayIteration: cartesian

A = [1 5 -5;
     0 3 2]
B = view(A, 1:2, 1:3)

@test each(index(A)) == eachindex(A) == 1:length(A)
@test each(index(B)) == CartesianIndices((1:2, 1:3))
@test each(index(A, :, 1:2)) == 1:4
@test each(index(B, :, 1:2)) == CartesianIndices((1:2, 1:2))
@test each(index(A, :, 2:3)) == 3:6
@test each(index(B, :, 2:3)) == CartesianIndices((1:2, 2:3))
@test each(index(A, 1, :)) == CartesianIndices((1, 1:3))

let k
k = 0
for v in each(A)
    @test v == A[k+=1]
end

k = 0
for I in each(index(A))
    @test A[I] == A[k+=1]
end

k = 0
for I in eachindex(A)
    @test A[I] == A[k+=1]
end

k = 0
for j in inds(A, 2)
    for v in each(A, :, j)
        @test v == A[k+=1]
    end
end

k = 0
for j in inds(A, 2)
    for I in each(index(A, :, j))
        @test A[I] == A[k+=1]
    end
end

k = 0
for j in inds(A, 2)
    for I in eachindex(stored(A, :, j))
        @test A[I] == A[k+=1]
    end
end

k = 0
for j in inds(A, 2)
    for v in each(stored(A, :, j))
        @test v == A[k+=1]
    end
end

end

A = Int[1 3; 2 4]
B = Array{Int}(undef, 2, 2)
C = PermutedDimsArray(Array{Int}(undef, 2, 2), [2,1])
R = reshape(view(Array{Int}(undef, 3, 2, 3), 1:2, 1:1, 1:2), (2, 2))

function badcopy!(dest, src)
    for (I, s) in zip(eachindex(dest), src)
        dest[I] = s
    end
    dest
end

fill!(B, -1)
@test badcopy!(B, A) == A
fill!(C, -1)
badcopy!(C, A)
@test_broken C[2,1] != A[2,1]   # oops!
@test_broken C[2,1] == A[1,2]

function goodcopy!(dest, src)
    for (I, s) in sync(index(dest), src)
        dest[I] = s
    end
    dest
end

fill!(B, -1)
@test goodcopy!(B, A) == A
fill!(C, -1)
goodcopy!(C, A)
@test C[2,1] == A[2,1]
@test C[1,2] == A[1,2]
fill!(R, -1)
@test goodcopy!(R, A) == A

D = ATs.OA(Array{Int}(undef, 2, 2), (-1, 2))
@test_throws DimensionMismatch goodcopy!(D, A)
E = ATs.OA(A, (-1,2))
goodcopy!(D, E)
@test D[0,3] == 1
@test D[1,3] == 2
@test D[0,4] == 3
@test D[1,4] == 4
D = ATs.OA(Array{Int}(undef,2,2), (-2,2))
@test_throws DimensionMismatch goodcopy!(D, E)
D = ATs.OA(Array{Int}(undef,2,2), (-1,1))
@test_throws DimensionMismatch goodcopy!(D, E)

function goodcopy2!(dest, src)
    for (Idest, Isrc) in sync(index(dest), index(src))
        dest[Idest] = src[Isrc]
    end
    dest
end

fill!(B, -1)
@test goodcopy2!(B, A) == A
fill!(C, -1)
goodcopy2!(C, A)
@test C[2,1] == A[2,1]
@test C[1,2] == A[1,2]
fill!(R, -1)
#@test goodcopy!(R, A) == A

#iter = sync(index(A), index(R))
#(IA, IR) = first(iter)
#@test isa(IR, Base.ReshapedIndex)

# TODO: uncomment when sync+stored is implemented
# function goodcopy3!(dest, src)
#     for (I, s) in sync(index(dest), stored(src))
#         dest[I] = s
#     end
#     dest
# end

# fill!(B, -1)
# @test goodcopy3!(B, A) == A
# fill!(C, -1)
# goodcopy3!(C, A)
# @test C[2,1] == A[2,1]
# @test C[1,2] == A[1,2]

# 1-argument case
function mysum1a(A)
    s = 0.0
    for (a,) in sync(A)
        s += a
    end
    s
end

function mysum1b(A)
    s = 0.0
    for (I,) in sync(index(A))
        s += A[I]
    end
    s
end

@test mysum1a(A) ≈ sum(A)
@test mysum1b(A) ≈ sum(A)

# 3-argument form
function mysum!(dest, A, B)
    for (Idest, a, b) in sync(index(dest), A, B)
        dest[Idest] = a + b
    end
    dest
end

C = PermutedDimsArray([10 30; 20 40], [2,1])
D = fill(-1, (2,2))
@test mysum!(D, A, C) == Int[11 23; 32 44]
D = fill(-1, (2,2))
@test mysum!(D, C, A) == Int[11 23; 32 44]
Cc = goodcopy!(similar(A), C)
@test isa(Cc, Array)
D = fill(-1, (2,2))
@test mysum!(D, Cc, A) == Int[11 23; 32 44]
goodcopy!(D, Cc)
mysum!(C, D, A)
@test C[1,2] == 23
@test C[2,1] == 32

# sync with dimension iterators
fill!(B, -1)
for (I, a) in sync(index(B, :, 2), value(A, :, 1))
    B[I] = a
end
@test B == [-1 A[1,1]; -1 A[2,1]]
fill!(B, -1)
for (IB, IA) in sync(index(B, :, 1), index(A, :, 2))
    B[IB] = A[IA]
end
@test B == [A[1,2] -1; A[2,2] -1]

## optimized ReshapedArray iterators
A = reshape(1:12, 4, 3)    # LinearFast
@test each(index(A, :, 2)) == 5:8
@test each(index(A, 2:3, 3)) == 10:11
A = view(copy(reshape(1:15, 5, 3)), 1:4, :)  # LinearSlow
a = reshape(A, 12)

let k

@test each(index(a, 1:4)) == CCI(CartesianIndices(size(A)), cartesian(CartesianIndex(1,1),CartesianIndex(4,1)))
k = 0
for I in each(index(a, 1:4))
    @test a[I] == a[k+=1]
end
@test each(index(a, 3:9)) == CCI(CartesianIndices(size(A)),
                                 cartesian(CartesianIndex(3,1),CartesianIndex(1,3)))
k = 2
for I in each(index(a, 3:9))
    @test a[I] == a[k+=1]
end
end

function sum_cols_slow!(S, A)  # slow for ReshapedArrays
    fill!(S, 0)
    @assert inds(S,2) == inds(A,2)
    for j in inds(A,2)
	tmp = S[1,j]
	@inbounds for i in inds(A, 1)
	    tmp += A[i,j]
	end
	S[1,j] = tmp
    end
    S
end

function sum_cols_fast!(S, A)
    fill!(S, 0)
    @assert inds(S,2) == inds(A,2)
    for j in inds(A,2)
	tmp = S[1,j]
	@inbounds for I in each(index(A, :, j))
	    tmp += A[I]
	end
	S[1,j] = tmp
    end
    S
end

const can_inline = Base.JLOptions().can_inline == 1
dim1 = can_inline ? 10^6 : 10
A = rand(dim1,2,5)
B = view(A, 1:size(A,1)-1, :, :)
R = reshape(B, (size(B,1),prod(size(B)[2:end])))
@test isa(each(index(R, :, 2)), CCI)
S1 = zeros(1,size(R,2))
S2 = similar(S1)
@test sum_cols_fast!(S1, R) == sum_cols_slow!(S2, R)

let nfailures
if can_inline
    nfailures = 0
    for i = 1:5
        ts = @elapsed sum_cols_slow!(S2, R)
        tf = @elapsed sum_cols_fast!(S2, R)
        nfailures += tf > ts
    end
    @test_broken nfailures <= 1
end
end


text:
module Models

using ..FastAI

using BSON
using Flux
using Zygote
using DataDeps

include("layers.jl")
include("blocks.jl")

include("xresnet.jl")
include("unet.jl")

export xresnet18, xresnet50, UNetDynamic

end


text:
module TestMethods

  using GBIF
  using Test

  set = occurrences()

  @test typeof(set[1]) == GBIFRecord
  @test length(set[1:4]) == 4

end


text:
##
## METIS
## Minimum interface to the Metis graph partitioner
##

const METIS_NOPTIONS = 40
const METIS_OK = Int32(1)
const metis_options = - ones(Cint, METIS_NOPTIONS)
const libmetis = "/Users/lcambier/.julia/packages/Homebrew/l8kUw/deps/usr/lib/libmetis.dylib"

function vertex_sep_fast!(n, colptr::Array{Cint,1}, rowval::Array{Cint,1}, part::Array{Cint,1})
    sepSize = zeros(Cint,1)
    n2 = Cint(n)
    err = ccall((:METIS_ComputeVertexSeparator,libmetis), Cint,
                (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},     Ptr{Cint}, Ptr{Cint}),
                Ref(n2),    colptr,    rowval,    C_NULL,    metis_options, sepSize,   part)
    err == METIS_OK || error("METIS_ComputeVertexSeparator returned error code $err")
    return
end

function vertex_sep_bissect_fast!(n, colptr::Array{Cint,1}, rowval::Array{Cint,1}, part::Array{Cint,1})
    nparts = Cint(2)
    one = Cint(1)
    two = Cint(2)
    n2 = Cint(n)
    objval = zeros(Cint,1)
    err = ccall((:METIS_PartGraphRecursive,libmetis), Int32,
                (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                 Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                 Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),
                Ref(n2), Ref(one), colptr, rowval, C_NULL, C_NULL, C_NULL, Ref(two),
                C_NULL, C_NULL, metis_options, objval, part)
    err == METIS_OK || error("METIS_PartGraphRecursive returned error code $err")
    for i1 = 1:n # update part: those with (0), connected to at least one (1), gets a (2)
        if part[i1] == 1
            continue
        end
        for k0 = colptr[i1]:colptr[i1+1]-1 # 0-based
            k1 = k0 + 1
            j0 = rowval[k1]
            j1 = j0 + 1
            if part[j1] == 1
                part[i1] = 2 # Becomes a sep
                break
            end
        end
    end
end


text:
using BOMoD
using Random
using StatsBase
using LinearAlgebra
using CSV
include("Binary_Quadratic.jl")
###################
#Set Design and make space
##################

# random
rng = MersenneTwister(1)
# set up BOMoD use
# make binary moduels
mod0 = Mod(Symbol(0))
mod1 = Mod(Symbol(1))

# no constrains model of lenght 10, not all constructs are generated
all_moduels = Group_Mod([mod0,mod1])
σ²_mean = 10^-6
# the design
design = construct_design(all_moduels,lenght_c)
# the space
space = BOMoD.getspace(design)

####################
#settings evaluation
###################
# evaluad differet BO
# number of new Q
n_Q = 1
# number of simimulations
const n_iter = 10
const n_cycles = 4
# number of BO cycles
####################
#toy_data settings
###################
const lenght_c = 10
const λ = 0
const Lc = 10

##################
#BO-fixed settings
const Test_kernel = LevStehnoexp(1) # Kernel
n_samples = 10 # number of sampels in batch
##################
# setting up the toydata
const Q = compute_decay(rng,lenght_c,Lc)
# xstar, optimal construct

#Mykernel = SumKernel(CosStehno(),LevStehnoexp(1),1)
Mykernel = LevStehnoexp(1)
# matrix to save diff  x* and the max predicted construct.
# row represents independet cycels
# collumes represent the number of the BO prosses is repeated

max_random = Matrix{Float64}(undef, n_iter, n_cycles)
max_model = Matrix{Float64}(undef, n_iter, n_cycles)

for i in 1:n_iter
    println("n =$i")
    #choose random n points, x are index of the row in dataframe to link proper y
    points = sample(rng,space,n_samples,with_index = true )
    # liste to save the BO outcome an compare to radom
    #modelpoints = Array{Float64}(undef,n_samples)
    #radompoints = Array{Float64}(undef,n_samples)
    # set starting point
    #index_train_BO = points[:,2]
    #index_train_random = points[:,2]
    construct_train_BO = points[:,1]
    construct_train_random = points[:,1]

    # do the BO in n_cycles
    for j in 1:n_cycles
        # sample random
        # filter sample constructs
        x_test_random = [construct for construct in space if sum([isequal(construct,train_con) for train_con in construct_train_random]) == 0 ]
        myrandompoints = StatsBase.sample(rng,x_test_random,n_samples)

        # construct generated y-values
        y_train = map(x -> unpackconstruct(x) |> (x -> quad(x,Q,λ)),construct_train_BO)
        # use BO
        x_test_BO = [construct for construct in space if sum([isequal(construct,train_con) for train_con in construct_train_BO]) == 0 ]
        Model = gp_optimised(construct_train_BO,y_train,Test_kernel)
        new_predictions = save_thompson_sampling(Model.GP_model,x_test_BO,10,10^-6)

        #save new constructs
        construct_train_random = [construct_train_BO;myrandompoints[:,1]]
        construct_train_BO = [construct_train_BO;new_predictions[:,1]]

        #max values
        max_random[i,j] =  maximum(map(x -> unpackconstruct(x) |> (x -> quad(x,Q,λ)),construct_train_random))
        max_model[i,j] = maximum(map(x -> unpackconstruct(x) |> (x -> quad(x,Q,λ)),construct_train_BO ))

    end
end

writedlm("C:\\Users\\samdi\\.julia\\dev\\BOMoD\\scripts\\model_max.txt",max_model)
writedlm("C:\\Users\\samdi\\.julia\\dev\\BOMoD\\scripts\\random_max.txt",max_random)
gr()
plt = plot(xaxis ="x",yaxis = "maxf(x)");
plot1 = makeplot(max_random,"random",plt,:red);
plot2 = makeplot(max_model,"BO",plt,:blue);

plot(plt)
function makeplot(Data_matrix,name,plt,c= :blue)
    μ = [mean(Data_matrix[:,i]) for i in 1:size(Data_matrix)[2]]
    sigma = [std(Data_matrix[:,i]) for i in 1:size(Data_matrix)[2]]
    x_plot = collect(1:size(Data_matrix)[2])
    return plot!(plt,x_plot,μ,ribbon = sigma,label=name,color = c,lw=1.5,fillalpha = 0.3)
end


text:
using LinearAlgebra, FastSolversForWeightedTV, Flux, TestImages, BenchmarkTools
using CUDA; CUDA.allowscalar(false)
using Random; Random.seed!(123)

# Dimension
# dim = 2
dim = 3

# Optimization options
flag_gpu = true
# flag_gpu = false
dim == 2 && (steplength = 1f0/8f0)
dim == 3 && (steplength = 1f0/12f0)
opt = opt_fista(steplength; niter=5, tol_x=nothing, nesterov=true)

# TV norm
n = tuple(256*ones(Int64, dim)...)
h = tuple(ones(Float32, dim)...)
g = gradient_norm(2,1,n,h; T=ComplexF32); flag_gpu && (g = g |> gpu)

# Numerical phantom
y = Float32.(TestImages.shepp_logan(n[1:2]...))
dim == 3 && (y = repeat(y; outer=(1,1,n[3])))
y = y/norm(y, Inf)
y += 0.1f0*randn(ComplexF32, n)
flag_gpu && (y = y |> gpu)

# Proxy
λ = 0.5f0*norm(y)^2/g(y)
@btime proxy(y, λ, g, opt);

# Projection test
ε = 0.9f0*g(y)
@btime project(y, ε, g, opt);

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
"""
$(TYPEDEF)

A `Shape` is used for visualization of a [`Body`](@ref). All `Shape`s have an `xoffset` and a `qoffset` parameter to specify the position and orientation offset between the actual [`Body`](@ref) and the visualized `Shape`.
# Important attributes
* `xoffset`: The relative position offset (vector) between `Body` and `Shape`.
* `qoffset`: The relative orientation offset (quaternion) between `Body` and `Shape`.
* `color`:   The color of the `Shape`.
"""
abstract type Shape{T} end

struct EmptyShape{T} <: Shape{T}
    EmptyShape() = new{Float64}()
end

"""
$(TYPEDEF)

A `Mesh` can be used to visualize arbitrary geometries.
# Important attributes
* `path`: The path to the geometry file.
"""
mutable struct Mesh{T} <: Shape{T}
    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    path::String
    scale::SVector{3,T}
    color::RGBA

    function Mesh(path::String;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((xoffset, qoffset))...)

        new{T}(xoffset, qoffset, path, scale, color)
    end

    function Mesh(path::String, m::Real, J::AbstractMatrix;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), name::String="", color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((m, J, xoffset, qoffset))...)

        return Body(m, J; name=name, shape=new{T}(xoffset, qoffset, path, scale, color))
    end
end

"""
$(TYPEDEF)

A `Box`.
# Important attributes
* `xyz`: The box size in x, y, and z direction (vector).
"""
mutable struct Box{T} <: Shape{T}
    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    xyz::SVector{3,T}
    scale::SVector{3,T}
    color::RGBA

    function Box(x::Real, y::Real, z::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((x, y, z, xoffset, qoffset))...)

        new{T}(xoffset, qoffset, [x;y;z], scale, color)
    end

    function Box(x::Real, y::Real, z::Real, m::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), name::String="", color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((x, y, z, m, xoffset, qoffset))...)
        J = 1 / 12 * m * diagm([y^2 + z^2;x^2 + z^2;x^2 + y^2])

        return Body(m, J; name=name, shape=new{T}(xoffset, qoffset, [x;y;z], scale, color))
    end
end

"""
$(TYPEDEF)

A `Cylinder` along the z-axis.
# Important attributes
* `rh`: The cylinder radius and height (vector).
"""
mutable struct Cylinder{T} <: Shape{T}
    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    rh::SVector{2,T}
    scale::SVector{3,T}
    color::RGBA

    # Cylinder points in the z direction
    function Cylinder(r::Real, h::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((r, h, xoffset, qoffset))...)

        new{T}(xoffset, qoffset, [r;h], scale, color)
    end

    function Cylinder(r::Real, h::Real, m::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), name::String="", color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((r, h, m, xoffset, qoffset))...)
        J = 1 / 2 * m * diagm([r^2 + 1 / 6 * h^2;r^2 + 1 / 6 * h^2;r^2])

        return Body(m, J; name=name, shape=new{T}(xoffset, qoffset, [r;h], scale, color))
    end
end

"""
$(TYPEDEF)

A `Sphere`.
# Important attributes
* `r`: The sphere radius.
"""
mutable struct Sphere{T} <: Shape{T}
    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    r::T
    scale::SVector{3,T}
    color::RGBA

    function Sphere(r::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((r, xoffset, qoffset))...)

        new{T}(xoffset, qoffset, r, scale, color)
    end

    function Sphere(r::Real, m::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), name::String="", color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((r, m, xoffset, qoffset))...)
        J = 2 / 5 * m * diagm([r^2 for i = 1:3])

        return Body(m, J; name=name, shape=new{T}(xoffset, qoffset, r, scale, color))
    end
end

"""
$(TYPEDEF)

A square `Pyramid` with a base in the x-y-plane.
* `wh`: The pyramid width and height (vector).
"""
mutable struct Pyramid{T} <: Shape{T}
    xoffset::SVector{3,T}
    qoffset::UnitQuaternion{T}

    wh::SVector{2,T}
    scale::SVector{3,T}
    color::RGBA

    # Pyramid points in the z direction, Center of mass at 1/4 h
    function Pyramid(w::Real, h::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((w, h, xoffset, qoffset))...)

        new{T}(xoffset, qoffset, [w;h], scale, color)
    end

    function Pyramid(w::Real, h::Real, m::Real;
            xoffset::AbstractVector = szeros(3), qoffset::UnitQuaternion = one(UnitQuaternion),
            scale::AbstractVector = sones(3), name::String="", color = RGBA(0.75, 0.75, 0.75)
        )
        T = promote_type(eltype.((w, h, m, xoffset, qoffset))...)
        J = 1/80 * m * diagm([4*w^2+3*h^2;4*w^2+3*h^2;8*w^2])

        return Body(m, J; name=name, shape=new{T}(xoffset, qoffset, [w;h], scale, color))
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, shape::Shape{T}) where {T}
    summary(io, shape)
    println(io,"")
    println(io," xoffset: "*string(shape.xoffset))
    println(io," qoffset: "*string(shape.qoffset))
    println(io," color:   "*string(shape.color))
    println(io," scale:   "*string(shape.scale))
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, shape::EmptyShape{T}) where {T}
    summary(io, shape)
    println(io,"")
end


text:
using JLD2
using DataFrames
using StatsPlots

@load "paper/jls.jld2" jls_df

p = @df jls_df[jls_df[:, :Type] .== "AT", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoPlot_AT.pdf")
p = @df jls_df[jls_df[:, :Type] .== "AS", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoPlot_AS.pdf")
p = @df jls_df[jls_df[:, :Type] .== "BS", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoPlot_BS.pdf")

# @df jls_df[jls_df[:, :Type] .== "AT", :] scatter(:n, :std, group=:N)
# @df jls_df[jls_df[:, :Type] .== "AS", :] scatter(:n, :std, group=:N)
# @df jls_df[jls_df[:, :Type] .== "BS", :] scatter(:n, :std, group=:N)

@load "paper/jlsInd.jld2" jlsInd_df

p = @df jlsInd_df[jlsInd_df[:, :Type] .== "AT", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoIndPlot_AT.pdf")
p = @df jlsInd_df[jlsInd_df[:, :Type] .== "AS", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoIndPlot_AS.pdf")
p = @df jlsInd_df[jlsInd_df[:, :Type] .== "BS", :] plot(:n, :mean, group=:N, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks = [50; 400; 800; 1600; 3200],
  size=(400,250))
Plots.savefig(p, "paper/lgDemoIndPlot_BS.pdf")


text:
@system IntercellularSpace(Weather) begin
    A_net ~ hold
    #TODO: interface between boundary/stomata/intercellular space (i.e. soil layers?)
    gvc ~ hold

    #FIXME: duplicate in Stomata
    Ca(CO2, P_air): co2_air => (CO2 * P_air) ~ track(u"μbar")

    #HACK: high temperature simulation requires higher upper bound
    Cimax(Ca): intercellular_co2_upper_limit => 2Ca ~ track(u"μbar")
    Cimin: intercellular_co2_lower_limit => 0 ~ preserve(u"μbar")
    Ci(Ca, Ci, A_net, gvc): intercellular_co2 => begin
        Ca - Ci ⩵ A_net / gvc
    end ~ bisect(min=Cimin, upper=Cimax, u"μbar")
end