text:
using DelimitedFiles

function move(pos, cmd, dir)
	w = cmd[1]
	n = parse(Int, cmd[2:end])
	if w == 'N'
		pos += [n,0]
	elseif w == 'S'
		pos -= [n,0]
	elseif w == 'E'
		pos += [0,n]
	elseif w == 'W'
		pos -= [0,n]
	elseif w == 'L'
		dir = [cosd(n) sind(n);-sind(n) cosd(n)]*dir
	elseif w == 'R'
		dir = [cosd(n) -sind(n);sind(n) cosd(n)]*dir
	elseif w == 'F'
		pos += n*dir
	end
	pos = round.(pos)
	dir = round.(dir)
	(Int.(pos),Int.(dir))
end

# returns solution of part 1
function x()
	t = readdlm("12.txt", String)
	pos = [0,0]
	dir = [0,1]
	for c in t
		(pos, dir) = move(pos, c, dir)
	end
	abs(pos[1])+abs(pos[2])
end

function move2(spos, wpos, cmd)
	w = cmd[1]
	n = parse(Int, cmd[2:end])
	if w == 'N'
		wpos += [n,0]
	elseif w == 'S'
		wpos -= [n,0]
	elseif w == 'E'
		wpos += [0,n]
	elseif w == 'W'
		wpos -= [0,n]
	elseif w == 'L'
		diff = wpos - spos
		diff = [cosd(n) sind(n);-sind(n) cosd(n)]*diff
		wpos = spos + diff
	elseif w == 'R'
		diff = wpos - spos
		diff = [cosd(n) -sind(n);sind(n) cosd(n)]*diff
		wpos = spos + diff
	elseif w == 'F'
		mv = wpos - spos
		spos += n*mv
		wpos += n*mv
	end
	spos = round.(spos)
	wpos = round.(wpos)
	(Int.(spos),Int.(wpos))
end

# returns solution of part 2
function y()
	t = readdlm("12.txt", String)
	spos = [0,0]
	wpos = [1,10]
	for c in t
		(spos, wpos) = move2(spos, wpos, c)
	end
	abs(spos[1])+abs(spos[2])
end


text:
using Printf

randN(N) = () -> rand(1:N) == 1 ? 1 : 0
function unbiased(biased::Function)
    this, that = biased(), biased()
    while this == that this, that = biased(), biased() end
    return this
end

@printf "%2s | %10s | %5s | %5s | %8s" "N" "bias./unb." "1s" "0s" "pct ratio"
const nrep = 10000
for N in 3:6
    biased = randN(N)

    v = collect(biased() for __ in 1:nrep)
    v1, v0 = count(v .== 1), count(v .== 0)
    @printf("%2i | %10s | %5i | %5i | %5.2f%%\n", N, "biased", v1, v0, 100 * v1 / nrep)

    v = collect(unbiased(biased) for __ in 1:nrep)
    v1, v0 = count(v .== 1), count(v .== 0)
    @printf("%2i | %10s | %5i | %5i | %5.2f%%\n", N, "unbiased", v1, v0, 100 * v1 / nrep)
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Simulate_MassWithDamper

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

# Damper
@forceElement PDamper(; d=1.0) begin
    v = ModiaMath.RealScalar("v", causality=ModiaMath.Input,  numericType=ModiaMath.WR)
    f = ModiaMath.RealScalar("f", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeTorque(damper::PDamper, sim::ModiaMath.SimulationState)
    damper.f.value = -damper.d*damper.v.value
end

@assembly MassWithDamper(;Lx = 0.1, Ly=0.1, Lz=0.1, m=1.0) begin
   world  = Object3D(visualizeFrame=true)

   body   = Object3D(Solid(SolidBox(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   prismatic = Prismatic(world, body, axis=2, v_start=0.1)

   damper         = PDamper(d=0.2)
   damper_adaptor = Modia3D.AdaptorForceElementToPFlange(v=damper.v, f=damper.f)
   Modia3D.connect(damper_adaptor, prismatic)
end

gravField      = UniformGravityField(g=9.81, n=[0,0,-1])
massWithDamper = MassWithDamper(sceneOptions=SceneOptions(gravityField=gravField, visualizeFrames=true, defaultFrameLength=0.3))
model  = Modia3D.SimulationModel( massWithDamper )
result = ModiaMath.simulate!(model, stopTime=5.0, tolerance=1e-4, log=true)

ModiaMath.plot(result, ["prismatic.s", "prismatic.v", "prismatic.a", "prismatic.f"])

println("... success of Simulate_MassWithDamper.jl!")
end


text:
# Strand
# ======
#
# Strand type for genomic ranges.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
# Outer constructors
* [`Strand(strand::Char)`](@ref)
* [`Strand(strand::UInt8)`](@ref)

[`Strand`](@ref) can take four kinds of values listed in the next table:

| Symbol | Constant              | Meaning                           |
| :----- | :-------------------- | :-------------------------------- |
| `'?'`  | [`STRAND_NA`](@ref)   | strand is unknown or inapplicable |
| `'+'`  | [`STRAND_POS`](@ref)  | positive strand                   |
| `'-'`  | [`STRAND_NEG`](@ref)  | negative strand                   |
| `'.'`  | [`STRAND_BOTH`](@ref) | non-strand-specific feature       |
"""
primitive type Strand 8 end

Base.convert(::Type{Strand}, strand::UInt8) = reinterpret(Strand, strand)

"""
    Strand(strand::UInt8)
"""
Strand(strand::UInt8) = convert(Strand, strand)
Base.convert(::Type{UInt8}, strand::Strand) = reinterpret(UInt8, strand)

Base.isless(a::Strand, b::Strand) = convert(UInt8, a) < convert(UInt8, b)

"Unknown strand ('?')"
const STRAND_NA   = convert(Strand, 0b000)

"Positive strand ('+')"
const STRAND_POS  = convert(Strand, 0b001)

"Negative strand ('-')"
const STRAND_NEG  = convert(Strand, 0b010)

"Both strand ('.')"
const STRAND_BOTH = convert(Strand, 0b011)

function Base.convert(::Type{Strand}, strand::Char)
    if strand == '+'
        return STRAND_POS
    end

    if strand == '-'
        return STRAND_NEG
    end

    if strand == '.'
        return STRAND_BOTH
    end

    if strand == '?'
        return STRAND_NA
    end

    error("'$(strand)' is not a valid strand")
end

"""
    Strand(strand::Char)
"""
Strand(strand::Char) = convert(Strand, strand)

function Base.convert(::Type{Char}, strand::Strand)
    if strand == STRAND_NA
        return '?'
    end

    if strand == STRAND_POS
        return '+'
    end

    if strand == STRAND_NEG
        return '-'
    end

    if strand == STRAND_BOTH
        return '.'
    end

    error("'$(strand)' is not a valid strand")
end

function Base.show(io::IO, strand::Strand)
    if strand == STRAND_NA
        print(io, "STRAND_NA")
    elseif strand == STRAND_POS
        print(io, "STRAND_POS")
    elseif strand == STRAND_NEG
        print(io, "STRAND_NEG")
    elseif strand == STRAND_BOTH
        print(io, "STRAND_BOTH")
    else
        print(io, "(undefined strand)")
    end
end

function Base.print(io::IO, strand::Strand)
    if strand == STRAND_NA
        print(io, "?")
    elseif strand == STRAND_POS
        print(io, "+")
    elseif strand == STRAND_NEG
        print(io, "-")
    elseif strand == STRAND_BOTH
        print(io, ".")
    else
        error("try to print an invalid strand")
    end
end

function Base.write(io::IO, strand::Strand)
    return write(io, convert(UInt8, strand))
end

function Base.read(io::IO, ::Type{Strand})
    return convert(Strand, read(io, UInt8))
end


text:
# tests the
using RecursiveArrayTools, DifferenceEquations, Test
@testset "matrix as AbstractVectorOfArray" begin
    A = rand(2, 10)
    vv = VectorOfArray([A[:, i] for i in 1:size(A, 2)])
    mv = MatrixVectorOfArray(A)
    @test mv[1] == A[:, 1]
    @test mv[1] == vv[1]
    @test size(mv) == (2, 10)
    @test mv[5][2] == A[2, 5]
    @test mv[5][2] == vv[5][2]
    @test mean(mv) == mean(vv)
    @test map(x -> 2 * x, mv) == map(x -> 2 * x, eachcol(mv.u))
    @test mapreduce(x -> 2 * x, +, mv) == mapreduce(x -> 2 * x, +, vv)
    @test [v for v in mv] == [v for v in vv]  # interator
    mv[5] = [1.0, 0.1]
    vv[5] = [1.0, 0.1]
    @test mv == vv
    mv[2, 6] = 5.9
    vv[2, 6] = 5.9
    @test mv == vv

    mv_sim = similar(mv)
    fill!(mv_sim, 0.0)

    @inferred MatrixVectorOfArray(A)
    @inferred mv[1]
    @inferred mv[2, 2]
    @inferred similar(mv)
end


text:
abstract abstract_newton_direction;
abstract abstract_residuals;

type class_theta
    # mu = 1.0 <=> no mu change
    # mu = 0.0 <=> aggresively decrease mu
    mu::Float64
    dual::Float64
    primal::Float64
end

include("class_homogeneous_residuals.jl")
include("class_homogeneous_newton.jl")

max_forks_repo_name: ohinder/Everest


text:
# Objective functions
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export log_likelihood, ∇log_likelihood, mse, ∇mse

###################################################################################################
# Mean squared error
"""
    f = mse(X, Y)

Mean squared error between arrays/tensor X and Y

See also: [`∇mse`](@ref)
"""
mse(X, Y) = .5f0/size(X, 4)*norm(X - Y, 2)^2

"""
    ∇f = ∇mse(X, Y)

Gradient of the MSE loss with respect to input tensors X and Y.

See also: [`mse`](@ref)
"""
∇mse(X, Y) = 1f0/size(X, 4)*(X - Y)

###################################################################################################
# Log-likelihood

"""
    f = log_likelihood(X; μ=0f0, σ=1f0)

Log-likelihood of X for a Gaussian distribution with given mean μ and variance σ. All elements of X are assumed to be iid.

See also: [`∇log_likelihood`](@ref)
"""
log_likelihood(X; μ=0f0, σ=1f0) = 1f0/size(X, 4)*sum(-.5f0*((X .- μ)/σ).^2)

"""
    ∇f = ∇log_likelihood(X; μ=0f0, σ=1f0)

Gradient of the log-likelihood function with respect to the input tensor X.

See also: [`log_likelihood`](@ref)
"""
∇log_likelihood(X; μ=0f0, σ=1f0) = -1f0/size(X, 4)*(X .- μ)/σ^2


text:
using Cxx

# Adding path to tetgen folder
const path_to_tetgen = "/home/olli/programming/forks/tetgen_src"

# Adding source files
tetgen_source = "/home/olli/programming/forks/tetgen_src/tetgen.cxx"
predicates_source = "/home/olli/programming/forks/tetgen_src/predicates.cxx"

# Adding header dir path
addHeaderDir(path_to_tetgen, kind=C_System)

# include source files
cxxinclude(predicates_source)
cxxinclude(tetgen_source)

# Creating tetgen objects
b_ = @cxxnew tetgenbehavior()
in_ = @cxxnew tetgenio()
out = @cxxnew tetgenio()
addin_ = @cxxnew tetgenio()
bgmin_ = @cxxnew tetgenio()

# including polynomial file
# It's necessary to procide 2 arguments for parse_commandline. "test" is only for that purpose
icxx"""
char *file_name[] = {"test", "example.poly"};
$b_->parse_commandline(2, file_name);
"""

# checking input
icxx"""
if ($b_->refine) {
    if (!$in_->load_tetmesh($b_->infilename)) {
      std::cout << "Something went wrong in refine section?"<< std::endl;
    }
} else {
if (!$in_->load_plc($b_->infilename, (int) $b_->object)) {
      std::cout << "Something went wrong in load plc?"<< std::endl;
    }
}
"""

# Adding points
icxx"""
  if ($b_->insertaddpoints) {
    if (!$addin_->load_node($b_->addinfilename)) {
      $addin_->numberofpoints = 0l;
    }
  }
  if ($b_->metric) {
    if (!$bgmin_->load_tetmesh($b_->bgmeshfilename)) {
      $bgmin_->numberoftetrahedra = 0l;
    }
  }

  if ($bgmin_->numberoftetrahedra > 0l) {
    tetrahedralize($b_, $in_, NULL, $addin_, $bgmin_);
  } else {
    tetrahedralize($b_, $in_, NULL, $addin_, NULL);
  }
"""


text:
function get_ρ(x_e::T, penalty, xmin) where {T}
    if PENALTY_BEFORE_INTERPOLATION
        return density(penalty(x_e), xmin)
    else
        return penalty(density(x_e, xmin))
    end
end

function get_ρ_dρ(x_e::T, penalty, xmin) where {T}
    d = ForwardDiff.Dual{T}(x_e, one(T))
    if PENALTY_BEFORE_INTERPOLATION
        p = density(penalty(d), xmin)
    else
        p = penalty(density(d, xmin))
    end
    g = p.partials[1]
    return p.value, g
end


text:
using MechanicalSketch
import MechanicalSketch: background, background, sethue, O, EM, finish, empty_figure
import MechanicalSketch: lumin, get_current_lumin, get_current_RGB, color_with_lumin
import MechanicalSketch: sethue, fontface, textpath, textextents, PT, fontsize, gsave, grestore
import MechanicalSketch: newpath, translate, m, text, fillpreserve, clip, rect, clipreset
import MechanicalSketch: setline, PALETTE, color_from_palette

let
BACKCOLOR = color_with_lumin(PALETTE[2], 70);
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_10.png"))
    background(BACKCOLOR)
    sethue(PALETTE[2])
end

# Electric power symbol, implemented in 'power.jl' `showpower`

# parameters
power = 3.0kW
maxpower = 5kW |> typeof(power)
fractionpower = power / maxpower
rounded = round(typeof(power), power, digits = 1)
fs = 4EM
symbsize = 1.5fs

restart()

begin
    # Prepare color selection
    luminback = lumin(BACKCOLOR)
    luminfront = get_current_lumin()
    deltalumin = luminback - luminfront
    avglumin = luminfront + 0.5deltalumin
    curcol = get_current_RGB()
    avgcol = color_with_lumin(curcol, luminfront + 0.5 * deltalumin)
    contrastcol = color_with_lumin(curcol, luminfront + 2 * deltalumin);
end;

gsave()
fontsize(symbsize)
fontface("DejaVu Sans")
str = "⚡"
WI, HE = textextents(str)[3:4]
p = O + (-11m, 6m)
begin
    # Outline symbol
    sethue(avgcol)
    fontface("DejaVu Sans Bold")
    fontsize(1.0989010989010988 * symbsize)
    text(str, p - fs.*(0.03, -0.04))
    # Shadow text
    fontface("DejaVu Sans")
    fontsize(fs)
    text(string(rounded), p + (1.5WI, 0) - fs.*(0.03, -0.04))
    # Inside symbol
    sethue(color_with_lumin(color_from_palette("red"), luminback))
    fontsize(symbsize)
    text(str, p)
    # Front text
    fontsize(fs)
    text(string(rounded), p + (1.5WI, 0))
end

# Make a mask
begin
    clipreset()
    newpath()
    translate(p)
    fontsize(symbsize)
    textpath(str)
    fillpreserve()
    clip()
    translate(-p)
end

sethue(curcol)
rect(p, WI , -fractionpower * HE, :fill)
sethue(contrastcol)
setline(PT)
rect(p, WI , -fractionpower * HE, :stroke)
clipreset()
grestore()
finish()
nothing
end


text:
""" Import a file. """
function importFile(name::String)
    dataDir = joinpath(dirname(pathof(gcSolver)), "..", "data")
    return DataFrame(CSV.File(joinpath(dataDir, name)))
end

""" Import all of the data into one dataframe. """
function importData(monomeric = false)
    yData = importFile("WTDimericMutSingleCellData.csv")
    monDF = importFile("MonomericMutSingleCellData.csv")

    append!(yData, monDF)
    if monomeric
        yData = filter(row -> row.Bivalent ∈ [0], yData)
    end

    affDF = importFile("WTmutAffData.csv")
    exprDF = importFile("RecQuantitation.csv")

    exprDF = stack(exprDF, [:Treg, :Thelper, :NK, :CD8]; variable_name = "Cell")
    exprDF = unstack(exprDF, :Receptor, :value)
    select!(exprDF, Not(:Moment))

    fullData = innerjoin(yData, affDF, on = :Ligand => :Mutein)
    fullData = leftjoin(fullData, exprDF, on = :Cell)

    return fullData
end

""" Import the saved fit. """
function importFit()
    return importFile("actualFit.csv")
end

""" Import Date pSTAT conversion factors. """
function importConvFrame()
    return importFile("DateConvFrame.csv")
end

"""Creates Sigma for Var Propagation"""
function getSigma(cellType)
    sigma = zeros(3, 3)

    momentDF = importFile("receptor_moments.csv")
    covDF = importFile("receptor_covariances.csv")

    momentDF = momentDF[!, ["Cell Type", "Receptor", "Variance"]]
    momentDF = groupby(momentDF, ["Cell Type", "Receptor"])
    momentDF = combine(momentDF, :Variance => mean)
    filter!(row -> row["Cell Type"] .== cellType, momentDF)
    for i = 1:3
        sigma[i, i] = momentDF[i, "Variance_mean"]
    end

    covDF = covDF[!, ["Cell Type", "CD25:Receptor", "Covariance"]]
    covDF = groupby(covDF, ["Cell Type", "CD25:Receptor"])
    covDF = combine(covDF, :Covariance => mean)
    filter!(row -> row["Cell Type"] .== "Treg", covDF)
    for i = 2:3
        sigma[1, i] = sigma[i, 1] = covDF[i - 1, "Covariance_mean"]
    end

    return sigma
end


text:
struct Ice_output_config
    name::String
    path::String
    nslices::Int
    nfiles::Int
    nechoes::Int
    nchannels::Int
    dtype::Type
    size::Tuple{Int, Int}
end

function Ice_output_config(name, path, nslices, nfiles; nechoes=1, nchannels=1, dtype=Int16)
    return Ice_output_config(name, path, nslices, nfiles, nechoes, nchannels, dtype, getsize(path))
end

function get_setting(T, lines, setting; offset=3, default=0)
    for iLine in 1:length(lines)
        if occursin(setting, lines[iLine])
            try
                return parse(T, lines[iLine + offset])
            catch
                return default
            end
        end
    end
    return default
end

function read_volume(cfg)
    volume = create_volume(cfg)

    for i in 1:cfg.nfiles
        num = lpad(i, 5, '0')
        imahead = joinpath(cfg.path, "MiniHead_ima_$num.IceHead")
        file = joinpath(cfg.path, "WriteToFile_$num.ima")

        if occursin(cfg.name, read(imahead, String))
            vol = Array{cfg.dtype}(undef, cfg.size...)
            read!(file, vol)
            lines = readlines(imahead)
            eco = get_setting(Int, lines, "EchoNumber"; default=1)
            slc = getslice(lines)
            rescale_slope = get_setting(Float32, lines, "RescaleSlope"; offset=4, default=1)
            rescale_intercept = get_setting(Float32, lines, "RescaleIntercept"; offset=4, default=0)
            volume[:,:,slc,eco] .= (vol .+ rescale_intercept) .* rescale_slope
        end
    end

    return volume
end

function getslice(lines)
    slc = get_setting(Int, lines, "Actual3DImaPartNumber"; default=nothing)
    if (isnothing(slc))
        slc = get_setting(Int, lines, "AnatomicalSliceNo")
    end
    return slc + 1
end

function getsize(path)
    lines = readlines(joinpath(path, "MiniHead_ima_00001.IceHead"))
    return (get_setting(Int, lines, "NoOfCols"), get_setting(Int, lines, "NoOfRows"))
end

function create_volume(cfg)
    return zeros(Float32, cfg.size..., cfg.nslices, cfg.nechoes)
end


text:
#=
    Contents: Non app-specific handlers.
=#

"""
Process login credentials.
Login request must be a POST request.
If login is successful, redirect to success_redirect, else return bad_request.
"""
function login!(req, res)
    acdata = config[:acdata]
    if req.method == "POST"                                           # Require that login requests are POST requests
        username, password = extract_username_password(req)
	if login_credentials_are_valid(username, password, acdata)    # Successful login: Redirect
	    # Create session
	    if using_clientside_sessions()
		session = session_create!(username)
		session_write!(res, session)
	    else
		session_create!(res, username)
	    end

	    # Redirect
	    res.data = get_redirect_location(username, :login, :success_redirect)
	else                                                          # Unsuccessful login: Return fail message
	    res.status = 400    # Bad request
	    res.data   = config[:login][:fail_msg]
	end
    else
	badrequest!(res)
    end
end

"User has clicked the logout button/link: Redirect to redirect_path."
function logout!(req, res)
    is_not_logged_in(req) && (notfound!(res); return)
    username = get_username(req)
    loc      = get_redirect_location(username, :logout, :redirect)
    redirect!(res, loc)
    if using_clientside_sessions()
	session_delete!(res)
    else
	session_id = read_sessionid(req)
	session_delete!(res, session_id)
    end
end

"""
User has submitted a new password.

Lockout occurs if the supplied original password is wrong too many times (defined by cfg["pwdreset_cfg"]).
"""
function process_pwdreset!(req, res)
    is_not_logged_in(req) && (notfound!(res); return)
    acdata = config[:acdata]
    if req.method == "POST"
        form_id, nonce, current_pwd, new_pwd, new_pwd2 = extract_currpwd_newpwd_newpwd2(req)

	# Check whether nonce is valid
	nonce_is_valid = false
	if using_clientside_sessions()
	    session = session_read(req)
	    dct     = session["forms"]
	    if haskey(dct, form_id) && dct[form_id] == nonce
		nonce_is_valid = true
	    end
	else
	    session_id = read_sessionid(req)
	    if nonce == session_get(session_id, "forms", form_id)
		nonce_is_valid = true
	    end
	end

	# If nonce is valid then check login credentials
	if nonce_is_valid
	    username  = get_username(req)
	    if login_credentials_are_valid(username, current_pwd, acdata) && new_pwd == new_pwd2    # Successful password reset: Redirect
		set_password!(username, new_pwd)
		res.data = get_redirect_location(username, :pwdreset, :success_redirect)
	    else                                                # Unsuccessful password reset: Return 400: Bad Request
		res.status = 400    # Bad request
		res.data   = config[:pwdreset][:fail_msg]
	    end
	else
	    badrequest!(res)
        end
    else
        badrequest!(res)
    end
end

"Page displayed when user clicks Reset Password link."
function reset_password!(req, res)
    is_not_logged_in(req) && (notfound!(res); return)
    if using_clientside_sessions()
	session  = session_read(req)
	res.data = pwdreset_form(session)
    else
	session_id = read_sessionid(req)
	res.data   = pwdreset_form(session_id)
    end
end

# EOF


text:
@testset "Construction" begin
    @test OneHotVector(10, 1) == OneHotVector(10, 1)
    @test_throws ArgumentError OneHotVector(1, 10)
    @test_throws ArgumentError OneHotVector(0, 1)

    arr = zeros(Bool, 4)
    @test_throws ArgumentError OneHotVector(arr)
    arr[2] = true
    @test OneHotVector(arr) == OneHotVector(4, 2)
    arr[4] = true
    @test_throws ArgumentError OneHotVector(arr)
end

@testset "Basics" begin
    @test eltype(OneHotVector) == Bool
    @test length(OneHotVector(1, 1)) == 1
    @test length(OneHotVector(23241, 9)) == 23241
end

@testset "Getindex" begin
    v = OneHotVector(11, 3)
    @test_throws BoundsError v[0]
    @test_throws BoundsError v[-5]
    @test_throws BoundsError v[15]

    @test !(v[11])
    @test v[3]

    @test first(OneHotVector(1, 1))
end

@testset "Addition/subtraction" begin
    for len in [1, 2, 9, 11, 115, 500]
        for i in 1:3
            v1 = OneHotVector(len, rand(1:len))
            v2 = rand(-10:10, len)
            @test v1 + v2 == [(i+j) for (i,j) in zip(v1, v2)]
            @test v1 - v2 == [(i-j) for (i,j) in zip(v1, v2)]
            @test v2 + v1 == [(i+j) for (i,j) in zip(v2, v1)]
            @test v2 - v1 == [(i-j) for (i,j) in zip(v2, v1)]
        end
    end
end

@testset "Misc" begin
    for len in [1, 2, 4, 9, 11, 22, 115, 500]
        for i in 1:5
            v = OneHotVector(len, rand(1:len))
            v´ = collect(v)

            @test v[:] === v
            @test count(identity, v) == count(v) == 1

            for F in [argmax, argmin, sum, count, allunique,
                      x -> findfirst(isodd, x), x -> findfirst(isequal(4), x),
                      x -> findfirst(y -> false, x), x -> findfirst(~, x)]
                @test F(v) == F(v´)
            end

            for F in [reverse, x -> circshift(x, 1), x -> circshift(x, -5),
                      x -> circshift(x, 500), x -> filter(iseven, x),
                      x -> filter(y -> true, x), x -> filter(isequal(3), x)]
                @test collect(F(v)) == F(v´)
            end
        end
    end
end


text:
module GaussianBasis

using Molecules
using Formatting
import Molecules: Atom, symbol

export BasisFunction, BasisSet

@doc raw"""
     BasisFunction

Basis object holding basis function information.

# Fields

| Name        | Type                               |   Description |
|:------------|:-----------------------------------|:-----------------------------------------------------------|
|`l`          | `Cint`                             | Angular momentum number |
|`coef`       | `Vector{Cdouble}`                  | Array holding expansion coefficients for the primitives |
|`exp`        | `Vector{Cdouble}`                  | Array holding exponents for primitives          |

# Examples

```julia
julia> bf = BasisFunction(1, [1/√2, 1/√2], [5.0, 1.2])
P shell with 3 basis built from 2 primitive gaussians

χ₁₋₁ =    0.7071067812⋅Y₁₋₁⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₋₁⋅r¹⋅exp(-1.2⋅r²)

χ₁₀  =    0.7071067812⋅Y₁₀⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₀⋅r¹⋅exp(-1.2⋅r²)

χ₁₁  =    0.7071067812⋅Y₁₁⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₁⋅r¹⋅exp(-1.2⋅r²)
```
"""
struct BasisFunction
    l::Cint
    coef::Vector{Cdouble}
    exp::Vector{Cdouble}
end

include("BasisParser.jl")
include("BasisSet.jl")
include("Libcint.jl")
include("Integrals.jl")
include("Gradients.jl")

end # module


text:
using CollisionDetection
using StaticArrays
using Base.Test
using JLD
CD = CollisionDetection

@test CD.childsector([1,1,1],[0,0,0]) == 7
@test CD.childsector([-1,1,-1],[0,0,0]) == 2

@test CD.childsector([-1,1,-1,1],[0,0,0,0]) == 10
@test CD.childsector([1,1],[0,0]) == 3

@test CD.childcentersize(SVector(0.,0.,0.), 1., 3) == (SVector(0.5,0.5,-0.5), 0.5)
@test CD.childcentersize(SVector(0.,0.,0.), 2., 7) == (SVector(1.0,1.0,1.0), 1.0)

@test CD.childcentersize(SVector(1.,1.,1.,1.), 2., 9) == (SVector(2.,0.,0.,2.), 1.)

d, n = 3, 200
data = SVector{d,Float64}[rand(SVector{d,Float64}) for i in 1:n]
radii = abs.(rand(n))
tree = CD.Octree(data, radii)

@test length(tree) == n

println("Testing allways true iterator")
sz = 0
for box in CD.boxes(tree)
    sz += length(box)
end
@test sz == n

println("Testing allways false iterator")
sz = 0
for box in CD.boxes(tree, (c,s)->false)
    sz += length(box.data)
end
@test sz == 0

println("Testing query for box in sector [7,0]")
sz = 0
for box in CD.boxes(tree, (c,s)->CD.fitsinbox([1,1,1]/8,0,c,s))
    sz += length(box)
end
@show sz

T = Float64
P = SVector{2,T}
n = 40000
data = P[rand(P) for i in 1:n]
radii = abs.(rand(T,n))
tree = CD.Octree(data, radii)

println("Testing allways true iterator [2D]")
sz = 0
for box in CD.boxes(tree)
    sz += length(box)
end
@test sz == n

println("Testing allways false iterator [2D]")
sz = 0
for box in CD.boxes(tree, (c,s)->false)
    sz += length(box.data)
end
@test sz == 0

println("Testing query for box in sector [7,0], [2D]")
sz = 0
for box in CD.boxes(tree, (c,s)->CD.fitsinbox([1,1]/8,0,c,s))
    sz += length(box)
end
@show sz

## Test the case of degenerate bounding boxes
fn = joinpath(dirname(@__FILE__),"assets","front.jld")
V,F = jldopen(fn,"r") do file
        read(file, "V"), read(file, "F")
end

fn = joinpath(dirname(@__FILE__),"assets","back.jld")
U,G = jldopen(fn,"r") do file
        read(file, "V"), read(file, "F")
end

radii = zeros(length(V))
tree = Octree(V, radii)

tree_ctr, tree_hs = tree.center, tree.halfsize

u = U[3]
atol = sqrt(eps())
@test norm(u-V[4]) < atol
pred(c,s) = fitsinbox(Array(u), 0.0, c, s+atol)
@test pred(tree_ctr, tree_hs)
it = boxes(tree, pred)
st = start(it)
@test !done(it, st)

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection


text:
using CodeTracking
using JuliaInterpreter: getfile, locals
import ..Atom: fullpath, strlimit

function datatip(word, path, row, column, state::DebuggerState = STATE)
  # frame validation
  state.frame === nothing && return nothing

  frame = active_frame(state)
  scope = frame.framecode.scope

  # only work for methods
  scope isa Module && return nothing

  # path identity check
  path != fullpath(getfile(frame)) && return nothing

  # line number check
  defstr, startline = definition(String, scope)
  endline = startline + sum(c === '\n' for c in defstr)
  startline <= row <= endline || return nothing

  # local bindings
  for v in locals(frame)
    str = string(v.name)
    if str == word
      valstr = repr(MIME("text/plain"), v.value, context = :limit => true)
      return [Dict(:type => :snippet, :value => strlimit(valstr, 1000))]
    elseif str ∈ split(word, '.') # interpret dot accessors
      valstr = repr(MIME("text/plain"), interpret(word), context = :limit => true)
      return [Dict(:type => :snippet, :value => strlimit(valstr, 1000))]
    end
  end

  return nothing # when no local binding exists
end


text:
###################################################
## networkviz.jl
## add point&click node information to visual (default visualizer)
###################################################

"""
    `NetworkViz`: Basic NetworkVisualizer. Shows node information in title bar after click
"""
mutable struct NetworkViz <: NetworkVisualizer
    # Mandatory attributes
    network::Network
    window::RenderWindow
    view::View
    nodes::Vector{CircleShape}
    roads::Dict{Tuple{Int,Int},Line}
    nodeRadius::Float64
    colors::VizColors

    "node positions KD-tree"
    tree::KDTree
    "selected node"
    selectedNode::Int

    "contructor: initialize kd-tree"
    function NetworkViz(n::Network; colors::VizColors=RoadTypeColors())
        obj = new()
        # KD-tree with positions
        dataPos = Array{Float64,2}(2,length(n.nodes))
        for (i,node) in enumerate(n.nodes)
           dataPos[1,i] = node.x
           dataPos[2,i] = node.y
        end

        obj.network = n
        obj.tree = KDTree(dataPos)
        obj.selectedNode = 1
        obj.colors = colors
        return obj
    end
end

function visualEvent(v::NetworkViz, event::Event)
    if get_type(event) == EventType.MOUSE_BUTTON_PRESSED && get_mousebutton(event).button == MouseButton.LEFT
        x,y = get_mousebutton(event).x, get_mousebutton(event).y
        coord = pixel2coords(v.window,Vector2i(x,y))

        id = knn(v.tree,[Float64(coord.x),-Float64(coord.y)],1)[1][1]
        set_fillcolor(v.nodes[v.selectedNode], nodeColor(v.colors, v.network.nodes[v.selectedNode]))
        set_fillcolor(v.nodes[id], SFML.red)
        v.selectedNode = id
        set_title(v.window, "Node : $id in: $(in_neighbors(v.network.graph,id)) out: $(out_neighbors(v.network.graph,id))")
    end
end

function visualRedraw(v::NetworkViz)
    set_fillcolor(v.nodes[v.selectedNode], SFML.red)
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module JuliaZH

import Base.Docs: DocStr

const keywords = Dict{Symbol, DocStr}()

function __init__()
    ENV["REPL_LOCALE"] = "zh_CN"
    if ccall(:jl_generating_output, Cint, ()) == 0
        include(joinpath(@__DIR__, "i18n_patch.jl"))
        include(joinpath(@__DIR__, "basedocs.jl"))
    end
end

end # module


text:
abstract type AbstractMonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis end

Base.length(basis::AbstractMonomialBasis) = length(basis.monomials)
Base.copy(basis::AbstractMonomialBasis) = typeof(basis)(copy(basis.monomials))

MP.nvariables(basis::AbstractMonomialBasis) = MP.nvariables(basis.monomials)
MP.variables(basis::AbstractMonomialBasis) = MP.variables(basis.monomials)
MP.monomialtype(::Type{<:AbstractMonomialBasis{MT}}) where MT = MT

empty_basis(MB::Type{<:AbstractMonomialBasis{MT}}) where {MT} = MB(MP.emptymonovec(MT))
function maxdegree_basis(B::Type{<:AbstractMonomialBasis}, variables, maxdegree::Int)
    return B(MP.monomials(variables, 0:maxdegree))
end
function basis_covering_monomials(B::Type{<:AbstractMonomialBasis}, monos::AbstractVector{<:AbstractMonomial})
    return B(monos)
end

# The `i`th index of output is the index of occurence of `x[i]` in `y`,
# or `0` if it does not occur.
function multi_findsorted(x, y)
    I = zeros(Int, length(x))
    j = 1
    for i in eachindex(x)
        while j ≤ length(y) && x[i] < y[j]
            j += 1
        end
        if j ≤ length(y) && x[i] == y[j]
            I[i] = j
        end
    end
    return I
end

function merge_bases(basis1::MB, basis2::MB) where MB<:AbstractMonomialBasis
    monos = MP.mergemonovec([basis1.monomials, basis2.monomials])
    I1 = multi_findsorted(monos, basis1.monomials)
    I2 = multi_findsorted(monos, basis2.monomials)
    return MB(monos), I1, I2
end

"""
    struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
        monomials::MV
    end

Monomial basis with the monomials of the vector `monomials`.
For instance, `MonomialBasis([1, x, y, x^2, x*y, y^2])` is the monomial basis
for the subspace of quadratic polynomials in the variables `x`, `y`.
"""
struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractMonomialBasis{MT, MV}
    monomials::MV
end
MonomialBasis(monomials::AbstractVector) = MonomialBasis(MP.monovec(monomials))

MP.polynomialtype(::Union{MonomialBasis{MT}, Type{<:MonomialBasis{MT}}}, T::Type) where MT = MP.polynomialtype(MT, T)
MP.polynomial(f::Function, mb::MonomialBasis) = MP.polynomial(f, mb.monomials)

MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type) = MP.polynomial(Q, mb.monomials, T)

function MP.coefficients(p, ::Type{<:MonomialBasis})
    return MP.coefficients(p)
end


text:
# ***** FLOATS and INTS *****
# x = 3.2     # float
# y = 3       # int
# z = 3.      # float

# println("x: ", typeof(x))
# println("y: ", typeof(y))
# println("z: ", typeof(z))
# println("sum: ", typeof(x + y), " ", x + y)

# ***** STRINGS and CHARS *****
# a = 'a'         # char
# word = "hello"  # string
# wrong = 'word'

# ***** BOOLEANS *****
# t = true
# f = false
# println(t, " ", typeof(t))
# println(f, " ", typeof(f))

# ***** control flow with if, elseif, else *****
x = -5
y = 0

if (x > 0) && (y > 0)
    println("x and y are positive")
elseif (x < 0) && (y < 0)
    println("x and y are negative")
elseif (x > 0) || (y > 0)
    println("x or y are positive")
else
    println("x or y are 0")
end

# ***** control flow with ternary operator *****
x > 0 ? println("x is positive") : println("x is not positive")

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# Copyright (c) 2019 Guido Kraemer
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

module BitcoinPrimitives

using BitConverter, Secp256k1, Ripemd, SHA

export
    CompactSizeUInt, Outpoint,
    TxIn, TxOut, Tx,
    Block, Header,
    Script, Witness,
    serialize, iscoinbase, coinbase_height,
    script, type,
    hash256,
    target, difficulty, check_pow

include("lib/CompactSizeUInt.jl")
include("transaction/Tx.jl")
include("block/Block.jl")
include("lib/hash.jl")

end # module


text:
# SDPA file format reader-writer

# Redefine the following method to branch on the filename extension of `optimizer` supports more formats
MOI.write(optimizer::SOItoMOIBridge, filename::String) = MOI.write(optimizer.sdoptimizer, filename)
MOI.write(optimizer::AbstractSDOptimizer, filename::String) = writesdpa(optimizer, filename)
MOI.read!(optimizer::SOItoMOIBridge, filename::String) = MOI.read!(optimizer.sdoptimizer, filename)
MOI.read!(optimizer::AbstractSDOptimizer, filename::String) = readsdpa!(optimizer, filename)

function writesdpa(optimizer::AbstractSDOptimizer, filename::String)
    endswith(filename, ".sdpa") || error("Filename must end with .sdpa")
    open(filename, "w") do io
        nconstrs = getnumberofconstraints(optimizer)
        println(io, nconstrs)
        nblocks = getnumberofblocks(optimizer)
        println(io, nblocks)
        for blk in 1:nblocks
            print(io, getblockdimension(optimizer, blk))
            if blk != nblocks
                print(io, ' ')
            end
        end
        println(io)
        for c in 1:nconstrs
            print(io, getconstraintconstant(optimizer, c))
            if c != nconstrs
                print(io, ' ')
            end
        end
        println(io)
        for (val, blk, i, j) in getobjectivecoefficients(optimizer)
            println(io, "0 $blk $i $j $val")
        end
        for c in 1:nconstrs
            for (val, blk, i, j) in getconstraintcoefficients(optimizer, c)
                println(io, "$c $blk $i $j $val")
            end
        end
    end
end

nextline(io::IO) = chomp(readline(io))

function readsdpa!(optimizer::AbstractSDOptimizer, filename::String)
    endswith(filename, ".sdpa") || error("Filename '$filename' must end with .sdpa")
    open(filename, "r") do io
        line = nextline(io)
        while line[1] == '"' || line[1] == '*' # Comment
            line = nextline(io)
        end
        nconstrs = parse(Int, line)
        nblocks = parse(Int, nextline(io))
        blkdims = parse.(Int, split(nextline(io)))
        init!(optimizer, blkdims, nconstrs)
        T = coefficienttype(optimizer)
        constraint_constants = parse.(T, split(nextline(io)))
        for c in 1:nconstrs
            setconstraintconstant!(optimizer, constraint_constants[c], c)
        end
        while !eof(io)
            line = nextline(io)
            isempty(line) && break
            s = split(line)
            c = parse(Int, s[1])
            0 ≤ c ≤ nconstrs || error("Invalid constraint index $c in '$filename', it should be an integer between 0 and $nconstrs")
            blk = parse(Int, s[2])
            1 ≤ blk ≤ nblocks || error("Invalid block index $blk in '$filename', it should be an integer between 0 and $nblocks")
            i = parse(Int, s[3])
            j = parse(Int, s[4])
            val = parse(T, s[5])
            if iszero(c)
                setobjectivecoefficient!(optimizer, val, blk, i, j)
            else
                setconstraintcoefficient!(optimizer, val, c, blk, i, j)
            end
        end
    end
end


text:
#posit_cgen_helpers.jl.  A bunch of helper functions that make the cgen
#process smoother.

#symbols for each number
const sym = Dict(1 => :a, 2 => :b, 3 => :c, 4 => :d)

float_for_posit_size(n) = (n < 32) ? 32 : 64
c_ftype(n) = Dict(32 => "float", 64 => "double")[n]

#C internal use conversion functions
to_f(n, es, v) = "p$(n)e$(es)_to_f($v)"
to_p(n, es, v) = "f_to_p$(n)e$(es)($v)"

#C datatype names
p(n,es) = "p$(n)e$(es)_t"
#C function names
p(n,es,op) = "p$(n)e$(es)_$(op)"
#C++  class names
c(n,es) = "P$(n)e$(es)"

function c_literal(n::Unsigned)
  if isa(n, UInt8)
    "((uint8_t) 0x$(hex(n,2)))"
  elseif isa(n, UInt16)
    "((uint16_t) 0x$(hex(n,4)))"
  elseif isa(n, UInt32)
    "0x$(hex(n,8))L"
  elseif isa(n, UInt64)
    "0x$(hex(n,16))LL"
  end
end

const bits_to_UInt = Dict(8=>UInt8, 16=>UInt16, 32=>UInt32, 64=>UInt64)

top_bits(n) = one(bits_to_UInt[n]) << (n - 1)
zero_bits(n) = zero(bits_to_UInt[n])


text:
struct ObservablePair{S, T} <: AbstractObservable{T}
    first::AbstractObservable{S}
    second::AbstractObservable{T}
    f
    g
    links::Tuple{ObserverFunction,ObserverFunction}

    function ObservablePair(first::AbstractObservable{S}, second::AbstractObservable{T}; f = identity, g = identity) where {S, T}

        # the two observables should trigger each other, but only in one direction
        # as otherwise there will be an infinite loop of updates
        done = Ref(false)
        link1 = on(first) do val
            if !done[]
                done[] = true
                second[] = f(val)
                done[] = false
            end
        end
        link2 = on(second) do val
            if !done[]
                done[] = true
                first[] = g(val)
                done[] = false
            end
        end

        new{S, T}(first, second, f, g, (link1, link2))
    end
end

ObservablePair(first::AbstractObservable; f = identity, g = identity) =
    ObservablePair(first, Observable{Any}(f(first[])); f = f, g = g)

observe(o::ObservablePair) = o.second

function off(o::ObservablePair)
    for i in 1:2
        off(o[i], o.links[i])
    end
end

Base.iterate(p::ObservablePair, i=1) = i > 2 ? nothing : (getfield(p, i), i + 1)
Base.getindex(p::ObservablePair, i::Int) = getfield(p, i)
Base.firstindex(p::ObservablePair) = 1
Base.lastindex(p::ObservablePair) = 2
Base.length(p::ObservablePair) = 2
Base.first(p::ObservablePair) = p.first
Base.last(p::ObservablePair) = p.second


text:
using LJH
using DataStructures
using Nullables
using Test

@testset "ljh3" begin

header_extra = OrderedDict(["a"=>"b","c"=>"d","e"=>Dict("ea"=>"eb")])
fname = tempname()
fw = LJH.create3(fname, 9.6e-6, header_extra)

traces = [rand(UInt16,rand(1:1000)) for i=1:10000];
first_rising_samples = [rand(1:length(trace)) for trace in traces];
frame1indexs = 1:length(traces);
timestamp_usecs = frame1indexs.*1000;

for i = 1:length(traces)
    write(fw, traces[i],first_rising_samples[i],
     frame1indexs[i], timestamp_usecs[i])
end
@test fw[1].data == traces[1]
@test fw[77].data == traces[77]
close(fw)

f = LJH3File(fname)
f2 = ljhopen(fname)
records = [record for record in f];
@test fw.index == f.index
@test traces == [record.data for record in records]
@test first_rising_samples == [record.first_rising_sample for record in records]
@test collect(frame1indexs) == [record.frame1index for record in records]
@test timestamp_usecs == [record.timestamp_usec for record in records]
@test records == collect(f)
@test f[1].data == traces[1]
@test f[77].data == traces[77]
@test all(f[key]==f.header[key] for key in keys(f))
@test all(f[key]==f.header[key] for key in keys(header_extra))
@test length(f2) == length(traces) # test that length works before collect is called on f2
@test f2.index == f.index
LJH.seekto(f,length(traces)-1)
@test get(LJH.tryread(f))==records[end-1]
@test get(LJH.tryread(f))==records[end]
@test isnull(LJH.tryread(f))
@test LJH.filename(f) == LJH.filename(fw) == fname
close(f)
close(f2)
rm(fname)
end

@testset "LJH3Record and LJHRecord joint API" begin
    FrameTime, PretrigNSamples, NRow, row = 9.6e-6, 100, 30, 1
    # rowcount of LJH2 files = framecount*NRow+row, framecount=frame1index
    r = LJH.LJHRecord{FrameTime, PretrigNSamples, NRow}(1:1000,10*NRow+row,20)
    r3 = LJH.LJH3Record{FrameTime}(1:1000,PretrigNSamples,10,20)
    @test LJH.data(r) == LJH.data(r3)
    @test LJH.frameperiod(r) == LJH.frameperiod(r3)
    @test LJH.first_rising_sample(r) == LJH.first_rising_sample(r3)
    @test LJH.frame1index(r) == LJH.frame1index(r3)
    @test LJH.timestamp_usec(r) == LJH.timestamp_usec(r3)
    @test length(r) == length(r3)
end

# using BenchmarkTools
# @benchmark read(seekstart(f.io)) setup=(f=LJH3File(fname)) teardown=close(f) evals=1
# @benchmark collect(f) setup=(f=LJH3File(fname)) teardown=close(f) evals=1


text:
# sample code for the future

function usage(args, core)
 println("julia> repl()")
 println(" > core c1 data/his001.cnf   # readcore from file to c1")
 println(" > core c2 data/his002.cnf")
 println(" > core c3 data/his001.cnf")
 println(" > prcore c1                 # printcore c1")
 println(" > prove c1 15 5             # simpleprove c1 with ns=15 and nc=5")
 println(" > prove c3 15 5")
 println(" > proof c1                  # show proofs of c1")
 println(" > mgu c1                    # show mgus of c1")
 println(" > mgu c1 t                  # show mgus of c1 with input vars")
 println()
 println("or script f1 or string f1 already made.")
 println("julia> io = open(f1, read=true) or ")
 println("julia> io = IOBuffer(\"core c1 data/his001.cnf\\nprove c1 5 5\")")
 println("julia> repl(io)")
 println("... and so on...")
 println("--- end of usage")
end

function help()
 println("help,?       : shows help")
 println("end          : finish this loop")
 println("core xx cnf  : read core")
 println("prove xx n m : simpleprove(xx, n, m)")
 println("ls   xx      : ls xx")
 println("proofs xx    : show proofs of core")
 println("mgu   xx     : show mgus of core")
 println("cores        : show loaded core names")
 println("prcore xx    : print core")
 println("usage        : show sample command sequence")
 println("--- end of help")
end

function separate(str)
 sps = findall(isspace, str)
 s = 1
 tokens = []
 for e in sps
   push!(tokens, strip(str[s:e]))
   s=e
 end
 push!(tokens, strip(str[s:end]))
 tokens
end

function proofsfn(args, cores)
 name=args[1]
 printproofs1(cores[name])
end

function mgusfn(args, cores)
 name = args[1]
 if length(args) == 2
   printmgus(cores[name], true)
 else
   printmgus(cores[name])
 end
end

function provefn(args, cores)
 if length(args) >= 3
   name,sns,snc = args
 else
   println("no max num of step and max num of contradictions")
   return
 end
 ns, nc = parse(Int, sns), parse(Int, snc)
 try
  rcore = simpleprovercore(cores[name], ns, nc)[2]
  cores[name] = rcore
 catch
  println("error in proving")
 end
end

function lsfn(args, cores)
 if 0 == length(args)
  run(`ls`)
 else
  run(`ls $(args[1])`)
 end
end

function corefn(args, cores)
  name=args[1]
  cnf=args[2]
  println("cnf=$cnf => $name")
  core = readcore(cnf)
  cores[name] = core
end

# direct called functions
function prcore(args, cores)
  printcore(cores[args[1]])
end

function cores(args, cores)
  println(keys(cores))
end

commandmap = Dict(
 "help"  => :help,
 "ls"    => :lsfn,
 "core"  => :corefn,
 "prove" => :provefn,
 "proof" => :proofsfn,
 "mgu"   => :mgusfn
)

function repl(io=stdin)
#function repl()
 println("Now repl is exerimental.")
 println("help show what commands are available")
 println("end exist the loop")
 println()

 print("> ")

 cores = Dict()
 for line in eachline(io)
 try
  cmdline = separate(line)
  if cmdline[1] == "end"; return
  elseif cmdline[1] == "";
  elseif cmdline[1] == "#";
  elseif cmdline[1] == "?"; help()
  elseif cmdline[1] == "help"; help()
  else
   args = cmdline[2:end]
   cmd = strip(cmdline[1])
   if cmd in keys(commandmap)
     eval(Expr(:call, commandmap[cmd], args, cores))
   else
     eval(Expr(:call, Symbol(cmd), args, cores))
   end
  end
  catch e
   println(e)
  end
 print("> ")
 end
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
include("tif/eval_tif.jl")
include("constraints/eval_constraints.jl")
include("fill_up/fill_up.jl")
include("maximin.jl")
include("minimax.jl")
include("no_obj.jl")
include("custom.jl")


text:
export StoppingCriterion
include("problem.jl")
include("options.jl")

include("debugOptions.jl")
include("recordOptions.jl")

include("costPlan.jl")
include("gradient_plan.jl")
include("hessianPlan.jl")
include("proximalPlan.jl")
include("subgradient_plan.jl")

include("stoppingCriterion.jl")
include("stepsize.jl")


text:
# test serialization of distributions

using IncrementalInference
using Test

##

@testset "Packing Categorical" begin
##

ctg = Categorical(5)

packed = packDistribution(ctg)

@test packed isa PackedSamplableBelief
@test packed isa IncrementalInference.PackedCategorical

upck = unpackDistribution(packed)

@test upck isa Categorical
@test isapprox(ctg, upck)

##
end

@testset "Packing Normal" begin
##

nor = Normal()

packed = packDistribution(nor)

@test packed isa PackedSamplableBelief
@test packed isa PackedNormal

upck = unpackDistribution(packed)

@test upck isa Normal
@test isapprox(nor, upck)

##
end

@testset "Packing for FullNormal" begin
##

mv = MvNormal([0.0;1.0], [1 0.1; 0.1 1])
packed = packDistribution(mv)

@test packed isa PackedSamplableBelief
@test packed isa PackedFullNormal

upck = unpackDistribution(packed)

@test upck isa FullNormal
@test isapprox(mv, upck)

##
end

@testset "Packing for DiagNormal" begin
##

mv = MvNormal([0.0; 1.0], [4.0; 4.0])
packed = packDistribution(mv)

@test packed isa PackedSamplableBelief
@test packed isa PackedDiagNormal

upck = unpackDistribution(packed)

@test upck isa DiagNormal
@test isapprox(mv, upck)

##
end

@testset "Packing for ZeroMeanFullNormal" begin
##

mv = MvNormal([1 0.1; 0.1 1])
packed = packDistribution(mv)

@test packed isa PackedSamplableBelief
@test packed isa PackedZeroMeanFullNormal

upck = unpackDistribution(packed)

@test upck isa ZeroMeanFullNormal
@test isapprox(mv, upck)

##
end

@testset "Packing for ZeroMeanDiagNormal" begin
##

mv = MvNormal([4.0;4.0])
packed = packDistribution(mv)

@test packed isa PackedSamplableBelief
@test packed isa PackedZeroMeanDiagNormal

upck = unpackDistribution(packed)

@test upck isa ZeroMeanDiagNormal
@test isapprox( mv.Σ.diag, upck.Σ.diag )
# @test isapprox(mv, upck)

##
end

@testset "Packing for AliasingScalarSampler" begin
##

bss = AliasingScalarSampler([1.0;2.0], [0.6;0.4])
packed = packDistribution(bss)

@test packed isa PackedSamplableBelief
@test packed isa IncrementalInference.PackedAliasingScalarSampler

upck = unpackDistribution(packed)

@test upck isa AliasingScalarSampler
@test isapprox( bss.domain, upck.domain )
@test isapprox( bss.weights.values, upck.weights.values )
@test isapprox( bss.weights.sum, upck.weights.sum )

##
end

@testset "Packing for ManifoldKernelDensity" begin
##

T = ContinuousEuclid{2}
pts = [randn(2) for _ in 1:50]

mkd = manikde!(T, pts)
packed = packDistribution(mkd)

@test packed isa PackedSamplableBelief
@test packed isa IncrementalInference.PackedManifoldKernelDensity

upck = unpackDistribution(packed)

@test upck isa ManifoldKernelDensity

@test isapprox(mkd, upck)

##
end

@testset "Packing of Rayleigh" begin
##

r = Rayleigh(1.1)
r_ = packDistribution(r)

@test r_ isa PackedSamplableBelief
@test r_ isa PackedRayleigh

r__ = unpackDistribution(r_)

@test r__ isa Rayleigh
@test isapprox(r.σ, r__.σ)

##
end

## Legacy tests

# @testset "hard-coded test of PackedPrior to Prior" begin
# ##

# pt = PackedPrior("KDE:100:[1.5015]:[-98.8276 -101.803 -98.1296 -96.1897 -99.3076 -99.1881 -101.721 -101.332 -100.431 -97.7293 -96.7652 -99.3806 -95.5593 -104.237 -94.9318 -101.691 -102.255 -98.9559 -99.3386 -99.2361 -102.483 -102.896 -97.0244 -98.9643 -99.4457 -101.298 -103.263 -2.75251 5.14065 0.327863 3.60042 -0.604114 -0.0564047 -0.804898 3.05469 1.4974 1.34657 2.22745 4.78117 1.89485 -0.48091 6.63068 0.63013 -3.11422 1.73705 5.22904 -1.73223 2.47453 1.10584 -0.0179944 3.65585 4.50016 -1.95042 98.2664 98.9983 103.748 100.789 98.4127 101.397 104.364 102.125 96.3685 103.59 99.0581 100.913 101.461 105.211 103.513 99.3325 101.201 98.05 103.508 99.9785 104.624 100.202 100.258 101.579 96.6931 95.4181 299.02 296.804 301.322 298.127 299.578 298.36 296.339 300.156 299.641 297.731 299.822 296.941 295.857 299.482 302.531 301.875 302.192 301.999 300.634 294.084 300.44]")

# tt = convert(BallTreeDensity, pt.Z)
# @test isa(tt, BallTreeDensity)

# F = Prior

# upt = convert(F, pt)

# @test isa(upt, Prior)

# ## TODO add more tests
# end

#


text:
# This code is from the Modia project and is licensed as follows:
# https://github.com/ModiaSim/Modia.jl/blob/b61daad643ef7edd0c1ccce6bf462c6acfb4ad1a/LICENSE

################################################
#
# Functions to tear systems of equations
#
# Author: Martin Otter, DLR-SR (first version: Jan. 14, 2017)
#
# Details are described in the paper:
#   Otter, Elmqvist (2017): Transformation of Differential Algebraic Array Equations to
#                           Index One Form. Modelica'2017 Conference.
#
################################################

"""
    (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs)

Equations es shall be solved with respect to variables vs. The function returns
the teared equation so that if vTear is given, vSolved can be computed from eSolved
in a forward sequence (so solving eSolved[1] for vSolved[1], eSolved[2] for vSolved[2],
and so on). vTear must be selected, so that the equations eResidues are fulfilled.
Equations es are the union of eSolved and eResidue.
Variables vs are the union of vSolved and vTear.

Gsolvable defines the variables that can be explicitly solved in every equation without influencing the solution space
(= rank preserving operation).
"""
function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}, vs::Vector{Int})
    G = ict.graph
    vActive = BitSet(vs)

    for eq in es  # iterate only over equations that are not in eSolvedFixed
        for vj in Gsolvable[eq]
            if G.matching[vj] === unassigned && (vj in vActive)
                r = add_edge_checked!(ict, Iterators.filter(!=(vj), 𝑠neighbors(G.graph, eq)), vj) do G
                    G.matching[vj] = eq
                    G.ne += length(𝑠neighbors(G.graph, eq)) - 1
                end
                r && break
            end
        end
    end

    vSolved = filter(v->G.matching[v] !== unassigned, topological_sort(ict))
    inv_matching = Union{Missing, Int}[missing for _ = 1:nv(G)]
    for (v, eq) in pairs(G.matching)
        eq === unassigned && continue
        inv_matching[v] = eq
    end
    eSolved = getindex.(Ref(inv_matching), vSolved)
    vTear = setdiff(vs, vSolved)
    eResidue = setdiff(es, eSolved)
    return (eSolved, vSolved, eResidue, vTear)
end


text:
using DashBootstrapComponents

standalone_radio_check = html_div([
    html_div([
        dbc_checkbox(
            id = "standalone-checkbox",
            label = "This is a checkbox",
            value = false,
        ),
        dbc_switch(
            id = "standalone-switch",
            label = "This is a toggle switch",
            value = false,
        ),
        dbc_radiobutton(
            id = "standalone-radio",
            label = "This is a radio button",
            value = false,
        ),
    ]),
    html_p(id = "standalone-radio-check-output"),
]);

callback!(
    app,
    Output("standalone-radio-check-output", "children"),
    Input("standalone-checkbox", "value"),
    Input("standalone-switch", "value"),
    Input("standalone-radio", "value"),
) do checkbox_checked, switch_checked, radio_checked

    output_string =
        """Selections:  """ *
        """Checkbox: $(checkbox_checked ? "True" : "False"), """ *
        """Toggle Switch: $(switch_checked ? "True" : "False"), """ *
        """Radio Button: $(radio_checked ? "True" : "False")"""
    return output_string
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
export Config
using YAML

function Config(;filename="cfg/base.yaml")
    YAML.load_file(filename)
end

function Config(filename::String)
    YAML.load_file(filename)
end

function Config(filenames::Array{String})
    merge([YAML.load_file(f) for f in filenames]...)
end

function Config(cfg::Dict)
    cfg
end

function Config(cfg::Dict, other::Dict)
    merge(cfg, other)
end

function Config(cfg::Dict, other::String)
    merge(cfg, Config(other))
end

function Config(cfg::Dict, other::Array{String})
    merge(cfg, Config(other))
end

function input_lengths(cfg::Dict)
    inputs = Array{Int64}(undef, 10)
    inputs[1] = cfg["n_cell_params"]
    inputs[2] = cfg["n_cell_params"]
    inputs[3] = cfg["n_cell_params"]
    inputs[4] = cfg["n_cell_params"] + cfg["n_cell_state"] + 1
    inputs[5] = cfg["n_cell_params"] + cfg["n_cell_state"]
    inputs[6] = 2*cfg["n_cell_params"]
    inputs[7] = 2*cfg["n_cell_params"]
    inputs[8] = 2*cfg["n_cell_params"] + cfg["n_conn_params"] + cfg["n_conn_state"]
    inputs[9] = cfg["n_conn_params"] + cfg["n_conn_state"] + 2
    inputs[10] = 2*cfg["n_cell_params"] + cfg["n_conn_params"] + cfg["n_conn_state"]
    inputs
end

function output_lengths(cfg::Dict)
    outputs = Array{Int64}(undef, 10)
    outputs[1] = 1
    outputs[2] = cfg["n_cell_params"]
    outputs[3] = 1
    outputs[4] = cfg["n_cell_state"] + 1
    outputs[5] = cfg["n_cell_params"]
    outputs[6] = 1
    outputs[7] = cfg["n_conn_params"]
    outputs[8] = 1
    outputs[9] = cfg["n_conn_state"] + 1
    outputs[10] = cfg["n_conn_params"]
    outputs
end


text:
function evo_sweep(MPSvec::Vector{<:Any}, Cpre::AbstractArray{<:Number}, MPOvec::Vector{<:MPOsparseT}, RBlocks::Vector{<:Any}, dτ::Union{Float64,ComplexF64})

	#
	#	foward sweep
	#
	Lenv = [Array{ComplexF64,2}(I,1,1)]
	LBlocks = [Lenv]
	for site in 1:length(MPSvec)-1

		#
		# evolve A_site forward in time
		#

		localMPO = MPOvec[site]
		Renv = RBlocks[end-(site-1)]
		Lenv = LBlocks[site]
		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)
		AR_site = MPSvec[site]
		@tensor MPSvec[site][α, d, β] := Cpre[α, γ]*AR_site[γ, d, β]

		#MPS = Heff(MPSvec[site])
		#@tensor energyTest[] := MPS[a,b,c]*conj(MPSvec)[site][a,b,c]
		#@show energyTest

		ACnew, info = exponentiate(x->Heff(x), -0.5*dτ, MPSvec[site]; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

		#
		# split A_new_site and left can it
		#

		ALnew, S, Vdagger = leftCanSite(ACnew; optSVD = true)
		S = S./sign(S[1,1])
		snorm = sqrt(sum(S.^2)) #normalize S
		S = S./snorm
		C = S*Vdagger
		MPSvec[site] = ALnew

		#
		# build up Lenv
		#
		Lenv = applyTM_MPO([MPSvec[site]], [localMPO], LBlocks[end]; left=true)
		push!(LBlocks, Lenv) # (int,1,....,N-1)

		#
		# evolve C back in time
		#
		Hceff(x) = applyHCeff(x, Lenv, Renv)
		Cpre, info = exponentiate(x->Hceff(x), 0.5*dτ, C; ishermitian = true, tol = 10e-16)

	end
	#@show "at N-1"
	#
	#	evolve site N
	#
	Lenv=LBlocks[end]
	Renv=RBlocks[1]
	localMPO = MPOvec[end]
	Heff(x) = applyHeff(x, localMPO, Lenv, Renv)
	AR_site = MPSvec[end]
	@tensor v_0[α, d, β] := Cpre[α, γ]*AR_site[γ, d, β]

	ACnew, info = exponentiate(Heff, -1.0*dτ, v_0; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

	MPSvec[end] = ACnew

	RBlocks = [Renv]

	@show  "energy last site ",applyTM_MPO(MPSvec[1:end], MPOvec[1:end], [Array{ComplexF64,2}(I,1,1)]; left = false)

	#@show "at N"
	#
	#	evolve back
	#
	Cvec = []
	for site in length(MPSvec):-1:2

		#here MPS hast to be reshaped

		AR, S, U = rightCanSite(MPSvec[site];optSVD=true)
		MPSvec[site] = AR
		S = S./sign(S[1,1])
		snorm = sqrt(sum(S.^2))
		S = S./snorm
		C = U*S
		push!(Cvec, C)

		#
		# add RBlock
		#
		#
		Renv = applyTM_MPO([MPSvec[site]],[MPOvec[site]], RBlocks[end]; left = false)
		push!(RBlocks, Renv)#[int,N,.....,2]

		#
		#	evolve C backwards
		#
		Lenv = LBlocks[site]
		Hceff(x) = applyHCeff(x, Lenv, Renv)
		Cnext, info = exponentiate(Hceff, 0.5*dτ, C)

		#
		#	evolve AC forward
		#
		@tensor v_0[α, d, β] := MPSvec[site-1][α, d, γ]*Cnext[γ, β]
		Lenv = LBlocks[site-1]
		Renv = RBlocks[end]
		localMPO = MPOvec[site-1]
		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)

		ACnew, info = exponentiate(Heff, -0.5*dτ, v_0; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

		MPSvec[site-1] = ACnew

	end

	AR, S, U = rightCanSite(MPSvec[1]; optSVD = true)
	MPSvec[1] = AR
	S = S./sign(S[1,1])
	snorm = sqrt(sum(S.^2))
	S = S./snorm
	Cpre = U*S
	#@show Cpre

	return Cpre, RBlocks, Cvec[end:-1:1]
end

function evo_sweep_2Site(MPSvec::Vector{<:Any}, Cpre::AbstractArray{<:Number}, MPOvec::Vector{<:MPOsparseT}, RBlocks::Vector{<:Any}, dτ::Union{Float64,ComplexF64}; maxDim::Int = 200)

	#
	#	foward sweep
	#
	Lenv = [Array{ComplexF64,2}(I,1,1)]
	LBlocks = [Lenv]
	AC = MPSvec[1]
	for site in 2:length(MPSvec)

		#
		# 	evolve 2 sites forward
		#

		# 	get MPOs
		localMPO1 = MPOvec[site-1]
		localMPO2 = MPOvec[site]

		#	get Environment
		Renv = RBlocks[end-(site-1)]
		Lenv = LBlocks[site-1]
		# change to pop!
		#Renv = pop!(RBlocks)
		#Lenv = LBlocks[end]

		Heff2(x) = applyHeff2(x, localMPO1, localMPO2, Lenv, Renv)
		@tensor AC2[α, d1, d2, β] := AC[α, d1, γ]*MPSvec[site][γ, d2, β]
		AC2new, info = exponentiate(x->Heff2(x), -0.5*dτ, AC2; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

		#
		#	split AC2new
		#
		(a,d1, d2, b) = size(AC2new)
		AC2 = reshape(AC2new, a*d1, d2*b)
		U, S, V = svd(AC2)

		if length(S) > maxDim
			S = S[1:maxDim] # truncate
		end
		sdim = length(S)

		S = diagm(0=>S)
		S = S./sign(S[1,1])
		snorm = sqrt(sum(S.^2)) #normalize S
		S = S./snorm

		AL = reshape(U[:,1:sdim], a, d1, size(S)[1])
		MPSvec[site-1] = AL
		AC = reshape(S*(V')[1:sdim,:], size(S)[2], d2, b)

		#
		#	evolve AC backwards
		#

		#
		#	build up Lenv
		#
		Lenv = applyTM_MPO([MPSvec[site-1]], [localMPO1], LBlocks[end]; left=true)
		push!(LBlocks, Lenv) # (int,1,....,N-1)

		#
		#	evolve C back in time
		#
		Heff(x) = applyHeff(x, localMPO2, Lenv, Renv)
		AC, info = exponentiate(x->Heff(x), 0.5*dτ, AC; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

	end

	#pop!(LBlokcs)
	Renv=RBlocks[1]
	RBlocks = [Renv]
	#
	#	evolve back
	#
	Cvec = []
	for site in length(MPSvec)-1:-1:1

		# 	get MPOs
		localMPO1 = MPOvec[site]
		localMPO2 = MPOvec[site+1]

		#	get Environment
		Renv = RBlocks[end]
		Lenv = LBlocks[site]
		# with pop!
		#Lenv = pop!(LBlocks)
		#Renv = RBlocks[end]

		#
		#	evolve AC forward
		#
		Heff2(x) = applyHeff2(x, localMPO1, localMPO2, Lenv, Renv)
		@tensor AC2[α, d1, d2, β] := AC[γ, d2, β]*MPSvec[site][α, d1, γ]
		AC2new, info = exponentiate(x->Heff2(x), -0.5*dτ, AC2; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

		#
		#	split AC2new
		#
		(a,d1, d2, b) = size(AC2new)
		AC2 = reshape(AC2new, a*d1, d2*b)
		U, S, V = svd(AC2)

		if length(S) > maxDim
			S = S[1:maxDim] # truncate
		end
		sdim = length(S)

		S = diagm(0=>S)
		S = S./sign(S[1,1])
		snorm = sqrt(sum(S.^2)) #normalize S
		S = S./snorm

		AC = reshape(U[:,1:sdim]*S, a, d1, size(S)[1])
		AR = reshape((V')[1:sdim, :], size(S)[2], d2, b)
		MPSvec[site+1] = AR

		#
		# add RBlock
		#
		#
		Renv = applyTM_MPO([MPSvec[site+1]],[localMPO2], RBlocks[end]; left = false)
		push!(RBlocks, Renv)#[int,N,.....,2]

		#
		#	evolve C backwards
		#
		Lenv = LBlocks[site]
		Heff(x) = applyHeff(x, localMPO1, Lenv, Renv)
		AC, info = exponentiate(x->Heff(x), 0.5*dτ, AC; ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)

	end

	MPSvec[1] = AC
	AR, S, U = rightCanSite(MPSvec[1]; optSVD = true)
	MPSvec[1] = AR
	S = S./sign(S[1,1])
	snorm = sqrt(sum(S.^2))
	S = S./snorm
	Cpre = U*S
	#@show Cpre

	return Cpre, RBlocks

end

function vmps_sweep(MPS::Vector{<:Any}, Cpre::AbstractArray, MPOvec::Vector{<:MPOsparseT}, RBlocks::Vector{<:Any})

	#
	#	foward sweep
	#
	L = Array{ComplexF64}(I,1,1)
	Lenv = [L]
	LBlocks = [Lenv]
	for site in 1:length(MPS)-1

		#
		# evolve A_site forward in time
		#

		localMPO = MPOvec[site]
		Renv = RBlocks[end-(site-1)]
		Lenv = LBlocks[site]
		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)
		AR_site = MPS[site]
		@tensor MPS[site][α, d, β] := Cpre[α, γ]*AR_site[γ, d, β]

		E, ACnew, info = eigsolve(x->Heff(x), MPS[site], 1, :SR, ishermitian = true, tol= 10e-14)
		@show E
		ACnew = ACnew[1]
		#@show info

		#
		# split A_new_site and left can it
		#

		ALnew, S, Vdagger = leftCanSite(ACnew; optSVD = true)

		S = S/sign(S[1,1])
		S = S/(sqrt(dot(S,S))) #normalize S
		Cpre = S*Vdagger

		MPS[site] = ALnew

		#
		# build up Lenv
		#
		Lenv = applyTM_MPO([MPS[site]], [localMPO], LBlocks[end]; left=true)
		push!(LBlocks, Lenv) # (int,1,....,N-1)

	end
	#@show "at N-1"
	#
	#	evolve site N
	#
	Lenv=LBlocks[end]
	Renv=RBlocks[1]
	localMPO = MPOvec[end]
	Heff(x) = applyHeff(x, localMPO, Lenv, Renv)
	AR_site = MPS[end]

	@tensor v_0[α, d, β] := Cpre[α, γ]*AR_site[γ, d, β]

	#ACnew, info = exponentiate(Heff, -1.0*dτ, SparseTensor(v_0); ishermitian = true, tol=10e-16, verbosity = 0, maxiter = 20)
	E, ACnew, info = eigsolve(x->Heff(x), v_0, 1, :SR, ishermitian = true, tol= 10e-14)
	@show E

	MPS[end] = ACnew[1]
	RBlocks = [Renv]

	#
	#	evolve back
	#
	Cvec = []
	for site in length(MPS):-1:2

		#here MPS hast to be reshaped

		AR, S, U = rightCanSite(MPS[site]; optSVD=true)
		MPS[site] = AR
		S = S/sign(S[1,1])
		snorm = sqrt(dot(S,S))
		S = S/snorm
		Cnext = U*S
		push!(Cvec, Cnext)

		#
		# add RBlock
		#
		Renv = applyTM_MPO([MPS[site]], [MPOvec[site]], RBlocks[end]; left = false)
		push!(RBlocks, Renv)#[int,N,.....,2]

		#
		#	evolve AC forward
		#
		@tensor v_0[α, d, β] := MPS[site-1][α, d, γ]*Cnext[γ, β]
		Lenv = LBlocks[site-1]
		Renv = RBlocks[end]
		localMPO = MPOvec[site-1]
		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)

		E, ACnew, info = eigsolve(x->Heff(x), v_0, 1, :SR, ishermitian = true, tol= 10e-14)
		@show E
		MPS[site-1] = ACnew[1]

	end

	AR, S, U = rightCanSite(MPS[1]; optSVD = true)
	MPS[1] = AR

	S = S/sign(S[1,1])
	snorm = sqrt(dot(S,S))
	S = S/snorm
	Cpre = U*S

	return Cpre, RBlocks, Cvec[end:-1:1]
end

function applyMPO_sweep(MPSvec, new_MPSvec, MPOvec, RBlocks)
	Lenv = [Array{ComplexF64,2}(I,1,1)]
	LBlocks = [Lenv]
	Cpre = Array{ComplexF64,2}(I,1,1)
	N = size(MPSvec)[1]

	for site in 1:N

		# multiply B(site) with Cpre to form M
		localMPO = MPOvec[site]
		Renv = RBlocks[end-(site-1)]
		Lenv = LBlocks[site]

		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)
		AR_site = MPSvec[site]
		@tensor MPSvec[site][α, d, β] := Cpre[α, γ]*AR_site[γ, d, β]

		#
		#	apply Heff
		#
		new_M = Heff(MPSvec[site])
		#move center
		A, S, V = leftCanSite(new_M, optSVD=true)

		# update site
		new_MPSvec[site] = A

		A, S, V = leftCanSite(MPSvec[site], optSVD=true)
		S = S./sign(S[1,1])
		S = S./(sqrt(sum(S.^2))) #normalize S
		Cpre = S*V

		# update site
		MPSvec[site] = A

		if site < N
			#
			# build up Lenv
			#
			Lenv = applyTM_MPO([MPSvec[site]],[new_MPSvec[site]], [localMPO], LBlocks[end]; left=true)
			push!(LBlocks, Lenv) # (int,1,....,N-1)
		end
	end

	Renv = [Array{ComplexF64,2}(I,1,1)]
	RBlocks = [Renv]

	Cnext = Cpre #Array{ComplexF64, 2}([[1]])
	for site in N:-1:1

		AL_site = MPSvec[site]
		@tensor MPSvec[site][α, d, β] := AL_site[α, d, γ]*Cnext[γ, β]
		Lenv = LBlocks[site]
		Renv = RBlocks[end]

		localMPO = MPOvec[site]
		bDim = typeof(localMPO.bDim) == Int ? (localMPO.bDim, localMPO.bDim) : localMPO.bDim
		Heff(x) = applyHeff(x, localMPO, Lenv, Renv)

		# multiply A(site) with Cpre to form M
		new_M = Heff(MPSvec[site])

		#move center
		B, __, __ = rightCanSite(new_M, optSVD=true)
		#B, S, U = rightCanSite(new_M, optSVD=true)
		#S = S./sign(S[1,1])
		#snorm = sqrt(sum(S.^2))
		#S = S./snorm
		#Cnext = U*S

		# update site
		new_MPSvec[site] = B

		B, S, U = rightCanSite(MPSvec[site], optSVD=true)
		S = S./sign(S[1,1])
		snorm = sqrt(sum(S.^2))
		S = S./snorm
		Cnext = U*S

		MPSvec[site] = B

		# update right environment
		if site > 1
			Renv = applyTM_MPO([MPSvec[site]], [new_MPSvec[site]], [MPOvec[site]], RBlocks[end]; left = false)
			push!(RBlocks, Renv)#[int,N,.....,2]
		end
	end

	return new_MPSvec, RBlocks
end

function iter_applyMPO(MPSvec::Vector{<:Any}, MPOvec::Vector{<:MPOsparseT}, new_MaxD::Int; Niter = 1)
	N = size(MPSvec)[1]
	d = size(MPSvec[1])[2]

	new_MPSvec, Cvec = randMPS(N, d, new_MaxD)

	# bring both MPS in right can form
	rightCanMPS(MPSvec)
	rightCanMPS(new_MPSvec)

	Renv = [Array{ComplexF64,2}(I,1,1)]
	RBlocks = [Renv]
	for site = 1:length(MPSvec)-1
		Renv = applyTM_MPO([MPSvec[end-(site-1)]], [new_MPSvec[end-(site-1)]], [MPOvec[end-(site-1)]], RBlocks[site]; left = false)
		push!(RBlocks, Renv)#[init,N,.....,2]
	end

	# sweep through the system and minimize || |ψ_old> - |ψ_new> ||^2
	for i in 1:Niter
		new_MPSvec, RBlocks = applyMPO_sweep(MPSvec, new_MPSvec, MPOvec, RBlocks)
	end

	return new_MPSvec
end

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
export ConcatCube, concatenateCubes
export mergeAxes
import ..Cubes: YAXArray, caxes, iscompressed, cubechunks, chunkoffset
using DiskArrayTools: diskstack, DiskArrayTools

function Base.map(op, incubes::YAXArray...)
    axlist = copy(caxes(incubes[1]))
    all(i -> caxes(i) == axlist, incubes) || error("All axes must match")
    props = merge(getattributes.(incubes)...)
    YAXArray(
        axlist,
        broadcast(op, map(c -> getdata(c), incubes)...),
        props,
        mapreduce(i -> i.cleaner, append!, incubes),
    )
end

"""
    function concatenateCubes(cubelist, cataxis::CategoricalAxis)

Concatenates a vector of datacubes that have identical axes to a new single cube along the new
axis `cataxis`
"""
function concatenatecubes(cl, cataxis::CubeAxis)
    length(cataxis.values) == length(cl) ||
        error("cataxis must have same length as cube list")
    axlist = axcopy.(caxes(cl[1]))
    T = eltype(cl[1])
    N = ndims(cl[1])
    cleaners = CleanMe[]
    append!(cleaners, cl[1].cleaner)
    for i = 2:length(cl)
        all(caxes(cl[i]) .== axlist) ||
            error("All cubes must have the same axes, cube number $i does not match")
        eltype(cl[i]) == T || error(
            "All cubes must have the same element type, cube number $i does not match",
        )
        ndims(cl[i]) == N || error("All cubes must have the same dimension")
        append!(cleaners, cl[i].cleaner)
    end
    props = mapreduce(getattributes, merge, cl, init = getattributes(cl[1]))
    YAXArray([axlist..., cataxis], diskstack([getdata(c) for c in cl]), props, cleaners)
end
function concatenatecubes(; kwargs...)
    cubenames = String[]
    for (n, c) in kwargs
        push!(cubenames, string(n))
    end
    cubes = map(i -> i[2], collect(kwargs))
    findAxis("Variable", cubes[1]) === nothing ||
        error("Input cubes must not contain a variable kwarg concatenation")
    concatenateCubes(cubes, CategoricalAxis("Variable", cubenames))
end

max_forks_repo_name: meggart/YAXArrays


text:
######
# This file is part of the BivariateCopulas.jl package.
#
#   Module definition
#
#           University of Liverpool, Institute for Risk and Uncertainty
#
#                                           Authors: Ander Gray & Liam T. Henry
#                                           Email:  ander.gray@liverpool.ac.uk
######

module BivariateCopulas

using Reexport

@reexport using Distributions             #   https://github.com/JuliaStats/Distributions.jl
using LinearAlgebra                       #   Main Library
using Interpolations                      #   https://github.com/JuliaMath/Interpolations.jl

import Distributions: cdf, quantile, rand

import Base: rand

global  n  = 200;                       # Number of descretizations
global bOt = 0.0001;                    # smallest quantile to use if left tail is unbounded. Should we use these to define the copula ranges?
global tOp = 0.9999;                    # largest quantile to use if right tail is unbounded
#global bOt = 0;
#global tOp = 1;

abstract type AbstractCopula <: Real end
abstract type AbstractJoint <: Real end
abstract type AbstractMarginal <: ContinuousUnivariateDistribution end

include("copulas.jl")
include("NormalDistribution.jl")

export
    copula, density, cdf, sample, rand, CholeskyGaussian,
    conditional, calcCdf, calcDensity, M, W, Pi, Frank, Clayton, Gaussian,
    norm_cdf, bivariate_cdf, mvNormCdf, Joint, invCdf, marginal, quantile,
    ρCopula, SpearmanCopula, KendalCopula, τCopula

end #BivariateCopulas


text:
using Test
using Statistics
using DelimitedFiles
using Dates

@testset "nse" begin
        for i in 1: 1000
                Observed_Discharge = rand(5:0.5:80, 1000)
                Modelled_Discharge = rand(5:0.5:80, 1000)
                Mean_Observed_Discharge = ones(1000) * mean(Observed_Discharge)
                Nominator = sum((Modelled_Discharge - Observed_Discharge).^2)
                Denominator = sum((Observed_Discharge - Mean_Observed_Discharge).^2)
                NSE = 1 - Nominator / Denominator
                @test NSE == nse(Observed_Discharge, Modelled_Discharge)
                # observed is the same as modelled
                Modelled_Discharge = Observed_Discharge
                @test nse(Observed_Discharge, Modelled_Discharge) == 1
        end
end

@testset "nselog" begin
        for i in 1: 1000
                Observed_Discharge = rand(5:0.5:80, 1000)
                Modelled_Discharge = rand(5:0.5:80, 1000)
                Mean_Observed_Discharge = ones(1000) * mean(Observed_Discharge)
                Nominator = sum((log.(Modelled_Discharge) - log.(Observed_Discharge)).^2)
                Denominator = sum((log.(Observed_Discharge) - log.(Mean_Observed_Discharge)).^2)
                NSE = 1 - Nominator / Denominator
                @test NSE == lognse(Observed_Discharge, Modelled_Discharge)
                # observed is the same as modelled
                Modelled_Discharge = Observed_Discharge
                @test lognse(Observed_Discharge, Modelled_Discharge) == 1
        end
end

@testset "VE" begin
        for i in 1: 10000
                Observed_Discharge = rand(5:0.5:80, 1000)
                Modelled_Discharge = rand(5:0.5:80, 1000)
                Nominator = sum(abs.(Modelled_Discharge - Observed_Discharge))
                Denominator = sum(Observed_Discharge)
                VE = 1 - Nominator / Denominator
                @test VE == volumetricefficiency(Observed_Discharge, Modelled_Discharge)
                # observed is the same as modelled
                Modelled_Discharge = Observed_Discharge
                @test volumetricefficiency(Observed_Discharge, Modelled_Discharge) == 1
        end
end

@testset "FDC" begin
        for i in 1:1000
                Discharge = rand(5:0.5:80, 1000)
                SortedQ = sort(Discharge, rev = true)
                Rank = collect(1 : length(Discharge))
                @test minimum(Discharge) == SortedQ[end]
                @test maximum(Discharge) == SortedQ[1]
                Exceedance = Rank / (length(Discharge) + 1)
                SortedQ_Function, Exceedance_Function = flowdurationcurve(Discharge)
                @test Exceedance == Exceedance_Function
                @test SortedQ == SortedQ_Function
        end
end

@testset "autocorrelation" begin
        for i in 1:1000
                Discharge = rand(5:0.5:80, 1000)
                Timelag = 1
                Discharge_1 = Discharge[1 + Timelag: end]
                Discharge_0 = Discharge[1: end - Timelag]
                Mean_Discharge = ones(length(Discharge_0)) * mean(Discharge_0)
                Nominator = sum((Discharge_0 - Mean_Discharge) .* (Discharge_1 - Mean_Discharge))
                Denominator = sum((Discharge_0 - Mean_Discharge).^2)
                Correlation = Nominator / Denominator
                @test round(Correlation, digits=12) == round(autocorrelation(Discharge, Timelag), digits=12)
                #@test round(Correlation, digits=12) == round(autocorrelation2(Discharge, Timelag), digits = 12)
                # if timelag is zero correlation should be one
                Timelag = 0
                @test 1 == autocorrelation(Discharge, Timelag)
                #@test 1 == autocorrelation2(Discharge, Timelag)
        end
end

@testset "autocorrelationcurve" begin
        for i in 1: 1000
                Discharge = rand(5:0.5:80, 1000)
                Timelag = 30
                Discharge_0 = Discharge[1 : end - Timelag]
                Discharge_0_Shift = Discharge[1 + Timelag : end]
                Mean_Discharge = ones(length(Discharge_0)) * mean(Discharge_0)
                AC = zeros(Timelag + 1)
                AC_cor = zeros(Timelag + 1)
                count = 0
                count_large = 0
                for i in 0 : Timelag
                        Discharge_1 = Discharge[1 + i : end - Timelag + i]
                        Discharge_1_shift = Discharge[1 + Timelag - i : end - i]
                        Nominator = sum((Discharge_0 - Mean_Discharge) .* (Discharge_1 - Mean_Discharge))
                        Denominator = sum((Discharge_0 - Mean_Discharge).^2)
                        AC[i + 1] = Nominator / Denominator
                        @test AC[1] >= AC[i + 1]
                        if i > 0 && AC[i + 1] <= AC[i]
                                count += 1
                        elseif i > 0 && AC[i + 1] > AC[i]
                                count_large += 1
                        end
                        # for other function
                        AC_cor[i + 1] = cor(Discharge_0_Shift, Discharge_1_shift)
                end
                Lags = collect(0:Timelag)
                #@test count >= count_large
                AC_Function, Lags_Function = autocorrelationcurve(Discharge, Timelag)
                @test AC == AC_Function
                @test Lags == Lags_Function
                #test other function
                # AC_Function, Lags_Function = autocorrelationcurve(Discharge, Timelag)
                # @test AC_cor == AC_Function
                # @test Lags == Lags_Function
        end
end

@testset "runoff" begin
        for i in 1 :100
                Area = rand(50.:5.:500.)[1]
                Timespan = 1096
                Precipitation = rand(5:0.5:80, Timespan)
                Discharge = rand(5:0.5:80, Timespan)
                Timeseries = readdlm("Defreggental/tas_model_timeseries.txt")[1 : Timespan]
                Timeseries_DateFormat = Date.(Timeseries, dateformat"y,m,d")
                Runoff_1 = sum(Discharge[1:31]) / Area * (1000 * 3600 * 24) / sum(Precipitation[1:31])
                Runoff_2 = sum(Discharge[32:59]) / Area * (1000 * 3600 * 24) / sum(Precipitation[32:59])
                Runoff_Function, Month = monthlyrunoff(Area, Precipitation, Discharge,Timeseries_DateFormat)
                @test round(Runoff_Function[1], digits= 6) == round(Runoff_1, digits = 6)
                @test round(Runoff_Function[2], digits= 6) == round(Runoff_2, digits = 6)

                Runoff_mean_yearly = averagemonthlyrunoff(Area, Precipitation, Discharge,Timeseries_DateFormat)
                Runoff_51 = sum(Discharge[366: 366+30]) / Area * (1000 * 3600 * 24) / sum(Precipitation[366 : 30 + 366])
                Runoff_52 = sum(Discharge[731: 731+30]) / Area * (1000 * 3600 * 24) / sum(Precipitation[731: 30 + 731])
                meanRunoff = mean([Runoff_1 Runoff_51 Runoff_52])
                @test round(meanRunoff, digits=5) == round(Runoff_mean_yearly[1], digits=5)

                # if Precipitation is equal to Discharge than Runoff = 1
                Timespan = 31
                Precipitation = rand(5:0.5:80, Timespan)
                Discharge = rand(5:0.5:80, Timespan)
                Timeseries = readdlm("Defreggental/tas_model_timeseries.txt")[1 : Timespan]
                Timeseries_DateFormat = Date.(Timeseries, dateformat"y,m,d")
                Discharge = Precipitation * Area / (1000 * 3600 * 24)
                Runoff_Function, Month = monthlyrunoff(Area, Precipitation, Discharge,Timeseries_DateFormat)
                @test Runoff_Function[1] == 1.0

                # calculate yearly mean runoff
        end

end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
function jlpeval(sexp)
    ts = typof(sexp)

    if ts <: JLPType
    elseif ts <: JLPAbstList
        jlpasbtlisteval(sexp)
    end
end

reservedwords = Dict("true" => (x -> JLPBool(parse(Bool, x))),
                     "false" => (x -> JLPBool(parse(Bool, x))))

@enum Token begin
    jlpListStart
    jlpListEnd
    jlpInt
    jlpFloat
    jlpSymbol
    jlpString
    jlpChar
end

struct Tokenval
    token::Token
    val::String
end

function isliststarted(tkn)
    isnothing(tkn) || tkn == jlpListStart
end

function tokenize(code::String)
    code = code
    state = Union{Nothing, Token}[jlpListStart]
    result = Vector[[]]
    tmpstr   = Char[]

    tmpreset!() = tmpstr = Char[]

    issttstarted() = !isempty(state) && isliststarted(state[end])

    addtmp!(c::Char) = push!(tmpstr, c)

    function addtmp2res!()
        push!(result[end], join(tmpstr))
        tmpreset!()
    end

    function addval2res!()
        if !isempty(state) && !issttstarted()
            st = pop!(state)
            ts = join(tmpstr)
            tmpreset!()
            an = if st == jlpInt
                     JLPInteger(parse(Int64, ts))
                 elseif st == jlpFloat
                     JLPFloat(parse(Float64, ts))
                 elseif st == jlpSymbol
                     get(reservedwords, ts, JLPSymbol)(ts)
                 end
            push!(result[end], an)
        end
    end

    for c in code
        if !isempty(state) && state[end] == jlpString
            if (length(tmpstr) != 0 && tmpstr[end] != '\\' && c == '"') || (length(tmpstr) == 0 && c == '"')
                res = pop!(result[end])
                push!(result[end], JLPString(JLPSymbol(res.val), join(tmpstr)))
                tmpreset!()
                pop!(state)
                continue
            end
            addtmp!(c)
            continue
        elseif c == ' ' || c == '\n'
            if !issttstarted()
                addval2res!()
                push!(state, nothing)
            end
            continue
        else
            if length(state) != 0 && isnothing(state[end])
                pop!(state)
            end
        end

        if '0' <= c <= '9'
            if issttstarted()
                push!(state, jlpInt)
            elseif length(tmpstr) == 1 && tmpstr[end] == '-'
                pop!(state)
                push!(state, jlpInt)
            end
            addtmp!(c)

        elseif c == '.'
            if state[end] == jlpInt
                pop!(state)
                push!(state, jlpFloat)
            end
            addtmp!(c)

        elseif c == '('
            if !isempty(state) && !issttstarted()
                pop!(state)
            end
            push!(state, jlpListStart)
            push!(result, [])
            push!(result[end], Tokenval(jlpListStart, join(tmpstr)))
            tmpreset!()

        elseif c == ')'
            if !issttstarted()
                addval2res!()
            end
            bef = JLPNil()
            while length(result[end]) != 1
                tmp = pop!(result[end])
                bef = JLPList(!(typeof(tmp) <: Tokenval) || (typeof(tmp) <: Tokenval && tmp.val == "") ?
                              JLPNil() : JLPSymbol(tmp.val), tmp, bef)
            end
            tmp = pop!(result[end])
            pop!(result)
            push!(result[end], JLPList(!(typeof(tmp) <: Tokenval) || (typeof(tmp) <: Tokenval && tmp.val == "") ?
                                       JLPNil() : JLPSymbol(tmp.val), typeof(bef) <: JLPNil ? bef : bef.car, typeof(bef) <: JLPNil ? bef : bef.cdr))
            pop!(state)

        elseif c == '"'
            if !isempty(state) && !issttstarted()
                pop!(state)
            end
            push!(state, jlpString)
            push!(result[end], Tokenval(jlpString, join(tmpstr)))
            tmpreset!()
        else
            if issttstarted()
                push!(state, jlpSymbol)
            end
            addtmp!(c)
        end
    end
    if !issttstarted()
        addval2res!()
        push!(state, nothing)
    end
    result[end]
end


text:
mat(ex::SubHam{<:Tuple, <:LocalOperator}) = ex.val

# convert config to indices
Base.to_indices(::Type{L}, inds...) where L <: SiteLabel = map(x->to_indices(L, x), inds)

Base.to_indices(::Type{Bit}, val::Number) = Int(val + 1)
Base.to_indices(::Type{Spin}, val::Number) = Int(div(val + 1, 2))
Base.to_indices(::Type{Half}, val::Number) = Int(val + 0.5)

@generated function Base.to_indices(::Type{Bit}, val::NTuple{N, T}) where {N, T}
    ex = :(val[1])
    for i = 2:N
        factor = 2^(i-1)
        ex = :(val[$i] * $factor + $ex)
    end
    :(ex)
end

@generated function Base.to_indices(::Type{Spin}, val::NTuple{N, T}) where {N, T}
    ex = :(div(val[1] + 1, 2))
    for i = 2:N
        factor = 2^(i-1)
        ex = :(div(val[$i] + 1, 2) * $factor + $ex)
    end
    :(Int($ex))
end

@generated function Base.to_indices(::Type{Half}, val::NTuple{N, T}) where {N, T}
    ex = :(val[1] + 0.5)
    for i = 2:L
        factor = 2^(i-1)
        ex = :((val[$i] + 0.5) * $factor + $ex)
    end
    :(Int($ex))
end

struct NumberOfBody{K} end
NumberOfBody(::NumberOfBody{K}) where K = NumberOfBody{K}()
NumberOfBody(::LocalOperator{T, N}) where {T, N} = NumberOfBody{N}()
NumberOfBody(x::SubHam) = NumberOfBody{length(x.placeholders)}()
NumberOfBody(x::Operator) = NumberOfBody(x.expr)

struct PlaceHolderType{T} end
PlaceHolderType(::SubHam{PT}) where PT = PlaceHolderType{PT}()

# NOTE: only operator acts on a single SubHam
PlaceHolderType(x::Operator) = PlaceHolderType(x.expr)
# PlaceHolderType(x::AddOrSub) =

"""
    LOpEigen{K}

Number of eigen states in local operator.
"""
struct LOpEigen{K} end

LOpEigen(::LocalOperator{T, N, O}) where {T, N, O} = LOpEigen{O}()
LOpEigen(x::Operator) = LOpEigen(x.expr)
LOpEigen(x::AddOrSub) = (LOpEigen(x.lhs); LOpEigen(x.rhs))
LOpEigen(x::SubHam{<:Tuple, <:LocalOperator}) = LOpEigen(x.val)

function LOpEigen(x::SubHam{<:Tuple, <:Tuple})
    n = LOpEigen(first(x.val))

    # check if all the local operator has the same number of eigens
    for each in x.val
        n === LOpEigen(each) || error("number of eigen states in given local operators are not the same")
    end
    n
end

struct LocalOperatorValType{T, OpType} end
LocalOperatorValType(::MT) where {T, MT <: AbstractMatrix{T}} = LocalOperatorValType{T, MT}()
LocalOperatorValType(x::LocalOperator) = LocalOperatorValType(x.val)
LocalOperatorValType(x::SubHam{<:Tuple, <:LocalOperator}) = LocalOperatorValType(x.val)
LocalOperatorValType(x::Operator) = LocalOperatorValType(x.expr)

# multiple operators
LocalOperatorValType(x::SubHam{<:Tuple, <:Tuple}) = LocalOperatorValType(x.val...)
LocalOperatorValType(a::SubHam, x::SubHam...) = LocalOperatorValType(LocalOperatorValType(a), LocalOperatorValType(x...))

LocalOperatorValType(x::LocalOperatorValType{T}...) where T = LocalOperatorValType{T, Tuple{_op_type(x...)...}}()
_op_type(::Tuple{}) = ()
_op_type(::LocalOperatorValType{T, Op}) where {T, Op} = Op
_op_type(::LocalOperatorValType{T, Op1}, ::LocalOperatorValType{T, Op2}, ops...) where {T, Op1, Op2} = (Op1, Op2, _op_type(ops)...)

function is_other_sites_the_same(::NumberOfBody{1}, site, region, lhs::Sites, rhs::Sites)
    @inbounds for other_site in Lattices.sites(region)
        if other_site != site
            s = region[other_site]
            lhs[s] == rhs[s] || return false
        end
    end
    true
end

function is_other_sites_the_same(::NumberOfBody, sites, region, lhs::Sites, rhs::Sites)
    for other_site in Lattices.sites(region)
        # NOTE: we seperate single site because of the following:
        #
        # # chain lattice
        # julia> 1 in 1
        #        true
        #
        # # square lattice
        # julia> (1, 1) in (1, 1)
        #        false
        if !(other_site in sites)
            s = region[other_site]
            lhs[s] == rhs[s] || return false
        end
    end
    true
end

Base.getindex(ex::OnRegion, lhs::Sites, rhs::Sites) = eval_expr(ex, lhs, rhs)

## eval Operator expressions
# unzip OnRegion
eval_expr(ex::OnRegion, lhs::Sites, rhs::Sites) =
    eval_expr(ex.expr, ex.region, lhs, rhs)

# 1. dispatch to specific local operator
eval_expr(ex::Operator, region, lhs::Sites, rhs::Sites) =
    eval_expr(ex, region, LocalOperatorValType(ex), PlaceHolderType(ex), lhs, rhs)

# 2. if there is no specific method for the operator, then fall back to general methods
eval_expr(ex::Operator, region, ::LocalOperatorValType, ::PlaceHolderType, lhs, rhs) =
    eval_expr(ex, region, LOpEigen(ex), NumberOfBody(ex), PlaceHolderType(ex), lhs, rhs)

# H1 + H2
eval_expr(ex::Add, region, lhs::Sites, rhs::Sites) =
    eval_expr(ex.lhs, region, lhs, rhs) + eval_expr(ex.rhs, region, lhs, rhs)
# H1 - H2
eval_expr(ex::Sub, region, lhs::Sites, rhs::Sites) =
    eval_expr(ex.lhs, region, lhs, rhs) - eval_expr(ex.rhs, region, lhs, rhs)

# α * H
eval_expr(ex::ScalarProd, region, lhs::Sites, rhs::Sites) =
    ex.val * eval_expr(ex.expr, region, lhs, rhs)

# sum(A[i]) on vertex
# NOTE:
# The summation looks like
#
# U x I x I + I x U x I + I x I x U
#
# which means the 2x2 local operator U
# defines which config should be the same/or not
#
# So we will check if position (excludes U) are the same in each loop
# this determine if this block is zero.
# TODO:
# 1. for diagnol local operators, we can directly check if all configurations
# are the same.
# 2. for local operators like X, we can directly check if all configurations
# are opposite (regarding to LOpEigen{2} operators)
function eval_expr(
    expr::Sum{<:SubHam}, region,
    ::LOpEigen{2}, n::NumberOfBody{1}, # 2x2 local operator
    ::PlaceHolderType{Tuple{Vertex}},
    lhs::Sites{L}, rhs::Sites{L}) where L

    Aij = zero(eltype(expr))

    @inbounds for site in Lattices.sites(region)
        s = region[site]
        lhs_ind = to_indices(L, lhs[s])
        rhs_ind = to_indices(L, rhs[s])

        # check identities
        Aij += is_other_sites_the_same(n, site, region, lhs, rhs) ? expr.expr[lhs_ind, rhs_ind] : 0
    end
    Aij
end

######################### specialized methods ##################################
# # sum(X[i]) over all vertex
# function eval_expr(
#     ex::Sum{<:SubHam}, region,
#     ::LocalOperatorValType{T, <:Pauli.X},
#     ::PlaceHolderType{Tuple{Vertex}},
#     lhs::Sites{L}, rhs::Sites{L}) where {T, L}
#
#     # check if all the configurations are opposite
#     @inbounds for site in Lattices.sites(region)
#         s = region[site]
#         # NOTE: since this is for Pauli.X
#         # we can assume there is only two
#         # states, so just use !=
#         lhs[s] != rhs[s] || return zero(T)
#     end
#
#     # T(length(lhs))
#     # TODO:
#     # NOT DONE YET
# end

# sum(Z[i]) over all vertex
function eval_expr(
    ex::Sum{<:SubHam}, region,
    ::LocalOperatorValType{T, <:Pauli.Z},
    ::PlaceHolderType{Tuple{Vertex}},
    lhs::Sites{L}, rhs::Sites{L}) where {T, L}

    # check if all the configurations are the same
    @inbounds for site in Lattices.sites(region)
        s = region[site]
        lhs[s] == rhs[s] || return zero(T)
    end

    # down down -> -1
    # up   up   -> +1
    downs, ups = n_downs_ups(lhs)
    T(ups - downs)
end

# sum(Z[i]Z[j]) on neighbors
function eval_expr(
    expr::Sum{<:SubHam}, region,
    ::LOpEigen{2}, n::NumberOfBody{2}, # 2x2 local operator
    ::PlaceHolderType{Tuple{Neighbor{K}, Neighbor{K}}},
    lhs::Sites{L}, rhs::Sites{L}) where {K, L}

    Aij = zero(eltype(expr))
    @inbounds for (site_a, site_b) in Lattices.edges(region, length=K)
        s_a = region[site_a]
        s_b = region[site_b]

        lhs_ind = to_indices(L, lhs[s_a], lhs[s_b])
        rhs_ind = to_indices(L, rhs[s_a], rhs[s_b])

        Aij += is_other_sites_the_same(n, (site_a, site_b), region, lhs, rhs) ? expr.expr[lhs_ind, rhs_ind] : 0
    end
    Aij
end

# partial indexing
# TODO: this should use partial eval
# TODO: add assertion
Base.getindex(ex::OnRegion, lhs::Sites, colon::Colon) = eval_expr(ex.expr, ex.region, lhs, colon)
Base.getindex(ex::OnRegion, colon::Colon, rhs::Sites) = eval_expr(ex.expr, ex.region, colon, rhs)

function eval_expr(
    ex, region,
    lhs::Sites{L, T, N}, ::Colon) where {L, T, N}

    # TODO:
    # wrap this as an iterator?
    # for rhs in IterName
    # blabla
    # end

    nelem = 1<<length(region)
    Ai = Vector{eltype(ex)}(undef, nelem)
    rhs = Sites{L, T, N}(fill(down(L), size(lhs)))
    for i = 1:nelem
        Ai[i] = eval_expr(ex, region, lhs, rhs)
        rhs << 1
    end

    Ai
end

function eval_expr(
    ex, region,
    ::Colon, rhs::Sites{L, T, N}) where {L, T, N}

    # TODO:
    # wrap this as an iterator?
    # for rhs in IterName
    # blabla
    # end

    nelem = 1<<length(region)
    Ai = Vector{eltype(ex)}(undef, nelem)
    lhs = Sites{L, T, N}(fill(down(L), size(rhs)))
    for i = 1:nelem
        Ai[i] = eval_expr(ex, region, lhs, rhs)
        lhs << 1
    end

    Ai
end

################################################################################

# TODO:
# sum(A1[i]A2[j]...) on vertex
LocalOperator(x::SubHam{<:Tuple, <:LocalOperator}) = x.val
LocalOperator{T, 1}(x::SubHam{<:Tuple, <:LocalOperator{T, 1}}) where T = x.val
LocalOperator{T, 1}(x::SubHam{<:Tuple, <:Tuple}) where T = prod(LocalOperator{T, 1}, x.val)

# TODO: use a repeated product instead
# 1. A1[i]A2[j] does not go through (1, 1), (2, 2), ..., this should use @subham A1[i]A2[i]
# 2. parse A1[i]A2[i] => (A1 * A2)[i] in make_subham, not in SubHamExpr, we keep A1[i]A2[i] for printing

function eval_expr(ex::Sum{<:SubHam}, region,
    ::LOpEigen{2}, n::NumberOfBody{N}, # 2^Nx2^N local operator
    ::PlaceHolderType{NTuple{N, Vertex}},
    lhs::Sites{L}, rhs::Sites{L}) where {N, L}

    error("Summation over multiple vertex is not supported yet. See issue #1")
    Aij = zero(eltype(ex))

    it = Lattices.sites(region)

    it = Iterators.product((it for i=1:N)...)

    local_h = LocalOperator{eltype(ex), 1}(ex.expr)

    @inbounds for _sites in it
        lhs_ind = Tuple(to_indices(L, lhs[region[site]]) for site in _sites)
        rhs_ind = Tuple(to_indices(L, rhs[region[site]]) for site in _sites)

        Aij += is_other_sites_the_same(n, _sites, region, lhs, rhs) ? ex.expr[lhs_ind, rhs_ind] : 0
    end
    Aij
end

# struct NonZeroIterator
# end
#
# function eval_expr(
#     ex::Sum{<:SubHam}, region,
#     ::LOpEigen{2}, ::NumberOfBody{2}, ::PlaceHolderType{Tuple{Vertex}},
#     lhs::Sites{L}) where L
#
#     Ai = []
#     @inbounds for site in Lattices.sites(region)
#         s = region[site]
#         lhs_ind = to_indices(L, lhs[s])
#         rhs_ind = to_indices(L, rhs[s])
#
#         Aij += expr.expr[lhs_ind, rhs_ind]
#     end
# end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
using Base: PkgId, loaded_modules, package_callbacks, @get!
using Base.Meta: isexpr

export @require

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

loaded(pkg) = haskey(Base.loaded_modules, pkg)

const _callbacks = Dict{PkgId, Vector{Function}}()
callbacks(pkg) = @get!(_callbacks, pkg, [])

listenpkg(f, pkg) =
  loaded(pkg) ? f() : push!(callbacks(pkg), f)

function loadpkg(pkg)
  if haskey(_callbacks, pkg)
    fs = _callbacks[pkg]
    delete!(_callbacks, pkg)
    foreach(Base.invokelatest, fs)
  end
end

function withpath(f, path)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end

function err(f, listener, mod)
  try
    f()
  catch e
    @warn """
      Error requiring $mod from $listener:
      $(sprint(showerror, e, catch_backtrace()))
      """
  end
end

function parsepkg(ex)
  isexpr(ex, :(=)) || @goto fail
  mod, id = ex.args
  (mod isa Symbol && id isa String) || @goto fail
  return id, String(mod)
  @label fail
  error("Requires syntax is: `@require Pkg=\"uuid\"`")
end

macro require(pkg, expr)
  pkg isa Symbol &&
    return Expr(:macrocall, Symbol("@warn"), __source__,
                "Requires now needs a UUID; please see the readme for changes in 0.7.")
  id, modname = parsepkg(pkg)
  pkg = :(Base.PkgId(Base.UUID($id), $modname))
  quote
    if !isprecompiling()
      listenpkg($pkg) do
        withpath($(string(__source__.file))) do
          err($__module__, $modname) do
            $(esc(:(eval($(Expr(:quote, Expr(:block,
                                            :(const $(Symbol(modname)) = Base.require($pkg)),
                                            expr)))))))
          end
        end
      end
    end
  end
end