text:
"""
    randn([rng=GLOBAL_RNG], T=Quaternion{Float64}, [dims...])

Generate a normally distributed random quaternion of type `T` with mean 0 and standard
deviation of norm 1.  Optionally generate an *array* of such quaternions.  This module
currently provides an implementation for the types `QuaternionF16`, `QuaternionF32`, and
`QuaternionF64` (the default).  The values are drawn from the spherically symmetric
quaternionic normal distribution of variance 1 (corresponding to each component having
independent normal distribution with mean zero and variance 1/4).

Note that this function works with `Quaternion{BigFloat}`, even though `Base.randn` does
not work with `BigFloat`; we just use the [Box-Muller
transform](https://en.wikipedia.org/wiki/Box–Muller_transform) to obtain the desired
result.

If the quaternion type passed in is `Rotor`, the result will be normalized correctly.
Because the distribution is spherically symmetric, the result is a truly random
rotation.

If the quaternion type is `QuatVec`, the result will have a 0 scalar component, and the
vector will have mean 0 standard deviation of norm 1.

# Examples
```julia
julia> randn(QuaternionF64)
0.4336736009756228 - 0.45087190792840853𝐢 - 0.24723937675211696𝐣 - 0.4514571469326208𝐤
julia> randn(QuaternionF16, 2, 2)
2×2 Matrix{QuaternionF16}:
   0.4321 + 1.105𝐢 + 0.2664𝐣 - 0.1359𝐤   0.064 + 0.9263𝐢 - 0.4138𝐣 + 0.05505𝐤
 0.2512 - 0.2585𝐢 - 0.2803𝐣 - 0.00964𝐤  -0.1256 + 0.1848𝐢 + 0.03607𝐣 - 0.752𝐤
```

"""
Base.randn(rng::AbstractRNG, QT::Type{<:AbstractQuaternion{T}}) where {T<:AbstractFloat} =
    QT(randn(rng, T)/2, randn(rng, T)/2, randn(rng, T)/2, randn(rng, T)/2)

function Base.randn(rng::AbstractRNG, QT::Type{<:AbstractQuaternion{BigFloat}})
    # Use the Box-Muller transform to get randn BigFloats from rand BigFloat
    c = rand(rng, BigFloat, 4)
    l1 = √(-log(c[1])/2)
    s2, c2 = sincospi(2*c[2])
    l3 = √(-log(c[3])/2)
    s4, c4 = sincospi(2*c[4])
    QT(l1*c2, l1*s2, l3*c4, l3*s4)
end

_q3v_factor(::Type{T}) where T = inv(√T(3))

Base.randn(rng::AbstractRNG, QT::Type{QuatVec{T}}) where {T<:AbstractFloat} =
    QT(0, randn(rng, T)*_q3v_factor(T), randn(rng, T)*_q3v_factor(T), randn(rng, T)*_q3v_factor(T))

function Base.randn(rng::AbstractRNG, QT::Type{QuatVec{BigFloat}})
    # Use the Box-Muller transform to get randn BigFloats from rand BigFloat
    c = rand(rng, BigFloat, 4)
    l1 = √(-log(c[1])/2)
    s2 = sinpi(2*c[2])
    l3 = √(-log(c[3])/2)
    s4, c4 = sincospi(2*c[4])
    QuatVec(0, l1*s2, l3*c4, l3*s4)
end


text:
using ModelingToolkit
using Test
MT = ModelingToolkit

@variables t x
struct MyNLS <: MT.AbstractSystem
    name
    systems
end
@test_logs (:warn,) tmp = independent_variables(MyNLS("sys", []))
tmp = independent_variables(MyNLS("sys", []))
@test tmp == []

struct MyTDS <:MT.AbstractSystem
    iv
    name
    systems
end
@test_logs (:warn,) iv = independent_variables(MyTDS(t, "sys", []))
iv = independent_variables(MyTDS(t, "sys", []))
@test all(isequal.(iv, [t]))

struct MyMVS <:MT.AbstractSystem
    ivs
    name
    systems
end
@test_logs (:warn,) ivs = independent_variables(MyMVS([t, x], "sys", []))
ivs = independent_variables(MyMVS([t, x], "sys", []))
@test all(isequal.(ivs, [t, x]))

ext: jl

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
# Test functions for parsing
function fixable1(x)
    return x
end
fixable2(x) = x
function fixable3(A)
    s = 0
    fi = firstindex(A)
    for i in eachindex(A)
        for j in fi:i-1
            s += A[j]
        end
    end
    return s
end
function unfixable1(A)
    s = 0
    fi = firstindex(A)
    for i in eachindex(A)
        for j in fi:i-1
            s += A[j]
        end end
    return s
end

# Generated functions
@generated function generated1(A::AbstractArray{T,N}, val) where {T,N}
    ex = Expr(:tuple)
    for i = 1:N
        push!(ex.args, :val)
    end
    return ex
end
call_generated1(ndims) = generated1(fill(0, ntuple(d->1, ndims)...), 7)

# getproperty is defined in sysimg.jl
getline(lnn) = lnn.line

function f63()
    x = 1 + 1
    @info "hello"
    y = 7
end


text:
using DataDeps

"""
    extract_targz(file)

Platform-independent file extraction
"""
function extract_targz(file)
  if Sys.iswindows()
    run(pipeline(`7z x -tgzip -so $file`, `7z x -si -ttar`))
  else
    run(`tar -xzf $file`)
  end
end

function data_folder_moist_thermo()
  register(DataDep("MoistThermoConstructorTestData",
                   "Data to test moist thermodynamic constructors",
                   "https://caltech.box.com/shared/static/0741fg18sav94jyt94j4lyymiao4eobw.gz",
                   "d826e1cf57ba8bfcdccf2fb2b3f264c508e7929032a2b6639b379a8d1b1d5dfd",
                   post_fetch_method=extract_targz))
  datafolder = datadep"MoistThermoConstructorTestData"
  return datafolder
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using Pkg
using CSV
using Statistics
using OnlinePCA
using Distributions
using DelimitedFiles
using Arpack
using HypergeometricFunctions
using Distances
using StatsBase

mat = CSV.read("/home/cbl06/Documents/Matrix.csv")
genes = mat[1]
cells = mat[1,:]
@time arr = convert(Matrix{Int32}, mat[2:2169])
arrlog = Float32.(log.(arr .+ 1))
X = arr
bool = var(X, dims = 2)[:] .!= 0
X = X[bool, :]
@time standardize(ZScoreTransform, Float32.(X), dims=2)
rmin = minimum(X, dims = 2)
rmax = maximum(X, dims = 2)
rangeMM = rmax - rmin
X = (X .- rmin) ./ rangeMM
StatsBase.standardize(UnitRangeTransform, [0.0 -0.5 0.5; 0.0 1.0 2.0], dims=2)
standardize(UnitRangeTransform, Array{Float64,2}(X), dims =1)
X = convert(Matrix{Float32}, X)
X = vcat(X, 1 .- X)
Acol = sum(X, dims = 1)
Arow = sum(X, dims = 2)
Acol_sqrt = convert(Matrix{Float32}, sqrt.(Acol))
Arow_sqrt = convert(Matrix{Float32}, sqrt.(Arow))
@time Z = X ./ Acol_sqrt ./ Arow_sqrt

U, S, V = (Arpack.svds(Z, nsv = 50))[1]
gene = (Z * V)[:, 1:end .!= 1]
cell = (size(Z)[2] * V)[:, 1:end .!= 1]
@time dist = Distances.pairwise(Euclidean(), cell, gene, dims = 1)

P = mapslices(x->sortperm(x; rev = false, alg = PartialQuickSort(200)), dist, dims = 1)
mapslices(sortperm, dist, dims = 2)

sil = sortperm(dist[:,1], rev = false)[1:200]
sortperm(eachcol(dist[:, 1:5]))
dist[sil,1]

sortperm(dist[:,32], rev = true, alg = PartialQuickSort(200))

function NormalizeData(X::AbstractArray; scale_factor::Int = 10000)::Array{Float32,2}
    return log1p.(scale_factor * (X ./ sum(X, dims = 1)))
end
@time NormalizeData(arr)
scale(arr, d::MvLogNormal)
StatsBase.standardize(ZScoreTransform, [0.0 -0.5 5.5; 0.0 1.0 2.0], dims=2)

JuliaSeurat.NormalizeData(arrlog)

max_stars_repo_name: Cortalak/Pissaro

max_issues_repo_name: Cortalak/Pissaro

max_forks_repo_name: Cortalak/Pissaro


text:
module Day07

using Statistics

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day07"), String)

export part1, part2

"""
The total fuel required to align all crab submarines to `target`
This function accepts an optional `fuel_cost_function` which calculates the fuel
required to move a crab submarine from one location to another. This function
should act on non-negative integer values (ie. absolute difference). The default
is the `identity` function, in which the distance between the two locations is
also the fuel cost.
"""
function fuel_to_align(
    positions::Vector{Int64},
    target::Integer,
    fuel_cost_function::Function = identity
)
    absolute_movements = abs.(positions .- target)
    fuel_costs = fuel_cost_function.(absolute_movements)
    sum(fuel_costs)
end

"""
I'm calling this `brute_force` since it inelegantly searches the entire solution
space for a minimum cost. If this function isn't good enough for part 2 I can
look for leaner search strategies (starting from median position maybe?).

Starting with the assumption that the optimal position must be somewhere between
the minimum and maximum positions, we can check each possible position between
the two to determine the minimum fuel consumption.
"""
function brute_force(positions::Vector{Int64}, fuel_cost_function::Function = identity)
    minimum_fuel = Inf
    for position = minimum(positions):maximum(positions)
        required_fuel = fuel_to_align(positions, position, fuel_cost_function)
        minimum_fuel = min(required_fuel, minimum_fuel)
    end
    Integer(minimum_fuel)
end

function part1(input = input)
    positions = [parse(Int64, position) for position in split(input, ",")]
    brute_force(positions)
end

"""
I think this is an appropriate name? We're summing an arithmetic sequence.
"""
arithmetic_fuel_cost(movement::Int64) = movement * (movement + 1) / 2

function part2(input = input)
    positions = [parse(Int64, position) for position in split(input, ",")]
    brute_force(positions, arithmetic_fuel_cost)
end

end # module


text:
using Documenter, WaterWaves

makedocs(;
    modules=[WaterWaves],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/TheoGerostathis/WaterWaves.jl/blob/{commit}{path}#L{line}",
    sitename="WaterWaves.jl",
    authors="Theodoros Gerostathis",
    assets=String[],
)

deploydocs(;
    repo="github.com/TheoGerostathis/WaterWaves.jl",
)


text:
using ITensors,
      Test

@testset "Site Types" begin

  N = 10

  @testset "Star in operator strings" begin
    sites = spinHalfSites(N)
    @test_throws ArgumentError op(sites, "Sp", 1)
    @test sites[1] isa Index
    Sz = op(sites,"Sz",2)
    SzSz = op(sites,"Sz * Sz",2)
    @test SzSz ≈ multSiteOps(Sz,Sz)
    Sy = op(sites,"Sy",2)
    SySy = op(sites,"Sy * Sy",2)
    @test SySy ≈ multSiteOps(Sy,Sy)

    sites = spinOneSites(N)
    @test_throws ArgumentError op(sites, "Sp", 1)
    Sz = op(sites,"Sz",2)
    SzSz = op(sites,"Sz * Sz",2)
    @test SzSz ≈ multSiteOps(Sz,Sz)
    Sy = op(sites,"Sy",2)
    SySy = op(sites,"Sy * Sy",2)
    @test SySy ≈ multSiteOps(Sy,Sy)
  end
end


text: println("good")


text:
export fit_tree!, fit_tree

using ...JLBoostTrees: is_left_child
using Tables

function viewtable(table, rows, cols)
    t = Tables.columns(table)
    v = map(cols) do c
      col = Tables.getcolumn(t, c)
      c => view(col, rows)
    end
    (; v...)
  end

function tree_diag_print(jlt)
    if jlt.parent===nothing
        parent_feature=nothing
        parent_split=nothing
        split_sign = ""
    else
        parent_feature=jlt.parent.splitfeature
        parent_split=jlt.parent.split
        split_sign = is_left_child(jlt) ? "<=" : ">"
    end
    jlt.splitfeature, jlt.split, parent_feature, split_sign, parent_split
end

"""
	_fit_tree(loss, df, target, features, warm_start, jlt, node_colsample_strategy)

Fit a tree by following a algorithm

Parameters:
* loss
    A SupervisedLoss from LossFunctions where `deriv` and `deriv2` are defined
* tbl
    A Tables.jl compatible table
* target
    The column name. Likely to be a `Symbol` or a `String`
* features
    The iterable of features name. The type of the elements should be Symbol or String
* warm_start
    A vector weight to serve as starting weights
* jlt
    The JLBoost tree to update
* stopping_criterion: Function
* tree_growth: Function
    A function to control where to grow the tree
* lambda = 0
    The L1 Norm regulization constnat
* gamma = 0
    The L2 Norm regulization constnat
* colsample_bynode = 1 (NOT IMPLEMENTED YET)
    What proportion of features to sample for each node
* colsample_bylevel = 1 (NOT IMPLEMENTED YET)
    What proportion of features to sample for each level
"""
function _fit_tree!(loss, tbl, target, features, warm_start,
    jlt::AbstractJLBoostTree = JLBoostTree(0.0),
    tree_growth = depth_wise,
    stopping_criterion = max_depth_stopping_criterion(6);
    lambda = 0, gamma = 0,
    verbose = false, #colsample_bynode = 1, colsample_bylevel = 1,
	kwargs...)

    @assert Tables.istable(tbl)

    tblc = Tables.columns(tbl)
    dfo = DataFrame(tbl)

    @assert nrow(tblc) >= 2 # seriously? you have so few records

	# make absolutely sure that target is not part of features
    features = setdiff(features, [target])

    if verbose
        println("`_fit_tree!`: Current state of tree $jlt")
    end

    no_more_gains_to_found = false

    while !no_more_gains_to_found && !stopping_criterion(jlt)
        if verbose
            println("Tree Depth: $(treedepth(jlt))")
        end
        # at the beginning there is only one leaf node which is the parent for all nodes
        # amongst the end nodes compute the best split and choose the best split based on the leaf
        # nodes keep only those where a split has not been decided
        # println("the nodes considered for expansion are $(get_leaf_nodes(jlt) |> length)")

        # println.(get_leaf_nodes(jlt) .|> tree_diag_print)

        leaf_nodes = filter(x->ismissing(x.splitfeature) && ismissing(x.gain), get_leaf_nodes(jlt))

        # println("the nodes considered for expansion are $(length(leaf_nodes))")
        # println.(leaf_nodes .|> tree_diag_print)
        # keep track of the best gains at each node as we do not want to store the gain in the tree

        best_split_dict = Dict()

        # for all nodes eligible for splitting
        # compute the best split feature and the best split point
        # set the split point
        # leaf_node = leaf_nodes[1]
        for leaf_node in leaf_nodes
            if leaf_node.parent === nothing
                # if the node is the parent
                tblc_filtered = tblc
                warm_start_filtered = warm_start
            else
                keeprow = keeprow_vec(tbl, leaf_node)
                #tblc_filtered = view(tblc, keeprow, :)
                tblc_filtered = Tables.columns(DataFrame(view(dfo, keeprow, :)))
                #println("TYPE")
                #println(typeof(tblc_filtered))
                #print(keeprow)
                #println(size(keeprow))
                #println(typeof(tblc))
                #tblc_filtered = tblc[keeprow, :]
                warm_start_filtered = view(warm_start, keeprow)
                #warm_start_filtered = warm_start[keeprow]
                if sum(keeprow) <= 2
                    # no rows are kept, so move to next node
                    leaf_node.split = typemin(Float64)
                    leaf_node.splitfeature = Symbol("too few records for split")
                    leaf_node.gain = typemin(Float64)
                    #println("this branch has too few records $(leaf_node)")
                    continue
                else
                    # println("detective")
                    # println(extrema(tblc_filtered.SepalLength))
                end
            end

            # compute the gain for all splits for all features
            split_with_best_gain = find_best_split(loss, tblc_filtered, features[1], target, warm_start_filtered,
                                lambda, gamma; verbose=verbose, kwargs...)

            for feature in Iterators.drop(features, 1)

                feature_split =find_best_split(loss, tblc_filtered, feature, target, warm_start_filtered,
                                    lambda, gamma; verbose=verbose, kwargs...)

                if feature_split.gain > split_with_best_gain.gain
                    split_with_best_gain = feature_split
                end
            end

            # remember the split but do not set children
            if verbose
                println("found a best split at $(split_with_best_gain.feature) $(split_with_best_gain.split_at) $(split_with_best_gain.gain) $(split_with_best_gain.further_split) for $(leaf_node)")
            end
            best_split_dict[leaf_node] = split_with_best_gain
            # set the parent tree node
            leaf_node.split = split_with_best_gain.split_at
            leaf_node.splitfeature = split_with_best_gain.feature
            leaf_node.gain = split_with_best_gain.gain

            key = collect(keys(best_split_dict))[1]
            # println("Type is $(typeof(key))")
            # println(key)
            #println(typeof(key))
            # reset the best split
            split_with_best_gain = ()
        end

        # tree_growth phase
        # select the node to grow based on growth function
        # the tree_growth function will return the list of
        # nodes_to_split = tree_growth(jlt)
        nodes_to_split::Vector{<:AbstractJLBoostTree} = tree_growth(jlt)
        # println("Nodes to split")
        # println.(nodes_to_split)
        # println("Dict nodes")
        # println.(keys(best_split_dict))

        #@assert haskey(best_split_dict, nodes_to_split[1])
        # println(best_split_dict)
        no_more_gains_to_found = true
        for node_to_split in nodes_to_split

            # there needs to be positive gain then apply split to the tree
            # print(best_split_dict)
            # println(best_split_dict)
            # println("node to split is next line")
            # # println(typeof(node_to_split))

            # println("mehmehmeh")

            if !haskey(best_split_dict, node_to_split)
                if verbose
                    println("Won't split on $(node_to_split)")
                end
                continue
            else
                if verbose
                    println("Will split on $(node_to_split)")
                end
            end
            split_with_best_gain = best_split_dict[node_to_split]

            if split_with_best_gain.further_split && (split_with_best_gain.gain > 0)
                no_more_gains_to_found = false
                left_treenode = JLBoostTree(split_with_best_gain.lweight; parent = node_to_split)
                right_treenode = JLBoostTree(split_with_best_gain.rweight; parent = node_to_split)
                node_to_split.children = [left_treenode, right_treenode]
            end
        end
    end # end !stopping_criterion(jlt)
    jlt
end


text:
"""
绘制相关的功能
"""
module Drawers

using Plots
using ..Basics

export @savepng
export draw_lines, draw_lines!
export draw_polygon, draw_polygon!
export draw_points, draw_points!
export draw_figure

"""
绘制散点图
"""
function draw_points(pts::Array{Point2D, 1}; kw...)
    fig = plot()
    draw_points!(fig, pts; kw...)
    return fig
end

"""
绘制散点图
"""
function draw_points!(plt, pts::Array{Point2D, 1}; kw...)
    scatter!(plt, [pt.x for pt in pts], [pt.y for pt in pts];
    legend=false, kw...)
    return plt
end

"""绘制线"""
function draw_lines(lns::Array{Segment, 1}; lcidx=nothing, kw...)
    fig = plot()
    draw_lines!(fig, lns; lcidx=lcidx, kw...)
    return fig
end

"""绘制线"""
function draw_lines!(plt, lns::Array{Segment, 1}; lcidx=nothing, kw...)
    plot!(plt, legend=false, kw...)
    colors = [:red, :orange, :yellow, :green, :grays, :blue, :purple]
    custom_color = !isnothing(lcidx)
    for (idx, seg) in enumerate(lns)
        cidx = custom_color ? lcidx[idx] : idx
        plot!(plt, [seg.pt1.x, seg.pt2.x], [seg.pt1.y, seg.pt2.y],
        color=colors[(cidx-1)%7+1])
    end
end

"""绘制多边形"""
function draw_polygon(pls::Array{<:AbstractPolygon, 1}; pcidx=nothing, kw...)
    fig = plot()
    draw_polygon!(fig, pls; pcidx=pcidx, kw...)
    return fig
end

"""绘制多边形"""
function draw_polygon!(plt, pls::Array{<:AbstractPolygon, 1}; pcidx=nothing, kw...)
    plot!(plt, legend=false, kw...)
    colors = [:red, :orange, :yellow, :green, :grays, :blue, :purple]
    custom_color = !isnothing(pcidx)
    for (idx, ply) in enumerate(pls)
        cidx = custom_color ? pcidx[idx] : idx
        ptxs = [pt.x for pt in ply.vertex]
        push!(ptxs, ply.vertex[1].x)
        ptys = [pt.y for pt in ply.vertex]
        push!(ptys, ply.vertex[1].y)
        plot!(plt,
            ptxs, ptys,
            lw=0,
            fill=(0, colors[(cidx-1)%7+1]),
            fillalpha=0.5
        )
    end
end

"""
绘制各个成分都有的图形
"""
function draw_figure(
    pts::Union{Nothing, Array{Point2D, 1}},
    lns::Union{Nothing, Array{Segment, 1}},
    pls::Union{Nothing, Array{<:AbstractPolygon, 1}};
    ptext=nothing, lcidx=nothing, pcidx=nothing
    )
    plt = plot()
    isnothing(pts) ? nothing : draw_points!(plt, pts, text=ptext)
    isnothing(lns) ? nothing : draw_lines!(plt, lns, lcidx=lcidx)
    isnothing(pls) ? nothing : draw_polygon!(plt, pls, pcidx=pcidx)
    return plt
end

"""
保存图片
"""
macro savepng(plt, name)
    esc(quote
        $png($plt, $name)
    end)
end

end


text:
##########################
# @@@ CORE FUNCTIONS @@@ #
##########################

function f̃(p::Proposition, vals::Vector{Float64})
    return sum(map(x -> p.rel(x, p.a), vals))/length(vals)
end

function ==̃ₕ(x::Int64, y::Int64; horizon::Int64=4)::Float64
    # |x-y| ≥ h
    if abs(x-y) ≥ horizon
        return float(0)
    # |x-y| < h
    else
        return (horizon - abs(x-y))/horizon
    end
end

function <̃ₕ(x::Int64, y::Int64; horizon::Int64=4)::Float64
    # y-x > h
    if y-x > horizon
        return float(1)
    # y ≦ x
    elseif y ≤ x
        return float(0)
    # y-x ≦ h
    elseif y-x ≤ horizon
        return (y-x)/horizon
    end
end

function R̃ₓ(rel::AbstractIntervalRelation,
        w₁::NTuple{2, Int64},
        w₂::NTuple{2, Int64};
        horizon::Int64=1)::Float64
    x, y, w, z = w₁[1], w₁[2], w₂[1], w₂[2]
    # relations
    if eltype(rel) == :L
        return <̃ₕ(y, w; horizon)
    elseif eltype(rel) == :A
        return ==̃ₕ(y, w; horizon)
    elseif eltype(rel) == :O
        return ⊓(<̃ₕ(x, w; horizon), <̃ₕ(w, y; horizon), <̃ₕ(y, z; horizon))
    elseif eltype(rel) == :E
        return ⊓(<̃ₕ(x, w; horizon), ==̃ₕ(y, z; horizon))
    elseif eltype(rel) == :D
        return ⊓(<̃ₕ(x, w; horizon), <̃ₕ(z, y; horizon))
    elseif eltype(rel) == :B
        return ⊓(==̃ₕ(x, w; horizon), <̃ₕ(z, y; horizon))
    # inverse relations
    elseif eltype(rel) == :InvL
        return <̃ₕ(z, x; horizon)
    elseif eltype(rel) == :InvA
        return ==̃ₕ(z, x; horizon)
    elseif eltype(rel) == :InvO
        return ⊓(<̃ₕ(w, x; horizon), <̃ₕ(x, z; horizon), <̃ₕ(z, y; horizon))
    elseif eltype(rel) == :InvE
        return ⊓(<̃ₕ(w, x; horizon), ==̃ₕ(y, z; horizon))
    elseif eltype(rel) == :InvD
        return ⊓(<̃ₕ(w, x; horizon), <̃ₕ(y, z; horizon))
    elseif eltype(rel) == :InvB
        return ⊓(==̃ₕ(x, w; horizon), <̃ₕ(y, z; horizon))
    end
end

# global H = Dict{UInt, Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}}()

function mc(I::ModalInstance, φ::ClassificationRule; minh::Int=1, maxh::Int=4, horizon::Int=4, memo::Bool=false)
    if !memo
        return _mc(I,φ; horizon)
    else
        return _mc3(I,φ,H; minh=minh, maxh=maxh, horizon=horizon)
    end
end

function _mc2(I::ModalInstance,
        φ::ClassificationRule,
        H::Dict{UInt, Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}};
        minh::Int=1,
        maxh::Int=4,
        horizon::Int=4)::Float64

    # @show φ

    # println("#1")
    if haskey(H, hash(_antecedent(φ))) && haskey(H[hash(_antecedent(φ))], hash(I)) && haskey(H[hash(_antecedent(φ))][hash(I)], horizon) && haskey(H[hash(_antecedent(φ))][hash(I)][horizon], (1,2))
        # println("##1")
        return H[hash(_antecedent(φ))][hash(I)][horizon][(1,2)]
    end

    sub = _collect_sorted_subformulas(φ)
    # L = Dict{Tuple{Union{Node{AbstractOperator}, Node{AbstractProposition}, Node{AbstractHeytingAlgebraValue}},NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y]
    for (i,ψ) in enumerate(sub)

        # @show ψ

        # println("#2")
        # if haskey(H, hash(ψ)) && haskey(H[hash(ψ)],hash(I)) && haskey(H[hash(ψ)][hash(I)],horizon)
        #     # println("##2")
        #     continue
        # end

        # println("#3")
        if !haskey(H, hash(ψ))
            # println("##3")
            H[hash(ψ)] = Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}()
            H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
            for h in minh:maxh
                # println("###3")
                H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
            end
        else
            # println("#4")
            if !haskey(H[hash(ψ)], hash(I))
                # println("##4")
                H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
                for h in minh:maxh
                    # println("###4")
                    H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
                end
            end
        end

        if typeof(ψ.data) <: AbstractProposition
            # println("#prop")
            for (x, y) ∈ intervals
                v = f̃(ψ.data, I[1][ψ.data.A][x:y])
                for h in minh:maxh
                    H[hash(ψ)][hash(I)][h][(x,y)] = v
                end
            end

        elseif typeof(ψ.data) <: Float64
            # println("#float")
            for (x, y) ∈ intervals
                for h in minh:maxh
                    H[hash(ψ)][hash(I)][h][(x,y)] = ψ.data
                end
            end

        elseif typeof(ψ.data) == ExistentialIntervalRelation{:G}
            # println("#<G>")
            for h in minh:maxh
                # println("###<G>")
                S = Float64[]
                for op in [exIntRel(:L),exIntRel(:A),exIntRel(:O),exIntRel(:E),exIntRel(:D),exIntRel(:B),exIntRel(:InvL),exIntRel(:InvL),exIntRel(:InvA),exIntRel(:InvO),exIntRel(:InvE),exIntRel(:InvD),exIntRel(:InvB)]
                    # println("####<G>")
                    for (x, y) ∈ intervals
                        # println("#####<G>")
                        s = 0.0
                        for (w, z) ∈ intervals
                            # println("#######<G>")
                            # @show H[ψ.right][I][h][(w,z)]
                            # if !haskey(H, ψ.right)
                            #     v = _mc(I,ψ.right,H; minh, maxh, horizon)
                            # else
                            #     v = H[ψ.right][I][h][(w,z)]
                            # end
                            # if !haskey(H, ψ.right)
                            #     println("@@@@@ missing ψ.right = $(ψ.right) in H -- i = $i")
                            #     sf = _collect_sorted_subformulas(φ)
                            #     for ξ in sf
                            #         @show (i,ξ)
                            #     end
                            # end
                            s = ⊔(s, ⊓(R̃ₓ(op, (x, y), (w, z); horizon=h), H[hash(ψ.right)][hash(I)][h][(w,z)]))
                        end

                        push!(S, s)
                    end
                end
                # @show S
                # if !haskey(H, hash(ψ))
                #     println("@@@@@ missing ψ = $(ψ) in H -- i = $i")
                #     sf = _collect_sorted_subformulas(φ)
                #     for ξ in sf
                #         @show ξ
                #     end
                # end
                H[hash(ψ)][hash(I)][h][(1,2)] = ⊔(S...)
                # @show H[ψ][I][h][(1,2)]
            end
            # println("##<G>")
            return H[hash(ψ)][hash(I)][horizon][(1,2)]

        elseif typeof(ψ.data) == UniversalIntervalRelation{:G}
            println("#[G]")
            for h in minh:maxh
                S = Float64[]
                for op in [univIntRel(:L),univIntRel(:A),univIntRel(:O),univIntRel(:E),univIntRel(:D),univIntRel(:B),univIntRel(:InvL),univIntRel(:InvA),univIntRel(:InvO),univIntRel(:InvE),univIntRel(:InvD),univIntRel(:InvB)]
                    for (x, y) ∈ intervals
                        s = 1.0
                        for (w, z) ∈ intervals
                            s = ⊓(s, ↦(R̃ₓ(op, (x, y), (w, z); horizon=h), H[hash(ψ.right)][hash(I)][h][(w,z)]))
                        end
                        push!(S, s)
                    end
                end
                H[hash(ψ)][hash(I)][h][(1,2)] = ⊓(S...)
            end
            println("##[G]")
            return H[hash(ψ)][hash(I)][horizon][(1,2)]

        elseif typeof(ψ.data) <: AbstractExistentialIntervalRelation
            # println("#ex")
            for h in minh:maxh
                # println("###ex")
                for (x, y) ∈ intervals
                    # println("####ex")
                    s = 0.0
                    for (w, z) ∈ intervals
                        # println("######ex")
                        s = ⊔(s, ⊓(R̃ₓ(ψ.data, (x, y), (w, z); horizon=h), ⊓(H[hash(ψ.left)][hash(I)][h][(w,z)], H[hash(ψ.right)][hash(I)][h][(w,z)])))
                    end
                    # println("####ex")
                    # @show s
                    # @show typeof(s)
                    H[hash(ψ)][hash(I)][h][(x,y)] = s
                end
                # println("###ex")
            end
            # println("##ex")

        elseif typeof(ψ.data) <: AbstractUniversalIntervalRelation
            println("#univ")
            for h in minh:maxh
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals
                        # @show H[ψ.left][I][h][(w,z)]
                        # @show typeof(H[ψ.left][I][h][(w,z)])
                        # @show H[ψ.right][I][h][(w,z)]
                        # @show typeof(H[ψ.right][I][h][(w,z)])
                        # @show R̃ₓ(ψ.data, (x, y), (w, z); horizon=h)
                        # @show typeof(R̃ₓ(ψ.data, (x, y), (w, z); horizon=h))
                        s = ⊓(s, ↦(R̃ₓ(ψ.data, (x, y), (w, z); horizon=h), ↦(H[hash(ψ.left)][hash(I)][h][(w,z)], H[hash(ψ.right)][hash(I)][h][(w,z)])))
                    end
                    H[hash(ψ)][hash(I)][h][(x,y)] = s
                end
            end
            println("##univ")
        elseif typeof(ψ.data) <: BinaryRelation{:∧}
            # println("#∧")
            for h in minh:maxh
                # println("###∧")
                for (x, y) in intervals
                    # println("####∧")
                    H[hash(ψ)][hash(I)][h][(x,y)] = ⊓(H[hash(ψ.left)][hash(I)][h][(x,y)], H[hash(ψ.right)][hash(I)][h][(x,y)])
                end
            end
            # println("##∧")
        elseif typeof(ψ.data) <: BinaryRelation{:→}
            # println("#→")
            for h in minh:maxh
                # println("###→")
                for (x, y) in intervals
                    # println("####→")
                    # @show H[ψ.left][I][h][(x,y)]
                    # @show typeof(H[ψ.left][I][h][(x,y)])
                    # @show H[ψ.right][I][h][(x,y)]
                    # @show typeof(H[ψ.right][I][h][(x,y)])
                    # error(1)
                    H[hash(ψ)][hash(I)][h][(x,y)] = ↦(H[hash(ψ.left)][hash(I)][h][(x,y)], H[hash(ψ.right)][hash(I)][h][(x,y)])
                end
            end
            # println("##→")
        end
    end
end

function _mc3(I::ModalInstance,
    φ::ClassificationRule,
    H::Dict{Tuple{UInt, UInt, UInt, Tuple{UInt, UInt}}, Float64};
    minh::Int=1,
    maxh::Int=4,
    horizon::Int=4)::Float64

    # @show φ

    if haskey(H, (hash(_antecedent(φ)), hash(I), hash(horizon), (hash(1), hash(2))))
        return H[(hash(_antecedent(φ)), hash(I), hash(horizon), (hash(1), hash(2)))]
    end
    # # println("#1")
    # if haskey(H, hash(_antecedent(φ))) && haskey(H[hash(_antecedent(φ))], hash(I)) && haskey(H[hash(_antecedent(φ))][hash(I)], horizon) && haskey(H[hash(_antecedent(φ))][hash(I)][horizon], (1,2))
    #     # println("##1")
    #     return H[hash(_antecedent(φ))][hash(I)][horizon][(1,2)]
    # end

    sub = _collect_sorted_subformulas(φ)
    L = Dict{Tuple{UInt,NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y && y-x+1 ≤ log2(N)]
    for ψ in sub

        # # Subformula already checked?
        # if haskey(H, (hash(ψ), hash(I), hash(horizon), (hash(1), hash(2))))
        #     continue
        # end

        # @show ψ

        # println("#2")
        # if haskey(H, hash(ψ)) && haskey(H[hash(ψ)],hash(I)) && haskey(H[hash(ψ)][hash(I)],horizon)
        #     # println("##2")
        #     continue
        # end

        # # println("#3")
        # if !haskey(H, hash(ψ))
        #     # println("##3")
        #     H[hash(ψ)] = Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}()
        #     H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
        #     for h in minh:maxh
        #         # println("###3")
        #         H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
        #     end
        # else
        #     # println("#4")
        #     if !haskey(H[hash(ψ)], hash(I))
        #         # println("##4")
        #         H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
        #         for h in minh:maxh
        #             # println("###4")
        #             H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
        #         end
        #     end
        # end

        if typeof(ψ.data) <: AbstractProposition
            # println("#prop")
            for (x, y) ∈ intervals
                v = f̃(ψ.data, I[1][ψ.data.A][x:y])
                # for h in minh:maxh
                L[(hash(ψ), (x, y))] = v
                # end
            end

        elseif typeof(ψ.data) <: Float64
            # println("#float")
            for (x, y) ∈ intervals
                # for h in minh:maxh
                L[(hash(ψ), (x, y))] = ψ.data
                # end
            end

        elseif typeof(ψ.data) == ExistentialIntervalRelation{:G}
            # println("#<G>")
            # for h in minh:maxh
                # println("###<G>")
            S = Float64[]
            for op in [exIntRel(:L),exIntRel(:A),exIntRel(:O),exIntRel(:E),exIntRel(:D),exIntRel(:B),exIntRel(:InvL),exIntRel(:InvL),exIntRel(:InvA),exIntRel(:InvO),exIntRel(:InvE),exIntRel(:InvD),exIntRel(:InvB)]
                # println("####<G>")
                for (x, y) ∈ intervals
                    # println("#####<G>")
                    s = 0.0
                    for (w, z) ∈ intervals
                        # println("#######<G>")
                        # @show H[ψ.right][I][h][(w,z)]
                        # if !haskey(H, ψ.right)
                        #     v = _mc(I,ψ.right,H; minh, maxh, horizon)
                        # else
                        #     v = H[ψ.right][I][h][(w,z)]
                        # end
                        # if !haskey(H, ψ.right)
                        #     println("@@@@@ missing ψ.right = $(ψ.right) in H -- i = $i")
                        #     sf = _collect_sorted_subformulas(φ)
                        #     for ξ in sf
                        #         @show (i,ξ)
                        #     end
                        # end
                        s = ⊔(s, ⊓(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end

                    push!(S, s)
                end
            end
            # @show S
            # if !haskey(H, hash(ψ))
            #     println("@@@@@ missing ψ = $(ψ) in H -- i = $i")
            #     sf = _collect_sorted_subformulas(φ)
            #     for ξ in sf
            #         @show ξ
            #     end
            # end
            # H[hash(ψ)][hash(I)][horizon][(1,2)] = ⊔(S...)
            v = ⊔(S...)
            H[(hash(ψ), hash(I), hash(horizon), (hash(1), hash(2)))] = v
                # @show H[ψ][I][h][(1,2)]
            # end
            # println("##<G>")
            return v

        elseif typeof(ψ.data) == UniversalIntervalRelation{:G}
            # println("#[G]")
            # for h in minh:maxh
            S = Float64[]
            for op in [univIntRel(:L),univIntRel(:A),univIntRel(:O),univIntRel(:E),univIntRel(:D),univIntRel(:B),univIntRel(:InvL),univIntRel(:InvA),univIntRel(:InvO),univIntRel(:InvE),univIntRel(:InvD),univIntRel(:InvB)]
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals
                        s = ⊓(s, ↦(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end
                    push!(S, s)
                end
            end
            v = ⊓(S...)
            # H[hash(ψ)][hash(I)][horizon][(1,2)] = ⊓(S...)
            H[(hash(ψ), hash(I), hash(horizon), (hash(1), hash(2)))] = v
            # end
            # println("##[G]")
            return v

        elseif typeof(ψ.data) <: AbstractExistentialIntervalRelation
            # println("#ex")
            # for h in minh:maxh
                # println("###ex")
            for (x, y) ∈ intervals
                # println("####ex")
                s = 0.0
                for (w, z) ∈ intervals
                    # println("######ex")
                    s = ⊔(s, ⊓(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ⊓(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                # println("####ex")
                # @show s
                # @show typeof(s)
                L[(hash(ψ), (x, y))] = s
            end
                # println("###ex")
            # end
            # println("##ex")

        elseif typeof(ψ.data) <: AbstractUniversalIntervalRelation
            # println("#univ")
            # for h in minh:maxh
            for (x, y) ∈ intervals
                s = 1.0
                for (w, z) ∈ intervals
                    # @show H[ψ.left][I][h][(w,z)]
                    # @show typeof(H[ψ.left][I][h][(w,z)])
                    # @show H[ψ.right][I][h][(w,z)]
                    # @show typeof(H[ψ.right][I][h][(w,z)])
                    # @show R̃ₓ(ψ.data, (x, y), (w, z); horizon=h)
                    # @show typeof(R̃ₓ(ψ.data, (x, y), (w, z); horizon=h))
                    s = ⊓(s, ↦(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ↦(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                L[(hash(ψ), (x, y))] = s
            end
            # end
            # println("##univ")
        elseif typeof(ψ.data) <: BinaryRelation{:∧}
            # println("#∧")
            # for h in minh:maxh
                # println("###∧")
            for (x, y) in intervals
                # println("####∧")
                L[(hash(ψ), (x, y))] = ⊓(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
            # end
            # println("##∧")
        elseif typeof(ψ.data) <: BinaryRelation{:→}
            # println("#→")
            # for h in minh:maxh
                # println("###→")
            for (x, y) in intervals
                # println("####→")
                # @show H[ψ.left][I][h][(x,y)]
                # @show typeof(H[ψ.left][I][h][(x,y)])
                # @show H[ψ.right][I][h][(x,y)]
                # @show typeof(H[ψ.right][I][h][(x,y)])
                # error(1)
                L[(hash(ψ), (x, y))] = ↦(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
            # end
            # println("##→")
        end
    end
    return H[hash(_antecedent(φ))][hash(I)][horizon][(1,2)]
end

function _mc(I::ModalInstance,
        φ::ClassificationRule,
        H::Dict{UInt, Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}};
        minh::Int=1,
        maxh::Int=4,
        horizon::Int=4)::Float64

    # @show φ

    # println("#1")
    if haskey(H, hash(_antecedent(φ))) && haskey(H[hash(_antecedent(φ))], hash(I)) && haskey(H[hash(_antecedent(φ))][hash(I)], horizon) && haskey(H[hash(_antecedent(φ))][hash(I)][horizon], (1,2))
        # println("##1")
        return H[hash(_antecedent(φ))][hash(I)][horizon][(1,2)]
    end

    sub = _collect_sorted_subformulas(φ)
    L = Dict{Tuple{UInt,NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y]
    for (i,ψ) in enumerate(sub)

        # @show ψ

        # println("#2")
        # if haskey(H, hash(ψ)) && haskey(H[hash(ψ)],hash(I)) && haskey(H[hash(ψ)][hash(I)],horizon)
        #     # println("##2")
        #     continue
        # end

        # println("#3")
        if !haskey(H, hash(ψ))
            # println("##3")
            H[hash(ψ)] = Dict{UInt, Dict{Int, Dict{NTuple{2,Int}, Float64}}}()
            H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
            for h in minh:maxh
                # println("###3")
                H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
            end
        else
            # println("#4")
            if !haskey(H[hash(ψ)], hash(I))
                # println("##4")
                H[hash(ψ)][hash(I)] = Dict{Int, Dict{NTuple{2,Int}, Float64}}()
                for h in minh:maxh
                    # println("###4")
                    H[hash(ψ)][hash(I)][h] = Dict{NTuple{2,Int}, Float64}()
                end
            end
        end

        if typeof(ψ.data) <: AbstractProposition
            # println("#prop")
            for (x, y) ∈ intervals
                v = f̃(ψ.data, I[1][ψ.data.A][x:y])
                # for h in minh:maxh
                L[(hash(ψ), (x, y))] = v
                # end
            end

        elseif typeof(ψ.data) <: Float64
            # println("#float")
            for (x, y) ∈ intervals
                # for h in minh:maxh
                L[(hash(ψ), (x, y))] = ψ.data
                # end
            end

        elseif typeof(ψ.data) == ExistentialIntervalRelation{:G}
            # println("#<G>")
            # for h in minh:maxh
                # println("###<G>")
            S = Float64[]
            for op in [exIntRel(:L),exIntRel(:A),exIntRel(:O),exIntRel(:E),exIntRel(:D),exIntRel(:B),exIntRel(:InvL),exIntRel(:InvL),exIntRel(:InvA),exIntRel(:InvO),exIntRel(:InvE),exIntRel(:InvD),exIntRel(:InvB)]
                # println("####<G>")
                for (x, y) ∈ intervals
                    # println("#####<G>")
                    s = 0.0
                    for (w, z) ∈ intervals
                        # println("#######<G>")
                        # @show H[ψ.right][I][h][(w,z)]
                        # if !haskey(H, ψ.right)
                        #     v = _mc(I,ψ.right,H; minh, maxh, horizon)
                        # else
                        #     v = H[ψ.right][I][h][(w,z)]
                        # end
                        # if !haskey(H, ψ.right)
                        #     println("@@@@@ missing ψ.right = $(ψ.right) in H -- i = $i")
                        #     sf = _collect_sorted_subformulas(φ)
                        #     for ξ in sf
                        #         @show (i,ξ)
                        #     end
                        # end
                        s = ⊔(s, ⊓(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end

                    push!(S, s)
                end
            end
            # @show S
            # if !haskey(H, hash(ψ))
            #     println("@@@@@ missing ψ = $(ψ) in H -- i = $i")
            #     sf = _collect_sorted_subformulas(φ)
            #     for ξ in sf
            #         @show ξ
            #     end
            # end
            H[hash(ψ)][hash(I)][horizon][(1,2)] = ⊔(S...)
                # @show H[ψ][I][h][(1,2)]
            # end
            # println("##<G>")
            return H[hash(ψ)][hash(I)][horizon][(1,2)]

        elseif typeof(ψ.data) == UniversalIntervalRelation{:G}
            # println("#[G]")
            # for h in minh:maxh
            S = Float64[]
            for op in [univIntRel(:L),univIntRel(:A),univIntRel(:O),univIntRel(:E),univIntRel(:D),univIntRel(:B),univIntRel(:InvL),univIntRel(:InvA),univIntRel(:InvO),univIntRel(:InvE),univIntRel(:InvD),univIntRel(:InvB)]
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals
                        s = ⊓(s, ↦(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end
                    push!(S, s)
                end
            end
            H[hash(ψ)][hash(I)][horizon][(1,2)] = ⊓(S...)
            # end
            # println("##[G]")
            return H[hash(ψ)][hash(I)][horizon][(1,2)]

        elseif typeof(ψ.data) <: AbstractExistentialIntervalRelation
            # println("#ex")
            # for h in minh:maxh
                # println("###ex")
            for (x, y) ∈ intervals
                # println("####ex")
                s = 0.0
                for (w, z) ∈ intervals
                    # println("######ex")
                    s = ⊔(s, ⊓(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ⊓(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                # println("####ex")
                # @show s
                # @show typeof(s)
                L[(hash(ψ), (x, y))] = s
            end
                # println("###ex")
            # end
            # println("##ex")

        elseif typeof(ψ.data) <: AbstractUniversalIntervalRelation
            # println("#univ")
            # for h in minh:maxh
            for (x, y) ∈ intervals
                s = 1.0
                for (w, z) ∈ intervals
                    # @show H[ψ.left][I][h][(w,z)]
                    # @show typeof(H[ψ.left][I][h][(w,z)])
                    # @show H[ψ.right][I][h][(w,z)]
                    # @show typeof(H[ψ.right][I][h][(w,z)])
                    # @show R̃ₓ(ψ.data, (x, y), (w, z); horizon=h)
                    # @show typeof(R̃ₓ(ψ.data, (x, y), (w, z); horizon=h))
                    s = ⊓(s, ↦(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ↦(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                L[(hash(ψ), (x, y))] = s
            end
            # end
            # println("##univ")
        elseif typeof(ψ.data) <: BinaryRelation{:∧}
            # println("#∧")
            # for h in minh:maxh
                # println("###∧")
            for (x, y) in intervals
                # println("####∧")
                L[(hash(ψ), (x, y))] = ⊓(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
            # end
            # println("##∧")
        elseif typeof(ψ.data) <: BinaryRelation{:→}
            # println("#→")
            # for h in minh:maxh
                # println("###→")
            for (x, y) in intervals
                # println("####→")
                # @show H[ψ.left][I][h][(x,y)]
                # @show typeof(H[ψ.left][I][h][(x,y)])
                # @show H[ψ.right][I][h][(x,y)]
                # @show typeof(H[ψ.right][I][h][(x,y)])
                # error(1)
                L[(hash(ψ), (x, y))] = ↦(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
            # end
            # println("##→")
        end
    end
    return H[hash(_antecedent(φ))][hash(I)][horizon][(1,2)]
end

function _mc(I::ModalInstance, φ::ClassificationRule; horizon::Int=4)
    sub = _collect_sorted_subformulas(φ)
    L = Dict{Tuple{UInt,NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y]
    for ψ ∈ sub
        if typeof(ψ.data) <: AbstractProposition
            for (x, y) ∈ intervals # setdiff(𝕀, [(-1, 0)])
                # L[(ψ, (x, y))] = ((x, y) != (-1, 0) ? f̃(ψ.data, ts[ψ.data.attribute][x:y]) : 0.0) # TODO is (-1, 0) needed for a proposition? Perhaps not
                L[(hash(ψ), (x ,y))] = f̃(ψ.data, I[1][ψ.data.A][x:y]) # TODO works only for time series
            end
        elseif typeof(ψ.data) <: Float64
            for (x, y) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                L[(hash(ψ), (x, y))] = ψ.data
            end
        elseif typeof(ψ.data) == ExistentialIntervalRelation{:G}
            S = Float64[]
            for op in [exIntRel(:L),exIntRel(:A),exIntRel(:O),exIntRel(:E),exIntRel(:D),exIntRel(:B),exIntRel(:InvL),exIntRel(:InvL),exIntRel(:InvA),exIntRel(:InvO),exIntRel(:InvE),exIntRel(:InvD),exIntRel(:InvB)]
                for (x, y) ∈ intervals
                    s = 0.0
                    for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                        s = ⊔(s, ⊓(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end
                    push!(S, s)
                end
            end
            return ⊔(S...)
        elseif typeof(ψ.data) == UniversalIntervalRelation{:G}
            S = Float64[]
            for op in [univIntRel(:L),univIntRel(:A),univIntRel(:O),univIntRel(:E),univIntRel(:D),univIntRel(:B),univIntRel(:InvL),univIntRel(:InvA),univIntRel(:InvO),univIntRel(:InvE),univIntRel(:InvD),univIntRel(:InvB)]
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                        s = ⊓(s, ↦(R̃ₓ(op, (x, y), (w, z); horizon=horizon), L[(hash(ψ.right), (w, z))]))
                    end
                    push!(S, s)
                end
            end
            return ⊓(S...)
        elseif typeof(ψ.data) <: AbstractExistentialIntervalRelation
            for (x, y) ∈ intervals # (typeof(ψ.data) === ExistentialLater ? 𝕀 : setdiff(𝕀,[(-1, 0)]))
                s = 0.0
                for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                    s = ⊔(s, ⊓(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ⊓(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                L[(hash(ψ), (x, y))] = s
            end
        elseif typeof(ψ.data) <: AbstractUniversalIntervalRelation
            for (x, y) ∈ intervals # (typeof(ψ.data) === UniversalLater ? 𝕀 : setdiff(𝕀,[(-1,0)]))
                s = 1.0
                for (w, z) ∈ intervals # setdiff(𝕀,[(-1,0)])
                    s = ⊓(s, ↦(R̃ₓ(ψ.data, (x, y), (w, z); horizon=horizon), ↦(L[(hash(ψ.left), (w, z))], L[(hash(ψ.right), (w, z))])))
                end
                L[(hash(ψ), (x, y))] = s
            end
        elseif typeof(ψ.data) <: BinaryRelation{:∧}
            for (x, y) in intervals
                # println("####∧")
                L[(hash(ψ), (x, y))] = ⊓(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
        elseif typeof(ψ.data) <: BinaryRelation{:→}
            for (x, y) in intervals
                # println("####∧")
                L[(hash(ψ), (x, y))] = ↦(L[(hash(ψ.left), (x, y))], L[(hash(ψ.right), (x, y))])
            end
        end
    end
end

function mmmc(ℐ::ClassificationDataset, Γ::ClassificationRules)
    accuracy = 0.0

    Threads.@threads for i in 1:length(ℐ.instances)
        I = ℐ.instances[i]
        class = ℐ.classes[i]
        L = Dict{ClassificationRule, Float64}() # ρ ↦ truth value

        for ρ in Γ
            # lock(lk) do
            L[ρ] = mc(I, ρ; minh=Γ.minh, maxh=Γ.maxh, horizon=Γ.horizon, memo=true)
            # end
        end

        maxv = argmax(L)
        maxv_rules = [k for (k,v) ∈ L if v == L[maxv]]
        pred = _consequent(maxv_rules[rand(1:length(maxv_rules))])

        if pred == class
            accuracy = accuracy + 1.0
        end
    end
    return accuracy/length(ℐ)
end

function z(Γ::ClassificationRules)
    round(mmmc(Γ.ds, Γ), digits=4), round(1 - _complexity(Γ)/_maxsize(Γ), digits=4) # 1 - sum(sizes)/(max_size*maxr) #1 - sum(sizes)/(max_size*100) # (mean_size(crules)/max_size) * (length(classes(td))/length(rules(crules))) # 1.0 - sum(sizes)/(max_size*12)
end

function CV(Γ::ClassificationRules)
    return abs(length(Γ.ds.unique_classes) - length(_classes(Γ)))
end

function init()
    ex_relations    = [exIntRel(s) for s in [:L,:A,:O,:E,:D,:B,:InvL,:InvA,:InvO,:InvE,:InvD,:InvB]]
    univ_relations  = [univIntRel(s) for s in [:L,:A,:O,:E,:D,:B,:InvL,:InvA,:InvO,:InvE,:InvD,:InvB]]
    init_relations  = [exIntRel(:G),univIntRel(:G)]
    prop_relations  = [conjunction, implication]
    relations       = vcat(ex_relations, univ_relations,prop_relations)
    _rand_rules(init_relations,relations,HeytingChainAlgebra,train,[≤,<,==,>,≥]; minh=1,maxh=3,maxmd=4,maxdepth=6,minnumrules=4,maxnumrules=48)
end

function hypervolume_indicator(y::Vector{Tuple{Float64,Float64}}; refpoint=(0,0))
    sort!(y, by=x->x[1], rev=true, alg=QuickSort)

    hi = 0.0
    # (index, tuple)
    for (i, t) in enumerate(y)
        # if it is the first, then
        if i == 1
            hi += abs(t[1] - refpoint[1]) * abs(t[2] - refpoint[2])
        else
            s = y[i-1]
            hi += abs(t[1] - refpoint[1]) * abs(t[2] - s[2])
        end
    end

    return round(hi, digits=4)
end

function _maxsize(Γ::ClassificationRules)
    maxmd = Γ.maxmd
    maxnumrules = Γ.maxnumrules
    maxdepth = Γ.maxdepth

    # 1 for the root
    # 2^(maxmd-1) for the nodes with modalities, -1 for the root; e.g., ◊(φ ∧ φ)
    #
    return (1 + 2*(2^(maxmd-1) - 1) + 2^(maxdepth-maxmd) + 2^(maxmd-1) + 2) * maxnumrules #(1 + 2 + (maxmd > 2 ? ((maxmd-1) * 2^(maxmd-2)) : 0 ) + 2^(maxmd-1))*maxnumrules + maxnumrules*2
end

#################################
# @@@ End of CORE FUNCTIONS @@@ #
#################################


text:
for (TT,init_f,clear_f) in [(:ArbReal,@libarb(_arb_vec_init),@libarb(_arb_vec_clear)),
                            (:ArbComplex,@libarb(_acb_vec_init),@libarb(_acb_vec_clear))]
    @eval begin
        init(::Type{T}, len::Integer) where {T<:$TT} =
            ccall($init_f, Ptr{T}, (Int32,), len)

        clear!(ptr::Ptr{T}, len::Integer) where {T<:$TT} =
            ccall($clear_f, Cvoid, (Ref{T}, Cint), ptr, len)
    end
end

"""
    ArblibVector{T}(ptr, len)

Vector of length `len` containing elements of type `T` with memory
handled by libarb.
"""
mutable struct ArblibVector{T} <: AbstractVector{T}
    ptr::Ptr{T}
    len::Int32

    """
    ArblibVector(::Type{T}, len)

Allocate uninitialized vector of length `len` containing elements of type `T`.
"""
    ArblibVector(::Type{T}, len::Int) where T = new{T}(init(T, len), len)
end

const ArbVector{P} = ArblibVector{ArbReal{P}}
const AcbVector{P} = ArblibVector{ArbComplex{P}}

Base.eltype(v::ArblibVector{T}) where T = T
Base.length(v::ArblibVector) = v.len
Base.size(v::ArblibVector) = (v.len,)

function Base.getindex(v::ArblibVector{T}, i::Integer) where T
    i ∈ 1:v.len || throw(BoundsError(v, i))
    unsafe_load(v.ptr, i)
end

"""
    ArblibVector(v::AbstractVector{T})

Construct `ArblibVector` and fill with data from `v`.
"""
function ArblibVector(v::AbstractVector{T}) where T
    vv = ArblibVector(T, length(v))
    @inbounds for (i,e) in enumerate(v)
        unsafe_store!(vv.ptr, e, i)
    end
    vv
end

# ArblibVectors of ArbFloats are stored as ArbReals, since there is no
# arf_ptr available.
ArblibVector(v::AbstractVector{ArbFloat{P}}) where P =
    ArblibVector(ArbReal.(v))

"""
    Vector(v::ArblibVector{T})

Construct new `Vector` with elements copied from `v`.
"""
function Base.Vector(v::ArblibVector{T}) where T
    vv = Vector{T}(undef, v.len)
    @inbounds for i = 1:v.len
        vv[i] = unsafe_load(v.ptr, i)
    end
    vv
end

"""
    free!(v::ArblibVector{T})

Free the memory allocated for `v` and null the pointer and length to
prevent future accesses to released memory.
"""
function free!(v::ArblibVector{T}) where T
    clear!(v.ptr, v.len)
    v.ptr = C_NULL
    v.len = 0
end

"""
    take!(v::ArblibVector)

Copy the data of `v` to a [`Vector`](@ref) and release the memory.
"""
function Base.take!(v::ArblibVector)
    vv = Vector(v)
    free!(v)
    vv
end


text:
# ---
# title: 530. Minimum Absolute Difference in BST
# id: problem530
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/minimum-absolute-difference-in-bst/description/>
# hidden: true
# ---
#
# Given a binary search tree with non-negative values, find the minimum
# [absolute difference](https://en.wikipedia.org/wiki/Absolute_difference)
# between values of any two nodes.
#
# **Example:**
#
#
#
#     Input:
#
#        1
#         \
#          3
#         /
#        2
#
#     Output:
#     1
#
#     Explanation:
#     The minimum absolute difference is 1, which is the difference between 2 and 1 (or between 2 and 3).
#
#
#
#
# **Note:**
#
#   * There are at least two nodes in this BST.
#   * This question is the same as 783: <https://leetcode.comhttps://leetcode.com/problems/minimum-distance-between-bst-nodes/>
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TTT
    include("../../examples/tictactoe.jl")
end

@testset "tictactoe" begin
    g = TTT.TicTacToe()
    reset!(g)
    while !terminated(g)
        @test sum(act!(g, rand(actions(g)))) == 0
        render(g)
    end
end


text:
# We get cube diameters from a cubemap. The cubemap of a 3 × 3 image is an array with the
# following structure,
#
#      1 2 3 4 5
#    1 ⋅ — ⋅ — ⋅
#    2 | □ | □ |
#    3 ⋅ — ⋅ — ⋅
#    4 | □ | □ |
#    5 ⋅ — ⋅ — ⋅
#
# where `⋅` are the vertices, `—` and `|` are 1-cubes and `□` are 2-cubes. Each value in the
# resulting array is equal to the birth time of a cube, which is equal to the `maximum` of
# the birth times of cofaces.
#
# See https://link.springer.com/chapter/10.1007%2F978-3-642-23175-9_7 for more info.
@generated function _cubemap(input::Array{T,N}) where {T,N}
    quote
        result = similar(input, size(input) .* 2 .- 1)
        result .= typemin(T)
        @inbounds @nloops $N i input begin
            index = CartesianIndex{$N}(@ntuple $N i)
            value = input[index]
            @nexprs $N j -> (range_j = max(2i_j - 2, 1):min(2i_j, size(result, j)))
            (@nref $N result range) .= max.((@nref $N result range), value)
        end
        return result
    end
end

_one_hot(i, ::Val{N}) where {N} = CartesianIndex{N}(ntuple(isequal(i), Val(N)))

@inline _map(_, ::Tuple{}) = ()
@inline _map(f, t) = tuple(f(t[1]), _map(f, TupleTools.tail(t))...)

# Convenience functions from converting cubemap index to vertices and back.
@inline function _from_cubemap(root::CartesianIndex{K}, ::Val{N}) where {K,N}
    2^count(iseven, Tuple(root)) == N || throw(ArgumentError("invalid N"))
    result = ntuple(_ -> root, Val(N))
    for (i, j) in enumerate(Tuple(root))
        if iseven(j)
            result = result .+ ntuple(Val(N)) do k
                ifelse(isodd(k), 1, -1) * _one_hot(i, Val(K))
            end
        end
        result = TupleTools.sort(result; by=Base.Fix2(getindex, i))
    end
    result = TupleTools.sort(result)
    return _map(result) do c
        CartesianIndex{K}((Tuple(c) .+ 1) .÷ 2)
    end
end

function _to_cubemap(vertices::NTuple{N}) where {N}
    K = length(first(vertices))
    result = ntuple(Val(K)) do i
        sum(2 .* getindex.(vertices, i) .- 1) ÷ N
    end
    return CartesianIndex{K}(result)
end

function _is_valid(vertices::Tuple{})
    return false
end
function _is_valid(vertices::NTuple{1})
    return true
end
function _is_valid(vertices::NTuple{N}) where {N}
    floor(log2(N)) == log2(N) || return false
    allunique(vertices) || return false
    K = length(first(vertices))
    for i in 1:K
        l, h = extrema(getindex.(vertices, i))
        h - l ≤ 1 || return false
    end
    return true
end

"""
    Cube{D, T, K} <: AbstractCell{D, T, CartesianIndex{K}}

A `Cube` is similar to a `Simplex`, but it has `2^D` vertices instead of `D+1`. The vertices
are encoded as the position in the CubeMap (see reference in [`Cubical`](@ref)). A `Cube`'s
vertices are of type `CartesianIndex{K}`.

# Example

```jldoctest
julia> Cube{1}(CartesianIndex(1, 2), 1.0)
Cube{1, Float64, 2}((1, 2), 1.0)

```
"""
struct Cube{D,T,K} <: AbstractCell{D,T,CartesianIndex{K}}
    root::NTuple{K,Int32}
    birth::T

    function Cube{D,T,K}(root::CartesianIndex{K}, birth) where {D,T,K}
        D < 0 && throw(DomainError("dimension must be a non-negative integer"))
        return new{D,T,K}(Int32.(Tuple(root)), convert(T, birth))
    end
end
function Cube{D}(root::CartesianIndex{K}, birth::T) where {D,T,K}
    return Cube{D,T,K}(root, birth)
end
function Cube{D}(vertices, birth) where {D}
    vs = tuple(CartesianIndex.(vertices)...)
    _is_valid(vs) || throw(ArgumentError("invalid vertices"))
    root = _to_cubemap(vs)
    return Cube{D}(root, birth)
end

birth(cube::Cube) = cube.birth
index(cube::Cube) = CartesianIndex(cube.root)
Base.sign(cube::Cube) = 1
Base.:-(cube::Cube) = cube
Base.abs(cube::Cube) = cube

@generated Base.length(::Type{<:Cube{D}}) where {D} = :($(2^D))
vertices(cube::Cube{D}) where {D} = _from_cubemap(index(cube), Val(length(cube)))

"""
    Cubical{T, K} <: AbstractFiltration{CartesianIndex{K}, T}

`Cubical` is used to compute sublevel persistent homology on `N`-dimensional images, which
are of type `AbstractArray{T, N}`.

This type uses the CubeMap structure to find birth times of cubes (see reference).

# Constructor

* `Cubical(image::AbstractArray{T, N}, threshold=maximum(image))`

# Reference

Wagner, H., Chen, C., & Vuçini, E. (2012). [Efficient computation of persistent homology for
cubical data.](https://link.springer.com/chapter/10.1007/978-3-642-23175-9_7) In Topological
methods in data analysis and visualization II (pp. 91-106). Springer, Berlin, Heidelberg.

# Example

```jldoctest
julia> image = [1 0 0; 0 2 0; 0 0 0];

julia> ripserer(Cubical(image))[1]
1-element 0-dimensional PersistenceDiagram:
 [0.0, ∞)

julia> ripserer(Cubical(image))[2]
1-element 1-dimensional PersistenceDiagram:
 [1.0, 2.0)

```
"""
struct Cubical{K,T,A<:AbstractArray{T,K}} <: AbstractFiltration{CartesianIndex{K},T}
    data::A
    cubemap::A
    threshold::T
end

function Cubical(
    data::AbstractArray{T,K}; threshold=maximum(data), verbose=false
) where {T,K}
    return Cubical{K,T,typeof(data)}(data, _cubemap(data), T(threshold))
end

nv(cf::Cubical) = length(cf.data)
threshold(cf::Cubical) = cf.threshold
simplex_type(::Type{<:Cubical{K,T}}, D) where {K,T} = Cube{D,T,K}

# when working with vertices, we don't cubemap them.
births(cf::Cubical) = cf.data
vertices(cf::Cubical) = CartesianIndices(cf.data)

function edges(cf::Cubical{K}) where {K}
    result = edge_type(cf)[]
    for i in CartesianIndices(cf.cubemap)
        if count(iseven, Tuple(i)) == 1
            edge = unsafe_simplex(cf, Val(1), i, 1)
            if !isnothing(edge)
                push!(result, edge)
            end
        end
    end
    return result
end

function simplex(cf::Cubical{N,T}, ::Val{D}, vertices, sign=1) where {D,T,N}
    if _is_valid(vertices)
        root = _to_cubemap(vertices)
        return unsafe_simplex(cf, Val(D), root, sign)
    else
        return nothing
    end
end

function unsafe_simplex(cf::Cubical{N,T}, ::Val{D}, new_root, _) where {D,T,N}
    birth = get(cf.cubemap, new_root, missing)
    if ismissing(birth) || birth > cf.threshold
        return nothing
    else
        return Cube{D,T,N}(new_root, birth)
    end
end

struct CubeCoboundary{K,D,F<:Cubical{K},C<:Cube{D}}
    filtration::F
    cube::C
end

function coboundary(filtration::Cubical{K}, cube::Cube{D}) where {K,D}
    return CubeCoboundary{K,D,typeof(filtration),typeof(cube)}(filtration, cube)
end

function Base.iterate(cob::CubeCoboundary{K,D}, (i, sign)=(K, 1)) where {K,D}
    orig_root = index(cob.cube)
    while i ≤ K
        # To get coboundary in decreasing root order, we have to traverse i from K to 1 and
        # then back to K with sign reversed.
        next_i = i - sign
        next_sign = sign
        if next_i ≤ 0
            next_sign = -1
            next_i = 1
        end

        if isodd(orig_root[i])
            new_root = orig_root + sign * _one_hot(i, Val(K))
            cofacet = unsafe_simplex(cob.filtration, Val(D + 1), new_root, sign)
            if !isnothing(cofacet)
                _cofacet::simplex_type(cob.filtration, D + 1) = cofacet
                return _cofacet, (next_i, next_sign)
            end
        end

        i, sign = next_i, next_sign
    end
    return nothing
end

struct CubeBoundary{K,D,F<:Cubical{K},C<:Cube{D}}
    filtration::F
    cube::C
end

function boundary(filtration::Cubical{K}, cube::Cube{D}) where {K,D}
    return CubeBoundary{K,D,typeof(filtration),typeof(cube)}(filtration, cube)
end

function Base.iterate(bnd::CubeBoundary{K,D}, (i, sign)=(K, -1)) where {K,D}
    orig_root = index(bnd.cube)
    while i ≤ K
        next_i = i + sign
        next_sign = sign
        if next_i ≤ 0
            next_sign = 1
            next_i = 1
        end

        if iseven(orig_root[i])
            new_root = orig_root + sign * _one_hot(i, Val(K))
            facet = unsafe_simplex(bnd.filtration, Val(D - 1), new_root, sign)
            if !isnothing(facet)
                _facet::simplex_type(bnd.filtration, D - 1) = facet
                return _facet, (next_i, next_sign)
            end
        end

        i, sign = next_i, next_sign
    end
    return nothing
end

struct CubicalColumnsToReduce{D,K,F<:Cubical{K}}
    filtration::F
end

function columns_to_reduce(filtration::Cubical{K}, itr) where {K}
    D = dim(eltype(itr)) + 1
    return CubicalColumnsToReduce{D,K,typeof(filtration)}(filtration)
end

Base.eltype(::CubicalColumnsToReduce{D,<:Any,F}) where {D,F} = simplex_type(F, D)

function Base.iterate(cols::CubicalColumnsToReduce{D,K}, i=1) where {D,K}
    last_i = lastindex(cols.filtration.cubemap)
    cart = CartesianIndices(cols.filtration.cubemap)
    while i < last_i
        root = cart[i]
        if count(iseven, Tuple(root)) == D
            simplex = unsafe_simplex(cols.filtration, Val(D), root, 1)
            if !isnothing(simplex)
                _simplex::simplex_type(cols.filtration, D) = simplex
                return simplex, i + 1
            end
        end
        i += 1
    end
    return nothing
end

struct CubicalDist{C} <: AbstractMatrix{Int}
    filtration::C
end

Base.size(cd::CubicalDist) = (nv(cd.filtration), nv(cd.filtration))
function Base.getindex(cd::CubicalDist, i::Integer, j::Integer)
    ci = CartesianIndices(cd.filtration.data)[i]
    cj = CartesianIndices(cd.filtration.data)[j]
    return sum(abs, Tuple(ci) .- Tuple(cj))
end

distance_matrix(cf::Cubical) = CubicalDist(cf)


text:
# this function shows the behavior of the
# likelihood, the log-likelihood and the average
# log-likelihood as the sample size increases
using Distributions
n = 10 # try increasing this to 50, 1000, etc.
p = 0.3
y = rand(Bernoulli(p),n)
f =  pdf.(Bernoulli(p), y)
@show n
@show ℒ = prod(f)
@show lnℒ = sum(log.(f))
@show avglnℒ = (1/n)*sum(log.(f))
nothing

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
using Compat.Test
using Compat.LinearAlgebra
using MathProgBase

function mixintprogtest(solver)
    @testset "Testing mixintprog with $solver" begin
        # integer knapsack problem
        sol = mixintprog(-[5,3,2,7,4],Float64[2 8 4 2 5],'<',10,:Int,0,1,solver)
        @test sol.status == :Optimal
        @test isapprox(sol.objval, -16.0, atol=1e-6)
        @test isapprox(norm(sol.sol - [1.0, 0.0, 0.0, 1.0, 1.0]), 0.0, atol=1e-4)
    end
end


text:
using Luxor, Contour, Images

function geostats(s)
  Δ = s/500 # scale factor, original design was 500units
  circle(O, 245Δ, :clip)
  sethue("white")
  circle(O, 245Δ, :fill)

  # generate contours from image
  pwd()
  i = load("contours.png")
  img = imresize(i, ratio=1.75)
  H, W = size(img)
  imgg = imfilter(img, Kernel.gaussian(15))
  x = -W/2.0:W/2.0 - 1  # W elements
  y = -H/2.0:H/2.0 - 1  # H elements
  z = Float64.(Gray.(imgg))
  # just two contours required
  clevels = Contour.contourlevels(z, 2)
  contourlines = contours(x, y, z, clevels)

  # draw the contours
  @layer begin
    scale(Δ)
    setline(8Δ)
    for (n, cl) in enumerate(levels(contourlines))
      for l in lines(cl)
        # convert contour points to Luxor coords
        sethue("black")
        pgon = map(p -> Point(p[2], p[1]), zip(coordinates(l)...))
        poly(pgon, :stroke)

        # optional coloring?
        setgray([.6, .8][n])
        poly(pgon, :fill, close=true)
      end
    end
  end

  clipreset()

  sethue("black")
  setline(8Δ)
  circle(O, 245Δ, :stroke)

  sethue(Luxor.julia_purple)
  circle(Point(190, 120)Δ, 55Δ, :fill)

  sethue(Luxor.julia_green)
  circle(Point(70, -55)Δ, 50Δ, :fill)

  sethue(Luxor.julia_red)
  circle(Point(-180, 12)Δ, 40Δ, :fill)
end

function logo(s, fname)
  Drawing(s, s, fname)
  origin()
  geostats(s)
  finish()
  preview()
end

function logotext(w, h, fname)
  Drawing(w, h, fname)
  #background("white")
  origin()
  table = Table([h], [h, w - h])
  @layer begin
    translate(table[2])
    background("white") # comment for transparent background
    sethue("black")
    # find all fonts available on Linux with `fc-list | -f 2 -d ":"`
    fontface("Julius Sans One")
    fontsize(h/2.5)
    text("GeoStats.jl", halign=:center, valign=:middle)
  end
  @layer begin
    translate(table[1])
    geostats(h)
  end
  finish()
  preview()
end

for ext in [".svg", ".png"]
  logo(120, "../assets/logo"*ext)
  logotext(350, 100, "../assets/logo-text"*ext)
end


text:
using CpuId

@testset "cpu_times" begin
    if current_system() isa Linux
        record = PSUtil.cpu_times()
        fields = (:user, :nice, :system, :idle, :iowait, :irq, :softirq, :steal, :guest, :guest_nice)
        @test length(record) == length(fields)
        @test keys(record) == fields
        @test eltype(record) == Int
    else
        @test_broken False
        @info "`cpu_times` is not tested on $(current_system())"
    end
end

@testset "cpu_percent" begin
    if current_system() isa Linux
        record = PSUtil.cpu_percent()

        @test 0<= PSUtil.cpu_percent() <= 100
        @test 0<= PSUtil.cpu_percent(1) <= 100
        @test 0<= first(PSUtil.cpu_percent(0, 0)) <= 100
        @test 0<= first(PSUtil.cpu_percent(0, 0, 1)) <= 100

        # test if it actaully blocks the process
        t = @elapsed PSUtil.cpu_percent(;interval=1)
        @test 1< t < 1.5

        # test if the concurrent version only blocks about 1 seconds
        t = @elapsed map(t->t.result,
            @sync map(0:PSUtil.cpu_count()-1) do i
                t = @async PSUtil.cpu_percent(i; interval=1)
            end
        )
        @test 1< t < 1.5

        # test if the non-blocking version doesn't block the routine
        PSUtil.cpu_percent(0, 0) # hit precompilation
        t = @elapsed PSUtil.cpu_percent(0, 0)
        @test t < 0.1

        err_msg = @capture_err PSUtil.cpu_percent(0, 0; interval=1)
        @test occursin("keyword `interval` is not used in non-block version.", err_msg)
    else
        # TODO: cpu_times is only supported on Linux now and thus we can't test this `cpu_percent` on other platforms
        @test_broken False
        @info "`cpu_percent` is not tested on $(current_system())"
    end
end

@testset "cpu_count" begin
    @test PSUtil.cpu_count() == PSUtil.cpu_count(logical=true)
    @test CpuId.cputhreads() == PSUtil.cpu_count(logical=true)
    @test CpuId.cpucores() ==  PSUtil.cpu_count(logical=false)
end


text:
# indicator of second-order cones

export IndSOC, IndRotatedSOC

"""
**Indicator of the second-order cone**

    IndSOC()

Returns the indicator of the second-order cone (also known as ice-cream cone or Lorentz cone), that is
```math
C = \\left\\{ (t, x) : \\|x\\| \\leq t \\right\\}.
```
"""
struct IndSOC <: ProximableFunction end

function (f::IndSOC)(x::AbstractVector{T}) where T <: Real
  # the tolerance in the following line should be customizable
  if norm(x[2:end]) - x[1] <= 1e-14
    return zero(T)
  end
  return +Inf
end

is_convex(f::IndSOC) = true
is_set(f::IndSOC) = true

function prox!(y::AbstractVector{T}, f::IndSOC, x::AbstractVector{T}, gamma::T=one(T)) where T <: Real
  @views nx = norm(x[2:end])
  t = x[1]
  if t <= -nx
    y .= zero(T)
  elseif t >= nx
    y .= x
  else
    r = T(0.5) * (one(T) + t / nx)
    y[1] = r * nx
    @views y[2:end] .= r .* x[2:end]
  end
  return zero(T)
end

fun_name(f::IndSOC) = "indicator of the second-order cone"
fun_dom(f::IndSOC) = "AbstractArray{Real,1}"
fun_expr(f::IndSOC) = "x ↦ 0 if x[1] >= ||x[2:end]||, +∞ otherwise"
fun_params(f::IndSOC) = "none"

function prox_naive(f::IndSOC, x::AbstractVector{T}, gamma=1.0) where T <: Real
  nx = norm(x[2:end])
  t = x[1]
  if t <= -nx
    y = zero(x)
  elseif t >= nx
    y = x
  else
    y = zero(x)
    r = 0.5 * (1 + t / nx)
    y[1] = r * nx
    y[2:end] .= r .* x[2:end]
  end
  return y, 0.0
end

# ########################
# ROTATED SOC
# ########################

"""
**Indicator of the rotated second-order cone**

    IndRotatedSOC()

Returns the indicator of the *rotated* second-order cone (also known as ice-cream cone or Lorentz cone), that is
```math
C = \\left\\{ (p, q, x) : \\|x\\|^2 \\leq 2\\cdot pq, p \\geq 0, q \\geq 0 \\right\\}.
```
"""
struct IndRotatedSOC <: ProximableFunction end

function (f::IndRotatedSOC)(x::AbstractVector{T}) where T <: Real
  if x[1] >= -1e-14 && x[2] >= -1e-14 && norm(x[3:end])^2 - 2*x[1]*x[2] <= 1e-14
    return zero(T)
  end
  return +Inf
end

is_convex(f::IndRotatedSOC) = true
is_set(f::IndRotatedSOC) = true

function prox!(y::AbstractVector{T}, f::IndRotatedSOC, x::AbstractVector{T}, gamma::T=one(T)) where T <: Real
  # sin(pi/4) = cos(pi/4) = 0.7071067811865475
  # rotate x ccw by pi/4
  x1 = 0.7071067811865475*x[1] + 0.7071067811865475*x[2]
  x2 = 0.7071067811865475*x[1] - 0.7071067811865475*x[2]
  # project rotated x onto SOC
  @views nx = sqrt(x2^2+norm(x[3:end])^2)
  t = x1
  if t <= -nx
    y .= zero(T)
  elseif t >= nx
    y[1] = x1
    y[2] = x2
    @views y[3:end] .= x[3:end]
  else
    r = T(0.5) * (one(T) + t / nx)
    y[1] = r * nx
    y[2] = r * x2
    @views y[3:end] = r .* x[3:end]
  end
  # rotate back y cw by pi/4
  y1 = 0.7071067811865475*y[1] + 0.7071067811865475*y[2]
  y2 = 0.7071067811865475*y[1] - 0.7071067811865475*y[2]
  y[1] = y1
  y[2] = y2
  return zero(T)
end

fun_name(f::IndRotatedSOC) = "indicator of the rotated second-order cone"
fun_dom(f::IndRotatedSOC) = "AbstractArray{Real,1}"
fun_expr(f::IndRotatedSOC) = "x ↦ 0 if x[1] ⩾ 0, x[2] ⩾ 0, norm(x[3:end])² ⩽ 2*x[1]*x[2], +∞ otherwise"
fun_params(f::IndRotatedSOC) = "none"

function prox_naive(f::IndRotatedSOC, x::AbstractVector{T}, gamma=1.0) where T <: Real
  g = IndSOC()
  z = copy(x)
  z[1] = 0.7071067811865475*x[1] + 0.7071067811865475*x[2]
  z[2] = 0.7071067811865475*x[1] - 0.7071067811865475*x[2]
  y, = prox_naive(g, z, gamma)
  y1 = 0.7071067811865475*y[1] + 0.7071067811865475*y[2]
  y2 = 0.7071067811865475*y[1] - 0.7071067811865475*y[2]
  y[1] = y1
  y[2] = y2
  return y, 0.0
end


text:
@testset "Subpopulation report" begin

@testset "function subpopulationPopReport" begin
    subpop1 = Subpopulation( "Subpop1" )
    addSubpopulationCondition!( subpop1, MPcondition( "branch", !=, "A" ) )
    subpop2 = Subpopulation( "Subpop2" )
    addSubpopulationCondition!( subpop2, MPcondition( "time in node", >=, 84 ) )
    subpop3 = Subpopulation( "Subpop3" )
    addSubpopulationCondition!( subpop3, MPcondition( "started as", ==,
        "Reserve junior" ) )
    subpop4 = Subpopulation( "Subpop4" )
    addSubpopulationCondition!( subpop4, MPcondition( "tenure", <, 12 ) )
    subpop5 = Subpopulation( "Subpop5" )
    addSubpopulationCondition!( subpop5, MPcondition( "was", !=, "B junior" ) )
    subpop6 = Subpopulation( "Subpop6" )
    subpop7 = Subpopulation( "Subpop7" )
    addSubpopulationCondition!( subpop7, MPcondition( "age", >=, 240 ) )

    setSubpopulationSourceNode!.( [subpop1, subpop2, subpop3, subpop4,
        subpop5, subpop6, subpop7], ["", "", "", "Reserve junior", "Career",
        "Branch A", "Foo"] )

    @test isempty( subpopulationPopReport( mpSim, [-12.0, -6.0], subpop1,
        subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 ) )
    @test isempty( subpopulationPopReport( mpSim, [12.0, 6.0], subpop7 ) )
    @test isempty( subpopulationPopReport( mpSim, -12, subpop1, subpop2,
        subpop3, subpop4, subpop5, subpop6, subpop7 ) )

    nodeReport = nodePopReport( mpSim, 12, "A junior", "B junior",
        "Reserve junior", "A senior", "B senior", "Master", "" )

    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3,
        subpop4, subpop5, subpop6, subpop7 )
    @test all( report[:, "Subpop1"] .==
        nodeReport[:, "B junior"] +
        nodeReport[:, "Reserve junior"] +
        nodeReport[:, "B senior"] +
        nodeReport[:, "Master"] )
    @test all( report[:, "Subpop2"] .==
        vcat( fill( 0, 7 ), 24, 48, 64, 56, fill( 48, 15 ) ) )
    @test all( report[:, "Subpop3"] .==
        vcat( 10, 20, 28, fill( 36, 7 ), 28, fill( 20, 15 ) ) )
    @test all( report[:, "Subpop4"] .== 10 )
    @test all( report[:, "Subpop5"] .==
        vcat( 0, 0, 16, 32:8:80, 72, fill( 64, 15 ) ) )
    @test all( report[:, "Subpop6"] .==
        nodeReport[:, "A junior"] +
        nodeReport[:, "A senior"] )
    @test "Subpop7" ∉ names( report )
end  # @testset "function subpopulationPopReport"

end  # @testset "Subpopulation report"


text:
module BouncingSphere_with_time

using ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

@time begin BouncingSphere = Model(
        boxHeigth = 0.1,
        groundMaterial = VisualMaterial(color="DarkGreen", transparency=0.5),
        gravField = UniformGravityField(g=9.81, n=[0, -1, 0]),
        world = Object3D(feature=Scene(gravityField=:gravField,
                                            visualizeFrames=false,
                                            defaultFrameLength=0.2,
                                            visualizeBoundingBox = true,
                                            enableContactDetection=true,
                                            visualizeContactPoints=false)),
        worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
        ground = Object3D(parent=:world,
                        translation=:[0.0,-boxHeigth/2,0.0],
                        feature=Solid(shape=Box(lengthX=4.0, lengthY=:boxHeigth, lengthZ=0.7),
                                        visualMaterial=:groundMaterial,
                                        solidMaterial="Steel",
                                        collision=true)),
        sphere = Object3D(feature=Solid(shape=Sphere(diameter=0.2),
                                        visualMaterial=VisualMaterial(color="Blue"),
                                        solidMaterial="Steel",
                                        massProperties=MassPropertiesFromShapeAndMass(mass=0.001),
                                        collision=true)),
        prism = Prismatic(obj1=:world, obj2=:sphere, axis=2, s=Var(init=1.0))
    )
end

@time bouncingSphere = @instantiateModel(buildModia3D(BouncingSphere), unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 2.2
dtmax = 0.1
tolerance = 1e-8
requiredFinalStates = [0.01948634300615831, 0.10092552754875266]
@time simulate!(bouncingSphere, stopTime=stopTime, tolerance=tolerance, dtmax=dtmax, log=false, logStates=false, logEvents=false, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
@time plot(bouncingSphere, ["prism.s", "prism.v"], figure=1)

end


text:
"""
Type: GenericState

Methods: update!, reinit!

A generic State to describe the state of a problem at a point x.

Tracked data include:
 - x : our current iterate
 - current_time  : time
 - current_score : score

Note: by default, unknown entries are set to nothing.
"""
mutable struct GenericState <: AbstractState

    x :: AbstractVector

    #Current time
    current_time  :: FloatVoid
    #Current score
    current_score :: FloatVoid

    function GenericState(x             :: AbstractVector;
                          current_time  :: FloatVoid = nothing,
                          current_score :: FloatVoid = nothing)

      return new(x, current_time, current_score)
   end
end

"""
update!: generic update function for the State

The function compare the kwargs and the entries of the State.
If the type of the kwargs is the same as the entry or the entry was nothing, then
it is updated.

Set convert to true, to update even incompatible types.
"""
function update!(stateatx :: AbstractState; convert = false, kwargs...)

 kwargs = Dict(kwargs)

 for k ∈ fieldnames(typeof(stateatx))
  if (k ∈ keys(kwargs)) && (convert || getfield(stateatx, k) == nothing || typeof(kwargs[k]) ∈ [typeof(getfield(stateatx, k)), Nothing])
   setfield!(stateatx, k, kwargs[k])
  end
 end

 return stateatx
end

"""
reinit!: function that set all the entries at void except the mandatory x

Note: If x is given as a keyword argument it will be prioritized over
the argument x.
"""
function reinit!(stateatx :: AbstractState, x :: Iterate; kwargs...)

 for k ∈ fieldnames(typeof(stateatx))
   if k != :x setfield!(stateatx, k, nothing) end
 end

 return update!(stateatx; x=x, kwargs...)
end

"""
reinit!: shorter version of reinit! reusing the x in the state

Note: If x is given as a keyword argument it will be prioritized over
the argument x.
"""
function reinit!(stateatx :: AbstractState; kwargs...)
 return reinit!(stateatx, stateatx.x; kwargs...)
end

"""
_domain_check: verifies is there is a NaN in the State entries

return true if a NaN has been found
"""
function _domain_check(stateatx :: AbstractState)
 domainerror = false

 for k ∈ fieldnames(typeof(stateatx))
   try domainerror = domainerror || (true in isnan.(getfield(stateatx, k))) catch end
 end

 return domainerror
end


text:
function Base.display(M::TensorMesh3D)

	println("$(M.dim)-dimensional tensor mesh of size $(M.n)")
	println("Number of cells:   $(M.nc)")
	println("Number of faces:   $(M.nf) = $(sum(M.nf))")
	println("Number of edges:   $(M.ne) = $(sum(M.ne))")
	println("Number of nodes:   $(prod(M.n.+1))")
	println("Coordinate origin: $(M.x0)")
	println("Domain size:       $(sum(M.h1))m x $(sum(M.h2))m x $(sum(M.h3))m")
	println("Minimum cell size: $(minimum(M.h1))m x $(minimum(M.h2))m x $(minimum(M.h3))m")
	println("Maximum cell size: $(maximum(M.h1))m x $(maximum(M.h2))m x $(maximum(M.h3))m")

end

function Base.display(M::RegularMesh)

	println("$(M.dim)-dimensional regular mesh of size $(M.n)")
	println("Number of cells:   $(M.nc)")
	println("Number of faces:   $(M.nf) = $(sum(M.nf))")
	println("Number of edges:   $(M.ne) = $(sum(M.ne))")
	println("Number of nodes:   $(prod(M.n.+1))")
	println("Coordinate origin: $(M.x0)")
	println("Domain size:       $(M.domain)")
	println("Cell size:         $(M.h)")

end


text:
# Packages
using Plots

# Generate time series with a doubling in volatility after 500 observations
Σ = cholesky([1 -1; -1 2])
T = 500
y = randn(T, 2) * Σ.U
Σ = cholesky([2 -2; -2 4])
T = 500
y = vcat(y, randn(T, 2) * Σ.U)

# Struct containing y and priors
s = LocalLevel(y, [5 -5], 10*ones(1, 1), Matrix(I, 2, 2), 0.95, 0.95)

# Estimate
m = estimate(s)

# Plot simulated data and estimated means
scatter(s.y, color = [:blue :red], markeralpha = 0.5, layout = (2, 1), label = "")
plot!(m.predicted.μ, color = [:blue :red], linewidth = 2; label = "")

# Plot true and estimated variances
plot([m.predicted.Σ[:, 1, 1] m.predicted.Σ[:, 2, 2]], color = [:blue :red], linewidth = 1; label = "", layout = (2, 1))
plot!([[ones(T);ones(T)*2] [ones(T)*2;ones(T)*4]], color = :grey, linewidth = 2)

# Plot estimated variances vs. variance of state
plot([m.Σ[:, 1, 1] m.Σ[:, 2, 2]], color = [:blue :red], linewidth = 1; label = "", layout = (2, 1))
plot!([m.Φ[:, 1, 1] m.Φ[:, 2, 2]], color = :grey, label = "")


text:
function load_knead(dir=ENV["ANALYSIS_FILES"])
    run(Cmd([
        "kneaddata_read_count_table",
        "--input", joinpath(dir, "kneaddata/"),
        "--output", joinpath(dir, "kneaddata", "kneaddata_read_counts.tsv")
    ]))

    df = CSV.read("/grace/sequencing/processed/mgx/kneaddata/kneaddata_read_counts.tsv", DataFrame; missingstring="NA")
    for col in names(df, r"Homo_sapiens")
        df[!, col] = map(p-> coalesce(p...), zip(df[!, col], df[!, replace(col, "Homo_sapiens"=>"hg37dec_v0.1")]))
    end
    df."decontaminated Homo_sapiens pair1" = map(x-> ismissing(x) ? missing :  x isa Real ? x : parse(Float64, x), df."decontaminated Homo_sapiens pair1")
    df."decontaminated Homo_sapiens orphan2" = map(x-> ismissing(x) ? missing :  x isa Real ? x : parse(Float64, x), df."decontaminated Homo_sapiens orphan2")

    df = select(df, Not(r"hg37"))
    CSV.write("data/read_counts.csv", df)
    return df
end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
# This files reimplements gmsh/tutorial/t3.geo in Julia.

import gmsh

model = gmsh.model
factory = model.geo

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)

model.add("t3")

# Copied from t1.jl...
lc = 1e-2
factory.addPoint(0, 0, 0, lc, 1)
factory.addPoint(.1, 0,  0, lc, 2)
factory.addPoint(.1, .3, 0, lc, 3)
factory.addPoint(0, .3, 0, lc, 4)
factory.addLine(1, 2, 1)
factory.addLine(3, 2, 2)
factory.addLine(3, 4, 3)
factory.addLine(4, 1, 4)
factory.addCurveLoop([4, 1, -2, 3], 1)
factory.addPlaneSurface([1], 1)
model.addPhysicalGroup(0, [1, 2], 1)
model.addPhysicalGroup(1, [1, 2], 2)
model.addPhysicalGroup(2, [1], 6)
model.setPhysicalName(2, 6, "My surface")
# ...end of copy

h = 0.1
angle = 90.

# Extruding the mesh in addition to the geometry works as in .geo files: the
# number of elements for each layer and the (end) height of each layer are
# specified in two vectors.
ov = factory.extrude([(2,1)], 0, 0, h, [8,2], [0.5,1])

#/ Rotational and twisted extrusions are available as well with the built-in CAD
# kernel. The last (optional) argument for the Extrude/Revolve/Twist commands
# specified whether the extruded mesh should be recombined or not.
ov = factory.revolve([(2,28)], -0.1,0,0.1, 0,1,0, -pi/2, [7])
ov = factory.twist([(2,50)], 0,0.15,0.25, -2*h,0,0, 1,0,0, angle*pi/180.,
                   [10], [], true)

model.addPhysicalGroup(3, [1, 2, ov[2][2]], 101)

factory.synchronize()
model.mesh.generate(3)
gmsh.write("t3.msh")
gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
function extensiblefunction!(@nospecialize(f::Function), varargs...)
    return extensiblefunction!(f, varargs)
end

function extensiblefunction!(@nospecialize(f::Function),
                             @nospecialize(varargs::Tuple))
    global _registry_extensibleunion_to_genericfunctions
    global _registry_genericfunctions_to_extensibleunions
    if !haskey(_registry_genericfunctions_to_extensibleunions, f)
        _registry_genericfunctions_to_extensibleunions[f] = Set{Any}()
    end
    for i = 1:length(varargs)
        if isextensibleunion(varargs[i])
            push!(_registry_extensibleunion_to_genericfunctions[varargs[i]],
                  f)
            push!(_registry_genericfunctions_to_extensibleunions[f],
                  varargs[i])
        else
            throw(ArgumentError(
                "Argument is not a registered extensible union."))
        end
    end
    _update_all_methods_for_extensiblefunction!(f)
    return f
end

function isextensiblefunction(@nospecialize(f::Function))
    global _registry_genericfunctions_to_extensibleunions
    return haskey(_registry_genericfunctions_to_extensibleunions, f)
end


text:
using Plots,LaTeXStrings
using LinearAlgebra,FundamentalsNumericalComputation

gr()
default(legend=nothing,linewidth=2)
plt = plot(aspect_ratio=1,showaxis=false,grid=false,xaxis=:none,yaxis=false)

A=[-2;1];
b=[3,2];

plot!([-3*A[1],3*A[1]],[-3,3]*A[2],l=(:dash,:black))
plot!([0,b[1]],[0,b[2]],l=:darkred)
x = A\b;y = A*x;
plot!([0,y[1]],[0,y[2]],l=:darkblue)
plot!([b[1], y[1]],[b[2], y[2]],l=:darkblue)

scatter!([0],[0],m=(5,:black,:black))

sq = [0.4im, -0.4+0.4im, -0.4];
sq = y[1]+1im*y[2] .+ sq*sign(y[1]+1im*y[2]);
plot!(real(sq),imag(sq),l=:darkgray)
xlims!(-3.5,3.5)
ylims!(-3.5,3.5)

annotate!(-2,0.25,text("range",:right))
annotate!(-1.93,0.25,text(L"(\!\mathbf{A})",:left))
annotate!(.7,0.9,text(L"\mathbf{b}"))
annotate!(1.,-0.8,text(L"\mathbf{A} \mathbf{x}",:right))
annotate!(2.4,0.4,text(L"\mathbf{b} - \mathbf{A} \mathbf{x}",:left))

savefig(plt,"normaleqns2d.svg")

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
using Test

@testset "MuttsDicts" begin
    @testset "MuttsDicts" begin
        include("MuttsDicts.jl")
    end
end


text:
using Documenter, StatFiles

makedocs(
	modules = [StatFiles],
	sitename = "StatFiles.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/StatFiles.jl.git"
)


text:
using OpenTelemetryAPI
using OpenTelemetryInstrumentationBase

OpenTelemetryInstrumentationBase.init()


text:
@testset "O(logn)" begin
    d = Drawing(400, 200, "test.png")
    background("white")
    sethue("black")
    latex(L"\mathcal{O}(\log{n})", 50)
    translate(50, 70)
    orangeblue = blend(O, Point(0, 50), "orange", "blue")
    setblend(orangeblue)
    latex(L"\mathcal{O}\left(\frac{\log{x}}{2}\right)", 20)
    finish()
    @test_reference "refs/ologn.png" load("test.png")
    rm("test.png")
end


text:
module Lib

using Revise

include("neat/Neat.jl")
using .Neat
export sigmoid
export Connection
export Node, NodeType
export Network

lol = function ()
    println("debug")

    conf = Dict(
        "num_inputs" => 2,
        "num_outputs" => 1
    )

    net = Network(conf)
    n = Node(5, INPUT::NodeType)
    @show n
    Neat.add_node!(net, n)

    println(net)
    #conn = Connection(1, 2)
    #add_connection!(n1, conn)
end

test = function ()

return "lib"

end

end # module


text:
# ---
# title: 770. Basic Calculator IV
# id: problem770
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, String, Stack
# link: <https://leetcode.com/problems/basic-calculator-iv/description/>
# hidden: true
# ---
#
# Given an `expression` such as `expression = "e + 8 - a + 5"` and an evaluation
# map such as `{"e": 1}` (given in terms of `evalvars = ["e"]` and `evalints =
# [1]`), return a list of tokens representing the simplified expression, such as
# `["-1*a","14"]`
#
#   * An expression alternates chunks and symbols, with a space separating each chunk and symbol.
#   * A chunk is either an expression in parentheses, a variable, or a non-negative integer.
#   * A variable is a string of lowercase letters (not including digits.) Note that variables can be multiple letters, and note that variables never have a leading coefficient or unary operator like `"2x"` or `"-x"`.
#
# Expressions are evaluated in the usual order: brackets first, then
# multiplication, then addition and subtraction. For example, `expression = "1 +
# 2 * 3"` has an answer of `["7"]`.
#
# The format of the output is as follows:
#
#   * For each term of free variables with non-zero coefficient, we write the free variables within a term in sorted order lexicographically. For example, we would never write a term like `"b*a*c"`, only `"a*b*c"`.
#   * Terms have degree equal to the number of free variables being multiplied, counting multiplicity. (For example, `"a*a*b*c"` has degree 4.) We write the largest degree terms of our answer first, breaking ties by lexicographic order ignoring the leading coefficient of the term.
#   * The leading coefficient of the term is placed directly to the left with an asterisk separating it from the variables (if they exist.)  A leading coefficient of 1 is still printed.
#   * An example of a well formatted answer is `["-2*a*a*a", "3*a*a*b", "3*b*b", "4*a", "5*c", "-6"]`
#   * Terms (including constant terms) with coefficient 0 are not included.  For example, an expression of "0" has an output of [].
#
# **Examples:**
#
#
#
#     Input: expression = "e + 8 - a + 5", evalvars = ["e"], evalints = [1]
#     Output: ["-1*a","14"]
#
#     Input: expression = "e - 8 + temperature - pressure",
#     evalvars = ["e", "temperature"], evalints = [1, 12]
#     Output: ["-1*pressure","5"]
#
#     Input: expression = "(e + 8) * (e - 8)", evalvars = [], evalints = []
#     Output: ["1*e*e","-64"]
#
#     Input: expression = "7 - 7", evalvars = [], evalints = []
#     Output: []
#
#     Input: expression = "a * b * c + b * a * c * 4", evalvars = [], evalints = []
#     Output: ["5*a*b*c"]
#
#     Input: expression = "((a - b) * (b - c) + (c - a)) * ((a - b) + (b - c) * (c - a))",
#     evalvars = [], evalints = []
#     Output: ["-1*a*a*b*b","2*a*a*b*c","-1*a*a*c*c","1*a*b*b*b","-1*a*b*b*c","-1*a*b*c*c","1*a*c*c*c","-1*b*b*b*c","2*b*b*c*c","-1*b*c*c*c","2*a*a*b","-2*a*a*c","-2*a*b*b","2*a*c*c","1*b*b*b","-1*b*b*c","1*b*c*c","-1*c*c*c","-1*a*a","1*a*b","1*a*c","-1*b*c"]
#
#
# **Note:**
#
#   1. `expression` will have length in range `[1, 250]`.
#   2. `evalvars, evalints` will have equal lengths in range `[0, 100]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end