text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ServiceDelegationPropertiesFormat2 <: SwaggerModel
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    actions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: actions
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ServiceDelegationPropertiesFormat2(;serviceName=nothing, actions=nothing, provisioningState=nothing)
        o = new()
        validate_property(ServiceDelegationPropertiesFormat2, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(ServiceDelegationPropertiesFormat2, Symbol("actions"), actions)
        setfield!(o, Symbol("actions"), actions)
        validate_property(ServiceDelegationPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ServiceDelegationPropertiesFormat2

const _property_map_ServiceDelegationPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("serviceName")=>Symbol("serviceName"), Symbol("actions")=>Symbol("actions"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ServiceDelegationPropertiesFormat2 = Dict{Symbol,String}(Symbol("serviceName")=>"String", Symbol("actions")=>"Vector{String}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ServiceDelegationPropertiesFormat2 }) = collect(keys(_property_map_ServiceDelegationPropertiesFormat2))
Swagger.property_type(::Type{ ServiceDelegationPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceDelegationPropertiesFormat2[name]))}
Swagger.field_name(::Type{ ServiceDelegationPropertiesFormat2 }, property_name::Symbol) =  _property_map_ServiceDelegationPropertiesFormat2[property_name]

function check_required(o::ServiceDelegationPropertiesFormat2)
    true
end

function validate_property(::Type{ ServiceDelegationPropertiesFormat2 }, name::Symbol, val)
end


text:
using Nernst
using Test

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

ext: jl


text:
function get_endpoints(N)
	inners = 2 .^ collect(0:N-1)
	outers = reverse(inners)
	combs = Float64[]
	for (inner, outer) in zip(inners, outers)
		endpoints = [1, 2]
		c = repeat(endpoints, inner=inner, outer=outer)
		if Base.isempty(combs)
			combs = c
		else
			combs = hcat(combs, c)
		end
	end
	return combs
end

get_sequential_endpoints(N, i) = digits(i-1, base=2, pad=N) .+ 1

function d_dsw(X⃗::Vector{Interval}, Y⃗::Vector{Interval}; squared::Bool=false)
	p = length(X⃗)
	num_endpoints = 2 * p

	d_min = Inf
	d_max = -Inf
	for idx_endpoint = 1:2^num_endpoints
		endpoints = get_sequential_endpoints(num_endpoints, idx_endpoint)
		dⱼᵢ = 0
		for i = 1:p
			dⱼᵢ += (X⃗[i][endpoints[i]] - Y⃗[i][endpoints[p + i]])^2
		end
		if (!squared) dⱼᵢ ^= 0.5 end

		d_min = min(d_min, dⱼᵢ)
		d_max = max(d_max, dⱼᵢ)
	end
	Interval(d_min, d_max)
end

function d_dsw(X⃗::FuzzyVector, Y⃗::FuzzyVector; squared::Bool=false)
	levels = X⃗[1].levels
	num_levels = length(levels)

	if X⃗ == Y⃗
		d = SingletonFuzzyNumber(levels, number=0)
	elseif length(X⃗) ≠ length(Y⃗)
		return false
	else
		grades = Vector{Interval}(undef, num_levels)
		for (lvl, α) in enumerate(levels)
			X⃗_cut = cut(X⃗, α)
			Y⃗_cut = cut(Y⃗, α)
			grades[lvl] = d_dsw(X⃗_cut, Y⃗_cut; squared=squared)
		end
		d = FuzzyNumber(levels, grades)
	end
	d
end

function u(X⃗::FuzzyVector, prototypes::Vector{FuzzyVector}; m::Real=2.0, method::String="lfcm")
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	levels = X⃗[1].levels
	num_levels = length(levels)
	c = length(prototypes)

	Iⱼ = [i for i = 1:c if X⃗ == prototypes[i]]
	uⱼ = Vector{FuzzyNumber}(undef, c)
	if isempty(Iⱼ) # Iⱼ = ∅
		for i = 1:c
			grades = Vector{Interval}(undef, num_levels)
			for (lvl, α) in enumerate(levels)
				X_cut = cut(X⃗, α)
				C_cut = cut(prototypes,  α)
				if method == "lfcm"
					grades[lvl] = u_lfcm(X_cut, C_cut, i, m=m)
				elseif method == "dsw"
					grades[lvl] = u_dsw(X_cut, C_cut, i, m=m)
				end
			end
			uⱼ[i] = FuzzyNumber(levels, grades)
		end
	else # Iⱼ ≠ ∅
		for i = 1:c
			if i ∉ Iⱼ
				uⱼ[i] = SingletonFuzzyNumber(levels, number=0)
			else
				uⱼ[i] = SingletonFuzzyNumber(levels, number=1 / length(Iⱼ))
			end
		end
	end
	uⱼ
end

function u_dsw(X⃗::Vector{Interval}, prototypes::Vector{Vector{Interval}}, i::Int; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	h = 1 / (1 - m)
	c = length(prototypes)
	p = length(X⃗)

	C⃗ᵢ = prototypes[i]
	num_endpoints = (c + 1) * p

	u_min = Inf
	u_max = -Inf
	for idx_endpoint = 1:2^num_endpoints
		endpoints = get_sequential_endpoints(num_endpoints, idx_endpoint)

		dⱼᵢ = 0
		for j = 1:p
			c_idx = p + (i-1)*p + j
			dⱼᵢ += (X⃗[j][endpoints[j]] - C⃗ᵢ[j][endpoints[c_idx]])^2
		end
		if dⱼᵢ != 0
			dⱼᵢ ^= h
		end

		∑ = 0
		for k = 1:c
			C⃗ₖ = prototypes[k]
			dⱼₖ = 0
			for j = 1:p
				c_idx = p + (k-1)*p + j
				dⱼₖ += (X⃗[j][endpoints[j]] - C⃗ₖ[j][endpoints[c_idx]])^2
			end
			if dⱼₖ != 0
				dⱼₖ ^= h
			end
			∑ += dⱼₖ
		end
		u = ∑ == 0 ? nothing : dⱼᵢ / ∑
		u_min = min(u_min, u)
		u_max = max(u_max, u)
	end
	Interval(u_min, u_max)
end

function u_lfcm(X⃗::Vector{Interval}, prototypes::Vector{Vector{Interval}}, i::Int; m::Real=2.0, d_method::String="interval")
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	c = length(prototypes)

	D_squared = Vector{Interval}(undef, c)

	for k = 1:c
		if d_method == "dsw"
			D_squared[k] = d_dsw(X⃗, prototypes[k], squared=true)
		else
			D_squared[k] = d_interval(X⃗, prototypes[k], squared=true)
		end
	end

	u_lfcm(D_squared, i; m=m)
end

function u_lfcm(D_squared::Vector{Interval}, i::Int; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	h = 1 / (1 - m)
	c = length(D_squared)

	u_a = -1
	u_b = -1

	# corollary 1
	if D_squared[i].right != 0
		for k = 1:c
			if k != i && D_squared[k].left == 0
				u_a = 0
				break
			end
		end
	end

	# corollary 2
	if D_squared[i].right == D_squared[i].left == 0
		for k = 1:c
			if k != i && D_squared[k].left != 0
				u_a = 1
				break
			end
		end
	end

	# corollary 3
	if D_squared[i].left != 0
		for k = 1:c
			if k != i && D_squared[k].right == 0
				u_b = 0
				break
			end
		end
	end

	# corollary 4
	if D_squared[i].left == 0
		for k = 1:c
			if k != i && D_squared[k].right != 0
				u_b = 1
				break
			end
		end
	end

	d2jih = (D_squared[i].right)^h
	d1jih = (D_squared[i].left)^h

	u_a_denom = d2jih
	u_b_denom = d1jih
	for k = 1:c
		if k != i
			u_a_denom += (D_squared[k].left)^h
			u_b_denom += (D_squared[k].right)^h
		end
	end
	if u_a == -1
		u_a = d2jih / u_a_denom
	end
	if u_b == -1
		u_b = d1jih / u_b_denom
	end

	Interval(u_a, u_b)
end

function c_dsw(X::Vector{FuzzyVector}, u::Matrix{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	levels = X[1][1].levels
	num_levels = length(levels)
	N, c = size(u)
	p = length(X[1])

	C = Vector{FuzzyVector}(undef, c)
	for i = 1:c
		C⃗ = Vector{FuzzyNumber}(undef, p)
		for j = 1:p
			grades = Vector{Interval}(undef, num_levels)
			for lvl = 1:num_levels
				num_endpoints = N + N

				c_min = Inf
				c_max = -Inf
				@floop for idx_endpoint = 1:2^num_endpoints
					endpoints = get_sequential_endpoints(num_endpoints, idx_endpoint)
					numerator = 0
					for k = 1:N
						u_i = N + k
						numerator += u[k, i][lvl][endpoints[u_i]]^m * X[k][j][lvl][endpoints[k]]
					end

					∑ = 0
					for k = 1:N
						u_i = N + k
						∑ += u[k, i][lvl][endpoints[u_i]]
					end
					cᵢⱼ = ∑ == 0 ? nothing : numerator / ∑
					@reduce() do (c_min; cᵢⱼ)
						if cᵢⱼ < c_min
							c_min = cᵢⱼ
						end
					end
					@reduce() do (c_max; cᵢⱼ)
						if cᵢⱼ > c_max
							c_max = cᵢⱼ
						end
					end
				end
				grades[lvl] = Interval(c_min, c_max)
			end
			C⃗[j] = FuzzyNumber(levels, grades)
		end
		C[i] = FuzzyVector(C⃗)
	end
	C
end

function km_iwa(X::Vector{Interval}, u::Vector{Interval}; bound::String, m::Real=2.0) # Karnik and Mendel's interval weighted average
    N = size(X, 1)
    x = Vector{Float64}(undef, N)
    for k = 1:N
		if bound == "lower"
			x[k] = round(X[k].left, digits=6) # a_i
		elseif bound == "upper"
			x[k] = round(X[k].right, digits=6) # b_i
		else
			throw(ArgumentError("bound must be \"lower\" or \"upper\""))
		end
    end
    x_sorted = sort(x)
    sorted_indices = sortperm(x)

    w = Vector{Float64}(undef, N)
    c = Vector{Float64}(undef, N)
    d = Vector{Float64}(undef, N)
    for k = 1:N
        w[k] = round(mid(u[k]), digits=6)
        c[k] = round(u[k].left, digits=6)
        d[k] = round(u[k].right, digits=6)
    end
    y′ = round(sum(w.^m .* x) / sum(w.^m), digits=6)

	for iteration = 1:N+1
		iteration <= N || error("KM IWA -> infinite loop")

		# println(x_sorted)
        _k = nothing
        for k = 1:N-1
            if x_sorted[k] <= y′ && y′ <= x_sorted[k+1]
				# println(x_sorted[k], " <= ", y′, " <= ", x_sorted[k+1], " | k=$k")
                _k = k
                break
            end
        end
		!isnothing(_k) || error("Can't find _k")

		# println("_k ", _k)
		left_indices = sorted_indices[1:_k]
		right_indices = sorted_indices[_k+1:end]
		if bound == "lower"
			w[left_indices] = d[left_indices]
			w[right_indices] = c[right_indices]
		elseif bound == "upper"
			w[left_indices] = c[left_indices]
			w[right_indices] = d[right_indices]
		end
        w = max.(w, 0) # handle the weird case that w is negative??
		# println(bound, " ", w)

        y_k = round(sum(w.^m .* x) / sum(w.^m), digits=6)
        if isapprox(y′, y_k, atol=1e-5)
            break
		# else
		# 	println("$y′ ≂̸ $y_k")
        end
        y′ = y_k
    end
    y′
end

function c_karnik(X::Vector{FuzzyVector}, u::Matrix{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	levels = X[1][1].levels
	num_levels = length(levels)
	N, c = size(u)
	p = length(X[1])

	C = Vector{FuzzyVector}(undef, c)
	for i = 1:c
		C⃗ = Vector{FuzzyNumber}(undef, p)
		for j = 1:p
			grades = Vector{Interval}(undef, num_levels)
			for (lvl, α) in enumerate(levels)
				X_cut = cut(X, j, α)
				u_cut = cut(FuzzyVector(u[:, i]), α)
				c_left = km_iwa(X_cut, u_cut; bound="lower", m=m)
				c_right = km_iwa(X_cut, u_cut; bound="upper", m=m)
				grades[lvl] = Interval(c_left, c_right)
			end
			C⃗[j] = FuzzyNumber(levels, grades)
		end
		C[i] = FuzzyVector(C⃗)
	end
	C
end

function c_karnik(X::Vector{FuzzyVector}, u::Vector{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	levels = X[1][1].levels
	num_levels = length(levels)
	c = 1
	p = length(X[1])

	C⃗ = Vector{FuzzyNumber}(undef, p)
	for j = 1:p
		grades = Vector{Interval}(undef, num_levels)
		for (lvl, α) in enumerate(levels)
			X_cut = cut(X, j, α)
			u_cut = cut(FuzzyVector(u), α)
			c_left = km_iwa(X_cut, u_cut; bound="lower", m=m)
			c_right = km_iwa(X_cut, u_cut; bound="upper", m=m)
			grades[lvl] = Interval(c_left, c_right)
		end
		C⃗[j] = FuzzyNumber(levels, grades)
	end
	FuzzyVector(C⃗)
end

function d_interval(A⃗::FuzzyVector, B⃗::FuzzyVector; squared::Bool=false)
	if A⃗ == B⃗
		𝟎 = SingletonFuzzyNumber(A⃗[1].levels, number=0)
		return 𝟎
	elseif length(A⃗) ≠ length(B⃗)
		return false
	end

	levels = A⃗[1].levels
	num_levels = length(levels)

	grades = Vector{Interval}(undef, num_levels)
	for (lvl, α) in enumerate(levels)
		A_cut = cut(A⃗, α)
		B_cut = cut(B⃗, α)
		grades[lvl] = d_interval(A_cut, B_cut; squared=squared)
	end
	FuzzyNumber(levels, grades)
end

function d_interval(A⃗::Vector{Interval}, B⃗::Vector{Interval}; squared::Bool=false)
	p = length(A⃗)

	d = Interval(0)
	for i = 1:p
		d += (A⃗[i] - B⃗[i])^2
	end
	if (!squared) d ^= 0.5 end

	d
end

function u_interval(squared_fuzzy_distances::Vector{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	c = length(squared_fuzzy_distances)
	levels = squared_fuzzy_distances[1].levels
	num_levels = length(levels)

	I = Int64[]
	𝟎 = SingletonFuzzyNumber(levels, number=0)
	for i = 1:c
		if squared_fuzzy_distances[i] == 𝟎
			push!(I, i)
		end
	end

	u⃗ = Vector{FuzzyNumber}(undef, c)
	if isempty(I)
		for i = 1:c
			grades = Vector{Interval}(undef, num_levels)
			for (lvl, α) in enumerate(levels)
				d²_cut = cut(squared_fuzzy_distances, α)
				grades[lvl] = u_interval(d²_cut, i; m=m)
			end
			u⃗[i] = FuzzyNumber(levels, grades)
		end
	else
		for i = 1:c
			if i ∉ I
				u⃗[i] = 𝟎
			else
				β = SingletonFuzzyNumber(levels, number=1 / length(I))
				u⃗[i] = β
			end
		end
	end

	u⃗
end

function u_interval(squared_fuzzy_distances::Vector{Interval}, i::Int; m::Real=2.0)
	h = 1 / (1 - m)
	squared_fuzzy_distances[i]^h / sum(squared_fuzzy_distances.^h)
end

function c_interval(X⃗::Vector{FuzzyVector}, u::Vector{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	levels = X⃗[1][1].levels
	num_levels = length(levels)
	N = length(X⃗)
	p = length(X⃗[1])

	C⃗ = Vector{FuzzyNumber}(undef, p)
	for j = 1:p # each dimension j
		membership_grades = Vector{Interval}(undef, num_levels)
		for lvl = 1:num_levels
			sum_numerator = Interval(0)
			sum_denominator = Interval(0)
			for k = 1:N
				sum_numerator += (u[k][lvl]^m * X⃗[k][j][lvl])
				sum_denominator += u[k][lvl]^m
			end
			membership_grades[lvl] = sum_numerator / sum_denominator
		end
		C⃗[j] = FuzzyNumber(levels, membership_grades)
	end
	C⃗ = FuzzyVector(C⃗)
	C⃗
end

function d2(X⃗::Vector{FuzzyVector}, C⃗::Vector{FuzzyVector})
	N = length(X⃗)
	c = length(C⃗)
	squared_fuzzy_distances = Matrix{FuzzyNumber}(undef, N, c)
	for i = 1:c
		for j = 1:N
			squared_fuzzy_distances[j, i] = FuzzySets.d_interval(X⃗[j], C⃗[i], squared=true)
		end
	end
	squared_fuzzy_distances
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0ad6e310-97e2-11eb-28bc-99ff80496ecb
using Revise, FourierTools, ImageCore, TestImages, PlutoUI, ImageShow

# ╔═╡ 19339dd6-97e2-11eb-2449-05f0a5ca2f9d
begin
	img = Float32.(testimage("fabio_512_gray"))
	z = zeros(Float32, (768, 768))
	FourierTools.center_set!(z, img)
end

# ╔═╡ 1e72338e-97e2-11eb-0ca1-374f5ab2d2da
@bind Δ Slider(-1024:1:1024)

# ╔═╡ 1aeea42c-97e2-11eb-3af3-d377949794ed
md"Δ=$Δ"

# ╔═╡ 3752a30c-97e2-11eb-3469-dbf857775870
Gray.(FourierTools.shear(z, Δ))

# ╔═╡ 39b3ba50-97e2-11eb-3153-9757657de483
Gray.(FourierTools.shear(img, Δ, 2, 1))

# ╔═╡ Cell order:
# ╠═0ad6e310-97e2-11eb-28bc-99ff80496ecb
# ╠═19339dd6-97e2-11eb-2449-05f0a5ca2f9d
# ╠═1aeea42c-97e2-11eb-3af3-d377949794ed
# ╠═1e72338e-97e2-11eb-0ca1-374f5ab2d2da
# ╠═3752a30c-97e2-11eb-3469-dbf857775870
# ╠═39b3ba50-97e2-11eb-3153-9757657de483


text:
function invariants_check(x)
prim=zeros(6,1);
pv=zeros(8,1);

v=zeros(12,1);
prim[1]=x[1] + x[2] + x[3]

prim[2]=x[4] + x[5] + x[6]

prim[3]=x[1]^2 + x[2]^2 + x[3]^2

prim[4]=x[1]*x[4] + x[2]*x[6] + x[3]*x[5]

prim[5]=x[4]^3 + x[5]^3 + x[6]^3

prim[6]=x[1]^3 + x[2]^3 + x[3]^3 + x[4]^2*x[5] + x[4]*x[6]^2 + x[5]^2*x[6]

pv[1]=x[1]*x[5] + x[2]*x[4] + x[3]*x[6]

pv[2]=x[4]^2 + x[5]^2 + x[6]^2

pv[3]=x[1]^2*x[3] + x[1]*x[2]^2 + x[2]*x[3]^2

pv[4]=x[1]*x[2]*x[4] + x[1]*x[3]*x[5] + x[2]*x[3]*x[6]

pv[5]=x[1]*x[4]^2 + x[2]*x[6]^2 + x[3]*x[5]^2

pv[6]=x[1]^2*x[5] + x[2]^2*x[4] + x[3]^2*x[6]

pv[7]=x[1]*x[4]*x[5] + x[2]*x[4]*x[6] + x[3]*x[5]*x[6]

pv[8]=x[4]^2*x[6] + x[4]*x[5]^2 + x[5]*x[6]^2

v[1] = 1
v[2] = pv[1]
v[3] = pv[2]
v[4] = pv[3]
v[5] = pv[4]
v[6] = pv[5]
v[7] = pv[6]
v[8] = pv[7]
v[9] = pv[8]
v[10] = pv[1]*pv[2]
v[11] = pv[2]^2
v[12] = pv[4]^2
return prim, v, pv

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using Turing
using Test

obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

@model newinterface(obs) = begin
  p ~ Beta(2,2)
  for i = 1:length(obs)
    obs[i] ~ Bernoulli(p)
  end
  p
end

# Turing.TURING[:modelex]

# newinterface(data)
#
# ga = VarInfo()
# sampler = HMCSampler{HMC}(HMC(100, 1.5, 3))
# ga = newinterface(data, ga, sampler)
# newinterface
#
# chain = sample(newinterface, HMC(100, 1.5, 3))

chain = sample(newinterface(obs), HMC(100, 0.75, 3, :p, :x); chunk_size=2)


text:
"""
    $(TYPEDEF)

An image file directory is a sorted collection of the tags representing this
plane in the TIFF file. They behave like dictionaries, so given an IFD called
`ifd`, we can add new tags as follows:

```jldoctest; setup = :(ifd = TiffImages.IFD(UInt32))
julia> ifd[TiffImages.IMAGEDESCRIPTION] = "Some details";

julia> ifd[TiffImages.IMAGEWIDTH] = 512;

julia> ifd
IFD, with tags:
	Tag(IMAGEWIDTH, 512)
	Tag(IMAGEDESCRIPTION, "Some details")
```
"""
struct IFD{O <: Unsigned}
    tags::OrderedDict{UInt16, Tag}
end

IFD(o::Type{O}) where {O <: Unsigned} = IFD{O}(OrderedDict{UInt16, Tag}())

Base.length(ifd::IFD) = length(ifd.tags)
Base.keys(ifd::IFD) = keys(ifd.tags)
Base.iterate(ifd::IFD) = iterate(ifd.tags)
Base.iterate(ifd::IFD, n::Integer) = iterate(ifd.tags, n)
Base.getindex(ifd::IFD, key::TiffTag) = getindex(ifd, UInt16(key))
Base.getindex(ifd::IFD{O}, key::UInt16) where {O} = getindex(ifd.tags, key)
Base.in(key::TiffTag, v::IFD) = in(UInt16(key), v)
Base.in(key::UInt16, v::IFD) = in(key, keys(v))
Base.delete!(ifd::IFD, key::TiffTag) = delete!(ifd, UInt16(key))
Base.delete!(ifd::IFD, key::UInt16) = delete!(ifd.tags, key)

Base.setindex!(ifd::IFD, value::Tag, key::UInt16) = setindex!(ifd.tags, value, key)
Base.setindex!(ifd::IFD, value::Tag, key::TiffTag) = setindex!(ifd.tags, value, UInt16(key))

Base.setindex!(ifd::IFD, value, key::TiffTag) = setindex!(ifd, value, UInt16(key))
function Base.setindex!(ifd::IFD{O}, value, key::UInt16) where {O <: Unsigned}
    setindex!(ifd, Tag(key, value), UInt16(key))
end

function load!(tf::TiffFile, ifd::IFD)
    for idx in keys(ifd)
        ifd[idx] = load(tf, ifd[idx])
    end
end

function Base.show(io::IO, ifd::IFD)
    print(io, "IFD, with tags: ")
    for tag in sort(collect(keys(ifd)))
        print(io, "\n\t", ifd[tag])
    end
end

function Base.read(tf::TiffFile{O}, ::Type{IFD}) where O <: Unsigned
    # Regular TIFF's use 16bits instead of 32 bits for entry data
    N = O == UInt32 ? read(tf, UInt16) : read(tf, O)

    entries = OrderedDict{UInt16, Tag}()

    for i in 1:N
        tag = read(tf, Tag)
        entries[tag.tag] = tag
    end

    next_ifd = Int(read(tf, O))
    IFD{O}(entries), next_ifd
end

function Base.iterate(file::TiffFile{O}) where {O}
    seek(file.io, file.first_offset)
    iterate(file, (read(file, IFD)))
end

"""
    iterate(file, state) -> IFD, Int

Advances the iterator to the next IFD.

**Output**
- `Vector{Int}`: Offsets within file for all strips corresponding to the current
   IFD
- `Int`: Offset of the next IFD
"""
function Base.iterate(file::TiffFile, state::Tuple{Union{IFD{O}, Nothing}, Int}) where {O}
    curr_ifd, next_ifd_offset = state
    # if current element doesn't exist, exit
    (curr_ifd == nothing) && return nothing
    (next_ifd_offset <= 0) && return (curr_ifd, (nothing, 0))

    seek(file.io, next_ifd_offset)
    next_ifd, next_ifd_offset = read(file, IFD)

    return (curr_ifd, (next_ifd, next_ifd_offset))
end

struct IFDLayout
    nsamples::Int
    nrows::Int
    ncols::Int
    nbytes::Int
    readtype::DataType
    rawtype::DataType
    mappedtype::DataType
    compression::CompressionType
    interpretation::PhotometricInterpretations
end

function output(ifd::IFD)
    nrows = Int(ifd[IMAGELENGTH].data)
    ncols = Int(ifd[IMAGEWIDTH].data)

    samplesperpixel = Int(ifd[SAMPLESPERPIXEL].data)
    sampleformats = fill(UInt16(0x01), samplesperpixel)
    if SAMPLEFORMAT in ifd
        sampleformats = ifd[SAMPLEFORMAT].data
    end

    interpretation = Int(ifd[PHOTOMETRIC].data)

    strip_nbytes = ifd[STRIPBYTECOUNTS].data
    nbytes = Int(sum(strip_nbytes))
    bitsperpixel = ifd[BITSPERSAMPLE].data
    rawtypes = Set{DataType}()
    mappedtypes = Set{DataType}()
    for i in 1:samplesperpixel
        rawtype = rawtype_mapping[(SampleFormats(sampleformats[i]), bitsperpixel[i])]
        push!(rawtypes, rawtype)
        if rawtype <: Unsigned
            push!(mappedtypes, Normed{rawtype, sizeof(rawtype)*8})
        elseif rawtype <: Signed
            push!(mappedtypes, Fixed{rawtype, sizeof(rawtype)*8-1})
        else
            push!(mappedtypes, rawtype)
        end
    end
    (length(rawtypes) > 1) && error("Variable per-pixel storage types are not yet supported")
    rawtype = first(rawtypes)
    readtype = rawtype

    compression = COMPRESSION_NONE
    if COMPRESSION in ifd
        compression = CompressionType(ifd[COMPRESSION].data)
    end

    if compression != COMPRESSION_NONE
        # recalculate nbytes if the data is compressed since the inflated data
        # is most likely larger than the bytes on disk
        nbytes = nrows*ncols*samplesperpixel*sizeof(rawtype)
        readtype = UInt8
    end

    IFDLayout(
        samplesperpixel, nrows, ncols,
        nbytes,
        readtype,
        rawtype,
        first(mappedtypes),
        compression,
        PhotometricInterpretations(interpretation))
end

function Base.read!(target::AbstractArray{T, N}, tf::TiffFile, ifd::IFD) where {T, N}
    layout = output(ifd)

    rowsperstrip = layout.nrows
    (ROWSPERSTRIP in keys(ifd)) && (rowsperstrip = ifd[ROWSPERSTRIP].data)
    nstrips = ceil(Int, layout.nrows / rowsperstrip)

    strip_nbytes = ifd[STRIPBYTECOUNTS].data

    if layout.compression != COMPRESSION_NONE
        # strip_nbytes is the number of bytes pre-inflation so we need to
        # calculate the expected size once decompressed and update the values
        strip_nbytes = fill(rowsperstrip*layout.ncols, length(strip_nbytes))
        strip_nbytes[end] = (layout.nrows - (rowsperstrip * (nstrips-1))) * layout.ncols
    end

    strip_offsets = ifd[STRIPOFFSETS].data

    if PLANARCONFIG in ifd
        planarconfig = ifd[PLANARCONFIG].data
        (planarconfig != 1) && error("Images with data stored in planar format not yet supported")
    end

    if nstrips > 1
        startbyte = 1
        for i in 1:nstrips
            seek(tf, strip_offsets[i])
            nbytes = Int(strip_nbytes[i] / sizeof(T))
            read!(tf, view(target, startbyte:(startbyte+nbytes-1)), layout.compression)
            startbyte += nbytes
        end
    else
        seek(tf, strip_offsets[1])
        read!(tf, target, layout.compression)
    end
end

function Base.write(tf::TiffFile{O}, ifd::IFD{O}) where {O <: Unsigned}
    N = length(ifd)
    O == UInt32 ? write(tf, UInt16(N)) : write(tf, UInt64(N))

    sort!(ifd.tags)

    # keep track of which tags are too large to fit in the IFD slot and need a
    # remote location for their data
    remotedata = OrderedDict{Tag, Vector{Int}}()
    for (key, tag) in ifd
        pos = position(tf.io)
        if !write(tf, tag)
            remotedata[tag] = [pos]
        end
    end

    # end position, write a zero by default, but this should be updated if any
    # more IFDs are written
    ifd_end_pos = position(tf.io)
    write(tf, O(0))

    for (tag, poses) in remotedata
        data_pos = position(tf.io)
        # add NUL terminator to the end of Strings that don't have it already
        data = (eltype(tag) == String && !endswith(tag.data, '\0')) ? tag.data * "\0" : tag.data
        write(tf, data)
        push!(poses, data_pos)
    end

    for (tag, poses) in remotedata
        orig_pos, data_pos = poses
        seek(tf, orig_pos)
        write(tf, tag, data_pos)
    end

    seek(tf, ifd_end_pos)

    return ifd_end_pos
end


text:
using StatsBase
using Statistics
using Test

function learn_MeanFunctions()
    data = [1:5...]

    @test 3.0 == mean(data)
end

function learn_ScalarStatistics()
    data = [1:5...]
    count = length(data)

    # var
    @test 2.0 == var(data, corrected=false) # population: /n from Statistics
    @test 2.0 == var(data, uweights(count)) # population: /n from StatsBase

    @test 2.5 == var(data)                                  # sample: /(n-1) from Statistics
    @test 2.5 == var(data, uweights(count), corrected=true) # sample: /(n-1) from StatsBase

    # std
    @test sqrt(2.0) == std(data, corrected=false) # population: /n from Statistics
    @test sqrt(2.0) == std(data, uweights(count)) # population: /n from Statistics

    @test sqrt(2.5) == std(data)
    @test sqrt(2.5) == std(data, uweights(count), corrected=true)

    # zscore = (data - mu)/sigma
    mu = 2.0          # mean, from above
    sigma = sqrt(2.0) # std, from above
    @test zscore(data, mu, sigma) == [-0.7071067811865475, 0.0, 0.7071067811865475, 1.414213562373095, 2.1213203435596424]
end

@testset "StatsBase" begin
    learn_MeanFunctions()
    learn_ScalarStatistics()
end

nothing


text:
L = (4,4,4,4)
β = 6
NTRACE = 3
#gparam = Setup_Gauge_action(β)
gparam =  GaugeActionParam_standard(β,NTRACE)

BoundaryCondition=[1,1,1,-1]
Nwing = 1
initial="cold"
NC =3

hop= 0.141139#Hopping parameter
r= 1#Wilson term
eps= 1e-19
Dirac_operator= "Wilson"
MaxCGstep= 3000

fparam = FermiActionParam_Wilson(hop,r,eps,Dirac_operator,MaxCGstep)


text:
## compose_matrix.jl : methods to compose example matrices

"""
    laplace1d(n)

1d Laplacian on [0,1] using n points
"""
laplace1d(n) = stencil2mat(SArray{Tuple{3}}([-n^2 2n^2 -n^2]), n)

"""
    laplace2d(n, m)

2d Laplacian on [0,1]^2 using n x m points
"""
laplace2d(n, m) = stencil2mat(SArray{Tuple{3,3}}([0 -m^2 0 -n^2 2(n^2+m^2) -n^2 0 -m^2 0]), n, m)

"""
    laplace3d(n, m, l)

3d Laplacian on [0,1]^3 using n x m x l points
"""
laplace3d(n, m, l) = stencil2mat(SArray{Tuple{3,3,3}}([0 0 0 0 -l^2 0 0 0 0 0 -m^2 0 -n^2 2(m^2+n^2+l^2) -n^2 0 -m^2 0 0 0 0 0 -l^2 0 0 0 0]), n, m, l)

"""
    elliptic1d(k)

Generate system matrix for general 1d elliptic PDE defined on [0,1].
The diffusion coefficient k must given on an equidistant grid with nodes {1,...,n-1}.
The step size is thus h = 1/n.
"""
function elliptic1d(k)
    n = length(k) + 1
    stencil2mat(SArray{Tuple{3}}([-n^2 2n^2 -n^2]), k)
end

"""
    elliptic2d(k)

Generate system matrix for general 2d elliptic PDE defined on [0,1]^2.
The diffusion coefficient k must given on an equidistant grid with nodes {1,...,n-1} x {1,..,m-1}. The step sizes are thus h_x = 1/n and h_y = 1/m.
"""
function elliptic2d(k)
    n, m = size(k) .+ 1
    stencil2mat(SArray{Tuple{3,3}}([0 -m^2 0 -n^2 2(n^2+m^2) -n^2 0 -m^2 0]), k)
end

"""
    elliptic3d(k)

Generate system matrix for general 3d elliptic PDE defined on [0,1]^3.
The diffsion coefficient k must given on an equidistant grid with nodes {1,...,n-1} x {1,..,m-1} x {0,...,l-1}. The step sizes are thus h_x = 1/n, h_y = 1/m and h_z = 1/l.
"""
function elliptic3d(k)
    n,m,l = size(k) .+ 1
    stencil2mat(SArray{Tuple{3,3,3}}([0 0 0 0 -l^2 0 0 0 0 0 -m^2 0 -n^2 2(n^2+m^2+l^2) -n^2 0 -m^2 0 0 0 0 0 -l^2 0 0 0 0]), k)
end

# construction for constant stencil problems
function stencil2mat(stencil::SArray, n::Int...)
    R = CartesianIndices(n.-1)
    I1, Iend = extrema(R)
    Is = Int[]
    Js = Int[]
    Vs = Float64[]
    @inbounds for I in R
        is,js,vs = _stencil2mat(stencil,I,I1,Iend)
        push!(Is,is...)
        push!(Js,js...)
        push!(Vs,vs...)
    end
    sparse(Is,Js,Vs)
end

# main driver code - constant stencil
@noinline function _stencil2mat(stencil::SArray,I,I1,Iend)
    R = CartesianIndices(UnitRange.(max(I1, I-I1).I, min(Iend, I+I1).I)) # TODO: #29440 replaces UnitRange with : (colon), see https://discourse.julialang.org/t/psa-replacement-of-ind2sub-sub2ind-in-julia-0-7/14666/9 [repeats 3 times]
    Is = fill(0,length(R))
    Js = fill(0,length(R))
    Vs = fill(0.,length(R))
    sz = LinearIndices(Iend.I)
    @inbounds for (i,J) in enumerate(R)
        idx = J-I+2I1
        Is[i] = sz[I]
        Js[i] = sz[J]
        Vs[i] = stencil[idx]
    end
    Is,Js,Vs
end

# main driver code for length 7 constant 1d stencils (used in `Cubic()` interpolation)
function _stencil2mat(stencil::SVector{7},I,I1,Iend)
    R = CartesianIndices(UnitRange.(max(I1, I-3I1).I, min(Iend, I+3I1).I))
    Is = fill(0,length(R))
    Js = fill(0,length(R))
    Vs = fill(0.,length(R))
    sz = LinearIndices(Iend.I)
    @inbounds for (i,J) in enumerate(R)
        idx = J-I+4I1
        Is[i] = sz[I]
        Js[i] = sz[J]
        Vs[i] = stencil[idx]
    end
    if I.I[1] == 2 # correction for boundary
        Vs[1] += 1/16
    elseif I.I[1] == Iend.I[1]-1
        Vs[end] += 1/16
    end
    Is,Js,Vs
end

# construction for non-constant stencil problems
function stencil2mat(stencil::SArray, k::AbstractArray)
    n = size(k)
    R = CartesianIndices(n)
    I1, Iend = extrema(R)
    Is = Int[]
    Js = Int[]
    Vs = Float64[]
    for I in R
        is,js,vs = _stencil2mat(k,stencil,I,I1,Iend)
        push!(Is,is...)
        push!(Js,js...)
        push!(Vs,vs...)
    end
    sparse(Is,Js,Vs)
end

# main driver code - non-constant stencil
function _stencil2mat(k::AbstractArray,stencil::SArray,I,I1,Iend)
    R = CartesianIndices(UnitRange.(max(I1, I-I1).I, min(Iend, I+I1).I))
    #R = max(I1, I-I1):min(Iend, I+I1)
    Is = fill(0, length(R))
    Js = fill(0, length(R))
    Vs = fill(0., length(R))
    sz = LinearIndices(Iend.I)
    for (i,J) in enumerate(R)
        idx = J-I+2I1
        Is[i] = sz[I]
        Js[i] = sz[J]
        Vs[i] = k[J]*stencil[idx]
    end
    Is,Js,Vs
end


text:
module RIPEMD

##############################################################################
##
## TODO
##
##############################################################################

# - Fix bug for more than one chunk
# - Refactor RIPEMD-160 and SHA-256 to share boilerplate (E.g., read/write, padding)
# - Digest -> init, update, finalize

##############################################################################
##
## Refereces
##
##############################################################################

# - https://maemo.gitorious.org/maemo-pkg/python-crypto/source/8651b0eace17916fe7ba14923dbe4054f255ec2a:lib/Crypto/Hash/RIPEMD160.py
# - https://github.com/bitcoin/bitcoin/blob/master/src/crypto/ripemd160.cpp

##############################################################################
##
## Exports
##
##############################################################################

export ripemd160

##############################################################################
##
## Notes
##
##############################################################################

# WARNING: This module is in need of some clean up

# Pseudocode: http://homes.esat.kuleuven.be/~bosselae/ripemd/rmd160.txt
# Additional notes: https://en.bitcoin.it/wiki/RIPEMD-160
#
# RIPEMD-160 is an iterative hash function that operates on 32-bit words.
# The round function takes as input a 5-word chaining variable and a 16-word
# message block and maps this to a new chaining variable. All operations are
# defined on 32-bit words. Padding is identical to that of MD4.

##############################################################################
##
## Definitions
##
##############################################################################

const BLOCK_SIZE = 64  # SHA256 operates on 64 (512 bit) blocks
const WORDS_PER_BLOCK = 16

# Nonlinear functions at bit level: exor, mux, -, mux, -
const f = [(x, y, z) -> x $ y $ z,
           (x, y, z) -> (x & y) | (~x & z),
           (x, y, z) -> (x | ~y) $ z,
           (x, y, z) -> (x & z) | (y & ~z),
           (x, y, z) -> x $ (y | ~z)]

# Added constants (hexadecimal)
# NOTE: First 32 bits of the fractional parts of the square (K0)/cube (K1)
# roots of the first 4 primes, 2 through 7.
const K0 = [0x00000000, 0x5A827999, 0x6ED9EBA1, 0x8F1BBCDC, 0xA953FD4E]
const K1 = [0x50A28BE6, 0x5C4DD124, 0x6D703EF3, 0x7A6D76E9, 0x00000000]

# Selection of message word
# NOTE: Julia is 1-indexed
const r0 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
            8, 5, 14, 2, 11, 7, 16, 4, 13, 1, 10, 6, 3, 15, 12, 9,
            4, 11, 15, 5, 10, 16, 9, 2, 3, 8, 1, 7, 14, 12, 6, 13,
            2, 10, 12, 11, 1, 9, 13, 5, 14, 4, 8, 16, 15, 6, 7, 3,
            5, 1, 6, 10, 8, 13, 3, 11, 15, 2, 4, 9, 12, 7, 16, 14]
const r1 = [5, 15, 8, 1, 10, 3, 12, 5, 14, 7, 16, 9, 2, 11, 4, 13,
            7, 12, 4, 8, 1, 14, 6, 11, 15, 16, 9, 13, 5, 10, 2, 3,
            16, 6, 2, 4, 8, 15, 7, 10, 12, 9, 13, 3, 11, 1, 5, 14,
            9, 7, 5, 2, 4, 12, 16, 1, 6, 13, 3, 14, 10, 8, 11, 15,
            13, 16, 11, 5, 2, 6, 9, 8, 7, 3, 14, 15, 1, 4, 10, 12]

# Amount for rotate left (rol)
# NOTE: Julia is 1-indexed
const s0 = [11, 14, 15, 12, 5, 8, 7, 9, 11, 13, 14, 15, 6, 7, 9, 8,
            7, 6, 8, 13, 11, 9, 7, 15, 7, 12, 15, 9, 11, 7, 13, 12,
            11, 13, 6, 7, 14, 9, 13, 15, 14, 8, 13, 6, 5, 12, 7, 5,
            11, 12, 14, 15, 14, 15, 9, 8, 9, 14, 5, 6, 8, 6, 5, 12,
            9, 15, 5, 11, 6, 8, 13, 12, 5, 12, 13, 14, 11, 8, 5, 6]
const s1 = [8, 9, 9, 11, 13, 15, 15, 5, 7, 7, 8, 11, 14, 14, 12, 6,
            9, 13, 15, 7, 12, 8, 9, 11, 7, 7, 12, 7, 6, 15, 13, 11,
            9, 7, 15, 11, 8, 6, 6, 14, 12, 13, 5, 14, 13, 13, 7, 5,
            15, 5, 8, 11, 14, 14, 6, 14, 6, 9, 12, 9, 12, 5, 15, 8,
            8, 5, 12, 9, 12, 5, 14, 6, 8, 13, 6, 5, 15, 13, 11, 11]

##############################################################################
##
## Initialize helper functions
##
##############################################################################

# Hard-coding for 32-bit unsigned ints
# Technically, we could use the '>>' operator for right shift because
# we're only using unsigned ints, but we explicitly use logical right shift
ROTRIGHT(num, shift) = (num >>> shift) | (num << (32 - shift))
ROTLEFT(num, shift) = (num << shift) | (num >>> (32 - shift))

##############################################################################
##
## Function definitions
##
##############################################################################

function transform!(state, block)
  A0, B0, C0, D0, E0 = state
  A1, B1, C1, D1, E1 = state

  X = zeros(Uint32, WORDS_PER_BLOCK)

  for i = 1:WORDS_PER_BLOCK
    X[i] = uint32(block[4 * (i - 1) + 4]) << 24 +
           uint32(block[4 * (i - 1) + 3]) << 16 +
           uint32(block[4 * (i - 1) + 2]) << 8 +
           uint32(block[4 * (i - 1) + 1])
  end

  for j = 1:80
    # TODO: Replace with cld(j, 16) in Julia v0.4
    div, rem = divrem(j, 16)
    i = div + (rem == 0 ? 0 : 1)
    k = 5 - i + 1

    T = uint32(A0 + f[i](B0, C0, D0) + X[r0[j]] + K0[i])
    T = uint32(ROTLEFT(T, s0[j]) + E0)

    A0 = E0
    E0 = D0
    D0 = ROTLEFT(C0, 10)
    C0 = B0
    B0 = T

    T = uint32(A1 + f[k](B1, C1, D1) + X[r1[j]] + K1[i])
    T = uint32(ROTLEFT(T, s1[j]) + E1)

    A1 = E1
    E1 = D1
    D1 = ROTLEFT(C1, 10)
    C1 = B1
    B1 = T
  end

  T = uint32(state[2] + C0 + D1)
  state[2] = uint32(state[3] + D0 + E1)
  state[3] = uint32(state[4] + E0 + A1)
  state[4] = uint32(state[5] + A0 + B1)
  state[5] = uint32(state[1] + B0 + C1)
  state[1] = uint32(T)

  return
end

function ripemd160(msg::ASCIIString; is_hex=false)

  if is_hex
    len = int(length(msg) / 2)
    result = zeros(Uint8, len)
    for i = 1:len
      result[i] = uint8(parseint(msg[2 * i - 1: 2 * i],16))
    end
    msg = result
  else
    # We only want byte array literal (i.e., character array)
    msg = msg.data
  end

  # Get original length and bit lengths
  len = length(msg)
  bitlen = len * 8

  # Append the bit '1' to the message.
  append!(msg, [0x80])

  state = [0x67452301, 0xEFCDAB89, 0x98BADCFE, 0x10325476, 0xC3D2E1F0]

  # Divide up message into blocks of BLOCK_SIZE = 512 bits
  # and run through transformation
  while length(msg) >= BLOCK_SIZE
    transform!(state, msg[1:BLOCK_SIZE])
    msg = msg[BLOCK_SIZE + 1:end]
  end

  # Get the number of characters untransformed
  rem = length(msg)

  # If there are any characters remaining
  if rem > 0

    # Append k bits '0', where k is the minimum number >= 0 such that the
    # resulting message length (modulo 512 in bits) is 448.
    if length(msg) > BLOCK_SIZE - 8
      msg = append!(msg, zeros(Uint8, BLOCK_SIZE - rem))
      transform!(state, msg)
      msg = zeros(Uint8, BLOCK_SIZE)
    else
      msg = append!(msg, zeros(Uint8, BLOCK_SIZE - rem))
    end

    # Append length of message (without the '1' bit or padding), in bits, as
    # 64-bit big-endian integer (this will make the entire post-processed
    # length a multiple of 512 bits)
    msg[64] = (bitlen >>> 56) & 0xff
    msg[63] = (bitlen >>> 48) & 0xff
    msg[62] = (bitlen >>> 40) & 0xff
    msg[61] = (bitlen >>> 32) & 0xff
    msg[60] = (bitlen >>> 24) & 0xff
    msg[59] = (bitlen >>> 16) & 0xff
    msg[58] = (bitlen >>> 8) & 0xff
    msg[57] = bitlen & 0xff

    # Process the last block
    transform!(state, msg)
  end

  # Assemble digest and return
  return reduce((x, y) -> string(x, y), map(x -> lpad(hex(ntoh(uint32(x))), 8, "0"), state))
end

end # module RIPEMD


text:
module JetPackWaveFD

using CRC32c, Distributed, DistributedJets, DistributedOperations, Formatting, JetPack, Jets, LinearAlgebra, Statistics, UUIDs, WaveFD

function isvalid_srcfieldfile(m::AbstractArray, hostname, srcfieldfile, chksum)
    _chksum = crc32c(m)
    validhost = hostname == gethostname()
    validfile = isfile(srcfieldfile)
    validmodel = _chksum == chksum
    (validhost && validfile && validmodel, _chksum)
end

function realizewavelet(wavelet, sz, sx, st, dtmod, ntmod_wav)
    local wavelet_realization
    if isa(wavelet, Array)
        if length(sz) == 1 && ndims(wavelet) == 1
            wavelet_realization = reshape(wavelet, length(wavelet), 1)
        else
            wavelet_realization = wavelet
        end
        if size(wavelet_realization) != (ntmod_wav, length(sx))
            throw(ArgumentError("expected size(wavelet_realization) = $ntmod_wav - got $(size(wavelet_realization)) -- dtmod=$(dtmod)"))
        end
    else
        wavelet_realization = zeros(Float32, ntmod_wav, length(sx))
        for is = 1:length(sx)
            wavelet_realization[:,is] .= get(wavelet, convert(Array{Float64,1}, st[is] .+ dtmod * [0:ntmod_wav-1;]))
        end
    end
    wavelet_realization
end

@inline function megacells_per_second(nz, nx, nt, time)
    nt > 1 ? nz*nx*nt/1_000_000/time : 0.0
end

@inline function megacells_per_second(nz, ny, nx, nt, time)
    nt > 1 ? nz*ny*nx*nt/1_000_000/time : 0.0
end

include("ginsu.jl")
include("illumination.jl")
include("jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl")
include("jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl")
include("jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl")
include("jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl")
include("jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl")
include("jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl")
include("jop_sinc_regular.jl")

export modelindex, srcillum, srcillum!

end


text:
################################################################################
# Functions as exported by Clang but removed parameter types.                  #
# This regular expression should help                                          #
# ([A-z_0-9]*)(\{([A-z,_,0-9])*\})*                                          #
################################################################################
function tox_dns3_new(server_public_key)
    ccall((:tox_dns3_new,libtoxcore),Ptr{Void},(Ptr{Uint8},),server_public_key)
end

function tox_dns3_kill(dns3_object)
    ccall((:tox_dns3_kill,libtoxcore),Void,(Ptr{Void},),dns3_object)
end

function tox_generate_dns3_string(dns3_object,string,string_max_len,request_id,name,name_len)
    ccall((:tox_generate_dns3_string,libtoxcore),Cint,(Ptr{Void},Ptr{Uint8},Uint16,Ptr{Uint32},Ptr{Uint8},Uint8),dns3_object,string,string_max_len,request_id,name,name_len)
end

function tox_decrypt_dns3_TXT(dns3_object,tox_id,id_record,id_record_len,request_id)
    ccall((:tox_decrypt_dns3_TXT,libtoxcore),Cint,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Uint32,Uint32),dns3_object,tox_id,id_record,id_record_len,request_id)
end


text:
ENV["PYTHON"] = ""
import Pkg; Pkg.build("PyCall")

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
import Pkg
import RetroCap
import Test
using UUIDs

Test.@testset "RetroCap.jl" begin
    Test.@testset "Compress submodule" begin
        Test.@testset "compress_main.jl" begin
            include("test_compress_main.jl")
        end
        Test.@testset "recompress.jl" begin
            include("test_recompress.jl")
        end
    end

    Test.@testset "assert.jl" begin
        Test.@test_nowarn RetroCap.always_assert(true)
        Test.@test RetroCap.always_assert(true) isa Nothing
        Test.@test RetroCap.always_assert(true) == nothing
        Test.@test Test.@test_nowarn RetroCap.always_assert(true) isa Nothing
        Test.@test Test.@test_nowarn RetroCap.always_assert(true) == nothing
        Test.@test_throws RetroCap.AlwaysAssertionError RetroCap.always_assert(false)
    end
    Test.@testset "compat_entries.jl" begin
        Test.@testset "_compute_cap_upper_bound" begin
            Test.@test RetroCap._compute_cap_upper_bound(v"1.2.3") == "1"
            Test.@test RetroCap._compute_cap_upper_bound(v"1.2.0") == "1"
            Test.@test RetroCap._compute_cap_upper_bound(v"1.0.3") == "1"
            Test.@test RetroCap._compute_cap_upper_bound(v"1.0.0") == "1"
            Test.@test RetroCap._compute_cap_upper_bound(v"0.2.3") == "0.2"
            Test.@test RetroCap._compute_cap_upper_bound(v"0.2.0") == "0.2"
            Test.@test RetroCap._compute_cap_upper_bound(v"0.0.3") == "0.0.3"
            Test.@test RetroCap._compute_cap_upper_bound(v"0.0.0") == "0.0.0"
        end
        Test.@testset "_extract_lower_bound" begin
            Test.@test RetroCap._extract_lower_bound("1.2.3 - foo") == v"1.2.3"
            Test.@test RetroCap._extract_lower_bound("1.2 - foo") == v"1.2"
            Test.@test RetroCap._extract_lower_bound("1 - foo") == v"1"
            Test.@test RetroCap._extract_lower_bound(" - foo") == v"0"
        end
        Test.@testset "generate_compat_entry" begin
            Test.@test "0.0.0 - 1" == RetroCap.generate_compat_entry(RetroCap.Package("PkgA"),
                                                                     RetroCap.Package("PkgB"),
                                                                     RetroCap.UpperBound(),
                                                                     "0",
                                                                     Pkg.Types.semver_spec("0"),
                                                                     v"1.2.3",
                                                                     nothing)
        end
        Test.@testset "too_many_breaking_nonzero_versions" begin
            Test.@test !RetroCap.too_many_breaking_nonzero_versions(Pkg.Types.semver_spec("0"))
            Test.@test RetroCap.too_many_breaking_nonzero_versions(Pkg.Types.semver_spec(">=1.2.3"))
        end
        Test.@testset "too_many_breaking_zero_versions" begin
            Test.@test RetroCap.too_many_breaking_zero_versions(Pkg.Types.semver_spec("0"))
            Test.@test !RetroCap.too_many_breaking_zero_versions(Pkg.Types.semver_spec(">=1.2.3"))
        end
        Test.@testset "next_breaking_release" begin
            Test.@test v"0.0.0" == RetroCap.next_breaking_release(nothing)
            Test.@test v"0.2.0" == RetroCap.next_breaking_release(v"0.1.0")
            Test.@test v"2.0.0" == RetroCap.next_breaking_release(v"1.0.0")
        end
    end
    Test.@testset "Monotonic" begin
        RetroCap.with_temp_dir() do tmp_dir
            cd(tmp_dir) do
                # Create a fake registry
                mkdir("Fake")
                pkga = joinpath("Fake", "A", "APkg")
                pkgb = joinpath("Fake", "B", "BPkg")
                pkgc = joinpath("Fake", "C", "CPkg")
                uuida, uuidb, uuidc = uuid4(), uuid4(), uuid4()
                uuidr = uuid4()
                open(joinpath("Fake", "Registry.toml"), "w") do io
                    println(io, """
                    name = "Fake"
                    uuid = "$uuidr"
                    repo = "fake"

                    [packages]
                    $uuida = { name = "APkg", path = "A/APkg" }
                    $uuidb = { name = "BPkg", path = "B/BPkg" }
                    $uuidc = { name = "CPkg", path = "C/CPkg" }
                    """)
                end
                map(mkpath, [pkga, pkgb, pkgc])
                open(joinpath(pkga, "Versions.toml"), "w") do io
                    println(io, """
                    ["1.0.0"]
                    git-tree-sha1 = "aa"

                    ["2.0.0"]
                    git-tree-sha1 = "bb"

                    ["3.0.0"]
                    git-tree-sha1 = "cc"
                    """)
                end
                open(joinpath(pkgb, "Versions.toml"), "w") do io
                    println(io, """
                    ["0.1.0"]
                    git-tree-sha1 = "aa"

                    ["0.2.0"]
                    git-tree-sha1 = "bb"
                    """)
                end
                open(joinpath(pkgc, "Deps.toml"), "w") do io
                    println(io, """
                    [0]
                    APkg = "$uuida"

                    ["0.1-0.2"]
                    BPkg = "$uuidb"
                    """)
                end
                open(joinpath(pkgc, "Versions.toml"), "w") do io
                    println(io, """
                    ["0.1.0"]
                    git-tree-sha1 = "aa"

                    ["0.2.0"]
                    git-tree-sha1 = "bb"

                    ["0.3.0"]
                    git-tree-sha1 = "cc"
                    """)
                end
                open(joinpath(pkgc, "Compat.toml"), "w") do io
                    println(io, """
                    ["0.1.0"]
                    APkg = "1-*"
                    BPkg = "0.1-0"

                    ["0.2.0"]
                    APkg = "2"
                    BPkg = "0.1-0.2"

                    ["0.3.0"]
                    APkg = "3"
                    """)
                end
                RetroCap.add_caps(RetroCap.MonotonicUpperBound(), RetroCap.CapLatestVersion(), "Fake")
                str = read(joinpath(pkgc, "Compat.toml"), String)
                Test.@test str == """
                ["0-0.1"]
                APkg = "1.0.0 - 2"
                BPkg = "0.1.0 - 0.2"

                ["0.2"]
                APkg = "2"
                BPkg = "0.1-0.2"

                ["0.3-0"]
                APkg = "3"
                """
            end
        end
    end
    Test.@testset "Run RetroCap on the General registry" begin
        RetroCap.with_temp_dir() do tmp_dir
            cd(tmp_dir)

            run(`git clone https://github.com/JuliaRegistries/General.git General`)
            run(`git clone https://github.com/BioJulia/BioJuliaRegistry.git BioJuliaRegistry`)

            RetroCap.add_caps(RetroCap.UpperBound(), RetroCap.ExcludeLatestVersion(), "General")
            RetroCap.add_caps(RetroCap.UpperBound(), RetroCap.ExcludeLatestVersion(), Any["General", "BioJuliaRegistry"])

            RetroCap.add_caps(RetroCap.UpperBound(), RetroCap.CapLatestVersion(), "General")
            RetroCap.add_caps(RetroCap.UpperBound(), RetroCap.CapLatestVersion(), Any["General", "BioJuliaRegistry"])

            RetroCap.add_caps(RetroCap.MonotonicUpperBound(), RetroCap.CapLatestVersion(), "General")
            RetroCap.add_caps(RetroCap.MonotonicUpperBound(), RetroCap.CapLatestVersion(), Any["General", "BioJuliaRegistry"])
        end
    end
end


text:
# main file for evaluating the residual

import PDESolver: updateMetricDependents, evalResidual

"""
### EulerEquationMod.init

  This function performs any operations that need to be performed before
  the first residual evaluation.
  Any operations that need to be performed at the beginning of *each*
  residual evaluation should go in dataPrep

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * pmesh (equal to `mesh` by default)

  **Keyword Arguments**

   * init_mesh: set the boundary condition functors, default true,
                should false if `eqn`is nested inside another physics
                module's [`AbstractSolutionData`](@ref) object.
"""
# high level functions
function init(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
              eqn::NSData{Tsol, Tres}, opts, pmesh=mesh;
              init_mesh=true) where {Tmsh, Tsol, Tres}

  # initialize the Euler equation object.  Tell it not to use a source term,
  # because the Navier Stokes module will handle the it
  src_orig = opts["SRCname"]
  opts["SRCname"] = "SRC0"
  # see `getBCFunctors` in this module for why init_mesh=true
  EulerEquationMod.init(mesh, sbp, eqn.euler_eqn, opts, init_mesh=true)
  opts["SRCname"] = src_orig

  if init_mesh
#    getBCFunctors(mesh, sbp, eqn, opts)
#    getBCFunctors(pmesh, sbp, eqn, opts)
  end

  getSRCFunctors(mesh, sbp, eqn, opts)
  if mesh.isDG
    getFluxFunctors(mesh, sbp, eqn, opts)
  end

  return nothing
end

function evalResidual(mesh::AbstractMesh, sbp::AbstractOperator, eqn::NSData,
                     opts::Dict, t=0.0)

  @assert eqn.commsize == 1

  params = eqn.params
  time = params.time
  time.t_dataprep += @elapsed dataPrep(mesh, sbp, eqn, opts)

  time.t_volume += @elapsed if opts["addVolumeIntegrals"]
    evalVolumeIntegrals(mesh, sbp, eqn, opts)
#    println("volume integral @time printed above")
  end

  time.t_bndry += @elapsed if opts["addBoundaryIntegrals"]
    # do in inviscid-type boundary integral
    # Because we only support non-precompute, all this function does is
    # the integration, so we can re-use the inviscid one
    EulerEquationMod.evalBoundaryIntegrals(mesh, sbp, eqn.euler_eqn, opts)
#   println("boundary integral @time printed above")
  end

  time.t_face += @elapsed if mesh.isDG && opts["addFaceIntegrals"]
    # invisicd face integrals
    EulerEquationMod.evalFaceIntegrals(mesh, sbp, eqn.euler_eqn, opts)
#    println("face integral @time printed above")
  end

  time.t_source += @elapsed evalSourceTerm(mesh, sbp, eqn, opts)
#  println("source integral @time printed above")

  if eqn.params.isViscous == true
    time.t_face += @elapsed evalFaceIntegrals_vector(mesh, sbp, eqn, opts)
    # do the non-inviscid-type boundary integral
    time.t_bndry += @elapsed evalBoundaryIntegrals_vector(mesh, sbp, eqn, opts)
  end

  # apply inverse mass matrix to eqn.res, necessary for CN
  if opts["use_Minv"]
    EulerEquationMod.applyMassMatrixInverse3D(mesh, sbp, eqn.euler_eqn, opts, eqn.res)
  end

  return nothing
end

"""
### NavierStokesMod.dataPrep

  This function calculates all the quantities that are stored over the entire
  mesh, and performs a few sanity check such as ensuring density and pressure
  are positive.

  This function also calls `dataPrep` for Euler.

  This is a high level function.
"""
function dataPrep(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::NSData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  # do the Euler part
  EulerEquationMod.dataPrep(mesh, sbp, eqn.euler_eqn, opts)

  # do the viscous part
  if eqn.params.isViscous == true
    # fill!(eqn.vecflux_face, 0.0)
    fill!(eqn.vecflux_faceL, 0.0)
    fill!(eqn.vecflux_faceR, 0.0)

    calcViscousFlux_interior(mesh, sbp, eqn, opts)

    fill!(eqn.vecflux_bndry, 0.0)
    calcViscousFlux_boundary(mesh, sbp, eqn, opts)
  end

  return nothing
end

"""
  Does the volume integrals

  Calls the Euler module do to the first part and then differentiates the
  result again.
"""
function evalVolumeIntegrals(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator, eqn::NSData{Tsol, Tres, Tdim},
                             opts) where {Tmsh,  Tsol, Tres, Tdim}

  EulerEquationMod.evalVolumeIntegrals(mesh, sbp, eqn.euler_eqn, opts)

  # now differentiate the volume flux agains for the viscous part
  if eqn.params.isViscous == true
    weakdifferentiate2!(mesh, sbp, eqn, eqn.res)
  end

  return nothing
end

"""
  Updates the residual with the source term contribution.

  Source terms are defined in this module (not Euler), although currently
  the source term functors are compatible.

  **Inputs**:

   * mesh : Abstract mesh type
   * sbp  : Summation-by-parts operator
   * eqn  : Euler equation object
   * opts : options dictonary

  Outputs: none

  Aliasing restrictions: none

"""
function evalSourceTerm(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::NSData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  if opts["use_src_term"]
    # the source term signature are compatible, so use the Euler function
    # to evaluate it
    EulerEquationMod.applySourceTerm(mesh, sbp, eqn, opts, eqn.src_func)
  end

  return nothing
end  # end function

function updateMetricDependents(mesh::AbstractMesh, sbp::AbstractOperator,
                                eqn::NSData, opts)

  # all the arrays are aliased, so use the Euler version
  EulerEquationMod.updateMetricDependents(mesh, sbp, eqn.euler_eqn, opts)

  return nothing
end

function majorIterationCallback(itr::Integer,
                                mesh::AbstractMesh,
                                sbp::AbstractOperator,
                                eqn::NSData, opts, f::IO)

  # the Euler callback does useful things and doesn't have anything incorrect
  # for viscous, so use it.  Eventually Navier Stokes should have its own
  EulerEquationMod.majorIterationCallback(itr, mesh, sbp, eqn.euler_eqn, opts, f)

  return nothing
end


text:
# Alignment Anchor
# ================
#
# Sequence alignment anchor type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Alignment operation with anchoring positions.
"""
struct AlignmentAnchor
    seqpos::Int
    refpos::Int
    op::Operation
end

function AlignmentAnchor(pos::Tuple{Int,Int}, op)
    return AlignmentAnchor(pos[1], pos[2], op)
end

function Base.show(io::IO, anc::AlignmentAnchor)
    print(io, "AlignmentAnchor(", anc.seqpos, ", ", anc.refpos, ", '", anc.op, "')")
end


text:
using Elefridge
using ZfpCompression
using PyPlot
using PyCall
using Blosc
xr = pyimport("xarray")
ccrs = pyimport("cartopy.crs")

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,"_go3.grib"),readdir(path))
grib = xr.open_dataarray(joinpath(path,filelist[end]),engine="cfgrib")
X = grib.data
lat = grib.latitude.data
lon = grib.longitude.data

o3clin24 = Array(LinQuant24Array(X))
o3clog16 = Array(LogQuant16Array(X))
o3clin16 = Array(LinQuant16Array(X))
o3clog8 = Array(LogQuant8Array(X))
o3clin8 = Array(LinQuant8Array(X))

for level in 10:10:137

    ## compression
    rbits = [10,7,5,3,1]
    cfs = [32/24,32/16,32/16,32/8,32/8]

    lat_div(n::Integer) = Array(-90:180/(n-1):90)
    lon_div(n::Integer) = Array(0:360/n:360)

    ## PLOT
    ioff()
    o3 = X[level,:,:]'
    vmin,vmax = extrema(o3)
    fig,ax1 = plt.subplots(1,1,figsize=(10,6),sharex=true,sharey=true,
                        subplot_kw=Dict("projection"=>ccrs.Robinson(central_longitude=180)))

    lond = lon_div(length(rbits)+1)

    # Float32
    ax1.pcolormesh(lon[lon .<= lond[2]],lat,o3[lon .<= lond[2],:]',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree())

    for (i,o3c) in zip(3:length(lond),[o3clin24,o3clog16,o3clin16,o3clog8,o3clin8])
        lonwhere = (lon .<= lond[i]) .& (lon .> lond[i-1]-0.5)

        o3cl = o3c[level,:,:]'

        ax1.pcolormesh(lon[lonwhere],lat,o3cl[lonwhere,:]',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree())
    end

    ax1.coastlines(color="k",lw=0.5)

    for longi in lond[2:end-1]
        ax1.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
        # ax1.plot(repeat([longi,],6),lat_div(31)[end-5:end],"w",transform=ccrs.PlateCarree(),lw=1.5)
    end

    ax1.text(0.127,0.105,"Compression\nfactor",transform=ax1.transAxes,fontweight="bold",color="white")

    for (i,x) in enumerate([0.23,0.33,0.43,0.54,0.64,0.76])
        c = i == 1 ? 1 : Int(round(cfs[i-1]))
        c = i == 2 ? "1.3" : c
        ax1.text(x,0.03,"$(c)x",transform=ax1.transAxes,fontweight="bold",fontsize=14,
                                color="white")
    end

    for (x,txt) in zip([0.24,0.34,0.43,0.53,0.62,0.72],
            ["32-bit\nFloat","24-bit\nLin","16-bit\nLog","16-bit\nLin",
            "8-bit\nLog","8-bit\nLin"])
        ax1.text(x,-0.07,txt,transform=ax1.transAxes)
    end

    ax1.set_title(L"O$_3$ compression: Lin/LogQuantization",fontweight="bold")

    tight_layout()
    savefig("/Users/milan/git/Elefridge.jl/maps/o3/linlog_o3_$level.png",dpi=200)
    close(fig)
end


text:
using KrylovMethods
using Base.Test
using MatrixDepot

println("=== Testing Lanczos helpers ===")
symOrth = KrylovMethods.symOrtho

v     = ([0;randn()],[randn();0],randn(2))
for k=1:3
	println("Givens rotation for v=[$(v[k][1]),$(v[k][2])]")
	c,s,r = symOrth(v[k][1],v[k][2])
	@test abs(c.^2 + s.^2 - 1) < 1e-14
	Q = [c s; -s c]
	w = [norm(v[k]);0]
	@test norm(Q*v[k]-w) < 1e-13
end

println("=== Testing Lanczos Bidiagonalization ===")
A = sprandn(10,10,.4) + speye(10)
p = randn(10)
k = 9

U,B,V = lanczosBidiag(A,p,k)

@test norm(full(U'*A*V-B))/norm(full(A)) < 1e-15

println("=== Testing Lanczos Tridiagonalization ===")

# test stopping for low-rank matrix
A = randn(100,10)
H = A*A'
T,V = lanczosTridiag(H,ones(100),100,doReorth=true)
@test size(T,1)==11
@test size(V,2)==11

A   = matrixdepot("rosser",4)
Aff  = x -> A*x
T1,V1 = lanczosTridiag(A,ones(size(A,1)),size(A,1),tol=1e-10,doReorth=true)
T2,V2 = lanczosTridiag(Aff,ones(size(A,1)),size(A,1),tol=1e-10,doReorth=false)
sT1  = svd(full(T1))[2]
sT2  = svd(full(T2))[2]
sA   = svd(full(A))[2]

@test abs(sT1[1]-sA[1])/sA[1] < 1e-4
@test abs(sT2[1]-sA[1])/sA[1] < 1e-4
B = V1'*V1
@test norm(B - diagm(diag(B))) < 1e-10
@test norm(diag(B)- ones(size(B,1))) < 1e-10
println("=== Lanczos: All tests passed. ===")


text:
export rotate_mirror

const dir_left_rotate = Dict(
    left => down,
    down => right,
    right => up,
    up => left
)

const dir_right_rotate = Dict(
    left => up,
    up => right,
    right => down,
    down => left
)

const dir_180_rotate = Dict(
    left => right,
    right => left,
    up => down,
    down => up
)

const dir_mirror = Dict(
    left => up,
    up=>left,
    down=>right,
    right=>down
)

function rotate_board_dir(board, dir)
    [
        (board, dir),
        (rotl90(board), dir_left_rotate[dir]),
        (rot180(board), dir_180_rotate[dir]),
        (rotr90(board), dir_right_rotate[dir])
    ]
end

function rotate_mirror(board, dir)
    vcat(
        rotate_board_dir(board, dir),
        rotate_board_dir(board |> transpose |> collect, dir_mirror[dir])
    )
end

function rotate_mirror(board)
    tboard = transpose(board) |> collect

    return board,  rotl90(board),  rot180(board),  rotr90(board),
           tboard, rotl90(tboard), rot180(tboard), rotr90(tboard)
end


text:
struct ConfusionMatrix{T<:Integer}
    c::Array{T,2}
end

size(c::ConfusionMatrix,kws...)=size(c.c,kws...)
sum(c::ConfusionMatrix,kw...)=sum(c.c,kw...)
+(a::ConfusionMatrix,b::ConfusionMatrix)=ConfusionMatrix(a.c.+b.c)
eltype(a::ConfusionMatrix)=eltype(a.c)

function confusion(labels, predicts, classes)
    predict_real = zeros(Integer, classes, classes)
    for (p, r) in zip(predicts, labels)
        predict_real[p + 1,r + 1] += 1
    end
    return ConfusionMatrix(predict_real)
end

tp(c::ConfusionMatrix,t) = c.c[t,t]
fp(c::ConfusionMatrix,t) = sum(c.c[t,1:t - 1]) + sum(c.c[t,t + 1:end])
tn(c::ConfusionMatrix,t) = sum(c.c[1:t - 1,1:t - 1]) + sum(c.c[1:t - 1,t + 1:end]) + sum(c.c[t + 1:end,1:t - 1]) + sum(c.c[1 + t:end,1 + t:end])
fn(c::ConfusionMatrix,t) = sum(c.c[1:t - 1,t]) + sum(c.c[t + 1:end,t])

function tpr(c::ConfusionMatrix, t=0)
    if t == 0
        cs = size(c, 1)
        TP = [tp(c, i) for i in 1:cs]
        FN = [fp(c, i) for i in 1:cs]
    else
        TP = tp(c, t)
        FN = fn(c, t)
    end
    return TP ./ (TP .+ FN)
end

function fpr(c::ConfusionMatrix, t=0)
    if t == 0
        cs = size(c, 1)
        FP = [fp(c, i) for i in 1:cs]
        TN = [tn(c, i) for i in 1:cs]
    else
        FP = fp(c, t)
        TN = tn(c, t)
    end
    return FP ./ (FP .+ TN)
end

recall = tpr
sensitivity = tpr

function precision(c::ConfusionMatrix, t=0)
    if t == 0
        cs = size(c, 1)
        TP = [tp(c, i) for i in 1:cs]
        FP = [fp(c, i) for i in 1:cs]
    else
        TP = tp(c, t)
        FP = fp(c, t)
    end
    return TP ./ (TP .+ FP)
end

function accuracy(c::ConfusionMatrix, t=0)
    if t == 0
        cs = size(c, 1)
        TP = [tp(c, i) for i in 1:cs]
        TN = [tn(c, i) for i in 1:cs]
    else
        TP = tp(c, t)
        TN = tn(c, t)
    end
    return (TP .+ TN) ./ sum(c)
end

function threshold_scan(labels, scores)
    len=length(labels)
    p = sortperm(scores)
    c = zeros(Int64, 2, 2)
    tn, fp, fn, tp = CartesianIndices(c)
    T = sum(labels)
    N = length(labels) - T
    c[tp] = T
    c[fp] = N
    result=Array{ConfusionMatrix,1}(undef,len-1)
    for (j,i) in enumerate(p)
        if labels[i] == true
            c[fn] += 1
            c[tp] -= 1
        else
            c[tn] += 1
            c[fp] -= 1
        end
        j == len && break
        result[j]=ConfusionMatrix(copy(c))
    end
    result
end

function auc(x::Array{<:Real}, y::Array{<:Real})
    n = length(x)
    length(y) ≠ n && throw(DimensionMismatch("`x` and `y` must have same length (`x` is $n, `y` is $(length(y)))."))
    issorted(x) || issorted(x, rev=true) || throw(ArgumentError("`x` must be sorted in increasing order."))
    rev = zero(promote_type(eltype(x), eltype(y)))
    p = sortperm(x)
    x₀=0
    y₀=y[p[1]]
    for i = 1:n
        Δx = x[p[i]] - x₀
        ȳ = (y[p[i]] + y₀) / 2
        if !isnan(Δx) && !isnan(ȳ)
            rev += Δx * ȳ
        end
        x₀=x[p[i]]
        y₀=y[p[i]]
    end
    rev
end

export threshold_scan,auc,tpr,fpr,recall,precision,confusion,ConfusionMatrix,accuracy

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

module Descriptors

using JuLIP

using ACE: SparseSHIP, PolyTransform, PolyCutoff1s, SHIPBasis

export SHIPDescriptor, descriptors

"""
`SHIPDescriptor(; deg=nothing, wY=1.5, rcut=nothing, r0=1.0, p=2)`

This returns a `desc::SHIPBasis` object which can be interpreted as a descriptor
map. Call `descriptors(desc, py_at)` where `py_at` is a Python object
(of ase Atoms type) to get the descriptors.

Call `Dict(desc)` to obtain a dictionary that fully describes the descriptor
and can be serialised to a JSON file / deserialised.

* `bodyorder` : specify body-order; `bodyorder = 3` corresponds to SOAP??;
`bodyorder = 5` to SNAP??;
* `deg, wY` : specify polynomial degree restriction; the basis will contain all
tensor products `Pk * Ylm` such that `k + wY * l ≦ deg`.
* `rcut` : cutoff radius
* `r0` : an estimate for nearest-neighbour distance - not crucial
* `p` : specifies distance transform, u = (r/r0)^(-p); i.e., polynomials
`Pk` are polynomials in `u` not in `r`. (e.g. p = 1 => Coulomb coordinates)
"""
function SHIPDescriptor(species = :X;
                        bodyorder=3, deg=nothing,
                        rcut=nothing, r0=2.5, p=2,
                        kwargs... )
   spec = SparseSHIP(species, bodyorder-1, deg; kwargs...)
   trans = PolyTransform(p, r0)
   fcut = PolyCutoff1s(2, rcut)
   return SHIPBasis(spec, trans, fcut)
end

function descriptors(basis::SHIPBasis, _at)
   at = Atoms(_at)
   B = zeros(Float64, length(basis), length(at))
   for i = 1:length(at)
      B[:, i] = site_energy(basis, at, i)
   end
   return B
end

end


text:
import GeoStatsBase: solve_single

# define a dummy solver for testing
@simsolver Dummy begin end
function solve_single(problem::SimulationProblem,
                      var::Symbol, solver::Dummy, preproc)
  npts = npoints(domain(problem))
  V = variables(problem)[var]
  vcat(fill(zero(V), npts÷2), fill(one(V), npts÷2))
end


text:
module MixedModelsPermutations

using BlockDiagonals
using LinearAlgebra
using MixedModels # we add several methods
using Random
using SparseArrays
using StaticArrays
using Statistics
using StatsBase
using StatsModels
using Tables

using MixedModels: MixedModelBootstrap
using MixedModels: fixef!, stderror!
using MixedModels: getθ!, updateL! # setθ! is imported for extending
using MixedModels: replicate

import MixedModels: allpars, coefpvalues, issingular, setθ!, tidyβ, tidyσs, residuals

export MixedModelPermutation,
       nonparametricbootstrap,
       olsranef,
       permutation,
       permutationtest,
       residuals

include("mixedmodelpermutation.jl")
include("inflation.jl")
include("nonparametricbootstrap.jl")
include("permutation.jl")
include("ols.jl")

end


text:
#! /usr/bin/julia

# Rosetta Code, Pernicious numbers

ispernicious(n::Int) = isprime(count_ones(n))

pcnt = 0
i = 0
print("    ")
while pcnt < 25
    i += 1
    ispernicious(i) || continue
    pcnt += 1
    print(i, " ")
end
println()

print("    ")
for i in 888888877:888888888
    ispernicious(i) || continue
    print(i, " ")
end
println()

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run subset of tests
=#

using Test
using Printf
using Hypatia

test_files = [
    "polyutils",
    "cone",
    "native",
    "moi",
    "examples",
    # require optional dependencies:
    # "pardiso",
    # "hsl",
    ]

println()
@info("starting all tests")
println()
timings = Dict{String, Float64}()

@testset "all tests" begin
all_test_time = @elapsed for t in test_files
    @info("starting $t tests")
    test_time = @elapsed include("run$(t)tests.jl")
    flush(stdout); flush(stderr)
    @info("finished $t tests in $(@sprintf("%8.2e seconds", test_time))")
    println()
    timings[t] = test_time
end

@info("finished all tests in $(@sprintf("%8.2e seconds", all_test_time))")
println("\ntest suite timings (seconds):")
display(timings)
println()
end

println()
;


text:
module YaoPlots

export plot

using Compose

plot(;kwargs...) = x->plot(x;kwargs...)

include("helperblock.jl")
include("vizcircuit.jl")
include("zx_plot.jl")

end


text:
using Base.Test
using Distributions
# using Gadfly
using Klara

# Target distribution shape, unscaled
g(x) = exp(-dot(x, x)/2)
log_g(x) = -dot(x, x)/2

# Candidate (starting) density, use Normal(0, 2)
g0(x) = 1/(sqrt(2*pi)*2)*exp(-dot(x, x)/8)
log_g0(x) = -dot(x, x)/8-log(2 * sqrt(2*pi))

# Compute M0 such that candidate dominates target
h0(x) = -g(x)/g0(x)

#res0 = optimize(h0, ones(L))
#M0 = -res0.f_minimum

M0 = 10.0
log_M0 = log(M0)

# Candidate (starting) Distribution, use Normal(0, 2)
g1(x) = exp(sum(logpdf(Normal(0, 2), x)))
log_g1(x) = sum(logpdf(Normal(0, 2), x))

# Compute M1 such that candidate dominates target
h1(x) = -g(x)/sum(g1(x))
#res1 = optimize(h1, ones(L))
#M1 = -res1.f_minimum
M1 = 10.0
log_M1 = log(M1)

model = likelihood_model(BasicContUnvParameter(:p, logtarget=log_g), false)

sampler = ARS(log_g0, proposalscale=log_M0)

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

v0 = Dict(:p=>1.)

job = BasicMCJob(model, sampler, mcrange, v0)

run(job)

chain1 = output(job)

sampler = ARS(log_g1, proposalscale=log_M1)

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

v0 = Dict(:p=>1.)

job = BasicMCJob(model, sampler, mcrange, v0)

run(job)

chain2= output(job)

println()
println("M: ", [M0 M1])
println()

println("Acceptance rate: ", acceptance(chain1, diagnostics=false))

println("Acceptance rate: ", acceptance(chain2, diagnostics=false))

theta = hcat(chain1.value, chain2.value);
for i in sample(1:length(theta), 200)
  @test isapprox(h0(theta[i])/M0, h1(theta[i])/M1)
end

# plot(x=chain1.value, Geom.histogram(bincount=20))

# plot(x=chain2.value, Geom.histogram(bincount=20))


text:
# PLOTS_DEFAULTS = Dict(:markersize => 7, :linewidth => 3, :legend => true, :warn_on_unsupported => true) # for presentation
fontsize = 16;
PLOTS_DEFAULTS = Dict(
     :legend => true,
     :warn_on_unsupported => true,
     :legendfontsize => fontsize,
     :xguidefontsize => fontsize,
     :yguidefontsize => fontsize,
     :guidefontsize => fontsize,
     :linewidth => 1,
     :markersize => 8,
     # :fontfamily => "serif",
     # :fontfamily => "IPAexGothic",
     :tickfontsize => fontsize
     )

max_stars_repo_name: uchiiii/dotfiles

max_issues_repo_name: uchiiii/dotfiles

max_forks_repo_name: uchiiii/dotfiles


text:
struct Bucket{K, V}
    items::Vector{Pair{K, V}}
end

function Base.iterate(b::Bucket)
    if isempty(b.items)
        return nothing
    end
    return (b.items[1], 2)
end

function Base.iterate(b::Bucket, state::Int)
    if length(b.items) < state
        return nothing
    end
    return (b.items[state], state + 1)
end

"""
    Base.setindex!(b::Bucket{K, V}, value::V, key::K) where {K, V}

Set `value` for `key`, returns a tuple `(value, replaced::Boolean)`.
"""
function Base.setindex!(b::Bucket{K, V}, value::V, key::K) where {K, V}
    for i in eachindex(b.items)
        (k, v) = b.items[i]
        if k == key
            b.items[i] = key => value
            return (value, true)
        end
    end
    push!(b.items, key => value)
    return (value, false)
end

function Base.getindex(b::Bucket{K, V}, key::K) where {K, V}
    for (k, v) in b.items
        if k == key
            return v
        end
    end
    throw(KeyError("Key $key not in collection."))
end

function Base.get(b::Bucket{K, V}, key::K, default::V) where {K, V}
    for (k, v) in b.items
        if k == key
            return v
        end
    end
    return default
end

function Base.haskey(b::Bucket{K, V}, key::K) where {K, V}
    for (k, v) in b.items
        if k == key
            return true
        end
    end
    return false
end

function Base.delete!(b::Bucket{K, V}, key::K) where {K, V}
    for (i, pair) in enumerate(b.items)
        (k, v) = pair
        if k == key
            deleteat!(b.items, i)
            return true
        end
    end
    return false
end

function Base.empty!(b::Bucket)
    empty!(b.items)
    return nothing
end


text:
module StructuredProximalOperators

using Parameters
import LinearAlgebra: norm, dot
import Manifolds:
    check_manifold_point,
    check_tangent_vector,
    distance,
    embed,
    embed!,
    exp!,
    inner,
    log!,
    manifold_dimension,
    project,
    project!,
    representation_size,
    retract,
    retract!,
    zero_tangent_vector
using LinearAlgebra
using SparseArrays

using ManifoldsBase
using Manifolds

using JuMP

using Random
using PolynomialRoots

using PGFPlotsX, LaTeXStrings, ForwardDiff

import Manifolds: show
import Base: show, ==, <, copy

## Manifolds exports
export Manifold
export SVDMPoint, UMVTVector
export Euclidean, ℝ
export manifold_dimension
export randomMPoint, randomTVector
export inner, norm
export project, project!
export retract, retract!
export egrad_to_rgrad!, egrad_to_rgrad, ehess_to_rhess!, ehess_to_rhess
export check_manifold_point, check_tangent_vector, show

export is_manifold_point, is_tangent_vector
export representation_size
export embed!, embed, zero_tangent_vector
export embedding_dimension

## Regularizers exports
export g, prox_αg, prox_αg!, ∇M_g, ∇M_g!, ∇²M_g_ξ, ∇²M_g_ξ!
export model_g_subgradient!, build_subgradient_from_normalcomp, build_normalcomp_from_subgradient

export l1Manifold, regularizer_l1
export FixedRankMatrices, regularizer_lnuclear
export PSphere, regularizer_distball
export ProductManifold, regularizer_group
export PlaneParabola, maxquadAL

export wholespace_manifold

export copy
export Regularizer

abstract type Regularizer end

<(::Euclidean, ::Euclidean) = false

function copy(M::Euclidean{N,𝔽}) where {N,𝔽}
    return Euclidean(representation_size(M)...; field = 𝔽)
end

function copy(M::ProductManifold{𝔽,TM}) where {𝔽,TM}
    return ProductManifold(M.manifolds...)
end

function name(M::Euclidean{dim,𝔽}; short = true) where {dim,𝔽}
    rep = representation_size(M)
    res = "$𝔽 ^" * string(rep)
    if length(rep) == 1
        res = "$𝔽 ^" * string(rep[1])
    end
    return res
end

## Helper
softthresh(x, α) = sign(x) * max(0, abs(x) - α)

## Interface

function g(::Regularizer, x)
    return error("g not implemented for regularizer $(typeof(g)) and point $(typeof(x)).")
end

function prox_αg(reg::T, x, α) where {T<:Regularizer}
    res = zero(x)
    M = prox_αg!(reg, res, x, α)
    return res, M
end
function prox_αg!(g, res, x, α)
    return error("prox_αg! not implemented for regularizer $(typeof(g)), point $(typeof(x)).")
end

function ∇M_g(g, M, x)
    return error("∇M_g not implemented for regularizer $(typeof(g)), manifold $M, point $(typeof(x)).")
end
function ∇M_g!(g, M, res, x)
    return error("∇M_g! not implemented for regularizer $(typeof(g)), manifold $M, result $(typeof(res)), point $(typeof(x)).")
end

function ∇²M_g_ξ(g, M, x, ξ)
    return error("∇²M_g_ξ not implemented for regularizer $(typeof(g)), manifold $M, point $(typeof(x)), vector $(typeof(ξ)).")
end
function ∇²M_g_ξ!(g, M, res, x, ξ)
    return error("∇²M_g_ξ not implemented for regularizer $(typeof(g)), manifold $M, result $(typeof(res)), point $(typeof(x)), vector $(typeof(ξ)).")
end

##
egrad_to_rgrad!(M::Manifold, gradf_x, x, ∇f_x) = project!(M, gradf_x, x, ∇f_x)
function egrad_to_rgrad(M::Manifold, x, ∇f_x)
    gradf = zero_tangent_vector(M, x)
    return project!(M, gradf, x, ∇f_x)
end
function ehess_to_rhess(M::Manifold, x, ∇f_x, ∇²f_ξ, ξ)
    Hessf_xξ = Manifolds.allocate(ξ)
    return ehess_to_rhess!(M, Hessf_xξ, x, ∇f_x, ∇²f_ξ, ξ)
end

function wholespace_manifold(::Tr, x) where {Tr<:Regularizer}
    return Euclidean(size(x)...)
end

## Manifolds
include("manifolds/l1subspace.jl")
include("manifolds/FixedRankMatrices.jl")
include("manifolds/PShpere.jl")
include("manifolds/planeParabola.jl")

## Regularizers
include("regularizers/regularizer_l1.jl")
include("regularizers/regularizer_lnuclear.jl")
include("regularizers/regularizer_distball.jl")
include("regularizers/regularizer_group.jl")
include("regularizers/maxquadAL.jl")

include("compare_smoothcurves.jl")

export compare_curves,
    display_curvescomparison,
    check_e2r_gradient_hessian,
    check_retraction,
    check_regularizer_gradient_hessian

end


text:
function update_psi!(state::State, data::Data, prior::Prior)
  m, s = params(prior.psi)
  vkernel = zero.(state.psi)
  mkernel = zero.(state.psi)

  # Update kernels
  for i in samplenames
    yfinite = ref_yfinite(data, i)
    lam = ref_lambda(state, i)
    zeta = ref_zeta(state, i)
    v = ref_v(state, i)
    for n in eachindex(v)
      k = lam[n]
      vkernel[k] += zeta[n]^2 * v[n] / state.omega[k]
      mkernel[k] += zeta[n] * (yfinite[n] - state.mu[k]) * v[n] / state.omega[k]
    end
  end

  vnew = 1 ./ (s^-2 .+ vkernel)
  mnew = vnew .* (m/s^2 .+ mkernel)

  state.psi .= randn(prior.K) .* sqrt.(vnew) + mnew

  update_sigma!(state)
  update_phi!(state)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
@testset "basic integration" begin
    m = SquareRooter(4)
    report = IC.train!(m, Step(2), InvalidValue(), NumberLimit(3); verbosity=0);
    @test report[1] == (Step(2), (new_iterations = 6,))
    @test report[2] == (InvalidValue(), (done=false, log=""))
    report[3] == (NumberLimit(3),
                  (done=true,
                   log="Stop triggered by NumberLimit(3) "*
                   "stopping criterion. "))

    m = SquareRooter(4)
    @test_logs((:info, r"final loss"),
               (:info, r"final training loss"),
               (:info, r"Stop triggered by Num"),
               IC.train!(m, Step(2),  InvalidValue(), NumberLimit(3)));
    @test_logs((:info, r"Using these controls"),
               (:info, r"Stepping model for 2 more iterations"),
               (:info, r"Stepping model for 2 more iterations"),
               (:info, r"Stepping model for 2 more iterations"),
               (:info, r"final loss"),
               (:info, r"final training loss"),
               (:info, r"total control cycles\: 3"),
               (:info, r"A total of 6 iterations added"),
               (:info, r"Stop triggered by NumberLimit"),
               IC.train!(m, Step(2),
                         InvalidValue(),
                         NumberLimit(3);
                         verbosity=2));
end

struct WrappedSquareRooter
    s::SquareRooter
end

@testset "expose and Callback(..., raw=true)" begin
    IC.train!(w::WrappedSquareRooter, n) = IC.train!(w.s, n)
    IC.loss(w::WrappedSquareRooter) = IC.loss(w.s)
    IC.train!(w::WrappedSquareRooter, n) = IC.train!(w.s, n)
    IC.expose(w::WrappedSquareRooter) = w.s

    wrapper = WrappedSquareRooter(SquareRooter(4))

    roots = Any[]
    wrappers = Any[]

    g(w) = push!(wrappers, deepcopy(w))

    c1 = Callback(s->push!(roots, s.root))
    c2 = Callback(g, raw=true)

    IC.train!(wrapper, Step(1), c1, c2, NumberLimit(2), verbosity=0)

    @test map(w -> w.s.root, wrappers) ≈ roots
end

@testset "skip integration" begin
    m = SquareRooter(4)
    numbers = Int[]
    IC.train!(m,
              Step(1),
              IterationControl.skip(
                  WithNumberDo(x->push!(numbers, x)), predicate=3),
              NumberLimit(10), verbosity=0)
    @test numbers == [3, 6, 9]
end

@testset "integration test related to #38" begin
    model = IterationControl.SquareRooter(4)
    @test_logs((:info, r"number"),
               (:info, r"number"),
               (:info, r"final loss"),
               (:info, r"final training loss"),
               (:info, r"Stop triggered by"),
               IC.train!(model,
                         Step(1),
                         Threshold(2.1),
                         WithNumberDo(),
                         IterationControl.skip(WithLossDo(), predicate=3)))
end

@testset "integration test with PQ" begin
    # get complete training losses:
    model = SquareRooter(1e40)
    IC.train!(model, 70)
    tlosses = IC.training_losses(model)

    # get complete losses:
    model = SquareRooter(1e40)
    losses = map(1:70) do _
        IC.train!(model, 1)
        IC.loss(model)
    end

    # train with PQ criterion:
    K = 5
    model = SquareRooter(1e40)
    pq = PQ(k=K)
    states = Any[]
    f(s) = push!(states, s)
    controls = [Step(1), NumberLimit(70), IC.with_state_do(pq, f=f)]
    IC.train!(model, controls..., verbosity=0)

    # check internal state of PQ makes sense, as far as `loss` and
    # `training_losses` are concerned:
    @test map(s->s.loss, states) == losses
    @test map(enumerate(states)) do (j, s)
        L = min(K, j) # what length of PQ training_losses should be
        t= s.loss == losses[j] &&
            s.training_losses == reverse(tlosses[j-L+1:j])
    end |> all
end

# https://github.com/JuliaAI/MLJIteration.jl/issues/36
@testset "integration test related to MLJIteration.jl #36" begin
    model = SquareRooter(NaN)
    report =
        IC.train!(model, Step(1), InvalidValue(), NumberLimit(3), verbosity=0)
    @test report[1][2].new_iterations == 1
    @test occursin("Stopping", report[2][2].log)
end


text:
struct QuadratureRule{dim,shape,T}
    weights::Vector{T}
    points::Vector{Vec{dim,T}}
end

getweights(qr::QuadratureRule) = qr.weights

getpoints(qr::QuadratureRule) = qr.points

#quadratures
struct GrundmannMoeller <: AbstractQuadratureRule end
struct Strang <: AbstractQuadratureRule end
struct DefaultQuad <: AbstractQuadratureRule end
struct GaussLegendre <: AbstractQuadratureRule end

# utils
function (::Type{QuadratureRule{2,RefTetrahedron}})(quad_type::DefaultQuad, order::Int)
    if order <= 6
        return QuadratureRule{2,RefTetrahedron}(Strang(),order)
    elseif order > 6 && isodd(order)
        s = Int((order-1)/2)
        return QuadratureRule{2,RefTetrahedron}(GrundmannMoeller(),s)
    else
        throw(ArgumentError("Quadrature rule of order $order not available"))
    end
end

# Get GaussLegendre weigths and points from FastGaussQuadrature
function (::Type{QuadratureRule{1,RefTetrahedron}})(quad_type::GaussLegendre, order::Int)
    points, weights = gausslegendre(order)
    # Shift interval from (-1,1) to (0,1)
    weights *= 0.5
    points = points .+ 1.0; points /= 2.0
    return QuadratureRule{1,RefTetrahedron,Float64}(weights, [Tensor{1,1}([x]) for x in points])
end

function (::Type{QuadratureRule{1,RefTetrahedron}})(quad_type::DefaultQuad, order::Int)
    return QuadratureRule{1,RefTetrahedron}(GaussLegendre(), order)
end

"""
Face Quadrature Rules
"""
function create_face_quad_rule(quad_rule::QuadratureRule{1,shape,T}, itp::Interpolation{2,shape}) where {T,shape}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    geom_face_interpol = get_default_geom_interpolator(shape, Val{1})
    face_quad_rule = QuadratureRule{2,shape,T}[]
    _populate_face_quad_rule!(face_quad_rule, geom_face_interpol, p, w)
end

function create_face_quad_rule(quad_rule::QuadratureRule{1,shape,T}) where {T,shape}
    w = getweights(quad_rule)
    p = getpoints(quad_rule)
    geom_face_interpol = get_default_geom_interpolator(shape, Val{1})
    face_quad_rule = QuadratureRule{2,shape,T}[]
    _populate_face_quad_rule!(face_quad_rule, geom_face_interpol, p, w)
end

function _populate_face_quad_rule!(face_quad_rule::Vector{QuadratureRule{2,shape,T}}, geom_face_interpol::Interpolation{1,shape,order}, p::Vector{Vec{1,T}}, w::Vector{T}) where {shape, order, T}
    face_coords = reference_edges(shape, Val{2})
    for j = 1:get_num_faces(shape, Val{2})
        new_points = Vec{2,T}[]
        for (qp, ξ) in enumerate(p)
            #Map from reference dim-1 shape to reference dim shape face/edge
            η = zero(Vec{2,T})
            for (k,x) in enumerate(face_coords[j])
                η += value(geom_face_interpol, k, ξ)*x
            end
            push!(new_points, η)
        end
        push!(face_quad_rule, QuadratureRule{2,shape,T}(w, new_points))
    end
    return face_quad_rule
end

#OTher Functions
"""
function integrate(f::function,qr::QuadratureRule)
integrate function f on reference shape of quadrature qr
"""
function integrate(f::Function,qr::QuadratureRule)
    p = getpoints(qr); w = getweights(qr)
    int_val = zero(typeof(f(p[1])))
    for (i,x) in enumerate(p)
        int_val += f(x)*w[i]
    end
    int_val
end


text:
export monitor

"""
Use: monitor(username::String;
                delay::Number = 0.05,
                iterations::Int = 10,
                command::String = "julia",
                logfile::String = Pkg.dir("CPUmonitor")*"/log/$(now()).txt",
                verbose::Bool = true)

Tracks CPU stats on local processor.

Saves the results in /myjulia/CPUmonitor/log/
"""
function monitor(username::String;
                    delay::Number = 0.05,
                    iterations::Int = 10,
                    command::String = "julia",
                    logfile::String = tempname()*".CPUmonitor.log.$(now()).txt",
                    verbose::Bool = true)
    verbose && println("Log: $logfile")
    logfile = track(username, delay, iterations, command, logfile, verbose)
    stat = parselog(logfile, verbose)
end


text:
function regulateheart!(system::CVSystem,n::Int64)
    # change activation function parameters based on new HR
    push!(system.heart.activation.th,1/system.cns.H[n+2]*minTos);
    system.heart.activation.tau1 = 0.269*system.heart.activation.th[end];
    system.heart.activation.tau2 = 0.452*system.heart.activation.th[end];
    t = linspace(0,system.heart.activation.th[end],10000);
    g1 = (t/system.heart.activation.tau1).^system.heart.activation.m1;
    g2 = (t/system.heart.activation.tau2).^system.heart.activation.m2;
    push!(system.heart.activation.k,maximum((g1./(1+g1)).*(1./(1+g2)))^-1)

    # change max. ventriclular elastance
    push!(system.heart.lv.Emax,system.cns.Emaxlv[n+2]);
    push!(system.heart.rv.Emax,system.cns.Emaxrv[n+2]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
struct Cartpole{T} <: AbstractModel
    mc::T
    mp::T
    l::T
    g::T
end

Cartpole() = Cartpole(1.0, 0.2, 0.5, 9.81)

function RobotDynamics.dynamics(model::Cartpole, x, u)
    mc = model.mc  # mass of the cart in kg (10)
    mp = model.mp   # mass of the pole (point mass at the end) in kg
    l = model.l   # length of the pole in m
    g = model.g  # gravity m/s^2

    q = x[ @SVector [1,2] ]
    qd = x[ @SVector [3,4] ]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp*g*l*s]
    B = @SVector [1, 0]

    qdd = -H\(C*qd + G - B*u[1])
    return [qd; qdd]
end

RobotDynamics.state_dim(::Cartpole) = 4
RobotDynamics.control_dim(::Cartpole) = 1
RobotDynamics.diffmethod(::Cartpole) = RobotDynamics.ForwardAD()


text:
"""
Calculating the input number times four

  timesFour(times_by::integer)

Return 'times_by*4'
"""
function timesFour(times_by::Integer)
    times_by*4
end

max_stars_repo_name: AnimatorOfSouls/Random-Projects

max_issues_repo_name: AnimatorOfSouls/Random-Projects

max_forks_repo_name: AnimatorOfSouls/Random-Projects


text:
# using Revise
using Test, PseudoArcLengthContinuation, LinearAlgebra, Setfield, SparseArrays, ForwardDiff
const PALC = PseudoArcLengthContinuation

n = 250*150
M = 30
par = nothing
sol0 = rand(2n)				# 585.977 KiB
orbitguess_f = rand(2n*M+1)	# 17.166MiB
pb = PeriodicOrbitTrapProblem(
			(x, p) -> x.^2,
			(x, p) -> (dx -> 2 .* dx),
			rand(2n),
			rand(2n),
			M)

pbg = PeriodicOrbitTrapProblem(
			(x, p) -> x.^2,
			(x, p) -> (dx -> 2 .* dx),
			pb.ϕ,
			pb.xπ,
			M ; ongpu = true)

pbi = PeriodicOrbitTrapProblem(
			(o, x, p) -> o .= x.^2 ,
			((o, x, p, dx) -> o .= 2 .* dx),
			pb.ϕ,
			pb.xπ,
			M ; isinplace = true)
@test PALC.isInplace(pb) == false
# @time PALC.POTrapFunctional(pb, res, orbitguess_f)
# @time PALC.POTrapFunctional(pbi, res, orbitguess_f)
res = @time pb(orbitguess_f, par)
resg = @time pbg(orbitguess_f, par)
resi = @time pbi(orbitguess_f, par)
@test res == resi
@test res == resg

res = @time pb(orbitguess_f, par, orbitguess_f)
resg = @time pbg(orbitguess_f, par, orbitguess_f)
resi = @time pbi(orbitguess_f, par, orbitguess_f)
@test res == resi
@test res == resg

@time PALC.POTrapFunctional!(pbi, resi, orbitguess_f, par)
@time PALC.POTrapFunctionalJac!(pbi, resi, orbitguess_f, par, orbitguess_f)
@test res == resi

# @code_warntype PALC.POTrapFunctional!(pbi, resi, orbitguess_f)

# using BenchmarkTools
# @btime pb($orbitguess_f) 					# 13.535 ms (188 allocations: 34.34 MiB)
# @btime pbi($orbitguess_f) 					# 7.869 ms (128 allocations: 17.17 MiB)
# @btime pb($orbitguess_f, $orbitguess_f) 	# 25.500 ms (373 allocations: 51.51 MiB)
# @btime pbi($orbitguess_f, $orbitguess_f)  	# 12.595 ms (253 allocations: 17.18 MiB)
# @btime PALC.POTrapFunctional!($pbi, $resi, $orbitguess_f) # 7.104 ms (126 allocations: 5.88 KiB)
# @btime PALC.POTrapFunctionalJac!($pbi, $resi, $orbitguess_f, $orbitguess_f) # 10.528 ms (251 allocations: 11.69 KiB)
#
#
# using IterativeSolvers, LinearMaps
#
# Jmap = LinearMap{Float64}(dv -> pbi(orbitguess_f, dv), 2n*M+1 ; ismutating = false)
# gmres(Jmap, orbitguess_f; verbose = false, maxiter = 1)
# @time gmres(Jmap, orbitguess_f; verbose = false, maxiter = 10)
#
# Jmap! = LinearMap{Float64}((o, dv) -> PALC.POTrapFunctionalJac!(pbi, o, orbitguess_f, dv), 2n*M+1 ; ismutating = true)
# gmres(Jmap!, orbitguess_f; verbose = false, maxiter = 1)
# @time gmres(Jmap!, orbitguess_f; verbose = false, maxiter = 10)
#
# # @code_warntype PALC.POTrapFunctional!(pbi, resi, orbitguess_f)
# # @profiler PALC.POTrapFunctionalJac!(pbi, resi, orbitguess_f, orbitguess_f)
#
# Jmap2! = LinearMap{Float64}((o, dv) -> pbi(o, orbitguess_f, dv), 2n*M+1 ; ismutating = true)
# gmres(Jmap2!, orbitguess_f; verbose = false, maxiter = 1)
# @time gmres(Jmap2!, orbitguess_f; verbose = false, maxiter = 10)
#
# Jmap3! = LinearMap{Float64}((o, dv) -> (o .= pbi( orbitguess_f, dv)), 2n*M+1 ; ismutating = true)
# gmres(Jmap3!, orbitguess_f; verbose = false, maxiter = 1)
# @time gmres(Jmap3!, orbitguess_f; verbose = false, maxiter = 10)
#
# using ProfileView, Profile
# @profview gmres!(res, Jmap!, orbitguess_f; verbose = false, maxiter = 1)
# @profview gmres!(res, Jmap!, orbitguess_f; verbose = false, maxiter = 10)

####################################################################################################
# test whether we did not make any mistake in the improved version of the PO functional
function _functional(poPb, u0, p)
	M = poPb.M
	N = poPb.N
	T = u0[end]
	h = T / M

	u0c = reshape(u0[1:end-1], N, M)
	outc = similar(u0c)
	outc[:, 1] .= (u0c[:, 1] .- u0c[:, M-1]) .- h/2 .* (poPb.F(u0c[:, 1], p) .+ poPb.F(u0c[:, M-1], p))

	for ii = 2:M-1
		outc[:, ii] .= (u0c[:, ii] .- u0c[:, ii-1]) .- h/2 .* (poPb.F(u0c[:, ii], p) .+ poPb.F(u0c[:, ii-1], p))

	end

	# closure condition ensuring a periodic orbit
	outc[:, M] .= u0c[:, M] .- u0c[:, 1]

	return vcat(vec(outc),
			dot(u0c[:, 1] .- poPb.xπ, poPb.ϕ)) # this is the phase condition
end

function _dfunctional(poPb, u0, p, du)
	# jacobian of the functional

	M = poPb.M
	N = poPb.N
	T = u0[end]
	dT = du[end]
	h = T / M
	dh = dT / M

	u0c = reshape(u0[1:end-1], N, M)
	duc = reshape(du[1:end-1], N, M)
	outc = similar(u0c)

	outc[:, 1] .= (duc[:, 1] .- duc[:, M-1]) .- h/2 .* (poPb.J(u0c[:, 1], p)(duc[:, 1]) .+ poPb.J(u0c[:, M-1], p)(duc[:, M-1]))

	for ii = 2:M-1
		outc[:, ii] .= (duc[:, ii] .- duc[:, ii-1]) .- h/2 .* (poPb.J(u0c[:, ii], p)(duc[:, ii]) .+ poPb.J(u0c[:, ii-1], p)(duc[:, ii-1]))
	end

	outc[:, 1] .-=  dh/2 .* (poPb.F(u0c[:, 1], p) .+ poPb.F(u0c[:, M-1], p))
	for ii = 2:M-1
		outc[:, ii] .-= dh/2 .* (poPb.F(u0c[:, ii], p) .+ poPb.F(u0c[:, ii-1], p))
	end

	# closure condition ensuring a periodic orbit
	outc[:, M] .= duc[:, M] .- duc[:, 1]

	return vcat(vec(outc),
			dot(duc[:, 1], poPb.ϕ)) # this is the phase condition

end

res = @time pb(orbitguess_f, par)
_res = _functional(pb, orbitguess_f, par)
@test res ≈ _res

_du = rand(length(orbitguess_f))
res = @time pb(orbitguess_f, par, _du)
_res = _dfunctional(pb, orbitguess_f, par, _du)
@test res ≈ _res

####################################################################################################
# test whether the analytical version of the Jacobian is right
n = 50
pbsp = PeriodicOrbitTrapProblem(
			(x, p) -> cos.(x),
			(x, p) -> spdiagm(0 => -sin.(x)),
			rand(2n),
			rand(2n),
			10)
orbitguess_f = rand(2n*10+1)
dorbit = rand(2n*10+1)
Jfd = sparse( ForwardDiff.jacobian(x->pbsp(x, par), orbitguess_f) )
Jan = pbsp(Val(:JacFullSparse), orbitguess_f, par)
@test norm(Jfd - Jan, Inf) < 1e-6
####################################################################################################
# test whether the inplace version of computation of the Jacobian is right
n = 1000
pbsp = PeriodicOrbitTrapProblem(
			(x, p) -> x.^2,
			(x, p) -> spdiagm(0 => 2 .* x),
			rand(2n),
			rand(2n),
			M)

sol0 = rand(2n)
orbitguess_f = rand(2n*M+1)
Jpo = pbsp(Val(:JacFullSparse), orbitguess_f, par)
Jpo2 = copy(Jpo)
pbsp(Val(:JacFullSparseInplace), Jpo2, orbitguess_f, par)
@test nnz(Jpo2 - Jpo) == 0
####################################################################################################
# test of the version with inhomogenous time discretisation
pbsp = PeriodicOrbitTrapProblem(
			(x, p) -> cos.(x),
			(x, p) -> spdiagm(0 => -sin.(x)),
			rand(2n),
			rand(2n),
			10)

pbspti = PeriodicOrbitTrapProblem(
			(x, p) -> cos.(x),
			(x, p) -> spdiagm(0 => -sin.(x)),
			pbsp.ϕ,
			pbsp.xπ,
			ones(9) ./ 10)
PALC.getM(pbspti)
orbitguess_f = rand(2n*10+1)
@test pbspti.xπ ≈ pbsp.xπ
@test pbspti.ϕ ≈ pbsp.ϕ
pbspti(orbitguess_f, par)
@test pbsp(orbitguess_f, par) ≈ pbspti(orbitguess_f, par)


text:
"""
Modia module with electric component models (inspired from Modelica Standard Library).

* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2021: Hilding Elmqvist
* License: MIT (expat)

"""
#module Electric

using ModiaLang
using Unitful

Pin = Model( v = potential, i = flow )

OnePort = Model( p = Pin, n = Pin, equations = :[
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i ]
    )

"""
  Resistor(R=1.0u"Ω")

Electrical resistor

`R` - Resistance Ω
"""
Resistor = OnePort | Model( R = 1.0u"Ω", equations = :[ R*i = v ] )
# @showModel(Resistor)

Capacitor = OnePort | Model( C = 1.0u"F", v=Var(init=0.0u"V"), equations = :[ C*der(v) = i ] )

Inductor = OnePort | Model( L = 1.0u"H", i=Var(init=0.0u"A"), equations = :[ L*der(i) = v ] )

ConstantVoltage = OnePort | Model( V = 1.0u"V", equations = :[ v = V ] )

Ground = Model( p = Pin, equations = :[ p.v = 0.0u"V" ] )

# Ideal operational amplifier (norator-nullator pair), but 3 pins
IdealOpAmp3Pin = Model(
    in_p = Pin,
    in_n = Pin,
    out = Pin,
    equations = :[
        in_p.v = in_n.v
        in_p.i = 0u"A"
        in_n.i = 0u"A" ]
)

# Partial generic voltage source using the input signal as source voltage
PartialSignalVoltage = Model(
    v = input,
    p = Pin,
    n = Pin
)

# Generic voltage source using the input signal (without and with unit) as source voltage
SignalVoltage = PartialSignalVoltage | Model(
    equations = :[
        p.v - n.v = v
                0 = p.i + n.i
                i = p.i ]
)
UnitlessSignalVoltage = PartialSignalVoltage | Model(
    equations = :[
        p.v - n.v = v*u"V"
                0 = p.i + n.i
                i = p.i ]
)

# Partial sensor to measure the current in a branch
PartialCurrentSensor = Model(
    i = output,
    p = Pin, # (info = "Positive pin")
    n = Pin, # (info = "Negative pin")
)

# Sensor to measure the current in a branch
CurrentSensor = PartialCurrentSensor | Model(
    equations = :[
        p.v = n.v
          0 = p.i + n.i
          i = p.i]
)
UnitlessCurrentSensor = PartialCurrentSensor | Model(
    equations = :[
        p.v = n.v
        0 = p.i + n.i
        i = p.i/u"A"]
)

#=
# Step voltage source
@model StepVoltage begin
    V = Parameter(1.0, start = 1.0, info = "Voltage") #, T = Unitful.V)
    startTime = Parameter(0.0, start = 0.0, info = "Start time") # , T = Unitful.s)
    OnePort()
equations
    v = if time < startTime; 0 else V end
end

@model VoltageSource begin
    OnePort()
    offset = Par(0.0) # Voltage offset
    startTime = Par(0.0) # Time offset
    signalSource = SignalSource(offset=offset, startTime=startTime)
equations
    v = signalSource.y
end

#=
@model SineVoltage1 begin
    # Sine voltage source
    V = Parameter() # Amplitude of sine wave
    phase = Par(0.0) # Phase of sine wave
    freqHz = Parameter() # Frequency of sine wave
    VoltageSource(signalSource=Sine(amplitude=V, freqHz=freqHz, phase=phase))
end
=#

# Sinusoidal voltage source
@model SineVoltage begin
    V = Parameter() # Amplitude of sine wave
    phase = Par(0.0) # Phase of sine wave
    freqHz = Parameter() # Frequency of sine wave
    VoltageSource()
equations
    v = V*sin(10*time)
end

# Ideal diode
@model IdealDiode begin # Ideal diode
    OnePort()
    Ron = Par(1.0E-2) # Forward state-on differential resistance (closed diode resistance)
    Goff = Par(1.0E-2) # Backward state-off conductance (opened diode conductance)
    Vknee = Par(0) # Forward threshold voltage
    # off = Variable(start=true) # Switching state
    s = Float(start=0.0) # Auxiliary variable for actual position on the ideal diode characteristic
    #=
    s = 0: knee point
    s < 0: below knee point, diode conducting
    s > 0: above knee point, diode locking
    =#
equations
    # off := s < 0
#        v = s * if !positive(s); 1 else Ron end + Vknee
#        i = s * if !positive(s); Goff else 1 end + Goff * Vknee
    v = s * if !(s>0); 1 else Ron end + Vknee
    i = s * if !(s>0); Goff else 1 end + Goff * Vknee
end

@model Diode begin
    OnePort()
    Ids=Par(1.e-6) # "Saturation current";
    Vt=Par(0.04) # "Voltage equivalent of temperature (kT/qn)";
    Maxexp = Par(15) # "Max. exponent for linear continuation";
    R=Par(1.e8) # "Parallel ohmic resistance";
equations
    i = if v/Vt > Maxexp; Ids*(exp(Maxexp)*(1 + v/Vt - Maxexp) - 1) else Ids*(exp(v/Vt) - 1) end + v/R
end
=#
#end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang