text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using JSON, Statistics, GitHub, Base64, SHA, Downloads, Dates, CairoMakie

include("benchmark-library.jl")

ctx = github_context()

projects = [
    "MakieDev" => "prio-obs-precompiles",
    "Latency" => "all",
    "MakieMesh" => "mesh-precompiles",
    "MakieCompiled" => "moah-precompiles]",
    "MakieMaster" => "master",
    "MakieAttributes" => "attribute-precompiles",
    "MakieBlock" => "blocks-precompile",
]

infos = map(projects) do (name, branch)
    return BenchInfo(
        project=joinpath(raw"C:\Users\sdani\SimiWorld\ProgrammerLife", name),
        branch=branch,
        commit=branch * "1"
    )
end

plot_url = run_benchmarks(ctx, [GitHub.branch(ctx.repo, "master"), GitHub.tag(ctx.repo, "v0.16.3"), infos...])

i1 = BenchInfo(commit="444eae5ce174d23d53c181144b357382bd57afa8", branch="julia1.9+this-pr")
i2 = BenchInfo(commit="7ccf35e789fccdda2429b39b6b10f4e91adcc5fd", branch="julia1.9+no-static-arrays")
i3 = BenchInfo(julia="julia-1-7-1", commit="f04ea93cde4384bad3f008bddf291a15ac35774a", branch="julia1.7+more-precompiles")
i4 = BenchInfo(julia="julia-1-7-1", commit="4a127bfd8dee589bbc3c2da1c91573a73752f36c", branch="julia1.7+v0.16.5")
i5 = BenchInfo(julia="julia-1-9-0-DEV-5", commit="7841b37e84921c4be8ec028b6aeb6cd12c036650", branch="master-yesterday")

fig = run_benchmarks(ctx, [i4, i3, i5, GitHub.branch(ctx.repo, "master"), i2, i1])

resize!(fig.scene, (1500, 600))
fig


text:
using FractalTools
using PyPlot
using Triangulation

close("all")

datapath = joinpath(@__DIR__,  "../../data/gmsh_data")
# datapath = joinpath(pathof(FractalTools), "data/gmsh_data")
filename = "t3.msh"

# Read mesh file.
mesh = read_mesh(joinpath(datapath, filename))
meshes = triangular_partition(mesh)
regions = gmsh2matplotlibtri.(meshes)
# region = regions[3]

# Compute z-values on mesh points.
f(x, y) = x.^2 .+ y.^2
z = map(region ->  f(region.x,region.y), regions)

# Calculate integral
I = map(i ->  integrate(regions[i], z[i], α=0.001), 1 : length(z))
@show I
#------------------------ Plots --------------------#

figure("interpolation")
plot_trisurf(region.x ,region.y, z)
triplot(region)
gcf()

ext: jl


text:
using Documenter, SMM

makedocs(
 pages = Any[
"Home" => "index.md"])


text:
using AbstractPlotting
using AbstractPlotting: textslider
using Colors
using Observables

function main()
    tslider1, _a = textslider(range(-10, stop = 10, length = 101), "a", start=0.)
    tslider2, _b = textslider(range(-10, stop = 10, length = 101), "b", start=0.)
    xs = -π:0.01:π
    ab = lift(tuple, _a,_b)
    pos = lift(ab) do (a, b)
        [Point2f0(x,sin(x)) for x in xs]
    end
    p = lines(pos)

    onany(pos) do pos
             AbstractPlotting.update!(p)
    end
    scene = hbox(p,vbox(tslider1,tslider2))
    RecordEvents(scene, "output")
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Plots

#points = rand(-5:5, 5).+im*rand(-5:5, 5)
#points = [0+3im, 1+3im, 2+0im, 1-3im, 0-3im, 1+0im, 0+3im]
#points = [1+im, 1-im, -1-im, -1+im, 1+im]
points = [784-413im,779-353im,576-354im,666-490im,573-604im,778-604im,790-534im,783-535im,768-562im,754-573im,738-577im,613-577im,691-479im,622-370im,738-368im,754-375im,771-396im,776-412im,784-413im]
#points="";document.onclick=()=>{points+=`${event.clientX}-${event.clientY}im,`;console.log(points)}
evolacc = 26
convacc = 50

function idoit()
    limx = [1.5*minimum(real.(points))-maximum(real.(points))/2, 1.5*maximum(real.(points))-minimum(real.(points))/2]
    limy = [1.5*minimum(imag.(points))-maximum(imag.(points))/2, 1.5*maximum(imag.(points))-minimum(imag.(points))/2]

    number = size(points, 1)
    times = [(i-1)/(number-1) for i in 1:number]

    function c(n)
        if n==0
            sum((points[i+1]-points[i])/(times[i+1]-times[i])*(times[i+1]*(times[i+1]/2-times[i])+(times[i]^2)/2)+points[i]*(times[i+1]-times[i]) for i in 1:(number-1))
        else
            sum(im*(points[i+1]-points[i])/(2*pi*n)*exp(-im*2*pi*n*times[i+1])+(im*points[i]/(2*pi*n)+(points[i+1]-points[i])/(times[i+1]-times[i])/(4*pi^2*n^2))*(exp(-im*2*pi*n*times[i+1])-exp(-im*2*pi*n*times[i])) for i in 1:(number-1))
        end
    end

    anim = @animate for N in convert.(Int, vcat(collect(1:convacc), convacc*ones(8,1)))
        tran = 0:0.005:1

        function fseries(x)
            sum(c(i)*exp(im*i*2*pi*x) for i in (-N):N)
        end

        iaxis = imag.(fseries.(tran))
        xaxis = real.(fseries.(tran))

        plot(xaxis, iaxis, color = :yellow, xlims = limx, ylims = limy, labels = "N = $N",
                background_color = :black, foreground_color = :white, grid = true, aspectratio = 1
            )
        scatter!(real.(points), imag.(points), color = :white, labels = false, marker = :cross)
    end

    gif(anim, "seriesconverge.gif", fps = 8)

    anim2 = @animate for time in vcat(collect(0:0.01:1), ones(30,1))
        evolve = 0:0.001:time

        function fseries2(x)
            sum(c(i)*exp(im*i*2*pi*x) for i in (-evolacc):evolacc)
        end

        ievol = imag.(fseries2.(evolve))
        xevol = real.(fseries2.(evolve))

        function arbitraryf(n)
            if iseven(n)==true
                a = -n/2
                b = n/2
            else
                a = 0.5-n/2
                b = 0.5+n/2
            end
            sum(c(i)*exp(im*i*2*pi*time) for i in a:b)
        end

        phs = vcat(0, [arbitraryf(i) for i in 0:16])

        plot(xevol, ievol, color = :yellow, xlims = limx, ylims = limy, labels = "N = $evolacc",
            background_color = :black, foreground_color = :white, grid = true, aspectratio = 1
            )
        #scatter!(real.(points), imag.(points), color = :white, labels = false, marker = :cross)
        plot!(real.(phs), imag.(phs), color = :white, labels = false, alpha = 0.65)
    end

    gif(anim2, "seriesevolution.gif", fps = 30)
end

@time idoit()

max_stars_repo_name: rxwop/fourierseries

max_issues_repo_name: rxwop/fourierseries

max_forks_repo_name: rxwop/fourierseries


text:
"""
Provides a definition for a SingleSwap mutation operator.
"""
type SingleSwapMutationDefinition <: MutationDefinition
  stage::AbstractString
  rate::Float

  SingleSwapMutationDefinition(stage::AbstractString, rate::Float) = new(stage, rate)
end

type SingleSwapMutation <: Mutation
  stage::AbstractString
  rate::Float
  representation::Representation

  SingleSwapMutation(stage::AbstractString, rate::Float, rep::Representation) =
    new(stage, rate, rep)
end

"""
Composes a single swap mutation operator from its definition.
"""
function compose!(d::SingleSwapMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  SingleSwapMutation(d.stage, d.rate, sp.stages[d.stage].representation)
end

"""
Swaps the position of two randomly selected genes within a provided chromosome
with one another.

**Parameters:**

* `stage::AbstractString`, the name of the developmental stage that this
  operator should be applied to. Defaults to the genotype if no stage is
  specified.
* `rate::Float`, the probability that this operator will perform a swap upon
   invocation.
"""
single_swap() = single_swap(0.5)
single_swap(rate::Float) = single_swap("", rate)
single_swap(stage::AbstractString) = single_swap(stage, 0.5)
single_swap(stage::AbstractString, rate::Float) =
  SingleSwapMutationDefinition(stage, rate)

"""
Possibly performs a single swap mutation operation on a provided chromosome,
with a probability given by the mutation rate of the underlying operator. If the
mutation is performed, the positions of two randomly selected genes are swapped
with one another.
"""
function mutate!{T}(o::SingleSwapMutation, input::Vector{T})
  rand() < o.rate && return input

  # Perform the swap between two randomly selected genes.
  i, j = sample(1:length(input), 2, replace=false)
  input[i], input[j] = input[j], input[i]

  return input
end


text:
"""
    step!(new, old, ::DistributionStrategy)

Perform one iteration of Conway's Game of Life.
"""
step!

"""
    prepare!(grid)

Prepare the grid to be used as the old/current world stage for `step!`.
Assuming that only inner cells, `grid[2:end-1,2:end-1]`, contain actual
world information, mirror the edges of that grid.
"""
function prepare!(grid)
    m, n = size(grid)
    @views begin
        grid[1,2:n-1] .= grid[m-1,2:n-1]
        grid[m,2:n-1] .= grid[2,2:n-1]
        # Now we are good to include the corners as well:
        grid[1:m,1] .= grid[1:m,n-1]
        grid[1:m,n] .= grid[1:m,2]
    end
    nothing
end

function life_rule!(new, old, newindices, oldindices)
    neighbours = [CartesianIndex(x,y) for x in -1:1, y in -1:1 if x!=0 || y!=0]
    for (i,j) in zip(CartesianIndices(oldindices), CartesianIndices(newindices))
        nc = sum(old[i+n] for n in neighbours)
        if old[i]
            new[j] = 2 <= nc <= 3
        else
            new[j] = nc == 3
        end
    end
    nothing
end

interior(r::OrdinalRange) = r[2:end-1]
interior(t::Tuple) = map(interior, t)
interior(m::AbstractMatrix) = @view m[interior(axes(m))...]


text:
using SafeTestsets

# The absolute most basic tests
@safetestset "Most Basic Tests" begin
    include("test_mostBasic.jl")
end

# Test MonteCarlo capabilities
@safetestset "MonteCarlo Tests" begin
    include("test_montecarlo.jl")
end

# The affine constraint tests
@safetestset "Affine Constraints Tests" begin
    include("test_affineConstraints.jl")
end

# The base constraint manager tests
@safetestset "Base Constraint Manager Tests" begin
    include("test_constraintManager.jl")
end

 # Test the angle cone
@safetestset "Angle SOCP Constraint Test" begin
    include("test_angleSOCP.jl")
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_volume

"""
    vtk_volume(;kwargs...)

A Volume component.
Volume is exposing a vtkImageData to a downstream filter
It takes the following set of properties:
  - state: { image: { ...imagedata-props }, field: { ...dataArray } }
Keyword arguments:
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `state` (Dict; optional): State of the volume
"""
function vtk_volume(; kwargs...)
        available_props = Symbol[:id, :port, :state]
        wild_props = Symbol[]
        return Component("vtk_volume", "Volume", "dash_vtk", available_props, wild_props; kwargs...)
end

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
function set_inactive!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer, tweet_list::AbstractArray
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    this_agent.active = false
    agent_edges = [e for e in edges(graph) if (src(e) == agent_idx || dst(e) == agent_idx)]
    for e in agent_edges
        rem_edge!(graph, e)
    end
    empty!(this_agent.feed)
    for t in tweet_list
        if t.source_agent == agent_idx
            t.weight = -1
        end
    end
    return state
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
#!/usr/bin/env julia

module VCF

export CheckChrmosome, records, header, version
export HeadInfos, chrange, getchr, fetch
export alts, uniqalts, info, infoAA
export E1000G, merge

include("read.jl")

"""
For checking the input chromosome arguement
"""
function CheckChrmosome(chr)
	match = 0
	possible = ["X"; "Y"; "M"; range(1,22)]
	for chrom in possible
		if string(chrom) == chr
			match += 1
		end
	end
	return match
end

"""
Function to extract only the vcf record containing lines
"""
function records(file)
	lines = read(file)
	record = []
	for ln in lines
		if ismatch(r"^[\d]|X|Y|M", ln)
			push!(record, "$ln")
		end
	end
	return record
end

"""
Function for extracting the header lines
"""
function header(lines)
	head = []
	for ln in lines
		if ismatch(r"^#", ln)
			push!(head, "$ln\n")
		end
	end
	if sizeof(head) >= 1
		return head
	else
		print("Header information not found in the input .vcf file!\n")
	end
end

"""
Function for getting Variant Call Format version information
"""
function version(lines)
	headers = header(lines)
	if sizeof(headers) < 1
		print("Header information not found in the input .vcf file!\n")
	elseif ismatch(r"=VCF", headers[1]) && sizeof(headers) >= 1
		head = headers[1]
		ln = split(head, "=")
		version = ln[2]
		return version
	else
		print("Version information not found in the input .vcf file!\n")
	end
end

"""
Function for extracting INFO provided in header lines
"""
function HeadInfos(lines)
	headers = header(lines)
	info = []
	if sizeof(headers) < 1
		print("Header information not found in the input .vcf file!\n")
	elseif sizeof(headers) >= 1
		for ln in headers
			if ismatch(r"^##INFO", ln)
				push!(info, "$ln\n")
			end
		end
		return info
	else
		print("INFO not found in the input .vcf file!\n")
	end
end

"""
Function for extracting start and end position of queried chromosome
"""
function chrange(lines, chr::AbstractString=chr)
	inpt_chr = CheckChrmosome(chr)
	if inpt_chr >= 1
		pos = []
		for rec in lines
			spl = split(rec, "\t")
			if spl[1] == chr
				push!(pos, spl[2])
				#println(spl[2])
			end
		end
		chr_strt = parse(Int, pos[1])
		chr_end   = parse(Int, pos[end])
		return chr_strt, chr_end
	else
		println("Please check the requested Chromosome name! \t *** $chr ***")
	end
end

function getchr(records)
	chrs = []
	for ln in records
		spl = split(ln, "\t")
		push!(chrs, spl[1])
	end
	list  = sort(unique(chrs))
	for ch in list
		print("$ch\n")
	end
end

"""
Function for extracting queried segment or entire chromosome
"""
function fetch(records, chr::AbstractString, strt=0, stop=0)
	inpt_chr = CheckChrmosome(chr)
	if inpt_chr >= 1
		pos1 = []
		# If specific range is not asked, then entire chomosome is selected
		if strt == 0 && stop == 0
			strt, stop = chrange(records, chr)
		end

		for rec2 in records
			spl2 = split(rec2, "\t")
			if spl2[1] == chr
				push!(pos1, spl2[2])
		#		println(spl2[2])
			end
		end

		pos_strt = parse(Int, pos1[1])
		pos_end   = parse(Int, pos1[end])

		for rec3 in records
			spl3 = split(rec3, "\t")
			cur_pos = parse(Int, spl3[2])
			region  = spl3[1]

			if (strt < pos_strt && stop > pos_end)|| strt < pos_strt || stop > pos_end || stop < pos_strt
				println("Request out of range!")
			elseif region == chr && strt <= cur_pos <= stop
				print(rec3)
			end
		end
	else
		println("Please check the requested Chromosome name! \t *** $chr ***")
	end
end

"""
Get unique pair of alteration records REF => ALT
"""
function alts(records)
	alt = []
	for line in records
		tmp = split(line, "\t")
		# if ID column present
		if ismatch(r"[\d]", tmp[3])
			#tmp[4] => REF ; $tmp[5] => ALT
			pair = tmp[4] *" => " * tmp[5]
			push!(alt, pair)
		else
			#tmp[3] => REF ; $tmp[4] => ALT
			pair = tmp[3] *" => " * tmp[4]
			push!(alt, pair)
		end
	end
	alt = unique(alt)
	return alt
end

"""
Get unique list of alteration (ALT) types
"""
function uniqalts(records)
	alt = []
	for line in records
		tmp = split(line, "\t")
		#
		 if ID column present
		if ismatch(r"[\d]", tmp[3])
			#$tmp[5] => ALT
			push!(alt, tmp[5])
		else
			#$tmp[4] => ALT
			push!(alt, tmp[4])
		end
	end
	alt = unique(alt)
	return alt
end

"""
INFO Fixed field from input VCF file
"""
function info(records)
	info = []
	for line in records
		tmp = split(line, "\t")
		#tmp[5] => ALT
		push!(info, tmp[8])
	end
	info = unique(info)
	return info
end

"""
AA : ancestral allele
"""
function infoAA(records)
	inf = info(records)
	aa = []
	for line in inf
		tmp = split(line, ";")
		for col in tmp
			#tmp[end] => AA
			if ismatch(r"AA=", col)
				push!(aa, col)
			end
		end
	end
	return aa
end

"""
1000G : membership in 1000 Genomes
"""
function E1000G(records)
	tg = []
	for line in records
		if ismatch(r"E_1000G", line)
			push!(tg, line)
		end
	end
	return tg
end

function merge(vcfiles)

end

end
#function contigs
#function filters
#function formats
#function metadata


text:
using EarthAlbedo
using Test

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


text:
function dual_model_and_map(primal_model::MOI.ModelLike)
    dual = Dualization.dualize(primal_model)
    return dual.dual_model, dual.primal_dual_map
end


text:
a = 3

b = 4

function f(x, y)
  a = 3
  return x + y
end

c = f(a, b)

d = 5

println(c)

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
function yLP(S, b, lb, ub, c, d;
        ϵ = 1e-5,
        sense = JuMP.MOI.MAX_SENSE,
        solver = Clp.Optimizer
    )

    lp_model = JuMP.Model(solver)
    JuMP.set_silent(lp_model)
    M, N = size(S)

    # Variables
    y = JuMP.@variable(lp_model, y[1:N])
    r = JuMP.@variable(lp_model, r)

    # Constraints
    JuMP.@constraint(lp_model, S * y - r .* b .== 0.0)
    JuMP.@constraint(lp_model, d' * y == 1.0)
    JuMP.@constraint(lp_model, y - r .* lb .>= 0.0)
    JuMP.@constraint(lp_model, y - r .* ub .<= 0.0)
    JuMP.@constraint(lp_model, r >= ϵ)

    # objective
    JuMP.@objective(lp_model, sense, c' * y)

    # optimize
    JuMP.optimize!(lp_model)

    yval = JuMP.value.(y)
    rval = JuMP.value.(r)
    v = yval ./ rval # sol
    y = (c' * v) / (d' * v) # yield

    status = JuMP.termination_status(lp_model)
    status, v, y
end
yLP(net, d; kwargs...) = yLP(net.S, net.b, net.lb, net.ub, net.c, d; kwargs...)


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using Random
using SparseArrays
using ADCMEKit

Random.seed!(118)

# trial number specifies a folder to store results
trialnum = 6

# geometry setup in domain [0,1]^2
solid_left = 0.4
solid_right = 0.6
solid_top = 0.5
solid_bottom = 0.6

chip_left = 0.46
chip_right = 0.54
chip_top = 0.5
chip_bottom = 0.52

# solid_left = 0.45
# solid_right = 0.55
# solid_top = 0.5
# solid_bottom = 0.52

# chip_left = 0.48
# chip_right = 0.52
# chip_top = 0.5
# chip_bottom = 0.505

k_mold = 0.014531
k_chipdie = 2.60475
k_air = 0.64357
nu = 0.47893 # equal to 1/Re
power_source = 82.46295 #82.46295 = 1.0e6 / air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

epsilon_k = 0.3

function k_exact(x, y)
    # 3.0 + 100000 * (x - 0.5)^3 / (1 + y^2)
    # 5 * exp((-(x-0.5)^2-(y-0.5)^2)/ 0.00002) +
    # 3 * exp((-(x-0.48)^2-(y-0.505)^2)/ 0.00005) +
    # 6 * exp((-(x-0.51)^2-(y-0.502)^2)/ 0.00001) + 2.604
    k_mold + (rand() < 0.5) * ((1 - epsilon_k) + rand() * 2 * epsilon_k) * k_chipdie

end

function k_nn(x, y)
    out = fc([x y], [20,20,20,1])^2 + 0.01 # N x 1
    squeeze(out)
end

m = 50
n = 50
h = 1/n
NT = 10    # number of iterations for Newton's method, 8 is good for m=400

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
chip_fem_idx = Array{Int64, 1}([])
chip_fvm_idx = Array{Int64, 1}([])
chip_fem_top_idx = Array{Int64, 1}([])

for i = 1:(m+1)
    for j = 1:(n+1)
        if (i-1)*h >= solid_left-1e-9 && (i-1)*h <= solid_right+1e-9 && (j-1)*h >= solid_top-1e-9 && (j-1)*h <= solid_bottom+1e-9
            # print(i, j)
            global solid_fem_idx = [solid_fem_idx; (j-1)*(m+1)+i]
            if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_bottom+1e-9
                global chip_fem_idx = [chip_fem_idx; (j-1)*(m+1)+i]
            end
            if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_top+1e-9
                global chip_fem_top_idx = [chip_fem_top_idx; (j-1)*(m+1)+i]
            end
        end
    end
end

for i = 1:m
    for j = 1:n
        if (i-1)*h + h/2 >= solid_left-1e-9 && (i-1)*h + h/2 <= solid_right+1e-9 &&
            (j-1)*h + h/2 >= solid_top-1e-9 && (j-1)*h + h/2 <= solid_bottom+1e-9
            global solid_fvm_idx = [solid_fvm_idx; (j-1)*m+i]
            if (i-1)*h + h/2 >= chip_left-1e-9 && (i-1)*h + h/2 <= chip_right+1e-9 && (j-1)*h + h/2 >= chip_top-1e-9 && (j-1)*h + h/2<= chip_bottom+1e-9
                global chip_fvm_idx = [chip_fvm_idx; (j-1)*m+i]
            end
        end
    end
end

# initialize space varying k and heat source

xy = fem_nodes(m, n, h)
chip_x, chip_y = xy[chip_fem_idx, 1], xy[chip_fem_idx, 2]
# k_chip = @. k_nn(chip_x, chip_y); k_chip=stack(k_chip)
k_chip = Variable(k_mold .* ones(length(chip_fem_idx)))
k_chip_exact = @. k_exact(chip_x, chip_y)

k_fem = k_air * constant(ones((m+1)*(n+1)))
k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
k_fem = scatter_update(k_fem, chip_fem_idx, k_chip)
kgauss = fem_to_gauss_points(k_fem, m, n, h)

heat_source_fem = zeros((m+1)*(n+1))
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = fem_to_gauss_points(heat_source_fem, m, n, h)

B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
heat_source = constant(compute_fem_source_term1(heat_source_gauss, m, n, h))

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))

bd = bcnode("all", m, n, h)

# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

bd = [bd; bd .+ (m+1)*(n+1);
     2*(m+1)*(n+1)+1; 2*(m+1)*(n+1)+m;
    #  (2*(m+1)*(n+1)+m*n )+1:(2*(m+1)*(n+1)+m*n )+m+1]
     bd .+ (2*(m+1)*(n+1)+m*n )]

# add solid region into boundary condition for u, v, p
bd = [bd; solid_fem_idx; solid_fem_idx .+ (m+1)*(n+1); solid_fvm_idx .+ 2(m+1)*(n+1)]

function compute_residual(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = fem_to_gauss_points(T, m, n, h)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, m, n, h)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, m, n, h) * T - heat_source
    R = [F;G;H0;T0]
    return R
end

function compute_jacobian(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]

    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, m, n, h))

    gradT = eval_grad_on_gauss_pts1(T, m, n, h)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = -buoyance_coef * constant(compute_fem_mass_matrix1(m, n, h))
    T_mat = [SparseTensor(spzeros((m+1)*(n+1), (m+1)*(n+1))); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(m*n, (m+1)*(n+1)))]

    N = 2*(m+1)*(n+1) + m*n
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros((m+1)*(n+1), m*n)) M]]
end

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

function plot_velo_pres_temp_cond(k)

    S_true = V_data
    S = run(sess, V_computed)

    figure(figsize=(15,4))
    subplot(131)
    title("exact x velocity")
    visualize_scalar_on_fem_points(S_true[1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    subplot(132)
    title("predicted x velocity")
    visualize_scalar_on_fem_points(S[1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    subplot(133)
    title("difference in x velocity")
    visualize_scalar_on_fem_points(S[1:(m+1)*(n+1)] .* u_std .- S_true[1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    tight_layout()
    savefig("xzchip_figures$trialnum/pix/xzchipv0_pix_velox$k.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact y velocity")
    visualize_scalar_on_fem_points(S_true[(m+1)*(n+1)+1: 2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    subplot(132)
    title("predicted y velocity")
    visualize_scalar_on_fem_points(S[(m+1)*(n+1)+1: 2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    subplot(133)
    title("difference in y velocity")
    visualize_scalar_on_fem_points(S[(m+1)*(n+1)+1: 2*(m+1)*(n+1)]  .* u_std .- S_true[(m+1)*(n+1)+1: 2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
    tight_layout()
    savefig("xzchip_figures$trialnum/pix/xzchipv0_pix_veloy$k.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact pressure")
    visualize_scalar_on_fvm_points(S_true[ 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std, m, n, h);gca().invert_yaxis()
    subplot(132)
    title("predicted pressure")
    visualize_scalar_on_fvm_points(S[ 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n]  .* p_std, m, n, h);gca().invert_yaxis()
    subplot(133)
    title("difference in pressure")
    visualize_scalar_on_fvm_points(S[ 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std .- S_true[ 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std,  m, n, h);gca().invert_yaxis()
    tight_layout()
    savefig("xzchip_figures$trialnum/pix/xzchipv0_pix_pres$k.png")

    figure(figsize=(15,4))
    subplot(131)
    title("exact temperature")
    visualize_scalar_on_fem_points(S_true[ 2*(m+1)*(n+1)+m*n+1:end] .* T_infty .+ T_infty, m, n, h);gca().invert_yaxis()
    subplot(132)
    title("predicted temperature")
    visualize_scalar_on_fem_points(S[ 2*(m+1)*(n+1)+m*n+1:end] .* T_infty .+ T_infty, m, n, h);gca().invert_yaxis()
    subplot(133)
    title("difference in temperature")
    visualize_scalar_on_fem_points(S[ 2*(m+1)*(n+1)+m*n+1:end]  .* T_infty .- S_true[2*(m+1)*(n+1)+m*n+1:end] .* T_infty, m, n, h);gca().invert_yaxis()
    tight_layout()
    savefig("xzchip_figures$trialnum/pix/xzchipv0_pix_temp$k.png")

    figure(figsize=(15,4))
    subplot(131)
    xx = chip_left : h : chip_right
    yy = chip_top : h : chip_bottom
    k_chip_exact_2d = reshape(k_chip_exact, length(yy), length(xx))
    pcolor(xx, yy, k_chip_exact_2d)
    gca().set_aspect(1)
    colorbar()
    title("exact chip conductivity")

    subplot(132)
    k_chip_ = run(sess, k_chip)
    k_chip_2d = reshape(k_chip_, length(yy), length(xx))
    pcolor(xx, yy, k_chip_2d)
    gca().set_aspect(1)
    colorbar()
    title("predicted chip conductivity")

    subplot(133)
    pcolor(xx, yy, k_chip_2d .- k_chip_exact_2d)
    gca().set_aspect(1)
    colorbar()
    title("difference in chip conductivity")

    tight_layout()
    savefig("xzchip_figures$trialnum/pix/xzchipv0_pix_cond$k.png")

end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(m*n+3*(m+1)*(n+1)))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n+(m+1)*(n+1)))

# sess = Session(); init(sess)
# output = run(sess, S)
V_computed = S[end, :]

V_data = matread("xzchip_figures$trialnum/xzchipv0_fn_data_m$(m)_n$n.mat")["V"]

# sample_size = 100
# idx = rand(1:(m+1)*(n+1), sample_size)
# idx = [idx; (m+1)*(n+1) .+ idx; 2*(m+1)*(n+1)+m*n .+ idx] # observe velocity and temperature
# observed_data = V_data[idx]

noise = false
noise_level = 0.05
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

# loss = mean((V_computed[idx] .- observed_data)^2)
loss = mean((V_computed .- V_data)^2)
loss = loss * 1e10
# ---------------------------------------------------
# create a session and run
max_iter = 1
sess = Session(); init(sess)

@info run(sess, loss, k_chip=>k_chip_exact)
lineview(sess, k_chip, loss, k_chip_exact, one.(k_chip_exact))
savefig("xzchip_figures$trialnum/lineview.png")

# for k = 1:100
#     loss_ = BFGS!(sess, loss, max_iter)
#     matwrite("xzchip_figures5/loss$k.mat", Dict("L"=>loss_))
#     close("all"); semilogy(loss_); title("loss vs. iteration")
#     savefig("xzchip_figures5/loss$k.png")
#     plot_velo_pres_temp_cond(k)
#     ADCME.save(sess, "xzchip_figures5/nn$k.mat")
# end


text:
export ft_settimeouts

function ft_settimeouts(ft_handle::Culong,
                         readtimeout::Integer,
                         writetimeout::Integer)
  ft_status = ccall((:FT_SetTimeouts, d2xx),
                     Cuint,
                     (Culong, Cuint, Cuint),
                     ft_handle, readtimeout, writetimeout)
  checkstatus(ft_status)
  return nothing
end


text:
import Term.consoles: console_height, console_width, console, err_console, Console
import Term: tprint, tprintln

println("\nTesting tprint, stdout temporarily disabled")

@suppress_out begin
    @testset "\e[34mTPRINT" begin
        @test_nothrow tprint("string")
        @test_nothrow tprint("[red]adasd[/red]")
        @test_nothrow tprint("[blue on_green]adasda")
        @test_nothrow tprint("[red]dadas[green]insdai[/green]outssdrse[blue]fsfsf[/blue]")

        @test_nothrow tprint(Panel("test"))
        @test_nothrow tprint(TextBox("test"))

        @test_nothrow tprint(1)
        @test_nothrow tprint(:x)

        @test_nothrow tprint(1, Panel("test"), "test")

        @test_nothrow tprintln("test")

        @test_nothrow tprintln(Panel("test"))

        @test_nothrow tprintln(1, Panel("test"), "test")
    end
end

@testset "\e[34mCONSOLE" begin
    @test console_height() == displaysize(stdout)[1]
    @test console_height(stdout) == displaysize(stdout)[1]

    @test console_width() == displaysize(stdout)[2]
    @test console_width(stdout) == displaysize(stdout)[2]

    @test Console(stdout) == Console()

end


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @load

# load training set -
full_training_data_frame = load_tga_training_data()

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> ((x == 2 || x == 3) && y == has_TM_flag), full_training_data_frame)
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:14]))
validation_data_array = convert.(Float32, Matrix(experimental_data_table[!, 15:19]))

# what is the size of the system?
(P, D) = size(experimental_data_table)

# setup sampling -
δ = 0.9
number_of_samples = 100
number_of_input_types = 12
input_range_array = Array{Float32,2}(undef, number_of_input_types, 2)
for i ∈ 1:number_of_input_types
    L = (1 - δ) * minimum(input_data[:, i])
    U = (1 + δ) * maximum(input_data[:, i])
    input_range_array[i, 1] = L
    input_range_array[i, 2] = U
end

# archive -
ensemble_archive = Array{Matrix{Float32},1}()

# build an input array -
sample_input_array_random = Array{Float32,2}(undef, number_of_samples, number_of_input_types)
for i ∈ 1:number_of_samples
    for j ∈ 1:number_of_input_types
        L = input_range_array[j, 1]
        U = input_range_array[j, 2]
        r = rand()
        v = L + r * (U - L)
        sample_input_array_random[i, j] = v

        if (j == number_of_input_types)
            sample_input_array_random[i, j] = has_TM_flag
        end
    end
end

# main training loop -
for i ∈ 1:P

    # load a model -
    model_name = "deep_coag_model-L$(i)O-TF-TM-$(has_TM_flag)-V2-V3.bson"
    model_file_path = joinpath(_PATH_TO_MODELS, "network-1-epoch-12k", model_name)
    @load model_file_path deep_coag_model

    # sample -
    simulated_output_data = sample(deep_coag_model, sample_input_array_random)
    push!(ensemble_archive, simulated_output_data)
end

S = vcat(ensemble_archive...)

# save -
# path_to_input = joinpath(_PATH_TO_DATA, "synthetic_input_array.dat")
# CSV.write(path_to_input, Tables.table(sample_input_array_random))

# path_to_output = joinpath(_PATH_TO_DATA, "synthetic_output_array.dat")
# CSV.write(path_to_output, Tables.table(S))

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
mutable struct BadAgent <: AbstractAgent
    useless::Int
    pos::Int
end
mutable struct BadAgentId <: AbstractAgent
    id::Float64
end
struct ImmutableAgent <: AbstractAgent
    id::Int
end
mutable struct DiscreteVelocity <: AbstractAgent
    id::Int
    pos::NTuple{2, Float64}
    vel::NTuple{2, Int}
    diameter::Float64
end
mutable struct ParametricAgent{T <: Integer} <: AbstractAgent
    id::T
    pos::NTuple{2, T}
    weight::T
    info::String
end

@testset "Model construction" begin
    # Shouldn't use ImmutableAgent since it cannot be edited
    agent = ImmutableAgent(1)
    @test_logs (:warn, "AgentType should be mutable. Try adding the `mutable` keyword infront of `struct` in your agent definition.") ABM(agent)
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(agent; warn=false)) <: AbstractAgent
    # Cannot use BadAgent since it has no `id` field
    @test_throws ArgumentError ABM(BadAgent)
    agent = BadAgent(1,1)
    @test_throws ArgumentError ABM(agent)
    # Cannot use BadAgent in a grid space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(BadAgent, GridSpace((1,1)))
    @test_throws ArgumentError ABM(agent, GridSpace((1,1)))
    # Cannot use BadAgentId since `id` has an invalid type
    @test_throws ArgumentError ABM(BadAgentId)
    agent = BadAgentId(1.0)
    @test_throws ArgumentError ABM(agent)
    # Cannot use Agent0 in a grid space context since it has no `pos` field
    @test_throws ArgumentError ABM(Agent0, GridSpace((1,1)))
    agent = Agent0(1)
    @test_throws ArgumentError ABM(agent, GridSpace((1,1)))
    # Cannot use Agent3 in a graph space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(Agent3, GraphSpace(Agents.Graph(1)))
    agent = Agent3(1, (1,1), 5.3)
    @test_throws ArgumentError ABM(agent, GraphSpace(Agents.Graph(1)))
    # Cannot use Agent3 in a continuous space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(Agent3, ContinuousSpace((1,1), 0.1))
    @test_throws ArgumentError ABM(agent, ContinuousSpace((1,1), 0.1))
    # Cannot use Agent4 in a continuous space context since it has no `vel` field
    @test_throws ArgumentError ABM(Agent4, ContinuousSpace((1,1), 0.1))
    agent = Agent4(1, (1,1), 5)
    @test_throws ArgumentError ABM(agent, ContinuousSpace((1,1), 0.1))
    # Shouldn't use DiscreteVelocity in a continuous space context since `vel` has an invalid type
    @test_logs (:warn, "`vel` field in Agent struct should be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace.") ABM(DiscreteVelocity, ContinuousSpace((1,1), 0.1))
    agent = DiscreteVelocity(1, (1,1), (2,3), 2.4)
    @test_logs (:warn, "`vel` field in Agent struct should be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace.") ABM(agent, ContinuousSpace((1,1), 0.1))
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(agent, ContinuousSpace((1,1), 0.1); warn=false)) <: AbstractAgent
    # Shouldn't use ParametricAgent since it is not a concrete type
    @test_logs (:warn, "AgentType is not concrete. If your agent is parametrically typed, you're probably seeing this warning because you gave `Agent` instead of `Agent{Float64}` (for example) to this function. You can also create an instance of your agent and pass it to this function. If you want to use `Union` types for mixed agent models, you can silence this warning.") ABM(ParametricAgent, GridSpace((1,1)))
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(ParametricAgent, GridSpace((1,1)); warn=false)) <: AbstractAgent
    # ParametricAgent{Int} is the correct way to use such an agent
    @test Agents.agenttype(ABM(ParametricAgent{Int}, GridSpace((1,1)))) <: AbstractAgent
    #Type inferance using an instance can help users here
    agent = ParametricAgent(1, (1,1), 5, "Info")
    @test Agents.agenttype(ABM(agent, GridSpace((1,1)))) <: AbstractAgent
    #Mixed agents
    @test Agents.agenttype(ABM(Union{Agent0,Agent1}; warn=false)) <: AbstractAgent
    @test_logs (:warn, "AgentType is not concrete. If your agent is parametrically typed, you're probably seeing this warning because you gave `Agent` instead of `Agent{Float64}` (for example) to this function. You can also create an instance of your agent and pass it to this function. If you want to use `Union` types for mixed agent models, you can silence this warning.") ABM(Union{Agent0,Agent1})
    @test_throws ArgumentError ABM(Union{Agent0,BadAgent}; warn=false)
end

@testset "sample!" begin
  Random.seed!(6459)
  model = ABM(Agent2)
  for i in 1:20; add_agent!(model, rand()/rand()); end
  allweights = [i.weight for i in values(model.agents)]
  mean_weights = sum(allweights)/length(allweights)
  sample!(model, 12, :weight)
  @test Agents.nagents(model) == 12
  allweights = [i.weight for i in values(model.agents)]
  mean_weights_new = sum(allweights)/length(allweights)
  @test mean_weights_new > mean_weights
  sample!(model, 40, :weight)
  @test Agents.nagents(model) == 40
  allweights = [i.weight for i in values(model.agents)]
  mean_weights_new = sum(allweights)/length(allweights)
  @test mean_weights_new > mean_weights

  Random.seed!(6459)
  model2 = ABM(Agent3, GridSpace((10, 10)))
  for i in 1:20; add_agent_single!(Agent3(i, (1,1), rand()/rand()), model2); end
  @test sample!(model2, 10, :weight) == nothing
  allweights = [i.weight for i in values(model2.agents)]
  mean_weights = sum(allweights)/length(allweights)
  sample!(model2, 12, :weight)
  @test Agents.nagents(model2) == 12
  allweights = [i.weight for i in values(model2.agents)]
  mean_weights_new = sum(allweights)/length(allweights)
  @test mean_weights_new > mean_weights

  sample!(model2, 40, :weight)
  @test Agents.nagents(model2) == 40

  Random.seed!(6459)
  #Guarantee all starting weights are unique
  model3 = ABM(Agent2)
  while true
    for i in 1:20; add_agent!(model3, rand()/rand()); end
    allweights = [i.weight for i in values(model3.agents)]
    allunique(allweights) && break
  end
  # Cannot draw 50 samples out of a pool of 20 without replacement
  @test_throws ErrorException sample!(model3, 50, :weight; replace=false)
  sample!(model3, 15, :weight; replace=false)
  allweights = [i.weight for i in values(model3.agents)]
  @test allunique(allweights)
  model3 = ABM(Agent2)
  while true
    for i in 1:20; add_agent!(model3, rand()/rand()); end
    allweights = [i.weight for i in values(model3.agents)]
    allunique(allweights) && break
  end
  sample!(model3, 100, :weight; replace=true)
  allweights = [i.weight for i in values(model3.agents)]
  @test !allunique(allweights)
end

@testset "add_agent! (discrete)" begin
  properties = Dict(:x1=>1)
  space = GraphSpace(complete_digraph(10))
  model = ABM(Agent7, space; properties=properties)
  attributes = (f1=true,f2=1)
  add_agent!(1, model, attributes...)
  attributes = (f2=1,f1=true)
  add_agent!(1, model; attributes...)
  @test model.agents[1].id != model.agents[2].id
  @test model.agents[1].pos == model.agents[2].pos
  @test model.agents[1].f1 == model.agents[2].f1
  @test model.agents[1].f2 == model.agents[2].f2
  @test add_agent_single!(model, attributes...).pos ∈ 1:10
  for id in 4:11
      agent = Agent7(id, 2, attributes...)
      add_agent_single!(agent, model)
  end
  @test !has_empty_positions(model)
  agent = Agent7(12,5, attributes...)
  add_agent_single!(agent, model)
  @test_throws KeyError model[12]
  add_agent!(agent, model)
  @test model[12].pos ∈ 1:10

  agent = Agent7(13,5, attributes...)
  @test add_agent!(agent, 3, model).pos == 3

  model = ABM(Agent1, GridSpace((10,10)))
  agent = Agent1(1,(3,6))
  @test add_agent!(agent, (7,8), model).pos == (7,8)
end

@testset "add_agent! (continuous)" begin
  properties = Dict(:x1=>1)
  space2d = ContinuousSpace((1,1), 0.1; periodic=true)
  model = ABM(Agent8, space2d; properties=properties)
  attributes = (f1=true,f2=1)
  add_agent!(model, attributes...)
  attributes = (f2=1,f1=true)
  add_agent!(model; attributes...)
  @test model.agents[1].id != model.agents[2].id
  @test model.agents[1].f1 == model.agents[2].f1
  @test model.agents[1].f2 == model.agents[2].f2
  agent = Agent8(3, (0,0), false, 6)
  @test 0 <= add_agent!(agent, model).pos[1] <= 1
  agent.id = 4
  @test add_agent!(agent, (0.5, 0.5), model).pos[1] ≈ 0.5 atol=1e-3
end

@testset "move_agent!" begin
  # GraphSpace
  model = ABM(Agent5, GraphSpace(path_graph(6)))
  agent = add_agent!(model, 5.3)
  init_pos = agent.pos
  # Checking specific indexing
  move_agent!(agent, rand([i for i in 1:6 if i != init_pos]), model)
  new_pos = agent.pos
  @test new_pos != init_pos
  # Checking a random move
  ni = 0; init_pos = agent.pos
  while agent.pos == init_pos
      move_agent!(agent, model)
  end
  @test ni < Inf

  # GridSpace
  model = ABM(Agent1, GridSpace((5,5)))
  agent = add_agent!((2,4), model)
  move_agent!(agent, (1,3), model)
  @test agent.pos == (1,3)
  ni = 0; init_pos = agent.pos
  while agent.pos == init_pos
      move_agent!(agent, model)
  end
  @test ni < Inf

  model = ABM(Agent1, GridSpace((2,1)))
  agent = add_agent!((1,1), model)
  move_agent_single!(agent, model)
  @test agent.pos == (2,1)
  agent2 = add_agent!((1,1), model)
  move_agent_single!(agent2, model)
  # Agent shouldn't move since the grid is saturated
  @test agent2.pos == (1,1)

  # ContinuousSpace
  model = ABM(Agent6, ContinuousSpace((1,1), 0.1))
  agent = add_agent!((0.0, 0.0), model, (0.5, 0.0), 1.0)
  move_agent!(agent, model)
  @test agent.pos == (0.5, 0.0)
end

@testset "kill_agent!" begin
  # No Space
  model = ABM(Agent0)
  add_agent!(model)
  agent = add_agent!(model)
  @test nagents(model) == 2
  kill_agent!(agent, model)
  @test nagents(model) == 1
  # GraphSpace
  model = ABM(Agent5, GraphSpace(path_graph(6)))
  add_agent!(model, 5.3)
  add_agent!(model, 2.7)
  @test nagents(model) == 2
  kill_agent!(model.agents[1], model)
  @test nagents(model) == 1
  kill_agent!(2, model)
  @test nagents(model) == 0
  # GridSpace
  model = ABM(Agent1, GridSpace((5,5)))
  add_agent!((1,3), model)
  add_agent!((1,3), model)
  add_agent!((5,2), model)
  @test nagents(model) == 3
  for id in copy(ids_in_position((1,3), model))
    kill_agent!(id, model)
  end
  @test nagents(model) == 1
  # ContinuousSpace
  model = ABM(Agent6, ContinuousSpace((1,1), 0.1))
  add_agent!((0.7,0.1), model, (15,20), 5.0)
  add_agent!((0.2,0.9), model, (8,35), 1.7)
  @test nagents(model) == 2
  kill_agent!(model[1], model)
  @test nagents(model) == 1
end

@testset "genocide!" begin
  # Testing no space
  model = ABM(Agent0)
  for i in 1:10
      a = Agent0(i)
      add_agent!(a, model)
  end
  genocide!(model)
  @test nagents(model) == 0
  for i in 1:10
      a = Agent0(i)
      add_agent!(a, model)
  end
  genocide!(model, 5)
  @test nagents(model) == 5
  genocide!(model, a -> a.id < 3)
  @test nagents(model) == 3

  model = ABM(Agent3, GridSpace((10, 10)))

  # Testing genocide!(model::ABM)
  for i in 1:20
    agent = Agent3(i, (1,1), rand())
    add_agent_single!(agent, model)
  end
  genocide!(model)
  @test nagents(model) == 0

  # Testing genocide!(model::ABM, n::Int)
  for i in 1:20
    # Explicitly override agents each time we replenish the population,
    # so we always start the genocide with 20 agents.
    agent = Agent3(i, (1,1), rand())
    add_agent_single!(agent, model)
  end
  genocide!(model, 10)
  @test nagents(model) == 10

  # Testing genocide!(model::ABM, f::Function) with an anonymous function
  for i in 1:20
    agent = Agent3(i, (1,1), rand())
    add_agent_single!(agent, model)
  end
  @test nagents(model) == 20
  genocide!(model, a -> a.id > 5)
  @test nagents(model) == 5

  Random.seed!(6465)
  # Testing genocide!(model::ABM, f::Function) when the function is invalid
  # (i.e. does not return a bool)
  for i in 1:20
    agent = Agent3(i, (rand(1:10), rand(1:10)), i*2)
    add_agent_pos!(agent, model)
  end
  @test_throws TypeError genocide!(model, a -> a.id)
  N = nagents(model)

  # Testing genocide!(model::ABM, f::Function) with a named function
  # No need to replenish population since the last test fails
  function complex_logic(agent::A) where A <: AbstractAgent
    if agent.pos[1] <= 5 && agent.weight > 25
      true
    else
      false
    end
  end
  genocide!(model, complex_logic)
  @test nagents(model) < N

  space2d = ContinuousSpace((1,1), 0.1; periodic=true)
  model = ABM(Agent8, space2d)
  attributes = (f1=true,f2=1)
  for _ in 1:10 add_agent!(model, attributes...) end
  genocide!(model)
  @test nagents(model) == 0
  for _ in 1:10 add_agent!(model, attributes...) end
  genocide!(model, 5)
  @test nagents(model) == 5
  genocide!(model, a -> a.id < 3)
  @test nagents(model) == 3
end

mutable struct Daisy <: AbstractAgent
  id::Int
  pos::Tuple{Int, Int}
  breed::String
end
mutable struct Land <: AbstractAgent
  id::Int
  pos::Tuple{Int, Int}
  temperature::Float64
end
@testset "fill space" begin
  space = GridSpace((10, 10))
  model = ABM(Land, space)
  fill_space!(model, 15)
  @test nagents(model) == 100
  for a in allagents(model)
    @test a isa Land
    @test a.temperature == 15
  end

  space = GridSpace((10, 10))
  model = ABM(Union{Daisy, Land}, space; warn = false)
  fill_space!(Daisy, model, "black")
  @test nagents(model) == 100
  for a in values(model.agents)
    @test a isa Daisy
    @test a.breed == "black"
  end

  space = GridSpace((10, 10), periodic = true)
  model = ABM(Union{Daisy, Land}, space; warn = false)
  temperature(pos) = (pos[1]/10, ) # make it Tuple!
  fill_space!(Land, model, temperature)
  @test nagents(model) == 100
  for a in values(model.agents)
    @test a.temperature == a.pos[1]/10
  end

end

@testset "model step order" begin
    function model_step!(model)
        for a in allagents(model)
            if a.weight > 1.0
                model.count += 1
            end
        end
    end
    function agent_step!(a, model)
        a.weight += 1
    end

    for bool in (true, false)
        model = ABM(Agent2; properties = Dict(:count => 0))
        for i in 1:100
            add_agent!(model, rand())
        end
        step!(model, agent_step!, model_step!, 1, bool)
        if bool
            @test model.count == 100
        else
            @test model.count == 0
        end
    end
end


text:
function get_tennisracket(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    T=Float64)

    origin = Origin{T}(name=:origin)
    mass = 1.0
    h = 1.0
    bodies = [Box(h / 25.0, h / 2.0, h, mass, color=RGBA(1.0, 0.0, 0.0), name=:box)]
    joints = [JointConstraint(Floating(origin, bodies[1]), name=:floating_joint)]
    mechanism = Mechanism(origin, bodies, joints,
        timestep=timestep,
        gravity=gravity)
    return mechanism
end

function initialize_tennisracket!(mechanism::Mechanism{T};
    x=zeros(3),
    q=one(Quaternion),
    v=zeros(3),
    ω=zeros(3)) where T

    joint = get_joint(mechanism, :floating_joint)
    zero_velocity!(mechanism)
    set_minimal_coordinates!(mechanism, joint, [x; rotation_vector(q)])
    set_minimal_velocities!(mechanism, joint, [v; ω])
end


text:
using Base.Test
using Lora

# println("    Testing BasicMCJob constructors...")

# Example 1: multivariate parameter

p = BasicContMuvParameter(
  :p, 1, signature=:low, logtarget=(state, states) -> state.logtarget = -dot(state.value, state.value)
)
model = likelihood_model([p])

sampler = MH([1., 1.])

tuner = VanillaMCTuner()
# tuner = VanillaMCTuner(verbose=true)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

vstate = [BasicContMuvParameterState([1.25, 3.11], [:value, :logtarget], [:accept])]

outopts = Dict(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

job = BasicMCJob(
  model,
  1,
  sampler,
  tuner,
  mcrange,
  vstate,
  outopts,
  true,
  true,
  false,
  false
)

@time run(job)

chain = output(job)

mean(chain)

acceptance(chain)

# Example 2: univariate parameter

# using Lora

p = BasicContUnvParameter(:p, 1, signature=:low, logtarget=(state, states) -> state.logtarget = -abs2(state.value))
model = likelihood_model([p])

sampler = MH()

tuner = VanillaMCTuner()

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

vstate = [BasicContUnvParameterState(5.1)]

outopts = Dict(:monitor=>[:value, :logtarget])

job = BasicMCJob(
  model,
  1,
  sampler,
  tuner,
  mcrange,
  vstate,
  outopts,
  true,
  true,
  false,
  false
)

@time run(job)

chain = output(job)

mean(chain)

acceptance(chain, diagnostics=false)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Basic goto tests

function goto_test1()
    @goto a
    return false
    @label a
    return true
end
@test goto_test1()

@test eval(:(@label a)) === nothing

@test Expr(:error, "label \"a\" referenced but not defined") ==
    Meta.lower(@__MODULE__, :(@goto a))

@test Expr(:error, "label \"a\" defined multiple times") ==
    Meta.lower(@__MODULE__, quote
        function goto_test2()
            @goto a
            @label a
            @label a
            return
        end
    end)

@test Expr(:error, "label \"a\" referenced but not defined") ==
    Meta.lower(@__MODULE__, quote
        function goto_test3()
            @goto a
            return
        end
    end)

@test Expr(:error, "misplaced label") ==
    Meta.lower(@__MODULE__, quote
        function goto_test4()
            @goto a
            try
                @label a
            catch
            end
        end
    end)

# test that labels in macros are reassigned if unescaped
macro goto_test5_macro1()
    return :(@label a)
end
macro goto_test5_macro2()
    return :(@goto a)
end
macro goto_test5_macro3()
    return esc(:(@label a))
end

@test Expr(:error, "label \"a\" referenced but not defined") ==
    Meta.lower(@__MODULE__, quote
        function goto_test5_1()
            @goto a
            @goto_test5_macro1
            return
        end
    end)

let e = Meta.lower(@__MODULE__, quote
        function goto_test5_2()
            @goto_test5_macro2
            @label a
            return
        end
    end)
    @test (e::Expr).head === :error
    @test contains(e.args[1], r"label \"#\d+#a\" referenced but not defined")
end

function goto_test5_3()
    @goto a
    return false
    @goto_test5_macro3
    return true
end
@test goto_test5_3()

@test Expr(:error, "goto from a try/finally block is not permitted") ==
    Meta.lower(@__MODULE__, quote
        function goto_test6()
            try
                @goto a
            finally
            end
            @label a
            return
        end
    end)

function goto_test6()
    @goto a
    @label a
end

@test goto_test6() === nothing

function goto_test7(x)
    @label a
    if x
        @goto a
    end
end

@test goto_test7(false) === nothing

module GotoMacroTest
    macro goto_test8_macro()
        quote
            function $(esc(:goto_test8))()
                @label a
                @goto a
            end
        end
    end
end

GotoMacroTest.@goto_test8_macro

# issue #15600
function t0_15600(flag)
    flag && @goto return2
    return 1
    @label return2
    return 2
end
@test t0_15600(true) == 2
@test t0_15600(false) == 1
function t1_15600(flag)
    flag || @goto return2
    return 1
    @label return2
    return 2
end
@test t1_15600(true) == 1
@test t1_15600(false) == 2

# issue #15561
function f15561()
    a = @goto crater
    @label crater
end
@test f15561() === nothing

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
function OpenDialog(filterlist::Cstring, defaultpath::Cstring)
    outpath = Ref(Ptr{UInt8}())
    status = @ccall libnfd.NFD_OpenDialog(filterlist::Cstring, defaultpath::Cstring, outpath::Ref{Ptr{UInt8}})::NFDResult
    return status, outpath[]
end

function OpenDialogMultiple(filterlist::Cstring, defaultpath::Cstring)
    outpathset = Ref(NFDPathSet(C_NULL, C_NULL, 0))
    status = @ccall libnfd.NFD_OpenDialogMultiple(filterlist::Cstring, defaultpath::Cstring, outpathset::Ref{NFDPathSet})::NFDResult
    return status, outpathset[]
end

function SaveDialog(filterlist::Cstring, defaultpath::Cstring)
    outpath = Ref(Ptr{UInt8}())
    status = @ccall libnfd.NFD_SaveDialog(filterlist::Cstring, defaultpath::Cstring, outpath::Ref{Ptr{UInt8}})::NFDResult
    return status, outpath[]
end

function PickFolder(defaultpath::Cstring)
    outpath = Ref(Ptr{UInt8}())
    status = @ccall libnfd.NFD_PickFolder(defaultpath::Cstring, outpath::Ref{Ptr{UInt8}})::NFDResult
    return status, outpath[]
end

function GetError()
    @ccall libnfd.NFD_GetError()::Cstring
end

function PathSetFree(pathset::NFDPathSet)
    refpathset = Ref(pathset)
    @ccall libnfd.NFD_PathSet_Free(refpathset::Ref{NFDPathSet})::Cvoid
end

function Free(ptr)
    @ccall libnfd.NFD_Free(ptr::Ptr{Cvoid})::Cvoid
end


text:
module PwDynamicDefinition
using ValidatedNumerics
using ..DynamicDefinition
using ..Contractors
using TaylorSeries: Taylor1

using ..DynamicDefinition: derivative

export PwMap, preim, nbranches, plottable

"""
Dynamic based on a piecewise defined map.

The map is defined as T(x) = Ts[k](x) if x ∈ [endpoints[k], endpoints[k+1]).

We set is_full[k]==true if we can prove/assume that the kth branch is full, i.e., Ts[k](endpoints[k])==0 and Ts[k](endpoints[k+1])==1 or vice versa.
"""
struct PwMap <: Dynamic
	Ts::Array{Function, 1}
	endpoints::Array{Interval, 1}
	is_full
	orientations # these will be filled in automatically, usually
end

import ..DynamicDefinition: domain

DynamicDefinition.domain(S::PwMap) = hull(S.endpoints[1], S.endpoints[end])

PwMap(Ts, endpoints) = PwMap(Ts, endpoints, fill(false, length(endpoints)-1))
PwMap(Ts, endpoints, is_full) = PwMap(Ts, map(Interval, endpoints), is_full, [sign(derivative(Ts[i], (endpoints[i]+endpoints[i+1])/2)) for i in 1:length(endpoints)-1])

DynamicDefinition.nbranches(D::PwMap)=length(D.endpoints)-1

DynamicDefinition.is_full_branch(D::PwMap) = all(D.is_full)

function DynamicDefinition.preim(D::PwMap, k, y, ϵ)
	@assert 1 <= k <= nbranches(D)
	domain = hull(D.endpoints[k], D.endpoints[k+1])
	root(x->D.Ts[k](x)-y, domain, ϵ)
end

function DynamicDefinition.plottable(D::PwMap, x)
	for k in nbranches(D)
		domain = hull(D.endpoints[k], D.endpoints[k+1])
		if x in domain
			return D.Ts[k](x)
		end
	end
end

"""
Computes the hull of an iterable of intervals
"""
common_hull(S) = interval(minimum(x.lo for x in S), maximum(x.hi for x in S))

function DynamicDefinition.derivative(n, D::PwMap, x)
	common_hull(derivative(n, D.Ts[i], x ∩ hull(D.endpoints[i],D.endpoints[i+1])) for i in 1:length(D.endpoints)-1)
end

function DynamicDefinition.distorsion(D::PwMap, x)
	common_hull(distorsion(D.Ts[i], x ∩ hull(D.endpoints[i],D.endpoints[i+1])) for i in 1:length(D.endpoints)-1)
end

end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
using Matter
using Matter: Vector

@testset "vector arithmetic" begin
    a = Vector(0., 1.)
    b = Vector(-1., 2.)
    @test a + b == Vector(-1., 3.)
    @test a - b == Vector(1., -1.)
    @test a * 2 == Vector(0., 2.)
    @test a / 2 == Vector(0., 1. / 2.)
    @test a == Vector(0., 1.)
end

@testset "in-place arithmetic" begin
    a = Vector(0., 1.)
    b = Vector(-1., 2.)
    @test add!(a, b) == Vector(-1., 3.)
    a = Vector(0., 1.)
    @test sub!(a, b) == Vector(1., -1.)
    a = Vector(0., 1.)
    @test mult!(a, 2) == Vector(0., 2.)
    a = Vector(0., 1.)
    @test div!(a, 2) == Vector(0., 1. / 2.)
    a = Vector(0., 1.)
    @test a == Vector(0., 1.)
end

@testset "vector operations" begin
    a = Vector(0., 1.)
    b = cross(a)
    @test det(a, b) == 0.
    @test magnitude(a) == 1.
    @test normalize!(a) == Vector(0., 1.)
    a = Vector(1., 0.)
    rotate!(a, π / 2)
    @test isapprox(a, Vector(0., 1.))
    a = Vector(2., 0.)
    rotate!(a, π / 2, Vector(1., 0.))
    @test isapprox(a, Vector(1., 1.))
end


text:
import Base.==
"""
This type represents a single syllable of a sound.

A single syllable is a set of points in a sonogram. These points are stored in two arrays of equal length, `times` and `harmonics`. Every index in the sonogram of the form (times[i], harmonics[i]) ∀ i within the length of these arrays is part of the syllable, and any index not of this form is not part of the syllable.

Remember: Syllable and SyllableBox types just say *where* the syllable is in the sonogram where they were found. They contain no values from said sonogram (See [syllable and syllable box types](@ref syllSyllBoxTypes) for details)

Fields
------
 - `harmonics`: harmonic (down-column) coordinates of the syllable
 - `times`: time (along-row) coordinates of the syllable. *Note*: This is not in standard units of time. The unit of time represented depends on the window step used in building the sonogram
"""
immutable Syllable{T<:Integer}
    harmonics::Array{T, 1}
    times::Array{T, 1}
end

function ==(s1::Syllable, s2::Syllable)
    return (s1.harmonics == s2.harmonics) && (s1.times == s2.times)
end

"""
This type represents a rectangle around a syllable that minimally contains it
Remember: Syllable and SyllableBox types just say *where* the syllable is in the sonogram where they were found. They contain no values from said sonogram. (See [syllable and syllable box types](@ref syllSyllBoxTypes) for details)

Let h be index down columns (harmonic) and t be index along rows (time slice).
Fields
------
 - `hmin`: Lowest harmonic in the syllable.
 - `tmin`: Start time of the syllable
 - `hmax`: Highest harmonic in the syllable.
 - `tmax`: End time of the syllable.
"""
immutable SyllableBox{T<:Integer}
    hmin::T
    tmin::T
    hmax::T
    tmax::T
end

function ==(sb1::SyllableBox, sb2::SyllableBox)
    return (sb1.tmin == sb2.tmin) && (sb1.hmin == sb2.hmin) && (sb1.tmax == sb2.tmax) && (sb1.hmax == sb2.hmax)
end

"""A helper function to turn a linear index into a 2-D Array into a row index (I.E., a position down in a column)"""
function index1(size1::Integer, index::Integer)
    @assert index > 0
    if (index%size1) != 0
	return (index%size1)
    else
	return size1
    end
end

"""A helper function to turn a linear index into a 2-D Array into a column index (I.E., a position along in a row)"""
function index2(size1::Integer, index::Integer)
    return div(index-1, size1)+1
end

"""Given dimension 1 of a 2-D array and two indecies, determine if index1 is adjacent to index2"""
function is_adjacent(size1::Integer, index1::Integer, index2::Integer)
    if abs(index1 - index2)  <=  1 || abs(abs(index1 - index2) - size1) <= 1
	return  true
    end

    return false
end

"""
 `find_syllables{T}(sonogram::AbstractArray{T, 2}, thresholdlow::T, thresholdhigh::T)`

Find syllables in a bird call, with an algorithm very similar to the one presented by [1]

Parameters
----------
 - `sonogram`: A sonogram of the bird call
 - `thresholdlow`: The threshold below which values in `sonogram` are judged to be noise
 - `thresholdhigh`: The threshold above which values in `sonogram` trigger the syllable search (see below)

Return value
------------
 - An `Array` of `Syllable`s in `sonogram`

Algorithm
---------
The algorithm is (roughly) as follows:
 1) Search through `sonogram` for values greater than `thresholdhigh`. For each such value `v`:
  1.1) If the value is already in a previously found syllable, skip it and go on to the next
  1.2) Note every index `i` in `sonogram` greater than `thresholdlow` and adjacent to the index of `v`
  1.3) Note every index `i_2` in `sonogram greater than `thresholdlow` and adjacent to `i`
  1.4) And so on, until no more adjacent values greater than `thresholdlow` are found. This clump of adjacent values constitutes one _syllable_
 Note: Steps 1.1-1.4 are what are referred to above as the "syllable search"
 Further note: The above describes what is, in effect done, but not how. How is below in the code.
"""
function find_uncombined_syllables{T}(sonogram::AbstractArray{T, 2}, thresholdlow::Real, thresholdhigh::Real)
    syllables = Array{Syllable{Int}, 1}()
    const size1 = size(sonogram, 1)

    indecies = find(x->x>thresholdlow, sonogram)

    while !isempty(indecies)
	seed_index = indecies[1]
	adjacent_indecies = find(x->is_adjacent(size1, x, seed_index), indecies)
	to_add_to_adjacent_indecies = adjacent_indecies # Note: this doesn't actually get added; It's just so that the loop below will work properly the first time around

	done_adding_indecies = false
	while !done_adding_indecies
	    nadj_inds = length(adjacent_indecies)
	    last_added = to_add_to_adjacent_indecies
	    to_add_to_adjacent_indecies = Int[]
	    for index in view(indecies, last_added)
		for (j, jindex) in enumerate(indecies)
		    if is_adjacent(size1, jindex, index) && !(j in adjacent_indecies)
			push!(adjacent_indecies, j)
			push!(to_add_to_adjacent_indecies, j)
		    end
		end
	    end
	    if length(adjacent_indecies) == nadj_inds
		done_adding_indecies = true
	    end
	end
	sort!(adjacent_indecies)

	syllables_temp = indecies[adjacent_indecies]
	syllables_temp1 = map(i->index1(size1, i), syllables_temp)
	syllables_temp2 = map(i->index2(size1, i), syllables_temp)
	@assert !any(x->x<=0, syllables_temp1)
	@assert !any(x->x>size1, syllables_temp1)
	@assert !any(x->x<=0, syllables_temp2)
	@assert !any(x->x>size(sonogram, 2), syllables_temp2)

	if any(x->x>thresholdhigh, view(sonogram, syllables_temp))
	   push!(syllables, Syllable(syllables_temp1, syllables_temp2))
	end

	deleteat!(indecies, adjacent_indecies)

    end

    return syllables
end

"""
Take a set of syllables (problably output from `find_uncombined_syllables`) and combine syllables seperated by small time-gaps

Parameters
----------
 - `Syllables`: A set of syllables as an array of arrays, each array containing a set of indecies that is one syllable
 - `maxTimeGap`: The maximum time gap (in seconds) allowed between syllables before they are considered really and truly seperate
 - `conversion_factor`: The multiplier to convert maxTimeGap from units of seconds into units of indecies (which is what syllables are in).
"""
function combine_syllables{T}(syllables_in::AbstractArray{Syllable{T}, 1}, maxTimeGap::Real, conversion_factor::Real=samplerate/windowstep)
    maxTimeGap *= conversion_factor
    syllables_in = deepcopy(syllables_in)
    syllables = Syllable{Int}[]
    nsyllables = length(syllables_in)

    while length(syllables_in) > 0
        syllable_temp = syllables_in[1] #[syllables_in[1][1], syllables_in[1][2]]
	todelete = Int[1]
	for i in 2:length(syllables_in)
	    max_col_syll1 = maximum(syllables_in[i-1].times)
	    min_col_syll2 = minimum(syllables_in[i].times)

	    if min_col_syll2 <= (max_col_syll1 + maxTimeGap + 1)
		append!(syllable_temp.times, syllables_in[i].times)
		append!(syllable_temp.harmonics, syllables_in[i].harmonics)
		push!(todelete, i)
	    else
		break
	    end
	end

	push!(syllables, syllable_temp)

	deleteat!(syllables_in, todelete)
    end

    return syllables
end

"""
`function syllable_box(syllable{Tuple{Array{Int, 1}, Array{Int, 1}}})`
Returns the upper left and lower right corners of a box containing the syllable `syllable`
"""
function syllable_box(syllable::Syllable)::SyllableBox
    @assert !any(x->x==0, syllable.times)
    @assert !any(x->x==0, syllable.harmonics)
    return SyllableBox(minimum(syllable.harmonics), minimum(syllable.times), maximum(syllable.harmonics), maximum(syllable.times))
end

"""
`function syllable_boxes(syllables::AbstractArray{Tuple{Array{Int, 1}, Array{Int, 1}}, 1})`
Returns the upper left and lower right corners of a box containing each syllable in `syllables`
"""
function syllable_boxes{T}(syllables::AbstractArray{Syllable{T}, 1})
    boxes = map(syllable_box, syllables)

    return boxes
end

export Syllable, SyllableBox
export find_uncombined_syllables
export combine_syllables
export syllable_box
export syllable_boxes


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

using Markdown
using InteractiveUtils

# ╔═╡ 53d9d2c0-0369-403e-8378-858663608a35
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Latexify"),
        Pkg.PackageSpec(name="LaTeXStrings"),
        Pkg.PackageSpec(name="Plots"),
        Pkg.PackageSpec(name="ControlSystems"),
		Pkg.PackageSpec(name="ColorSchemes"),
		Pkg.PackageSpec(name="PlotThemes")
    ])
    using Latexify
    using LaTeXStrings
	using Plots
	using ControlSystems
	using ColorSchemes
	using PlotThemes
end

# ╔═╡ 451411e9-6eb1-4ad1-97cd-946300511c0d
begin
	Pkg.add(url="https://github.com/Pocket-titan/DarkMode")
    import DarkMode
end

# ╔═╡ 550a6097-7c37-41dc-aa8c-0da7b523f31a
begin
	md"# Teoría de Circuitos II - Filtro Pasa Todo

Tomamos el siguiente circuito, con el objetivo de analizar y simular su transferencia."
end

# ╔═╡ 67c80a93-838d-4f72-b917-75ff417e1d71
md"Si se usa el modelo OpAmp ideal $:\, V^+=V^-=V_x, \quad Z_{in}\rightarrow\infty$ se encuentra la siguiente ecuación matricial aplicando ley de nodos a $V^+$y $V^-$:"

# ╔═╡ 70bb620a-b3d3-4d48-a591-f2e3abbb5f59
begin
	ex1 = :([1/R_1+1/R_2 -1/R_2; 1/R_1+jwC 0])
	ex2 = :([V_x; V_2])
	ex3 = :([1/R_1; jwC])
	latexstring(latexraw(ex1),latexraw(ex2),"=",latexraw(ex3),"V_1")
end

# ╔═╡ b9353662-c3ca-460d-8724-ff01d8692b32
md"La transferencia $H(jw) = \frac{V_2}{V_1}$ se puede obtener despejando la ecuación y resolviendo para $V_2$."

# ╔═╡ ec139584-2a0a-4583-8392-ee5bcd7096f7
begin
	ex4 = :([1/R_1+1/R_2 -1/R_2; 1/R_1+jwC 0]^(-1))
	latexstring(latexraw(ex2),"=",latexraw(ex4),latexraw(ex3),"V_1")
end

# ╔═╡ 14dfa10b-ab0d-4b6b-b3b9-aad4838c157b
md"$\vdots$ Contemplando $R_1 = R_2$:"

# ╔═╡ 8eaa07ab-f6f1-45aa-a929-f25958bfef73
begin
	ex5 = :(H(jw)=(jCw-1/(R_3))/(jCw+1/(R_3)))
	latexify(ex5)
end

# ╔═╡ 67d13c95-9f5d-4753-90ed-b15612b89945
md"La simulación numérica la vamos a realizar en el lenguaje Julia, usando el módulo ControlSystems.jl"

# ╔═╡ 97467fe3-f873-40a0-a402-f4c2d398b956
begin
	    R1 = 1e3;
	    R2 = 1e3;
	    C  = 1e-6;
	    R3 = 1e3;
	md"Una vez definidas las variables R1, R2, C, R3, definimos los coeficientes del numerador y denominador."
end

# ╔═╡ 3f47cbfb-ccb5-4f19-99fb-387ccdc26554
num = [C/R2, -1/(R1*R3)]

# ╔═╡ 536d941c-24e3-454d-b8af-1da615cee757
den = [C/R2,  1/(R2*R3)]

# ╔═╡ 578b5413-6af8-429b-8e6f-4581704a21e6
md"ControlSystems.tf retorna H como un objeto TransferFunction:"

# ╔═╡ 5d68641a-ac0b-48c8-b7aa-9f4798f1f729
H = tf(num, den)

# ╔═╡ 0b356ae9-a7a0-4e19-add1-1d7498c10869
md"ControlSystems.bodeplot toma el obtjeto TransferFunction y obtiene la respuesta en frecuencia,que efectivamente corresponde a un filtro pasa todo."

# ╔═╡ 006447ea-7984-4eb2-916d-e85181b639c3
begin
	theme(:juno)
	bodeplot(H, palette=:seaborn_muted6, thm=:juno)
end

# ╔═╡ 77224aa9-e2cc-493c-a810-4affe4b33428
md"ControlSystems.rlocusplot toma el objeto TransferFunction y obtiene el diagrama de polos y ceros, que evidencia la estabilidad del sistema y hace visible la condición de filtro pasa todo: simetría de polos y ceros respecto al eje jw."

# ╔═╡ 2162cb78-efdf-4b65-bed0-e87b2b28d829
begin
	theme(:juno)
	rlocusplot(H; aspect_ratio=:equal, framestyle=:origin, palette=:seaborn_muted6)
end

# ╔═╡ cd8440fe-d58d-4c2f-a233-2f1aaec2c647
md"Finalmente, volviendo al modelo en LTSpice comparamos los resultados de la simulación numérica con una simulación circuital"

# ╔═╡ c5d02935-1c7c-47bd-8b52-d0f238c8cd9b
md"La simulación de LTSpice, a diferencia del modelo ideal, implementa el modelo integrador del amplificador operacional. El resultado es el siguiente:"

# ╔═╡ 26701700-7abf-40a6-a11f-776684f48ba7
md"Es un comportamiento equivalente a la simulación numérica a bajas frecuencias, filtro pasa todo, pero el comportamiento difiere a altas frecuencias donde el efecto del capacitor del modelo integrador es apreciable, y la transferencia real disminuye."

# ╔═╡ b69f9cbc-b116-4206-9ede-2a3a2d6bc135
begin
	struct Wow
		filename
	end

	function Base.show(io::IO, ::MIME"image/png", w::Wow)
		write(io, read(w.filename))
	end
end

# ╔═╡ 2d83491c-b238-40d2-903a-e6368089bdd3
Wow("Images/circuito.png")

# ╔═╡ ecd5ba2f-0337-4e8d-bdb0-4e8e2c15607b
Wow("Images/circuito.png")

# ╔═╡ 0193b0ca-26ff-45ce-bdaa-650a5fc807dd
Wow("Images/sim.png")

# ╔═╡ f99891ef-ab72-4bcb-8e14-10a70cdc5750
DarkMode.enable(theme="nord")

# ╔═╡ Cell order:
# ╟─550a6097-7c37-41dc-aa8c-0da7b523f31a
# ╟─2d83491c-b238-40d2-903a-e6368089bdd3
# ╟─67c80a93-838d-4f72-b917-75ff417e1d71
# ╟─70bb620a-b3d3-4d48-a591-f2e3abbb5f59
# ╟─b9353662-c3ca-460d-8724-ff01d8692b32
# ╟─ec139584-2a0a-4583-8392-ee5bcd7096f7
# ╟─14dfa10b-ab0d-4b6b-b3b9-aad4838c157b
# ╟─8eaa07ab-f6f1-45aa-a929-f25958bfef73
# ╟─67d13c95-9f5d-4753-90ed-b15612b89945
# ╟─97467fe3-f873-40a0-a402-f4c2d398b956
# ╠═3f47cbfb-ccb5-4f19-99fb-387ccdc26554
# ╠═536d941c-24e3-454d-b8af-1da615cee757
# ╟─578b5413-6af8-429b-8e6f-4581704a21e6
# ╠═5d68641a-ac0b-48c8-b7aa-9f4798f1f729
# ╟─0b356ae9-a7a0-4e19-add1-1d7498c10869
# ╠═006447ea-7984-4eb2-916d-e85181b639c3
# ╟─77224aa9-e2cc-493c-a810-4affe4b33428
# ╠═2162cb78-efdf-4b65-bed0-e87b2b28d829
# ╟─cd8440fe-d58d-4c2f-a233-2f1aaec2c647
# ╟─ecd5ba2f-0337-4e8d-bdb0-4e8e2c15607b
# ╟─c5d02935-1c7c-47bd-8b52-d0f238c8cd9b
# ╟─0193b0ca-26ff-45ce-bdaa-650a5fc807dd
# ╟─26701700-7abf-40a6-a11f-776684f48ba7
# ╟─b69f9cbc-b116-4206-9ede-2a3a2d6bc135
# ╟─53d9d2c0-0369-403e-8378-858663608a35
# ╟─451411e9-6eb1-4ad1-97cd-946300511c0d
# ╟─f99891ef-ab72-4bcb-8e14-10a70cdc5750

max_stars_repo_path: 1 - Filtro Pasa Todo/notebook.jl

max_stars_repo_name: Matrsz/Teoria-Circuitos-II-UTN

max_issues_repo_path: 1 - Filtro Pasa Todo/notebook.jl

max_issues_repo_name: Matrsz/Teoria-Circuitos-II-UTN

max_forks_repo_path: 1 - Filtro Pasa Todo/notebook.jl

max_forks_repo_name: Matrsz/Teoria-Circuitos-II-UTN


text:
using LinearAlgebra
using SparseArrays
using NetworkLearning
using LightGraphs: Graph, add_edge!
using SimpleWeightedGraphs: SimpleWeightedGraph
using LearnBase: ObsDim, ObsDimension
using Statistics: mean
using Test
using Future: copy!
using DelimitedFiles: readdlm

# Test components
include("t_components.jl")
Test.@testset "Network Learning (various components)" begin
	t_components();
end

# Test observation-based learning
include("t_observation_networklearner.jl")
Test.@testset "Network Learning (observation-based)" begin
	t_observation_networklearner();
end

# Test entity-based learning
include("t_entity_networklearner.jl")
Test.@testset "Network Learning (entity-based)" begin
	t_entity_networklearner();
end


text:
# Test the multivariate distribution support for HMC sampler

using Turing, Distributions

# Define helper functions
function sigmoid(t)
  return 1 / (1 + exp.(-t))
end

# Define NN flow
function nn(x, b1, w11, w12, w13, bo, wo)
  h = tanh.([w11' * x + b1[1]; w12' * x + b1[2]; w13' * x + b1[3]])
  return sigmoid((wo' * h)[1] + bo)
end

# Generating training data
N = 20
M = round(Int64, N / 4)
x1s = rand(M) * 5
x2s = rand(M) * 5
xt1s = Array([[x1s[i]; x2s[i]] for i = 1:M])
append!(xt1s, Array([[x1s[i] - 6; x2s[i] - 6] for i = 1:M]))
xt0s = Array([[x1s[i]; x2s[i] - 6] for i = 1:M])
append!(xt0s, Array([[x1s[i] - 6; x2s[i]] for i = 1:M]))

xs = [xt1s; xt0s]
ts = [ones(M); ones(M); zeros(M); zeros(M)]

# Define model

alpha = 0.16            # regularizatin term
var = sqrt(1.0 / alpha) # variance of the Gaussian prior

@model bnn(ts) = begin
  b1 ~ MvNormal([0 ;0; 0], [var 0 0; 0 var 0; 0 0 var])
  w11 ~ MvNormal([0; 0], [var 0; 0 var])
  w12 ~ MvNormal([0; 0], [var 0; 0 var])
  w13 ~ MvNormal([0; 0], [var 0; 0 var])
  bo ~ Normal(0, var)

  wo ~ MvNormal([0; 0; 0], [var 0 0; 0 var 0; 0 0 var])
  for i = rand(1:N, 10)
    y = nn(xs[i], b1, w11, w12, w13, bo, wo)
    ts[i] ~ Bernoulli(y)
  end
  b1, w11, w12, w13, bo, wo
end

# Sampling
chain = sample(bnn(ts), HMC(10, 0.1, 5))


text:
#=
	This script generates multiple simulated gibbs runs in parallel, all on the same simulated data.
	There are two ways to run this file, both which requires launching Julia in the top diretory.
	The first way is to run the following command from the terminal
	```
	$ julia -p 10 scripts/scriptParallel.jl
	```
	In the above command julia is started with 10 independent workers and simulations are made in parallel.
	You can change this to the number of cores/workers you want.

	The second way is the launch the script from within julia:
	```
	$ julia
	julia> addprocs(10)   # start 10 parallel workers
	julia> include("scripts/scriptParallel.jl")
	```
=#

#--- set the seed or generate a random one
# const seed = Base.Random.RANDOM_SEED
const seed = Uint32[2657077506, 531413685, 532641760, 1678774672] # big bump
srand(seed)

# ---- set directory name for the simulation
const scriptname = "scriptParallel"
const savepath = joinpath("simulations", "$(scriptname)_$(seed[1])") #<--change the directory name here

# ---- parameters of the simulation run
const maxiter = 2500  # sets the maximum number of gibbs iterations for each worker
const percentNyqForC = 0.5   # used for T l_max
const numofparsForP  = 1500  # used for P l_max
const hrfactor = 2.0
const pixel_size_arcmin = 2.0
const n = 2.0^9
const beamFWHM = 0.0
const nugget_at_each_pixel = (4.0)^2
begin
	local deltx =  pixel_size_arcmin * π / (180 * 60) # rads
	local period = deltx * n # side length in rads
	local deltk =  2 * π / period
	local nyq = (2 * π) / (2 * deltx)
	const maskupP  = √(deltk^2 * numofparsForP / π)  #l_max for for phi
	const maskupC  = min(9000.0, percentNyqForC * (2 * π) / (2 * pixel_size_arcmin * π / (180*60))) #l_max for for phi
	println("muK_per_arcmin = $(sqrt(nugget_at_each_pixel * (pixel_size_arcmin^2)))") # muK per arcmin
	println("maskupP = $maskupP") # muK per arcmin
	println("maskupC = $maskupC") # muK per arcmin
	println("$(seed[1])") # muK per arcmin
end
const scale_grad =  2.0e-3
const scale_hmc  =  2.0e-3

# ------------ load modules and functions
@everywhere push!(LOAD_PATH, pwd()*"/src")
@everywhere using Interp
@everywhere FFTW.set_num_threads(2)
require("cmb.jl")
require("fft.jl")
require("funcs.jl") # use reload after editing funcs.jl
require("Interp.jl")

# --------- generate cmb spectrum class for high res and low res
parlr = setpar(
	pixel_size_arcmin, n, beamFWHM, nugget_at_each_pixel,
	maskupC, maskupP
)
parhr = setpar(
	pixel_size_arcmin./hrfactor, hrfactor*n, beamFWHM, nugget_at_each_pixel,
	maskupC, maskupP
)

# -------- Simulate data: ytx, maskvarx, phix, tildetx
ytk_nomask, tildetk, phix, tx_hr = simulate_start(parlr);
phik = fft2(phix, parlr)
tmpdo = maximum(parlr.grd.x)*0.3
tmpup = maximum(parlr.grd.x)*0.4
maskboolx = tmpdo .<= parlr.grd.x .<= tmpup
maskvarx  = parlr.nugget_at_each_pixel .* ones(size(parlr.grd.x))
maskvarx[maskboolx] = Inf
ytx = ifft2r(ytk_nomask, parlr)
ytx[maskboolx] = 0.0
ytk = fft2(ytx, parlr)

# -----------  set the save directory and save stuff
isdir(savepath) && run(`rm -rf $savepath`)
run(`mkdir $savepath`)
for r in workers()
	run(`mkdir $savepath/job$r`)
end
run(`cp src/funcs.jl $savepath/funcs_save.jl`)
run(`cp scripts/$scriptname.jl $savepath/$(scriptname)_save.jl`)
writecsv("$savepath/seed.csv", seed)
writecsv("$savepath/x.csv", parlr.grd.x)
writecsv("$savepath/y.csv", parlr.grd.y)
writecsv("$savepath/k1.csv", parlr.grd.k1)
writecsv("$savepath/k2.csv", parlr.grd.k2)
writecsv("$savepath/ytx.csv", ytx)
writecsv("$savepath/phix.csv", phix)
writecsv("$savepath/tildetx_lr.csv", ifft2r(tildetk, parlr))
writecsv("$savepath/tx_lr.csv", tx_hr[1:3:end, 1:3:end] )
writecsv("$savepath/maskvarx.csv", maskvarx)
# print out the quadratic estimate for comparision
al_lr = ttk_al(parlr)
qex_lr = ifft2r(ttk_est(ytk_nomask, parlr) .* (parlr.cPP ./ (2.0 .* al_lr + parlr.cPP)) , parlr)
writecsv("$savepath/qex_lr.csv", qex_lr)

# ------------------ initalized and run the gibbs
@everywhere function gibbsloop(its, parhr, parlr, ytx, maskvarx, maskupC, savepath, r, scale_grad, scale_hmc, hrfactor)
	acceptclk   = [1] #initialize acceptance record
	tx_hr_curr  = zero(parhr.grd.x)
	ttx         = zero(parhr.grd.x)
	p1hr, p2hr  = zero(parhr.grd.x), zero(parhr.grd.x)
	phik_curr   = zero(fft2(ytx, parlr))
	tildetx_hr_curr = zero(parhr.grd.x)

	for bglp = 1:its
		tic()
		# ----- update tildetx_hr_curr
		if bglp % 100 == 1
			p1hr[:], p2hr[:] = gibbspass_t!(tx_hr_curr, ttx, phik_curr, ytx,
				maskvarx, parlr, parhr, repstretch(2.0, 1.2maskupC, 75, 1000)
			)
		end
		p1hr[:], p2hr[:] = gibbspass_t!(tx_hr_curr, ttx, phik_curr, ytx,
			maskvarx, parlr, parhr, fill(Inf, 400)
		)

		tildetx_hr_curr[:] = spline_interp2(
			parhr.grd.x, parhr.grd.y, tx_hr_curr,
			parhr.grd.x + p1hr, parhr.grd.y + p2hr
		)

		# ----- update phik_curr
		if bglp <= 10
			gradupdate!(phik_curr, tildetx_hr_curr, parlr, parhr, scale_grad);
		else
			push!(acceptclk, hmc!(phik_curr, tildetx_hr_curr, parlr, parhr, scale_hmc))
		end
		if bglp % 50 == 1
			writecsv("$savepath/job$r/tildetx_lr_curr_$bglp.csv", tildetx_hr_curr[1:int(hrfactor):end,1:int(hrfactor):end])
			writecsv("$savepath/job$r/tx_lr_curr_$bglp.csv", tx_hr_curr[1:int(hrfactor):end,1:int(hrfactor):end])
			writecsv("$savepath/job$r/phix_curr_$bglp.csv", ifft2r(phik_curr, parlr))
			writecsv("$savepath/job$r/acceptclk.csv", acceptclk)
		end
		toc()
	end # for
end # function

# ------- run parallel jobs
@sync @parallel for r in workers()
	gibbsloop(maxiter, parhr, parlr, ytx, maskvarx, maskupC, savepath, r, scale_grad, scale_hmc, hrfactor)
end

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
"""
    @declare_input mymap(rt, x::String, y::Int) :: String

Macro to create Salsa Input functions, which are the interface to store and retrieve
"input values" on a `Salsa.Runtime`.

This macro creates a function with the same name to retrieve input values from a Runtime,
as well as a setter function and deleter function to set values of the specified return
type on the Runtime, stored for a given key.

Salsa input functions can be called from derived functions, and will return the currently
stored value. Whenever the user sets/deletes a new value for the given input, it will
invalidate all dependent derived functions, so that they will be rerun and can access the
new value.

# Example
```julia-repl
julia> @declare_input student_grade(rt, name::String)::Float64
(student_grade, set_student_grade!, delete_student_grade!)

julia> rt = Runtime();

julia> set_student_grade!(rt, "Nathan", 3.5)

julia> student_grade(rt, "Nathan")
3.5

julia> delete_student_grade!(rt, "Nathan")

julia> student_grade(rt, "Nathan")
ERROR: KeyError: key DependencyKey(key=InputKey{typeof(student_grade),Tuple{String}}(), args=("Nathan",)) not found
```
"""
macro declare_input(e::Expr)
    # We require the inputs to have a type specified so that Salsa can store them in a
    # typed dictionary per-input.
    @assert e.head === :(::) "Missing type annotation. Expected: @declare_input i(...)::Type"
    @assert length(e.args) >= 1 && e.args[1] isa Expr && e.args[1].head === :call
    value_t = e.args[end]

    callexpr = e.args[1]
    inputname = callexpr.args[1]

    args = callexpr.args[2:end]

    # _argnames and _argtypes fill in anonymous names for unnamed args (`::Int`) and `Any`
    # for untyped args. `fullargs` will have all args w/ names and types.
    argnames = _argnames(args)
    argtypes = _argtypes(args)
    # NOTE: We don't support default values in inputs!
    # TODO: Do we want to support this? We'd probably have to / want to change the order of
    # the fields in `set_*!()`, so that the value comes first? Or maybe we could fake this
    # by generating multiple methods, but frankly that seems a bit weird.
    fullargs = [Expr(:(::), argnames[i], argtypes[i]) for i = 1:length(args)]
    # Put the filled-out args back into the original expression.
    args = callexpr.args[2:end] = fullargs

    # TODO: Use the returntype to strongly type the DefaultStorage dictionaries!
    # TODO: Use base's deduced return type (can be more specific) for DefaultStorage.

    runtime_arg = args[1]
    implicit_value_arg = gensym("value")

    # Eval all the argument types, to make sure they're valid.
    args_typetuple = Tuple(Core.eval(__module__, t) for t in argtypes)
    TT = Tuple{args_typetuple[2:end]...}  # Will be Tuple{} if not 2+ args.

    # Assert that the @declare_input starts with a runtime arg. The arg can be untyped, or
    # unnamed, it just has to be at least able to hold a Runtime.
    if !(
        length(args_typetuple) >= 1 &&
        # This allows any input like these: `_`, `rt`, `::Runtime`, `::MyRuntime`, `r::Any`
        # And disallows inputs like these: `::Int`, `x::String`.
        (args_typetuple[1] <: Runtime || args_typetuple[1] >: Runtime)
    )
        err_str = "@declare_input functions must accept a `Runtime` as the first argument."
        if length(args_typetuple) >= 1
            err_str *= " Got unexpected $(args_typetuple[1]) instead."
        end
        throw(ArgumentError(err_str))
    end

    # If we've made it here, there shouldn't be any reason the definitions below will fail,
    # so it's okay to pre-declare the getter function type, which we use in the InputKey,
    # below.
    getter_f = @eval __module__ function $inputname end

    # Build the Key type here, at macro parse time, since it's expensive to construct at runtime.
    # (Use type of function, not obj, because closures are not isbits)
    input_key_t = InputKey{typeof(getter_f),TT}
    dependency_key_expr = :($input_key_t(($(argnames[2:end]...),)))
    getter_body = quote
        $memoized_lookup_unwrapped($runtime_arg, $dependency_key_expr)::$value_t
    end

    setter_body = quote
        $Salsa._safe_setter_body($runtime_arg, $dependency_key_expr, $implicit_value_arg)
    end

    deleter_body = quote
        $Salsa._safe_deleter_body($runtime_arg, $dependency_key_expr)
    end

    # Construct the Getter function from the provided call expr.
    # (deepcopy the callexpr so that the subsequent calls don't modify it)
    getter = Expr(:function, deepcopy(callexpr), getter_body)

    # Construct the Deleter next, since it has the same interface as the Getter.
    deleter_name = callexpr.args[1] = Symbol("delete_$(inputname)!")
    # (deepcopy the callexpr so that the setter doesn't modify it)
    deleter = Expr(:function, deepcopy(callexpr), deleter_body)

    # Construct the Setter
    # Inject the value to be set & rename the function.
    setter_name = callexpr.args[1] = Symbol("set_$(inputname)!")
    push!(callexpr.args, :($implicit_value_arg::$value_t))
    setter = Expr(:function, callexpr, setter_body)

    esc(quote
        Core.@__doc__ $getter
        $setter
        $deleter

        # For type stability
        @inline function $Salsa.get_user_function(
            $(fullargs[1]),
            ::$input_key_t,
        )
            return $getter
        end

        # Return all the generated functions as a hint to REPL users what we're generating.
        ($inputname, $setter_name, $deleter_name)
    end)
end

function _safe_setter_body(runtime::Salsa._TopLevelRuntime, key, value)
    return set_input!(runtime, key, value)
end
function _safe_deleter_body(runtime::Salsa._TopLevelRuntime, key)
    return delete_input!(runtime, key)
end

function _safe_setter_body(::Salsa._TracingRuntime, key, value)
    error("Attempted impure operation in a derived function!: set_input!(rt, $key, $value)")
end
function _safe_deleter_body(::Salsa._TracingRuntime, key)
    error("Attempted impure operation in a derived function!: delete_input!(rt, $key)")
end


text:
#functions for calling the matrix vector products
#i.e. function A(y_in,Nx,Ny,h) computes the product A*u and stores it in y

#Right now only for  the second-order accurate SBP operators for constant coefficients from
#Mattsson and Nordstrom, 2004.

#Written for the 2D domain (x,y) \in (a, b) \times (c, d),
#stacking the grid function in the vertical direction,
#so that, e.g., \partial u / \partial x \approx D \kron I and
#\partial u / \partial y \approx I \kron D
#the faces of the 2D domain are label bottom to top, left to right, i.e.
#side 1 is y = c, side 2 is y = d
#side 3 is x = a, side 4 is x = b

#D2x and D2y compute approximations to
#\partial^2 u / \partial x^2 and
#\partial^2 u / \partial y^2, respectively

#FACEtoVOL(face,u_face,Nx,Ny) maps the face value u_face to a full length vector at the
#nodes corresponding to face

#BxSx \approx the traction on faces 3 and 4
#BySy \approx the traction on faces 1 and 2

#BxSx_tran and BySy_tran are the transposes of BxSx and BySy

# variables order: y_in, Nx, Ny, N, hx, hy, y_out(container)

##
# y_in: input Array
# y_out: output Array
##

function D2x_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1})
	#N = Nx*Ny
	#hx = Float64(1/(Nx-1))
	#hy = Float64(1/(Ny-1))
	#y = similar(y_in)
	#y = similar(y_in)
	@inbounds  for idx = 1:Ny
		y_out[idx] = (y_in[idx] - 2*y_in[Ny + idx] + y_in[2*Ny + idx]) / hx^2
	end

	@inbounds for idx1 = Ny+1:N-Ny
		y_out[idx1] = (y_in[idx1 - Ny] - 2 * y_in[idx1] + y_in[idx1 + Ny]) / hx^2
	end

	@inbounds for idx2 = N-Ny+1:N
		y_out[idx2] = (y_in[idx2 - 2*Ny] -2 * y_in[idx2 - Ny] + y_in[idx2]) / hx^2
	end
	# return y_out
end

function D2y_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1})
	#N = Nx*Ny
	#hx = Float64(1/(Nx-1))
	#hy = Float64(1/(Ny-1))
	@inbounds for idx = 1:Ny:N-Ny+1
		y_out[idx] = (y_in[idx] - 2 * y_in[idx + 1] + y_in[idx + 2]) / hy^2
	end

	@inbounds for idx1 = Ny:Ny:N
		y_out[idx1] = (y_in[idx1 - 2] - 2 * y_in[idx1 .- 1] + y_in[idx1]) / hy^2
	end

	@inbounds for j = 1:Nx
		@inbounds for idx = 2+(j-1)*Ny:j*Ny-1
			y_out[idx] = (y_in[idx - 1] - 2 * y_in[idx] + y_in[idx + 1]) / hy^2
		end
	end

	# return y_out

end

function Dx_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1})

	@inbounds for idx = 1:Ny
		y_out[idx] = (y_in[idx + Ny] - y_in[idx]) / hx
	end

	@inbounds for idx1 = Ny+1:N-Ny
		y_out[idx1] = (y_in[idx1 + Ny]-y_in[idx1 - Ny]) / (2*hx)
	end

	@inbounds for idx2 = N-Ny+1:N
		y_out[idx2] = (y_in[idx2]-y_in[idx2 .- Ny]) ./ hx
	end
	# return y_out
end

function Dy_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1})

	@inbounds for idx = 1:Ny:N-Ny+1
		y_out[idx] = (y_in[idx + 1] - y_in[idx]) / h
	end

	@inbounds for idx = Ny:Ny:N
		y_out[idx] = (y_in[idx] - y_in[idx - 1]) /h
	end

	@inbounds for j = 1:Nx
		@inbounds for idx = 2+(j-1)*Ny:j*Ny-1
		 	y_out[idx] = (y_in[idx + 1] - y_in[idx - 1]) / (2*h)
		end
	end
	# return y_out
end

function Hxinv_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, coef::Float64, y_out::Array{Float64,1})
	#N = Nx*Ny
	#y = similar(y_in)
	@inbounds for idx = 1:Ny
		y_out[idx] = coef * (2*y_in[idx]) * (1/hx)
	end

	@inbounds for idx1 = Ny+1:N-Ny
		y_out[idx1] = coef * (1*y_in[idx1]) * (1/hx)
	end

	@inbounds for idx2 = N-Ny+1:N
		y_out[idx2] = coef * (2*y_in[idx2]) * (1/hx)
	end

	# return y_out
end

function Hyinv_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, coef::Float64, y_out::Array{Float64,1})
	#N = Nx*Ny
	#y = similar(y_in)
	#hx = Float64(1/(Nx-1))
	#hy = Float64(1/(Ny-1))
	@inbounds for idx = 1:Ny:N-Ny+1
		y_out[idx] = coef * (2*y_in[idx]) * (1/hy)
	end

	@inbounds for idx1 = Ny:Ny:N
		y_out[idx1] = coef * (2*y_in[idx1]) * (1/hy)
	end

	@inbounds for i = 1:Nx
		@inbounds for idx2 = 2+(i-1)*Ny:i*Ny-1
			y_out[idx2] = coef * (y_in[idx2]) * (1/hy)
		end
	end
	# return y_out
end

function Hx_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64,coef::Float64 ,y_out::Array{Float64,1})
	@inbounds for idx = 1:Ny
		y_out[idx] = coef * hx*y_in[idx]/2
	end

	@inbounds for idx1 = Ny+1:N-Ny
		y_out[idx1] = coef* hx*y_in[idx1]
	end

	@inbounds for idx2 = N-Ny+1:N
		y_out[idx2] = coef * hx*y_in[idx2]/2
	end
	# return y_out
end

function Hy_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1})
	#N = Nx*Ny
	#y = similar(y_in)

	@inbounds for idx = 1:Ny:N-Ny+1
		y_out[idx] = hy*y_in[idx]/2
	end

	@inbounds for idx1 = Ny:Ny:N
		y_out[idx1] = hy*y_in[idx1]/2
	end

	@inbounds for i = 1:Nx
		@inbounds for idx2 = 2 + (i-1)*Ny:i*Ny-1
			y_out[idx2] = hy*y_in[idx2]
		end
	end
	# return y_out
end

function FACEtoVOL_beta!(y_in_face::Array{Float64,1},  face::Int64, Nx::Int64, Ny::Int64, N::Int64, y_outs::Array{Array{Float64,1},1})
	if face == 1
		idx = 1:Ny:N-Ny+1
	elseif face==2
		idx = Ny:Ny:N
	elseif face==3
		idx = 1:Ny
	elseif face == 4
		idx = N-Ny+1:N
	else
	end
	y_out = y_outs[face]
	y_out[idx] = y_in_face
	# return y_out
end

function VOLtoFACE_beta!(y_in::Array{Float64,1}, face::Int64,Nx::Int64, Ny::Int64, N::Int64, y_outs::Array{Array{Float64,1},1}) ## Has some issue
	if face == 1
			idx = 1:Ny:N-Ny+1
	elseif face == 2
			idx = Ny:Ny:N
	elseif face == 3
			idx = 1:Ny
	elseif face == 4
			idx = N-Ny+1:N
	else
	end
    y_out = y_outs[face]
	y_out[idx] = y_in[idx]

	# return y_out
end

function Bx_beta!(Nx::Int64,Ny::Int64, N::Int64, y_out::Array{Float64,1})

	@inbounds for idx=1:Ny
		y_out[idx] = -1
	end

	@inbounds for idx = N-Ny+1:N
		y_out[idx] = 1
	end
	# return y_out
end

function By_beta!(Nx::Int64,Ny::Int64,N::Int64,y_By::Array{Float64,1})
	@inbounds for idx = 1:Ny:N-Ny+1
		y_By[idx] = -1
	end

	@inbounds for idx1 = Ny:Ny:N
		y_By[idx1] = 1
	end
	# return y_By
end

function BxSx_beta!(y_in::Array{Float64,1},Nx::Int64,Ny::Int64,N::Int64,hx::Float64,hy::Float64,y_out::Array{Float64,1})
	@inbounds for idx = 1:Ny
		y_out[idx] = (1/hx) * (1.5 * y_in[idx] - 2 * y_in[idx + Ny] + 0.5 * y_in[idx + 2*Ny])
		y_out[N-Ny + idx] = (1/hx) * (0.5 * y_in[N-3*Ny + idx] - 2 * y_in[N-2*Ny + idx] + 1.5 * y_in[N-Ny + idx])
	end
	# return y_out
end

function BySy_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64,y_out::Array{Float64,1})
	@inbounds for idx = 1:Ny:N-Ny+1
		y_out[idx] = (1/hy) * (1.5 * y_in[idx] - 2 * y_in[idx .+ 1] + 0.5 * y_in[idx .+ 2])
	end

	@inbounds for idx = Ny:Ny:N
		y_out[idx] = (1/hy) * (0.5 * y_in[idx - 2] - 2 * y_in[idx - 1] + 1.5 * y_in[idx])
	end

	# return y_out
end

function BxSx_tran_beta!(y_in::Array{Float64,1},Nx::Int64,Ny::Int64,N::Int64,hx::Float64,hy::Float64,y_out::Array{Float64,1}) # be careful with += expression
	#hx = Float64(1/(Nx-1))
   	#hy = Float64(1/(Ny-1))
	@inbounds for idx1 = 1:Ny
		y_out[idx1] = (1.5 * y_in[idx1]) * (1/hx)

		# for idx = Ny+1:2*Ny
		y_out[idx1+Ny] = (-2 * y_in[idx1]) * (1/hx)
		# end

		# for idx  = 2*Ny+1:3*Ny
		y_out[idx1+2Ny] = (0.5 * y_in[idx1]) * (1/hx)
		# end
	end

	@inbounds for idxN = N-Ny+1:N
		y_out[idxN] = (1.5 * y_in[idxN]) * (1/hx)

		# for idx = N-2*Ny+1:N-Ny
		y_out[idxN-Ny] = (-2 * y_in[idxN]) * (1/hx)
		# end

		# for idx = N-3*Ny+1:N-2*Ny
		y_out[idxN-2Ny] = (0.5 * y_in[idxN]) * (1/hx)
		# end
	end
	# return y_out
end

function BySy_tran_beta!(y_in::Array{Float64,1}, Nx::Int64, Ny::Int64, N::Int64, hx::Float64, hy::Float64, y_out::Array{Float64,1}) # Be Careful about double foor loops
	@inbounds for idx1 = 1:Ny:N-Ny+1
		y_out[idx1] = (1.5 * y_in[idx1]) * (1/hy)

		# for idx = Ny+1:2*Ny
		y_out[idx1+1] = (-2 * y_in[idx1]) * (1/hy)
		# end

		# for idx  = 2*Ny+1:3*Ny
		y_out[idx1+2] = (0.5 * y_in[idx1]) * (1/hy)
		# end
	end

	@inbounds for idxN = Ny:Ny:N
		y_out[idxN] = (1.5 * y_in[idxN]) * (1/hy)

		# for idx = N-2*Ny+1:N-Ny
		y_out[idxN-1] = (-2 * y_in[idxN]) * (1/hy)
		# end

		# for idx = N-3*Ny+1:N-2*Ny
		y_out[idxN-2] = (0.5 * y_in[idxN]) * (1/hy)
		# end
	end

	# return y_out
end

# Previous function from deriv_ops.jl

# function FACEtoVOL(y_in_face, face, Nx, Ny)
# 	N = Nx*Ny
# 	y = zeros(N)

# 	if face == 1
# 		idx = 1:Ny:N-Ny+1
# 	elseif face == 2
# 		idx = Ny:Ny:N
# 	elseif face == 3
# 		idx = 1:Ny
# 	elseif face == 4
# 		idx = N-Ny+1:N
# 	else
# 	end

# 	y[idx] = u_face

# 	# return y

# end

# function VOLtoFACE(y_in, face, Nx, Ny)
# 	N = Nx*Ny
#         y = zeros(N)

#         if face == 1
#                 idx = 1:Ny:N-Ny+1
#         elseif face == 2
#                 idx = Ny:Ny:N
#         elseif face == 3
#                 idx = 1:Ny
#         elseif face == 4
#                 idx = N-Ny+1:N
#         else
#         end

# 	y[idx] = y_in[idx]
#         # return y
# end

# function Hxinv(y_in, Nx, Ny, h)
# 	N = Nx*Ny
# 	y = zeros(N)

# 	idx = 1:Ny
# 	y[idx] = (2*y_in[idx]) .* (1/h)

# 	idx = Ny+1:N-Ny
# 	y[idx] = (1*y_in[idx]) .* (1/h)

# 	idx = N-Ny+1:N
# 	y[idx] = (2*y_in[idx]) .* (1/h)

# 	# return y
# end

# function Hyinv(y_in, Nx, Ny, h)
# 	N = Nx*Ny
# 	y = zeros(N)

# 	idx = 1:Ny:N-Ny+1
# 	y[idx] = (2*y_in[idx]) .* (1/h)

# 	idx = Ny:Ny:N
# 	y[idx] = (2*y_in[idx]) .* (1/h)

# 	for i = 1:Nx
# 		idx = 2+(i-1).*Ny:i*Ny-1
# 		y[idx] = y_in[idx] .* (1/h)
# 	end

# 	# return y

# end

# function Hx(y_in, Nx, Ny, h)
# 	N = Nx*Ny
#         y = zeros(N)

#         idx = 1:Ny
# 	y[idx] = h .* (1/2)*y_in[idx]

#         idx = Ny+1:N-Ny
#         y[idx] = h .* 1*y_in[idx]

#         idx = N-Ny+1:N
# 	y[idx] = h .* (1/2)*y_in[idx]

#         # return y

# end

# function Hy(y_in, Nx, Ny, h)
# 	N = Nx*Ny
#         y = zeros(N)

#         idx = 1:Ny:N-Ny+1
# 	y[idx] = h .* (1/2)*y_in[idx]

#         idx = Ny:Ny:N
# 	y[idx] = h .* (1/2)*y_in[idx]

#         for i = 1:Nx
#                 idx = 2+(i-1).*Ny:i*Ny-1
#                 y[idx] = h .* y_in[idx]
#         end

#         # return y

# end

# function BxSx_tran(y_in, Nx, Ny, h)
# 	N = Nx*Ny
# 	y = zeros(N)

# 	idx1 = 1:Ny
# 	y[idx1] += (1.5 .* y_in[idx1]) .* (1/h)
# 	idx = Ny+1:2*Ny
# 	y[idx] += (-2 .* y_in[idx1]) .* (1/h)
# 	idx  = 2*Ny+1:3*Ny
# 	y[idx] += (0.5 .* y_in[idx1]) .* (1/h)

# 	idxN = N-Ny+1:N
# 	y[idxN] += (1.5 .* y_in[idxN]) .* (1/h)
# 	idx = N-2*Ny+1:N-Ny
# 	y[idx] += (-2 .* y_in[idxN]) .* (1/h)
# 	idx = N-3*Ny+1:N-2*Ny
# 	y[idx] += (0.5 .* y_in[idxN]) .* (1/h)

# 	# return y
# end


text:
export DynamicDistributionModel

"""
    DynamicDistributionModel(f::Tf, ns::Int, na::Int) -> DynamicDistributionModel{Tf}

Use a general function `f` to store the transformations. `ns` and `na` are the number of states and actions.
"""
struct DynamicDistributionModel{Tf<:Function} <: AbstractEnvironmentModel
    f::Tf
    ns::Int
    na::Int
end

RLBase.get_observation_space(m::DynamicDistributionModel) = DiscreteSpace(m.ns)
RLBase.get_action_space(m::DynamicDistributionModel) = DiscreteSpace(m.na)

(m::DynamicDistributionModel)(s, a) = m.f(s, a)


text:
function checkfortype(dataobject::InfoType, datatype::Symbol)
    if !dataobject.hydro && datatype==:hydro
        error("[Mera]: Simulation has no hydro files!")
    elseif !dataobject.amr && datatype==:amr
        error("[Mera]: Simulation has no amr files!")
    elseif !dataobject.gravity && datatype==:gravity
        error("[Mera]: Simulation has no gravity files!")
    elseif !dataobject.rt && datatype==:rt
        error("[Mera]: Simulation has no rt files!")
    elseif !dataobject.particles && datatype==:particles
        error("[Mera]: Simulation has no particle files!")
    elseif !dataobject.clumps && datatype==:clumps
        error("[Mera]: Simulation has no clump files!")
    elseif !dataobject.sinks && datatype==:sinks
        error("[Mera]: Simulation has no sink files!")
    end
end

function checklevelmax(dataobject::InfoType, lmax::Real)
    if dataobject.levelmax < lmax
        error("[Mera]: Simulation lmax=$(dataobject.levelmax) < your lmax=$lmax")
    elseif lmax < dataobject.levelmin
        error("[Mera]: Simulation lmin=$(dataobject.levelmin) > your lmin=$lmax")
    end
end

# use lmax in case user forces to load a uniform grid from amr data (lmax=levelmin)
function checkuniformgrid(dataobject::InfoType, lmax::Real)
    isamr = true
    if lmax == dataobject.levelmin
        isamr = false
    end
    return isamr
end

function checkuniformgrid(dataobject::DataSetType, lmax::Real)
    isamr = true
    if lmax == dataobject.info.levelmin
        isamr = false
    end
    return isamr
end