text:
"""
    The structure provides convenient computation of
    a Groebner basis of given ideal at a point.

    See evaluate function for more details
"""

### Probably we want to dispatch on coefficients type
mutable struct GroebnerEvaluator
    ideal::IdealContext
end

"""
    Convenience ctor 1

    ?? convenience
"""
function GroebnerEvaluator(ideal)
    GroebnerEvaluator(ideal)
end

#=
"""
    Convenience ctor 2
"""
function GroebnerEvaluator(
                           ideal,
                           eval_ring,
                           coeff_ring,
                           ground;
                           saturated=true)

    GroebnerEvaluator(ideal,
                       eval_ring,
                       coeff_ring,
                       ground,
                       saturated
    )
end
=#

"""
    Evaluates all coeffs of the underlying_ideal of G at the given point p,
    then computing the Groebner basis of the resulting ideal
    Returns the computed Groebner basis

    Groebner basis generators are designed to be the elements
    of G.eval_ring over G.ground

    Underlying ideal polynomial coefficients live in G.coeff_ring and
    must agree with substituting the given point p

    Dependent of the saturating variable "t" polynomials are
    erased from the resulting Groebner basis
"""
function AbstractAlgebra.evaluate(G::GroebnerEvaluator, p)
    context = G.ideal

    I = context.I
    ground = context.ground
    evalring = context.evalring

    singular_ground = base_ring(evalring)

    # TODO: change
    lift = typeof(p[1]) <: Nemo.gfp_elem ? x -> x.data : x -> x

    Is = [
            map_coefficients(c -> singular_ground(lift(evaluate(c, p))), f)
            for f in I
         ]

    Is = [
            change_base_ring(base_ring(evalring), f, parent=evalring)
            for f in Is
    ]

    ideal = Singular.Ideal(evalring, Is)
    # @info "I am gpoing to compute GB of $ideal"
    gb = GroebnerBasis.f4(ideal, reducegb=0, monorder=:lex)
    # this should never happen ideally (but it does!!)
    # @info gb
    if length(gens(gb)) == 1
        @info gb
        @warn "F4 failed. Switching to Singular.std"
        gb = Singular.std(ideal, complete_reduction=true)
    end

    gb = collect(gens(gb))

    # ideal = Singular.Ideal(eval_ring, gb)
    # gb = collect(gens(Singular.std(ideal, complete_reduction=true)))

    # @info "After reduction $gb"

    if context.saturated
        t = last(gens(evalring))
        gb = filter(f -> degree(f, t) == 0, gb)
    end

    gb = sort(gb, by=collect ∘ AbstractAlgebra.exponent_vectors)

    normalize(f) = !isconstant(f) ? f * inv(leading_coefficient(f)) : f
    gb = map(normalize, gb)

    gb
end

"""
    Short-cut for evaluate
"""
function (G::GroebnerEvaluator)(p)
    evaluate(G, p)
end

function AbstractAlgebra.nvars(G::GroebnerEvaluator)
    return length(gens(G.ideal.basepolyring))
end

"""
    Returns a random point suitable for evaluation
"""
function generate_point(G::GroebnerEvaluator; M=Inf)
    if M == Inf
        [ rand(G.ideal.ground) for _ in 1:nvars(G) ]
    else
        [ G.ideal.ground(rand(0:M)) for _ in 1:nvars(G) ]
    end
end

function AbstractAlgebra.base_ring(G::GroebnerEvaluator)
    G.ideal.ground
end


text:
using Documenter
import PoseComposition

pages_in_order = [
    "index.md",
    "operations.md",
    "further_api_reference.md",
]

makedocs(
    sitename="PoseComposition.jl",
    pages = pages_in_order,
    expandfirst = pages_in_order,
)

deploydocs(
    repo = "github.com/probcomp/PoseComposition.jl.git",
    devbranch = "main",
)

ext: jl


text:
println("**********Testing Distribution**********")

facts("Creating a Distribution") do
  fv = FeatureVector()
  c = Cluster()
  ds = DataSet()
  d1 = Distribution(fv)
  d2 = Distribution(c)
  d3 = Distribution(ds)

  @fact typeof(d1) => Distribution{FeatureVector}
  @fact typeof(d2) => Distribution{Cluster}
  @fact typeof(d3) => Distribution{DataSet}
end

facts("If sent in with a FeatureVector, d.total is set") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact d1.total => 7
end

facts("Get probability of seeing feature in Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  value = probability(d1,"word")

  @fact value => 4/7
end

facts("Get features of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact features(d1) => keys(d1.space)
end

facts("Check isempty on a Distribution") do
  d1 = Distribution{FeatureVector}()
  d2 = Distribution{Cluster}()
  d3 = Distribution{DataSet}()

  @fact isempty(d1) => isempty(d1.space)
  @fact isempty(d2) => isempty(d2.space)
  @fact isempty(d3) => isempty(d3.space)
end

facts("entropy(Distribution) returns 0 if empty Distribution") do
  d1 = Distribution{FeatureVector}()

  @fact entropy(d1) => 0
end

facts("entropy(Distribution) returns entropy of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  @fact entropy(d1) => 0.9852281360342515
end

facts("entropy of different Distribution types") do
  fv1 = FeatureVector(["c"=>1,"b"=>1,"a"=>1,"d"=>1])
  fv2 = FeatureVector(["f"=>1,"g"=>1,"e"=>1,"h"=>1])
  fva = FeatureVector(["1"=>1,"2"=>1,"3"=>1,"4"=>1])
  fvb = FeatureVector(["8"=>1,"7"=>1,"6"=>1,"5"=>1])

  c1 = Cluster()
  c2 = Cluster()
  c1["fv1"] = fv1
  c1["fv2"] = fv2
  c2["fva"] = fva
  c2["fvb"] = fvb

  ds = DataSet()
  ds["c1"] = c1
  ds["c2"] = c2

  fv1_dist = Distribution(fv1)
  fv2_dist = Distribution(fv2)
  fva_dist = Distribution(fva)
  fvb_dist = Distribution(fvb)
  c1_dist = Distribution(c1)
  c2_dist = Distribution(c2)
  ds_dist = Distribution(ds)

  @fact entropy(fv1_dist) => 2
  @fact entropy(fv2_dist) => 2
  @fact entropy(fva_dist) => 2
  @fact entropy(fvb_dist) => 2
  @fact entropy(c1_dist) => 3
  @fact entropy(c2_dist) => 3
  @fact entropy(ds_dist) => 4
end

facts("info_gain(Distribution,Distribution) returns 0 if empty Distributions") do
  d1 = Distribution{FeatureVector}()
  d2 = Distribution{FeatureVector}()

  @fact info_gain(d1,d2) => 0
end

facts("info_gain(Distribution,Distribution) returns info_gain of Distribution") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  dict2 = ["∂" => .1, "happy" => .3]
  fv2 = FeatureVector(dict2)
  d2 = Distribution(fv2)

  @fact info_gain(d1,d2) => 0.17395001157511858
end

facts("Smoothing FeatureVector has correct probabilties before and after") do
  dict1 = ["word" => 7, "another" => 13]
  fv1 = FeatureVector(dict1)
  d1 = Distribution(fv1)

  #no smoothing
  @fact probability(d1,"word") => (7/20)
  @fact probability(d1,"another") => (13/20)
  @fact probability(d1,"unk") => (0/20)

  #delta smoothing
  delta_smoothing!(d1)
  @fact probability(d1,"word") => (8/(3+20))
  @fact probability(d1,"another") => (14/(3+20))
  @fact probability(d1,"unk") => (1/(3+20))

  #removing smoothing
  remove_smoothing!(d1)
  @fact probability(d1,"word") => (7/20)
  @fact probability(d1,"another") => (13/20)
  @fact probability(d1,"unk") => (0/20)
end

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
function get_expected_buses(::Type{T}, sys::System) where {T <: AggregationTopology}
    expected_buses = Dict{String, Vector{String}}()
    for bus in get_components(Bus, sys)
        agg = get_aggregation_topology_accessor(T)(bus)
        name = get_name(agg)
        if !haskey(expected_buses, name)
            expected_buses[name] = Vector{String}()
        end
        push!(expected_buses[name], get_name(bus))
    end

    return expected_buses
end

function test_aggregation_topologies(sys::System, expected_areas, expected_zones)
    expected_buses_by_area = get_expected_buses(Area, sys)
    expected_buses_by_zone = get_expected_buses(LoadZone, sys)

    areas = collect(get_components(Area, sys))
    @test length(areas) == expected_areas
    area_mapping = get_aggregation_topology_mapping(Area, sys)
    for area in areas
        area_name = get_name(area)
        buses = sort!([get_name(x) for x in get_buses(sys, area)])
        @test buses == sort(expected_buses_by_area[area_name])
        @test buses == sort!([get_name(x) for x in area_mapping[area_name]])
    end

    zones = collect(get_components(LoadZone, sys))
    zone_mapping = get_aggregation_topology_mapping(LoadZone, sys)
    @test length(zones) == expected_zones
    for zone in zones
        zone_name = get_name(zone)
        buses = sort!([get_name(x) for x in get_buses(sys, zone)])
        @test buses == sort(expected_buses_by_zone[zone_name])
        @test buses == sort!([get_name(x) for x in zone_mapping[zone_name]])
    end
end

@testset "Test topology mappings" begin
    sys = create_rts_system()
    test_aggregation_topologies(sys, 3, 21)
end

@testset "Test PM areas and load zones" begin
    path = joinpath(DATA_DIR, "matpower", "RTS_GMLC.m")
    sys = System(PowerSystems.PowerModelsData(path))
    test_aggregation_topologies(sys, 3, 21)
end

@testset "Test get_components_in_aggregation_topology" begin
    sys = create_rts_system()
    areas = collect(get_components(Area, sys))
    @test !isempty(areas)
    area = areas[1]
    generators = get_components_in_aggregation_topology(ThermalStandard, sys, area)

    for gen in generators
        bus = get_bus(gen)
        @test IS.get_uuid(get_area(bus)) == IS.get_uuid(area)
    end
end


text:
using GPUArrays
using Base.Test, GPUArrays.TestSuite

function run_broadcasting(Typ)
    @testset "broadcast" begin
        test_broadcast(Typ)
        test_vec3(Typ)
    end
end

test_idx(idx, A::AbstractArray{T}) where T = A[idx] * T(2)

function testv3_1(a, b)
    x = a .* b
    y =  x .+ b
    @inbounds return y[1]
end

function testv3_2(a, b)
    x = a .* b
    y =  x .+ b
    return y
end

# more complex function for broadcast
function test_kernel(a::T, b) where T
    c = a
    for i=1:5
        c += T(b)
        c = a * T(2)
    end
    return c
end

function test_broadcast(Typ)
    for ET in supported_eltypes()
        N = 10
        T = Typ{ET}
        @testset "broadcast $ET" begin
            @testset "RefValue" begin
                cidx = rand(UInt32(1):UInt32(N), 2*N)
                gidx = Typ(cidx)
                cy = TestSuite.toarray(ET, (2*N,))
                gy = Typ(cy)
                cres = zeros(ET, size(cidx))
                gres = Typ(cres)
                gres .= test_idx.(gidx, Base.RefValue(gy))
                cres .= test_idx.(cidx, Base.RefValue(cy))
                @test Array(gres) == cres
            end
            @testset "Tuple" begin
                against_base(T, (3, N), (3, N), (N,), (N,), (N,)) do out, arr, a, b, c
                    res2 = broadcast!(out, arr, (a, b, c)) do xx, yy
                        xx + sum(yy)
                    end
                end
            end
            ############
            # issue #27
            against_base((a, b)-> a .+ b, T, (4, 5, 3), (1, 5, 3))
            against_base((a, b)-> a .+ b, T, (4, 5, 3), (1, 5, 1))

            ############
            # issue #22
            dim = (32, 32)
            against_base(T, dim, dim, dim) do tmp, a1, a2
                tmp .=  a1 .+ a2 .* ET(2)
            end

            ############
            # issue #21
            if ET in (Float32, Float64)
                against_base((a1, a2)-> muladd.(ET(2), a1, a2), T, dim, dim)
                #########
                # issue #41
                # The first issue is likely https://github.com/JuliaLang/julia/issues/22255
                # since GPUArrays adds some arguments to the function, it becomes longer longer, hitting the 12
                # so this wont fix for now
                against_base(T, dim, dim, dim, dim, dim, dim) do a1, a2, a3, a4, a5, a6
                    @. a1 = a2 + (1.2) *((1.3)*a3 + (1.4)*a4 + (1.5)*a5 + (1.6)*a6)
                end

                against_base(T, dim, dim, dim, dim) do u, uprev, duprev, ku
                    fract = ET(1//2)
                    dt = ET(1.4)
                    dt2 = dt^2
                    @. u = uprev + dt*duprev + dt2*(fract*ku)
                end
                against_base((x) -> (-).(x), T, (2, 3))

                against_base(T, dim, dim, dim, dim, dim, dim) do utilde, gA, k1, k2, k3, k4
                    btilde1 = ET(1)
                    btilde2 = ET(1)
                    btilde3 = ET(1)
                    btilde4 = ET(1)
                    dt = ET(1)
                    @. utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
                end
            end

            against_base((x) -> fill!(x, 1), T, (3,3))
            against_base((x, y) -> map(+, x, y), T, (2, 3), (2, 3))

            against_base((x) -> 2x, T, (2, 3))
            against_base((x, y) -> x .+ y, T, (2, 3), (1, 3))
            against_base((z, x, y) -> z .= x .+ y, T, (2, 3), (2, 3), (2,))

            T = Typ{ET}
            against_base(A -> A .= identity.(ET(10)), T, (40, 40))
            against_base(A -> test_kernel.(A, ET(10)), T, (40, 40))
            against_base(A -> A .* ET(10), T, (40, 40))
            against_base((A, B) -> A .* B, T, (40, 40), (40, 40))
            against_base((A, B) -> A .* B .+ ET(10), T, (40, 40), (40, 40))
        end
    end
end

function test_vec3(Typ)
    @testset "vec 3" begin
        N = 20

        xc = map(x-> ntuple(i-> rand(Float32), Val{3}), 1:N)
        yc = map(x-> ntuple(i-> rand(Float32), Val{3}), 1:N)

        x = Typ(xc)
        y = Typ(yc)

        res1c = zeros(Float32, N)
        res2c = similar(xc)

        res1 = Typ(res1c)
        res2 = Typ(res2c)

        res1 .= testv3_1.(x, y)
        res1c .= testv3_1.(xc, yc)
        @test Array(res1) ≈ res1c

        res2 .= testv3_2.(x, y)
        res2c .= testv3_2.(xc, yc)
        @test all(map((a,b)-> all((1,2,3) .≈ (1,2,3)), Array(res2), res2c))
    end
end


text:
module REISE

using CSV: CSV
using DataFrames: DataFrames
using Dates: Dates
using JuMP: JuMP
using MAT: MAT
using Requires: Requires
import SparseArrays: sparse, SparseMatrixCSC

include("types.jl")         # Defines Case, Results, Storage, DemandFlexibility,
#     VariablesOfInterest
include("read.jl")          # Defines read_case, read_storage, read_demand_flexibility
include("prepare.jl")       # Defines reise_data_mods
include("model.jl")         # Defines _build_model (used in interval_loop)
include("loop.jl")          # Defines interval_loop
include("query.jl")         # Defines get_results (used in interval_loop)
include("save.jl")          # Defines save_input_mat, save_results

function __init__()
    Requires.@require Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" begin
        include(joinpath("solver_specific", "gurobi.jl"))
    end
end

"""
    REISE.run_scenario(;
        interval=24, n_interval=3, start_index=1, outputfolder="output",
        inputfolder=pwd())

Run a scenario consisting of several intervals.
'interval' specifies the length of each interval in hours.
'n_interval' specifies the number of intervals in a scenario.
'start_index' specifies the starting hour of the first interval, to determine
    which time-series data should be loaded into each intervals.
'inputfolder' specifies where to load the relevant data from. Required files
    are 'case.mat', 'demand.csv', 'hydro.csv', 'solar.csv', and 'wind.csv'.
'outputfolder' specifies where to store the results. Defaults to an `output`
    subdirectory of inputfolder. This folder will be created if it does not exist at
    runtime.
'optimizer_factory' is the solver used for optimization. If not specified, Gurobi is
    used by default.
"""
function run_scenario(;
    num_segments::Int=1,
    interval::Int,
    n_interval::Int,
    start_index::Int,
    inputfolder::String,
    outputfolder::Union{String,Nothing}=nothing,
    threads::Union{Int,Nothing}=nothing,
    optimizer_factory=nothing,
    solver_kwargs::Union{Dict,Nothing}=nothing,
    model_kwargs::Union{Dict,Nothing}=nothing,
)
    isnothing(optimizer_factory) && error("optimizer_factory must be specified")
    # Setup things that build once
    # If no solver kwargs passed, instantiate an empty dict
    solver_kwargs = something(solver_kwargs, Dict())
    # If outputfolder not given, by default assign it inside inputfolder
    isnothing(outputfolder) && (outputfolder = joinpath(inputfolder, "output"))
    # If outputfolder doesn't exist (isdir evaluates false) create it (mkdir)
    isdir(outputfolder) || mkdir(outputfolder)
    stdout_filepath = joinpath(outputfolder, "stdout.log")
    stderr_filepath = joinpath(outputfolder, "stderr.err")
    case = read_case(inputfolder)
    storage = read_storage(inputfolder)
    demand_flexibility = read_demand_flexibility(inputfolder, interval)
    println("All scenario files loaded!")
    case = reise_data_mods(case; num_segments=num_segments)
    save_input_mat(case, storage, inputfolder, outputfolder)
    # Create final model kwargs by merging defaults with user-specified
    default_model_kwargs = Dict(
        "case" => case,
        "storage" => storage,
        "demand_flexibility" => demand_flexibility,
        "interval_length" => interval,
    )
    if isnothing(model_kwargs)
        model_kwargs = default_model_kwargs
    else
        model_kwargs = merge(default_model_kwargs, model_kwargs)
    end
    # If a number of threads is specified, add to solver settings dict
    isnothing(threads) || (solver_kwargs["Threads"] = threads)
    println("All preparation complete!")
    # While redirecting stdout and stderr...
    println("Redirecting outputs, see stdout.log & stderr.err in outputfolder")
    redirect_stdout_stderr(stdout_filepath, stderr_filepath) do
        # Loop through intervals
        m = interval_loop(
            optimizer_factory,
            model_kwargs,
            solver_kwargs,
            interval,
            n_interval,
            start_index,
            inputfolder,
            outputfolder,
        )
    end
    return m
end

# Module end
end


text:
## gcmfaces type definition + methods

## type definition

abstract type AbstractGcmfaces{T, N} <: AbstractArray{T, N} end

struct gcmfaces{T, N} <: AbstractGcmfaces{T, N}
   nFaces::Int
   grTopo::String
   f::Array{Array{T,N},1}
end

# additional constructors

function gcmfaces(nFaces::Int,grTopo::String,v1::AbstractArray)
  tmp_eltype=eltype(v1[1])
  tmp_ndims=ndims(v1[1])
  gcmfaces{tmp_eltype,tmp_ndims}(nFaces,grTopo,v1)
end

function gcmfaces(nFaces::Int,grTopo::String)
  tmp0=Array{Float64}(undef,1,1)
  tmp1=fill(tmp0,nFaces)
  gcmfaces{Float64,2}(nFaces,grTopo,tmp1)
end

gcmfaces() = gcmfaces(MeshArrays.nFaces,MeshArrays.grTopo);

## Convenience functions

function fsize(A::AbstractGcmfaces{T, N}) where {T,N}
  fs=Array{NTuple{N, Int}}(undef,A.nFaces)
  for i=1:A.nFaces
    fs[i]=size(A.f[i]);
  end
  return fs
end

function fsize(A::AbstractGcmfaces{T, N},i::Integer) where {T,N}
  fs=size(A.f[i]);
end

## Interface Methods

function Base.size(A::AbstractGcmfaces{T, N}) where {T,N}
  fs=fsize(A)
  tmp1=0
  for i=1:A.nFaces
    tmp1=tmp1+fs[i][1]*fs[i][2]
  end
  #N>2 ? s=(tmp1,1,fs[1][3]) : s=(tmp1,1)
  if N==2; s=(tmp1,1);
  elseif N==3; s=(tmp1,1,fs[1][3]);
  elseif N==4; s=(tmp1,1,fs[1][3],fs[1][4]);
  else; error("N>4 case not implemented yet");
  end
  return s
end

function Base.size(A::AbstractGcmfaces{T, N},dim::Integer) where {T,N}
  tmp1=size(A)
  s=tmp1[dim]
end

function Base.getindex(A::AbstractGcmfaces{T, N}, i::Vararg{Int, N}) where {T,N}

  f=0
  j=0
  k=0
  tmp1=0
  for iFace=1:A.nFaces
    tmpsize=fsize(A,iFace)
    tmp11=tmpsize[1]*tmpsize[2]
    tmp2=tmp1+tmp11
    if tmp1<i[1]<=tmp2
      f=iFace;
      tmp3=(i[1]-tmp1);
      k=Int(ceil(tmp3/tmpsize[1]))
      j=Int(tmp3-tmpsize[1]*(k-1))
    end;
    tmp1=tmp1+tmp11
  end
  if ndims(A)==2; val=A.f[f][j,k];
  elseif ndims(A)==3; val=A.f[f][j,k,i[3]];
  elseif ndims(A)==4; val=A.f[f][j,k,i[3],i[4]];
  else; error("ndims>4 case not implemented yet")
  end
  return val
end

## view

function Base.view(a::AbstractGcmfaces{T, N}, I::Vararg{Any, N}) where {T,N}
  nFaces=a.nFaces;
  grTopo=a.grTopo;
  v1=Array{Any}(undef,nFaces);
  for iFace=1:nFaces
    if N==2;
      v1[iFace]=view(a.f[iFace],:,:);
    elseif N==3;
      v1[iFace]=view(a.f[iFace],:,:,I[3]);
    elseif N==4;
      v1[iFace]=view(a.f[iFace],:,:,I[3],I[4]);
    else;
      error("N>4 case not implemented yet")
    end
  end
  c=gcmfaces(nFaces,grTopo,v1);
  return c;
end

# Custom pretty-printing

function Base.show(io::IO, z::gcmfaces)

#    @printf io " MeshArrays instance with \n"
    printstyled(io, " gcmfaces array \n",color=:normal)
    printstyled(io, "  grid type   = ",color=:normal)
    printstyled(io, "$(z.grTopo)\n",color=:blue)
    printstyled(io, "  # of faces  = ",color=:normal)
    printstyled(io, "$(z.nFaces)\n",color=:blue)
    if ~isassigned(z.f);
      printstyled(io, "  data type   = ",color=:normal)
      printstyled(io, "unassigned\n",color=:green)
      printstyled(io, "  face sizes  = ",color=:normal)
      printstyled(io, "unassigned\n",color=:green)
    else
      printstyled(io, "  data type   = ",color=:normal)
      printstyled(io, "$(typeof(z.f[1][1]))\n",color=:blue)
      printstyled(io, "  face sizes  = ",color=:normal)
      printstyled(io, "$(size(z.f[1]))\n",color=:blue)
      for iFace=2:z.nFaces
        printstyled(io, "                ",color=:normal)
        printstyled(io, "$(size(z.f[iFace]))\n",color=:blue)
      end
    end

    return
end

# basic operations

import Base: +, -, *, /, getindex
import Base: isnan, isinf, isfinite
import Base: maximum, minimum, sum, fill

function +(a::gcmfaces)
  c=gcmfaces(a.nFaces,a.grTopo)
  for iFace=1:a.nFaces
    c.f[iFace]=+a.f[iFace];
  end
  return c
end

function +(a::gcmfaces,b::gcmfaces)
  c=gcmfaces(a.nFaces,a.grTopo)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]+b.f[iFace];
  end
  return c
  #the following modifies a:
  #  c=a
  #  c.f=a.f .+ b.f
  #the following fails immutability:
  #  c=gcmfaces(a.nFaces,a.grTopo)
  #  c.f=a.f .+ b.f
end

function +(a::Number,b::gcmfaces)
  c=gcmfaces(b.nFaces,b.grTopo)
  for iFace=1:b.nFaces
    c.f[iFace]=a.+b.f[iFace];
  end
  return c
  #the following is deprecated synthax as of v0.7:
  #  c=gcmfaces(b.nFaces,b.grTopo)
  #  c.f=a .+ b.f
end

function +(a::gcmfaces,b::Number)
  c=b+a
  return c
end

function -(a::gcmfaces)
  c=gcmfaces(a.nFaces,a.grTopo)
  for iFace=1:a.nFaces
    c.f[iFace]=-a.f[iFace];
  end
  return c
end

function -(a::gcmfaces,b::gcmfaces)
  c=gcmfaces(a.nFaces,a.grTopo)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]-b.f[iFace];
  end
  return c
end

function -(a::Number,b::gcmfaces)
  c=gcmfaces(b.nFaces,b.grTopo)
  for iFace=1:b.nFaces
    c.f[iFace]=a.-b.f[iFace];
  end
  return c
end

function -(a::gcmfaces,b::Number)
  c=gcmfaces(a.nFaces,a.grTopo)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace].-b;
  end
  return c
end

function *(a::gcmfaces,b::gcmfaces)
  nFaces=a.nFaces;
  grTopo=a.grTopo;
  c=gcmfaces(nFaces,grTopo);
  for iFace=1:nFaces
    c.f[iFace]=a.f[iFace].*b.f[iFace];
  end
  return c
end

function *(a::Number,b::gcmfaces)
  c=gcmfaces(b.nFaces,b.grTopo);
  for iFace=1:b.nFaces
    c.f[iFace]=a*b.f[iFace];
  end
  return c
end

function *(a::gcmfaces,b::Number)
  c=b*a
  return c
end

function /(a::gcmfaces,b::gcmfaces)
  nFaces=a.nFaces;
  grTopo=a.grTopo;
  c=gcmfaces(nFaces,grTopo);
  for iFace=1:nFaces
    c.f[iFace]=a.f[iFace]./b.f[iFace];
  end
  return c
end

function /(a::Number,b::gcmfaces)
  nFaces=b.nFaces;
  grTopo=b.grTopo;
  c=gcmfaces(nFaces,grTopo);
  for iFace=1:nFaces
    c.f[iFace]=a./b.f[iFace];
  end
  return c
end

function /(a::gcmfaces,b::Number)
  nFaces=a.nFaces;
  grTopo=a.grTopo;
  c=gcmfaces(nFaces,grTopo);
  for iFace=1:nFaces
    c.f[iFace]=a.f[iFace]./b;
  end
  return c
end

##

function isnan(a::gcmfaces)
    nFaces=a.nFaces;
    grTopo=a.grTopo;
    c=gcmfaces(a.nFaces,a.grTopo);
    for iFace=1:a.nFaces
      c.f[iFace]=isnan.(a.f[iFace]);
    end
    return c
end

function isinf(a::gcmfaces)
    nFaces=a.nFaces;
    grTopo=a.grTopo;
    c=gcmfaces(a.nFaces,a.grTopo);
    for iFace=1:a.nFaces
      c.f[iFace]=isinf.(a.f[iFace]);
    end
    return c
end

function isfinite(a::gcmfaces)
    nFaces=a.nFaces;
    grTopo=a.grTopo;
    c=gcmfaces(a.nFaces,a.grTopo);
    for iFace=1:a.nFaces
      c.f[iFace]=isfinite.(a.f[iFace]);
    end
    return c
end

function sum(a::gcmfaces)
    c=0.0;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=c+sum(tmp1);
    end
    return c
end

function maximum(a::gcmfaces)
    c=-Inf;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=max(c,maximum(tmp1));
    end
    return c
end

function minimum(a::gcmfaces)
    c=Inf;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=min(c,minimum(tmp1));
    end
    return c
end

function fill(val::Any,a::gcmfaces)
    c=gcmfaces(a.nFaces,a.grTopo);
    for iFace=1:a.nFaces
      c.f[iFace]=fill(val,fsize(a,iFace));
    end
    return c
end


text:
"""
Non-negative CANDECOMP tensor decomposition.
"""
function nncp(tnsr::AbstractArray{T,N},
              r::Integer;
              tol::Float64=1e-4,
              maxiter::Integer=100,
              compute_error::Bool=false,
              verbose::Bool=true) where {T,N}

    all(x -> x >= 0, tnsr) || throw(ArgumentError("Input tensor must be nonnegative."))
    num_modes = _check_tensor(tnsr, r)
    T_norm = norm(tnsr)

    factors = [abs.(F) .* (T_norm ^ (1/num_modes) / norm(F)) for F::Matrix{T} in _random_factors(size(tnsr), r)]
    factors_old = deepcopy(factors)
    factors_exp = deepcopy(factors)
    gram = [F'F for F in factors]

    niters = 0
    converged = false
    obj_old = 0.0
    t = 1.0
    t_old = 1.0

    LB = fill!(similar(tnsr, num_modes), 1.0)
    LB_old = fill!(similar(tnsr, num_modes), 1.0)

    U = similar(tnsr, r, r)

    pb = Progress(maxiter, "NNCP iterations ")
    while !converged && niters < maxiter
        update!(pb, niters)
        LB_old .= LB
        local M::typeof(U)
        for i in 1:num_modes
            idx = [num_modes:-1:i+1; i-1:-1:1] # num_modes except i
            U .= reduce((x, y) -> x .* y, gram[idx]) # element-wise product of gram matrices
            LB[i] = norm(U)
            M = _row_unfold(tnsr, i) * khatrirao(factors[idx])
            factors[i] = max.(0, factors_exp[i] .- (factors_exp[i] * U .- M) .* (1 / LB[i]))
            mul!(gram[i], factors[i]', factors[i])
        end

        obj = dot(gram[num_modes], U) - 4 * dot(factors[num_modes], M)
        if obj > obj_old
            for i in eachindex(factors_old)
                copyto!(factors_exp[i], factors_old[i])
            end
        else
            t = 1 + sqrt(1 + 4 * t_old^2)
            min_w = (t_old - 2) / t
            for i in eachindex(factors)
                w = min(min_w, sqrt(LB_old[i] / LB[i]))
                factors_exp[i] .= factors[i] .+ w .* (factors[i] .- factors_old[i])
                copyto!(factors_old[i], factors[i])
            end
            converged = abs(obj_old - obj) < tol * abs(obj_old)
            t_old = t
            obj_old = obj
        end
        niters += 1
    end
    finish!(pb)
    verbose && _iter_status(converged, niters, maxiter)

    res = CANDECOMP(tuple(factors...), fill!(similar(tnsr, r), 1.0))
    compute_error && _set_rel_residue(res, tnsr)
    return res
end


text:
# safe log(x) function, log 0 = 0
safelog(x) = x > zero(x) ? log2(x) : zero(x)
safediv(x, y) = y == zero(y) ? zero(y) : x / y
# general functions for entropy

"""
    make_joint_distribution(N::NT) where {NT<:AbstractEcologicalNetwork}

Returns a double stochastic matrix from the adjacency or incidence matrix.
Raises an error if the matrix contains negative values. Output in bits.
"""
function make_joint_distribution(N::NT) where {NT<:AbstractEcologicalNetwork}
    !(typeof(N) <: QuantitativeNetwork) || any(N.A .≥ 0) || throw(DomainError("Information only for nonnegative interaction values"))
    return N.A / sum(N.A)
end

"""
   entropy(P::AbstractArray)

Computes the joint entropy of a double stochastic matrix. Does not perform any
checks whether the matrix is normalized. Output in bits.
"""
function entropy(P::AbstractArray)
    return - sum(P .* safelog.(P))
end

"""
   entropy(P::AbstractArray, dims::I)

Computes the marginal entropy of a double stochastic matrix. `dims` indicates
whether to compute the entropy for the rows (`dims`=1) or columns (`dims`=2).
Does not perform any checks whether the matrix is normalized. Output in bits.
"""
function entropy(P::AbstractArray, dims::I) where I <: Int
    return entropy(sum(P', dims=dims))
end

"""
   entropy(N::AbstractEcologicalNetwork)

Computes the joint entropy of an ecological network. Output in bits.
"""
function entropy(N::NT) where {NT<:AbstractEcologicalNetwork}
    P = make_joint_distribution(N)
    return entropy(P)
end

"""
   entropy(N::AbstractEcologicalNetwork, dims::I)

Computes the marginal entropy of an ecological network. `dims` indicates
whether to compute the entropy for the rows (`dims`=1) or columns (`dims`=2).
Output in bits.
"""
function entropy(N::NT, dims::I) where {NT<:AbstractEcologicalNetwork, I <: Int}
    P = make_joint_distribution(N)
    return entropy(P, dims)
end

"""
   conditional_entropy(P::AbstractArray, given::I)

Computes the conditional entropy of double stochastic matrix. If `given = 1`,
it is the entropy of the columns, and visa versa when `given = 2`. Output in bits.
"""
function conditional_entropy(P::AbstractArray, given::I) where I <: Int
    dims = (given % 2) + 1
    return - sum(P .* safelog.(safediv.(P, sum(P, dims=dims))))
end

"""
   conditional_entropy(N::AbstractEcologicalNetwork, given::I)

Computes the conditional entropy of an ecological network. If `given = 1`,
it is the entropy of the columns, and visa versa when `given = 2`.
"""
function conditional_entropy(N::NT, given::I) where {NT<:AbstractEcologicalNetwork, I <: Int}
    P = make_joint_distribution(N)
    return conditional_entropy(P, given)
end

"""
    mutual_information(P::AbstractArray)

Computes the mutual information of a double stochastic matrix. Output in bits.
"""
function mutual_information(P::AbstractArray)
    I = entropy(P, 1) - conditional_entropy(P, 2)
    return max(I, zero(I))  # might give small negative value
end

"""
    mutual_information(P::AbstractArray)

Computes the mutual information of an ecological network. Output in bits.
"""
function mutual_information(N::NT) where {NT<:AbstractEcologicalNetwork}
    P = make_joint_distribution(N)
    return mutual_information(P)
end

"""
    variation_information(P::AbstractArray)

Computes the variation of information of a double stochastic matrix. Output in
bits.
"""
function variation_information(P::AbstractArray)
    return conditional_entropy(P, 1) + conditional_entropy(P, 2)
end

"""
    variation_information(N::AbstractEcologicalNetwork)

Computes the variation of information of an ecological network. Output in
bits.
"""
function variation_information(N::NT) where {NT <: AbstractEcologicalNetwork}
    P = make_joint_distribution(N)
    return variation_information(P)
end

"""
    diff_entropy_uniform(P::AbstractArray)

Computes the difference in entropy of the marginals compared to the entropy of
an uniform distribution. The parameter `dims` indicates which marginals are used,
with both if no value is provided. Output in bits.
"""
function diff_entropy_uniform(P::AbstractArray)
    D = log2(length(P)) - entropy(P, 1) - entropy(P, 2)
    return max(zero(D), D)
end

"""
    diff_entropy_uniform(P::AbstractArray, dims::I)

Computes the difference in entropy of the marginals compared to the entropy of
an uniform distribution. The parameter `dims` indicates which marginals are used,
with both if no value is provided. Output in bits.
"""
function diff_entropy_uniform(P::AbstractArray, dims::I) where {I <: Int}
    D = log2(size(P, dims)) - entropy(P, dims)
    return max(zero(D), D)
end

"""
    diff_entropy_uniform(N::AbstractEcologicalNetwork, dims::I=nothing)

Computes the difference in entropy of the marginals compared to the entropy of
an uniform distribution. The parameter `dims` indicates which marginals are used,
with both if no value is provided. Output in bits.
"""
function diff_entropy_uniform(N::NT, dims=nothing) where {NT <: AbstractEcologicalNetwork}
    if dims == nothing
        P = make_joint_distribution(N)
        return diff_entropy_uniform(P)
    else  # compute marginals
        typeof(dims) <: Int || throw(ArgumentError("dims should be an integer (1 or 2)"))
        P = make_joint_distribution(N)
        return diff_entropy_uniform(P, dims)
    end
end

"""
    information_decomposition(N::AbstractEcologicalNetwork; norm::Bool=false, dims::I=nothing)

Performs an information theory decomposition of a given ecological network, i.e.
the information content in the normalized adjacency matrix is split in:

- `:D` : difference in entropy of marginals compared to an uniform distribition
- `:I` : mutual information
- `:V` : variation of information / conditional entropy

If `norm=true`, the components are normalized such that their sum is equal to 1.
One can optinally give the dimision, indicating whether to compute the indices
for the rows (`dims=1`), columns (`dims=2`) or the whole matrix (default).

Result is returned in a Dict. Outputs in bits.
"""
function information_decomposition(N::NT; norm::Bool=false, dims::I=nothing) where {NT <: AbstractEcologicalNetwork, I <: Union{Int, Nothing}}
    decomposition = Dict{Symbol, Float64}()
    P = make_joint_distribution(N)
    if dims == nothing
        # difference marginal entropy
        decomposition[:D] = diff_entropy_uniform(P)
        # mutual information
        decomposition[:I] = 2mutual_information(P)
        # variance of information
        decomposition[:V] = variation_information(P)
    else
        decomposition[:D] = diff_entropy_uniform(P, dims)
        decomposition[:I] = mutual_information(P)
        decomposition[:V] = conditional_entropy(P, (dims % 2) + 1)
    end
    if norm
        potential_information = sum(values(decomposition))
        for key in keys(decomposition)
            decomposition[key] /= potential_information
        end
    end
    return decomposition
end

"""
    convert2effective(indice::Real)

Convert an information theory indices in an effective number (i.e. number of
corresponding interactions).
"""
function convert2effective(indice::R) where R <: Real
    return 2.0^indice
end

"""
    potential_information(N::NT)

Computes the maximal potential information in a network, corresponding to
every species interacting with every other species. Compute result for the
marginals using the optional parameter `dims`. Output in bits.
"""
function potential_information(N::NT) where NT <: AbstractEcologicalNetwork
    m, n = size(N)
    return log2(m) + log2(n)
end

"""
    potential_information(N::NT, dims::I)

Computes the maximal potential information in a network, corresponding to
every species interacting with every other species. Compute result for the
marginals using the optional parameter `dims`. Output in bits.
"""
function potential_information(N::NT, dims::I) where {NT <: AbstractEcologicalNetwork, I <: Int}
    return log2(size(N)[dims])
end


text:
# Choose number of agents
num_agents = 7

# Define the world polygon (counter-clockwise)
polyOutersp = [
             0.0 0.0;
             3.0 0.0;
             3.5 1.2;
             3.5 2.0;
             1.5 3.0;
            -0.5 2.0;
            -0.5 1.2;
            ]'

# Define the obstacle (clockwise)
polyInnersp = [
              0.7 1.0;
              0.8 1.5;
              1.5 2.0;
              2.2 1.5;
              2.3 1.0;
            ]'

poly = [polyOutersp, polyInnersp]

# Define the virtual gates
sp = [
     -0.5  1.5;
      1.5  3.0;
      2.2  1.5;
      2.3  1.0;
      1.2  1.0;
     ]'

ep = [
      0.8  1.5;
      1.5  2.0;
      3.5  2.0;
      3.0  0.0;
      1.2  0.0;
     ]'
clr = [
        :blue,
        :red,
        :green,
        :red,
        :green,
      ]

# Define initialization region
init_center = [0.2, 0.7]
init_width = 0.8

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics


text:
#!/usr/bin/env julia

# Import libraries.
using Turing, StatsPlots, Random

# Set the true probability of heads in a coin.
p_true = 0.5

# Iterate from having seen 0 observations to 100 observations.
Ns = 0:100

# Draw data from a Bernoulli distribution, i.e. draw heads or tails.
Random.seed!(12)
data = rand(Bernoulli(p_true), last(Ns))

# Declare our Turing model.
@model function coinflip(y)
    # Our prior belief about the probability of heads in a coin.
    p ~ Beta(1, 1)

    # The number of observations.
    N = length(y)
    for n in 1:N
        # Heads or tails of a coin are drawn from a Bernoulli distribution.
        y[n] ~ Bernoulli(p)
    end
end

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
iterations = 1000
ϵ = 0.05
τ = 10

# Start sampling.
#@time chain = sample(coinflip(data), HMC(ϵ, τ), iterations)
@time chain = sample(coinflip(data), HMC(ϵ, τ),MCMCThreads(),iterations,1)

# Plot a summary of the sampling process for the parameter p, i.e. the probability of heads in a coin.
#histogram(chain[:p])

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
# ---
# title: 563. Binary Tree Tilt
# id: problem563
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/binary-tree-tilt/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the sum of every tree node 's
# **tilt**._
#
# The **tilt** of a tree node is the **absolute difference** between the sum of
# all left subtree node **values** and all right subtree node **values**. If a
# node does not have a left child, then the sum of the left subtree node
# **values** is treated as `0`. The rule is similar if there the node does not
# have a right child.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/20/tilt1.jpg)
#
#
#
#     Input: root = [1,2,3]
#     Output: 1
#     Explanation:
#     Tilt of node 2 : |0-0| = 0 (no children)
#     Tilt of node 3 : |0-0| = 0 (no children)
#     Tile of node 1 : |2-3| = 1 (left subtree is just left child, so sum is 2; right subtree is just right child, so sum is 3)
#     Sum of every tilt : 0 + 0 + 1 = 1
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/20/tilt2.jpg)
#
#
#
#     Input: root = [4,2,9,3,5,null,7]
#     Output: 15
#     Explanation:
#     Tilt of node 3 : |0-0| = 0 (no children)
#     Tilt of node 5 : |0-0| = 0 (no children)
#     Tilt of node 7 : |0-0| = 0 (no children)
#     Tilt of node 2 : |3-5| = 2 (left subtree is just left child, so sum is 3; right subtree is just right child, so sum is 5)
#     Tilt of node 9 : |0-7| = 7 (no left child, so sum is 0; right subtree is just right child, so sum is 7)
#     Tilt of node 4 : |(3+5+2)-(9+7)| = |10-16| = 6 (left subtree values are 3, 5, and 2, which sums to 10; right subtree values are 9 and 7, which sums to 16)
#     Sum of every tilt : 0 + 0 + 0 + 2 + 7 + 6 = 15
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/20/tilt3.jpg)
#
#
#
#     Input: root = [21,7,14,1,1,2,2,3,3]
#     Output: 9
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-1000 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
include("simplepolysos.jl")
include("sosdemo1.jl")


text:
using FoodCachingPlotting
using Test

@testset "FoodCachingPlotting.jl" begin
    # Write your tests here.
end


text:
using Wavelets
using Base.Test

print("util ...\n")

@test dyadicdetailn(0) == 1
@test dyadicdetailn(1) == 2
@test dyadicdetailindex(0,1) == 2

j = 5
@test dyadicdetailn(j)+1 == dyadicdetailrange(j)[1]
@test dyadicdetailindex(j,3) == dyadicdetailrange(j)[3]
J = 7
@test ndyadicscales(2^J) == J
@test dyadicdetailn(maxdyadiclevel(2^J)) == 2^(J-1)
@test ndyadicscales(rand(2^J)) == J
@test maxdyadiclevel(rand(2^J)) == J-1

@test tl2dyadiclevel(512,1) == 8
@test dyadiclevel2tl(512,9) == 0
@test tl2dyadiclevel(rand(2^9),1) == 8
@test dyadiclevel2tl(rand(2^9),9) == 0

n = 64
L = 2
@test maxtransformlevels(n) == ndyadicscales(n)
@test detailindex(n, L, 3) == dyadicdetailindex(tl2dyadiclevel(n,L), 3)
@test detailrange(n, L) == dyadicdetailrange(tl2dyadiclevel(n,L))
@test detailn(n, L) == dyadicdetailn(tl2dyadiclevel(n,L))

# UTILITY FUNCTIONS

@test mirror([1]) == [1]
@test mirror([2,3]) == [2,-3]
@test mirror([2,3,4]) == [2,-3,4]
@test mirror([4.9,5,6,7]) == [4.9,-5,6,-7]

@test upsample([1]) == [1,0]
@test upsample([1],1) == [0,1]
@test upsample([1,2]) == [1,0,2,0]
@test upsample([1,2],1) == [0,1,0,2]
@test downsample([1,2]) == [1]
@test downsample([1,2],1) == [2]
@test downsample([1,2,3,4]) == [1,3]
@test downsample([1,2,3,4],1) == [2,4]

@test iscube(rand(4)) == true
@test iscube(rand(4,4,4,4)) == true
@test iscube(rand(4,5)) == false
@test isdyadic(rand(4)) == true
@test isdyadic(rand(32,32)) == true
@test isdyadic(rand(6)) == false

n = 64
@test wcount(randn(n)) == n
@test wcount(rand(n),0) == n
@test wcount(rand(n),1.01) == 0
@test wcount(rand(n).-2.5,1.01) == n
@test wcount([0,2,6,7,8],6.6) == 2
@test wcount([10,-11,6,7,8.0,-5,-8,0],3,level=2) == 3
@test wcount(randn(n,n)) == n*n
@test wcount(rand(n,n),1.01) == 0
@test wcount([-1 2;3 4],2.5) == 2

n = 32
x = randn(n)
y = randn(n)
for sh in (-43, -5, 0, 1, 28, 41)
    @test circshift(x,sh) == Util.circshift!(copy(x),sh)
    @test circshift(x,-sh) == Util.circshift!(copy(x),-sh)
    @test x == Util.circshift!(Util.circshift!(copy(x),sh),-sh)
    @test circshift(x,sh) == Util.circshift!(similar(x),x,sh)
    @test circshift(x,-sh) == Util.circshift!(similar(x),x,-sh)
end

a0 = [1,2]
exp = [1,2]
n = 2
@test Util.split!(copy(a0)) == exp
@test Util.split!(copy(a0), n, similar(a0)) == exp
@test Util.split!(similar(a0), copy(a0), n) == exp
@test Util.merge!(copy(a0)) == exp
@test Util.merge!(copy(a0), n, similar(a0)) == exp
@test Util.merge!(similar(a0), copy(a0), n) == exp

a0 = [-1,2,3.3,4]
exp = [-1,3.3,2,4]
n = 4
@test Util.split!(copy(a0)) == exp
@test Util.split!(copy(a0), n, similar(a0)) == exp
@test Util.split!(similar(a0), copy(a0), n) == exp
a0 = [-1,3.3,2,4]
exp = [-1,2,3.3,4]
@test Util.merge!(copy(a0)) == exp
@test Util.merge!(copy(a0), n, similar(a0)) == exp
@test Util.merge!(similar(a0), copy(a0), n) == exp

a0 = randn(128)
@test a0 == Util.merge!(Util.split!(a0))

# with strides
a0 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
ia = 3
inca = 2
@test Util.split!(similar(a0), copy(a0), ia, inca, 2)[1:2] == [3,5]
@test Util.split!(similar(a0), copy(a0), ia, inca, 4)[1:4] == [3,7,5,9]
@test Util.split!(similar(a0), copy(a0), 1, 3, 4)[1:4] == [1,7,4,10]
@test Util.split!(similar(a0), copy(a0), 1, 1, 8)[1:8] ==  Util.split!(similar(a0), copy(a0), 8)[1:8]

@test Util.merge!(similar(a0), ia, inca, copy(a0), 2)[ia:inca:ia+(2-1)*inca] == [1,2]
@test Util.merge!(similar(a0), ia, inca, copy(a0), 4)[ia:inca:ia+(4-1)*inca] == [1,3,2,4]
ia = 1
inca = 3
@test Util.merge!(similar(a0), ia, inca, copy(a0), 4)[ia:inca:ia+(4-1)*inca] == [1,3,2,4]
@test Util.merge!(similar(a0), 1, 1, copy(a0), 8)[1:8] ==  Util.merge!(similar(a0), copy(a0), 8)[1:8]

n = 128
x = randn(n)
for L = 0:ndyadicscales(n)
    for st in (:full, :dwt)
        @test isvalidtree(x, maketree(n, L, st))
    end
end
tree = maketree(n, 4)
tree[5] = false
@test !isvalidtree(x, tree)
tree = maketree(n, 4)
tree[18] = true
@test isvalidtree(x, tree)
tree[7] = false
@test !isvalidtree(x, tree)
tree[7] = true
tree[9] = false
@test !isvalidtree(x, tree)

EE = Exception
@test_throws EE maketree(n, 4, :foo)

makewavelet(wavelet(WT.db2, WT.Filter))
for tf in ("Blocks", "HeaviSine")
    testfunction(8, tf)
end
EE = Exception
@test_throws EE testfunction(8, "sksksks")


text:
mutable struct ImageFile
    f
    caption
end
ImageFile(f) = ImageFile(f, "")

"""

Take an image file and encode it

The \\ keeps linebreaks when run through Base.Markdown.parse

"""
gif_to_data_tpl = """

{{{:caption}}}

"""

function gif_to_data(imgfile::AbstractString, caption="")
    data = base64encode(read(imgfile, String))
    Mustache.render(gif_to_data_tpl, data=data, caption=md(caption))
end

function gif_to_data(img, caption="")
    imgfile = tempname() * ".gif"
    io = open(imgfile, "w")
    show(io, "image/png", img)
    close(io)
    gif_to_data(imgfile, caption)
end


text:
# This code presents an example to illustrate how SmolyakApprox can be used

using SmolyakApprox

function test_smolyak_approx()

  d  = 5  # Set the number of dimensions
  mu = 3  # Set the level of approximation

  g2, m2 =  smolyak_grid(clenshaw_curtis_equidistant,d,mu)
  grid, multi_ind = smolyak_grid(chebyshev_gauss_lobatto,d,mu)  # Construct the Smolyak grid and the multi index

  # An arbitrary test function

  function test(grid)

    y_value = (grid[:,1].+1).^0.1.*exp.(grid[:,2]).*log.(grid[:,3].+2).^0.2.*(grid[:,4].+2).^0.8.*(grid[:,5].+7).^0.1

    return y_value

  end

  y = test(grid)  # Evaluate the test function on the Smolyak grid
  y_pl = test(g2)

  point = [0.75, 0.45, 0.82, -0.15, -0.95]  # Choose a point to evaluate the approximated function

  # One way of computing the weights and evaluating the approximated function

  weights = smolyak_weights(y,grid,multi_ind)        # Compute the Smolyak weights
  y_hat = smolyak_evaluate(weights,point,multi_ind)  # Evaluate the approximated function

  #= A second way of computing the weights and evaluating the approximated function that
  computes the interpolation matrix just once. =#

  interp_mat = smolyak_inverse_interpolation_matrix(grid,multi_ind)  # Compute the interpolation matrix
  w = smolyak_weights(y,interp_mat)             # Compute the Smolyak weights
  y_hatt = smolyak_evaluate(w,point,multi_ind)  # Evaluate the approximated function

  # Piecewise linear

  w_pl = smolyak_pl_weights(y_pl,g2,m2)
  y_pl_hat = smolyak_pl_evaluate(w_pl,point,g2,m2)

  # Evaluate the exact function at point

  y_actual = test(point')

  # Now consider the ansiotropic case

  mu = [3, 2, 2, 2, 3]
  grid, multi_ind = smolyak_grid(chebyshev_gauss_lobatto,d,mu)  # Construct the Smolyak grid and the multi index
  y = test(grid)
  weights = smolyak_weights(y,grid,multi_ind)                   # Compute the Smolyak weights
  y_hat_ansio = smolyak_evaluate(weights,point,multi_ind)       # Evaluate the approximated function
  weights_th = smolyak_weights_threaded(y,grid,multi_ind)                   # Compute the Smolyak weights

  g3, m3 =  smolyak_grid(clenshaw_curtis_equidistant,d,mu)
  y_pl = test(g3)
  w_pl_ansio = smolyak_pl_weights(y_pl,g3,m3)
  y_pl_hat_ansio = smolyak_pl_evaluate(w_pl,point,g3,m3)
  w_pl_ansio_th = smolyak_pl_weights_threaded(y_pl,g3,m3)

  # Now test the full grid results

  mu = 3
  grid_full, multi_ind_full = smolyak_grid_full(chebyshev_gauss_lobatto,d,mu) # Construct the Smolyak grid and the multi index
  y_full = test(grid_full)
  weights_full = smolyak_weights_full(y_full,grid_full,multi_ind_full)        # Compute the Smolyak weights
  y_hat_full = smolyak_evaluate_full(weights_full,point,multi_ind_full)       # Evaluate the approximated function

  return y_hat, y_hatt, y_pl_hat, y_actual, y_hat_ansio, y_pl_hat_ansio, y_hat_full

end

test_smolyak_approx()


text:
function find_direction(A, ∇f, scheme::QuasiNewton{<:DirectApprox})
   return -(A\∇f)
end
function find_direction(A, ∇f, scheme::QuasiNewton{<:InverseApprox})
   return -A*∇f
end

function find_direction!(d, B, ∇f, scheme::QuasiNewton{<:DirectApprox})
   d .= -(B\∇f)
   d
end
function find_direction!(d, A, ∇f, scheme::QuasiNewton{<:InverseApprox})
   rmul!(mul!(d, A, ∇f), -1)
   d
end


text:
function directorySearch(directory, searchString)
    outfile  = open("../data/fileList.txt","w")
    fileList = filter(x->occursin(searchString, x), readdir(directory))

    for file in fileList
        println(outfile, file)
    end
    close(outfile)
end

directorySearch(pwd(),".jl")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module ClimaCorePlots

import RecipesBase
import TriplotBase

import ClimaCore:
    ClimaCore,
    DataLayouts,
    Fields,
    Geometry,
    Meshes,
    Operators,
    Topologies,
    Spaces

RecipesBase.@recipe function f(space::Spaces.SpectralElementSpace2D;)
    quad = Spaces.quadrature_style(space)
    quad_name = Base.typename(typeof(quad)).name
    dof = Spaces.Quadratures.degrees_of_freedom(quad)

    topology = Spaces.topology(space)
    mesh = topology.mesh
    @assert mesh isa Meshes.RectilinearMesh "plotting only defined for RectilinearMesh"

    n1 = Meshes.nelements(mesh.intervalmesh1)
    n2 = Meshes.nelements(mesh.intervalmesh2)

    coord_field = Fields.coordinate_field(space)
    x1coord = vec(parent(coord_field)[:, :, 1, :])
    x2coord = vec(parent(coord_field)[:, :, 2, :])

    coord_symbols = propertynames(coord_field)

    seriestype := :scatter
    title --> "$n1 × $n2 $quad_name{$dof} element space"
    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    marker --> :cross
    markersize --> 1
    seriescolor --> :blue
    legend --> false

    (x1coord, x2coord)
end

RecipesBase.@recipe function f(space::Spaces.ExtrudedFiniteDifferenceSpace)
    coord_field = Fields.coordinate_field(space)
    data = Fields.field_values(coord_field)
    Ni, Nj, _, Nv, Nh = size(data)

    #TODO: assumes VIFH layout
    @assert Nj == 1 "plotting only defined for 1D extruded fields"

    hspace = space.horizontal_space

    quad = Spaces.quadrature_style(hspace)
    quad_name = Base.typename(typeof(quad)).name
    dof = Spaces.Quadratures.degrees_of_freedom(quad)

    coord_symbols = propertynames(coord_field)
    hcoord = vec(parent(coord_field)[:, :, 1, :])
    vcoord = vec(parent(coord_field)[:, :, 2, :])

    stagger = space.staggering isa Spaces.CellCenter ? :center : :face

    seriestype := :scatter
    title --> "$Nh $quad_name{$dof} element × $Nv level $stagger space"
    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    marker --> :cross
    markersize --> 1
    seriescolor --> :blue
    legend --> false

    (hcoord, vcoord)
end

RecipesBase.@recipe function f(field::Fields.FiniteDifferenceField)
    # unwrap the data to plot
    space = axes(field)
    coord_field = Fields.coordinate_field(space)

    xdata = parent(field)[:, 1]
    ydata = parent(Spaces.coordinates_data(space))[:, 1]

    coord_symbols = propertynames(coord_field)

    # set the plot attributes
    xguide --> "value"
    yguide --> (
        if field isa Spaces.FaceFiniteDifferenceSpace
            "$(coord_symbols[1]) faces"
        else
            "$(coord_symbols[1]) centers"
        end
    )

    # fix the ylim to the column space (domain)
    ylims --> extrema(ydata)
    linewidth --> 2
    legend --> false

    (xdata, ydata)
end

RecipesBase.@recipe function f(
    field::Fields.SpectralElementField2D;
    interpolate = 10,
)
    @assert interpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    # compute the interpolated data to plot
    space = axes(field)
    topology = Spaces.topology(space)
    mesh = topology.mesh

    Nu = interpolate
    coord_field = Fields.coordinate_field(space)

    M_coords = Operators.matrix_interpolate(coord_field, Nu)
    M = Operators.matrix_interpolate(field, Nu)

    domain = Meshes.domain(mesh)
    x1min = Geometry.component(domain.interval1.coord_min, 1)
    x2min = Geometry.component(domain.interval2.coord_min, 1)
    x1max = Geometry.component(domain.interval1.coord_max, 1)
    x2max = Geometry.component(domain.interval2.coord_max, 1)

    # our interpolated field is transposed
    x1coord = [Geometry.component(pt, 1) for pt in M_coords[:, 1]]
    x2coord = [Geometry.component(pt, 2) for pt in M_coords[1, :]]

    coord_symbols = propertynames(coord_field)

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (x1coord, x2coord, M')
end

function _slice_triplot(field, hinterpolate, ncolors)
    data = Fields.field_values(field)
    Ni, Nj, _, Nv, Nh = size(data)

    space = axes(field)
    @assert Nj == 1

    hcoord_field = getproperty(Fields.coordinate_field(space), 1)
    vcoord_field = getproperty(Fields.coordinate_field(space), 2)

    if hinterpolate ≥ 1
        Nu = hinterpolate
        uquad = Spaces.Quadratures.ClosedUniform{Nu}()
        M_hcoord = Operators.matrix_interpolate(hcoord_field, uquad)
        M_vcoord = Operators.matrix_interpolate(vcoord_field, uquad)
        M_data = Operators.matrix_interpolate(field, uquad)

        hcoord_data = vec(M_hcoord)
        vcoord_data = vec(M_vcoord)
        data = vec(M_data)
        triangles = hcat(Spaces.triangles(Nv, Nu, Nh)...)'
    else
        hcoord_data = vec(parent(hcoord_field))
        vcoord_data = vec(parent(vcoord_field))
        data = vec(parent(data))
        triangles = hcat(Spaces.triangles(Nv, Ni, Nh)...)'
    end
    cmap = range(extrema(data)..., length = ncolors)

    z = TriplotBase.tripcolor(
        hcoord_data,
        vcoord_data,
        data,
        triangles,
        cmap;
        bg = NaN,
        # unique number of nodal element coords
        # (number of nodal values minus number of shared faces)
        px = (Ni * Nh) - (Nh - 1),
        py = Nv,
    )
    # some plots backends need coordinates in sorted order
    return (sort(unique(hcoord_data)), sort(unique(vcoord_data)), z')
end

# 2D hybrid plot
RecipesBase.@recipe function f(
    field::Fields.Field{<:Any, S};
    hinterpolate = 0,
    ncolors = 256,
) where {
    S <: Spaces.ExtrudedFiniteDifferenceSpace{
        <:Any,
        <:Spaces.IntervalSpectralElementSpace1D,
    },
}
    hcoord, vcoord, data = _slice_triplot(field, hinterpolate, ncolors)
    coord_symbols = propertynames(Fields.coordinate_field(axes(field)))

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (hcoord, vcoord, data)
end

function _slice_along(field, coord)
    # which axis are we slicing?
    axis = if coord isa ClimaCore.Geometry.XPoint
        1
    elseif coord isa ClimaCore.Geometry.YPoint
        2
    else
        error(
            "coord must be a ClimaCore.Geometry.XPoint or ClimaCore.Geometry.YPoint",
        )
    end
    space = axes(field)
    hspace = space.horizontal_space
    htopo = ClimaCore.Spaces.topology(hspace)
    hmesh = htopo.mesh
    linear_idx = LinearIndices(ClimaCore.Meshes.elements(hmesh))

    # get the axis interval mesh we are slicing along
    hmesh_slice, hmesh_ortho = if axis == 1
        (htopo.mesh.intervalmesh1, htopo.mesh.intervalmesh2)
    else
        (htopo.mesh.intervalmesh2, htopo.mesh.intervalmesh1)
    end

    # check that the coordinate is in the domain
    domain = hmesh_slice.domain
    if coord < domain.coord_min || coord > domain.coord_max
        error("$coord is outside of the field axis domain: $domain")
    end

    # get the elment offset given the axis slice
    slice_h, _ = ClimaCore.Meshes.containing_element(hmesh_slice, coord)
    hidx = axis == 1 ? linear_idx[slice_h, 1] : linear_idx[1, slice_h]

    # find the node idx we want to slice along the given axis element
    hcoord_data = hspace.local_geometry.coordinates
    hdata = ClimaCore.slab(hcoord_data, hidx)
    hnode_idx = 1
    for i in axes(hdata)[axis]
        pt = axis == 1 ? hdata[i, 1] : hdata[1, i]
        axis_value = Geometry.component(pt, axis)
        coord_value = Geometry.component(coord, 1)
        if axis_value > coord_value
            break
        end
        hnode_idx = i
    end

    # construct the slice field space
    htopo_ortho = ClimaCore.Topologies.IntervalTopology(hmesh_ortho)
    hspace_ortho = ClimaCore.Spaces.SpectralElementSpace1D(
        htopo_ortho,
        ClimaCore.Spaces.quadrature_style(hspace),
    )
    vspace_ortho =
        ClimaCore.Spaces.CenterFiniteDifferenceSpace(space.vertical_topology)

    if space.staggering === ClimaCore.Spaces.CellFace
        vspace_ortho = ClimaCore.Spaces.FaceFiniteDifferenceSpace(slice_vspace)
    end
    space_ortho = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
        hspace_ortho,
        vspace_ortho,
    )

    # construct the slice field ortho to the slice coordinate
    ortho_field = Fields.Field(eltype(field), space_ortho)

    # copy the data into this field for the given node idx for all elements intersecting slice coord
    field_data = ClimaCore.Fields.todata(field)
    ortho_data = ClimaCore.Fields.todata(ortho_field)

    for i in 1:size(linear_idx, axis == 1 ? 2 : 1)
        for v in 1:ClimaCore.Spaces.nlevels(space_ortho)
            hidx = axis == 1 ? linear_idx[slice_h, i] : linear_idx[i, slice_h]
            ijslab = ClimaCore.slab(field_data, v, hidx)
            islab = ClimaCore.slab(ortho_data, v, i)
            # copy the nodal data
            for ni in 1:size(islab)[1]
                islab[ni] =
                    axis == 1 ? ijslab[hnode_idx, ni] : ijslab[ni, hnode_idx]
            end
        end
    end
    return ortho_field
end

# 3D hybrid plot
RecipesBase.@recipe function f(
    field::Fields.Field{<:Any, S};
    slice = nothing,
    hinterpolate = 0,
    ncolors = 256,
) where {
    S <: Spaces.ExtrudedFiniteDifferenceSpace{
        <:Any,
        <:Spaces.RectilinearSpectralElementSpace2D,
    },
}
    if slice === nothing
        error("must specify coordinate axis slice for 3D hybrid plots")
    end
    if length(slice) != 3
        error("must specify a length 3 slice for 3D hybrid plots")
    end
    if count(c -> c isa Colon, slice) != 2
        error("can only specify one axis to slice")
    end
    if !(slice[3] isa Colon)
        error("slicing only along the X and Y axis is supported")
    end
    slice_coord = if slice[1] isa Colon
        ClimaCore.Geometry.YPoint{eltype(slice[2])}(slice[2])
    else
        ClimaCore.Geometry.XPoint{eltype(slice[1])}(slice[1])
    end
    slice_field = _slice_along(field, slice_coord)
    hcoord, vcoord, data = _slice_triplot(slice_field, hinterpolate, ncolors)

    coord_symbols = propertynames(Fields.coordinate_field(axes(slice_field)))

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (hcoord, vcoord, data)
end

function _unfolded_pannel_matrix(field, interpolate)
    space = axes(field)
    topology = Spaces.topology(space)
    mesh = topology.mesh
    nelem = Meshes.nelements(mesh)
    panel_size = mesh.ne

    quad_from = Spaces.quadrature_style(space)
    quad_to = Spaces.Quadratures.Uniform{interpolate}()
    Imat = Spaces.Quadratures.interpolation_matrix(Float64, quad_to, quad_from)

    dof = interpolate

    pannel_range(i) =
        ((panel_size * dof) * (i - 1) + 1):((panel_size * dof) * i)

    # construct a matrix to fill in the rotated / flipped pannel data
    unfolded_panels =
        fill(NaN, ((panel_size * dof) * 3, (panel_size * dof) * 4))

    # temporary pannels as we have to rotate / flip some and not all operators are in place
    # TODO: inefficient memory wise, but good enough for now
    panels = [fill(NaN, (panel_size * dof, panel_size * dof)) for _ in 1:6]

    interpolated_data =
        DataLayouts.IJFH{Float64, interpolate}(Array{Float64}, nelem)
    field_data = Fields.field_values(field)

    Operators.tensor_product!(interpolated_data, field_data, Imat)

    # element index ordering defined by a specific layout
    eidx = 1
    for panel_idx in 1:6
        panel_data = panels[panel_idx]
        # elements are ordered along fastest axis defined in sphere box mesh
        for ex2 in 1:panel_size, ex1 in 1:panel_size
            # compute the nodal extent index range for this element
            x1_nodal_range = (dof * (ex1 - 1) + 1):(dof * ex1)
            x2_nodal_range = (dof * (ex2 - 1) + 1):(dof * ex2)
            # transpose the data as our plotting axis order is
            # reverse nodal element order (x1 axis varies fastest)
            data_element = permutedims(parent(interpolated_data)[:, :, 1, eidx])
            panel_data[x2_nodal_range, x1_nodal_range] = data_element
            eidx += 1
        end
    end

    # (px, py, rot) for each panel
    # px, py are the locations of the panel
    # rot is number of clockwise rotations (0:3)
    # https://extranet.gfdl.noaa.gov/~atw/ferret/cubed_sphere/
    # Equatorial strip, poles connected through the Americas
    panel_locations =
        [(4, 2, 0), (1, 2, 0), (3, 3, 2), (2, 2, 1), (3, 2, 1), (3, 1, 1)]
    for (i, (px, py, rot)) in enumerate(panel_locations)
        unfolded_panels[pannel_range(py), pannel_range(px)] .= if rot == 0
            panels[i]
        elseif rot == 1
            reverse(transpose(panels[i]), dims = 1)
        elseif rot == 2
            reverse(panels[i], dims = (1, 2))
        else
            reverse(transpose(panels[i]), dims = 2)
        end
    end
    #=
    unfolded_panels[pannel_range(1), pannel_range(2)] =
        reverse(panels[5], dims = 1)
    unfolded_panels[pannel_range(2), pannel_range(1)] =
        reverse(panels[4], dims = 2)
    unfolded_panels[pannel_range(2), pannel_range(2)] = transpose(panels[1])
    unfolded_panels[pannel_range(2), pannel_range(3)] = transpose(panels[2])
    unfolded_panels[pannel_range(2), pannel_range(4)] =
        reverse(panels[6], dims = 2)
    unfolded_panels[pannel_range(3), pannel_range(2)] =
        reverse(panels[3], dims = 2)
    =#
    return unfolded_panels
end

RecipesBase.@recipe function f(
    field::Fields.CubedSphereSpectralElementField2D;
    interpolate = 10,
)
    @assert interpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    unfolded_panels = _unfolded_pannel_matrix(field, interpolate)

    # construct the title for info about the field space
    space = axes(field)
    nelem = Topologies.nlocalelems(Spaces.topology(space))
    quad_from = Spaces.quadrature_style(space)
    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
    quad_from_name = Base.typename(typeof(quad_from)).name

    # set the plot attributes
    seriestype := :heatmap
    title --> "$nelem $quad_from_name{$dof_in} element space"
    xguide --> "panel x1"
    yguide --> "panel x2"
    seriescolor --> :balance

    (unfolded_panels)
end

RecipesBase.@recipe function f(
    field::Fields.Field{<:Any, S};
    level = nothing,
    hinterpolate = 10,
) where {
    S <: Spaces.ExtrudedFiniteDifferenceSpace{
        <:Any,
        <:Spaces.CubedSphereSpectralElementSpace2D,
    },
}
    @assert hinterpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    space = axes(field)
    if level === nothing
        if space isa Spaces.FaceExtrudedFiniteDifferenceSpace
            level = Utilities.PlusHalf(0)
        else
            level = 1
        end
    end
    level_field = Fields.level(field, level)
    unfolded_panels = _unfolded_pannel_matrix(level_field, hinterpolate)

    # construct the title for info about the field space
    nlevel = Spaces.nlevels(space)
    nelem = Topologies.nlocalelems(Spaces.topology(space))
    quad_from = Spaces.quadrature_style(space)
    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
    quad_from_name = Base.typename(typeof(quad_from)).name

    # set the plot attributes
    seriestype := :heatmap
    title -->
    "level $level of $nlevel × $nelem $quad_from_name{$dof_in} element space"
    xguide --> "panel x1"
    yguide --> "panel x2"
    seriescolor --> :balance

    (unfolded_panels)
end

end # module


text:
import Base.randperm

"The path graph on n vertices"
function pathGraph(n::Int64)
  x = append!(collect(1:(n-1)), collect(2:n))
  y = append!(collect(2:n), collect(1:(n-1)))
  w = ones(2 * (n - 1))
  return sparse(x, y, w)
end # pathGraph

"The complete graph"
function completeGraph(n::Int64)
  return sparse(ones(n,n) - eye(n))
end # completeGraph

"""The simple ring on n vertices"""
function ringGraph(n::Int64)
    a = spdiagm(ones(n-1),1,n,n)
    a[1,n] = 1
    a = a + a'
end

"""A generalization of a ring graph.
The vertices are integers modulo n.
Two are connected if their difference is in gens.
For example,

```
generalizedRing(17, [1 5])
```
"""
function generalizedRing(n::Int64, gens)
    k = length(gens)
    m = 2*n*k
    ai = zeros(Int64,m)
    aj = zeros(Int64,m)
    ind = 1
    for i in 0:(n-1)
        for j in 1:k
            ai[ind] = i
            aj[ind] = mod(i+gens[j],n)
            ind = ind + 1
            ai[ind] = i
            aj[ind] = mod(i-gens[j],n)
            ind = ind + 1
        end
    end
    return sparse(1+ai,1+aj,ones(m),n,n)
    #return ai, aj
end

"""A random generalized ring graph of degree k.
Gens always contains 1, and the other k-1 edge types
are chosen from an exponential distribution"""
function randGenRing(n::Int64, k::Integer)
    # if any of n, 2n, 3n etc. is in gens we will have self loops
    gens = [0]
    while 0 in (gens % n)
        gens = [1; 1 + ceil(Integer,exp(rand(k-1)*log(n-1)))]
    end

    return generalizedRing(n, gens)
end

"""
~~~
generalizedNecklace{Tv, Ti}(A::SparseMatrixCSC{Tv, Ti}, H::SparseMatrixCSC, k::Int64)
~~~

 Constructs a generalized necklace graph starting with two graphs A and H. The
resulting new graph will be constructed by expanding each vertex in H to an
instance of A. k random edges will be generated between components. Thus, the
resulting graph may have weighted edges.
"""

function generalizedNecklace{Tv, Ti}(A::SparseMatrixCSC{Tv, Ti}, H::SparseMatrixCSC, k::Int64)
  a = findnz(A)
  h = findnz(H)

  # these are square matrices
  n = A.n
  m = H.n

  newI = Ti[]
  newJ = Ti[]
  newW = Tv[]

  # duplicate the vertices in A so that each vertex in H corresponds to a copy of A
  for i in 1:m
    newI = append!(newI, a[1] + n * (i - 1))
    newJ = append!(newJ, a[2] + n * (i - 1))
    newW = append!(newW, a[3])
  end

  # for each edge in H, add k random edges between two corresponding components
  # multiedges will be concatenated to a single edge with higher cost
  for i in 1:length(h[1])
    u = h[1][i]
    v = h[2][i]

    if (u < v)
      #component x is from 1 + (x - 1) * n to n + (x - 1) * n
      for edgeToAdd in 1:k
        newU = rand(1:n) + n * (u - 1)
        newV = rand(1:n) + n * (v - 1)
        append!(newI, [newU, newV])
        append!(newJ, [newV, newU])
        append!(newW, [1, 1])
      end
    end
  end

  return sparse(newI, newJ, newW)
end # generalizedNecklace

"""The d dimensional hypercube.  Has 2^d vertices"""
function hyperCube(d::Int64)
  a = sparse([0 1; 1 0])

  for i = 1:(d-1)
    k = 2^i
    D = speye(k)
    a = [a D; D a]
  end

  return a
end # hyperCube

"""The complete binary tree on n vertices"""
function completeBinaryTree(n::Int64)

  k = div(n-1,2)
  a = sparse(collect(1:k),2*collect(1:k),1,n,n) + sparse(collect(1:k),2*collect(1:k)+1,1,n,n)

  if 2*k+1 < n
    a[n-1,n] = 1
  end

  a = a + a'

  return a
end # completeBinaryTree

""" An n by n grid with random weights. User can specify the weighting scheme. """
function wGrid2(n::Int64; weightGen::Function=rand)
    gr2 = sparse(grid2(n));

    gr2.nzval = Float64[weightGen() for i in 1:nnz(gr2)]

    # symmetrize
    gr2 = tril(gr2) + tril(gr2)'

    return gr2
end

""" An n^3 grid with random weights. User can specify the weighting scheme. """
function wGrid3(n::Int64; weightGen::Function=rand)
    gr2 = grid2(n);

    a = kron(speye(n), gr2);
    b = kron(gr2, speye(n));

    gr3 = sparse(a + b);
    gr3.nzval = Float64[weightGen() for i in 1:nnz(gr3)]

    # symmetrize
    gr3 = tril(gr3) + tril(gr3)'

    return gr3
end

"""An n-by-m grid graph.  iostropy is the weighting on edges in one direction."""
function grid2(n::Int64, m::Int64; isotropy=1)
  a = kron(speye(n),spdiagm(ones(m-1),1,m,m))
  a = a + isotropy*kron(spdiagm(ones(n-1),1,n,n), speye(m))
  a = a + a'
  return a
end # grid2

grid2(n::Int64) = grid2(n,n)

"""Coordinates for plotting the vertices of the n-by-m grid graph"""
function grid2coords(n::Int64, m::Int64)
  x = kron(collect(1:n),ones(m))
  y = kron(ones(n),collect(1:m))
  return x, y
end # grid2coords

grid2coords(n) = grid2coords(n, n)

"""A random matching on n vertices"""
function randMatching(n::Int64)

  p = randperm(n)
  n1 = convert(Int64,floor(n/2))
  n2 = 2*n1
  a = sparse(p[1:n1],p[(n1+1):n2],ones(n1),n,n)

  a = a + a'

  return a

end # randMatching

"""A sum of k random matchings on n vertices"""
function randRegular(n::Int64, k::Int64)
  a = randMatching(n)
  for i in 2:k
    a = a + randMatching(n)
  end

  return a
end # randRegular

"""Create a graph on n vertices.
For each vertex, give it k edges to randomly chosen prior
vertices.
This is a variety of a preferential attachment graph.
"""
function grownGraph(n::Int64, k::Int64)
  a = spzeros(n,n)

  for i = 1:k
    a = a + sparse(2:n,ceil(Integer,collect(1:n-1).*rand(n-1)),1,n,n)
  end

  a = a + a'
end # grownGraph

# used in grownGraphD
function randSet(n::Integer,k::Integer)
    if n == k
        return collect(1:n)
    elseif n < k
        error("n must be at least k")
    else

        s = sort(ceil(Integer,n*rand(k)))
        good = (minimum(s[2:end]-s[1:(end-1)]) > 0)
        while good == false
            s = sort(ceil(Integer,n*rand(k)))
            good = (minimum(s[2:end]-s[1:(end-1)]) > 0)
        end

        return s

    end
end

"""Like a grownGraph, but it forces the edges to all be distinct.
It starts out with a k+1 clique on the first k vertices"""
function grownGraphD(n::Int64, k::Int64)
    a = spzeros(n,n)

    u = zeros(Int64, k*(n-k-1))
    v = zeros(Int64, k*(n-k-1))

    for i in (k+2):n
        nb = randSet(i-1,k)
        u[(i-k-2)*k + collect(1:k)] = i
        v[(i-k-2)*k + collect(1:k)] = nb
    end

    a = sparse(u,v,1,n,n)

    (ai,aj) = findnz(triu(ones(k+1,k+1),1))
    a = a + sparse(ai,aj,1,n,n)
    a = a + a'

end # grownGraphD

"""A preferential attachment graph in which each vertex has k edges to those
that come before.  These are chosen with probability p to be from a random vertex,
and with probability 1-p to come from the endpoint of a random edge.
It begins with a k-clique on the first k+1 vertices."""
function prefAttach(n::Int64, k::Int64, p::Float64)
    if n == (k+1)
        return sparse(ones(Float64,n,n) - eye(Float64,n))
    elseif n <= k
        error("n must be more than k")
    else

        u = zeros(Int64,n*k)
        v = zeros(Int64,n*k)

        # fill in the initial clique
        # this will accidentally double every edge in the clique
        # we clean it up at the end
        ind = 1
        for i in 1:(k+1)
            for j in 1:(k+1)
                if i != j
                    u[ind] = i
                    v[ind] = j
                    ind += 1
                end
            end
        end

        s = zeros(Int64,k)
        for i in (k+2):n
            distinct = false
            while distinct == false
                for j in 1:k
                    if rand(Float64) < p
                        s[j] = rand(1:(i-1))
                    else
                        s[j] = v[rand(1:(k*(i-1)))]
                    end
                end
                s = sort(s)
                distinct = true
                for ii in 1:(k-1)
                    if s[ii] == s[ii+1]
                        distinct = false
                    end
                end
                # distinct = (minimum(s[2:end]-s[1:(end-1)]) > 0)

            end

            for j in 1:k
                u[ind] = i
                v[ind] = s[j]
                ind += 1
            end

        end # for i

        w = ones(Float64,n*k)

        w[1:(k*(k+1))] = 1/2

        a = sparse(u,v,w,n,n)
        a = a + a'
        return a
    end
end

"""Randomly permutes the vertex indices"""
function randperm(mat::AbstractMatrix)
    perm = randperm(mat.n)
    return mat[perm,perm]
end

randperm(f::Expr) = randperm(eval(f))

"""Generate a random graph on n vertices with m edges.
The actual number of edges will probably be smaller, as we sample
with replacement"""
function ErdosRenyi(n::Integer, m::Integer)
    ai = rand(1:n, m)
    aj = rand(1:n, m)
    ind = (ai .!= aj)
    mat = sparse(ai[ind],aj[ind],1,n,n)
    uniformWeight!(mat)
    return mat
end

"""Generate an ER graph with average degree k,
and then return the largest component.
Will probably have fewer than n vertices.
If you want to add a tree to bring it back to n,
try ErdosRenyiClusterFix."""
function ErdosRenyiCluster(n::Integer, k::Integer)
    m = ceil(Integer,n*k/2)
    ai = rand(1:n, m)
    aj = rand(1:n, m)
    ind = (ai .!= aj)
    mat = sparse(ai[ind],aj[ind],1,n,n)
    mat = mat + mat'

    return biggestComp(mat)
end

"""Like an Erdos-Renyi cluster, but add back a tree so
it has n vertices"""
function ErdosRenyiClusterFix(n::Integer, k::Integer)
    m1 = ErdosRenyiCluster(n, k)
    n2 = n - size(m1)[1]
    if (n2 > 0)
        m2 = completeBinaryTree(n2)
        return joinGraphs(m1,m2,1)
    else
        return m1
    end
end

"""Generate a random graph with n vertices from one of our natural distributions"""
function pureRandomGraph(n::Integer)

    gr = []
    wt = []

    push!(gr,:(pathGraph($n)))
    push!(wt,1)

    push!(gr,:(ringGraph($n)))
    push!(wt,3)

    push!(gr,:(completeBinaryTree($n)))
    push!(wt,3)

    push!(gr,:(grownGraph($n,2)))
    push!(wt,6)

    push!(gr,:(grid2(ceil(Integer,sqrt($n)))[1:$n,1:$n]))
    push!(wt,6)

    push!(gr,:(randRegular($n,3)))
    push!(wt,6)

    push!(gr,:(ErdosRenyiClusterFix($n,2)))
    push!(wt,6)

    if n >= 4
        push!(gr,:(randGenRing($n,4)))
        push!(wt,6)
    end

    i = sampleByWeight(wt)

    # make sure get a connected graph
    its = 0
    mat = eval(gr[i])

    while (~isConnected(mat)) && (its < 100)
        i = sampleByWeight(wt)
        mat = eval(gr[i])
        its += 1
    end
    if its == 100
        error("Getting a disconnected graph from $(gr[i])")
    end

    if (sum(diag(mat)) > 0)
        error("nonzero diag from $(gr[i])")
    end

    return floatGraph(mat)

end

"""sample an index with probability proportional to its weight given here"""
function sampleByWeight(wt)
    r = rand(1)*sum(wt)
    find(cumsum(wt) .> r)[1]
end

"""A Chimera graph with some weights.  The weights just appear when graphs are combined.
For more interesting weights, use `wtedChimera`"""
function semiWtedChimera(n::Integer)

    if (n < 2)
        gr = sparse([0.0])

        return randperm(gr)
    end

    r = rand()^2

    if (n < 30) || (rand() < .2)

        gr = pureRandomGraph(n)

        return randperm(gr)
    end

    if (n < 200)
        # just join disjoint copies of graphs

        n1 = 10 + floor(Integer,(n-20)*rand())
        n2 = n - n1
        k = ceil(Integer,exp(rand()*log(min(n1,n2)/2)))

        gr = joinGraphs(r*chimera(n1),chimera(n2),k)

        return randperm(gr)
    end

    # split with probability .7

    if (rand() < .7)
        n1 = ceil(Integer,10*exp(rand()*log(n/20)))

        n2 = n - n1
        k = floor(Integer,1+exp(rand()*log(min(n1,n2)/2)))

        gr = joinGraphs(r*chimera(n1),chimera(n2),k)

        return randperm(gr)

    else
        n1 = floor(Integer,10*exp(rand()*log(n/100)))

        n2 = floor(Integer, n / n1)

        if (rand() < .5)

            gr = productGraph(r*chimera(n1),chimera(n2))

        else

            k = floor(Integer,1+exp(rand()*log(min(n1,n2)/10)))
            gr = generalizedNecklace(r*chimera(n1),chimera(n2),k)

        end

        n3 = n - size(gr)[1]
        if (n3 > 0)
            gr = joinGraphs(gr,chimera(n3),2)

        end

        return randperm(gr)

    end
end

"""Builds a chimeric graph on n vertices.
The components come from pureRandomGraph,
connected by joinGraphs, productGraph and generalizedNecklace"""
function chimera(n::Integer)

    gr = semiWtedChimera(n)
    unweight!(gr)

    return gr

end

"""Builds the kth chimeric graph on n vertices.
It does this by resetting the random number generator seed.
It should captute the state of the generator before that and then
return it, but it does not yet."""
function chimera(n::Integer, k::Integer)
    srand(100*n+k)
    g = chimera(n)
    return g
end

"""Applies one of a number of random weighting schemes to the edges of the graph"""
function randWeight(a)

    if (rand() < .2)
        return a
    end

    n = a.n
    (ai,aj) = findnz(a)
    m = length(ai)

    # potentials or edge-based

    if (rand() < .3)
        w = rand(m)

    else
        v = randn(a.n)

        # mult by matrix ?
        if (rand() < .5)

            invdeg = spdiagm(1./(a*ones(size(a)[1])))
            if (rand() < .5)
                for i in 1:10
                    v = a * (invdeg * v)
                    v = v - mean(v)
                end
            else
                for i in 1:10
                    v = v - a * (invdeg * v)
                    v = v - mean(v)
                end
            end
        end

        w = abs(v[ai]-v[aj])

    end

    # reciprocate or not?

    w[w.==0] = 1
    w[isnan(w)] = 1

    if (rand() < .5)
        w = 1./w
    end

    w = w / mean(w)

    ar = sparse(ai,aj,w,n,n)
    ar = ar + ar';
    return ar
end

"""Builds the kth wted chimeric graph on n vertices.
It does this by resetting the random number generator seed.
It should captute the state of the generator before that and then
return it, but it does not yet."""
function wtedChimera(n::Integer, k::Integer)
    srand(100*n+k)
    g = wtedChimera(n)
    return g
end

function semiWtedChimera(n::Integer, k::Integer)
    srand(100*n+k)
    g = semiWtedChimera(n)
    return g
end

"""Generate a chimera, and then apply a random weighting scheme"""
function wtedChimera(n::Integer)
    return randWeight(semiWtedChimera(n))
end


text:
using Thrift
using Test

import Thrift: process, meta

global srvr

# include the generated module, which in-turn includes our implementation code in `proto_tests_impl.jl`
# prepend with @everywhere when running a TProcessPoolServer on multiple Julia processors
@isdefined(srvcctrl) || include("gen-jl/srvcctrl/srvcctrl.jl");
@isdefined(proto_tests) || include("gen-jl/proto_tests/proto_tests.jl");

import .proto_tests: ProtoTestsProcessor, ProtoTestsClient, InvalidOperation, AllTypes, AllTypesDefault, TestEnum
import .proto_tests: test_hello, test_exception, test_oneway, ping, test_enum, test_types, test_types_default
import .srvcctrl: start_service, stop_service

transport_factory(x) = x
protocol_factory(x) = TBinaryProtocol(x)
#protocol_factory(x) = TCompactProtocol(x)

function make_server()
    # create a server instance with our choice of protocol and transport
    srvr_processor = ProtoTestsProcessor()
    srvr_transport = TServerSocket(19999)

    #srvr = TProcessPoolServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)
    #srvr = TTaskServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)
    srvr = TSimpleServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)

    @info("server transport: $(typeof(srvr_transport))")
    @info("server protocol : $(typeof(protocol_factory(srvr_transport)))")
    @info("server processor: $(typeof(srvr_processor))")
    @info("server type: $(typeof(srvr))")
    srvr
end

srvr = make_server()
# start serving client requests
@debug("server starting to serve requests...")

if haskey(ENV, "TEST_SRVR_ASYNC")
    @async try
        start_service(srvr)
    catch ex
        isa(ex, Base.IOError) || @error("Server stopped with $ex")
    finally
        @info("server stopped serving requests.")
    end
else
    start_service(srvr)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d1b3d67a-2134-11eb-37c0-31d07bf1dec8
begin
	using PlutoUI
	using Plots
	using LinearAlgebra
end

# ╔═╡ 442454c8-2135-11eb-0531-fb9566ef35a1
@bind a Slider(1:20; default=8, show_value=true)

# ╔═╡ f0f3a8be-2135-11eb-2707-7168315ff140
A = [2.1 0.1;0.3 a*0.1]

# ╔═╡ 5d6030c8-2136-11eb-0778-1702370cdde4
begin

			y = [det(A)^a for a in 1:20]
			plot(x = 1:20,y,size=(400,300))

end

# ╔═╡ Cell order:
# ╠═d1b3d67a-2134-11eb-37c0-31d07bf1dec8
# ╠═442454c8-2135-11eb-0531-fb9566ef35a1
# ╠═f0f3a8be-2135-11eb-2707-7168315ff140
# ╠═5d6030c8-2136-11eb-0778-1702370cdde4


text:
"""
# Description

Rewrite an expression written in terms of column names (represented as symbols)
into an expression written in terms of indexing into a tuple using numeric
indices.

See docs/expression_rules.md for details.

# Arguments

1. `e::Any`: An expression.
2. `index::Dict{Symbol, Int}`: A dictionary mapping symbols to indices.

# Return Values

1. `e::Any`: An expression in which column name symbols are replaced by tuple
    indexing expressions.

# Examples

```
julia> rewrite_column_names(
           :(a + b),
           :t,
           Dict((:a, false) => 1, (:b, false) => 2),
       )
:(t[1] + t[2])
```
"""
function rewrite_column_names(
    @nospecialize(e::Any),
    tuple_name::Any,
    index::Dict{ColumnName, Int},
)::Any
    if isa(e, Expr) && e.head == :(=)
        Expr(
            :(=),
            e.args[1],
            rewrite_column_names(e.args[2], tuple_name, index),
        )
    elseif isa(e, Expr) && e.head == :kw
        Expr(
            :kw,
            e.args[1],
            rewrite_column_names(e.args[2], tuple_name, index),
        )
    elseif isa(e, Expr) && e.head == :call
        rewritten_args = Any[]
        for i in 2:length(e.args)
            push!(
                rewritten_args,
                rewrite_column_names(e.args[i], tuple_name, index),
            )
        end
        Expr(
            :call,
            e.args[1],
            rewritten_args...,
        )
    elseif isa(e, Expr) && e.head == :$
        e
    elseif isa(e, Expr) && e.head == :comparison
        rewritten_args = Any[]
        for i in 1:length(e.args)
            if i % 2 == 1
                push!(
                    rewritten_args,
                    rewrite_column_names(e.args[i], tuple_name, index),
                )
            else
                push!(rewritten_args, e.args[i])
            end
        end
        Expr(e.head, rewritten_args..., )
    elseif isa(e, Expr) && e.head == :macrocall && isa(e.args[1], GlobalRef) && e.args[1].name == Symbol("@cmd")
        if startswith(e.args[3], '$')
            Expr(:ref, tuple_name, index[ColumnName(Symbol(strip(e.args[3], '$')), true)])
        else
            Expr(:ref, tuple_name, index[ColumnName(Symbol(e.args[3]), false)])
        end
    elseif isa(e, Expr)
        rewritten_args = Any[]
        for i in 1:length(e.args)
            push!(
                rewritten_args,
                rewrite_column_names(e.args[i], tuple_name, index),
            )
        end
        Expr(e.head, rewritten_args...,)
    elseif isa(e, Symbol)
        Expr(:ref, tuple_name, index[ColumnName(e, false)])
    else
        e
    end
end

"""
# Description

Rewrite an expression written in terms of column names (represented as symbols)
into an expression written in terms of safe column names (symbols from gensym
calls).

See docs/expression_rules.md for details.

# Arguments

1. `e::Any`: An expression.
2. `mapping::Dict{Symbol, Symbol}`: A dictionary mapping symbols to symbols.

# Return Values

1. `e::Any`: An expression in which column name symbols are replaced by tuple
    indexing expressions.

# Examples

```
julia> rewrite_column_names_broadcast(
    :(a + b),
    Dict((:a, false) => gensym(), (:b, false) => gensym()),
)
:(var"##256" + var"##257")
```
"""
function rewrite_column_names_broadcast(
    @nospecialize(e::Any),
    mapping::Dict{ColumnName, Symbol},
)
    if isa(e, Expr) && e.head == :(=)
        Expr(
            :(=),
            e.args[1],
            rewrite_column_names_broadcast(e.args[2], mapping),
        )
    elseif isa(e, Expr) && e.head == :kw
        Expr(
            :kw,
            e.args[1],
            rewrite_column_names_broadcast(e.args[2], mapping),
        )
    elseif isa(e, Expr) && e.head == :call
        rewritten_args = Any[]
        for i in 2:length(e.args)
            push!(
                rewritten_args,
                rewrite_column_names_broadcast(e.args[i], mapping),
            )
        end
        Expr(
            :call,
            e.args[1],
            rewritten_args...,
        )
    elseif isa(e, Expr) && e.head == :$
        e
    elseif isa(e, Expr) && e.head == :comparison
        rewritten_args = Any[]
        for i in 1:length(e.args)
            if i % 2 == 1
                push!(
                    rewritten_args,
                    rewrite_column_names_broadcast(e.args[i], mapping),
                )
            else
                push!(rewritten_args, e.args[i])
            end
        end
        Expr(e.head, rewritten_args..., )
    elseif isa(e, Expr) && e.head == :macrocall && isa(e.args[1], GlobalRef) && e.args[1].name == Symbol("@cmd")
        if startswith(e.args[3], '$')
            mapping[ColumnName(Symbol(strip(e.args[3], '$')), true)]
        else
            mapping[ColumnName(Symbol(e.args[3]), false)]
        end
    elseif isa(e, Expr)
        rewritten_args = Any[]
        for i in 1:length(e.args)
            push!(
                rewritten_args,
                rewrite_column_names_broadcast(e.args[i], mapping),
            )
        end
        Expr(e.head, rewritten_args...,)
    elseif isa(e, Symbol)
        mapping[ColumnName(e, false)]
    else
        e
    end
end


text:
#==============================================================================#
# HTML2MD.jl
#
# Julia wrapper for the NodeJS to-markdown module.
#
# Copyright OC Technology Pty Ltd 2017 - All rights reserved
#==============================================================================#

module HTML2MD

export html2md

using NodeJS
using Requests.post
using Requests.URI

const tcp_port = 34562

server_process = nothing

function start_node_server()

    global server_process

    @assert server_process == nothing

    if !isdir(joinpath(@__DIR__, "node_modules"))
        run(setenv(`$(npm_cmd()) install to-markdown`, dir=@__DIR__))
    end

    server_process = spawn(`$(nodejs_cmd()) -e """
        const http = require('http')
        const port = $tcp_port
        const h2m = require('to-markdown')

        const requestHandler = (request, response) => {
            if (request.method == 'POST') {
                var body = '';
                request.on('data', function (data) { body += data; });
            }

            request.on('end', function () {
                response.end(h2m(body, {converters: [{

                filter: function (node) {
                    return node.nodeName === 'A' && !node.getAttribute('href')
                },
                replacement: function(content) {
                    return '[' + content + '](@ref)';
                }
            },{
                filter: ['fullname'],
                replacement: function(content) {
                    return '\n\n# ' + content + '\n\n'
                }
            },{
                filter: ['important'],
                replacement: function (content) {
                  return '\n**Important**\n> ' + content + '\n\n'
                }
            },{
                filter: ['note'],
                replacement: function (content) {
                  return '\n**Note**\n> ' + content + '\n\n'
                }
            },{
                filter: ['em', 'i'],
                replacement: function (content) {
                  return '*' + content + '*'
                }
            }]}))
            });
        }

        http.createServer(requestHandler).listen(port, (err) => {
          if (err) {
            return console.log('HTML2MD Error:', err)
          }
        })
    """`)

    atexit(() -> kill(server_process))

    sleep(2)
    server_started = true
end

function html2md(html)

    global server_process

    if server_process == nothing
        start_node_server()
    end

    html = replace(html, "\\", "\\\\")
    html = replace(html, "\$", "\\\$")
    html = replace(html, "\"\"\"", "\\\"\\\"\\\"")

    md = String(post(URI("http://localhost:$tcp_port"), html).data)

    # Work around for https://github.com/domchristie/to-markdown/issues/181
    md = replace(md, r"([0-9])\\[.]", s"\1.")

    return md
end

end # module HTML2MD

#==============================================================================#
# End of file.
#==============================================================================#


text:
println("**********Testing FeatureVector**********")

facts("Creating a FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector()

  @fact fv1.total => 7
  @fact fv2.total => 0
  @fact fv1.map => dict1
  @fact typeof(fv1.map) => typeof(dict1)
  @fact typeof(fv2.map) => typeof(Dict{Any,Number}())
end

facts("Get value of FeatureVector given key") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)

  @fact fv1.total => 7
  @fact fv1["word"] => fv1.map["word"]
end

facts("Sanitize Dictionary values of 0") do
  dict = ["word" => 4, "another" => 3, "help"=>0]
  fv = FeatureVector(dict)
  @fact fv.total => 7
  fv["the"] = 1
  @fact fv.total => 8
  fv["the"] = 0
  @fact fv.total => 7

  @fact fv.map => ["word"=>4,"another"=>3]
  add!(fv,(-1*fv))
  @fact isempty(fv) => true
  @fact fv.total => 0
end

facts("Set value of FeatureVector given key") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  value = fv1["word"]

  fv1["word"] = 7

  @fact (fv1["word"] != dict1["word"]) => true
  @fact fv1["word"] => 7
  @fact dict1["word"] => 4
end

facts("Check to see if a key exists in a FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)

  @fact haskey(fv1, "word") => true
  @fact haskey(fv1, "another") => true
  @fact haskey(fv1, "nope") => false
end

facts("Get keys of FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)

  @fact keys(fv1) => Base.keys(fv1.map)
end

facts("Get values of FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)

  @fact values(fv1) => Base.values(fv1.map)
end

facts("Check isempty on a FeatureVector") do
  fv1 = FeatureVector()

  @fact isempty(fv1) => true

  fv1["the"] = 3

  @fact isempty(fv1) => false

end

facts("Make copy of a FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = copy(fv1)

  for key in keys(fv1)
          @fact fv1[key] => fv2[key]
  end
  @fact fv2 == fv1 => false
  @fact typeof(fv1) => typeof(fv2)
  @fact fv2.total => fv1.total
end

facts("Find common type of a FeatureVector and an empty FeatureVector") do
  dict1 = ["word" => 4, "∂" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector()

  @fact find_common_type(fv1,fv2) => (UTF8String,Int64)
  @fact find_common_type(fv2,fv1) => (UTF8String,Int64)
end

facts("Find common type of two FeatureVector") do
  dict1 = ["word" => 4, "another" => 3]
  dict2 = ["word" => .2, "∂" => .8]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector(dict2)

  @fact find_common_type(fv1,fv2) => (UTF8String,Float64)
end

facts("Adding two FeatureVectors") do
  dict1 = ["word" => 4, "another" => 3]
  dict2 = ["word" => 2, "∂" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector(dict2)
  dict3 = ["word"=>6,"∂"=>3,"another"=>3]
  fv3 = (fv1 + fv2)
  fv4 = (fv2 + fv1)
  for key in Base.keys(dict3)
    @fact fv3[key] => dict3[key]
  end
  for key in Base.keys(dict3)
    @fact fv4[key] => dict3[key]
  end
  @fact fv3.total => (fv1.total + fv2.total)
  @fact fv4.total => (fv1.total + fv2.total)
end

facts("Subtracting two FeatureVectors") do
  dict1 = ["word" => 4, "another" => 3]
  dict2 = ["word" => 2, "∂" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector(dict2)
  dict3 = ["word"=>2,"∂"=>-3,"another"=>3]
  fv3 = (fv1 - fv2)
  dict4 = ["word"=>-2,"∂"=>3,"another"=>-3]
  fv4 = (fv2 - fv1)
  for key in Base.keys(dict3)
    @fact fv3[key] => dict3[key]
  end
  for key in Base.keys(dict4)
    @fact fv4[key] => dict4[key]
  end
  @fact fv3.total => (fv1.total - fv2.total)
  @fact fv4.total => (fv2.total - fv1.total)
end

facts("Multiply a FeatureVector by a scalar") do
  dict1 = ["word" => 4, "another" => 3]
  dict2 = ["word" => .4, "∂" => .3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector(dict2)
  fv3 = FeatureVector()

  fv1 = fv1 * 2
  fv2 = fv2 * 2
  fv4 = fv3 * 2

  @fact fv1["word"] => 8
  @fact fv1["another"] => 6
  @fact fv1.total => 14
  @fact fv2["word"] => .8
  @fact fv2["∂"] => .6
  @fact fv2.total => 1.4
  @fact fv3 => fv4
  @fact fv3.total => 0
  @fact fv4.total => 0
end

facts("Divide a FeatureVector by a scalar") do
  dict1 = ["word" => 4, "another" => 3]
  dict2 = ["word" => .4, "∂" => .3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector(dict2)
  fv3 = FeatureVector()

  fv1 = fv1 / 2
  fv2 = fv2 / 2
  fv4 = fv3 / 2

  @fact fv1["word"] => 2
  @fact fv1["another"] => 1.5
  @fact fv1.total => 3.5
  @fact fv2["word"] => .2
  @fact fv2["∂"] => .15
  @fact fv2.total => .35
  @fact fv3 => fv4
  @fact fv3.total => 0
  @fact fv4.total => 0
end

facts("Rationalize a FeatureVector by a scalar") do
  dict1 = ["word" => 4, "another" => 3]
  fv1 = FeatureVector(dict1)
  fv2 = FeatureVector()

  fv1 = fv1 // 2
  fv3 = fv2 // 2

  @fact fv1["word"] => 2//1
  @fact fv1["another"] => 3//2
  @fact fv1.total => 7//2
  @fact fv2 => fv3
  @fact fv2.total => 0
  @fact fv3.total => 0
end

facts("Find dist_cos between two FeatureVectors") do
  fv = FeatureVector()
  fw = FeatureVector(Dict(["x"=>0, "y"=>0]))
  fv1 = FeatureVector(Dict(["x"=>1, "y"=>0]))
  fv2 = FeatureVector(Dict(["x"=>0, "y"=>1]))
  fv3 = FeatureVector(Dict(["x"=>1, "y"=>1]))
  fv4 = FeatureVector(Dict(["x"=>1.0, "y"=>sqrt(3)]))

  @fact dist_cos(fv3,fv3) => 0
  @fact isnan(dist_cos(fw,fv3)) => true
  @fact isnan(dist_cos(fv,fv1)) => true
  @fact 1 - cos(90*pi/180) - 1e-15 <= dist_cos(fv1,fv2) <= 1 - cos(90*pi/180) + 1e-15 => true
  @fact 1 - cos(45*pi/180) - 1e-15 <= dist_cos(fv2,fv3) <= 1 - cos(45*pi/180) + 1e-15 => true
  @fact 1 - cos(45*pi/180) - 1e-15 <= dist_cos(fv1,fv3) <= 1 - cos(45*pi/180) + 1e-15 => true
  @fact 1 - cos(60*pi/180) - 1e-15 <= dist_cos(fv1,fv4) <= 1 - cos(60*pi/180) + 1e-15 => true
  @fact 1 - cos(30*pi/180) - 1e-15 <= dist_cos(fv2,fv4) <= 1 - cos(30*pi/180) + 1e-15 => true
  @fact 1 - cos(0*pi/180) - 1e-15 <= dist_cos(fv1,fv1) <= 1 - cos(0*pi/180) + 1e-15 => true
end

facts("Find dist_zero between two FeatureVectors") do
  fv = FeatureVector()
  fv1 = FeatureVector(Dict(["x"=>2,"y"=>7]))
  fv2 = FeatureVector(Dict(["x"=>0,"y"=>4]))
  fv3 = FeatureVector(Dict(["x"=>3,"y"=>0]))
  fv4 = FeatureVector(Dict(["x"=>0,"y"=>0]))

  @fact dist_zero(fv,fv4) => 0
  @fact dist_zero(fv1,fv1) => 0
  @fact dist_zero(fv1,fv) == dist_zero(fv1,fv4) == 2 => true
  @fact dist_zero(fv2,fv) == dist_zero(fv3,fv4) == 1 => true
  @fact dist_zero(fv1,fv2) == dist_zero(fv1,fv3) == 1 => true
  @fact dist_zero(fv2,fv3) => 2
end

facts("Find dist_taxicab between two FeatureVectors") do
  fv = FeatureVector()
  fv1 = FeatureVector(Dict(["x"=>2,"y"=>2]))
  fv2 = FeatureVector(Dict(["x"=>0,"y"=>2]))
  fv3 = FeatureVector(Dict(["x"=>2,"y"=>0]))
  fv4 = FeatureVector(Dict(["x"=>0,"y"=>0]))

  @fact dist_taxicab(fv,fv4,false) => 0
  @fact dist_taxicab(fv1,fv1,false) => 0
  @fact dist_taxicab(fv1,fv,false) == dist_taxicab(fv1,fv4,false) == 4 => true
  @fact dist_taxicab(fv2,fv,false) == dist_taxicab(fv3,fv4,false) == 2 => true
  @fact dist_taxicab(fv1,fv2,false) == dist_taxicab(fv1,fv3,false) == 2 => true
  @fact dist_taxicab(fv2,fv3,false) => 4
end

facts("Find dist_euclidean between two FeatureVectors") do
  fv = FeatureVector()
  fv1 = FeatureVector(Dict(["x"=>2,"y"=>2]))
  fv2 = FeatureVector(Dict(["x"=>0,"y"=>2]))
  fv3 = FeatureVector(Dict(["x"=>2,"y"=>0]))
  fv4 = FeatureVector(Dict(["x"=>0,"y"=>0]))

  @fact dist_euclidean(fv,fv4,false) => 0
  @fact dist_euclidean(fv1,fv1,false) => 0
  @fact dist_euclidean(fv1,fv,false) == dist_euclidean(fv1,fv4,false) == sqrt(8) => true
  @fact dist_euclidean(fv2,fv,false) == dist_euclidean(fv3,fv4,false) == 2 => true
  @fact dist_euclidean(fv1,fv2,false) == dist_euclidean(fv1,fv3,false) == 2 => true
  @fact dist_euclidean(fv2,fv3,false) => sqrt(8)
end

facts("Find dist_infinite between two FeatureVectors") do
  fv = FeatureVector()
  fv1 = FeatureVector(Dict(["x"=>10,"y"=>5]))
  fv2 = FeatureVector(Dict(["x"=>0,"y"=>5]))
  fv3 = FeatureVector(Dict(["x"=>10,"y"=>0]))
  fv4 = FeatureVector(Dict(["x"=>0,"y"=>0]))

  @fact dist_infinite(fv,fv4,false) => 0
  @fact dist_infinite(fv1,fv,false) == dist_infinite(fv1,fv4,false) == 10 => true
  @fact dist_infinite(fv1,fv1,false) => 0
  @fact dist_infinite(fv2,fv,false) => 5
  @fact dist_infinite(fv3,fv4,false) => 10
  @fact dist_infinite(fv1,fv2,false) => 10
  @fact dist_infinite(fv1,fv3,false) => 5
  @fact dist_infinite(fv2,fv3,false) => 10
end