text:
using PyFenicsAD
using FiniteDifferences
using Zygote
using PyCall
using Test

firedrake = pyimport("firedrake")
firedrake_adjoint = pyimport("firedrake_adjoint")
ufl = pyimport("ufl")

mesh = firedrake.UnitSquareMesh(3, 2)
V = firedrake.FunctionSpace(mesh, "P", 1)

function assemble_firedrake(u, kappa0, kappa1)

    x = firedrake.SpatialCoordinate(mesh)
    f = x[1]

    inner, grad, dx = ufl.inner, ufl.grad, ufl.dx
    J_form = 0.5 * inner(kappa0 * grad(u), grad(u)) * dx - kappa1 * f * u * dx
    # J = firedrake.assemble(J_form)
    # Avoid PyCall's automatic conversion to Julia types
    J = pycall(firedrake.assemble, PyObject, J_form)
    return J
end

templates = (firedrake.Function(V), firedrake.Constant(0.0), firedrake.Constant(0.0))
inputs = (ones(V.dim()), ones(1) * 0.5, ones(1) * 0.6)

zygote_assemble_firedrake(inputs...) =
    evaluate_primal(assemble_firedrake, templates, inputs...)[1]

Zygote.@adjoint function zygote_assemble_firedrake(inputs...)
    pyout = pycall(evaluate_primal, PyObject, assemble_firedrake, templates, inputs...)
    numpy_output, firedrake_output, firedrake_inputs, tape =
        [get(pyout, PyObject, i) for i = 0:3]

    function pullback_fun(g)
        vjp_out = evaluate_pullback(firedrake_output, firedrake_inputs, tape, g)
    end

    return get(pyout, 0), pullback_fun
end

hh(inputs...) = zygote_assemble_firedrake(inputs...)
hh0(x) = hh(x, inputs[2], inputs[3])
hh1(y) = hh(inputs[1], y, inputs[3])
hh2(z) = hh(inputs[1], inputs[2], z)

@testset "zygote_assemble_forward" begin
    out = zygote_assemble_firedrake(inputs...)

    u1 = firedrake.interpolate(firedrake.Constant(1.0), V)
    J = assemble_firedrake(u1, firedrake.Constant(0.5), firedrake.Constant(0.6))
    @test out[1] == PyCall.convert(Float64, J)
end

@testset "zygote_assemble_vjp" begin
    out, pullback_fun = Zygote.pullback(hh, inputs...)

    g = ones(size(out))
    vjp_out = pullback_fun(g)

    fdm = FiniteDifferences.central_fdm(2, 1)

    fdm_jac1 = FiniteDifferences.jacobian(fdm, hh0, inputs[1])
    fdm_jac2 = FiniteDifferences.jacobian(fdm, hh1, inputs[2])
    fdm_jac3 = FiniteDifferences.jacobian(fdm, hh2, inputs[3])

    @test isapprox(vjp_out[1], fdm_jac1[1]', atol = 1e-5)
    @test isapprox(vjp_out[2], fdm_jac2[1]', atol = 1e-5)
    @test isapprox(vjp_out[3], fdm_jac3[1]', atol = 1e-5)
end

mesh = firedrake.UnitSquareMesh(6, 5)
V = firedrake.FunctionSpace(mesh, "P", 1)

function solve_firedrake(kappa0, kappa1)

    x = firedrake.SpatialCoordinate(mesh)
    f = x[1]

    u = firedrake.Function(V)
    bcs = [firedrake.DirichletBC(V, firedrake.Constant(0.0), "on_boundary")]

    inner, grad, dx = ufl.inner, ufl.grad, ufl.dx
    JJ = 0.5 * inner(kappa0 * grad(u), grad(u)) * dx - kappa1 * f * u * dx
    v = firedrake.TestFunction(V)
    F = firedrake.derivative(JJ, u, v)
    firedrake.solve(F == 0, u, bcs = bcs)
    return u
end

templates = (firedrake.Constant(0.0), firedrake.Constant(0.0))
inputs = (ones(1) * 0.5, ones(1) * 0.6)

zygote_solve_firedrake(inputs...) =
    evaluate_primal(solve_firedrake, templates, inputs...)[1]

Zygote.@adjoint function zygote_solve_firedrake(inputs...)
    pyout = pycall(evaluate_primal, PyObject, solve_firedrake, templates, inputs...)
    numpy_output, firedrake_output, firedrake_inputs, tape =
        [get(pyout, PyObject, i) for i = 0:3]

    function pullback_fun(g)
        vjp_out = vjp_out = evaluate_pullback(firedrake_output, firedrake_inputs, tape, g)
    end

    return get(pyout, 0), pullback_fun
end

@testset "zygote_solve_forward" begin
    out = zygote_solve_firedrake(inputs...)
    u = solve_firedrake(firedrake.Constant(0.5), firedrake.Constant(0.6))
    @test isapprox(out, to_numpy(u))
end

@testset "zygote_solve_pullback" begin
    out, pullback_fun = Zygote.pullback(zygote_solve_firedrake, inputs...)
    # g = np.ones_like(numpy_output)
    g = ones(size(out))
    vjp_out = pullback_fun(g)
    @test isapprox(vjp_out[1], [-1.13533303])
    @test isapprox(vjp_out[2], [0.94611086])
end


text:
using DataFrames
using CSV
using DelimitedFiles
using Glob
using Statistics
using LinearAlgebra
include("../src/StochasticBlockModel.jl")

println("Start execution")

println(pwd())
time_limit = 20.0
output_path = "experiments/experiment-S1-milp.txt"
verbose = false
file_num_start = 0

### -------------------------------------------------------------------
all_datasets = glob("instances/S1/in/*.in")

println(string("Datasets found: ", length(all_datasets)))

symm_break_const_list = [true]
symm_w_const_list = [true]

### -------------------------------------------------------------------

method = "milp"

df = DataFrame(instance_path = String[],
               dataset_name = String[],
               n = Int[],
               m = Int[],
               q = Int[],
               avg_degree = Float64[],
               method = String[],
               LB = Float64[],
               UB = Float64[],
               status = String[],
               solvetime = Float64[],
               nodecount = Int[],
               lazycount = Int[],
               symm_break_const = Bool[],
               symm_w_const = Bool[],
               strong_assort = Bool[],
               weak_assort = Bool[],
               z = Matrix{Int}[],
               w = Matrix{Float64}[]
               )

# Write column names of DataFrame to output file
df_columns = string(join(string.(names(df)), "\t"),"\n")
open(output_path, "a+") do io
   write(io, df_columns)
end

for (num_file, filepath) in enumerate(all_datasets)
    if num_file < file_num_start
        continue
    end

    dataset_path = split(filepath, "/")[end]
    println("File number: $num_file")
    println(filepath)

    # Load Dataset
    dataset = StochasticBlockModel.Dataset(filepath)
    println(string("Size: ", dataset.n))

    # enforce symmetry on w
    for symm_w in symm_w_const_list
        # symmetry breaking constraints
        for symm_break in symm_break_const_list
            # Estimate SBM
            sbm, z, opt_results = StochasticBlockModel.MILP(dataset, time_limit=time_limit, symm_break_const=symm_break, symm_w_const=symm_w)

            avg_degree = mean(dataset.k)
            nodecount = opt_results.nodecount
            if isnothing(nodecount)
                nodecount = 0
            end
            lazycount = opt_results.lazycount
            if isnothing(lazycount)
                lazycount = 0
            end

            # Strong assortativity
            strong_assort = (minimum(diag(sbm.w)) >= maximum([sbm.w[i,j]
                            for i=1:dataset.n_communities, j=1:dataset.n_communities if j>i]))

            # Weak assortativity
            weak_assort = all([sbm.w[i,i] >= maximum([sbm.w[i,j] for j=1:dataset.n_communities if i!=j])
                                                    for i=1:dataset.n_communities])

            row = (filepath, dataset_path, dataset.n, dataset.m, dataset.n_communities, avg_degree, method,
                    opt_results.LB, opt_results.UB, string(opt_results.status), opt_results.solvetime,
                    nodecount, lazycount, symm_break, symm_w, strong_assort,
                    weak_assort, z, sbm.w)

            push!(df, row)

            # Write new line to output file
            CSV.write(output_path, df[end:end,:], delim="\t", append=true)

            # Save prob matrix w and assignments x
            res_folder = joinpath(join(split(output_path,"/")[1:end-1], "/"), string("n",dataset.n))
            if ~isdir(res_folder)
                mkdir(res_folder)
            end

            # Save w
            w_savepath = replace(joinpath(res_folder, dataset_path),".in"=>"-$method.w")
            writedlm(w_savepath, sbm.w, '\t')

            # Save z
            z_savepath = replace(joinpath(res_folder, dataset_path),".in"=>"-$method.z")
            writedlm(z_savepath, z, '\t')
        end
    end
end

println("End of execution.")

ext: jl

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
function testLogistic(p::Float64, C::Float64, m::Int64) where F<:Function
  testFunction(BernoulliFactory.logistic, m, C*p/(1.0+C*p), makef(p), C,
  GLOBAL_RNG)
end

Random.seed!(192837465)

testLogistic(0.4, 0.5, ntrials)
testLogistic(0.2, 1.0, ntrials)
testLogistic(0.1, 2.0, ntrials)
testLogistic(0.4, 2.0, ntrials)
testLogistic(0.4, 5.0, ntrials)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#* $Id$
#*--------------------------------------------------------------------
type Resistor
	Resistor()=begin
		new(
			resistance(),
			wire(Dict{Symbol,Any}(
				:Brief=>"Inlet",
				:PosX=>0.3831,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			wire(Dict{Symbol,Any}(
				:Brief=>"Outlet",
				:PosX=>0.3529,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			[
				:(inlet.V - outlet.V = R * outlet.i),
				:(outlet.i = inlet.i),
			],
			[
				"","",
			],
			[:R,],
			[:inlet,:outlet,]
		)
	end
	R::resistance
	inlet::wire
	outlet::wire
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Resistor
function setEquationFlow(in::Resistor)
	addEquation(1)
	addEquation(2)
end
function atributes(in::Resistor,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Resistor"
	fields[:Brief]="Electrical Resistor."
	drive!(fields,_)
	return fields
end
Resistor(_::Dict{Symbol,Any})=begin
	newModel=Resistor()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
########### Defines the parameters & plot for Fig 2 b,d,f,h,j ###########
c = 4.0 # set the host contact rate as a constant 4.0
# import required variables/modules common to all simulations
include("../Model_scripts/fig2params.jl");
include("../Model_scripts/Functions.jl");
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N'
########### Figure 2b ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "c=4.0\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_2b.png"
)
########### Figure 2d ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculate weighted H0
H0_w = calculate_H0()
# calculate weighted V0
V0_w = calculate_V0()
# plot the weighted R0
plot_reproductive_rate(5, size(Np)[1], "$figure_directory/graph_2d.png")
########### Figure 2f ###########
# plot the average V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_2f.png")
########### Figure 2h ###########
# calculate the weighted virulence (for noise reduction)
vir_w_avg = calculate_average_virulence()
# plot the average horizontal transmission & weighted virulence
plot_horizontal_and_virulence(1, size(Np)[1], "$figure_directory/graph_2h.png")
########### Figure 2j ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_2j.png")


text:
"""
    MSD!(data::Dict, alpha=10)
This will add a variable "MSD" to `data` which is created by `read_dump`. "MSD" is the abbreviation of Mean Square Displacement (MSD), which could be used to calculate diffusion coefficients with the Einstein relationship. `alpha` is the ration between simulation time and the maximum of MSD calculating times.
"""
function MSD!(data::Dict, alpha=10)
    try
        data["coord_corr"]
    catch
        PBC!(data)
    end
    # Reading Input
    num_steps = data["num_steps"]
    num_atoms = data["num_atoms"]
    coord_corr = copy_array(data["coord_corr"])

    # MSD
    num_msd_steps = convert(Int64, round(num_steps / alpha))
    msd = zeros(num_msd_steps, num_atoms, 4)
    for step = 0 : num_msd_steps-1
        msd[step+1, :, 1:3] = mean((coord_corr[1+step:end, :, :] .- coord_corr[1:end-step, :, :]).^2, dims=1)
        msd[step+1, :, 4] = sum(msd[step+1, :, 1:3], dims=2)
    end
    data["MSD"] = msd
    println("\"MSD\"\t\thas been added to `data`")
end

"""
    diffusion(data::Dict, alpha=10)
This will return the diffusion coefficients of each atoms in `data`.

Notes: time_step shoulb be added by `time_step!()` before calling this function.
"""
function diffusion(data::Dict, alpha=10)
    try
        data["MSD"]
    catch
        MSD!(data, alpha)
    end

    try
        data["time_step"]
    catch
        error("\"time_step\" is not contained in `data`. Please call `time_step!()` before calling this function")
    end

    # Reading Input
    MSD = copy_array(data["MSD"])
    time_step = prod(data["time_step"])
    num_atoms = data["num_atoms"]

    # Diffusion
    num_step_msd = size(MSD, 1)
    t = [i*time_step for i=0:num_step_msd-1]
    res = zeros(num_atoms)

    for atom = 1:num_atoms
        res[atom] = leastsq(t, MSD[:, atom, 4])[2]
    end
    return res
end


text:
module TestParse

using Test
using OpenQASM.Types
using YaoCompiler
using YaoCompiler.Intrinsics
using YaoCompiler.QASM
using RBNF: Token

qasm"""OPENQASM 2.0;
include "qelib1.inc";

gate post q {x q;}
"""

@testset "qasm gate codegen" begin
    ast = @code_qasm gate = true cu3(0.1, 0.2, 0.3)
    @test ast isa Gate
    @test length(ast.decl.cargs) == 3
    @test length(ast.decl.qargs) == 2
    @test length(ast.body) == 5
    @test ast.decl.cargs[1].str == "theta"
    @test ast.decl.cargs[2].str == "phi"
    @test ast.decl.cargs[3].str == "lambda"

    inst = ast.body[1]
    @test inst isa Instruction
    @test inst.name == "u1"
    @test length(inst.cargs) == 1
    @test length(inst.qargs) == 1
    @test length(inst.cargs[1]) == 3
    @test length(inst.cargs[1][1]) == 3
    carg = inst.cargs[1][1]
    @test carg[1].str == "lambda"
    @test carg[2].str == "-"
    @test carg[3].str == "phi"
    @test inst.cargs[1][2].str == "/"
    @test inst.cargs[1][3].str == "2"
end

@testset "@code_qasm optimize=true passes=:julia gate=true" begin
    ast = @code_qasm optimize = true passes = :julia gate = true cu3(0.1, 0.2, 0.3)
    @test ast.body[1] isa Instruction
    @test ast.body[1].name == "rz"
    @test ast.body[1].cargs[1] isa Token{:int}
    @test ast.body[1].cargs[1].str == "0"

    @test ast.body[2] isa Instruction
    @test ast.body[2].name == "ry"
    @test ast.body[2].cargs[1] isa Token{:int}
    @test ast.body[2].cargs[1].str == "0"

    @test ast.body[3] isa Instruction
    @test ast.body[3].name == "rz"
    @test ast.body[3].cargs[1] isa Tuple
    @test ast.body[3].cargs[1][1][1].str == "lambda"
    @test ast.body[3].cargs[1][1][2].str == "-"
    @test ast.body[3].cargs[1][1][3].str == "phi"
    @test ast.body[3].cargs[1][2].str == "/"
    @test ast.body[3].cargs[1][3].str == "2.0"

    @test ast.body[4] isa Instruction
    @test ast.body[4].name == "CX"
    @test isempty(ast.body[4].cargs)
    @test length(ast.body[4].qargs) == 2

    @test ast.body[5] isa Instruction
    @test ast.body[5].name == "rz"
    @test ast.body[5].cargs[1].fn == :(-)
    @test ast.body[5].cargs[1].arg isa Tuple
    @test ast.body[5].cargs[1].arg[1].str == "theta"
    @test ast.body[5].cargs[1].arg[2].str == "/"
    @test ast.body[5].cargs[1].arg[3].str == "2.0"

    @test ast.body[6] isa Instruction
    @test ast.body[6].name == "ry"
    @test ast.body[6].cargs[1].str == "0"
end

@device function circuit()
    1:3 => ccx()
    (2, 3) => cx()
    2 => h()
    c = @measure 1
    if c == 1
        3 => post()
    end
    3 => rx(1.0)
    return (c = c,)
end

@testset "branches & const inline" begin
    ast = @code_qasm optimize = true circuit()
    @test ast isa MainProgram
    # @test ast.prog[3] isa Instruction
    # @test ast.prog[3].name == "h"
end

end # TestParse


text:
using LinearAlgebra

mutable struct NesterovAcceleration{R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}}
    k::Int
    s::T
end

function NesterovAcceleration(
    x::T,
) where {R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}}
    NesterovAcceleration{R,C,T}(0, zero(x))
end

function update!(L::NesterovAcceleration{R,C,T}, s, y) where {R,C,T}
    L.s .= s
    L.k += 1
    return L
end

function reset!(L::NesterovAcceleration{R,C,T}) where {R,C,T}
    L.k = 0
end

import Base: *

function (*)(L::NesterovAcceleration, v)
    w = similar(v)
    mul!(w, L, v)
end

import LinearAlgebra: mul!

function mul!(d::T, L::NesterovAcceleration{R,C,T}, v::T) where {R,C,T}
    if L.k == 0
        d .= 0
    elseif L.k == 1
        d .= v
    else
        d .= v .- (L.k - 1) / (L.k + 2) .* L.s
    end
end


text:
@testset "Particle Group 2D2V" begin

    n_particles = 10
    charge = -1.0
    mass = 1.0
    n_weights = 1

    particle_group =
        ParticleGroup{2,2}(n_particles, charge = 1.0,
                           mass = 1.0, n_weights = 1)

    for ipart = 1:n_particles
        GEMPIC.set_x!(particle_group, ipart, [ipart, 0.0])
        GEMPIC.set_v!(particle_group, ipart, [ipart^2, 0.0])
        GEMPIC.set_weights!(particle_group, ipart, ipart / n_particles)
    end

    alpha = 0.1
    kx = 0.5

    mesh = TwoDGrid( 0, 2π/kx, 128, 0, 1, 16)

    sampler = LandauDamping(alpha, kx)

    ParticleInCell.sample!(particle_group, mesh, sampler)

    @test true

end


text:
#
# Generate random numbers with a normal distribution
#
# mean: the mean
# sigma: the standard deviation
#

normal(mean :: Float64, sigma :: Float64) = sigma*sqrt(-2*log(rand()))*cos(2*pi*rand()) + mean


text:
using Images, TestImages
using Base.Test

img = testimage("autumn_leaves")
@assert colorspace(img) == "RGBA"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint16

img = testimage("cameraman")
@assert colorspace(img) == "Gray"
@assert ndims(img) == 2
@assert colordim(img) == 0
@assert eltype(img) == Uint8

img = testimage("earth_apollo17")
@assert colorspace(img) == "RGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8

img = testimage("fabio")
@assert colorspace(img) == "Gray"
@assert ndims(img) == 2
@assert colordim(img) == 0
@assert eltype(img) == Uint8

img = testimage("house")
@assert colorspace(img) == "GrayAlpha"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8

img = testimage("jetplane")
@assert colorspace(img) == "GrayAlpha"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8

img = testimage("lighthouse")
@assert colorspace(img) == "RGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8

img = testimage("mandrill")
@assert colorspace(img) == "RGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8

img = testimage("moonsurface")
@assert colorspace(img) == "Gray"
@assert ndims(img) == 2
@assert colordim(img) == 0
@assert eltype(img) == Uint8

img = testimage("mountainstream")
@assert colorspace(img) == "RGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8


text:
using StanModels, Distributions

# Define the Stan language model

m2_1s = "
// Inferring a Rate
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
";

# Define the Stanmodel.

sm = SampleModel("m2.1s", m2_1s);

# Use 16 observations

N = 15
d = Binomial(9, 0.66)
k = rand(d, N)
n = repeat([9], N)

# Input data for cmdstan

m2_1_data = Dict("N" => N, "n" => n, "k" => k);

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m2_1_data);

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
using Test
using LinearAlgebra: norm
using NonNegLeastSquares: nonneg_lsq
using SparseArrays

#test specific
using Random
using PyCall
const pyopt = pyimport_conda("scipy.optimize", "scipy")

function test_case1()
    A = [ 0.53879488  0.65816267
          0.12873446  0.98669198
          0.24555042  0.00598804
          0.80491791  0.32793762 ]

    b = [0.888, 0.562, 0.255, 0.077]
    x = [0.15512102, 0.69328985] # approx solution from scipy
    return A, b, x
end

function test_case2()
    A = [ -0.24  -0.82   1.35   0.36   0.35
          -0.53  -0.20  -0.76   0.98  -0.54
           0.22   1.25  -1.60  -1.37  -1.94
          -0.51  -0.56  -0.08   0.96   0.46
           0.48  -2.25   0.38   0.06  -1.29 ]
    b = [-1.6, 0.19, 0.17, 0.31, -1.27]
    x = [2.2010416, 1.19009924, 0.0, 1.55001345, 0.0]
    return A, b, x
end

function test_case3() # non-float
    A = ones(Int, 4, 3)
    b = 2*ones(Int, 4)
    x = 2*ones(Int, 3)
    return A, b, x
end

function test_algorithm(fh::Function, ε::Real=1e-5)
    # Solve A*x = b for x, subject to x >=0
    A, b, x = test_case1()
    @test norm(fh(A,b) - x) < ε

    A, b, x = test_case2()
    @test norm(fh(A,b) - x) < ε

    # Test a bunch of random cases
    for i = 1:100
        m,n = rand(1:10),rand(1:10)
        A3 = randn(m,n)
        b3 = randn(m)
        x3,resid = pyopt.nnls(A3,b3)
        if resid > ε
            @test norm(fh(A3,b3) - x3) < ε
        else
            @test norm(A3*fh(A3,b3) - b3) < ε
        end
    end
end

nnls(A,b) = nonneg_lsq(A, b; alg=:nnls)
nnls_gram(A,b) = nonneg_lsq(A'*A, A'*b; alg=:nnls, gram=true)
fnnls(A,b) = nonneg_lsq(A, b; alg=:fnnls)
fnnls_gram(A,b) = nonneg_lsq(A'*A, A'*b; alg=:fnnls, gram=true)
pivot(A,b) = nonneg_lsq(A, b; alg=:pivot)
pivot_comb(A,b) = nonneg_lsq(A, b; alg=:pivot, variant=:comb)
pivot_cache(A,b) = nonneg_lsq(A, b; alg=:pivot, variant=:cache)

algs = [nnls, nnls_gram, fnnls, fnnls_gram, pivot, pivot_comb, pivot_cache]
errs = fill(1e-5, length(algs))

for (f, ε) in zip(algs, errs)
    print("testing ")
    @show f
    test_algorithm(f, ε)
    println("done")
end

#= non-float test fails, so revisit later
@testset "pivot_cache-non-float" begin
    A, b, x = test_case3()
    xi = pivot_cache(A, b)
    xf = pivot_cache(Float32.(A), Float32.(b))
    @test xi ≈ xf
end
=#

@testset "comb" begin
    A, b, x = test_case2()
    x0 = pivot_comb(A, b)
    P! = falses(size(A,2),1)
    x1 = nonneg_lsq(A, b; alg=:pivot, variant=:comb, P! = P!)
    @test x0 == x1
    @test P! != falses(size(A,2),1)
    @test all(x0[(!).(P!)] .== 0)
end

@testset "NNLS" begin include("nnls_test.jl") end
@testset "FNNLS" begin include("fnnls_test.jl") end
@testset "Pivot" begin include("pivot_test.jl") end
@testset "Sparse" begin include("sparse_test.jl") end


text:
#
# Qubit site type
#

# Define Qubit space in terms of
# Qubit/2 space, but use different
# defaults for QN names

"""
    space(::SiteType"Qubit";
          conserve_qns = false,
          conserve_parity = conserve_qns,
          conserve_number = false,
          qnname_parity = "Parity",
          qnname_number = "Number")

Create the Hilbert space for a site of type "Qubit".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"Qubit";
  conserve_qns=false,
  conserve_parity=conserve_qns,
  conserve_number=false,
  qnname_parity="Parity",
  qnname_number="Number",
)
  if conserve_number && conserve_parity
    return [
      QN((qnname_number, 0), (qnname_parity, 0, 2)) => 1,
      QN((qnname_number, 1), (qnname_parity, 1, 2)) => 1,
    ]
  elseif conserve_number
    return [QN(qnname_number, 0) => 1, QN(qnname_number, 1) => 1]
  elseif conserve_parity
    return [QN(qnname_parity, 0, 2) => 1, QN(qnname_parity, 1, 2) => 1]
  end
  return 2
end

ITensors.val(::ValName"0", ::SiteType"Qubit") = 1
ITensors.val(::ValName"1", ::SiteType"Qubit") = 2
ITensors.val(::ValName"Up", ::SiteType"Qubit") = 1
ITensors.val(::ValName"Dn", ::SiteType"Qubit") = 2
ITensors.val(::ValName"↑", ::SiteType"Qubit") = 1
ITensors.val(::ValName"↓", ::SiteType"Qubit") = 2

ITensors.state(::StateName"0", ::SiteType"Qubit") = [1.0, 0.0]
ITensors.state(::StateName"1", ::SiteType"Qubit") = [0.0, 1.0]
ITensors.state(::StateName"+", ::SiteType"Qubit") = [1.0, 1.0] / √2
ITensors.state(::StateName"-", ::SiteType"Qubit") = [1.0, -1.0] / √2
ITensors.state(::StateName"i", ::SiteType"Qubit") = [1.0, im] / √2
ITensors.state(::StateName"-i", ::SiteType"Qubit") = [1.0, -im] / √2
ITensors.state(::StateName"Up", t::SiteType"Qubit") = state(StateName("0"), t)
ITensors.state(::StateName"Dn", t::SiteType"Qubit") = state(StateName("1"), t)
ITensors.state(::StateName"↑", t::SiteType"Qubit") = state(StateName("0"), t)
ITensors.state(::StateName"↓", t::SiteType"Qubit") = state(StateName("1"), t)

# Pauli eingenstates
ITensors.state(::StateName"X+", t::SiteType"Qubit") = state(StateName("+"), t)
ITensors.state(::StateName"X-", t::SiteType"Qubit") = state(StateName("-"), t)

ITensors.state(::StateName"Y+", t::SiteType"Qubit") = state(StateName("i"), t)
ITensors.state(::StateName"Y-", t::SiteType"Qubit") = state(StateName("-i"), t)

ITensors.state(::StateName"Z+", t::SiteType"Qubit") = state(StateName("0"), t)
ITensors.state(::StateName"Z-", t::SiteType"Qubit") = state(StateName("1"), t)

# SIC-POVMs
state(::StateName"Tetra1", t::SiteType"Qubit") = state(StateName("Z+"), t)
state(::StateName"Tetra2", t::SiteType"Qubit") = [
  1 / √3
  √2 / √3
]
state(::StateName"Tetra3", t::SiteType"Qubit") = [
  1 / √3
  √2 / √3 * exp(im * 2π / 3)
]
state(::StateName"Tetra4", t::SiteType"Qubit") = [
  1 / √3
  √2 / √3 * exp(im * 4π / 3)
]

#
# 1-Qubit gates
#
ITensors.op(::OpName"X", ::SiteType"Qubit") = [
  0 1
  1 0
]

ITensors.op(::OpName"σx", t::SiteType"Qubit") = op("X", t)

ITensors.op(::OpName"σ1", t::SiteType"Qubit") = op("X", t)

ITensors.op(::OpName"Y", ::SiteType"Qubit") = [
  0.0 -1.0im
  1.0im 0.0
]

ITensors.op(::OpName"σy", t::SiteType"Qubit") = op("Y", t)

ITensors.op(::OpName"σ2", t::SiteType"Qubit") = op("Y", t)

ITensors.op(::OpName"iY", ::SiteType"Qubit") = [
  0 1
  -1 0
]
ITensors.op(::OpName"iσy", t::SiteType"Qubit") = op("iY", t)

ITensors.op(::OpName"iσ2", t::SiteType"Qubit") = op("iY", t)

ITensors.op(::OpName"Z", ::SiteType"Qubit") = [
  1 0
  0 -1
]

ITensors.op(::OpName"σz", t::SiteType"Qubit") = op("Z", t)

ITensors.op(::OpName"σ3", t::SiteType"Qubit") = op("Z", t)

function ITensors.op(::OpName"√NOT", ::SiteType"Qubit")
  return [
    (1 + im)/2 (1 - im)/2
    (1 - im)/2 (1 + im)/2
  ]
end

ITensors.op(::OpName"√X", t::SiteType"Qubit") = op("√NOT", t)

ITensors.op(::OpName"H", ::SiteType"Qubit") = [
  1/sqrt(2) 1/sqrt(2)
  1/sqrt(2) -1/sqrt(2)
]

# Rϕ with ϕ = π/2
ITensors.op(::OpName"Phase", ::SiteType"Qubit"; ϕ::Number=π / 2) = [
  1 0
  0 exp(im * ϕ)
]

ITensors.op(::OpName"P", t::SiteType"Qubit"; kwargs...) = op("Phase", t; kwargs...)

ITensors.op(::OpName"S", t::SiteType"Qubit") = op("Phase", t; ϕ=π / 2)

## Rϕ with ϕ = π/4
ITensors.op(::OpName"π/8", ::SiteType"Qubit") = [
  1 0
  0 1 / sqrt(2)+im / sqrt(2)
]

ITensors.op(::OpName"T", t::SiteType"Qubit") = op("π/8", t)

# Rotation around X-axis
function ITensors.op(::OpName"Rx", ::SiteType"Qubit"; θ::Number)
  return [
    cos(θ / 2) -im*sin(θ / 2)
    -im*sin(θ / 2) cos(θ / 2)
  ]
end

# Rotation around Y-axis
function ITensors.op(::OpName"Ry", ::SiteType"Qubit"; θ::Number)
  return [
    cos(θ / 2) -sin(θ / 2)
    sin(θ / 2) cos(θ / 2)
  ]
end

# Rotation around Z-axis
function ITensors.op(::OpName"Rz", ::SiteType"Qubit"; ϕ::Number)
  return [
    exp(-im * ϕ / 2) 0
    0 exp(im * ϕ / 2)
  ]
end

# Rotation around generic axis n̂
function ITensors.op(::OpName"Rn", ::SiteType"Qubit"; θ::Real, ϕ::Real, λ::Real)
  return [
    cos(θ / 2) -exp(im * λ)*sin(θ / 2)
    exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2)
  ]
end

function ITensors.op(::OpName"Rn̂", t::SiteType"Qubit"; kwargs...)
  return ITensors.op(OpName("Rn"), t; kwargs...)
end

#
# 2-Qubit gates
#

ITensors.op(::OpName"CNOT", ::SiteType"Qubit") = [
  1 0 0 0
  0 1 0 0
  0 0 0 1
  0 0 1 0
]

ITensors.op(::OpName"CX", t::SiteType"Qubit") = op("CNOT", t)

ITensors.op(::OpName"CY", ::SiteType"Qubit") = [
  1 0 0 0
  0 1 0 0
  0 0 0 -im
  0 0 im 0
]

ITensors.op(::OpName"CZ", ::SiteType"Qubit") = [
  1 0 0 0
  0 1 0 0
  0 0 1 0
  0 0 0 -1
]

function ITensors.op(::OpName"CPHASE", ::SiteType"Qubit"; ϕ::Number)
  return [
    1 0 0 0
    0 1 0 0
    0 0 1 0
    0 0 0 exp(im * ϕ)
  ]
end
ITensors.op(::OpName"Cphase", t::SiteType"Qubit"; kwargs...) = op("CPHASE", t; kwargs...)

function ITensors.op(::OpName"CRx", ::SiteType"Qubit"; θ::Number)
  return [
    1 0 0 0
    0 1 0 0
    0 0 cos(θ / 2) -im*sin(θ / 2)
    0 0 -im*sin(θ / 2) cos(θ / 2)
  ]
end
ITensors.op(::OpName"CRX", t::SiteType"Qubit"; kwargs...) = ITensors.op("CRx", t; kwargs...)

function ITensors.op(::OpName"CRy", ::SiteType"Qubit"; θ::Number)
  return [
    1 0 0 0
    0 1 0 0
    0 0 cos(θ / 2) -sin(θ / 2)
    0 0 sin(θ / 2) cos(θ / 2)
  ]
end
ITensors.op(::OpName"CRY", t::SiteType"Qubit"; kwargs...) = ITensors.op("CRy", t; kwargs...)

function ITensors.op(::OpName"CRz", ::SiteType"Qubit"; ϕ::Number)
  return [
    1 0 0 0
    0 1 0 0
    0 0 exp(-im * ϕ / 2) 0
    0 0 0 exp(im * ϕ / 2)
  ]
end
ITensors.op(::OpName"CRZ", t::SiteType"Qubit"; kwargs...) = ITensors.op("CRz", t; kwargs...)

function ITensors.op(::OpName"CRn", ::SiteType"Qubit"; θ::Number, ϕ::Number, λ::Number)
  return [
    1 0 0 0
    0 1 0 0
    0 0 cos(θ / 2) -exp(im * λ)*sin(θ / 2)
    0 0 exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2)
  ]
end
function ITensors.op(::OpName"CRn̂", t::SiteType"Qubit"; kwargs...)
  return ITensors.op("CRn", t; kwargs...)
end

ITensors.op(::OpName"SWAP", ::SiteType"Qubit") = [
  1 0 0 0
  0 0 1 0
  0 1 0 0
  0 0 0 1
]
ITensors.op(::OpName"Swap", t::SiteType"Qubit") = op("SWAP", t)

function ITensors.op(::OpName"√SWAP", ::SiteType"Qubit")
  return [
    1 0 0 0
    0 (1 + im)/2 (1 - im)/2 0
    0 (1 - im)/2 (1 + im)/2 0
    0 0 0 1
  ]
end
ITensors.op(::OpName"√Swap", t::SiteType"Qubit") = op("√SWAP", t)

ITensors.op(::OpName"iSWAP", t::SiteType"Qubit") = [
  1 0 0 0
  0 0 im 0
  0 im 0 0
  0 0 0 1
]
ITensors.op(::OpName"iSwap", t::SiteType"Qubit") = op("iSWAP", t)

function ITensors.op(::OpName"√iSWAP", t::SiteType"Qubit")
  return [
    1 0 0 0
    0 1/√2 im/√2 0
    0 im/2 1/√2 0
    0 0 0 1
  ]
end
ITensors.op(::OpName"√iSwap", t::SiteType"Qubit") = op("√iSWAP", t)

# Ising (XX) coupling gate
function ITensors.op(::OpName"Rxx", t::SiteType"Qubit"; ϕ::Number)
  return [
    cos(ϕ) 0 0 -im*sin(ϕ)
    0 cos(ϕ) -im*sin(ϕ) 0
    0 -im*sin(ϕ) cos(ϕ) 0
    -im*sin(ϕ) 0 0 cos(ϕ)
  ]
end
ITensors.op(::OpName"RXX", t::SiteType"Qubit"; kwargs...) = op("Rxx", t; kwargs...)

# Ising (YY) coupling gate
function ITensors.op(::OpName"Ryy", ::SiteType"Qubit"; ϕ::Number)
  return [
    cos(ϕ) 0 0 im*sin(ϕ)
    0 cos(ϕ) -im*sin(ϕ) 0
    0 -im*sin(ϕ) cos(ϕ) 0
    im*sin(ϕ) 0 0 cos(ϕ)
  ]
end
ITensors.op(::OpName"RYY", t::SiteType"Qubit"; kwargs...) = op("Ryy", t; kwargs...)

# Ising (ZZ) coupling gate
function ITensors.op(::OpName"Rzz", ::SiteType"Qubit"; ϕ::Number)
  return [
    exp(-im * ϕ) 0 0 0
    0 exp(im * ϕ) 0 0
    0 0 exp(im * ϕ) 0
    0 0 0 exp(-im * ϕ)
  ]
end
ITensors.op(::OpName"RZZ", t::SiteType"Qubit"; kwargs...) = op("Rzz", t; kwargs...)

#
# 3-Qubit gates
#

function ITensors.op(::OpName"Toffoli", ::SiteType"Qubit")
  return [
    1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0
    0 0 1 0 0 0 0 0
    0 0 0 1 0 0 0 0
    0 0 0 0 1 0 0 0
    0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 1
    0 0 0 0 0 0 1 0
  ]
end

ITensors.op(::OpName"CCNOT", t::SiteType"Qubit") = op("Toffoli", t)

ITensors.op(::OpName"CCX", t::SiteType"Qubit") = op("Toffoli", t)

ITensors.op(::OpName"TOFF", t::SiteType"Qubit") = op("Toffoli", t)

function ITensors.op(::OpName"Fredkin", ::SiteType"Qubit")
  return [
    1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0
    0 0 1 0 0 0 0 0
    0 0 0 1 0 0 0 0
    0 0 0 0 1 0 0 0
    0 0 0 0 0 0 1 0
    0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 1
  ]
end

ITensors.op(::OpName"CSWAP", t::SiteType"Qubit") = op("Fredkin", t)
ITensors.op(::OpName"CSwap", t::SiteType"Qubit") = op("Fredkin", t)

ITensors.op(::OpName"CS", t::SiteType"Qubit") = op("Fredkin", t)

#
# 4-Qubit gates
#

function ITensors.op(::OpName"CCCNOT", ::SiteType"Qubit")
  return [
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
  ]
end

# spin-full operators
ITensors.op(::OpName"Sz", ::SiteType"Qubit") = [
  0.5 0.0
  0.0 -0.5
]

ITensors.op(::OpName"Sᶻ", t::SiteType"Qubit") = op(OpName("Sz"), t)

ITensors.op(::OpName"S+", ::SiteType"Qubit") = [
  0 1
  0 0
]

ITensors.op(::OpName"S⁺", t::SiteType"Qubit") = op(OpName("S+"), t)

ITensors.op(::OpName"Splus", t::SiteType"Qubit") = op(OpName("S+"), t)

ITensors.op(::OpName"S-", ::SiteType"Qubit") = [
  0 0
  1 0
]

ITensors.op(::OpName"S⁻", t::SiteType"Qubit") = op(OpName("S-"), t)

ITensors.op(::OpName"Sminus", t::SiteType"Qubit") = op(OpName("S-"), t)

ITensors.op(::OpName"Sx", ::SiteType"Qubit") = [
  0.0 0.5
  0.5 0.0
]

ITensors.op(::OpName"Sˣ", t::SiteType"Qubit") = op(OpName("Sx"), t)

ITensors.op(::OpName"iSy", ::SiteType"Qubit") = [
  0.0 0.5
  -0.5 0.0
]

ITensors.op(::OpName"iSʸ", t::SiteType"Qubit") = op(OpName("iSy"), t)

ITensors.op(::OpName"Sy", ::SiteType"Qubit") = [
  0.0 -0.5im
  0.5im 0.0
]

ITensors.op(::OpName"Sʸ", t::SiteType"Qubit") = op(OpName("Sy"), t)

ITensors.op(::OpName"S2", ::SiteType"Qubit") = [
  0.75 0.0
  0.0 0.75
]

ITensors.op(::OpName"S²", t::SiteType"Qubit") = op(OpName("S2"), t)

ITensors.op(::OpName"ProjUp", ::SiteType"Qubit") = [
  1 0
  0 0
]

ITensors.op(::OpName"projUp", t::SiteType"Qubit") = op(OpName("ProjUp"), t)

ITensors.op(::OpName"ProjDn", ::SiteType"Qubit") = [
  0 0
  0 1
]

ITensors.op(::OpName"projDn", t::SiteType"Qubit") = op(OpName("ProjDn"), t)


text:
using Soss

p = @model r,s begin
    p_bad ~ Beta(1,3) |> iid(s)
    bad ~ For(s) do j
            Bernoulli(p_bad[j])
        end

    fpr ~ Beta(1,3) |> iid(r)
    tpr ~ Beta(3,1) |> iid(r)

    pos_rate = hcat(fpr, tpr)

    y ~ For(r,s) do i,j
        Bernoulli(pos_rate[i, bad[j] + 1])
    end
    return y
end;

markovBlanket(p, :bad)

+

# ASVI
q = @model r,s,λα begin
    λ = λα.λ
    α = λα.α
    p_bad ~ For(s) do j
        a = 1 * λ.p_bad[j,1] + α.p_bad[j,1] * (1 - λ.p_bad[j,1])
        b = 3 * λ.p_bad[j,2] + α.p_bad[j,2] * (1 - λ.p_bad[j,2])
        Beta(a,b)
    end
    bad ~ For(s) do j
            Bernoulli(p_bad[j] * λ.bad[j] + α.bad[j] * (1 - λ.bad[j]))
        end
    tpr ~ For(r) do i
            a = 3 * λ.tpr[i,1] + α.tpr[i,1] * (1 - λ.tpr[i,1])
            b = 1 * λ.tpr[i,2] + α.tpr[i,2] * (1 - λ.tpr[i,2])
            Beta(a,b)
        end
    fpr ~ For(r) do i
        a = 1 * λ.fpr[i,1] + α.fpr[i,1] * (1 - λ.fpr[i,1])
        b = 3 * λ.fpr[i,2] + α.fpr[i,2] * (1 - λ.fpr[i,2])
        Beta(a,b)
    end
end

tr = as((
    λ=as((
        p_bad = as(Array, as𝕀, 14, 2),
        bad = as(Array, as𝕀, 14),
        tpr = as(Array, as𝕀, 5, 2),
        fpr = as(Array, as𝕀, 5, 2))
        ),
    α = as((
        p_bad = as(Array, asℝ₊, 14, 2),
        bad = as(Array, as𝕀, 14),
        tpr = as(Array, asℝ₊, 5, 2),
        fpr = as(Array, asℝ₊, 5, 2))
        )
    ))

x = zeros(124)
λα = tr(x)

y = rand(p(r=5, s=14));
posterior = p(r=5, s=14) | (;y)

function elbo(x)
    λ = tr(x)

    qλ = q(r=5,s=14, λα=λα)
    r = rand(qλ)
    logdensity_def(posterior, r) - logdensity_def(qλ, r)
end

elbo(x)

using Random
Random.seed!(1);

y_obs = truth.y

obs = (y=y,)

particles(p(r=5,s=14))

entropy(p(r=5,s=14))

lik = predictive(p, :p_bad, :fpr, :tpr)

import TransformVariables
_lik = as((
           p_bad = as(Array, as𝕀, 14),
           fpr = as(Array, as𝕀, 5),
           tpr = as(Array, as𝕀, 5)))

ϕ = Optim.minimizer(result)
t_lik(ϕ)

q = @model r,s,λ begin
    p_bad ~ For(j -> Beta(1+λ.p_bad[j].α, 1+λ.p_bad[j].β), s)

    bad ~ For(j -> Bernoulli(p_bad[j]), s)

    fpr ~ For(i -> Beta(1+λ.fpr[i].α, 1+λ.fpr[i].β), r)
    tpr ~ For(i -> Beta(1+λ.tpr[i].α, 1+λ.tpr[i].β), r)
end;

import TransformVariables: as, asℝ₊
trλ = as((λ=as((
    p_bad = as(Array, as((α=asℝ₊, β=asℝ₊)), 14),
    fpr =   as(Array, as((α=asℝ₊, β=asℝ₊)), 5),
    tpr =   as(Array, as((α=asℝ₊, β=asℝ₊)), 5),
)),
))

trλ.dimension

x0 = zeros(48)
λ = trλ(x0)

logdensity_def(p(r=5, s=14) | (;y), rand(q(r=5,s=14, λ=λ.λ)))

function elbo(p, args, obs, q, λ)
    qargs = merge(args, λ)
    mean(logdensityof(p(args), merge(particles(q(qargs)), obs))) + entropy(q(args))
end

ϕ0 = zeros(48);
elbo(p,(r=5,s=14), obs, q, trλ(ϕ0))

using Optim
result = optimize(
      ϕ -> -elbo(p,(r=5,s=14),obs, q, trλ(ϕ))
    , ϕ0
    , Optim.Options(
        show_trace=true,
        show_every=100,
        iterations=1000
    )

)

ϕ0 = Optim.minimizer(result)
λ = trλ(ϕ0).λ

particles(q(r=5,s=14, λ=λ)) |> pairs

elbo(p,q,(r=5,s=14),trλ(result))

###################
# ADVI

using AdvancedVI
# using DistsAD

bad = rand(Bernoulli(0.5) |> iid(14))

tr = as(p(r=5,s=14), (y=y,));
d = tr.dimension

getq(λ) = MvNormal(λ[1:d], exp.(λ[d .+ (1:d)]))

advi = ADVI(10, 10_000)

function logπ(θ)
    logdensityof(p(r=5,s=14),  merge((y=y,), tr(θ)))
end

logπ(randn(24))

λ = vi(logπ, advi, getq, zeros(48)) |> particles |> tr

pairs(λ)

q = @model r,s,λ begin
    p_bad ~ Uniform() |> iid(s)

    bad ~ For(s) do j
        Bernoulli(λ.bad[j].p)
    end

    fpr ~ For(r) do i Beta(λ.fpr[i].α, λ.fpr[i].β, check_args=false) end
    tpr ~ For(r) do i Beta(λ.tpr[i].α, λ.tpr[i].β, check_args=false) end
end;

function logweightedmean(ℓ1,x1,ℓ2,x2)
    if isinf(ℓ1)
        return (x1+x2)/2
    else
        w = exp(ℓ2 - ℓ1)/2
        target = (x1 + w*x2)/(1 + w)
        return 0.7*x1 + 0.3*target
    end
end

λ = ( bad = [(p=0.5,) for j in 1:14]
    , fpr   = [(α=1.0, β=1.0) for i in 1:5]
    , tpr   = [(α=1.0, β=1.0) for i in 1:5]
)

importanceSample(p(r=5,s=14), q(r=5,s=14, λ=λ), (y=truth.y,))

function f(λ, N=100, warmup=10)
    λ0 = deepcopy(λ)

    for iternum in 1:N

        if iternum < warmup
            ℓx = importanceSample(p(r=5,s=14), q(r=5,s=14, λ=λ0), (y=truth.y,))
        else
            ℓx = importanceSample(p(r=5,s=14), q(r=5,s=14, λ=λ), (y=truth.y,))
        end

        x = ℓx.val
        ℓ = ℓx.ℓ

        for j in eachindex(λ.bad)
            bad = x.bad[j]
            p = λ.bad[j].p

            p = logweightedmean(λ._ℓ, p, ℓ, bad)
            λ = @set λ.bad[j].p = p
        end

        for j in eachindex(λ.p_bad)
            p_bad = x.p_bad[j]

            α = λ.p_bad[j].α
            β = λ.p_bad[j].β
            α_plus_β = α + β

            Ex = α / α_plus_β
            Ex² = α * β / (α_plus_β^2 * (α_plus_β + 1)) + Ex^2

            # Update moments
            Ex = logweightedmean(λ._ℓ, Ex, ℓ, p_bad)
            Ex² = logweightedmean(λ._ℓ, Ex², ℓ, p_bad^2)

            # Parameter estimates by method of moments
            V = Ex² - Ex^2
            n = Ex * (1 - Ex) / V - 1
            α = max(floatmin(Float32), n * Ex)
            β = max(floatmin(Float32), n * (1 - Ex))

            λ = @set λ.p_bad[j].α = α
            λ = @set λ.p_bad[j].β = β
        end

        for i in eachindex(λ.tpr)
            tpr = x.tpr[i]

            α = λ.tpr[i].α
            β = λ.tpr[i].β
            α_plus_β = α + β

            Ex = α / α_plus_β
            Ex² = α * β / (α_plus_β^2 * (α_plus_β + 1)) + Ex^2

            # Update moments
            Ex = logweightedmean(λ._ℓ, Ex, ℓ, tpr)
            Ex² = logweightedmean(λ._ℓ, Ex², ℓ, tpr^2)

            # Parameter estimates by method of moments
            V = Ex² - Ex^2
            n = Ex * (1 - Ex) / V - 1
            α = max(floatmin(Float32), n * Ex)
            β = max(floatmin(Float32), n * (1 - Ex))

            λ = @set λ.tpr[i].α = α
            λ = @set λ.tpr[i].β = β
        end

        for i in eachindex(λ.fpr)
            fpr = x.fpr[i]

            α = λ.fpr[i].α
            β = λ.fpr[i].β
            α_plus_β = α + β

            Ex = α / α_plus_β
            Ex² = α * β / (α_plus_β^2 * (α_plus_β + 1)) + Ex^2

            # Update moments
            Ex = logweightedmean(λ._ℓ, Ex, ℓ, fpr)
            Ex² = logweightedmean(λ._ℓ, Ex², ℓ, fpr^2)

            # Parameter estimates by method of moments
            V = Ex² - Ex^2
            n = Ex * (1 - Ex) / V - 1
            α = max(floatmin(Float32), n * Ex)
            β = max(floatmin(Float32), n * (1 - Ex))

            λ = @set λ.fpr[i].α = α
            λ = @set λ.fpr[i].β = β
        end

        λ = @set λ._ℓ = logaddexp(λ._ℓ, ℓ)
    end
    return λ
end

λ = ( _ℓ = -Inf
    , p_bad = [(α=1.0, β=1.0) for j in 1:14]
    , bad   = [(p=0.5,) for j in 1:14]
    , fpr   = [(α=1.0, β=1.0) for i in 1:5]
    , tpr   = [(α=1.0, β=1.0) for i in 1:5]
);

λ = f(λ,20000,100);
λ.bad

sourceParticles()(m, Val{1000})

particles(m(r=5,s=14))
post=dynamicHMC(m2(r=5,s=14), (y=truth.y,))

dynamicHMC(m2(r=5,s=14), (y=truth.y,))

markovBlanket(m, :bad)

m = @model s begin
    p_bad ~ Beta(1,3) |> iid(s)
    bad ~ Normal(p_bad[1],1) |> iid(s)
    end

particles(m(s=3))


text:
using JSON
using InteractiveUtils: subtypes

function savefb(filename::AbstractString, nsolt::AbstractNsolt, mode::AbstractString="w")
    if match(r".*\.json$", filename) === nothing
        filename = string(filename, ".json")
    end
    open(io->savefb(io, nsolt), filename, mode)
end

savefb(io::IOStream, nsolt::AbstractNsolt) = print(io, serialize(nsolt))

function serialize(nsolt::NS) where {NS<:AbstractNsolt}
    # convert to Dict
    dict_nsolt = nsolt |> JSON.json |> JSON.parse

    nsolt_configs, nsolt_params = separate_nsolt_configs_params(NS, dict_nsolt)

    delete!(nsolt_configs, "nStages")
    push!(nsolt_configs, "polyphaseOrder" => collect(orders(nsolt)))

    if NS <: Cnsolt
        nsolt_params["Φ"] = map(diag(nsolt.Φ)) do p
            Dict( "re" => real(p), "im" => imag(p))
        end
    end

    outdict = Dict(
        "NsoltType" => string(nameof(NS)),
        "ElementType" => string(eltype(NS)),
        "Configs" => nsolt_configs,
        "Params" => nsolt_params,
    )
    JSON.json(outdict)
end

function separate_nsolt_configs_params(::Type{NS}, dict::Dict) where {NS<:AbstractNsolt}
    keynames = [ "decimationFactor", "polyphaseOrder", "nChannels", "perm" ]
    configs = filter(dict) do d
        any(d.first .== keynames)
    end
    params = filter(dict) do d
        !any(d.first .== keynames)
    end

    (configs, params)
end

loadfb(filename::AbstractString) = open(io->MDCDL.loadfb(io), filename)

function loadfb(io::IOStream)
    dic = JSON.parse(read(io, String))

    deserialize(dic)
end

function deserialize(dict::Dict;
        nsolttypeset=[ subtypes(Rnsolt)..., subtypes(Cnsolt)... ],
        eltypeset=subtypes(AbstractFloat)
    )

    strNS = dict["NsoltType"]
    strT = dict["ElementType"]
    dic_nsolt_configs = dict["Configs"]
    dic_nsolt_params = dict["Params"]

    NS = Dict( string.(nsolttypeset) .=> nsolttypeset )[strNS]
    T = Dict( string.(eltypeset) .=> eltypeset )[strT]

    nsolt = create_nsolt_by_dict(NS, T, dic_nsolt_configs)

    return set_nsolt_params_by_dict!(nsolt, dic_nsolt_params)
end

function create_nsolt_by_dict(::Type{NS}, T::Type, dict::Dict) where {NS<:AbstractNsolt}
    df = get_nsolt_config(NS, dict, "decimationFactor")
    ord = get_nsolt_config(NS, dict, "polyphaseOrder")
    nch = get_nsolt_config(NS, dict, "nChannels")
    perm = if haskey(dict, "perm")
        get_nsolt_config(NS, dict, "perm")
    else
        (collect(1:length(df))...,)
    end

    NS(T, df, ord, nch, perm=perm)
end

function set_nsolt_params_by_dict!(nsolt::NS, dict::Dict) where {NS<:RnsoltTypeI}
    initpm = [ nsolt.W0 => "W0", nsolt.U0 => "U0", nsolt.CJ => "CJ"]
    proppm = [ nsolt.Udks => "Udks" ]
    set_nsolt_params_by_dict!(NS, dict, initpm, proppm)

    return nsolt
end

function set_nsolt_params_by_dict!(nsolt::NS, dict::Dict) where {NS<:RnsoltTypeII}
    initpm = [ nsolt.W0 => "W0", nsolt.U0 => "U0", nsolt.CJ => "CJ"]
    proppm = [ nsolt.Wdks => "Wdks", nsolt.Udks => "Udks" ]
    set_nsolt_params_by_dict!(NS, dict, initpm, proppm)

    return nsolt
end

function set_nsolt_params_by_dict!(nsolt::NS, dict::Dict) where {NS<:CnsoltTypeI}
    initpm = [ nsolt.V0 => "V0", nsolt.FJ => "FJ"]
    proppm = [ nsolt.Wdks => "Wdks", nsolt.Udks => "Udks", nsolt.θdks => "θdks" ]
    set_nsolt_params_by_dict!(NS, dict, initpm, proppm)

    return nsolt
end

function set_nsolt_params_by_dict!(nsolt::NS, dict::Dict) where {NS<:CnsoltTypeII}
    initpm = [ nsolt.V0 => "V0", nsolt.FJ => "FJ"]
    proppm = [ nsolt.Wdks => "Wdks", nsolt.Udks => "Udks", nsolt.θ1dks => "θ1dks", nsolt.Ŵdks => "Ŵdks", nsolt.Ûdks => "Ûdks", nsolt.θ2dks => "θ2dks" ]
    set_nsolt_params_by_dict!(NS, dict, initpm, proppm)

    return nsolt
end

function set_nsolt_params_by_dict!(NS::Type, dict::Dict, initpm::Array, proppm::Array)
    foreach(initpm) do pm
        pm.first .= get_nsolt_config(NS, dict, pm.second)
    end

    foreach(proppm) do pm
        a = get_nsolt_config(NS, dict, pm.second)
        setparams_dk!(pm.first, a)
    end
end

function setparams_dk!(dst, src)
    foreach(dst, src) do dstk, srck
        foreach(dstk, srck) do dstkd, srckd
            dstkd .= srckd
        end
    end
end

get_nsolt_config(::Type{NS}, dict::Dict, akey::AbstractString) where {NS} = get_nsolt_config(NS, Val(Symbol(akey)), dict[akey])

function get_nsolt_config(::Type{NS}, ::Val{:decimationFactor}, data) where {NS<:AbstractNsolt}
    (Int.(data)...,)
end

function get_nsolt_config(::Type{NS}, ::Val{:polyphaseOrder}, data) where {NS<:AbstractNsolt}
    (Int.(data)...,)
end

function get_nsolt_config(::Type{NS}, ::Val{:nChannels}, data) where {NS<:Rnsolt}
    (Int.(data)...,)
end

function get_nsolt_config(::Type{NS}, ::Val{:perm}, data) where {NS<:AbstractNsolt}
    (Int.(data)...,)
end

function get_nsolt_config(::Type{NS}, ::Val{:nChannels}, data) where {NS<:Cnsolt}
    Int(sum(data))
end

function get_nsolt_config(::Type{NS}, ::Val{:CJ}, data) where {T,NS<:Rnsolt{T}}
    T.(hcat(data...))
end

function get_nsolt_config(::Type{NS}, ::Val{:FJ}, data) where {T,NS<:Cnsolt{T}}
    map(hcat(data...)) do FJpm
        Complex{T}(FJpm["re"] + FJpm["im"]*im)
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:W0}, data) where {T,NS<:Rnsolt{T}}
    T.(hcat(data...))
end

function get_nsolt_config(::Type{NS}, ::Val{:U0}, data) where {T,NS<:Rnsolt{T}}
    T.(hcat(data...))
end

function get_nsolt_config(::Type{NS}, ::Val{:V0}, data) where {T,NS<:Cnsolt{T}}
    T.(hcat(data...))
end

function get_nsolt_config(::Type{NS}, ::Val{:Udks}, data) where {T,NS<:AbstractNsolt{T}}
    map(data) do Ud
        map(Ud) do Udk
            T.(hcat(Udk...))
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:Wdks}, data) where {T,NS<:RnsoltTypeII{T}}
    map(data) do Wd
        map(Wd) do Wdk
            T.(hcat(Wdk...))
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:Wdks}, data) where {T,NS<:Cnsolt{T}}
    map(data) do Wd
        map(Wd) do Wdk
            T.(hcat(Wdk...))
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:θdks}, data) where {T,NS<:CnsoltTypeI{T}}
    map(data) do θd
        map(θd) do θdk
            T.(θdk)
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:θ1dks}, data) where {T,NS<:CnsoltTypeII{T}}
    map(data) do θd
        map(θd) do θdk
            T.(θdk)
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:Ŵdks}, data) where {T,NS<:CnsoltTypeII{T}}
    map(data) do Wd
        map(Wd) do Wdk
            T.(hcat(Wdk...))
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:Ûdks}, data) where {T,NS<:CnsoltTypeII{T}}
    map(data) do Wd
        map(Wd) do Wdk
            T.(hcat(Wdk...))
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:θ2dks}, data) where {T,NS<:CnsoltTypeII{T}}
    map(data) do θd
        map(θd) do θdk
            T.(θdk)
        end
    end
end

function get_nsolt_config(::Type{NS}, ::Val{:Φ}, data) where {T,NS<:Cnsolt{T}}
    Diagonal(map(p->Complex{T}(p["re"] + p["im"]*im), data))
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
module convnet

using mnist

# tipuri
typealias Neurons{T} Dict{Int, Array{T}}
typealias Errors{T}  Dict{Int, Array{T}}
typealias Weights{T} Dict{Int, Array{T}}
typealias Gradients{T} Dict{Int, Array{T}}

typealias NetworkConfig Array{Dict{String,Any}, 1}

type FwdComputation{T}
  batchSize::Int64
  A::Neurons{T}
  Z::Neurons{T}
end

type BackComputation{T}
  batchSize::Int64
  E::Errors{T}
  G::Gradients{T}
end

function init_weights(T::DataType, config::NetworkConfig,
                      minValue::Float64 = -0.5, maxValue::Float64 = 0.5)
  # Dictionary of weights
  weights = Weights{T}();
  # Size of inputs and number of maps from previous layer
  width = config[1]["width"];
  maps_no = 1;
  # Init weights
  for l in 2:(length(config))
    layer = config[l];
    if layer["type"] == "CONV"
      weights[l] = rand(T, layer["maps"], maps_no, layer["r"] * layer["r"] + 1);
      weights[l] = minValue + weights[l] * maxValue;
      maps_no = layer["maps"];
      width = width - layer["r"] + 1;
    elseif layer["type"] == "MAXPOOL"
      width = width / layer["r"];
    elseif layer["type"] == "FULL"
      weights[l] = rand(layer["m"], convert(Int64,maps_no * width * width + 1));
      weights[l] = minValue + weights[l] * maxValue;
      maps_no = layer["m"];
      width = 1;
    end
  end
  return weights
end

function prepare_network(T::DataType, config::NetworkConfig, batchSize)
  A = Neurons{T}();
  Z = Neurons{T}();
  E = Neurons{T}();
  G = init_weights(T, config, 0.0, 0.0);
  # Size of inputs and number of maps from previous layer
  width = config[1]["width"];
  maps_no = 1;
  # Init weights
  for l in 2:(length(config))
    layer = config[l];
    if layer["type"] == "CONV"
      width = width - layer["r"] + 1;
      maps_no = layer["maps"];
      A[l] = zeros(T, batchSize, maps_no, width, width);
      Z[l] = zeros(T, batchSize, maps_no, width, width);
      E[l] = zeros(T, batchSize, maps_no, width, width);
    elseif layer["type"] == "MAXPOOL"
      width = convert(Int64, width / layer["r"]);
      Z[l] = zeros(T, batchSize, maps_no, width, width);
      E[l] = zeros(T, batchSize, maps_no, width, width);
    elseif layer["type"] == "FULL"
      maps_no = layer["m"];
      width = 1;
      A[l] = zeros(T, batchSize, maps_no);
      Z[l] = zeros(T, batchSize, maps_no);
      E[l] = zeros(T, batchSize, maps_no);
    end
  end
  return FwdComputation(batchSize, A, Z), BackComputation(batchSize, E, G);
end

function describe_weights{T}(W::Weights{T})
  for w in W
    print("Weights on layer #", w[1], " have size ", size(w[2]));
    println(" (", prod(size(w[2])), " parameters)");
  end
end

layers = Dict{String, Any}[{"type" => "INPUTS", "width" => 28},
                           {"type" => "CONV", "r" => 5, "maps" => 6,
                            "mapping" => Bool[1 1 1 1 1 1],
                            "activation" => "sigmoid"},
                           {"type" => "MAXPOOL", "r" => 2},
                           {"type" => "CONV", "r" => 5, "maps" => 13,
                            "mapping" => Bool[1 1 0 1 1 1 1 1 0 0 1 1 1;
                                              1 1 1 0 1 1 1 1 1 0 0 1 1;
                                              1 1 1 1 0 1 1 1 1 1 0 0 1;
                                              1 1 1 1 1 0 1 1 1 1 1 0 0;
                                              1 1 1 1 1 1 0 0 1 1 1 1 0;
                                              1 0 1 1 1 1 1 0 0 1 1 1 1],
                            "activation" => "sigmoid"},
                           {"type" => "MAXPOOL", "r" => 2},
                           {"type" => "FULL", "m" => 400, "activation" => "sigmoid"},
                           {"type" => "FULL", "m" => 10, "activation" => "softmax"}
                           ];

sigmoid(x) = 1 ./ (1 .+ exp(-x));

function softmax{T}(a::Array{T,1})
  y = 1 ./ 1 .+ exp(a);
  y = y / sum(y);
end

function forward{T}(config::NetworkConfig, X::Array{T,2},
                    W::Weights{T}, nn::FwdComputation{T})
  maps_no = 1;
  width = config[1]["width"];

  nn.Z[1] = reshape(X, nn.batchSize, 1, width, width);

  for l in 2:(length(config)-1)

    println("Gigi va da bani: ", l);

    layer = config[l];

    if layer["type"] == "CONV"
      println("Gigi ", l, " este convolutiv si recalcitrant");
      new_width = width - layer["r"] + 1;
      new_maps_no = layer["maps"];

      r = layer["r"];
      for n in 1:nn.batchSize
        for m in 1:new_maps_no
          for i in 1:new_width
            for j in 1:new_width
              nn.A[l][n, m, i, j] = zero(T);
              for prev_m in 1:maps_no
                input = [one(T) reshape(nn.Z[l-1][n, prev_m, i:(i+r-1), j:(j+r-1)], 1, r * r)];
                nn.A[l][n, m, i, j] += layer["mapping"][prev_m,m] * (input * reshape(W[l][m,prev_m,:], r*r+1,1))[1];
              end
            end
          end
        end
      end
      nn.Z[l] = sigmoid(nn.A[l]);

      width = new_width;
      maps_no = new_maps_no;
    elseif layer["type"] == "MAXPOOL"
      r = layer["r"];
      new_width = convert(Int64, width / r);
      for n in 1:nn.batchSize
        for (m in 1:maps_no)
          for (i in 1:new_width)
            for (j in 1:new_width)
              nn.Z[l][n, m, i, j] =
                maximum(nn.Z[l-1][n, m, (r*(i-1)+1):(r*i), (r*(i-1)+1):(r*i)]);
            end
          end
        end
      end
      width = new_width;
    elseif layer["type"] == "FULL"
      for (n in 1:nn.batchSize)
        for (j in 1:layer["m"])
          nn.A[l][n,j] =
            ([1 reshape(nn.Z[l-1][n,1:end,1:end,1:end], 1, width * maps_no * width)]
             *
             vec(W[l][j,:]))[1];
        end
      end
      if layer["activation"] == "sigmoid"
        nn.Z[l] = sigmoid(nn.A[l]);
      elseif layer["activation"] == "softmax"
        nn.Z[l] = softmax(nn.A[l]);
      end
      width = 1;
      maps_no = layer["m"];
    end
  end
end

function backprop{T}(config::NetworkConfig, X::Array{T, 2}, L::Array{T, 2},
                     W::Weights{T}, nn::FwdComputation{T},
                     err::BackComputation{T})
  # Se calculează erorile și gradienții.
  # Ultimul strat se calculează înaintea celorlalte.
  # Într-un layer maxpul facem așa:
  #   1. calculăm eroarea (ca sumă, plm, după toate cele de pe stratul viitor)
  #   2. calculăm gradient (plm, super simplu, o înmulțire)
  #   3. punem eroarea și la elementul maxim de la convoluta de mai devreme
  # La o convolută facem după cum se scurge mai departe:
  #   1. ne-a pus maxpul erorile în erori
  #   2. calculăm doar gradienții
  #   3. gradienții să fie medie peste toate erorile puse în aceeași mapă (hartă n.r.)
  # La o full e cel mai simplu
  #   1. se adună și se înmulțesc
  # încă ceva, presupune că ultimul strat este full conectat

  forward(config, X, W, nn);  # forward
  # Last layer
  println("dam la iesire intai");

  l = length(config);

  layer = config[l];

  if layer["activation"] == "softmax" || layer["activation"] == "sigmoid"
    err.E[l][:,:] = nn.Z[l] .- L';
  else
    println("Exceptie: nu s-a implementat inca");
  end

  for j in 1:layer["m"]
    if config[l-1]["type"] == "FULL"
      err.G[l][j, 2:end] = sum(err.E[l][j] * nn.Z[l-1],1);
      err.G[l][j, 1] = err.E[l][j];
    elseif config[l-1]["type"] == "MAXPOOL"
      i = 1;
      err.G[l][j, i] = err.E[l][j];
      for m in 1:config[l-1]["maps"]
        for l1 in 1:config[l-1]["width"]
          for l2 in 1:config[l-1]["width"]
            i = i+1;
            err.G[l][j, i] = sum(err.E[l][j] * nn.Z[l-1][m, l1, l2]);
          end
        end
      end
    else
      println("Daca ai ajuns aici esti un chinez!");
    end
  end

  # All the other layers
  for l in (length(config)-1):-1:2
    println("facem stratul ", l);
    layer = config[l];
    if layer["type"] == "CONV"
      println("Gigi e convolut");
      # presupunem ca maxpul ne-a pus eroarea
    elseif layer["type"] == "MAXPOOL"
      # maxpul -> full
      if config[l+1]["type"] == "FULL"
        for m in 1:layer["maps"]
          width = size(nn.Z[l],3);
          for i in 1:width
            for j in 1:width
              err.E[l+1] * W[l+1][:,(m-1) * width*width + (i-1)*width + j + 1];
            end
          end
        end
      else
        #it should be conv
        for m in 1:layer["maps"]
          width = size(nn.Z[l],3);
          for i in 1:width
            for j in 1:width
              sum(err.E[l+1][,:m,:,:] * [:,:,2:])
            end
          end
        end
      end

    elseif layer["type"] == "FULL"
      println("Gigi e complet conectat");
    else
      println("Daca ai ajuns aici esti un oier!");
    end
  end
end

function test()
  d = read_mnist_data();
  w = init_weights(Float64, layers, -0.1, 0.1);
  nn, err = prepare_network(Float64, layers, 100);
  describe_weights(w);
  backprop(layers, d["train_img"][:,1:100], d["train_label"][:,1:100], w, nn, err);
end

end

max_stars_repo_name: tudor-berariu/ml-hacks

max_issues_repo_name: tudor-berariu/ml-hacks

max_forks_repo_name: tudor-berariu/ml-hacks


text:
@testset "Timing Accuracy" begin
    n = 7
    t0 = 7257.0
    h = 0.01
    tmax = 10.0
    intr = Integrator(ahl21!,h,t0,tmax)

    # Read in initial conditions:
    H = [n,ones(Int64, n - 1)...]
    elements = readdlm("elements.txt", ',')[1:n,:]
    ic = ElementsIC(t0,H,elements)
    s = State(ic)

    # Holds transit times, derivatives, etc.
    tt = TransitTiming(tmax,ic);

    # Run integrator and calculate times
    intr(s,tt,grad=false)

    # Check that we get back the initial transit time (up to tolerance)
    tol = 1e-6 # Percent error
    for i in 1:n-1
        @test abs((ic.elements[i+1,3] - tt.tt[i+1,1])/ic.elements[i+1,3]) < tol
    end
end


text:
using FinancialSymbology
using Test, Documenter

import FinancialSymbology: figiidtype, makeurl, makejob, splitjobs

idstrings = [
    "FR0000121014", "US88160R1014", "AAPL US Equity", "5505072", "CH0210483332",
    "BYVZLD1", "BDDXSM4", "42751Q105", "FR0000120693", "654106103", "GB0002374006"
]

idchecks = [
    Isin, Isin, Ticker, Sedol, Isin, Sedol, Sedol, Cusip, Isin, Cusip, Isin
]

figitypes = [
    "ID_ISIN","ID_ISIN","TICKER","ID_SEDOL","ID_ISIN",
    "ID_SEDOL","ID_SEDOL","ID_CUSIP","ID_ISIN","ID_CUSIP","ID_ISIN"
]

sedoljob = Dict([
    :idType=>"ID_SEDOL",
    :idValue=>"5505072"
])

tickerjob = Dict([
    :idType=>"TICKER",
    :idValue=>"AAPL",
    :exchCode=>"US",
    :marketSecDes=>"Equity"
])

jobs = [sedoljob, tickerjob]

splitjob = [[sedoljob],[tickerjob]]

ids = makeidentifier.(idstrings)

@testset "IdChecks.jl" begin
    @test all(map((id, t) -> id isa t, ids, idchecks))
    doctest(FinancialSymbology; manual=false)
end

@testset "OpenFigiApi.jl" begin
    @test all(map((id, t) -> figiidtype(id) == t, ids, figitypes))
    @test makeurl(OpenFigiAPI()) == "https://api.openfigi.com/v3/mapping"
    job = makejob(ids[4])
    @test all(map(k -> job[k] == sedoljob[k], sedoljob |> keys |> collect))
    job = makejob(ids[3])
    @test all(map(k -> job[k] == tickerjob[k], tickerjob |> keys |> collect))
    resp = fetchsecuritydata(ids)
    @test length(resp) == length(ids)
    fetchsecuritydata(ids[1])
    fetchsecuritydata([ids[1]])
    fetchsecuritydata(ids[2])
end

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
#=
Test resynthesizer in 3D using commonly available test data.
=#

# The tested item
include("../../src/inpaint.jl")

using ImageView, Images

# Grab the open source test data
using TestImages
# "mri" is short name
fullImage = testimage("mri");

println("fullImage is:")
println(summary(fullImage))

#=
See JuliaImages.
The image has named axes.
We only want to convert the data.
=#
image = fullImage.data
println("image is:")
println(summary(image))

# Create a mask
mask = falses(size(image))
mask[3:4, 2:3, 5:6] .= true

println("mask is:")
println(summary(mask))

function test(data, mask)
    result = inpaint(data, mask)

    # Assert result is 3D array

    #=
    Chunk it back into the Image for display.
    =#
    #fullImage.data = result
    imshow(result)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# set voltage limit to max value
# SUM
export psp_setvlimitmax

const VLIMIT_MAX = [0x53, 0x55, 0x4D, 0x0D]

"""
    psp_setvlimitmax(io)

Set voltage limit to maximum value.
"""
function psp_setvlimitmax(io_psp::IO)
  write(io_psp, VLIMIT_MAX)
  return nothing
end


text:
module Revise

using FileWatching, REPL, Distributed, UUIDs
import LibGit2
using Base: PkgId

using OrderedCollections: OrderedDict

export revise, includet, MethodSummary

"""
    Revise.watching_files[]

Returns `true` if we watch files rather than their containing directory.
FreeBSD and NFS-mounted systems should watch files, otherwise we prefer to watch
directories.
"""
const watching_files = Ref(Sys.KERNEL == :FreeBSD)

"""
    Revise.polling_files[]

Returns `true` if we should poll the filesystem for changes to the files that define
loaded code. It is preferable to avoid polling, instead relying on operating system
notifications via `FileWatching.watch_file`. However, NFS-mounted
filesystems (and perhaps others) do not support file-watching, so for code stored
on such filesystems you should turn polling on.

See the documentation for the `JULIA_REVISE_POLL` environment variable.
"""
const polling_files = Ref(false)
function wait_changed(file)
    try
        polling_files[] ? poll_file(file) : watch_file(file)
    catch err
        if Sys.islinux() && err isa SystemError && err.errnum == 28  # ENOSPC
            @warn """Your operating system has run out of inotify capacity.
            Check the current value with `cat /proc/sys/fs/inotify/max_user_watches`.
            Set it to a higher level with, e.g., `echo 65536 | sudo tee -a /proc/sys/fs/inotify/max_user_watches`.
            This requires having administrative privileges on your machine (or talk to your sysadmin).
            See https://github.com/timholy/Revise.jl/issues/26 for more information."""
        end
        rethrow(err)
    end
    return nothing
end

"""
    Revise.tracking_Main_includes[]

Returns `true` if files directly included from the REPL should be tracked.
The default is `false`. See the documentation regarding the `JULIA_REVISE_INCLUDE`
environment variable to customize it.
"""
const tracking_Main_includes = Ref(false)

include("relocatable_exprs.jl")
include("types.jl")
include("parsing.jl")
include("exprutils.jl")
include("pkgs.jl")
include("git.jl")
include("recipes.jl")
include("logging.jl")

### Globals to keep track of state

"""
    Revise.watched_files

Global variable, `watched_files[dirname]` returns the collection of files in `dirname`
that we're monitoring for changes. The returned value has type [`WatchList`](@ref).

This variable allows us to watch directories rather than files, reducing the burden on
the OS.
"""
const watched_files = Dict{String,WatchList}()

"""
    Revise.revision_queue

Global variable, `revision_queue` holds the names of files that we need to revise, meaning
that these files have changed since we last processed a revision.
This list gets populated by callbacks that watch directories for updates.
"""
const revision_queue = Set{String}()

"""
    Revise.fileinfos

`fileinfos` is the core information that tracks the relationship between source code
and julia objects, and allows re-evaluation of code in the proper module scope.
It is a dictionary indexed by absolute paths of files;
`fileinfos[filename]` returns a value of type [`Revise.FileInfo`](@ref).
"""
const fileinfos = Dict{String,FileInfo}()

"""
    Revise.module2files

`module2files` holds the list of filenames used to define a particular
module. This is only used by `revise(MyModule)` to "refresh" all the definitions in a module.
"""
const module2files = Dict{Symbol,Vector{String}}()

"""
    Revise.included_files

Global variable, `included_files` gets populated by callbacks we register with `include`.
It's used to track non-precompiled packages and, optionally, user scripts (see docs on
`JULIA_REVISE_INCLUDE`).
"""
const included_files = Tuple{Module,String}[]  # (module, filename)

"""
    Revise.basesrccache

Full path to the running Julia's cache of source code defining `Base`.
"""
const basesrccache = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia", "base.cache")

"""
    Revise.basebuilddir

Julia's top-level directory when Julia was built, as recorded by the entries in
`Base._included_files`.
"""
const basebuilddir = begin
    sysimg = filter(x->endswith(x[2], "sysimg.jl"), Base._included_files)[1][2]
    dirname(dirname(sysimg))
end

"""
    Revise.juliadir

Constant specifying full path to julia top-level source directory.
This should be reliable even for local builds, cross-builds, and binary installs.
"""
const juliadir = begin
    local jldir = basebuilddir
    if !isdir(joinpath(jldir, "base"))
        # Binaries probably end up here. We fall back on Sys.BINDIR
        jldir = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia")
        if !isdir(joinpath(jldir, "base"))
            while true
                trydir = joinpath(jldir, "base")
                isdir(trydir) && break
                trydir = joinpath(jldir, "share", "julia", "base")
                if isdir(trydir)
                    jldir = joinpath(jldir, "share", "julia")
                    break
                end
                jldirnext = dirname(jldir)
                jldirnext == jldir && break
                jldir = jldirnext
            end
        end
    end
    jldir
end
const cache_file_key = Dict{String,String}() # corrected=>uncorrected filenames

"""
    Revise.dont_watch_pkgs

Global variable, use `push!(Revise.dont_watch_pkgs, :MyPackage)` to prevent Revise
from tracking changes to `MyPackage`. You can do this from the REPL or from your
`.julia/config/startup.jl` file.

See also [`Revise.silence`](@ref).
"""
const dont_watch_pkgs = Set{Symbol}()
const silence_pkgs = Set{Symbol}()
const depsdir = joinpath(dirname(@__DIR__), "deps")
const silencefile = Ref(joinpath(depsdir, "silence.txt"))  # Ref so that tests don't clobber

function use_compiled_modules()
    return Base.JLOptions().use_compiled_modules != 0
end

"""
    fmrep = eval_revised(fmnew::FileModules, fmref::FileModules)

Implement the changes from `fmref` to `fmnew`, returning a replacement [`FileModules`](@ref)
`fmrep`.
"""
function eval_revised(fmnew::FileModules, fmref::FileModules)
    fmrep = FileModules(first(keys(fmref)))  # replacement for fmref
    for (mod, fmmnew) in fmnew
        fmrep[mod] = fmm = FMMaps()
        if haskey(fmref, mod)
            eval_revised!(fmm, mod, fmmnew, fmref[mod])
            for p in workers()
                p == myid() && continue
                remotecall(eval_revised_dummy!, p, mod, fmmnew, fmref[mod])
            end
        else  # a new submodule (see #43)
            eval_and_insert_all!(fmm, mod, fmmnew.defmap)
            for p in workers()
                p == myid() && continue
                remotecall(eval_and_insert_all_dummy!, p, mod, fmmnew.defmap)
            end
        end
    end
    fmrep
end

function eval_revised!(fmmrep::FMMaps, mod::Module,
                       fmmnew::FMMaps, fmmref::FMMaps)
    # Update to the state of fmmnew, preventing any unnecessary evaluation
    with_logger(_debug_logger) do
        @debug "Diff" _group="Parsing" activemodule=fullname(mod) newexprs=setdiff(keys(fmmnew.defmap), keys(fmmref.defmap)) oldexprs=setdiff(keys(fmmref.defmap), keys(fmmnew.defmap))
        # Delete any methods missing in fmmnew
        # Because we haven't yet evaled and gotten signature types, let's do this "from scratch" using method signatures.
        # (It might be more efficient to eval first and delete later, but this proves to be safer.)
        oldsigs, newsigs = filter_signatures(mod, keys(fmmref.defmap)), filter_signatures(mod, keys(fmmnew.defmap))
        for sig in oldsigs
            sig ∈ newsigs && continue
            @debug "DiffSig" _group="Parsing" activemodule=fullname(mod) missingsig=sig newsigs=newsigs
            for sigt in sigex2sigts(mod, sig)
                m = get_method(sigt)
                if isa(m, Method)
                    Base.delete_method(m)
                    @debug "DeleteMethod" _group="Action" time=time() deltainfo=(sigt, MethodSummary(m))
                else
                    mths = Base._methods_by_ftype(sigt, -1, typemax(UInt))
                    io = IOBuffer()
                    println(io, "Extracted method table:")
                    println(io, mths)
                    info = String(take!(io))
                    @warn "Revise failed to find any methods for signature $sigt\n  Perhaps it was already deleted.\n$info"
                end
            end
        end
        # Eval any expressions unique to fmmnew
        for (def,val) in fmmnew.defmap
            @assert def != nothing
            defref = getkey(fmmref.defmap, def, nothing)
            if defref != nothing
                # The same expression is found in both, only update the lineoffset
                if val !== nothing
                    sigtref, oldoffset = fmmref.defmap[defref]
                    lnref = firstlineno(defref)
                    lnnew = firstlineno(def)
                    lineoffset = (isa(lnref, Integer) && isa(lnnew, Integer)) ? lnnew-lnref : 0
                    fmmrep.defmap[defref] = (sigtref, lineoffset)
                    for sigt in sigtref
                        fmmrep.sigtmap[sigt] = defref
                    end
                    oldoffset != lineoffset && @debug "LineOffset" _group="Action" time=time() deltainfo=(sigtref, lnnew, oldoffset=>lineoffset)
                else
                    fmmrep.defmap[defref] = nothing
                end
            else
                eval_and_insert!(fmmrep, mod, def=>val)
            end
        end
    end
    return fmmrep
end

function filter_signatures(mod::Module, defs)
    sigs = Set{RelocatableExpr}()
    for def in defs
        while is_trivial_block_wrapper(def)
            def = def.args[end]
        end
        def isa ExLike || continue
        if def.head == :macrocall
            _, def = macexpand(mod, convert(Expr, def))
            def isa ExLike || continue
        end
        def.head == :function || def.head == :(=) || continue
        sig = get_signature(def)
        sig isa ExLike || continue
        push!(sigs, relocatable!(sig))
    end
    return sigs
end

eval_revised_dummy!(mod::Module, fmmnew::FMMaps, fmmref::FMMaps) =
    eval_revised!(FMMaps(), mod, fmmnew, fmmref)

function eval_and_insert!(fmm::FMMaps, mod::Module, pr::Pair)
    def, val = pr.first, pr.second
    ex = convert(Expr, def)
    try
        if isdocexpr(ex) && mod == Base.__toplevel__
            Core.eval(Main, ex)
        else
            Core.eval(mod, ex)
        end
        @debug "Eval" _group="Action" time=time() deltainfo=(mod, relocatable!(ex))
        if val isa RelocatableExpr
            instantiate_sigs!(fmm, def, val, mod)
        else
            fmm.defmap[def] = val
            if val !== nothing
                for sigt in val[1]
                    fmm.sigtmap[sigt] = def
                end
            end
        end
    catch err
        @error "failure to evaluate changes in $mod"
        showerror(stderr, err)
        println_maxsize(stderr, "\n", ex; maxlines=20)
    end
    return fmm
end

function eval_and_insert_all!(fmm::FMMaps, mod::Module, defmap::DefMap)
    for pr in defmap
        eval_and_insert!(fmm, mod, pr)
    end
end

eval_and_insert_all_dummy!(mod::Module, defmap::DefMap) =
    eval_and_insert_all!(FMMaps(), mod, defmap)

"""
    Revise.init_watching(files)

For every filename in `files`, monitor the filesystem for updates. When the file is
updated, either [`revise_dir_queued`](@ref) or [`revise_file_queued`](@ref) will
be called.
"""
function init_watching(files)
    udirs = Set{String}()
    for file in files
        dir, basename = splitdir(file)
        haskey(watched_files, dir) || (watched_files[dir] = WatchList())
        push!(watched_files[dir], basename)
        if watching_files[]
            @async revise_file_queued(file)
        else
            push!(udirs, dir)
        end
    end
    for dir in udirs
        updatetime!(watched_files[dir])
        @async revise_dir_queued(dir)
    end
    return nothing
end

"""
    revise_dir_queued(dirname)

Wait for one or more of the files registered in `Revise.watched_files[dirname]` to be
modified, and then queue the corresponding files on [`Revise.revision_queue`](@ref).
This is generally called within an `@async`.
"""
@noinline function revise_dir_queued(dirname)
    if !isdir(dirname)
        sleep(0.1)   # in case git has done a delete/replace cycle
        if !isfile(dirname)
            with_logger(SimpleLogger(stderr)) do
                @warn "$dirname is not an existing directory, Revise is not watching"
            end
            return nothing
        end
    end
    latestfiles = watch_files_via_dir(dirname)  # will block here until file(s) change
    for file in latestfiles
        push!(revision_queue, file)
    end
    @async revise_dir_queued(dirname)
end

# See #66.
"""
    revise_file_queued(filename)

Wait for modifications to `filename`, and then queue the corresponding files on [`Revise.revision_queue`](@ref).
This is generally called within an `@async`.

This is used only on platforms (like BSD) which cannot use [`revise_dir_queued`](@ref).
"""
function revise_file_queued(file)
    if !isfile(file)
        sleep(0.1)  # in case git has done a delete/replace cycle
        if !isfile(file)
            with_logger(SimpleLogger(stderr)) do
                @error "$file is not an existing file, Revise is not watching"
            end
            return nothing
        end
    end

    wait_changed(file)  # will block here until the file changes
    push!(revision_queue, file)
    @async revise_file_queued(file)
end

"""
    Revise.revise_file_now(file)

Process revisions to `file`. This parses `file` and computes an expression-level diff
between the current state of the file and its most recently evaluated state.
It then deletes any removed methods and re-evaluates any changed expressions.

`file` must be a key in [`Revise.fileinfos`](@ref)
"""
function revise_file_now(file)
    if !haskey(fileinfos, file)
        println("Revise is currently tracking the following files: ", keys(fileinfos))
        error(file, " is not currently being tracked.")
    end
    fi = fileinfos[file]
    maybe_parse_from_cache!(fi, file)
    fmref = fi.fm
    topmod = first(keys(fi.fm))
    fmnew = parse_source(file, topmod)
    if fmnew != nothing
        fmrep = eval_revised(fmnew, fmref)
        fileinfos[file] = FileInfo(fmrep, fi)
    end
    nothing
end

function instantiate_sigs!(fm::FileModules)
    for (mod, fmm) in fm
        instantiate_sigs!(fmm, mod)
    end
    return fm
end

function instantiate_sigs!(fmm::FMMaps, mod::Module)
    for (def, sig) in fmm.defmap
        if sig isa RelocatableExpr
            instantiate_sigs!(fmm, def, sig, mod)
        end
    end
    return fmm
end

function instantiate_sigs!(fmm::FMMaps, def::RelocatableExpr, sig::RelocatableExpr, mod::Module)
    sigts = sigex2sigts(mod, sig, def)
    # Insert into the maps
    fmm.defmap[def] = (sigts, 0)
    for sigt in sigts
        fmm.sigtmap[sigt] = def
    end
    return def
end

"""
    revise()

`eval` any changes in the revision queue. See [`Revise.revision_queue`](@ref).
"""
function revise()
    sleep(0.01)  # in case the file system isn't quite done writing out the new files
    for file in revision_queue
        revise_file_now(file)
    end
    empty!(revision_queue)
    tracking_Main_includes[] && queue_includes(Main)
    nothing
end

# This variant avoids unhandled task failures
function revise(backend::REPL.REPLBackend)
    sleep(0.01)  # in case the file system isn't quite done writing out the new files
    for file in revision_queue
        try
            revise_file_now(file)
        catch err
            put!(backend.response_channel, (err, catch_backtrace()))
        end
    end
    empty!(revision_queue)
    tracking_Main_includes[] && queue_includes(Main)
    nothing
end

"""
    revise(mod::Module)

Reevaluate every definition in `mod`, whether it was changed or not. This is useful
to propagate an updated macro definition, or to force recompiling generated functions.
"""
function revise(mod::Module)
    for file in module2files[Symbol(mod)]
        for (mod,fmm) in fileinfos[file].fm
            for def in keys(fmm.defmap)
                Core.eval(mod, convert(Expr, def))
            end
        end
    end
    return true  # fixme try/catch?
end

"""
    Revise.track(mod::Module, file::AbstractString)
    Revise.track(file::AbstractString)

Watch `file` for updates and [`revise`](@ref) loaded code with any
changes. `mod` is the module into which `file` is evaluated; if omitted,
it defaults to `Main`.
"""
function track(mod::Module, file::AbstractString)
    isfile(file) || error(file, " is not a file")
    file = normpath(abspath(file))
    fm = parse_source(file, mod)
    if fm != nothing
        instantiate_sigs!(fm)
        fileinfos[file] = FileInfo(fm)
    end
    init_watching((file,))
end

track(file::AbstractString) = track(Main, file)

"""
    includet(filename)

Load `filename` and track any future changes to it.
`includet` is deliberately non-recursive, so if `filename` loads any other files,
they will not be automatically tracked.
(See [`Revise.track`](@ref) to set it up manually.)

`includet` is intended for "user scripts," e.g., a file you use locally for a specific
purpose such as loading a specific data set or performing some kind of analysis.
Do not use `includet` for packages, as those should be handled by `using` or `import`.
If `using` and `import` aren't working, you may have packages in a non-standard location;
try fixing it with something like `push!(LOAD_PATH, "/path/to/my/private/repos")`.
"""
includet(file::AbstractString) = (Base.include(Main, file); track(Main, file))

"""
    Revise.silence(pkg)

Silence warnings about not tracking changes to package `pkg`.
"""
function silence(pkg::Symbol)
    push!(silence_pkgs, pkg)
    if !isdir(depsdir)
        mkpath(depsdir)
    end
    open(silencefile[], "w") do io
        for p in silence_pkgs
            println(io, p)
        end
    end
    nothing
end
silence(pkg::AbstractString) = silence(Symbol(pkg))

## Utilities

"""
    method = get_method(sigt)

Get the method `method` with signature-type `sigt`. This is used to provide
the method to `Base.delete_method`. See also [`get_signature`](@ref).

If `sigt` does not correspond to a method, returns `nothing`.

# Examples

```jldoctest; setup = :(using Revise), filter = r"in Main at.*"
julia> mymethod(::Int) = 1
mymethod (generic function with 1 method)

julia> mymethod(::AbstractFloat) = 2
mymethod (generic function with 2 methods)

julia> Revise.get_method(Tuple{typeof(mymethod), Int})
mymethod(::Int64) in Main at REPL[0]:1

julia> Revise.get_method(Tuple{typeof(mymethod), Float64})
mymethod(::AbstractFloat) in Main at REPL[1]:1

julia> Revise.get_method(Tuple{typeof(mymethod), Number})

```
"""
function get_method(@nospecialize(sigt))
    mths = Base._methods_by_ftype(sigt, -1, typemax(UInt))
    length(mths) == 1 && return mths[1][3]
    if !isempty(mths)
        # There might be many methods, but the one that should match should be the
        # last one, since methods are ordered by specificity
        i = lastindex(mths)
        while i > 0
            m = mths[i][3]
            m.sig == sigt && return m
            i -= 1
        end
    end
    return nothing
end

"""
    rex = get_def(method::Method)

Return the RelocatableExpr defining `method`.
The source-file defining `method` must be tracked.
If it is in Base, this will execute `track(Base)` if necessary.
"""
function get_def(method::Method; modified_files=revision_queue)
    filename = String(method.file)
    yield()   # magic bug fix for the OSX test failures. TODO: figure out why this works (prob. Julia bug)
    startswith(filename, "REPL") && error("methods defined at the REPL are not yet supported")
    # First we check for whether `filename` is consistent with the module.
    # If not, most likely this is a method that was defined via a macro defined in a different module.
    # E.g., https://github.com/timholy/Rebugger.jl/issues/3
    modfiles = get(module2files, nameof(method.module), nothing)
    if modfiles != nothing && filename ∉ modfiles
        # Try to find the actual file
        for file in modfiles
            fi = fileinfos[file]
            maybe_parse_from_cache!(fi, file)
            haskey(fi.fm, method.module) || continue
            map = fi.fm[method.module].sigtmap
            def = get(map, method.sig, nothing)
            def == nothing || return copy(def)
        end
    end
    recipemod = method.module
    if !haskey(fileinfos, filename)
        ret = find_file(filename, recipemod)
        if ret == nothing
            @warn "file $filename is not tracked by Revise"
            return nothing
        end
        filename, recipemod = ret
    end
    filename isa AbstractString || return nothing
    if !haskey(fileinfos, filename)
        @info "tracking $recipemod"
        track(recipemod; modified_files=modified_files)
    end
    fi = fileinfos[filename]
    maybe_parse_from_cache!(fi, filename)
    map = fi.fm[method.module].sigtmap
    haskey(map, method.sig) && return copy(map[method.sig])
    @warn "$(method.sig) was not found, perhaps it was generated by code"
    nothing
end

function find_file(filename, mod)
    # For files not in fileinfos, see if we can identify their origin in Base or the stdlibs
    handledmods = (Base, Core.Compiler)
    newmod = parentmodule(mod)
    while mod ∉ handledmods && newmod != mod
        mod = newmod
        newmod = parentmodule(mod)
    end
    if mod ∈ handledmods
        basefile = Base.find_source_file(filename)  # check base
        basefile isa AbstractString && return fixpath(realpath(basefile)), mod
    elseif occursin("stdlib", filename)
        return fixpath(realpath(filename)), mod
    end
    return nothing
end

function printf_maxsize(f::Function, io::IO, args...; maxchars::Integer=500, maxlines::Integer=20)
    # This is dumb but certain to work
    iotmp = IOBuffer()
    for a in args
        print(iotmp, a)
    end
    print(iotmp, '\n')
    seek(iotmp, 0)
    str = read(iotmp, String)
    if length(str) > maxchars
        str = first(str, (maxchars+1)÷2) * "…" * last(str, maxchars - (maxchars+1)÷2)
    end
    lines = split(str, '\n')
    if length(lines) <= maxlines
        for line in lines
            f(io, line)
        end
        return
    end
    half = (maxlines+1) ÷ 2
    for i = 1:half
        f(io, lines[i])
    end
    maxlines > 1 && f(io, ⋮)
    for i = length(lines) - (maxlines-half) + 1:length(lines)
        f(io, lines[i])
    end
end
println_maxsize(args...; kwargs...) = println_maxsize(stdout, args...; kwargs...)
println_maxsize(io::IO, args...; kwargs...) = printf_maxsize(println, stdout, args...; kwargs...)

function fix_line_statements!(ex::Expr, file::Symbol, line_offset::Int=0)
    if ex.head == :line
        ex.args[1] += line_offset
        ex.args[2] = file
    else
        for (i, a) in enumerate(ex.args)
            if isa(a, Expr)
                fix_line_statements!(a::Expr, file, line_offset)
            elseif isa(a, LineNumberNode)
                ex.args[i] = file_line_statement(a::LineNumberNode, file, line_offset)
            end
        end
    end
    ex
end

file_line_statement(lnn::LineNumberNode, file::Symbol, line_offset) =
    LineNumberNode(lnn.line + line_offset, file)

function update_stacktrace_lineno!(trace)
    for i = 1:length(trace)
        t, n = trace[i]
        if t.linfo isa Core.MethodInstance
            sigt = t.linfo.def.sig
            file = String(t.file)
            if haskey(fileinfos, file)
                fm = fileinfos[file].fm
                for (mod, fmm) in fm
                    if haskey(fmm.sigtmap, sigt)
                        def = fmm.sigtmap[sigt]
                        lineoffset = fmm.defmap[def][2]
                        if lineoffset != 0
                            t = StackTraces.StackFrame(t.func, t.file, t.line+lineoffset, t.linfo, t.from_c, t.inlined, t.pointer)
                            trace[i] = (t, n)
                            break
                        end
                    end
                end
            end
        end
    end
    return trace
end

function macroreplace!(ex::Expr, filename)
    for i = 1:length(ex.args)
        ex.args[i] = macroreplace!(ex.args[i], filename)
    end
    if ex.head == :macrocall
        m = ex.args[1]
        if m == Symbol("@__FILE__")
            return String(filename)
        elseif m == Symbol("@__DIR__")
            return dirname(String(filename))
        end
    end
    return ex
end
macroreplace!(s, filename) = s

@noinline function run_backend(backend)
    while true
        tls = task_local_storage()
        tls[:SOURCE_PATH] = nothing
        ast, show_value = take!(backend.repl_channel)
        if show_value == -1
            # exit flag
            break
        end
        # Process revisions
        revise(backend)
        # Now eval the input
        REPL.eval_user_input(ast, backend)
    end
    nothing
end

"""
    steal_repl_backend(backend = Base.active_repl_backend)

Replace the REPL's normal backend with one that calls [`revise`](@ref) before executing
any REPL input.
"""
function steal_repl_backend(backend = Base.active_repl_backend)
    @async begin
        # terminate the current backend
        put!(backend.repl_channel, (nothing, -1))
        fetch(backend.backend_task)
        # restart a new backend that differs only by processing the
        # revision queue before evaluating each user input
        backend.backend_task = @async run_backend(backend)
    end
    nothing
end

function wait_steal_repl_backend()
    iter = 0
    # wait for active_repl_backend to exist
    while !isdefined(Base, :active_repl_backend) && iter < 20
        sleep(0.05)
        iter += 1
    end
    if isdefined(Base, :active_repl_backend)
        steal_repl_backend(Base.active_repl_backend)
    else
        @warn "REPL initialization failed, Revise is not in automatic mode. Call `revise()` manually."
    end
end

"""
    Revise.async_steal_repl_backend()

Wait for the REPL to complete its initialization, and then call [`steal_repl_backend`](@ref).
This is necessary because code registered with `atreplinit` runs before the REPL is
initialized, and there is no corresponding way to register code to run after it is complete.
"""
function async_steal_repl_backend()
    mode = get(ENV, "JULIA_REVISE", "auto")
    if mode == "auto"
        atreplinit() do repl
            @async wait_steal_repl_backend()
        end
    end
    return nothing
end

function __init__()
    # Base.register_root_module(Base.__toplevel__)
    if isfile(silencefile[])
        pkgs = readlines(silencefile[])
        for pkg in pkgs
            push!(silence_pkgs, Symbol(pkg))
        end
    end
    push!(Base.package_callbacks, watch_package)
    push!(Base.include_callbacks,
        (mod::Module, fn::AbstractString) -> push!(included_files, (mod, normpath(abspath(fn)))))
    mode = get(ENV, "JULIA_REVISE", "auto")
    if mode == "auto"
        if isdefined(Base, :active_repl_backend)
            steal_repl_backend(Base.active_repl_backend)
        elseif isdefined(Main, :IJulia)
            Main.IJulia.push_preexecute_hook(revise)
        end
        if isdefined(Main, :Atom)
            for x in ["eval", "evalall", "evalrepl"]
                old = Main.Atom.handlers[x]
                Main.Atom.handle(x) do data
                    revise()
                    old(data)
                end
            end
        end
    end
    polling = get(ENV, "JULIA_REVISE_POLL", "0")
    if polling == "1"
        polling_files[] = watching_files[] = true
    end
    rev_include = get(ENV, "JULIA_REVISE_INCLUDE", "0")
    if rev_include == "1"
        tracking_Main_includes[] = true
    end
    # Correct line numbers for code moving around
    Base.update_stackframes_callback[] = update_stacktrace_lineno!
end

## WatchList utilities
function systime()
    tv = Libc.TimeVal()
    tv.sec + tv.usec/10^6
end
function updatetime!(wl::WatchList)
    wl.timestamp = systime()
end
Base.push!(wl::WatchList, filename) = push!(wl.trackedfiles, filename)
WatchList() = WatchList(systime(), Set{String}())
Base.in(file, wl::WatchList) = in(file, wl.trackedfiles)

include("precompile.jl")
_precompile_()

end # module


text:
using DynamicIterators
using Distributions

D = MvNormal([1.0, 0.5], [1.0 0.5; 0.5 1.5] )
struct Move{T}
    x::T
    σ::Float64
    i::Int
end
m1(x) = Move(x, 0.1, 1)
m2(x) = Move(x, 0.1, 2)
Base.rand(M::Move) = M.x + M.σ*randn()*[M.i-1, 2-M.i]
Distributions.logpdf(M::Move, x) = logpdf(Normal(M.x[M.i], M.σ), x[M.i])
MH1 = MetropolisHastings(D, m1, logpdf)
MH2 = MetropolisHastings(D, m2, logpdf)

I = Evolve(i->rand(1:2))

MH = mixture(I, (MH1, MH2))

X = values(trace(MH, 1=>(1, [0.0, 0.0]), endtime(2000)))


text:
export circuit_eqs, sheet_resistance, currents

"""
```julia
circuit_eqs(nwn::NWN_circuit{M,T,N})
```
Creates the system of equations which solves the
electrical problem of the nanowire network under
the given model `M`. Removes ill conditioned nodes (those not
connected through the circuit to a voltage source
or ground).
"""
function circuit_eqs(nwn::NWN_circuit{M,T,N}) where {M,T,N}
    nₗ, nₑ, n₀ = length(nwn.lines), length(nwn.elecs), length(nwn.grnds)
    nₗₙ = nwn.ns[nₗ] # number of line related nodes
    nₑₙ = nwn.ns[nₗ + nₑ] - nₗₙ
    # n₀ₙ = nwn.ns[nₗ + nₑ + n₀] - nₗₙ - nₑₙ

    # include only nodes part of the actual circuit
    con_comps = connected_components(nwn.graph)
    connected = vcat(con_comps[.!(con_comps .⊆ [1:nₗₙ])]...)

    # remove grounds from equations as MNA requires
    inds = setdiff(connected, (nwn.ns[nₗ+nₑ]+1):nwn.ns[nₗ+nₑ+n₀])

    # view of Laplacian of graph with indices needed
    L = -laplacian_matrix(nwn.graph)[inds,inds]

    # generates B matrix
    B = spzeros(T,length(inds),nₑₙ)
    vs = spzeros(T,nₑₙ)
    for i ∈ 1:nₑ
        js = (nwn.ns[nₗ+i]-nwn.ls[nₗ+i]+1):nwn.ns[nₗ+i]
        B[CartesianIndex.(findall(inds.∈[js]),js.-nₗₙ)] .= one(T)
        vs[js.-nₗₙ] .= nwn.volts[i]
    end

    # generates system of equations in the form Ax=z
    A = [L B;
        transpose(B) spzeros(T,nₑₙ,nₑₙ)]
    z = [spzeros(T,length(inds));
        vs]
    A,z,inds
end

"""
```julia
sheet_resistance(nwn::NWN_circuit)
```
Calculates the sheet resistance of a nanowire network object.
Solves the network using MNA and returns the ratio of the voltage
to the current. Requires the network to have only one electrode.
"""
function sheet_resistance(nwn::NWN_circuit)
    @assert length(nwn.elecs)==1 "Sheet resistance is only defined for one electrode and one ground."
    A, z, inds = circuit_eqs(nwn)
    sol = A\collect(z)
    nwn.volts[1]/sum(sol[length(inds)+1:end])
end

function sheet_resistance(nwn::NWN_circuit{M,BigFloat,N}) where {M,N}
    @assert length(nwn.elecs)==1 "Sheet resistance is only defined for one electrode and one ground."
    A, z, inds = circuit_eqs(nwn)
    # For BigFloat the matrix must be dense
    sol = collect(A)\collect(z)
    nwn.volts[1]/sol[end]
end

"""
```julia
currents(nwn::NWN_circuit{M,T,N})
```
Returns a list of currents as iterated by `edges(nwn.graph)` and using thhe node labeling as defined by the model `M`.
"""
function currents(nwn::NWN_circuit{M,T,N}) where {M,T,N}
	A,z,inds = circuit_eqs(nwn)
	sol = A\collect(z)
	Vs = zeros(T,nv(nwn.graph))
	Vs[inds] = sol[1:length(inds)]
	es = edges(nwn.graph)
	# (Vₒᵤₜ-Vᵢₙ)/R
	(Vs[dst.(es)].-Vs[src.(es)]).*weight.(es)
end


text:
# Modelo de turbina com sangria
type Turbina_sangra
	Turbina_sangra()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Entalpia(),
			Eficiencia(Dict{Symbol,Any}(
				:Brief=>"Eficiencia da turbina"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia da turbina",
				:PosX=>0.9,
				:PosY=>0.45
			)),
			Fracao(Dict{Symbol,Any}(
				:Brief=>"Fracao massica da sangria"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>0,
				:PosY=>0.25
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>1,
				:PosY=>0.85
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{outx}",
				:PosX=>0.85,
				:PosY=>1
			)),
			[
				:(H_IS = PP2.propPS(Fout.P,Fin.S)),
				:(Fout.H = (H_IS - Fin.H) * EF_T + Fin.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
				:(Fin.F * (Fin.H - Fout.H) = POT_TURB),
				:(Fouts.F = Fin.F * y),
				:(Fout.F = Fin.F - Fouts.F),
				:(Fouts.P = Fout.P),
				:(Fouts.T = Fout.T),
				:(Fouts.S = Fout.S),
				:(Fouts.H = Fout.H),
			],
			[
				"","","","","","","","","","",
			],
			[:PP2,],
			[:H_IS,:EF_T,:POT_TURB,:y,:Fin,:Fout,:Fouts,]
		)
	end
	PP2::DanaPlugin
	H_IS::Entalpia
	EF_T::Eficiencia
	POT_TURB::Potencia
	y::Fracao
	Fin::Corrente
	Fout::Corrente
	Fouts::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Turbina_sangra
function setEquationFlow(in::Turbina_sangra)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
end
function atributes(in::Turbina_sangra,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/turbina_sa"
	drive!(fields,_)
	return fields
end
Turbina_sangra(_::Dict{Symbol,Any})=begin
	newModel=Turbina_sangra()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#=
    IndexLike

There are two components to any index, the keys and the values. These differ
from dictionaries in that the values must always be a subtype of
`AbstractUnitRange`. This facilitates a number of convenient traits and
assumptions concerning and index's behavior. When this behavior is overly
restrictive it makes sense use an index as component of some other object.
For example, a dictionary could potentially utilize the keys of an index
to map to the location of the dictionaries values.

# Keys

The only restriction to a collection of keys is that they all be unique to the
other keys in the same collection. This behavior is similar to what's seen if
you use `Set` or the keys of most `AbstractDict` subtypes. Keys of an index
allow the same customization and optimization that is commonly seen throughout
Julia via variety of dynamic, static, or immutable subtypes of `AbstractVector`.
Therefore, the only unique quality of keys for an index is the dichotomy offered
by the `ContinuousTrait` or `DiscreteTrait`.

* Continuous keys : Continuous keys imply an order the inherit to it's type.
  Therefore, all keys that are continuous are assumed to be some sort of range.
  This has the potential to speed up many operations but limits the ability to
  use keys with semantic meaning.
* Discrete keys: All vectors that are not ranges are assumed to be discrete,
  offering increased flexibility at the cost of performance.

It's worth noting that given some customization may help convert a system that
uses discrete keys into one of continuous keys. For example, using the
`Unitful.jl` package can facilitate the use of semantic units of measurement
as a well defined range.

# Values
=#