text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Poppler"
version = v"0.87.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://poppler.freedesktop.org/poppler-0.87.0.tar.xz", "6f602b9c24c2d05780be93e7306201012e41459f289b8279a27a79431ad4150e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/poppler-*/

export CXXFLAGS="-I${prefix}/include/openjpeg-2.3"

if [[ "${target}" == "${MACHTYPE}" ]]; then
    # When building for the host platform, the system libexpat is picked up
    rm /usr/lib/libexpat.so*
fi

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_GTK_TESTS=OFF \
    -DENABLE_CMS=lcms2 \
    -DENABLE_GLIB=ON \
    -DENABLE_QT5=OFF \
    -DENABLE_UNSTABLE_API_ABI_HEADERS=ON \
    -DWITH_GObjectIntrospection=OFF \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libpoppler-cpp", :libpoppler_cpp),
    LibraryProduct("libpoppler-glib", :libpoppler_glib),
    LibraryProduct("libpoppler", :libpoppler),
    ExecutableProduct("pdfattach", :pdfattach),
    ExecutableProduct("pdfdetach", :pdfdetach),
    ExecutableProduct("pdffonts", :pdffonts),
    ExecutableProduct("pdfimages", :pdfimages),
    ExecutableProduct("pdfinfo", :pdfinfo),
    ExecutableProduct("pdfseparate", :pdfseparate),
    ExecutableProduct("pdftocairo", :pdftocairo),
    ExecutableProduct("pdftohtml", :pdftohtml),
    ExecutableProduct("pdftoppm", :pdftoppm),
    ExecutableProduct("pdftops", :pdftops),
    ExecutableProduct("pdftotext", :pdftotext),
    ExecutableProduct("pdfunite", :pdfunite),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Cairo_jll"),
    Dependency("Fontconfig_jll"),
    # Dependency("GTK3_jll"),
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("JpegTurbo_jll"),
    # TODO: v4.3.0 is available, use that next time
    Dependency("Libtiff_jll"; compat="4.1.0"),
    Dependency("OpenJpeg_jll"),
    # Dependency("gdk_pixbuf_jll"),
    Dependency("libpng_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5")

ext: jl

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
#--------------------------------------------------------------------#
#     Input data used for heat and moisture transfer calculation     #
#--------------------------------------------------------------------#

#DODATI RUBNE UVJETE

#########
#Timestep if necessary#
##########

#dt=1        #Timestep in h
#T=8760      #Time domain

###########
#Mesh size#
###########

dl=0.02#defined in m

########################
#Geometry and materials#
########################

#Elements -- matrix of vectors E_1=[material_key,x_0,y_0,x_1,y_1] to E_n

ne=23 #number of elements

E=zeros(Float64,ne,9) #Geometry and material matrix + boundary conditions

#Element 1 (Concrete)
E[1,1]=12
E[1,2]=0
E[1,3]=0
E[1,4]=0.76
E[1,5]=0.08
E[1,6]=1
E[1,7]=1.65
E[1,8]=2200
E[1,9]=1000

#Element 2 (XPS)
E[2,1]=14
E[2,2]=0
E[2,3]=0.08
E[2,4]=0.08
E[2,5]=1
E[2,6]=2
E[2,7]=0.035
E[2,8]=50
E[2,9]=1015

#Element 3 (Reinforced concrete)
E[3,1]=13
E[3,2]=0.08
E[3,3]=0.08
E[3,4]=0.68
E[3,5]=1
E[3,6]=3
E[3,7]=2.6
E[3,8]=2500
E[3,9]=1000

#Element 4 (Reinforced concrete)
E[4,1]=13
E[4,2]=0.16
E[4,3]=1
E[4,4]=0.68
E[4,5]=1.08
E[4,6]=4
E[4,7]=2.6
E[4,8]=2500
E[4,9]=1000

#Element 5 (XPS)
E[5,1]=14
E[5,2]=0
E[5,3]=1
E[5,4]=0.08
E[5,5]=1.08
E[5,6]=5
E[5,7]=0.035
E[5,8]=50
E[5,9]=1015

#Element 6 (XPS)
E[6,1]=14
E[6,2]=0.08
E[6,3]=1
E[6,4]=0.16
E[6,5]=1.08
E[6,6]=6
E[6,7]=0.035
E[6,8]=50
E[6,9]=1015

#Element 7 (PVC hydroinsulation)
#E[7,1]=15
#E[7,2]=0.08
#E[7,3]=1
#E[7,4]=0.16
#E[7,5]=1.01
#E[7,6]=7
#E[7,7]=0.14
#E[7,8]=1200
#E[7,9]=1000

#Element 8 (PVC hydroinsulation)
#E[8,1]=15
#E[8,2]=0.15
#E[8,3]=1.01
#E[8,4]=0.16
#E[8,5]=1.77
#E[8,6]=8
#E[8,7]=0.14
#E[8,8]=1200
#E[8,9]=1000

#Element 9 (XPS)
E[7,1]=14
E[7,2]=0.68
E[7,3]=0.08
E[7,4]=0.76
E[7,5]=1.08
E[7,6]=7
E[7,7]=0.035
E[7,8]=50
E[7,9]=1015

#Element 10 (Concrete)
E[8,1]=12
E[8,2]=0.76
E[8,3]=1
E[8,4]=1.4
E[8,5]=1.08
E[8,6]=8
E[8,7]=1.65
E[8,8]=2200
E[8,9]=1000

#Element 11 (XPS)
E[9,1]=14
E[9,2]=0.08
E[9,3]=1.08
E[9,4]=0.16
E[9,5]=1.86
E[9,6]=9
E[9,7]=0.035
E[9,8]=50
E[9,9]=1015

#Element 12 (Reinforced concrete)
E[10,1]=13
E[10,2]=0.16
E[10,3]=1.08
E[10,4]=0.36
E[10,5]=1.86
E[10,6]=10
E[10,7]=2.6
E[10,8]=2500
E[10,9]=1000

#Element 13 (PVC hydroinsulation)
#E[13,1]=15
#E[13,2]=0.36
#E[13,3]=1.08
#E[13,4]=0.37
#E[13,5]=1.27
#E[13,6]=13
#E[13,7]=0.14
#E[13,8]=1200
#E[13,9]=1000

#Element 14 (PVC hydroinsulation)
#E[14,1]=15
#E[14,2]=0.37
#E[14,3]=1.08
#E[14,4]=1.45
#E[14,5]=1.09
#E[14,6]=14
#E[14,7]=0.14
#E[14,8]=1200
#E[14,9]=1000

#Element 15 (Mineral wool)
E[11,1]=11
E[11,2]=0.36
E[11,3]=1.08
E[11,4]=1.4
E[11,5]=1.18
E[11,6]=11
E[11,7]=0.035
E[11,8]=150
E[11,9]=1030

#Element 16 (Cement screed)
E[12,1]=16
E[12,2]=0.36
E[12,3]=1.18
E[12,4]=1.4
E[12,5]=1.26
E[12,6]=12
E[12,7]=1.6
E[12,8]=2000
E[12,9]=1100

##############################
#Boundary condition materials#
##############################

#Element 17 (Constant temperature soil)
E[13,1]=100
E[13,2]=0
E[13,3]=-dl
E[13,4]=0.76
E[13,5]=0
E[13,6]=13
E[13,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[13,8]=5
E[13,9]=1

#Element 18 (Variable temperature soil 0-108 cm)
E[14,1]=101
E[14,2]=-dl
E[14,3]=0
E[14,4]=0
E[14,5]=1.08
E[14,6]=18
E[14,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[14,8]=5
E[14,9]=1

#Element 19 (Variable temperature soil 108-124 cm)
E[15,1]=102
E[15,2]=0.08-dl
E[15,3]=1.08
E[15,4]=0.08
E[15,5]=1.24
E[15,6]=15
E[15,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[15,8]=5
E[15,9]=1

#Element 20 (Environmental conditions)
E[16,1]=103
E[16,2]=0.08-dl
E[16,3]=1.24
E[16,4]=0.08
E[16,5]=1.86
E[16,6]=16
E[16,7]=2/dl*(25)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[16,8]=5
E[16,9]=1

#Element 21 (Adiabatic boundary conditions)
E[17,1]=104
E[17,2]=0.08
E[17,3]=1.86
E[17,4]=0.36
E[17,5]=1.86+dl
E[17,6]=17
E[17,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[17,8]=1
E[17,9]=-1 #Adiabatic BC key --> -1

#Element 22 (Internal conditions)
E[18,1]=105
E[18,2]=0.36
E[18,3]=1.26
E[18,4]=0.36+dl
E[18,5]=1.86
E[18,6]=18
E[18,7]=2/dl*(8)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[18,8]=22
E[18,9]=1

#Element 23 (Internal conditions)
E[19,1]=106
E[19,2]=0.36+dl
E[19,3]=1.26
E[19,4]=1.4
E[19,5]=1.26+dl
E[19,6]=19
E[19,7]=2/dl*(8)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[19,8]=22
E[19,9]=1

#Element 24 (Adiabatic boundary conditions)
E[20,1]=107
E[20,2]=1.4
E[20,3]=1
E[20,4]=1.4+dl
E[20,5]=1.26
E[20,6]=20
E[20,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[20,8]=1
E[20,9]=-1 #Adiabatic BC key --> -1

#Element 25 (Constant temperature soil)
E[21,1]=108
E[21,2]=0.76+dl
E[21,3]=1-dl
E[21,4]=1.4
E[21,5]=1
E[21,6]=21
E[21,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[21,8]=5
E[21,9]=1

#Element 26 (Variable temperature soil 0-100 cm)
E[22,1]=109
E[22,2]=0.76
E[22,3]=0
E[22,4]=0.76+dl
E[22,5]=1
E[22,6]=22
E[22,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[22,8]=5
E[22,9]=1

#Element 27 (Constant temperature soil)
E[23,1]=110
E[23,2]=0
E[23,3]=1.08
E[23,4]=0.08-dl
E[23,5]=1.08+dl
E[23,6]=23
E[23,7]=2/dl*(1000)#(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[23,8]=5
E[23,9]=1


text:
# Based on https://www.amazon.com/Quantum-Monte-Carlo-Methods-Algorithms/dp/1107006422

"""
    std_error(F::FullBinner)

Estimate of the standard error of the time series mean by performing a full binning analysis.

"Full" binning means that we bin the data multiple times, considering all bin sizes compatible
with the length of the time series.
"""
function std_error end

std_error(F::FullBinner{<:Number}) = _std_error(F)
function std_error(F::FullBinner{<:AbstractArray})
    v = VectorOfArray(F)
    r = similar(F[1], real(eltype(v))) # std error is always real
    @inbounds for i in eachindex(F[1])
        r[i] = _std_error(v[i,:]) # TODO: maybe a view might be faster here
    end
    return r
end

"""
    tau(F::FullBinner[, std_error])

Calculate the autocorrelation time tau.

If `std_error` is provided, no binning is necessary (much faster).
"""
function tau(F::FullBinner) end

tau(F::FullBinner{<:AbstractArray}, stderr = std_error(F)) = 0.5*(stderr.^2 * length(F) ./ var(F) .- 1)
tau(F::FullBinner{<:Number}, stderr = std_error(F)) = 0.5*(stderr^2 * length(F) / var(F) - 1)
_tau(Rvalue::Float64) = (Rvalue - 1)/2

@inline _std_error(X::AbstractVector{<:Real}) = _std_error_from_R_function(X)
@inline function _std_error(X::AbstractVector{<:Complex})
    std_err_real = _std_error_from_R_function(real(X))
    std_err_imag = _std_error_from_R_function(imag(X))
    sqrt(std_err_real^2 + std_err_imag^2)
end

"""
Estimates the error coefficient `R` by considering `<R>(bss)` (averaged R values)
and taking the largest value compatible with at least 32 bins.
"""
@inline function _std_error_from_R_function(X::AbstractVector{T}) where T<:Real
    bss, Rs, means = R_function(X; min_nbins=32)
    R = length(means) > 0 ? means[end] : 1.0
    # conv = isconverged(X, means)
    # bs = findfirst(r -> r >= R, Rs)
    return _R2std_error(X, R)
end

"""
    all_binning_errors(X) -> bs, stds, cum_stds

Returns all compatible bin sizes `bs`, the corresponding standard errors `stds`,
and the cumulative standard errors `cum_stds` (reduced fluctuations).

Only bin size that lead to at least 32 bins are considered.
"""
function all_binning_errors end

all_binning_errors(F::FullBinner{<:Number}) = _all_binning_errors(F)

function _all_binning_errors(X::AbstractVector{<:Real})
    bss, Rs, means = R_function(X; min_nbins=32)
    return bss, _R2std_error.(Ref(X), Rs), _R2std_error.(Ref(X), means)
end
function _all_binning_errors(X::AbstractVector{<:Complex})
    bss_real, stds_real, cum_stds_real = _all_binning_errors(real(X))
    bss_imag, stds_imag, cum_stds_imag = _all_binning_errors(imag(X))
    # @assert bss_real == bss_imag
    stds = sqrt.(stds_real.^2 .+ stds_imag.^2)
    cum_stds = sqrt.(cum_stds_real.^2 .+ cum_stds_imag.^2)
    return bss_real, stds, cum_stds
end

# TODO: Improve this (until then, don't export)
"""
  isconverged(X)

Checks if the estimation of the one sigma error is converged.

Returns `true` once the mean `R` value is converged up to 0.1% accuracy.
This corresponds to convergence of the error itself up to ~3% (sqrt).
"""
function isconverged(X::AbstractVector{T}) where T<:Real
  bss, Rs, means = R_function(X)
  return isconverged(X, means)
end
function isconverged(X::AbstractVector{T}, means::AbstractVector) where T<:Real
    len = length(means)
    len < 10 && (return false)
    lastn = min(len, 200)
    start = len-(lastn-1)
    return @views maximum(abs.(diff(means[start:end])))<1e-3 # convergence condition
end
isconverged(X::AbstractVector{<:Complex}) = isconverged(real(X)) && isconverged(imag(X)) # check if that's really what we want here


text:
using Plots
plotlyjs()
using PlutoUI
using FFTW
using MRIReco
using LinearAlgebra
using IterativeSolvers

#p = scatter(rand(10), rand(10));
#display(p)
println("Warm-up done.")


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

using Markdown
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 878a4750-b15e-11ea-2584-8feba490699f
using Test

# ╔═╡ fd0763a0-b163-11ea-23b4-a7bae7052e19
md"# File picker"

# ╔═╡ 592d3ac0-b15e-11ea-313c-db05bc7303cc
md"# Notebook interaction"

# ╔═╡ 6dde0352-b15e-11ea-2fa8-7327cc366c1a
md"## Running multiple cells"

# ╔═╡ 7370dcc0-b15e-11ea-234b-23584c864b61
ma = 1

# ╔═╡ 75b21a30-b15e-11ea-3046-2170ec097e63
mb = 2

# ╔═╡ 7b74dd40-b15e-11ea-291a-d7e10a185718
@test ma + mb == 3

# ╔═╡ f0b821b0-b15f-11ea-1f64-dd33aa85b54e
md"## Moving cells"

# ╔═╡ 4980fc10-b163-11ea-081b-c1335699a8f6
md"## Stopping cells"

# ╔═╡ 4e0c5090-b163-11ea-3cf9-a1b17dd2f053

# ╔═╡ d2c1d090-b162-11ea-0c17-2b234c098cf9
md"# CodeMirror"

# ╔═╡ d890a190-b162-11ea-31dd-8d603787e5c5
md"## Autocomplete"

# ╔═╡ e141f910-b162-11ea-039b-3ba1414cbd07

# ╔═╡ 39724bae-b160-11ea-0af5-737755a57e48
md"# Network"

# ╔═╡ f33dc5b0-b160-11ea-15e3-1b54697b5872

# ╔═╡ 95492630-b15f-11ea-3f16-410fe87cd32c
md"# Errors"

# ╔═╡ c1c894f0-b166-11ea-31e0-dd1cddb7b7ca
md"## Multiple definitions"

# ╔═╡ d05d7df0-b166-11ea-2d8f-39b31a42890d
1
2

# ╔═╡ bdeba8e0-b166-11ea-3913-3138ab428083
md"## Stack traces"

# ╔═╡ 84888e20-b160-11ea-1d61-c5934251d6dd
html"

"

# ╔═╡ 9dc4a0a0-b15f-11ea-361c-87742cf3f2a2
function ef(x)

	sqrt(-x)
end

# ╔═╡ aab109c0-b15f-11ea-275d-31e21fcda8c4
ef(1)

# ╔═╡ 976bc2a0-b160-11ea-3e7a-9f033b0f2daf
function eg(x)

	sqrt(-x)
end

# ╔═╡ 9c74f9b2-b160-11ea-35fb-917cb1120f5b
eg(1)

# ╔═╡ ea3f77f0-b166-11ea-046e-ef39bfc57d0f
md"## Bad errors"

# ╔═╡ ee92a570-b166-11ea-38fa-0f3563b3f790
TODO: add that error showing error to the Julia tests

# ╔═╡ 4f32c16e-b15e-11ea-314e-01700baeb92d
md"# Bonds"

# ╔═╡ 3a14b3f0-b165-11ea-153d-796416ee5ccc
md"## Lossy"

# ╔═╡ 41a75500-b165-11ea-2519-bbd0feaef6cf
@bind bl1 html"<input type='range' max='100000'>"

# ╔═╡ 4ccbf670-b165-11ea-1951-c17ffb8a58cf
sleep(.5); bl1

# ╔═╡ 8bb26902-b165-11ea-048c-d7f7a72006ee
@assert bl1 isa Int64

# ╔═╡ e559eaf0-b165-11ea-0d81-ffc480afe8f3
@bind bl2 html"<input type='range' max='100000'>"

# ╔═╡ e63be680-b165-11ea-0fd3-bd4e0bf92eb8
bl2

# ╔═╡ 59966a90-b163-11ea-1786-e56e45f06dd0
md"## Recursive"

# ╔═╡ 55bade10-b163-11ea-327f-132c93471713
md"## Scrolling"

# ╔═╡ f2c0bb90-b162-11ea-24a1-3f864a09e5ee
@bind bw1 html"<input type='range' value='0'>"

# ╔═╡ 56e6f440-b15e-11ea-1327-09932af5b5bd
HTML("

")

# ╔═╡ 2296ac80-b163-11ea-3d00-ed366fa9ce3e
@bind bw2 html"<input type='range' value='0'>"

# ╔═╡ 20d72230-b163-11ea-39c2-69bf2c422d50
HTML("

")

# ╔═╡ 54ccad80-b163-11ea-1777-c5364408948e

# ╔═╡ 32b5edc0-b15d-11ea-09d6-3b889f6d397a
md"# Rich display"

# ╔═╡ 3be84600-b166-11ea-1d24-59212363543f
md"## `text/plain`"

# ╔═╡ 42f0a872-b166-11ea-0c71-355d62f67fca
ra = 1:100

# ╔═╡ 794bc212-b166-11ea-0840-fddb29190841
1:13

# ╔═╡ 95898660-b166-11ea-1db1-df7f3c4f1353
"I am not bold"

# ╔═╡ 2859a890-b161-11ea-14e9-b7ddaf08195a
md"## Tree view"

# ╔═╡ 23f41dd2-b15c-11ea-17d2-45b3e83093ba
Ref(Dict(:a => [1,md"![](https://upload.wikimedia.org/wikipedia/commons/thumb/1/1c/PDS_70.jpg/567px-PDS_70.jpg)", md"# Hello"], [3,4] => :b))

# ╔═╡ 88bd7aae-b15f-11ea-270e-ab00e6a01203
["asdf", "I am not bold"]

# ╔═╡ 479df620-b161-11ea-0d2b-293fac4e46bb
md"## Markdown"

# ╔═╡ 57a82490-b162-11ea-0fdb-73c57b08830f
md"Writing a notebook is not just about writing the final document — Pluto empowers the experiments and discoveries that are essential to getting there.

Writing a notebook is not just about writing the final document — Pluto empowers the experiments and discoveries that are essential to getting there."

# ╔═╡ 5dbbc6c0-b162-11ea-0890-7f0e994b2691
md"Writing a notebook is not just about writing the final document — Pluto empowers the experiments and discoveries that are essential to getting there."

# ╔═╡ 1bb05fc0-b15d-11ea-3dae-7734f66a0c56
md"# Testing machinery"

# ╔═╡ 9ac925d0-b15d-11ea-2abd-7db360900be0
html"""

"""

# ╔═╡ 7e2cc6c0-b15d-11ea-32b0-15394cdebd35
function ask(kind, str::Markdown.MD)
	HTML("
" * sprint(show, MIME"text/html"(), str) * "
")
end

# ╔═╡ ce7f9d90-b163-11ea-0ff7-cf7378644741
ask("visual", md"During launch: 👆👆 Header should be BLUE and say LOADING")

# ╔═╡ 11f3ae90-b164-11ea-1027-0d2e6e7048dd
ask("visual", md"Hover shows path outlined and RENAME button")

# ╔═╡ 26d2b310-b164-11ea-0029-f7b04ee5ba73
ask("manual", md"Autocomplete shows on every keystroke")

# ╔═╡ 1c6229f0-b165-11ea-0f1d-674022c43971
ask("manual", md"Autocompleting a **directory** using ENTER does what you expect")

# ╔═╡ 427e4980-b164-11ea-3604-69e5af4ad390
ask("manual", md"Delete `jl` from the end, it gives **one** suggestion: the original path, with a sheet-of-paper icon")

# ╔═╡ 6e04dd80-b164-11ea-0b33-7dc165574a35
ask("manual", md"Delete `?.jl` (4 charachters) from the end, it gives **two** suggestions:
- the original path, with a sheet-of-paper icon
- a new suggested path, which ends with `.jl (new)` and has a different icon

1. Click RENAME
1. _There is a confirmation pop-up_
1. Click CANCEL
1. _The path changes back to the original_
1. Change the name again
1. Press `ENTER`
1. _There is a confirmation pop-up_
1. Click OK
1. _Path remains the same_
1. _The notebook is LOADING for a second_
")

# ╔═╡ 5fd66630-b15e-11ea-29b3-5fdf839d105b
ask("visual", md"Fold symbol should follow page scroll for long cells")

# ╔═╡ 9b26e0c0-b15e-11ea-1fbb-a34d46c520fa
ask("manual", md"Swap names `ma` and `mb`, then use `Ctrl+S` to run both cells. Then again, now use 'Submit all changes' button at top of notebook")

# ╔═╡ f4015940-b15f-11ea-0f3a-1714c79023c3
ask("manual", md"Move this cell:
- before itself (should not change pos)
- after itself (should not change pos)
- somewhere else

Drag to end of the screen, **should autoscroll**")

# ╔═╡ 1b569b72-b167-11ea-1462-63674f7e13db
ask("manual", md"Test autocomplete")

# ╔═╡ 4902abb2-b160-11ea-1120-f501bf151fc2
ask("manual", md"Refresh quickly & slowly a couple of times")

# ╔═╡ 58721630-b160-11ea-30b4-e9667b4df6da
ask("manual", md"Change UUID in URL to something invalid - should redirect to `/`")

# ╔═╡ c8595660-b166-11ea-0b6f-31d05ec9f51f
ask("manual", md"Click on the error")

# ╔═╡ e4a7c0b0-b15f-11ea-2eff-ad120d8cb859
ask("manual", md"Click both yellow items in the stack trace, the relevant line should be selected and scrolled into view")

# ╔═╡ ba734f30-b15f-11ea-0490-a55c3739e4bf
ask("visual", md"Stack trace should have **5** frames and they should look nice

The last **2** should have yellow 'filenames'")

# ╔═╡ a4f07420-b160-11ea-2c62-9115046e4e24
ask("manual broken", md"The `eg` cell should get unfolded when the stack frame is clicked")

# ╔═╡ 5f37b100-b165-11ea-2f65-158f1392aea9
ask("manual", md"Move the slider, the second cell should **not** execute on intermediate values. It should **not** be executing 500ms after releasing the slider.")

# ╔═╡ b3047b10-b165-11ea-014e-25d14d54558c
ask("manual", md"1. Move the slider
1. Uncomment the second cell and run
1. _The second cell has the slider value_
1. Move the slider
1. _Second cell updates_")

# ╔═╡ 0d1ca080-b163-11ea-2927-7b46d602f837
ask("manual", md"Moving slider should **not scroll** window")

# ╔═╡ 0c499910-b163-11ea-3fee-c5d006a97ff7
ask("manual", md"Moving slider should **scroll** window to **keep this message at a fixed position**")

# ╔═╡ 486c7770-b166-11ea-22df-b38b69fb51ad
ask("visual", md"Shows assignee `ra`")

# ╔═╡ 5f28c770-b166-11ea-3099-afb5ec07119b
ask("visual", md"Output can be scrolled horizontally")

# ╔═╡ 80e7b1f0-b166-11ea-1434-ff74883f762d
ask("visual", md"No scrollbar")

# ╔═╡ 893d7bf2-b166-11ea-1580-63694795488b
ask("visual", md"HTML tags are visible")

# ╔═╡ de56f5d2-b15c-11ea-320b-31b4da177c3a
ask("manual", md"I must be an expandable tree")

# ╔═╡ f3e429be-b15e-11ea-2a93-4516283aa4c0
ask("visual", md"Element alignment should be on baseline")

# ╔═╡ 4096de22-b15f-11ea-0178-17a2ac624e96
ask("visual broken", md"HTML should be escaped")

# ╔═╡ 9df64d10-b161-11ea-2c8c-ff6f0b978cc2
ask("visual", md"Tests visibile inside the following Markdown")

# ╔═╡ 78ef6e12-b162-11ea-1474-b3e749a7b705
ask("visual", md"Assignment to `低调又牛逼的类型系统` must be visible")

# ╔═╡ e16977de-b15a-11ea-3cbd-35af9445d31d
低调又牛逼的类型系统 = md"""
 $(ask("visual", md"Heading in simplified Chinese:"))

#### 低调又牛逼的类型系统

 $(ask("visual", md"Code and LaTeX math inline"))

 比如，`sin.(A)` 会计算 A 中每个元素的 $\sin(A)$ 值

 $(ask("visual", md"Bigger and centered LaTeX math"))

$$\frac{\pi^2}{2}$$

 $(ask("visual broken", md"Heading & text in Arabic displayed right-to-left"))

## الكَوكَب

يظن بأن كوكب الأرض هو مركز الكون

 $(ask("visual", md"Greek"))

Πλθτο ισ εεν πλανεετ! Some maths:

 $(ask("visual broken", md"Identical symbols:"))

$$\pi = π$$

 $(ask("visual", md"Image **should not** fill page when widescreen, image **should** fill page when narrow screen"))

![](https://upload.wikimedia.org/wikipedia/commons/thumb/1/1c/PDS_70.jpg/567px-PDS_70.jpg)"""

# ╔═╡ 21bd9950-b15d-11ea-2632-41b1c66563bd
ask("visual", md"These three paragraphs must have equal spacing between them")

# ╔═╡ Cell order:
# ╟─ce7f9d90-b163-11ea-0ff7-cf7378644741
# ╠═878a4750-b15e-11ea-2584-8feba490699f
# ╟─fd0763a0-b163-11ea-23b4-a7bae7052e19
# ╟─11f3ae90-b164-11ea-1027-0d2e6e7048dd
# ╟─26d2b310-b164-11ea-0029-f7b04ee5ba73
# ╟─1c6229f0-b165-11ea-0f1d-674022c43971
# ╟─427e4980-b164-11ea-3604-69e5af4ad390
# ╟─6e04dd80-b164-11ea-0b33-7dc165574a35
# ╟─592d3ac0-b15e-11ea-313c-db05bc7303cc
# ╟─5fd66630-b15e-11ea-29b3-5fdf839d105b
# ╟─6dde0352-b15e-11ea-2fa8-7327cc366c1a
# ╟─9b26e0c0-b15e-11ea-1fbb-a34d46c520fa
# ╠═7370dcc0-b15e-11ea-234b-23584c864b61
# ╠═75b21a30-b15e-11ea-3046-2170ec097e63
# ╟─7b74dd40-b15e-11ea-291a-d7e10a185718
# ╟─f0b821b0-b15f-11ea-1f64-dd33aa85b54e
# ╟─f4015940-b15f-11ea-0f3a-1714c79023c3
# ╟─4980fc10-b163-11ea-081b-c1335699a8f6
# ╠═4e0c5090-b163-11ea-3cf9-a1b17dd2f053
# ╟─d2c1d090-b162-11ea-0c17-2b234c098cf9
# ╟─d890a190-b162-11ea-31dd-8d603787e5c5
# ╟─1b569b72-b167-11ea-1462-63674f7e13db
# ╠═e141f910-b162-11ea-039b-3ba1414cbd07
# ╟─39724bae-b160-11ea-0af5-737755a57e48
# ╟─4902abb2-b160-11ea-1120-f501bf151fc2
# ╟─58721630-b160-11ea-30b4-e9667b4df6da
# ╠═f33dc5b0-b160-11ea-15e3-1b54697b5872
# ╟─95492630-b15f-11ea-3f16-410fe87cd32c
# ╟─c1c894f0-b166-11ea-31e0-dd1cddb7b7ca
# ╟─c8595660-b166-11ea-0b6f-31d05ec9f51f
# ╠═d05d7df0-b166-11ea-2d8f-39b31a42890d
# ╟─bdeba8e0-b166-11ea-3913-3138ab428083
# ╟─e4a7c0b0-b15f-11ea-2eff-ad120d8cb859
# ╟─ba734f30-b15f-11ea-0490-a55c3739e4bf
# ╠═aab109c0-b15f-11ea-275d-31e21fcda8c4
# ╟─84888e20-b160-11ea-1d61-c5934251d6dd
# ╠═9dc4a0a0-b15f-11ea-361c-87742cf3f2a2
# ╟─a4f07420-b160-11ea-2c62-9115046e4e24
# ╟─976bc2a0-b160-11ea-3e7a-9f033b0f2daf
# ╠═9c74f9b2-b160-11ea-35fb-917cb1120f5b
# ╟─ea3f77f0-b166-11ea-046e-ef39bfc57d0f
# ╠═ee92a570-b166-11ea-38fa-0f3563b3f790
# ╟─4f32c16e-b15e-11ea-314e-01700baeb92d
# ╟─3a14b3f0-b165-11ea-153d-796416ee5ccc
# ╟─5f37b100-b165-11ea-2f65-158f1392aea9
# ╠═41a75500-b165-11ea-2519-bbd0feaef6cf
# ╠═4ccbf670-b165-11ea-1951-c17ffb8a58cf
# ╠═8bb26902-b165-11ea-048c-d7f7a72006ee
# ╟─b3047b10-b165-11ea-014e-25d14d54558c
# ╠═e559eaf0-b165-11ea-0d81-ffc480afe8f3
# ╠═e63be680-b165-11ea-0fd3-bd4e0bf92eb8
# ╟─59966a90-b163-11ea-1786-e56e45f06dd0
# ╟─55bade10-b163-11ea-327f-132c93471713
# ╟─0d1ca080-b163-11ea-2927-7b46d602f837
# ╠═f2c0bb90-b162-11ea-24a1-3f864a09e5ee
# ╠═56e6f440-b15e-11ea-1327-09932af5b5bd
# ╠═20d72230-b163-11ea-39c2-69bf2c422d50
# ╟─0c499910-b163-11ea-3fee-c5d006a97ff7
# ╠═2296ac80-b163-11ea-3d00-ed366fa9ce3e
# ╠═54ccad80-b163-11ea-1777-c5364408948e
# ╟─32b5edc0-b15d-11ea-09d6-3b889f6d397a
# ╟─3be84600-b166-11ea-1d24-59212363543f
# ╟─486c7770-b166-11ea-22df-b38b69fb51ad
# ╟─5f28c770-b166-11ea-3099-afb5ec07119b
# ╠═42f0a872-b166-11ea-0c71-355d62f67fca
# ╟─80e7b1f0-b166-11ea-1434-ff74883f762d
# ╠═794bc212-b166-11ea-0840-fddb29190841
# ╟─893d7bf2-b166-11ea-1580-63694795488b
# ╟─95898660-b166-11ea-1db1-df7f3c4f1353
# ╟─2859a890-b161-11ea-14e9-b7ddaf08195a
# ╟─de56f5d2-b15c-11ea-320b-31b4da177c3a
# ╟─f3e429be-b15e-11ea-2a93-4516283aa4c0
# ╟─23f41dd2-b15c-11ea-17d2-45b3e83093ba
# ╟─4096de22-b15f-11ea-0178-17a2ac624e96
# ╠═88bd7aae-b15f-11ea-270e-ab00e6a01203
# ╟─479df620-b161-11ea-0d2b-293fac4e46bb
# ╟─9df64d10-b161-11ea-2c8c-ff6f0b978cc2
# ╟─78ef6e12-b162-11ea-1474-b3e749a7b705
# ╟─e16977de-b15a-11ea-3cbd-35af9445d31d
# ╟─21bd9950-b15d-11ea-2632-41b1c66563bd
# ╟─57a82490-b162-11ea-0fdb-73c57b08830f
# ╟─5dbbc6c0-b162-11ea-0890-7f0e994b2691
# ╟─1bb05fc0-b15d-11ea-3dae-7734f66a0c56
# ╠═9ac925d0-b15d-11ea-2abd-7db360900be0
# ╠═7e2cc6c0-b15d-11ea-32b0-15394cdebd35


text: #!/usr/bin/env julia

max_stars_repo_name: xijiang/Refactered-NSG

max_issues_repo_name: xijiang/Refactered-NSG

max_forks_repo_name: xijiang/Refactered-NSG


text:
"""
# LinuxTouchEvents.jl

Read touch screen events from a
[Linux Input Event](https://www.kernel.org/doc/html/v4.14/input/input.html)
device.

Tested with the [Raspberry Pi Touch Display](https://www.raspberrypi.org/documentation/hardware/display/README.md).

Use this module to get touch input in terminal-mode or frame-buffer apps.

(Not intended for apps that use a high-level API like GTK/X11 or SDL.)
"""
module LinuxTouchEvents

export TouchEventChannel

using ReadmeDocs
using UnixIO
using UnixIO: C
using UnixIO.Debug

README"## Interface"

README"""
    TouchEventChannel([dev="/dev/input/event0"])

Open a channel to read touch events from a Linux Input Event device.

e.g.

    c = TouchEventChannel()
    while true
        x, y = take!(c)
        println("Touch at (\$x, \$y)")
    end
"""
struct TouchEventChannel
    io::UnixIO.FD
    width::Int
    height::Int
    function TouchEventChannel(dev="/dev/input/event0")
        size = read("/sys/class/graphics/fb0/virtual_size", String)
        new(UnixIO.open(dev, C.O_RDONLY),
            (parse(Int, x) for x in split(size, ","))...)
    end
end

# https://www.kernel.org/doc/html/v4.14/input/input.html#event-interface
struct input_event
    tv_sec::Culong
    tv_usec::Culong
    type::Cushort
    code::Cushort
    value::Cuint
end

read_event(io) = reinterpret(input_event, read(io, sizeof(input_event)))[1]

@db function wait_for_event(io, type, code, value=nothing)
    while true
        e = read_event(io)                                             ;@db 4 e
        if e.type == type &&
           e.code == code &&
           (e.value == value || value == nothing)
            @db return e.value
        end
    end
end

# https://www.kernel.org/doc/html/v4.14/input/event-codes.html#event-types
# <linux/input-event-codes.h> https://git.io/JInkG
const EV_KEY = 0x01
const BTN_TOUCH = 0x14a

const EV_ABS = 0x03
const ABS_X  = 0x00
const ABS_Y  = 0x01

@db function Base.take!(t::TouchEventChannel)
    wait_for_event(t.io, EV_KEY, BTN_TOUCH, 1)
    x = wait_for_event(t.io, EV_ABS, ABS_X)
    y = wait_for_event(t.io, EV_ABS, ABS_Y)
    @db return x/t.width, y/t.height
end

Base.close(t::TouchEventChannel) = close(t.io)

end # module


text:
# To speed up the development workflow, use `TestEnv`.
# For example:
# ```
# $ julia --project
#
# julia> ENV["TEST_MLJBASE"] = "true"
#
# julia> using TestEnv; TestEnv.activate()
#
# julia> include("test/preliminaries.jl")
# [...]
#
# julia> include("test/resampling.jl")
# [...]
# ```
include("preliminaries.jl")

@conditional_testset "misc" begin
    @test include("utilities.jl")
    @test include("static.jl")
end

@conditional_testset "interface" begin
    @test include("interface/interface.jl")
    @test include("interface/data_utils.jl")
end

@conditional_testset "measures" begin
    @test include("measures/measures.jl")
    @test include("measures/measure_search.jl")
    @test include("measures/doc_strings.jl")
end

@conditional_testset "resampling" begin
    @test include("resampling.jl")
end

@conditional_testset "data" begin
    @test include("data/data.jl")
    @test include("data/datasets.jl")
    @test include("data/datasets_synthetic.jl")
end

@conditional_testset "sources" begin
    @test include("sources.jl")
end

@conditional_testset "machines" begin
    @test include("machines.jl")
end

@conditional_testset "composition_learning_networks" begin
    @test include("composition/learning_networks/nodes.jl")
    @test include("composition/learning_networks/inspection.jl")
    @test include("composition/learning_networks/machines.jl")
end

@conditional_testset "composition_models" begin
    @test include("composition/models/methods.jl")
    @test include("composition/models/from_network.jl")
    @test include("composition/models/inspection.jl")
    @test include("composition/models/deprecated.jl")
    @test include("composition/models/pipelines.jl")
    @test include("composition/models/transformed_target_model.jl")
    @test include("composition/models/stacking.jl")
    @test include("composition/models/static_transformers.jl")
end

@conditional_testset "operations" begin
    @test include("operations.jl")
end

@conditional_testset "hyperparam" begin
    @test include("hyperparam/one_dimensional_ranges.jl")
    @test include("hyperparam/one_dimensional_range_methods.jl")
end


text:
module GShock
using Random, DataFrames, XLSX, LinearAlgebra, Statistics,
	StatsBase, Distributions, Plots, CSV, RCall;

################################################################################
#  [1] STRUCTURES FOR CONTAINERS
################################################################################
include("GSDataLoader.jl")
include("mainsup.jl")
# ------------------------------------------------------
# [1.1] Estimated Parameters
struct Param
    B::Array   # Estimated βs
    Σ::Array   # Variance Covariance matrix
    Γ::Array   # Median of identificacion matrix
    ξ::Array   # Median of Eigenvectors
end
# ------------------------------------------------------
# [1.2] Post Estimation
struct Post
    Mean
    Qntls
end
# ------------------------------------------------------
# [1.3] Output structure
struct GSsolve
    Par	      # Estimated parameters
    IrfGS     # Impulse response function
    FevGS     # Forecast Error Variances decomposition
	IrfNF     # Impulse response function Non fundamental
    FevNF     # Forecast Error Variances decomposition Non fundamental
end

################################################################################
#  [2] MAIN FUNCTIONS
################################################################################

# ------------------------------------------------------
# [2.1] Estimation of Global Shocks
# ------------------------------------------------------
function GSstimation(y, p, h; xblock = false, GOS = true, nx = 0, VarGS = 1:3, nmodls = 5000,
		Lτ = 1, Uτ = 5, quint = [0.16 0.50 0.84], NF = true)
    # --------------------------------------------------------
    # [2.1.1] VAR in reduce form
    m = size(y)[2]
    Y, X = VARData(y, p)
    B, Φ = OLSbetas(Y, X, m, p, Xblock = xblock, Nx = nx)

    # --------------------------------------------------------
    #  Stability Checking
    if maximum(abs.(eigvals(Φ))) >= 1
        println("Not stationary model, we report only stationary draws")
    end

    # --------------------------------------------------------
    # [2.1.2] Residual of the reduce form
    T, mn = size(X)
    E = Y - X * B
    Σ = (E' * E) / (T - mn)
    SS = (E' * E)

    # --------------------------------------------------------
    # [2.1.3] Global Shock identification
    display("Model is being estimated...")
    Γ     = Array{Float64,2}(undef, m, nmodls) # All the identification vectors
    aux1  = (xblock * nx) + (1 - xblock) * m
    ξ 	  = Array{Float64,3}(undef, aux1, aux1, nmodls)
    IRFGS = Array{Float64,3}(undef, m, h, nmodls) # Container for IRF of GS
    FEVGS = Array{Float64,3}(undef, m, h, nmodls) # Container for IRF of GS
	IRFNF = Array{Float64,3}(undef, m, h, nmodls) # Container for IRF of GS
    FEVNF = Array{Float64,3}(undef, m, h, nmodls) # Container for IRF of GS
	U     = Array{Float64,2}(undef, T, nmodls) # Identified shocks
    for i = 1:nmodls
		β, Φ, γ, ξ[:, :, i] =
    	GSsimulation(Y, X, B, SS, p, Lτ, Uτ, Xblock = xblock, nx = nx, varGS = VarGS, nonfun = NF);
        # β, Φ, Γ, ξ
		Uaux   = γ\(E')
		U[:,i] = vec(Uaux[1,:])
        IRF, FEV = irf_fevd(Φ, γ, h, m)
        IRFGS[:, :, i] = IRF[:, 1, :]
        FEVGS[:, :, i] = FEV[:, 1, :]
		IRFNF[:, :, i] = IRF[:, 2, :]
		FEVNF[:, :, i] = FEV[:, 2, :]
        Γ[:, i] = γ[:, 1]
    end

    # --------------------------------------------------------
    # [2.1.4] Percentiles : Julia does not have a quintile
    # command that work in RN arrays
	# -------- Means --------
    irf1m = dropdims(mean(IRFGS, dims = 3), dims = 3)
    fev1m = dropdims(mean(FEVGS, dims = 3), dims = 3)
	irf2m = dropdims(mean(IRFNF, dims = 3), dims = 3)
    fev2m = dropdims(mean(FEVNF, dims = 3), dims = 3)

	# -------- Percentiles --------
    irf1q = Qntls(IRFGS, nmodls, quint, m, h)
    fev1q = Qntls(FEVGS, nmodls, quint, m, h)
	irf2q = Qntls(IRFNF, nmodls, quint, m, h)
	fev2q = Qntls(FEVNF, nmodls, quint, m, h)

	# -------------------------
    Γ = median(Γ, dims = 2)
    ξ = median(ξ, dims = 3)
    display("It ends!")

	# ------------------------
	# Outputs
	mytup =(irfgs = IRFGS, fevgs= FEVGS,irfnf = IRFNF,fevnf=FEVNF);
	params = Param(B, Σ, Γ, ξ)
	post1  = Post(irf1m, irf1q)
	post2  = Post(fev1m, fev1q)
	post3  = Post(irf2m, irf2q)
	post4  = Post(fev2m, fev2q)
	gsolve = GSsolve(params, post1, post2, post3, post4)
    return gsolve, mean(U,dims=2) , mytup
end

# ------------------------------------------------------
# [2.2] Bayesian Simulation
# ------------------------------------------------------

function GSsimulation(Y, X, B, SS, p, Lτ, Uτ; Xblock = false, nx = 3, varGS = 1:3, nonfun = true)
    # --------------------------------------------------------
    # [2.2.1] Drawing βₙ
    T, m   = size(Y)
    σdist  = InverseWishart(T, SS)
    b      = vec(B)
    rc, cc = size(B)
    Σ      = rand(σdist)                   # Sampling from an Inverse  Wishart
    XX     = Symmetric(inv(X'X), :L)
    kr     = kron(Σ, XX)
    disti  = MvNormal(b, kr)
    β, Φ   = BetaDraw(disti, rc, cc, m, p)
    # --------------------------------------------------------
    # [2.2.2] Cholesky Identification
    E1 = Y - X * β
    Σ1 = (E1' * E1) / (T - size(X, 2))
    C1 = cholesky(Σ1).L

    # --------------------------------------------------------
    # [2.2.3] Impulse-Response Function :
    #   The IRF[i,j,k] is the response k-period ahead of
    #   the variable i when is hit by the shock j;
    aux_IRF = Array{Float64,3}(undef, m * p, m * p, Uτ)
    IRF 	= Array{Float64,3}(undef, m, m, Uτ)
    aux_IRF[:, :, 1] = Φ
    IRF[:, :, 1]   	 = aux_IRF[1:m, 1:m, 1] * C1

    for i = 2:Uτ
        aux_IRF[:, :, i] = Φ * aux_IRF[:, :, i-1]
        IRF[:, :, i]     = aux_IRF[1:m, 1:m, i] * C1
    end

    # --------------------------------------------------------
    # [2.2.4] Global Shocks Identification
    # [Step 1] Matrix Λ and their weights
    tvars = length(varGS)
	naux  = (Xblock * nx) + (1 - Xblock) * m;
    Λ = zeros(naux, naux,  tvars)
    λ = zeros(tvars)
    for i = 1:tvars
        Λ[:, :, i] = Λmatrix(IRF, varGS[i], m; Lτ = Lτ, Uτ = Uτ, Xblock = Xblock, nx = nx)
        λ[i] 	   = tr(Λ[:, :, i])
    end
    λ = prod(λ) ./ λ

    # [Step 2] Matrix ξ
    Ξ = dropdims(sum(reshape(λ, 1, 1, :) .* Λ, dims = 3), dims = 3)
    # [Step 3] Eigenvector sorted by the max eigenvalue
    ξ = eigen(Ξ).vectors[:, end:-1:1]
    # [Step 4] NON FUNDAMENTAL SHOCK
	if nonfun
		nξ  = size(ξ)[1];
	    L   = [-ξ[2:end,1] ./ ξ[1,1] I(nξ-1)];
	    BL  = L*L';
	    BL  = Array(BL);
	    pos_tot=2;
	    Λaux= L*Λ[:,:,pos_tot]*L';
	    ψ   = eigen(Λaux, BL).vectors[:,end];
	    Ψ   = convert(Array,(ψ'*L)');
	    Ξ   = [ξ[:,1] Ψ];
	    # New Identification matrix
	    nsp = nullspace(Array(Ξ'))
	    ξ   = [Ξ nsp[:,1:(nξ-2)]];
	end
    # [Step 5] New Identification matrix
    if Xblock
        Γ = [C1[:, 1:nx] * ξ C1[:, nx+1:end]]
    else
        Γ = C1 * ξ
    end
    Γ = sign.(diag(Γ))' .* Γ  # Normalization to improve identification
    return β, Φ, Γ, ξ # betas, Companion Form, Identification, Eigenvectors
end
end

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
@taylorize function lotka_volterra!(du, u, p, t)
    u1u2 = u[1] * u[2]
    du[1] = 3.0 * (u[1] - u1u2)
    du[2] = u1u2 - u[2]
    return du
end

function lotka_volterra_hybrid(; nsplit=1,
                                 ε = 0.008,
                                 ε_ext=1e-4, # threshold for the outer approximation
                                 n_int=50)   # number of directions for the inner approximation

    # generate external / internal polytopic approximations of the guard
    B = Ball2([1.0, 1.0], 0.15) # "exact"
    B_ext = overapproximate(B, ε_ext) # outer approximation
    B_int = underapproximate(B, PolarDirections(n_int)) # inner approximation
    B_int = tohrep(convert(VPolygon, B_int)) # cast to Hrep
    B_intᶜ = complement(B_int)

    # define modes
    aut = LightAutomaton(3)
    outside = @system(x' = lotka_volterra!(x), dim: 2, x ∈ B_intᶜ)
    inside = @system(x' = lotka_volterra!(x), dim: 2, x ∈ B_ext)
    outside_unconstrained = @system(x' = lotka_volterra!(x), dim: 2, x ∈ Universe(2))

    # define the transition graph
    add_transition!(aut, 1, 2, 1)
    add_transition!(aut, 2, 3, 2)
    T_out_in = @map(x -> x, dim:2, x ∈ B_ext)
    T_in_out = @map(x -> x, dim:2, x ∈ B_intᶜ)

    # initial-value problem
    H = HybridSystem(automaton=aut, modes=[outside, inside, outside_unconstrained],
                                           resetmaps=[T_out_in, T_in_out])

    # initial states with splitting
    X0 = Hyperrectangle(low=[1.3-ε, 1.], high=[1.3+ε, 1.])
    X0s = split(X0, [nsplit, 1])
    X0st = [(X0s_i, 1) for X0s_i in X0s]
    return InitialValueProblem(H, X0st)
end

#=
# TEST RUN:

prob = lotka_volterra_hybrid(nsplit=10, ε_ext=1e-4, n_int = 30, ε = 0.008);
@time sol = solve(prob,
                  tspan=(0.0, 3.64),
                  alg=TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=true, disjointness=RA.ZonotopeEnclosure()),
                  max_jumps=2,
                  intersect_source_invariant=false,
                  intersection_method=RA.BoxIntersection(),
                  clustering_method=RA.BoxClustering(),
                  disjointness_method=RA.BoxEnclosure());
solz = overapproximate(sol, Zonotope);

B = Ball2([1.0, 1.0], 0.15) # "exact"
B_ext = overapproximate(B, 1e-6) # outer approximation
plot(solz, vars=(1, 2))
plot!(B_ext, ratio=1)
=#


text:
module Heterogeneous

end # module


text:
# ClimaCore hacks
Base.isless(zp1::CC.Geometry.ZPoint, zp2::CC.Geometry.ZPoint) = Base.isless(zp1.z, zp2.z)
Base.isless(zp::CC.Geometry.ZPoint, val::Number) = Base.isless(zp.z, val)
Base.isless(val::Number, zp::CC.Geometry.ZPoint) = Base.isless(val, zp.z)

Base.:+(zp::CC.Geometry.ZPoint, val) = Base.:+(zp.z, val)
Base.:+(val, zp::CC.Geometry.ZPoint) = Base.:+(val, zp.z)
Base.:/(val, zp::CC.Geometry.ZPoint) = Base.:/(val, zp.z)
Base.:+(zp1::CC.Geometry.ZPoint, zp2::CC.Geometry.ZPoint) = Base.:+(zp1.z, zp2.z)

Base.:*(zp1::CC.Geometry.ZPoint, zp2::CC.Geometry.ZPoint) = Base.:*(zp1.z, zp2.z)

Base.log(zp::CC.Geometry.ZPoint) = Base.log(zp.z)
Base.:-(zp::CC.Geometry.ZPoint) = Base.:-(zp.z)
Base.:-(zp::CC.Geometry.ZPoint, val) = Base.:-(zp.z, val)
Base.:-(val, zp::CC.Geometry.ZPoint) = Base.:-(val, zp.z)
Base.:-(zp1::CC.Geometry.ZPoint, zp2::CC.Geometry.ZPoint) = Base.:-(zp1.z, zp2.z)

Base.convert(::Type{Float64}, zp::CC.Geometry.ZPoint) = zp.z

struct Grid{FT, CS, FS, SC, SF}
    zmin::FT
    zmax::FT
    Δz::FT
    Δzi::FT
    nz::Int
    cs::CS
    fs::FS
    zc::SC
    zf::SF
    function Grid(Δz::FT, nz::Int) where {FT <: AbstractFloat}
        z₀, z₁ = FT(0), FT(nz * Δz)

        domain = CC.Domains.IntervalDomain(
            CC.Geometry.ZPoint{FT}(z₀),
            CC.Geometry.ZPoint{FT}(z₁),
            boundary_tags = (:bottom, :top),
        )

        mesh = CC.Meshes.IntervalMesh(domain, nelems = nz)

        cs = CC.Spaces.CenterFiniteDifferenceSpace(mesh)
        fs = CC.Spaces.FaceFiniteDifferenceSpace(cs)
        zc = CC.Fields.coordinate_field(cs)
        zf = CC.Fields.coordinate_field(fs)

        #Set the inverse grid spacing
        Δzi = 1 / Δz

        zmin = minimum(parent(zf))
        zmax = maximum(parent(zf))
        CS = typeof(cs)
        FS = typeof(fs)
        SC = typeof(zc)
        SF = typeof(zf)
        return new{FT, CS, FS, SC, SF}(zmin, zmax, Δz, Δzi, nz, cs, fs, zc, zf)
    end
end

# Index of the first interior cell above the surface
kc_surface(grid::Grid) = Cent(1)
kf_surface(grid::Grid) = CCO.PlusHalf(1)
kc_top_of_atmos(grid::Grid) = Cent(grid.nz)
kf_top_of_atmos(grid::Grid) = CCO.PlusHalf(grid.nz + 1)

is_surface_center(grid::Grid, k) = k == kc_surface(grid)
is_toa_center(grid::Grid, k) = k == kc_top_of_atmos(grid)
is_surface_face(grid::Grid, k) = k == kf_surface(grid)
is_toa_face(grid::Grid, k) = k == kf_top_of_atmos(grid)

zc_surface(grid::Grid) = grid.zc[kc_surface(grid)]
zf_surface(grid::Grid) = grid.zf[kf_surface(grid)]
zc_toa(grid::Grid) = grid.zc[kc_top_of_atmos(grid)]
zf_toa(grid::Grid) = grid.zf[kf_top_of_atmos(grid)]

real_center_indices(grid::Grid) = CenterIndices(grid)
real_face_indices(grid::Grid) = FaceIndices(grid)

struct FaceIndices{Nstart, Nstop, G}
    grid::G
    function FaceIndices(grid::G) where {G <: Grid}
        Nstart, Nstop = kf_surface(grid).i, kf_top_of_atmos(grid).i
        new{Nstart, Nstop, G}(grid)
    end
end

struct CenterIndices{Nstart, Nstop, G}
    grid::G
    function CenterIndices(grid::G) where {G <: Grid}
        Nstart, Nstop = kc_surface(grid).i, kc_top_of_atmos(grid).i
        new{Nstart, Nstop, G}(grid)
    end
end

Base.keys(ci::CenterIndices) = 1:length(ci)
Base.keys(fi::FaceIndices) = 1:length(fi)

n_start(::CenterIndices{Nstart}) where {Nstart} = Nstart
n_start(::FaceIndices{Nstart}) where {Nstart} = Nstart
n_stop(::CenterIndices{Nstart, Nstop}) where {Nstart, Nstop} = Nstop
n_stop(::FaceIndices{Nstart, Nstop}) where {Nstart, Nstop} = Nstop

Base.getindex(ci::CenterIndices, i::Int) = Cent(Base.getindex(n_start(ci):n_stop(ci), i))
Base.getindex(fi::FaceIndices, i::Int) = CCO.PlusHalf(Base.getindex(n_start(fi):n_stop(fi), i))

Base.length(::FaceIndices{Nstart, Nstop}) where {Nstart, Nstop} = Nstop - Nstart + 1
Base.length(::CenterIndices{Nstart, Nstop}) where {Nstart, Nstop} = Nstop - Nstart + 1

Base.iterate(fi::CenterIndices{Nstart, Nstop}, state = Nstart) where {Nstart, Nstop} =
    state > Nstop ? nothing : (Cent(state), state + 1)

Base.iterate(fi::FaceIndices{Nstart, Nstop}, state = Nstart) where {Nstart, Nstop} =
    state > Nstop ? nothing : (CCO.PlusHalf(state), state + 1)

Base.iterate(fi::Base.Iterators.Reverse{T}, state = Nstop) where {Nstart, Nstop, T <: CenterIndices{Nstart, Nstop}} =
    state < Nstart ? nothing : (Cent(state), state - 1)

Base.iterate(fi::Base.Iterators.Reverse{T}, state = Nstop) where {Nstart, Nstop, T <: FaceIndices{Nstart, Nstop}} =
    state < Nstart ? nothing : (CCO.PlusHalf(state), state - 1)

face_space(grid::Grid) = grid.fs
center_space(grid::Grid) = grid.cs

Base.eltype(::Grid{FT}) where {FT} = FT


text:
function projection2D(model::MetNet, ider1::IDER_TYPE, ider2::IDER_TYPE;
        l::Int = 20
    )

    ider1_L, ider1_U = bounds(model, ider1)
    ider1_range = range(ider1_L, ider1_U; length = l)
    ider2_L, ider2_U = bounds(model, ider2)
    ider2_range = range(ider2_L, ider2_U; length = l)

    proj = Dict{Float64, Vector{Float64}}()
    for flx1 in ider1_range
        try
            fixxing(model, ider1, flx1) do
                flx2_lb, flx2_ub = fva(model, ider2) .|> first
                proj[flx1] = [flx2_lb, flx2_ub]
            end
        catch err; end
    end
    return proj
end


text:
include("./engine.jl")

const DEBUG = false

function execute(e::Engine, source, args)
    ops = replace(source, r"[^\>\<\+\-\.\,\[\]\:\;\#\$]" => "")
    jump = 0
    aptr, cptr = 1, 1
    buf = []

    while cptr <= length(ops)
        if ops[cptr] == '>'
            DEBUG && print("inc_ptr\n")
            inc_ptr(e)
        elseif ops[cptr] == '<'
            DEBUG && print("dec_ptr\n")
            dec_ptr(e)
        elseif ops[cptr] == '+'
            DEBUG && print("inc_dat\n")
            inc_dat(e)
        elseif ops[cptr] == '-'
            DEBUG && print("dec_dat\n")
            dec_dat(e)
        elseif ops[cptr] == '.'
            DEBUG && print("output\n")
            push!(buf, output(e))
        elseif ops[cptr] == ','
            DEBUG && print("input\n")
            input(e, args[aptr])
            aptr += 1
        elseif ops[cptr] == ':'
            DEBUG && print("ioutput\n")
            push!(buf, ioutput(e))
        elseif ops[cptr] == ';'
            DEBUG && print("iinput\n")
            iinput(e, args[aptr])
            aptr += 1
        elseif ops[cptr] == '#'
            DEBUG && print("copy_to_top\n")
            copy_to_top(e)
        elseif ops[cptr] == '$'
            DEBUG && print("copy_to_cell\n")
            copy_to_cell(e)
        elseif ops[cptr] == '['
            DEBUG && print("loop\n")
            if !loop(e)
                cptr += findfirst(c -> c==']', ops[cptr:end])-1
            else
                jump = cptr-1
            end
        else
            DEBUG && print("loop end $(jump)\n")
            cptr = jump
        end
        cptr += 1
    end

    return buf
end


text:
include("helpers.jl")

laff_dot(x::T, y::T) where T <: Matrix{Number} = begin

    my, ny = verifydimensions(x, y)

    x = reshape(x, my, ny)

    sum(x .* y)
end

## Testar ##

x = Matrix([1,2,3], 1, 3)
y = Matrix([4,5,6], 3, 1)

laff_dot(x, y)

max_stars_repo_path: LAFF/Semana 1/laff_dot.jl

max_issues_repo_path: LAFF/Semana 1/laff_dot.jl

max_forks_repo_path: LAFF/Semana 1/laff_dot.jl


text:
include("io.jl")
include("line.jl")
include("plot_utils.jl")
include("functions.jl")
include("atmosphere.jl")
include("lambda_iteration.jl")

global nλ_bb = 50
global nλ_bf = 20

function recover_regular(ϵ::AbstractFloat,
                         maxiter::Integer,
                         quadrature::String,
                         DATA::String)

    println("---Recovering simulation---")
    atmos, S_new, populations = read_quantities(DATA; periodic=true)
    line = HydrogenicLine(test_atom(nλ_bb, nλ_bf)..., atmos.temperature)

    # LTE populations
    LTE_pops = LTE_populations(line, atmos)

    # Find continuum extinction and absorption extinction (without Thomson and Rayleigh)
    α_cont = α_absorption.(line.λ0,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(line.λ0,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])

    # destruction probability (Should I include line???)
    ελ = destruction(LTE_pops, atmos.electron_density, atmos.temperature, line)
    println("Minimum $(minimum(ελ)) destruction probability")

    # Start with the source function as the Planck function
    B_0 = Array{Float64, 4}(undef, (length(line.λ), size(α_cont)...))u"kW*m^-2*nm^-1"
    for l in eachindex(line.λ)
        B_0[l, :, :, :] = B_λ.(line.λ[l], atmos.temperature)
    end

    S_old = zero(B_0)

    C = calculate_C(atmos, LTE_pops)

    local convergence
    h5open(DATA, "r") do file
        convergence = read(file, "convergence")[:]
    end

    local i
    for j in 1:length(convergence)
        if convergence[j] == 0
            i = j
            break
        end
    end

    local J_new
    # check where ε < 5e-3, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, DATA)
        #############################
        # Calculate radiation field #
        #############################
        S_old = copy(S_new)
        J_new, damping_λ = J_λ_regular(S_old, α_cont, populations,
                                       atmos, line, quadrature)

        for l in eachindex(line.λ)
            S_new[l,:,:,:] = (1 .- ελ).*J_new[l,:,:,:] .+ ελ.*B_0[l,:,:,:]
        end

        #############################
        #      Calculate rates      #
        #############################
        R = calculate_R(atmos, line, J_new, damping_λ, LTE_pops)

        #############################
        #    Update populations     #
        #############################
        populations = get_revised_populations(R, C, atmos.hydrogen_populations*1.0)

        ############################
        #Write current data to file#
        ############################
        write_to_file(populations[:, 2:end-1, 2:end-1, :], DATA)
        write_to_file(S_new[:, :, 2:end-1, 2:end-1], DATA)

        ############################
        #     Iteration is done    #
        ############################
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont, populations
    end

    println("Converged in $i iterations")

    h5open(DATA, "r+") do file
       file["time"][:] = 0.8888888888
    end

    return J_new, S_new, α_cont, populations
end

function recover_voronoi(ϵ::AbstractFloat,
                         maxiter::Integer,
                         quadrature::String,
                         DATA::String)

    println("---Recovering simulation---")
    sites, S_new, populations = read_sites(DATA)
    line = HydrogenicLine(test_atom(nλ_bb, nλ_bf)..., sites.temperature)

    println(size(S_new))
    println(length(line.λ))

    return

    # Start in LTE
    LTE_pops = LTE_populations(line, sites)

    # Find continuum extinction and absorption extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(line.λ0,
                           sites.temperature,
                           sites.electron_density*1.0,
                           LTE_pops[:,1].+LTE_pops[:,2],
                           LTE_pops[:,3]) .+
             α_scattering.(line.λ0,
                           sites.electron_density,
                           LTE_pops[:,1])

    # Start with the source function as the Planck function
    B_0 = Array{Float64, 2}(undef, (length(line.λ), sites.n))u"kW*m^-2*nm^-1"

    for l in eachindex(line.λ)
        B_0[l, :] = B_λ.(line.λ[l], sites.temperature)
    end

    S_old = zero(B_0)

    # destruction probability (Should I include line???)
    ελ = destruction(LTE_pops, sites.electron_density, sites.temperature, line)
    println("Minimum $(minimum(ελ)) destruction probability")

    C = calculate_C(sites, LTE_pops)

    local convergence
    h5open(DATA, "r") do file
        convergence = read(file, "convergence")[:]
    end

    local i
    for j in 1:length(convergence)
        if convergence[j] == 0
            i = j
            break
        end
    end

    local J_new
    # check where ε < 1e-2, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, DATA)
        #############################
        # Calculate radiation field #
        #############################
        S_old = copy(S_new)
        J_new, damping_λ = J_λ_voronoi(S_old, α_cont, populations,
                                       sites, line, quadrature)

        for l in eachindex(line.λ)
            S_new[l,:] = (1 .- ελ).*J_new[l,:] .+ ελ.*B_0[l,:]
        end

        #############################
        #      Calculate rates      #
        #############################
        R = calculate_R(sites, line, J_new, damping_λ, LTE_pops)

        #############################
        #    Update populations     #
        #############################
        populations = get_revised_populations(R, C, sites.hydrogen_populations)

        ############################
        #Write current data to file#
        ############################
        write_to_file(populations, DATA)
        write_to_file(S_new, DATA)

        ############################
        #     Iteration is done    #
        ############################
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont, populations
    end

    println("Converged in $i iterations")

    h5open(DATA, "r+") do file
       file["time"][:] = 0.8888888888
    end

    return J_new, S_new, α_cont, populations
end

"""
    read_sites(DATA::String)

read quantities from irregular grid simulation from a hdf5 file
"""
function read_sites(DATA::String)
    local positions, temperature, electron_density, hydrogen_populations, S_λ, populations
    local velocity_z, velocity_x, velocity_y, boundaries
    h5open(DATA, "r") do file
        positions = read(file, "positions")[:, :]*u"m"
        boundaries = read(file, "boundaries")[:]u"m"

        temperature = read(file, "temperature")[:]*u"K"

        electron_density = read(file, "electron_density")[:]*u"m^-3"
        hydrogen_populations = read(file, "hydrogen_populations")[:]*u"m^-3"

        velocity_z = read(file, "velocity_z")[:]u"m*s^-1"
        velocity_x = read(file, "velocity_x")[:]u"m*s^-1"
        velocity_y = read(file, "velocity_y")[:]u"m*s^-1"

        S_λ = read(file, "source_function")[:, :]*u"kW*m^-2*nm^-1"
        populations = read(file, "populations")[:, :]*u"m^-3"
    end

    sites_file = "../data/sites_recover.txt"
    neighbours_file = "../data/neighbours_recover.txt"
    # write sites to file
    write_arrays(positions[2, :],
                 positions[3, :],
                 positions[1, :],
                 sites_file)

    x_min = ustrip(boundaries[3])
    x_max = ustrip(boundaries[4])
    y_min = ustrip(boundaries[5])
    y_max = ustrip(boundaries[6])
    z_min = ustrip(boundaries[1])
    z_max = ustrip(boundaries[2])

    # export sites to voro++, and compute grid information
    println("---Preprocessing grid---")

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
                   $(x_min) $(x_max)
                   $(y_min) $(y_max)
                   $(z_min) $(z_max)`)

    n_sites = length(positions[1,:])

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         temperature,
                         electron_density,
                         hydrogen_populations,
                         velocity_z,
                         velocity_x,
                         velocity_y,
                         z_min*1u"m", z_max*1u"m",
                         x_min*1u"m", x_max*1u"m",
                         y_min*1u"m", y_max*1u"m",
                         n_sites)

    return sites, S_λ, populations
end

recover_voronoi(1e-3, 100, "../quadratures/ul7n12.dat",
                "../data/voronoi_ul7n12_extinction_3e5_recover.h5")

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
export BlockShuffle, CycleShuffle, CircShift

#########################################################################
# BlockSuffle
#########################################################################
"""
    BlockShuffle(n::Int) <: Surrogate

A block shuffle surrogate constructed by dividing the time series
into `n` blocks of roughly equal width at random indices (end
blocks are wrapped around to the start of the time series).

Block shuffle surrogates roughly preserve short-range temporal properties
in the time series (e.g. correlations at lags less than the block length),
but break any long-term dynamical information (e.g. correlations beyond
the block length).

Hence, these surrogates can be used to test any null hypothesis aimed at
comparing short-range dynamical properties versus long-range dynamical
properties of the signal.
"""
struct BlockShuffle <: Surrogate
    n::Int
end

Base.show(io::IO, bs::BlockShuffle) = show(io, "BlockShuffle(n=$(bs.n))")

# Split time series in two by default.
BlockShuffle() = BlockShuffle(2)

function get_uniform_blocklengths(L::Int, n::Int)
     # Compute block lengths
    N = floor(Int, L/n)
    R = L % n
    blocklengths = [N for i = 1:n]
    for i = 1:R
        blocklengths[i] += 1
    end

    return blocklengths
end

function surrogenerator(x::AbstractVector, bs::BlockShuffle)
    L = length(x)
    bs.n < L || error("The number of blocks exceeds number of available points")
    Ls = get_uniform_blocklengths(L, bs.n)
    cs = cumsum(Ls)
    # will hold a rotation version of x
    xrot = similar(x)
    T = eltype(xrot)
    init = NamedTuple{(:L, :Ls, :cs, :xrot),Tuple{Int, Vector{Int}, Vector{Int}, Vector{T}}}((L, Ls, cs, xrot))
    return SurrogateGenerator(bs, x, init)
end

function (bs::SurrogateGenerator{<:BlockShuffle})()
    # TODO: A circular custom array implementation would be much more elegant here
    L = bs.init.L
    Ls = bs.init.Ls
    cs = bs.init.cs
    xrot = bs.init.xrot
    n = bs.method.n
    x = bs.x

    # Just create a temporarily randomly shifted array, so we don't need to mess
    # with indexing twice.
    circshift!(xrot, x, rand(1:L))

    # Block always must be shuffled (so ordered samples are not permitted)
    draw_order = zeros(Int, n)
    while any(draw_order .== 0) || all(draw_order .== 1:n)
       StatsBase.sample!(1:n, draw_order, replace = false)
    end

    # The surrogate.
    # TODO: It would be faster to re-allocate, but blocks may
    # be of different sizes and are shifted, so indexing gets messy.
    # Just append for now.
    T = eltype(x)
    s = Vector{T}(undef, 0)
    sizehint!(s, L)

    startinds = [1; cs .+ 1]
    @inbounds for i in draw_order
        inds = startinds[i]:startinds[i]+Ls[i]-1
        append!(s, xrot[inds])
    end

    return s
end

#########################################################################
# CycleShuffle
#########################################################################
"""
    CycleShuffle(n::Int = 7, σ = 0.5) <: Surrogate

Cycle shuffled surrogates[^Theiler1995] that identify successive local peaks in the data and shuffle the
cycles in-between the peaks. Similar to [`BlockShuffle`](@ref), but here
the "blocks" are defined as follows:
1. The timeseries is smoothened via convolution with a Gaussian (`DSP.gaussian(n, σ)`).
2. Local maxima of the smoothened signal define the peaks, and thus the blocks in between them.
3. The first and last index of timeseries can never be peaks and thus signals that
   should have peaks very close to start or end of the timeseries may not perform well. In addition,
   points before the first or after the last peak are never shuffled.
3. The defined blocks are randomly shuffled as in [`BlockShuffle`](@ref).

CSS are used to test the null hypothesis that the signal is generated by a periodic
oscillator with no dynamical correlation between cycles,
i.e. the evolution of cycles is not deterministic.

See also [`PseudoPeriodic`](@ref).

[^Theiler1995]: J. Theiler, On the evidence for low-dimensional chaos in an epileptic electroencephalogram, [Phys. Lett. A 196](https://doi.org/10.1016/0375-9601(94)00856-K)
"""
struct CycleShuffle{T <: AbstractFloat} <: Surrogate
    n::Int
    σ::T
end
CycleShuffle(n = 7, σ = 0.5) = CycleShuffle{typeof(σ)}(n, σ)

function surrogenerator(x::AbstractVector, cs::CycleShuffle)
    n, N = cs.n, length(x)
    g = DSP.gaussian(n, cs.σ)
    smooth = DSP.conv(x, g)
    r = length(smooth) - N
    smooth = iseven(r) ? smooth[r÷2+1:end-r÷2] : smooth[r÷2+1:end-r÷2-1]
    peaks = findall(i -> smooth[i-1] < smooth[i] && smooth[i] > smooth[i+1], 2:N-1)
    blocks = [collect(peaks[i]:peaks[i+1]-1) for i in 1:length(peaks)-1]
    init =  (blocks = blocks, s = copy(x), peak1 = peaks[1])
    SurrogateGenerator(cs, x, init)
end

function (sg::SurrogateGenerator{<:CycleShuffle})()
    blocks, s, peak1 = sg.init
    x = sg.x
    shuffle!(blocks)
    i = peak1
    for b in blocks
        s[(0:length(b)-1) .+ i] .= @view x[b]
        i += length(b)
    end
    return s
end

#########################################################################
# Timeshift
#########################################################################
"""
    CircShift(n) <: Surrogate
Surrogates that are circularly shifted versions of the original timeseries.

`n` can be an integer (meaning to shift for `n` indices), or any vector of integers,
which which means that each surrogate is shifted by an integer,
selected randomly among the entries in `n`.
"""
struct CircShift{N} <: Surrogate
    n::N
end

function surrogenerator(x, sd::CircShift)
    return SurrogateGenerator(sd, x, nothing)
end

function (sg::SurrogateGenerator{<:CircShift})()
    s = random_shift(sg.method.n)
    return circshift(sg.x, s)
end

random_shift(n::Integer) = n
random_shift(n::AbstractVector{<:Integer}) = rand(n)


text:
struct MissingProperty end
const MProperty = MissingProperty()

check_type(x, ptype, ::MissingProperty, default) = _to_type(x, ptype, default)
check_type(x, ptype, val, default) = _to_type(x, ptype, val)

_to_type(x, ptype::Function, val::Function) = __to_type(ptype(x), val(x))
_to_type(x, ptype::Function, val) = __to_type(ptype(x), val)
_to_type(x, ::Type{T}, val::Function) where {T} = __to_type(T, val(x))
_to_type(x, ::Type{T}, val) where {T} = __to_type(T, val)
__to_type(::Type{T}, val) where {T} = val isa T ? val : T(val)::T

function setter_error(x, pname)
    error("$(typeof(x)) does not have `properties` method. Cannot set $(pname) property.")
end

"""
    getter(x, pname, ptype, pdef)

Provides a type stable interface to retreiving properties from the properties of `x`.

* `x::Any`: Any object that may or may not have properties as defined in `ImageCore.jl`.
* `pname::Symbol`: The property name that should be searched for in `x`'s properties.
* `ptype::Union{Function,Type}`: If `pytype` is a `Type` then whatever is returned
  is first converted to this type. If `ptype` is a `Function` then `ptype(x)` is called
  to determine the appropriate return type.
* `pdef::Union{Function,Any}`: If `pdef` determines the defaul that is returned if
  `pname` is not found in the properties of `x` or if `x` has no properties. If `pdef`
  is a function then `pdef(x)` is called to reach the default returned value. Otherwise,
  `pdef` is returned. The default value should be consistent with what is defined by
  `ptype`.

# Examples
```
julia> new_property(x) = getter(x, "new_property", String, x -> "nothing")

julia> p = Dict{Symbol,Any}()

julia> new_property(p)
"nothing"

julia> p["new_property"] = "not nothing"

julia> new_property(p)
"not nothing"
```
"""
@inline getter(x, pname, ptype, pdef) = _pre_getter(HasProperties(x), x, pname, ptype, pdef)
@inline getter(x::AbstractDict, pname, ptype, pdef) = check_type(x, ptype, _getter(x, pname), pdef)
_pre_getter(::HasProperties{true}, x, pname, ptype, pdef) = check_type(x, ptype, _getter(properties(x), pname), pdef)
_pre_getter(::HasProperties{false}, x, pname, ptype, pdef) = check_type(x, ptype, MProperty, pdef)
_getter(x, pname) = get(x, pname, MProperty)

"""
    setter!(x, pname, ptype, val)

Provides a type stable interface to setting the properties the properties of `x`.

* `x::Any`: Any object that may or may not have properties as defined in `ImageCore.jl`.
* `pname::Symbol`: The property name that should be searched for in `x`'s properties.
* `ptype::Union{Function,Type}`: If `pytype` is a `Type` then `val` is ensured to
  be a subtype of this. If `ptype` is a `Function` then `ptype(x)` is called
  to determine the appropriate `Type` the `val` should be a subtype of.
* `val`: The value that the property `pname` is set to.

# Examples
```
julia> new_property!(x, val) = setter!(x, "new_property", String, val)

julia> p = Dict{Symbol,Any}()

julia> new_property(p, "not nothing")

julia> p["new_property"]
"not nothing"
```
"""
@inline setter!(x, pname, ptype, val) = _pre_setter!(HasProperties(x), x, pname, ptype, val)
setter!(x::AbstractDict, pname, ptype, val) = setindex!(x, _to_type(x, ptype, val), pname)
_pre_setter!(::HasProperties{true}, x, pname, ptype, val) = setindex!(properties(x), _to_type(x, ptype, val), pname)
_pre_setter!(::HasProperties{false}, x, pname, ptype, pdef) = setter_error(x, pname)


text:
"""
        grey(decisionMat, weights, fs, zeta)

    Perform GRA (Grey Relational Analysis) for a given decision matrix and weights.

# Arguments:
 - `decisionMat::Array{Float,2}`: n × m matrix of decision matrix.
 - `weights::Array{Float64, 1}`: m-vector of weights for criteria.
 - `fs::Array{Function, 1}`: m-vector of functions that are either maximize or minimize for each single criterion.
 - `zeta::Float64`: zeta parameter for the algorithm. The default is 0.5.

# Description
Applies GRA (Grey Relational Analysis).

# Output
- `::GameResult`: GreyResult object that holds many values including ordering of strategies or candidates and best index.

# Examples
```julia-repl
julia> # Decision matrix of 5 candidates and 6 criteria.
julia> df = DataFrame(
               :K1 => [105000.0, 120000, 150000, 115000, 135000],
               :K2 => [105.0, 110, 120, 105, 115],
               :K3 => [10.0, 15, 12, 20, 15],
               :K4 => [4.0, 4, 3, 4, 5],
               :K5 => [300.0, 500, 550, 600, 400],
               :K6 => [10.0, 8, 12, 9, 9]
        )
5×6 DataFrame
 Row │ K1        K2       K3       K4       K5       K6
     │ Float64   Float64  Float64  Float64  Float64  Float64
─────┼───────────────────────────────────────────────────────
   1 │ 105000.0    105.0     10.0      4.0    300.0     10.0
   2 │ 120000.0    110.0     15.0      4.0    500.0      8.0
   3 │ 150000.0    120.0     12.0      3.0    550.0     12.0
   4 │ 115000.0    105.0     20.0      4.0    600.0      9.0
   5 │ 135000.0    115.0     15.0      5.0    400.0      9.0

julia> # Direction of optimization for each single criterion
julia> functionlist = [minimum, maximum, minimum, maximum, maximum, minimum];

julia> # Weights
julia> w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40];

julia> result = grey(df, w, functionlist);

julia> result.scores
5-element Array{Float64,1}:
 0.525
 0.7007142857142857
 0.5464285714285715
 0.5762820512820512
 0.650952380952381

julia> result.bestIndex
2
```

# References
İşletmeciler, Mühendisler ve Yöneticiler için Operasyonel, Yönetsel ve Stratejik Problemlerin
Çözümünde Çok Kriterli Karar verme Yöntemleri, Editörler: Bahadır Fatih Yıldırım ve Emrah Önder,
Dora, 2. Basım, 2015, ISBN: 978-605-9929-44-8
"""
function grey(decisionMat::DataFrame, weights::Array{Float64,1}, fs::Array{Function,1}, zeta::Float64=0.5)::GreyResult

    #mat = convert(Matrix, decisionMat)
    mat = Matrix(decisionMat)

    nrows, ncols = size(mat)

    w = unitize(weights)

    referenceRow = apply_columns(fs, mat)

    normalizedReferenceRow = zeros(Float64, ncols)

    normalizedMat = similar(mat)
    for col in 1:ncols
        mmax = maximum(mat[:, col])
        mmin = minimum(mat[:, col])
        mrange = mmax - mmin
        for row in 1:nrows
            if fs[col] == maximum
                normalizedMat[row, col] = (mat[row, col] - mmin) / mrange
            elseif fs[col] == minimum
                normalizedMat[row, col] = (mmax - mat[row, col]) / mrange
            else
                @error fs[col]
                error("Function not defined")
            end
        end
        if fs[col] == maximum
            normalizedReferenceRow[col] = (referenceRow[col] - mmin) / mrange
        elseif fs[col] == minimum
            normalizedReferenceRow[col] = (mmax - referenceRow[col]) / mrange
        end
    end

    absoluteValueMat = similar(mat)
    for row in 1:nrows
        absoluteValueMat[row,:] = normalizedReferenceRow .- normalizedMat[row, :]  .|> abs
    end

    deltamin = absoluteValueMat |> minimum
    deltamax = absoluteValueMat |> maximum
    zetadeltamax = zeta * deltamax

    greytable = similar(absoluteValueMat)
    for row in 1:nrows
        for col in 1:ncols
            greytable[row, col] = (deltamin + zetadeltamax) / (absoluteValueMat[row, col] + zetadeltamax)
        end
    end

    weightedsums = zeros(Float64, nrows)
    for i in 1:nrows
        weightedsums[i] = w .* greytable[i, :] |> sum
    end

    orderings = sortperm(weightedsums)
    bestIndex = orderings |> last

    result = GreyResult(
        referenceRow,
        normalizedMat,
        absoluteValueMat,
        greytable,
        weightedsums,
        orderings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


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.

# =================================================================================================
"""
    mutable struct Pix
        ptr::Ptr{Cvoid}
    end

A wrapper for the PIX structure in the leptonica library.

__Values:__

| Name | Description
| :--- | :----------
| ptr  | The pointer to the Pix object allocated by the C library.

__Details:__

When the garbage collector collects this object the associated PIX object will be freed in the
library.
"""
mutable struct Pix
    ptr::Ptr{Cvoid}
    function Pix(ptr::Ptr{Cvoid})
        local retval = new(ptr)

        finalizer(retval) do obj
            pix_delete!(obj)
        end

        retval
    end
end

# =================================================================================================
"""
    show(
        io::IO,
        pix::Pix
    )::Nothing

Display summary information about the pix image.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | io   |         | The stream to write the information to.
| R | pix  |         | The image to display information about.
"""
function Base.show(
            io::IO,
            pix::Pix
        )::Nothing
    if pix.ptr != C_NULL
        local dims = pix_get_dimensions(pix)
        print(io, "Image ($(dims.w), $(dims.h)) at $(dims.d)ppi")
    else
        print(io, "Freed image.")
    end
    nothing
end

# =================================================================================================
"""
    is_valid(
        pix::Pix
    )::Bool

Check if the image has been freed or if it's still valid.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | pix  |         | The image to check.
"""
function is_valid(
            pix::Pix
        )::Bool
    return pix.ptr != C_NULL
end

# =================================================================================================
"""
    unsafe_convert(
        ::Type{Ptr{Cvoid}},
        pix::Pix
    )::Ptr{Cvoid}

"Convert" the image into a the pointer used by the Leptonica library.

__Arguments:__

| T | Name               | Default | Description
|---| :----------------- | :------ | :----------
| R | ::Type{Ptr{Cvoid}} |         | The type to convert into.
| R | pix                |         | The image to return the Leptonica pointer for.
"""
Base.unsafe_convert(::Type{Ptr{Cvoid}}, pix::Pix) = pix.ptr

# =================================================================================================
"""
    pix_delete!(
        pix::Pix
    )::Nothing

Release the Pix object so it can be freed.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | pix  |         | The image to release.

__Details:__

This method is called automatically by the garbage collector but can be called manually to release
the object early.  This method can be called multiple times without any negative effects.

Calling this method will free the object unless a reference is held by an external library.  Once
that library releases it's reference the Pix object should be fully freed.  However once
`pix_delete!()` is called on an object passing the object to any other library call will result in
an error.

Note: This method is not thread safe.
"""
function pix_delete!(pix::Pix)::Nothing
    if pix.ptr != C_NULL
        ccall((:pixDestroy, LEPTONICA), Cvoid, (Ptr{Ptr{Cvoid}},), Ref(pix.ptr))
        pix.ptr = C_NULL
    end
    nothing
end


text:
using Test
using TypedTables
using SplitApplyCombine
using Tables

include("properties.jl")
include("Table.jl")
include("FlexTable.jl")


text:
const SAC_WORD_SIZE = 4

"""
    decodesacbytes(T::Type, bs::Vector{UInt8}, needswap=false)

Decode the bytes `bs` into the type `T` swapping the byte order if
`needswap=true`. Returns an array of decoded values in the order they were
decoded.

    decodesacbytes(T::Type{String}, bs::Vector{UInt8}, stringsize=8)

Decode the bytes `bs` into `Strings` with each string being `stringsize`
bytes long. Returns an array of strings in the order they were decoded.
"""
decodesacbytes(T::Type{HeaderEnum}, bs::Vector{UInt8}, needswap=false) =
    reinterpret(HeaderEnum, decodesacbytes(Int32, bs, needswap))
decodesacbytes(T::Type{Bool}, bs::Vector{UInt8}, needswap=false) =
    decodesacbytes(Int32, bs, needswap) .!= 0 # SAC treats anything that isn't 0 as true
function decodesacbytes(T::Type, bs::Vector{UInt8}, needswap=false)
    @assert(length(bs) % 4 == 0, "Length of byte array is not a multiple of 4")
    decdata = reinterpret(T, bs)
    return needswap ? map(bswap, decdata) : decdata
end
function decodesacbytes(T::Type{String}, bs::Vector{UInt8}, stringsize::Integer=8)
    @assert(length(bs) % stringsize == 0,
            "Length of byte array is not a multiple of $stringsize")
    nstr = div(length(bs), stringsize)
    return [ascii(String(bs[(stringsize * n) + 1 : (stringsize * n) + stringsize]))
            for n in 0:1:nstr-1]
end

"Determine if a SAC file is in non-native endianness by checking the header
version number is between 1 and `SAC_HDR_VERSION`. Returns `true` if the file is
alien."
function isalienend(f::IOStream)
    p = position(f)
    seek(f, 76 * SAC_WORD_SIZE)
    nvhdr = reinterpret(Int32, read(f, 4))[1]
    seek(f,p)

    return !(1 <= nvhdr <= SAC_HDR_VERSION)
end

"""
    parsetext(T::Type, text::String, colwidth::Integer)

Parse values of type `T` from `text` where each value is `colwidth` columns wide
in the text. Returns an array of `T`s. Newline characters are stripped from
`text` before parsing.

Raises an error if the length of `text` is not a multiple of `colwidth` after
stripping newlines.
"""
parsetext(T::Type{Int32}, text::String) = parsetext(T, text, 10)
parsetext(T::Type{Float32}, text::String) = parsetext(T, text, 15)
parsetext(T::Type{HeaderEnum}, text::String) = reinterpret(HeaderEnum, parsetext(Int32, text))
parsetext(T::Type{Bool}, text::String) = parsetext(Int32, text) .!= 0
parsetext(T::Type{String}, text::String) = parsetext(T, text, 8)
function parsetext(T::Type, text::String, colwidth::Integer)
    cleantext = replace(text, '\n', "")
    @assert(length(cleantext) % colwidth == 0,
            "Text is not a multiple of $colwidth characters wide after stripping newline characters.")
    readvals = Array(T, div(length(cleantext), colwidth))

    for i in eachindex(readvals)
        field = cleantext[(i-1) * colwidth + 1 : (i-1) * colwidth + colwidth]
        readvals[i] = T <: Number ? parse(T, field) : strip(field)
    end

    return readvals
end

################################################################################
#                                  HEADER IO
################################################################################

"""
    readsachdr(fname::AbstractString; ascii=false)
    readsachdr(f::IOStream; ascii=false)

Read the header data (first 158 words/632 bytes) from the file named `fname` or
the stream `f`. Returns a `Header` instance constructed from the data.

Pass `ascii=true` to read the header of a SAC alpha file.
"""
readsachdr(fname::AbstractString; kwargs...) = open((f) -> readsachdr(f; kwargs...), fname)
function readsachdr(f::IOStream; ascii=false)
    ascii ? _readsachdr_ascii(f) : _readsachdr_binary(f)
end

function _readsachdr_ascii(f)
    seekstart(f)
    fhdrcontents = ascii(String(read(f, SAC_HDR_ASCII_END)))
    # For the alphanumeric section of the header, we have to offset from the
    # start of this section to account for the different length of the kevnm
    # field.
    hdrarr = [parsetext(Float32, fhdrcontents[SAC_HDR_ASCII_FLOAT_START:SAC_HDR_ASCII_INT_START-1]);
              parsetext(Int32, fhdrcontents[SAC_HDR_ASCII_INT_START:SAC_HDR_ASCII_ENUM_START-1]);
              parsetext(HeaderEnum, fhdrcontents[SAC_HDR_ASCII_ENUM_START:SAC_HDR_ASCII_BOOL_START-1]);
              parsetext(Bool, fhdrcontents[SAC_HDR_ASCII_BOOL_START:SAC_HDR_ASCII_ALPHA_START-1]);
              parsetext(String, fhdrcontents[SAC_HDR_ASCII_ALPHA_START:SAC_HDR_ASCII_ALPHA_START+7]);
              parsetext(String, fhdrcontents[SAC_HDR_ASCII_ALPHA_START+8:SAC_HDR_ASCII_ALPHA_START+23], 16);
              parsetext(String, fhdrcontents[SAC_HDR_ASCII_ALPHA_START+25:SAC_HDR_ASCII_END])]

    hdr = Header()
    for (val, field) in zip(hdrarr, fieldnames(Header))
        setfield!(hdr, field, val)
    end

    return hdr
end

function _readsachdr_binary(f::IOStream)
    seekstart(f)
    needswap = isalienend(f)
    bs = read(f, SAC_WORD_SIZE * SAC_HDR_NWORDS)
    hdr = Header()

    hdrvals = vcat(decodesacbytes(Float32, bs[1:SAC_WORD_SIZE * 70], needswap),
                   decodesacbytes(Int32, bs[(SAC_WORD_SIZE * 70) + 1 : SAC_WORD_SIZE * 85], needswap),
                   decodesacbytes(HeaderEnum, bs[(SAC_WORD_SIZE * 85) + 1 : SAC_WORD_SIZE * 105], needswap),
                   decodesacbytes(Bool, bs[(SAC_WORD_SIZE * 105) + 1 : SAC_WORD_SIZE * 110], needswap),
                   decodesacbytes(String, bs[(SAC_WORD_SIZE * 110) + 1 : SAC_WORD_SIZE * 112]),
                   decodesacbytes(String, bs[(SAC_WORD_SIZE * 112) + 1 : SAC_WORD_SIZE * 116], 16),
                   decodesacbytes(String, bs[(SAC_WORD_SIZE * 116) + 1 : SAC_WORD_SIZE * 158]))

    for (field, val) in zip(fieldnames(hdr), hdrvals)
        setfield!(hdr, field, val)
    end

    cleanhdr!(hdr)
    return hdr
end

################################################################################
#                                   DATA IO
################################################################################

"""
    readsac(fname::AbstractString; ascii=false)
    readsac(f::IOStream; ascii=false)

Read the SAC file at path `fname` or from the stream `f`. The return type is
determined using the file type declared in the file's header. *Note* this isn't
type stable.

    readsac{S<:AbstractSACData}(T::Type{S}, f::IOStream; ascii=false)

Read the SAC file of type `T` from the stream `f`. Raises an error if the type
`T` does not match the file type declared in the file's header. The returned
type is always of type `T`.

Pass `ascii=true` to read a SAC alpha file. *Note* this uses *significantly*
more memory than reading a binary SAC file.
"""
readsac(fname::AbstractString; kwargs...) = open((f) -> readsac(f; kwargs...), fname)
function readsac(f::IOStream; kwargs...)
    hdr = readsachdr(f; kwargs...)

    if hdr.iftype == itime && hdr.leven
        readsac(EvenTimeSeries, f, hdr; kwargs...)
    elseif hdr.iftype == itime && !hdr.leven
        readsac(UnevenTimeSeries, f, hdr; kwargs...)
    elseif hdr.iftype == irlim
        readsac(ComplexSpectrum, f, hdr; kwargs...)
    elseif hdr.iftype == iamph
        readsac(AmplitudeSpectrum, f, hdr; kwargs...)
    elseif hdr.iftype == ixy
        readsac(GeneralXY, f, hdr; kwargs...)
    end
end

readsac{S<:AbstractSACData}(T::Type{S}, f::IOStream; kwargs...) = readsac(T, f, readsachdr(f; kwargs...); kwargs...)
function readsac{S<:AbstractSACData}(T::Type{S}, f::IOStream, hdr::Header; kwargs...)
    if hdr.iftype != FILE_TYPE_ENUMS[T]
        error("File's header indicates it is not a $(T)")
    end
    if T == EvenTimeSeries
        T(hdr, readsac_data(f, hdr.npts; kwargs...)[1])
    elseif T == ComplexSpectrum
        T(hdr, complex(readsac_data(f, hdr.npts; kwargs...)...))
    else
        T(hdr, readsac_data(f, hdr.npts; kwargs...)...)
    end
end

"""
    readsac_data(f::IOStream, npts::Int32; ascii=false)

Read `npts` of the data section from the stream `f`. Will automatically swap the
byte order of the data if it isn't the host's byte order. Returns a 2-tuple
containing the first data section and, if present, the second data section.
"""
function readsac_data(f::IOStream, npts::Int32; ascii=false)
    ascii ? _readsac_data_ascii(f, npts) : _readsac_data_binary(f, npts)
end

function _readsac_data_binary(f::IOStream, npts::Int32)
    needswap = isalienend(f)
    seek(f, DATA_START)

    data1 = decodesacbytes(Float32, read(f, SAC_WORD_SIZE * npts), needswap)
    if eof(f)
        return (data1, Float32[])
    end

    data2 = decodesacbytes(Float32, read(f, SAC_WORD_SIZE * npts), needswap)
    return (data1, data2)
end

function _readsac_data_ascii(f::IOStream, npts::Int32)
    seek(f, ASCII_DATA_START)

    data = readstring(f)
    data = parsetext(Float32, data)

    if length(data) == 2 * npts
        return (data[1:npts], data[npts+1:end])
    elseif length(data) == npts
        return (data, Float32[])
    else
        error("Number of data points in file ($(length(data)) does not match
       the"*"number of points declared in the header ($(hdr.npts)).")
    end
end


text:
function compute_grad_ratio(grads)
    nsamps = length(grads)
    return grads[2:nsamps-1,:]./grads[1:nsamps-2,:]
end


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
#using PyCall
using Logging
using Printf
using Suppressor

#include("../../phantoms/simple_phantoms.jl")
#include("../tv_primaldual_flow.jl")
#include("../optical_flow.jl")
include("../../phantoms/sperm_phantom.jl")

H,W = 609, 609
detmin, detmax = -38.0, 38.0
function radon_operator(height, width, detcount, angles)
    proj_geom = ProjGeom(0.5, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, height, width, -38.0,38.0, -38.0,38.0)
    return A
end

#Create phantom at smaller gridsize than reconstruction
grid = collect(detmin:0.1:detmax)
r(s) = 1.0
frames, tracks = get_sperm_phantom(2,r,grid)
detcount = Int(floor(H*1.4))
#match size of input image (generating data)
nangles = 4
a = range(0.0,2*pi+pi/4, length=nangles*size(frames)[3])#rand(0.0:0.001:π, nangles*size(frames)[3])
angles = map(t -> a[(t-1)*nangles+1:(t-1)*nangles+nangles], 1:size(frames)[3])
As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles[t]),1:size(frames)[3])

bsex = zeros(size(As[1])[1],size(frames)[3])

map(t -> bsex[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
@info "adding gaussion noise at level 0.01"
rho = 0.01
e = randn(size(bsex));
e = rho*norm(bsex)*e/norm(e);
bs = bsex + e;

u0s = zeros(H,W,size(frames)[3])

As  = map(t -> radon_operator(H,W, detcount, angles[t]),1:size(u0s)[3])
cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

u0s = zeros(H,W,size(frames)[3])
mask = ones(size(u0s)...)

c=1.1
w_tvs = [0.0001,0.001,0.01,0.05,0.1,0.2,0.3,0.5]
w_flows = [0.0001,0.001,0.01,0.05,0.1,0.2,0.3,0.5]
prog = Progress(size(w_tvs,1)*size(w_flows,1),60, "Testing parameters")

for w_tv = w_tvs
    niter=450
    us_tv = zeros(H,W,size(frames)[3])
    for t = 1:size(frames)[3]
        A = As[t]
        p = bs[:,t]
        u0 = u0s[:,:,t]
        @suppress begin
            us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
        end
    end
    for w_flow = w_flows
        niter=50
        @suppress begin
            us_flow = recon2d_tv_primaldual_flow(As, bs, deepcopy(us_tv), niter, 200, w_tv, w_flow, c, mask, 0.012)
            filenamepng = @sprintf "flow_2ms_%f_%f__3.png" w_tv w_flow;
            #PLOT
            plot(Gray.(us_flow[:,:,1]), aspect_ratio=:equal, framestyle=:none, legend=false, yflip=false)
            plot!(tracks[1][:,1], tracks[1][:,2], aspect_ratio=:equal, linewidth=5)
            savefig(filenamepng)
            next!(prog)
        end
    end
end

# angles = rand(1:0.001:π,3)
# #Create phantom at smaller gridsize than reconstruction
# frames, tracks = get_sperm_phantom(2,r, grid)
# detcount = Int(floor(H*1.4))
# #match size of input image (generating data)
# As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles[t]),1:size(frames)[3])
# bsex = zeros(size(As[1])[1],size(frames)[3])
#
# map(t -> bsex[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
# @info "adding gaussion noise at level 0.01"
# rho = 0.01
# e = randn(size(bsex));
# e = rho*norm(bsex)*e/norm(e);
# bs = bsex + e;
#
# As  = map(t -> radon_operator(H,W, detcount, angles[t]),1:size(u0s)[3])
#
# for w_tv = w_tvs
#     niter=450
#     us_tv = zeros(H,W,size(frames)[3])
#     for t = 1:size(frames)[3]
#         A = As[t]
#         p = bs[:,t]
#         u0 = u0s[:,:,t]
#         @suppress begin
#             us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
#         end
#     end
#     for w_flow = w_flows
#         niter=50
#         @suppress begin
#             us_flow = recon2d_tv_primaldual_flow(As, bs, deepcopy(us_tv), niter, 20, w_tv, w_flow, c)
#             filenamepng = @sprintf "flow_radn_ang_2ms_%f_%f_.png" w_tv w_flow;
#             #PLOT
#             plot(Gray.(us_flow[:,:,1]), aspect_ratio=:equal, framestyle=:none, legend=false, yflip=false)
#             savefig(filenamepng)
#             next!(p)
#         end
#     end
# end
# cd(cwd)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
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

const NCURSES_ATTR_SHIFT =8
function NCURSES_BITS( m, shf)
    m << (shf + NCURSES_ATTR_SHIFT)
end

COLOR_PAIR( n ) = NCURSES_BITS(n, 0)

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

# not really used. See getmouse() hack in ccall.jl
NCURSES_MOUSE_MASK(b,m) = m<<((b-1)*6) # NCURSES_MOUSE_VERSION=1
const NCURSES_BUTTON_PRESSED = @compat 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 ) # NCURSES_MOUSE_VERSION=1


text:
"""visualization functions"""

using Plots

COLORS = ["royalblue", "orange", "lime", "gray", "coral", "brown4", "aquamarine4"]

"""get color given agent's index"""
function get_color(i::Int64)::String
    return COLORS[mod1(i, length(COLORS))]
end

"""plot circle in 2D space"""
function plot_circle!(
    x::Float64,
    y::Float64,
    r::Float64,
    color::String = "black",
    lw::Float64 = 1.0;
    fillalpha::Float64 = 0.0,
)
    Xs = cos.(collect(0:0.1:pi+0.1)) * r
    Ys = sin.(collect(0:0.1:pi+0.1)) * r
    offset_x = fill(x, length(Xs))
    offset_y = fill(y, length(Ys))
    plot!(
        Xs + offset_x,
        Ys + offset_y,
        fillrange = -Ys + offset_y,
        color = color,
        label = nothing,
        lw = lw,
        fillalpha = fillalpha,
    )
    plot!(Xs + offset_x, -Ys + offset_y, color = color, lw = lw, label = nothing)
end

"""plot sphere in 3D space"""
function plot_sphere!(
    x::Float64,
    y::Float64,
    z::Float64,
    r::Float64,
    color = "black",
    step::Int64 = 20,
)
    for a = -r:r/step:r
        b = sqrt(r^2 - a^2)
        Xs = cos.(collect(0:0.1:2π+0.1)) * b
        Ys = sin.(collect(0:0.1:2π+0.1)) * b
        plot!(
            Xs + fill(x, size(Xs)),
            Ys + fill(y, size(Ys)),
            fill(z + a, size(Xs)),
            legend = false,
            color = color,
        )
    end
end

"""plot obstacles"""
function plot_obs!(obstacles::Vector{Obs}) where {Obs<:Obstacle}
    for o in obstacles
        if Obs == CircleObstacle2D
            plot_circle!(o.x, o.y, o.r; fillalpha = 1.0)
        elseif Obs == CircleObstacle3D
            plot_sphere!(o.x, o.y, o.z, o.r, "black")
        end
    end
    return plot!()
end

"""plot roadmaps"""
function plot_roadmap!(
    V::Union{Nothing,Vector{Vector{Node{State}}}},
    ins_params...;
    lw::Float64 = 0.2,
) where {State<:AbsState}
    isnothing(V) && return nothing

    for i = 1:length(V)
        params = Dict(
            :lw => lw,
            :markersize => 2,
            :markershape => :circle,
            :legend => nothing,
            :color => get_color(i),
            :markerstrokecolor => get_color(i),
        )
        for v in V[i]
            for u_id in v.neighbors
                u = V[i][u_id]
                plot_motion!(v.q, u.q, map(arr -> arr[i], ins_params)..., params)
            end
            if isempty(v.neighbors)
                plot_motion!(v.q, v.q, map(arr -> arr[i], ins_params)..., params)
            end
        end
    end

    return plot!()
end

"""plot start and goal states for all agents"""
function plot_start_goal!(
    config_init::Vector{State},
    config_goal::Vector{State},
    ins_params...;
) where {State<:AbsState}

    for (i, (q_init, q_goal)) in enumerate(zip(config_init, config_goal))
        params = Dict(
            :seriestype => :scatter,
            :markersize => 5,
            :label => nothing,
            :color => get_color(i),
        )
        plot_start_goal!(q_init, q_goal, map(arr -> arr[i], ins_params)..., params)
    end
    return plot!()
end

"""plot trajectory of all agents"""
function plot_traj!(
    solution::Union{Nothing,Vector{Vector{Node{State}}}},
    ins_params...;
    lw::Float64 = 3.0,
    VERBOSE::Int64 = 0,
) where {State<:AbsState}

    isnothing(solution) && return

    N = length(solution[1])
    T = length(solution)
    for (t, Q_to) in enumerate(solution[2:end])
        VERBOSE > 0 && @printf("\rplotting t=%d/%d", t, T)
        Q_from = solution[t]
        for i = 1:N
            params = Dict(:color => get_color(i), :lw => lw, :label => nothing)
            plot_motion!(Q_from[i].q, Q_to[i].q, map(arr -> arr[i], ins_params)..., params)
        end
    end
    VERBOSE > 0 && @printf("\n")
    return plot!()
end

"""setup artboard"""
function plot_init!(State::DataType)
    if State in [StatePoint3D, StateArm33, StateCapsule3D]
        plot3d(
            size = (400, 400),
            xlim = (0, 1),
            ylim = (0, 1),
            zlim = (0, 1),
            xticks = nothing,
            yticks = nothing,
            zticks = nothing,
        )
    else
        plot(
            size = (400, 400),
            xlim = (0, 1),
            ylim = (0, 1),
            framestyle = :box,
            yflip = true,
            xmirror = true,
            xticks = nothing,
            yticks = nothing,
        )
    end
end

"""save figure"""
function safe_savefig!(filename::Union{Nothing,String} = nothing)
    isnothing(filename) && return
    dirname = join(split(filename, "/")[1:end-1], "/")
    !isdir(dirname) && mkpath(dirname)
    savefig(filename)
end

"""
    plot_instance!(
        config_init::Vector{State},
        config_goal::Vector{State},
        obstacles::Vector{Obs} where {Obs<:Obstacle},
        ins_params...;                                   # positions, radius
        filename::Union{Nothing,String} = nothing,       # save plot if given
    ) where {State<:AbsState}

visualize instance
"""
function plot_instance!(
    config_init::Vector{State},
    config_goal::Vector{State},
    obstacles::Vector{Obs} where {Obs<:Obstacle},
    ins_params...;
    filename::Union{Nothing,String} = nothing,
) where {State<:AbsState}

    N = length(config_init)
    plot_init!(State)
    plot_obs!(obstacles)
    foreach(
        i -> plot_agent!(config_init[i], map(arr -> arr[i], ins_params)..., get_color(i)),
        1:N,
    )
    plot_start_goal!(config_init, config_goal, ins_params...)
    safe_savefig!(filename)
    return plot!()
end

"""
    plot_res!(
        config_init::Vector{State},
        config_goal::Vector{State},
        obstacles::Vector{Obs} where {Obs<:Obstacle},
        ins_params...;                                                    # positions, radius
        roadmaps::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
        solution::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
        filename::Union{Nothing,String} = nothing,                        # save plot if given
    ) where {State<:AbsState}

visualize planning result
"""
function plot_res!(
    config_init::Vector{State},
    config_goal::Vector{State},
    obstacles::Vector{Obs} where {Obs<:Obstacle},
    ins_params...;
    roadmaps::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
    solution::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
    filename::Union{Nothing,String} = nothing,
    VERBOSE::Int64 = 0,
) where {State<:AbsState}

    plot_init!(State)
    plot_obs!(obstacles)
    plot_roadmap!(roadmaps, ins_params...)
    plot_traj!(
        solution,
        ins_params...;
        lw = isnothing(roadmaps) ? 1.0 : 3.0,
        VERBOSE = VERBOSE,
    )
    plot_start_goal!(config_init, config_goal, ins_params...)
    foreach(
        i -> plot_agent!(config_init[i], map(arr -> arr[i], ins_params)..., get_color(i)),
        1:length(config_init),
    )
    safe_savefig!(filename)
    return plot!()
end

"""
    plot_anim!(
        config_init::Vector{State},
        config_goal::Vector{State},
        obstacles::Vector{Obs} where {Obs<:Obstacle},
        ins_params...;
        solution::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
        filename::String = "tmp.gif",
        fps::Int64 = 10,
        interpolate_depth::Union{Nothing,Int64} = nothing,  # plots intermediate states
        VERBOSE::Int64 = 0,
    ) where {State<:AbsState}

visualize animation
"""
function plot_anim!(
    config_init::Vector{State},
    config_goal::Vector{State},
    obstacles::Vector{Obs} where {Obs<:Obstacle},
    ins_params...;
    solution::Union{Nothing,Vector{Vector{Node{State}}}} = nothing,
    filename::String = "tmp.gif",
    fps::Int64 = 10,
    interpolate_depth::Union{Nothing,Int64} = nothing,
    VERBOSE::Int64 = 0,
) where {State<:AbsState}

    if isnothing(solution)
        @warn "solution has not computed yet!"
        return
    end

    T = length(solution)
    anim = @animate for (t, Q) in enumerate(vcat(solution, [solution[end]]))
        VERBOSE > 0 && @printf("\rplotting t = %d / %d", t, T)
        plot_init!(State)
        plot_obs!(obstacles)
        plot_traj!(solution, ins_params...; lw = 1.0)
        plot_start_goal!(config_init, config_goal, ins_params...)

        if !isnothing(interpolate_depth) && interpolate_depth > 0 && 1 < t <= T
            # compute intermediate states
            C_arr = Array{Any}(undef, 2 + sum(map(k -> 2^k, 0:interpolate_depth-1)))
            add_mid_config(d, ind1, ind2) = begin
                if d > 0
                    C1 = C_arr[ind1]
                    C2 = C_arr[ind2]
                    C = map(e -> get_mid_status(e...), zip(C1, C2))
                    ind = Int64((ind1 + ind2) / 2)
                    C_arr[ind] = C
                    add_mid_config(d - 1, ind1, ind)
                    add_mid_config(d - 1, ind, ind2)
                end
            end

            C_arr[1] = map(v -> v.q, solution[t-1])
            C_arr[end] = map(v -> v.q, solution[t])
            add_mid_config(interpolate_depth, 1, length(C_arr))
            for Q_tmp in C_arr
                for (i, q) in enumerate(Q_tmp)
                    plot_agent!(q, map(arr -> arr[i], ins_params)..., get_color(i))
                end
            end
        else
            for (i, v) in enumerate(Q)
                plot_agent!(v.q, map(arr -> arr[i], ins_params)..., get_color(i))
            end
        end
    end
    VERBOSE > 0 && println()

    dirname = join(split(filename, "/")[1:end-1], "/")
    !isdir(dirname) && mkpath(dirname)

    return gif(anim, filename, fps = fps)
end

"""
    plot_tpg!(
        TPG::Vector{Vector{Action{State}}};
        init::Bool = true,
        offset::Int64 = 0,
        filename::Union{Nothing,String} = nothing,
    )

visualize temporal plan graph
"""
function plot_tpg!(
    TPG::Vector{Vector{Action{State}}};
    init::Bool = true,
    offset::Int64 = 0,
    filename::Union{Nothing,String} = nothing,
) where {State<:AbsState}

    N = length(TPG)
    init && plot(size = (200, 400), xlim = (0.5, N + 0.5), xticks = 1:N)
    params =
        Dict(:label => nothing, :color => :black, :markersize => 5, :markershape => :circle)
    for (i, actions) in enumerate(TPG)
        for a in actions
            for (j, b_id) in a.successors
                k = findfirst(a -> a.id == b_id, TPG[j])
                if k != nothing
                    b = TPG[j][k]
                    plot!([i, j], [a.t + offset, b.t + offset]; params...)
                end
            end
            annotate!(i - 0.1, a.t + offset, text(a.id, :black, :right, 3))
        end
        length(actions) == 1 && scatter!([i], [actions[1].t + offset]; params...)
    end
    safe_savefig!(filename)
    return plot!()
end