text:
module UCTransition
#
# This module models the transition from Legacy Means Tested benefits
# to Universal Credit.
#
# Currently it's based crudely on Scotland-wide transition figures from HoC Research
# https://commonslibrary.parliament.uk/constituency-data-universal-credit-roll-out/#caseload
#

using ScottishTaxBenefitModel
using .RunSettings  # : Settings, MT_Routing
using .Intermediate: MTIntermediate, HHIntermed
using .ModelHousehold
using .Definitions
using .Randoms: testp
using .Results: tozero!, HouseholdResult, BenefitUnitResult, LMTResults, UCResults
using .STBIncomes

export
    route_to_uc_or_legacy,
    route_to_uc_or_legacy!

@enum ClaimantType trans_all trans_housing trans_w_kids trans_incapacity trans_jobseekers

const PROPS_ON_UC = Dict(
    trans_all => 0.59,
    trans_housing => 0.61,
    trans_w_kids => 0.55,
    trans_incapacity => 0.32,
    trans_jobseekers => 0.95 )

"""
Allocate a benefit unit to ether UC or Legacy Benefits. Very crudely.
"""
function route_to_uc_or_legacy(
    settings :: Settings,
    bu       :: BenefitUnit,
    intermed :: MTIntermediate ) :: LegacyOrUC
    if settings.means_tested_routing != modelled_phase_in
        return settings.means_tested_routing == uc_full ? uc_bens : legacy_bens
    end
    prob = 0.0
    if intermed.num_job_seekers > 0
        prob = PROPS_ON_UC[trans_jobseekers]
    elseif intermed.limited_capacity_for_work
        prob = PROPS_ON_UC[trans_incapacity]
    elseif intermed.num_children > 0
        prob = PROPS_ON_UC[trans_w_kids]
    elseif intermed.benefit_unit_number == 1
        prob = PROPS_ON_UC[trans_housing]
    else
        prob = PROPS_ON_UC[trans_all]
    end
    head = get_head( bu )
    switch = testp( head.onerand, prob, Randoms.UC_TRANSITION )
    return switch ? uc_bens : legacy_bens
end

function route_to_uc_or_legacy!(
    results  :: HouseholdResult,
    settings :: Settings,
    hh       :: Household,
    intermed :: HHIntermed )
    bus = get_benefit_units( hh )
    RT = eltype( results.income ) # FIXME whole where RT thing

    for bno in eachindex( bus )
        im = intermed.buint[bno]
        bres = results.bus[bno]
        if bres.uc.basic_conditions_satisfied # FIXME This condition needs some thought.
            bres.route = route_to_uc_or_legacy( settings, bus[bno], im )
            if bres.route == legacy_bens
                tozero!( bres, UNIVERSAL_CREDIT )
                bres.uc = UCResults{RT}()
                if( bno == 1 ) && ( bres.legacy_mtbens.ctr_recipient > 0 )# ctr assignment hack
                    bres.pers[bres.legacy_mtbens.ctr_recipient].income[COUNCIL_TAX_BENEFIT] = bres.legacy_mtbens.ctr
                else
                    tozero!( bres, COUNCIL_TAX_BENEFIT )
                end
            elseif bres.route == uc_bens
                # nuke every old thing for everyone who's in scope for UC
                tozero!( bres, LEGACY_MTBS...)
                bres.legacy_mtbens = LMTResults{RT}()
                if( bno == 1 ) && ( bres.uc.recipient > 0 )# ctr assignment hack
                    bres.pers[bres.uc.recipient].income[COUNCIL_TAX_BENEFIT] = bres.uc.ctr
                else
                    tozero!( bres, COUNCIL_TAX_BENEFIT )
                end
                ## FIXME TRANSITIONAL PAYMENTS
            end
        end
        # .. so some futher nuking ..
        if settings.means_tested_routing == uc_full
            # these cease to exist, even for pensioners and students
            tozero!( bres,
                NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE,
                NON_CONTRIB_JOBSEEKERS_ALLOWANCE,
                INCOME_SUPPORT,
                WORKING_TAX_CREDIT,
                CHILD_TAX_CREDIT )
            bres.legacy_mtbens = LMTResults{RT}()
            # this last is entirely arbitrary, but
            # takes out essentially all ft students, etc.
            if im.age_oldest_adult < 50
                tozero!( bres, LEGACY_MTBS...)
            end
        end
    end
end

end # Module UCTransition

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset begin
    ga = rast(rand(4, 4, 10))

    @test size(ga[:,:, 1]) == (4, 4, 1)
    @test size(ga[2:4, 1:4, :]) == (3, 4, 10)
end

ext: jl


text:
#
# PauliString struct
#

struct PauliString
	N::Int  # length
	sites::Vector{Int} # To do: change this to bit number for "used sites" 000101 = operator at site 4 and 6
	# vector which contains site indicator where the string has non-trivial operators
	baseIdx::Vector{Int} # indicates the operator (σx, σy, σz)
	coef::Base.RefValue{ComplexF64}

	function PauliString(N::Int, site::Vector{Int}, baseIdx::Vector{Int}, coef::Number)
		length(baseIdx[baseIdx.>3]) == 0 || throw(ArgumentError("baseIdx $(baseIdx) containes unvalid arguments >3.
									Note the convention is: (1,2,3) = (σx, σy, σz)"))
		return new(N, site, baseIdx, Ref(convert(ComplexF64, coef)))
	end
end

coef(p::PauliString) = p.coef[]
set_coef!(p::PauliString, a::Number) = p.coef[] = a

import Base.==
import Base.+
import Base.*
import Base.conj
(==)(a::PauliString, b::PauliString) = a.sites == b.sites && a.baseIdx == b.baseIdx && coef(a) == coef(b) ? true : false
(∥)(a::PauliString, b::PauliString) = a.sites == b.sites && a.baseIdx == b.baseIdx ? true : false
(+)(a::PauliString, b::PauliString) = a∥b ? PauliString(a.N, a.sites, a.baseIdx, coef(a)+coef(b)) : [a,b]
(*)(a::PauliString, c::Number) = PauliString(a.N, a.sites, a.baseIdx, coef(a)*c)
(*)(c::Number, a::PauliString) = PauliString(a.N, a.sites, a.baseIdx, coef(a)*c)

function (*)(a::PauliString, b::PauliString)

	p1_sites = a.sites
	p2_sites = b.sites
	p1_idx = a.baseIdx.+1 #  in order that 0 (id) are not get cut out later on
	p2_idx = b.baseIdx.+1

	pArray = sparse([[1 for i in p1_sites]...,
			[2 for i in p2_sites]...],
			[p1_sites...,p2_sites...],[p1_idx..., p2_idx...])
	new_σ = []
	sites = [Set([x[2] for x in  findall(x->x!=0,pArray)])...][end:-1:1]
	#@show sites
	for col in sites
		#@show pArray[1,col], pArray[2,col]
		i = pArray[1,col] == 0 ? 1 : pArray[1,col]
		j = pArray[2,col] == 0 ? 1 : pArray[2,col]
		push!(new_σ, σσ(i, j))  # +1 to match the convention
     					   # 0 = σ⁰, 1 = σˣ , ...
	end

	new_idx = [Int(new_σ[i][1])-1 for i in 1:size(sites)[1]]
	c = *([new_σ[i][2] for i in 1:size(sites)[1]]...)

	idx = sortperm(sites)
	sites  = sites[idx]
	new_idx = new_idx[idx]

	return PauliString(a.N, sites, new_idx, c*coef(a)*coef(b))
end

function Base.conj!(a::PauliString)
 	set_coef!(a, (coef(a))')
	return a
end

function Base.conj(a::PauliString)
	b = deepcopy(a)
	set_coef!(b, (coef(b))')
	return b
end

function Base.show(io::IO, p::PauliString)
	print_string = fill("𝟙-",p.N)
	op_symbols = ["σˣ-","σʸ-","σᶻ-"]
	print_string[p.sites] .= op_symbols[p.baseIdx]

	output_string = chop(string("($(coef(p)))⋅",print_string...))
	print(io, output_string)
end

function rev_sites(a::PauliString)
	return PauliString(a.N, [a.N.-(a.sites.-1)...], a.baseIdx, coef(a))
end

# get cycles of a permutation
function getcycles(p::Vector{T}) where {T<:Any}
	isperm(p) || throw(DomainError("p is not a permutation"))

	n = length(p)
	used = falses(n)
	cycles = []

	for k = 1:n
		if used[k]; continue; end
		used[k] = true
		j = p[k]

		cycle = [j]
		# go through one cycle
		while !used[j]
			used[j] = true
			j = p[j]
			push!(cycle, j)
		end

		if length(cycle)>1
			cycle[1:end] = [cycle[end], cycle[1:end-1]...]
		end
		push!(cycles, cycle)

	end
	return cycles
end

# see
# https://github.com/toivoh/ConicHulls.jl/blob/888109ab6fe7d334ffa48e649947f0c293a967ae/src/Common.jl#L37
function isevenperm(p::Vector{T}) where {T<:Any}
	isperm(p) || throw(DomainError("p is not a permutation"))

	n = length(p)
	used = falses(n)
	even = true

	for k = 1:n
		if used[k]; continue; end

		# Each even cycle flips even (an odd number of times)
		used[k] = true
		j = p[k]

		# go through one cycle
		while !used[j]
			used[j] = true
			j = p[j]
			even = !even
		end
	end

	return even
end
isoddperm(p) = !isevenperm(p)

"""
    σσ(i,j)

calculates σᵢσⱼ = δᵢⱼ I   + im*ϵᵢⱼₗσₗ

# Arguments:
    - i:
    - j:

return:

"""
function σσ(i::Int, j::Int)
	coef = 1.0
	#@show i, j
	# input is σ⁰ = 1, σˣ = 2, ...
	i-=1; j-=1 # change to σ⁰ = 0, σˣ = 1, ...

	if i == 0 || j == 0 #either i or j is σ⁰
		return i+j+1, coef
	elseif i == j
		return 1, coef
	else

		# determin σ
		σ_res = filter(x->x ∉ [i,j],[1,2,3])[1]

		# determin coeff
		coef *= isoddperm([i,j,σ_res]) ? -1.0im : 1.0im

		return σ_res+1, coef
	end
end

"""
	[σᵢ , σⱼ] = 2*im*ϵᵢⱼₖ σₖ

"""
function com_σσ(i::Int, j::Int)
	i-=1; j-=1
	coef = 1
	if i == 0 || j == 0
		return 0,0
	elseif i == j
		return 0,0
	else
		# determin σ
		σ_res = filter(x->x ∉ [i,j],[1,2,3])[1]

		# determin coeff
		coef *= isoddperm([i,j,σ_res]) ? -1.0im : 1.0im

		return σ_res+1, 2*coef
	end
end

function commutePauli(p1::PauliString, p2::PauliString)
	p1_sites = p1.sites
	p2_sites = p2.sites
	p1_idx = p1.baseIdx
	p2_idx = p2.baseIdx

	pArray = sparse([[1 for i in p1_sites]...,
			[2 for i in p2_sites]...],
			[p1_sites...,p2_sites...],[p1_idx..., p2_idx...])

	# get unique sites to avoide double counting
	sites = [Set([x[2] for x in  findall(x->x!=0,pArray)])...][end:-1:1]
	new_σσ = []
	for x in sites

		# +1 for convention: label 1 is used for σ0 in σσ routine
		o1 = pArray[1,x]+1
		o2 = pArray[2,x]+1

		o1_o2 = σσ(o1, o2)
		o2_o1 = σσ(o2, o1)
		push!(new_σσ,(o1_o2, o2_o1))
	end

	coef_op = *([x[1][2] for x in new_σσ]...)-*([x[2][2] for x in new_σσ]...)

	if !iszero(coef_op)
		new_p = []
		new_site = []
		for (i,x) in enumerate(new_σσ)
			new_o = x[1][1]-1
			if !iszero(new_o)
				push!(new_p, new_o)
				push!(new_site, [sites...][i])
			end
		end
		res_p = PauliString(p1.N, Vector{Int}(new_site), Vector{Int}(new_p), coef_op*coef(p1)*coef(p2))

		return res_p

	end
end

"""
    adH(H, p, order)
calculates all elements of adH up to order n

# Arguments:

return:
!!! To do put output in a dict specifying each order.
"""
function adH(H, p, order)
	pauliList = Vector{PauliString}([deepcopy(p)])
	new_pauliList = []
	for n in 1:order
		for pauli_string in pauliList

			# get operator strings in H with support overlaping with pauli_string
			sup_string = getSupport(pauli_string)
			H_accesible = filter(x-> getSupport(x)[1]<=sup_string[2] && getSupport(x)[2]>=sup_string[1], H)
			res = filter(x->x!=nothing, [commutePauli(x,pauli_string)  for x in H_accesible])

			#res = filter(x->x!=nothing, [commutePauli(x,y)  for x in H])

			push!(new_pauliList, res...)

		end
		pauliList = Vector{PauliString}(new_pauliList)
	end

	if new_pauliList != []
		return pauliList
	else
		return nothing
	end
end

"""
    simplify_stringList(L; tol = nothing)

# Arguments:

return:
"""
function simplify_stringList(L::Vector{PauliString}; tol = nothing)

	# remove operator strings with coef = 0
	filter!(x->coef(x) != 0, L)
	# sort by number of sites wiht non-trivial operator
	sort!(L, by = x -> size(x.sites)[1])

	number_op = 1
	new_L = Vector{PauliString}([])
	while L != []

		# retrieve all operator with number_op non-trivial operator
		# and store them temporary in tmp
		idx = findlast(x->size(x.sites)[1] == number_op, L)
		tmp = idx != nothing ? [popfirst!(L) for i in 1:idx] : []

		# while tmp is not empty check for equivalent operator strings
		while tmp != []

			idx = findall(x->x∥tmp[1], tmp) # find all similar PauliStrings
			# if they sum to coef != 0  add them to the new_L
			if sum([coef.(tmp[idx])...]) != 0

				# erase zeros from string
				non_zero_idx = findall(x->x!=0,tmp[idx][1].baseIdx)
				if non_zero_idx != []
					new_base_idx = tmp[idx][1].baseIdx[non_zero_idx]
					new_sites = tmp[idx][1].sites[non_zero_idx]

					sort_sites = sortperm(new_sites)
					new_sites = new_sites[sort_sites]
					new_base_idx = new_base_idx[sort_sites]

				else

					new_base_idx = [0]
					new_sites = [1]
				end

				push!(new_L, PauliString(tmp[idx][1].N, new_sites, new_base_idx, sum([coef.(tmp[idx])...])))
			end
			deleteat!(tmp, idx) # remove visited strings from tmp

		end

		number_op += 1 # go to next number_op sector

	end

	# filter strings which coef are below tol if tol!=nothing
	if tol !=nothing new_L = filter(x->abs(coef(x))>=tol,new_L) end

	return new_L
end

getSupport(p::PauliString)  = (minimum(p.sites),maximum(p.sites))

function estimateLoc(pvec::Vector{PauliString})
	norm = 1/sum(abs.(coef.(pvec)))

	loc = norm*sum(abs.(coef.(pvec)).*[x[2]-x[1] for x in getSupport.(pvec)])

	return loc
end

max_stars_repo_name: Gawatz/PauliStrings

max_issues_repo_name: Gawatz/PauliStrings

max_forks_repo_name: Gawatz/PauliStrings


text:
module Lockyv2
using Dates

const config = Dict{Int8,Dict}(
1=> Dict{String,Any}(
"seed"=> 7077,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
2=> Dict{String,Any}(
"seed"=> 5566,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
3=> Dict{String,Any}(
"seed"=> 9106,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
4=> Dict{String,Any}(
"seed"=> 5579,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
5=> Dict{String,Any}(
"seed"=> 111,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
6=> Dict{String,Any}(
"seed"=> 9044,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
7=> Dict{String,Any}(
"seed"=> 9099,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
8=> Dict{String,Any}(
"seed"=> 9047,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
9=> Dict{String,Any}(
"seed"=> 9133,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
10=> Dict{String,Any}(
"seed"=> 9199,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
),
11=> Dict{String,Any}(
"seed"=> 9190,
"shift"=> 7,
"tlds"=> ["ru", "info", "biz", "click", "su", "work", "pl", "org", "pw", "xyz"]
)
)

function ror32(v, s)
    v &= 0xFFFFFFFF
    return ((v >> s) | (v << (32-s))) & 0xFFFFFFFF
end

function rol32(v, s)
    v &= 0xFFFFFFFF
    return ((v << s) | (v >> (32-s))) & 0xFFFFFFFF
end

function gen_one(date=missing, config_nr=1, domain_nr=missing, add_tld=false)
    domain = ""
    if date === missing
        date = Dates.Date(Dates.now())
    end
    if domain_nr === missing
        domain_nr = rand(Int16)
    end
    c = config[config_nr]
    seed_shifted = rol32(c["seed"], 17)
    dnr_shifted = rol32(domain_nr, 21)

    k = 0
    year,month,day = Dates.year(date),Dates.month(date),Dates.day(date)
    for _ in 1:7
        t_0 = ror32(0xB11924E1 * (year + k + 0x1BF5), c["shift"]) & 0xFFFFFFFF
        t_1 = ((t_0 + 0x27100001) ⊻ k) & 0xFFFFFFFF
        t_2 = (ror32(0xB11924E1 * (t_1 + c["seed"]), c["shift"])) & 0xFFFFFFFF
        t_3 = ((t_2 + 0x27100001) ⊻ t_1) & 0xFFFFFFFF
        t_4 = (ror32(0xB11924E1 * (trunc(Int64,day/2) + t_3), c["shift"])) & 0xFFFFFFFF
        t_5 = (0xD8EFFFFF - t_4 + t_3) & 0xFFFFFFFF
        t_6 = (ror32(0xB11924E1 * (month + t_5 - 0x65CAD), c["shift"])) & 0xFFFFFFFF
        t_7 = (t_5 + t_6 + 0x27100001) & 0xFFFFFFFF
        t_8 = (ror32(0xB11924E1 * (t_7 + seed_shifted + dnr_shifted), c["shift"])) & 0xFFFFFFFF
        k = ((t_8 + 0x27100001) ⊻ t_7) & 0xFFFFFFFF
        year += 1
    end

    length = (k % 11) + 7
    for i in 1:length
        k = (ror32(0xB11924E1*rol32(k, i), c["shift"]) + 0x27100001) & 0xFFFFFFFF
        domain *= 'a' + k % 25
    end

    k = ror32(k*0xB11924E1, c["shift"])
    tlds = c["tlds"]
    tld_i = ((k + 0x27100001) & 0xFFFFFFFF)%10 + 1

    if add_tld
        domain *= '.' + tlds[tld_i]
    end

    return domain
end

function gen_many(cnt=100)
    res = []
    date = Dates.Date(Dates.now())
    each_domain_items = trunc(Int64, cnt/11)
    remain_items = cnt%11
    for i in 1:11
        for j in 1:each_domain_items
            res = [res..., gen_one(date,i,j)]
        end
    end
    for i in 1:remain_items
        res = [res..., gen_one(date,rand(1:11),rand(1:each_domain_items))]
    end
    return res
end

end  # module Lockyv2


text:
function frqtr!(wc::AbstractVector, c::AbstractVector, α,
                prev::Vector=Array{eltype(wc),1}(length(wc)))
    fill!(wc, zero(eltype(wc)))
    dst_order = length(wc) - 1

    m1 = length(c)-1
    for i in -m1:0
        copy!(prev, wc)
        if dst_order >= 0
            @inbounds wc[1] = c[-i+1]
        end
        for m=2:dst_order+1
            @inbounds wc[m] = prev[m-1] + α*(prev[m] - wc[m-1])
        end
    end

    wc
end

function frqtr(c::AbstractVector, order=25, α=0.35)
    frqtr!(Array{eltype(c),1}(order+1), c, α)
end


text:
using LibSerialPort

v_uint8_eoc = [0xff, 0xff, 0xff]  # Array{UInt8,1}
v_char_eoc = Char.(v_uint8_eoc)  # Array{Char,1} ['ÿ', 'ÿ', 'ÿ']
s_eoc = String(copy(v_uint8_eoc))  # "\xff\xff\xff"

function _format_command(cmd)
    cmd * s_eoc
end

port = "/dev/ttyUSB0"
baudrate = 9600

sp = open(port, baudrate)

cmd = "sendme"
cmd = _format_command(cmd)
write(sp, cmd)

function my_readuntil(sp::SerialPort, delim::AbstractString, timeout_ms::Real)
   return String(take!(readuntil(sp, Vector{Char}(delim), timeout_ms)))
end

function my_readuntil(sp::SerialPort, delim::Vector{T}, timeout_ms::Real) where {T}
    start_time = time_ns()
    out = Vector{T}()
    lastchars = T[0 for i=1:length(delim)]
    while !eof(sp)
        if (time_ns() - start_time)/1e6 > timeout_ms
            break
        end
        if bytesavailable(sp) > 0
            c = read(sp, T)
            push!(out, c)
            lastchars = circshift(lastchars,-1)
            lastchars[end] = c
            if lastchars == delim
                break
            end
        end
    end
    return out
end

timeout_ms = 4000

#s = readuntil(sp, s_eoc, timeout_ms)
#println(s)
#r = transcode(UInt8, s)
#println(r)

#r = my_readuntil(sp, v_char_eoc, timeout_ms)
#println(r)

r = my_readuntil(sp, v_uint8_eoc, timeout_ms)
println(r)


text:
using PyPlot
using DelimitedFiles

if !isdir("figures_summary")
  mkdir("figures_summary")
end

m = 15
n = 30
h = 30.0
dz = 3.0 # meters
dx = 3.0
nz = Int64(round((m * h) / dz)) + 1
nx = Int64(round((n * h) / dx)) + 1
z_src = (collect(5:10:nz-5) .- 1 ) .* dz .+ dz/2.0
x_src = (5-1)ones(Int64, size(z_src)) .* dx .+ dx/2.0
z_rec = (collect(5:1:nz-5) .- 1) .* dz .+ dz/2.0
x_rec = (nx-5-1) .* ones(Int64, size(z_rec)) .*dx .+ dx/2.0

z_inj = (9-1)*h + h/2.0
x_inj = (3-1)*h + h/2.0
z_prod = (9-1)*h + h/2.0
x_prod = (28-1)*h + h/2.0

iter = 100
Prj_names = ["CO2_patchy_pgs"]
K_name = "/K$iter.txt"

rc("axes", titlesize=20)
rc("axes", labelsize=18)
rc("xtick", labelsize=18)
rc("ytick", labelsize=18)
rc("legend", fontsize=20)

for indStage = 2:11
  figure()
  iPrj = 1
  K = readdlm(Prj_names[iPrj] * "/Stage$indStage/" * K_name)
  imshow(K, extent=[0,n*h,m*h,0]);
  xlabel("Distance (m)")
  ylabel("Depth (m)")
  cb = colorbar()
  clim([20, 120])
  cb.set_label("Permeability (md)")
  shot_inds = collect(1:length(z_src))
  scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
  scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
  scatter(x_inj, z_inj, c="r", marker=">")
  scatter(x_prod, z_prod, c="r", marker="<")
  savefig("figures_summary/K_$(Prj_names[iPrj])_stage_$indStage.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);
end


text:
export QuantumCircuitModel

"""
    QCModelOptions
The composite mutable struct, `QCModelOptions`, holds various optimization model options for enhancements
with defualt options set to the values provided by `get_default_options` function.
"""
mutable struct QCModelOptions
    model_type                         :: String

    all_valid_constraints              :: Int64
    commute_gate_constraints           :: Bool
    involutory_gate_constraints        :: Bool
    redundant_gate_pair_constraints    :: Bool
    identity_gate_symmetry_constraints :: Bool
    visualize_solution                 :: Bool
    idempotent_gate_constraints        :: Bool
    convex_hull_gate_constraints       :: Bool
    unit_magnitude_constraints         :: Bool

    time_limit                         :: Float64
    relax_integrality                  :: Bool
    optimizer_log                      :: Bool
    objective_slack_penalty            :: Float64
end

"""
    get_default_options()
This function returns the default options for building the struct `QCModelOptions`.
"""
function get_default_options()
    model_type                         = "compact_formulation" # compact_formulation, balas_formulation

    all_valid_constraints              = 0    # -1, 0, 1
    commute_gate_constraints           = true # true, false
    involutory_gate_constraints        = true # true, false
    redundant_gate_pair_constraints    = true # true, false
    identity_gate_symmetry_constraints = true # true, false
    visualize_solution                 = true # true, false

    idempotent_gate_constraints        = false # true, false
    convex_hull_gate_constraints       = false # true, false
    unit_magnitude_constraints         = false # true, false

    time_limit                         = 10800 # float value
    relax_integrality                  = false # true, false
    optimizer_log                      = true  # true, false
    objective_slack_penalty            = 1E3   # > 0 value

    return QCModelOptions(model_type,
                          all_valid_constraints,
                          commute_gate_constraints,
                          involutory_gate_constraints,
                          redundant_gate_pair_constraints,
                          identity_gate_symmetry_constraints,
                          visualize_solution,
                          idempotent_gate_constraints,
                          convex_hull_gate_constraints,
                          unit_magnitude_constraints,
                          time_limit,
                          relax_integrality,
                          optimizer_log,
                          objective_slack_penalty)
end

"""
    QuantumCircuitModel
The composite mutable struct, `QuantumCircuitModel`, holds dictionaries for input data, abstract JuMP model for optimization,
variable references and result from solving the JuMP model.
"""
mutable struct QuantumCircuitModel
    data      :: Dict{String,Any}
    model     :: JuMP.Model
    options   :: QCModelOptions
    variables :: Dict{Symbol,Any}
    result    :: Dict{String,Any}

    "Constructor for struct `QuantumCircuitModel`"
    function QuantumCircuitModel(data::Dict{String,Any})
        data      = data
        model     = JuMP.Model()
        options   = QCO.get_default_options()
        variables = Dict{Symbol,Any}()
        result    = Dict{String,Any}()

        qcm = new(data, model, options, variables, result)

        return qcm
    end

end

"""
    GateData
The composite mutable struct, `GateData`, type of the gate, the complex matrix form
of the gate, full sized real form of the gate, inverse of the gate and a boolean which
states if the gate has all real entries.
"""
mutable struct GateData
    type    :: String
    complex :: Array{Complex{Float64},2}
    real    :: Array{Float64,2}
    inverse :: Array{Float64,2}
    isreal  :: Bool

    "Constructor for struct `GateData`"
    function GateData(gate_type::String, num_qubits::Int64)
        type    = gate_type
        complex = QCO.get_full_sized_gate(type, num_qubits)
        real    = QCO.complex_to_real_gate(complex)
        inverse = inv(real)
        isreal  = iszero(imag(complex))

        gate = new(type, complex, real, inverse, isreal)

        return gate
    end

end


text:
import ColorSchemes: viridis, inferno, twilight, deep, matter, ice, algae,  balance, curl
import Colors.hex
function plotly_cs(colorscheme; n_entries=11)
    scale = LinRange(0, 1, n_entries)
    colors = [get(colorscheme, s) for s in scale]
    return [[s, hex(color)] for (s, color) in zip(scale, colors)]
end

function normalize(x, a, b) #maps  the interval [a,b]  to [0,1] and finds the val
    (a < b) & (a <= x <= b) || error("($a, $b) is not an interval or $x not in ($a, $b)")
    return (x-a)/(b-a)
end

function asymmetric_div_scheme(data:: AbstractArray,  div_cs; ref_point=0.0, N=13)
    # defines an asymettric colorscheme adaptated to data, from a diverging one, and a reference point
    # data: abstract array
    # div_cmap: symmetric diverging colorscheme
    # ref_point:  reference point, i.e. a value in data range to be colormapped to the
    #             color of highest lightness in the diverging colorscheme
    # N:  number of entries in the asymmetric Plotly asymmetric colorscale

    dataMin, dataMax = extrema(filter(x->!isnan(x), data))
    (dataMin < ref_point < dataMax) ||
        error("The reference point, $ref_point is not in ($dataMin, $dataMax)")

    if dataMax-ref_point > ref_point - dataMin
        left_val = 2*ref_point-dataMax
        right_val = dataMax

        s = normalize(dataMin, left_val,right_val)
        refP_norm = normalize(ref_point, left_val, right_val)
        n = Int64(round(N*(refP_norm-s)))

        t = LinRange(s, refP_norm, n)
        step = t[2]-t[1]
        t = vcat(t, LinRange(refP_norm+step, 1.0, N-n))

    else
        left_val = dataMin
        right_val = 2*ref_point-dataMin
        s = normalize(dataMax, left_val,right_val)
        refP_norm = normalize(ref_point, left_val, right_val)
        n = Int64(round(N*refP_norm))
        t = LinRange(0, refP_norm, n)
        step = t[2]-t[1]
        t= vcat(t, LinRange(refP_norm+step, s, N-n) )
    end

    tn = [normalize(tv, t[1], t[end]) for tv in t]
    colors = [get(div_cs, tv) for tv in t]
    return [[tn[k], hex(color)] for (k, color) in enumerate(colors)]
end

pl_hsv=[[0.0, "rgb(0, 242, 242)"], #Saturation=1, Value(Brightness)=0.95,
 [0.083, "rgb(0, 121, 242)"],
 [0.166, "rgb(0, 0, 242)"],
 [0.25, "rgb(121, 0, 242)"],
 [0.333, "rgb(242, 0, 242)"],
 [0.416, "rgb(242, 0, 121)"],
 [0.5, "rgb(242, 0, 0)"],
 [0.583, "rgb(242, 121, 0)"],
 [0.666, "rgb(242, 242, 0)"],
 [0.75, "rgb(121, 242, 0)"],
 [0.833, "rgb(0, 242, 0)"],
 [0.916, "rgb(0, 242, 121)"],
 [1.0, "rgb(0, 242, 242)"]];

pl_ehtcmap = [[0.0, "#342724"], #plotly version of ehtcmap from ehtplot
 [0.1, "#682018"],
 [0.2, "#8e1d15"],
 [0.3, "#a82c21"],
 [0.4, "#c33d26"],
 [0.5, "#da560f"],
 [0.6, "#e87500"],
 [0.7, "#f49600"],
 [0.8, "#fbb71c"],
 [0.9, "#fadb74"],
 [1.0, "#fef7e7"]];

pl_BrBG=[[0.0,"#543005"],
[0.1,"#8c510a"],
[0.2, "#bf812d"],
[0.3, "#dfc27d"],
[0.4, "#f6e8c3"],
[0.5, "#f5f5f5"],
[0.6, "#c7eae5"],
[0.7, "#80cdc1"],
[0.8, "#35978f"],
[0.9, "#01665e"],
[1.0, "#003c30"]];

pl_viridis = plotly_cs(viridis, n_entries=13)
pl_inferno = plotly_cs(inferno)
pl_twilight = plotly_cs(twilight, n_entries=21)
pl_deep = plotly_cs(deep)
pl_matter = plotly_cs(matter)
pl_algae = plotly_cs(algae)
pl_ice = plotly_cs(ice)
pl_balance = plotly_cs(balance) #diverging cs
pl_curl = plotly_cs(curl); #diverging cs


text:
using Documenter, TraitSimulation

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = Documenter.HTML(),
    sitename = "TraitSimulation",
    modules = [TraitSimulation]
)

deploydocs(
    repo   = "github.com/sarah-ji/TraitSimulation.jl.git",
    target = "build"
)

max_stars_repo_name: sarah-ji/Tutorials

max_issues_repo_name: sarah-ji/Tutorials

max_forks_repo_name: sarah-ji/Tutorials


text:
using Cat
using Random

@category Model

@arrow Model Normal :: (Float64, Float64) --> Float64
@arrow Model Uniform :: Nothing --> Float64

# A functor Model => Set that composes as a state monad
@interpretation Sample (=>) Model begin
    rng::AbstractRNG
    N::Int64
    samples::Dict{Model.Arrow, Any}
    Sample(N) = new(Random.GLOBAL_RNG, N, Dict())
end

"Hook for common state update rules; called when @interpret is used"
function interp_state_hook(s::Sample, m::Model.Arrow, value_expr)
  if !haskey(s.samples, m)
      s.samples[m] = value_expr()
  end
  s.samples[m]
end

"Sample from a normal"
@interpret function (s::Sample)(m::Normal, μ, σ)
    μ .+ σ.*randn(s.rng, s.N)
end
"Sample from a uniform"
@interpret function (s::Sample)(m::Uniform, _)
    rand(s.rng, s.N)
end

z = Normal(1.0, Uniform())
y = Normal(z, 2.0)
x = Normal(z, y)

samples = Sample(5000)(x, nothing)

μ = sum(samples)/length(samples)
σ = sqrt( sum( (samples .- μ).^2 ) / length(samples) )

println("Sample μ/σ: $μ, $σ")

@interpretation Filter (<=) Model begin
    rng::AbstractRNG
    N::Int64
end

@interpret function (s::Sample)(m::Normal, z)

end


text:
function draw_gantt_text(ctx; delay=2, width=40, height=20, window_length=20)
    root_time = time_ns()
    window_logs = []
    if window_length !== nothing
        window_start = -window_length
    else
        window_start = 0
    end
    window_finish = 0
    procs = Set()
    for gp in ctx.procs
        for proc in Dagger.get_processors(gp)
            push!(procs, proc)
        end
    end

    while continue_rendering[]
        sleep(delay)
        logs = Dagger.get_logs!(ctx.log_sink, true) # get raw events
        update_window_logs!(window_logs, logs; root_time=root_time, window_start=window_start)

        #= TODO: Concatenate and render profile data
        if ctx.profile
            prof_data, prof_lidict = logs_to_stackframes(window_logs)
            prof_to_svg(prof_path, prof_data, prof_lidict, image_idx; width=width)
        end
        =#

        for proc in unique(map(x->x[1].timeline[2], filter(x->x[1].category==:compute, window_logs)))
            push!(procs, proc)
        end
        procs_len = length(procs)
        print("\e[2J\e[1;1H")
        for (proc_idx, proc) in enumerate(sort(collect(procs); lt=proclt))
            proc_row = ((proc_idx-1)*2)+1
            print("\e[$(proc_row)H")
            print("$(getname(proc)), $(window_start) s - $(window_finish) s")
            for log in filter(x->x[1].timeline[2]==proc, filter(x->x[1].category==:compute, window_logs))
                log_start_s = (log[1].timestamp-root_time)/(1000^3)
                log_finish_s = if length(log) == 2
                    log_finish_s = (log[2].timestamp-root_time)/(1000^3)
                else
                    window_length+window_start
                end
                log_start_s > (window_length+window_start) && continue # skip premature entries
                log_finish_s < window_start && continue # skip finished entries
                xstart = floor(Int, ((log_start_s-window_start)/(window_finish-window_start))*width)
                xfinish = ceil(Int, ((log_finish_s-window_start)/(window_finish-window_start))*width)
                xstart = max(xstart, 1)
                xfinish = min(xfinish, width)
                print("\e[$(proc_row+1);$(xstart)H")
                print(repeat("X", xfinish-xstart))
            end
            println()
        end
        #= TODO: Scheduler init
        for log in filter(x->x[1].category==:scheduler_init, window_logs)
            log_start_s = (log[1].timestamp-root_time)/(1000^3)
            log_finish_s = if length(log) == 2
                log_finish_s = (log[2].timestamp-root_time)/(1000^3)
            else
                window_finish+1
            end
            xstart = ((log_start_s-window_start)/(window_finish-window_start))*width
            xfinish = ((log_finish_s-window_start)/(window_finish-window_start))*width
            sethue("red")
            rect(Point(xstart,0),#=xfinish-xstart=#1,height)
        end
        =#
        window_finish = (time_ns()-root_time)/(1000^3)
        if window_length !== nothing
            window_start = window_finish-window_length
        end
    end
end


text:
using Documenter, GeoMakie

# use the README as the home page
readme = read(joinpath(dirname(@__DIR__), "README.md"), String)
# HTML not supported in the Markdown parser used by Documenter
href = """"""
homepage = replace(readme, href => "Makie.jl")
write(joinpath(@__DIR__, "src", "index.md"), homepage)

makedocs(;
    modules=[GeoMakie],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Example" => "example.md",
        "Recipes" => [
            "stock.md"
        ],
        "API Reference" => "api.md"
    ],
    sitename="GeoMakie.jl",
    authors="Anshul Singhvi",
)

deploydocs(;
    repo="github.com/JuliaPlots/GeoMakie.jl",
    push_preview = true,
)


text: for i in 1:10 println("Door $(i^2) is open.") end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Test VoidStopping" begin
  stp = VoidStopping()

  @test typeof(stp) <: AbstractStopping

  io = IOBuffer()
  show(io, stp)
end


text:
for num in 2:1499
	prime = true
	for i in 2:num-1
		if num % i == 0
			prime = false
		end
	end
	if prime

	end
end

max_stars_repo_name: bichanna/slap

max_issues_repo_name: bichanna/slap

max_forks_repo_name: bichanna/slap


text:
include("ekfslam-sim.jl")
include("gr-draw.jl")

# Global object containing simulation state, updated in real time.
# It is exposed for monitoring, debugging, and flexible visualization.
simdata = ekfsim_setup(10, "course1.txt")

init_plot_window(simdata.scene.boundaries)
draw_map(simdata.scene.landmarks, simdata.scene.waypoints)

let ellipses = []

    global monitor

    function monitor(simdata::SimData)
        scene, vehicle, state = simdata.scene, simdata.vehicle, simdata.state
        draw_scene(scene, state, vehicle)
        if simdata.state_updated
            draw_laser_lines(simdata.z, state.x[1:3])
            ellipses = compute_landmark_ellipses(state.x, state.cov)
        end
        draw_landmark_ellipses(ellipses)
        GR.updatews()
        return
    end

end

sim!(simdata, monitor, [simdata])


text:
@testset "Crazy Types as Loop Invariants" begin
    a = 0
    @test_throws TypeError begin
        @loopinvariant 2 for i = 1:10
            a += 1
        end
    end

    @test_throws TypeError begin
        @loopinvariant :sym for i = 1:10
            a += 1
        end
    end

    @test_throws TypeError begin
        @loopinvariant (x -> 10 + x) for i = 1:10
            a += 1
        end
    end

    @test begin
        @loopinvariant true for i = 1:10
            a += 1
        end
        a == 10
    end

    a = 0
    @test begin
        @loopinvariant all((x -> 10 + x < 100).([a, 10, 20, 30, a])) for i = 1:10
            a += 1
        end
        a == 10
    end
end


text:
module FEIterators

using StaticArrays
using MeshCore
using MeshCore: nshapes, indextype, nrelations, nentities, IncRel, shapedesc
using MeshSteward: Mesh, baseincrel, increl
using ..RefShapes: manifdim, manifdimv
using ..FElements: refshape, nfeatofdim
import ..FElements: jacjac
using ..FESpaces.FEFields: FEField, ndofsperterm
using ..FESpaces: FESpace, doftype
import ..FESpaces: ndofsperel
using ..QPIterators: QPIterator, bfungradpar

#= TODO is it more natural to have access to the geometry from the font element space or from the iterator? =#

"""
    FEIterator{FES, IR, G, IT, T, V, IR0, IR1, IR2, IR3, F0, F1, F2, F3}

Type of finite element iterator. Parameterized with the types of
- `FES`: finite element space,
- `IR`: base incidence relation of the mesh,
- `G`: type of the geometry attribute,
- `IT`: type of integer indices, such as the  numbers of nodes and degrees of freedom,
- `T`: type of the degree of freedom value (real double, complex float, ... ),
- `V`: `Val` representation of the manifold dimension of the base relation elements,
- `IR0`, `IR1`, `IR2`, `IR3`: types of incidence relations with which degrees
  of freedom are associated in the finite element space, for each of the
  manifolds dimensions 0, 1, 2, 3,
- `F0`, `F1`, `F2`, `F3`: types of fields with which degrees
  of freedom are associated in the finite element space, for each of the
  manifolds dimensions 0, 1, 2, 3.
"""
struct FEIterator{FES, IR, G, IT, T, V, IR0, IR1, IR2, IR3, F0, F1, F2, F3}
    fesp::FES
    _bir::IR
    _geom::G
    _dofs::Vector{IT}
    _dofvals::Vector{T}
    _entmdim::Vector{IT}
    _dofcomp::Vector{IT}
    _nodes::Vector{IT}
    _irs0::Union{Nothing, IR0}
    _irs1::Union{Nothing, IR1}
    _irs2::Union{Nothing, IR2}
    _irs3::Union{Nothing, IR3}
    _fld0::Union{Nothing, F0}
    _fld1::Union{Nothing, F1}
    _fld2::Union{Nothing, F2}
    _fld3::Union{Nothing, F3}
    _manifdimv::V
    _state::Vector{IT}

    function FEIterator(fesp::FES) where {FES}
        _bir = baseincrel(fesp.mesh)
        _geom = MeshCore.attribute(_bir.right, "geom")
        _dofs = zeros(Int64, ndofsperel(fesp))
        _dofvals = zeros(doftype(fesp), ndofsperel(fesp))
        _entmdim = zeros(Int64, ndofsperel(fesp))
        _dofcomp = zeros(Int64, ndofsperel(fesp))
        _nodes = zeros(Int64, nfeatofdim(fesp.fe, 0))
        _fld0 = nothing
        _fld1 = nothing
        _fld2 = nothing
        _fld3 = nothing
        _irs0 = nothing
        _irs1 = nothing
        _irs2 = nothing
        _irs3 = nothing
        0 in keys(fesp._irsfields) && (_irs0 = fesp._irsfields[0][1]; _fld0 = fesp._irsfields[0][2])
        1 in keys(fesp._irsfields) && (_irs1 = fesp._irsfields[1][1]; _fld1 = fesp._irsfields[1][2])
        2 in keys(fesp._irsfields) && (_irs2 = fesp._irsfields[2][1]; _fld2 = fesp._irsfields[2][2])
        3 in keys(fesp._irsfields) && (_irs3 = fesp._irsfields[3][1]; _fld3 = fesp._irsfields[3][2])
        _manifdimv = Val(manifdim(refshape(fesp.fe)))
        p = 1
        _irs0 != nothing && (p = _init_e_d!(_entmdim, _dofcomp, p, _irs0, _fld0))
        _irs1 != nothing && (p = _init_e_d!(_entmdim, _dofcomp, p, _irs1, _fld1))
        _irs2 != nothing && (p = _init_e_d!(_entmdim, _dofcomp, p, _irs2, _fld2))
        _irs3 != nothing && (p = _init_e_d!(_entmdim, _dofcomp, p, _irs3, _fld3))
        return new{FES, typeof(_bir), typeof(_geom), eltype(_dofs), doftype(fesp), typeof(_manifdimv), typeof(_irs0), typeof(_irs1), typeof(_irs2), typeof(_irs3), typeof(_fld0), typeof(_fld1), typeof(_fld2), typeof(_fld3)}(fesp, _bir, _geom, _dofs, _dofvals, _entmdim, _dofcomp, _nodes, _irs0, _irs1, _irs2, _irs3, _fld0, _fld1, _fld2, _fld3, _manifdimv, [0])
    end
end

"""
    Base.iterate(it::FEIterator, state = 1)

Advance the iterator to the next entity.

The nodes of the finite element are cached, as is a vector of all the degrees
of freedom represented on the element.
"""
function Base.iterate(it::FEIterator, state = 1)
    if state > nrelations(it._bir)
        return nothing
    else
        return (_update!(it, state), state+1)
    end
end

"""
    Base.length(it::FEIterator)

Number of elements represented by this iterator.
"""
Base.length(it::FEIterator)  = nrelations(it._bir)

"""
    ndofsperel(it::FEIterator)

Retrieve the number of degrees of freedom per element.
"""
ndofsperel(it::FEIterator) = ndofsperel(it.fesp)

"""
    eldofs(it::FEIterator)

Retrieve the vector of the element degrees of freedom
"""
eldofs(it::FEIterator) = it._dofs

"""
    elnodes(it::FEIterator)

Retrieve the vector of the nodes of the element.
"""
elnodes(it::FEIterator) = it._nodes

"""
    eldofentmdims(it::FEIterator)

Retrieve the vector of the entity dimensions for each element degree of freedom.

Each degree of freedom is associated with some entity of the finite element:
vertices, edges, faces, and so on. This vector records the dimension of the
manifold entity with which each degree of freedom is associated.
"""
eldofentmdims(it::FEIterator) = it._entmdim

"""
    eldofcomps(it::FEIterator)

Retrieve the vector of the component numbers for each element degree of freedom.

If multiple copies of the finite element are referenced in the finite element
space, each copy is referred to as component.
"""
eldofcomps(it::FEIterator) = it._dofcomp

function _init_e_d!(mdim, dofcomp, p, ir, fl)
    ndpt = ndofsperterm(fl)
    for k in 1:nentities(ir, 1)
        for i in 1:ndpt
            mdim[p] = MeshCore.manifdim(shapedesc(ir.right))
            dofcomp[p] = i
            p = p + 1
        end
    end
    return p
end

function _storedofs!(d, p, e, ir, fl)
    ndpt = ndofsperterm(fl)
    for k in 1:nentities(ir, e)
        gk = ir[e, k]
        for i in 1:ndpt
            d[p] = fl.dofnums[gk][i]
            p = p + 1
        end
    end
    return p
end

function _storedofvals!(d, p, e, ir, fl)
    ndpt = ndofsperterm(fl)
    for k in 1:nentities(ir, e)
        gk = ir[e, k]
        for i in 1:ndpt
            d[p] = fl.dofvals[gk][i]
            p = p + 1
        end
    end
    return p
end

function _update!(it::FEIterator, state)
    it._state[1] = state
    copyto!(it._nodes, it._bir[state])
    p = 1
    it._irs0 != nothing && (p = _storedofs!(it._dofs, p, state, it._irs0, it._fld0))
    it._irs1 != nothing && (p = _storedofs!(it._dofs, p, state, it._irs1, it._fld1))
    it._irs2 != nothing && (p = _storedofs!(it._dofs, p, state, it._irs2, it._fld2))
    it._irs3 != nothing && (p = _storedofs!(it._dofs, p, state, it._irs3, it._fld3))
    return it
end

"""
    eldofvals(it::FEIterator)

Provide access to vector of element degrees of freedom.
"""
function eldofvals(it::FEIterator)
    p = 1
    it._irs0 != nothing && (p = _storedofvals!(it._dofvals, p, it._state[1], it._irs0, it._fld0))
    it._irs1 != nothing && (p = _storedofvals!(it._dofvals, p, it._state[1], it._irs1, it._fld1))
    it._irs2 != nothing && (p = _storedofvals!(it._dofvals, p, it._state[1], it._irs2, it._fld2))
    it._irs3 != nothing && (p = _storedofvals!(it._dofvals, p, it._state[1], it._irs3, it._fld3))
    return it._dofvals
end

"""
    jacjac(it::FEIterator, qpit::QPIterator)

Compute the Jacobian matrix and the Jacobian determinant.

The finite element iterator cooperates with the quadrature point iterator here
to compute the Jacobian at the current integration point.
"""
function jacjac(it::FEIterator, qpit::QPIterator)
    return jacjac(it.fesp.fe, it._geom, it._nodes, qpit._geomscalbfungrad_ps[qpit._pt])
end

"""
    location(it::FEIterator, qpit::QPIterator)

Calculate the location of the quadrature point.
"""
function location(it::FEIterator, qpit::QPIterator)
    n = it._nodes[1]
    loc = it._geom[n] * qpit._geomscalbfuns[qpit._pt][1]
    for i in 2:length(it._nodes)
        n = it._nodes[i]
        loc += it._geom[n] * qpit._geomscalbfuns[qpit._pt][i]
    end
    return loc
end

end


text:
# 	This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith

export string,
       show

import Base: string,
             show

string(m::MExpr) = convert(String, m)
show(io::IO, m::MExpr) = print(io, convert(String, m))

function show(io::IO, ::MIME"text/plain", m::MExpr)
    input = "'(" * replace(convert(String, m), r";" => ");\n'(") * ")"
    write(ms.input, "$(replace(input, r";" => "\$"))\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        print(io, out)
    else
        mcall(m"display2d: true")
        write(ms, replace(input, r";" => "; print(ascii(3))\$ "))
        str = read(ms)
        str = rstrip(str, '\n')
        mcall(m"display2d: false")
        sp = split(str, "\x03")
        for k in 1:length(sp)
            print(io, sp[k])
        end
    end
end

function show(io::IO, ::MIME"text/latex", m::MExpr)
    check = "'(" * replace(convert(String, m), r";" => ")\$\n'(") * ")"
    write(ms.input, "$check\$\n print(ascii(4))\$")
    out = (readuntil(ms.output, EOT) |> String
                                     |> str -> rstrip(str, EOT))
    if occursin(synerr, out) || occursin(maxerr, out)
        @warn "Invalid Maxima expression"
        print(io, out)
    else
        write(ms, "tex('("*replace(convert(String, m), r";" => "))\$\ntex('(") * "))")
        texstr = read(ms)
        print(io, replace(texstr, r"\nfalse\n" => ""))
    end
end


text:
#include("./WebSockets.jl/src/WebSockets.jl")
using WebSockets
#import HTTP.Servers: MbedTLS

#import HTTP.Servers: MbedTLS                # For further imports
#import HTTP.Servers.MbedTLS:SSLConfig

#using HTTP, MbedTLS
#client = HTTP.Client(tlsconfig=MbedTLS.SSLConfig(false))
#HTTP.get(client, url; options...)

#sslconf = WebSockets.SSLConfig(cert, key)
#ServerWS(h,w, sslconfig = sslconf)

#tlsconfig = MbedTLS.SSLConfig("cert.pem", "key.pem")
#tlsconfig = WebSocket.SSLConfig("cert.pem", "cert_key.pem")

#require_ssl_verification = false
#HTTP.request("GET", "https://localhost:8080", require_ssl_verification = true)

#print(WebSocket)
#print(WebSCoket)

#WebSockets.default_options().sslconfig = tlsconfig

Base.exit_on_sigint(false)

try
  WebSockets.open("wss://127.0.0.1:8080/ws", require_ssl_verification=false) do ws
    writeguarded(ws, "hej")
    msg, stillopen = readguarded(ws)
    println("Received:", String(msg))
    if stillopen
      println("The connection is active, but we leave. WebSockets.jl will close properly.")
    else
      println("Disconnect during reading.")
    end
  end
catch err
  typeof(err) == InterruptException && rethrow(err)
  print(err)
end
#  data, success = readguarded(ws_client)
#  if success
#    println("received:", String(data))
#  end


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

mutable struct ApplicationGatewayRequestRoutingRulePropertiesFormat <: SwaggerModel
    ruleType::Any # spec type: Union{ Nothing, String } # spec name: ruleType
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    backendHttpSettings::Any # spec type: Union{ Nothing, SubResource } # spec name: backendHttpSettings
    httpListener::Any # spec type: Union{ Nothing, SubResource } # spec name: httpListener
    urlPathMap::Any # spec type: Union{ Nothing, SubResource } # spec name: urlPathMap
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewayRequestRoutingRulePropertiesFormat(;ruleType=nothing, backendAddressPool=nothing, backendHttpSettings=nothing, httpListener=nothing, urlPathMap=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("ruleType"), ruleType)
        setfield!(o, Symbol("ruleType"), ruleType)
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("backendHttpSettings"), backendHttpSettings)
        setfield!(o, Symbol("backendHttpSettings"), backendHttpSettings)
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("httpListener"), httpListener)
        setfield!(o, Symbol("httpListener"), httpListener)
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("urlPathMap"), urlPathMap)
        setfield!(o, Symbol("urlPathMap"), urlPathMap)
        validate_property(ApplicationGatewayRequestRoutingRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayRequestRoutingRulePropertiesFormat

const _property_map_ApplicationGatewayRequestRoutingRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("ruleType")=>Symbol("ruleType"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("backendHttpSettings")=>Symbol("backendHttpSettings"), Symbol("httpListener")=>Symbol("httpListener"), Symbol("urlPathMap")=>Symbol("urlPathMap"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayRequestRoutingRulePropertiesFormat = Dict{Symbol,String}(Symbol("ruleType")=>"String", Symbol("backendAddressPool")=>"SubResource", Symbol("backendHttpSettings")=>"SubResource", Symbol("httpListener")=>"SubResource", Symbol("urlPathMap")=>"SubResource", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewayRequestRoutingRulePropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayRequestRoutingRulePropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayRequestRoutingRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayRequestRoutingRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayRequestRoutingRulePropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayRequestRoutingRulePropertiesFormat[property_name]

const _allowed_ApplicationGatewayRequestRoutingRulePropertiesFormat_ruleType = ["Basic", "PathBasedRouting"]

function check_required(o::ApplicationGatewayRequestRoutingRulePropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayRequestRoutingRulePropertiesFormat }, name::Symbol, val)
    if name === Symbol("ruleType")
        Swagger.validate_param(name, "ApplicationGatewayRequestRoutingRulePropertiesFormat", :enum, val, _allowed_ApplicationGatewayRequestRoutingRulePropertiesFormat_ruleType)
    end
end


text:
mutable struct MovingWindowSolver{𝒫_T, τ_T, Δτ_T} <: InitialValueSolver
    𝒫::𝒫_T
    τ::τ_T
    Δτ::Δτ_T
end

# function MovingWindowSolver(parallelsolver::TimeParallelSolver, τ, Δτ)
#     if Δτ < τ / parallelsolver.P
#         error("Select Δτ ≥ τ / P!")
#     elseif Δτ > τ
#         error("Select Δτ ≤ τ!")
#     end
#     new(parallelsolver, τ, Δτ)
# end

MovingWindowSolver(𝒫; τ, Δτ) = MovingWindowSolver(𝒫, τ, Δτ)
@doc (@doc MovingWindowSolver) MoWA(args...; kwargs...) = MovingWindowSolver(args...; kwargs...)


text:
#! /usr/bin/julia

# Rosetta Code, Check input device is a terminal

if isa(STDIN, Base.TTY)
    println("This program sees STDIN as a TTY.")
else
    println("This program does not see STDIN as a TTY.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# utility functions

"""
`update_bin_weights!`: Update the bin weights from the ensemble

### Arguments
* `B` - bin data structure
* `E` - particle ensemble
"""
function update_bin_weights!(B::TB, E::TE) where {TB<:Bins, TE<:Ensemble}

   n_bins = length(B);

   # this loops over bins
   for i in 1:n_bins
      particle_ids = findall(isequal(i), E.b);

      B.ν[i] = sum(E.ω[particle_ids]);
      B.n[i] = length(particle_ids);
   end
   B
end

"""
`Points_to_Bins`: Convenience function for constructing a bin structure using a
sequence of points

### Arguments
`points` - An array of points defining the Voronoi cells
"""
function Points_to_Bins(points; d::TD = Nothing()) where {TD}

   B = Bins{typeof(points[1]), Float64, Int, Int, TD}([],[],[],[],[]);
   for point in points
      push!(B, point , 0, 0, 0, deepcopy(d));
   end
   return B
end

"""
`Voronoi_to_Bins`: Convenience function for constructing a bin structure using a
sequence of points as the Voronoi sites

### Arguments
`sites` - An array of points defining the Voronoi cells
"""
function Voronoi_to_Bins(sites; d::TD = Nothing()) where {TD}

   B = Bins{typeof(sites[1]), Float64, Int, Int, TD}([],[],[],[],[]);
   for site in sites
      push!(B, site , 0, 0, 0, deepcopy(d));
   end
   return B
end

"""
`Voronoi_bin_id`: Convenience function for bin id in Voronoi based bins

### Arguments
`X` - Point thats bin is to be determined
`tree` - A nearest neighbors tree structure constructed with `KDTree`
"""
function Voronoi_bin_id(X, tree)
   return knn(tree, X, 1)[1][1]
end

"""
`setup_Voronoi_bins` - Convenience function for constructing bins, a bin id function,
and a rebinning function based on a set of Voronoi points

### Arguments
`voronoi_pts` - User specified Voronoi cell centers
"""
function setup_Voronoi_bins(voronoi_pts; d::TD = Nothing()) where {TD}
   # construct the bins
   B = Voronoi_to_Bins(voronoi_pts, d = d);
   # construct the bin id function
   tree = KDTree(hcat(voronoi_pts...));
   bin_id = x-> Voronoi_bin_id(x,tree);
   # define the rebinning function
   function rebin!(E, B, t)
       @. E.b = bin_id(E.ξ);
      update_bin_weights!(B, E);
      E, B
   end

   return B, bin_id, rebin!

end

"""
`Dirac_to_Ensemble`: Convenience function for construction an ensemble from a
single initial walker.  This hard codes the weights to be Float64.

### Arguments
`X` - Starting state of all walkers
`n_particles` - Number of walkers in the ensemble
"""
function Dirac_to_Ensemble(X::TP, n_particles::TI; d::TD = Nothing()) where {TP, TI<:Integer, TD}
   ω = 1.0/n_particles
   E = Ensemble{TP, Float64, TI, TD}([deepcopy(X) for _ in 1:n_particles],
                                    [deepcopy(X) for _ in 1:n_particles],
                                    ω * ones(n_particles),ω * ones(n_particles),
                                    zeros(TI, n_particles),zeros(TI, n_particles),
                                    zeros(TI, n_particles),
                                    [deepcopy(d) for _ in 1:n_particles],
                                    [deepcopy(d) for _ in 1:n_particles]);
   return E
end

function trivial_analysis!(E, B, t)
   E, B
end


text:
using PyCall: python

ENV["TA_LIBRARY_PATH"] = joinpath(DEPOT_PATH[1], "talib/lib")
ENV["TA_INCLUDE_PATH"] = joinpath(DEPOT_PATH[1], "talib/include")
ENV["JULIA_DEPOT_PATH"] = DEPOT_PATH[1]

run(`$python -m pip install --upgrade pip`)
run(`$python -m pip install pandas tsfresh pyarrow`)
if !Sys.iswindows()
    run(`bash talib-install.sh`)
    run(`$python -m pip --no-cache-dir install TA-Lib`)
end


text:
# artist.jl
## https://developer.spotify.com/documentation/web-api/reference/artists/

## https://developer.spotify.com/documentation/web-api/reference/artists/get-artist/
@doc """
# Get an Artist
**Summary**: Get Spotify catalog information for a single artist identified by their unique Spotify ID.

`artist_id` _Required_: The Spotify artist ID. Up to 50 artist ID's can be passed by comma delimiting the ID's

[Reference](https://developer.spotify.com/documentation/web-api/reference/artists/get-artist/)
""" ->
function artist_get(artist_id)
    return spotify_request("artists/$artist_id")
end

## https://developer.spotify.com/documentation/web-api/reference/artists/get-artists-albums/
@doc """
# Get an Artists Albums
**Summary**: Get Spotify catalog information about an artist's albums.

`artist_id` _Required_: The Spotify artist ID. Up to 50 artist ID's can be passed by comma delimiting the ID's\n
`include_groups` _Optional_: A comma-separated list of keywords that will be used to filter the response. If not supoplied, all album types will be returned.
Valid values:
* `album`
* `single`
* `compilation`
* `appears_on`

`country` _Optional_: An ISO 3166-1 alpha-2 country code string. Use this to limit the response to one particular geographical market. Default "US"\n
`limit` _Optional_: The number of album objects to return. Default 20\n
`offset` _Optional_: The index of the first album to return. Default 0\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/artists/get-artists-albums/)
""" ->
function artist_get_albums(artist_id, include_groups="None", country="US", limit=20, offset=0)
    if include_groups == "None"
        return spotify_request("artists/$artist_id?include_groups=$include_groups&country=$country&limit=$limit&offset=$offset")
    else
        return spotify_request("artists/$artist_id?country=$country&limit=$limit&offset=$offset")
    end
end

## https://developer.spotify.com/documentation/web-api/reference/artists/get-artists-top-tracks/
@doc """
# Get an Artist's Top Tracks
**Summary**: Get Spotify catalog information about an artist's top tracks by country.\n

`artist_id` _Required_: The Spotify artist ID.\n
`country` _Optional_: An ISO 3166-1 alpha-2 country code string. Use this to limit the response to one particular geographical market. Default "US"\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/artists/get-artists-top-tracks/)
""" ->
function artist_top_tracks(artist_id, country="US")
    return spotify_request("artists/$artist_id/top-tracks?country=$country")
end

## https://developer.spotify.com/documentation/web-api/reference/artists/get-related-artists/
@doc """
# Get an Artist's Related Artists
**Summary**: Get spotify catalog information about artists similar to a given artist. Similarity is based on analysis of the Spotify community's listening history.

`artist_id` _Required_: The Spotify artist ID.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/artists/get-related-artists/)
""" ->
function artist_get_related_artists(artist_id)
    return spotify_request("artists/$artist_id/related-artists")
end


text:
abstract type Aircraft end

get_mass_props(ac::Aircraft) = ac.mass_props
get_pfm(ac::Aircraft) = ac.pfm
get_aerodynamics(ac::Aircraft) = ac.aerodynamics
get_propulsion(ac::Aircraft) = ac.propulsion

get_name(ac::Aircraft) = "Generic Aricraft"
get_wing_area(ac::Aircraft) = 0.0  # m²
get_wing_span(ac::Aircraft) = 0.0  # m
get_chord(ac::Aircraft) = 0.0  # m
get_arp(ac::Aircraft) = [0, 0, 0]  # m

get_empty_mass_props(ac::Aircraft) = RigidSolid(0.0, zeros(3), zeros(3, 3))
get_payload_mass_props(ac::Aircraft) = RigidSolid(0.0, zeros(3), zeros(3, 3))

function calculate_aircraft(ac::Aircraft, fcs::FCS, aerostate::AeroState,
                            state::State, grav::Gravity; consume_fuel=false)
    aero = get_aerodynamics(ac)
    prop = get_propulsion(ac)
    mass_props = get_mass_props(ac)

    # Calculate propulsion
    prop = calculate_propulsion(prop, fcs, aerostate, state; consume_fuel=consume_fuel)
    # Calculate aerodynamics
    aero = calculate_aerodynamics(ac, aero, fcs, aerostate, state)
    # Calculate gravity forces
    grav_force = get_gravity_body(grav, get_attitude(state)) * mass_props.mass
    grav_pfm = PointForcesMoments(mass_props.cg, grav_force, [0, 0, 0])

    β = get_beta(aerostate)
    α = get_alpha(aerostate)
    pfm = grav_pfm + get_pfm(prop) + get_pfm(aero)

    # TODO: where to discount fuel consumption?
    # Here and produce an aircraft with a mass and grav forces "not coherent"
    # or before calculating grav forces?
    # Calculate mass properties
    mass_props = (
        get_fuel_mass_props(prop) +
        get_empty_mass_props(ac) +
        get_payload_mass_props(ac)
        )
    # Return aircraft object
    typeof(ac)(mass_props, pfm, aero, prop)
end


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 = "wget2"
version = v"2.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/wget/wget2-$(version).tar.gz", "4fe2fba0abb653ecc1cc180bea7f04212c17e8fe05c85aaac8baeac4cd241544")
]

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

./configure \
--prefix=${prefix} \
--build=${MACHTYPE} \
--host=${target} \
--enable-shared=yes \
--enable-static=no \
--without-libpsl \
--without-libhsts \
--without-libnghttp2 \
--without-bzip2 \
--without-gpgme  \
--without-lzma \
--without-brotlidec \
--without-lzip \
--without-libidn2 \
--without-libidn \
--without-libpcre2 \
--without-libpcre \
--without-libmicrohttpd \
--without-plugin-support \
--disable-doc \
--disable-manylibs

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 = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("wget2", :wget2),
    #ExecutableProduct("wget2_noinstall", :wget2_noinstall),
    LibraryProduct("libwget", :libwget)
]

# Dependencies that must be installed before this package can be built
#Nettle and OpenSSL needed for mingw builds
dependencies = [
    Dependency(PackageSpec(name="GnuTLS_jll", uuid="0951126a-58fd-58f1-b5b3-b08c7c4a876d"))
    Dependency("Gettext_jll"; compat="=0.21.0")
    Dependency("Nettle_jll"; compat="~3.7.2")
    Dependency("OpenSSL_jll")
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    limit_lengthscale(lap_eig, λ_max)

Only keep modes with length scales larger than minimum
"""
function limit_lengthscale(lap_eig, λ_max)
    (; values, funcs, moments, massrelax) = lap_eig
    inds = values .≤ λ_max
    values = values[inds]
    funcs = funcs[:, inds]
    moments = [m[inds, inds] for m ∈ moments]
    massrelax = massrelax[inds, inds]
    (; values, funcs, moments, massrelax)
end


text:
using Pkg
using Caesar

## Read all test lines from file
fid = open("/tmp/caesar/dbgCSMHexDelayMatrix2.log","r")
lines = readlines(fid)
close(fid)

## split the lines into usable parts
sstr = split.(lines, r" -- ")
resultFolder = String[]
passedStr = String[]
injDelayStr = String[]
for x in sstr
  pic = split(x[2], r", ")
  push!(passedStr, pic[1])
  push!(injDelayStr, pic[2])
  push!(resultFolder, x[1])
end

## which tests were true and false
maskTrue = passedStr .|> px->split(px, '=')[2] == "true"
maskFalse = passedStr .|> px->split(px, '=')[2] == "false"
@assert sum(maskTrue) + sum(maskFalse) == length(lines) "True/False masks must reflect all lines"

# see latest IIF for functions

## True case

csmCounterTrue, trxTrue = calcCSMOccurancesFolders(resultFolder[maskTrue])
maxOccuranTrue = calcCSMOccuranceMax(csmCounterTrue, percentage=false)

open("/tmp/caesar/maxOccTrueTEST.txt","w") do io
  printCSMOccuranceMax(maxOccuranTrue, fid=io, percentage=false)
end

## False case

csmCounterFalse, trxFalse = calcCSMOccurancesFolders(resultFolder[maskFalse])
maxOccuranFalse = calcCSMOccuranceMax(csmCounterFalse)

open("/tmp/caesar/maxOccFalseTEST.txt","w") do io
  printCSMOccuranceMax(maxOccuranFalse, fid=io)
end

using JSON

open("/tmp/caesar/trxTrue2.json", "w") do io
  println(io, JSON.json(trxTrue))
end

open("/tmp/caesar/trxFalse2.json", "w") do io
  println(io, JSON.json(trxFalse))
end

## dev

csmCounterTrue[1][1]

#


text:
using CSV, DataFrames
using Distributions
using Turing, MCMCChains
using StatsPlots

"""
    load_therapeutic_touch_data(;plot = false)

Loads the TherapeuticTouchData into a DataFrame. If plot is true, it also
displays the histogram of the observed accuracy of practitioners (figure 9.9).
"""
function load_therapeutic_touch_data(;plot = false)
    df = CSV.read("./data/TherapeuticTouchData.csv")
    if plot
        freqtable = combine(groupby(df, :s), :y => mean) # compute the accuracy of each practitioner
        fig = histogram(freqtable[:y_mean], nbins = 8, xlabel = "Observed accuracy", ylabel = "# Practitioners", legend = false)
        display(fig)
    end
    f = x -> parse(Int,x[2:end]) # parse the string col
    df[:subject_n] = f.(df[:s])
    return df
end

"""
    tt_sample_basic_model(;iterations = 1000, nchains = 3)

Estimates the accuracy separately for each subject of the TherapeuticTouchData with ``\\beta(1,1)`` prior. Returns an MCMCChain.
"""
function tt_sample_basic_model(;iterations = 1000, nchains = 3)
    df = load_therapeutic_touch_data()
    # model definition
    @model tt_unpooled(obs, subj) = begin
        # Assumptions: one θ for each subject
        n_subjects = length(unique(subj))
        θ = Vector{Real}(undef, n_subjects)
        for i=1:n_subjects
            θ[i] ~ Beta(1,1) # unifor priors
        end
        # Observations
        for j=1:length(obs)
            obs[j] ~ Bernoulli(θ[subj[j]])
        end
    end
    obs, subj = (df[:y], df[:subject_n])
    ϵ = 0.05
    τ = 10
    chains = mapreduce(c -> sample(tt_unpooled(obs, subj),HMC(ϵ, τ), iterations), chainscat, 1:nchains)
    return chains
end

"""
    tt_sample_hierarchical_model(;iterations = 1000, nchains = 3)

Estimates the accuracy for each subject of the TherapeuticTouchData in a hierarchical model.
This model is shown in the book in figure 9.13.
Returns an MCMCChain.
"""
function tt_sample_hierarchical_model(;iterations = 1000, nchains = 3)
    df = load_therapeutic_touch_data()
    # model definition
    @model therapeutic_touch(obs, subj) = begin
        # Assumptions: one θ for each subject
        n_subjects = length(unique(subj))
        ω ~ Beta(1,1)
        κ ~ Gamma(0.01, 100) # notice the different parametrization wrt PYMC3
        a = ω * κ + 1
        b = (1-ω) * κ + 1
        θ = Vector{Real}(undef, n_subjects)
        for i=1:n_subjects
            θ[i] ~ Beta(a,b)
        end
        # Observations
        for j=1:length(obs)
            obs[j] ~ Bernoulli(θ[subj[j]])
        end
    end
    obs, subj = (df[:y], df[:subject_n])
    ϵ = 0.05
    τ = 10
    chains = mapreduce(c -> sample(therapeutic_touch(obs, subj),HMC(ϵ, τ), iterations), chainscat, 1:nchains)
    return chains
end

"""
    plot_tt_shrinkage_effect(basic, hierarchical)

Demonstrates the effect of shrinkage obtained using a hierarchical model instead of the basic model for the TherapeuticTouchData.
Usage:
```julia-repl
julia> basic = tt_sample_basic_model(nchains=1);
julia> hierarchical = tt_sample_hierarchical_model(nchains=1);
julia> plot_shrinkage_effect(basic, hierarchical)
```
"""
function plot_tt_shrinkage_effect(basic, hierarchical)
    basic_means = mean(Array(basic),dims=1)
    hierarchical_means = reshape(mean(Array(hierarchical),dims=1)[1:end-2], (1,28));
    labels = reshape(["θ[$(i)]" for i=1:28],(1,28))
    p = plot(["basic", "hierarchical"],
        [basic_means; hierarchical_means],
        markershape = :circle,
        legend = :outertopright,
        labels = labels,
        size = (700, 500),
        xlabel = "Model",
        ylabel = "Estimated accuracy"
        )
    display(p)
end


text:
Base.String(p::NeuroPath) = joinpath(String(dirname(p)), basename(p))

# repath
repath(p::StudyPath) = joinpath(dirname(p), basename(p))
repath(p::Study) = basename(p)

repath(p::DataPath) = joinpath(repath(p), basename(p))
repath(p::Data) = basename(p)

repath(p::SubjectPath) = joinpath(repath(p), basename(p))
repath(p::Subject) = basename(p)

repath(p::SessionPath) = joinpath(repath(p), basename(p))
repath(p::Session) = basename(p)

repath(p::ModalityPath) = joinpath(repath(p), basename(p))
repath(p::Modality) = basename(p)

repath(p::FilePath) = joinpath(repath(p), basename(p))
repath(p::File) = basename(p)

Base.basename(p::StudyPath) = getfield(p, :study)
Base.basename(p::DataPath) = getfield(p, :basename)
Base.basename(p::SubjectPath) = getfield(p, :subject)
Base.basename(p::SessionPath) = getfield(p, :session)
Base.basename(p::ModalityPath) = getfield(p, :modality)
Base.basename(p::FilePath) = getfield(p, :basename)
Base.basename(p::DerivativePath) = getfield(p, :pipeline)

Base.dirname(p::NeuroPath) = getfield(p, :dirname)
Base.dirname(p::NeuroPath{Nothing}) = ""
Base.dirname(p::PathVector) = getfield(p, :dirname)

Base.splitpath(p::StudyPath) = (dirname(p), Study(basename(p)))
Base.splitpath(p::Study) = (p,)

Base.splitpath(p::DataPath) = (splitpath(dirname(p))..., Data(basename(p)))
Base.splitpath(p::Data) = (p,)

Base.splitpath(p::SubjectPath) = (splitpath(dirname(p))..., Subject(basename(p)))
Base.splitpath(p::Subject) = (p,)

Base.splitpath(p::SessionPath) = (splitpath(dirname(p))..., basename(p))
Base.splitpath(p::Session) = (p,)

Base.splitpath(p::ModalityPath) = (splitpath(dirname(p))..., Modality(basename(p)))
Base.splitpath(p::Modality) = (p,)

Base.joinpath(x::AbstractString, y::Study, z...) = joinpath(joinpath(x, y), z...)
Base.joinpath(x::NeuroPath, y, z...) = joinpath(joinpath(x, y), z...)
Base.joinpath(x::NeuroPath, y::AbstractVector{NeuroPath}) = PathVector(x, y)
Base.joinpath(x::AbstractString, y::Study) = StudyPath(x, basename(y))
Base.joinpath(x::StudyPath, y::Data) = DataPath(x, basename(y))
Base.joinpath(x::StudyPath, y::Subject) = SubjectPath(x, basename(y))
Base.joinpath(x::Union{StudyPath,DataPath}, y::Subject) = SubjectPath(x, basename(y))
Base.joinpath(x::DataPath, y::Derivative) = DerivativePath(x, basename(y))
Base.joinpath(x::SubjectPath, y::Session) = SessionPath(x, basename(y))
Base.joinpath(x::Union{SubjectPath,SessionPath}, y::Modality) = ModalityPath(x, basename(y))
Base.joinpath(x::NeuroPath, y::File) = FilePath(x, basename(y))
@inline function Base.joinpath(x::ModalityPath, y::File)
    me = basename(x)
    if me === "anat"
        return _anat_file(x, y)
    elseif me === "dwi"
        return _dwi_file(x, y)
    elseif s === "eeg"
        return _eeg_file(x, y)
    elseif s === "fmap"
        return _fmap_file(x, y)
    elseif s === "func"
        return _func_file(x, y)
    elseif s === "ieeg"
        return _ieeg_file(x, y)
    elseif s === "meg"
        return _meg_file(x, y)
    elseif s === "perf"
        return _perf_file(x, y)
    elseif s === "pet"
        return _pet_file(x, y)
    else # s === :beh
        return _beh_file(x, y)
    end
end

Base.splitext(p::NeuroPath) = splitext(basename(p))

Base.isdir(p::FilePath) = false
Base.isdir(p::NeuroPath) = isdir(repath(p))

Base.isfile(p::FilePath) = isfile(repath(p))
Base.isfile(p::NeuroPath) = false

Base.ispath(p::NeuroPath) = isdir(p)
Base.ispath(p::FilePath) = isfile(p)

Base.mkdir(p::NeuroPath) = mkdir(repath(p))

Base.touch(p::FilePath) = touch(repath(p))
Base.cp(x::NeuroPath, y::NeuroPath) = cp(repath(x), repath(y))
Base.mv(x::NeuroPath, y::NeuroPath) = mv(repath(x), repath(y))
Base.rm(p::NeuroPath) = rm((p))
Base.open(p::FilePath) = open(repath(p))
Base.open(p::FilePath, mode) = open(repath(p), mode)
Base.readdir(p::NeuroPath) = readdir(repath(p))
Base.read(p::FilePath) = read(repath(p))
Base.write(p::FilePath, data) = write(repath(p), data)

is_derivative(::Nothing) = false
is_derivative(p::StudyPath) = false
is_derivative(p::ModalityPath) = basename(p) === "derivatives"
is_derivative(p::NeuroPath) = is_derivative(dirname(p))

"""
    session(p::NeuroPath)
    session(s::AbstractString)

Return the session from a `NeuroPath` or compose a new instance of `SessionPath`.
"""
session(p::FilePath) = session(dirname(p))
session(p::ModalityPath) = session(dirname(p))
session(p::SessionPath) = p
session(p::AbstractString) = SessionPath(p)

"""
    subject(p::NeuroPath)
    subject(s::AbstractString)

Return the subject from a `NeuroPath` or compose a new instance of `SubjectPath`.
"""
subject(p::FilePath) = subject(dirname(p))
subject(p::ModalityPath) = subject(dirname(p))
subject(p::SessionPath) = subject(dirname(p))
subject(p::SubjectPath) = p
subject(p::AbstractString) = SubjectPath(p)

"""
    modality(p::NeuroPath)
    modality(s::AbstractString)

Return the modality from a `NeuroPath` or compose a new instance of `ModalityPath`.
"""
modality(p::ModalityPath) = p
modality(p::FilePath) = modality(dirname(p))
modality(p::AbstractString) = ModalityPath(p)

"""
    study(p::NeuroPath)
    study(s::AbstractString)

Return the study from a `NeuroPath` or compose a new instance of `StudyPath`.
"""
study(p::StudyPath) = p
study(p::NeuroPath) = studydir(dirname(p))
study(p::AbstractString) = StudyPath(p)

"""
    derivative(p::NeuroPath)
    derivative(s::AbstractString)

Return the derivative subpath from a `NeuroPath` or compose a new instance of `DerivativePath`.
"""
derivative(p::NeuroPath) = derivative(dirname(p))
derivative(p::DerivativePath) = p
derivative(p::AbstractString) = DerivativePath(p)

file(p::AbstractString) = FilePath(p)

## Derivative
function Base.relpath(p::DataPath{StudyPath{D}}, r::Derivative) where {D}
    DerivativePath(DataPath(dirname(p), "derivatives"), basename(r))
end
function Base.relpath(p::SubjectPath{StudyPath{D}}, r::Derivative) where {D}
    SubjectPath(DerivativePath(DataPath(dirname(p), "derivatives"), basename(r)), basename(p))
end
Base.relpath(p::SessionPath, r::Derivative) = SessionPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::ModalityPath, r::Derivative) = ModalityPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::FilePath, r::Derivative) = ModalityPath(relpath(dirname(p), r), basename(p))

## Subject
Base.relpath(p::SubjectPath, r::Subject) = SubjectPath(dirname(p), basename(r))
Base.relpath(p::SessionPath, r::Subject) = SessionPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::ModalityPath, r::Subject) = ModalityPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::FilePath, r::Subject) = FilePath(relpath(dirname(p), r), basename(p))

## Session
Base.relpath(p::SessionPath, r::Session) = SessionPath(dirname(p), basename(r))
Base.relpath(p::ModalityPath, r::Session) = ModalityPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::FilePath, r::Session) = FilePath(relpath(dirname(p), r), basename(p))

## Modality
Base.relpath(p::ModalityPath, r::Modality) = ModalityPath(relpath(dirname(p), r), basename(p))
Base.relpath(p::FilePath, r::Modality) = ModalityPath(relpath(dirname(p), r), basename(p))

Base.mkpath(p::NeuroPath; mode::Integer = 0o777) = _mkpath(p, mode)
function _mkpath(x, mode)
    if !isdir(x)
        mkdir(x, mode)
    end
    return x
end
function _mkpath(x::NeuroPath, mode)
    pd = _mkpath(dirname(x), mode)
    pp = joinpath(p, basename(x))
    if !isdir(pp)
        mkdir(pp, mode)
    end
    return pp
end


text:
module JsonIO
using JSON

export loadDataset, saveDataset

function loadDataset(jsonFile::String)
    # Load a json dset
    f = open(jsonFile)
    dset = JSON.parse(f)
    close(f)
    return dset
end

function saveDataset(dset::Dict, jsonFile::String)
    # Save a Dict as json file
    f = open(jsonFile, "w")
    write(f, JSON.json(dset))
    close(f)
end

end#module


text:
using TensorDecompositions

function tensordecomp2(Λ::AbstractArray{T,3}, p::Int=5;
    scaleout::Bool=true) where T
    Λ=convert(Array{Float64}, Λ)
    size(Λ)
    𝒯=TensorDecompositions.candecomp(Λ, p,  (randn(size(Λ,1), k), randn(size(Λ,2), k), randn(size(Λ,3), k)), compute_error=true, method=:ALS);
    # 𝒯=TensorToolbox.cp_als(Λ,p)
    𝕏=TensorDecompositions.compose(𝒯)
    A=mean(𝕏[:,:,i] for i=1:size(𝕏,3))
    if scaleout
    A ./=maximum(abs.(A))
    end
    return A
end

#=
u = randn(10); v = randn(20); w = randn(30)
T = cat(map(x -> x * u * v', w)..., dims=3) + 0.2 * randn(10, 20, 30)
k=2
F = candecomp(T, k, (randn(10, k), randn(20, k), randn(30, k)), compute_error=true, method=:ALS);

# https://juliapackages.com/p/ntfk

import NTFk
import TensorDecompositions

csize = (2, 3, 4)
tsize = (5, 10, 15)
tucker_orig = NTFk.rand_tucker(csize, tsize; factors_nonneg=true, core_nonneg=true)

T_orig = TensorDecompositions.compose(tucker_orig)
T_orig .*= 1000

sizes = [csize, (1,3,4), (3,3,4), (2,2,4), (2,4,4), (2,3,3), (2,3,5)]
tucker_estimated, csize_estimated = NTFk.analysis(T_orig, sizes, 3; eigmethod=[false,false,false], progressbar=false, tol=1e-16, max_iter=100000, lambda=0.);
=#

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
# Simple recursive descent parser

module Parser

using Base.Iterators
using ..Lists
using ..Keywords
import Base: convert

export parse, SExpression

"""
Union type of all types allowed in primitive s-expressions. The output of `parse` is
guaranteed to be an `SExpression`.
"""
const SExpression =
    Union{BigInt,
          Float64,
          Rational{BigInt},
          Complex{BigInt},
          Complex{Float64},
          Complex{Rational{BigInt}},
          Nothing,
          Bool,
          String,
          Symbol,
          List,
          Keyword}

convert(::Type{SExpression}, x::Integer) = BigInt(x)
convert(::Type{SExpression}, x::AbstractFloat) = Float64(x)
convert(::Type{SExpression}, x::Rational) = Rational{BigInt}(x)
convert(::Type{SExpression}, x::Complex{<:Integer}) = Complex{BigInt}(x)
convert(::Type{SExpression}, x::Complex{<:AbstractFloat}) = Complex{Float64}(x)
convert(::Type{SExpression}, x::Complex{<:Rational}) = Complex{Rational{BigInt}}(x)
convert(::Type{SExpression}, x::AbstractString) = String(x)
convert(::Type{SExpression}, xs::Tuple) = List(convert.(SExpression, xs)...)
convert(::Type{SExpression}, xs::AbstractVector) = List(convert.(SExpression, xs)...)
SExpression(x) = convert(SExpression, x)

include("numeric.jl")
include("util.jl")

const _DELIMITERS = collect("()[]{}\",'`;")

const _NEWLINE = (
        '\u000D', '\u000A', '\u0085', '\u000B', '\u000C', '\u2028', '\u2029')
isnewline(c) = c ∈ _NEWLINE

const _CLOSE = Dict('(' => ')', '[' => ']', '{' => '}', '"' => '"')
closer(c) = _CLOSE[c]

const _READER_MACROS = Dict(
    '\'' => :quote,
    '`' => :quasiquote,
    ',' => :unquote)

const _ESCAPED = Dict('a' => '\x07', 'b' => '\x08', 't' => '\x09',
                      'n' => '\x0a', 'v' => '\x0b', 'f' => '\x0c',
                      'r' => '\x0d', 'e' => '\x1b', '"' => '"',
                      '\'' => '\'', '\\' => '\\')

isdelimiter(c) = isspace(c) || c in _DELIMITERS

macro syntaxcheck(ex)
    quote
        if !$(esc(ex))
            # TODO improve the error messages
            error("invalid syntax")
        end
    end
end

"""
    skipws(io::IO)

Skip all whitespace characters at beginning of the given `io` object.
"""
function skipws(io::IO)
    while !eof(io)
        c = read(io, Char)
        if !isspace(c)
            skip(io, -readsize(io, c))
            return
        end
    end
end

"""
    readsymbol(io::IO)

Read up until the next delimiter.
"""
function readsymbol(io::IO)
    buf = IOBuffer()
    while !eof(io)
        c = read(io, Char)
        if isdelimiter(c)
            skip(io, -readsize(io, c))
            return String(take!(buf))
        end
        write(buf, c)
    end
    String(take!(buf))
end

"""
A dot as in (a . b).
"""
struct Dot end

"""
Build an improper list (a b . c).
"""
function improperlist(prefix, terminal)
    if isempty(prefix)
        terminal
    else
        Cons(prefix[1], improperlist(prefix[2:end], terminal))
    end
end

"""
Check that there are no . in objs, and return the list if so. Otherwise process the dots as
they should be.
"""
function listify(objs) :: List
    objs_ = collect(objs)
    dots = count(isequal(Dot()), objs_)
    if dots == 0
        objs_
    elseif dots == 1
        if length(objs_) ≤ 2 || objs_[end - 1] !== Dot()
            error("invalid use of .")
        else
            improperlist(objs_[1:end-2], objs_[end])
        end
    elseif dots == 2
        if length(objs_) == 5 && objs_[2] === objs_[4] === Dot()
            [objs_[3], objs_[1], objs_[5]]
        else
            error("invalid use of .")
        end
    else
        error("invalid use of .")
    end
end

"""
    nextobject(io::IO)

Obtain the next Racket object from the given `io` object, wrapped in `Some`, or `nothing` if
there is no more input left.
"""
function nextobject(io::IO) :: Union{Some{<:Union{Dot, SExpression}}, Nothing}
    skipws(io)
    if eof(io)
        return nothing
    end

    c = peek(io, Char)
    if c == ';'
        skipcomment(io)
        nextobject(io)
    elseif c in "([{"
        read(io, Char)
        cl = closer(c)
        objects = readobjectsuntil(io, cl)
        Some(listify(objects))
    elseif c in ")]}"
        error("mismatched superfluous $c")
    elseif c in keys(_READER_MACROS)
        read(io, Char)
        nextobj = nextobject(io)
        if nextobj === nothing
            error("lonely reader macro $c")
        else
            Some(List(_READER_MACROS[c], something(nextobj)))
        end
    elseif c == '"'
        read(io, Char)
        Some(parsestring(io))
    elseif c == '#'
        read(io, Char)
        Some(readhash(io))
    else
        str = readsymbol(io)
        if str == "."
            Some(Dot())
        else
            asnumber = tryparse(Number, str)
            Some(something(asnumber, Symbol(str)))
        end
    end
end

"""
    readobjectsuntil(io::IO, cl::Char)

Read as many Racket objects are are available, then read the closer character.
Return a vector of the read objects.
"""
function readobjectsuntil(io::IO, cl::Char)
    objects = []
    while !eof(io)
        skipws(io)
        c = peek(io, Char)
        if c == cl
            read(io, Char)
            return objects
        else
            obj = nextobject(io)
            if obj === nothing
                break
            else
                push!(objects, something(obj))
            end
        end
    end
    error("unexpected end of file before expected closing $cl character")
end

"""
    readhash(io::IO)

Read the `io` as if a `#` symbol was just seen, and return the read object.
"""
function readhash(io::IO)
    x = read(io, Char)
    if x == 't'
        @syntaxcheck eof(io) || isdelimiter(peek(io, Char)) || begin
            read(io, Char) == 'r'
            read(io, Char) == 'u'
            read(io, Char) == 'e'
        end
        true
    elseif x == 'T'
        @syntaxcheck eof(io) || isdelimiter(peek(io, Char))
        true
    elseif x == 'f'
        @syntaxcheck eof(io) || isdelimiter(peek(io, Char)) || begin
            read(io, Char) == 'a'
            read(io, Char) == 'l'
            read(io, Char) == 's'
            read(io, Char) == 'e'
        end
        false
    elseif x == 'F'
        @syntaxcheck eof(io) || isdelimiter(peek(io, Char))
        false
    elseif x == ':'
        @syntaxcheck !eof(io)
        Keyword(readsymbol(io))
    else
        error("Unsupported syntax: #$x")
    end
end

"""
    skipcomment(io::IO)

Skip all characters until the end of line or end of file.
"""
function skipcomment(io::IO)
    while !eof(io) && !isnewline(read(io, Char)) end
end

"""
    parsestring(io::IO)

Read a string from the given `io` stream, using Racket parsing rules.
"""
function parsestring(io::IO)
    buf = IOBuffer()
    while !eof(io)
        c = read(io, Char)
        if c == '\\'
            x = read(io, Char)
            if x in '0':'7'
                val = x - '0'
                for _ in 1:2
                    if peek(io, Char) in '0':'7'
                        val *= 8
                        val += read(io, Char) - '0'
                    end
                end
                Char(val)
            elseif x in "xuU"
                # TODO
                error("escape sequence \\$x not yet supported")
            elseif x == '\r'
                if !eof(io) && peek(io, Char) == '\n'
                    read(io, Char)
                end
            elseif x == '\n'  # no-op
            elseif x in keys(_ESCAPED)
                write(buf, _ESCAPED[x])
            else
                error("invalid escape sequence \\$x")
            end
        elseif c == '"'
            return String(take!(buf))
        else
            write(buf, c)
        end
    end
    error("unexpected EOF while waiting for terminating \"")
end

"""
    parse(s::AbstractString)

Read the given string `s` as a single s-expression. The alternative vocabulary
`Base.parse(SExpression, s)` is also supported and may be preferred.
"""
function parse(s::AbstractString)
    buf = IOBuffer(s)
    obj = nextobject(buf)
    obj2 = nextobject(buf)
    if obj === nothing
        error("no object to read")
    elseif obj2 === nothing
        something(obj) :: SExpression
    else
        error("extra content after end of expression")
    end
end

"""
    parse(io::IO)

Read a single object from the given `io` stream.
"""
parse(io::IO) = let x = nextobject(io)
    if x === nothing
        error("no object to read")
    else
        something(x) :: SExpression
    end
end

# allow `parse(SExpression, o)`, which reads nicer than `SExpression.parse(o)`
Base.parse(::Type{SExpression}, o) = parse(o) :: SExpression

"""
    parsefile(filename::AbstractString)

Parse a file into a single list.
"""
parsefile(filename::AbstractString) = open(parseall, filename)

"""
    parseall(io::IO)
    parseall(s::AbstractString)

Parse all objects from the given stream or string into a single list.
"""
function parseall(io::IO)
    result = []
    while (obj = nextobject(io); obj !== nothing)
        push!(result, something(obj) :: SExpression)
    end
    convert(List, result)
end
parseall(s::AbstractString) = parseall(IOBuffer(s))

end


text:
"""
    read_opf(file, attr_type = Dict, mtg_type = MutableNodeMTG)

Read an OPF file, and returns an MTG.

# Arguments

- `file::String`: The path to the opf file.
- `attr_type::DataType = Dict`: the type used to hold the attribute values for each node.
- `mtg_type = MutableNodeMTG`: the type used to hold the mtg encoding for each node (*i.e.*
link, symbol, index, scale). See details section below.

# Details

`attr_type` should be:

- `NamedTuple` if you don't plan to modify the attributes of the mtg, *e.g.* to use them for
plotting or computing statistics...
- `MutableNamedTuple` if you plan to modify the attributes values but not adding new attributes
very often, *e.g.* recompute an attribute value...
- `Dict` or similar (*e.g.* `OrderedDict`) if you plan to heavily modify the attributes, *e.g.*
adding/removing attributes a lot

The `MultiScaleTreeGraph` package provides two types for `mtg_type`, one immutable
([`MultiScaleTreeGraph.NodeMTG`](@ref)), and one mutable
([`MultiScaleTreeGraph.MutableNodeMTG`](@ref)). If you're planning on modifying the mtg
encoding of some of your nodes, you should use [`MultiScaleTreeGraph.MutableNodeMTG`](@ref),
and if you don't want to modify anything, use [`MultiScaleTreeGraph.NodeMTG`](@ref) instead
as it should be faster.

# Note

See the documentation of the MTG format from the MTG package documentation for further details,
*e.g.* [The MTG concept](https://vezy.github.io/MultiScaleTreeGraph.jl/stable/the_mtg/mtg_concept/).

# Returns

The MTG root node.

# Examples

```julia
using PlantGeom
file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
# file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","coffee.opf")
opf = read_opf(file)
```
"""
function read_opf(
    file,
    attr_type = Dict,
    mtg_type = MultiScaleTreeGraph.MutableNodeMTG,
    geom_type = Meshes.SimpleMesh
)

    doc = readxml(file)
    xroot = root(doc)
    line = [2]

    if xroot.name != "opf"
        error("The file is not an OPF")
    end

    if xroot["version"] != "2.0"
        error("Cannot read OPF files version other than 2.0")
    end

    editable = parse(Bool, xroot["editable"])

    opf_attr = Dict{Symbol,Any}()
    # node = elements(xroot)[5]
    for node in eachelement(xroot)
        if node.name == "meshBDD"
            push!(opf_attr, :meshBDD => parse_meshBDD!(node))
        end

        if node.name == "materialBDD"
            push!(
                opf_attr,
                :materialBDD => parse_materialBDD!(node)
            )
        end

        if node.name == "shapeBDD"
            shapeBDD = parse_opf_elements!(
                node,
                [String, Int64, Int64]
            )

            # Increment index by 1 because Julia is 1-based indexing (OPF is 0-based):
            for (key, value) in shapeBDD
                value["materialIndex"] += 1
                value["meshIndex"] += 1
            end

            push!(opf_attr, :shapeBDD => shapeBDD)
        end

        if node.name == "attributeBDD"
            push!(opf_attr, :attributeBDD => parse_opf_attributeBDD!(node))
        end

        if node.name == "topology"
            ref_meshes = parse_ref_meshes(opf_attr)

            global mtg = parse_opf_topology!(
                node,
                nothing,
                get_attr_type(opf_attr[:attributeBDD]),
                attr_type,
                mtg_type,
                ref_meshes
            )

            append!(
                mtg,
                MultiScaleTreeGraph.node_attributes(
                    attr_type,
                    Dict(:ref_meshes => ref_meshes)
                )
            )

            return mtg
        end
    end
end

"""
Parse an array of values from the OPF into a Julia array (Arrays in OPFs
are not following XML recommendations)
"""
function parse_opf_array(elem, type = Float64)
    if type == String
        strip(elem)
    else
        parsed = map(x -> x == "NA" ? nothing : parse(type, x), split(elem))
        if length(parsed) == 1
            return parsed[1]
        else
            return parsed
        end
    end
end

"""
Parse the meshBDD using [`parse_opf_array`](@ref)
"""
function parse_meshBDD!(node)
    # MeshBDD:
    meshes = Dict{Int,Dict{String,Any}}()

    for m in eachelement(node)
        m.name != "mesh" ? @warn("Unknown node element in meshBDD: $(m.name)") : nothing
        mesh = Dict{String,Any}()
        mesh["name"] = m["name"]
        mesh["enableScale"] = parse(Bool, m["enableScale"])
        for i in eachelement(m)
            if i.name == "faces"
                push!(mesh, i.name => parse_opf_array(i.content, Int) .+ 1)
                # NB: adding 1 to the faces because the opf is 0-based but Julia is 1-based
            else
                push!(mesh, i.name => parse_opf_array(i.content))
            end
        end
        push!(meshes, parse(Int, m["Id"]) + 1 => mesh)
    end

    return meshes
end

"""
Parse the materialBDD using [`parse_opf_elements!`](@ref)
"""
function parse_materialBDD!(node)
    metBDDraw = parse_opf_elements!(
        node,
        [Float64, Float64, Float64, Float64, Float64]
    )
    metBDD = Dict{Int,Phong}()
    for (key, value) in metBDDraw
        # key = 1
        # value = metBDDraw[key]
        push!(metBDD, key => materialBDD_to_material(value))
    end

    return metBDD
end

"""
Generic parser for OPF elements.

# Arguments

- `opf::OrderedDict`: the opf Dict (using [XMLDict.xml_dict])
- `elem_types::Array`: the target types of the element (e.g. "[String, Int64]")

# Details

`elem_types` should be of the same length as the number of elements found in each
item of the subchild.
elem_types = [Float64, Float64, Float64, Float64, Float64, Float64]
"""
function parse_opf_elements!(node, elem_types)
    elem_dict = Dict()
    for m in eachelement(node)
        elems_ = Dict()
        for (i, el) in enumerate(elements(m))
            content = parse_opf_array(el.content, elem_types[i])
            if length(content) > 0
                push!(elems_, el.name => content)
            end
        end
        push!(elem_dict, parse(Int, m["Id"]) + 1 => elems_)
    end
    return elem_dict
end

"""
 Parse the opf attributes as a Dict.
"""
function parse_opf_attributeBDD!(node)
    elem_dict = Dict()
    for m in eachelement(node)
        push!(elem_dict, m["name"] => m["class"])
    end
    return elem_dict
end

"""
Get the attributes types in Julia `DataType`.
"""
function get_attr_type(attr)
    attr_Type = Dict{String,DataType}()
    for i in keys(attr)
        if attr[i] in ["Object", "String", "Color", "Image"]
            push!(attr_Type, i => String)
        elseif attr[i] == "Integer"
            push!(attr_Type, i => Int32)
        elseif attr[i] in ["Double", "Metre", "Centimetre", "Millimetre", "10E-5 Metre", "Metre_100"]
            push!(attr_Type, i => Float32)
        elseif attr[i] == "Boolean"
            push!(attr_Type, i => Bool)
        else
            error("Attribute type `$(attr[i])` not recognised in attributeBDD.")
        end
    end
    return attr_Type
end

"""
Parse the geometry element of the OPF.

# Note
The transformation matrix is 3*4.
elem = elem.content
"""
function parse_geometry(elem)
    geom = Dict{Symbol,Union{Int,Float64,SMatrix{3,4}}}()
    for i in eachelement(elem)
        if i.name == "shapeIndex"
            push!(geom, :shapeIndex => parse(Int, i.content) + 1)
        elseif i.name == "mat"
            push!(geom, :mat => SMatrix{3,4}(reshape(parse_opf_array(i.content), 4, 3)'))
        elseif i.name == "dUp"
            push!(geom, :dUp => parse(Float64, i.content))
        elseif i.name == "dDwn"
            push!(geom, :dDwn => parse(Float64, i.content))
        end
    end
    return geom
end

"""
Parser for OPF topology.

# Note

The transformation matrices in `geometry` are 3*4.
parse_opf_topology!(elem, node_i, features)
node = elem
mtg = node_i
features = get_attr_type(opf_attr[:attributeBDD])

# Debugging:
mtg = nothing

node = elem
mtg = node_i
parse_opf_topology!(
                node,
                nothing,
                get_attr_type(opf_attr[:attributeBDD]),
                attr_type,
                mtg_type,
                ref_meshes
            )
"""
function parse_opf_topology!(node, mtg, features, attr_type, mtg_type, ref_meshes)

    link = "/" # default, for "topology" and "decomp"

    if node.name == "branch"
        link = "+"
    elseif node.name == "follow"
        link = "<"
    end

    id = parse(Int, node["id"])

    MTG = mtg_type(
        link,
        node["class"],
        id,
        parse(Int, node["scale"])
    )

    if mtg !== nothing
        node_i = Node(
            id,
            mtg,
            MTG,
            MultiScaleTreeGraph.init_empty_attr(attr_type)
        )
    else
        # First node:
        node_i = Node(id, MTG, MultiScaleTreeGraph.init_empty_attr(attr_type))
    end

    # node_i.children
    attrs = Dict{Symbol,Any}()

    # Array with the last node computed. This is used for the "follow" nodes, because they always
    # refer to the last computed node...
    last_node = Node[node_i]

    # Handle the children, can be attributes of children nodes:
    # elem = elements(node)[1]
    for elem in eachelement(node)
        # If an element is an attribute, add it to the attributes of the Node:
        if elem.name in keys(features)
            push!(attrs, Symbol(elem.name) => parse_opf_array(elem.content, features[elem.name]))
        elseif elem.name == "geometry"

            geom = parse_geometry(elem)

            # Parse the geometry (transformation, reference mesh + index and dUp and dDwn):
            if haskey(geom, :shapeIndex)
                # Rotation + Scaling. No need to decouple them here, but in case we need to
                # in the future, see: https://stackoverflow.com/a/29618569/6947799
                # See also this for decomposition: https://colab.research.google.com/drive/1ImBB-N6P9zlNMCBH9evHD6tjk0dzvy1_

                #! OK what I could do is use my own transformation function that adds w (=1)
                #! to the Point3 when transforming it with the 4x4 matrix?

                transformation = Translation(geom[:mat][1:3, 4]) ∘ LinearMap(geom[:mat][1:3, 1:3])
                # NB: We read an homogeneous transformation matrix from the OPF, but we work
                # with cartesian coordinates in PlantGeom by design. So we deconstruct our
                # homogeneous matrix into the two corresponding rotation and translation
                # matrices, and create a transformation from them using transformation
                # composition from CoordinateTransformations.jl. Note that our homogeneous
                # matrix is rotate first, then translate (hence the order of transformation)

                push!(
                    attrs,
                    Symbol(elem.name) => geometry(
                        ref_meshes.meshes[geom[:shapeIndex]],
                        geom[:shapeIndex],
                        transformation,
                        geom[:dUp],
                        geom[:dDwn],
                        nothing
                    )
                )
            end
        elseif elem.name in ["decomp", "branch"]
            last_node[1] = parse_opf_topology!(elem, node_i, features, attr_type, mtg_type, ref_meshes)
        elseif elem.name == "follow"
            # We use last_node here instead of node_i because a following node refers to the last computed node
            last_node[1] = parse_opf_topology!(elem, last_node[1], features, attr_type, mtg_type, ref_meshes)
        else
            error("Attribute $(elem.name) not found in attributeBDD (or badly written?)")
        end
    end

    node_i.attributes = MultiScaleTreeGraph.node_attributes(attr_type, attrs)

    return node_i
end


text:
using DrWatson
@quickactivate "GRAPETests"
import Downloads

DOWNLOADS = Dict(
    "https://raw.githubusercontent.com/JuliaQuantumControl/GRAPE.jl/data-dump/TLS/opt_result_LBFGSB.jld2" =>
        joinpath(datadir("TLS"), "opt_result_LBFGSB.jld2"),
    "https://raw.githubusercontent.com/JuliaQuantumControl/GRAPE.jl/data-dump/TLS/opt_result_OptimLBFGS.jld2" =>
        joinpath(datadir("TLS"), "opt_result_OptimLBFGS.jld2"),
    "https://raw.githubusercontent.com/JuliaQuantumControl/GRAPE.jl/data-dump/PE_OCT.jld2" =>
        joinpath(datadir(), "PE_OCT.jld2"),
    "https://raw.githubusercontent.com/JuliaQuantumControl/GRAPE.jl/data-dump/PE_OCT_direct.jld2" =>
        joinpath(datadir(), "PE_OCT_direct.jld2"),
)

function download_dump(url, destination; force=false, verbose=true)
    if !isfile(destination) || force
        verbose && (@info "Downloading $url => $destination")
        mkpath(dirname(destination))
        Downloads.download(url, destination)
    else
        verbose && (@info "$destination OK")
    end
end

@info "Download Dumps"
for (url, destination) in DOWNLOADS
    download_dump(url, destination)
end