text:
"""
Loads in PowerModels network data given the name of a network case file,
then starts creating samples
"""
function create_samples(net::String, K=Inf; U=0.0, S=0.0, V=0.0, max_iter=Inf, T=Inf, discard=false, variance=false,
						input_vars=DEFAULT_INPUTS, output_vars=DEFAULT_OUTPUTS, dual_vars=DEFAULT_DUALS,
						sampler=sample_polytope_cprnd, sampler_opts=Dict{Symbol,Any}()::Dict{Symbol}, A=[]::Array, b=[]::Array,
						pd_max=nothing, pd_min=nothing, pf_min=0.7071, pf_lagging=true, save_certs=false, save_max_load=false,
						print_level=0, stat_track=false, save_while=false, save_infeasible=false, save_path="", net_path="",
						model_type=PM.QCLSPowerModel, r_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL),
						opf_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	net, net_path = load_net(net, net_path, print_level)

	return create_samples(net, K; U=U, S=S, V=V, max_iter=max_iter, T=T, discard=discard, variance=variance,
							  input_vars=input_vars, output_vars=output_vars, dual_vars=dual_vars, A=A, b=b,
							  sampler=sampler, sampler_opts=sampler_opts, save_max_load=save_max_load,
							  pd_max=pd_max, pd_min=pd_min, pf_min=pf_min, pf_lagging=pf_lagging, save_certs=save_certs,
							  print_level=print_level, stat_track=stat_track, save_while=save_while,
							  save_infeasible=save_infeasible, save_path=save_path, net_path=net_path,
							  model_type=model_type, r_solver=r_solver, opf_solver=opf_solver)
end

"""
    create_samples(net, K; <keyword arguments>)

Creates an AC OPF dataset for the given PowerModels network dictionary. Generates samples until one of the given stopping criteria is met.
Takes options to determine how to sample points, what information to save, and what information is printed.

# Examples
```
julia> results = create_samples("case5.m", 100; T=1000, net_path="data")
```

# Arguments
- 'net::Dict': network information stored in a PowerModels.jl format specified dictionary
- 'K::Integer': the maximum number of samples before stopping sampling
- 'U::Float': the minimum % of unique active sets sampled in the previous 1 / U samples to continue sampling
- 'S::Float': the minimum % of saved samples in the previous 1 / L samples to continue sampling
- 'V::Float': the minimum % of feasible samples that increase the variance of the dataset in the previous 1 / L samples to continue sampling
- 'T::Integer': the maximum time for the sampler to run in seconds.
- 'max_iter::Integer': maximum number of iterations for the sampler to run for.
- 'sampler::Function': the sampling function to use. This function must take arguements A and b, and can take optional arguments.
- 'sampler_opts::Dict': a dictionary of optional arguments to pass to the sampler function.
- 'A::Array': defines the initial sampling space polytope Ax<=b. If not provided, initializes to a default.
- 'b::Array': defines the initial sampling space polytope Ax<=b. If not provided, initializes to a default.
- 'pd_max::Array': the maximum active load values to use when initializing the sampling space and constraining the loads. If nothing, finds the maximum load at each bus with the given relaxed model type.
- 'pd_min::Array': the minimum active load values to use when initializing the sampling space and constraining the loads. If nothing, this is set to 0 for all loads.
- 'pf_min::Array/Float:' the minimum power factor for all loads in the system (Number) or an array of minimum power factors for each load in the system.
- 'pf_lagging::Bool': indicating if load power factors can be only lagging (True), or both lagging or leading (False).
- 'reset_level::Integer': determines how to reset the load point to be inside the polytope before sampling. 2: Reset closer to nominal load & chebyshev center, 1: Reset closer to chebyshev center, 0: Reset at chebyshev center.
- 'save_certs::Bool': specifies whether the sampling space, Ax<=b (A & b matrices) are saved to the results dictionary.
- 'save\\_max_load::Bool': specifies whether the max active load demands used are saved to the results dictionary.
- 'model_type::Type': an abstract PowerModels type indicating the network model to use for the relaxed AC-OPF formulations (Max Load & Nearest Feasible)
- 'r_solver': an optimizer constructor used for solving the relaxed AC-OPF optimization problems.
- 'opf_solver': an optimizer constructor used to find the AC-OPF optimal solution for each sample.
- 'print_level::Integer': from 0 to 3 indicating the level of info to print to console, with 0 indicating minimum printing.
- 'stat_track::Integer': from 0 to 3 indicating the level of stats info saved during each iteration	0: No information saved, 1: Feasibility, New Certificate, Added Sample, Iteration Time, 2: Variance for all input & output variables
- 'save_while::Bool': indicates whether results and stats information is saved to a csv file during processing.
- 'save_infeasible::Bool': indicates if infeasible samples are saved. If true saves infeasible samples in a seperate file from feasible samples.
- 'save_path::String:' a string with the file path to the desired result save location.
- 'net_path::String': a string with the file path to the network file.
- 'variance::Bool': indicates if dataset variance information is tracked for each unique active set.
- 'discard::Bool': indicates if samples that do not increase the variance within a unique active set are discarded.

See 'OPF-Learn: An Open-Source Framework for Creating Representative AC Optimal Power Flow Datasets'
for more information on how the AC OPF datasets are created.

Modified from AgenerateACOPFsamples.m written by Ahmed S. Zamzam
"""
function create_samples(net::Dict, K=Inf; U=0.0, S=0.0, V=0.0, max_iter=Inf, T=Inf, discard=false, variance=false,
							input_vars=DEFAULT_INPUTS, output_vars=DEFAULT_OUTPUTS, dual_vars=DEFAULT_DUALS,
							sampler=sample_polytope_cprnd, sampler_opts=Dict{Symbol,Any}()::Dict{Symbol}, A=[]::Array, b=[]::Array,
							pd_max=nothing, pd_min=nothing, pf_min=0.7071, pf_lagging=true, reset_level=0, save_certs=false, save_max_load=false,
							print_level=0, stat_track=false, save_while=false, save_infeasible=false, save_path="", net_path="",
							model_type=PM.QCLSPowerModel, r_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL),
							opf_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	now_str = Dates.format(Dates.now(), "dd-mm-yyy_HH.MM.SS")  # Get date & time for result file names
	net_name = net["name"]
	save_order = vcat(input_vars, output_vars, dual_vars)
    if net_path == ""  # Set net path to save path if not specified for saving network max loads
		net_path = save_path
	end

	# Gather network information used during processing
	A, b, x, results, fnfp_model, base_load_feasible, net_r = initialize(net, pf_min, pf_lagging, pd_max, pd_min,
																		 input_vars, output_vars, dual_vars,
																		 save_certs, save_infeasible, save_while,
																		 stat_track, save_max_load, A, b,
																		 sampler, sampler_opts, net_name,
																		 net_path, model_type, r_solver,
																		 reset_level, print_level)

	AC_inputs = results["inputs"]
	AC_outputs = results["outputs"]
	duals = results["duals"]

	save_infeasible && (infeasible_AC_inputs = results["infeasible_inputs"])
	save_certs && (store_cert(results, A, b))
	stat_track > 0 && (stats = results["stats"])

	# Set up sample processing while loop
	k = 0  # Count of feasible samples collected
    i = 1  # Count of iterations
	s = 0  # Count of samples since last saved (not discarded) sample
	u = 0  # Count of samples since last unique active set found
	v = 0  # Count of samples since last increase in variance seen
	start_time = time()  # Start time in seconds
	m = size(A, 1)  # Number of polytope planes
    while (k < K) & (u < (1 / U)) & (s < (1 / S)) & (v < 1 / V) &
		  (i < max_iter) & ((time() - start_time) < T)
        iter_start_time = time()

		if print_level > 0
			println("Samples: $(k) / $(K),\t Iter: $(i)")
		elseif i % (max_iter / 10) == 0
			println("Iter: $(i)")
		end

		# Set stats to defaults
		iter_stats = Dict(
						  :new_cert	=> false,
						  :new_set 	=> false,
						  :var_inc 	=> false,
						  :added 	=> false,
						  :feasible	=> false,
						  :iter_time => -1.0,
						  :active_set => [],
						  )

        # Generate sample of load profile
        m_ = size(A, 1)
        if m < m_  # New infeasibility certificate was added to the polytope
			center, radius = chebyshev_center(A, b)
			if reset_level > 1
				x = (base_load_feasible * 0.1 + x * 0.9) * 0.9 + center' * 0.1
			elseif reset_level > 0
				x = x * 0.9 + center' * 0.1
			else
				x = center'
			end
			m = m_
		end

        # Sample uniformly from the interior of the convex polytope space
		x = sampler(A, b, x, 1; sampler_opts...)

        # Set network loads to sampled values
        set_network_load(net, x, scale_load=false)

        # Solve OPF for the load sample
		result, feasible = run_ac_opf(net, print_level=print_level, solver=opf_solver)
		print_level > 0 && println("OPF SUCCESS: " * string(feasible))

        if feasible
			s, u, v, k, iter_stats = store_feasible_sample(s, u, v, k, i, K, iter_stats,
										               AC_inputs, AC_outputs, duals, dual_vars,
													   x, result, discard, variance, net_name,
													   now_str, save_path, save_while, save_order,
													   print_level)
			print_level > 0 && println("Samples: $(k) / $(K),\t Iter: $(i)")
        else
			save_infeasible && store_infeasible_sample(infeasible_AC_inputs, x, result,
							save_while, net_name, now_str, save_order, dual_vars, save_path)
            px = x[1:Integer(length(x)/2)]
            qx = x[Integer(length(x)/2) + 1:end]

            r, pd, qd, solved = find_nearest_feasible(fnfp_model, px, qx, print_level=print_level)
			r = sum((pd .- px).^2 + (qd .- qx).^2)

			print_level > 0 && println("R: $(r)")
            if (r > R_TOLERANCE) & solved
                iter_stats[:new_cert] = true

				# Solve OPF for the random sample:
                # try to solve the non-convex AC OPF at this solution
                # note that: while this solution is feasible for the relaxation
                # it can be infeasible for the original problem.
                xp_ = pd .- INFEASIBILITY_CERT_SHIFT
                xq_ = qd .- INFEASIBILITY_CERT_SHIFT

				xp_[xp_ .< 0] .= 0 # Ensure loads do not go negative if shifted
				xq_[xq_ .< 0] .= 0

                x_star = vcat(xp_, xq_)
                x_x_star = (x - x_star)'

                # Add infeasibility certificate
                A = vcat(A, x_x_star)
                b = vcat(b, x_x_star * x_star)
				save_while && save_cert(A, b, net_name*"_"*now_str, dir=save_path)
				save_certs && store_cert(results, A, b)

				# Slightly shift FNFP to make it more likely to be feasible
                xp_ = pd .- FNFP_SHIFT
                xq_ = qd .- FNFP_SHIFT

				# Removes any negative values introduced by the shift
				xp_[xp_ .< 0] .= 0
				xq_[xq_ .< 0] .= 0

				x_star = vcat(xp_, xq_)
                x = x_star
                set_network_load(net, x, scale_load=false)

                # Solve OPF for the relaxation feasible sample
				result, feasible = run_ac_opf(net, print_level=print_level, solver=opf_solver)
				print_level > 0 && println("FNFP OPF SUCCESS: " * string(feasible))

                if feasible
					s, u, v, k, iter_stats = store_feasible_sample(s, u, v, k, i, K, iter_stats,
												  AC_inputs, AC_outputs, duals, dual_vars,
												  x, result, discard, variance, net_name,
												  now_str, save_path, save_while, save_order,
												  print_level)
					println("Samples: $(k) / $(K),\t Iter: $(i)")
                else
					save_infeasible && store_infeasible_sample(infeasible_AC_inputs, x, result,
								    save_while, net_name, now_str, save_order, dual_vars, save_path)
				end
			end
		end
		iter_elapsed_time = time() - iter_start_time
		iter_stats[:iter_time] = iter_elapsed_time

		if stat_track > 0
			update_stats!(stats, duals, iter_stats, save_level=stat_track)
			save_while && (save_stats(iter_stats, net_name*"_"*now_str*"_stats", dir=save_path))
		end

        i = i + 1
	end

	# Need to convert the unique active sets Set to and Array, due to a bug with PyJulia
	results["duals"]["unique_active_sets"] = collect(results["duals"]["unique_active_sets"])
    return results
end

"""
Gathers the network information needed during the dataset creation process, then initializes
objects used to create and process load samples.
"""
function initialize(net, pf_min, pf_lagging, pd_max, pd_min,
					input_vars, output_vars, dual_vars,
					save_certs, save_infeasible, save_while, stat_track, save_max_load,
					A, b, sampler, sampler_opts,
					net_name, net_path, model_type, r_solver, reset_level, print_level)
	# Get system information for relaxed AC OPF runs
	net_info = gather_net_info(net)
	num_loads = net_info["num_loads"]
	pg_max = net_info["pg_max"]
	base_load = net_info["base_load"]

	# Set pf_min & base_load in relaxed network info dict
	net_r = create_net_r(net, pf_min, pf_lagging, base_load)

	# Construct objects for storing feasible point data
    results = initialize_results(input_vars, output_vars, dual_vars, net_info,
								 save_certs, save_infeasible, stat_track)

    # Find the maximum servable load for each load bus in the system
	if !isnothing(pd_min)
		net_r["pd_min"] = pd_min
		if save_max_load
			!haskey(results, "load_constraints") && (results["load_constraints"]=Dict())
			results["load_constraints"]["pd_min"] = pd_min
		end
	end
	pd_max = find_max_loads(pd_max, net_r, net_path, net_name, model_type, r_solver,
							save_max_load, results, save_while, print_level)

    # Initialize polytope, Ax <= b, parameter data structures
	A, b = initialize_polytope(num_loads, pd_max, pd_min, pf_min, pf_lagging, pg_max,
							   A, b, sampler, sampler_opts)

	# Create find nearest feasible point model without the objective
	pm = PM.instantiate_model(net_r, model_type, build_opf_var_load)
	fnfp_model = find_nearest_feasible_model(pm, print_level=print_level, solver=r_solver)

	# Initialize sampler origin
	x, base_load_feasible = initialize_sample_origin(A, b, base_load, num_loads, fnfp_model,
													 reset_level, print_level)

	return A, b, x, results, fnfp_model, base_load_feasible, net_r
end

"""
Creates A & b matrices defining a polytope, Ax<=b, if A & b are not given.
If A & b are given, checks their validity.

The initialized polytope has load constraints:
constrains the maximum active demand to the given maximum active demand values,
constrains the minimum active demand to the given minimum active demand values,
constrains the minimum reactive demand to zero,
constrains the reactive demand with the given minimum power factor,
constrains the maximum total active load to the sum of generator active power ratings.
"""
function initialize_polytope(num_loads, pd_max, pd_min, pf_min, pf_lagging, pg_max, A=[], b=[],
							 sampler=sample_polytope_cprnd, sampler_opts=Dict())
	if isempty(A) & isempty(b)
		isnothing(pd_min) && (pd_min=zeros((num_loads, 1)))
		pf_min isa Number && (pf_min = ones(num_loads) .* pf_min)
		d = tan.(acos.(pf_min))

		# Max active demand constaint
		max_pd_constr = hcat(I(num_loads), zeros((num_loads, num_loads)))
		# Min active demand constraint
		min_pd_constr = hcat(-I(num_loads), zeros((num_loads, num_loads)))
		# Max reactive demand constraint
		max_qd_constr = hcat(-I(num_loads).*d, I(num_loads))
		# Min reactive demand constraint
		if pf_lagging
			min_qd_constr = hcat(zeros((num_loads, num_loads)), -I(num_loads))
		else
			min_qd_constr = hcat(I(num_loads).*d, -I(num_loads))
		end
		# Total active demand constraint
		tot_pd_constr = hcat(ones((1, num_loads)), zeros((1, num_loads)))

		A = vcat(max_pd_constr,
				 min_pd_constr,
				 max_qd_constr,
				 min_qd_constr,
				 tot_pd_constr,
				 )
		b = vcat(pd_max,
				 pd_min,
				 zeros((num_loads, 1)),
				 zeros((num_loads, 1)),
				 [sum(pg_max)],
				 )
	else  # Check if given polytope definition matrices are valid
		@assert size(A,1) == size(b,1) "A & b dimensions do not align. Must align with equation, Ax<=b"
		@assert size(A,2) == (2 * num_loads) "A must be N x (2*num_loads)"
		@assert size(b,2) == 1 "b must be N x 1"
		try  # to gather a sample with the given polytope definition
			center, radius = chebyshev_center(A, b)
			x = center'
			x = sampler(A, b, x, 1; sampler_opts...)
		catch
			@assert false "The sampler could not find a sample from the given polytope, Ax<=b. The given space must be enclosed."
		end
	end
	return A, b
end

"Finds a point to initialize the sampler from based on the network, polytope, and origin reset option"
function initialize_sample_origin(A, b, base_load, num_loads, fnfp_model, reset_level, print_level)
	# Initialize load sample, x, to chebyshev center
	center, radius = chebyshev_center(A, b)
    x = center'

	# Ensure that the base load is feasible for the relaxation & within the polytope
	if reset_level > 1
		base_pd = base_load[1:num_loads]
		base_qd = base_load[num_loads+1:end]
		r, pd, qd = find_nearest_feasible(fnfp_model, base_pd, base_qd, print_level=print_level)
		base_load_feasible = vcat(pd, qd)
		x = base_load_feasible * 0.9 + center' * 0.1
	else
		base_load_feasible = center'
	end
	return x, base_load_feasible
end


text:
# MATLAB session variable
global s1 = nothing

"""
Start MATLAB session
"""
function init_MATLAB()
    println("Starting MATLAB session...")
    global s1 = MSession()
    # cwd = joinpath(Pkg.dir("ReconstructionEvaluations"), "src")
    cwd = joinpath(pwd(), "ReconstructionEvaluations", "src")
    println("Changing MATLAB userpath to $cwd")
    eval_string(s1, "userpath('$cwd')")
end

function check_MATLAB()
    if s1 == nothing
        init_MATLAB()
    end
end

ext: jl


text:
using CollectiveSpins
using PyCall
@pyimport matplotlib.pyplot as plt

const Γeff_range = Float64[-1:1:1]
const Ωeff_range = Float64[-1:0.5:1]
const T = Float64[0:0.01:5]
const phi0 = 0.
const theta0 = pi/2.

for Ωeff=Ωeff_range
    for Γeff=Γeff_range
        state0 = CollectiveSpins.meanfield.blochstate(phi0, theta0)
        tout, states_t = CollectiveSpins.meanfield.timeevolution_symmetric(T, state0, Ωeff, Γeff)
        sx = [CollectiveSpins.meanfield.sx(state)[1] for state in states_t]
        sy = [CollectiveSpins.meanfield.sy(state)[1] for state in states_t]
        sz = [CollectiveSpins.meanfield.sz(state)[1] for state in states_t]
        plt.figure(1)
        plt.plot(tout, sx)
        plt.ylabel("\$\\sigma_x\$")
        plt.figure(2)
        plt.plot(tout, sy)
        plt.ylabel("\$\\sigma_y\$")
        plt.figure(3)
        plt.plot(tout, sz)
        plt.ylabel("\$\\sigma_z\$")
    end
    plt.figure(1)
    pycall(plt.gca()["set_color_cycle"], PyAny, nothing)
    plt.figure(2)
    pycall(plt.gca()["set_color_cycle"], PyAny, nothing)
    plt.figure(3)
    pycall(plt.gca()["set_color_cycle"], PyAny, nothing)
end

plt.show()


text:
module QuickActors

export QuickActor, QuickAddr, QuickScheduler, run!, send!, spawn!

using ActorInterfaces.Classic
using ActorInterfaces.Implementation
using ActorInterfaces.Implementation.Tick

struct QuickAddr <: Addr
    id::UInt64
end
# Classic.SendStyle(::Type{QuickAddr}) = Sendable()

struct Envelope{TMsg <: Any}
    to::QuickAddr
    msg::TMsg
end

struct QuickScheduler <: Implementation.Scheduler
    actorcache::Dict{QuickAddr, Any}
    msgs::Array{Envelope}
    QuickScheduler() = new(Dict(), [])
end

struct QuickActor{TBehavior}
    addr::QuickAddr
    behavior::TBehavior
end

function QuickActor(behavior)
    return QuickActor{typeof(behavior)}(gen_addr(), behavior)
end
gen_addr() = QuickAddr(rand(UInt64))

struct ActorContext{TActor}
    scheduler::QuickScheduler
    actor::TActor
end

Classic.self(;ctx::ActorContext) = ctx.actor.addr

function Tick.tick!(sdl::QuickScheduler)::Bool
    isempty(sdl.msgs) && return false
    envelope = popfirst!(sdl.msgs)
    ctx = get!(sdl.actorcache, envelope.to, nothing)
    isnothing(ctx) && error("Actor with address $(msg.to) not scheduled.")
    deliver!(ctx, envelope.msg)
    return true
end

function deliver!(ctx, msg)
    ctx.actor.behavior(msg...; ctx)
end

function send!(sdl::QuickScheduler, target::QuickAddr, msg...)
    push!(sdl.msgs, Envelope(target, msg))
    return nothing
end
# send!(sdl::QuickScheduler, target::Addr, msg) = send!(sdl, target, msg, SendStyle(typeof(msg)))

function spawn!(sdl::QuickScheduler, behavior)
    actor = QuickActor(behavior)
    sdl.actorcache[actor.addr] = ActorContext(sdl, actor)
    return actor.addr
end

function spawn!(sdl::QuickScheduler, behavior, aquintance1, aquintances...)
    stateful_bhv = (msg...; ctx) -> behavior(aquintance1, aquintances..., msg...; ctx)
    return spawn!(sdl, stateful_bhv)
end

function Classic.send(target::QuickAddr, msg...; ctx::ActorContext)
    send!(ctx.scheduler, target, msg...)
end

Classic.spawn(behavior, aquintances...; ctx::ActorContext) = spawn!(ctx.scheduler, behavior, aquintances...)

function Classic.become(newbhv; ctx::ActorContext)
    sdl = ctx.scheduler
    actor = QuickActor(self(;ctx), newbhv)
    sdl.actorcache[self(;ctx)] = ActorContext(sdl, actor)
end

function Classic.become(newbhv, aquintance1, aquintances...; ctx::ActorContext)
    stateful_bhv = (msg...; ctx) -> newbhv(aquintance1, aquintances..., msg...; ctx)
    Classic.become(stateful_bhv; ctx)
end

end # module QuickActors


text:
module Faker
using Dates
using Random

import YAML

data=Dict()
data["en"]=Dict()
data["en"]["faker"]=Dict()

function language(file::String) #"es-MX"
	original= pwd()
	cd(abspath(joinpath(dirname(Base.find_package("Faker")), "../src/locales")))
	newdata = YAML.load(open("$(file).yml"))[file]["faker"]
	for (k, v) in newdata
        if (haskey(data["en"]["faker"], k))
			data["en"]["faker"][k]=merge(data["en"]["faker"][k],newdata[k])
		end
	end
	cd(original)
end

function loadfiles(nivel::String)
	original= pwd()
	cd(abspath(joinpath(dirname(Base.find_package("Faker")), "../src/locales/en")))
	if (nivel=="base")
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("address.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("business.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("color.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("company.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("file.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("internet.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("lorem.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("name.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("nation.yml"))["en"]["faker"])
		data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open("phone_number.yml"))["en"]["faker"])
	else
		ls= readdir()
		for i=1:length(ls)
			if ( (occursin(Regex("((.)*\\.(yml))"),ls[i])) && (!occursin(r"^(\.)",ls[i])))
				data["en"]["faker"]=merge(data["en"]["faker"],YAML.load(open(ls[i]))["en"]["faker"])
			end
		end
	end
	cd(original)
end

loadfiles("base")

function executor(data :: Array)

 fn= data[rand(1:length(data))]

if occursin(r"#{.+}",fn)
	elements = collect((m.match for m = eachmatch(r"#{[a-zA-Z\_\.]+}", fn)))
	for i in elements
      fn= replace(fn, i => eval(Meta.parse(i[3:end-1]*"()")),count=1)
    end
end
return fn

end

include("base.jl")
include("persona.jl")
include("creditcard.jl")
include("datatime.jl")
include("misc.jl")
include("phone_number.jl")
include("color.jl")
include("company.jl")
include("file.jl")
include("internet.jl")
include("lorem.jl")
include("direccion.jl")
include("user_agent.jl")

end # module


text:
## QUESTION 2

include("md.jl");

fname = "../../swp/1069_auto.mdv"

# Create a new simulation object
q = MolecularDynamicsTrial(
    1000,               # number of bodies
    0.75,               # particle density
    1.069,              # desired temperature; also use 1.304
    0.032,              # time step dx, in units of τ
    20000,              # number of iterations
    3.405e-10,          # length units, from Verlet, in m
    119.8,              # energy scale for potential, in Kelvin
    39.948 / 6.0221e26  # argon atomic mass, in kg
);

# Place the particles
mdplace!(q);

# Run for 199 steps to thermalize
verlet!(q, 199)

# Save the trial
save(fname, q)

# Relax
relax!(q)

# Save the relaxed trial
save(fname, q)

# Thermalize and relax 19 more times
thermalize(fname, q, 200, 19)

# Run for another 16000 steps
for i in 1:32
    verlet!(q, 500)
    save(fname, q)
end

# Extract information
extract("../../swp/set02", q, 10)


text:
export PyPtr, PyObject, Py_NULL, METH_O, METH_NOARGS, METH_FASTCALL
export Py_ssize_t, Addr
export PyMethodDef
export Py_LE, Py_LT, Py_NE, Py_EQ, Py_GE, Py_GT

const Py_ssize_t = Cssize_t
const Py_hash_t = Py_ssize_t
const Addr = UInt64

# Ref: https://github.com/python/cpython/blob/3.8/Include/object.h#L104
struct PyObject
    ob_refcnt::Py_ssize_t
    ob_type::Ptr{PyObject}
end

const PyPtr = Ptr{PyObject}
const Py_NULL = reinterpret(PyPtr, C_NULL)

const METH_VARARGS = 0x0001
const METH_KEYWORDS = 0x0002
const METH_NOARGS = 0x0004
const METH_O = 0x0008
const METH_FASTCALL = 0x0080

# Ref: https://github.com/python/cpython/blob/3.8/Include/methodobject.h#L51
mutable struct PyMethodDef
    ml_name::Cstring
    ml_meth::Ptr{Nothing}
    ml_flags::Cint
    ml_doc::Cstring # may be NULL
end

# Ref: https://github.com/python/cpython/blob/3.8/Include/descrobject.h#L11
struct PyGetSetDef
    name::Ptr{UInt8}
    get::Ptr{Cvoid}
    set::Ptr{Cvoid} # may be NULL for read-only members
    doc::Ptr{UInt8} # may be NULL
    closure::Ptr{Cvoid} # pass-through thunk, may be NULL
end

# (o, o) -> p
# struct binaryfunc
#     unbox::Ptr{Nothing}
# end
# function (f::binaryfunc)(o1::PyPtr, o2::PyPtr)
#     ccall(f.unbox, PyPtr, (PyPtr, PyPtr), o1, o2)
# end
# DIO_ExceptCode(::binaryfunc) = Py_NULL
# (o) -> p
# struct unaryfunc
#     unbox::Ptr{Nothing}
# end
# function (f::unaryfunc)(o)
#     ccall(f.unbox, PyPtr, (PyPtr,), o)
# end
# DIO_ExceptCode(::unaryfunc) = Py_NULL
# # (o, o, o) -> p
# struct ternaryfunc
#     unbox::Ptr{Nothing}
# end
# function (f::ternaryfunc)(o1::PyPtr, o2::PyPtr, o3::PyPtr)
#     ccall(f.unbox, PyPtr, (PyPtr, PyPtr, PyPtr), o1, o2, o3)
# end
# DIO_ExceptCode(::ternaryfunc) = Py_NULL

const binaryfunc = Ptr{Nothing}
const unaryfunc = Ptr{Nothing}
const ternaryfunc = Ptr{Nothing}

# Ref: https://github.com/python/cpython/blob/3.8/Include/cpython/object.h#L95
struct PyNumberMethods
    nb_add::binaryfunc
    nb_subtract::binaryfunc
    nb_multiply::binaryfunc
    nb_remainder::binaryfunc
    nb_divmod::binaryfunc
    nb_power::ternaryfunc
    nb_negative::unaryfunc
    nb_positive::unaryfunc
    nb_absolute::unaryfunc
    nb_bool::Ptr{Nothing}
    nb_invert::unaryfunc
    nb_lshift::binaryfunc
    nb_rshift::binaryfunc
    nb_and::binaryfunc
    nb_xor::binaryfunc
    nb_or::binaryfunc
    nb_int::unaryfunc
    nb_reserved::Ptr{Nothing}
    nb_float::unaryfunc
    nb_inplace_add::binaryfunc
    nb_inplace_subtract::binaryfunc
    nb_inplace_multiply::binaryfunc
    nb_inplace_remainder::binaryfunc
    nb_inplace_power::ternaryfunc
    nb_inplace_lshift::binaryfunc
    nb_inplace_rshift::binaryfunc
    nb_inplace_and::binaryfunc
    nb_inplace_xor::binaryfunc
    nb_inplace_or::binaryfunc
    nb_floor_divide::binaryfunc
    nb_true_divide::binaryfunc
    nb_inplace_floor_divide::binaryfunc
    nb_inplace_true_divide::binaryfunc
    nb_index::unaryfunc
    nb_matrix_multiply::binaryfunc
    nb_inplace_matrix_multiply::binaryfunc
end

struct PyMemberDef
    name::Ptr{UInt8}
    typ::Cint
    offset::Int # warning: was Cint for Python <= 2.4
    flags::Cint
    doc::Ptr{UInt8}
end

# Ref: https://github.com/python/cpython/blob/3.8/Include/cpython/object.h#L177
mutable struct PyTypeObject
    # PyObject_HEAD (for non-Py_TRACE_REFS build):
    ob_refcnt::Py_ssize_t
    ob_type::PyPtr
    ob_size::Py_ssize_t # PyObject_VAR_HEAD

    # PyTypeObject fields:
    tp_name::Ptr{UInt8} # required, should be in format "<module>.<name>"

    # warning: these two were Cint for Python <= 2.4
    tp_basicsize::Int # required, = sizeof(instance)
    tp_itemsize::Int

    tp_dealloc::Ptr{Cvoid}
    tp_vectorcall_offset :: Addr
    tp_getattr::Ptr{Cvoid}
    tp_setattr::Ptr{Cvoid}
    tp_as_async::Ptr{Cvoid}
    tp_repr::Ptr{Cvoid}

    tp_as_number::Ptr{PyNumberMethods}
    tp_as_sequence::Ptr{Cvoid}
    tp_as_mapping::Ptr{Cvoid}

    tp_hash::Ptr{Cvoid}
    tp_call::Ptr{Cvoid}
    tp_str::Ptr{Cvoid}
    tp_getattro::Ptr{Cvoid}
    tp_setattro::Ptr{Cvoid}

    tp_as_buffer::Ptr{Cvoid}

    tp_flags::Clong # Required, should default to Py_TPFLAGS_DEFAULT

    tp_doc::Ptr{UInt8} # normally set in example code, but may be NULL

    tp_traverse::Ptr{Cvoid}

    tp_clear::Ptr{Cvoid}

    tp_richcompare::Ptr{Cvoid}

    tp_weaklistoffset::Int

    # added in Python 2.2:
    tp_iter::Ptr{Cvoid}
    tp_iternext::Ptr{Cvoid}

    tp_methods::Ptr{PyMethodDef}
    tp_members::Ptr{PyMemberDef}
    tp_getset::Ptr{PyGetSetDef}
    tp_base::Ptr{Cvoid}

    tp_dict::PyPtr
    tp_descr_get::Ptr{Cvoid}
    tp_descr_set::Ptr{Cvoid}
    tp_dictoffset::Int

    tp_init::Ptr{Cvoid}
    tp_alloc::Ptr{Cvoid}
    tp_new::Ptr{Cvoid}
    tp_free::Ptr{Cvoid}
    tp_is_gc::Ptr{Cvoid}

    tp_bases::PyPtr
    tp_mro::PyPtr
    tp_cache::PyPtr
    tp_subclasses::PyPtr
    tp_weaklist::PyPtr
    tp_del::Ptr{Cvoid}

    # added in Python 2.6:
    tp_version_tag::Cuint

    # only used for COUNT_ALLOCS builds of Python
    tp_allocs::Int
    tp_frees::Int
    tp_maxalloc::Int
    tp_prev::Ptr{Cvoid}
    tp_next::Ptr{Cvoid}
end

const Py_LT = Cint(0)
const Py_LE = Cint(1)
const Py_EQ = Cint(2)
const Py_NE = Cint(3)
const Py_GT = Cint(4)
const Py_GE = Cint(5)


text:
using NumDiff
using Base.Test

"""Векторозначная неавтономная функция от многих переменных
f(t, x_1, x_2,...,x_n) -> [f1, f2,...,fn]
"""
function f1(t::Float64, x::Vector{Float64})::Vector{Float64}
  return x .^ 2 + t*x
end

"""Частные производные по всем векторным аргументам функции f1"""
function dxf1(t::Float64, x::Vector{Float64})::Matrix{Float64}
  n = length(x)
  return diagm(2.0 .* x) + t * eye(n,n)
end

"""Векторозначная автономная функция от многих переменных
f(x_1, x_2,...,x_n) -> [f1, f2,...,fn]
"""
function f2(x::Vector{Float64})::Vector{Float64}
  return x .^2
end

"""Частные производные по всем векторным аргументам функции f2"""
function dxf2(x::Vector{Float64})::Matrix{Float64}
  return diagm(2.0 .* x)
end

"""Скалярная автономная функция от многих переменных
f(x_1, x_2,...,x_n) -> f::Float64
"""
function f3(x::Vector{Float64})::Float64
  return dot(x, x)
end

"""Частные производные по всем векторным аргументам функции f3"""
function dxf3(x::Vector{Float64})::Vector{Float64}
  return 2.0 .* x
end

"""Скалярная неавтономная функция от многих переменных
f(t, x_1, x_2,...,x_n) -> f::Float64
"""
function f4(t::Float64, x::Vector{Float64})::Float64
  return dot(x, x) + t
end

"""∂f4/∂x1, ∂f4/∂x2,...,∂f4/∂xn"""
function dxf4(t::Float64, x::Vector{Float64})::Vector{Float64}
  return 2.0 .* x
end

"""Векторознаяная функция от одной переменной"""
function f5(x::Float64)::Vector{Float64}
  return [1.0, 2.0, 3.0] * x
end

"""∂f5/∂x1, ∂f5/∂x2,...,∂f5/∂xn"""
function dxf5(x::Float64)::Vector{Float64}
  return [1.0, 2.0, 3.0]
end

# Тестирование
t_0 = 0.5
x_0 = [1.0, 2.0, 3.0]
# для всех порядков
@testset for od in 2:2:8
  # для всех переменных
  @testset for i in 1:length(x_0)
    @test NumDiff.derivative(f1, t_0, x_0; k=i, order=od) ≈ dxf1(t_0, x_0)[:, i]
    @test NumDiff.derivative(f2, x_0; k=i, order=od) ≈ dxf2(x_0)[:, i]
    @test NumDiff.derivative(f3, x_0; k=i, order=od) ≈ dxf3(x_0)[i]
    @test NumDiff.derivative(f4, t_0, x_0; k=i, order=od) ≈ dxf4(t_0, x_0)[i]
  end
  @test NumDiff.derivative(f5, 2.0; order=8) ≈ dxf5(2.0)
end


text:
using Test
using QuantumOptics, CollectiveSpins

@testset "mpc" begin

cs = CollectiveSpins

# System parameters
function test_2spin()
    a = 0.54
    γ = 1
    e_dipole = [0, 0, 1]
    T = [0:0.05:5.;]
    phi = [0.3, 0.8]
    theta = [1.0/2*pi, 3.0/2*pi]

    system = SpinCollection(Vector{Float64}[[0, 0, 0], [a, 0, 0]], e_dipole; gammas=γ)
    N = length(system.spins)

    # Meanfield + Correlations
    state0 = cs.mpc.blochstate(phi, theta)
    tout, state_mpc_t = cs.mpc.timeevolution(T, system, state0)

    # Quantum: master equation
    function fout(t, rho::AbstractOperator)
        i = something(findfirst(isequal(t), T), 0)
        rho_mpc = cs.mpc.densityoperator(state_mpc_t[i])
        @test tracedistance(rho, rho_mpc)<1e-5
    end

    Ψ₀ = cs.quantum.blochstate(phi, theta)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    cs.quantum.timeevolution(T, system, ρ₀, fout=fout)
end

function test_3spin()
    a = 0.54
    γ = 1.
    e_dipole = [0, 0, 1]
    T = [0:0.05:5.;]
    phi = [0.3, 0.8, 1.6]
    theta = [1.0/2*pi, 3.0/2*pi, 1.2*pi]

    system = SpinCollection(cs.geometry.triangle(a), e_dipole; gammas=γ)
    N = length(system.spins)

    # Meanfield
    state0 = cs.meanfield.blochstate(phi, theta)
    tout, state_mf_t = cs.meanfield.timeevolution(T, system, state0)

    # Meanfield + Correlations
    state0 = cs.mpc.blochstate(phi, theta)
    tout, state_mpc_t = cs.mpc.timeevolution(T, system, state0)

    # Quantum: master equation
    function fout(t, rho::AbstractOperator)
        i = something(findfirst(isequal(t), T), 0)
        rho_mf = cs.meanfield.densityoperator(state_mf_t[i])
        rho_mpc = cs.mpc.densityoperator(state_mpc_t[i])
        td_mf = tracedistance(rho, rho_mf)
        td_mpc = tracedistance(rho, rho_mpc)
        if (i>1)
            @test td_mpc < td_mf
        end
        @test td_mpc<0.05
    end

    Ψ₀ = cs.quantum.blochstate(phi, theta)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    cs.quantum.timeevolution(T, system, ρ₀, fout=fout)
end

test_2spin()
test_3spin()

end # testset


text:
include("read.jl")
include("write.jl")


text:
include("common.jl")

mmesh = Mesh(joinpath(PDATA,"rectangle.stl"))
c1=(1.0, 0.0)
c2=(1.0, 1.0)
fault_idx, fault_coords, fault_normal_vec = get_fault_info(mmesh, c1, c2)

mmesh, conn = fix_mesh(mmesh, fault_idx, rectangle_redblack)
close("all")
visualize_mesh(mmesh)
v = 0.05*(rand(mmesh.nnode).-0.5)
u = 0.05*(rand(mmesh.nnode).-0.5)
scatter(mmesh.nodes[:,1] + u, mmesh.nodes[:,2] + v, s = 2)
savefig("mesh.png")

bdnode = bcnode((x, y)->(x<1e-5 || x>2-1e-5), mmesh)

Conn = compute_slip_boundary_condition( conn, mmesh)
NT = 100
Δt = 0.004
λ = 1.0
μ = 1.0
invη = ones(get_ngauss(mmesh))
sn = size(conn,1)
slipvec = [zeros(NT, sn) ones(NT, sn)] # acceleration
gravity = 0.0

U, Sigma, Varepsilon = subduction_solver(invη, gravity, Conn, slipvec, bdnode, mmesh)

sess = Session(); init(sess)
Uval, Stress, Strain = run(sess, [U, Sigma, Varepsilon])
p = visualize_displacement(Uval[1:5:end,:], mmesh)
saveanim(p, "subduction_rectangle.gif")

close("all")
u = Uval[end,:]
visualize_scalar_on_fem_points(u[mmesh.ndof+1:2mmesh.ndof],mmesh)
savefig("y_displacement.png")


text:
export make_lattice, makelattice
export Lattice
export dimension

"""
    Lattice{O}

Represent a lattice.

# Arguments
* `unitcell::UnitCell{O}`
* `orthocube::OrthoCube`
* `bravais_coordinates::Vector{Vector{Int}}`
* `supercell::UnitCell{Tuple{O, Vector{Int}}}`
"""
struct Lattice{O}
    unitcell::UnitCell{O}
    orthocube::OrthoCube
    bravais_coordinates::Vector{Vector{Int}}
    supercell::UnitCell{Tuple{O, Vector{Int}}}
end

function Base.:(==)(lhs::Lattice{O}, rhs::Lattice{O}) where O
    return (
        lhs.unitcell == rhs.unitcell
        && lhs.orthocube == rhs.orthocube
        && lhs.bravais_coordinates == rhs.bravais_coordinates
    )
end

function makelattice(
    unitcell::UnitCell{O},
    shape_matrix::AbstractMatrix{<:Integer}
) where O
    dim = dimension(unitcell)
    if size(shape_matrix) != (dim, dim)
        throw(DimensionMismatch("unitcell and shape_matrix should have the same dimension"))
    end

    orthocube = OrthoCube(shape_matrix)
    generator_translations = find_generators(orthocube)
    unitcell_coordinates = generate_coordinates(orthocube, generator_translations)

    new_latticevectors = unitcell.latticevectors * orthocube.shape_matrix

    new_unitcell = makeunitcell(new_latticevectors; SiteType=Tuple{O, Vector{Int}})
    for uc in unitcell_coordinates, (orbname, orbcoord) in unitcell.sites
        cc = fract2carte(unitcell, orbcoord)
        new_cc = cc + unitcell.latticevectors * uc
        new_orbcoord = carte2fract(new_unitcell, new_cc)
        new_orbname = (orbname, uc)
        addsite!(new_unitcell, new_orbname, new_orbcoord)
    end

    return Lattice{O}(unitcell, orthocube, unitcell_coordinates, new_unitcell)
end
#TODO unit testing for lattice with wrong dimensions

"""
    makelattice(unitcell, shape, [generators])

Create a lattice with periodic boundary condition, using the unitcell, shape, and
translation generators.

# Arguments
* `unitcell::UnitCell{O}`
* `shape_matrix::AbstractMatrix{<:Integer}`: shape of the Bravais lattice. This can also be
a single integer, which is equivalent to a identity matrix times the number.

# Optional Arguments
* `generator_translations::AbstractMatrix{<:Integer}`
"""
function makelattice(
    unitcell::UnitCell{O},
    shape_matrix::AbstractMatrix{<:Integer},
    generator_translations::AbstractMatrix{<:Integer},
) where O
    dim = dimension(unitcell)
    if size(shape_matrix) != (dim, dim)
        throw(DimensionMismatch("unitcell and shape_matrix should have the same dimension"))
    end
    orthocube = OrthoCube(shape_matrix)
    unitcell_coordinates = generate_coordinates(orthocube, generator_translations)

    new_latticevectors = unitcell.latticevectors * orthocube.shape_matrix

    new_unitcell = makeunitcell(new_latticevectors; SiteType=Tuple{O, Vector{Int}})
    for uc in unitcell_coordinates, (orbname, orbcoord) in unitcell.sites
        cc = fract2carte(unitcell, orbcoord)
        new_cc = cc + unitcell.latticevectors * uc
        new_orbcoord = carte2fract(new_unitcell, new_cc)
        new_orbname = (orbname, uc)
        addsite!(new_unitcell, new_orbname, new_orbcoord)
    end

    return Lattice{O}(unitcell, orthocube, unitcell_coordinates, new_unitcell)
end

function makelattice(unitcell::UnitCell, scale::Integer)
    return makelattice(unitcell, scale*ones(Int, (1,1)))
end

make_lattice = makelattice

"""
    dimension(lattice)

Spatial dimension of the lattice.

# Arguments
- `lattice::Lattice`
"""
dimension(lattice::Lattice) = dimension(lattice.unitcell)


text:
using FactCheck, Images, Color, FixedPointNumbers
@osx_only using TestImages

facts("Read stream") do
    @osx_only begin
        # OS X only for the moment, because TestImages is only
        # required for OS X (see test/REQUIRE).  That seems to cause
        # a testing failure on Travis for Linux.
        # I assume that at some point, streaming will expand and other tests
        # will be needed, at which point we can open up this test
        workdir = joinpath(tempdir(), "ImagesStreamTest")
        testfile = joinpath(workdir, "jigsaw_tmpl.png")
        if !isfile(testfile)
            mkpath(workdir)
            download("http://www.imagemagick.org/Usage/images/jigsaw_tmpl.png", testfile)
        end
        open(testfile, "r") do io
            img = imread(io, Images.OSXNative)
        end
        @pending something --> true
        rm(workdir, recursive=true)
    end
end


text:
using GLFW, Test

window = GLFW.CreateWindow(800, 600, "InexactError")
@test !GLFW.WindowShouldClose(window)
# If you call the julia function GLFW.SetWindowShouldClose(window,
# truefalse), then only 0 or 1 can be passed. However, at least on
# some platforms (Linux), clicking on the window's close icon can send
# other values (like 189). Make sure such values don't cause trouble for
# GLFW.WindowShouldClose.
ccall((:glfwSetWindowShouldClose, GLFW.libglfw), Cvoid, (GLFW.Window, Cint), window, 189)
@test GLFW.WindowShouldClose(window)
GLFW.DestroyWindow(window)


text:
"""
    Reaction{C,U}

Contain a "calculate" function and an "update" function.
The "calculate" function calculates the a "rate" of reaction with the system state,
which determines the probability weight of the reaction be selected,
and the "update" function updates the system state,
when this reaction were selected randomly.
"""
struct Reaction{C,U}
    c::C # function to calculate "rate"
    u::U # function to update state
end

# generated type stable map for tuple,
# Base.map for large Tuple is not type stable
# however, this package must process large Tuples
@generated function gmap(f, tp::Tuple)
    N = length(tp.parameters)
    return Expr(:tuple, (:(f(tp[$i])) for i in 1:N)...)
end

# generated type stable accumulate for tuple
# Base.accumulate for tuple return a tuple
# this method return a Vector
@generated function gaccumulate(tp::Tuple)
    ret_type = promote_type(tp.parameters...)
    N = length(tp.parameters)
    return quote
        ret = Vector{$ret_type}(undef, $N)
        ret[1] = tp[1]
        $((:(@inbounds ret[$i] = ret[$(i - 1)] + tp[$i]) for i in 2:N)...)
        return ret
    end
end

# apply the given reaction to update system state
# generate if-elseif branch for each reaction for type stable code
# as[i], rs[i] are type unstable without @generated
@generated function applyreaction(hook!, rng::AbstractRNG, t, rs, ps, as, acc, ind, rn)
    N = length(rs.parameters)  # number of reactions
    sub = Symbol("sub_$N")
    ex = Expr(
        :elseif,
        :(ind == $N),
        quote
            $sub = sample(as[$N], rn - acc[$(N - 1)])
            (rs[$N].u)(rng, t, $sub, ps)
            return hook!(rng, t, $sub, ps)
        end,
        :(return :finnish) # this else is never executed but for type stable
    ) # last elseif block with the last reaction
    for i in N-1:-1:2
        sub = Symbol("sub_$i")
        ex = Expr(
            :elseif,
            :(ind == $i),
            quote
                $sub = sample(as[$i], rn - acc[$(i - 1)])
                (rs[$i].u)(rng, t, $sub, ps)
                return hook!(rng, t, $sub, ps)
            end,
            ex,
        ) # append elseif block form last to the second
    end
    return Expr(:if, :(ind == 1), quote
        sub_1 = sample(as[1], rn)
        (rs[1].u)(rng, t, sub_1, ps)
        return hook!(rng, t, sub_1, ps)
    end, ex) # append if block at first
end

"""
    gillespie!(hook!, rng::AbstractRNG, c::ContinuousClock, ps::NamedTuple, rs::Tuple)

Simulate the system using the Gillespie algorithm with the given parameters,
and return a tuple of end time
and the terminate state `:finnish`, `:zero` or any other state returns by the `hook!`.
The terminate state `:finnish` means that simulation reach to the end time,
and `:zero` means the simulation break because the total "reaction rate" is zero,
besides, `hook!` should return a symbol terminate state like `:break`.
if the return value of `hook!` is not `:finnish`, the simulation will be terminated.
The clock `c` and parameters `ps` will be updated during the simulation.

# Arguments

- `hook!`: a function with similar arguments to "update" functions,
   and it's recommended to create `hook!` with [`@ufunc`](@ref) macro;
   unlike "update" functions, `hook` will be called after each reaction
   and should return a terminate state used to terminate the simulation if it is not `:finnish`.
- `rng`: a random number generator for generate random numbers;
- `c`: a clock recording time, which must be the reference clock of recorded variables in `ps`;
- `ps`: a NamedTuple contains state, parameters even args used by `hook!` of the system;
- `rs`: a tuple contains reactions, all parameters required by reactions must be in `ps` with same name.
"""
function gillespie!(hook!, rng::AbstractRNG, c::ContinuousClock, ps::NamedTuple, rs::Tuple)
    term_state = :finnish # terminate state
    t = currenttime(c)
    for _ in c
        # use t before update it may cause type un-stability
        as = gmap(r -> (r.c)(currenttime(c), ps), rs)   # calculate "rate" for each reaction
        as_sum = gmap(sum, as)             # sum of "rate" for each reaction
        as_sum_acc = gaccumulate(as_sum)   # accumulate of the sum of "rate" for all reactions
        as_sum_sum = as_sum_acc[end]       # sum of "rate" for all reactions
        if iszero(as_sum_sum)              # if all "rate" are zero
            term_state = :zero             # mark terminate state as break
            break                          # break the loop
        end
        τ = -log(rand(rng, typeof(as_sum_sum))) / as_sum_sum   # calculate τ
        t = increase!(c, τ) # increase clock time
        rn = rand(rng, typeof(as_sum_sum)) * as_sum_sum # generate random number, use this rand number twice is OK
        ind = searchsortedfirst(as_sum_acc, rn)::Int # sample reaction index, use _sample instead of sample to return Int
        term_state = applyreaction(hook!, rng, t, rs, ps, as, as_sum_acc, ind, rn) # apply reaction of index ind to update system state
        term_state == :finnish || break # break if state is not :finnish
    end
    return t, term_state
end

"""
    gillespie([hook!, rng::AbstractRNG,] c, ps::NamedTuple, rs::Tuple)

Simulate the system using the Gillespie algorithm with the given parameters,
and return a tuple of updated `ps`, end time and terminate state.
More about terminate state, see [`gillespie!`](@ref).

# Arguments

- `hook!`: a function with similar arguments to "update" functions
   and recommended to created with [`@ufunc`](@ref) macro;
   unlike "update" functions, `hook` will be called after each reaction
   and should return a terminate state used to terminate the simulation if it is not `:finnish`.
- `rng`: a random number generator for generate random numbers;
- `c`: a `ContinuousClock`, a end time or a tuple of a begin and a end time;
- `ps`: a NamedTuple contains state, parameters even args used by `hook!` of the system;
- `rs`: a tuple contains reactions, all parameters required by reactions must be in `ps` with same name.

!!! info
    In order to update single value state,
    it must be converted to a [`ScalarType`](@ref Scalar.ScalarType) with [`scalar`](@ref Scalar.scalar).
"""
function gillespie(hook!, rng::AbstractRNG, c, ps::NamedTuple, rs::Tuple)
    # process c and ps
    # 1. deepcopy to avoid change the original values
    # 2. if c is a real or a tuple, convert it to ContinuousClock
    c′, ps′ = _process_args(c, ps)
    t, term_state = gillespie!(hook!, rng, c′, ps′, rs)
    return ps′, t, term_state
end
gillespie(c, ps::NamedTuple, rs::Tuple) =
    gillespie((_...) -> :finnish, Random.GLOBAL_RNG, c, ps, rs)
gillespie(rng::AbstractRNG, c, ps::NamedTuple, rs::Tuple) =
    gillespie((_...) -> :finnish, rng, c, ps, rs)
gillespie(hook!, c, ps::NamedTuple, rs::Tuple) =
    gillespie(hook!, Random.GLOBAL_RNG, c, ps, rs)

_process_args(c::ContinuousClock, ps) = deepcopy((c, ps))
_process_args(c::Real, ps) = ContinuousClock(c), deepcopy(ps)
_process_args((s, e)::Tuple{Real,Real}, ps) = ContinuousClock(e, s), deepcopy(ps)


text:
module AdminController

using Genie, Genie.Renderer, Genie.Requests, Genie.Exceptions, Genie.Renderer.Html, Genie.Flash
using SearchLight, SearchLight.Relationships

using CRUDModels, CRUDOperations
using ViewHelper

using Contacts, Clients, Species, Tissuetypes
using Markers

export create, read, update, delete, list
export sidebar_links

import Base: parse
parse(T::Type{String}, str::String) = str
export parse

################################################################################
sidebar_links = [:client, :contact, :specie, :tissuetype]

################################################################################
################################################################################
function index()
    html(:admin, :index; context = @__MODULE__, sidebar_links = sidebar_links)
end

function create()
    model = get_model(payload(:model))
    field_types = model_field_types(model)
    attrs = Dict(field => parse(field_types[field], postpayload(field))
                 for field in intersect(keys(field_types), fieldnames(model)))
    object = model(; attrs...) |> save!
    flash("Object created successfully.")
    redirect(:admin_model_read, model = nameof(model), id = object.id)
end

function read()
    model = get_model(payload(:model); context = @__MODULE__)
    object = findone(model, id = payload(:id))
    if isnothing(object)
        flash("Object not found.")
        redirect(:admin_model_list, model = nameof(model))
    else
        html(:admin, :read; context = @__MODULE__, model = model, object = object, sidebar_links = sidebar_links)
    end
end

function update()
    #model = get_model(payload(:model))
    #object = findone(model, id = payload(:id))
    #if !isnothing(object)
    #    field_types = model_field_types(model)
    #    attrs = Dict(field => parse(field_types[field], postpayload(field))
    #                 for field in intersect(keys(field_types), fieldnames(model)))
    #    object = model_update(model, object; attrs)
    #    flash("Object updated successfully.")
    #else
    #    flash("Object not found.")
    #end
    redirect(:admin_model_read, model = model_symbol(model), id = object.id)
end

function delete()
    model = get_model(payload(:model))
    object = findone(model, id = payload(:id))
    if !isnothing(object)
        SearchLight.delete!(object)
        flash("Object deleted successfully.")
    else
        flash("Object not found.")
    end
    redirect(:admin_model_list, model = nameof(model))
end

function list()
    model = get_model(payload(:model); context = @__MODULE__)
    objects = all(model)
    html(:admin, :list; context = @__MODULE__, model = model, objects = objects, sidebar_links = sidebar_links)
end

end


text:
module makeSlideWindows

using BioSequences, DataFrames

export getGCcontentInSlideWindowView

function getGCcontentInSlideWindowView(seq::BioSequence{DNAAlphabet{4}}, windowSize::Int, asProb::Bool, seqName::String)
	inputSeqLen = length(seq)
	@assert inputSeqLen >= windowSize
	restSeqLen = inputSeqLen % windowSize
	windowsNum = restSeqLen == 0 ? inputSeqLen ÷ windowSize : inputSeqLen ÷ windowSize + 1
	GCcontent = asProb ? zeros(Float32, windowsNum) : zeros(Int32, windowsNum)
	StartPos = zeros(Int32, windowsNum)
	StopPos = zeros(Int32, windowsNum)
	Chromosome = fill(seqName, windowsNum)
	for i = 0:(windowsNum-1)
		idx = i + 1
		start = i * windowSize + 1
		stop = restSeqLen == 0 ? (i + 1) * windowSize :
		                         ( i == (windowsNum - 1) ? i * windowSize + restSeqLen : (i+1) * windowSize)
		StartPos[idx] = start
		StopPos[idx] = stop
		if i == (windowsNum - 1)
			@assert stop == inputSeqLen
		end
		subseq = seq[start:stop]
		actgComposition = BioSequences.Composition(subseq).counts
		numG = haskey(actgComposition, DNA_G) ? actgComposition[DNA_G] : 0
		numC = haskey(actgComposition, DNA_C) ? actgComposition[DNA_C] : 0
		numGC = numG + numC
		GCcontent[idx] = asProb ? round(numGC/windowSize, digits = 4) : numGC
	end
	output = DataFrame(Chrom = Chromosome, Start = StartPos, Stop = StopPos, GC = GCcontent)
	return output
end

end


text:
"""
    VREP{N, AM, S, D} <: AbstractContinuousPost

Implementation of a linear reachability solver using the vertex representation.
"""
struct VREP{N, AM, S, D} <: AbstractContinuousPost
    δ::N
    approx_model::AM
    static::S
    dim::D
end

# convenience constructor using symbols
function VREP(; δ::N,
                dim::Union{Int, Missing}=missing,
                static::Bool=!ismissing(dim),
                backend=nothing,
                approx_model=Forward(sih=:concrete, exp=:base, setops=:vrep, backend=backend)) where {N}

    n = ismissing(dim) ? missing : Val(dim)
    return VREP(δ, approx_model, Val(static), n)
end

step_size(alg::VREP) = alg.δ
numtype(::VREP{N}) where {N} = N

function setrep(alg::VREP{N, AM, Val{false}, Missing}) where {N, AM}
    return VPolytope{N, Vector{N}}
end

function setrep(alg::VREP{N, AM, Val{false}, D}) where {N, AM, D}
    return VPolygon{N, Vector{N}}
end

function setrep(alg::VREP{N, AM, Val{true}, Missing}) where {N, AM}
    error("the set representation of this algorithm requires the dimension to be specified, but it is $(alg.dim)")
end

function setrep(alg::VREP{N, AM, Val{true}, Val{2}}) where {N, AM}
    return VPolygon{N, SVector{2, N}}
end

function setrep(alg::VREP{N, AM, Val{true}, Val{n}}) where {N, AM, n}
    return VPolytope{N, SVector{n, N}}
end

function rsetrep(alg::VREP{N}) where {N}
    ST = setrep(alg)
    RT = ReachSet{N, ST}
end

include("post.jl")
include("reach_homog.jl")


text:
# https://github.com/JuliaNLSolvers/Optim.jl
# https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/

time_start = time()

using PowerModels
using Optim
using JuMP, Ipopt

pkg_load_time = time() - time_start

time_start = time()

file_name = "../data/pglib_opf_case5_pjm.m"
#file_name = "../data/pglib_opf_case118_ieee.m"

data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

data_load_time = time() - time_start

time_start = time()

model = Model(Ipopt.Optimizer)
#set_optimizer_attribute(model, "print_level", 0)

@variable(model, jmp_va[i in keys(ref[:bus])])
@variable(model, ref[:bus][i]["vmin"] <= jmp_vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)

@variable(model, ref[:gen][i]["pmin"] <= jmp_pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
@variable(model, ref[:gen][i]["qmin"] <= jmp_qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])

@variable(model, -ref[:branch][l]["rate_a"] <= jmp_p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
@variable(model, -ref[:branch][l]["rate_a"] <= jmp_q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])

@objective(model, Min, sum(gen["cost"][1]*jmp_pg[i]^2 + gen["cost"][2]*jmp_pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]))

for (i,bus) in ref[:ref_buses]
    @constraint(model, jmp_va[i] == 0)
end

for (i,bus) in ref[:bus]
    bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
    bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]

    @constraint(model,
        sum(jmp_p[a] for a in ref[:bus_arcs][i]) ==
        sum(jmp_pg[g] for g in ref[:bus_gens][i]) -
        sum(load["pd"] for load in bus_loads) -
        sum(shunt["gs"] for shunt in bus_shunts)*jmp_vm[i]^2
    )

    @constraint(model,
        sum(jmp_q[a] for a in ref[:bus_arcs][i]) ==
        sum(jmp_qg[g] for g in ref[:bus_gens][i]) -
        sum(load["qd"] for load in bus_loads) +
        sum(shunt["bs"] for shunt in bus_shunts)*jmp_vm[i]^2
    )
end

# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
    f_idx = (i, branch["f_bus"], branch["t_bus"])
    t_idx = (i, branch["t_bus"], branch["f_bus"])

    p_fr = jmp_p[f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
    q_fr = jmp_q[f_idx]                     # q_fr is a reference to the optimization variable q[f_idx]
    p_to = jmp_p[t_idx]                     # p_to is a reference to the optimization variable p[t_idx]
    q_to = jmp_q[t_idx]                     # q_to is a reference to the optimization variable q[t_idx]

    vm_fr = jmp_vm[branch["f_bus"]]         # vm_fr is a reference to the optimization variable vm on the from side of the branch
    vm_to = jmp_vm[branch["t_bus"]]         # vm_to is a reference to the optimization variable vm on the to side of the branch
    va_fr = jmp_va[branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
    va_to = jmp_va[branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch

    g, b = PowerModels.calc_branch_y(branch)
    tr, ti = PowerModels.calc_branch_t(branch)
    tm = tr^2 + ti^2
    g_fr = branch["g_fr"]
    b_fr = branch["b_fr"]
    g_to = branch["g_to"]
    b_to = branch["b_to"]

    # From side of the branch flow
    @NLconstraint(model, p_fr ==  (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
    @NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )

    # To side of the branch flow
    @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
    @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )

    # Voltage angle difference limit
    @constraint(model, va_fr - va_to <= branch["angmax"])
    @constraint(model, va_fr - va_to >= branch["angmin"])

    # Apparent power limit, from side and to side
    @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
    @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
end

optimize!(model)
#cost = objective_value(model)

bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])

bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])

for (i,bus) in ref[:bus]
    if length(ref[:bus_loads][i]) > 0
        bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
        bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
    end

    if length(ref[:bus_shunts][i]) > 0
        bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
        bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
    end
end

br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])

br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_tm = Dict(i => 0.0 for (i,branch) in ref[:branch])

br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])

for (i,branch) in ref[:branch]
    g, b = PowerModels.calc_branch_y(branch)
    tr, ti = PowerModels.calc_branch_t(branch)

    br_g[i] = g
    br_b[i] = b

    br_tr[i] = tr
    br_ti[i] = ti
    br_tm[i] = tr^2 + ti^2

    br_g_fr[i] = branch["g_fr"]
    br_b_fr[i] = branch["b_fr"]
    br_g_to[i] = branch["g_to"]
    br_b_to[i] = branch["b_to"]
end

var_lookup = Dict{String,Int}()

var_init = Float64[]
var_lb = Float64[]
var_ub = Float64[]

var_idx = 1
for (i,bus) in ref[:bus]
    #push!(var_init, 0.0) #va
    push!(var_init, value(jmp_va[i]))
    push!(var_lb, -Inf)
    push!(var_ub, Inf)
    var_lookup["va_$(i)"] = var_idx
    global var_idx += 1

    #push!(var_init, 1.0) #vm
    push!(var_init, value(jmp_vm[i]))
    push!(var_lb, bus["vmin"])
    push!(var_ub, bus["vmax"])
    var_lookup["vm_$(i)"] = var_idx
    global var_idx += 1
end

for (i,gen) in ref[:gen]
    #push!(var_init, 0.0) #pg
    push!(var_init, value(jmp_pg[i])) # non-standard start
    push!(var_lb, gen["pmin"])
    push!(var_ub, gen["pmax"])
    var_lookup["pg_$(i)"] = var_idx
    global var_idx += 1

    #push!(var_init, 0.0) #qg
    push!(var_init, value(jmp_qg[i])) # non-standard start
    push!(var_lb, gen["qmin"])
    push!(var_ub, gen["qmax"])
    var_lookup["qg_$(i)"] = var_idx
    global var_idx += 1
end

for (l,i,j) in ref[:arcs]
    branch = ref[:branch][l]

    #push!(var_init, 0.0) #p
    push!(var_init, value(jmp_p[(l,i,j)]))
    push!(var_lb, -branch["rate_a"])
    push!(var_ub,  branch["rate_a"])
    var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
    global var_idx += 1

    #push!(var_init, 0.0) #q
    push!(var_init, value(jmp_q[(l,i,j)]))
    push!(var_lb, -branch["rate_a"])
    push!(var_ub,  branch["rate_a"])
    var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
    global var_idx += 1
end

@assert var_idx == length(var_init)+1

function opf_objective(x)
    cost = 0.0
    for (i,gen) in ref[:gen]
        pg = x[var_lookup["pg_$(i)"]]
        cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
    end
    return cost
end

function opf_constraints(c,x)
    va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
    vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])

    pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
    qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])

    p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
    q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])

    vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch])
    vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch])
    va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
    va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])

    va_con = [va[i] for (i,bus) in ref[:ref_buses]]

    #     @constraint(model,
    #         sum(p[a] for a in ref[:bus_arcs][i]) ==
    #         sum(pg[g] for g in ref[:bus_gens][i]) -
    #         sum(load["pd"] for load in bus_loads) -
    #         sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
    #     )
    power_balance_p_con = [
       sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) -
       bus_pd[i] -
       bus_gs[i]*vm[i]^2 -
       sum(p[a] for a in ref[:bus_arcs][i])
       for (i,bus) in ref[:bus]
    ]

    #     @constraint(model,
    #         sum(q[a] for a in ref[:bus_arcs][i]) ==
    #         sum(qg[g] for g in ref[:bus_gens][i]) -
    #         sum(load["qd"] for load in bus_loads) +
    #         sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
    #     )
    power_balance_q_con = [
       sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) -
       bus_qd[i] +
       bus_bs[i]*vm[i]^2 -
       sum(q[a] for a in ref[:bus_arcs][i])
       for (i,bus) in ref[:bus]
    ]

    # @NLconstraint(model, p_fr ==  (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
    power_flow_p_from_con = [
       (br_g[l]+br_g_fr[l])/br_tm[l]*vm_fr[l]^2 +
       (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
       (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
       p[(l,i,j)]
       for (l,i,j) in ref[:arcs_from]
    ]

    # @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
    power_flow_p_to_con = [
       (br_g[l]+br_g_to[l])*vm_to[l]^2 +
       (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
       (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
       p[(l,i,j)]
       for (l,i,j) in ref[:arcs_to]
    ]

    # @NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
    power_flow_q_from_con = [
       -(br_b[l]+br_b_fr[l])/br_tm[l]*vm_fr[l]^2 -
       (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
       (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
       q[(l,i,j)]
       for (l,i,j) in ref[:arcs_from]
    ]

    # @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
    power_flow_q_to_con = [
       -(br_b[l]+br_b_to[l])*vm_to[l]^2 -
       (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*cos(va_fr[l]-va_to[l])) +
       (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
       q[(l,i,j)]
       for (l,i,j) in ref[:arcs_to]
    ]

    # @constraint(model, va_fr - va_to <= branch["angmax"])
    # @constraint(model, va_fr - va_to >= branch["angmin"])
    power_flow_vad_con = [
       va_fr[l] - va_to[l]
       for (l,i,j) in ref[:arcs_from]
    ]

    # @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
    power_flow_mva_from_con = [
       p[(l,i,j)]^2 + q[(l,i,j)]^2
       for (l,i,j) in ref[:arcs_from]
    ]

    # @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
    power_flow_mva_to_con = [
       p[(l,i,j)]^2 + q[(l,i,j)]^2
       for (l,i,j) in ref[:arcs_to]
    ]

    c .= [
        va_con...,
        power_balance_p_con...,
        power_balance_q_con...,
        power_flow_p_from_con...,
        power_flow_p_to_con...,
        power_flow_q_from_con...,
        power_flow_q_to_con...,
        power_flow_vad_con...,
        power_flow_mva_from_con...,
        power_flow_mva_to_con...,
    ]

    return c
end

con_lbs = Float64[]
con_ubs = Float64[]

#@constraint(model, va[i] == 0)
for (i,bus) in ref[:ref_buses]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
end

#power_balance_p_con
for (i,bus) in ref[:bus]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
    #push!(con_lbs, -Inf)
    #push!(con_ubs, Inf)
end

#power_balance_q_con
for (i,bus) in ref[:bus]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
    #push!(con_lbs, -Inf)
    #push!(con_ubs, Inf)
end

#power_flow_p_from_con
for (l,i,j) in ref[:arcs_from]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
end

#power_flow_p_to_con
for (l,i,j) in ref[:arcs_to]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
end

#power_flow_q_from_con
for (l,i,j) in ref[:arcs_from]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
end

#power_flow_q_to_con
for (l,i,j) in ref[:arcs_to]
    push!(con_lbs, 0.0)
    push!(con_ubs, 0.0)
end

#power_flow_vad_con
for (l,i,j) in ref[:arcs_from]
    branch = ref[:branch][l]
    push!(con_lbs, branch["angmin"])
    push!(con_ubs, branch["angmax"])
end

#power_flow_mva_from_con
for (l,i,j) in ref[:arcs_from]
    branch = ref[:branch][l]
    push!(con_lbs, -Inf)
    push!(con_ubs, branch["rate_a"]^2)
end

#power_flow_mva_to_con
for (l,i,j) in ref[:arcs_to]
    branch = ref[:branch][l]
    push!(con_lbs, -Inf)
    push!(con_ubs, branch["rate_a"]^2)
end

println("variables: $(length(var_init)), $(length(var_lb)), $(length(var_ub))")
println("constraints: $(length(opf_constraints(zeros(length(con_lbs)), var_init))), $(length(con_lbs)), $(length(con_ubs))")

df = TwiceDifferentiable(opf_objective, var_init)
dfc = TwiceDifferentiableConstraints(opf_constraints, var_lb, var_ub, con_lbs, con_ubs)

model_build_time = time() - time_start

time_start = time()

options = Optim.Options(show_trace=true)

# NOTE: had to change initial guess to be an interior point
res = optimize(df, dfc, var_init, IPNewton(), options)
#res = optimize(df, dfc, var_init, LBFGS(), options) #  StackOverflowError:
#res = optimize(df, dfc, var_init, NelderMead(), options) #  StackOverflowError:
display(res)

sol = res.minimizer
cost = res.minimum

#println(cost)
#println(sol)

sol_eval = opf_constraints(zeros(length(con_lbs)), sol)
vio_lb = [max(v,0) for v in (con_lbs .- sol_eval)]
vio_ub = [max(v,0) for v in (sol_eval .- con_ubs)]
const_vio = vio_lb .+ vio_ub
#println(const_vio)
println("total constraint violation: $(sum(const_vio))")

solve_time = time() - time_start

println("")
println("\033[1mSummary\033[0m")
println("   case......: $(file_name)")
println("   cost......: $(round(Int, cost))")
println("   pkg time..: $(pkg_load_time)")
println("   data time.: $(data_load_time)")
println("   build time: $(model_build_time)")
println("   solve time: $(solve_time)")
println("")

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
@testset "restrict" begin
    imgcol = rand(RGB{Float64}, 5, 6)
    imgcolax = AxisArray(imgcol, :y, :x)
    imgr = restrict(imgcolax, (1,2))
    @test pixelspacing(imgr) == (2,2)
    @test pixelspacing(imgcolax) == (1,1)  # issue https://github.com/JuliaImages/Images.jl/issues/347

    imgr = @inferred(restrict(imgcolax, Axis{:y}))
    @test pixelspacing(imgr) == (2, 1)
    imgr = @inferred(restrict(imgcolax, Axis{:x}))
    @test pixelspacing(imgr) == (1, 2)

    @testset "OffsetArray" begin
        A = ones(4, 5)
        Ao = OffsetArray(A, 0:3, -2:2)
        Aor = restrict(Ao)
        @test axes(Aor) == axes(OffsetArray(restrict(A), 0:2, -1:1))

        # restricting AxisArrays with offset axes
        AA = AxisArray(Ao, Axis{:y}(0:3), Axis{:x}(-2:2))
        @test axes(AA) === axes(Ao)
        AAr = restrict(AA)
        axs = axisvalues(AAr)
        @test axes(axs[1])[1] == axes(Aor)[1]
        @test axes(axs[2])[1] == axes(Aor)[2]
        AAA = AxisArray(Aor, axs)  # just test that it's constructable (another way of enforcing agreement)
    end
end


text:
abstract DimensionKmethod

typealias TupleRanges Tuple{UnitRange{Int64},Vararg{UnitRange{Int64}}}

immutable OnePointLR <: DimensionKmethod end
immutable OnePointR <: DimensionKmethod end
immutable OnePointL <: DimensionKmethod end
export DimensionKmethod, OnePointR, OnePointLR, OnePointL

function dimDiff( u :: Array{Float64}, varidx :: Int64, m :: OnePointLR)

  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx, Ridx) = map((2:(dimlength-1), 1:(dimlength -2), 3:(dimlength))) do  dimll
    pp=map(enumerate(size(u))) do xx
      (xx[1] == varidx) ? dimll : 1:xx[2]
    end
    (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])/2

  (midx, Lidx, Ridx) = map((1:1, 1:1, 2:2)) do  dimll
    pp =map(enumerate(size(u))) do xx
      (xx[1] == varidx) ? dimll : 1:xx[2]
    end
    (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  (midx, Lidx, Ridx) =
      map((dimlength:dimlength, (dimlength-1):(dimlength-1),dimlength:dimlength)) do dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  return df
end

function dimDiff( u :: Array{Float64}, varidx :: Int64, m :: OnePointL)

  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx, Ridx) =
      map((2:dimlength, 1:(dimlength -1), 2:dimlength)) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  (midx, Lidx, Ridx) = map((1:1, 1:1, 2:2)) do  dimll
    pp=map(enumerate(size(u))) do xx
      (xx[1] == varidx) ? dimll : 1:xx[2]
    end
    (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  return df
end

function dimDiff( u :: Array{Float64}, varidx :: Int64, m :: OnePointR)
  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx, Ridx) =
      map((1:(dimlength-1), 1:(dimlength -1), 2:(dimlength))) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : (1:xx[2])
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  (midx, Lidx, Ridx) =
      map((dimlength:dimlength, (dimlength-1):(dimlength-1),dimlength:dimlength)) do dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] - u[Lidx...])

  return df
end

export dimDiff
function dimAverage( u :: Array{Float64}, varidx :: Int64, m :: OnePointR)
  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx, Ridx) =
      map((1:(dimlength-1), 1:(dimlength -1), 2:(dimlength))) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : (1:xx[2])
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] + u[Lidx...])/2.0

  (midx, Lidx, Ridx) =
      map((dimlength:dimlength, (dimlength-1):(dimlength-1),dimlength:dimlength)) do dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] + u[Lidx...])/2.0

  return df
end

function dimAverage( u :: Array{Float64}, varidx :: Int64, m :: OnePointL)
  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx, Ridx) =
      map((2:dimlength, 1:(dimlength -1), 2:dimlength)) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : (1:xx[2])
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] + u[Lidx...])/2.0

  (midx, Lidx, Ridx) = map((1:1, 1:1, 2:2)) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = (u[Ridx...] + u[Lidx...])/2.0

  return df
end
export dimAverage

function dimMoves( u :: Array{Float64}, varidx :: Int64, m:: OnePointL)
  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Lidx) =
      map((2:dimlength, 1:(dimlength -1))) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : (1:xx[2])
        end
        (pp...)
  end
  df[midx...] = u[Lidx...]

  (midx, Lidx) = map((1:1, 1:1)) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = u[Lidx...]

  return df
end

function dimMoves( u :: Array{Float64}, varidx :: Int64, m:: OnePointR)
  df = zeros(size(u))
  dimlength =size(u,varidx)
  local midx :: TupleRanges
  local Lidx :: TupleRanges
  local Ridx :: TupleRanges

  (midx, Ridx) =
      map((1:(dimlength-1), 2:(dimlength))) do  dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : (1:xx[2])
        end
        (pp...)
  end
  df[midx...] = u[Ridx...]

  (midx, Ridx) =
      map((dimlength:dimlength ,dimlength:dimlength)) do dimll
        pp = map(enumerate(size(u))) do xx
          (xx[1] == varidx) ? dimll : 1:xx[2]
        end
        (pp...)
  end
  df[midx...] = u[Ridx...]
  return df
end


text:
using Test
using Random
using UnitTestDesign

# ensure these little functions are mutually exclusive.
micro_set = [
    [0, 0],
    [1, 0],
    [0, 1],
    [1, 1],
    [1, 2]
]
function micro_test(micro_func)
    [micro_func(a, b) for (a, b) in micro_set]
end
@test micro_test(UnitTestDesign.ignores) == [1, 0, 0, 0, 0]
@test micro_test(UnitTestDesign.skips) == [0, 1, 0, 0, 0]
@test micro_test(UnitTestDesign.misses) == [0, 0, 1, 0, 0]
@test micro_test(UnitTestDesign.matches) == [0, 0, 0, 1, 0]
@test micro_test(UnitTestDesign.mismatch) == [0, 0, 0, 0, 1]

### Core comparisons.
# They will all use the same test suite but will give different answers.
compare_suite = [
    # [case, tuple]
    [[1, 0, 0, 0], [2, 0, 0, 0]],  # mismatch
    [[3, 0, 4, 0], [3, 0, 4, 0]],  # exact match
    [[1, 1, 0, 2], [1, 0, 0, 2]],  # cover, partial case
    [[1, 0, 0, 3], [1, 1, 0, 3]],  # incomplete / partial cover
    [[0, 1, 0, 2], [1, 0, 3, 0]],  # crossed
    [[1, 3, 1, 2], [0, 3, 0, 2]]   # cover, complete case
]
function cs_test(compare_func)
    cs_res = zeros(Bool, length(compare_suite))
    code = 0
    for cs_idx in 1:length(compare_suite)
        case = compare_suite[cs_idx][1]
        tuple = compare_suite[cs_idx][2]
        cs_res[cs_idx] = compare_func(case, tuple)
        code <<= 1
        code |= cs_res[cs_idx]
    end
    (cs_res, code)
end
# match everything
@test cs_test(UnitTestDesign.case_compatible_with_tuple)[2] == 31
cpc_res = cs_test(UnitTestDesign.case_partial_cover)
@test cpc_res[2] == 29
cpt_res = cs_test(UnitTestDesign.case_covers_tuple)
@test cpt_res[2] == 25

### coverage_by_parameter
cbp_cases = [
    [[0 0 0; 0 0 0; 0 1 0]', 3, [0, 1, 0]],
    [[0 0 0; 0 0 0; 0 1 0]', 2, [0, 0, 0]],
    [[0 0 7; 0 0 4; 0 1 0]', 3, [0, 1, 2]]
]
for cbp_case in cbp_cases
    mc = UnitTestDesign.MatrixCoverage(collect(cbp_case[1]), cbp_case[2], [2, 2, 7])
    res = UnitTestDesign.coverage_by_parameter(mc)
    @test res == cbp_case[3]
end

mc_opcm = UnitTestDesign.one_parameter_combinations_matrix([2, 2, 3], 2)
@test mc_opcm.remain == 4 * 3

### coverage_by_value

# (coverage matrix, remaining uncovered, arity,
#  parameter, coverage of that parameter)
cbv_cases = [
    [[0 0 0; 0 0 0; 0 0 0]', 3, [2,2,2], 1, [0, 0]],
    [[0 0 0; 0 0 0; 0 0 0]', 3, [2,2,2], 2, [0, 0]],
    [[0 0 0; 0 0 0; 0 0 0]', 3, [2,3,2], 2, [0, 0, 0]], # sees arity of choice
    [[0 1 0; 0 1 0; 0 0 0]', 3, [2,3,2], 2, [2, 0, 0]], # does multiple rows
    [[0 1 0; 0 1 0; 0 3 0]', 3, [2,3,2], 2, [2, 0, 1]], # does multiple columns
    [[0 1 0; 0 1 0; 0 3 0]', 2, [2,3,2], 2, [2, 0, 0]]  # excludes extra rows
]
for cbv_case in cbv_cases
    mc = UnitTestDesign.MatrixCoverage(collect(cbv_case[1]), cbv_case[2], cbv_case[3])
    res = UnitTestDesign.coverage_by_value(mc, cbv_case[4])
    @test res == cbv_case[5]
end

### most_matches_existing(allc, row_cnt, arity, existing, param_idx)
mmm_allc = [1 1 1 0 0; 1 2 1 0 0; 0 0 1 2 2]'
mmm_arity = [2, 2, 2, 2, 2]
mmm_cases = [
    [[1, 0, 0, 0, 0], 2, [1, 1]],
    [[1, 0, 0, 0, 0], 4, [0, 0]],
    [[0, 0, 0, 0, 1], 2, [0, 0]],
    [[0, 0, 0, 0, 1], 3, [0, 0]],
    [[0, 0, 0, 0, 2], 3, [1, 0]],
    [[1, 0, 0, 0, 0], 3, [2, 0]],
    [[0, 0, 1, 0, 0], 1, [2, 0]],
    [[1, 0, 1, 0, 0], 2, [1, 1]],
    [[1, 1, 0, 0, 0], 3, [1, 0]],
    [[1, 2, 0, 0, 0], 3, [1, 0]],
    [[0, 0, 1, 0, 2], 4, [0, 1]]
]
for mmm_case in mmm_cases
    mc = UnitTestDesign.MatrixCoverage(collect(mmm_allc), 3, mmm_arity)
    res = UnitTestDesign.most_matches_existing(mc, mmm_case[1], mmm_case[2])
    @test res == mmm_case[3]
end

# Let's understand this most_matches_existing function by assuming,
# at this point, that it obeys the condition that the given parameter
# index must be nonzero. Let's always put that last. And let's just make
# one row to cover or not cover. Test for yes or no.
mm2_arity = [2, 2, 2, 2]
# test is a) a tuple to cover, b) an existing array. The param_idx is the last entry,
# so that's always zero for the "existing" array.
mm2_cases = [
    [[1; 1; 0; 1], [1, 1, 0, 0], 1],
    [[1; 1; 0; 1], [1, 0, 0, 0], 1],
    [[1; 1; 0; 1], [0, 0, 0, 0], 1],
    [[1; 1; 0; 1], [1, 0, 2, 0], 0],  # This could be a 1. This is subtle. Must match all existing.
    [[1; 1; 0; 1], [1, 1, 2, 0], 0]   # again, compatible but not exact.
]
for mm2_case in mm2_cases
    mc = UnitTestDesign.MatrixCoverage(reshape(mm2_case[1], length(mm2_arity), 1), 1, mm2_arity)
    res = UnitTestDesign.most_matches_existing(mc, mm2_case[2], length(mm2_arity))
    @test res == [mm2_case[3], 0]
end

## insert_tuple_into_tests
test_set = [
    1 1 1;
    1 0 3;
    1 3 2;
    1 2 2
]
allc = UnitTestDesign.MatrixCoverage(
    [
        1 2 0
        1 0 1
        0 0 3
        0 0 0
    ],
    3,
    [2, 3, 3]
)
UnitTestDesign.insert_tuple_into_tests(test_set, allc)

wider = [
    1 1 0; 1 2 0; 2 1 0;
    2 2 0; 3 1 0; 3 2 0
]
arity4 = [3, 2, 2]
allc4 = [
    1 0 1; 2 0 1; 3 0 1; 0 1 1; 0 2 1;
    1 0 2; 2 0 2; 3 0 2; 0 1 2; 0 2 2
]'
mc4 = UnitTestDesign.MatrixCoverage(collect(allc4), size(allc4, 2), arity4)
mm4 = UnitTestDesign.matches_from_missing(mc4, wider[1, :], 3)
@test mm4 == [2, 2]

### first_match_for_parameter(mc, param_idx)
arity = [3, 2, 2, 2]
mat = [0 1 2 0; 2 2 1 0]'
mcfm = UnitTestDesign.MatrixCoverage(collect(mat), size(mat, 2), arity)
fm1 = UnitTestDesign.first_match_for_parameter(mcfm, 3)
@test fm1 == [0, 1, 2, 0]
blank = UnitTestDesign.first_match_for_parameter(mcfm, 1)
@test blank == [2, 2, 1, 0]
lack = UnitTestDesign.first_match_for_parameter(mcfm, 4)
@test lack == [0, 0, 0, 0]

### fill_consistent_matches
arity = [2, 4, 4]
mat = [
    0 0 0;
    1 0 0;
    0 2 0;
    0 0 3;
    0 0 4
]'
mc = UnitTestDesign.MatrixCoverage(collect(mat), size(mat, 2), arity)
res1 = UnitTestDesign.fill_consistent_matches(mc, [0, 2, 0])
@test res1 == [1, 2, 3]

arity = [2, 4, 4]
mat2 = [
    0 0 0;
    1 2 0;
    0 2 0;
    1 0 4;
    0 0 3
]'
mc2 = UnitTestDesign.MatrixCoverage(collect(mat2), size(mat2, 2), arity)
res2 = UnitTestDesign.fill_consistent_matches(mc2, [0, 2, 0])
@test res2 == [1, 2, 4]

arity = [2, 4, 4]
mat3 = [
    0 0 0;
    1 2 0;
    0 2 0;
    0 0 3;
    1 0 4
]'
mc3 = UnitTestDesign.MatrixCoverage(collect(mat3), size(mat3, 2), arity)
res3 = UnitTestDesign.fill_consistent_matches(mc3, [0, 2, 0])
@test res3 == [1, 2, 3]

### add_coverage!(allc, row_cnt, entry)
# This checks that the rows of the matrix are reordered
# when a new test covers values. It expects a particular
# reordering due to swapping. Very much an implementation test.
ac_cases = [
    [[1 1 0; 1 2 0; 1 0 1; 1 0 2]', 4, 2, [1, 1, 1], 2, [1 0 2; 1 2 0; 1 1 0; 1 0 1]'],
    [[1 1 0; 1 2 0; 1 0 1; 1 0 2]', 3, 2, [1, 1, 1], 1, [1 2 0; 1 1 0; 1 0 1; 1 0 2]'],
    [[1 0 2; 0 2 2; 0 1 2]', 1, 2, [1, 0, 2], 0, [1 0 2; 0 2 2; 0 1 2]'],
    [[1 0 2; 0 2 2; 0 1 2]', 1, 2, [1, 0, 1], 1, [1 0 2; 0 2 2; 0 1 2]']
]
for ac_case in ac_cases
    mcac = UnitTestDesign.MatrixCoverage(collect(ac_case[1]), ac_case[2], [2, 3, 2])
    UnitTestDesign.add_coverage!(mcac, ac_case[4])
    @test mcac.remain == ac_case[5]
    @test mcac.allc == ac_case[6]
end

### match_score
ms_cases = [
    [[1 1 0; 1 2 0; 0 1 3]', 3, 2, [4,4,4], 0],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3]', 4, 2, [1,1,1], 1],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3]', 4, 2, [1,2,1], 2],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3]', 4, 2, [1,1,3], 2],
    [[1 1 0; 1 2 0; 0 2 1; 0 1 3]', 3, 2, [1,1,3], 1]
]
for ms_case in ms_cases
    mc_ms = UnitTestDesign.MatrixCoverage(collect(ms_case[1]), ms_case[2], [4, 4, 4])
    cnt = UnitTestDesign.match_score(mc_ms, ms_case[4])
    @test cnt == ms_case[5]
end

## remove_combinations
rc_cases = [
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2]', (x -> x[2] == 1 && x[3] == 1), [1 1 2; 1 2 1; 1 2 2]'],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2]', (x -> x[2] == 1), [1 2 1; 1 2 2]'],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2]', (x -> x[1] == 1 && x[3] == 1), [1 1 2; 1 2 2]'],
    [[1 1 1; 1 1 2; 1 2 1; 1 2 2]', (x -> x[2] < x[3]), [1 1 1; 1 2 1; 1 2 2]']
]
for rc_case in rc_cases
    mc_rc = UnitTestDesign.MatrixCoverage(collect(rc_case[1]), 4, [2, 3, 2])
    UnitTestDesign.remove_combinations!(mc_rc, rc_case[2])
    @test mc_rc.allc == rc_case[3]
end

## multi_way_coverage
mwc = UnitTestDesign.multi_way_coverage([2,3,4,2,2,3], Dict(3 => [[1,3,4,5]]), 2)
@test maximum(sum(mwc .!= 0, dims = 1)) == 3


text:
# issue 505
@testset "mdd+indent" begin
    h = """
        @def x = Dict(
            :a => 5,
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "5", h)
    h = """
        @def x = Dict(
            :a => (1, 2, 3),
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a][1]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "1", h)
    h = """
        @def x = Dict(
            :a => (1,
                   2,
                   3),
            :b => 7
            )
        A
        ```julia:ex
        locvar(:x)[:a][2]
        ```
        \\show{ex}
        """ |> fd2html
    @test occursin(
        "2", h)
end

# Blocks of definitions
@testset "mddefblock" begin
    s = """
        +++
        a = 5
        b = "hello"
        +++
        {{b}} {{a}}
        """ |> fd2html
    # danger of the stuff
    @test isapproxstr(s, "
hello 5
")
    s = """
        +++
        a = 5
        b = "hello"
        +++
        {{b}} {{a}}
        +++
        a = 7
        +++
        {{a}}
        """ |> fd2html
    @test isapproxstr(s, "
hello 7 7
")
    # more things
    s = """
       +++
       foo(x, y, z) = x^2 + 2y + z
       out = foo(2,3,4)
       +++
       {{out}}
       """ |> fd2html
    @test isapproxstr(s, "
14
")
    # errors
    s = """
        +++
        s = sqrt(-1)
        +++
        {{s}}
        """
    @test_throws ErrorException s |> fd2html
end

# blocks of definition with date
@testset "mddefblock+date" begin
    s = """
        +++
        a = 5
        pubdate = Date(2013, 9, 4)
        +++
        {{a}} {{pubdate}}
        """ |> fd2html
    @test isapproxstr(s, "
 5 2013-09-04
")
end

# windows return
@testset "mddefwindows" begin
    s = "+++\r\na = \"hello\"\r\n+++\r\n# hi\r\n{{a}}" |> fd2html
    @test isapproxstr(s,
        """

hi

        hello
        """
    )
end


text:
using Distributed

@everywhere include("../RandomWalk.jl")

numtraj = 1000000
time = vcat(0, [2^(i/3) for i=-15:50])
# Nvals = [2, 4, 6, 10, 12, 15, 20, 24, 30, 40, 50, 60, 80, 100]
Nvals = [4, 10, 20]
biases = [-2.0, 2.0, -1.5, 1.5, -1.0, 1.0,-0.6, 0.6]

for N in Nvals, b in biases
    println("Begin N = $N, b = $b")
    path = "/data/niels/sim/rw/biased/N_$(N)_b_$b"
    mkpath(path)
    traj = RandomWalk.biased_random_walk(numtraj, 2*N, 300000, time, bias=b)
    RandomWalk.savedata("$path/traj.h5", time, traj, 2*N)
end

max_stars_repo_name: nielsmde/RandomWalk

max_issues_repo_name: nielsmde/RandomWalk

max_forks_repo_name: nielsmde/RandomWalk


text:
using MCMCBenchmarks, Test, Random

@testset "SDT Tests " begin
    Random.seed!(6505)
    path = pathof(MCMCBenchmarks)
    include(joinpath(path,"../../Models/SDT/SDT_Functions.jl"))
    include(joinpath(path,"../../Models/SDT/SDT.jl"))
    Nd = 10^4
    d = 2.0
    c = 0.0
    Nreps = 1
    ProjDir = @__DIR__
    cd(ProjDir)
    samplers=(CmdStanNUTS(CmdStanConfig, ProjDir),
        AHMCNUTS(AHMC_SDT, AHMCconfig),
        DHMCNUTS(sampleDHMC)
        )
    options = (Nsamples=2000, Nadapt=1000, delta=.8, Nd=Nd)
    results = benchmark(samplers, simulateSDT, Nreps;options...)
    @test results[!,:d_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ d atol = .05
    @test results[!,:d_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ d atol = .05
    @test results[!,:d_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ d atol = .05
    @test results[!,:c_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ c atol = .05
    @test results[!,:c_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ c atol = .05
    @test results[!,:c_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ c atol = .05
    @test results[!,:d_sampler_rhat][1] ≈ 1 atol = .03
    @test results[!,:c_sampler_rhat][1] ≈ 1 atol = .03
    isdir("tmp") && rm("tmp", recursive=true)
end


text:
module Feature_Test

using SpikeSorting

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

d=DetectNeg()
cluster=ClusterOSort();
align=AlignMax();
reduce=ReductionNone();
thres=ThresholdMeanN();
num_channels=1;

(buf,nums)=output_buffer(num_channels);

function make_voltage(num_channels)

    v=rand(1:1000, 1000, num_channels);
    count = 1
    for i=525:550
        v[i,1] = v[i,1] + count
        count = count + 10
    end
    for i=551:775
        v[i,1]+=1
        count -= 10
    end
    v
end

v=make_voltage(num_channels)
#=
Time
=#

f=FeatureTime()
s1=create_multi(d,cluster,align,f,reduce,thres,num_channels);

cal!(s1,v,buf,nums,0)
cal!(s1,v,buf,nums,1)
cal!(s1,v,buf,nums,2)
onlinesort!(s1,v,buf,nums)

#=
Discrete Derivatives
=#
#=
f=FeatureDD()
s1=create_multi(d,cluster,align,f,reduce,thres,num_channels);

cal!(s1,v,buf,nums,0)
cal!(s1,v,buf,nums,1)
cal!(s1,v,buf,nums,2)
onlinesort!(s1,v,buf,nums)

=#
#=
Curvature
=#
#=
f=FeatureCurv()
s1=create_multi(d,cluster,align,f,reduce,thres,num_channels);

cal!(s1,v,buf,nums,0)
cal!(s1,v,buf,nums,1)
cal!(s1,v,buf,nums,2)
onlinesort!(s1,v,buf,nums)

facts() do

end
=#
end


text:
struct Ket <: AbstractState
    data::Vector{Scalar}
end

struct Bra <: AbstractState
    data::Vector{Scalar}
end

adjoint(ψ::Ket) = Bra(conj(data(ψ)))
adjoint(ϕ::Bra) = Ket(conj(data(ϕ)))

*(ϕ::Bra, ψ::Ket) = check_dimensions_match(ϕ, ψ) && data(ϕ) ⋅ data(ψ)
*(ψ::Ket, ϕ::Bra) = check_dimensions_match(ψ, ϕ) && Operator(data(ψ) * data(ϕ)')


text:
@node MvNormalMeanCovariance Stochastic [ out, (μ, aliases = [ mean ]), (Σ, aliases = [ cov ]) ]

conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :out } }) = MvNormalMeanCovariance
conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :μ } })   = MvNormalMeanCovariance
conjugate_type(::Type{ <: MvNormalMeanCovariance }, ::Type{ Val{ :Σ } })   = InverseWishart

# default method for mean-field assumption
@average_energy MvNormalMeanCovariance (q_out::Any, q_μ::Any, q_Σ::Any) = begin
    # naive: (ndims(q_out) * log2π + mean(logdet, q_Σ) + tr(cholinv(mean(q_Σ))*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)'))) / 2
    dim = ndims(q_out)

    m_mean, v_mean = mean_cov(q_μ)
    m_out, v_out   = mean_cov(q_out)
    inv_m_Σ        = mean(inv, q_Σ)

    result = zero(promote_type(eltype(q_out), eltype(q_μ), eltype(q_Σ)))
    result += mean(logdet, q_Σ)
    result += dim * log2π
    @inbounds for k1 in 1:dim, k2 in 1:dim   # optimize trace operation (indices can be interchanges because of symmetry)
        result += inv_m_Σ[k1,k2] * (v_out[k1,k2] + v_mean[k1,k2] + (m_out[k2]-m_mean[k2]) * (m_out[k1] - m_mean[k1]))
    end
    result /= 2

    return result
end

# default method for structured mean-field assumption
@average_energy MvNormalMeanCovariance (q_out_μ::Any, q_Σ::Any) = begin
    # naive: @views (d*log2π + mean(logdet, q_Σ) + tr(cholinv(mean(q_Σ))*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ))) / 2
    dim = div(ndims(q_out_μ), 2)

    m, V = mean_cov(q_out_μ)
    inv_m_Σ = mean(inv, q_Σ)

    result = zero(promote_type(eltype(q_out_μ), eltype(q_Σ)))
    result += mean(logdet, q_Σ)
    result += dim * log2π
    @inbounds for k1 in 1:dim, k2 in 1:dim   # optimize trace operation (indices can be interchanges because of symmetry)
        result += inv_m_Σ[k1,k2] * (V[k1,k2] + V[dim+k1,dim+k2] - V[dim+k1,k2] - V[k1,dim+k2] + (m[k1] - m[dim+k1]) * (m[k2] - m[dim+k2]))
    end
    result /= 2

    return result
end


text:
# WITH clause.

mutable struct WithClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    recursive::Bool
    list::Vector{SQLClause}

    WithClause(;
               over = nothing,
               recursive = false,
               list) =
        new(over, recursive, list)
end

WithClause(list...; over = nothing, recursive = false) =
    WithClause(over = over, recursive = recursive, list = SQLClause[list...])

"""
    WITH(; over = nothing, recursive = false, list)
    WITH(list...; over = nothing, recursive = false)

A `WITH` clause.

# Examples

```jldoctest
julia> c = FROM(:person) |>
           WHERE(OP("IN", :person_id,
                          FROM(:essential_hypertension) |>
                          SELECT(:person_id))) |>
           SELECT(:person_id, :year_of_birth) |>
           WITH(FROM(:condition_occurrence) |>
                WHERE(OP("=", :condition_concept_id, 320128)) |>
                SELECT(:person_id) |>
                CTE(:essential_hypertension));

julia> print(render(c))
WITH "essential_hypertension" AS (
  SELECT "person_id"
  FROM "condition_occurrence"
  WHERE ("condition_concept_id" = 320128)
)
SELECT "person_id", "year_of_birth"
FROM "person"
WHERE ("person_id" IN (
  SELECT "person_id"
  FROM "essential_hypertension"
))
```

```jldoctest
julia> c = FROM(:essential_hypertension) |>
           SELECT(OP("*")) |>
           WITH(recursive = true,
                FROM(:concept) |>
                WHERE(OP("=", :concept_id, 320128)) |>
                SELECT(:concept_id, :concept_name) |>
                UNION(all = true,
                      FROM(:eh => :essential_hypertension) |>
                      JOIN(:cr => :concept_relationship,
                           OP("=", (:eh, :concept_id), (:cr, :concept_id_1))) |>
                      JOIN(:c => :concept,
                           OP("=", (:cr, :concept_id_2), (:c, :concept_id))) |>
                      WHERE(OP("=", (:cr, :relationship_id), "Subsumes")) |>
                      SELECT((:c, :concept_id), (:c, :concept_name))) |>
                CTE(:essential_hypertension, columns = [:concept_id, :concept_name]));

julia> print(render(c))
WITH RECURSIVE "essential_hypertension" ("concept_id", "concept_name") AS (
  SELECT "concept_id", "concept_name"
  FROM "concept"
  WHERE ("concept_id" = 320128)
  UNION ALL
  SELECT "c"."concept_id", "c"."concept_name"
  FROM "essential_hypertension" AS "eh"
  JOIN "concept_relationship" AS "cr" ON ("eh"."concept_id" = "cr"."concept_id_1")
  JOIN "concept" AS "c" ON ("cr"."concept_id_2" = "c"."concept_id")
  WHERE ("cr"."relationship_id" = 'Subsumes')
)
SELECT *
FROM "essential_hypertension"
```
"""
WITH(args...; kws...) =
    WithClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(WITH), pats::Vector{Any}) =
    dissect(scr, WithClause, pats)

function PrettyPrinting.quoteof(c::WithClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(WITH))
    if c.recursive !== false
        push!(ex.args, Expr(:kw, :recursive, c.recursive))
    end
    if isempty(c.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(c.list, ctx))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::WithClause, c′) =
    WithClause(over = rebase(c.over, c′), list = c.list, recursive = c.recursive)


text:
function EulerSDE(traceT,drift,D,Y0)
	b=sqrt(2*D);
    dts = diff(traceT);
    steps = length(dts);
    sigma=sqrt.(dts);  # variance=sigma^2
    dW = sigma.*randn((steps,)) .+ 0

    traceY = fill(NaN,size(traceT))
    traceY[1] = Y0;
    y = Y0;
    for i = 2:length(traceT)
        y = y + drift(y)*dts[i-1] + F_ext*dts[i-1] + b*dW[i-1]
        traceY[i] = y;
    end
    return traceY
end


text:
using Test
using ClimaCore.Geometry, ClimaCore.DataLayouts
using LinearAlgebra, StaticArrays

@testset "AxisTensors" begin
    x = Geometry.Covariant12Vector(1.0, 2.0)
    y = Geometry.Contravariant12Vector(1.0, 4.0)

    M = Geometry.Axis2Tensor(
        (Geometry.Cartesian12Axis(), Geometry.Covariant12Axis()),
        [1.0 0.0; 0.5 2.0],
    )

    @test dot(x, y) == x' * y == 9.0
    @test dot(y, x) == y' * x == 9.0

    @test x == x
    @test x != Geometry.components(x)
    @test x != Geometry.Contravariant12Vector(Geometry.components(x)...)

    @test x[1] == 1.0
    @test y[2] == 4.0
    @test M[2] == 0.5
    @test M[2, 1] == 0.5
    @test M[:, 1] == Geometry.Cartesian12Vector(1.0, 0.5)
    @test M[1, :] == Geometry.Covariant12Vector(1.0, 0.0)

    @test x * y' ==
          x ⊗ y ==
          Geometry.AxisTensor(
              (axes(x, 1), axes(y, 1)),
              Geometry.components(x) * Geometry.components(y)',
          )

    @test Geometry.components(M * inv(M)) == @SMatrix [1.0 0.0; 0.0 1.0]
    @test Geometry.components(inv(M) * M) == @SMatrix [1.0 0.0; 0.0 1.0]

    @test x ⊗ 3 == Geometry.Covariant12Vector(3.0, 6.0)
    @test x ⊗ (1, (a = 2, b = 3)) == (
        Geometry.Covariant12Vector(1.0, 2.0),
        (
            a = Geometry.Covariant12Vector(2.0, 4.0),
            b = Geometry.Covariant12Vector(3.0, 6.0),
        ),
    )

    @test Geometry.components(M * inv(M)) == @SMatrix [1.0 0.0; 0.0 1.0]
    @test Geometry.components(inv(M) * M) == @SMatrix [1.0 0.0; 0.0 1.0]

    @test M * y == Geometry.Cartesian12Vector(1.0, 8.5)
    @test M \ Geometry.Cartesian12Vector(1.0, 8.5) == y

    @test_throws DimensionMismatch dot(x, x)
    @test_throws DimensionMismatch M * x
    @test_throws DimensionMismatch M \ x

    @test DataLayouts.basetype(typeof(x)) == Float64
    @test DataLayouts.typesize(Float64, typeof(x)) == 2
end

@testset "transform" begin
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)
    @test_throws InexactError Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0),
    )

    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test_throws InexactError Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0) * Geometry.Cartesian1Vector(1.0)',
    )
end