text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using Random
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
	using ViewerGL
	GL = ViewerGL
	include("../examples/3d/views.jl")
end

function build_copFE(FV::Lar.Cells, EV::Lar.Cells)
	copFE = Lar.u_coboundary_1(FV, EV) # unsigned
	faceedges = [findnz(copFE[f,:])[1] for f=1:size(copFE,1)]

	f_edgepairs = Array{Array{Int64,1}}[]
	for f=1:size(copFE,1)
		edgepairs = Array{Int64,1}[]
		for v in FV[f]
			push!(edgepairs, [e for e in faceedges[f] if v in EV[e]])
		end
		push!(f_edgepairs, edgepairs)
	end
	copFE.=0
	for f=1:size(copFE,1)
		facepairs = sort(sort.(f_edgepairs[f]))
		# setting a reference
		copFE[f, facepairs[1][1]] = 1
		for (e1,e2) in facepairs
			@assert copFE[f, e1] != 0
			if copFE[f, e2] == 0
				copFE[f, e2] =
					-((EV[e1][1]==EV[e2][1])||(EV[e1][2]==EV[e2][2])) * copFE[f,e1]+
					((EV[e1][1]==EV[e2][2])||(EV[e1][2]==EV[e2][1])) * copFE[f,e1]
			end
		end
	end
	return copFE
end

Random.seed!(1234);

scaling = 1.2
no_sphere = 10
expl = 2.0
pt1, pt2 = 5, 10
#pt1, pt2 = 10, 10

carry = Array{Tuple{Array{Float64,2},Array{Array{Int64,1},1}}}(undef, no_sphere)
for i = 1 : no_sphere
	carry[i] = Lar.struct2lar(Lar.Struct([
		Lar.t(rand(0.5:1e-2:5, 3)...),
		Lar.sphere(rand(1:1e-2:3) * scaling)([pt1, pt2])
	]))
end
V, FV = Lar.struct2lar(Lar.Struct(carry))
FV = sort(sort.(FV))

model = CAGD.Model(V);
function getFaceEdges(fV::Array{Int,1})
    return [fV[1], fV[2]], [fV[1], fV[3]], [fV[2], fV[3]]
end
EV = unique([(map(f -> getFaceEdges(f), FV)...)...])
FE = [[3*i+1, 3*i+2, 3*i+3] for i = 0 : length(FV)-1]
CAGD.addModelCells!(model, 1, EV, signed = true)
CAGD.addModelCells!(model, 2, build_copFE(FV, EV))
model = CAGD.mergeModelVertices(model, signed_merge=true)

if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end


text:
# Use modeling and optimizer package
using JuMP
using GLPK

# Define the model
model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
@constraint(model, 6x + 8y >= 100)
@constraint(model, 7x + 12y >= 120)
@objective(model, Min, 12x + 20y)

# Optimize the model
optimize!(model)

# Show results
@show value(x);
@show value(y);
@show objective_value(model);

ext: jl


text:
module Covid19Modelling

using CSV
using Measures
using Flux
using DifferentialEquations
using DiffEqFlux
using LaTeXStrings
using Zygote
using BSON
using LinearAlgebra
using NNlib

include("utils.jl")
include("data.jl")
include("model.jl")
include("train.jl")

export train_all_districts, prediction

end


text:
module PlaymatSimulator

using Reexport

@reexport using Images
@reexport using GameOne
@reexport using DataStructures

export Animations
export in_bounds, copy_actor, kill_actor!

const SDL2 = GameOne.SDL2
Animations = include("animations.jl")
include("auth.jl")
include("IO.jl")
include("logging.jl")

function copy_actor(a::Actor)
    c = occursin(".", a.label) ? Actors.Image(a.label) : Actors.Text(a.label, a.data[:font_path])
    c.w = a.w
    c.h = a.h
    c.x = a.x + 20
    c.y = a.y + 20
    c.data = a.data

    c
end

function download_card_by_name(card_name::String, card_data::Array, found_card=false)
	for (i,c) in enumerate(card_data)
		if card_data[i]["name"] == c["name"]
			download_mtg_imgs([ c ])
			found_card = true
			break
		end
	end

	if !found_card
		@warn "$card_name not found in the database!"
	else
		card_name
	end
end

function kill_actor!(a::Actor)
    for (i,sf) in enumerate(a.surfaces)
        a.surfaces[i] = SDL2.C_NULL
    end

    for (i,tx) in enumerate(a.textures)
        a.textures[i] = SDL2.C_NULL
    end

    SDL2.FreeSurface.(a.surfaces)
    SDL2.DestroyTexture.(a.textures)
end

function download_mtg_imgs(deck_cards::Vector{Any}, save_dir="images", throttle=0.1)
	for card in deck_cards
		println("""Processing MtG card $(card["name"])""")

		try
			if "card_faces" in keys(card)

				if "image_uris" in keys(card["card_faces"][begin])
					uri_front = card["card_faces"][1]["image_uris"]["border_crop"]
					uri_back = card["card_faces"][2]["image_uris"]["border_crop"]

					img_front = get_card_img(uri_front)
					img_back = get_card_img(uri_back)

					safesave("$save_dir/$(card["name"]).png", img_front)
					safesave("$save_dir/$(card["name"]).png", img_back)
				else
					img = get_card_img(card["image_uris"]["border_crop"])
					safesave("$save_dir/$(card["name"]).png", img)
				end

			else
				img = get_card_img(card["image_uris"]["border_crop"])
				safesave("images/$(card["name"]).png", img)
			end

		catch e
			@warn e
		end

		sleep(throttle)
	end
end

end # module

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com

"

module SpinVectorO3Dynamics
    include("..\\samplers\\tools.jl")

    using ..SamplingTools: hj_info, add_noise, add_cross_talk, discrete_aggregation
    using DifferentialEquations
    using Random
    using LinearAlgebra
    using Interpolations
    using RecursiveArrayTools

    function O3_Bloch_model(theta_vector::Array{Float64, 1},
                             phi_vector::Array{Float64, 1},
                             biases::Array{Float64, 1},
                             JM::Array{Float64, 2},
                             run_time::Float64,
                             open_system::Bool, # stochastic only
                             temp::Float64, # stochastic only
                             damping::Float64, # stochastic only
                             trajectories::Int64, # stochastic only
                             at_fn,
                             bt_fn,
                             ts,
                             abstol::Float64,
                             reltol::Float64)
        nspins = size(theta_vector, 1)
        t_span = (0.0, run_time)
        i1 = collect(1:nspins)
        i2 = collect((nspins + 1):(2*nspins))
        i3 = collect((2*nspins + 1):(3*nspins))
        # convert to the Bloch equations
        Mx = sin.(theta_vector) .* cos.(phi_vector)
        My = sin.(theta_vector) .* sin.(phi_vector)
        Mz = cos.(theta_vector)
        u_0 = vcat(Mx, My, Mz)

        model! = function (du, u, p, t)
            mx, my, mz = u[i1], u[i2], u[i3]
            hx = 2.0 * at_fn(t)
            hz = 2.0 .* bt_fn(t) .* (biases  .+ (JM * mz))
            du[i1] = -hz .* my
            du[i2] = (hx .* mz) + (hz .* mx)
            du[i3] = -hx .* my
            nothing
        end

        noise_model! = function (du, u, p, t)
            mx, my, mz = u[i1], u[i2], u[i3]
            hx = 2.0 * at_fn(t)
            hz = 2.0 .* bt_fn(t) .* (biases  .+ (JM * mz))
            du[i1] = (-hz .* (damping .* mx .* mz + my)) -
                       (damping .* hx .* (mz.^2 + my.^2))

            du[i2] = (hz .* (mx - damping .* my .* mz)) +
                       (hx .* (damping .* mx .* my + mz))

            du[i3] = (-damping .* hz .* (mx.^2 + my.^2)) +
                       (hx .* (damping .* mx .* mz - my))
            nothing
        end

        σ_model! = function (du, u, p, t)
            # constants infront of wiener process
            mx, my, mz = u[i1], u[i2], u[i3]
            c = sqrt(2*damping*temp)
            du[i1] = c .* (my + mz)
            du[i2] = c .* (mx + mz)
            du[i3] = c .* (mx + my)
            nothing
        end

        # # callback to coonserve system properties
        # function g(resid,u,p,t)
        #     mx, my, mz = u[i1], u[i2], u[i3]
        #     dr = mx.^2 + my.^2 + mz.^2 .- 1.0
        #     resid[i1] = dr
        #     resid[i2] = dr
        #     resid[i3] = dr
        #
        # end
        # cb = DifferentialEquations.ManifoldProjection(g)

        # supplying jacobian can provide speed increases

        if open_system
            # f =  DifferentialEquations.ODEFunction(noise_model!, jac=noise_model_jac!)
            f = noise_model!
            prob = DifferentialEquations.SDEProblem(f, σ_model!, u_0, t_span)
            ensembleprob = EnsembleProblem(prob)
            # can be upgraded to use GPU
            sol = DifferentialEquations.solve(ensembleprob, EnsembleThreads(),
                                              trajectories=trajectories,
                                              saveat=ts, abstol=abstol,
                                              reltol=reltol, maxiters=1e8)
        else
            # f =  DifferentialEquations.ODEFunction(model!, jac=model_jac!)
            f = model!
            prob = DifferentialEquations.ODEProblem(f, u_0, t_span)
            sol = DifferentialEquations.solve(prob, saveat=ts, # callback=cb,
                                              abstol=abstol, reltol=reltol,
                                              maxiters=1e8)
        end

        return sol
    end

    function O3_model(theta_vector::Array{Float64, 1},
                         phi_vector::Array{Float64, 1},
                         biases::Array{Float64, 1},
                         JM::Array{Float64, 2},
                         run_time::Float64,
                         open_system::Bool, # stochastic only
                         temp::Float64, # stochastic only
                         damping::Float64, # stochastic only
                         trajectories::Int64, # stochastic only
                         at_fn,
                         bt_fn,
                         ts,
                         abstol::Float64,
                         reltol::Float64)
        if open_system
            @error ErrorException("Use Bloch equation for open system evaluation")
        end
        nspins = size(theta_vector, 1)
        t_span = (0.0, run_time)
        ω_θ = zeros(Float64, nspins)
        ω_ϕ = zeros(Float64, nspins)
        u_0 = ArrayPartition(theta_vector, ω_θ, phi_vector, ω_ϕ)

        model! = function (du, u, p, t)
            # θ contribution
            θ, ω, ϕ, ρ = u.x[1], u.x[2], u.x[3], u.x[4]
            hx = at_fn(t)
            hz = bt_fn(t) .* (biases + (JM * cos.(θ)))
            du.x[1] .= ω
            du.x[2] .= (-0.5 .* ρ .* sin.(θ)) + ((hx .* cos.(θ) .* cos.(ϕ)) +
                       (sin.(θ) .* hz))

            # ϕ contributions
            du.x[3] .= ρ
            du.x[4] .= (0.5 .* ω .* sin.(θ)) - (hx .* sin.(θ) .* sin.(ϕ))
            nothing
        end

        # f =  DifferentialEquations.ODEFunction(model!, jac=model_jac!)
        f = model!
        prob = DifferentialEquations.ODEProblem(f, u_0, t_span)
        sol = DifferentialEquations.solve(prob, saveat=ts, abstol=abstol,
                                          reltol=reltol, maxiters=1e8)
        return sol
    end

    function solve(h::Dict, J::Dict; anneal_time=100.0, run_time=0.0, # in nanoseconds
                    temp=0.1, friction_constant=1e-3, trajectories=1000,
                    # for individual qubit schedules, input a
                    # dictionary of interpolations
                    # NOTE if interpolations aren't continuous then
                    # instabilities arise in the ODE integrator
                    a_fn=nothing, b_fn=nothing, Bloch=true,
                    add_xtalk=false, noisy=false, open_system=false,
                    initial_theta_state=Array{Float64, 1}(),
                    initial_phi_state=Array{Float64, 1}(),
                    saveat=collect(0:0.01:1),  # points in normalized time
                    abstol=1e-9, reltol=1e-6, kwargs...)

        nspins = length(h)
        # check to see if we have individual qubit schedules
        independent_qubit_sechdules = false
        if isa(a_fn, Dict) | isa(b_fn, Dict)
            # if one input is not a dict but one if, convert to dict of lists
            if !isa(a_fn, Dict)
                interp = deepcopy(a_fn)
                a_fn = Dict([(i, interp) for i = 1:nspins])
            elseif !isa(b_fn, Dict)
                interp = deepcopy(b_fn)
                b_fn = Dict([(i, interp) for i = 1:nspins])
            end

            independent_qubit_sechdules = true
        elseif a_fn == nothing
            # linear schedules
            a_fn = LinearInterpolation(saveat, 1.0 .- saveat, extrapolation_bc=Flat())
        elseif b_fn == nothing
            b_fn = LinearInterpolation(saveat, saveat, extrapolation_bc=Flat())
        end

        if run_time == 0
            run_time = anneal_time
        end

        # define coefficient interpolation functions
        if independent_qubit_sechdules
            # vectorise the independent_schedules
            at_fn = function (t)
                a_vec = zeros(Float64, nspins)
                for i = 1:nspins
                    a_vec[i] = a_fn[i](t/anneal_time)
                end
                return a_vec
            end
            bt_fn = function (t)
                b_vec = zeros(Float64, nspins)
                for i = 1:nspins
                    b_vec[i] = b_fn[i](t/anneal_time)
                end
                return b_vec
            end
        else
            at_fn = (t) -> a_fn(t/anneal_time)
            bt_fn = (t) -> b_fn(t/anneal_time)
        end

        # prepare initial state if defined
        if length(initial_theta_state) == nspins
            if length(initial_phi_state) != nspins
                initial_phi_state = zeros(Float64, nspins)
            end
            @assert ndims(initial_theta_state) == 1
            @assert ndims(initial_phi_state) == 1
            theta_vector = deepcopy(initial_theta_state)
            phi_vector = deepcopy(initial_phi_state)
        else
            # starting in the transverse field direction as default
            theta_vector = fill(0.5 * pi, nspins)
            phi_vector = zeros(Float64, nspins)
        end

        # add noise effects
        if noisy
            h, J = add_noise(h, J)
        end
        # add cross-talk effects
        if add_xtalk
            h, J = add_xtalk(h, J)
        end

        # prepare problem info
        JM, M, degs = hj_info(h, J)
        JS = deepcopy(JM)
        biases = diag(JM)
        JM -= diagm(biases)
        JM += JM'

        if Bloch
            _fn = O3_Bloch_model
        else
            _fn = O3_model
        end

        # run the protocol
        sol = _fn(theta_vector, phi_vector, biases, JM, run_time,
                      open_system, # stochastic only
                      temp, # stochastic only
                      friction_constant, # stochastic only
                      trajectories, # stochastic only
                      at_fn, bt_fn, run_time.*saveat, abstol,
                      reltol)

        z_idxs = (2*nspins + 1):(3*nspins)
        if open_system
            spin_vectors = [Array{Float64, 1}(sign.(sol[i].u[end][z_idxs])) for i =1:trajectories]
            results = discrete_aggregation(spin_vectors, trajectories, collect(UpperTriangular(JS)))
        else
            if !Bloch
                final_sol = sol.u[end].x[1]
                spins = Array{Float64, 1}(sign.(cos.(final_sol)))
            else
                final_sol = sol.u[end][z_idxs]
                spins = Array{Float64, 1}(sign.(final_sol))
            end

            energy = biases' * spins + 0.5 * (spins' * (JM * spins))
            results = [spins, energy]
        end
        return results, sol
    end

    export solve, O3_model, O3_Bloch_model

end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
#import Pkg;
#Pkg.add("HypothesisTests")
#Pkg.add("Random")

using HypothesisTests, Random

Random.seed!(0); x= randn(100); y= randn(100);

println(OneSampleTTest(x,y))


text:
# RANDOM PARTITIONS

export RandomPartition

# This is a work in process. We'll use the Chinese Restaurant
# Process to generate a random partition of an n-set. We
# hope that the probability we generate a given partition
# is precisely 1/bell(n). We're using ideas from this website:
# http://djalil.chafai.net/blog/2012/05/03/generating-uniform-random-partitions/

# using SimpleTools

"""
`RandomPartition(n)` creates a random partition of {1,2,...,n}. The probability
we see any particular partition is (nearly?) `1/bell(n)`.

THIS IS A WORK IN PROGRESS. It doesn't seem to be working to give the same
probability to all `n`-element partitions.
"""
function RandomPartition(n::Int)
  nn = 10n
  wts = zeros(nn)
  for j=1:nn
    J = BigInt(j)
    wts[j] = Float64(BigInt(J)^n) / Float64(factorial(J))
  end
  k = random_choice(wts)
  d = Dict{Int,Int}()
  for i=1:n
    c = mod(rand(Int),n)+1
    d[i] = c
  end
  return Partition(d)
end

function RP_test(n::Int,reps::Int)
  d = Dict{Partition{Int},Int}()
  for i=1:reps
    P = RandomPartition(n)
    if haskey(d,P)
      d[P]+=1
    else
      d[P]=1
    end
  end
  for P in keys(d)
    println(d[P]/reps,"\t",P)
  end
  println("Goal = ", Float64(1/bell(n)))
  return d
end


text:
function buildNTwoSiteHorizontal(Aj::Vector{ITensor},
                                 L::Environments,
                                 R::Environments,
                                 IEnvs,
                                 row::Int,
                                 col::Int,
                                 ncol::Int,
                                 ϕ::ITensor)::ITensor
    Ny       = length(Aj)
    is_cu    = is_gpu(Aj)
    N        = is_cu ? cuITensor(1.0) : ITensor(1.0)
    for si in inds(Aj[row], "Site")
        N   *= spinI(si; is_gpu=is_cu)
    end
    workingN = N
    if row > 1
        workingN *= IEnvs[:below][row - 1]
    end
    workingN *= L.I[row]
    workingN *= ϕ
    workingN *= R.I[row]
    if row < Ny
        workingN *= IEnvs[:above][end - row]
    end
    return workingN
end

function buildHIedgeTwoSiteHorizontal(Aj::Vector{ITensor},
                                      E::Environments,
                                      row::Int,
                                      col::Int,
                                      ncol::Int,
                                      side::Symbol,
                                      ϕ::ITensor)
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    IH_a   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IH_b   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IH     = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for work_row in 1:row-1
        IH_b *= Aj[work_row] * E.H[work_row]
        IH_b *= dag(prime(Aj[work_row], "Link"))
    end
    @inbounds for work_row in row+1:Ny
        IH_a *= Aj[work_row] * E.H[work_row]
        IH_a *= dag(prime(Aj[work_row], "Link"))
    end
    IH *= E.H[row]
    IH *= IH_b
    op  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    for si in inds(Aj[row], "Site")
        op *= spinI(si; is_gpu=is_cu)
    end
    IH *= ϕ
    IH *= op
    IH *= IH_a
    AAinds = inds(prime(ϕ))
    @assert hasinds(inds(IH), AAinds) "inds(IH): $(inds(IH)),\n inds(AA): $AAinds"
    @assert hasinds(AAinds, inds(IH)) "inds(IH): $(inds(IH)),\n inds(AA): $AAinds"
    return (IH,)
end

function buildHIsTwoSiteHorizontal(Aj::Vector{ITensor},
                                   L::Environments,
                                   R::Environments,
                                   row::Int,
                                   col::Int,
                                   ncol::Int,
                                   ϕ::ITensor;
                                   edge::Symbol=:none)
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    @timeit "buildHIedge" begin
        if edge == :left
            return buildHIedgeTwoSiteHorizontal(Aj, R, row, col, ncol, :left, ϕ)
        end
        if edge == :right
            return buildHIedgeTwoSiteHorizontal(Aj, L, row, col, ncol, :right, ϕ)
        end
    end
    HLI_a = is_cu ? cuITensor(1.0) : ITensor(1.0)
    HLI_b = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR_a = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR_b = is_cu ? cuITensor(1.0) : ITensor(1.0)
    HLI   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @timeit "build HI lower" begin
        @inbounds for work_row in 1:row-1
            tL     = Aj[work_row]*L.H[work_row]
            HLI_b *= tL
            HLI_b *= R.I[work_row]
            HLI_b *= dag(prime(Aj[work_row], "Link"))
            tR     = Aj[work_row]*R.H[work_row]
            IHR_b *= tR
            IHR_b *= L.I[work_row]
            IHR_b *= dag(prime(Aj[work_row], "Link"))
        end
    end
    @timeit "build HI upper" begin
        @inbounds for work_row in reverse(row+1:Ny)
            tL     = Aj[work_row]*L.H[work_row]
            HLI_a *= tL
            HLI_a *= R.I[work_row]
            HLI_a *= dag(prime(Aj[work_row], "Link"))
            tR     = Aj[work_row]*R.H[work_row]
            IHR_a *= tR
            IHR_a *= L.I[work_row]
            IHR_a *= dag(prime(Aj[work_row], "Link"))
        end
    end
    @timeit "join up HIs" begin
        HLI  *= L.H[row]
        HLI  *= HLI_b
        HLI  *= ϕ
        HLI  *= R.I[row]
        IHR  *= L.I[row]
        IHR  *= IHR_b
        IHR  *= ϕ
        IHR  *= R.H[row]
        HLI  *= HLI_a
        IHR  *= IHR_a
        op    = is_cu ? cuITensor(1.0) : ITensor(1.0)
        for si in inds(Aj[row], "Site")
            op   *= spinI(si; is_gpu=is_cu)
        end
        HLI  *= op
        IHR  *= op
    end
    AAinds = inds(prime(ϕ))
    @assert hasinds(inds(IHR), AAinds)
    @assert hasinds(inds(HLI), AAinds)
    @assert hasinds(AAinds, inds(IHR))
    @assert hasinds(AAinds, inds(HLI))
    return (HLI, IHR)
end

function verticalTermsTwoSiteHorizontal(Aj::Vector{ITensor},
                                        L::Environments,
                                        R::Environments,
                                        AI,
                                        AV,
                                        H,
                                        row::Int,
                                        col::Int,
                                        ncol::Int,
                                        ϕ::ITensor)::Vector{ITensor}
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    vTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for opcode in 1:length(H)
        thisVert = dummy
        op_row_a = H[opcode].sites[1][1]
        op_col_a = H[opcode].sites[1][2]
        op_row_b = H[opcode].sites[2][1]
        op_col_b = H[opcode].sites[2][2]
        if op_row_b < row || row < op_row_a
            local V, I
            if op_row_a > row
                V = AV[:above][opcode][end - row]
                I = row > 1 ? AI[:below][row - 1] : dummy
            elseif op_row_b < row
                V = AV[:below][opcode][row - op_row_b]
                I = row < Ny ? AI[:above][end - row] : dummy
            end
            thisVert *= V
            thisVert *= L.I[row]
            thisVert *= ϕ
            thisVert *= R.I[row]
            thisVert *= I
            for si in inds(Aj[row], "Site")
                thisVert *= spinI(si; is_gpu=is_cu)
            end
        elseif row == op_row_a
            low_row  = op_row_a - 1
            high_row = op_row_b
            AIL = low_row > 0   ? AI[:below][low_row]        : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            if col == op_col_a
                sA   = firstind(Aj[op_row_a], "Site,c$col")
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[op_row_b], "Site,c$col")
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op_a *= spinI(firstind(Aj[op_row_a], "Site,c$ncol"); is_gpu=is_cu)
                op_b *= spinI(firstind(Aj[op_row_b], "Site,c$ncol"); is_gpu=is_cu)
            else
                sA   = firstind(Aj[op_row_a], "Site,c$ncol")
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[op_row_b], "Site,c$ncol")
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op_a *= spinI(firstind(Aj[op_row_a], "Site,c$col"); is_gpu=is_cu)
                op_b *= spinI(firstind(Aj[op_row_b], "Site,c$col"); is_gpu=is_cu)
            end
            thisVert  = AIH
            thisVert *= L.I[op_row_b]
            thisVert *= Aj[op_row_b] * op_b
            thisVert *= R.I[op_row_b]
            thisVert *= dag(Aj[op_row_b])'
            thisVert *= L.I[op_row_a]
            thisVert *= ϕ
            thisVert *= R.I[op_row_a]
            thisVert *= AIL
            thisVert *= op_a
        elseif row == op_row_b
            low_row  = op_row_a - 1
            high_row = op_row_b
            AIL = low_row > 0 ? AI[:below][low_row] : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            if op_col_a == col
                sA   = firstind(Aj[op_row_a], "Site,c$col")
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[op_row_b], "Site,c$col")
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op_a *= spinI(firstind(Aj[op_row_a], "Site,c$ncol"); is_gpu=is_cu)
                op_b *= spinI(firstind(Aj[op_row_b], "Site,c$ncol"); is_gpu=is_cu)
            else
                sA   = firstind(Aj[op_row_a], "Site,c$ncol")
                op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[op_row_b], "Site,c$ncol")
                op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op_a *= spinI(firstind(Aj[op_row_a], "Site,c$col"); is_gpu=is_cu)
                op_b *= spinI(firstind(Aj[op_row_b], "Site,c$col"); is_gpu=is_cu)
            end
            thisVert  = AIL
            thisVert *= L.I[op_row_a]
            thisVert *= Aj[op_row_a] * op_a
            thisVert *= R.I[op_row_a]
            thisVert *= dag(Aj[op_row_a])'
            thisVert *= L.I[op_row_b]
            thisVert *= ϕ
            thisVert *= R.I[op_row_b]
            thisVert *= AIH
            thisVert *= op_b
        end
        @assert hasinds(inds(thisVert), AAinds) "inds of thisVert and AAinds differ!\n$(inds(thisVert))\n$AAinds\n"
        @assert hasinds(AAinds, inds(thisVert)) "inds of thisVert and AAinds differ!\n$(inds(thisVert))\n$AAinds\n"
        if hasinds(AAinds, inds(thisVert)) && hasinds(inds(thisVert), AAinds)
            vTerms[opcode] = thisVert
        end
    end
    return vTerms
end

function interiorTermsTwoSiteHorizontal(Aj::Vector{ITensor},
                                        L::Environments,
                                        R::Environments,
                                        AI,
                                        AInt,
                                        H,
                                        row::Int,
                                        col::Int,
                                        ncol::Int,
                                        ϕ::ITensor)::Vector{ITensor}
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    intTerms = Vector{ITensor}(undef, length(H))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    AAinds = inds(prime(ϕ))
    @inbounds for opcode in 1:length(H)
        thisInter = copy(dummy)
        op_row    = H[opcode].sites[1][1]
        op_col_a  = H[opcode].sites[1][2]
        op_col_b  = H[opcode].sites[2][2]
        if op_row != row
            local IT, I
            if op_row > row
                IT = AInt[:above][opcode][end - row]
                I  = row > 1 ? AI[:below][row - 1] : dummy
            else
                IT = AInt[:below][opcode][row - 1]
                I = row < Ny ? AI[:above][end - row] : dummy
            end
            thisInter *= IT
            thisInter *= L.I[row]
            thisInter *= ϕ
            thisInter *= R.I[row]
            thisInter *= I
            for si in inds(Aj[row], "Site")
                thisInter *= spinI(si; is_gpu=is_cu)
            end
            @assert hasinds(inds(thisInter), AAinds)
            @assert hasinds(AAinds, inds(thisInter))
        else
            low_row  = op_row - 1
            high_row = op_row
            AIL = low_row > 0   ? AI[:below][low_row]        : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            thisInter  = AIL
            thisInter *= L.I[row]
            thisInter *= ϕ
            thisInter *= R.I[row]
            thisInter *= AIH
            if op_col_a == col
                sA   = firstind(Aj[row], "Site,c$col")
                op_a = copy(H[opcode].ops[1])
                op_a = replaceind!(op_a, H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[row], "Site,c$ncol")
                op_b = copy(H[opcode].ops[2])
                op_b = replaceind!(op_b, H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op = op_a * op_b
            else
                sA   = firstind(Aj[row], "Site,c$col")
                op_a = copy(H[opcode].ops[2])
                op_a = replaceind!(op_a, H[opcode].site_ind, sA)
                op_a = replaceind!(op_a, H[opcode].site_ind', sA')
                sB   = firstind(Aj[row], "Site,c$ncol")
                op_b = copy(H[opcode].ops[1])
                op_b = replaceind!(op_b, H[opcode].site_ind, sB)
                op_b = replaceind!(op_b, H[opcode].site_ind', sB')
                op = op_a * op_b
            end
            thisInter *= op
            @assert hasinds(inds(thisInter), AAinds)
            @assert hasinds(AAinds, inds(thisInter))
        end
        intTerms[opcode] = thisInter
    end
    return intTerms
end

function fieldTermsTwoSiteHorizontal(Aj::Vector{ITensor},
                                     L::Environments,
                                     R::Environments,
                                     AI,
                                     AF,
                                     H,
                                     row::Int,
                                     col::Int,
                                     ncol::Int,
                                     ϕ::ITensor)::Vector{ITensor}
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    fTerms = Vector{ITensor}(undef, length(H))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    AAinds = inds(prime(ϕ))
    @inbounds for opcode in 1:length(H)
        thisField = copy(dummy)
        op_row    = H[opcode].sites[1][1]
        op_col    = H[opcode].sites[1][2]
        if op_row != row
            local F, I
            if op_row > row
                F = AF[:above][opcode][end - row]
                I = row > 1 ? AI[:below][row - 1] : dummy
            else
                F = AF[:below][opcode][row - 1]
                I = row < Ny ? AI[:above][end - row] : dummy
            end
            thisField *= F
            thisField *= L.I[row]
            thisField *= ϕ
            thisField *= R.I[row]
            thisField *= I
            for si in inds(Aj[row], "Site")
                thisField *= spinI(si; is_gpu=is_cu)
            end
        else
            low_row  = op_row - 1
            high_row = op_row
            AIL = low_row > 0   ? AI[:below][low_row]        : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            thisField  = AIL
            thisField *= L.I[row]
            thisField *= ϕ
            thisField *= R.I[row]
            thisField *= AIH
            if op_col == col
                sA = firstind(Aj[row], "Site,c$col")
                op = copy(H[opcode].ops[1])
                op = replaceind!(op, H[opcode].site_ind, sA)
                op = replaceind!(op, H[opcode].site_ind', sA')
                op *= spinI(firstind(Aj[row], "Site,c$ncol"); is_gpu=is_cu)
            else
                sA = firstind(Aj[row], "Site,c$ncol")
                op = copy(H[opcode].ops[1])
                op = replaceind!(op, H[opcode].site_ind, sA)
                op = replaceind!(op, H[opcode].site_ind', sA')
                op *= spinI(firstind(Aj[row], "Site,c$col"); is_gpu=is_cu)
            end
            thisField *= op
        end
        @assert hasinds(inds(thisField), AAinds)
        @assert hasinds(AAinds, inds(thisField))
        fTerms[opcode] = thisField
    end
    return fTerms
end

function connectLeftTermsTwoSiteHorizontal(Aj::Vector{ITensor},
                                           L::Environments,
                                           R::Environments,
                                           AI, AL, H,
                                           row::Int, col::Int, ncol::Int,
                                           ϕ::ITensor)::Vector{ITensor}
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    lTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    lcol, rcol = extrema([col, ncol])
    @inbounds for opcode in 1:length(H)
        op_row_b = H[opcode].sites[2][1]
        op_b     = copy(H[opcode].ops[2])
        as       = firstind(Aj[op_row_b], "Site,c$lcol")
        op_b     = replaceind!(op_b, H[opcode].site_ind, as)
        op_b     = replaceind!(op_b, H[opcode].site_ind', as')
        thisHori = is_cu ? cuITensor(1.0) : ITensor(1.0)
        if op_row_b != row
            local ancL, I
            if op_row_b > row
                ancL = AL[:above][opcode][end - row]
                I = row > 1 ? AL[:below][opcode][row - 1] : dummy
            else
                ancL = AL[:below][opcode][row - 1]
                I = row < Ny ? AL[:above][opcode][end - row] : dummy
            end
            thisHori  = ancL
            thisHori *= L.InProgress[row, opcode]
            thisHori *= ϕ
            thisHori *= R.I[row]
            thisHori *= I
            for si in inds(Aj[row], "Site")
                thisHori *= spinI(si; is_gpu=is_cu)
            end
        else
            low_row  = (op_row_b <= row) ? op_row_b - 1 : row - 1;
            high_row = (op_row_b >= row) ? op_row_b + 1 : row + 1;
            if low_row >= 1
                thisHori *= AL[:below][opcode][low_row]
            end
            thisHori *= R.I[row]
            thisHori *= ϕ
            thisHori *= L.InProgress[row, opcode]
            thisHori *= op_b
            thisHori *= spinI(firstind(Aj[row], "Site,c$rcol"); is_gpu=is_cu)
            if high_row <= Ny
                thisHori *= AL[:above][opcode][end - high_row + 1]
            end
        end
        @assert hasinds(inds(thisHori), AAinds)
        @assert hasinds(AAinds, inds(thisHori))
        lTerms[opcode] = thisHori
    end
    return lTerms
end

function connectRightTermsTwoSiteHorizontal(Aj::Vector{ITensor},
                                            L::Environments,
                                            R::Environments,
                                            AI, AR, H,
                                            row::Int, col::Int, ncol::Int,
                                            ϕ::ITensor)::Vector{ITensor}
    Ny     = length(Aj)
    is_cu  = is_gpu(Aj)
    rTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    lcol, rcol = extrema([col, ncol])
    @inbounds for opcode in 1:length(H)
        op_row_a = H[opcode].sites[1][1]
        op_a     = copy(H[opcode].ops[1])
        as       = firstind(Aj[op_row_a], "Site,c$rcol")
        op_a     = replaceind!(op_a, H[opcode].site_ind, as)
        op_a     = replaceind!(op_a, H[opcode].site_ind', as')
        thisHori = is_cu ? cuITensor(1.0) : ITensor(1.0)
        if op_row_a != row
            local ancR, I
            if op_row_a > row
                ancR = AR[:above][opcode][end - row]
                I = row > 1 ? AR[:below][opcode][row - 1] : dummy
            else
                ancR = AR[:below][opcode][row - 1]
                I = row < Ny ? AR[:above][opcode][end - row] : dummy
            end
            thisHori = ancR
            thisHori *= R.InProgress[row, opcode]
            thisHori *= ϕ
            thisHori *= L.I[row]
            thisHori *= I
            for si in inds(Aj[row], "Site")
                thisHori *= spinI(si; is_gpu=is_cu)
            end
        else
            low_row  = (op_row_a <= row) ? op_row_a - 1 : row - 1;
            high_row = (op_row_a >= row) ? op_row_a + 1 : row + 1;
            if low_row >= 1
                thisHori *= AR[:below][opcode][low_row]
            end
            thisHori *= L.I[row]
            thisHori *= R.InProgress[row, opcode]
            thisHori *= ϕ
            thisHori *= op_a
            if high_row <= Ny
                thisHori *= AR[:above][opcode][end - high_row + 1]
            end
            thisHori *= spinI(firstind(Aj[row], "Site,c$lcol"); is_gpu=is_cu)
        end
        @assert hasinds(inds(thisHori), AAinds)
        @assert hasinds(AAinds, inds(thisHori))
        rTerms[opcode] = thisHori
    end
    return rTerms
end

function buildLocalHTwoSiteHorizontal(A::Vector{ITensor},
                                      L::Environments, R::Environments,
                                      AncEnvs, H,
                                      row::Int, col::Int, ncol::Int,
                                      ϕ::ITensor;
                                      verbose::Bool=false)
    field_H_terms = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Vertical)
    term_count    = 1 + length(field_H_terms) + length(vert_H_terms)
    Ny            = length(A)
    Nx            = size(H, 2)
    @timeit "build Ns" begin
        N   = buildNTwoSiteHorizontal(A, L, R, AncEnvs[:I], row, col, ncol, ϕ)
    end
    den = scalar(collect(N * dag(ϕ)'))
    if verbose
        println("--- NORM ---")
        println(den)
    end
    local left_H_terms, right_H_terms
    mincol, maxcol    = extrema([col, ncol])
    interior_H_terms  = getDirectional(vcat(H[:, mincol]...), Horizontal)
    term_count       += length(interior_H_terms)
    if mincol > 1
        left_H_terms  = getDirectional(vcat(H[:, mincol - 1]...), Horizontal)
        term_count   += length(left_H_terms)
    end
    if maxcol < Nx
        right_H_terms = getDirectional(vcat(H[:, maxcol]...), Horizontal)
        term_count   += length(right_H_terms)
    end
    if 1 < col < Nx && 1 < ncol < Nx
        term_count += 1
    end
    Hs = Vector{ITensor}(undef, term_count)
    term_counter = 1
    @debug "\t\tBuilding I*H and H*I row $row col $col"
    @timeit "build HIs" begin
        edge = (1 < col < Nx && 1 < ncol < Nx) ? :none : ((col == 1 || ncol == 1) ? :left : :right)
        HIs = buildHIsTwoSiteHorizontal(A, L, R, row, col, ncol, ϕ; edge=edge)
        Hs[term_counter] = HIs[1]
        if length(HIs) == 2
            Hs[term_counter+1] = HIs[2]
        end
        term_counter += length(HIs)
        if verbose
            println("--- HI TERMS ---")
            for HI in HIs
                println(scalar(HI * dag(ϕ)')/den)
            end
        end
    end
    @debug "\t\tBuilding vertical H terms row $row col $col"
    @timeit "build vertical terms" begin
        vTs = verticalTermsTwoSiteHorizontal(A, L, R, AncEnvs[:I], AncEnvs[:V], vert_H_terms, row, col, ncol, ϕ)
        Hs[term_counter:term_counter+length(vTs) - 1] = vTs
        term_counter += length(vTs)
        if verbose
            println( "--- vT TERMS ---")
            for vT in vTs
                println(scalar(vT * dag(ϕ)')/den)
            end
        end
    end
    @debug "\t\tBuilding field H terms row $row col $col"
    @timeit "build field terms" begin
        fTs = fieldTermsTwoSiteHorizontal(A, L, R, AncEnvs[:I], AncEnvs[:F], field_H_terms, row, col, ncol, ϕ)
        Hs[term_counter:term_counter+length(fTs) - 1] = fTs[:]
        term_counter += length(fTs)
        if verbose
            println( "--- fT TERMS ---")
            for fT in fTs
                println(scalar(fT * dag(ϕ)')/den)
            end
        end
    end
    @debug "\t\tBuilding interior H terms row $row col $col"
    @timeit "build interior terms" begin
        iTs = interiorTermsTwoSiteHorizontal(A, L, R, AncEnvs[:I], AncEnvs[:IT], interior_H_terms, row, col, ncol, ϕ)
        Hs[term_counter:term_counter+length(iTs) - 1] = iTs[:]
        term_counter += length(iTs)
        if verbose
            println( "--- interiorT TERMS ---")
            for iT in iTs
                println(scalar(iT * dag(ϕ)')/den)
            end
        end
    end
    if mincol > 1
        @debug "\t\tBuilding left H terms row $row col $col"
        @timeit "build left terms" begin
            lTs = connectLeftTermsTwoSiteHorizontal(A, L, R, AncEnvs[:I], AncEnvs[:L], left_H_terms, row, col, ncol, ϕ)
            Hs[term_counter:term_counter+length(lTs) - 1] = lTs[:]
            term_counter += length(lTs)
            if verbose
                println( "--- lT TERMS ---")
                for lT in lTs
                    println(scalar(lT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilt left terms"
    end
    if maxcol < Nx
        @debug "\t\tBuilding right H terms row $row col $col"
        @timeit "build right terms" begin
            rTs = connectRightTermsTwoSiteHorizontal(A, L, R, AncEnvs[:I], AncEnvs[:R], right_H_terms, row, col, ncol, ϕ)
            Hs[term_counter:term_counter+length(rTs) - 1] = rTs[:]
            term_counter += length(rTs)
            if verbose
                println( "--- rT TERMS ---")
                for rT in rTs
                    println(scalar(rT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilt right terms"
    end
    return Hs, N
end

function intraColumnGaugeHorizontal(A::Vector{ITensor}; kwargs...)::Vector{ITensor}
    Ny = length(A)
    @inbounds for row in reverse(2:Ny)
        @debug "\tBeginning intraColumnGauge for col $col row $row"
        cmb_is     = IndexSet(inds(A[row], "Site"))
        cmb_is     = IndexSet(cmb_is..., inds(A[row], "Link,r")...)
        cmb = combiner(cmb_is, tags="CMB")
        Lis = IndexSet(combinedind(cmb)) #cmb_is
        if row < Ny
            Lis = IndexSet(Lis..., commonind(A[row], A[row + 1]))
        end
        old_ci    = commonind(A[row], A[row-1])
        Ac        = A[row]*cmb
        U, S, V   = svd(Ac, Lis; mindim=dim(old_ci), maxdim=dim(old_ci), cutoff=0.0, utags=tags(old_ci))
        A[row]    = U*cmb
        A[row-1] *= (S*V)
        #new_ci    = commonind(A[row], A[row-1])
        #A[row]    = replaceind!(A[row], new_ci, old_ci)
        #A[row-1]  = replaceind!(A[row-1], new_ci, old_ci)
    end
    return A
end

struct ITensorMapTwoSiteHorizontal
  A::Vector{ITensor}
  H::Matrix{Vector{Operator}}
  L::Environments
  R::Environments
  AncEnvs::NamedTuple
  row::Int
  col::Int
  ncol::Int
end
Base.eltype(M::ITensorMapTwoSiteHorizontal)  = eltype(M.A[M.row])
Base.size(M::ITensorMapTwoSiteHorizontal)    = ITensors.dim(M.A[M.row])
function (M::ITensorMapTwoSiteHorizontal)(v::ITensor)
    Hs, N  = buildLocalHTwoSiteHorizontal(M.A, M.L, M.R, M.AncEnvs, M.H, M.row, M.col, M.ncol, v)
    localH = sum(Hs)
    return noprime(localH)
end

function optimizeLocalHTwoSiteHorizontal(A::Vector{ITensor},
                                         L::Environments, R::Environments,
                                         AncEnvs, H,
                                         row::Int, col::Int, ncol::Int;
                                         kwargs...)
    Ny            = length(A)
    Nx            = size(H, 2)
    is_cu         = is_gpu(A)
    field_H_terms = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]...), Vertical)
    @debug "\tBuilding H for col $col row $row"
    @timeit "build H" begin
        Hs, N = buildLocalHTwoSiteHorizontal(A, L, R, AncEnvs, H, row, col, ncol, A[row])
    end
    initial_N = real(scalar(collect(N * dag(A[row])')))
    @timeit "sum H terms" begin
        localH = sum(Hs)
    end
    initial_E = real(scalar(collect(deepcopy(localH) * dag(A[row])')))
    @info "Initial energy at row $row col $col : $(initial_E/(initial_N*Nx*Ny)) and norm : $initial_N"
    #println("Initial energy at row $row col $col : $(initial_E/(initial_N*Nx*Ny)) and norm : $initial_N")
    @debug "\tBeginning davidson for col $col row $row"
    mapper   = ITensorMapTwoSiteHorizontal(A, H, L, R, AncEnvs, row, col, ncol)
    λ, new_A = davidson(mapper, A[row]; maxiter=2, kwargs...)
    new_E    = λ
    N        = buildNTwoSiteHorizontal(A, L, R, AncEnvs[:I], row, col, ncol, new_A)
    new_N    = real(scalar(collect(N * dag(new_A)')))
    @info "Optimized energy at row $row col $col : $(new_E/(new_N*Nx*Ny)) and norm : $new_N"
    #println("Optimized energy at row $row col $col : $(new_E/(new_N*Nx*Ny)) and norm : $new_N")
    @timeit "restore intraColumnGauge" begin
        if row < Ny
            @debug "\tRestoring intraColumnGauge for col $col row $row"
            cmb_is  = IndexSet(inds(A[row], "Site"))
            cmb_is  = IndexSet(cmb_is..., inds(A[row], "Link,r")...)
            cmb     = combiner(cmb_is, tags="CMB")
            ci      = combinedind(cmb)
            Lis     = IndexSet(ci)
            if row > 1
                Lis = IndexSet(Lis..., commonind(A[row], A[row - 1]))
            end
            old_ci  = commonind(A[row], A[row+1])
            svdA    = new_A*cmb
            Ris     = uniqueinds(inds(svdA), Lis)
            U, S, V = svd(svdA, Ris; mindim=dim(old_ci), maxdim=dim(old_ci), vtags=tags(old_ci))
            new_ci  = commonind(V, S)
            #replaceind!(V, new_ci, old_ci)
            A[row]  = V*cmb
            newU    = S*U*A[row+1]
            #replaceind!(newU, new_ci, old_ci)
            A[row+1] = newU
            if row < Ny - 1
                nI    = is_cu ? cuITensor(1.0) : ITensor(1.0)
                for si in inds(A[row+1], "Site")
                    nI    *= spinI(si; is_gpu=is_cu)
                end
                newAA  = AncEnvs[:I][:above][end - row - 1]
                newAA *= L.I[row+1]
                newAA *= A[row+1] * nI
                newAA *= R.I[row+1]
                newAA *= dag(A[row+1])'
                AncEnvs[:I][:above][end - row] = newAA
            end
        else
            A[row] = initial_N > 0 ? new_A/√new_N : new_A
        end
    end
    return A, AncEnvs
end

function buildAncsHorizontal(A::fPEPS, Aj::Vector{ITensor}, L::Environments, R::Environments, H, col::Int, ncol::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @timeit "makeAncIs" begin
        Ia = makeAncillaryIsHorizontal(A, Aj, L, R, col, ncol)
        Ib = Vector{ITensor}(undef, Ny)
        Is = (above=Ia, below=Ib)
    end
    @timeit "makeAncVs" begin
        vH = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Vertical)
        Va = makeAncillaryVsHorizontal(A, Aj, L, R, vH, col, ncol)
        Vb = [Vector{ITensor}() for ii in 1:length(Va)]
        Vs = (above=Va, below=Vb)
    end
    @timeit "makeAncFs" begin
        fH = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Field)
        Fa = makeAncillaryFsHorizontal(A, Aj, L, R, fH, col, ncol)
        Fb = [Vector{ITensor}() for ii in 1:length(Fa)]
        Fs = (above=Fa, below=Fb)
    end
    mincol, maxcol = extrema([col, ncol])
    @timeit "makeAncInts" begin
        iH = getDirectional(vcat(H[:, mincol]...), Horizontal)
        ITa = makeAncillaryIntsHorizontal(A, Aj, L, R, iH, col, ncol)
        ITb = [Vector{ITensor}() for ii in 1:length(ITa)]
        ITs = (above=ITa, below=ITb)
    end
    Ls = (above=Vector{ITensor}(), below=Vector{ITensor}())
    Rs = (above=Vector{ITensor}(), below=Vector{ITensor}())
    if mincol > 1
        @timeit "makeAncLs" begin
            lH = getDirectional(vcat(H[:, mincol-1]...), Horizontal)
            La = makeAncillarySideHorizontal(A, Aj, L, R, lH, col, ncol, :left)
            Lb = [Vector{ITensor}() for ii in  1:length(La)]
            Ls = (above=La, below=Lb)
        end
    end
    if maxcol < Nx
        @timeit "makeAncRs" begin
            rH = getDirectional(vcat(H[:, maxcol]...), Horizontal)
            Ra = makeAncillarySideHorizontal(A, Aj, R, L, rH, col, ncol, :right)
            Rb = [Vector{ITensor}() for ii in  1:length(Ra)]
            Rs = (above=Ra, below=Rb)
        end
    end
    Ancs = (I=Is, V=Vs, F=Fs, L=Ls, R=Rs, IT=ITs)
    return Ancs
end

function updateAncsHorizontal(A::fPEPS, Aj::Vector{ITensor},
                              L::Environments, R::Environments,
                              AncEnvs, H,
                              row::Int, col::Int, ncol::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)

    Is, Vs, Fs, Ls, Rs, ITs = AncEnvs
    @debug "\tUpdating ancillary identity terms for col $col row $row"
    Ib = updateAncillaryIsHorizontal(A, Aj, Is[:below], L, R, row, col, ncol)
    Is = (above=Is[:above], below=Ib)

    @debug "\tUpdating ancillary vertical terms for col $col row $row"
    vH = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Vertical)
    Vb = updateAncillaryVsHorizontal(A, Aj, Vs[:below], Ib, L, R, vH, row, col, ncol)
    Vs = (above=Vs[:above], below=Vb)

    @debug "\tUpdating ancillary field terms for col $col row $row"
    fH = getDirectional(vcat(H[:, col]..., H[:, ncol]...), Field)
    Fb = updateAncillaryFsHorizontal(A, Aj, Fs[:below], Ib, L, R, fH, row, col, ncol)
    Fs = (above=Fs[:above], below=Fb)

    mincol, maxcol    = extrema([col, ncol])
    @debug "\tUpdating interior field terms for col $col row $row"
    itH = getDirectional(vcat(H[:, mincol]...), Horizontal)
    ITb = updateAncillaryIntsHorizontal(A, Aj, ITs[:below], Ib, L, R, itH, row, col, ncol)
    ITs = (above=ITs[:above], below=ITb)

    if mincol > 1
        @debug "\tUpdating ancillary left terms for col $col row $row"
        lH = getDirectional(vcat(H[:, mincol-1]...), Horizontal)
        Lb = updateAncillarySideHorizontal(A, Aj, Ls[:below], Ib, L, R, lH, row, col, ncol, :left)
        Ls = (above=Ls[:above], below=Lb)
    end
    if maxcol < Nx
        @debug "\tUpdating ancillary right terms for col $col row $row"
        rH = getDirectional(vcat(H[:, maxcol]...), Horizontal)
        Rb = updateAncillarySideHorizontal(A, Aj, Rs[:below], Ib, R, L, rH, row, col, ncol, :right)
        Rs = (above=Rs[:above], below=Rb)
    end
    Ancs = (I=Is, V=Vs, F=Fs, L=Ls, R=Rs, IT=ITs)
    return Ancs
end

function joincols(A::fPEPS, col::Int, ncol::Int)::Tuple{Vector{ITensor}, Vector{ITensor}}
    Ny, Nx = size(A)
    Aj = Vector{ITensor}(undef, Ny)
    cis = Vector{ITensor}(undef, Ny-1)
    for row in 1:Ny
        Aj[row] = A[row, col]*A[row, ncol]
        if row > 1
            Aj[row] *= cis[row-1]
        end
        if row < Ny
            ci = combiner(commonindex(A[row, col], A[row+1,col]), commonindex(A[row, ncol], A[row+1, ncol]), tags="Link,u,ji$row")
            Aj[row] *= ci
            cis[row] = ci
        end
    end
    return Aj, cis
end

function unjoincols(A::fPEPS, Aj::Vector{ITensor}, cis::Vector{ITensor}, col::Int, ncol::Int)::fPEPS
    Ny, Nx = size(A)
    for row in 1:Ny
        if row > 1
            Aj[row] *= cis[row-1]
        end
        if row < Ny
            Aj[row] *= cis[row]
        end
    end
    for row in 1:Ny
        col_is       = commoninds(Aj[row], A[row, col])
        cnc_ci       = commonind(A[row, col], A[row, ncol])
        U, S, V      = svd(Aj[row], col_is; mindim=dim(cnc_ci), maxdim=dim(cnc_ci), lefttags=tags(cnc_ci), cutoff=0.)
        A[row, col]  = U
        A[row, ncol] = S*V
    end
    D = dim(commonind(A[1, col], A[1, ncol]))
    return A
end

function sweepColumnHorizontal(A::fPEPS,
                               L::Environments, R::Environments,
                               H, col::Int, ncol::Int;
                               kwargs...)
    Ny, Nx = size(A)
    #println("Sweeping column $col with next column $ncol")
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    @debug "Beginning buildAncs for col $col"
    Aj, jis = joincols(A, col, ncol)
    @timeit "intraColumnGauge Aj" begin
        Aj = intraColumnGaugeHorizontal(Aj; kwargs...)
    end
    AncEnvs = buildAncsHorizontal(A, Aj, L, R, H, col, ncol)
    @inbounds for row in 1:Ny
        if row > 1
            @timeit "updateAncs" begin
                AncEnvs = updateAncsHorizontal(A, Aj, L, R, AncEnvs, H, row-1, col, ncol)
            end
        end
        @debug "Beginning optimizing H for col $col"
        Aj, AncEnvs = optimizeLocalHTwoSiteHorizontal(Aj, L, R, AncEnvs, H, row, col, ncol; kwargs...)
    end
    #A = unjoincols(A, Aj, jis, col, ncol)
    #return A
    return A, Aj
end


text:
"""
    mmread(filename|io)

Read `Matrixmarket` format file (extension `.mtx`) and return sparse or dense matrix.
Symmetric and Hermitian matrices use the corresponding wrapper types.
Patterns result in sparse matrics with element type `Bool`.
They may be converted to numerical types by multiplying with a number.
"""
function mmread(filename::AbstractString)
    open(filename, "r") do file
        if stat(file).size == 0
            throw(DataError("matrix file $filename is empty"))
        end
        mmread(file)
    end
end

using Mmap

const COORD = "coordinate"
const ARRAY = "array"
const MATRIX = "matrix"
const MATRIXM = "%%matrixmarket"
const COMPLEX = "complex"
const REAL = "real"
const INTEGER = "integer"
const PATTERN = "pattern"
const GENERAL = "general"
const SYMMETRIC = "symmetric"
const HERMITIAN = "hermitian"
const SKEW_SYMMETRIC = "skew-symmetric"

function mmread(file::IO)
    line = lowercase(readline(file))
    tokens = split(line)
    if tokens[1] != MATRIXM
        parserr(string("Matrixmarket: invalid header:", line))
    end
    line = readline(file)
    while length(strip(line)) == 0 || line[1] == '%'
        line = readline(file)
    end
    if tokens[2] == MATRIX
        mmread_matrix(file, line, tokens[3:end]...)
    else
        parserr(string("Matrixmarket: unsupported type: ", line))
    end
end

# mmap for regular files - else read
function getbytes(io::IOStream)
    isfile(io) ? Mmap.mmap(io, grow=false, shared=false) : read(io)
end
getbytes(io::IO) = read(io)

function mmread_matrix(file::IO, line, form, field, symm)
    FMAP = Dict(REAL => (3, Float64),
                COMPLEX => (4, ComplexF64),
                INTEGER => (3, Int64),
                PATTERN => (2, Bool))

    ty, T = get(FMAP, field) do
            parserr("Matrixmarket: unsupported field $field (only real/complex/pattern)")
    end

    SMAP = Dict(GENERAL => (1, 0, Any),
                SYMMETRIC => (0, 1, Symmetric),
                SKEW_SYMMETRIC => (1, 1, Array),
                HERMITIAN => (0, 1, Hermitian))

    p1, pc, wrapper = get(SMAP, symm) do
        parserr("Matrixmarket: unsupported symmetry $symm (general/symmetric/hermitian,skew-symmetric)")
    end

    if form == COORD
        m, n, nz = parseint(line)
        b = getbytes(file)
        rv = Vector{Int}(undef, nz)
        cv = Vector{Int}(undef, nz)
        vv = Vector{T}(undef, nz)
        parseloop!(Val(ty), b, rv, cv, vv)
        result = mksparse!(m, n, rv, cv, vv)
    elseif form == ARRAY
        m, n = parseint(line)
        b = getbytes(file)
        p = 1
        result = zeros(T, m, n)
        for c = 1:n
            for r = (pc*c+p1):m
                p, v = parsenext(T, b, p)
                result[r,c] = v
            end
        end
    else
        parserr("Matrixmarket: unsupported format '$form'")
    end

    wrap(result, wrapper)
end

wrap(result, ::Type{T}) where T<:Union{Symmetric,Hermitian} = T(result, :L)
wrap(result, ::Type{Array}) = result - mtranspose(result)
wrap(result, ::Type{Any}) = result

function parseloop!(::Val{4}, c::Vector{UInt8}, rv, cv, vv::Vector{T}) where T<:Complex
    nz = length(rv)
    R = real(T)
    p = 1
    for i = 1:nz
        p, rv[i] = parsenext(Int, c, p)
        p, cv[i] = parsenext(Int, c, p)
        p, r = parsenext(R, c, p)
        p, s = parsenext(R, c, p)
        vv[i] = r + s*im
    end
end

function parseloop!(::Val{3}, c::Vector{UInt8}, rv, cv, vv::Vector{T}) where T <:Real
    nz = length(rv)
    p = 1
    for i = 1:nz
        p, rv[i] = parsenext(Int, c, p)
        p, cv[i] = parsenext(Int, c, p)
        p, vv[i] = parsenext(T, c, p)
    end
end

function parseloop!(::Val{2}, c::Vector{UInt8}, rv, cv, vv::Vector{T}) where T<:Number

    nz = length(rv)
    p = 1
    for i = 1:nz
        p, rv[i] = parsenext(Int, c, p)
        p, cv[i] = parsenext(Int, c, p)
    end
    fill!(vv, T(1))
end

function parseint(line::AbstractString)
    tokens = split(line)
    parse.(Int, tokens)
end

"""
    mksparse(m, n, rowval, colval, nzval)
    mksparse!(m, n, rowval, colval, nzval)

Construct a `SparseMatrixCSC` of dimensions `(m,n)` from the data given in the
three input vectors of equal lengths.

`mksparse!` destoys the content of `colval`.

`A[rowval[i],colval[i]] == nzval[i] for i ∈ 1:length(nzval)`. All other entries are zero.
"""
mksparse(m, n, rv, cv, vv) = mksparse!(m, n, rv, copy(cv), vv)

function mksparse!(m::Integer, n::Integer, rv::AbstractVector{Ti}, cv::AbstractVector{Ti},
                  vv::AbstractVector{Tv}) where {Ti<:Integer,Tv<:Number}

    nz = length(rv)
    length(cv) == nz == length(vv) || argerr("all vectors need same length")
    micv, mcv = extrema(cv)
    mirv, mrv = extrema(rv)
    micv > 0 && mcv <= n || daterr("all column indices must be >= 1 and <= $n")
    mirv > 0 && mrv <= m || daterr("all row indices must be >= 1 and <= $m")
    sizeof(Ti) <= 8 || argerr("Index type greater 64 bits not supported")
    sr = count_ones(Ti(nextpow(2, mrv + 1) - 1))
    sh = count_zeros(Ti(nextpow(2, mcv + 1) - 1))
    sr < sh || argerr("combined size of column and row indices exceeds $Ti size")
    # Ti must be able to keep (nz + 1)
    (nz + 1) % Ti != nz + 1 && argerr("Ti($Ti) cannot store nz($nz)")
    # compress row, col into Ti
    # t = UInt64[]
    # push!(t, time_ns())
    colptr = zeros(Ti, n+1)
    colptr[1] = 1
    # push!(t, time_ns())
    @inbounds for i = 1:nz
        cvi = cv[i]
        cv[i] = cvi << sr | rv[i]
        colptr[cvi+1] += 1
    end
    # push!(t, time_ns())
    cumsum!(colptr, colptr)
    # push!(t, time_ns())
    p = specialsort(cv, sr)
    # push!(t, time_ns())
    # println("times: $(diff(t) ./ 1e6) ms")
    SparseMatrixCSC{Tv,Ti}(m, n, colptr, rv[p], vv[p])
end

function specialsort(cv::Vector{Int}, sr::Int)
    nz = length(cv)
    if nz > 10000
        x = isqrt(nz) << sr
        p = sortperm(cv .÷ x)
        sortperm!(p, cv, initialized=true)
    else
        sortperm(cv)
    end
end

"""
    colval(A)
reconstruct colum-indices from colptr of `SparseMatrixCSC`.
"""
function colval(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    nz = nnz(A)
    cv = Vector{Ti}(undef, nz)
    colptr = A.colptr
    for j = 1:A.n
        for i = colptr[j]:colptr[j+1]-1
            cv[i] = j
        end
    end
    cv
end

"""
    mtranspose(A)
Materialized transpose of a matrix
"""
mtranspose(A::SparseMatrixCSC) = mksparse!(A.n, A.m, colval(A), copy(A.rowval), copy(A.nzval))
mtranspose(A::Matrix) = Matrix(transpose(A))
mtranspose(A) = transpose(A)

"""
    madjoint(A)
Materialized adjoint of sparse Matrix
"""
madjoint(A) = conj!(mtranspose(A))

"""
    mmreadcomment(filename)
return info comment strings for MatrixMarket format files
"""
function mmreadcomment(filename::AbstractString)
    io = IOBuffer()
    mark(io)
    open(filename,"r") do mmfile
        skip = 0
        while !eof(mmfile)
            s = readline(mmfile)
            skip = isempty(strip(s)) || s[1] == '%' ? 0 : skip + 1
            skip <= 1 && println(io, s)
            if skip == 1
                length(split(s)) == 3 && break
                reset(io)
                mark(io)
            end
        end
    end
    String(take!(io))
end

"""
    mmreadheader(filename)
Read header information from mtx file.
"""
function mmreadheader(file::AbstractString)
    if isfile(file)
        open(file) do io
          if stat(io).size == 0
            return nothing
          end
          line = lowercase(readline(io))
          while true
            token = split(line)
            if length(token) >= 4 &&
                token[1] == MATRIXM &&
                token[2] == MATRIX &&
                token[3] in [COORD, ARRAY]

                hdr = Dict{Symbol,Any}()
                field = :none
                while startswith(line, '%') || isempty(strip(line))
                    field = push_hdr!(hdr, line, field)
                    line = readline(io)
                end
                res = try parseint(line) catch; [] end
                if length(res) != (token[3] == COORD ? 3 : 2)
                    daterr("MatrixMarket file '$file' invalid sizes: '$line'")
                end
                hdr[:m] = res[1]
                hdr[:n] = res[2]
                length(res) >= 3 && (hdr[:nz] = res[3])
                hdr[:format] = token[3]
                hdr[:field] = token[4]
                hdr[:symmetry] = token[5]
                if haskey(hdr, :notes)
                    hdr[:notes] = join(wordlist(String(take!(hdr[:notes]))), ' ')
                end
                if haskey(hdr, :date)
                    val = hdr[:date]
                    hdr[:date] = isempty(val) ? 0 : parse(Int, hdr[:date])
                end
                if length(hdr) >= 6 && get(hdr, :nz, 0) > 0
                    return hdr
                else
                    while !eof(io) && !startswith(line, '%')
                        line = readline(io)
                    end
                    if eof(io)
                        return hdr
                    end
                    line = lowercase(line)
                end
            else
                daterr("file '$file' is not a MatrixMarket file")
            end
          end
        end
    else
        nothing
    end
end

"""
    wordlist(string)

Separate words is string by spaces and delimiters.
Return list of unique words, which can be used as keywords.
"""
function wordlist(s::AbstractString)
    list = unique!(split(s, r"[][\s(){}`\"'*]", keepempty = false))
    list = replace.(list, Ref(r"[.:,;']$" => ""))
    # remove all lowercase words with less than ... chars
    list = filter!(x->!(length(x)<4 && all(islowercase.(collect(x))) || length(x) < 2), list)
    unique!(list)
end

function push_hdr!(hdr, line::AbstractString, field::Symbol)
    isempty(strip(line)) && return field
    if field == :notes
        if !startswith(line, "%---")
            println(get!(hdr, field) do; IOBuffer() end, strip(line[2:end]))
        end
        return field
    end
    reg = r"^% *([^:[]+): *(.*)$"
    regtitle = r"^% *\[([^]]*)]"
    if (m = match(reg, line)) !== nothing
        s = Symbol(m[1])
        if s in (:name, :kind, :ed, :fields, :author, :date)
            field = s
            value = strip(m[2])
            hdr[field] = value
        elseif s == :notes
            field = s
        end
    elseif (m = match(regtitle, line)) !== nothing
        field = :title
        value = m[1]
        hdr[field] = value
    end
    field
end

### parsing decimal integers and floats
function _parsenext(v::Vector{UInt8}, p1::Int)
    iaccu = Unsigned(0)
    daccu = Unsigned(0)
    eaccu = 0
    df = 0
    sig = 0
    esig = 0
    i = p1
    n = length(v)
    c = v[i]
    while c == 0x20 || c == 0x0a || c == 0x0d || c == 0x09
        c = v[i += 1]
    end
    n0 = i
    if c == UInt8('+')
        sig = 1
        c = v[i += 1]
    elseif c == UInt8('-')
        sig = -1
        c = v[i += 1]
    end
    ne = i
    di = 0
    while 0x30 <= c <= 0x39
        c -= 0x30
        di += di > 0 || c > 0
        iaccu = iaccu * 10 + c
        c = v[i += 1]
    end
    if c == UInt8('.')
        c = v[i += 1]
        d0 = i
        while 0x30 <= c <= 0x39
            c -= 0x30
            di += di > 0 || c > 0
            daccu = daccu * 10 + c
            c = v[i += 1]
        end
        df = i - d0
        ne += 1
    end
    i > ne || parserr("Invalid decimal number: '$(String(v[n0:min(end,n0+5)]))'")
    if i > ne && ( c == UInt8('e') || c == UInt8('E') )
        c = v[i += 1]
        if c == UInt8('+')
            esig = 1
            c = v[i += 1]
        elseif c == UInt8('-')
            esig = -1
            c = v[i += 1]
        end
        while 0x30 <= c <= 0x39
            eaccu = eaccu * 10 + c - 0x30
            c = v[i += 1]
        end
        if esig < 0
            eaccu = -eaccu
        end
    end
    i == n0 && parserr("No decimal number found: '$(String(v[n0:min(end,n0+5)]))'")
    i, iaccu, daccu, eaccu, sig, df, di
end

function parsenext(T::Type{<:Signed}, v::Vector{UInt8}, p1::Int)
    i, iaccu, daccu, eaccu, sig, df = _parsenext(v, p1)
    daccu == 0 && eaccu == 0 && df == 0 || error("1")
    i, T(iaccu) * ifelse(sig < 0, T(-1), T(1))
end
function parsenext(T::Type{<:Unsigned}, v::Vector{UInt8}, p1::Int)
    i, iaccu, daccu, eaccu, sig, df = _parsenext(v, p1)
    daccu == 0 && eaccu == 0 && sig >= 0 && df == 0 || error("2")
    i, T(iaccu)
end
function parsenext(T::Type{<:AbstractFloat}, v::Vector{UInt8}, p1::Int)
    i, iaccu, daccu, eaccu, sig, df, di = _parsenext(v, p1)
    eaccu -= df
    if di <= 16 && iaccu != 0
        daccu += 10^df * iaccu
    end
    f = if di <= 16 && daccu < UInt(1)<<53
        if 0 <= eaccu < 23
            T(daccu) * exp10(eaccu)
        elseif 0 < -eaccu < 23
            T(daccu) / exp10(-eaccu)
        else
            parse(T, String(view(v, p1:i-1)))
        end
    else
        parse(T, String(view(v, p1:i-1)))
    end
    i, (sig < 0 ? -f : f)
end
function parsenext(T::Type{<:Complex}, c, p)
    R = real(T)
    p, r = parsenext(R, c, p)
    p, s = parsenext(R, c, p)
    p, r + s*im
end


text:
module WakeData

# package code goes here

include("reformat.jl")
include("join.jl")
include("stat.jl")
include("model.jl")

end # module


text:
#Circuit analysis tools
#-------------------------------------------------------------------------------
module CircuitAnalysis

import MDDatasets: DS #To dispatch on a symbol & minimize namespace pollution

include("base.jl")
include("functions.jl")
include("unitconvert.jl")
include("show.jl")

#==Exported symbols
===============================================================================#
#Objects
#-------------------------------------------------------------------------------

#Functions
#-------------------------------------------------------------------------------
export dB20, dB10, dB, dBm, dBW, Vpk, Ipk, VRMS, IRMS

#Value tags:
#-------------------------------------------------------------------------------
export impedance
export admittance
export capacitance
export inductance

#==Other interface tools (symbols not exported to avoid collisions):
================================================================================
#Already in base:
==#

#==Symbols *not* exported (useful for implementing custom functions)
================================================================================
#Tagged data types:
	TImpedance{T}
	TAdmittance{T}
	TInductance{T}
	TCapacitance{T}
=#

end #module CircuitAnalysis


text:
#####################
# Circuit line types to interface between file parsers and circuit compilers
#####################

"""Circuit and vtree node ids used for IO"""
const ID = UInt32

"""A parsed circuit file format line"""
abstract type CircuitFormatLine <: FormatLine end

"""A file consisting for circuit formal lines"""
const CircuitFormatLines = AbstractVector{<:CircuitFormatLine}

"""A string comment line for circuit files"""
struct CircuitCommentLine <: CircuitFormatLine
    comment::AbstractString
end

"""A header line for circuit files"""
abstract type CircuitHeaderLine <: CircuitFormatLine end

"""A circuit format line without child IDs"""
abstract type LeafCircuitLine <: CircuitFormatLine end

"""A circuit format line with child IDs"""
abstract type InnerCircuitLine <: CircuitFormatLine end

"""A line that represents a logical literal"""
abstract type LiteralLine <: LeafCircuitLine end

struct SddHeaderLine <: CircuitHeaderLine
    num_nodes::Int
end

struct PsddHeaderLine <: CircuitHeaderLine
    num_nodes::Int
end

struct LcHeaderLine <: CircuitHeaderLine end

"""
A line representing a weighted single literal (for example a logistic circuit literal).
"""
struct WeightedLiteralLine <: LiteralLine
    node_id::ID
    vtree_id::ID
    literal::Lit
    normalized::Bool
    weights::Vector{Float64}
end

"""
A line representing a single literal without parameters.
"""
struct UnweightedLiteralLine <: LiteralLine
    node_id::ID
    vtree_id::ID
    literal::Lit
    normalized::Bool
end

is_normalized(l::LiteralLine) = l.normalized

import ...LogicCircuits: literal, variable # extend
literal(l::LiteralLine) = l.literal
variable(l::LiteralLine) = lit2var(literal(l))

"""
A line representing either the true or false logical constants
"""
abstract type ConstantLine <: LeafCircuitLine end

"""
A weighted constant line for a known variable
"""
struct WeightedNamedConstantLine <: ConstantLine
    node_id::ID
    vtree_id::ID
    variable::Var
    # constant::Bool always true because one cannot associate a weight with the models of false
    # normalized::Bool always true because we have a single variable field given above, which only exists if the node is normalized
    weight::Float64
end

struct AnonymousConstantLine <: ConstantLine
    node_id::ID
    constant::Bool
    normalized::Bool
end

is_normalized(l::AnonymousConstantLine) = l.normalized
is_normalized(::WeightedNamedConstantLine) = true
variable(l::WeightedNamedConstantLine) = l.variable
lnconstant(::WeightedNamedConstantLine) = true
lnconstant(l::AnonymousConstantLine) = l.constant

"""
Paired boxes, or elements, are conjunctions
in a larger decision node line
"""
abstract type Element end
abstract type NormalizedElement <: Element end

struct LCElement <: NormalizedElement
    prime_id::ID
    sub_id::ID
    weights::Vector{Float64}
end

struct PSDDElement <: NormalizedElement
    prime_id::ID
    sub_id::ID
    weight::Float64
end

abstract type TrimmedElement <: Element end

struct SDDElement <: TrimmedElement
    prime_id::ID
    sub_id::ID
end

is_normalized(::NormalizedElement) = true
is_normalized(::TrimmedElement) = false

"""
A line representing a decision node in the circuit (an OR of AND elements)
"""
struct DecisionLine{ET<:Element} <: InnerCircuitLine
    node_id::ID
    vtree_id::ID
    num_elements::UInt32
    elements:: Vector{ET}
end

"""
A line representing a bias node in the circuit (an OR with one child)
"""
struct BiasLine <: InnerCircuitLine
    node_id::ID
    weights::Vector{Float64}
    BiasLine(weights) = new(typemax(ID), weights)
end

"""
String representations of each line type
"""
# TODO rename this IO module so as to not conflict with Core.IO below
Base.show(io::Core.IO, ln::SddHeaderLine) =  print(io, "sdd $(ln.num_nodes)")
Base.show(io::Core.IO, ln::PsddHeaderLine) =  print(io, "psdd $(ln.num_nodes)")
Base.show(io::Core.IO, ln::LcHeaderLine) =  print(io, "Logistic Circuit")
Base.show(io::Core.IO, ln::CircuitCommentLine) =  print(io, "c $(ln.comment)")
Base.show(io::Core.IO, ln::DecisionLine) =  print(io, "D $(ln.node_id) $(ln.vtree_id) $(ln.num_elements) $(ln.elements)")
Base.show(io::Core.IO, ln::BiasLine) =  print(io, "B $(ln.weights)")
Base.show(io::Core.IO, ln::WeightedNamedConstantLine) =  print(io, "T $(ln.node_id) $(ln.vtree_id) $(ln.variable) $(ln.weight)")
Base.show(io::Core.IO, ln::UnweightedLiteralLine) =  print(io, "L $(ln.node_id) $(ln.vtree_id) $(ln.literal)")

Base.show(io::Core.IO, ln::WeightedLiteralLine) = begin
    @assert ln.normalized
    if ln.literal > 0
        print(io, "T $(ln.node_id) $(ln.vtree_id) $(ln.literal) $(ln.weights)")
    else
        print(io, "F $(ln.node_id) $(ln.vtree_id) $(-ln.literal) $(ln.weights)")
    end
end

Base.show(io::Core.IO, ln::AnonymousConstantLine) = begin
    @assert !ln.normalized
    if ln.constant
        print(io, "T $(ln.node_id)")
    else
        print(io, "F $(ln.node_id)")
    end
end

Base.show(io::Core.IO, e::SDDElement) =  print(io, "$(e.prime_id) $(e.sub_id)")
Base.show(io::Core.IO, e::PSDDElement) =  print(io, "$(e.prime_id) $(e.sub_id) $(e.weight)")
Base.show(io::Core.IO, e::LCElement) =  print(io, "($(e.prime_id) $(e.sub_id) $(e.weights))")

Base.show(io::Core.IO, v::Vector{<:Union{Element, AbstractFloat}}) =  join(io, v, " ")


text:
@doc doc"""
    Euclidean{T<:Tuple} <: Manifold

Euclidean vector space $\mathbb R^n$.

# Constructor

    Euclidean(n)

generates the $n$-dimensional vector space $\mathbb R^n$.

   Euclidean(m, n)

generates the $mn$-dimensional vector space $\mathbb R^{m \times n}$, whose
elements are interpreted as $m \times n$ matrices.
"""
struct Euclidean{T<:Tuple} <: Manifold where {T} end

Euclidean(n::Int) = Euclidean{Tuple{n}}()
Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}()

function representation_size(::Euclidean{Tuple{n}}, ::Type{T}) where {n,T<:Union{MPoint, TVector, CoTVector}}
    return (n,)
end

function representation_size(::Euclidean{Tuple{m,n}}, ::Type{T}) where {m,n,T<:Union{MPoint, TVector, CoTVector}}
    return (m,n)
end

@generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...)

struct EuclideanMetric <: RiemannianMetric end

@traitimpl HasMetric{Euclidean,EuclideanMetric}

function local_metric(::MetricManifold{<:Manifold,EuclideanMetric}, x)
    return Diagonal(ones(SVector{size(x, 1),eltype(x)}))
end

function inverse_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x)
    return local_metric(M, x)
end

det_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = one(eltype(x))

log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = zero(eltype(x))

@inline inner(::Euclidean, x, v, w) = dot(v, w)
@inline inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)

distance(::Euclidean, x, y) = norm(x-y)
norm(::Euclidean, x, v) = norm(v)
norm(::MetricManifold{<:Manifold,EuclideanMetric}, x, v) = norm(v)

exp!(M::Euclidean, y, x, v) = (y .= x + v)

log!(M::Euclidean, v, x, y) = (v .= y - x)

function zero_tangent_vector!(M::Euclidean, v, x)
    fill!(v, 0)
    return v
end

project_point!(M::Euclidean, x) = x

function project_tangent!(M::Euclidean, w, x, v)
    w .= v
    return w
end

"""
    projected_distribution(M::Euclidean, d, [x])

Wraps standard distribution `d` into a manifold-valued distribution. Generated
points will be of similar type to `x`. By default, the type is not changed.
"""
function projected_distribution(M::Euclidean, d, x)
    return ProjectedPointDistribution(M, d, project_point!, x)
end

function projected_distribution(M::Euclidean, d)
    return ProjectedPointDistribution(M, d, project_point!, rand(d))
end

"""
    normal_tvector_distribution(S::Euclidean, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(M::Euclidean{Tuple{N}}, x, σ) where N
    d = Distributions.MvNormal(zero(x), σ)
    return ProjectedTVectorDistribution(M, x, d, project_tangent!, x)
end


text:
#=
the online version of QAPLIB – A Quadratic Assignment Problem Library by R.E. Burkard, S.E. Karisch and F. Rendl, (Journal of Global Optimization 10:391-403, 1997.)
https://coral.ise.lehigh.edu/data-sets/qaplib/
=#

D = Vector{Matrix{Float64}}()
F = Vector{Matrix{Float64}}()
push!(D, [0 5 6; 5 0 3.6; 6 3.6 0])
push!(F, [0 10 3; 10 0 6.5; 3 6.5 0])


text:
# Training script for h5 files using Julia to coordinate batches
using go

if length(ARGS) < 4
    println(STDERR, "usage: train_hf5.jl input_hf5file MODEL_NAME num_epochs examples_per_epoch")
    exit()
end

INPUT = ARGS[1]
MODEL_NAME = ARGS[2]
EPOCHS = parse(Int, ARGS[3])
NUM_EXAMPLES = parse(Int, ARGS[4])

# Set the feature extractors
feats = go.LIBERTIES

println(STDERR, "Initializing model...")
clf = go.CLF_DCNN(feats)

println(STDERR, "Training...")
tm = time()
checkpoint_path = "../models/$(MODEL_NAME)_checkpoint.hf5"
# Save the model structure
go.save_model(clf, "../models/", MODEL_NAME, save_weights=false, save_yaml=true)
go.julia_train_h5(clf, INPUT,
                  epochs=EPOCHS,
                  examples_per_epoch=NUM_EXAMPLES,
                  n_validation_examples=10000,
                  checkpoint_path=checkpoint_path)
println("Took $(time() - tm) seconds")

println(STDERR, "Saving model...")
go.save_model(clf, "../models/", MODEL_NAME, save_weights=true, save_yaml=false)


text:
"""
	interface.jl

Defines the interface to each UncProgram in the UncNLPrograms enviroment.
The interface should feel familiar to interacting with NLPModels.jl,
apart of the JuliaSmoothOptimizers orginization.
"""

# macro from NLPModels.jl
macro lencheck(l, vars...)
    exprs = Expr[]
    for var in vars
      varname = string(var)
      push!(exprs, :(
        if length($(esc(var))) != $(esc(l))
          throw(DimensionError($varname, $(esc(l)), length($(esc(var)))))
        end
      ))
    end
    Expr(:block, exprs...)
end

"""
    obj

```math
f(x)
````
Evaluate the objective function at the point x.
"""
function obj(nlp::UncProgram, x::AbstractArray{<:Real})
    @lencheck nlp.n x
    return nlp.f(x)
end

"""
    grad

```math
∇f(x)
```
Evaluate the gradient of the objective function at a point x.
"""
function grad(nlp::UncProgram, x::AbstractArray{<:Real})
    @lencheck nlp.n x
    g = zeros(length(x))
    return nlp.g!(g, x)
end

"""
    objgrad

```math
f(x), ∇f(x)
```
Evaluate the gradient and it's objective function at a point x,
by a iterating over the dimesions in an economical fashian.
"""
function objgrad(nlp::UncProgram, x::AbstractArray{<:Real})
    @lencheck nlp.n x
    g = zeros(length(x))
    return nlp.fg!(x, g)
end

"""
    hessAD

```math
∇²f(x)
```
Evaluate the full Hessian matrix at a point x by computing
the gradient of the Jacobian formula using ForwardDiff.jl,
a forward automatic differentation package.
"""
function hessAD(nlp::UncProgram, x::AbstractArray{<:Real})
    @lencheck nlp.n x
    g = zeros(length(x))
    return ForwardDiff.jacobian(nlp.g!, g, x)
end

"""
    adjdim!

```math
f(x), ∇f(x)
```
Change the dimensions of a specified problem in the TestSet
"""
function adjdim!(nlp::UncProgram, n::Int)
    #@warn "adjdim!: This operation has not been tested to gaurentee similair output"
    n, x0 = nlp.init(n)
    nlp = UncProgram(nlp, n, x0)
    TestSet[nlp.name] = nlp
    return nlp
end

"""
    Programs

```math
f(x), ∇f(x)
```
A list of unconstrained nonlinear programming problems in the current testing enviroment.
"""
function Programs()
    for nlp in values(TestSet)
        @printf "%s with dimension %d\n" nlp.name nlp.n
    end
    return keys(TestSet)
end # changes must be propogated to gHSTest

"""
    SelectProgram

```math
f(x), ∇f(x)
```
Returns an instance of UncProgram.
"""
function SelectProgram(key::String; n=nothing)
    if key ∈ keys(TestSet)
        !isa(n, Nothing) && adjdim!(TestSet[key], n)
        return TestSet[key]
    end
    @warn "The program $key is not in the testing enviroment"
end


text:
module FolderStorage

using AbstractStorage, Random, Serialization

struct Folder <: Container
    foldername::String
    nretry::Int
end

"""
    Folder(name[; nretry=10])

Return a representation of a POSIX folder.
"""
Folder(foldername; nretry=10) = Folder(isabspath(foldername) ? foldername : normpath(joinpath(pwd(), foldername)), nretry)

"""
    Container(Folder, d::Dict)

Return a representation of a POSIX folder where, for example, `d=Dict("foldername"=>"name", "nretry"=>10)`.
"""
AbstractStorage.Container(::Type{Folder}, d::Dict, session=nothing; nretry=10) =
    Folder(d["foldername"]; nretry = get(d, "nretry", nretry))

"""
    mkpath(c::Folder)

equivalent to mkpath(c.foldername).
"""
Base.mkpath(c::Folder) = mkpath(c.foldername)
Base.mkpath(c::Folder, o::AbstractString) = mkpath(joinpath(c.foldername, splitpath(o)[1:end-1]...))

function writebytes(c::Folder, o::AbstractString, data::AbstractArray{UInt8})
    mkpath(c, o)
    filename = joinpath(c.foldername, o)
    for i = 1:c.nretry
        write(filename, data) == length(data) && return nothing
        @warn "problem writing to $c/$o, attempt $i."
        sleep(0.1*2^(i-1))
    end
    error("problem writing to $c/$o in 10 attempts.")
end

"""
    write(c::Folder, filename, data)

Equivalent to write(joinpath(c.foldername, filename), data)
"""
function Base.write(c::Folder, o::AbstractString, data::AbstractString)
    mkpath(c, o)
    write(joinpath(c.foldername, o), data)
end

function Base.write(c::Folder, o::AbstractString, data::DenseArray{T}) where {T<:Number}
    databytes = unsafe_wrap(Array, convert(Ptr{UInt8}, pointer(data)), sizeof(data), own=false)
    writebytes(c, o, databytes)
end

function Base.write(c::Folder, o::AbstractString, data::DenseArray)
    io = IOBuffer()
    serialize(io, data)
    databytes = take!(io)
    writebytes(c, o, databytes)
end

function byterange(iblock, block_size, block_remainder)
    isremainder = iblock <= block_remainder
    firstbyte = (iblock - 1)*block_size + (isremainder ? iblock : block_remainder + 1)
    lastbyte = firstbyte + (isremainder ? block_size : block_size - 1)
    firstbyte,lastbyte
end

function readbytes_thread(c, o, data, threadid, thread_size, thread_remainder)
    firstbyte,lastbyte = byterange(threadid, thread_size, thread_remainder)
    io = open(joinpath(c.foldername, o))
    seek(io, firstbyte-1)
    read!(io, view(data, firstbyte:lastbyte))
    close(io)
    nothing
end

function readbytes!(c::Folder, o::String, data::Vector{UInt8})
    _nthreads = clamp(Threads.nthreads(), 1, length(data))
    thread_size, thread_remainder = divrem(length(data), _nthreads)
    @sync for threadid = 1:_nthreads
        @async Threads.@spawn readbytes_thread(c, o, data, threadid, thread_size, thread_remainder)
    end
    data
end

Base.read(c::Folder, o::AbstractString, ::Type{String}) = read(joinpath(c.foldername, o), String)

function Base.read(c::Folder, o::AbstractString)
    databytes = readbytes!(c, o, Vector{UInt8}(undef, filesize(c, o)))
    io = IOBuffer(databytes)
    deserialize(io)
end

"""
    Base.read!(c::Folder, filename, String, data)

Equivalent to `read!(joinpath(c.foldername, filename), String, data).`
"""
function Base.read!(c::Folder, o::AbstractString, data::DenseArray{T}) where {T<:Number}
    databytes = unsafe_wrap(Array, convert(Ptr{UInt8}, pointer(data)), (sizeof(data),))
    readbytes!(c, o, databytes)
    data
end

function Base.read!(c::Folder, o::AbstractString, data::DenseArray{T}) where {T}
    databytes = Vector{UInt8}(undef, filesize(c, o))
    readbytes!(c, o, databytes)
    io = IOBuffer(databytes)
    data .= deserialize(io)
    data
end

function Base.deepcopy(src::Folder)
    dst = Folder(src.foldername*"-copy-"*randstring(4))
    isdir(src.foldername) && cp(src, dst)
    dst
end

Base.filesize(c::Folder, o::AbstractString) = filesize(joinpath(c.foldername, o))
Base.cp(src::Folder, dst::Folder) = cp(src.foldername, dst.foldername, force=true)
Base.copy(src::Folder) = Folder(src.foldername*"-copy-"*randstring(4))
Base.readdir(src::Folder) = readdir(src.foldername)
Base.isdir(src::Folder) = isdir(src.foldername)

AbstractStorage.scrubsession(c::Folder) = c

AbstractStorage.minimaldict(c::Folder) = Dict("foldername"=>c.foldername)

function Base.isfile(c::Folder, object::AbstractString)
    # TODO: this feels like a kludge.  the trouble is that the file can
    # be written to jointpath(c,object) or a set of files, one for each
    # thread.
    isfile(joinpath(c.foldername, object)) || isfile(string(joinpath(c.foldername, object), "-1"))
end

"""
    rm(c::Folder)

Equivalent to `rm(c.foldername)`.
"""
function Base.rm(c::Folder)
    for itry = 1:c.nretry
        try
            rm(c.foldername, recursive=true, force=true)
            break
        catch
            @warn "unable to remove $(c.foldername), trial $itry"
            itry == c.nretry && rethrow()
        end
        sleep(0.1*2^(itry-1))
    end
    nothing
end

"""
    rm(c::Folder, filename)

Equivalent to `rm(joinpath(c.foldername, filename))`.
"""
function Base.rm(c::Folder, filename)
    for itry = 1:c.nretry
        try
            rm(joinpath(c.foldername, filename), recursive=true, force=true)
            break
        catch
            @warn "unable to remove $(joinpath(c.foldername, filename)), trial $itry"
            itry == c.nretry && rethrow()
        end
        sleep(0.1*2^(itry-1))
    end
    nothing
end

export Folder

end


text:
using DelimitedFiles
using StatsBase

function solve(input)
    boxesIds = readdlm(input)
    vectorIds = map(vectorize, boxesIds)
    # part 1
    counts = map(letterCounts,vectorIds)
    groups = sum(counts)
    checksum = groups[1] * groups[2]
    println("checksum: $checksum")
    # part 2
    N = length(vectorIds)
    for i in 1:N
        for j in i+1:N
            diff = findall(!iszero,vectorIds[i]-vectorIds[j])
            if length(diff) == 1
                println(boxesIds[i])
                println(boxesIds[j])
                println("i $i j $j diff $diff")
                result = string(boxesIds[i][1:diff[1]-1],boxesIds[i][diff[1]+1:end])
                println("result: $result")
                break
            end
        end
    end

end

function vectorize(s)
    return Vector{UInt8}(s)
end

function letterCounts(v)
    c=counts(v)
    return [in(2,c),in(3,c)]
end

solve("input.txt")


text:
@info "Train baseline lstm language model using enwik8 dataset..."

using Knet

include("../src/data.jl")
include("../src/baseline.jl")
include("../src/train.jl")

datadir = "../data/enwik8"
jld2dir = "../jld2/enwik8.jld2"
BATCHSIZE = 32

if !isfile(jld2dir)
    println("Reading data from directory: $datadir")
    println("Setting batch size to $BATCHSIZE")
    vocab = Vocab("$datadir/train.txt")
    trainfile = TextReader("$datadir/train.txt", vocab)
    validfile = TextReader("$datadir/valid.txt", vocab)
    testfile = TextReader("$datadir/test.txt", vocab)
    dtrn = TextData(trainfile, batchsize=BATCHSIZE)
    ddev = TextData(validfile, batchsize=BATCHSIZE)
    dtst = TextData(testfile, batchsize=BATCHSIZE)
    println("Saving data from $jld2dir")
    Knet.save(jld2dir, "dtrn", dtrn, "dtst", dtst, "ddev", ddev)
else
    println("Loading data from $jld2dir")
    (dtrn, dtst, ddev) = Knet.load(jld2dir, "dtrn", "dtst", "ddev")
    vocab = dtrn.src.vocab
    if dtrn.batchsize != BATCHSIZE
        changebatchsize!(dtrn, BATCHSIZE)
        changebatchsize!(ddev, BATCHSIZE)
        changebatchsize!(dtst, BATCHSIZE)
    end;
end;

println()
@info "Initializing the model and collecting training data..."
epochs, em_size, hidden_size, layers = 30, 256, 256, 2;

println("embedding size: ", em_size);
println("hidden size: ", hidden_size);
println("layers: ", layers);
println("epochs: ", epochs)
ctrn = collect(dtrn);
trn = collect(flatten(shuffle!(ctrn) for i in 1:epochs));
trnmini = ctrn[1:20];
dev = collect(ddev);
model = SimpleLSTMModel(em_size, hidden_size, vocab; layers=layers, dropout=0.2);
# initopt!(model, length(trn))

println()
@info "Starting training, total iteration no: $(length(trn))"
bestmodel = trainadam!(model, trn, dev)

# println()
# @info "Finished training, Starting evaluation ..."
# trnloss = loss(bestmodel, dtrn);
# println("Training set scores:       ", report(trnloss))
devloss = loss(bestmodel, ddev);
println("Development set scores:    ", report(devloss))
testloss = loss(bestmodel, dtst);
println("Test set scores:           ", report(testloss))

@info "Generate text using the trained model"
print(generate(bestmodel, start="United Nations ", maxlength=1024))

# @info "Saving the model as model.jld2"
# Knet.save("model.jld2", "model", model); # uncomment for saving


text:
using RealContinuedFractions
using Test
using TestSetExtensions

@testset ExtendedTestSet "RealContinuedFractions.jl" begin
    @testset "constructor" begin
        @test ContinuedFraction(Float64[1,2,3]) == ContinuedFraction{Vector{Float64}}([1,2,3])
        @test ContinuedFraction(Float64[1,2,3]) != ContinuedFraction{Vector{Float64}}([1,2,3,4])
    end
    @testset "contfrac" begin
        @testset "0" begin
            @test contfrac(0) == ContinuedFraction{Vector{Int}}([0])
            @test contfrac(0.) == ContinuedFraction{Vector{Int}}([0])
            @test contfrac(0//1) == ContinuedFraction{Vector{Int}}([0])
            @test contfrac(Int8, 0) == ContinuedFraction{Vector{Int8}}([0])
            @test contfrac(Float64, 0) == ContinuedFraction{Vector{Float64}}([0])
        end
        @testset "π" begin
            @test contfrac(π, 5) == ContinuedFraction([3, 7, 15, 1, 292])
            @test_throws InexactError contfrac(Int8, π, 5)
            @test contfrac(Float64, π, 5) == ContinuedFraction([3., 7, 15, 1, 292])
            @test contfrac(big(π), 75) == ContinuedFraction{Vector{Int}}([
                3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2,
                2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2,
                2, 6, 3, 5, 1, 1, 6, 8, 1, 7, 1, 2, 3, 7, 1, 2, 1, 1, 12,
                1, 1, 1, 3, 1, 1, 8, 1, 1, 2, 1, 6, 1, 1, 5, 2, 2, 3, 1, 2])
        end
    end
    @testset "fromcontfrac" begin
        @test fromcontfrac(contfrac(π, 4)) == 355//113
        @test fromcontfrac(Int32, contfrac(big(π), 15)) == 245850922//78256779
        @test_throws OverflowError fromcontfrac(Int32, contfrac(big(π), 25))
        @test fromcontfrac(contfrac(big(π), 25)) == 8958937768937//2851718461558
        @test fromcontfrac(Int, contfrac(big(π), 25)) == 8958937768937//2851718461558
        @test_throws OverflowError fromcontfrac(Int, contfrac(big(π), 35))
        @test_throws OverflowError fromcontfrac(Int128, contfrac(big(π), 80))
        @test fromcontfrac(Float64, contfrac(big(π), 25)) == 3.141592653589793
    end
    @testset "roundtrip" begin
        @testset "contfrac ∘ fromcontfrac" begin
            for T in (UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64), len in 1:50
                q = filter(r -> r > 0, rand(T, len))
                isempty(q) && (q = T[0])
                q[end] == 1 && (q[end] = 2)
                cf = ContinuedFraction{Vector{T}}(q)
                @test contfrac(T, fromcontfrac(BigInt, cf)) == cf
            end
        end
        @testset "fromcontfrac ∘ contfrac" begin
            for _ in 1:100
                x = randn()
                @test fromcontfrac(Float64, contfrac(x, 1024)) ≈ x
                r = rationalize(x)
                @test fromcontfrac(contfrac(r)) == r
            end
        end
    end
    @testset "convergents" begin
        @test convergents(contfrac(π, 12)) == [
            3//1, 22//7, 333//106, 355//113, 103993//33102, 104348//33215, 208341//66317,
            312689//99532, 833719//265381, 1146408//364913, 4272943//1360120, 5419351//1725033]
        @test convergents(contfrac(big(π), 25))[end-5:end] == [
            14885392687//4738167652, 21053343141//6701487259, 1783366216531//567663097408,
            3587785776203//1142027682075, 5371151992734//1709690779483, 8958937768937//2851718461558]
    end
    @testset "convergent" begin
        cf = contfrac(big(π), 35)
        @test convergent(cf) == convergents(cf)[end]
        @test convergent(cf) == -6844906428587781119//7849984890429250040
        @test convergent(BigInt, cf) == 265099323460521503743//84383735478118508040
    end
end

max_stars_repo_name: FedericoStra/RealContinuedFractions

max_issues_repo_name: FedericoStra/RealContinuedFractions

max_forks_repo_name: FedericoStra/RealContinuedFractions


text:
# Conjugates for
#
#   Beta - Bernoulli
#   Beta - Binomial
#

posterior_canon(pri::Beta, ss::BernoulliStats) = Beta(pri.α + ss.cnt1, pri.β + ss.cnt0)
posterior_canon(pri::Beta, ss::BinomialStats) = Beta(pri.α + ss.ns, pri.β + (ss.ne * ss.n - ss.ns))

complete(G::Type{Bernoulli}, pri::Beta, p::Float64) = Bernoulli(p)

# specialized fit_map and posterior_randmodel methods for Binomial
#
# n is needed to create a Binomial distribution (which can not be provided through complete)
#

fit_map(pri::Beta, ss::BinomialStats) = Binomial(ss.n, posterior_mode(pri, ss))
fit_map(pri::Beta, G::Type{Binomial}, data::(Int, Array)) = fit_map(pri, suffstats(G, data))
fit_map(pri::Beta, G::Type{Binomial}, data::(Int, Array), w::Array) = fit_map(pri, suffstats(G, data, w))

posterior_randmodel(pri::Beta, ss::BinomialStats) = Binomial(ss.n, posterior_rand(pri, ss))

function posterior_randmodel(pri::Beta, G::Type{Binomial}, data::(Int, Array))
	posterior_randmodel(pri, suffstats(G, data))
end

function posterior_randmodel(pri::Beta, G::Type{Binomial}, data::(Int, Array), w::Array)
	posterior_randmodel(pri, suffstats(G, data, w))
end


text:
using FactCheck

include("tux.jl")

exitstatus()

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo