text:
import Base: copy, reset

abstract type Model end

"""
$(TYPEDEF)

Dynamics model

Holds all information required to uniquely describe a dynamic system, including
a general nonlinear dynamics function of the form `ẋ = f(x,u)`, where x ∈ ℜⁿ are
the states and u ∈ ℜᵐ are the controls.

Dynamics function `Model.f` should be in the following forms:
    'f!(ẋ,x,u)' and modify ẋ in place
"""
struct AnalyticalModel <:Model
    f::Function # continuous dynamics (ie, differential equation)
    n::Int # number of states
    m::Int # number of controls
    state_error::Function #how to find the error for the rollout
    exp::Function #expansion function for finding directional derivatives
    evals::Vector{Int}
    info::Dict{Symbol,Any}

    # Construct a model from an explicit differential equation
    function AnalyticalModel(f::Function, n::Int64, m::Int64,
         state_error::Function, exp::Function, d::Dict{Symbol,Any}=Dict{Symbol,Any}())
        d[:evals] = 0
        evals = [0,]
        # Make dynamics inplace
        if is_inplace_dynamics(f,n,m)
            f! = f
        else
            f! = wrap_inplace(f)
        end
        new(f!,n,m,state_error,exp,evals,d)
    end
end

""" $(SIGNATURES) Create a Model given an inplace analytical function for the continuous dynamics with n states and m controls"""
Model(f::Function, n::Int64, m::Int64, state_error::Function, exp::Function, d::Dict{Symbol,Any}=Dict{Symbol,Any}()) = AnalyticalModel(f,n,m,state_error,exp,d)

evals(model::Model) = model.evals[1]
reset(model::Model) = begin model.evals[1] = 0; return nothing end

function dynamics(model::Model,xdot,x,u)
    model.f(xdot,x,u)
    model.evals[1] += 1
end

"""
$(TYPEDEF)
RigidBodyDynamics model. Wrapper for a RigidBodyDynamics Mechanism
"""
struct RBDModel <: Model
    f::Function # continuous dynamics (ie, differential equation)
    n::Int # number of states
    m::Int # number of controls
    mech::Mechanism  # RigidBodyDynamics Mechanism
    evals::Vector{Int}
    info::Dict{Symbol,Any}
end

"""
$(SIGNATURES)
 Construct model from a `Mechanism` type from `RigidBodyDynamics`
 """
function Model(mech::Mechanism)
    m = length(joints(mech))  # subtract off joint to world
    Model(mech,ones(m))
end

"""
$(SIGNATURES)
Model(mech::Mechanism, torques::Array{Bool, 1}) Constructor for an underactuated mechanism, where torques is a binary array
that specifies whether a joint is actuated.
"""
function Model(mech::Mechanism, torques::Array)

    # construct a model using robot dynamics equation assembed from URDF file
    n = num_positions(mech) + num_velocities(mech) + num_additional_states(mech)
    num_joints = length(joints(mech))  # subtract off joint to world

    if length(torques) != num_joints
        error("Torque underactuation specified does not match mechanism dimensions")
    end

    m = convert(Int,sum(torques)) # number of actuated (ie, controllable) joints
    torque_matrix = 1.0*Matrix(I,num_joints,num_joints)[:,torques.== 1] # matrix to convert from control inputs to mechanism joints

    statecache = StateCache(mech)
    dynamicsresultscache = DynamicsResultCache(mech)

    function f(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T
        state = statecache[T]
        dyn = dynamicsresultscache[T]
        dynamics!(view(ẋ,1:n), dyn, state, x, torque_matrix*u)
        return nothing
    end
    d = Dict{Symbol,Any}()

    evals = [0,]
    RBDModel(f, n, m, mech, evals, d)
end

"""$(SIGNATURES) Construct a fully actuated model from a string to a urdf file"""
function Model(urdf::String)
    # construct model using string to urdf file
    mech = parse_urdf(urdf)
    Model(mech)
end

"""$(SIGNATURES) Construct a partially actuated model from a string to a urdf file, where torques is a binary array that specifies whether a joint is actuated."""
function Model(urdf::String,torques::Array{Float64,1})
    # underactuated system (potentially)
    mech = parse_urdf(urdf)
    Model(mech,torques)
end

"""
$(SIGNATURES)
Determine if the dynamics in model are in place. i.e. the function call is of
the form `f!(xdot,x,u)`, where `xdot` is modified in place. Returns a boolean.
"""
function is_inplace_dynamics(model::Model)::Bool
    x = rand(model.n)
    u = rand(model.m)
    xdot = rand(model.n)
    try
        model.f(xdot,x,u)
    catch x
        if x isa MethodError
            return false
        end
    end
    return true
end

function is_inplace_dynamics(f::Function,n::Int64,m::Int64)::Bool
    x = rand(n)
    u = rand(m)
    xdot = rand(n)
    try
        f(xdot,x,u)
    catch x
        if x isa MethodError
            return false
        end
    end
    return true
end

"""
$(SIGNATURES)
Makes the dynamics function `f(x,u)` appear to operate as an inplace operation of the
form `f!(xdot,x,u)`.
"""
function wrap_inplace(f::Function)
    f!(xdot,x,u) = copyto!(xdot, f(x,u))
    f!(xdot,x) = copyto!(xdot, f(x))
end


text:
using OceanTurb, Printf

function diffusivity!(K, KU, model)
    for i in eachindex(K)
        @inbounds K[i] = KU(model, i)
    end
    return nothing
end

@use_pyplot_utils

usecmbright()

c = Constants(f=1e-4, β=0.0)

modelsetup = (N=128, L=128, stepper=:BackwardEuler, constants=c)

Fb = 0.0 #-1e-7
Fu = -1e-4
Fe = 0.0 #-1e-9
N² = 5e-6
Δt = 1
times = (0, 2, 8, 32) .* hour

K₀ = 1e-5

kpp = ModularKPP.Model(; modelsetup...,
    diffusivity = ModularKPP.LMDDiffusivity(KU₀=K₀, KT₀=K₀, KS₀=K₀)
    )

tke = KPP_TKE.Model(; modelsetup...,
    diffusivity = ModularKPP.LMDDiffusivity(KU₀=K₀, KT₀=K₀, KS₀=K₀),
            tke = KPP_TKE.TKEParameters(Cτ=Inf, CDe=0.5, KU₀=K₀, KT₀=K₀, KS₀=K₀, Ke₀=1e-5)
    )

tkeK = FaceField(tke.grid)
kppK = FaceField(kpp.grid)

# Initial condition and fluxes
Fθ = Fb / (c.α * c.g)
dTdz = N² / (c.α * c.g)
T₀(z) = 20 + dTdz*z
S₀(z) = 20 + dTdz*z

models = (kpp, tke)

for model in models
    model.solution.T = T₀
    model.solution.S = S₀

    model.bcs.T.top = FluxBoundaryCondition(Fθ)
    model.bcs.T.bottom = GradientBoundaryCondition(dTdz)

    model.bcs.U.top = FluxBoundaryCondition(Fu)
end

tke.bcs.e.top = FluxBoundaryCondition(Fe)

fig, axs = subplots(ncols=4, figsize=(8, 6))

sca(axs[1])
removespines("top", "right")
xlabel(L"T")
ylabel(L"z \, \mathrm{(m)}")

sca(axs[2])
removespines("top", "left", "right")
axs[2].tick_params(left=false, labelleft=false)
xlabel(L"U")

#=
sca(axs[3])
removespines("top", "left", "right")
axs[3].tick_params(left=false, labelleft=false)
xlabel(L"V")
=#

sca(axs[3])
removespines("top", "left", "right")
axs[3].tick_params(left=false, labelleft=false)
xlabel(L"e")

sca(axs[4])
removespines("top", "left")
axs[4].tick_params(left=false, labelleft=false, right=true, labelright=true)
xlabel(L"K")
ylabel(L"z \, \mathrm{(m)}")

for i = 1:length(times)
    for model in models
        run_until!(model, Δt, times[i])
    end

    diffusivity!(kppK, ModularKPP.KU, kpp)
    diffusivity!(tkeK, KPP_TKE.K_mixing_time, tke)

    @printf("""
        t : %.1f hours

          surface temperature
          =============

               kpp T(z=0) : %.6f
               tke T(z=0) : %.6f
    \n""", time(kpp)/hour,
    kpp.solution.T[end-1],
    tke.solution.T[end-1],
    )

    if i == 1
        vlabel = (kpp="CVMix KPP", tke="TKE")
    else
        vlabel = (kpp="", tke="")
    end

    #=
    if i == 1
        tlabel = text(kpp.solution.T[end], 0.5,
            @sprintf("\$ t = %.0f \$ hours", time(kpp)/hour),
            verticalalignment="bottom", horizontalalignment="center", color=defaultcolors[i])
    else
        tlabel = text(maximum(kpp.solution.T.data)-0.003, -kpp.state.h,
            @sprintf("\$ t = %.0f \$ hours", time(kpp)/hour),
            verticalalignment="bottom", horizontalalignment="left", color=defaultcolors[i])
    end
    =#

    sca(axs[1])
    plot(kpp.solution.T, "-",  color=defaultcolors[i], label=vlabel.kpp, alpha=0.8, markersize=1.5)
    plot(tke.solution.T, "--", color=defaultcolors[i], label=vlabel.tke, alpha=0.8, markersize=1.5)

    sca(axs[2])
    plot(kpp.solution.U, "-",  color=defaultcolors[i], label=vlabel.kpp, alpha=0.8, markersize=1.5)
    plot(tke.solution.U, "--", color=defaultcolors[i], label=vlabel.tke, alpha=0.8, markersize=1.5)

    #=
    sca(axs[3])
    plot(kpp.solution.V, "-",  color=defaultcolors[i], label=vlabel.kpp, alpha=0.8, markersize=1.5)
    plot(tke.solution.V, "--", color=defaultcolors[i], label=vlabel.tke, alpha=0.8, markersize=1.5)
    =#

    sca(axs[3])
    plot(tke.solution.e, "-", color=defaultcolors[i], label=vlabel.tke, alpha=0.8, markersize=1.5)

    sca(axs[4])
    plot(kppK, "-", color=defaultcolors[i], label=vlabel.kpp, alpha=0.8, markersize=1.5)
    plot(tkeK, "--", color=defaultcolors[i], label=vlabel.tke, alpha=0.8, markersize=1.5)
end

#=
for ax in axs
    sca(ax)
    legend(fontsize=10)
end
=#

gcf()

ext: jl


text:
module test_24_Reuse

using ModiaResult
using Unitful
using OrderedCollections
ModiaResult.@usingModiaPlot

t = range(0.0, stop=10.0, length=100)

result1 = OrderedDict{String,Any}()

result1["time"] = t*u"s"
result1["phi"]  = sin.(t)u"rad"
result1["w"]    = cos.(t)u"rad/s"

result2 = OrderedDict{String,Any}()
result2["time"] = t*u"s"
result2["phi"]  = 1.2*sin.(t)u"rad"
result2["w"]    = 0.8*cos.(t)u"rad/s"

println("\n... test_24_Reuse:")
ModiaResult.printResultInfo(result1)
ModiaResult.printResultInfo(result2)

plot(result1, ("phi", "w"), prefix="Sim 1:", heading="Test reuse")
plot(result2, ("phi", "w"), prefix="Sim 2:", reuse=true)

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
############################################################################################
# Include files
include("generated.jl")
include("helper.jl")
include("data.jl")
include("utility.jl")
include("chains.jl")
include("jitter.jl")
include("printing.jl")

############################################################################################
# Export


text:
#! /usr/bin/julia

# Rosetta Code, Handle a signal

ccall(:jl_exit_on_sigint, Void, (Cint,), 0)
tic()
ticks = 0
try
    while true
        sleep(0.5)
        ticks += 1
        println(ticks)
    end
end
println()
toc()

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
@testset "PLOTUtils" begin

#=
#functions to test
genotype_heatmap2
dp_heatmap2
=#

@testset "genotype_heatmap2(x,title)" begin

plot = genotype_heatmap2(x,"test_genotype_title")
end

@testset "dp_heatmap2(x,title)" begin

plot = dp_heatmap2(x,"test_read_depth_title")
end

end


text:
# Kamil Jarosz

function simplify(ex)
	# jeżeli nie znamy wyrażenia to go nie upraszczamy
	ex
end

function simplify(ex::Expr)
	if ex.head != :call
		return ex
	end

	op = ex.args[1]
	# argumenty po uproszczeniu
	args = map(simplify, ex.args[2:end])

	if op == :+
		# nie sumujemy zer
		args = filter(a -> a != 0, args)

		# suma z 0 składników to 0
		if length(args) == 0
			return 0
		end

		# stałe sumujemy w jedną
		constant = 0
		args2 = []

		for arg in args
			if isa(arg, Number)
				constant += arg
			else
				push!(args2, arg)
			end
		end

		# jeśli stała nie jest 0 to ją dodajemy do wynikowych argumentów
		if constant != 0
			push!(args2, constant)
		end

		if length(args2) == 1
			return args2[1]
		else
			return Expr(:call, :+, args2...)
		end
	elseif op == :-
		a = args[1]
		b = args[2]

		if isa(a, Number) && isa(b, Number)
			# liczba - liczba
			return a - b
		elseif a == 0
			# 0 - b
			# tutaj b nie jest stałą bo warunek
			#   wyżej byłby prawdziwy
			return :(-$(b))
		elseif b == 0
			# a - 0
			return a
		else
			return :($(a) - $(b))
		end
	elseif op == :*
		# jedynki ignorujemy
		args = filter(v -> v != 1, args)

		if length(args) == 0
			# iloczyn 0 składników to 1
			return 1
		elseif any(v -> v == 0, args)
			# iloczyn 0*[...] to 0
			return 0
		end

		# stałe mnożymy w jedną
		constant = 1
		args2 = []

		for arg in args
			if isa(arg, Number)
				constant *= arg
			else
				push!(args2, arg)
			end
		end

		# jeśli stała nie jest 1 to ją dodajemy do wynikowych argumentów
		if constant != 1
			push!(args2, constant)
		end

		if length(args2) == 1
			return args2[1]
		else
			return Expr(:call, :*, args2...)
		end
	elseif op == :/
		a = args[1]
		b = args[2]

		# z upraszczaniem dzielenia jest taki problem, że
		#   możemy przez przypadek usunąć błąd dzielenia
		#   przez zero, dlatego generalnie nie skracamy

		if isa(a, Number) && isa(b, Number)
			# liczba / liczba
			return a/b
		else
			# to tutaj zachowamy błąd
			return :($(a)/$(b))
		end
	else
		# nieznany operator to nie upraszczamy
		return ex
	end
end

function autodiff(ex::Number, sym::Symbol)::Number
	# pochodna z liczby po symbolu to 0
	0
end

function autodiff(ex::Symbol, sym::Symbol)::Number
	if ex == sym
		# pochodna z x po x to 1
		return 1
	else
		# pochodna z symbolu po x to 0
		return 0
	end
end

function autodiff(ex::Expr, sym::Symbol)::Expr
	if ex.head != :call
		throw(ArgumentError("Unknown expression: " * string(ex)))
	end

	op = ex.args[1]
	args = ex.args[2:end]

	if op == :+ || op == :-
		# pochodna sumy to suma pochodnych
		return Expr(:call, op, map(a -> autodiff(a, sym), args)...)
	elseif op == :*
		# pochodna iloczynu
		ret = []

		for (i, a) in enumerate(args)
			# args[i] = a
			toMultiply = [autodiff(a, sym), args[1:end .!= i]...]
			push!(ret, Expr(:call, :*, toMultiply...))
		end

		return Expr(:call, :+, ret...)
	elseif op == :/
		# pochodna ilorazu
		a = args[1]
		ap = autodiff(a, sym)
		b = args[2]
		bp = autodiff(b, sym)

		numerator = :($(ap)*$(b)-$(a)*$(bp))
		denominator = :($(a)*$(b)*$(a)*$(b))
		return :($(numerator)/$(denominator))
	elseif op == :^
		# pochodna potęgi
		a = args[1]
		ap = autodiff(a, sym)
		b = args[2]
		bp = autodiff(b, sym)

		return :($(a)^$(b) * ($(b)*$(ap)/$(a) + $(bp)*log($(a))))
	else
		throw(ArgumentError("Unknown operator: " * string(op)))
	end
end

function test_diff(expr, sym)
	dexpr = autodiff(expr, sym)
	#println(dexpr)

	dexpr = simplify(dexpr)
	println(dexpr)

	println()
	#dump(dexpr)
end

test_diff(:(1/x*5+(2+2)), :x)
test_diff(:(1/3*x^8 + 12*x^2 + 7), :x)
test_diff(:(2*x*y*z + 2*y + x^x), :x)


text:
import MaximumLikelihoodProblems

import GLM
import Statistics
import Test

σ_true = 0.5
β_true = [1.0, 2.0, -1.0]

# @show σ_true
# @show σ_hat
# @show β_true
# @show β_hat

Test.@test typeof(σ_hat) == typeof(σ_true)
Test.@test σ_hat isa Real
Test.@test σ_true isa Real
Test.@test isfinite(σ_hat)
Test.@test σ_hat > 0.1
Test.@test abs(σ_hat - σ_true) < 0.03
Test.@test abs2(σ_hat - σ_true) < 1e-3

Test.@test typeof(β_hat) == typeof(β_true)
Test.@test ndims(β_hat) == ndims(β_true)
Test.@test size(β_hat) == size(β_true)
absolute_error = abs.(β_hat - β_true)
square_error = abs2.(β_hat - β_true)
absolute_error_proportional = absolute_error ./ abs.(β_true)
square_error_proportional = square_error ./ abs.(β_true)
Test.@test sum(absolute_error) < 0.030
Test.@test sum(square_error) < 1e-3
Test.@test sum(absolute_error_proportional) < 0.030
Test.@test sum(square_error_proportional) < 1e-3
Test.@test maximum(absolute_error) < 0.013
Test.@test maximum(square_error) < 1e-3
Test.@test maximum(absolute_error_proportional) < 0.013
Test.@test maximum(square_error_proportional) < 1e-3
Test.@test Statistics.mean(absolute_error) < 0.010
Test.@test Statistics.mean(square_error) < 1e-3
Test.@test Statistics.mean(absolute_error_proportional) < 0.007
Test.@test Statistics.mean(square_error_proportional) < 1e-3

gradient_vector_at_θ_hat = MaximumLikelihoodProblems.Internal.gradient_vector(transformed_gradient_problem,
                                                                              θ_hat)
hessian_matrix_at_θ_hat = MaximumLikelihoodProblems.Internal.hessian_matrix(transformed_gradient_problem,
                                                                            θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_zero(gradient_vector_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_hermitian(hessian_matrix_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_negative_definite(hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_positive_definite(-hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)

n, p = size(X)
beta_hat_ols = ( X' * X )\( X' * y )
beta_hat_mle = beta_hat_ols
y_hat_ols = X * beta_hat_ols
epsilon_hat_ols = y - y_hat_ols
S_ols = epsilon_hat_ols' * epsilon_hat_ols
sigma_2_hat_ols = (S_ols) / (n - p) # OLS estimator for σ²
sigma_2_hat_mle = (n - p) / (n) * sigma_2_hat_ols # MLE estimator for σ²
sigma_hat_ols = sqrt(sigma_2_hat_ols)
sigma_hat_mle = sqrt(sigma_2_hat_mle)
Test.@test isapprox(σ_hat, sigma_hat_ols; atol = 1e-4)
Test.@test isapprox(σ_hat, sigma_hat_mle; atol = 1e-4)
Test.@test isapprox(β_hat, beta_hat_ols)
Test.@test isapprox(β_hat, beta_hat_mle)

external_model_glm = GLM.lm(X, y)
sigma_hat_external_model_glm = GLM.dispersion(external_model_glm)
Test.@test isapprox(σ_hat, sigma_hat_external_model_glm; atol = 1e-4)
beta_hat_external_model_glm = GLM.coef(external_model_glm)
Test.@test isapprox(β_hat, beta_hat_external_model_glm)


text:
#
#  CalculateNoise.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 30/03/2020.
#
#

module CalculateNoise

# Import Julia packages
using Random
using Distributions
using LinearAlgebra
using StaticArrays

# Import local program modules

function calculateNoise!(W,nParticles,threadRNG)

    # Loop over all monomers
    for ii=1:nParticles
        W[ii] = SVector{3}(rand(threadRNG,Uniform(-1.0,1.0),3))
        magW = norm(W[ii])
        W[ii] = W[ii]*rand(threadRNG,Normal(0.0,1.0))/magW
    end

    return nothing
end

export calculateNoise!

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
using Printf: @printf, @sprintf
using Statistics: mean, std

####################
## main
####################
const N = 50 # the number of vector
x = 100(rand(N)) #
x_mean = mean(x)
x_std = std(x)
ind = sortperm(x)[end:-1:end-4]
####################

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()
#plotlyjs()

# bar
bar(x, xlims=(0.0, Float64(N+1)), ylims=(0.0,105.), legend=false, tickfont=font(12,"sans-serif"),
    title="N=$(N), mean="*@sprintf("%.2f",x_mean)*", σ="*@sprintf("%0.2f",x_std), # GRバックエンドの場合はσが正しく表示されない
    size=(800,600))

# horizontal line
xm = [0.0,Float64(N+1)]
ym = [x_mean, x_mean]
plot!(xm, ym, line=:solid, c=:green, lab="mean")

# 1st - 5th large
scatter!(ind, x[ind], ms=8., color=:red)

# save
savefig(joinpath(figdir,"rand_bar.png"))

####################


text:
#Plot a 2D-array as an image to the terminal
#Needs xterm-256color to run properly
#try to add export TERM=xterm-256color to your .bashrc
module ImageTerm
using Colors
using Images
export imageterm,sethd
const use_hd=Bool[false]

sethd(hd::Bool)=use_hd[1]=hd

function colstr(r::Integer,g::Integer,b::Integer)
  #r,g,b should be between 1 and 6
  r=min(r,6);r=max(r,1)
  g=min(g,6);g=max(g,1)
  b=min(b,6);b=max(b,1)
  return(36*(r-1)+6*(g-1)+b+15)
end
function colstr(c::AbstractRGB)
  return(colstr(ceil(Int,red(c)*6),ceil(Int,green(c)*6),ceil(Int,blue(c)*6)))
end

function stringfull(c)
  return(string("\e[48;5;",c,"m  "))
end

function stringtd(ct,cb)
  return(string("\e[48;5;",ct,"m\e[38;5;",cb,"m\u2584"))
end

abstract Normalization{T}
immutable LinearNormalization{T} <: Normalization{T}
  missval::T
  minval::T
  maxval::T
end

type ColorMapMiss{T,K}
  colors::Vector{RGB{T}}
  misscolor::RGB{T}
  missval::K
end
ColorMapMiss{T,V}(colors::Vector{RGB{T}},misscolor::RGB{V},missval)=ColorMapMiss(colors,convert(RGB{T},misscolor),missval)

getColor(c::ColorMapMiss,x::Number)=(x==c.missval) ? c.misscolor : c.colors[round(Int,x*(length(c.colors)-1))+1]
normalize(n::LinearNormalization,x::Number)= x==n.missval ? x : (x-n.minval)/(n.maxval-n.minval)
normalize!(n::LinearNormalization,x::Array)=for i in eachindex(x) x[i]=normalize(n,x[i]) end

function getPlotArray(a,npc,npl,transpose)
  mout=zeros(elemtype(a),npl,npc)
  if transpose
    f1=(size(a,1)-1)/(npc-1)
    f2=(size(a,2)-1)/(npl-1)
    for i=1:npc,j=1:npl
      mout[j,i]=readval(a,min(round(Int,(i-1)*f1)+1,size(a,1)),min(round(Int,(j-1)*f2)+1,size(a,2)))
    end
  else
    f1=(size(a,1)-1)/(npl-1)
    f2=(size(a,2)-1)/(npc-1)
    for i=1:npc,j=1:npl
      mout[j,i]=readval(a,min(round(Int,(j-1)*f1)+1,size(a,1)),min(round(Int,(i-1)*f2)+1,size(a,2)))
    end
  end

  mout
end

function getMinMax(x,missval)
  mi=typemax(eltype(x))
  ma=typemin(eltype(x))
  for ix in x
    if ix!=missval
      if ix<mi mi=ix end
      if ix>ma ma=ix end
    end
  end
  mi,ma
end

Base.one(x::UTF8String)=""

function imageterm(io,a::AbstractMatrix;col::Vector=colormap("Blues"),missval::Number=1.0e32,misscol=colorant"black",hd::Bool=use_hd[1],title::AbstractString="",subtitle::AbstractString="",transpose=false)
  ldim=size(a)
  l,c=Base.tty_size()
  legendwidth = 15
  titleheight = title=="" ? 0 : 1
  # Determine size available for the plot
  npc = div(c-legendwidth,hd ? 1 : 2)
  npl = (l-titleheight)*(hd ? 2 : 1)
  #Determine final plotsize
  npl,npc=getPlotsize(ldim[1],ldim[2],npl,npc,transpose)
  #Get the data necessary for plotting from input
  pa=getPlotArray(a,npc,npl,transpose)
  #Normalize plotting data to 0..1
  no=LinearNormalization(pa,convert(eltype(pa),missval))
  normalize!(no,pa)
  #
  cm=ColorMapMiss(col,misscol,missval)
  #
  rgbar = pa2rgb(pa,cm)
  #
  sa=rgbar2sa(rgbar,hd=hd)
  #Make the legend
  sl=get_legend(cm,size(sa,1),no)
  #Concatenate with plot
  s=reducedim(string,hcat(sa,sl),2,UTF8String(""))
  #Join title and subtitle
  vcat(get_title(title,npc),s,get_subtitle(subtitle,cm))
  print(io,join(s),"\n")
end
imageterm(a::AbstractMatrix;col::Vector=colormap("Blues"),missval::Number=1.0e32,misscol=colorant"black",hd::Bool=use_hd[1],title::AbstractString="",subtitle::AbstractString="",transpose=false)=
imageterm(STDOUT,a,col=col,missval=missval,misscol=misscol,hd=hd,title=title,subtitle=subtitle,transpose=transpose)

function imageterm(io,img::AbstractImageDirect;hd::Bool=use_hd[1])
  title=string(colorspace(img), " ", typeof(img).name)
  subtitle=string("data: ", summary(img.data))
  transpose=img.properties["spatialorder"][1]=="x" ? true : false
  imageterm(io,img.data,hd=hd,title=title,transpose=transpose)
end
imageterm(img::AbstractImageDirect;hd::Bool=use_hd[1])=imageterm(STDOUT,img,hd=hd)

function imageterm{T<:AbstractRGB}(io,a::AbstractMatrix{T};hd::Bool=use_hd[1],title::AbstractString="",subtitle::AbstractString="",transpose=false)
  ldim=size(a)
  l,c=Base.tty_size()
  titleheight = title=="" ? 0 : 1
  # Determine size available for the plot
  npc = div(c,hd ? 1 : 2)
  npl = (l-titleheight)*(hd ? 2 : 1)
  #Determine final plotsize
  npl,npc=getPlotsize(ldim[1],ldim[2],npl,npc,transpose)
  #Get the data necessary for plotting from input
  pa=getPlotArray(a,npc,npl,transpose)
  # Convert to STring representation
  sa=rgbar2sa(pa,hd=hd)
  s=reducedim(string,hcat(sa,UTF8String["\e[0m\n" for i=1:size(sa,1), j=1]),2,UTF8String(""))
  #Join title and subtitle
  vcat(get_title(title,npc),s,subtitle)
  print(io,join(s),"\n")
end
imageterm{T<:AbstractRGB}(a::AbstractMatrix{T};hd::Bool=false,title::AbstractString="",subtitle::AbstractString="",transpose=false)=
imageterm{T<:AbstractRGB}(STDOUT,a;hd=hd,title=title,subtitle=subtitle,transpose=transpose)

pa2rgb(pa,cm)=[getColor(cm,pa[irowpa,icol]) for irowpa=1:size(pa,1), icol=1:size(pa,2)]

#Convert values to String Array
function rgbar2sa(rgbar;hd=false)
  if hd==false
    sa=Array(UTF8String,size(rgbar))
    for i in eachindex(sa)
      sa[i]=stringfull(colstr(rgbar[i]))
    end
  else
    sa=Array(UTF8String,ceil(Int,size(rgbar,1)/2),size(rgbar,2))
    irowpa=1
    for irow=1:size(sa,1)
      for icol=1:size(sa,2)
        if irowpa<size(rgbar,2)
          sa[irow,icol]=stringtd(colstr(rgbar[irowpa,icol]),colstr(rgbar[irowpa+1,icol]))
        else
          sa[irow,icol]=stringtd(colstr(rgbar[irowpa,icol]),colstr(rgbar[irowpa,icol]))
        end
      end
      irowpa+=2
    end
  end
  sa
end

LinearNormalization(pa,missval)=LinearNormalization(missval,getMinMax(pa,missval)...)
get_legend(cm,npl,no::Normalization)=UTF8String[string("\e[0m ",repeat(stringfull(colstr(getColor(cm,(npl-iline)/(npl-1)))),2),"\e[0m",@sprintf("%10.2e",no.minval+(npl-iline)/(npl-1)*(no.maxval-no.minval)),"\n") for iline=1:npl]

get_title(title,npc)= title!="" ? string("\n\e[1m",^(" ",round(Int64,(npc-length(title))/2)),title,"\e[0m\n\n") : title
get_subtitle(subtitle,cm)= subtitle!="" ? "\e[0m; min=$(cm.norm.minval); max=$(cm.norm.maxval)" : subtitle

function getPlotsize(slines,scol,l,c,transpose)
  if transpose slines,scol = scol,slines end
  asprat=slines/scol
  if slines>l
    slines=l
    scol=round(Int,slines/asprat)
  end
  if scol>c
    scol=c
    slines=round(Int,scol*asprat)
  end
  return slines,scol
end

readval(a::AbstractArray,i,j)=getindex(a,i,j)
elemtype(a::AbstractArray)=eltype(a)
# Overwrite writemime for Images
Base.writemime(io::IO, ::MIME"text/plain", img::AbstractImageDirect) = imageterm(io, img)

#haskey(Pkg.installed(),"NetCDF") ? include("ncplot.jl") : nothing
end


text:
@testset "diffusion_matrix.jl" begin
    a = 0.0
    b = 2π
    N = 100

    # Test with ClosedIntervalDomain
    domain = ClosedIntervalDomain(a, b)
    D = diffusion_matrix(domain, N)
    @test D isa SparseMatrixCSC
    @test size(D) == (N + 1, N + 1)
    @test nnz(D) == 3(N + 1) - 2
    @test all(diag(D) .< 0)

    # Test with PeriodicIntervalDomain
    domain = PeriodicIntervalDomain(a, b)
    D = diffusion_matrix(domain, N)
    @test D isa SparseMatrixCSC
    @test size(D) == (N, N)
    @test nnz(D) == 3N
    @test all(diag(D, -1) .> 0)
    @test all(diag(D) .< 0)
    @test all(diag(D, +1) .> 0)
    @test D[1, end] > 0
    @test D[end, 1] > 0

    # Concrete unknown domain type
    domain = UnknownDomain()
    @test_throws Exception diffusion_matrix(domain, N)
end


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

mutable struct LogAnalyticsInputBase <: SwaggerModel
    blobContainerSasUri::Any # spec type: Union{ Nothing, String } # spec name: blobContainerSasUri
    fromTime::Any # spec type: Union{ Nothing, DateTime } # spec name: fromTime
    toTime::Any # spec type: Union{ Nothing, DateTime } # spec name: toTime
    groupByThrottlePolicy::Any # spec type: Union{ Nothing, Bool } # spec name: groupByThrottlePolicy
    groupByOperationName::Any # spec type: Union{ Nothing, Bool } # spec name: groupByOperationName
    groupByResourceName::Any # spec type: Union{ Nothing, Bool } # spec name: groupByResourceName

    function LogAnalyticsInputBase(;blobContainerSasUri=nothing, fromTime=nothing, toTime=nothing, groupByThrottlePolicy=nothing, groupByOperationName=nothing, groupByResourceName=nothing)
        o = new()
        validate_property(LogAnalyticsInputBase, Symbol("blobContainerSasUri"), blobContainerSasUri)
        setfield!(o, Symbol("blobContainerSasUri"), blobContainerSasUri)
        validate_property(LogAnalyticsInputBase, Symbol("fromTime"), fromTime)
        setfield!(o, Symbol("fromTime"), fromTime)
        validate_property(LogAnalyticsInputBase, Symbol("toTime"), toTime)
        setfield!(o, Symbol("toTime"), toTime)
        validate_property(LogAnalyticsInputBase, Symbol("groupByThrottlePolicy"), groupByThrottlePolicy)
        setfield!(o, Symbol("groupByThrottlePolicy"), groupByThrottlePolicy)
        validate_property(LogAnalyticsInputBase, Symbol("groupByOperationName"), groupByOperationName)
        setfield!(o, Symbol("groupByOperationName"), groupByOperationName)
        validate_property(LogAnalyticsInputBase, Symbol("groupByResourceName"), groupByResourceName)
        setfield!(o, Symbol("groupByResourceName"), groupByResourceName)
        o
    end
end # type LogAnalyticsInputBase

const _property_map_LogAnalyticsInputBase = Dict{Symbol,Symbol}(Symbol("blobContainerSasUri")=>Symbol("blobContainerSasUri"), Symbol("fromTime")=>Symbol("fromTime"), Symbol("toTime")=>Symbol("toTime"), Symbol("groupByThrottlePolicy")=>Symbol("groupByThrottlePolicy"), Symbol("groupByOperationName")=>Symbol("groupByOperationName"), Symbol("groupByResourceName")=>Symbol("groupByResourceName"))
const _property_types_LogAnalyticsInputBase = Dict{Symbol,String}(Symbol("blobContainerSasUri")=>"String", Symbol("fromTime")=>"DateTime", Symbol("toTime")=>"DateTime", Symbol("groupByThrottlePolicy")=>"Bool", Symbol("groupByOperationName")=>"Bool", Symbol("groupByResourceName")=>"Bool")
Base.propertynames(::Type{ LogAnalyticsInputBase }) = collect(keys(_property_map_LogAnalyticsInputBase))
Swagger.property_type(::Type{ LogAnalyticsInputBase }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LogAnalyticsInputBase[name]))}
Swagger.field_name(::Type{ LogAnalyticsInputBase }, property_name::Symbol) =  _property_map_LogAnalyticsInputBase[property_name]

function check_required(o::LogAnalyticsInputBase)
    (getproperty(o, Symbol("blobContainerSasUri")) === nothing) && (return false)
    (getproperty(o, Symbol("fromTime")) === nothing) && (return false)
    (getproperty(o, Symbol("toTime")) === nothing) && (return false)
    true
end

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


text:
module FluidScripts

end # module


text:
using Pkg
Pkg.activate(".")
using Optim
using ForwardDiff
using JLD
using LineSearches
using AdvancedMH
using AdaptiveMCMC
using MCMCChains

using Distributions
using StatsBase
using Base.Threads
nthreads()
using Printf
using Plots
using Plots.PlotMeasures
using LaTeXStrings
using Distributions
using StatsPlots

@time begin

#########################################################################################################################
# User inputs
#########################################################################################################################

# Goniometer settings (must be set before including set_up.jl file as they enter into model definitions)
TwoDeltaTheta = 0
phi = 2.5E-3
chi = 0.015*pi/180

include("set_up.jl")

# Define lattice vector
miller_vec = [0;0;1]

# Define bounds of prior prior_bounds[char,:] = (width of prior, height of prior) in nm in the lab frame
prior_bounds = [800 800; 500 1000; 500 1000; 800 800;
                1500 1500; 1500 1500; 1000 1500; 1000 1500;
                500 1200; 1000 1000; 1000 1000; 500 1200
]

# Set character of dislocation
char = 2  #parse(Int64,ARGS[1])
dis_lab = get_dis_lab(char, phi, chi)

# Set noise parameters
background_percent = 0.05 #parse(Float64,ARGS[2])  # magnitude of background noise as percentage of image intensity
noise_level = 0.01 #parse(Float64,ARGS[3])  # magnitude of electic readout noise as percentage of image intensity

# Set prior center in image frame
loc_center = [0.;0.]

# Define the prior based on the character
lower = loc_center - prior_bounds[char,:]/2
upper = loc_center + prior_bounds[char,:]/2
prior = Product(Uniform.(lower, upper))

# Define name for result files (plot generation at bottom of script)
save_str = "example_data//inference_example_char_$(char)_background_$(background_percent)_noise_$(noise_level)"
save_str = replace(save_str, "."=>"p")

# Grid sizes for computing the log-posterior contour plots
gridsize_1 = 50
gridsize_2 = 5

#########################################################################################################################
# Generate a noisy image to act as an evaluation image
#########################################################################################################################

# Draw the true location from the prior and transform it to the dislocation frame
true_loc_im = rand(prior)
true_loc_dis = image_to_dis(true_loc_im, dis_lab)

# Compute the noise-free image
im_true = compute_image(char, true_loc_dis, phi, chi, TwoDeltaTheta, miller_vec)

# Define noise model based on noise parameters
noise_std = noise_level*maximum(im_true)
noise_var = noise_std^2

background = background_percent*maximum(im_true)

pixel_var = im_true .+ background .+ noise_var

# Compute and save noisy evaluation image
im_obs = im_true .+ background + sqrt.(pixel_var).*randn(Npixels,Npixels)
save(save_str*"_im_obs.jld","im_obs", im_obs)

#########################################################################################################################
# MAP optimization and Laplace approximation
#########################################################################################################################

# Define negative log-likelihood function of dislocation position (in dislocation frame)
function nlogl(locs)

        -1*compute_loglikelihood(im_obs, noise_var, char, locs, phi, chi, TwoDeltaTheta, miller_vec, background)

end

# Define (unnormalized) log posterior through log-likelihood (the prior is uniform and so does not contribution)
logpi(locs) = -nlogl(locs)

# optimization through Optim routine
loc_opt = optimize(nlogl,
                   loc_center,
                   GradientDescent(linesearch=LineSearches.BackTracking(order=2)),
                   Optim.Options(extended_trace=true,
                                        store_trace=true,
                                        allow_f_increases=true,
                                        iterations=50,
                                        show_trace=true,
                                        g_abstol = 1E-5,
                                        x_abstol = 0.01)
)

# Unpack optimization results we want to save
num_iter = size(hcat(Optim.f_trace(loc_opt)...),2)
map_point = Optim.minimizer(loc_opt)
f_trace = hcat(Optim.f_trace(loc_opt)...)
loc_trace = hcat(Optim.x_trace(loc_opt)...)
runtime = Optim.time_run(loc_opt)

# Compute the Laplace approximation covariance (the negative inverse hessian of the posterior log-density at the MAP point)
map_cov = -inv(ForwardDiff.hessian(logpi, map_point))

# Save the data
save(save_str*"_opt_data.jld","true_loc_dis", true_loc_dis, "loc_trace", loc_trace, "f_trace", f_trace, "map_point", map_point, "num_iter", num_iter, "runtime", runtime, "map_cov", map_cov)

#########################################################################################################################
# MCMC sampling
#########################################################################################################################

# Define a distribution for drawing chain initialization points centered at the MAP point (alt. one could start chains at the MAP point)
mcmc_init = Product([Uniform(map_point[1] - 5, map_point[1] + 5),Uniform(map_point[2] - 5, map_point[2] + 5)])

# Define number and length of chains to compute [1000 length chain takes my laptop ~7-8 minutes to run]
nchains = 1
n_samples = 100
adaptive_rwhm_chains = zeros(n_samples, 2, nchains)
times = zeros(nchains)
Threads.@threads for i in 1:nchains
        θ_init = rand(mcmc_init)
        samples = @timed adaptive_rwm(θ_init, logpi, n_samples; algorithm=:asm, b=1)
        adaptive_rwhm_chains[:,:,i] = Array(samples.value.X)'
        times[i] = samples.time
end

save(save_str*"_mcmc.jld","chains", adaptive_rwhm_chains)

adaptive_rwhm_chains = load(save_str*"_mcmc.jld","chains")
all_chains = Chains(adaptive_rwhm_chains)
all_samples = Array(all_chains)

#########################################################################################################################
# Computing log-posterior contour for plotting
#########################################################################################################################

#####################################################################
# A course sampling of a region larger than the prior in image plane
#####################################################################

xs_im = 2*lower[1]:gridsize_1:2*upper[1]
ys_im = 2*lower[2]:gridsize_1:2*upper[2]
density_matrix_im = zeros(length(ys_im),length(xs_im))
dis_xyz = []
@threads for i = 1:length(ys_im)
        y = ys_im[i]

        for j = 1:length(xs_im)
                x = xs_im[j]
                loc = [x;y]

                append!(dis_xyz, image_to_dis(loc, dis_lab, dims=3))
                density_matrix_im[i,j] = nlogl(image_to_dis(loc,dis_lab))
        end
end

save(save_str*"_mat_im.jld","mat", density_matrix_im, "xs_im", xs_im, "ys_im", ys_im)

dis_xyz = reshape(dis_xyz, (3, Int(length(dis_xyz)/3)))

xs_dis = minimum(dis_xyz[1,:]):gridsize_1:maximum(dis_xyz[1,:])
ys_dis = minimum(dis_xyz[2,:]):gridsize_1:maximum(dis_xyz[2,:])

density_matrix_dis = zeros(length(ys_dis),length(xs_dis))
@threads for i = 1:length(ys_dis)
        y = ys_dis[i]

        for j = 1:length(xs_dis)
                x = xs_dis[j]
                loc = [x;y]

                density_matrix_dis[i,j] = nlogl(loc)
        end
end

save(save_str*"_mat_dis.jld","mat", density_matrix_dis, "xs_dis", xs_dis, "ys_dis", ys_dis)
#########################################################################
# Refined density data near the MAP point in the dislocation frame
#########################################################################

xs_0 = minimum(all_samples[:,1])
xs_0 = minimum([xs_0, true_loc_dis[1]]) - 10
xs_1 = maximum(all_samples[:,1])
xs_1 = maximum([xs_1, true_loc_dis[1]]) + 10
xs_mid = 0.5*(xs_1 + xs_0)
xs_length = xs_1 - xs_0

ys_0 = minimum(all_samples[:,2])
ys_0 = minimum([ys_0, true_loc_dis[2]]) - 10
ys_1 = maximum(all_samples[:,2])
ys_1 = maximum([ys_1, true_loc_dis[2]]) + 10
ys_mid = 0.5*(ys_1 + ys_0)
ys_length = ys_1 - ys_0
max_length = maximum([xs_length,ys_length])

xs_f_dis = xs_mid - max_length/2:gridsize_2:xs_mid + max_length/2
ys_f_dis = ys_mid - max_length/2:gridsize_2:ys_mid + max_length/2
density_matrix_f_dis = zeros(length(ys_f_dis),length(xs_f_dis))
@threads for i = 1:length(ys_f_dis)
        y = ys_f_dis[i]

        for j = 1:length(xs_f_dis)
                x = xs_f_dis[j]
                loc = [x;y]
                density_matrix_f_dis[i,j] = nlogl(loc)
        end
end

save(save_str*"_mat_f_dis.jld","mat", density_matrix_f_dis, "xs_f_dis", xs_f_dis, "ys_f_dis", ys_f_dis)

##############################################################
# Refined density data near the MAP point in the image plane
##############################################################

all_samples_im = reshape(dis_to_image(all_samples',dis_lab),size(all_samples'))'

xs_0 = minimum(all_samples_im[:,1])
xs_0 = minimum([xs_0, true_loc_im[1]]) - 10
xs_1 = maximum(all_samples_im[:,1])
xs_1 = maximum([xs_1, true_loc_im[1]]) + 10
xs_mid = 0.5*(xs_1 + xs_0)
xs_length = xs_1 - xs_0

ys_0 = minimum(all_samples_im[:,2])
ys_0 = minimum([ys_0, true_loc_im[2]]) - 10
ys_1 = maximum(all_samples_im[:,2])
ys_1 = maximum([ys_1, true_loc_im[2]]) + 10
ys_mid = 0.5*(ys_1 + ys_0)
ys_length = ys_1 - ys_0
max_length = maximum([xs_length,ys_length])

xs_f_im = xs_mid - max_length/2:gridsize_2:xs_mid + max_length/2
ys_f_im = ys_mid - max_length/2:gridsize_2:ys_mid + max_length/2
density_matrix_f_im = zeros(length(ys_f_im),length(xs_f_im))
@threads for i = 1:length(ys_f_im)
        y = ys_f_im[i]

        for j = 1:length(xs_f_im)
                x = xs_f_im[j]
                loc = [x;y]
                density_matrix_f_im[i,j] = nlogl(image_to_dis(loc,dis_lab))
        end
end

save(save_str*"_mat_f_im.jld","mat", density_matrix_f_im, "xs_f_im", xs_f_im, "ys_f_im", ys_f_im)

# ############################################################################################################################
# # Plotting
# ############################################################################################################################

include("plotting.jl")
load_name = save_str

figure_dir = "example_figures"

# Load data
im_obs = load(load_name*"_im_obs.jld","im_obs")
d = load(load_name*"_opt_data.jld")
true_loc_dis = d["true_loc_dis"]
true_loc_im = dis_to_image(true_loc_dis, dis_lab)
loc_trace_dis =d["loc_trace"]
loc_trace_im = reshape(dis_to_image(loc_trace_dis,dis_lab),size(loc_trace_dis))
f_trace =d["f_trace"]
map_point_dis =d["map_point"]
map_point_im = dis_to_image(map_point_dis,dis_lab)
num_iter = d["num_iter"]
runtime = d["runtime"]
map_cov = d["map_cov"]

adaptive_rwhm_chains = load(load_name*"_mcmc.jld","chains")
all_chains = Chains(adaptive_rwhm_chains)
all_samples = Array(all_chains)
all_samples_im = reshape(dis_to_image(all_samples',dis_lab),size(all_samples'))'

d = load(load_name*"_mat_im.jld")
xs_im = d["xs_im"]
ys_im = d["ys_im"]
density_matrix_im = d["mat"]

d = load(load_name*"_mat_dis.jld")
xs_dis = d["xs_dis"]
ys_dis = d["ys_dis"]
density_matrix_dis = d["mat"]

d = load(load_name*"_mat_f_dis.jld")
xs_f_dis = d["xs_f_dis"]
ys_f_dis = d["ys_f_dis"]
density_matrix_f_dis = d["mat"]

d = load(load_name*"_mat_f_im.jld")
xs_f_im = d["xs_f_im"]
ys_f_im = d["ys_f_im"]
density_matrix_f_im = d["mat"]

# Plot evaluation image
dfxm_image_plot_base(im_obs, char, true_loc_dis, figure_dir*"\\obs", prior=true, width=500)

# Figure 5c
contourf(xs_im, ys_im, density_matrix_im, levels=50, aspect_ratio=:equal, color=:reds)
rect = prior_bounds[char,:]
plot!(rectangle(rect...), opacity=.5, c=:gray, label=false, legend=:bottomright)
for t = 1:num_iter - 1
        plot!([loc_trace_im[1,t], loc_trace_im[1,t+1]], [loc_trace_im[2,t], loc_trace_im[2,t+1]],label=false,linewidth=3,c=:black)
end
scatter!(loc_trace_im[1,:],loc_trace_im[2,:], c=:orange, markershape = :utriangle, markersize = 4, label="iterates")
scatter!([map_point_im[1]],[map_point_im[2]],c=:lightgreen, markershape = :star5, markersize = 8, label=L"\xi^{\mathrm{MAP}}")
scatter!([true_loc_im[1]],[true_loc_im[2]], c=:cyan, markershape = :star7, markersize = 8, label=L"\xi^{\mathrm{true}}", legend=false)# outertop)
xlims!(-500,500)
ylims!(-600,600)
xticks!([-500,0,500])
yticks!([-500,0,500])
xlabel!(L"y_{\ell} \ (\mathrm{nm})", xtickfontsize=12, xguidefontsize=18, margin=5mm)
ylabel!(L"x_{\ell} \ (\mathrm{nm})", ytickfontsize=12, yguidefontsize=18)
plot!(size=(450,450),right_margin=15mm)
Plots.savefig(figure_dir*"\\iters.pdf")

# Figure 5b
prior_corners = [250; 500; -250; 500; -250; -500; 250; -500]
prior_corners_dis = reshape(image_to_dis(prior_corners, dis_lab), (2,4))
s_prior = Shape(prior_corners_dis[1,:],prior_corners_dis[2,:])

contourf(xs_dis, ys_dis, density_matrix_dis, levels=10, aspect_ratio=:equal, color=cgrad(:grays,rev=true))
plot!(s_prior, opacity=.5, c=:gray, linecolor=:white, label=false)
#plot!(s_region, c=false, linecolor=:white, label=false)
for t = 1:num_iter - 1
        plot!([loc_trace_dis[1,t], loc_trace_dis[1,t+1]], [loc_trace_dis[2,t], loc_trace_dis[2,t+1]],label=false,linewidth=3,c=:black)
end
scatter!(loc_trace_dis[1,:],loc_trace_dis[2,:], c=:orange, markershape = :utriangle, markersize = 4, label="iterates")
scatter!([map_point_dis[1]],[map_point_dis[2]],c=:lightgreen, markershape = :star5, markersize = 8, label=L"\xi^{\mathrm{MAP}}")
scatter!([true_loc_dis[1]],[true_loc_dis[2]], c=:cyan, markershape = :star7, markersize = 8, label=L"\xi^{\mathrm{true}}", legend=false)# outertop)
xlabel!(L"x_{d} \mathrm{ \ (nm)}", xtickfontsize=12, xguidefontsize=18, margin=5mm)
ylabel!(L"y_{d} \mathrm{ \ (nm)}", ytickfontsize=12, yguidefontsize=18)
xlims!(xs_dis[1],xs_dis[end])
ylims!(ys_dis[1],ys_dis[end])
xticks!([-500,0,500])
yticks!([-400,0,400])
plot!(size=(450,450),right_margin=15mm)
Plots.savefig(figure_dir*"\\iters_dis.pdf")

# Figure 6a
contourf(xs_f_dis,ys_f_dis,density_matrix_f_dis, levels=50, aspect_ratio=:equal, color=cgrad(:grays,rev=true))
scatter!(all_samples[:,1],all_samples[:,2],c=:plum1, alpha=.1, label="MCMC samples")
covellipse!(vec(map_point_dis), map_cov, n_std=2, aspect_ratio=1, color=(:yellow), label=false, alpha=.5)
scatter!([map_point_dis[1]],[map_point_dis[2]],c=:lightgreen, markershape = :star5, markersize = 8, label=L"\xi^{\mathrm{MAP}}")
scatter!([true_loc_dis[1]],[true_loc_dis[2]], c=:cyan, markershape = :star7, markersize = 8, label=L"\xi^{\mathrm{true}}", legend=false)# outertop)
xlabel!(L"x_{d} \mathrm{ \ (nm)}", xtickfontsize=14, xguidefontsize=18, margin=5mm)
ylabel!(L"y_{d} \mathrm{ \ (nm)}", ytickfontsize=14, yguidefontsize=18)
plot!(size=(450,450), right_margin=15mm)
Plots.savefig(figure_dir*"\\uq_dis.pdf")

# Figure 6b
vals,vecs = eigen(map_cov)
vecs_im = reshape(dis_to_image(vecs,dis_lab),(2,2))
map_cov_im = vecs_im*diagm(vals)*vecs_im'

contourf(xs_f_im,ys_f_im,density_matrix_f_im, levels=50,aspect_ratio=:equal, color=:reds)
scatter!(all_samples_im[:,1],all_samples_im[:,2], alpha=.2, c=:plum1, label="MCMC samples")
covellipse!(vec(map_point_im), map_cov_im, n_std=2, aspect_ratio=1, color=(:yellow), label=false, alpha=.5)
scatter!([map_point_im[1]],[map_point_im[2]],c=:lightgreen, markershape = :star5, markersize = 8, label=L"\xi^{\mathrm{MAP}}")
scatter!([true_loc_im[1]],[true_loc_im[2]], c=:cyan, markershape = :star7, markersize = 8, label=L"\xi^{\mathrm{true}}", legend=false)
xlabel!(L"y_{\ell} \mathrm{ \ (nm)}", xtickfontsize=14, xguidefontsize=18, margin=5mm)
ylabel!(L"x_{\ell} \mathrm{ \ (nm)}", ytickfontsize=14, yguidefontsize=18)
plot!(size=(450,450), right_margin=15mm)
Plots.savefig(figure_dir*"\\uq_im.pdf")

end

max_stars_repo_name: michael-c-brennan/DFXMTools

max_issues_repo_name: michael-c-brennan/DFXMTools

max_forks_repo_name: michael-c-brennan/DFXMTools


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

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec <: SwaggerModel
    conversion::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceConversion } # spec name: conversion
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    names::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames } # spec name: names
    preserveUnknownFields::Any # spec type: Union{ Nothing, Bool } # spec name: preserveUnknownFields
    scope::Any # spec type: Union{ Nothing, String } # spec name: scope
    versions::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion} } # spec name: versions

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec(;conversion=nothing, group=nothing, names=nothing, preserveUnknownFields=nothing, scope=nothing, versions=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("conversion"), conversion)
        setfield!(o, Symbol("conversion"), conversion)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("names"), names)
        setfield!(o, Symbol("names"), names)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("preserveUnknownFields"), preserveUnknownFields)
        setfield!(o, Symbol("preserveUnknownFields"), preserveUnknownFields)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("scope"), scope)
        setfield!(o, Symbol("scope"), scope)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec, Symbol("versions"), versions)
        setfield!(o, Symbol("versions"), versions)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec = Dict{Symbol,Symbol}(Symbol("conversion")=>Symbol("conversion"), Symbol("group")=>Symbol("group"), Symbol("names")=>Symbol("names"), Symbol("preserveUnknownFields")=>Symbol("preserveUnknownFields"), Symbol("scope")=>Symbol("scope"), Symbol("versions")=>Symbol("versions"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec = Dict{Symbol,String}(Symbol("conversion")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceConversion", Symbol("group")=>"String", Symbol("names")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionNames", Symbol("preserveUnknownFields")=>"Bool", Symbol("scope")=>"String", Symbol("versions")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion}")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionSpec)
    (getproperty(o, Symbol("group")) === nothing) && (return false)
    (getproperty(o, Symbol("names")) === nothing) && (return false)
    (getproperty(o, Symbol("scope")) === nothing) && (return false)
    (getproperty(o, Symbol("versions")) === nothing) && (return false)
    true
end

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


text:
using PerfChecker
using Documenter

DocMeta.setdocmeta!(PerfChecker, :DocTestSetup, :(using PerfChecker); recursive=true)

makedocs(;
    modules=[PerfChecker],
    authors="Azzaare <jf@baffier.fr>",
    repo="https://github.com/JuliaConstraints/PerfChecker.jl/blob/{commit}{path}#{line}",
    sitename="PerfChecker.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaConstraints.github.io/PerfChecker.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/PerfChecker.jl",
    devbranch="main",
)


text:
# # Saving SCF results on disk and SCF checkpoints
#
# For longer DFT calculations it is pretty standard to run them on a cluster
# in advance and to perform postprocessing (band structure calculation,
# plotting of density, etc.) at a later point and potentially on a different
# machine.
#
# To support such workflows DFTK offers the two functions [`save_scfres`](@ref)
# and [`load_scfres`](@ref), which allow to save the data structure returned
# by [`self_consistent_field`](@ref) on disk or retrieve it back into memory,
# respectively. For this purpose DFTK uses the
# [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) file format and Julia package.
# For the moment this process is considered an experimental feature and
# has a number of caveats, see the warnings below.
#
# !!! warning "Saving `scfres` is experimental"
#     The [`load_scfres`](@ref) and [`save_scfres`](@ref) pair of functions
#     are experimental features. This means:
#
#     - The interface of these functions
#       as well as the format in which the data is stored on disk can
#       change incompatibly in the future. At this point we make no promises ...
#     - JLD2 is not yet completely matured
#       and it is recommended to only use it for short-term storage
#       and **not** to archive scientific results.
#     - If you are using the functions to transfer data between different
#       machines ensure that you use the **same version of Julia, JLD2 and DFTK**
#       for saving and loading data.
#
# To illustrate the use of the functions in practice we will compute
# the total energy of the O₂ molecule at PBE level. To get the triplet
# ground state we use a collinear spin polarisation
# (see [Collinear spin and magnetic systems](@ref) for details)
# and a bit of temperature to ease convergence:

using DFTK
using LinearAlgebra
using JLD2

d = 2.079  # oxygen-oxygen bondlength
a = 9.0    # size of the simulation box
lattice = diagm(a * ones(3))
O = ElementPsp(:O, psp=load_psp("hgh/pbe/O-q6.hgh"))
atoms = [O => d / 2a * [[0, 0, 1], [0, 0, -1]]]
magnetic_moments = [O => [1., 1.]]

Ecut  = 10  # Far too small to be converged
model = model_PBE(lattice, atoms, temperature=0.02, smearing=smearing=Smearing.Gaussian(),
                  magnetic_moments=magnetic_moments)
basis = PlaneWaveBasis(model; Ecut, kgrid=[1, 1, 1])

scfres = self_consistent_field(basis, tol=1e-2, ρ=guess_density(basis, magnetic_moments))
save_scfres("scfres.jld2", scfres);
#-
scfres.energies

# The `scfres.jld2` file could now be transfered to a different computer,
# Where one could fire up a REPL to inspect the results of the above
# calculation:

using DFTK
using JLD2
loaded = load_scfres("scfres.jld2")
propertynames(loaded)
#-
loaded.energies

# Since the loaded data contains exactly the same data as the `scfres` returned by the
# SCF calculation one could use it to plot a band structure, e.g.
# `plot_bandstructure(load_scfres("scfres.jld2"))` directly from the stored data.

# ## Checkpointing of SCF calculations
# A related feature, which is very useful especially for longer calculations with DFTK
# is automatic checkpointing, where the state of the SCF is periodically written to disk.
# The advantage is that in case the calculation errors or gets aborted due
# to overrunning the walltime limit one does not need to start from scratch,
# but can continue the calculation from the last checkpoint.
#
# To enable automatic checkpointing in DFTK one needs to pass the `ScfSaveCheckpoints`
# callback to [`self_consistent_field`](@ref), for example:

callback = DFTK.ScfSaveCheckpoints()
scfres = self_consistent_field(basis; ρ=guess_density(basis, magnetic_moments),
                               tol=1e-2, callback=callback);

# Notice that using this callback makes the SCF go silent since the passed
# callback parameter overwrites the default value (namely `DefaultScfCallback()`)
# which exactly gives the familiar printing of the SCF convergence.
# If you want to have both (printing and checkpointing) you need to chain
# both callbacks:
callback = DFTK.ScfDefaultCallback() ∘ DFTK.ScfSaveCheckpoints(keep=true)
scfres = self_consistent_field(basis; ρ=guess_density(basis, magnetic_moments),
                               tol=1e-2, callback=callback);

# For more details on using callbacks with DFTK's `self_consistent_field` function
# see [Monitoring self-consistent field calculations](@ref).

# By default checkpoint is saved in the file `dftk_scf_checkpoint.jld2`, which is
# deleted automatically once the SCF completes successfully. If one wants to keep
# the file one needs to specify `keep=true` as has been done in the ultimate SCF
# for demonstration purposes: now we can continue the previous calculation
# from the last checkpoint as if the SCF had been aborted.
# For this one just loads the checkpoint with [`load_scfres`](@ref):

oldstate = load_scfres("dftk_scf_checkpoint.jld2")
scfres   = self_consistent_field(oldstate.basis, ρ=oldstate.ρ,
                                 ψ=oldstate.ψ, tol=1e-3);

# !!! note "Availability of `load_scfres`, `save_scfres` and `ScfSaveCheckpoints`"
#     As JLD2 is an optional dependency of DFTK these three functions are only
#     available once one has *both* imported DFTK and JLD2 (`using DFTK`
#     and `using JLD2`).

# (Cleanup files generated by this notebook)
rm("dftk_scf_checkpoint.jld2")
rm("scfres.jld2")


text:
struct LinearChainConditionalRandomField{Tv <: AbstractFloat}
    featurenames::Vector{String}
    scores::Vector{Tv}
    matrix::Matrix{Any}
end

function LinearChainConditionalRandomField([
    ["Y-1", "Y", "W"],
    ["Y-1", "Y", "P"],
    ["Y-1", "Y", "P-1"],
    ["Y", "W"],
    ["Y", "D"],
    ])
    fnames = ["Y-1", "Y", "P-1", "P", "P-2", "W-1", "W", "W+1", "D"]
    scores = fill(0., 0)

end

function display(model::LinearChainConditionalRandomField)
    nr = size(model.matrix, 1)
    nc = size(model.matrix, 2)
    rows = Vector{Any}(undef, nr)
    for i = 1:nr
        rows[i] = map(x -> ismissing(x) ? "" : string(x), model.matrix[i, :])
        pushfirst!(rows[i], model.featurenames[i])
    end
    pushfirst!(rows, Any["FN\\Score", model.scores...])
    return display(Markdown.MD(Markdown.Table(rows, [:r, fill(:c, nc)...])))
end

function LinearChainConditionalRandomField(
    edgeobservations::Vector{Function},
    nodeobservations::Vector{Function},
    observationfuncs::Vector{Function}
    )
    body
end
function inference(model::LinearChainConditionalRandomField)

end

function estimation(model::LinearChainConditionalRandomField)

end

function edgeobs1(tag::Vector{Int}, poswords::Matrix{String}, i::Int)
    return (tag[i - 1], tag[i], poswords[i])
end


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

# List of normalised generators for the six one-dimensional
# spaces of cusp forms of level 1.

using Nemo

R, q = PowerSeriesRing(ZZ, 20, "q")

GENERATOR = Dict(
  12 => q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + 401856*q^14 + 1217160*q^15 + 987136*q^16 - 6905934*q^17 + 2727432*q^18 + 10661420*q^19 + O(q^20),
  16 => q + 216*q^2 - 3348*q^3 + 13888*q^4 + 52110*q^5 - 723168*q^6 + 2822456*q^7 - 4078080*q^8 - 3139803*q^9 + 11255760*q^10 + 20586852*q^11 - 46497024*q^12 - 190073338*q^13 + 609650496*q^14 - 174464280*q^15 - 1335947264*q^16 + 1646527986*q^17 - 678197448*q^18 + 1563257180*q^19 + O(q^20),
  18 => q - 528*q^2 - 4284*q^3 + 147712*q^4 - 1025850*q^5 + 2261952*q^6 + 3225992*q^7 - 8785920*q^8 - 110787507*q^9 + 541648800*q^10 - 753618228*q^11 - 632798208*q^12 + 2541064526*q^13 - 1703323776*q^14 + 4394741400*q^15 - 14721941504*q^16 - 5429742318*q^17 + 58495803696*q^18 + 1487499860*q^19 + O(q^20),
  20 => q + 456*q^2 + 50652*q^3 - 316352*q^4 - 2377410*q^5 + 23097312*q^6 - 16917544*q^7 - 383331840*q^8 + 1403363637*q^9 - 1084098960*q^10 - 16212108*q^11 - 16023861504*q^12 + 50421615062*q^13 - 7714400064*q^14 - 120420571320*q^15 - 8939761664*q^16 + 225070099506*q^17 + 639933818472*q^18 - 1710278572660*q^19 + O(q^20),
  22 => q - 288*q^2 - 128844*q^3 - 2014208*q^4 + 21640950*q^5 + 37107072*q^6 - 768078808*q^7 + 1184071680*q^8 + 6140423133*q^9 - 6232593600*q^10 - 94724929188*q^11 + 259518615552*q^12 - 80621789794*q^13 + 221206696704*q^14 - 2788306561800*q^15 + 3883087691776*q^16 + 3052282930002*q^17 - 1768441862304*q^18 - 7920788351740*q^19 + O(q^20),
  26 => q - 48*q^2 - 195804*q^3 - 33552128*q^4 - 741989850*q^5 + 9398592*q^6 + 39080597192*q^7 + 3221114880*q^8 - 808949403027*q^9 + 35615512800*q^10 + 8419515299052*q^11 + 6569640870912*q^12 - 81651045335314*q^13 - 1875868665216*q^14 + 145284580589400*q^15 + 1125667983917056*q^16 - 2519900028948078*q^17 + 38829571345296*q^18 - 6082056370308940*q^19 + O(q^20))


text:
#define custom Distances Metric...
#see discussion of req methods here: https://github.com/JuliaStats/Distances.jl/issues/95

struct CorrectedKmer <: Metric
    k::Int
end

correctedkmer(a::AbstractArray, b::AbstractArray, k::Int) = CorrectedKmer(k)(a, b)
correctedkmer(a::Number, b::Number, k::Int) = CorrectedKmer(k)(a, b)

#looking back at RAD paper not sure about this... check!
function evaluate(dist::CorrectedKmer,a,b)
    D = sqeuclidean(a,b)
    #dist.N*log(1 -  minimum([D / (dist.N * 2),1])) / - dist.k
    D / (dist.k * (sum(a) + sum(b)))
end

@eval @inline (dist::CorrectedKmer)(a::AbstractArray, b::AbstractArray) = evaluate(dist, a, b)
@eval @inline (dist::CorrectedKmer)(a::Number, b::Number) = evaluate(dist, a, b)

#top level function
function sequmap(seqs::Array{String,1}, ndim::Int;
    k = 5,
    lookup_dic = NT_DICT,
    pca = true,
    pca_maxoutdim = 5,
    n_neighbors = 12,
    min_dist = 0.7,
    repulsion_strength = 0.1,
    metric = SqEuclidean(),
    umap_kwargs = Pair{Symbol,Any}[] #can also put all umap args together
    )
    vecs = []
    missing_chars = Set{Char}()
    for seq in seqs
        vec, missing_chars = kmer_embed(seq, k, kmer_count!; lookup_dic = lookup_dic, missing_chars = missing_chars);
        push!(vecs, vec)
    end
    if length(missing_chars) > 0
        @warn "Ignored the following characters missing from lookup: $(unique(missing_chars)). Check your sequence type!"
    end
    X = hcat(vecs...);
    X = convert(Array{Float64,2}, X)
    if pca
        M = fit(PCA, X; maxoutdim = pca_maxoutdim)
        PCAembedding = transform(M, X);
        proj = umap(PCAembedding, ndim; n_neighbors = n_neighbors, min_dist = min_dist, repulsion_strength = repulsion_strength, metric = metric, umap_kwargs...)
    else
        proj = umap(X, ndim; n_neighbors = n_neighbors, min_dist = min_dist, repulsion_strength = repulsion_strength, metric = metric, umap_kwargs...)
    end
    return proj
end

#top level function
function seqpca(seqs::Array{String,1}, ncomponents::Int;
    k = 5,
    lookup_dic = NT_DICT)
    vecs = []
    missing_chars = Set{Char}()
    for seq in seqs
        vec, missing_chars = kmer_embed(seq, k, kmer_count!; lookup_dic = lookup_dic, missing_chars = missing_chars);
        push!(vecs, vec)
    end
    if length(missing_chars) > 0
        @warn "Ignored the following characters missing from lookup: $(unique(missing_chars)). Check your sequence type!"
    end
    X = hcat(vecs...);
    X = convert(Array{Float64,2}, X)
    M = fit(PCA, X; maxoutdim = ncomponents)
    return transform(M, X)
end


text:
module MaxHelpingHandStrawberryCollectionField

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/StrawberryCollectionField" StrawberryCollectionField(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    delayBetweenBerries::Bool=true, includeGoldens::Bool=false)

const placements = Ahorn.PlacementDict(
    "Strawberry Collection Field (max480's Helping Hand)" => Ahorn.EntityPlacement(
        StrawberryCollectionField,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# Set up Problem
model = Dynamics.car
n,m = model.n, model.m
N = 51
Q = Diagonal(I,n)*0.01
R = Diagonal(I,m)*0.01
Qf = Diagonal(I,n)*0.01
xf = [0,1,0]
obj = LQRObjective(Q,R,Qf,xf,N)
n,m = model.n, model.m
bnd = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf], u_min=[0.1,-2], u_max=1.5)
bnd1 = BoundConstraint(n,m, u_min=bnd.u_min)
bnd_x = BoundConstraint(n,m, x_min=[-0.5, -0.01, -Inf], x_max=[0.5, 1.01, Inf])
goal = goal_constraint(xf)
obs = (([0.2, 0.6], 0.25),
       ([-0.5, 0.5], 0.4))
obs1 = TO.planar_obstacle_constraint(n,m, obs[1]..., :obstacle1)
obs2 = TO.planar_obstacle_constraint(n,m, obs[2]..., :obstacle2)
con = Constraints(N)
con[1] += bnd1
for k = 2:N-1
    con[k] += bnd  + obs1 + obs2
end
con[N] += goal
prob = Problem(rk4(model), obj, constraints=con, tf=3)

# Solve with ALTRO
initial_controls!(prob, ones(m,N-1))
ilqr = iLQRSolverOptions()
al = AugmentedLagrangianSolverOptions(opts_uncon=ilqr)
al.constraint_tolerance = 1e-3
al.constraint_tolerance_intermediate = 1e-1
solve!(prob, al)
max_violation(prob)

# Test Primal Dual variable
opts = ProjectedNewtonSolverOptions{Float64}(verbose=true)
solver = ProjectedNewtonSolver(prob,opts)
NN = length(solver.V.Z)
p = num_constraints(prob)
P = sum(p) + N*n

V = copy(solver.V)
Z = TO.primals(V)
@test V.X[1] == prob.X[1]
V.X[1][1] = 100
@test V[1] == 100
Z[2] = 150
@test V.X[1][2] == 150
@test V[2] == 150
Z .+= 1
@test V[2] == 151
Y = TO.duals(V)
@test length(Y) == P

# Create PN Solver
solver = ProjectedNewtonSolver(prob,opts)
NN = N*n + (N-1)*m
p = num_constraints(prob)
P = N*n + sum(p)

# Test functions
TO.dynamics_constraints!(prob, solver)
TO.update_constraints!(prob, solver)
TO.active_set!(prob, solver)
@test all(solver.a.primals)
TO.dynamics_jacobian!(prob, solver)
@test solver.∇F[1].xx == solver.Y[1:n,1:n]
@test solver.∇F[2].xx == solver.Y[n .+ (1:n),1:n]

TO.constraint_jacobian!(prob, solver)
@test solver.∇C[1] == solver.Y[2n .+ (1:p[1]), 1:n+m]

TO.cost_expansion!(prob, solver)

# Check cost gradient and hessian
solver = ProjectedNewtonSolver(prob,opts)
V0 = copy(solver.V.V)
TO.PrimalDual(V0,prob)
TO.cost_expansion!(prob, solver)

function evalcost(V)
    cost(prob, TO.PrimalDual(V,prob))
end
@test evalcost(V0) == cost(prob, solver.V)
@test ForwardDiff.gradient(evalcost, V0)[1:NN] ≈ solver.g
@test solver.H[1:n,1:n] == Q*prob.dt
@test solver.H[n .+ (1:m), n .+ (1:m)] == R*prob.dt
@test solver.H[end-n+1:end, end-n+1:end] == Qf

# Test Constraint Violation
solver = ProjectedNewtonSolver(prob,opts)
solver.opts.active_set_tolerance = 0.0
TO.dynamics_constraints!(prob, solver)
TO.update_constraints!(prob, solver)
TO.dynamics_jacobian!(prob, solver)
TO.constraint_jacobian!(prob, solver)
TO.active_set!(prob, solver)
Y,y = TO.active_constraints(prob, solver)
viol = TO.calc_violations(solver)
@test maximum(maximum.(viol)) == norm(y,Inf)
@test norm(y,Inf) == max_violation(prob)
@test max_violation(solver) == max_violation(prob)

# Test Projection
solver = ProjectedNewtonSolver(prob,opts)
solver.opts.feasibility_tolerance = 1e-10
solver.opts.active_set_tolerance = 1e-3
TO.update!(prob, solver)
Y,y = TO.active_constraints(prob, solver)
TO.projection!(prob, solver)
TO.update!(prob, solver, solver.V)
@test TO.max_violation(solver) < solver.opts.feasibility_tolerance
res0 = norm(TO.residual(prob, solver))
res,δλ = TO.multiplier_projection!(prob, solver)
@test res < res0

solver = ProjectedNewtonSolver(prob,opts)
λ0 = copy(solver.V.Y)
V0 = copy(solver.V.V)
solver.opts.feasibility_tolerance = 1e-10
solver.opts.active_set_tolerance = 1e-3
TO.update!(prob, solver)
Y,y = TO.active_constraints(prob, solver)
TO.primaldual_projection!(prob, solver)
solver.V.Y

TO.update!(prob, solver, solver.V)
@test TO.max_violation(solver) < solver.opts.feasibility_tolerance
res0 = norm(TO.residual(prob, solver))
res, = TO.multiplier_projection!(prob, solver)
@test res < res0
TO.duals(solver.V)[solver.a.duals]
typeof(cholesky(Symmetric(Y*Hinv*Y')))

# Build KKT
Hinv = inv(Diagonal(Array(solver.H)))
TO.cost_expansion!(prob, solver)
J0 = cost(prob, V)
res0 = norm(TO.residual(prob, solver))
@test res0 == res
viol0 = max_violation(solver)
TO.cost_expansion!(prob, solver)
δV, = TO.solveKKT_Shur(prob, solver, Hinv)
V_ = solver.V + 0.5*δV
TO.projection!(prob, solver, V_)
TO.cost_expansion!(prob, solver, V_)
TO.multiplier_projection!(prob, solver, V_)

V_ = TO.line_search(prob, solver, δV)
res = norm(TO.residual(prob, solver, V_))
@test res < res0

# Test Newton Step
solver = ProjectedNewtonSolver(prob,opts)
solver.opts.feasibility_tolerance = 1e-10
solver.opts.verbose = true
V_ = TO.newton_step!(prob, solver)
TO.update!(prob, solver, V_)
@test max_violation(solver) < 1e-10
@test cost(prob, V_) < cost(prob, solver.V)
cost(prob, V_) - cost(prob, solver.V)
@test norm(TO.residual(prob, solver, V_)) < 1e-4

using Plots

plot()
TO.plot_circle!(obs[1]...)
TO.plot_circle!(obs[2]...)
plot_trajectory!(prob.X)
plot_trajectory!(V_.X)

plot(prob.U)
plot!(V_.U)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
function position_calc(x...)
    _position_calc(filter(x->!isa(x, Void), x)...)
end
function glsllinspace(position::Grid, gi, index)
    """
    (((float(position.dims[$gi])-($(index)+1)) *
        position.minimum[$gi] + $(index)*position.maximum[$gi]) *
        position.multiplicator[$gi])
    """
end
function glsllinspace(grid::Grid{1}, gi, index)
    """
    (((float(position.dims)-($(index)+1)) *
        position.minimum + $(s)*position.maximum) *
        position.multiplicator)
    """
end
function grid_pos(grid::Grid{1})
    "$(glsllinspace(grid, 0, "index"))"
end
function grid_pos(grid::Grid{2})
    "vec2($(glsllinspace(grid, 0, "index2D.x")), $(glsllinspace(grid, 1, "index2D.y")))"
end
function grid_pos(grid::Grid{3})
    "vec3(
        $(glsllinspace(grid, 0, "index2D.x")),
        $(glsllinspace(grid, 1, "index2D.y")),
        $(glsllinspace(grid, 2, "index2D.z"))
    )"
end

function _position_calc{T<:AbstractFloat}(
        grid::Grid{2}, position_z::MatTypes{T}, target::Type{Texture}
    )
    """
    ivec2 index2D = ind2sub(position.dims, index);
    vec2 normalized_index = vec2(index2D) / vec2(position.dims);
    float height = texture(position_z, normalized_index+(offset/vec2(position.dims))).x;
    pos = vec3($(grid_pos(grid)), height);
    """
end

function _position_calc{T<:AbstractFloat}(
        position_x::MatTypes{T}, position_y::MatTypes{T}, position_z::MatTypes{T}, target::Type{Texture}
    )
"""
    ivec2 index2D = ind2sub(dims, index);
    vec2 normalized_index = vec2(index2D) / vec2(dims);
    vec2 offsetted_index = normalized_index + (offset/vec2(dims));
    pos = vec3(
        texture(position_x, offsetted_index).x,
        texture(position_y, offsetted_index).x,
        texture(position_z, offsetted_index).x
    );
"""
end

function _position_calc{T<:AbstractFloat}(
        position_x::VecTypes{T}, position_y::T, position_z::T, target::Type{TextureBuffer}
    )
    "pos = vec3(texelFetch(position_x, index).x, position_y, position_z);"
end
function _position_calc{T<:AbstractFloat}(
        position_x::VecTypes{T}, position_y::T, position_z::T, target::Type{GLBuffer}
    )
    "pos = vec3(position_x, position_y, position_z);"
end
function _position_calc{T<:FixedVector}(
        position_xyz::VecTypes{T}, target::Type{TextureBuffer}
    )
    "pos = texelFetch(position, index).xyz;"
end
function _position_calc{T<:FixedVector}(
        position_xyz::VecTypes{T}, target::Type{GLBuffer}
    )
    len = length(T)
    filler = join(ntuple(x->0, 3-len), ", ")
    needs_comma = len != 3 ? ", " : ""
    "pos = vec3(position $needs_comma $filler);"
end
function _position_calc{T<:AbstractFloat}(
        position_x::VecTypes{T}, position_y::VecTypes{T}, position_z::VecTypes{T},
        target::Type{TextureBuffer}
    )
    "pos = vec3(
        texelFetch(position_x, index).x,
        texelFetch(position_y, index).x,
        texelFetch(position_z, index).x
    );"
end
function _position_calc{T<:AbstractFloat}(
        position_x::VecTypes{T}, position_y::VecTypes{T}, position_z::VecTypes{T},
        target::Type{GLBuffer}
    )
    "pos = vec3(
        position_x,
        position_y,
        position_z
    );"
end
function _position_calc(
        position::Grid{1}, target
    )
    "
    pos = vec3($(grid_pos(position)), 0, 0);
    "
end
function _position_calc(
        position::Grid{2}, target
    )
    "
    ivec2 index2D = ind2sub(position.dims, index);
    pos = vec3($(grid_pos(position)), 0);
    "
end
function _position_calc{T}(
        position::Grid{2}, ::VecTypes{T}, target::Type{GLBuffer}
    )
    "
    ivec2 index2D = ind2sub(position.dims, index);
    pos = vec3($(grid_pos(position)), position_z);
    "
end
function _position_calc(
        position::Grid{3}, target
    )
    "
    ivec3 index2D = ind2sub(position.dims, index);
    pos = $(grid_pos(position));
    "
end


text:
include("LocalField/pAdic.jl")
include("LocalField/qAdic.jl")
include("LocalField/Ring.jl")
include("LocalField/Conjugates.jl")
include("LocalField/Poly.jl")


text:
# Returns Magnus and Neudecker's duplication matrix of size n
function duplication(n::Int64)
  a = tril(ones(n,n))
  i = find(a)
  a[i] = 1:size(i,1)
  a = a + tril(a,-1)'
  j = trunc(Integer, vec(a))
  m = (n*(n+1)/2)
  m = trunc(Integer,m)
  d = zeros(n*n,m)
  for r = 1:size(d,1)
      d[r, j[r]] = 1
  end
  return d
end


text:
using Test, LinearAlgebra
using DiffEqOperators, OrdinaryDiffEq

@testset "Matrix Free Operator" begin
  f = (du, u, p) -> begin
    cumsum!(du, u)
    @. du = p*du
  end
  p = 1.
  A = MatrixFreeOperator(f, (p,), size=(5,5), opnorm=(p)->5)
  b = rand(5)
  @test is_constant(A)
  prob = ODEProblem(A, b, (0,1.), p)
  @test_nowarn solve(prob, Tsit5())

  f = (du, u, p, t) -> cumsum!(du, u)
  args = (1, 1)
  A = MatrixFreeOperator(f,args, size=(5,5), opnorm=(p)->5)
  @test is_constant(A) == false
  b = rand(5)
  @test_nowarn solve(prob, LawsonEuler(krylov=true, m=5), dt=0.1)

  A1 = [[4 -6]; [1 -1]]
  A2 = [[1 2]; [3 2]]
  THRESHOLD = 1
  # analytic solution for the DE
  u1(t) = t < THRESHOLD ? exp(t)*(-1+3exp(t)) : exp(-3-t)*(exp(5)*(-2+7exp(1))+exp(5t)*(-3+8exp(1)))/5
  u2(t) = t < THRESHOLD ? -exp(t)/2+exp(2t)   : (2-7exp(1))*exp(2-t)/5+exp(-3+4t)*(-3+8exp(1))*(3/10)
  sol_analytic(t) = [u1(t), u2(t)]
  # setups for DE solvers
  function Q!(df, f, p, t)
    i = t >= THRESHOLD ? 2 : 1
    A = p[i]
    mul!(df, A, f)
  end
  p = (A1,A2)
  O = MatrixFreeOperator(Q!, (p, 0.), size=(2,2), opnorm=(p)->10)
  # solve DE numerically
  T = 2
  f_0 = [2.0; 1/2]
  prob = ODEProblem(O,f_0,(0.0,T),p)
  sol = solve(prob, LinearExponential(krylov=:simple, m=3), tstops=THRESHOLD)
  @test sol(THRESHOLD) ≈ sol_analytic(THRESHOLD) atol = 1e-3
  @test sol(T) ≈ sol_analytic(T) atol = 1e-3
end


text:
module Sheets

export RWGSheet, read_sheet_data, write_sheet_data, find_unique_periods
export rotate!, translate!, combine, recttri, SV2, MV2

using StaticArrays: SVector, MVector, SMatrix
using ..PSSFSSLen
using JLD2
using LinearAlgebra: norm
using RecipesBase

const MV2 = MVector{2,Float64}
const SV2 = SVector{2,Float64}

abstract type Sheet end

mutable struct RWGSheet <: Sheet
  style::String
  units::PSSFSSLength # Length unit
  s₁::SV2 # Direct lattice vector (specified units)
  s₂::SV2 # Direct lattice vector (specified units)
  β₁::SV2 # Reciprocal lattice vector (1/(specified units))
  β₂::SV2 # Reciprocal lattice vector (1/(specified units))
  dx::Float64 # Unit cell displacment in x (in specified units)
  dy::Float64 # Unit cell displacment in y (in specified units)
  rot::Float64 # Rotation angle for unit cell (deg)
  ρ::Vector{SV2} # Node coordinates
  e1::Vector{Int} # Edge connect. list. e1[i] is the initial node of edge i
  e2::Vector{Int} # Edge connect. list. e2[i] is the terminal node of edge i
  fv::Array{Int,2} #  Face/vertex list. fv[:,i] lists vertices of face i
  fe::Array{Int,2} #  Face/edge list. fe[:,i] lists edges of face i
  fr::Vector{Float64} # Face resistance list. fr[i] is the sheet resistance of face i (Ω/□)
  # The following fields are storage for face/face integrals:
  J::Vector{ComplexF64}
  J_ξ::Vector{ComplexF64}
  J_η::Vector{ComplexF64}
  K::Vector{ComplexF64}
  K_ξ::Vector{ComplexF64}
  K_η::Vector{ComplexF64}
  ρ_r::Vector{SV2}
  rinv::Vector{Float64}
  # Parameters that the face/face integrals depend on:
  ψ₁::Float64  # Incremental phase shift (radians)
  ψ₂::Float64  # Incremental phase shift (radians)
  u::Float64   # Smoothing parameter (1/(specified units))

  class::Char # Type of sheet. 'J' for electric current, 'M' for magnetic current
  info::String  # Informational comment
  # The following flag tells rwg_setup whether (.true.) or not (.false.)
  # to check for consistent edges at xi or eta = 0 and 1.  The default
  # value (.TRUE.) means that the check should be performed.
  ξη_check::Bool
  # The following flag tells rwg_setup whether (.true.) or not (.false.)
  # to Find Unique Face Pairs.
  fufp::Bool
end # struct
import Base.==
==(sh1::RWGSheet, sh2::RWGSheet) = all((getfield(sh1,f)==getfield(sh2,f) for
                                                    f in fieldnames(RWGSheet)))

# Add a zero-argument constructor:
RWGSheet() = RWGSheet("", u"mm",            # style, units
                      SV2([0.0,0.0]),       # s₁
                      SV2([0.0,0.0]),       # s₂
                      SV2([0.0,0.0]),       # β₁
                      SV2([0.0,0.0]),       # β₂
                      0.0, 0.0, 0.0,        # dx, dy, rot
                      SV2[],                # ρ
                      Int[], Int[],         # e1, e2
                      Array{Int}(undef,0,0), # fv
                      Array{Int}(undef,0,0), # fe
                      Float64[],            # fr
                      ComplexF64[],         # J
                      ComplexF64[],         # J_ξ
                      ComplexF64[],         # J_η
                      ComplexF64[],         # K
                      ComplexF64[],         # K_ξ
                      ComplexF64[],         # K_η
                      Array{SV2}(undef,0),  # ρ_r
                      Float64[],            # rinv
                      0.0, 0.0, 0.0,        # ψ₁, ψ₂, u
                      ' ', "",              # class, info
                      true, false)          # ξη_check, fufp

Base.show(io::IO, ::MIME"text/plain", s::RWGSheet) =
    print(io, "RWGSheet: style=", s.style, ", class=", s.class, ", ", length(s.ρ), " nodes, ", length(s.e1),
            " edges, ", size(s.fv,2), " faces")

"""
    read_sheet_data(filename::AbstractString)::RWGSheet

Read the sheet geometry data from a `JLD2` file named in `filename`.
"""
function read_sheet_data(filename::AbstractString)::RWGSheet
    jldopen(filename, "r") do file
        try
            return file["sheet"]
        catch
            @error "$(filename) does not contain sheet data"
        end
    end
end # function

"""
    write_sheet_data(filename::AbstractString, sheet::RWGSheet)

Write the sheet geometry data to a `JLD2` file named in `filename`.
"""
function write_sheet_data(filename::AbstractString, sheet::RWGSheet)
    jldopen(filename, "w") do file
        file["sheet"] = sheet
    end
end

"""
    find_unique_periods(junction::Vector{Int}, sheets)

Find the unique unit cells for the sheets used in the FSS analysis.

# Arguments
- `junction`:  An integer array of length `(Nlayer-1)` containing
               in location `i` the index of the FSS sheet located
               at the interface of dielectric layers `i` and `i+1`. If
               no sheet is present there, the value is 0.
- `sheets`:     An iterable that contains the FSS sheets.
# Return Value
- `upa`            (Unique Periodicity Array) An integer array of the
                   same length of junction, containing zeros in the same
                   locations.  The nonzero entries correspond to sheet
                   locations, and are numbered consecutively according
                   to the equivalence class of the sheet at that location.
                   Two sheets are equivalent if they have the same unit cell.
"""
function find_unique_periods(junction::Vector{Int}, sheets)
    all(t isa Sheet for t in sheets) || error("Elements of sheets must be of type Sheet")
    one_meter = map(x -> ustrip(Float64, x.units, 1.0u"m"), sheets)
    s1s2 = vcat(map(x -> hcat(x.s₁..., x.s₂...), sheets)...) # Each row is s1x s1y s2x s2y
    s1s2 = s1s2 ./ one_meter # All rows now are comparable (in meters)
    s1s2 = round.(s1s2, sigdigits=8)

    upa = zeros(Int, length(junction))
    Nup = 0  # Initialize Number of Unique Periodicities.
    for i in 1:length(upa)  # Step through each junction.
        isht = junction[i] # Sheet index.
        ((isht == 0) || (sheets[isht].style == "NULL")) && continue
        # Compare s1 and s2 of current (isht) sheet with previous sheets.
        for n in 1:Nup  # compare to one member of each equivalence class.
            # Find a sheet that is in equivalence class n:
            nsht = junction[findfirst(isequal(n), upa)]  # Index of sheet to be compared.
            # Compare unit cell of sheet nsht with that of sheet isht:
            if view(s1s2, isht, :) == view(s1s2, nsht, :)
                # Sheets are in the same equivalence class
                upa[i] = n  # Store equiv. class number.
                @goto NextOuterFor
            end
        end
        # If execution fell through to here, we found a sheet that is not
        # in an existing equivalence class.
        Nup += 1 # Bump count of equiv. classes.
        upa[i] = Nup  # Store equiv. class number.
        @label NextOuterFor
    end
    return upa
end

"""
    rotate!(sh::RWGSheet, rot::Real)

Rotate a sheet by rot degrees (counter-clockwise).
"""
function rotate!(sh::RWGSheet, rot::Real)
    rot == 0 && return
    s,c = sincosd(rot)
    rotmat = SMatrix{2,2}([c -s;s c])
    sh.s₁ = rotmat * sh.s₁
    sh.s₂ = rotmat * sh.s₂
    sh.β₁ = rotmat * sh.β₁
    sh.β₂ = rotmat * sh.β₂
    for n in eachindex(sh.ρ)
        sh.ρ[n] = rotmat * sh.ρ[n]
    end
    sh.rot = rot
    return sh
end

"""
    translate!(sh::RWGSheet, dx, dy)

Translate a sheet by dx in x and dy in y.
"""
function translate!(sh::RWGSheet, dx::Real, dy::Real)
    dx == dy == 0 && (return sh)
    tvec = [dx,dy]
    for n in eachindex(sh.ρ)
        sh.ρ[n] = tvec + sh.ρ[n]
    end
    return sh
end

"""
    combine(sh1::RWGSheet, sh2::RWGSheet, dup_coor::Char, dup_coor_value::Real)

Combine the triangulations stored in sheets `sh1` and `sh2`.

# Arguments:
- `sh1`, `sh2`:  sheets having
              initialized values for fields `units`, `ρ`,
              `e1`, `e2` `fv`, `fe`, and possibly `fr`.  It is assumed
              that the two triangulations do not overlap except possibly
              along a line defined by `dup_coor` and `dup_coor_value`,
              as discussed below.  If they do coincide along such
              a line, then they must share the same set of vertices
              and edges along this line.  These duplicate vertices
              and edges will be removed by this routine.
- `dup_coor`    Either 'x' or 'y' to indicate at which coordinate
              constant line the two triangulations may overlap,
              requiring redundant edges and nodes to be removed,
              or ' ' indicating that no search for duplicate nodes
              is required.
- `dup_coor_value`  The value of the coordinate at which the two
              input triangulations overlap.
# Return value
- `sh3`         A `RWGSheet` instance with the following member arrays
                initialized: `units`, `ρ`, `ec`, `fv`, `fe`, and `fr`.
"""
function combine(sh1::RWGSheet, sh2::RWGSheet, dup_coor::Char, dup_coor_value::Real)
    sh1.units == sh2.units || error("Inconsistent units for sh1 and sh2")
    # Count number of vertices located at the duplicate coordinate.
    # Save vertex indices of matching points in vcen1 and vcen2
    #
    Nvcen = 0
    Necen = 0
    vcen1 = Int[]
    vcen2 = Int[]
    if dup_coor ≠ ' '
        tol = 0.5e-4 * norm(sh1.ρ[sh1.e1[1]] - sh1.ρ[sh1.e2[1]])
        for i in 1:length(sh2.ρ)
            if dup_coor == 'x'
                test = sh2.ρ[i][1]
            elseif dup_coor == 'y'
                test = sh2.ρ[i][2]
            else
                test = typemax(typeof(test))
            end
            if abs(test - dup_coor_value) < tol
                # Test to see if there is a sh1 vertex at same coordinate:
                n1match = 0
                for n1 in 1:length(sh1.ρ)
                    if norm(sh1.ρ[n1] - sh2.ρ[i]) < tol
                        n1match = n1
                        break
                    end
                end
                if n1match ≠ 0
                    Nvcen += 1
                    push!(vcen2, i)
                    push!(vcen1, n1match)
                end
            end
        end
        Necen = Nvcen - 1 # Number of shared edges.
        ecen1 = zeros(Int, Necen); ecen2 = zeros(Int, Necen)
        # Locate and save indices in sh1 and sh2 of edges along the center line:
        Necen = 0
        for i in 1:length(sh2.e1) # Loop over sh2 edges
            if (sh2.e1[i] in vcen2) && (sh2.e2[i] in vcen2)
                # Edge i of sh2 is a shared edge.
                Necen += 1
                ecen2[Necen] = i
                i1 = findfirst(x->x==sh2.e1[i], vcen2)
                i2 = findfirst(x->x==sh2.e2[i], vcen2)
                # Now find matching sh1 edge:
                for j1 in 1:length(sh1.e1)
                    if  ((sh1.e1[j1] == vcen1[i1]) && (sh1.e2[j1] == vcen1[i2])) ||
                        ((sh1.e2[j1] == vcen1[i1]) && (sh1.e1[j1] == vcen1[i2]))
                        ecen1[Necen] = j1
                        @goto sh2edges
                    end
                end
            else
                continue
            end
            error("Unable to find matching duplicate edge in combine_sheet")
            @label sh2edges
        end
    end

    # Allocate triangulation arrays in new sheet:
    sh3 = RWGSheet()
    sh3.e1 = zeros(Int, length(sh1.e1) + length(sh2.e1) - Necen)
    sh3.e2 = zeros(Int, length(sh1.e2) + length(sh2.e2) - Necen)
    sh3.ρ = Vector{SV2}(undef, length(sh1.ρ) + length(sh2.ρ) - Nvcen)
    sh3.fe = zeros(Int, 3, size(sh1.fe,2) + size(sh2.fe,2))
    sh3.fv = zeros(Int, 3, size(sh1.fv,2) + size(sh2.fv,2))

    # Copy vertex locations:
    sh3.ρ[1:length(sh1.ρ)] = sh1.ρ
    if dup_coor == ' '
        sh3.ρ[(1+length(sh1.ρ)):end] = sh2.ρ
    else
        i2 = 1 + length(sh1.ρ)  # Node counter
        for i in 1:length(sh2.ρ)
            if !(i in vcen2[1:Nvcen])
                sh3.ρ[i2] = sh2.ρ[i]
                i2 += 1
            end
        end
    end
    # Copy edge-node matrices
    eoffset = length(sh1.e1)
    voffset = length(sh1.ρ)
    sh3.e1[1:length(sh1.e1)] = sh1.e1
    sh3.e2[1:length(sh1.e2)] = sh1.e2
    if dup_coor == ' '
        sh3.e1[(eoffset+1):end] = sh2.e1 .+ voffset
        sh3.e2[(eoffset+1):end] = sh2.e2 .+ voffset
    else
        for i in 1:length(sh2.e1)
            if i in ecen2
                # Edge i of sh2 is located on duplication line.
                # Do not include this duplicate edge in sh3, but
                # decrement the edge index offset:
                eoffset -= 1
            else
                if sh2.e1[i] in vcen2
                    # The initial point of edge i of sh2 is a duplicate vertex.
                    i2 = findfirst(x->x==sh2.e1[i], vcen2)
                    sh3.e1[i+eoffset] = vcen1[i2]
                else
                    # Ordinary point
                    sh3.e1[i + eoffset] = sh2.e1[i] + voffset
                end
                #
                if sh2.e2[i] in vcen2
                    # The terminal point of edge i of sh2 is on the duplication edge.
                    i2 = findfirst(x->x==sh2.e2[i], vcen2)
                    sh3.e2[i+eoffset] = vcen1[i2]
                else
                    # Ordinary point:
                    sh3.e2[i+eoffset] = sh2.e2[i] + voffset
                end
            end
        end
        # Correct vertex indices:
        for i in Nvcen:-1:1
            sh3.e1[sh3.e1 .> (length(sh1.ρ) + vcen2[i])] .-= 1
            sh3.e2[sh3.e2 .> (length(sh1.ρ) + vcen2[i])] .-= 1
        end
    end
    # Copy face/vertex matrix
    sh3.fv[:, 1:size(sh1.fv,2)] = sh1.fv
    sh3.fv[:, 1+size(sh1.fv,2):end] = sh2.fv .+ voffset # offset will be corrected later
    #  Correct duplicate vertices from sh2:
    for n2 in 1:Nvcen # Examine each duplicate vertex in sh2
        n1 = vcen1[n2] # Initialize index of matching vertex in sh1
        # Replace all references to vertex vcen2[n2] with ref to vcen1[n1]
        sh3.fv[sh3.fv .== vcen2[n2] + voffset] .= n1
    end
    # Correct vertex indices:
    for i in Nvcen:-1:1
        sh3.fv[sh3.fv .> voffset + vcen2[i]] .-= 1
    end
    # Copy face/edge matrix:
    foffset = size(sh1.fv,2)
    eoffset = length(sh1.e1)
    sh3.fe[:,1:size(sh1.fe,2)] = sh1.fe
    sh3.fe[:,1+size(sh1.fe,2):end] = sh2.fe .+ eoffset # offset will be corrected later
    #  Correct duplicate edges from sh2:
    if dup_coor ≠ ' '
        for n2 in 1:length(ecen2) # Examine each duplicate edge in sh2
            e2 = ecen2[n2]  # Index of duplicate edge in sh2.
            e1 = ecen1[n2]  # Index of duplicate edge in sh1.
            sh3.fe[sh3.fe .== (e2 + eoffset)] .= e1
        end
        # Correct edge indices:
        for i in length(ecen2):-1:1
            sh3.fe[sh3.fe .> eoffset + ecen2[i]] .-= 1
        end
    end
    return sh3
end

"""
    recttri(rhobl::SVector{2,Float64}, rhotr::SVector{2,Float64}, nx::Int, ny::Int)

Create a variable of type `RWGSheet` that contains the triangulation for
a rectangular strip.  The fields `ρ`, `e1`, `e2`, `fv`, and `fe` properly initialized.
"""
function recttri(rhobl::SV2, rhotr::SV2, nx::Int, ny::Int)
    nodecount = (nx+1) * (ny+1)  # Number of nodes.
    edgecount = 3*nx*ny + nx + ny  # Number of edges.
    facecount = 2*nx*ny  # Number of faces.

    sh = RWGSheet()
    sh.ρ = Vector{SV2}(undef, nodecount)

    # Set the node coordinates:
    drho = (rhotr - rhobl) ./ [nx, ny]
    n = 0  # Initialize node index.
    for j in 0:ny
        yj = j * drho[2]
        for i in 0:nx
            n += 1
            sh.ρ[n] = rhobl + SV2([i*drho[1], yj])
        end
    end

    sh.e1 = zeros(Int, edgecount)
    sh.e2 = zeros(Int, edgecount)
    e = 0  # Initialize edge index.
    # Do the horizontal edges:
    for j in 0:ny
        kadd = j * (nx+1)
        for i in 1:nx
            e += 1
            sh.e1[e] = i + kadd
            sh.e2[e] = sh.e1[e] + 1
        end
    end
    # Do the vertical edges:
    for j in 1:ny
        kadd = (j-1) * (nx+1) + 1
        for i in 0:nx
            e += 1
            sh.e1[e] = i + kadd
            sh.e2[e] = sh.e1[e] + (nx+1)
        end
    end
    # Do the diagonal edges:
    for j in 1:ny
        kadd1 = (j-1) * (nx+1)
        kadd2 = 1 + j * (nx+1)
        for i in 1:nx
            e += 1
            sh.e1[e] = i + kadd1
            sh.e2[e] = i + kadd2
        end
    end

    # Done with edges.  Begin setting up faces.
    # Allocate arrays whose length depends only on the number of faces:
    sh.fv = zeros(Int, 3, facecount)
    sh.fe = zeros(Int, 3, facecount)
    # Set up the face/vertex and face/edge matrices:
    nhe = nx*ny + nx  # Number of horizontal edges.
    nve = nx*ny + ny  # Number of vertical edges.
    nde = nx*ny       # Number of diagonal edges
    f = 0  # Initialize face index.
    for j in 1:ny
        nadd1 = (j-1) * (nx+1)
        nadd2 = 1 + j * (nx+1)
        for i in 1:nx
            f += 1  # Bump face index (upper left face).
            sh.fv[1,f] = i + nadd1  # Lower Left vertex.
            sh.fv[2,f] = i + nadd2  # Upper right vertex.
            sh.fv[3,f] = i + nadd2 - 1 # Upper left vertex.
            sh.fe[1,f] = i + j*nx # Upper edge.
            sh.fe[2,f] = i + nhe + nadd1 # Left edge
            sh.fe[3,f] = i + nhe + nve + (j-1)*nx # Diagonal edge
            f += 1  # Bump face index (lower right face).
            sh.fv[1,f] = sh.fv[1,f-1] # Lower Left vertex.
            sh.fv[2,f] = 1 + sh.fv[1,f-1] # Lower right vertex.
            sh.fv[3,f] = sh.fv[2,f-1] # Upper right vertex.
            sh.fe[1,f] = 1 + sh.fe[2,f-1] # Right edge.
            sh.fe[2,f] = sh.fe[3,f-1] # Diagonal edge.
            sh.fe[3,f] = sh.fe[1,f-1] - nx # Bottom edge
        end
    end

    return sh
end

"Plot recipe for RWGSheet"
@recipe function f(sh::RWGSheet; edges=true, faces=false, nodes=false,
                    edgenumbers=false, facenumbers=false, nodenumbers=false,
                    unitcell=false, rep=(1,1), fontsize=9)
    # set a default value for an attribute with `-->`.  Force it with `:=`.
    xguide --> "x ($(sh.units))"
    yguide --> "y ($(sh.units))"
    aspect_ratio := :equal

    if isa(rep[1], Int)
        mrange = 1:rep[1]
    elseif isa(rep[1], UnitRange)
        mrange = rep[1]
    else
        error("Illegal type for rep[1]")
    end
    if isa(rep[2], Int)
        nrange = 1:rep[2]
    elseif isa(rep[2], UnitRange)
        nrange = rep[2]
    else
        error("Illegal type for rep[2]")
    end

    for m in mrange, n in nrange
        x0, y0 = (m-1)*sh.s₁ + (n-1)*sh.s₂

        # Add series for faces
        if faces
            for i in 1:size(sh.fv,2)
                points = sh.ρ[sh.fv[:,i]]
                x = x0 .+ [point[1] for point in points]
                y = y0 .+ [point[2] for point in points]
                @series begin
                    seriestype := :shape
                    # ignore series in legend and color cycling
                    primary := false
                    linecolor := nothing
                    fillcolor --> :blue
                    fillalpha --> 0.8
                    markershape := :none
                    x, y
                end
            end
        end

        # Add series for edges
        if edges
            x = Float64[]
            y = Float64[]
            for i in 1:length(sh.e1)
                points = sh.ρ[[sh.e1[i], sh.e2[i]]]
                push!(x,NaN)
                push!(y,NaN)
                append!(x, [p[1] for p in points])
                append!(y, [p[2] for p in points])
            end
            x .+= x0
            y .+= y0
            @series begin
                seriestype := :path
                # ignore series in legend and color cycling
                primary := false
                linecolor --> :black
                linestyle := :solid
                fillcolor := nothing
                fillalpha := 0
                markershape := :none
                x, y
            end
        end

        # Add series for unit cell
        if unitcell
            points = [0*sh.s₁, sh.s₁, sh.s₁+sh.s₂, sh.s₂]
            x = x0 .+ [point[1] for point in points]; push!(x, x[1])
            y = y0 .+ [point[2] for point in points]; push!(y, y[1])
            @series begin
                seriestype := :path
                # ignore series in legend and color cycling
                primary := false
                linecolor := :blue
                linestyle := :dot
                fillcolor := nothing
                fillalpha := 0
                markershape := :none
                x, y
            end
        end

        # Add series for nodes
        if nodes
            x = x0 .+ [p[1] for p in sh.ρ]
            y = y0 .+ [p[2] for p in sh.ρ]
            @series begin
                seriestype := :scatter
                # ignore series in legend and color cycling
                primary := false
                linecolor := nothing
                markercolor --> :black
                markershape --> :circle
                markersize -->  1
                x, y
            end
        end

        # Add series for node numbers
        if nodenumbers
            x = x0 .+ [p[1] for p in sh.ρ]
            y = y0 .+ [p[2] for p in sh.ρ]
            @series begin
                seriestype := :scatter
                # ignore series in legend and color cycling
                primary := false
                linecolor := nothing
                markersize := 0
                markeralpha := 0
                markercolor := nothing
                markershape := :none
                annotations := [(x[i], y[i], string(i), fontsize) for i in 1:length(x)]
                x,y
            end
        end

        # Add series for edge numbers
        if edgenumbers
            x = zeros(Float64, length(sh.e1))
            y = zeros(Float64, length(sh.e1))
            for i in 1:length(sh.e1)
                x[i], y[i] = 0.5*sum(sh.ρ[[sh.e1[i], sh.e2[i]]])
            end
            x .+= x0
            y .+= y0
            @series begin
                seriestype := :scatter
                # ignore series in legend and color cycling
                primary := false
                linecolor := nothing
                markersize := 0
                markeralpha := 0
                markercolor := nothing
                markershape := :none
                annotations := [(x[i], y[i], string(i), fontsize) for i in 1:length(x)]
                x,y
            end
        end

        # Add series for face numbers
        if facenumbers
            x = zeros(Float64, size(sh.fv,2))
            y = zeros(Float64, size(sh.fv,2))
            for i in 1:size(sh.fv,2)
                x[i], y[i] = (1/3)*sum(sh.ρ[sh.fv[:,i]])
            end
            x .+= x0
            y .+= y0
            @series begin
                seriestype := :scatter
                # ignore series in legend and color cycling
                primary := false
                linecolor := nothing
                markersize := 0
                markeralpha := 0
                markercolor := nothing
                markershape := :none
                annotations := [(x[i], y[i], string(i), fontsize) for i in 1:length(x)]
                x,y
            end
        end
    end
end

end # module


text:
println("Hello world")

# something.
#new change on the main branch

#something again.

#just this one other thing. something added to the temp branch.

max_stars_repo_name: RodolfoBreves/testGit

max_issues_repo_name: RodolfoBreves/testGit

max_forks_repo_name: RodolfoBreves/testGit