text:
#
#	MPO structure
#
struct MPO{DT<:Union{localOp, AbstractArray}}
	bDim::Union{Int, Tuple{Int,Int}}
	Operator::Vector{DT}
	Op_index::Vector{Int}

	function MPO(bDim::Union{Tuple{Int,Int}, Int}, Operator::Union{Vector{<:localOp}, Vector{<:AbstractArray}}, Op_index::Vector{Int})
		allunique(Op_index) || throw(ArgumentError("operator indices contains duplicate"))
		dim = typeof(bDim) == Int ? bDim*bDim : *(bDim...)
		dim >= maximum(Op_index) || throw(ArgumentError("operator indices exceed MPO bond dimension"))

		return new{eltype(Operator)}(bDim, Operator, Op_index)
	end
end

# define MPOsparseT
const MPOsparseT = MPO{<:localOp}

import Base.eltype
eltype(mpo::MPOsparseT) = promote_type(eltype.(mpo.Operator)...)

getBondDim(mpo::MPO{DT}) where {DT<:Union{localOp, AbstractArray}} = typeof(mpo.bDim) == Int ? (mpo.bDim, mpo.bDim) : mpo.bDim

physDim(Op::localOp{DT, D}) where {DT<:Number, D} = D
physDim(mpo::MPO) = physDim(mpo.Operator[1])

import Base.getindex
function getindex(mpo::MPO{DT}, i::Int) where {DT<:Union{localOp, AbstractArray}}
	idx = findfirst(x->x == i, mpo.Op_index)

	if idx == nothing
		return nothing
	else
		return mpo.Operator[idx]
	end
end

getindex(mpo::MPO{DT}, i::Int, j::Int) where {DT<:Union{localOp, AbstractArray}} = getindex(mpo, LinearIndices((getBondDim(mpo)[1], getBondDim(mpo)[2]))[i,j])
getindex(mpo::MPO{DT}, i::Int, j::Colon) where {DT<:Union{localOp, AbstractArray}} = [getindex(mpo, x) for x in LinearIndices((getBondDim(mpo)[1], getBondDim(mpo)[2]))[i,:]]
getindex(mpo::MPO{DT}, i::Colon, j::Int) where {DT<:Union{localOp, AbstractArray}} = [getindex(mpo, x) for x in LinearIndices((getBondDim(mpo)[1], getBondDim(mpo)[2]))[:,j]]

"""
    eachLocalOp(mpo)

Iterates over all localOp and returns them in a tuple with the index position of the MPO.

# Arguments:
    - mpo:

return:

	!!!ADD example!!!

"""
function eachLocalOp(mpo::MPO{DT}) where {DT<:Union{localOp, AbstractArray}}
	idx = CartesianIndices(getBondDim(mpo))[mpo.Op_index]
	idx = [(x[1],x[2]) for x in idx]

	return zip(mpo.Operator, idx)
end

#
#	this should be changed
#
"""
	iteration over all terms in MPO
"""
function eachoperation(mpo::MPOsparseT)
	tmp = []
	@inbounds for x in eachLocalOp(mpo)
		for y in eachoperation(x[1])
			push!(tmp, [y[1], Tuple(y[2:3]), x[2]])
		end

	end
	return tmp
end

import Base.show
function Base.show(io::IO, mpo::MPOsparseT)

	bDim = getBondDim(mpo)
	idx = [CartesianIndices(bDim)[x] for x in mpo.Op_index]
	B = zeros(Int, bDim)
	B[idx] .= [1:length(mpo.Op_index)...]

	print(io, "operators in MPO:\n $B \n")

	for (i,x) in enumerate(mpo.Operator)

		print(io, " operator: $i --> $(Array(x))\n")

	end
end

function getMPOTensor(mpo::MPO{DT}) where {DT<:Union{localOp, AbstractArray}}
	dims = getBondDim(mpo)
	phys_dim = physDim(mpo)
	data_type = eltype(mpo)
	MPOTensor = zeros(data_type, (phys_dim, dims..., phys_dim))

	for (op, (i,j)) in eachLocalOp(mpo)
		MPOTensor[:, i, j, :] = Array(op)
	end

	return MPOTensor
end

function getMPOFromTensor(mpoTensor::Array{DT, 4}) where {DT<:Number}
	bond_dim = size(mpoTensor)[2:3]

	operator = localOp[]
	idx = Int[]
	for col in 1:bond_dim[2]
		for row in 1:bond_dim[1]

			op = mpoTensor[:,row,col,:]

			if iszero(op)!=true
				push!(operator, localOp(op))
				push!(idx, LinearIndices(bond_dim)[row, col])
			end
		end
	end

	mpo = MPO(bond_dim, operator, idx)
	return mpo
end

#
#	some MPO arithmetic
#
import Base.(*)
import Base.(+)

(*)(a::Number, MPOvec::Vector{<:MPO{DT}}) where {DT<:Union{localOp, AbstractArray}} = Vector{MPO{DT}}([MPO(MPOvec[1].bDim, a .* MPOvec[1].Operator, MPOvec[1].Op_index),MPOvec[2:end]...])
(+)(mpo_a::Vector{<:MPO{DT}}, mpo_b::Vector{<:MPO{DT}}) where {DT<:Union{localOp, AbstractArray{<:Number}}} = add_MPO(mpo_a, mpo_b)

function add_MPO(MPOvecA::Vector{<:MPO{DT}}, MPOvecB::Vector{<:MPO{DT}}) where {DT<:Union{localOp, AbstractArray{<:Number}}}
	@assert length(MPOvecA) == length(MPOvecB) "you are trying to add
					to MPO with different length"

	newMPOvec = Vector{MPO{DT}}([])

	for i in 1:size(MPOvecA)[1]

		MPOA = MPOvecA[i]
		MPOB = MPOvecB[i]

		bDim_A = MPOA.bDim
		bDim_B = MPOB.bDim

		bDim_A = typeof(bDim_A) == Int ? (bDim_A, bDim_A) : bDim_A
		bDim_B = typeof(bDim_B) == Int ? (bDim_B, bDim_B) : bDim_B

		# take care of last and first site
		new_bin = i != 1 ? bDim_A[1]+bDim_B[1] : 1
		new_bout = i != size(MPOvecA)[1] ? bDim_A[2]+bDim_B[2] : 1

		OP_vec_A = MPOA.Operator
		OP_idx_A = MPOA.Op_index

		OP_vec_B = MPOB.Operator
		OP_idx_B = MPOB.Op_index

		#
		#	bring into block diagonal form
		#

		# adjust OP_idx_A
		idx = CartesianIndices((bDim_A[2],bDim_A[1]))[OP_idx_A]
		OP_idx_A = LinearIndices((new_bout,bDim_A[1]))[idx]

		# adjust Op_idx_B
		idx = CartesianIndices((bDim_B[2],bDim_B[1]))[OP_idx_B]
		idx = [Array([Tuple(x)...]).+[new_bout-bDim_B[2],new_bin-bDim_B[1]] for x in idx]
		OP_idx_B = [LinearIndices((new_bout, new_bin))[x...] for x in idx]
		new_Op_vec = Vector{DT}([OP_vec_A..., OP_vec_B...])
		new_Op_idx = Vector{Int}([OP_idx_A..., OP_idx_B...])

		newMPO = MPO((new_bin, new_bout), new_Op_vec, new_Op_idx)
		push!(newMPOvec, newMPO)
	end

	return newMPOvec
end

function (*)(MPOvecA::Vector{<:MPOsparseT}, MPOvecB::Vector{<:MPOsparseT})

	@assert length(MPOvecA) == length(MPOvecB) "you are trying to mulitply
						    to MPO with different length"

	new_MPOvec = MPO{localOp}[]
	for i in 1:length(MPOvecA)

		mpoA = getMPOTensor(MPOvecA[i])
		mpoB = getMPOTensor(MPOvecB[i])

		@tensor new_mpo[α, α′, d, d′, β, β′] := mpoA[d, α, β, γ]*mpoB[γ, α′, β′, d′]

		new_mpo = reshape(new_mpo, *(size(new_mpo)[1:2]...), size(new_mpo)[3], size(new_mpo)[4], *(size(new_mpo)[5:end]...))
		new_mpo = permutedims(new_mpo, [2,1,4,3])
		new_mpo = getMPOFromTensor(new_mpo)
		push!(new_MPOvec, new_mpo)
	end

	return new_MPOvec
end

#
#	multiplication of MPO with Array
#	(multiplication acts as matrix multiplication only on the bond dimension)
#

function (*)(mpo::MPO{localOp{DT,D}}, a::AbstractArray{DT, 2}) where {DT<:Number, D}
	tmp = []

	#multiplication from the right
	@inbounds for col in 1:size(a)[2]
		for row in 1:getBondDim(mpo)[1]
			array_col = @view a[:,col]  # get column
			idx = findall(x->x!=0, array_col) # find nz elements in column
			if !isempty(idx)
				mpo_row = mpo[row,:]  # eventually replace with a view
				idx_nothing = findall(x->typeof(x) == Nothing, mpo_row) #check if there is a localOp
				idx = setdiff!(idx,idx_nothing) #set difference of idx: make sure that element of idx is not in idx_nothing
				if !isempty(idx)
					mpo_row = mpo_row[idx] #[mpo[row,i] for i in idx] # get elements which are at nz column index
					coef = array_col[idx] # [array_col[i] for i in idx]
					new_localOp = sum(coef.*mpo_row) #obtain actual new localOp
					push!(tmp, (new_localOp, (row,col)))
				end
			end
		end

	end

	# construct new OP
	Operator = [x[1] for x in tmp]
	Op_index = [LinearIndices((getBondDim(mpo)[1], size(a)[2]))[x[2]...] for x in tmp]

	idx = findall(x->x!=nothing, Operator)
	Operator = Vector{localOp{DT,D}}(Operator[idx])
	Op_index = Op_index[idx]
	# this would just be necessary if we would loop first over row
	#idx = sortperm(Op_index)
	#Op_index = Op_index[idx]
	#Operator = Operator[idx]

	return MPO((getBondDim(mpo)[1], size(a)[2]), Operator, Op_index)
end

function (*)(a::AbstractArray{DT, 2}, mpo::MPO{localOp{DT,D}}) where {DT<:Number, D}
	tmp = []

	#multiplication from the left
	@inbounds for col in 1:getBondDim(mpo)[2]
		for row in 1:size(a)[1]

			array_row = @view a[row,:]  # get row
			idx = findall(x->x!=0, array_row) # find nz elements in row
			if !isempty(idx)
				mpo_col = mpo[:,col]
				idx_nothing = findall(x->typeof(x) == Nothing, mpo_col) #check if there is a localOp
				idx = setdiff!(idx,idx_nothing)
				if !isempty(idx)
					mpo_col = mpo_col[idx] #[mpo[i, col] for i in idx] # get elements which are at nz column index
					coef = array_row[idx] #[array_row[i] for i in idx]
					new_localOp = sum(coef.*mpo_col) #obtain actual new localOp
					push!(tmp, (new_localOp, (row,col)))
				end
			end
		end

	end

	# construct new OP
	Operator = [x[1] for x in tmp]
	Op_index = [LinearIndices((size(a)[1], getBondDim(mpo)[2]))[x[2]...] for x in tmp]

	idx = findall(x->x!=nothing, Operator)
	Operator = Vector{localOp{DT,D}}(Operator[idx])
	Op_index = Op_index[idx]

	# this would just be necessary if we would loop first over row
	#idx = sortperm(Op_index)
	#Op_index = Op_index[idx]
	#Operator = Operator[idx]

	return MPO((size(a)[1], getBondDim(mpo)[2]), Operator, Op_index)
end

#
#	function for deparalization of MPO (shrinks complexity)
#
#	so far only for localOp structure not defined for Array operator
#

function checkPara(a, b)
	@assert length(a) == length(b)

	parallel = false
	c = 0
	if typeof.(a) != typeof.(b)
		return (false, 0.0)
	end

	idx = findall(x->typeof(x) != Nothing, a)
	#tmp = [parallelOp(a[x], b[x]) for x in idx]
	tmp = [a[x]∥b[x] for x in idx]  # check if operator are parallel
	if length(unique([x[2] for x in tmp])) != 1 # check if the scale between them are identical
		return (false, 0.0)
	end

	return (tmp[1][1], tmp[1][2])
end

function findParaRow(mpo::MPO{localOp{DT,D}}) where {DT<:Number, D}
	a = zeros(ComplexF64, getBondDim(mpo)[1])
	a[1] = 1.0
	T = [a]

	mpo_row = mpo[1,:]
	K = Any[mpo_row]

	for row_i in length(K)+1:getBondDim(mpo)[1]
		para = false

		for k_idx in 1:length(K)

			mpo_row = mpo[row_i, :]
			k = K[k_idx]

			para, scale = checkPara(mpo_row, k)

			if para == true

				T[k_idx][row_i] = scale
				break
			end
		end

		if para != true
			push!(K, mpo_row)
			a = zeros(ComplexF64, getBondDim(mpo)[1])
			a[row_i] = 1.0
			push!(T, a)
		end

	end

	# recat K into mpo

	K = hcat(K...)
	K = K[Transpose(LinearIndices(size(K)))]
	idx = findall(x->typeof(x) != Nothing, K)
	op_idx = LinearIndices(size(K))[idx]
	K = MPO(size(K), Vector{localOp{DT,D}}(K[idx]), op_idx)

	return K, hcat(T...)
end

function findParaCol(mpo::MPO{localOp{DT,D}}) where {DT<:Number, D}
	a = zeros(ComplexF64, getBondDim(mpo)[2])
	a[1] = 1.0
	T = [a]
	mpo_col = mpo[:,1]
	K = Any[mpo_col]

	for col_i in length(K)+1:getBondDim(mpo)[2]
		para = false

		for k_idx in 1:length(K)
			mpo_col = mpo[:,col_i]
			k = K[k_idx]

			para, scale = checkPara(mpo_col, k)
			if para == true

				T[k_idx][col_i] = scale
				break
			end
		end

		if para != true

			push!(K, mpo_col)
			a = zeros(ComplexF64, getBondDim(mpo)[2])
			a[col_i] = 1.0
			push!(T,a)
		end
	end

	# recat K into mpo
	K = hcat(K...)
	idx = findall(x->typeof(x) != Nothing, K)
	op_idx = LinearIndices(size(K))[idx]
	K = MPO(size(K), Vector{localOp{DT,D}}(K[idx]), op_idx)

	return K, Transpose(hcat(T...))
end

#
#	deparallizes mpo-chain
#
function depara!(MPOvec::Vector{MPO{localOp{DT,D}}}) where {DT<:Number, D}
	for site in 1:length(MPOvec)
		#@show "dePara col", site
		mpo = MPOvec[site]

		mpo, T = findParaCol(mpo)
		MPOvec[site] = mpo

		if site+1 <= length(MPOvec)
			mpo_next = MPOvec[site+1]
			mpo_next = T*mpo_next
			MPOvec[site+1] = mpo_next
		end
	end

	for site in length(MPOvec):-1:1
		#@show "dePara row", site
		mpo = MPOvec[site]
		mpo, T = findParaRow(mpo)
		MPOvec[site] = mpo
		if site-1 > 1
			mpo_pre = MPOvec[site-1]
			mpo_pre = mpo_pre*T
			MPOvec[site-1] = mpo_pre
		end
	end

	for site in 1:length(MPOvec)
		#@show "dePara col", site
		mpo = MPOvec[site]

		mpo, T = findParaCol(mpo)
		MPOvec[site] = mpo

		if site+1 <= length(MPOvec)
			mpo_next = MPOvec[site+1]
			mpo_next = T*mpo_next
			MPOvec[site+1] = mpo_next
		end
	end

	return MPOvec
end

max_stars_repo_name: Gawatz/MPOmodule

max_issues_repo_name: Gawatz/MPOmodule

max_forks_repo_name: Gawatz/MPOmodule


text:
# Copyright (c) 2021 Leonardo Banchi <leonardo.banchi@unifi.it>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using LinearAlgebra, BlockDiagonals

" symplectic matrix "
omega(n::Int) = BlockDiagonal([[0 1.; -1. 0] for k=1:n])
omega(A::AbstractMatrix) = omega(Int(size(A,1)/2))

" arrange diagonal elements in a block form "
wdiag(D) = BlockDiagonal([ Diagonal([d,d]) for d in D])

" check if a matrix is symplectic "
issymplectic(S) = S*omega(S)*S'≈omega(S)

" generate a random symplectic matrix "
function rand_symplectic(n)
	A = randn(2n,2n)
	exp(Hermitian(A)*omega(n))
end

" generate a random covariance matrix "
function rand_covmat(n,σ=1)
	v = 1/2 .+ σ*rand(n)
	S = rand_symplectic(n)
	Symmetric(S*wdiag(v)*S')
end

function _will_init(V)
	@assert	size(V,1) == size(V,2) && iseven(size(V,1))
	Int(size(V,1)/2), omega(V)
end

"""
	D,S = williamson(V)

	find the Williamson decomposition such that
	S'*wdiag(D)*S ≈ V
"""
function williamson(V)
	n, Ω = _will_init(V)
	iv12 = inv(sqrt(Symmetric(V)))
	J = iv12*Ω*iv12
	f = schur((J-J')/2)
	T, Q = f.Schur, f.vectors
	P = BlockDiagonal([ sign(T[2k-1,2k]) > 0 ?
					   [1. 0; 0 1.] : [0 1.; 1. 0] for k=1:n])
	D = inv.(diag(P*T*P,1)[1:2:2n])
	iS = wdiag(sqrt.(D))*P*Q'*iv12
	D, Ω*iS*Ω
end

"""
	D,S = williamsonNew(V)

	find the Williamson decomposition such that
	S'*wdiag(D)*S ≈ V
	new algorithm.
"""
function williamsonNew(V)
	n, Ω = _will_init(V)
	λ = sort(real(eigvals(V*Ω*1im)))[n+1:2n]
	S = zeros(2n,2n)
    for m=1:n
	    denom = λ[m]*prod(λ[k]^2-λ[m]^2 for k=1:n if k!=m)
        Vλ = complex.(V)-(1im*λ[m])*Ω
	    minorden(k) = det(view(Vλ,setdiff(1:2n,2m),setdiff(1:2n,k)))/denom
        fam = real(minorden(2m))
	    for k=1:n
            faeven = k==m ? sqrt(fam) : minorden(2k)/sqrt(fam)
            faodd = -minorden(2k-1)/sqrt(fam)
            S[2m-1,2k] 	= -real(faodd)
            S[2m,2k] = imag(faodd)
            S[2m-1,2k-1] = real(faeven)
            S[2m,2k-1] = -imag(faeven)
        end
	end
	λ, S
end

ext: jl


text:
using FileIO, GeoMakie, ColorSchemes, GLMakie

export createmaps, plotmap

using GeoMakie.GeoJSON, GeoMakie.GeoInterface
function geotest()
    source = LonLat()
    dest = WinkelTripel()
    states = download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")
    states_geo = GeoJSON.read(read(states, String))
    lons = LinRange(-179.5, 179.5, 360 ÷ 2)
    lats = LinRange(-89.5, 89.5, 180 ÷ 2)
    field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]
    points = map(CartesianIndices(size(field))) do xy
        x, y = Tuple(xy)
        Point2f0(lons[x], lats[y])
    end

    proj = GeoMakie.proj(source, dest)
    limits = FRect2D(Makie.apply_transform(proj, points))
    fig = Figure()

    ax = Makie.Axis(fig[1,1])
    limits!(limits)
    ax.scene.transformation.transform_func[] = proj

    # wireframe!(ax, lons, lats, field, color=(:gray, 0.2), transparency=true)
    # wireframe!(-180:15.0:180, -90:15.0:90, ones(25,13), color=(:black, 1.0), transparency=true)
    # wireframe!(repeat(-180:15.0:180, inner=15), repeat(-90:15.0:90, inner=15), ones(25*15,13*15), color=(:black, 1.0), transparency=true)
    # xgrid, ygrid = xygrid(repeat(-180:15.0:180, inner=15), repeat(-90:15.0:90, inner=15))
    # Proj4.transform!(source, dest, vec(xgrid), vec(ygrid))
    # wireframe!(xgrid, ygrid, ones(25*15,13*15), color=(:black, 1.0), transparency=true)

    n = length(GeoInterface.features(states_geo))
    lines!(ax, GeoMakie.coastlines())
    poly!(ax, states_geo, color= 1:n, strokecolor = :blue, strokewidth = 1)
    # display(fig)
    # nothing
    Makie.save("geomakie_figtest.png", ax.scene, resolution=(1000,600))
    GLMakie.destroy!(GLMakie.global_gl_screen())
end

function createmap(gisregion, regions, regionlist, lons, lats, colors, source, dest, xs, ys,
                    landcenters, popcenters, connected, connectedoffshore;
                    lines=false, labels=false, resolutionscale=1, textscale=1, legend=false)
    nreg = length(regionlist)
    scale = maximum(size(regions))/6500

    regions[regions.==NOREGION] .= nreg + 1

    xmin, xmax = extrema(xs)
    ymin, ymax = extrema(ys)
    aspect_ratio = (ymax - ymin) / (xmax - xmin)
    pngwidth = round(Int, resolutionscale*1.02*size(regions,1)) # allow for margins (+1% on both sides)
    pngsize = pngwidth, round(Int, pngwidth * aspect_ratio)     # use aspect ratio after projection transformation

    println("...constructing map...")
    proj = GeoMakie.proj(source, dest)
    fig = Figure()
    ax = Makie.Axis(fig[1,1], xtickcolor=RGBA(colorant"black", 0.4), ytickcolor=RGBA(colorant"black", 0.4),
                        xtickwidth=scale*2, ytickwidth=scale*2)
    # ax.scene.transformation.transform_func[] = proj
    # surface!(xs, ys; color=float.(regions), colormap=colors, shading=false, transparency=true)
    # xgrid, ygrid = xygrid(-180:15.0:180, -90:15.0:90)
    # display(xgrid)
    # Proj4.transform!(source, dest, vec(xgrid), vec(ygrid))
    # wireframe!(xgrid, ygrid, ones(25,13), color=(:black, 1.0), transparency=true)

    if lines
        println("...drawing transmission lines...")
        for (i,conn) in enumerate([connected, connectedoffshore])
            for reg1 = 1:nreg, reg2 = 1:nreg
                reg1 <= reg2 && continue
                !conn[reg1, reg2] && continue
                # line = [lonlatpoint(popcenters[reg1,:]), lonlatpoint(popcenters[reg2,:])]     # straight line on projected map
                line = greatcircletrack(popcenters[reg1,:], popcenters[reg2,:], 20)             # great circle segments
                projectedline = Makie.apply_transform(proj, Point2f0.(line))
                color = (i == 1) ? :black : :white
                lines!(projectedline, color=color, linewidth=resolutionscale*scale*textscale*5)
            end
        end
    end
    if labels
        for reg = 1:nreg
            pos = Makie.apply_transform(proj, lonlatpoint(popcenters[reg,:]))
            text!(string(regionlist[reg]); position=pos, align=(:center,:center), textsize=textscale*scale*resolutionscale*100)
        end
    end

    # draw the high resolution surface last, so the above calls to lines! go 20x faster
    surface!(xs, ys; color=float.(regions), colormap=colors, shading=false)

    println("...saving...")
    mkpath(in_datafolder("output"))
    filename = in_datafolder("output", "$gisregion.png")
    isfile(filename) && rm(filename)
    Makie.save(filename, ax.scene, resolution=pngsize)
    if legend
        makelegend(string.(regionlist), colors[2:end-1], scale=(4*scale)^0.7)
        img = load(filename)
        legend = autocrop(load(in_datafolder("output", "legend.png")))
        height = max(size(img,1), size(legend,1))
        legendwidth = size(legend, 2)
        combined = [ypad(img, height) ypad(legend, height) fill(RGB{N0f8}(1.0,1.0,1.0), (height, legendwidth÷5))]
        FileIO.save(filename, combined)
        rm(in_datafolder("output", "legend.png"))
    end
    GLMakie.destroy!(GLMakie.global_gl_screen())
end

function createmaps(gisregion; scenarioyear="ssp2_2050", lines=true, labels=true, resolutionscale=1, textscale=1, randseed=1, downsample=1)
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    regions = regions[1:downsample:end, 1:downsample:end]
    offshoreregions = offshoreregions[1:downsample:end, 1:downsample:end]
    lonrange = lonrange[1:downsample:end]
    latrange = latrange[1:downsample:end]
    textscale *= downsample
    nreg = length(regionlist)

    println("Mapping colors to regions (avoid same color in adjacent regions)...")
    mergeregions = copy(regions)
    mergeregions[regions.==0] .= offshoreregions[regions.==0]
    mergeconnected = connectedregions(mergeregions, nreg)
    colorindices = greedycolor(mergeconnected, 1:7, 1:nreg, randseed=randseed)
    onshorecolors = [RGB(0.2,0.3,0.4); colorschemes[:Set2_7].colors[colorindices]; RGB(0.4,0.4,0.4)]
    offshorecolors = [RGB(0.4,0.4,0.4); colorschemes[:Set2_7].colors[colorindices]; RGB(0.2,0.3,0.4)]
    # onshorecolors = RGBA.([RGB(0.2,0.3,0.4); colorschemes[:Set2_7].colors[colorindices]; RGB(0.4,0.4,0.4)], 0.8)
    # offshorecolors = RGBA.([RGB(0.4,0.4,0.4); colorschemes[:Set2_7].colors[colorindices]; RGB(0.2,0.3,0.4)], 0.8)

    println("\nProjecting coordinates (Mollweide)...")
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    source = LonLat()
    dest = Projection("+proj=moll +lon_0=$(mean(lons)) +ellps=WGS84")
    xs, ys = xygrid(lons, lats)
    Proj4.transform!(source, dest, vec(xs), vec(ys))

    println("\nFinding interregional transmission lines...")
    geocenters, popcenters = getregioncenters(regions, nreg, lonrange, latrange, res, scenarioyear)   # column order (lat,lon)
    landcenters = find_landarea_near_popcenter(regions, nreg, popcenters, lonrange, latrange, res)
    connected = connectedregions(regions, nreg)
    connectedoffshore = connectedregions(offshoreregions, nreg)
    connectedoffshore[connected] .= false

    println("\nOnshore map...")
    createmap(gisregion, regions, regionlist, lons, lats, onshorecolors, source, dest, xs, ys,
        landcenters, popcenters, connected, connectedoffshore; lines, labels, resolutionscale, textscale)
    println("\nOffshore map...")
    createmap("$(gisregion)_offshore", offshoreregions, regionlist, lons, lats, offshorecolors, source, dest, xs, ys,
        landcenters, popcenters, connected, connectedoffshore, lines=false, labels=false, resolutionscale=resolutionscale, textscale=textscale)
    # exit()
    return nothing
end

xygrid(lons, lats) = [lon for lon in lons, lat in lats], [lat for lon in lons, lat in lats]

# ColorBrewer Set2_7:  https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#colorbrewer-1
function maskmap(mapname, regions, regionlist, lonrange, latrange;
                    resolutionscale=1, textscale=1, randseed=1, legend=false, downsample=1)
    regions = regions[1:downsample:end, 1:downsample:end]
    lonrange = lonrange[1:downsample:end]
    latrange = latrange[1:downsample:end]
    nreg = length(regionlist)

    colors = [
        RGB([174,140,114]/255...),  # bad land type
        RGB([230,0,0]/255...),      # high population
        RGB([255,100,255]/255...),  # protected area
        RGB([120,170,80]/255...),   # no grid
        RGB([255,217,47]/255...),   # solar plant A
        RGB([255,150,0]/255...),    # solar plant B
        RGB([140,156,209]/255...),  # wind plant A
        RGB([110,136,169]/255...),  # wind plant B
    ]

    println("Mapping colors to regions (avoid same color in adjacent regions)...")
    connected = connectedregions(regions, nreg)
    # colorindices = (nreg > 7) ? greedycolor(connected, 1:7, 1:nreg, randseed=randseed) : collect(1:nreg)
    # colors = colorschemes[Symbol("Set2_$nreg")].colors[colorindices]
    onshorecolors = [RGB(0.3,0.3,0.45); colors; RGB(0.4,0.4,0.4)]

    println("\nProjecting coordinates (Mollweide)...")
    res = 0.01
    res2 = res/2
    lons = (-180+res2:res:180-res2)[lonrange]         # longitude values (pixel center)
    lats = (90-res2:-res:-90+res2)[latrange]          # latitude values (pixel center)
    source = LonLat()
    dest = Projection("+proj=moll +lon_0=$(mean(lons)) +ellps=WGS84")
    xs, ys = xygrid(lons, lats)
    Proj4.transform!(source, dest, vec(xs), vec(ys))

    println("\nOnshore map...")
    createmap(mapname, regions, regionlist, lons, lats, onshorecolors, source, dest, xs, ys,
        [], [], connected, connected, lines=false, labels=false, resolutionscale=resolutionscale, textscale=textscale, legend=legend)
    return nothing
end

function autocrop(img, padding::Int=0)
    rowrange = dataindexes_lat(vec(any(img .!= RGB(1.0,1.0,1.0), dims=2)), padding)
    colrange = dataindexes_lat(vec(any(img .!= RGB(1.0,1.0,1.0), dims=1)), padding)   # use _lat both times (_lon does it circularly)
    return img[rowrange, colrange]
end

function ypad(img, newheight)
    height, width = size(img)
    height == newheight && return img
    height > newheight && error("Can't pad an image to a smaller size.")
    firstrow = (newheight - height) ÷ 2
    newimg = fill(RGB{N0f8}(1.0,1.0,1.0), (newheight, width))
    newimg[firstrow:(firstrow + height - 1), :] = img
    return newimg
end

function makelegend(labels, colors; scale=2)
    i = .!isempty.(labels)
    labels = labels[i]     # don't plot legend entries with empty labels
    colors = colors[i]
    len = length(labels)
    markerpositions = Point2f0.(0, len:-1:1) .* scale
    textpositions = Point2f0.(0.4, (len:-1:1) .+ 0) .* scale
    fig = scatter(
        markerpositions,
        color = colors,
        marker = :rect,
        markersize = 64*scale,
        axis = (show_axis = false,)
    )
    annotations!(
        labels,
        textpositions,
        align = (:left, :center),
        textsize = 36*scale
    )
    scene = fig.axis.scene
    limits!(FRect(-1, 0.25, 4, len+0.5)*scale)     # allow for 10 lines (to get constant text size for any number of labels)
    filename = in_datafolder("output", "legend.png")
    isfile(filename) && rm(filename)
    Makie.save(filename, scene, resolution = scene.resolution.val .* scale)
    GLMakie.destroy!(GLMakie.global_gl_screen())
end

# first color in colorlist that is not in columncolors (which may have many repeated elements)
function firstcolor(columncolors, colorlist)
    colorcounts = zeros(Int, length(colorlist))
    for c in columncolors
        if c > 0
            colorcounts[colorlist[c]] += 1
        end
    end
    goodcolors = colorlist[colorcounts.==0]
    return isempty(goodcolors) ? 0 : rand(goodcolors)
end

function greedycolor(connected, colorlist, order; randseed=0)
    randseed > 0 && Random.seed!(randseed)
    n = size(connected, 1)
    colors = zeros(Int, n)
    for i in order
        columncolors = colors[connected[:,i]]
        col = firstcolor(columncolors, colorlist)
        if col == 0
            println("All colors used, retrying with randseed=$(randseed+1)...")
            return greedycolor(connected, colorlist, order; randseed=randseed+1)
        end
        colors[i] = col
    end
    # display(countmap(colors))
    return colors
end

# reverse coordinate order
lonlatpoint(latlon) = Point2f0(latlon[2], latlon[1])

# returns great circle track between points given as Point2f0 objects (lon, lat order, in degrees).
greatcircletrack(point1::AbstractArray, point2::AbstractArray, segments) = greatcircletrack(Tuple(point1), Tuple(point2), segments)

# returns great circle track between points given as (lon, lat) tuples (in degrees).
greatcircletrack(point1::Tuple, point2::Tuple, segments) = [reverse(intermediatepoint(point1, point2, f)) for f in LinRange(0, 1, segments)]

# intermediatepoint() ported from Javascript to Julia
# Source: https://www.movable-type.co.uk/scripts/latlong.html, original name 'intermediatePointTo()'
#
# Returns the point on a great circle at given fraction between point1 and point2.
#
# @param   {Tuple} point1 - (latitude,longitude) of starting point in degrees.
# @param   {Tuple} point2 - (latitude,longitude) of destination point in degrees.
# @param   {number} fraction - Fraction between the two points (0 = starting point, 1 = destination point).
# @returns {Tuple} Intermediate point between this point and destination point (latitude,longitude).
#
# @example
#   p1 = (52.205, 0.119)
#   p2 = (48.857, 2.351)
#   pInt = intermediatePointTo(p1, p2, 0.25)    # 51.3721°N, 0.7073°E
function intermediatepoint(point1::Tuple, point2::Tuple, fraction)
    φ1, λ1 = point1
    φ2, λ2 = point2

    Δφ = φ2 .- φ1   # distance between points
    Δλ = λ2 .- λ1   # distance between points
    a = sind(Δφ/2)*sind(Δφ/2) + cosd(φ1)*cosd(φ2)*sind(Δλ/2)*sind(Δλ/2)
    δ = 2 * atan(sqrt(a), sqrt(1-a))   # radians

    A = sin((1-fraction)*δ) / sin(δ)
    B = sin(fraction*δ) / sin(δ)

    x = A*cosd(φ1)*cosd(λ1) + B*cosd(φ2)*cosd(λ2)
    y = A*cosd(φ1)*sind(λ1) + B*cosd(φ2)*sind(λ2)
    z = A*sind(φ1) + B*sind(φ2)

    lat = atand(z, sqrt(x^2 + y^2))
    lon = atand(y, x)
    return (lat, lon)
end

plotmap(x; args...) = heatmap(reverse(x, dims=2); scale_plot=false, args...)

function plottraining(countryname)
    df, _ = loadtrainingdata()
    cc = df[:, :country] .== countryname
    normdemand = loaddemanddata()[cc, :normdemand]
    # normalize temp1 to the range [0.5, 1.5] for the 5%-95% quantiles
    normtemp1 = (df[cc,:temp1] .- df[cc,:temp1_qlow]) ./ (df[cc,:temp1_qhigh] .- df[cc,:temp1_qlow]) .+ 0.5
    t = 1:8760
    plot(t, normdemand, color=:black, limits=FRect(0,0,9000,50))
    plot!(t, normtemp1, color=:blue)
end

function plottimezones()
    tzindices, tznames = loadtimezones(1:36000, 1:18000)
    nreg = length(tznames)
    cc = connectedregions(tzindices, nreg)
    colorindices = greedycolor(cc, 1:12, 1:nreg, randseed=494)
    cs = colorschemes[:Set3_12].colors[colorindices]
    plotmap(tzindices[1:5:end, 1:5:end], colormap=cs)
end

function plottimeoffsets()
    tzindices, tznames = loadtimezones(1:36000, 1:18000)
    tz = [n[1:3] == "Etc" ? TimeZone(n, TimeZones.Class(:LEGACY)) : TimeZone(n) for n in tznames]
    offsetlist = tzoffset.(tz)
    offsetvalues = sort(unique(offsetlist))
    offsetdata = [offsetlist[tzi] for tzi in tzindices]
    offsetindices = indexin(offsetdata, offsetvalues)
    nreg = length(offsetvalues)
    cc = connectedregions(offsetindices, nreg)
    colorindices = greedycolor(cc, 1:12, 1:nreg, randseed=1)
    cs = colorschemes[:Set3_12].colors[colorindices]
    plotmap(offsetindices[1:5:end, 1:5:end], colormap=cs)
end

function testmap()
    lons = LinRange(1, 30, 360)
    lats = LinRange(51, 74, 180)

    field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

    source = LonLat()
    dest = WinkelTripel()

    xs, ys = xygrid(lons, lats)
    Proj4.transform!(source, dest, vec(xs), vec(ys))

    scene = surface(xs, ys; color = field, shading = false, show_axis = false, scale_plot = false)

    geoaxis!(scene, extrema(lons), extrema(lats); crs = (src = source, dest = dest,))

    # coastlines!(scene; crs = (src = source, dest = dest,))
end

max_stars_repo_name: euronion/GlobalEnergyGIS

max_issues_repo_name: euronion/GlobalEnergyGIS

max_forks_repo_name: euronion/GlobalEnergyGIS


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

module ServiceCommunityClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_BgpServiceCommunitiesApi.jl")

# export models
export convert, BGPCommunity
export convert, BgpServiceCommunity
export convert, BgpServiceCommunityListResult
export convert, BgpServiceCommunityPropertiesFormat
export convert, CloudError
export convert, CloudErrorBody
export convert, Resource

# export operations
export convert, BgpServiceCommunitiesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using FixedSizeArrays, ImmutableArrays

immutable BenchData <: FixedVector{Float64, 5}
	gc_bytes::Float64
	time_ns::Float64
	gc_time_ns::Float64
	gc_num_pause::Float64
	gc_num_full_sweep::Float64
end

function Base.show(io::IO, b::BenchData)
	print(io, """
		bytes allocated 	: $(b.gc_bytes)
		time_ns 			: $(b.time_ns)
		gc_time_ns 			: $(b.gc_time_ns)
		gc_num_pause 		: $(b.gc_num_pause)
		gc_num_full_sweep 	: $(b.gc_num_full_sweep)
		""")
end

macro timing(expr)
	tms 	= gensym()
	bts 	= gensym()
	gctms 	= gensym()
	gcps 	= gensym()
	gcswps 	= gensym()

	tms1 	= gensym()
	gcswps1 = gensym()
	gcps1 	= gensym()
	gctms1 	= gensym()
	bts1 	= gensym()
	esc(quote
		$bts 	 = Base.gc_bytes()
		$gctms 	 = Base.gc_time_ns()
		$gcps 	 = Base.gc_num_pause()
		$gcswps  = Base.gc_num_full_sweep()
		$tms 	 = Base.time_ns()

		$expr

		$tms1 	 = Base.time_ns()
		$gcswps1 = Base.gc_num_full_sweep()
		$gcps1 	 = Base.gc_num_pause()
		$gctms1  = Base.gc_time_ns()
		$bts1 	 = Base.gc_bytes()

		BenchData($bts1-$bts, $tms1-$tms, $gctms1-$gctms, $gcps1-$gcps, $gcswps1-$gcswps)
	end)
end
type Vec4{T} <: FixedVector{T, 4}
	x::T
	y::T
	z::T
	w::T
end

function test(N)
	bench = Dict(
		(:julia, :vector_creation) 		=> BenchData[],
		(:julia, :matrix_creation) 		=> BenchData[],
		(:julia, :vector_creation_rand) => BenchData[],
		(:julia, :matrix_creation_rand) => BenchData[],
		(:fsa, 	 :vector_creation) 		=> BenchData[],
		(:fsa, 	 :matrix_creation) 		=> BenchData[],

		(:fsa, 	 		:vector_creation_fromj) => BenchData[],
		(:fsa, 	 		:matrix_creation_fromj) => BenchData[],
		(:fsawrapper, 	:vector_creation) 		=> BenchData[],
		(:immutablearrays, 	:vector_creation) 		=> BenchData[],
	)
	for i=1:N
		push!(bench[(:julia, :vector_creation)], (@timing vecdata = [1,2,3,4]))
		push!(bench[(:julia, :matrix_creation)], (@timing matdata = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4]))

		push!(bench[(:julia, :vector_creation_rand)], (@timing vecrand = rand(4)))
		push!(bench[(:julia, :matrix_creation_rand)], (@timing matrand = rand(4,4)))
		nvec(1,2,3,4)
		nvec(1,2,3,4)

		push!(bench[(:fsa, :vector_creation)], (@timing vecrand = nvec(1,2,3,4)))
		push!(bench[(:fsa, :matrix_creation)], (@timing matrand = nvec((4,4), 1,2,3,4, 1,2,3,4, 1,2,3,4, 1,2,3,4)))

		push!(bench[(:fsa, :vector_creation_fromj)], (@timing vecrand = nvec(vecrand)))
		push!(bench[(:fsa, :matrix_creation_fromj)], (@timing matrand = nvec(matrand)))

		push!(bench[(:fsawrapper, :vector_creation)], (@timing vecrand = Point(1,2,3,4)))

		push!(bench[(:immutablearrays, :vector_creation)], (@timing vecrand = Vector4(1,2,3,4)))
	end
	bench
end
result = test(1000)

println(mean(result[(:fsa, :vector_creation)][100:end]))
println(mean(result[(:fsawrapper, :vector_creation)][100:end]))
println(mean(result[(:immutablearrays, :vector_creation)][100:end]))

#=
function test2(N)
	a,b,c,d = rand(Float64, 4)
	result =  LolVec(a,b,c,d)
	bench = 0.0
	for i=1:N
		a,b,c,d = rand(Float64, 4)
		tic()
		result = LolVec(a,b,c,d)
		b = toq()
		bench += b
	end

	result,bench
end

@time t= test2(10)
@time t,b= test2(10^4)
println(b)
#0.00622337 seconds (1 MB allocated)
#0.001912617 seconds (470 kB allocated)

#0.0036441109999999725 2mb
#0.0010228019999999729 1mb
#0.0012369930000000022 1mb
=#


text:
# Typical usage: run `generate_artifacts.jl . --deploy` from within the directory that contains the files
using Pkg.Artifacts, LibGit2, ghr_jll

deploy = true
if deploy && !haskey(ENV, "GITHUB_TOKEN")
    error("For automatic github deployment, export GITHUB_TOKEN!")
end

data_dir = @__DIR__

if deploy
    # Where we will put our tarballs
    tempdir = mktempdir()

    function get_git_remote_url(repo_path::String)
        repo = LibGit2.GitRepo(repo_path)
        origin = LibGit2.get(LibGit2.GitRemote, repo, "github")
        return LibGit2.url(origin)
    end

    # Try to detect where we should upload these weights to (or just override
    # as shown in the commented-out line)
    origin_url = get_git_remote_url(dirname(@__DIR__))
    deploy_repo = "$(basename(dirname(origin_url)))/$(splitext(basename(origin_url))[1])"

    # deploy_repo = "staticfloat/ObjectDetector.jl"
    tag = "v0.0.2"
end

jsonfiles = filter(x -> endswith(x, ".json"), readdir(data_dir))

for jsonfile in jsonfiles
    name = splitext(jsonfile)[1]
    @info("Generating artifact for $(jsonfile)")
    # Create a local artifact
    hash = create_artifact() do artifact_dir
        # Copy in weights
        cp(joinpath(data_dir, jsonfile), joinpath(artifact_dir, jsonfile))
    end

    # Spit tarballs to be hosted out to local temporary directory:
    if deploy
        tarball_hash = archive_artifact(hash, joinpath(tempdir, "$(name).tar.gz"))

        # Calculate tarball url
        tarball_url = "https://github.com/$(deploy_repo)/releases/download/$(tag)/$(name).tar.gz"

        # Bind this to an Artifacts.toml file
        @info("Binding $(name) in Artifacts.toml...")
        bind_artifact!(joinpath(@__DIR__, "..", "Artifacts.toml"), name, hash; download_info=[(tarball_url, tarball_hash)], lazy=true, force=true)
    end
end

if deploy
    # Upload tarballs to a special github release
    @info("Uploading tarballs to $(deploy_repo) tag `$(tag)`")
    ghr() do ghr_exe
        # This may need the -u flag for the github user.
        run(`$ghr_exe  -r $(basename(deploy_repo)) $(tag) $(tempdir)`)
    end

    @info("Artifacts.toml file now contains all bound artifact names")
end

max_stars_repo_name: rustyconover/muid-julia

max_issues_repo_name: rustyconover/muid-julia

max_forks_repo_name: rustyconover/MemorableUniqueIdentifier


text:
#Credits: this code, and its embedded documentation, was ported to Julia from the original Matlab implementation of [EpiGen](https://github.com/LevineLab/EpiGen).
#Please refer to [their paper](https://www.biorxiv.org/content/10.1101/637272v1) if you use this model. And if in doubt then please use the original code.

using StatsBase, Distributions, DataFrames, Dates, Random

using CairoMakie

"""
    one_plot(results)

Plot some of the output from `EpiGen`.
"""
function one_plot(results)
    f = Figure()
    ax = Axis(f[1,1], xlabel="time",ylabel="various units")
    lines!(ax,results[:,2],label="fitness")
    lines!(ax,results[:,5],label="genetic distance")
    lines!(ax,results[:,6],label="epigenetic distance")
    f[1, 2] = Legend(f, ax, framevisible = false)
    f
end

"""
    EpiGen(N=500,time=200)

Run the EpiGen model from [their paper](https://www.biorxiv.org/content/10.1101/637272v1)

```
results=EpiGen(500,200)
one_plot(results)
```
"""
function EpiGen(N=500,time=200)

#N=1000 #the number of individuals in the population
#time=1000 #the number of generations to run the simulation
t_split = 40 #if patch='TRUE', t_split is the number of generations (40 in this case) the population stays in one environment or the other
radius=1 #distance from the hypersphere origin
nugen = 10 #number of genetic (HT) mutations for each generation
nuepi = 90 #number of epigenetic (LT) mutations for each generation
epiback = 0.0205 #epigenetic reversion rate
mlimitgen=2 #maximum effect a genetic mutation can have
mlimitepi = 0.3 #maximum effect an epigenetic mutation can have
epigenetics = true #switch to toggle effects of epigenetics on or off
patch = true #switch to determine if population switches between one environment and another. If TRUE, population goes in and out of selection environment at a frequency of t_split generations. If FALSE, only stays in the selection environment.
mode = "negative" #When patch = 'TRUE', this will run the model in two modes that determines the selection rules for environment 2: random or negative. 'negative' turns on randomly sampling individuals weighted by the reciprocal number of genetic mutations when in environment 2 (stabilizing selection). 'random' turns on randomly sampling indviduals in environment 2 with no weighting. Random sampling weighted by fitness always occurrs in environment 1
runnumber = 1 #replicate number of model run

patch_intervals = fill(0.0, t_split) #Create cell array of matrices to store interval runs
patch_removed_intervals = fill(0.0, t_split) #Create cell array of matrices to store interval runs of removed individuals. These are individuals that were not sampled to go into the next generation.
results_mat = fill(0.0, t_split) #Create cell array of matrices to store results matrix runs
results_removed_mat = fill(0.0, t_split) #Create cell array of matrices to store results matrix runs of removed individuals
epi_mat = fill(0.0, t_split) #Create cell array of matrices to store epigenetic distribution run
gen_mat = fill(0.0, t_split) #Create cell array of matrices to store genetic distribution run
epi_removed_mat = fill(0.0, t_split) #Create cell array of matrices to store epigenetic distribution run
gen_removed_mat = fill(0.0, t_split) #Create cell array of matrices to store genetic distribution run

##

"""
Driver is a vector of 0's and 1's and has length = time in generations. The loop below will
step through every element of Driver. When patch = 'TRUE' and Driver == 1, then select
= 1, and the population is randomly sampled to the next generation weighted
by fitness, i.e. is in environment 1. When patch = 'TRUE', mode = 'negative', and
Driver == 0, then the population is randomly sampled to the next generation weighted
by the reciprocal of the number of genetic mutations (stabilizing selection),
i.e. is in environment 2.
"""

nt=time+t_split-1
Driver=fill(1,nt)
for q in t_split:t_split*2:time
    Driver[q:q+t_split-1].=0
end

##

#Population information
population_co=fill(0.0,(N,25)) #create populations
population_co[:,1] .= 1:N #Index for individuals
population_co[:,2] .=1 #initial distance from the optimum is radius
population_co[:,3] .=exp.((-population_co[:,2].^2)./2) #Fitness - Using a Gaussian fitness function
#Column 4 is for distance traveled by genetic mutation
#Column 5 is for number of genetic mutations
#Column 6 is last time point of genetic mutations
#Column 7 is for distance traveled by epigenetic mutations
#Column 8 is for number of epigenetic mutations
#Column 9 is last time point of epigenetic mutations

#Information of individuals in the population that were removed from the main reproducing population
population_removed=zeros(N,25) ##create populations
#Column 1 #Index for individuals
#Column 2 #initial distance from the optimum is radius
#Column 3 #fitness function
#Column 4 is for distance traveled by genetic mutation
#Column 5 is for number of genetic mutations
#Column 6 is last time point of genetic mutations
#Column 7 is for distance traveled by epigenetic mutations
#Column 8 is for number of epigenetic mutations
#Column 9 is last time point of epigenetic mutations

#matrix to populate in loop
population_disgen=fill(0.0,(N,nt)) #distribution of genetic mutations of kept population (i.e. sampled to go onto next generation)
population_disepi=fill(0.0,(N,nt)) #distribution of epigenetic mutations of kept population
population_removed_disgen = fill(0.0,(N,nt)) #distribution of genetic mutations of removed population (i.e. not sampled to go onto next generation)
population_removed_disepi = fill(0.0,(N,nt)) #distribution of epigenetic mutations of removed population

##

#Initialize geometric model
#dimensions of spherical phenotypic space
#Z-values have been solved manually, dimensions solved for are 5, 10, 15,
#20, 25, 30, 35, 40
#If some other value of n is needed, need to solve the integral for Z and add the value to Z.mat
zmat=[5 0.75; 10 1.1641; 15 1.4663; 20 1.7162; 25 1.9342; 30 2.1299; 35 2.3092; 40 2.4755];

#Store Z-value for 5 dimension. We use 5 dimensions for all model runs
Z5=zmat[1,2]

#the angle between the mutational vector and the vector running from the current phenotype to the origin
#Angles can take on values between -pi/2 to pi/2
phi_values=collect(range(-pi/2,pi/2,50000))

#Calculate the n = 5 probability density function for the random angle
prob_density5=Z5*cos.(phi_values) .^(zmat[1,1] - 2)

#Initialize Results of kept population
#Column 1 generations
#Column 2 mean population fitness
#Column 3 mean number of genetic mutations
#Column 4 mean number of epigenetic mutations
#Column 5 mean distance traveled by genetic mutations
#Column 6 mean distance traveled by epigenetic mutations
results_co=zeros(nt,28)
results_co[:,1] .= 1:nt
#colname={'time','meanfitness','nogenmut','noepimut','distgen','distepi','net_fi_carbon_gen','net_fi_p_gen','net_fi_i_gen','net_fi_carbon_epi','net_fi_p_epi','net_fi_i_epi'};
#colindex={1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};

#Initialize Results of removed population (i.e. individuals not sampled to
#next generation
results_removed=zeros(nt,28)
results_removed[:,1] .=1:nt

##

#select = 1 #Switch to randomly resample population either not weighted (0) or weighted (1) by fitness

####
#Looping over time
for i in 1:length(Driver)

    #Check if patch is TRUE or FALSE. If patch_tf == 1, run model with only
    #1 type of selection (i.e. 1 environment) where sampling is weighted by
    #fitness. If patch_ft == 0, run the model for two types of selection (i.e. two
    #environments). When Driver = 1, sample weighted by fitness. When
    #Driver = 0, sample weighted by the reciprocal number of genetic
    #mutations
    patch_tf = !patch

    if patch_tf
        select = 1 #only sample with weighting as fitness
    elseif (!patch_tf) && (Driver[i] == 1)
        select = 1 #resample population with weight by fitness
    elseif (!patch_tf) && (Driver[i] == 0)
        select = 0 #resample population with other weighting parameters (e.g., by the reciprocal number of genetic mutations)
    end

    ####
    #Make genetic mutations, mutational supply is set to nugen/1 unit time
    #Distances from optimum are calculated individually

    seed=1000*Dates.second(now())+Dates.millisecond(now())
    Random.seed!(seed)

    #Vector of genetic mutations with each value being the magnitude of
    #their effect
    mut_effect=0 .+(mlimitgen*radius).*rand(1,nugen)

    seed=1000*Dates.second(now())+Dates.millisecond(now())
    Random.seed!(2*seed)

    #Randomly sample 10 angles
    phi5=sample(phi_values, ProbabilityWeights(prob_density5), nugen; replace=true)

    ###Assign mutations to individuals###

    index_mut=1:N

    seed=1000*Dates.second(now())+Dates.millisecond(now())
    Random.seed!(3*seed)

    #Randomly sample 10 individuals to genetically mutate without replacement (multiple hits are not allowed)
    index_mutated=sample(index_mut,nugen; replace=false)

    #retrieve mutants from population
    mutated_orgs =population_co[index_mutated,:]
    #mutated_orgs =dropdims(mutated_orgs,dims=1)

    ###calculate distance traveled after genetic mutation
    dist_travelled5 = mutated_orgs[:,2] - sqrt.(mutated_orgs[:,2].^2 + mut_effect[:].^2 .+ 2.0 *mutated_orgs[:,2].*mut_effect[:].*sin.(phi5[:]))

    #Calculate distance from the optimum phenotype
    dist5 = sqrt.(mutated_orgs[:,2].^2 + mut_effect[:].^2 + 2.0*sin.(phi5[:]).*mutated_orgs[:,2].*mut_effect[:])

    ####
    #Distance from the optimum phenotype based on oligo conditions;
    #oligo == 1 ; n = 5 dimensions of the hypersphere

    #store distance from optimum for mutated individuals
    mutated_orgs[:,2] .= dist5

    #rename distance_travelled5 vector to dist_travelled_gen
    dist_travelled_gen = dist_travelled5

    #NOT USED IN THIS VERSION OF MODEL CALCULATIONS - IGNORE
    #Using law of sines, calculate alternative angle.
    fi5 = asin.((sin.(phi5[:] .+ (pi/2)).*(mut_effect[:]))./dist5)

    #Insert fi5
    mutated_orgs[:,10] .= mutated_orgs[:,10] .+ fi5

    #Save mean distance traveled toward optimum by genetic mutation
    results_co[i,25] = mean(mutated_orgs[:,2])
    #Save std of distance traveled toward optimum by genetic mutation
    results_co[i,26] = std(mutated_orgs[:,2])
    #Calculates fitness for the mutants after genetic mutation
    mutated_orgs[:,3] .= exp.((-mutated_orgs[:,2].^2)/2)
    ###Cumulative distance travelled by genetic mutations
    mutated_orgs[:,4] .= mutated_orgs[:,4] + dist_travelled_gen
    #Increase the number of genetic mutations by one
    mutated_orgs[:,5] .= mutated_orgs[:,5] .+ 1
    #save distance travelled after genetic mutations
    mutated_orgs[:,6] .= dist_travelled_gen
    #Storing the mutated individuals
    population_co[index_mutated,:].=mutated_orgs[:,:]
    #Storing distribution of genetic mutations
    for g=1:length(dist_travelled_gen)
        population_disgen[index_mutated[g], i] = dist_travelled_gen[g]
    end

    #If the model run includes epigenetic effects, run the following
    tf = epigenetics

    if tf == 1

        #Make epigenetic mutations, mutational supply is set to Nu.epi / 1 unit of time

        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(4*seed)

        #Mutational effects of epigenetics
        mut_effect_epi=(mlimitepi*radius).*rand(1, nuepi)

        #Random angle for epigenetic mutations

        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(5*seed)

        #Randomly draw random angles for epigenetic mutations
        phi_epi5=sample(phi_values, ProbabilityWeights(prob_density5), nuepi; replace=true)

        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(6*seed)

        #retrieve epimutants from population w/out replacement
        index_mutated_epi = sample(index_mut, nuepi, replace=false)
        mutated_orgs_epi = population_co[index_mutated_epi, :]

        ###distance traveled after epimutations
        dist_travelled_epi5 = mutated_orgs_epi[:,2] - sqrt.(mutated_orgs_epi[:,2].^2 + mut_effect_epi[:].^2 + 2.0*mutated_orgs_epi[:,2].*mut_effect_epi[:].*sin.(phi_epi5[:]))

        #Distance from the optimum phenotype epimutations
        dist_epi5 = sqrt.(mutated_orgs_epi[:,2].^2 + mut_effect_epi[:].^2 + 2.0*sin.(phi_epi5[:]).*mutated_orgs_epi[:,2].*mut_effect_epi[:])

        #Distance from the optimum phenotype based on oligo conditions;

        #Store distance from the optimum phenotype
        mutated_orgs_epi[:,2] .= dist_epi5

        #rename dist_travelled_epi5 to dist_travelled_e
        dist_travelled_e = dist_travelled_epi5;

        #Save mean distance travelled from optimum by epigenetic mutation
        results_co[i,27] = mean(mutated_orgs_epi[:,2])

        #Save std of distance travelled from optimum by epigenetic mutation
        results_co[i,28] = std(mutated_orgs_epi[:,2])

        #Calculate fitness for mutants
        mutated_orgs_epi[:,3] .= exp.((-mutated_orgs_epi[:,2].^2)/2)

        #Cumulative distance travelled by epigenetic mutations
        mutated_orgs_epi[:,7] .= mutated_orgs_epi[:,7] + dist_travelled_e

        #Increase the number of epigenetic mutations by one
        mutated_orgs_epi[:,8] .= mutated_orgs_epi[:,8] .+ 1

        #Store distance travelled
        mutated_orgs_epi[:,9] .= dist_travelled_e

        #CALCULATE ALTERNATIVE ANGLE - IGNORE THIS METRIC - NOT USED IN THIS MODEL VERSION
        #Using law of sines
        #fi_epi5 = asin.((sin.(phi_epi5[:] /+ pi/2).*(mut_effect_epi[:]))./dist_epi5[:])

        #Insert fi5 for epi
        #mutated_orgs_epi[:,18] = mutated_orgs_epi[:,18] + fi_epi5

        #Storing mutated individuals
        population_co[index_mutated_epi,:] .= mutated_orgs_epi

        #Distribution of epigenetic mutations
        for e=1:length(dist_travelled_e)
            population_disepi[index_mutated_epi[e], i] = dist_travelled_e[e]
        end

    end

    #As of now, multiple hits of either genetic and epigenetic mutations are not possible.
    #However, a single individual can have two hits if one is genetic and one epigenetic mutation.

    ######
    #Population regulation and reproduction

    #Soft selection is acting -> population size stays the same
    #Sampling individuals to the next time point with replacement

    ofsp_index = 1:N;

    if select == 0 #if true, population is in environment 2

        mode_type = (mode=="random") #If true, sample without weighting
        #Sampling the next generation, no weighting
        if mode_type == 1

        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(7*seed)

        sampled_index = sample(ofsp_index, N, replace=true)
        else
        #Sampling the next generation, negative weighting (i.e. the
        #reciprocal number of genetic mutations)
        reciprocal_weights = 1 ./ (0.1 .+ population_co[:,5])

        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(8*seed)

        sampled_index = sample(ofsp_index, ProbabilityWeights(reciprocal_weights), N; replace=true)

        end

    elseif select == 1 #if true, population is in environment 1

        #Sampling the next generation, weighting with fitness
        seed=1000*Dates.second(now())+Dates.millisecond(now())
        Random.seed!(9*seed)

        sampled_index = sample(ofsp_index, ProbabilityWeights(population_co[:,3]), N; replace=true)

    end

    #Population after reproduction and regulation
    population_co .= population_co[sampled_index,:]
    population_disgen .= population_disgen[sampled_index,:]
    population_disepi .= population_disepi[sampled_index,:]

    #Removed individuals from population
    pop_removed_index = setdiff(ofsp_index, sampled_index)
    population_removed = population_co[pop_removed_index,:]

    population_removed_disgen = population_disgen[pop_removed_index,:]
    population_removed_disepi = population_disepi[pop_removed_index,:]

    ######Store results removed population

    #Mean fitness of the population
    results_removed[i,2]= mean(population_removed[:,3])
    #Store the mean number of genetic mutations present in an individual
    results_removed[i,3]= mean(population_removed[:,5])
    #Store the standard deviation of mean genetic mutations present in an individual
    results_removed[i,23]= std(population_removed[:,5])
    #Store the mean number of epigenetic mutations present in an individual
    results_removed[i,4]= mean(population_removed[:,8])
    #Store the standard deviation of mean epigenetic mutations present in an individual
    results_removed[i,24]= std(population_removed[:,8])
    #Store the cumulative distance travelled towards the optimum by genetic mutations
    results_removed[i,5]= mean(population_removed[:,4])
    #Store the cumulative distance travelled towards the optimum by epigenetic mutations
    results_removed[i,6]= mean(population_removed[:,7])

    ######Store results sampled (non-removed) population

    #Mean fitness of the population
    results_co[i,2]= mean(population_co[:,3])
    #Store the mean number of genetic mutations present in an individual
    results_co[i,3]= mean(population_co[:,5])
    #Store the standard deviation of mean genetic mutations present in an individual
    results_co[i,23]= std(population_co[:,5])
    #Store the mean number of epigenetic mutations present in an individual
    results_co[i,4]= mean(population_co[:,8])
    #Store the standard deviation of mean epigenetic mutations present in an individual
    results_co[i,24]= std(population_co[:,8])
    #Store the cumulative distance travelled towards the optimum by genetic mutations
    results_co[i,5]= mean(population_co[:,4])
    #Store the cumulative distance travelled towards the optimum by epigenetic mutations
    results_co[i,6]= mean(population_co[:,7])
    #Store mean genetic mutation fi5
    results_co[i,7] = mean(population_co[:,10])
    #Store mean epigenetic mutation fi5
    results_co[i,15] = mean(population_co[:,18])

    ######
    #Loss of epigenetic effects
    #Next epigenetic changes are lost at a given rate
    #(This needs to be slower than the forward mutation rate of epigenetic effects)
    #Epigenetic changes revert independently (always at a probability of: epi.backmutation in each generation)

    #Number of epigenetic mutations present
    no_epi_mut = sum(population_co[:,8])

    if (no_epi_mut>0)
        #Check which individuals have epigenetic mutation
        epi_p_index = findall(population_co[:,8] .> 0)
        epi_present = population_co[epi_p_index,:]
        epi_present_gen = population_disgen[epi_p_index,:];
        epi_present_epi = population_disepi[epi_p_index,:];

        #Generate backmutations probabilistically
        n=Int(sum(epi_present[:,8]))
        if n > 0
            p=Binomial(1,epiback)
            epi_backmuts = rand(p,n)
        else
            epi_backmuts = fill(0.0,n)
        end

        #Check for each backmutation whether reversion occurs or not
        ###Mapping from individual mutations to individuals

        #Construct a helping vector with as many factors as there are inds with the equal
        #length of mutation vector

        help_vec = vcat(fill.( 1:length(epi_present[:,8]) ,  Int.(epi_present[:,8]) )...)

        ####
        #Sum by each factor
        ##Same function as "tapply" in R

        df = DataFrame(a=epi_backmuts, b=help_vec)
        gdf = groupby(df, :b)
        tmp = combine(gdf, :a => sum)
        list_of_reversions = tmp.a_sum

        #Logical vector if backmutation occurs in a given individual (TRUE) or not
        epi_back_mut_index = findall(list_of_reversions .> 0)

        #Take only those individuals that had backmutations
        epi_back_mut_ind = epi_present[epi_back_mut_index,:]
        epi_back_mut_ind_gen = epi_present_gen[epi_back_mut_index,:]
        epi_back_mut_ind_epi = epi_present_epi[epi_back_mut_index,:]

        #How many backmutations occurred in each individual that had backmutations
        number_of_reversions = list_of_reversions[epi_back_mut_index]

        #if any backmutations occurred
        if sum(list_of_reversions.>0) > 0

            seed=1000*Dates.second(now())+Dates.millisecond(now())
            Random.seed!(10*seed)

            for b=1:length(number_of_reversions)

                lost_mutation=zeros(1,number_of_reversions[b])

                current_org = [epi_back_mut_ind_epi[b]]
                current_org = current_org[findall((!).(current_org .== 0))]
                #randomly sample epi effects from the distribution
                if length(current_org) >= number_of_reversions[b]
                    lost_mutation[1,:]=sample(current_org, number_of_reversions[b],replace=false)
                end

                #Now individuals need to lose mutation
                if !isempty(lost_mutation)

                    #loop through all epimutations
                    for count=1:length(lost_mutation)
                        #ensure current org contains epimutations to be
                        #lost
                        logic = sum(current_org .== lost_mutation[1,count])

                        if logic == 1
                            #gather epimutation values to be lost
                            r = epi_back_mut_ind_epi[b,:]

                            #find indicies of lost mutations
                            idx2 = findall(r.==lost_mutation[1,count])

                            #if only 1 mutation is found,
                            #set to zero
                            if length(idx2) == 1
                                epi_back_mut_ind_epi[b,idx2] .= 0

                            #if more than 1 epimutation, loop through the
                            #number of epimutations to lose
                            elseif length(idx2) > 1
                                [epi_back_mut_ind_epi[b,idx2[l]] = 0 for l=1:length(idx2)]
                            end
                        end
                    end
                end

               #Subtract the effect of the mutations from distance travelled by epigenetics
               epi_back_mut_ind[b,7]= epi_back_mut_ind[b,7]-sum(lost_mutation[1,:])
               #Add the effect of the mutation from distance to optimum
               epi_back_mut_ind[b,2]= epi_back_mut_ind[b,2]+sum(lost_mutation[1,:])
               #Subtract the number of epigenetic mutations
               epi_back_mut_ind[b,8]=epi_back_mut_ind[b,8]-number_of_reversions[b]
               #Recalculate fitness for backmutated individuals
               epi_back_mut_ind[b,3] = exp((-epi_back_mut_ind[b,2].^2)/2)
            end

            #Replace backmutated individuals in the population_co dataframe
            epi_present[epi_back_mut_index,:]=epi_back_mut_ind
            population_co[epi_p_index,:].=epi_present

            #Replace the epigenetic mutational distributions of the
            #backmutated individuals in the population_disepi dataframe
            epi_present_epi[epi_back_mut_index,:]=epi_back_mut_ind_epi
            population_disepi[epi_p_index,:]=epi_present_epi
        end
    end
end

return results_co
end #function EpiGen()

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
using Test

using AIMASamples.Vacuum
@testset "Vacuum" begin
    @test RunTableDrivenVacuumAgent()  == RunTableDrivenVacuumAgentResult
    @test RunSimpleReflexVacuumAgent() == RunSimpleReflexVacuumAgentResult
    @test RunReflexVacuumAgent()       == RunReflexVacuumAgentResult
    @test RunModelBasedVacuumAgent()   == RunModelBasedVacuumAgentResult
end

using AIMASamples.Romania
@testset "Romania" begin
   @test solveRomanianMapProblemBFS()  == solveRomanianMapProblemResultBFS
   @test solveRomanianMapProblemUCS()  == solveRomanianMapProblemResultMinCost
   @test solveRomanianMapProblemDLS()  == solveRomanianMapProblemResultDLS
   @test solveRomanianMapProblemIDS()  == solveRomanianMapProblemResultIDS
   @test solveRomanianMapProblemGSB()  == solveRomanianMapProblemResultGSB
   @test solveRomanianMapProblemGSD()  == solveRomanianMapProblemResultGSD
   @test solveRomanianMapProblemGSU()  == solveRomanianMapProblemResultGSU
   @test solveRomanianMapProblemGSAS() == solveRomanianMapProblemResultGSAS
   @test solveRomanianMapProblemGSBF() == solveRomanianMapProblemResultGSBF
   @test solveRomanianMapProblemRBF()  == solveRomanianMapProblemResultRBF
end

using AIMASamples.SlideBlock
@testset "SlideBlock" begin
    @test begin
        ret = solveSlideBlockProblemRBFS()
        println("$(ret[1]) to\n\n$(ret[end])")   #Added for coverage only
        length(ret) == 27
    end
    @test length(solveSlideBlockProblemGSU())  == 27
    @test length(solveSlideBlockProblemGSBF()) != 27
    @test length(solveSlideBlockProblemGSBF2())!= 27
    @test length(solveSlideBlockProblemGSAS()) == 27
    @test length(solveSlideBlockProblemGSAS2())== 27
end

using AIMASamples.Queens
@testset "Queens" begin
    @test begin
        ret = solveNQueensProblemGSU()
        println("$(ret[1]) to\n\n$(ret[end])")   #Added for coverage only
        length(ret) == 9
    end
    @test length(solveNQueensProblemGSBF()) == 9
    @test length(solveNQueensProblemGSAS()) == 9
    @test length(solveNQueensProblemRBFS()) == 9
    @test begin #This is carried out for coverage only
        solveNQueensProblemHCS()
        solveNQueensProblemSAS()
        solveNQueensProblemGAS()
        solveNQueensProblemGAS9()
        1 == 1
    end
end


text:
# NOTE: This simulation fits two models, and computes the posterior probability
# of β | y.

ENV["GKSwstype"] = "nul"  # For StatsPlots
println("pid: ", getpid())

# Read command line args.
if length(ARGS) > 1
  resultsdir = ARGS[1]
  awsbucket = ARGS[2]
  snum = parse(Int, ARGS[3])
  println("ARGS: ", ARGS)
else
  resultsdir = "results/test/"
  awsbucket = nothing
  snum = 3
end
flush(stdout)

# Create results/images dir if needed
imgdir = "$(resultsdir)/img"
mkpath(imgdir)

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))

using CytofDensityEstimation
const CDE = CytofDensityEstimation
using LaTeXStrings
using Distributions
using StatsPlots
using HypothesisTests
import Random
using RCall
using BSON
import CytofDensityEstimation.Model.default_ygrid
using StatsFuns

include("scenarios.jl")

plotsize = (400, 400)

# Simulate data.
Random.seed!(5);
simdata = scenarios(snum)
yC, yT = simdata[:yC], simdata[:yT]

# Plot true data density.
function plot_trude_data_density()
  plot(default_ygrid(), pdf.(simdata[:mmC], default_ygrid()), lw=3,
       label=L"y_C", color=:blue)
  plot!(default_ygrid(), pdf.(simdata[:mmT], default_ygrid()), lw=3,
        label=L"y_T", color=:red)
end
plot_trude_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-true-density.pdf"))
closeall()

# Plot observed data density.
function plot_observed_data_density()
  density(yC[isfinite.(yC)],  lw=3, label=L"y_C", color=:blue)
  density!(yT[isfinite.(yT)], lw=3, label=L"y_T", color=:red)
end
plot_observed_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-kde.pdf"))
closeall()

# Define data, prior, and initial state.
Random.seed!(7)
K = 5 # 3
data = CDE.Model.Data(yC, yT)
prior_mu = let
  yfinite = [data.yC_finite; data.yT_finite]
  Normal(mean(yfinite), std(yfinite))
end
prior = CDE.Model.Prior(K, mu=prior_mu, nu=LogNormal(1.6, 0.4),
                        p=Beta(1, 99), omega=InverseGamma(.1, .1),
                        eta=Dirichlet(K, 1 / K), psi=Normal(-1, 1))
state = CDE.Model.State(data, prior)
tuners = CDE.Model.Tuners(K, 0.1)

println("Priors:")
for fn in fieldnames(CDE.Model.Prior)
  println(fn, " => ", getfield(prior, fn))
end

function run(; p, beta)
  # Fix these parameters.
  state.p = p
  state.beta = beta

  # Run chain.
  @time chain, laststate, summarystats = CDE.fit(
      state, data, prior, deepcopy(tuners), nsamps=[3000], nburn=3000, thin=1,
      fix=[:p, :beta], rep_aux=10)

  # Make directoary for results if needed.
  resdir = "$(resultsdir)/m$(beta)"
  mkpath(resdir)

  # Save results
  BSON.bson("$(resdir)/results.bson",
            Dict(:chain=>deepcopy(chain), :laststate=>deepcopy(laststate),
                 :summarystats=>deepcopy(summarystats),
                 :simdata=>deepcopy(simdata)))

  return chain, laststate, summarystats
end

# Load rstats.
@rimport stats as rstats

function postprocess(beta)
  # Load results.
  out = BSON.load("$(resultsdir)/m$(beta)/results.bson")

  # Print summary statistics.
  CDE.Model.printsummary(out[:chain], out[:summarystats])
  flush(stdout)

  # Make imgdir if needed.
  _imgdir = joinpath(resultsdir, "m$(beta)/img")
  mkpath(_imgdir)

  # Plot results.
  @time CDE.Model.plotpostsummary(out[:chain], out[:summarystats], yC, yT,
                                  _imgdir, simdata=out[:simdata],
                                  bw_postpred=.15, xlims_=(-6, 6))
  flush(stdout)
end

# Run chains.
chain0, laststate0, summarystats0 = run(p=0.5, beta=0)
chain1, laststate1, summarystats1 = run(p=0.5, beta=1)

# Post process.
postprocess(0)
postprocess(1)

# Print KS statistic.
ks_fit = rstats.ks_test(yC, yT)
println("KS p-value: ", ks_fit["p.value"][1])
flush(stdout)

# Bayes Factor
@time lbf = CDE.Model.log_bayes_factor(data, chain0[1], chain1[1])
pm1 = CDE.Model.posterior_prob1(lbf, logpriorodds=logit(state.p))
println("Log Bayes Factor: $(lbf) | P(M1|y): $(pm1)")

# Plot marginal posterior density
CDE.Model.plot_posterior_predictive(yC, yT, chain0, chain1, pm1, imgdir,
                                    bw_postpred=0.15, simdata=simdata, lw=.5,
                                    ls=:solid, xlims_=(-6, 6), digits=5,
                                    fontsize=7)
CDE.Model.plot_gamma(yC, yT, chain0, chain1, pm1, imgdir)

# DIC
dic0, dic1 = CDE.dic(chain0[1], data), CDE.dic(chain1[1], data)
println("(DIC0, DIC1): ($(round(dic0, digits=3)), $(round(dic1, digits=3)))")
flush(stdout)

# Send results to S3.
if awsbucket != nothing
  CDE.Util.s3sync(from=resultsdir, to=awsbucket, tags=`--exclude '*.nfs'`)
end

println("Done!")
flush(stdout)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module KoalaRidge

# new:
export RidgeRegressor # eg, RidgeRegressor

# needed in this module:
import Koala: Regressor, BaseType, SupervisedMachine, Transformer
import Koala: params, keys_ordered_by_values
import DataFrames: AbstractDataFrame, DataFrame
import KoalaTransforms
import MultivariateStats
import UnicodePlots

# to be extended (but not explicitly rexported):
import Koala: setup, fit, predict
import Koala: default_transformer_X, default_transformer_y, transform, inverse_transform

## Model type definitions

struct LinearPredictor <: BaseType
    coefficients::Vector{Float64}
    bias::Float64
end

# Following returns a `DataFrame` with three columns:
#
# column name | description
# :-----------|:-------------------------------------------------
# `:index`    | index of a feature used to train `predictor`
# `:feature`  | corresponding feature label provided by `features`
# `:coef`     | coefficient for that feature in the predictor
#
# The rows are ordered by the absolute value of the coefficients. If
# `rgs` is unfitted, an error is returned.
function coef_info(predictor::LinearPredictor, features)
    coef_given_index = Dict{Int, Float64}()
    abs_coef_given_index = Dict{Int, Float64}()
    v = predictor.coefficients
    for k in eachindex(v)
        coef_given_index[k] = v[k]
        abs_coef_given_index[k] = abs(v[k])
    end
    df = DataFrame()
    df[:index] = reverse(keys_ordered_by_values(abs_coef_given_index))
    df[:feature] = map(df[:index]) do index
        features[index]
    end
    df[:coef] = map(df[:index]) do index
        coef_given_index[index]
    end
    return df
end

mutable struct RidgeRegressor <: Regressor{LinearPredictor}
    lambda::Float64
end

# lazy keywork constructor
RidgeRegressor(; lambda=0.0) = RidgeRegressor(lambda)

default_transformer_X(model::RidgeRegressor) =
    KoalaTransforms.DataFrameToArrayTransformer()
default_transformer_y(model::RidgeRegressor) =
    KoalaTransforms.RegressionTargetTransformer()

struct Cache <: BaseType
    X::Matrix{Float64}
    y::Vector{Float64}
    features::Vector{Symbol}
end

setup(model::RidgeRegressor, X, y, scheme_X, parallel, verbosity) =
    Cache(X, y, scheme_X.spawned_features)

function fit(model::RidgeRegressor,
             cache, add, parallel, verbosity; args...)

    weights = MultivariateStats.ridge(cache.X, cache.y, model.lambda)

    coefficients = weights[1:end-1]
    bias = weights[end]

    predictor = LinearPredictor(coefficients, bias)

    # report on the relative strength of each feature in the predictor:
    report = Dict{Symbol, Any}()
    cinfo = coef_info(predictor, cache.features) # a DataFrame object
    u = Symbol[]
    v = Float64[]
    for i in 1:size(cinfo, 1)
        feature, coef = (cinfo[i, :feature], cinfo[i, :coef])
        coef = floor(1000*coef)/1000
        if coef < 0
            label = string(feature, " (-)")
        else
            label = string(feature, " (+)")
        end
        push!(u, label)
        push!(v, abs(coef))
    end
    report[:feature_importance_curve] = u, v

    return predictor, report, cache

end

function predict(predictor::LinearPredictor, pattern::Vector{Float64})

    ret = predictor.bias
    for i in eachindex(predictor.coefficients)
        ret += predictor.coefficients[i]*pattern[i]
    end

    return ret

end

predict(model::RidgeRegressor, predictor::LinearPredictor,
        X::Matrix{Float64}, parallel, verbosity) =
            Float64[predict(predictor, X[i,:]) for i in 1:size(X,1)]

end # of module


text:
#== import modules here ==#
using PyCall
bse = pyimport("basis_set_exchange")

using HDF5

#== initialization ==#
upper_to_lower = Dict(
                      "A" => "a",
                      "B" => "b",
                      "C" => "c",
                      "D" => "d",
                      "E" => "e",
                      "F" => "f",
                      "G" => "g",
                      "H" => "h",
                      "I" => "i",
                      "J" => "j",
                      "K" => "k",
                      "L" => "l",
                      "M" => "m",
                      "N" => "n",
                      "O" => "o",
                      "P" => "p",
                      "Q" => "q",
                      "R" => "r",
                      "S" => "s",
                      "T" => "t",
                      "U" => "u",
                      "V" => "v",
                      "W" => "w",
                      "X" => "x",
                      "Y" => "y",
                      "Z" => "z"
                     )
#==============================#
#== get atom/basis set pairs ==#
#==============================#

basis_list = Vector{String}([])
status_bsed = h5open("../records/bsed.h5","r") do bsed
  #== loop over atoms ==#
  for symbol in names(bsed)
    #== loop over basis sets in atom ==#
    for basis in names(bsed["$symbol"])
      #== add basis set to pair list if unique ==#
      if !in(basis, basis_list)
        push!(basis_list, basis)
      end
    end
  end
end

for basis in basis_list

  bs_psi4 = bse.get_basis(basis,fmt="psi4", header=false)
  bs_psi4 = replace(bs_psi4, "D+" => "E+")
  bs_psi4 = replace(bs_psi4, "D-" => "E-")

  bs_name = basis
  bs_name = replace(bs_name, "*" => "s")
  bs_name = replace(bs_name, "+" => "p")
  bs_name = replace(bs_name, "(" => "_")
  bs_name = replace(bs_name, ")" => "_")
  bs_name = replace(bs_name, "," => "_")
  for letter in upper_to_lower
    bs_name = replace(bs_name, letter[1] => letter[2])
  end

  println(bs_name)

  bs_file = open("$bs_name.gbs", "w") do file
    write(file, bs_psi4)
  end
end


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

## dummy stub for https://github.com/JuliaBinaryWrappers/GMP_jll.jl
module GMP_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libgmp, libgmpxx

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libgmp_handle = C_NULL
libgmp_path = ""
libgmpxx_handle = C_NULL
libgmpxx_path = ""

if Sys.iswindows()
    const libgmp = "libgmp-10.dll"
    const libgmpxx = "libgmpxx-4.dll"
elseif Sys.isapple()
    const libgmp = "@rpath/libgmp.10.dylib"
    const libgmpxx = "@rpath/libgmpxx.4.dylib"
else
    const libgmp = "libgmp.so.10"
    const libgmpxx = "libgmpxx.so.4"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libgmp_handle = dlopen(libgmp)
    global libgmp_path = dlpath(libgmp_handle)
    global libgmpxx_handle = dlopen(libgmpxx)
    global libgmpxx_path = dlpath(libgmpxx_handle)
end

is_available() = true

end  # module GMP_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
using DirectedRoundings
using Test

# arithmetic function examples
#   the directed rounding functions work with 1,2,3,4 params

function sqrt5x(a)
    sqrt(5*a)
end

function reldiff(a, b)
    return (a - b) / a
end

function affine_map(x, shift, scale)
    return (x + shift) * scale
end

function right_triangles_ratio(a_adjacent, a_opposite,
                               b_adjacent, b_opposite)
    a_hypot = hypot(a_adjacent, a_opposite)
    b_hypot = hypot(b_adjacent, b_opposite)
    return a_hypot / b_hypot
end

const FP = Float32

a = FP(3); b = FP(4096); c = inv(FP(17)); d = eps(a)/2;

round_nearest(inv,a)
round_up(inv, a)
round_down(inv, a)
round_tozero(inv, a)
round_fromzero(inv, a)

round_up(+, a, b)
round_down(hypot, a, b)

round_up(affine_map, a, b, c)
round_down(affine_map, a, b, c)

round_nearest(right_triangles_ratio, a, b, c, d)
round_up(right_triangles_ratio, a, b, c, d)
round_down(right_triangles_ratio, a, b, c, d)


text:
julia> using IntervalArithmetic

julia> n = 2.0
2.0

julia> @interval 2n/3 + 1
[2.33333, 2.33334]

julia> showall(ans)
Interval(2.333333333333333, 2.3333333333333335)

julia> a = @interval(0.1, 0.3)
[0.0999999, 0.300001]

julia> b = @interval(0.3, 0.6)
[0.299999, 0.600001]

julia> a + b
[0.399999, 0.900001]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function alphavalue(c)
	int(c) - int('A') +1
end

function totalalphavalue(s)
	mapreduce(alphavalue, +, collect(s))
end

function main()
	filepath = pwd() * "/p022_names.txt"
	text = open(readall, filepath)
	text = replace(text, ['"','\\'], "")
	names = split(text,',')

	sort!(names)

	namevalues = map(totalalphavalue, names)
	sum(namevalues .* (1:length(namevalues)))

end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL
using QHull, Plasm

V,CV = Lar.rod()([4,1])
ch = QHull.chull(convert(Lar.Points,V'))
FV = ch.simplices
EV = Lar.simplexFacets(FV)
GL.VIEW([ GL.GLGrid(V,EV,GL.COLORS[1],1) ]);

cyl = Lar.Struct([ Lar.t(0,0,-1.5), (V,FV,EV) ])
cyl = Lar.struct2lar(cyl)

V,FV = ([1.0 1.0 -0.5 -0.5 -0.5 -0.5; 0.0 0.0 0.86603 0.86603 -0.86603 -0.86603; 0.0 3.0 0.0 3.0 0.0 3.0],
Array{Int64,1}[[3, 4, 1, 2], [3, 4, 5, 6], [5, 6, 1, 2], [3, 1, 5], [2, 4, 6]])
EV = [[1,3],[3,5],[5,1], [2,4],[4,6],[6,2], [1,2],[3,4],[5,6]]
cyl = (V,FV,EV)

tris = Lar.Struct([ cyl,
                    Lar.Struct([Lar.r(pi/2,0,0), cyl ]) ,
                    Lar.Struct([Lar.r(0,pi/2,0), cyl ]) ,
                    Lar.Struct([Lar.r(0,pi,0), cyl ]) ,
                    ])
V,FV,EV = Lar.struct2lar(tris)
Plasm.view(V,EV) #GL.VIEW([ GL.GLGrid(V,FV) ]);

#VV = [[k] for k=1:size(V,2)]
#GL.VIEW( GL.numbering(.2)((V,[VV,EV,FV]), GL.COLORS[4]) );

cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
#cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
cop_FE = Lar.coboundary_1(FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EV, cop_FE)

W = convert(Lar.Points, V')
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)

GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.,1.,1.,99,1));
meshes = GL.GLExplode(V,CVs[1:end],4,4,4,99,0.5);
GL.VIEW( push!( meshes, GL.GLFrame) );


text:
# expand_deps.jl - collect & expand all dependencies of a variable in  AbstractExGraph

## collect_deps

function collect_deps!(g::AbstractExGraph, nd::ExNode, depth::Int, result::Set{Symbol})
    if depth > 0
        for dep in dependencies(nd)
            if haskey(g, dep)
                collect_deps!(g, g[dep], depth - 1, result)
            end
        end
    end
    push!(result, varname(nd))
end

function collect_deps!(g::AbstractExGraph, ex::Expr, depth::Int, result::Set{Symbol})
    vnames = get_var_names(ex, rec=true)
    for vname in vnames
        collect_deps!(g, vname, depth, result)
    end
end

function collect_deps!(g::AbstractExGraph, x::Symbol, depth::Int, result::Set{Symbol})
    if haskey(g, x)
        collect_deps!(g, g[x], depth, result)
    end
end

function collect_deps!(g::AbstractExGraph, x, depth::Int, result::Set{Symbol})
    # do nothing
end

function collect_deps(g::AbstractExGraph, nd::ExNode, depth::Int=typemax(Int))
    result = Set{Symbol}()
    collect_deps!(g, nd, depth, result)
    return result
end

function collect_deps(g::AbstractExGraph, ex::Expr, depth::Int=typemax(Int))
    result = Set{Symbol}()
    collect_deps!(g, ex, depth, result)
    return result
end

function collect_deps(g::AbstractExGraph, x::Symbol, depth::Int=typemax(Int))
    result = Set{Symbol}()
    collect_deps!(g, x, depth, result)
    return result
end

function collect_deps(g::AbstractExGraph, xs::Vector{Symbol}, depth::Int=typemax(Int))
    result = Set{Symbol}()
    for x in xs
        collect_deps!(g, x, depth, result)
    end
    return result
end

collect_deps(g::AbstractExGraph, x, depth::Int=typemax(Int)) = Set{Symbol}()

## expand_deps

function expand_deps!(g::AbstractExGraph, nd::ExNode{:input}, depth::Int, result::Vector{Expr})
    # do nothing
end

function expand_deps!(g::AbstractExGraph, nd::ExNode{:constat}, depth::Int, result::Vector{Expr})
    push!(result, to_expr(nd))
end

function expand_deps!(g::AbstractExGraph, nd::ExNode{:(=)}, depth::Int, result::Vector{Expr})
    if depth > 0
        expand_deps!(g, [g[var] for var in dependencies(nd)], depth - 1, result)
        push!(result, to_expr(nd))
    end
end

function expand_deps!(g::AbstractExGraph, nd::ExNode{:call}, depth::Int, result::Vector{Expr})
    if depth > 0
        for dep in dependencies(nd)
            expand_deps!(g, g[end], depth - 1, result)
        end
        push!(result, to_expr(nd))
    end
end

function expand_deps(g::AbstractExGraph, nd::ExNode, depth::Int=typemax(Int))
    deps = collect_deps(g, nd, depth)
    ex = Expr(:block)
    for nd in g.tape
        if !isa(nd, ExNode{:input}) && in(varname(nd), deps)
            push!(ex.args, to_expr(nd))
        end
    end
    return ex
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#frag-mult-test.jl
#test fragmented multiplication

function frag_mult(a::Array{UInt64,1}, b::Array{UInt64,1})
  #note that frag_mult fails for absurdly high length integer arrays.
  length(a) != length(b) && throw(ArgumentError("mismatched arrays"))

  #take these two UInt64 arrays and reinterpret them as UInt32 arrays
  a_32 = reinterpret(UInt32, a)
  b_32 = reinterpret(UInt32, b)

  scratchpad = zeros(UInt32, length(a) * 4)
  carries    = zeros(UInt32, length(a) * 4)

  for idx = 1:length(a_32)
    for jdx = 1:length(b_32)
      #multiply the two values.
      multres::UInt64 = a_32[idx] * b_32[jdx]
      #bin these into a high word and a lwo word
      multresbins = reinterpret(UInt32, [multres])
      #add the low word to the scratchpad
      scratchpad[idx + jdx - 1] += multresbins[1]
      #save the carry, if need be.
      (scratchpad[idx + jdx - 1] < multresbins[1]) && (carries[idx + jdx] += 1)
      #add the high word to the scratchpad
      scratchpad[idx + jdx] += multresbins[2]
      #save the carry, if need be.
      (scratchpad[idx + jdx] < multresbins[2]) && (carries[idx + jdx + 1] += 1)
    end
  end
  #go through and resolve the carries.
  for idx = 1:length(carries)
    scratchpad[idx] += carries[idx]
    (scratchpad[idx] < carries[idx]) && (carries[idx + 1] += 1)
  end
  reinterpret(UInt64, scratchpad)
end

function calculate(a::Array{UInt64,1})
  sum = big(0)
  for idx = 1:length(a)
    sum += big(a[idx]) << (64 * (idx - 1))
  end
  sum
end

<<<<<<< HEAD:oddsandends/frag-mult-test.jl
function test_mults(cells::Integer)
=======
function test_mults(cells::Int)
>>>>>>> directory-reorg:src/oddsandends/frag-mult-test.jl
  x = rand(UInt64, cells)
  y = rand(UInt64, cells)

  x_big = calculate(x)
  y_big = calculate(y)

  #calculate using fragment multiplication
  res_frag = calculate(frag_mult(x,y))
  res_big = x_big * y_big

  println(res_frag)
  println(res_big)
end

#test that fragment multiplication in general works.
test_mults(1)
test_mults(2)
test_mults(4)
test_mults(8)

#next up:  truncated fragment multiplication.
function trunc_frag_mult(a::Array{UInt64,1}, b::Array{UInt64,1})
  #note that frag_mult fails for absurdly high length integer arrays.
  length(a) != length(b) && throw(ArgumentError("mismatched arrays"))

  #take these two UInt64 arrays and reinterpret them as UInt32 arrays
  a_32 = reinterpret(UInt32, a)
  b_32 = reinterpret(UInt32, b)
  l = length(a_32)

  scratchpad = zeros(UInt32, l + 1)
  carries    = zeros(UInt32, l)

  #first indexsum is length(a_32)
  indexsum = l
  for (aidx = 1:(indexsum - 1))
    temp_res::UInt64 = a_32[aidx] * b_32[indexsum - aidx]
    temp_res_high::UInt32 = (temp_res >> 32)
    scratchpad[1] += temp_res_high
    (scratchpad[1] < temp_res_high) && (carries[1] += 1)
  end
  #now proceed with the rest of the additions.
  for aidx = 1:l
    for bidx = (l + 1 - aidx):l
      temp_res = a_32[aidx] * b_32[bidx]
      temp_res_bins = reinterpret(UInt32, [temp_res])
      temp_res_low::UInt32 = temp_res
      temp_res_high = (temp_res >> 32)

      scratchindex = aidx + bidx - l

      scratchpad[scratchindex] += temp_res_low
      (temp_res_low > scratchpad[scratchindex]) && (carries[scratchindex] += 1)

      scratchpad[scratchindex + 1] += temp_res_high
      (temp_res_high > scratchpad[scratchindex + 1]) && (carries[scratchindex + 1] += 1)
    end
  end

  #go through and resolve the carries.
  for idx = 1:length(carries) - 1
    scratchpad[idx + 1] += carries[idx]
    (scratchpad[idx + 1] < carries[idx]) && (carries[idx + 1] += 1)
  end
  reinterpret(UInt64, scratchpad[2:end])
end

<<<<<<< HEAD:oddsandends/frag-mult-test.jl
function test_truncs(cells::Integer)
=======
function test_truncs(cells::Int)
>>>>>>> directory-reorg:src/oddsandends/frag-mult-test.jl
  x = rand(UInt64, cells)
  y = rand(UInt64, cells)

  z = frag_mult(x, y)[cells + 1:2 *cells]
  w = trunc_frag_mult(x, y)

  println(z)
  println(w)
end

test_truncs(1)
test_truncs(2)
test_truncs(4)
test_truncs(8)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Five, Hydrothermal Venture."

include("Utils.jl")

using LinearAlgebra
using .Utils

const Line = Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}

"Calculate the endpoints of the `line`."
function calculate_endpoints(line :: Line)
    x₀, x₁ = sort([first.(line)...]) .+ 1
    y₀, y₁ = sort([last.(line)...]) .+ 1

    return ((x₀, y₀), (x₁, y₁))
end

"Calculate the slope of the `line`."
function calculate_slope(line :: Line)
    x₀, x₁ = first.(line)
    y₀, y₁ = last.(line)

    return (y₁ - y₀)/(x₁ - x₀)
end

"Superimpose the `lines` over the `plane` and update it (mutating)."
function update_lines!(lines :: Array{Line}, plane :: Array{Int64, 2})
    for l in lines
        ((x₀, y₀), (x₁, y₁)) = calculate_endpoints(l)
        slope = calculate_slope(l)

        if slope in [-Inf, 0, Inf]
            plane[y₀:y₁, x₀:x₁] += ones(Int, (y₁ - y₀ + 1, x₁ - x₀ + 1))
        else
            # Number of rotations.
            rot = Int64(div(-slope + 1, 2))
            plane[y₀:y₁, x₀:x₁] += rotr90(Diagonal(ones(Int, x₁ - x₀ + 1)), rot)
        end
    end
end

"""Calculate the number of intersections of at least two `lines` in a ℤ
cartesian plane.  Return both the axial and the total count."""
function solve(lines :: Array{Line})
    xₘₐₓ = max([max(first.(l)...) for l in lines]...)
    yₘₐₓ = max([max(last.(l)...) for l in lines]...)

    plane = zeros(Int, (yₘₐₓ + 1, xₘₐₓ + 1))

    #= In the context of this day, I used the term “axial” as referring strictly
     = to vertical and horizontal lines only.
     =#
    axials = filter(P -> P[1][1] == P[2][1] || P[1][2] == P[2][2], lines)
    update_lines!(axials, plane)
    axial_count = count(>=(2), plane)

    diagonals = filter(P -> P[1][1] != P[2][1] && P[1][2] != P[2][2], lines)
    update_lines!(diagonals, plane)
    total_count = count(>=(2), plane)

    return (axial_count, total_count)
end

const POINTS_REX = r"(\d+),(\d+) -> (\d+),(\d+)"

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = readlines(ARGS[1])
    matched_data = match.(POINTS_REX, input_data)
    parsed_data = [parse.(Int, line.captures) for line in matched_data]
    lines = [((x₀, y₀), (x₁, y₁)) for (x₀, y₀, x₁, y₁) in parsed_data]

    S₁, S₂ = solve(lines)

    @assert S₁ == 6710
    @assert S₂ == 20121

    println(S₁)
    println(S₂)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# Demonstrate/test functions in cornerfeatures

using ImageProjectiveGeometry
using FileIO
using PyPlot

rad = 6

#img = PyPlot.imread("../sampleimages/lena.tif");
img = load("../sampleimages/lena.tif");

# Convert image from Gray Fixed Point Numbers to a simple array of floats
img = float(img);
figure(10), PyPlot.imshow(img);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Input image")

println("Harris corner detector")
(cimg, r, c) = ImageProjectiveGeometry.harris(img, 1,k= 0.04, radius=rad, N=100, subpixel=false, img=img, fig=12)
figure(12); PyPlot.title("Harris corners")
cimg[cimg .< 0] .= 0;   # eliminate -ve values
figure(2); PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Harris corner strength image")
ImageProjectiveGeometry.keypause()

println("Noble corner detector")
(cimg, r, c) = ImageProjectiveGeometry.noble(img, 1, radius=rad, N=100, subpixel=false, img=img, fig=12);
figure(12); PyPlot.title("Noble corners")
figure(2); PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Noble corner strength image")
ImageProjectiveGeometry.keypause()

println("Shi-Tomasi corner detector")
(cimg, r, c) = ImageProjectiveGeometry.shi_tomasi(img, 1, radius=rad, N=100, subpixel=false, img=img, fig=12);
figure(12); PyPlot.title("Shi_Tomasi corners")
figure(2); PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Shi-Tomasi corner strength image")
ImageProjectiveGeometry.keypause()

println("-ve Hessian features")
hdet = ImageProjectiveGeometry.hessianfeatures(img, 1)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(-hdet, N=100, img=img, fig=12)
figure(12); PyPlot.title("-ve Hessian features")
figure(2); PyPlot.imshow(-hdet);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("-ve Hessian image")
ImageProjectiveGeometry.keypause()

println("+ve Hessian features")
hdet = ImageProjectiveGeometry.hessianfeatures(img, 1)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(+hdet, N=100, img=img, fig=12)
figure(12); PyPlot.title("+ve Hessian features")
figure(2); PyPlot.imshow(hdet);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("+ve Hessian image")
ImageProjectiveGeometry.keypause()

println("Coherence")
cimg = ImageProjectiveGeometry.coherence(img, 2)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(cimg, N=10, img=img, fig=12)
figure(2); PyPlot.imshow(cimg);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Coherence image")
ImageProjectiveGeometry.keypause()

println("Fast radial")
(S, So) = ImageProjectiveGeometry.fastradial(img, [1, 3, 5, 7])
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(-S, N=10, img=img, fig=12)
figure(12); PyPlot.title("Fast radial features")
figure(2); PyPlot.imshow(-S);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Fast radial image: orientation and gradient magnitude image")

figure(3); PyPlot.imshow(-So);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Fast radial image: orientation only image")

nothing


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

mutable struct VirtualMachineScaleSetExtensionProperties <: SwaggerModel
    forceUpdateTag::Any # spec type: Union{ Nothing, String } # spec name: forceUpdateTag
    publisher::Any # spec type: Union{ Nothing, String } # spec name: publisher
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    typeHandlerVersion::Any # spec type: Union{ Nothing, String } # spec name: typeHandlerVersion
    autoUpgradeMinorVersion::Any # spec type: Union{ Nothing, Bool } # spec name: autoUpgradeMinorVersion
    settings::Any # spec type: Union{ Nothing, Any } # spec name: settings
    protectedSettings::Any # spec type: Union{ Nothing, Any } # spec name: protectedSettings
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    provisionAfterExtensions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: provisionAfterExtensions

    function VirtualMachineScaleSetExtensionProperties(;forceUpdateTag=nothing, publisher=nothing, type=nothing, typeHandlerVersion=nothing, autoUpgradeMinorVersion=nothing, settings=nothing, protectedSettings=nothing, provisioningState=nothing, provisionAfterExtensions=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("forceUpdateTag"), forceUpdateTag)
        setfield!(o, Symbol("forceUpdateTag"), forceUpdateTag)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("publisher"), publisher)
        setfield!(o, Symbol("publisher"), publisher)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("typeHandlerVersion"), typeHandlerVersion)
        setfield!(o, Symbol("typeHandlerVersion"), typeHandlerVersion)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("autoUpgradeMinorVersion"), autoUpgradeMinorVersion)
        setfield!(o, Symbol("autoUpgradeMinorVersion"), autoUpgradeMinorVersion)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("settings"), settings)
        setfield!(o, Symbol("settings"), settings)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("protectedSettings"), protectedSettings)
        setfield!(o, Symbol("protectedSettings"), protectedSettings)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualMachineScaleSetExtensionProperties, Symbol("provisionAfterExtensions"), provisionAfterExtensions)
        setfield!(o, Symbol("provisionAfterExtensions"), provisionAfterExtensions)
        o
    end
end # type VirtualMachineScaleSetExtensionProperties

const _property_map_VirtualMachineScaleSetExtensionProperties = Dict{Symbol,Symbol}(Symbol("forceUpdateTag")=>Symbol("forceUpdateTag"), Symbol("publisher")=>Symbol("publisher"), Symbol("type")=>Symbol("type"), Symbol("typeHandlerVersion")=>Symbol("typeHandlerVersion"), Symbol("autoUpgradeMinorVersion")=>Symbol("autoUpgradeMinorVersion"), Symbol("settings")=>Symbol("settings"), Symbol("protectedSettings")=>Symbol("protectedSettings"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("provisionAfterExtensions")=>Symbol("provisionAfterExtensions"))
const _property_types_VirtualMachineScaleSetExtensionProperties = Dict{Symbol,String}(Symbol("forceUpdateTag")=>"String", Symbol("publisher")=>"String", Symbol("type")=>"String", Symbol("typeHandlerVersion")=>"String", Symbol("autoUpgradeMinorVersion")=>"Bool", Symbol("settings")=>"Any", Symbol("protectedSettings")=>"Any", Symbol("provisioningState")=>"String", Symbol("provisionAfterExtensions")=>"Vector{String}")
Base.propertynames(::Type{ VirtualMachineScaleSetExtensionProperties }) = collect(keys(_property_map_VirtualMachineScaleSetExtensionProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetExtensionProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetExtensionProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetExtensionProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetExtensionProperties[property_name]

function check_required(o::VirtualMachineScaleSetExtensionProperties)
    true
end

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


text:
# rotate2.jl
# time 2D image rotation

using SPECTrecon: plan_rotate
using SPECTrecon: imrotate!, imrotate_adj!
using BenchmarkTools: @btime

function imrotate1_time()
    T = Float32
    N = 100
    img = zeros(T, N, N)
    img[30:50, 20:60] .= 1
    output = similar(img)
    θ = 3*π/16
    plan = plan_rotate(N; nthread = 1, T = T, method = :one)[1]
    println("imrotate1-forw")
    @btime imrotate!($output, $img, $θ, $plan)
    # 533.316 μs (0 allocations: 0 bytes)
    println("imrotate1-back")
    @btime imrotate_adj!($output, $img, $θ, $plan)
    # 517.367 μs (0 allocations: 0 bytes)
    nothing
end

function imrotate2_time()
    T = Float32
    N = 100
    img = zeros(T, N, N)
    img[30:50, 20:60] .= 1
    output = similar(img)
    θ = 3*π/16
    plan = plan_rotate(N; nthread = 1, T = T, method = :two)[1]
    println("imrotate2-forw")
    @btime imrotate!($output, $img, $θ, $plan)
    # 162.171 μs (0 allocations: 0 bytes)
    println("imrotate2-back")
    @btime imrotate_adj!($output, $img, $θ, $plan)
    # 168.518 μs (0 allocations: 0 bytes)
    nothing
end

# run all functions, time may vary on different machines, but should be all zero allocation.
imrotate1_time()
imrotate2_time()


text:
module Simplifications

using Sylvia

export simplify

const base_ctx = Sylvia.Context(Sylvia.@__context__)
const ctx = Sylvia.Context(base_ctx)

function simplify end
Sylvia.@register ctx simplify 1

@sym [Number] a b c d

@scope base_ctx let
    @sym [Wild{Number}] x y z
    @! set simplify(x) --> x
end

@scope ctx let
    @sym [Wild{Number}] x y z
    _0 = @! 0::Number
    _1 = @! 1::Number

    @! set simplify(x*y) --> simplify(x)*simplify(y) |> simplify
    @! set simplify(x*y*z) --> simplify(x*y)*simplify(z) |> simplify
    @! set simplify(x + y) --> simplify(x) + simplify(y) |> simplify
    @! set simplify(x + y + z) --> simplify(x + y) + simplify(z) |> simplify
    @! set simplify(x*y + y) --> simplify(x + 1)*simplify(y) |> simplify
    @! set simplify(x*z + y*z) --> simplify(x + y)*simplify(z) |> simplify

    @! set simplify(x + x) --> 2simplify(x) |> simplify
    @! set simplify(x*x) --> simplify(x)^2 |> simplify
    @! set simplify(x*x*x) --> simplify(x*x)*simplify(x) |> simplify

    @! set simplify(x/x) --> _1
    @! set simplify(x/_1) --> simplify(x)

    @! set simplify(-x) --> -simplify(x)
    @! set simplify(-(-x)) --> simplify(x)

    @! set simplify(_0 + x) --> simplify(x)
    @! set simplify(x + _0) --> simplify(x)

    @! set simplify(_0*x) --> _0
    @! set simplify(x*_0) --> _0
    @! set simplify(_1*x) --> simplify(x)
    @! set simplify(x*_1) --> simplify(x)

    @! set simplify(x^_1) --> simplify(x)
end

end


text:
using DeltaRobot
using Test

delta_robot_params = DeltaRobotParams(2.5, 3.0, 4.0, 0.312, 0.410, 9.8, 0.114, 0.051, 0.0106, 0.0106, 0.0106); # Actual robot parameters with random d-params

@testset "Inverse Kinematics" begin

end

@testset "Forward Kinematics" begin
    p = DeltaRobotArmParams(0.3, 0.2, 0.1, 0.1, 0.3, 9.8, 0.075, 0.024)

end

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
module Lanczos

function lanczos(A,v1,steps)

    size(A)[1] == size(A)[2] || error("matrix must be square")

    ∘ = dot

    v = zeros(size(v1))
    r = v1
    b = 1

    T = zeros(Float64,steps+1,steps+1)      # Lanczos matrix
    K = zeros(Float64,size(A)[1],steps+1)   # vectors in Krylov space

    # First round
    vold = v
    v = r/b
    K[:,1] = v
    Av = A*v
    a = v∘Av
    T[1,1] = a
    r = Av - v*a - vold*b

    # Loop
    for i in 1:steps
        b = norm(r)
        T[i,i+1] = T[i+1,i] = b
        vold = v
        v = r/b                     # next v
        K[:,i+1] = v                # store vector
        Av = A*v
        a = v∘Av
        T[i+1,i+1] = a
        r = Av - v*a - vold*b       # orthogonalize
    end

    return T, K

end

lanczos(A,v1) = lanczos(A,v1,size(A)[1])

end