text:
using MacroTools: @q
import PyCall
using MLStyle
using Lazy

import SymPy
using SymPy: Sym, sympy, symbols

# stats = PyCall.pyimport_conda("sympy.stats", "sympy")
# import_from(stats)

export sym
sym(s::Symbol) = SymPy.symbols(s)
sym(s) = Base.convert(Sym, s)
function sym(expr::Expr)
    @match expr begin
        Expr(:call, f, args...) => :($f($(map(sym,args)...)))
        :($x[$j]) => begin
            j = symbols(:j, cls=sympy.Idx)
            x = sympy.IndexedBase(x)
            return x[j]
        end
        _ => begin
                 @show expr
                 error("sym: Argument type not implemented")
             end
    end
end

export symlogpdf
function symlogpdf(m::Model)
    result = @q begin
        ctx = Dict()

        ℓ = zero(Soss.Sym)
    end

    exprs = @as x m begin
        canonical(x)
        dropLines(x)
        x.body
        symlogpdf.(x)
    end

    append!(result.args, exprs)

    # result

    push!(result.args, :(ctx,ℓ))

    result

    # expandSums(ℓ)
end

function symlogpdf(st::Soss.Follows)
    d = st.rhs
    x = st.x
    :(ℓ += $(symlogpdf(d,x)))
end

function symlogpdf(st::Soss.Let)
    val = st.rhs
    x = st.x
    :(ctx[$(QuoteNode(x))] = Expr(:$, $val))
    # :($x = $val)
end

function symlogpdf(d::Expr, x::Symbol)
    @match d begin
        :(iid($n,$dist)) => begin
            j = symbols(:j, cls=sympy.Idx)
            dist = sym(dist)
            x = sympy.IndexedBase(x)
            n = sym(n)
            :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
        end

        :(For($f, 1:$n)) => begin
            n = sym(n)
            @match f begin
                :(($j,) -> begin $lineno; $dist end) => begin
                            j = symbols(j) # , cls=sympy.Idx)
                            # @show j
                            dist = sym(dist)
                            # @show dist
                            x = sympy.IndexedBase(x)
                            return :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
                end

                f => begin
                    @show f
                    error("symlogpdf: bad argument")
                end
            end

        end

        _ => :(logpdf($(sym(d)), $(sym(x))))
    end
end

function expandSums(s::Sym)
    s.args==() && return s
    if s.func==sympy.Sum
        s2 = expandSum(s)
        # Maybe it didn't simplify and we're done
        s2 == s && return s2
        # Or not, and we recurse
        s = s2
    end
    newargs = [expandSums(t) for t in s.args]
    s.func(newargs...)
end

function expandSum(s::Sym)
    # println("expandSum")
    # @show s
    # println()
    @assert s.func == sympy.Sum
    sfunc = s.args[1].func
    sargs = s.args[1].args
    limits = s.args[2]
    ix = limits.args[1]
    if sfunc == sympy.Add
        return sfunc([maybesum(t, limits) for t in sargs]...)
    elseif sfunc == sympy.Mul
        factors = sargs
        constants = [fac for fac in factors if !(ix in fac)]
        newconst = foldl(*,constants)
        newfacs = foldl(*,setdiff(factors, constants))
        newsum = sympy.Sum(newfacs, limits)
        return newconst * newsum
    else
        return s
    end

end

import Base.in
function Base.in(j::Sym, s::Sym)
    for t in s.args
        if j==t || in(j,t)
            return true
        end
    end
    return false
end

function maybesum(t::Sym, limits::Sym)
    # println("maybeSum")
    # @show t,limits
    # println()
    j = limits.args[1]
    thesum = sympy.Sum(t, limits)
    ifelse(j in t, thesum, thesum.doit())
end

# integrate(exp(ℓ), (sym(:μ), -oo, oo), (sym(:logσ),-oo,oo))
export marginal
function marginal(ℓ,v)
    f = ℓ.func
    f == sympy.Add || return ℓ
    newargs = filter(t -> sym(v) in t, collect(ℓ.args))
    foldl(+,newargs)
end

marginal(m::Model, v) = marginal(m |> symlogpdf, v)

# We should be able to reason about a marginal from its derivative
export dmarginal
function dmarginal(ℓ, v)
    @as x ℓ begin
        marginal(x,sym(v))
        diff(x, sym(v))
        expand(x)
        sympy.collect(x, sym(v))
    end
end

dmarginal(m::Model, v) = dmarginal(m |> symlogpdf, v)

# https://discourse.julialang.org/t/pyobjects-as-keys/26521/2
const symfuncs = Dict()
function __init__()
    merge!(symfuncs, Dict(
          sympy.log => log
        , sympy.Pow => :^
        , sympy.Abs => abs
        , sympy.Indexed => getindex
    ))
end

export codegen
function codegen(s::Sym)
    s.func == sympy.Add && begin
        @gensym add
        ex = @q begin
            $add = 0.0
        end
        for arg in s.args
            t = codegen(arg)
            push!(ex.args, :($add += $t))
        end
        push!(ex.args, add)
        # @show ex
        return ex
    end

    s.func == sympy.Mul && begin
        @gensym mul
        ex = @q begin
            $mul = 1.0
        end
        for arg in s.args
            t = codegen(arg)
            push!(ex.args, :($mul *= $t))
        end
        push!(ex.args, mul)
        return ex
    end

    # s.func == sympy.Sum && begin

    # end

    s.func ∈ keys(symfuncs) && begin
        # @show s
        @gensym symfunc
        argnames = gensym.("arg" .* string.(1:length(s.args)))
        argvals = codegen.(s.args)
        ex = @q begin end
        for (k,v) in zip(argnames, argvals)
            push!(ex.args, :($k = $v))
        end
        f = symfuncs[s.func]
        push!(ex.args, :($symfunc = $f($(argnames...))))
        push!(ex.args, symfunc)
        return ex
    end

    s.func == sympy.Sum && begin
        @gensym sum
        @gensym Δsum
        @gensym lo
        @gensym hi

        summand = codegen(s.args[1])
        (ix, ixlo, ixhi) = s.args[2].args

        ex = @q begin
            $sum = 0.0
            $lo = $(codegen(ixlo))
            $hi = $(codegen(ixhi))
            @inbounds @simd for $(codegen(ix)) = $lo:$hi
                $Δsum = $summand
                $sum += $Δsum
            end
            $sum
        end

        return ex
    end

    s.func == sympy.Symbol && return Symbol(string(s))
    s.func == sympy.Idx && return Symbol(string(s))
    s.func == sympy.IndexedBase && return Symbol(string(s))

    # @show s
    SymPy.is_real(s) && begin
        return SymPy.N(s)
    end

    @show s.func
    error("codegen")
end

codegen(s::AbstractFloat) = s

# export codegen
function codegen(m::Model)
    code = codegen(symlogpdf(m))
    unknowns = parameters(m) ∪ arguments(m)
    unkExpr = Expr(:tuple,unknowns...)
    @gensym logdensity
    result = @q begin
        function $logdensity(pars)
            @unpack $(unkExpr) = pars
            $code
        end
    end

    flatten(result)
end

# s = symlogpdf(normalModel).args[7].args[3]

# export fexpr
# fexpr = quote
#     f = function(μ,σ,x)
#         a = $(codegen(s))
#         return a
#     end
# end


text:
using WignerSymbols

function compute_all(max_angular)
    for l1 = 0:max_angular
        for l2 = 0:max_angular
            for l3 = 0:max_angular
                for m1 = -l1:l1
                    for m2 = -l2:l2
                        for m3 = -l3:l3
                            wigner3j(Float64, l1, l2, l3, m1, m2, m3)
                        end
                    end
                end
            end
        end
    end
end

for max_angular in [4, 8, 12, 16, 20]
    println("max_angular = $max_angular")
    # warmup & compile
    compute_all(max_angular)
    GC.gc()

    empty!(WignerSymbols.Wigner3j)
    @time compute_all(max_angular)
end

ext: jl

max_stars_repo_name: Luthaf/wigner

max_issues_repo_name: Luthaf/wigner

max_forks_repo_name: Luthaf/wigner


text:
VTKClasses = Dict(0=>vtk.vtkEmptyCell,
        1=>vtk.vtkVertex,
        2=>vtk.vtkPolyVertex,
        3=>vtk.vtkLine,
        4=>vtk.vtkPolyLine,
        5=>vtk.vtkTriangle,
        6=>vtk.vtkTriangleStrip,
        7=>vtk.vtkPolygon,
        8=>vtk.vtkPixel,
        9=>vtk.vtkQuad,
        10=>vtk.vtkTetra,
        11=>vtk.vtkVoxel,
        12=>vtk.vtkHexahedron,
        13=>vtk.vtkWedge,
        14=>vtk.vtkPyramid,
        15=>vtk.vtkPentagonalPrism,
        16=>vtk.vtkHexagonalPrism,
        21=>vtk.vtkQuadraticEdge,
        22=>vtk.vtkQuadraticTriangle,
        23=>vtk.vtkQuadraticQuad,
        24=>vtk.vtkQuadraticTetra,
        25=>vtk.vtkQuadraticHexahedron,
        26=>vtk.vtkQuadraticWedge,
        27=>vtk.vtkQuadraticPyramid,
        28=>vtk.vtkBiQuadraticQuad,
        29=>vtk.vtkTriQuadraticHexahedron,
        30=>vtk.vtkQuadraticLinearQuad,
        31=>vtk.vtkQuadraticLinearWedge,
        32=>vtk.vtkBiQuadraticQuadraticWedge,
        33=>vtk.vtkBiQuadraticQuadraticHexahedron,
        34=>vtk.vtkBiQuadraticTriangle,
        35=>vtk.vtkCubicLine,
        41=>vtk.vtkConvexPointSet,
        42=>vtk.vtkPolyhedron)

function PyVTK(dataset::VTKUnstructuredData)
    grid = vtk.vtkUnstructuredGrid()
    points = vtk.vtkPoints()

    if dim(dataset) == 2
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims([dataset.point_coords; zeros(num_of_points(dataset))'])))
    else
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims(dataset.point_coords)))
    end
    grid[:SetPoints](points)

    jl_cells = Int[]
    jl_cell_locs = Int[]
    k = 0
    for i in 1:num_of_cells(dataset)
        push!(jl_cell_locs, k)
        push!(jl_cells, length(dataset.cell_connectivity[i]))
        k += 1

        if dataset.cell_types[i] == 8
            pixel_cc = dataset.cell_connectivity[i]
            push!(jl_cells, pixel_cc[1]-1)
            push!(jl_cells, pixel_cc[2]-1)
            push!(jl_cells, pixel_cc[4]-1)
            push!(jl_cells, pixel_cc[3]-1)
            k += 4
        elseif dataset.cell_types[i] == 11
            voxel_cc = dataset.cell_connectivity[i]
            push!(jl_cells, voxel_cc[1]-1)
            push!(jl_cells, voxel_cc[2]-1)
            push!(jl_cells, voxel_cc[4]-1)
            push!(jl_cells, voxel_cc[3]-1)
            push!(jl_cells, voxel_cc[5]-1)
            push!(jl_cells, voxel_cc[6]-1)
            push!(jl_cells, voxel_cc[7]-1)
            push!(jl_cells, voxel_cc[8]-1)
            k += 8
        else
            for j in dataset.cell_connectivity[i]
                push!(jl_cells, j-1)
                k += 1
            end
        end
    end
    _cells = vtkns.numpy_to_vtk(np.array(jl_cells), 0, 12) #vtkTypeIdArray() of serialized cell connectivity
    ncells = num_of_cells(dataset)
    vtkCells = vtk.vtkCellArray()
    vtkCells[:SetCells](ncells, _cells)
    vtk_cell_types = vtkns.numpy_to_vtk(np.array(UInt8.(dataset.cell_types)), 0, 3) #vtkUnsignedCharArray() of cell types
    vtk_cell_locs = vtkns.numpy_to_vtk(np.array(jl_cell_locs), 0, 12) #vtkTypeIdArray() of cell locations
    grid[:SetCells](vtk_cell_types, vtk_cell_locs, vtkCells)

    for m in keys(dataset.point_data)
        vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_data[m]))
        vtk_data_array[:SetName](m)
        grid[:GetPointData]()[:AddArray](vtk_data_array)
    end

    for m in keys(dataset.cell_data)
        vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(dataset.cell_data[m]))
        vtk_data_array[:SetName](m)
        grid[:GetCellData]()[:AddArray](vtk_data_array)
    end

    return grid
end

function PyVTK(dataset::VTKPolyData)
    poly_pyclasses = [5,7,8,9]
    vert_pyclasses = [1,2]
    line_pyclasses = [3,4]
    strip_pyclasses = [6]

    new_inds = Int[]
    polydata = vtk.vtkPolyData()

    points = vtk.vtkPoints()
    if dim(dataset) == 2
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims([dataset.point_coords; zeros(num_of_points(dataset))'])))
    else
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims(dataset.point_coords)))
    end
    polydata[:SetPoints](points)

    jl_verts = Int[]
    jl_lines = Int[]
    jl_polys = Int[]
    jl_strips = Int[]
    vertinds = Int[]
    lineinds = Int[]
    polyinds = Int[]
    stripinds = Int[]

    for i in 1:num_of_cells(dataset)
        _cell_type = dataset.cell_types[i]
        if _cell_type ∈ vert_pyclasses
            push!(vertinds, i)
            push!(jl_verts, length(dataset.cell_connectivity[i]))
            for j in dataset.cell_connectivity[i]
                push!(jl_verts, j-1)
            end
        elseif _cell_type ∈ line_pyclasses
            push!(lineinds, i)
            push!(jl_lines, length(dataset.cell_connectivity[i]))
            for j in dataset.cell_connectivity[i]
                push!(jl_lines, j-1)
            end
        elseif _cell_type ∈ poly_pyclasses
            push!(polyinds, i)
            push!(jl_polys, length(dataset.cell_connectivity[i]))
            if _cell_type == 8
                pixel_cc = dataset.cell_connectivity[i]
                push!(jl_polys, pixel_cc[1]-1)
                push!(jl_polys, pixel_cc[2]-1)
                push!(jl_polys, pixel_cc[4]-1)
                push!(jl_polys, pixel_cc[3]-1)
            else
                for j in dataset.cell_connectivity[i]
                    push!(jl_polys, j-1)
                end
            end
        elseif _cell_type ∈ strip_pyclasses
            push!(stripinds, i)
            push!(jl_strips, length(dataset.cell_connectivity[i]))
            for j in dataset.cell_connectivity[i]
                push!(jl_strips, j-1)
            end
        else
            throw("Cell $i is not supported.")
        end
    end

    _verts = vtkns.numpy_to_vtk(np.array(jl_verts), 0, 12) #vtkTypeIdArray() of serialized cell connectivity
    nverts = length(vertinds)
    verts = vtk.vtkCellArray()
    verts[:SetCells](nverts, _verts)
    polydata[:SetVerts](verts)

    _lines = vtkns.numpy_to_vtk(np.array(jl_lines), 0, 12) #vtkTypeIdArray() of serialized cell connectivity
    nlines = length(lineinds)
    lines = vtk.vtkCellArray()
    lines[:SetCells](nlines, _lines)
    polydata[:SetLines](lines)

    _polys = vtkns.numpy_to_vtk(np.array(jl_polys), 0, 12) #vtkTypeIdArray() of serialized cell connectivity
    npolys = length(polyinds)
    polys = vtk.vtkCellArray()
    polys[:SetCells](npolys, _polys)
    polydata[:SetPolys](polys)

    _strips = vtkns.numpy_to_vtk(np.array(jl_strips), 0, 12) #vtkTypeIdArray() of serialized cell connectivity
    nstrips = length(stripinds)
    strips = vtk.vtkCellArray()
    strips[:SetCells](nstrips, _strips)
    polydata[:SetStrips](strips)

    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_data[m]))
        else
            vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_data[m]))
        end
        vtk_data_array[:SetName](m)
        polydata[:GetPointData]()[:AddArray](vtk_data_array)
    end

    new_cell_inds = chain(vertinds, lineinds, polyinds, stripinds)
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            new_cell_data = zeros(num_of_cells(dataset))
            i = 1
            for j in new_cell_inds
                new_cell_data[i] = dataset.cell_data[m][j]
                i += 1
            end
            vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(new_cell_data))
        else
            _size = size(dataset.cell_data[m])
            new_cell_data = zeros(_size...)
            i = 1
            for j in new_cell_inds
                new_cell_data[:,i] = dataset.cell_data[m][:,j]
                i += 1
            end
            vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(new_cell_data))
        end
        vtk_data_array[:SetName](m)
        polydata[:GetCellData]()[:AddArray](vtk_data_array)
    end

    return polydata
end

function PyVTK(dataset::VTKStructuredData)
    grid = vtk.vtkStructuredGrid()
    points = vtk.vtkPoints()

    _point_coords = reshape(dataset.point_coords, (dim(dataset), num_of_points(dataset)))
    if dim(dataset) == 2
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims([_point_coords; zeros(num_of_points(dataset))'])))
    else
        points[:SetData](vtkns.numpy_to_vtk(PyReverseDims(_point_coords)))
    end
    grid[:SetPoints](points)

    if _dim == 2
        grid[:SetDimensions](_point_extents..., 1)
    else
        grid[:SetDimensions](_point_extents...)
    end

    add_structured_data_to_grid!(grid, dataset)
    return grid
end

function PyVTK(dataset::VTKRectilinearData)
    grid = vtk.vtkRectilinearGrid()

    _dim = dim(dataset)
    _point_extents = extents(dataset)

    xCoords = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_coords[1]))
    yCoords = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_coords[2]))
    if _dim == 2
        zCoords = vtkns.numpy_to_vtk(PyReverseDims([0.]))
    else
        zCoords = vtkns.numpy_to_vtk(PyReverseDims(dataset.point_coords[3]))
    end

    if _dim == 2
        grid[:SetDimensions](_point_extents..., 1)
    else
        grid[:SetDimensions](_point_extents...)
    end

    grid[:SetXCoordinates](xCoords)
    grid[:SetYCoordinates](yCoords)
    grid[:SetZCoordinates](zCoords)

    add_structured_data_to_grid!(grid, dataset)
    return grid
end

function PyVTK(dataset::VTKUniformRectilinearData)
    grid = vtk.vtkImageData()

    _dim = dim(dataset)
    _point_extents = extents(dataset)

    if _dim == 2
        grid[:SetOrigin](dataset.origin..., 0)
        grid[:SetSpacing](dataset.spacing..., 0)
        grid[:SetDimensions](_point_extents..., 1)
    else
        grid[:SetOrigin](dataset.origin...)
        grid[:SetSpacing](dataset.spacing...)
        grid[:SetDimensions](_point_extents...)
    end

    add_structured_data_to_grid!(grid, dataset)
    return grid
end

function add_structured_data_to_grid!(grid::PyObject, dataset::AbstractVTKStructuredData)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            data_array = reshape(dataset.point_data[m], (1, num_of_points(dataset)))
        else
            data_array = reshape(dataset.point_data[m], (_var_dim, num_of_points(dataset)))
        end
        vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(data_array))
        vtk_data_array[:SetName](m)
        grid[:GetPointData]()[:AddArray](vtk_data_array)
    end

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            data_array = reshape(dataset.cell_data[m], (1, num_of_cells(dataset)))
        else
            data_array = reshape(dataset.cell_data[m], (_var_dim, num_of_cells(dataset)))
        end
        vtk_data_array = vtkns.numpy_to_vtk(PyReverseDims(data_array))
        vtk_data_array[:SetName](m)
        grid[:GetCellData]()[:AddArray](vtk_data_array)
    end
end

function PyVTK(dataset::VTKMultiblockData)
    mb = vtk.vtkMultiBlockDataSet()
    for i in 1:length(dataset)
        block = dataset[i]
        mb[:SetBlock](i-1, PyVTK(block))
    end
    return mb
end

function _VTKDataTypes(pyvtk::PyObject)
    if pyvtk[:GetClassName]() == "vtkMultiBlockDataSet"
        return extract_blocked_data(pyvtk)
    else
        return extract_simple_data(pyvtk)
    end
end

function write_vtk(pyvtk::PyObject, filepath_noext::String)
    classname = pyvtk[:GetClassName]()
    if classname == "vtkUnstructuredGrid"
        writer = vtk.vtkXMLUnstructuredGridWriter()
        writer[:SetFileName](string(filepath_noext, ".vtu"))
    elseif classname == "vtkStructuredGrid"
        writer = vtk.vtkXMLStructuredGridWriter()
        writer[:SetFileName](string(filepath_noext, ".vts"))
    elseif classname == "vtkRectilinearGrid"
        writer = vtk.vtkXMLRectilinearGridWriter()
        writer[:SetFileName](string(filepath_noext, ".vtr"))
    elseif classname == "vtkImageData"
        writer = vtk.vtkXMLImageDataWriter()
        writer[:SetFileName](string(filepath_noext, ".vti"))
    elseif classname == "vtkMultiBlockDataSet"
        writer = vtk.vtkXMLMultiBlockDataWriter()
        writer[:SetFileName](string(filepath_noext, ".vtm"))
    else
        throw("Type of VTK object is not supported.")
    end

    writer[:SetInputData](pyvtk)
    writer[:Update]()

    return filepath_noext
end


text:
# This file was generated, do not modify it.

using Pkg # hideall
Pkg.activate("_literate/A-composing-models/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

using MLJ
using PrettyPrinting
MLJ.color_off() # hide

KNNRegressor = @load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

scitype(X.age)

pipe = Pipeline(
    coercer = X -> coerce(X, :age=>Continuous),
    one_hot_encoder = OneHotEncoder(),
    transformed_target_model = TransformedTargetModel(
        model = KNNRegressor(K=3);
        target=UnivariateStandardizer()
    )
)

pipe.transformed_target_model.model.K = 2
pipe.one_hot_encoder.drop_last = true;

evaluate(
    pipe,
    X,
    height,
    resampling=Holdout(),
    measure=rms
) |> pprint

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using CuArrays
using CUDNN
# using Milk
include("../../src/core.jl")

using MLDatasets
using MLDataUtils
using PyPlot

include("model.jl")

function main()
    images, labels = CIFAR10.traindata()
    X_train = images
    # X_train = normalize(X_train)
    y_train = one_hot(maximum(labels) + 1, labels .+ 1)
    m = CIFAR10Model()
    @time m = fit!(m, X_train, y_train; cuda=true, n_epochs=100)

    # sample accuracy
    X = X_train[:,:,:,1:100] |> to_cuda
    y = y_train[:,1:100] |> to_cuda
    accuracy(predict(m, X), y)
end

# function normalize(X)
#     μ = mean(X, (1,2))
#     # σ2 = std(X, (1,2))
#     σ2 = 0.5
#     return (X .- μ) ./ σ2
# end


text:
#    Problem : GROUP A
#    *********
#	 A nondiagonal quartic test problem.
#
#	 This problem has an arrow-head type Hessian with a tridiagonal
#    central part and a border of width 1.
#    The Hessian is singular at the solution.
#
#    Origonal SIF Source:
#    A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
#    "Performance of a multi-frontal scheme for partially separable
#    optimization"
#    Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
#
#    NONDQUAR.SIF classification OUR2-AN-V-0
#
#    Number of variables is variable
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	n  = lastindex(x)
	fx =  (x[1]-x[2])^2 + (x[n-1] - x[n])^2
	for i in 1:n-2
      	fx += (x[i]+x[i+1]+x[n])^4
	end
    return fx
end

g! = (g, x) -> begin
	n    = lastindex(x)
	g[1] = 2(x[1] - x[2]) + 4(x[1] + x[2] + x[n])^3
	g[2] = -2(x[1] - x[2]) + 4(x[1] + x[2] + x[n])^3 + 4(x[2] + x[3] + x[n])^3
	gn   = 4(x[1] + x[2] + x[n])^3 +  4(x[2] + x[3] + x[n])^3
	for i in 3:n-2
		g[i] += 4(x[i] + x[i+1] + x[n])^3 + 4(x[i+1] + x[i+2] + x[n])^3
		gn   += 4(x[i] + x[i+1] + x[n])^3
	end
	g[n-1] = 2(x[n-1] - x[n]) + 4(x[n-2] + x[n-1] + x[n])^3
	g[n]   = gn - 2(x[n-1] - x[n])
	return g
end

fg! = (g, x) -> begin
	n    = lastindex(x)
	fx   =  (x[1]-x[2])^2 + (x[n-1] - x[n])^2 + (x[1]+x[2]+x[n])^4 + (x[2]+x[3]+x[n])^4
	g[1] = 2(x[1] - x[2]) + 4(x[1] + x[2] + x[n])^3
	g[2] = -2(x[1] - x[2]) + 4(x[1] + x[2] + x[n])^3 + 4(x[2] + x[3] + x[n])^3
	gn   = 4(x[1] + x[2] + x[n])^3 +  4(x[2] + x[3] + x[n])^3
	for i in 3:n-2
		g[i] += 4(x[i] + x[i+1] + x[n])^3 + 4(x[i+1] + x[i+2] + x[n])^3
		gn   += 4(x[i] + x[i+1] + x[n])^3
		fx += (x[i]+x[i+1]+x[n])^4
	end
	g[n-1] = 2(x[n-1] - x[n]) + 4(x[n-2] + x[n-1] + x[n])^3
	g[n]   = gn - 2(x[n-1] - x[n])
	return fx, g
end

init = (n::Int=5000) -> begin
	mod(n, 2) > 0 && @warn "NONDQUAR: number of variables must be even"
	q = max(1, div(n, 2))
	n = 2q
	x0 = [j % 2 == 0 ? -1.0 : 1.0 for j in 1:n]
    return n, x0
end

#TestSet["NONDQUAR"] = UncProgram("NONDQUAR", f, g!, fg!, init)


text:
using LinearAlgebra

# k_local(Elastic Modulus,thickness,height,length,number of nodes,row one, row two, column one, column two)

#creates local stiffness matrix and maps to global
function k_local(k,l,n,r_1,r_2,c_1,c_2)
    k_global = zeros(n,n)
    k_global[r_1,c_1] = k
    k_global[r_1,c_2] = -k
    k_global[r_2,c_1] = -k
    k_global[r_2,c_2] = k
    return k_global
end

# constants from problem statement
E = 29*10^6
l = 2
k1 = 50
k2 = 60
k3 = 55
# combining all globally mapped matrices
kG = k_local(k1,l,4,1,2,1,2)+
     k_local(k2,l,4,2,3,2,3)+
     k_local(k3,l,4,3,4,3,4)
kG[4,:] = [0 0 0 1] #adjustment for Reaction from fixed end

# Force Vector
f = Float64[60 0 0 0]
F = transpose(f)

# Displacement vector
u = kG\F
println("Nodal displacement vector (inches):")
display(u)

# Strain Function
ϵ(x2,x1,l) = (x2 - x1)/l

#Strain Vector
ϵM = [ϵ(u[2],u[1],l)
      ϵ(u[3],u[2],l)
      ϵ(u[4],u[3],l) ]
println("Element strain vector (inch/inch):")
display(ϵM)

# Stress Function
σ(E,ϵ) = E*ϵ

# Stress Vector
σM = [σ(E,ϵM[1])
      σ(E,ϵM[2])
      σ(E,ϵM[3])]
println("Element Stress vector (psi):")
display(σM)


text:
export get_opts, get_tol_norm
import Base.Iterators: Pairs

const SPQR_DEFAULTS = pairs((
                             get_details = 0,
                             repeatable = true,
                             k = 1,
                             tol = -1.0,
                             atol = 0.0,
                             rtol = 0.0,
                             tol_norm_type = 2,
                             nsvals_large = 1,
                             nsvals_small = 1,
                             implicit_null_space_basis = true,
                             start_with_A_transpose = false,
                             ssi_min_block = 3,
                             ssi_max_block = 10,
                             ssi_nblock_increment = 5,
                             ssi_min_iters = 3,
                             ssi_max_iters = 100,
                             ssi_convergence_factor = 0.1, # also 0.25 is often safe
                             ssp_min_iters = 4,
                             ssp_max_iters = 10,
                             ssp_convergence_factor = 0.1,
));

function get_opts(opts::Pairs, symbols::Symbol...)
    res = []
    for sym in symbols
        @assert sym in keys(SPQR_DEFAULTS)
        push!(res, get(opts, sym) do; getindex(SPQR_DEFAULTS, sym); end)
    end
    pushfirst!(res, merge_opts(opts, zip_named(symbols, res)))
    length(res) == 1 ? res[1] : Tuple(res)
end

function get_tol_norm(opts::Pairs, A::AbstractMatrix{T}) where T
    opts, tol = get_opts(opts, :tol)
    normest_A = -1.0
    if tol == -1.0
        opts, atol, rtol, tol_norm_type = get_opts(opts, :atol, :rtol, :tol_norm_type)
        if ( atol <= 0.0 && rtol <= 0.0 )
            rtol = max(size(A)...) * eps(real(T))
        end
        if rtol > 0.0
            normest_A = tol_norm_type != 2 ? opnorm(A, tol_norm_type) : normest(A, tol=0.01)
            xtol = rtol * normest_A
        end
        tol = max(atol, xtol)
        opts = merge_opts(opts, tol=tol, rtol=rtol)
    end
    opts, tol, normest_A
end

zip_named(vs, vv) = pairs(NamedTuple{Tuple(vs)}(Tuple(vv)))
kwlist(;kw...) = kw
merge_opts(x::Pairs; kw...) = kwlist(;x..., kw...)
merge_opts(x::Pairs, kw) = kwlist(;x..., kw...)


text:
function handle_msg(w::InputWidget, msg)
    if msg.content["data"]["method"] == "backbone"
        recv(w, msg.content["data"]["sync_data"]["value"])
    end
end

function handle_msg{T}(w::Button{T}, msg)
    try
        if msg.content["data"]["method"] == "custom" &&
            msg.content["data"]["content"]["event"] == "click"
            # click event occured
            recv(w, convert(T, w.value))
        end
    catch e
        warn(string("Couldn't handle Button message ", e))
    end
end

function handle_msg{view}(w::Options{view}, msg)
    try
        if msg.content["data"]["method"] == "backbone"
            key = msg.content["data"]["sync_data"]["value_name"]
            if haskey(w.options, key)
                recv(w, w.options[key])
            end
        end
    catch e
        warn(string("Couldn't handle ", view, " message ", e))
    end
end


text:
using Revise, Pkg, DataFrames

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using BloombergData

#need this since we are uing generated functions]
if @isdefined(FIRST_RUN)
  revise(BloombergData)
end

const FIRST_RUN = true

function BloombergDataScript()
  BloombergData.bbsession(:geoffpull, livepull=true)
end

@time BloombergDataScript()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#################### ModelChains ####################

#################### Constructors ####################

function ModelChains(c::Chains, m::Model)
  ModelChains(c.value, c.range, c.names, c.chains, m)
end

#################### Conversions ####################

Base.convert(::Type{Chains}, mc::ModelChains) =
  Chains(mc.value, mc.range, mc.names, mc.chains)

#################### Indexing ####################

function Base.getindex(mc::ModelChains, window, names, chains)
  c = getindex(convert(Chains, mc), window, names2inds(mc, names), chains)
  ModelChains(c, mc.model)
end

names2inds(mc::ModelChains, nodekey::Symbol) = names2inds(mc, [nodekey])

function names2inds(mc::ModelChains, nodekeys::Vector{Symbol})
  inds = Int[]
  missing = Symbol[]
  for key in nodekeys
    keyinds = indexin(names(mc.model, all=key), mc.names)
    0 in keyinds ? push!(missing, key) : append!(inds, keyinds)
  end
  if !isempty(missing)
    throw(ArgumentError(string(
      "chain values are missing for nodes : ",
      join(map(string, missing), ", ")
    )))
  end
  inds
end

function Base.keys(mc::ModelChains, ntype::Symbol, at...)
  values = Symbol[]
  m = mc.model
  nodekeys = ntype == :dependent ?
    keys(m, :dependent) :
    intersect(keys(m, ntype, at...), keys(m, :dependent))
  for key in nodekeys
    all(name -> name in mc.names, names(m, all=key)) && push!(values, key)
  end
  values
end

#################### Auxilliary Functions ####################

function link(c::ModelChains)
  cc = copy(c.value)
  inds_queue = 1:length(c.names)
  for key in intersect(keys(c.model, :monitor), keys(c.model, :stochastic))
    node = c.model[key]
    inds = findall(in(names(node)), c.names)
    if !isempty(inds)
      f(x) = unlist(node, relist(node, x), true)
      cc[:, inds, :] = mapslices(f, cc[:, inds, :], 2)
      inds_queue = setdiff(inds_queue, inds)
    end
  end
  for j in inds_queue
    x = cc[:, j, :]
    if minimum(x) > 0.0
      cc[:, j, :] = maximum(x) < 1.0 ? StatsFuns.logit.(x) : log.(x)
    end
  end
  cc
end


text:
################################################################################
# forestGrowthFitting Problem
#
# How fast are French forests growing? How much timber can they provide in one hectare?
# Surely forests provide multiple ecosystem services, but in this (simplified) exercise we focus on the above questions by looking at the so-called "raw data" from the French National Forest Inventory, both in terms of individual trees and in terms of inventoried plots, to fit a generic growth model of the forest stands in terms of volumes with respect to the age of the trees.
#

# ### Environment set-up and data loading

# 1) Start by setting the working directory to the directory of this file and activate it. If you have the provided `Manifest.toml` file in the directory, just run `Pkg.instantiate()`, otherwise manually add the packages Pipe, HTTP, CSV, DataFrames, LsqFit, StatsPlots.

# 2) Load the packages Pipe, HTTP, CSV, DataFrames, LsqFit, StatsPlots

# 3) Load from internet or from local files the following data:

ltURL     = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/blob/main/lessonsMaterial/02_JULIA2/forestExercise/data/arbres_foret_2012.csv?raw=true" # live individual trees data
dtURL     = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/blob/main/lessonsMaterial/02_JULIA2/forestExercise/data/arbres_morts_foret_2012.csv?raw=true" # dead individual trees data
pointsURL = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/blob/main/lessonsMaterial/02_JULIA2/forestExercise/data/placettes_foret_2012.csv?raw=true" # plot level data
docURL    = "https://github.com/sylvaticus/IntroSPMLJuliaCourse/blob/main/lessonsMaterial/02_JULIA2/forestExercise/data/documentation_2012.csv?raw=true" # optional, needed for the species label

# If you choosen to download the data from internet, you can make for each of the dataset a `@pipe` macro starting with `HTTP.get(URL).body`, continuing the pipe with `CSV.File(_)` and end the pipe with a DataFrame object.

# 4) These datasets have many variable we are not using in this exercise.
# Out of all the variables, select only for the `lt` and `dt` dataframes the columns "idp" (pixel id), "c13" (circumference at 1.30 meters) and "v" (tree's volume). Then vertical concatenate the two dataset in an overall `trees` dataset.
# For the `points` dataset, select only the variables "idp" (pixel id), "esspre" (code of the main forest species in the stand) and "cac" (age class).

# 5) As the French inventory system is based on a concentric sample method (small trees are sampled on a small area (6 metres radius), intermediate trees on a concentric area of 9 metres and only large trees (with a circonference larger than 117.5 cm) are sampled on a concentric area of 15 metres of radius), define the following function to compute the contribution of each tree to the volume per hectare:

"""
    vHaContribution(volume,circonference)

Return the contribution in terms of m³/ha of the tree.

The French inventory system is based on a concentric sample method: small trees are sampled on a small area (6 metres radius), intermediate trees on a concentric area of 9 metres and only large trees (with a circonference larger than 117.5 cm) are sampled on a concentric area of 15 metres of radius.
This function normalise the contribution of each tree to m³/ha.
"""
function vHaContribution(v,c13)
    if c13 < 70.5
        return v/(6^2*pi/(100*100))
    elseif c13 < 117.5
        return v/(9^2*pi/(100*100))
    else
        return v/(15^2*pi/(100*100))
    end
end
# Use the above function to compute `trees.vHa` based on `trees.v` and `trees.c13`.

# 6) Aggregate the `trees` dataframe by the `idp` column to retrieve the sum of `vHa` and the number of trees for each point, calling these two columns `vHa` and `ntrees`.

# 7) Join the output of step 6 (the trees dataframe aggregated "by point") with the original points dataframe using the column `idp`.

# 8) Use boolean selection to apply the following filters:
filter_nTrees           = points.ntrees .> 5
filter_IHaveAgeClass    = .! in.(points.cac,Ref(["AA","NR"]))
filter_IHaveMainSpecies = .! ismissing.(points.esspre)
filter_overall          = filter_nTrees .&& filter_IHaveAgeClass .&& filter_IHaveMainSpecies

# 9) Run the following line to parse the age class (originally as a string indicating the 5-ages group) to an integer and compute the mid-range of the class in years. For example, class "02" will become 7.5 years.
points.cac              = (parse.(Int64,points.cac) .- 1 ) .* 5 .+ 2.5

# 10) Define the following logistic model of the growth relation with respect to the age with 3 parameters and make its vectorised form:
logisticModel(age,parameters) = parameters[1]/(1+exp(-parameters[2] * (age-parameters[3]) ))
logisticModelVec(age,parameters) = # .... complete

# 11) Set `initialParameters` to 1000,0.05 and 50 respectivelly

# 12) Perform the fittin of the model using the function `curve_fit(model,X,Y,initial parameters)` and obtain the fitted parameter fitobject.param

# 13) Compute the standard error for each estimated parameter and the confidence interval at 10% significance level

# 14) Plot a chart of fitted (y) by stand age (x) (i.e. the logisticModel with the given parameters)

# 15) Add to the plot a scatter chart of the actual observed VHa

# 16) [OPTIONAL] Look at the growth curves of individual species.
# Try to perform the above analysis for individual species, for example plot the fitted curves for the 5 most common species

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using Documenter
using Calendars

makedocs(
    modules = [Calendars],
    sitename = "Calendars.jl",
    clean = true,
    checkdocs = :none,
    doctest = false,
    pages = [
        "Functions" => "index.md",
        "References" => "References.md"
    ]
)

# "Introduction" => "Calendars.md",
#deploydocs(
#    repo = "github.com/PeterLuschny/Calendars.jl.git";
#    push_preview = true
#)


text:
using ProgressMeter

"""
Canonicalize a GitHub repository URL
"""
function canonicalize_source_url(url)
    repo_regex = r"(https:\/\/)?github.com\/([^\/]+)\/([^\/]+)\/?$"

    m = match(repo_regex, url)
    if m !== nothing
        _, user, repo = m.captures
        if !endswith(repo, ".git")
            return "https://github.com/$user/$repo.git"
        end
    end
    url
end

struct GitTransferProgress
    total_objects::Cuint
    indexed_objects::Cuint
    received_objects::Cuint
    local_objects::Cuint
    total_deltas::Cuint
    indexed_deltas::Cuint
    received_bytes::Csize_t
end

function transfer_progress(progress::Ptr{GitTransferProgress}, p::Any)
    progress = unsafe_load(progress)
    p.n = progress.total_objects
    if progress.total_deltas != 0
        p.desc = "Resolving Deltas: "
        p.n = progress.total_deltas
        update!(p, Int(max(1, progress.indexed_deltas)))
    else
        update!(p, Int(max(1, progress.received_objects)))
    end
    return Cint(0)
end

macro compat_gc_preserve(args...)
    esc(Expr(:macrocall, Expr(:., :GC, QuoteNode(Symbol("@preserve"))), args...))
end

function clone(url, source_path)
    p = Progress(0, 1, "Cloning: ")
    @compat_gc_preserve p begin
        @static if VERSION >= v"0.7-"
            callbacks = LibGit2.RemoteCallbacks(transfer_progress=@cfunction(transfer_progress, Cint, (Ptr{GitTransferProgress}, Any)),
                payload = p)
        else
            callbacks = LibGit2.RemoteCallbacks(transfer_progress=cfunction(transfer_progress, Cint, Tuple{Ptr{GitTransferProgress}, Any}),
                payload = pointer_from_objref(p))
        end
        fetch_opts = LibGit2.FetchOptions(callbacks = callbacks)
        clone_opts = LibGit2.CloneOptions(fetch_opts=fetch_opts, bare = Cint(true))
        return LibGit2.clone(url, source_path, clone_opts)
    end
end

"""
    download_source(state::WizardState)

Ask the user where the source code is coming from, then download and record the
relevant parameters, returning the source `url`, the local `path` it is stored
at after download, and a `hash` identifying the version of the code. In the
case of a `git` source URL, the `hash` will be a git treeish identifying the
exact commit used to build the code, in the case of a tarball, it is the
`sha256` hash of the tarball itself.
"""
function download_source(state::WizardState)
    # First, ask the user where this is all coming from
    msg = replace(strip("""
    Please enter a URL (git repository or compressed archive) containing the
    source code to build:
    """), "\n" => " ")
    print(state.outs, msg, "\n> ")
    entered_url = readline(state.ins)
    println(state.outs)

    url = canonicalize_source_url(entered_url)
    if url != entered_url
        print(state.outs, "The entered URL has been canonicalized to\n")
        printstyled(state.outs, url, bold=true)
        println(state.outs)
        println(state.outs)
    end

    # Record the source path and the source hash
    source_path = joinpath(state.workspace, basename(url))
    local source_hash

    if endswith(url, ".git")
        # Clone the URL, record the current gitsha for the given branch
        repo = clone(url, source_path)

        msg = "You have selected a git repository. Please enter a branch, commit or tag to use.\n" *
        "Please note that for reproducability, the exact commit will be recorded, \n" *
        "so updates to the remote resource will not be used automatically; \n" *
        "you will have to manually update the recorded commit."
        print(state.outs, msg, "\n> ")
        treeish = readline(state.ins)
        println(state.outs)

        obj = try
            LibGit2.GitObject(repo, treeish)
        catch
            LibGit2.GitObject(repo, "origin/$treeish")
        end
        source_hash = LibGit2.string(LibGit2.GitHash(obj))

        # Tell the user what we recorded the current commit as
        print(state.outs, "Recorded as ")
        printstyled(state.outs, source_hash, bold=true)
        println(state.outs)
        close(repo)
    else
        # Download the source tarball
        source_path = joinpath(state.workspace, basename(url))

        if isfile(source_path)
            name, ext = splitext(basename(source_path))
            n = 1
            while isfile(joinpath(state.workspace, "$(name)_$n$ext"))
                n += 1
            end
            source_path = joinpath(state.workspace, "$(name)_$n$ext")
        end

        download_cmd = gen_download_cmd(url, source_path)
        oc = OutputCollector(download_cmd; verbose=true, tee_stream=state.outs)
        try
            if !wait(oc)
                error()
            end
        catch
            error("Could not download $(url) to $(state.workspace)")
        end

        # Save the source hash
        open(source_path) do file
            source_hash = bytes2hex(BinaryProvider.sha256(file))
        end
    end

    # Spit back the url, local path and source hash
    return url, source_path, source_hash
end

"""
    step1(state::WizardState)

It all starts with a single step, the unabashed ambition to leave your current
stability and engage with the universe on a quest to create something new, and
beautiful and unforseen.  It all ends with compiler errors.

This step selets the relevant platform(s) for the built binaries.
"""
function step1(state::WizardState)
    # Select a platform
    msg = "\t\t\t# Step 1: Select your platforms\n\n"
    printstyled(state.outs, msg, bold=true)
    terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
    platform_select = request(terminal,
        "Make a platform selection",
        RadioMenu([
            "All supported architectures",
            "Specific operating system",
            "Specific architecture",
            "Custom",
        ])
    )
    println(state.outs)

    # Set `state.platforms` accordingly
    result = nothing
    if platform_select == 1
        state.platforms = supported_platforms()
    elseif platform_select == 2
        oses = sort(unique(map(typeof, supported_platforms())), by = repr)
        while true
            result = request(terminal,
                "Select operating systems",
                MultiSelectMenu(map(repr, oses))
            )
            result = map(x->oses[x], collect(result))
            if isempty(result)
                println("Must select at least one operating system (CTRL-C to cancel)")
            else
                break
            end
        end
        state.platforms = collect(filter(x->typeof(x) in result, supported_platforms()))
    elseif platform_select == 3
        arches = sort(unique(map(arch, supported_platforms())), by = repr)
        while true
            result = request(terminal,
                "Select architectures",
                MultiSelectMenu(map(repr, arches))
            )
            result = map(x->arches[x], collect(result))
            if isempty(result)
                println("Must select at least one architecture (CTRL-C to cancel)")
            else
                break
            end
        end
        state.platforms = collect(filter(x->arch(x) in result, supported_platforms()))
    elseif platform_select == 4
        platfs = supported_platforms()
        while true
            result = request(terminal,
                "Select platforms",
                MultiSelectMenu(map(repr, platfs))
            )
            if isempty(result)
                println("Must select at least one platform (CTRL-C to cancel)")
            else
                break
            end
        end
        state.platforms = collect(map(x->platfs[x], collect(result)))
    else
        error("Somehow platform_select was not a valid choice!")
    end
    state.platforms = [typeof(p)(arch(p); libc=libc(p), call_abi=call_abi(p), compiler_abi=CompilerABI(:gcc8)) for p in state.platforms]
    println(state.outs)
end

"""
Canonicalize URL to a file within a GitHub repo
"""
function canonicalize_file_url(url)
    blob_regex = r"(https:\/\/)?github.com\/([^\/]+)\/([^\/]+)\/blob\/([^\/]+)\/(.+)"
    m = match(blob_regex, url)
    if m !== nothing
        _, user, repo, ref, filepath = m.captures
        if length(ref) != 40 || !all(c->isnumeric(c) || c in 'a':'f' || c in 'A':'F', ref)
            # Ask github to resolve this ref for us
            ref = GitHub.reference(GitHub.Repo("$user/$repo"), "heads/$ref").object["sha"]
        end
        return "https://raw.githubusercontent.com/$user/$repo/$ref/$(filepath)"
    end
    url
end

function obtain_binary_deps(state::WizardState)
    msg = "\t\t\t# Step 2b: Obtain binary dependencies (if any)\n\n"
    printstyled(state.outs, msg, bold=true)

    q = "Do you require any (binary) dependencies? "
    if yn_prompt(state, q, :n) == :y
        empty!(state.dependencies)
        terminal = TTYTerminal("xterm", state.ins, state.outs, state.outs)
        while true
            bindep_select = request(terminal,
                "How would you like to specify this dependency?",
                RadioMenu([
                    "Provide remote build.jl file",
                    "Paste in a build.jl file",
                    "Never mind",
                ])
            )
            println(state.outs)

            if bindep_select == 1
                println(state.outs, "Enter the URL to use: ")
                print(state.outs, "> ")
                url = readline(state.ins)
                println(state.outs)
                canon_url = canonicalize_file_url(url)
                if url != canon_url
                    print(state.outs, "The entered URL has been canonicalized to\n")
                    printstyled(state.outs, canon_url, bold=true)
                    println(state.outs)
                    println(state.outs)
                end
                push!(state.dependencies, RemoteBuildDependency(canon_url,
                    String(HTTP.get(canon_url).body)))
            elseif bindep_select == 2
                println(state.outs, "Please provide the build.jl file. Press ^D when you're done")
                script = String(read(state.ins))
                Base.reseteof(terminal)
                push!(state.dependencies, InlineBuildDependency(script))
            elseif bindep_select == 3
                break
            end

            q = "Would you like to provide additional dependencies? "
            if yn_prompt(state, q, :n) != :y
                println(state.outs)
                break
            end
        end
    end
end

function obtain_source(state::WizardState)
    msg = "\t\t\t# Step 2a: Obtain the source code\n\n"
    printstyled(state.outs, msg, bold=true)

    # Create the workspace that we'll stash everything within
    state.workspace = tempname()
    mkpath(state.workspace)

    # These are the metadata we need to know about all the sources we'll be
    # building over the course of this journey we're on together.
    state.source_urls = String[]
    state.source_files = String[]
    state.source_hashes = String[]

    while true
        url, file, hash = download_source(state)
        push!(state.source_urls, url)
        push!(state.source_files, file)
        push!(state.source_hashes, hash)
        println(state.outs)

        q = "Would you like to download additional sources? "
        if yn_prompt(state, q, :n) != :y
            break
        end
    end
    println(state.outs)
end

"""
    step2(state::WizardState)

This step obtains the source code to be built and required binary dependencies.
"""
function step2(state::WizardState)
    obtain_source(state)
    obtain_binary_deps(state)
end


text:
"""
    clusterkcenters(X::AbstractMatrix, kcluster::Int; nReplicates::Int=10) -> F

Perform clustering with K-center algorithm. Input data `X` should belong to AbstractMatrix type
and its columns corresponds to variables, rows are frames. Also, the number of clusters `kcluster` should be specified.

Returns a NamedTuple object `F` which contains cluster index for each sample in `F.indexOfCluster`,
the coordinates of cluster centers in `F.center`, the indices of cluster centers in `F.indexOfCenter`,
and distances of samples from the nearest centers in `F.distanceFromCenter`.

#  Example
```julia-repl
julia> using MDToolbox, Plots
julia> X = rand(1000, 2)
julia> F = clusterkcenters(X, 3)
julia> scatter(X[:, 1], X[:, 2], c=F.indexOfCluster)
```
# References
```
This function uses the method described in
[1] S. Dasgupta and P. M. Long, J. Comput. Syst. Sci. 70, 555 (2005).
[2] J. Sun, Y. Yao, X. Huang, V. Pande, G. Carlsson, and L. J. Guibas, Learning 24, 2 (2009).
```
"""
function clusterkcenters(t::AbstractMatrix, kcluster::Int; nReplicates::Int=10)
    nframe = size(t, 1)
    #ta2, com = decenter(ta)
    dim = size(t, 2)

    indexOfCenter = zeros(Int64, kcluster)
    indexOfCluster = zeros(Int64, nframe)
    distanceFromCenter = zeros(Float64, nframe, 1)

    indexOfCenter_out = similar(indexOfCenter)
    indexOfCluster_out = similar(indexOfCluster)
    distanceFromCenter_out = similar(distanceFromCenter)
    center_out = similar(t, kcluster, dim)
    distanceMax_out = Inf64

    for ireplica = 1:nReplicates
        # create the center of the 1st cluster
        indexOfCenter[1] = rand(1:nframe)
        # first, all points belong to the 1st cluster
        indexOfCluster .= ones(Int64, nframe)
        # distances of the points from the 1st center
        ref = t[indexOfCenter[1]:indexOfCenter[1], :]
        distanceFromCenter .= sum((ref .- t).^2, dims=2)

        for i = 2:kcluster
            indexOfCenter[i] = argmax(distanceFromCenter[:])
            ref = t[indexOfCenter[i]:indexOfCenter[i], :]
            dist = sum((ref .- t).^2, dims=2)
            index = dist[:] .< distanceFromCenter[:]
            if any(index)
                # updated if the dist to a new cluster is smaller than the previous one
                distanceFromCenter[index] .= dist[index]
                indexOfCluster[index] .= i
            end
        end

        distanceMax = maximum(distanceFromCenter)
        if (ireplica == 1) | (distanceMax < distanceMax_out)
          distanceMax_out = distanceMax
          indexOfCenter_out = indexOfCenter
          indexOfCluster_out = indexOfCluster
          center_out = t[indexOfCenter, :]
          distanceFromCenter_out .= distanceFromCenter
        end
        Printf.@printf("%d iteration  distance_max = %f  kcluster = %d\n", ireplica, sqrt(distanceMax_out), kcluster)
    end
    distanceFromCenter_out .= sqrt.(distanceFromCenter_out)

    return (indexOfCluster=indexOfCluster_out, center=center_out, distanceFromCenter=distanceFromCenter, indexOfCenter=indexOfCenter_out)
end

"""
    clusterkcenters(ta::AbstractMatrix, kcluster::Int; nReplicates::Int=10) -> F

Perform clustering with K-center algorithm for a TrjArray variable `ta`.
"""
function clusterkcenters(ta::TrjArray, kcluster::Int; nReplicates::Int=10)
    X = ta.xyz
    F = clusterkcenters(X, kcluster, nReplicates=nReplicates)
    center_ta = ta[F.indexOfCenter, :]
    return (indexOfCluster=F.indexOfCluster, center=center_ta, distanceFromCenter=F.distanceFromCenter, indexOfCenter=F.indexOfCenter)
end

"""
    clustercutoff(X::AbstractMatrix, rcut=10.0) -> F

Perform clustering using the given cutoff radius.

#  Example
```julia-repl
julia> using MDToolbox, Plots
julia> X = rand(1000, 2)
julia> F = clustercutoff(X, 0.5)
julia> scatter(X[:, 1], X[:, 2], c=F.indexOfCluster)
```
# References
```
This function uses the method described in
[1] X. Daura, K. Gademann, B. Jaun, D. Seebach, W. F. van Gunsteren, and A. E. Mark, Angewandte 38, 236–240 (1999).
```
"""
function clustercutoff(t::AbstractMatrix, rcut=10.0)
    nframe = size(t, 1)
    rcut2 = rcut^2

    indexOfCluster = zeros(Int64, nframe)
    id = indexOfCluster .== 0

    center = similar(t)
    k = 1

    while any(id)
        x = view(t, id, :)
        index = view(indexOfCluster, id)

        n_max = 0
        iframe_max = 0
        for iframe = 1:size(x, 1)
            dist = sum((x[iframe:iframe, :] .- x).^2, dims=2)
            n = sum(dist .< rcut2)
            if n_max < n
                n_max = n
                iframe_max = iframe
            end
        end

        center[k:k, :] .= x[iframe_max:iframe_max, :]
        dist = sum((center[k:k, :] .- x).^2, dims=2)
        index[dist[:] .< rcut2] .= k
        k += 1

        id = indexOfCluster .== 0
    end

    return (indexOfCluster=indexOfCluster, center=center[1:(k-1), :])
end

"""
    compute_cov(ta::AbstractMatrix, lagtime::Int=0) -> cov

Compute a variance-covariance or time-lagged covariance matrix from input data `X`
Input data `X` should belong to AbstractMatrix type and its columns corresponds to variables, rows are frames.
Optional input is the lagtime=`lagtime` for the calculation of the covariance matrix (default is `lagtime=0`).

#  Example
```julia-repl
julia> X = rand(1000, 100)
julia> cov = compute_cov(X)
```
"""
function compute_cov(X::AbstractMatrix; lagtime::Int=0)
    nframe = size(X, 1)
    X_centerized = X .- mean(X, dims=1)
    cov = X_centerized[1:(end-lagtime), :]' * X_centerized[(1+lagtime):end, :] ./ (nframe - 1)
end

"""
    rsvd(X::AbstractMatrix; k::Number=10) -> F

Perform the randomized SVD for input data `X`.
Input data `X` should belong to AbstractMatrix type and its columns corresponds to variables, rows are frames.
Users can specify the dimension `k` of subspace onto which the data is randomly projected (by default `k=10`).

Returns a NamedTuple object `F` whose members are same as usual SVD, `F.V`, `F.S`, `F.U`.

#  Example
```julia-repl
julia> X = randn(rand(1000, 10))
julia> F = rsvd(X)
```

# References
```
Halko, Nathan, Per-Gunnar Martinsson, and Joel A. Tropp.
"Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions."
SIAM review 53.2 (2011): 217-288.
```
"""
function rsvd(A::AbstractMatrix, k::Number=10)
    T = eltype(A)
    m, n = size(A)
    l = k + 2
    @assert l < n
    G = randn(T, n, l)
    #G .= G ./ sqrt.(sum(G.^2, dims=1))
    H = A * G
    @assert m > l
    Q = Matrix(qr!(H).Q)
    @assert m == size(Q, 1)
    @assert l == size(Q, 2)
    T = A' * Q
    V, σ, W = svd(T)
    U = Q * W
    λ = σ .* σ ./ m
    return (V=V[:,1:k], S=λ[1:k], U=U[:,1:k])
end

"""
    pca(X::AbstractMatrix; k=dimension) -> F

Perform principal component analysis (PCA). PCA captures degrees of freedom which have the largest variances.
Input data `X` should belong to AbstractMatrix type and its columns corresponds to variables, rows are frames.

Returns a NamedTuple object `F` which contains the principal components in `F.projection`,
the prncipal modes in the columns of the matrix `F.mode`, and the variances of principal components in `F.variance`.

If `k=dimension` is specified, the randomized SVD is used for reducing memory.
This algorithm first project the data into a randomly selected k+2-dimensional space,
then PCA is performed in the projected data. See the references for details.
Note that if the dimension of `X` is larger than 5000, the randomized SVD is forcibly used with `k=1000`.

#  Example
```julia-repl
julia> using MDToolbox, Plots, Statistics
julia> X = cumsum(rand(1000, 10))
julia> F = pca(X)
julia> plot(F.projection[:, 1], F.projection[:, 2])
```

# References
```
Halko, N., Martinsson, P.-G., Shkolnisky, Y. & Tygert, M.
An Algorithm for the Principal Component Analysis of Large Data Sets.
SIAM J. Sci. Comput. 33, 2580–2594 (2011).
```
"""
function pca(X::AbstractMatrix; k=nothing)
    is_randomized_pca = false
    if !isnothing(k)
        is_randomized_pca = true
        @printf "Warning: a randomized SVD approximation is used for PCA with reduced dimension k = %d\n" k
    end
    dim = size(X, 2)
    if !is_randomized_pca & (dim > 5000)
        k = 1000
        is_randomized_pca = true
        @printf "Warning: the dimension of input data is too large (dim > 5000)\n"
        @printf "Warning: randomized SVD approximation is used for PCA with reduced dimension k = %d\n" k
    end
    if is_randomized_pca
        nframe = size(X, 1)
        X_centerized = X .- mean(X, dims=1)
        F = rsvd(X_centerized, k)
        variance = F.S.^2 ./ (nframe - 1)
        mode = F.V
        projection = X_centerized * mode
    else
        # eigendecomposition of covariance matrix
        covar = compute_cov(X)
        F = eigen(covar, sortby = x -> -x)
        #lambda = F.values[end:-1:1]
        #W = F.vectors[:, end:-1:1]
        variance = F.values
        mode = F.vectors
        # projection
        X_centerized = X .- mean(X, dims=1)
        projection = X_centerized * mode
    end

    return (projection=projection, mode=mode, variance=variance)
end

"""
    tica(X::AbstractMatrix, lagtime::Int=1) -> F

This routine performs time-structure based Independent Component Analysis (tICA).
tICA captures degrees of freedom which are most important in the sense that their motions are very slow.
`X` belongs to AbstractMatrix type and its columns corresponds to variables, rows are frames.
User should specify the lagtime for the calculation of the time-lagged covariance matrix.
Returns a NamedTuple object `F` which contains independent components in `F.projection`,
the independent modes in the columns of the matrix `F.mode`, and the eigenvalues in `F.variance`.

If the dimension of `X` is larger than 5000, the randomized SVD approximation is forcibly used with `k=1000` to reduce memory.

#  Example
```julia-repl
julia> using MDToolbox, Plots, Statistics
julia> X = cumsum(rand(1000, 10))
julia> F = tica(X, 10)
julia> plot(F.projection[:, 1], F.projection[:, 2])
```

# References
```
Naritomi, Y. & Fuchigami, S.
Slow dynamics in protein fluctuations revealed by time-structure based independent component analysis:
The case of domain motions.
The Journal of Chemical Physics 134, 065101 (2011).
```
"""
function tica(X::AbstractMatrix, lagtime::Int=1)
    # standard and time-lagged covariance matrices
    covar0 = compute_cov(X, lagtime=0)
    covar  = compute_cov(X, lagtime=lagtime)

    # symmetrize the time-lagged covariance matrix
    covar .= 0.5 .* (covar .+ covar')

    # calc pseudo-inverse of covar0
    #covar0_inv = pinv(covar0)

    # remove degeneracy for solving the generalized eigenvalue problem
    #F = eigen(covar0, sortby = x -> -x)
    F = pca(X)
    pvariance = F.variance
    pmode = F.mode
    index = pvariance .> 10^(-6)
    pmode = pmode[:, index]
    covar0 = pmode'*covar0*pmode
    covar = pmode'*covar*pmode

    # solve the generalized eigenvalue problem
    F = eigen(covar, covar0, sortby = x -> -x)
    variance = F.values
    mode = F.vectors

    # projection
    mode = pmode * mode
    X_centerized = X .- mean(X, dims=1)
    projection = X_centerized * mode

    # normalize mode vectors
    fac = sqrt.(sum(mode.^2, dims=1))
    mode .= mode ./ fac

    return (projection=projection, mode=mode, variance=variance)
end


text:
Vtr = (ℂ^2,
        (ℂ^3)',
        ℂ^4,
        ℂ^5,
        (ℂ^6)')
Vℤ₂ = (ℂ[ℤ₂](0=>1, 1=>1),
        ℂ[ℤ₂](0=>1, 1=>2)',
        ℂ[ℤ₂](0=>3, 1=>2)',
        ℂ[ℤ₂](0=>2, 1=>3),
        ℂ[ℤ₂](0=>2, 1=>5))
Vℤ₃ = (ℂ[ℤ₃](0=>1, 1=>1, 2=>2),
        ℂ[ℤ₃](0=>3, 1=>2, 2=>1),
        ℂ[ℤ₃](0=>2, 1=>3, 2=>1),
        ℂ[ℤ₃](0=>1, 1=>2, 2=>3)',
        ℂ[ℤ₃](0=>2, 1=>4, 2=>3)')
VU₁ = (ℂ[U₁](0=>1, 1=>1, -1=>2),
        ℂ[U₁](0=>3, 1=>2, -1=>1),
        ℂ[U₁](0=>2, 1=>3, -1=>1),
        ℂ[U₁](0=>1, 1=>2, -1=>3)',
        ℂ[U₁](0=>2, 1=>4, -1=>3)')
VCU₁ = (ℂ[CU₁]((0,0)=>1, (0,1)=>1, 1=>2),
        ℂ[CU₁]((0,0)=>3, (0,1)=>2, 1=>1),
        ℂ[CU₁]((0,0)=>2, (0,1)=>3, 1=>1),
        ℂ[CU₁]((0,0)=>1, (0,1)=>2, 1=>3)',
        ℂ[CU₁]((0,0)=>2, (0,1)=>4, 1=>3)')
VSU₂ = (ℂ[SU₂](0=>1, 1//2=>1),
        ℂ[SU₂](0=>2, 1//2=>1, 1=>1),
        ℂ[SU₂](1//2=>2, 1=>1),
        ℂ[SU₂](0=>1, 1//2=>3)',
        ℂ[SU₂](0=>2, 3//2=>1)')

for (G,V) in ((Trivial, Vtr), (ℤ₂, Vℤ₂), (ℤ₃, Vℤ₃), (U₁, VU₁), (CU₁, VCU₁), (SU₂, VSU₂))
    println("------------------------------------")
    println("Tensors with symmetry: $G")
    println("------------------------------------")
    V1, V2, V3, V4, V5 = V
    @testset "Basic tensor properties" begin
        W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5
        for T in (Int, Float32, Float64, ComplexF32, ComplexF64, BigFloat)
            t = Tensor(zeros, T, W)
            @test eltype(t) == T
            @test norm(t) == 0
            @test codomain(t) == W
            @test space(t) == W
            @test domain(t) == one(W)
        end
    end
    @testset "Basic linear algebra" begin
        W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t = TensorMap(rand, T, W)
            @test eltype(t) == T
            @test codomain(t) == W.second
            @test domain(t) == W.first
            @test isa(@inferred(norm(t)), real(T))
            @test norm(t)^2 ≈ dot(t,t)
            α = rand(T)
            @test norm(α*t) ≈ abs(α)*norm(t)
            @test norm(t+t, 2) ≈ 2*norm(t, 2)
            @test norm(t+t, 1) ≈ 2*norm(t, 1)
            @test norm(t+t, Inf) ≈ 2*norm(t, Inf)
            p = 3*rand(Float64)
            @test norm(t+t, p) ≈ 2*norm(t, p)
            @test norm(t) ≈ norm(t')

            t2 = TensorMap(rand, T, W)
            β = rand(T)
            @test dot(β*t2,α*t) ≈ conj(β)*α*conj(dot(t,t2))
            @test dot(t2,t) ≈ conj(dot(t, t2))
            @test dot(t2,t) ≈ conj(dot(t2', t'))
            @test dot(t2,t) ≈ dot(t', t2')
        end
    end
    @testset "Basic linear algebra: test via conversion" begin
        W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t = TensorMap(rand, T, W)
            t2 = TensorMap(rand, T, W)
            @test norm(t, 2) ≈ norm(convert(Array,t), 2)
            @test dot(t2,t) ≈ dot(convert(Array,t2), convert(Array, t))
            α = rand(T)
            @test convert(Array, α*t) ≈ α*convert(Array,t)
            @test convert(Array, t+t) ≈ 2*convert(Array,t)
        end
    end
    @testset "Permutations: test via inner product invariance" begin
        W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5
        t = Tensor(rand, ComplexF64, W);
        t′ = Tensor(rand, ComplexF64, W);
        for k = 0:5
            for p in permutations(1:5)
                p1 = ntuple(n->p[n], StaticLength(k))
                p2 = ntuple(n->p[k+n], StaticLength(5-k))
                t2 = @inferred permuteind(t, p1, p2)
                @test norm(t2) ≈ norm(t)
                t2′= permuteind(t′, p1, p2)
                @test dot(t2′,t2) ≈ dot(t′,t)
            end
        end
    end
    @testset "Permutations: test via conversion" begin
        W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5
        t = Tensor(rand, ComplexF64, W);
        for k = 0:5
            for p in permutations(1:5)
                p1 = ntuple(n->p[n], StaticLength(k))
                p2 = ntuple(n->p[k+n], StaticLength(5-k))
                @test convert(Array, permuteind(t, p1, p2)) ≈ permutedims(convert(Array, t), (p1...,p2...))
            end
        end
    end
    @testset "Full trace: test self-consistency" begin
        t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1');
        t2 = permuteind(t, (1,2), (4,3))
        s = @inferred tr(t2)
        @test conj(s) ≈ tr(t2')
        @tensor s2 = t[a,b,b,a]
        @tensor t3[a,b] := t[a,c,c,b]
        @tensor s3 = t3[a,a]
        @test s ≈ s2
        @test s ≈ s3
    end
    @testset "Partial trace: test self-consistency" begin
        t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V1);
        @tensor t2[a,b] := t[c,d,b,d,c,a]
        @tensor t4[a,b,c,d] := t[d,e,b,e,c,a]
        @tensor t5[a,b] := t4[a,b,c,c]
        @test t2 ≈ t5
    end
    @testset "Trace: test via conversion" begin
        t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V1);
        @tensor t2[a,b] := t[c,d,b,d,c,a]
        @tensor t3[a,b] := convert(Array, t)[c,d,b,d,c,a]
        @test t3 ≈ convert(Array, t2)
    end
    @testset "Trace and contraction" begin
        t1 = Tensor(rand, ComplexF64, V1 ⊗ V2 ⊗ V3);
        t2 = Tensor(rand, ComplexF64, V2' ⊗ V4 ⊗ V1');
        t3 = t1 ⊗ t2
        @tensor ta[a,b] := t1[x,y,a]*t2[y,b,x]
        @tensor tb[a,b] := t3[x,y,a,y,b,x]
        @test ta ≈ tb
    end
    @testset "Tensor contraction: test via conversion" begin
        A1 = TensorMap(randn, ComplexF64, V1'*V2', V3')
        A2 = TensorMap(randn, ComplexF64, V3*V4, V5)
        rhoL = TensorMap(randn, ComplexF64, V1, V1)
        rhoR = TensorMap(randn, ComplexF64, V5, V5)' # test adjoint tensor
        H = TensorMap(randn, ComplexF64, V2*V4, V2*V4)
        @tensor HrA12[a, s1, s2, c] := rhoL[a, a'] * conj(A1[a', t1, b]) *
            A2[b, t2, c'] * rhoR[c', c] * H[s1, s2, t1, t2]

        @tensor HrA12array[a, s1, s2, c] := convert(Array, rhoL)[a, a'] *
            conj(convert(Array, A1)[a', t1, b]) *
            convert(Array, A2)[b, t2, c'] *
            convert(Array, rhoR)[c', c] *
            convert(Array, H)[s1, s2, t1, t2]

        @test HrA12array ≈ convert(Array, HrA12)
    end
    @testset "Multiplication and inverse: test compatibility" begin
        W1 = V1 ⊗ V2 ⊗ V3
        W2 = V4 ⊗ V5
        for T in (Float64, ComplexF64)
            t1 = TensorMap(rand, T, W1, W1)
            t2 = TensorMap(rand, T, W2, W2)
            t = TensorMap(rand, T, W1, W2)
            @test t1*(t1\t) ≈ t
            @test (t/t2)*t2 ≈ t
            @test t1\one(t1) ≈ inv(t1)
            @test one(t1)/t1 ≈ pinv(t1)
            @test_throws SpaceMismatch inv(t)
            @test_throws SpaceMismatch t2\t
            @test_throws SpaceMismatch t/t1
            tp = pinv(t)*t
            @test tp ≈ tp*tp
        end
    end
    @testset "Multiplication and inverse: test via conversion" begin
        W1 = V1 ⊗ V2 ⊗ V3
        W2 = V4 ⊗ V5
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t1 = TensorMap(rand, T, W1, W1)
            t2 = TensorMap(rand, T, W2, W2)
            t = TensorMap(rand, T, W1, W2)
            d1 = dim(W1)
            d2 = dim(W2)
            At1 = reshape(convert(Array, t1), d1, d1)
            At2 = reshape(convert(Array, t2), d2, d2)
            At = reshape(convert(Array, t), d1, d2)
            @test reshape(convert(Array, t1*t), d1, d2) ≈ At1*At
            @test reshape(convert(Array, t1'*t), d1, d2) ≈ At1'*At
            @test reshape(convert(Array, t2*t'), d2, d1) ≈ At2*At'
            @test reshape(convert(Array, t2'*t'), d2, d1) ≈ At2'*At'

            @test reshape(convert(Array, inv(t1)), d1, d1) ≈ inv(At1)
            @test reshape(convert(Array, pinv(t)), d2, d1) ≈ pinv(At)

            if T == Float32 || T == ComplexF32
                continue
            end

            @test reshape(convert(Array, t1\t), d1, d2) ≈ At1\At
            @test reshape(convert(Array, t1'\t), d1, d2) ≈ At1'\At
            @test reshape(convert(Array, t2\t'), d2, d1) ≈ At2\At'
            @test reshape(convert(Array, t2'\t'), d2, d1) ≈ At2'\At'

            @test reshape(convert(Array, t2/t), d2, d1) ≈ At2/At
            @test reshape(convert(Array, t2'/t), d2, d1) ≈ At2'/At
            @test reshape(convert(Array, t1/t'), d1, d2) ≈ At1/At'
            @test reshape(convert(Array, t1'/t'), d1, d2) ≈ At1'/At'
        end
    end
    @testset "Factorization" begin
        W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t = Tensor(rand, T, W)
            @testset "leftorth with $alg" for alg in (TensorKit.QR(), TensorKit.QRpos(), TensorKit.QL(), TensorKit.QLpos(), TensorKit.Polar(), TensorKit.SVD(), TensorKit.SDD())
                Q, R = @inferred leftorth(t, (3,4,2),(1,5); alg = alg)
                QdQ = Q'*Q
                @test QdQ ≈ one(QdQ)
                @test Q*R ≈ permuteind(t, (3,4,2),(1,5))
            end
            @testset "leftnull with $alg" for alg in (TensorKit.QR(), TensorKit.QRpos())
                N = @inferred leftnull(t, (3,4,2),(1,5); alg = alg)
                NdN = N'*N
                @test NdN ≈ one(NdN)
                @test norm(N'*permuteind(t, (3,4,2),(1,5))) < 100*eps(norm(t))
            end
            @testset "rightorth with $alg" for alg in (TensorKit.RQ(), TensorKit.RQpos(), TensorKit.LQ(), TensorKit.LQpos(), TensorKit.Polar(), TensorKit.SVD(), TensorKit.SDD())
                L, Q = @inferred rightorth(t, (3,4),(2,1,5); alg = alg)
                QQd = Q*Q'
                @test QQd ≈ one(QQd)
                @test L*Q ≈ permuteind(t, (3,4),(2,1,5))
            end
            @testset "rightnull with $alg" for alg in (TensorKit.LQ(), TensorKit.LQpos())
                M = @inferred rightnull(t, (3,4),(2,1,5); alg = alg)
                MMd = M*M'
                @test MMd ≈ one(MMd)
                @test norm(permuteind(t, (3,4),(2,1,5))*M') < 100*eps(norm(t))
            end
            @testset "svd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD())
                U, S, V = @inferred svd(t, (3,4,2),(1,5); alg = alg)
                UdU = U'*U
                @test UdU ≈ one(UdU)
                VVd = V*V'
                @test VVd ≈ one(VVd)
                @test U*S*V ≈ permuteind(t, (3,4,2),(1,5))
            end

            t = Tensor(rand, T, V1 ⊗ V1' ⊗ V2 ⊗ V2')
            @testset "eig and isposdef" begin
                D, V = eigen(t, (1,3), (2,4))
                VdV = V'*V
                VdV = (VdV + VdV')/2
                @test isposdef(VdV)
                t2 = permuteind(t, (1,3), (2,4))
                @test t2*V ≈ V*D
                @test !isposdef(t2) # unlikely for non-hermitian map
                t2 = (t2 + t2');
                D, V = eigen(t2)
                VdV = V'*V
                @test VdV ≈ one(VdV)
                λ = minimum(minimum(real(diag(b))) for (c,b) in blocks(D))
                @test isposdef(t2) == isposdef(λ)
                @test isposdef(t2 - λ*one(t2) + 0.1*one(t2))
                @test !isposdef(t2 - λ*one(t2) - 0.1*one(t2))
            end
        end
    end
    @testset "Tensor truncation" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            for p in (1, 2, 3, Inf)
                t = TensorMap(randn, T, V1 ⊗ V2 ⊗ V3 , V4 ⊗ V5)
                U₀, S₀, V₀, = svd(t)
                t = rmul!(t, 1/norm(S₀, p))
                U, S, V, ϵ = @inferred svd(t; trunc = truncerr(5e-1), p = p)
                # @show p, ϵ
                # @show domain(S)
                # @test min(space(S,1), space(S₀,1)) != space(S₀,1)
                U′, S′, V′, ϵ′ = svd(t; trunc = truncerr(nextfloat(ϵ)), p = p)
                @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
                U′, S′, V′, ϵ′ = svd(t; trunc = truncdim(dim(domain(S))), p = p)
                @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
                U′, S′, V′, ϵ′ = svd(t; trunc = truncspace(space(S,1)), p = p)
                @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
                # results with truncationcutoff cannot be compared because they don't take degeneracy into account, and thus truncate differently
                U, S, V, ϵ = svd(t; trunc = truncbelow(1/dim(domain(S₀))), p = p)
                # @show p, ϵ
                # @show domain(S)
                # @test min(space(S,1), space(S₀,1)) != space(S₀,1)
                U′, S′, V′, ϵ′ = svd(t; trunc = truncspace(space(S,1)), p = p)
                @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
            end
        end
    end
    @testset "Exponentiation" begin
        W = V1 ⊗ V2 ⊗ V3
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t = TensorMap(rand, T, W, W)
            s = dim(W)
            expt = @inferred exp(t)
            @test reshape(convert(Array, expt), (s,s)) ≈
                    exp(reshape(convert(Array, t), (s,s)))
        end
    end
    @testset "Tensor product: test via norm preservation" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t1 = TensorMap(rand, T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2)
            t2 = TensorMap(rand, T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1)
            t = @inferred (t1 ⊗ t2)
            @test norm(t) ≈ norm(t1) * norm(t2)
        end
    end
    @testset "Tensor product: test via conversion" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t1 = TensorMap(rand, T, V2 ⊗ V3 ⊗ V1, V1)
            t2 = TensorMap(rand, T, V2 ⊗ V1 ⊗ V3, V2)
            t = @inferred (t1 ⊗ t2)
            d1 = dim(codomain(t1))
            d2 = dim(codomain(t2))
            d3 = dim(domain(t1))
            d4 = dim(domain(t2))
            At = convert(Array, t)
            @show sizeof(At)
            @test reshape(At, (d1, d2, d3, d4)) ≈
                    reshape(convert(Array, t1), (d1, 1, d3, 1)) .*
                    reshape(convert(Array, t2), (1, d2, 1, d4))
        end
    end
    @testset "Tensor product: test via tensor contraction" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            t1 = Tensor(rand, T, V2 ⊗ V3 ⊗ V1)
            t2 = Tensor(rand, T, V2 ⊗ V1 ⊗ V3)
            t = @inferred (t1 ⊗ t2)
            @tensor t′[1, 2, 3, 4, 5, 6] := t1[1,2,3]*t2[4,5,6]
            @test t ≈ t′
        end
    end
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

# This implements the functions needed by the waiting time method
# for the choice of the next move

module WaitingTimes

using ExtractMacro
using DataStructures
using ..Common
using ..Interface

export THeap, pick_next, update_heap!

τΔE(X::AbstractGraph, C::Config, β::Real, ΔE) = max(1.0, Float64(exp(β * ΔE)))
τ(X::AbstractGraph, C::Config, β::Real, i::Integer) = τΔE(X, C, β, delta_energy(X, C, i))
gen_wt(τ::Float64) = begin
    wt = -τ * log1p(-rand())
    #@assert wt ≥ 0 wt, τ, log1p(-rand())
    return wt
end

const THeap = MutableBinaryMinHeap{Float64}

function THeap(X::AbstractGraph, C::Config, β::Real)
    N = getN(X)
    @assert C.N == N
    τs = [τ(X, C, β, i) for i = 1:N]
    theap = THeap()
    for i = 1:N
        j = push!(theap, gen_wt(τs[i]))
        @assert j == i # implementation detail check
    end
    return theap
end

pick_next(theap::THeap) = top_with_handle(theap)

function update_heap!(theap::THeap, X::AbstractGraph, C::Config, β::Real, i::Integer, t::Float64)
    ΔE = delta_energy(X, C, i)
    spinflip!(X, C, i)
    τi = τΔE(X, C, β, -ΔE)
    ti = t + gen_wt(τi)
    update!(theap, i, ti)
    for j in neighbors(X, i)
        τj = τ(X, C, β, j)
        tj = t + gen_wt(τj)
        update!(theap, j, tj)
    end
    return ΔE
end

end # module


text:
push!(LOAD_PATH, ".")

module Debug

const ga_debug = false
const ga_plot = false

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
reversedims(A::Array) = permutedims(A, ndims(A):-1:1) # @generated probably the "right" way

for MAType in (Float32MultiArray, Float64MultiArray, Int32MultiArray, Int64MultiArray)
    @eval function Base.convert(::Type{$MAType}, A::Array{<:Real})
        Arev = reversedims(A)
        msg = $MAType()
        msg.layout.dim = [MultiArrayDimension("",d,s) for (d,s) in zip(size(A), reverse(cumprod([size(Arev)...])))]
        msg.data = Arev[:]
        msg
    end

    @eval function Base.convert(::Type{Array{T}}, msg::$MAType) where {T}
        if isempty(msg.data)
            T[]
        else
            reversedims(reshape(msg.data, (Int(d.size) for d in msg.layout.dim[end:-1:1])...))
        end
    end
end

mutable struct PredictionOutputHolder
    stale::Bool
    y
    z
    r
    b
end

function prediction_callback(msg::PredictionOutput, pred_container::PredictionOutputHolder)
    pred_container.stale = false
    pred_container.y = Array{Float32}(msg.y)
    pred_container.z = Array{Float32}(msg.z)
    pred_container.r = Array{Float32}(msg.r)
    pred_container.b = Array{Float32}(msg.b)
end

const pred_container = PredictionOutputHolder(true,nothing,nothing,nothing,nothing)
const publisher = Publisher{PredictionInput}("prediction_input", queue_size=10)
const subscriber = Subscriber{PredictionOutput}("prediction_output", prediction_callback, (pred_container,), queue_size=10)

function predict!(car1, car2, extras, traj_lengths, car1_future, sample_ct,
                  output::PredictionOutputHolder=pred_container, timeout=1000)
    output.stale = true
    publish(publisher, PredictionInput(Float32MultiArray(car1),
                                       Float32MultiArray(car2),
                                       Float32MultiArray(extras),
                                       Int32MultiArray(traj_lengths),
                                       Float32MultiArray(car1_future),
                                       Float32MultiArray(Float32[]),
                                       Float32MultiArray(Float32[]),
                                       Int32(sample_ct)))
    for i in 1:timeout
        rossleep(Duration(0.001))
        output.stale || break
    end
    output
end

function predict!(car1, car2, extras, traj_lengths, car1_future_x, car1_future_y, sample_ct,
                  output::PredictionOutputHolder=pred_container, timeout=1000)
    output.stale = true
    publish(publisher, PredictionInput(Float32MultiArray(car1),
                                       Float32MultiArray(car2),
                                       Float32MultiArray(extras),
                                       Int32MultiArray(traj_lengths),
                                       Float32MultiArray(Float32[]),
                                       Float32MultiArray(car1_future_x),
                                       Float32MultiArray(car1_future_y),
                                       Int32(sample_ct)))
    for i in 1:timeout
        rossleep(Duration(0.001))
        output.stale || break
    end
    output
end

######
# const y_top = -1.83284f0
# const y_bot = -6.09485f0
const y_top = -1.5f0
const y_bot = -6.5f0
const xdd_choices = Float32[0, -3, 4, -6]
const y_targets = [y_top, y_bot]
const tibvp = SteeringBVP(TripleIntegratorDynamics(1, Float32), TimePlusQuadraticControl(SMatrix{1,1}(Float32(1e-3))))

trapz(x::Vector, dt, x0=0) = cumsum(vcat(0, 0.5*dt*(x[1:end-1] + x[2:end]))) + x0
trapz(x::Matrix, dt, x0=0) = cumsum(hcat(zeros(size(x,1)), 0.5*dt*(x[:,1:end-1] + x[:,2:end])), 2) .+ x0
vel(xd, yd) = sqrt.(xd.^2 + yd.^2)
acc(xd, xdd, yd, ydd, v=vel(xd,yd)) = (xd.*xdd + yd.*ydd)./v
curv(xd, xdd, yd, ydd, v=vel(xd,yd)) = (xd.*ydd - xdd.*yd)./v.^3

@inbounds @views function all_future_xy(xdd_choices, y_targets, state0, tibvp,
                                        a0 = [(1,1)], ply=4, t_hold=3, dt=Float32(0.1);
                                        des_length=length(a0)+ply*t_hold)
    Nx = length(xdd_choices)
    Ny = length(y_targets)
    X = zeros(Float32, Nx^ply, 1 + des_length, 3)
    Y = zeros(Float32, Ny^ply, 1 + des_length, 3)
    if a0 isa Vector{Tuple{Int,Int}}
        X[:,1,:] .= state0[1:2:end]'
        Y[:,1,:] .= state0[2:2:end]'
        t0 = 1 + length(a0)
        for (t, (xi, yi)) in enumerate(a0)
            X[:,t+1,3] .= xdd_choices[xi]
            q0 = SVector(Y[1,t,1],
                         Y[1,t,2],
                         Y[1,t,3])
            qf = SVector(y_targets[yi], 0f0, 0f0)
            c, u = tibvp(q0, qf, 100f0)
            if dt > duration(u)
                q = qf
            else
                q = u.x(q0, u.y, u.t, dt)
            end
            Y[:,t+1,1] .= q[1]
            Y[:,t+1,2] .= q[2]
            Y[:,t+1,3] .= q[3]
        end
    elseif a0 isa Tuple
        Ẍ0, Ÿ0 = a0
        X[:,1,:] .= state0[1:2:end]'
        Y[:,1,:] .= state0[2:2:end]'
        X[:,2:1+length(Ẍ0),3] .= Ẍ0'
        Y[:,2:1+length(Ÿ0),3] .= Ÿ0'
        t0 = 1 + length(Ÿ0)
        cumsum!(Y[:,2:t0,2], Y[:,1:t0-1,3]*dt, 2)    # X integrated below (joint with the case above)
        Y[:,2:t0,2] .+= Y[:,1,2]
        cumsum!(Y[:,2:t0,1], (Y[:,1:t0-1,2] .+ Y[:,2:t0,2])*dt./2, 2)
        Y[:,2:t0,1] .+= Y[:,1,1]
    else
        error("Incompatible a0 type.")
    end

    # xdd, xd, x
    for (r, ci) in enumerate(CartesianRange(Tuple(Nx for i in 1:ply)))
        for (p, v) in enumerate(reverse(ci.I))    # reverse so that later indices (actions for later timesteps) iterate fastest
            for t in 1:t_hold
                X[r,t0+t_hold*(p-1)+t,3] = xdd_choices[v]
            end
        end
    end
    for t in (t0+ply*t_hold+1):1+des_length
        X[:,t,3] = X[:,t-1,3]
    end
    cumsum!(X[:,2:end,2], X[:,1:end-1,3]*dt, 2)
    X[:,2:end,2] .+= X[:,1,2]
    cumsum!(X[:,2:end,1], (X[:,1:end-1,2] .+ X[:,2:end,2])*dt./2, 2)
    X[:,2:end,1] .+= X[:,1,1]

    # ydd, yd, y
    for (r, ci) in enumerate(CartesianRange(Tuple(Ny for i in 1:ply)))
        for (p, v) in enumerate(reverse(ci.I))    # reverse so that later indices (actions for later timesteps) iterate fastest
            curr_step = t0 + (p-1)*t_hold
            q0 = SVector(Y[r,curr_step,1],
                         Y[r,curr_step,2],
                         Y[r,curr_step,3])
            qf = SVector(y_targets[v], 0f0, 0f0)
            c, u = tibvp(q0, qf, 100f0)
            for t in 1:t_hold
                if dt*t > duration(u)
                    q = qf
                else
                    q = u.x(q0, u.y, u.t, t*dt)
                end
                Y[r,curr_step+t,1] = q[1]
                Y[r,curr_step+t,2] = q[2]
                Y[r,curr_step+t,3] = q[3]
            end
        end
    end
    for t in (t0+ply*t_hold+1):1+des_length
        Y[:,t,3] = Y[:,t-1,3]
        Y[:,t,2] = Y[:,t-1,2] + Y[:,t-1,3]*dt
        Y[:,t,1] = Y[:,t-1,1] + Y[:,t-1,2]*dt
    end
    X, Y
end

function propagate(xdd_choices, y_targets, car1_curr_state, tibvp, action_sequence)
    pfX, pfY = all_future_xy(xdd_choices, y_targets, car1_curr_state, tibvp, action_sequence, 0)
    [SVector{6,Float32}(pfX[1,i,1], pfY[1,i,1],
                        pfX[1,i,2], pfY[1,i,2],
                        pfX[1,i,3], pfY[1,i,3]) for i in 1:size(pfX, 2)]
end


text:
"""
Functor to update JS part when an observable changes.
We make this a Functor, so we can clearly identify it and don't sent
any updates, if the JS side requires to update an Observable
(so we don't get an endless update cycle)
"""
struct JSUpdateObservable
    session::Session
    id::String
end

function (x::JSUpdateObservable)(value)
    # Sent an update event
    send(x.session, payload=value, id=x.id, msg_type=UpdateObservable)
end

"""
Update the value of an observable, without sending changes to the JS frontend.
This will be used to update updates from the forntend.
"""
function update_nocycle!(obs::Observable, @nospecialize(value))
    obs.val = value
    for f in Observables.listeners(obs)
        if !(f isa JSUpdateObservable)
            Base.invokelatest(f, value)
        end
    end
    return
end

function jsrender(session::Session, obs::Observable)
    html = Observable{Any}(
        repr_richest(jsrender(session, obs[]))
    )
    dom = DOM.m_unesc("span", html[])
    on(session, obs) do data
        html[] = by_value(render_sub_session(session, data))
    end
    onjs(session, html, js"(html)=> {
        const new_node = JSServe.materialize_node(html)
        const dom = $(dom)
        dom.replaceChild(new_node, dom.childNodes[0])
    }")
    return dom
end

# on & map versions that deregister when session closes!
function Observables.on(f, session::Session, observable::Observable)
    to_deregister = on(f, observable)
    push!(session.deregister_callbacks, to_deregister)
    return to_deregister
end

function Observables.onany(f, session::Session, observables::Observable...)
    to_deregister = onany(f, observables...)
    append!(session.deregister_callbacks, to_deregister)
    return to_deregister
end

function Base.map(f, session::Session, observables::Observable...; result=Observable{Any}())
    # map guarantees to be run upfront!
    result[] = f(Observables.to_value.(observables)...)
    onany(session, observables...) do newvals...
        result[] = f(newvals...)
    end
    return result
end


text:
const M=100_000_000; const K=100;
srand(1);
svec1 = rand([string(rand(Char.(32:126), 16)...) for k in 1:M÷K], M);

using FastGroupBy
@time x= FastGroupBy.load_bits.(UInt128, svec1);
@time grouptwo!(x, svec1)
@time x = ntoh.(x);
@time sorttwo!(x, svec1);
issorted(svec1)


text:
InputFolder = './Images/Kidney/';
OutputFolder = './Results/Images/Kidney/';
Fill = 1;
Dilate = 5;
Lbl = 1;

@iA = '*_C0000*.tif';		% Image filter
@iS = '*_C0001*.tif';		% Image annotations filter

@fxg_lBlockClassify[iA, iS] > [S];
params.BlckSize = 32;
params.Feat = 'mnstdbifs';
params.ClassifierType = 'RF';
params.ClassifierFile = './Classifiers/ClassifierKidney_mnstdbifsRF.mat';
params.ExportAnnotations = './Classifiers/Annotation_Kidney.tif';
/endf

/show iA > S;
/keep S > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
module PandorasBoxTileGlitcher

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/tileGlitcher" TileGlitcher(x::Integer, y::Integer, allowAir::Bool=false, transformAir::Bool=false, target::String="Both", threshold::Number=0.1, rate::Number=0.05, flag::String="", customFgTiles::String="", customBgTiles::String="")

const placements = Ahorn.PlacementDict(
    "Tile Glitcher (Pandora's Box)" => Ahorn.EntityPlacement(
        TileGlitcher,
        "rectangle"
    ),
)

Ahorn.editingOptions(entity::TileGlitcher) = Dict{String, Any}(
    "target" => String[
        "None",
        "FG",
        "BG",
        "Both"
    ]
)

Ahorn.minimumSize(entity::TileGlitcher) = 8, 8
Ahorn.resizable(entity::TileGlitcher) = true, true

Ahorn.selection(entity::TileGlitcher) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::TileGlitcher, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.5, 0.3, 1.0, 0.4), (0.5, 0.3, 1.0, 1.0))
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
using ITensors, Test, Random

@testset "Basic DMRG" begin

  @testset "Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1",N)

    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    H = MPO(ampo,sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy,psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "QN-conserving Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1",N; conserve_qns=true)

    ampo = AutoMPO()
    for j=1:N-1
      ampo += "Sz",j,"Sz",j+1
      ampo += 0.5,"S+",j,"S-",j+1
      ampo += 0.5,"S-",j,"S+",j+1
    end
    H = MPO(ampo,sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites,state,4)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy,psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "Transverse field Ising" begin
    N = 32
    sites = siteinds("S=1/2",N)
    Random.seed!(432)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"Z",j,"Z",j+1)
      add!(ampo,-1.0,"X",j)
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    noise!(sweeps,1E-10)
    energy,psi = dmrg(H,psi0,sweeps,outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 1.0 - 1.0/sin(π/(4*N+2))
    @test abs((energy-energy_exact)/energy_exact) < 1e-4
  end

  @testset "Transverse field Ising, conserve Sz parity" begin
    N = 32
    sites = siteinds("S=1/2", N, conserve_szparity = true)
    Random.seed!(432)

    state = [isodd(j) ? "↑" : "↓" for j in 1:N]
    psi0 = randomMPS(sites, state)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"X",j,"X",j+1)
      add!(ampo,-1.0,"Z",j)
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    noise!(sweeps,1E-10)
    energy,psi = dmrg(H,psi0,sweeps,outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 1.0 - 1.0/sin(π/(4*N+2))
    @test abs((energy-energy_exact)/energy_exact) < 1e-4
  end

  @testset "DMRGObserver" begin

    # Test that basic constructors work
    observer = DMRGObserver()
    observer = DMRGObserver(;minsweeps=2,energy_tol=1E-4)

    # Test in a DMRG calculation
    N = 10
    sites = siteinds("S=1/2",N)
    Random.seed!(42)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N-1
      ampo += -1,"Sz",j,"Sz",j+1
    end
    for j = 1:N
      ampo += -0.2,"Sx",j
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10)
    cutoff!(sweeps,1E-12)

    observer = DMRGObserver(["Sz","Sx"], sites)

    E,psi = dmrg(H,psi0,sweeps,observer=observer,outputlevel=0)
    @test length(measurements(observer)["Sz"])==3
    @test length(measurements(observer)["Sx"])==3
    @test all(length.(measurements(observer)["Sz"]) .== N)
    @test all(length.(measurements(observer)["Sx"]) .== N)
    @test length(energies(observer))==3
    @test length(truncerrors(observer))==3
    @test energies(observer)[end]==E
    @test all(truncerrors(observer) .< 1E-9)

    orthogonalize!(psi,1)
    m = scalar(dag(psi[1])*noprime(op(sites, "Sz", 1)*psi[1]))
    @test measurements(observer)["Sz"][end][1] ≈ m
  end

  @testset "Sum of MPOs (ProjMPOSum)" begin
    N = 10
    sites = siteinds("S=1",N)

    ampoZ = AutoMPO()
    for j=1:N-1
      ampoZ += "Sz",j,"Sz",j+1
    end
    HZ = MPO(ampoZ,sites)

    ampoXY = AutoMPO()
    for j=1:N-1
      ampoXY += 0.5,"S+",j,"S-",j+1
      ampoXY += 0.5,"S-",j,"S+",j+1
    end
    HXY = MPO(ampoXY,sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    energy,psi = dmrg([HZ,HXY], psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "Excited-state DMRG" begin
    N = 10
    weight = 15.0

    sites = siteinds("S=1",N)
    sites[1] = Index(2,"S=1/2,n=1,Site")
    sites[N] = Index(2,"S=1/2,n=$N,Site")

    ampo = AutoMPO()
    for j=1:N-1
      ampo += "Sz",j,"Sz",j+1
      ampo += 0.5,"S+",j,"S-",j+1
      ampo += 0.5,"S-",j,"S+",j+1
    end
    H = MPO(ampo,sites)

    psi0i = randomMPS(sites,10)

    sweeps = Sweeps(4)
    maxdim!(sweeps, 10,20,100,100)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps,1E-10)

    energy0, psi0 = dmrg(H,psi0i, sweeps; outputlevel=0)
    @test energy0 < -11.5

    psi1i = randomMPS(sites,10)
    energy1,psi1 = dmrg(H,[psi0],psi1i,sweeps;outputlevel=0,weight=weight)

    @test energy1 > energy0
    @test energy1 < -11.1

    @test inner(psi1,psi0) < 1E-5
  end

  @testset "Fermionic Hamiltonian" begin
    N = 10
    t1 = 1.0
    t2 = 0.5
    V = 0.2
    s = siteinds("Fermion", N; conserve_qns = true)

    state = fill(1,N)
    state[1] = 2
    state[3] = 2
    state[5] = 2
    state[7] = 2
    psi0 = productMPS(s,state)

    ampo = AutoMPO()
    for j=1:N-1
      ampo += -t1,"Cdag",j,"C",j+1
      ampo += -t1,"Cdag",j+1,"C",j
      ampo += V,"N",j,"N",j+1
    end
    for j=1:N-2
      ampo += -t2,"Cdag",j,"C",j+2
      ampo += -t2,"Cdag",j+2,"C",j
    end
    H = MPO(ampo, s)

    sweeps = Sweeps(5)
    maxdim!(sweeps, 10, 20, 100, 100, 200)
    cutoff!(sweeps, 1E-8)
    noise!(sweeps, 1E-10)

    energy, psi = dmrg(H, psi0, sweeps; outputlevel = 0)
    @test (-6.5 < energy < -6.4)
  end

  @testset "Hubbard model" begin
    N = 10
    Npart = 8
    t1 = 1.0
    U  = 1.0
    V1 = 0.5
    sites = siteinds("Electron",N; conserve_qns=true)
    ampo = AutoMPO()
    for i=1:N
      ampo += (U,"Nupdn",i)
    end
    for b=1:N-1
      ampo += -t1,"Cdagup",b,"Cup",b+1
      ampo += -t1,"Cdagup",b+1,"Cup",b
      ampo += -t1,"Cdagdn",b,"Cdn",b+1
      ampo += -t1,"Cdagdn",b+1,"Cdn",b
      ampo += V1,"Ntot",b,"Ntot",b+1
    end
    H = MPO(ampo,sites)
    sweeps = Sweeps(6)
    maxdim!(sweeps,50,100,200,400,800,800)
    cutoff!(sweeps,1E-10)
    state = ["Up",
             "Dn",
             "Dn",
             "Up",
             "Emp",
             "Up",
             "Up",
             "Emp",
             "Dn",
             "Dn"
            ]
    psi0 = randomMPS(sites,state,10)
    energy,psi = dmrg(H,psi0,sweeps; outputlevel = 0)
    @test (-8.02 < energy < -8.01)
  end

  @testset "Input Without Ortho Center or Not at 1" begin
    N = 6
    sites = siteinds("S=1",N)

    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(1)
    maxdim!(sweeps,10)
    cutoff!(sweeps,1E-11)

    psi0 = randomMPS(sites,4)

    # Test that input works with wrong ortho center:
    orthogonalize!(psi0,5)
    energy,psi = dmrg(H, psi0, sweeps; outputlevel=0)

    # Test that input works with no ortho center:
    for j=1:N
      psi0[j] = randomITensor(inds(psi0[j]))
    end
    energy,psi = dmrg(H, psi0, sweeps; outputlevel=0)
  end

  @testset "Bug fixed in threaded block sparse" begin
    maxdim = 10
    nsweeps = 2
    outputlevel = 0
    cutoff = 0.0
    Nx = 4
    Ny = 2
    U = 4.0
    t = 1.0
    N = Nx * Ny
    sweeps = Sweeps(nsweeps)
    maxdims = min.(maxdim, [100, 200, 400, 800, 2000, 3000, maxdim])
    maxdim!(sweeps, maxdims...)
    cutoff!(sweeps, cutoff)
    noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
    sites = siteinds("Electron", N; conserve_qns = true)
    lattice = square_lattice(Nx, Ny; yperiodic = true)
    ampo = AutoMPO()
    for b in lattice
      ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
      ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
      ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
      ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
    end
    for n in 1:N
      ampo .+= U, "Nupdn", n
    end
    H = MPO(ampo,sites)
    Hsplit = splitblocks(linkinds, H)
    state = [isodd(n) ? "↑" : "↓" for n in 1:N]
    ψ0 = productMPS(sites, state)
    using_threaded_blocksparse = ITensors.enable_threaded_blocksparse()
    energy, _ = dmrg(H, ψ0, sweeps; outputlevel = outputlevel)
    energy_split, _ = dmrg(Hsplit, ψ0, sweeps; outputlevel = outputlevel)
    @test energy_split ≈ energy
    if !using_threaded_blocksparse
      ITensors.disable_threaded_blocksparse()
    end
  end
end

nothing


text:
using Test

# Run test suites
include("testdecimalplaces.jl")
include("testsigdigits.jl")


text:
import GraphIdx.Grid: incmat

incmat(2, 2, 2)
@time a = incmat(1000, 1000, 2)

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx