text:
export Combiner

struct Combiner <: TensorStorage{Number}
end

data(::Combiner) = error("Combiner storage has no data")

Base.eltype(::Type{<:Combiner}) = Nothing
Base.eltype(::StoreT) where {StoreT<:Combiner} = eltype(StoreT)

Base.promote_rule(::Type{<:Combiner},StorageT::Type{<:Dense}) = StorageT

#
# CombinerTensor (Tensor using Combiner storage)
#

const CombinerTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Combiner}

function contraction_output(::TensorT1,
                            ::TensorT2,
                            indsR::IndsR) where {TensorT1<:CombinerTensor,
                                                 TensorT2<:DenseTensor,
                                                 IndsR}
  TensorR = contraction_output_type(TensorT1,TensorT2,IndsR)
  return similar(TensorR,indsR)
end

function contraction_output(T1::TensorT1,
                            T2::TensorT2,
                            indsR) where {TensorT1<:DenseTensor,
                                          TensorT2<:CombinerTensor}
  return contraction_output(T2,T1,indsR)
end

function contract!!(R::Tensor{<:Number,NR},
                    labelsR::NTuple{NR},
                    T1::CombinerTensor{<:Number,N1},
                    labelsT1::NTuple{N1},
                    T2::Tensor{<:Number,N2},
                    labelsT2::NTuple{N2}) where {NR,N1,N2}
  if N1 ≤ 1
    #println("identity")
    return R
  elseif N1 + N2 == NR
    error("Cannot perform outer product involving a combiner")
  elseif count_common(labelsT1,labelsT2) == 1 && length(inds(T1)) == 2
    ci = commonindex(inds(T1), inds(T2))
    ui = uniqueindex(inds(T1), inds(T2))
    inds2        = [inds(T2)...]
    cpos1,cpos2  = intersect_positions(labelsT1,labelsT2)
    inds2[cpos2] = ui
    return Tensor(copy(store(T2)), IndexSet(inds2...))
  elseif count_common(labelsT1,labelsT2) == 1 && length(inds(T1)) != 2
    # This is the case of Index replacement or uncombining
    T2data      = data(store(T2))
    cpos1,cpos2 = intersect_positions(labelsT1,labelsT2)
    indsC = inds(T1)
    indsT = inds(T2)

    newlength = (length(indsC)-1) + (length(indsT)-1)
    newinds = Vector{Index}(undef,newlength)
    n = 1
    # Copy existing indices before one we are uncombining
    for i in 1:cpos2-1
      newinds[n] = indsT[i]
      n += 1
    end
    # Replace uncombined index with indices of combiner
    for j in 2:length(indsC)
      newinds[n] = indsC[j]
      n += 1
    end
    # Copy existing indices after one we are uncombining
    for i in cpos2+1:length(indsT)
      newinds[n] = indsT[i]
      n += 1
    end

    return Tensor(Dense(copy(T2data)), IndexSet(newinds...))
  elseif is_combiner(labelsT1,labelsT2)
    # This is the case of combining
    Alabels,Blabels = compute_contraction_labels(inds(T2),inds(T1))
    final_labels    = contract_labels(Blabels, Alabels)
    final_labels_n  = contract_labels(labelsT1,labelsT2)
    indsR = inds(R)
    if final_labels != final_labels_n
      perm  = getperm(final_labels_n, final_labels)
      indsR = permute(inds(R), perm)
      labelsR = permute(labelsR, perm)
    end
    cpos1,cposR = intersect_positions(labelsT1,labelsR)
    labels_comb = deleteat(labelsT1,cpos1)
    vlR = [labelsR...]
    for (ii, li) in enumerate(labels_comb)
      insert!(vlR, cposR+ii, li)
    end
    deleteat!(vlR, cposR)
    labels_perm = tuple(vlR...)
    perm = getperm(labels_perm,labelsT2)
    T2p = reshape(R,permute(inds(T2),perm))
    permutedims!(T2p,T2,perm)
    R = reshape(T2p,indsR)
  end
  return R
end

function contract!!(R::Tensor{<:Number,NR},
                    labelsR::NTuple{NR},
                    T1::Tensor{<:Number,N1},
                    labelsT1::NTuple{N1},
                    T2::CombinerTensor{<:Number,N2},
                    labelsT2::NTuple{N2}) where {NR,N1,N2}
  return contract!!(R,labelsR,T2,labelsT2,T1,labelsT1)
end


text:
using Molly

temp = 298u"K"
atoms, atoms_data, specific_inter_lists, general_inters, neighbor_finder, coords, box_size = readinputs(
            joinpath(dirname(pathof(Molly)), "..", "data", "5XER", "gmx_top_ff.top"),
            joinpath(dirname(pathof(Molly)), "..", "data", "5XER", "gmx_coords.gro"))
simulator = VelocityVerlet(dt=0.0002u"ps", coupling=AndersenThermostat(temp, 1.0u"ps"))

sys = System(
    atoms=atoms,
    atoms_data=atoms_data,
    general_inters=general_inters,
    specific_inter_lists=specific_inter_lists,
    coords=coords,
    velocities=[velocity(a.mass, temp) for a in atoms],
    box_size=box_size,
    neighbor_finder=neighbor_finder,
    loggers=Dict("temp" => TemperatureLogger(10),
                    "writer" => StructureWriter(10, "traj_5XER_1ps.pdb")),
)

simulate!(sys, simulator, 5_000)

ext: jl


text:
using Test, FinancialMonteCarlo, DifferentialEquations;
@show "Differential Equation Junctor"
Nsim = 10000;
Nstep = 30;
toll = 0.8;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
S0 = 100.0
K = 100.0;
D = 90.0;
r = 0.02
sigma = 0.2
T = 1.0;
d = 0.01;
u0 = 0.0;

p1 = 0.3;
lam = 5.0;
lamp = 30.0;
lamm = 20.0;

#Drift
dr_(u, p, t) = (r - d - sigma^2 / 2.0 - lam * (p1 / (lamp - 1) - (1 - p1) / (lamm + 1)))
#Diffusion
g_1(u, p, t) = sigma
#Time Window
tspan = (0.0, T)
rate(u, p, t) = lam * T
affect!(integrator) = (integrator.u = integrator.u + ((rand(mc.rng) < p1) ? FinancialMonteCarlo.quantile_exp(lamp, rand(mc.rng)) : -FinancialMonteCarlo.quantile_exp(lamm, rand(mc.rng))))
jump = ConstantRateJump(rate, affect!)

#Definition of the SDE
prob = SDEProblem(dr_, g_1, u0, tspan, rng = mc.rng)
jump_prob = JumpProblem(prob, Direct(), jump)
monte_prob = EnsembleProblem(jump_prob)
rfCurve = ZeroRate(r);
func(x) = S0 * exp(x);
model = FinancialMonteCarlo.MonteCarloDiffEqModel(monte_prob, func, Underlying(S0, d))

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)

@show FwdPrice = pricer(model, rfCurve, mc, FwdData);
@show EuPrice = pricer(model, rfCurve, mc, EUData);
@show AmPrice = pricer(model, rfCurve, mc, AMData);
@show BarrierPrice = pricer(model, rfCurve, mc, BarrierData);
@show AsianPrice1 = pricer(model, rfCurve, mc, AsianFloatingStrikeData);
@show AsianPrice2 = pricer(model, rfCurve, mc, AsianFixedStrikeData);

@test abs(FwdPrice - 98.8436678850961) < toll

Tneg = -T;
tspanNeg = (0.0, Tneg)
probNeg = SDEProblem(dr_, g_1, u0, tspanNeg, rng = mc.rng)
monte_probNeg = EnsembleProblem(probNeg)


text:
# Unit Tests
using Base.Test
using MLKernels

FloatingPointTypes = (Float32, Float64)
IntegerTypes = (Int32, UInt32, Int64, UInt64)
MOD = MLKernels

module MLKTest

    using MLKernels

    import Base.convert

    immutable TestKernel{T<:AbstractFloat} <: StandardKernel{T}
        a::MLKernels.HyperParameter{T}
        TestKernel(a::Variable{T}) = new(HyperParameter(a, unbounded(T)))
    end
    TestKernel{T<:AbstractFloat}(a::Argument{T}=1.0) = TestKernel{T}(Variable(a))
    MLKernels.phi(κ::TestKernel, x::Real, y::Real) = x*y
    MLKernels.phi(κ::TestKernel, x::Vector, y::Vector) = product(x)*product(y)

    immutable PairwiseTestKernel{T<:AbstractFloat} <: PairwiseKernel{T}
        a::MLKernels.HyperParameter{T}
        PairwiseTestKernel(a::Variable{T}) = new(HyperParameter(a, unbounded(T)))
    end
    eval(MLKernels.generate_outer_constructor(PairwiseTestKernel, (1,)))
    eval(MLKernels.generate_conversions(PairwiseTestKernel))

    PairwiseTestKernel{T<:AbstractFloat}(a::Argument{T}=1.0) = PairwiseTestKernel{T}(Variable(a))
    MLKernels.pairwise(κ::PairwiseTestKernel, x::Real, y::Real) = x*y
    function MLKernels.unsafe_pairwise(κ::PairwiseTestKernel, x::AbstractArray, y::AbstractArray)
        prod(x)*prod(y)
    end

    immutable TestClass{T<:AbstractFloat} <: CompositionClass{T}
        a::MLKernels.HyperParameter{T}
        TestClass(a::Variable{T}) = new(HyperParameter(a, unbounded(T)))
    end
    TestClass{T<:AbstractFloat}(a::Argument{T}=1.0) = TestClass{T}(Variable(a))
end

include("reference.jl")

include("hyperparameter.jl")
include("pairwisekernel.jl")
include("kernelcomposition.jl")
include("kerneloperation.jl")
include("pairwisefunction.jl")
include("kernelfunction.jl")


text:
@propagate_inbounds function Base.filter(f::Function, r::UnionRange)
    return isempty(r) ? empty(r) : @inbounds(r[find_all(f, r)])
end

@propagate_inbounds Base.filter(f::BitAnd, r::UnionRange) = _fltr_and(f, r, order(r))
@propagate_inbounds function _fltr_and(f, r, ro)
    _fltr_and(r, ro, find_all(f.f1, r, ro), find_all(f.f2, r, ro))
end
function _fltr_and(r, ro, inds1, inds2)
    if isempty(inds1)
        return empty(r)
    else
        return isempty(inds2) ? empty(r) : @inbounds(r[intersect(inds1, inds2)])
    end
end

@propagate_inbounds Base.filter(f::BitOr, r::UnionRange) = _fltr_or(f, r, order(r))
@propagate_inbounds function _fltr_or(f, r, ro)
    return _fltr_or(r, ro, find_all(f.f1, r, ro), find_all(f.f2, r, ro))
end
function _fltr_or(r, ro, inds1, inds2)
    if isempty(inds1)
        return isempty(inds2) ? empty(r) : @inbounds(r[inds2])
    else
        if isempty(inds2)
            return @inbounds(r[inds1])
        else
            if is_after(inds1, inds2)
                return vcat(@inbounds(inds2), @inbounds(inds1))
            elseif is_before(inds1, inds2)
                return vcat(@inbounds(r[inds1]), @inbounds(r[inds2]))
            else
                return @inbounds(r[_group_min(inds1, Forward, inds2, Forward):_group_max(inds1, Forward, inds2, Forward)])
            end

        end
    end
end

@propagate_inbounds function Base.filter(f::Fix2{typeof(!=)}, r::UnionRange)
    return _fltr_not(r, find_all(<(f.x), x), find_all(>(f.x), x))
end

function _fltr_not(r, inds1, inds2)
    if isempty(inds1)
        return isempty(inds2) ? empty(r) : @inbounds(r[inds2])
    else
        if isempty(inds2)
            return @inbounds(r[inds1])
        else
            return vcat(@inbounds(r[inds1]), @inbounds(r[inds2]))
        end
    end
end


text:
# "Generating gamma variates by a modified rejection technique"
# J.H. Ahrens, U. Dieter
# Communications of the ACM, Vol 25(1), 1982, pp 47-54
# doi:10.1145/358315.358390

# suitable for shape >= 1.0

struct GammaGDSampler{T<:Real} <: Sampleable{Univariate,Continuous}
    a::T
    s2::T
    s::T
    i2s::T
    d::T
    q0::T
    b::T
    σ::T
    c::T
    scale::T
end

function GammaGDSampler(g::Gamma{T}) where {T}
    a = shape(g)
    # Step 1
    s2 = a - 0.5
    s = sqrt(s2)
    i2s = 0.5/s
    d = 5.656854249492381 - 12.0s # 4*sqrt(2) - 12s

    # Step 4
    ia = 1.0/a
    q0 = ia*@horner(ia,
                    0.0416666664,
                    0.0208333723,
                    0.0079849875,
                    0.0015746717,
                    -0.0003349403,
                    0.0003340332,
                    0.0006053049,
                    -0.0004701849,
                    0.0001710320)

    if a <= 3.686
        b = 0.463 + s + 0.178s2
        σ = 1.235
        c = 0.195/s - 0.079 + 0.16s
    elseif a <= 13.022
        b = 1.654 + 0.0076s2
        σ = 1.68/s + 0.275
        c = 0.062/s + 0.024
    else
        b = 1.77
        σ = 0.75
        c = 0.1515/s
    end

    GammaGDSampler(T(a), T(s2), T(s), T(i2s), T(d), T(q0), T(b), T(σ), T(c), scale(g))
end

function calc_q(s::GammaGDSampler, t)
    v = t*s.i2s
    if abs(v) > 0.25
        return s.q0 - s.s*t + 0.25*t*t + 2.0*s.s2*log1p(v)
    else
        return s.q0 + 0.5*t*t*(v*@horner(v,
                                         0.333333333,
                                         -0.249999949,
                                         0.199999867,
                                         -0.1666774828,
                                         0.142873973,
                                         -0.124385581,
                                         0.110368310,
                                         -0.112750886,
                                         0.10408986))
    end
end

function rand(rng::AbstractRNG, s::GammaGDSampler)
    # Step 2
    t = randn(rng)
    x = s.s + 0.5t
    t >= 0.0 && return x*x*s.scale

    # Step 3
    u = rand(rng)
    s.d*u <= t*t*t && return x*x*s.scale

    # Step 5
    if x > 0.0
        # Step 6
        q = calc_q(s, t)
        # Step 7
        log1p(-u) <= q && return x*x*s.scale
    end

    # Step 8
    t = 0.0
    while true
        e = 0.0
        u = 0.0
        while true
            e = randexp(rng)
            u = 2.0rand(rng) - 1.0
            t = s.b + e*s.σ*sign(u)
            # Step 9
            t ≥ -0.718_744_837_717_19 && break
        end

        # Step 10
        q = calc_q(s, t)

        # Step 11
        (q > 0.0) && (s.c*abs(u) ≤ expm1(q)*exp(e-0.5t*t)) && break
    end

    # Step 12
    x = s.s+0.5t
    return x*x*s.scale
end

# "Computer methods for sampling from gamma, beta, poisson and bionomial distributions"
# J.H. Ahrens and U. Dieter
# Computing, 1974, Volume 12(3), pp 223-246
# doi:10.1007/BF02293108

# valid for 0 < shape <= 1
struct GammaGSSampler <: Sampleable{Univariate,Continuous}
    a::AbstractFloat
    ia::AbstractFloat
    b::AbstractFloat
    scale::AbstractFloat
end

function GammaGSSampler(d::Gamma)
    a = shape(d)
    ia = 1.0 / a
    b = 1.0+0.36787944117144233 * a
    GammaGSSampler(a, ia, b, scale(d))
end

function rand(rng::AbstractRNG, s::GammaGSSampler)
    while true
        # step 1
        p = s.b*rand(rng)
        e = randexp(rng)
        if p <= 1.0
            # step 2
            x = exp(log(p)*s.ia)
            e < x || return s.scale*x
        else
            # step 3
            x = -log(s.ia*(s.b-p))
            e < log(x)*(1.0-s.a) || return s.scale*x
        end
    end
end

# "A simple method for generating gamma variables"
# G. Marsaglia and W.W. Tsang
# ACM Transactions on Mathematical Software (TOMS), 2000, Volume 26(3), pp. 363-372
# doi:10.1145/358407.358414
# http://www.cparity.com/projects/AcmClassification/samples/358414.pdf

struct GammaMTSampler <: Sampleable{Univariate,Continuous}
    d::AbstractFloat
    c::AbstractFloat
    κ::AbstractFloat
end

function GammaMTSampler(g::Gamma)
    d = shape(g) - 1/3
    c = 1.0 / sqrt(9.0 * d)
    κ = d * scale(g)
    GammaMTSampler(d, c, κ)
end

function rand(rng::AbstractRNG, s::GammaMTSampler)
    while true
        x = randn(rng)
        v = 1.0 + s.c * x
        while v <= 0.0
            x = randn(rng)
            v = 1.0 + s.c * x
        end
        v *= (v * v)
        u = rand(rng)
        x2 = x * x
        if u < 1.0 - 0.331 * abs2(x2) || log(u) < 0.5 * x2 + s.d * (1.0 - v + log(v))
            return v*s.κ
        end
    end
end

# Inverse Power sampler
# uses the x*u^(1/a) trick from Marsaglia and Tsang (2000) for when shape < 1
struct GammaIPSampler{S<:Sampleable{Univariate,Continuous},T<:Real} <: Sampleable{Univariate,Continuous}
    s::S #sampler for Gamma(1+shape,scale)
    nia::T #-1/scale
end

function GammaIPSampler(d::Gamma,::Type{S}) where S<:Sampleable
    GammaIPSampler(Gamma(1.0 + shape(d), scale(d)), -1.0 / shape(d))
end
GammaIPSampler(d::Gamma) = GammaIPSampler(d,GammaMTSampler)

function rand(rng::AbstractRNG, s::GammaIPSampler)
    x = rand(rng, s.s)
    e = randexp(rng)
    x*exp(s.nia*e)
end

# function sampler(d::Gamma)
#     if d.shape < 1.0
#         # TODO: d.shape = 0.5 : use scaled chisq
#         GammaIPSampler(d)
#     elseif d.shape == 1.0
#         Exponential(d.scale)
#     else
#         GammaGDSampler(d)
#     end
# end


text:
# It is suggested that test code for MLJ.jl include files be placed in
# a file of the same name under "test/" (and included below) and that
# this test code be wrapped in a module. Any new module name will do -
# eg, `module TestDatasets` for code testing `datasets.jl`.

using Distributed
addprocs(2)

@everywhere begin
using MLJ
using MLJBase
using Test
using Random
end

@testset "utilities" begin
  @test include("utilities.jl")
end

@testset "parameters" begin
  @test include("parameters.jl")
end

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

@testset "tuning" begin
    @test include("tuning.jl")
end

@testset "learning_curves" begin
    @test include("learning_curves.jl")
end

@testset "ensembles" begin
    @test include("ensembles.jl")
end

@testset "matching models to data" begin
    @test include("model_matching.jl")
end

@testset "tasks" begin
  @test include("tasks.jl")
end

@testset "scitypes" begin
    @test include("scitypes.jl")
end


text:
using DynamicWalking2018
using Base.Test

const module_tempdir = joinpath(Base.tempdir(), string(module_name(DynamicWalking2018)))

using Conda
if is_windows()
    Conda.add("nbconvert") # jupyter doesn't seem to come with nbconvert on windows?
end

@testset "example notebooks" begin
    using IJulia

    outputdir = module_tempdir
    if !isdir(outputdir)
        mkpath(outputdir)
    end
    jupyter = IJulia.jupyter
    for f in filter(x -> endswith(x, "ipynb"), readdir("../notebooks"))
        if contains(f, "Julia is Fast") && !is_linux()
            # Can't run the gcc demo on osx or windows
            continue
        end
        notebook = joinpath("..", "notebooks", f)
        output = joinpath(outputdir, f)
        @test begin run(`$jupyter nbconvert --to notebook --execute $notebook --output $output --ExecutePreprocessor.timeout=500`); true end
    end
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to import and transform raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# --------------- Import data ---------------------------------------------------------------------------------------------------------------------------------

"""
    import_R(dir::String, year::Integer)

The function loads the WIOD Input-Output Table of specified year as DataFrame into the environment and
    extracts the raw intermediate, final demand matrices (Z, F) and the inventory adjustment matrix IV.

# Arguments
- `dir::String`: directory of raw data.
- `year::Integer`: specifies the year of the WIOD table which should be imported.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.

# Examples
```julia-repl
julia> Z, F, IV = import_R(dir, 2012)
```
"""

function import_R(dir::String, year::Integer, N::Integer, S::Integer)

    # load data for specified year and transform into a DataFrame (automatically)
    path = dir * "WIOT" * string(year) * "_October16_ROW.RData"
    df = RData.load(path)["wiot"] # 2472×2690 DataFrame
    transform!(df, [:Year, :RNr] .=> ByRow(Int64) .=> [:Year, :RNr], renamecols=false) # Years, row-industry-identifier as type: Int64

    # Intermediate demand
    Z = df[1:N*S, 6:N*S+5] # NS×NS
    Z = Matrix(convert.(Float64, Z))

    # Final demand
    F = df[1:N*S, N*S+6:end-1] # NS×5*N

    # Subset final demand for inventory adjustments
    inventory_columns = 5:5:5*N
    IV = F[:, inventory_columns] # NS×N
    F = F[:, Not(inventory_columns)] # NS×4*N

    IV = Matrix(convert.(Float64, IV)) # NS×N
    F = Matrix(convert.(Float64, F)) # NS×4*N

    return Z, F, IV
end

"""
    import_excel(dir::String, year::Integer)

The function loads the WIOD Input-Output Table of specified year as DataFrame into the environment and
    extracts the raw intermediate, final demand matrices (Z, F) and the inventory adjustment matrix IV.

# Arguments
- `dir::String`: directory of raw data.
- `year::Integer`: specifies the year of the WIOD table which should be imported.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.

# Examples
```julia-repl
julia> Z, F, IV = import_excel(dir, 2012)
```
"""

function import_excel(dir::String, year::Integer, N::Integer, S::Integer)

    # load data for specified year and transform into a DataFrame (automatically)
    path = dir * "WIOT" * string(year)[end-1:end] * "_ROW_Apr12.xlsx"
    df = DataFrame(XLSX.readxlsx(path)["WIOT_$year"][:], :auto)
    df = df[7:end,5:end]

    # Intermediate demand
    Z = df[1:N*S, 1:N*S] # NS×NS
    Z = Matrix(convert.(Float64, Z))

    # Final demand
    F = df[1:N*S, N*S+1:end-1] # NS×5*N

    # Subset final demand for inventory adjustments
    inventory_columns = 5:5:5*N
    IV = F[:, inventory_columns] # NS×N
    F = F[:, Not(inventory_columns)] # NS×4*N

    IV = Matrix(convert.(Float64, IV)) # NS×N
    F = Matrix(convert.(Float64, F)) # NS×4*N

    return Z, F, IV
end

# --------------- Transformations -----------------------------------------------------------------------------------------------------------------------------

"""
    create_matrices(df::DataFrame, N::Integer, S::Integer)

The function creates the origin country-industry destination country-industry intermediate demand matrix Z and
    the country-industry VA vector.

# Arguments
- `Z::Matrix`: NS×NS, raw origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, raw origin country-industry destination country final demand matrix F.
- `IV::Matrix`: NS×N, raw origin country-industry destination country inventory adjustment matrix IV.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, origin country-industry destination country final demand matrix F.
- `Y::Vector{Float64}`: NS×1, inventory adjusted country-industry gross output vector Y.
- `VA::Vector{Float64}`: NS×1, country-industry value added vector VA.

# Examples
```julia-repl
julia> Z, F, Y, VA = create_matrices(WIOD_2012, N, S)
```
"""
function create_matrices(Z::Matrix, F::Matrix, IV::Matrix, N::Integer, S::Integer)

    # remove zero and negative entries (page 37 Antras and Chor, 2018)
    #Z = ifelse.(Z .<= 0.0, 1e-18, Z)
    #F = ifelse.(F .<= 0.0, 1e-18, F)

    IV = [sum(IV[i,:]) for i in 1:N*S] # NS×1

    # Compute gross output from intermediate and final demand
    Y = [sum(Z[i,:]) + sum(F[i,:]) + IV[i] for i in 1:N*S] # NS×1
    Y = ifelse.(abs.(Y) .< 1e-18, 0.0, Y)

    # Inventory adjustment (spread IV regardless of origin over all columns)
    adj = Y .- IV # NS×1
    adj = ifelse.(adj .< 0.0, 0.0, adj) # in case negative

    Z = Z .* (repeat(Y, 1, N*S) ./ repeat(adj, 1, N*S)) # NS×NS
    Z = ifelse.(isnan.(Z), 0.0, Z)
    Z = ifelse.(isinf.(Z), 0.0, Z)
    Z = ifelse.(Z .< 0.0, 0.0, Z)

    F = F .* (repeat(Y, 1, N*4) ./ repeat(adj, 1, N*4)) # NS×N*4, inventory adjustment
    F = ifelse.(isnan.(F), 0.0, F)
    F = ifelse.(isinf.(F), 0.0, F)
    F = ifelse.(F .< 0.0, 0.0, F)
    F = [sum(F[i,j:j+3]) for i in 1:N*S, j in 1:4:N*4] # NS×N

    # Use computed values to have internal consistency
    Y = [sum(Z[i,:]) + sum(F[i,:]) for i in 1:N*S] # NS×1
    VA = [Y[i] - sum(Z[:,i]) for i in 1:N*S] # NS×1, compute value added (gross output minus inputs)
    VA = ifelse.(VA .< 0.0, 1e-18, VA) # in case one is negative

    return Z, F, Y, VA
end

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

"""
    create_trade_shares(Z::Matrix, F::Matrix, Y::Vector, N::Integer, S::Integer)

The function computes the origin country-industry destination country intermediate goods trade share matrix π_Z and
    origin country-industry destination country final goods trade share matrix π_F.

# Arguments
- `Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, origin country-industry destination country final demand matrix F.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `π_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix{Float64}`: NS×N, origin country-industry destination country foods goods trade share matrix π_Z.

# Examples
```julia-repl
julia> π_Z, π_F = create_trade_shares(Z, F, N, S)
```
"""
function create_trade_shares(Z::Matrix, F::Matrix, N::Integer, S::Integer)

    Z_agg = [sum(Z[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N*S] # S×NS
    π_Z = Z ./ repeat(Z_agg, N) # NS×NS
    π_Z = ifelse.(isnan.(π_Z), 0.0, π_Z)

    F_agg = [sum(F[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N] # S×N
    π_F = F ./ repeat(F_agg, N) # NS×N
    π_F = ifelse.(isnan.(π_F), 0.0, π_F)

    return π_Z, π_F
end

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

"""
    create_expenditure_shares(Z::Matrix, F::Matrix, Y::Vector, N::Integer, S::Integer)

The function computes the country value added vector VA_ctry respecting inventory adjustments, intermediate goods expenditure share matrix γ
    and final goods import expenditure share matrix α. Additionally, country-industry value added coeffiecient matrix VA_coeff as well as
    country trade balance vector TB_ctry and country final goods import vector F_ctry are produced.

# Arguments
- `Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, origin country-industry destination country final demand matrix F.
- `Y::Vector`: NS×1, origin country-industry gross output vector Y.
- `VA::Vector`: NS×1, country-industry value added vector VA.
- `π_Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix`: NS×N, origin country-industry destination country foods goods trade share matrix π_Z.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `γ::Matrix{Float64}`: S×NS, country-industry intermediate input expenditure share matrix γ.
- `α::Matrix{Float64}`: S×N, country-industry final expenditure expenditure share matrix α.
- `VA_coeff::Vector{Float64}`: NS×1, country-industry value added coefficients vector VA_coeff.
- `TB_ctry::Vector{Float64}`: NS×1, country trade balance (X-M) vector TB.
- `VA_ctry::Vector{Float64}`: N×1, country value added vector VA_ctry.
- `F_ctry::Vector{Float64}`: N×1, country final import demand vector VA_ctry.

# Examples
```julia-repl
julia> γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry = create_expenditure_shares(Z, F, Y, π_Z, π_F, VA, N, S)
```
"""
function create_expenditure_shares(Z::Matrix, F::Matrix, Y::Vector, VA::Vector, π_Z::Matrix, π_F::Matrix, N::Integer, S::Integer)

    # Country-industry level input expenditure shares
    # sum over reporting country only (eq. 37) => sum over column gives total intermediate consumption per country-industry (+ VA = Y)
    Z_agg = [sum(Z[i:S:(N-1)*S+i, j]) for i in 1:S, j in 1:N*S]
    Z_agg_VA = [Z_agg; VA']
    γ_VA = [Z_agg_VA[i,j]/sum(Z_agg_VA[:,j]) for i in 1:S+1, j in 1:N*S]
    γ = γ_VA[1:S, 1:N*S] # S×NS
    γ = ifelse.(isnan.(γ), 0.0, γ) # some entries in potentially equal to zero => NaN => assume 0
    γ = ifelse.(γ .> 1.0, 1.0, γ) # if value explode attribute reduce to 1.0 (for country-sector where VA is negative)

    # Country-industry level value added coefficients
    VA_coeff = γ_VA[S+1, 1:N*S] # NS×1
    VA_coeff = ifelse.(isnan.(VA_coeff), 1.0, VA_coeff) # some entries in potentially equal to zero => NaN, assume VA coefficient to be 1
    VA_coeff = ifelse.(VA_coeff .< 0.0, 0.0, VA_coeff)

    # Country level value added
    VA_ctry = [sum(VA[i:i+S-1]) for i in 1:S:N*S] # N×1

    # Final import demand at country level (sum over rows not columns - would give export demand)
    F_ctry = [sum(F[:,j]) for j in 1:N] # N×1

    # Country-industry final consumption expenditure (import) shares => columns sum to 1
    α = [sum(F[i:S:(N-1)*S+i, j])/F_ctry[j] for i in 1:S, j in 1:N] # S×N

    # Country level trade balance (exports - imports)
    # used computation in Antras and Chor (2018) quite different
    # need to calculate country level exports/imports (row/columns) (i.e. aggregate and remove intra-country trade)
    E = [sum(Y[i:i+S-1]) - sum(Z[i:i+S-1,i:i+S-1]) - sum(F[i:i+S-1,ceil(Int, i/S)]) for i in 1:S:N*S] # N×1

    M_Z = [sum(Z[:,j:j+S-1]) - sum(Z[j:j+S-1,j:j+S-1]) for j in 1:S:N*S]
    M_F = [sum(F[:,ceil(Int, j/S)]) - sum(F[j:j+S-1,ceil(Int, j/S)]) for j in 1:S:N*S]
    M = M_Z .+ M_F # N×1

    TB_ctry = E .- M # N×1

    # #### Antras and Chor (2018) calculation
    # E_A = zeros(N)
    # M_A = zeros(N)
    # for j in 1:N
    #     for s in 1:N
    #         for i in 1:N
    #             iiss = (i-1)*S+s
    #             jjss = (j-1)*S+s
    #             for r in 1:S
    #                 jjrr = (j-1)*S+r
    #                 iirr = (i-1)*S+r
    #                 M_A[j] += π_Z[iiss,jjrr]*γ[s,jjrr]*Y[jjrr]
    #                 E_A[j] += π_Z[jjss,iirr]*γ[s,iirr]*Y[iirr]
    #             end
    #             M_A[j] += π_F[iiss,j]*α[s,j]*F_ctry[j]
    #             E_A[j] += π_F[jjss,i]*α[s,i]*F_ctry[i]
    #         end
    #     end
    # end

    # TB_ctry = E_A .- M_A # N×1

    #TB_ctry == VA_ctry .- F_ctry # must hold, holds much better with own calculation?

    return γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Summary function
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    transform_WIOD_2016(dir::String, year::Integer)

The function imports and transforms the raw data for further use in the model.

# Arguments
- `dir::String`: directory of raw data.
- `year::Integer`: specifies the year of the WIOD table which should be imported.

# Output
- `Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix{Float64}`: NS×N, origin country-industry destination country final demand matrix F.
- `Y::Vector{Float64}`: NS×1, origin country-industry gross output vector Y.
- `F_ctry::Vector{Float64}`: N×1, country final import demand vector VA_ctry.
- `TB_ctry::Vector{Float64}`: NS×1, country trade balance (X-M) vector TB.
- `VA_ctry::Vector{Float64}`: N×1, country value added vector VA_ctry.
- `VA_coeff::Vector{Float64}`: NS×1, country-industry value added coefficients vector VA_coeff.
- `γ::Matrix{Float64}`: S×NS, country-industry intermediate input expenditure share matrix γ.
- `α::Matrix{Float64}`: S×N, country-industry final expenditure expenditure share matrix α.
- `π_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate goods trade share matrix π_Z.
- `π_F::Matrix{Float64}`: NS×N, origin country-industry destination country final goods trade share matrix π_F.

# Examples
```julia-repl
julia> Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F = transform_WIOD_2016(dir, 2014)
```
"""
function transform_WIOD_2016(dir::String, format::String, year::Integer, N::Integer, S::Integer)

    if format == "excel"
        Z, F, IV = import_excel(dir, year, N, S)
    else
        Z, F, IV = import_R(dir, year, N, S)
    end

    Z, F, Y, VA = create_matrices(Z, F, IV, N, S)
    π_Z, π_F = create_trade_shares(Z, F, N, S)
    γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry = create_expenditure_shares(Z, F, Y, VA, π_Z, π_F, N, S)

    return Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F
end

max_stars_repo_path: old stuff/transform_WIOD_2016_3.jl

max_issues_repo_path: old stuff/transform_WIOD_2016_3.jl

max_forks_repo_path: old stuff/transform_WIOD_2016_3.jl


text:
using Base.Cartesian

function mul_suboperator(a, b, index::Int)
    result = Ket(b.basis)
    result_data = reshape(result.data, reverse(result.shape)...)
    b_data
    for m=1:size(a.data,2)

    end
end

# println(macroexpand(:(@nloops 3 i d->(1:b.shape[d]) d->(if d==index && i_d!=m continue end) begin
#             for n=1:size(b.data,1)
#                 (@nref 3 result_data d->(4-d==index ? n : i_{4-d})) += a[n,m] * (@nref 3 b_date d->(i_{4-d}))
#             end
#         end)))

# # println(macroexpand(:(@nloops 3 i b begin
# #                 (@nref 3 result_data i) += a[n,m] * (@nref 3 b_date i)
# #         end)))
# println(macroexpand(:(function mul_suboperator(a::Operator, b::Ket, index::Int)
#     result = Ket(b.basis)
#     result_data = reshape(result.data, reverse(result.basis.shape)...)
#     b_data = reshape(b.data, reverse(result.basis.shape)...)
#     for m=1:size(a.data,2)
#         @nloops 3 i d->(1:b.basis.shape[d]) d->(if d==index && i_d!=m continue end) begin
#             for n=1:size(b.data,1)
#                 (@nref 3 result_data d->(4-d==index ? n : i_{4-d})) += a.data[n,m] * (@nref 3 b_data d->(i_{4-d}))
#             end
#         end
#     end
# end)))

# println(macroexpand(:(@ngenerate N Int function gemm_r!{N}(alpha, a::Array{Complex128, N}, b::Array{Complex128, 2},
#                                         index::Int, beta, result::Array{Complex128, N})
#     for n=1:size(b,1)
#         @nloops N i d->(1:size(a,N+1-d)) d->(if d==index && i_d!=n continue end) begin
#             for m=1:size(b,2)
#                  (@nref N result d->(N+1-d==index ? m : i_{N+1-d})) = beta*(@nref N result d->(N+1-d==index ? m : i_{N+1-d})) + alpha*(@nref N a d->(i_{N+1-d})) * b[n,m]
#             end
#         end
#     end
#     return 0
# end)))
# println(macroexpand(:(@ngenerate N Int function mul_suboperator_r{N}(a::Array{Complex128, N}, b::Array{Complex128, 2}, index::Int, result::Array{Complex128, N})
#     for n=1:size(b,1)
#         @nloops N i d->(1:size(a,N+1-d)) d->(if d==index && i_d!=n continue end) begin
#                 (@nref N result (d->(i_{N+1-d}))) = alpha * (@nref N psi (d->(i_{N+1-d}))) * op[m,n]
#         end
#     end
#     return 0
# end)))

type Index{I}
    data::Type{Array{Int,I}}
end

Index(i::Int) = Index(Array{Int,i})

#println(macroexpand(:(
@ngenerate RANK Nothing function t{RANK}(x::Type{Type{RANK}}, shape::Vector{Int}, strides::Vector{Int}, index_N::Int, stride_index::Int, a::Matrix, b::Vector, result::Vector)
    @nloops RANK i (d->1:shape[d])  begin
        for m =1:index_N
            for n=1:index_N
                N_n = (n-1)*stride_index + 1
                N_m = (m-1)*stride_index + 1
                @nexprs RANK (d->(s = (i_d-1)*strides[d]; N_n+=s; N_m+=s))
                result[N_m] += a[m,n]*b[N_n]
            end
        end
    end
    return nothing
end
#)))

function strides(shape::Vector{Int})
    N = length(shape)
    S = zeros(Int, N)
    S[N] = 1
    for i=N-1:-1:1
        S[i] = S[i+1]*shape[i+1]
    end
    return S
end

a = [-1  4 3; -3 2 1; 6 8 5]
v1 = [3,5]
v2 = [2,7,1]
v3 = [1:10000000]

v = kron(kron(v1,v2),v3)
shape = [length(v1), length(v2), length(v3)]

S = strides(shape)
index = 2

shape_reduced = [shape[1:index-1], shape[index+1:end]]
S_reduced = [S[1:index-1], S[index+1:end]]
stride_index = S[index]
index_N = shape[index]

println("shape", shape)
println("shape_reduced", shape_reduced)
println("stride", S)
println("stride_reduced", S_reduced)
println("index_N", index_N)
println("stride_index", stride_index)

result = zeros(Int,prod(shape))
@time t(Type{length(S)-1}, shape_reduced, S_reduced, index_N, stride_index, a, v, result)
result = zeros(Int,prod(shape))
@time t(Type{length(S)-1}, shape_reduced, S_reduced, index_N, stride_index, a, v, result)

@time kron(kron(v1,a*v2), v3)
@time kron(kron(v1,a*v2), v3)
#Index(3)
#t4(Index{3})

# @ngenerate N Nothing function t{N}(a::Array{Complex128, 2}, b::Array{Complex128, N}, index::Int, result::Array{Complex128, N})
#     for m=1:size(a,2)
#         @nloops N i d->(1:size(b,N+1-d)) d->(if d==index && i_d!=m continue end) begin
#             for n=1:size(a,1)
#                 #1+1
#                 #(@nref N result d->(N+1-d==index ? n : i_{N+1-d})) += 1#a[n,m] * (@nref N b d->(i_{N+1-d}))
#                 @time (@nref 2 a i) = 0
#                 @time (@nref 2 a i) = 0
#                 error()
#             end
#         end
#     end
#     return nothing
# end

# a = rand(Complex128, 2,2)
# b = rand(Complex128, 2,2,2,2,2,2,2,2,2)
# result = rand(Complex128, 2,2,2,2,2,2,2,2,2)

# f{N}(b::Array{Complex128,N}) = (b[2,2,2,2,2]=1; nothing)
# #@time f(b)
# #@time f(b)
# @time t(a,b,8, result)
# @time t(a,b,8, result)


text:
# get and set current directory

function pwd()
    b = Array(Uint8,1024)
    len = Csize_t[length(b),]
    uv_error(:getcwd, ccall(:uv_cwd, Cint, (Ptr{Uint8}, Ptr{Csize_t}), b, len))
    bytestring(b[1:len[1]-1])
end

function cd(dir::String)
    uv_error("chdir $dir", ccall(:uv_chdir, Cint, (Ptr{Uint8},), bytestring(dir)))
end
cd() = cd(homedir())

@unix_only function cd(f::Function, dir::String, args...)
    fd = ccall(:open,Int32,(Ptr{Uint8},Int32),".",0)
    systemerror(:open, fd == -1)
    try
        cd(dir)
        f(args...)
    finally
        systemerror(:fchdir, ccall(:fchdir,Int32,(Int32,),fd) != 0)
        systemerror(:close, ccall(:close,Int32,(Int32,),fd) != 0)
    end
end
@windows_only function cd(f::Function, dir::String, args...)
    old = pwd()
    try
        cd(dir)
        f(args...)
   finally
        cd(old)
    end
end
cd(f::Function) = cd(f, homedir())

function mkdir(path::String, mode::Unsigned=0o777)
    @unix_only ret = ccall(:mkdir, Int32, (Ptr{Uint8},Uint32), bytestring(path), mode)
    @windows_only ret = ccall(:_wmkdir, Int32, (Ptr{Uint16},), utf16(path))
    systemerror(:mkdir, ret != 0)
end

function mkpath(path::String, mode::Unsigned=0o777)
    isdirpath(path) && (path = dirname(path))
    dir = dirname(path)
    (path == dir || isdir(path)) && return
    mkpath(dir, mode)
    mkdir(path)
end

mkdir(path::String, mode::Signed) = error("mode must be an unsigned integer; try 0o$mode")
mkpath(path::String, mode::Signed) = error("mode must be an unsigned integer; try 0o$mode")

function rmdir(path::String)
    @unix_only ret = ccall(:rmdir, Int32, (Ptr{Uint8},), bytestring(path))
    @windows_only ret = ccall(:_wrmdir, Int32, (Ptr{Uint16},), utf16(path))
    systemerror(:rmdir, ret != 0)
end

# The following use Unix command line facilites

rm(path::String) = FS.unlink(path)
cp(src::String, dst::String) = FS.sendfile(src, dst)
mv(src::String, dst::String) = FS.rename(src, dst)
touch(path::String) = run(`touch $path`)

# Obtain a temporary filename.
@unix_only function tempname()
    d = get(ENV, "TMPDIR", C_NULL) # tempnam ignores TMPDIR on darwin
    p = ccall(:tempnam, Ptr{Uint8}, (Ptr{Uint8},Ptr{Uint8}), d, "julia")
    systemerror(:tempnam, p == C_NULL)
    s = bytestring(p)
    c_free(p)
    return s
end

# Obtain a temporary directory's path.
@unix_only tempdir() = dirname(tempname())

# Create and return the name of a temporary file along with an IOStream
@unix_only function mktemp()
    b = joinpath(tempdir(), "tmpXXXXXX")
    p = ccall(:mkstemp, Int32, (Ptr{Uint8}, ), b) # modifies b
    return (b, fdio(p, true))
end

@windows_only begin
function tempdir()
    temppath = Array(Uint16,32767)
    lentemppath = ccall(:GetTempPathW,stdcall,Uint32,(Uint32,Ptr{Uint16}),length(temppath),temppath)
    if lentemppath >= length(temppath) || lentemppath == 0
        error("GetTempPath failed: $(FormatMessage())")
    end
    resize!(temppath,lentemppath+1)
    return utf8(UTF16String(temppath))
end
tempname(uunique::Uint32=uint32(0)) = tempname(tempdir(), uunique)
function tempname(temppath::String,uunique::Uint32)
    tname = Array(Uint16,32767)
    uunique = ccall(:GetTempFileNameW,stdcall,Uint32,(Ptr{Uint16},Ptr{Uint16},Uint32,Ptr{Uint16}),
        utf16(temppath),utf16("jul"),uunique,tname)
    lentname = findfirst(tname,0)-1
    if uunique == 0 || lentname <= 0
        error("GetTempFileName failed: $(FormatMessage())")
    end
    resize!(tname,lentname+1)
    return utf8(UTF16String(tname))
end
function mktemp()
    filename = tempname()
    return (filename, open(filename,"r+"))
end
end

# Create and return the name of a temporary directory
@unix_only function mktempdir()
    b = joinpath(tempdir(), "tmpXXXXXX")
    p = ccall(:mkdtemp, Ptr{Uint8}, (Ptr{Uint8}, ), b)
    return bytestring(p)
end

@windows_only function mktempdir()
    seed::Uint32 = rand(Uint32)
    dir = tempdir()
    while true
        if uint16(seed) == 0
            seed += 1
        end
        filename = tempname(dir, seed)
        ret = ccall(:_wmkdir, Int32, (Ptr{Uint16},), utf16(filename))
        if ret == 0
            return filename
        end
        systemerror(:mktempdir, errno()!=EEXIST)
        seed += 1
    end
end

function readdir(path::String)
    # Allocate space for uv_fs_t struct
    uv_readdir_req = zeros(Uint8, ccall(:jl_sizeof_uv_fs_t, Int32, ()))

    # defined in sys.c, to call uv_fs_readdir, which sets errno on error.
    file_count = ccall(:jl_readdir, Int32, (Ptr{Uint8}, Ptr{Uint8}),
                       bytestring(path), uv_readdir_req)
    systemerror("unable to read directory $path", file_count < 0)

    # The list of dir entries is returned as a contiguous sequence of null-terminated
    # strings, the first of which is pointed to by ptr in uv_readdir_req.
    # The following lines extracts those strings into dirent
    entries = ByteString[]
    offset = 0

    for i = 1:file_count
        entry = bytestring(ccall(:jl_uv_fs_t_ptr_offset, Ptr{Uint8},
                                 (Ptr{Uint8}, Int32), uv_readdir_req, offset))
        push!(entries, entry)
        offset += sizeof(entry) + 1   # offset to the next entry
    end

    # Clean up the request string
    ccall(:jl_uv_fs_req_cleanup, Void, (Ptr{Uint8},), uv_readdir_req)

    entries
end

readdir() = readdir(".")

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
using SVD_Projection
using Test, Random
Random.seed!(1234)
discovery = rand([0., 1., 2.], 30, 20);
target= rand([0., 1., 2.], 30, 15);

V_PC = project_PC(discovery, target, N_PC=3)
@test V_PC[1,:] == [-0.12869065723376472; -0.05776618956235076; -0.12278055691277556]


text:
# GEMM
#

macro blis_linalg_gemm(typename, targetfunc, bliname)

    # Get method for the typed API backend.
    blifuncname = Symbol("bli_", cblas_typechar[ eval(typename) ], bliname)
    blifunc = getproperty(TypedBackend, blifuncname)

    return quote
        """
        BLIS-based GEMM with approximately the same interface as `LinearAlgebra.BLAS.gemm!`.
        Strided arrays are directly supported without further Julia-level dispatching.
        """
        $(esc(targetfunc))(tA::AbstractChar,
                           tB::AbstractChar,
                           alpha::$typename,
                           A::StridedMatrix{$typename},
                           B::StridedMatrix{$typename},
                           beta::$typename,
                           C::StridedMatrix{$typename}) = begin

            bli_tA = char_to_trans[tA]
            bli_tB = char_to_trans[tB]

            m,  k = size(A)
            k_, n = size(B)
            m_, n_= size(C)

            # Transpose sizes.
            if bli_tA == BLIS_TRANSPOSE || bli_tA == BLIS_CONJ_TRANSPOSE
                m, k = k, m
            end
            if bli_tB == BLIS_TRANSPOSE || bli_tB == BLIS_CONJ_TRANSPOSE
                k_, n = n, k_
            end

            (k == k_) || throw(DimensionMismatch("GEMM contracted size mismatch."))
            (m == m_ &&
             n == n_) || throw(DimensionMismatch("GEMM target buffer size mismatch."))

            $blifunc(bli_tA,
                     bli_tB,
                     m, n, k,
                     [alpha],
                     A, strides(A)...,
                     B, strides(B)...,
                     [beta],
                     C, strides(C)...)
            C

        end
    end
end

@blis_linalg_gemm(Float32, gemm!, gemm!)
@blis_linalg_gemm(Float64, gemm!, gemm!)
@blis_linalg_gemm(ComplexF32, gemm!, gemm!)
@blis_linalg_gemm(ComplexF64, gemm!, gemm!)


text:
Dict{Float32,String}(
19=>"North Europe",
8=>"Central Europe",
15=>"South Europe",
33=>"West South America",
2=>"Amazon",
18=>"Northeast Brazil",
26=>"Southeastern South America",
5=>"Central America/Mexico",
1=>"Alaska/N.W. Canada",
13=>"East North America",
32=>"West North America",
10=>"Central North America",
9=>"Canada/Greenland/Iceland" ,
11=>"East Africa",
30=>"West Asia" ,
29=>"West Africa" ,
22=>"Sahara",
21=>"Southern Africa",
7=>"Central Asia",
12=>"East Asia",
28=>"Tibetan Plateau",
16=>"North Asia",
25=>"Southeast Asia",
23=>"South Asia",
17=>"Northern Australia",
24=>"Southern Australia")


text:
using DelimitedFiles

function solve(input)
    raw = readdlm(input,'\n')
    coords = vcat(map(parseCoords,raw)...)
    best_time = 0
    smaller_height = height(coords,0)
    for time in 1:10460
        current_height = height(coords,time)
        if current_height < smaller_height
            smaller_height=current_height
            best_time = time
        end
    end
    @show best_time
    current = copy(coords[:,1:2])+best_time*coords[:,3:4]
    x_min,x_max = extrema(current[:,1])
    y_min, y_max = extrema(current[:,2])
    w = x_max-x_min+1
    h = y_max-y_min+1
    grid = fill(' ',(h,w))

    for point in 1:size(current)[1]
        grid[current[point,:][2]-y_min+1,current[point,:][1]-x_min+1]='#'
    end

    for line in 1:size(grid)[1]
        println(join(grid[line,:]))
    end

end

function parseCoords(s)
    m = match(r"position\=\<(?<x>.+),(?<y>.+)\> velocity\=\<(?<vx>.+),(?<vy>.+)\>", s)
    x = tryparse(Int64,m[:x])
    y = tryparse(Int64,m[:y])
    vx = tryparse(Int64,m[:vx])
    vy = tryparse(Int64,m[:vy])
    return [x y vx vy]
end

function height(coords,time)
    ymin,ymax = extrema(coords[:,2] + coords[:,4]*time)
    return ymax-ymin
end

solve("input.txt")


text:
# Configuration for running spatial simulation
export simtype
@everywhere simtype = 2
const num_trials = 1000
@everywhere const N_list = [6400,3200,1600,800,400,200,100,50,25]        # Meta-population size
N_mut_list = [0.125, 0.25, 0.5, 1.0]
#mutation_stddev_list = [0.005,0.01,0.02,0.04]
const num_subpops = 1                    # Number of subpopulations
const num_attributes_list = [1]        # number attributes for quantitative representation
const ngens = 1       # Generations after burn-in
const burn_in= 3.0    # generations of burn_in as a multiple of N
const ideal = 1.0
const neutral = true
const fit_max=1.0
const fit_slope=1.0
const additive_error=false
const wrap_attributes=false

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
import LinearAlgebra: kron, eigmin

const σx = [0.0 + 0.0im 1; 1 0]
const σy = [0.0 + 0.0im -1.0im; 1.0im 0]
const σz = [1.0 + 0.0im 0; 0 -1]
const σi = [1.0 + 0.0im 0; 0 1]
const σ₋ = [0.0im 0;1 0]
const σ₊ = [0.0im 1;0 0]
const σ = [σx, σy, σz, σi]
xvec = [[1.0 + 0.0im, 1.0] / sqrt(2), [1.0 + 0.0im, -1.0] / sqrt(2)]
yvec = [[1.0, 1.0im] / sqrt(2), [1.0, -1.0im] / sqrt(2)]
zvec = [[1.0 + 0.0im, 0], [0, 1.0 + 0.0im]]

const spσz = sparse(σz)
const spσx = sparse(σx)
const spσy = sparse(σy)
const spσi = sparse(I, 2, 2)

"""
    PauliVec

A constant that holds the eigenvectors of the Pauli matrices. Indices 1, 2 and 3 corresponds to the eigenvectors of ``σ_x``, ``σ_y`` and ``σ_z``.

# Examples
```julia-repl
julia> σx*PauliVec[1][1] == PauliVec[1][1]
true
```
"""
const PauliVec = [xvec, yvec, zvec]

"""
    ⊗(A, B)

Calculate the tensor product of `A` and `B`.

# Examples
```julia-repl
julia> σx⊗σz
4×4 Array{Complex{Float64},2}:
 0.0+0.0im   0.0+0.0im  1.0+0.0im   0.0+0.0im
 0.0+0.0im  -0.0+0.0im  0.0+0.0im  -1.0+0.0im
 1.0+0.0im   0.0+0.0im  0.0+0.0im   0.0+0.0im
 0.0+0.0im  -1.0+0.0im  0.0+0.0im  -0.0+0.0im
```
"""
⊗ = kron

"""
$(SIGNATURES)

Decompse matrix `mat` onto matrix basis `basis`

# Examples
```julia-repl
julia> matrix_decompose(1.0*σx+2.0*σy+3.0*σz, [σx,σy,σz])
3-element Array{Complex{Float64},1}:
 1.0 + 0.0im
 2.0 + 0.0im
 3.0 + 0.0im
```
"""
function matrix_decompose(mat::AbstractMatrix, basis::Vector{<:AbstractMatrix})
    dim = size(basis[1])[1]
    [tr(mat * b) / dim for b in basis]
end

"""
    check_positivity(m)

Check if matrix `m` is positive. Return `true` is the minimum eigenvalue of `m` is greater than or equal to 0.
"""
function check_positivity(m::AbstractMatrix)
    if !ishermitian(m)
        @warn "Input fails the numerical test for Hermitian matrix. Use the upper triangle to construct a new Hermitian matrix."
        d = Hermitian(m)
    else
        d = m
    end
    eigmin(d) >= 0
end

"""
    gibbs_state(h, β)

Calculate the Gibbs state of the matrix `h` at temperature `T` (mK).

# Examples
```julia-repl
julia> gibbs_state(σz, 10)
2×2 Array{Complex{Float64},2}:
 0.178338+0.0im       0.0+0.0im
      0.0+0.0im  0.821662+0.0im
```
"""
function gibbs_state(h, T)
    β = temperature_2_β(T)
    Z = 0.0
    w, v = eigen(Hermitian(h))
    res = zeros(eltype(v), size(h))
    for (i, E) in enumerate(w)
        t = exp(-β * E)
        vi = @view v[:, i]
        res += t * vi * vi'
        Z += t
    end
    res / Z
end

"""
    low_level_matrix(M, lvl)

Calculate the matrix `M` projected to lower energy subspace containing `lvl` energy lvl.

# Examples
```julia-repl
julia> low_level_matrix(σx⊗σx, 2)
4×4 Array{Complex{Float64},2}:
 -0.5+0.0im   0.0+0.0im   0.0+0.0im   0.5+0.0im
  0.0+0.0im  -0.5+0.0im   0.5+0.0im   0.0+0.0im
  0.0+0.0im   0.5+0.0im  -0.5+0.0im   0.0+0.0im
  0.5+0.0im   0.0+0.0im   0.0+0.0im  -0.5+0.0im
```
"""
function low_level_matrix(M, lvl)
    if lvl > size(M, 1)
        @warn "Subspace dimension bigger than total dimension."
        return M
    else
        w, v = eigen(Hermitian(M))
        res = zeros(eltype(v), size(M))
        for i in range(1, stop = lvl)
            vi = @view v[:, i]
            res += w[i] * vi * vi'
        end
        return res
    end
end

"""
$(SIGNATURES)

Test if `𝐔` is a unitary matrix. The function checks how close both `` 𝐔𝐔^† `` and `` 𝐔^†𝐔 `` are to `` I ``, with relative and absolute error given by `rtol`, `atol`.

# Examples
```julia-repl
julia> check_unitary(exp(-1.0im*5*0.5*σx))
true
```
"""
function check_unitary(𝐔::AbstractMatrix; rtol = 1e-6, atol = 1e-8)
    a1 = isapprox(
        𝐔 * 𝐔',
        Matrix{eltype(𝐔)}(I, size(𝐔)),
        rtol = rtol,
        atol = atol,
    )
    a2 = isapprox(
        𝐔' * 𝐔,
        Matrix{eltype(𝐔)}(I, size(𝐔)),
        rtol = rtol,
        atol = atol,
    )
    a1 && a2
end

"""
$(SIGNATURES)

Generate a log-uniformly distributed array with `num` elements between `a` and `b`. The base of log is `base` with default value 10.
"""
function log_uniform(a, b, num; base = 10)
    loga = log(base, a)
    logb = log(base, b)
    10 .^ range(loga, logb, length = num)
end

"""
$(SIGNATURES)

Calculate the partial trace of the density matrix `ρ`, assuming all the subsystems are qubits. `qubit_2_keep` denotes the indices whose corresponding qubits are not traced out.

# Examples
```julia-repl
julia> ρ1 = [0.4 0.2; 0.2 0.6]; ρ2 = [0.5 0; 0 0.5];
julia> partial_trace(ρ1⊗ρ2, [1])
2×2 Array{Float64,2}:
 0.4  0.2
 0.2  0.6
```
"""
function partial_trace(ρ::Matrix, qubit_2_keep::AbstractVector{Int})
    num = Int(log2(size(ρ, 1)))
    partial_trace(ρ, [2 for i in 1:num], qubit_2_keep)
end

"""
$(SIGNATURES)

Calculate the partial trace of the density matrix `ρ`. Assume `ρ` is in the tensor space of ``ℋ₁⊗ℋ₂⊗…``, `sys_dim` is an array of the corresponding sub-system dimensions. `dim_2_keep` is an array of the indices whose corresponding subsystems are not traced out.

# Examples
```julia-repl
julia> ρ1 = [0.4 0.2; 0.2 0.6]; ρ2 = [0.5 0; 0 0.5];
julia> partial_trace(ρ1⊗ρ2, [2, 2], [1])
2×2 Array{Float64,2}:
 0.4  0.2
 0.2  0.6
```
"""
function partial_trace(ρ::Matrix, sys_dim::AbstractVector{Int}, dim_2_keep::AbstractVector{Int})
    size(ρ, 1) != size(ρ, 2) && throw(ArgumentError("ρ is not a square matrix."))
    prod(sys_dim) != size(ρ, 1) && throw(ArgumentError("System dimensions do not multiply to density matrix dimension."))

    N = length(sys_dim)
    sys_dim = reverse(sys_dim)
    ρ = reshape(ρ, sys_dim..., sys_dim...)
    dim_2_keep = N .+ 1 .- dim_2_keep
    dim_2_trace = [i for i in 1:N if !(i in dim_2_keep)]
    traction_idx = collect(1:2*N)
    traction_idx[dim_2_trace .+ N] .= traction_idx[dim_2_trace]
    tensortrace(ρ, traction_idx)
end

"""
$(SIGNATURES)

Calculate the indices of the degenerate values in list `w`. `w` must be sorted in ascending order.

# Examples
```julia-repl
julia> w = [1,1,1,2,3,4,4,5]
julia> find_degenerate(w)
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [6, 7]
```
"""
function find_degenerate(w; atol::Real=1e-6, rtol::Real=0)
    dw = w[2:end] - w[1:end-1]
    inds = findall((x)->isapprox(x, 0.0, atol=atol, rtol=rtol), dw)
    res = Array{Array{Int, 1}, 1}()
    temp = Int[]
    for i in eachindex(inds)
        push!(temp, inds[i])
        if (i==length(inds)) || !(inds[i]+1 == inds[i+1])
            push!(temp, inds[i]+1)
            push!(res, temp)
            temp = Int[]
        end
    end
    res
end

"""
$(SIGNATURES)

Calculate the fidelity between two density matrices `ρ` and `σ`: ``Tr[\\sqrt{\\sqrt{ρ}σ\\sqrt{ρ}}]²``.

# Examples
```julia-repl
julia> ρ = PauliVec[1][1]*PauliVec[1][1]'
julia> σ = PauliVec[3][1]*PauliVec[3][1]'
julia> fidelity(ρ, σ)
0.49999999999999944
```
"""
function fidelity(ρ, σ)
    if !check_density_matrix(ρ) || !check_density_matrix(σ)
        throw(ArgumentError("Inputs are not valid density matrices."))
    end
    temp = sqrt(ρ)
    real(tr(sqrt(temp * σ * temp))^2)
end

"""
$(SIGNATURES)

Check whether matrix `ρ` is a valid density matrix. `atol` and `rtol` is the absolute and relative error tolerance to use when checking `ρ` ≈ 1.

# Examples
```julia-repl
julia> check_density_matrix(σx)
false
julia> check_density_matrix(σi/2)
true
```
"""
function check_density_matrix(ρ; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps())
    ishermitian(ρ) && (eigmin(ρ) >= 0) && isapprox(tr(ρ), 1, atol=atol, rtol=rtol)
end

"""
$(SIGNATURES)

Calculate the purity of density matrix `ρ`: ``tr(ρ²)``. This function does not check whether `ρ` is a valid density matrix.
"""
function purity(ρ)
    tr(ρ*ρ)
end

"""
$(SIGNATURES)

Check whether the input `ρ` is a pure state: ``tr(ρ²)≈1``. This function will first check if `ρ` is a valid density matrix. `atol` and `rtol` are the corresponding absolute and relative error tolerance used for float point number comparison.
"""
function check_pure_state(ρ; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps())
    check_density_matrix(ρ) && isapprox(purity(ρ), 1, atol=atol, rtol=rtol)
end


text:
# ----------------------------------------------------
# --- Metadata
# ----------------------------------------------------
mutable struct uhd_rx_metadata_t end

mutable struct uhd_tx_metadata_t end

mutable struct uhd_async_metadata_t end

const uhd_rx_metadata_handle = Ptr{uhd_rx_metadata_t}

const uhd_tx_metadata_handle = Ptr{uhd_tx_metadata_t}

const uhd_async_metadata_handle = Ptr{uhd_async_metadata_t}

# ----------------------------------------------------
# --- Ranges
# ----------------------------------------------------
struct uhd_range_t
    start::Cdouble
    stop::Cdouble
    step::Cdouble
end
mutable struct uhd_meta_range_t end
const uhd_meta_range_handle = Ptr{uhd_meta_range_t}

mutable struct uhd_sensor_value_t end
const uhd_sensor_value_handle = Ptr{uhd_sensor_value_t}

mutable struct uhd_string_vector_t end
const uhd_string_vector_handle = Ptr{uhd_string_vector_t}

# ----------------------------------------------------
# --- Tune request (board managment)
# ----------------------------------------------------
struct uhd_tune_request_t
    target_freq::Cdouble
    rf_freq_policy::uhd_tune_request_policy_t
    rf_freq::Cdouble
    dsp_freq_policy::uhd_tune_request_policy_t
    dsp_freq::Cdouble
    args::Ptr{Cchar}
end

struct uhd_tune_result_t
    clipped_rf_freq::Cdouble
    target_rf_freq::Cdouble
    actual_rf_freq::Cdouble
    target_dsp_freq::Cdouble
    actual_dsp_freq::Cdouble
end

struct uhd_usrp_rx_info_t
    mboard_id::Ptr{Cchar}
    mboard_name::Ptr{Cchar}
    mboard_serial::Ptr{Cchar}
    rx_id::Ptr{Cchar}
    rx_subdev_name::Ptr{Cchar}
    rx_subdev_spec::Ptr{Cchar}
    rx_serial::Ptr{Cchar}
    rx_antenna::Ptr{Cchar}
end

struct uhd_usrp_tx_info_t
    mboard_id::Ptr{Cchar}
    mboard_name::Ptr{Cchar}
    mboard_serial::Ptr{Cchar}
    tx_id::Ptr{Cchar}
    tx_subdev_name::Ptr{Cchar}
    tx_subdev_spec::Ptr{Cchar}
    tx_serial::Ptr{Cchar}
    tx_antenna::Ptr{Cchar}
end

# ----------------------------------------------------
# --- EEPROM
# ----------------------------------------------------
mutable struct uhd_dboard_eeprom_t end
const uhd_dboard_eeprom_handle = Ptr{uhd_dboard_eeprom_t}

mutable struct uhd_mboard_eeprom_t end
const uhd_mboard_eeprom_handle = Ptr{uhd_mboard_eeprom_t}

# ----------------------------------------------------
# --- Subdev
# ----------------------------------------------------
struct uhd_subdev_spec_pair_t
    db_name::Ptr{Cchar}
    sd_name::Ptr{Cchar}
end

mutable struct uhd_subdev_spec_t end
const uhd_subdev_spec_handle = Ptr{uhd_subdev_spec_t}

struct uhd_usrp_register_info_t
    bitwidth::Csize_t
    readable::Bool
    writable::Bool
end

# ----------------------------------------------------
# --- Streamers
# ----------------------------------------------------
struct uhd_stream_args_t
    cpu_format::Ptr{Cchar}
    otw_format::Ptr{Cchar}
    args::Ptr{Cchar}
    channel_list::Ptr{Csize_t}
    n_channels::Cint
end

struct uhd_stream_cmd_t
    stream_mode::uhd_stream_mode_t
    num_samps::Csize_t
    stream_now::Bool
    time_spec_full_secs::Int64
    time_spec_frac_secs::Cdouble
end

mutable struct uhd_rx_streamer end
const uhd_rx_streamer_handle = Ptr{uhd_rx_streamer}

mutable struct uhd_tx_streamer end
const uhd_tx_streamer_handle = Ptr{uhd_tx_streamer}

# ----------------------------------------------------
# --- Core USRP
# ----------------------------------------------------
mutable struct uhd_usrp end
const uhd_usrp_handle = Ptr{uhd_usrp}

mutable struct uhd_usrp_clock end
const uhd_usrp_clock_handle = Ptr{uhd_usrp_clock}

max_stars_repo_name: RGerzaguet/UHDBindings

max_issues_repo_name: RGerzaguet/UHDBindings


text:
# Show how  google drive can be used for file upload
# Make an account for the class
# Suggest making a  course based google  account
# Make a form (forms.google.com)
# Some useful form details can  be pre-filled.  You have access to (at  least):
# * studentName
# * studentLogin
# * studentID
# * setNumber
# * probNum
#
# Replace the pre-filled value in the URL with a PERL reference to a variable  from above. E.g.:
# https://docs.google.com/forms/d/e/1FAIpQLScZXmLG0N4GmbLn1GXMbFkoOXVQb2d37NQug9DkFte48fiEbw/viewform?usp=pp_url&entry.941198495=$studentName&entry.52654376=$setNumber&entry.2147005713=$probNumber
#
# Include a link using  [tag](link) markdown from a question, or a  label, as below

# Variables that

q = label(raw"""

Upload your file: [upload file](https://docs.google.com/forms/d/e/1FAIpQLScZXmLG0N4GmbLn1GXMbFkoOXVQb2d37NQug9DkFte48fiEbw/viewform?usp=pp_url&entry.941198495=$studentName&entry.52654376=$setNumber&entry.2147005713=$probNum)

""")

Page("File upload example", (q,))


text:
using Documenter

try
    using SimpleHypergraphs
catch
    if !("../src/" in LOAD_PATH)
	   push!(LOAD_PATH,"../src/")
	   @info "Added \"../src/\"to the path: $LOAD_PATH "
	   using SimpleHypergraphs
    end
end

makedocs(
    sitename = "SimpleHypergraphs",
    format = Documenter.HTML(),
    modules = [SimpleHypergraphs],
	pages = ["index.md", "reference.md"],
	doctest = true
)

deploydocs(
    repo ="github.com/pszufe/SimpleHypergraphs.jl.git",
	target="build"
)


text:
using Documenter, Collatz
# https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#Overview
# https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
const _devurl = "dev"
@time deploydocs(
    repo = "github.com/Skenvy/Collatz.git",
    deploy_config = Documenter.GitHubActions(),
    # A special branch for hosting the built docs
    branch = "gh-pages-julia",
    dirname = "julia",
    # Where to deploy to from a non-tagged commit?
    devurl = _devurl,
    devbranch = "main",
    versions = ["stable" => "v^", "v#.#", _devurl => _devurl],
    push_preview = true,
)

max_stars_repo_name: Skenvy/Collatz

max_issues_repo_name: Skenvy/Collatz

max_forks_repo_name: Skenvy/Collatz


text:
#  Thu Dec 12 09:49:16 CET 2019

#
# Simulation: cluster formation for cells with polarized movements
#

# auxiliary codes
include("aux.jl")
using Random
#
# PARAMETERS
#

# NOTE:: add a function to collect command line parameters
const _L = [64,64]
const _steps   = 30000
const _samples = 100
const _icoordination = 4
_rates   = Dict()
_rates["depolarization"]= 0.05e0
_rates["polarization"]  = 0.1e0  # total polariz rate ( 4 times larger )
_rates["displacement"]  = 1e0

#
# DATA INIT
#
iskip=1
data = zeros(Float64,Int(_steps/iskip),2)
L_tot = prod(_L)
#
# MONTE CARLO SAMPLES
#
for sample=1:_samples
    n0 = 32 #Int(prod(_L)*0.1)
    cell,lattice = init(_L,n0);
    idata_idx = 0
    for step=0:_steps-1
        N = size(cell)[1]
        #
        # data collection
        #

        if (mod(step,iskip) == 0)
            idata_idx += 1
            x = sum( lattice .== _nonpolarized)
            y = sum( lattice .>  _nonpolarized)
            data[idata_idx,1] += (x)/_samples
            data[idata_idx,2] += (y)/_samples
        end

        # shuffle the cell vector
        # shuffle!(cell);
        # update
        iflag = 0
        k = 1
        prob=0e0
        rng = rand() # use MersenneTwister
        while ((iflag < 1) & (k < N+1))
            iflag,prob = update(k,cell,lattice,_rates,prob,rng)
            k += 1
        end
    end
end


text:
#=
dictionary:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module Dictionary
    # dict
    d = Dict{String, Int}()

    # methods
    function add_key_val(key::String, val::Int)
        d[key] = val
        return d
    end

    function get_by_key(key::String)
        return d[key]
    end

    function iterate_for(d::Dict)
        for (key, val) in d
            println(key, val)
        end
    end

    function main()
        println("Add 'a', 1: ", add_key_val("a", 1))
        println("Add 'b', 2: ", add_key_val("b", 2))
        println("Add 'c', 3: ", add_key_val("c", 3))
        println("Get val of 'a': ", get_by_key("a"))
        println("Has key, 'b': ", haskey(d, "b"))
        iterate_for(d)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Dictionary
    Dictionary.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Test
using AdventOfCode.TestUtils

TestUtils.init()

@testset "AdventOfCode" begin
    include("day01.jl")
    include("day02.jl")
    include("day03.jl")
    include("day04.jl")
    include("day05.jl")
    include("day06.jl")
    include("day07.jl")
    include("day08.jl")
    include("day09.jl")
end

TestUtils.output("perfs.md")


text:
import UUIDs: UUID, uuid1
import .ExpressionExplorer: SymbolsState

"The building block of a `Notebook`. Contains code, output, reactivity data, mitochondria and ribosomes."
mutable struct Cell
    "Because Cells can be reordered, they get a UUID. The JavaScript frontend indexes cells using the UUID."
    cell_id::UUID
    code::String

    output_repr::Union{Vector{UInt8},String,Nothing}
    repr_mime::MIME
    errored::Bool
    runtime::Union{Missing,UInt64}
    code_folded::Bool
    queued::Bool
    running::Bool

    parsedcode::Union{Nothing,Expr}
    module_usings::Set{Expr}
    rootassignee::Union{Nothing,Symbol}
end

Cell(cell_id, code) = Cell(cell_id, code, nothing, MIME("text/plain"), false, missing, false, false, false, nothing, Set{Expr}(), nothing)
Cell(code) = Cell(uuid1(), code)
Cell() = Cell("")


text:
"""
Mean mutual information of column a (Dunn et. al. 2008).
Summation is over j=1 to N, j ≠ a. Total is N-1.
"""
function _mean_column(mi::Matrix{T}) where {T}
    (dropdims(sum(mi, dims=1), dims=1) .- diag(mi)) ./ (size(mi, 1)-1)
end

"""
Mean mutual information of column a (Dunn et. al. 2008).
Summation is over j=1 to N, j ≠ a. Total is N-1.

Overall mean mutual information (Dunn et. al. 2008).
2/(N*(N-1)) by the sum of MI where the indices run i=1 to N-1, j=i+1 to N (triu).
"""
function _mean_total(mi::Matrix{T}) where T
    values = matrix2list(mi)
    sum(values) / length(values)
end

"APC (Dunn et. al. 2008)"
function APC!(MI::Matrix{T}) where T
    nrow, ncol = size(MI)
    MI_mean = _mean_total(MI)
    #   if MI_mean ==  0.0
    #     return(fill!(MI, zero(T)))
    #   end
    MI_res_mean = _mean_column(MI)
    @inbounds for j in 1:ncol
        for i in 1:nrow
            if i != j
                MI[i,j] -= ((MI_res_mean[i] * MI_res_mean[j]) /  MI_mean)
            else
                MI[i,j] = NaN
            end
        end
    end
    MI
end

function APC!(MI::PairwiseListMatrix{T,true,Vector{T}}) where T <: AbstractFloat
    ncol = MI.nelements
    MI_col_sum = vec(sum_nodiag(MI, dims=1))
    MI_mean = sum(MI_col_sum) / (length(MI) - ncol)
    MI_col_mean = MI_col_sum ./ (ncol-1)
    k = 0
    list = MI.list
    @inbounds for i in 1:ncol
        for j in i:ncol
            k += 1
            if i != j
                list[k] -= ((MI_col_mean[i] * MI_col_mean[j]) /  MI_mean)
            else
                list[k] = NaN
            end
        end
    end
    MI
end

function APC!(MI::PairwiseListMatrix{T,false,Vector{T}}) where T <: AbstractFloat
    ncol = MI.nelements
    MI_col_sum = vec(sum_nodiag(MI, dims=1))
    MI_mean = sum(MI_col_sum) / (length(MI) - ncol)
    MI_col_mean = MI_col_sum ./ (ncol-1)
    k = 0
    list = MI.list
    @inbounds for i in 1:(ncol-1)
        for j in (i+1):ncol
            k += 1
            list[k] -= ((MI_col_mean[i] * MI_col_mean[j]) /  MI_mean)
        end
    end
    MI.diag[:] .= NaN
    MI
end

function APC!(MI::NamedArray{T,2,PairwiseListMatrix{T,D,Vector{T}},NTuple{2,OrderedDict{String,Int}}}) where {T,D}
    plm = getarray(MI)::PairwiseListMatrix{T,D,Vector{T}}
    APC!(plm)
    MI
end