text:
#2 -echelon system with no cost/financial considerations

using Distributions
using InventoryManagement

#define network connectivity
net = MetaDiGraph(path_digraph(2)) # 1 (distribution center) -> 2 (market)
materials = [:A]
bom = [0]
set_prop!(net, :materials, materials)
set_prop!(net, :bill_of_materials, bom)

#specify parameters, holding costs and capacity, market demands and penalty for unfilfilled demand
set_props!(net, 1, Dict(:initial_inventory => Dict(:A => Inf)))

set_props!(net, 2, Dict(:initial_inventory => Dict(:A => 100),
                        :demand_distribution => Dict(:A => Normal(5,0.5)),
                        :demand_frequency => Dict(:A => 2),
                        :supplier_priority => Dict(:A => [1])))

#specify sales prices, transportation costs, lead time
set_props!(net, 1, 2, Dict(:lead_time => Dict(:A => Poisson(5))))

#create environment
num_periods = 100
env = SupplyChainEnv(net, num_periods, backlog = true)

#define reorder policy parameters
policy_type = :rQ #(s, S) policy
review_period = 25 #review every 25 periods
r = Dict((2,:A) => 20) #lower bound on inventory
Q = Dict((2,:A) => 80) #base stock level

#run simulation with reorder policy
simulate_policy!(env, r, Q; policy_type, review_period)

#make plots
using StatsPlots
#unfulfilled
fig1 = plot(env.demand.period, env.demand.unfulfilled, linetype=:steppost, lab="backlog",
                    xlabel="period", ylabel="level", title="Node 2, Material A")
#add inventory position
inv_on_hand = filter(i -> i.level < Inf, env.inv_on_hand)
plot!(inv_on_hand.period, inv_on_hand.level, linetype=:steppost, lab = "on-hand inventory")


text:
# This file is designed to test the visibility interpolation algorithm.
import PyPlot
import PyPlot.plt
using LaTeXStrings

using DiskJockey.constants
using DiskJockey.gridding
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gauss

# Because of the flipped nature of the sky (but not flipped nature of the UV plane)
# there are some tricky conventions about how to pack the array.

# Test to see if we get the convolutional interpolation correct by using a 2D
# Elliptical Gaussian.
# Realistic Gaussian will have scale dimensions (fatter in x direction)
const mu_RA = 2.0 # [arcsec]
const mu_DEC = -0.5 # [arcsec]
const s_x = 1.0 # [arcsec]
const s_y = 2.0 # [arcsec]
const rho = 0.0
const theta = 20.0 # deg
const p0 = [mu_RA, mu_DEC, s_x, s_y, rho] # [arcsec]
const p_center = [0.0, 0.0, s_x, s_y, rho] # [arcsec]

# If the synthesized beam is 0.3", to oversample this by a factor of 10 would require something like 512x512 pixels
nx = 128
ny = 128

# full span of the image
ra = fftspace(10., nx) # [arcsec]
dec = fftspace(10., ny) # [arcsec]

# convert ra and dec in [arcsec] to radians, and then take the sin to convert to ll, mm
ll = sin.(ra * arcsec) # direction cosines
mm = sin.(dec * arcsec)

# The natural, rotated and shifted Gaussian image
img = imageGauss(ll, mm, p0, 1, theta=theta)

# A centered Gaussian image which we will then rotate and shift via the phase theorem
img_plain1 = imageGauss(ll, mm, p_center, 1) #theta = 0

lam0 = lam0s["12CO2-1"] # [microns] 12CO 2-1

skim = SkyImage(img, ra, dec, lam0)

skim_plain1 = SkyImage(img_plain1, ra, dec, lam0)

# Apply a centered correction function on the rotated and offset image and Fourier transform the image
corrfun!(skim, 0., 0.)
shift_fft = transform(skim)

# Now, create images that we will use as test cases for doing the interpolation before or after
# sampling

# Apply a centered correction function on the centered image (to rotate and phase shift downsampled
# visibilities later)
corrfun!(skim_plain1, 0, 0)
vis_fft_center = transform(skim_plain1)
# Apply Rotation
# Apply Phase shift

# Return a normalized instance that is symmetric about 0
function scale(data)
    s = maximum(abs.(data))
    return norm = plt[:Normalize](vmin=-s, vmax=s, clip=false)
end

# Basic plot of the shifted image

# Because the sky convention is different than the way the SkyImage is stored,
# we need to flip the array for plotting
fig, ax = plt[:subplots](nrows=2, figsize=(5, 8))

ext = (skim.ra[end], skim.ra[1], skim.dec[1], skim.dec[end])
ax[1][:imshow](flipdim(skim.data[:,:,1], 2), interpolation="none", origin="lower", cmap=plt[:get_cmap]("Greys"), extent=ext)
ax[1][:contour](flipdim(skim.data[:,:,1], 2), origin="lower", extent=ext)
ax[1][:set_title]("Sky Projection")
ax[1][:set_xlabel](L"$\alpha$ [arcsec]")
ax[1][:set_ylabel](L"$\delta$ [arcsec]")

ext = (ll[1], ll[end], mm[1], mm[end])
ax[2][:imshow](skim.data[:,:,1], interpolation="none", origin="lower", cmap=plt[:get_cmap]("Greys"), extent=ext)
ax[2][:contour](skim.data[:,:,1], origin="lower", extent=ext)
ax[2][:set_title]("Raw Array")
ax[2][:set_xlabel](L"$ll$")
ax[2][:set_ylabel](L"$mm$")

fig[:subplots_adjust](left=0.15, right=0.85, hspace=0.25)
plt[:savefig]("gaussian_img_array_rot.png")

# Basic plot of the *centered* Gaussian image.
fig, ax = plt[:subplots](nrows=1, figsize=(5, 5))
# Real, analytic Gaussian
ext = (skim.ra[end], skim.ra[1], skim.dec[1], skim.dec[end])
aximg = ax[:imshow](flipdim(skim_plain1.data[:,:,1], 2), interpolation="none", origin="lower", cmap=plt[:get_cmap]("Greys"), extent=ext) #, norm = scale(img))
ax[:set_title]("image")
ax[:set_xlabel](L"$\alpha$ [arcsec]")
ax[:set_ylabel](L"$\delta$ [arcsec]")
#[left, bottom, width, height]
cax = fig[:add_axes]([0.84, 0.25, 0.03, 0.45])
cb = fig[:colorbar](aximg, cax=cax)

fig[:subplots_adjust](left=0.15, right=0.85, hspace=0.25)
plt[:savefig]("gaussian_img_center_rot.png")

# It's nice to see these plots and what they look like, but it doesn't really make sense
# to compare them, unless we are comparing *interpolated* visibilities to the analytic formula.
# Just comparing the image is WRONG.

# THIS (following) IS THE ONLY THING WE SHOULD CARE ABOUT.

# Next, choose some (u,v) points within the bounds and see how the
# interpolated values compare to what the FTGauss would return.

nu = length(shift_fft.uu)
zer = zeros(nu)

function plot_1d(analytic, approx, fname)
    fig, ax = plt[:subplots](nrows=2, figsize=(5, 8))
    ax[1][:plot](shift_fft.uu, zer, ".k", label="Grid spacing")
    ax[1][:plot](uu, real(approx), "ob", label="Interp")
    ax[1][:plot](uu, real(analytic), ".r", label="Analytic")
    # ax[1][:plot](vis_fft.uu, real(analytic_u), "or", label="Analytic")
    ax[1][:set_xlim](-100, 100)
    ax[1][:set_title]("Real")
    ax[1][:set_xlabel](L"u [k $\lambda$]")
    ax[1][:legend]()

    ax[2][:plot](shift_fft.uu, zer, ".k", label="Grid spacing")
    ax[2][:plot](uu, imag(approx), "ob", label="Interp")
    ax[2][:plot](uu, imag(analytic), ".r", label="Analytic")
    # ax[2][:plot](vis_fft.uu, imag(analytic_u), "or", label="Analytic")
    ax[2][:set_xlim](-100, 100)
    ax[2][:set_title]("Imag")
    ax[2][:set_xlabel](L"u [k $\lambda$]")

    plt[:savefig](fname)
end

n = 100
uu = linspace(-100, 100, n) # [kλ]
approx = Array{Complex128}(n)
analytic = Array{Complex128}(n)

# First, let's see how the interpolated points, with *no shift*, correspond to the analytic form
# with no shift.

for i=1:n
    u = uu[i]
    v = 0.0
    #println("Interpolating at $u, $v")
    approx[i] = interpolate_uv(u, v, vis_fft_center)
    analytic[i] = FTGauss(u, v, [0.0, 0.0, s_x, s_y, rho], 1)
end

plot_1d(analytic, approx, "interpolation_norot_noshift.png")

# Next, let's see how the interpolated points, then rotated, correspond to the analytic form with rotation.

for i=1:n
    u = uu[i]
    v = 0.0

    up = cosd(theta) * u - sind(theta) * v
    vp = sind(theta) * u - cosd(theta) * v

    approx[i] = interpolate_uv(up, vp, vis_fft_center)
    analytic[i] = FTGauss(u, v, [0.0, 0.0, s_x, s_y, rho], 1, theta=theta)
end

plot_1d(analytic, approx, "interpolation_rot_noshift.png")

# Finally, let's see how the interpolated points, rotated, then phase-shifted, compare to the analytic form with rotation and phase shift.
mu = Float64[mu_RA, mu_DEC] * arcsec # [radians]

for i=1:n
    u = uu[i]
    v = 0.0

    up = cosd(theta) * u - sind(theta) * v
    vp = sind(theta) * u - cosd(theta) * v

    # Interpolate points at rotated u', v' points, but then shift according to
    # direction and amplitude of mu offset using u,v coordinates.

    # Now shift the interpolated points according to the phase theorem
    R = Float64[u, v] * 1e3 #[λ]
    # Not actually in polar phase form
    shift = exp(-2pi * 1.0im * (R' * mu)[1])

    #println("Interpolating at $u, $v")
    approx[i] = shift * interpolate_uv(up, vp, vis_fft_center)
    analytic[i] = FTGauss(u, v, p0, 1, theta=theta)

end

plot_1d(analytic, approx, "interpolation_rot_shift.png")

ext: jl


text:
function digitalroot(n)
    ns = string(n)
    sum = 0
    for i in ns
      sum += (Int(i) - Int('0'))
    end
    if(sum > 9)
      return digitalroot(sum)
    else
      return sum
    end
 end

max_stars_repo_name: brunopasseti/CodeWarsKatas

max_issues_repo_name: brunopasseti/CodeWarsKatas

max_forks_repo_name: brunopasseti/CodeWarsKatas


text:
tab(n::Int) = " "^n
no_indent(io::IO) = IOContext(io, :indent=>0)
no_indent_first_line(io::IO) = IOContext(io, :no_indent_first_line=>true)

# 1.0 compatibility
function indent_print(io::IO, ::Nothing)
    indent = get(io, :indent, 0)
    print(io, tab(indent), "nothing")
end

function indent_print(io::IO, xs...)
    indent = get(io, :indent, 0)
    Base.print(io, tab(indent), xs...)
end

function indent_println(io::IO, xs...)
    if get(io, :no_indent_first_line, false)
        indent_print(no_indent(io), xs..., "\n")
    else
        indent_print(io, xs..., "\n")
    end
end

function within_indent(f, io)
    f(indent(io))
end

function indent(io)
    IOContext(io, :indent => get(io, :indent, 0) + 4)
end

function show_option(io::IO, m::MIME"text/plain", x)
    if !get(io, :no_indent_first_line, false)
        indent_print(io)
    end

    print(io, GREEN_FG(string(typeof(x))), "(;")
    println(io)

    show_option_fields(io, m, x)
    indent_print(io, ")")
end

function show_option(io::IO, m::MIME"text/html", x)
    topmost = get(io, :configurations_style_sheet, true)
    if topmost
        print(io, """

        <dev class="configurations-option-type">
        """)
    end

    print(io, """

$(summary(x))
    """)
    show_option_fields(IOContext(io, :configurations_style_sheet=>false), m, x)
    print(io, """

    """)
    topmost && print(io, "</dev>")
end

function show_option_fields(io::IO, m::MIME"text/plain", x)
    within_indent(io) do io
        for name in fieldnames(typeof(x))
            value = getfield(x, name)
            if value != field_default(typeof(x), name)
                indent_print(io, LIGHT_BLUE_FG(string(name)), " = ")
                show_option_value(no_indent_first_line(io), m, value)
                println(io, ",")
            end
        end
    end
end

function show_option_fields(io::IO, m::MIME"text/html", x)
    print(io, "
")
    for name in fieldnames(typeof(x))
        value = getfield(x, name)
        if value != field_default(typeof(x), name)
            print(io, "
")
            print(io, "", string(name), "=")
            show_option_value(io, m, value)
            print(io, "
")
        end
    end
    print(io, "
")
end

function show_option_value(io::IO, m::MIME, x)
    if is_option(x)
        show_option(io, m, x)
    else
        show(io, x)
    end
end

function show_option_value(io::IO, m::MIME"text/html", x::AbstractDict)
    print(io, "", typeof(x), "(")

    print(io, "
")
    for (k, v) in x
        print(io, "
")
        print(io, "&nbsp&nbsp", repr(k), "=>")
        print(io, "")
        show_option_value(io, m, v)
        print(io, "
")
    end
    print(io, ")
")
end

function show_option_value(io::IO, m::MIME, x::AbstractDict)
    if !get(io, :no_indent_first_line, false)
        indent_print(io)
    end

    print(io, typeof(x))
    print(io, "(")

    println(io)
    within_indent(io) do io
        for (k, v) in x
            indent_print(io, repr(k), " => ")
            show_option_value(no_indent_first_line(io), m, v)
            println(io, ",")
        end
    end
    indent_print(io, ")")
end

function show_option_value(io::IO, m::MIME, list::Vector)
    if !(any(is_option, list) || length(list) > 4)
        return show(io, list)
    end

    println(io, "[")
    within_indent(io) do io
        for each in list
            indent_print(io)
            show_option_value(io, m, each)
            if length(list) > 1
                println(io, ",")
            end
        end
    end
    indent_print(io, "]")
end

function show_option_value(io::IO, m::MIME"text/html", list::Vector)
    if !(any(is_option, list) || length(list) > 4)
        return print(io, string(list))
    end

    print(io, "[")
    print(io, "
")
    for (i, each) in enumerate(list)
        print(io, "
")
        print(io, "", i, ".", "")
        show_option_value(io, m, each)
        print(io, "
")
    end
    print(io, "]
")
end


text:
import java.io.ByteArrayInputStream;
import java.io.ByteArrayOutputStream;
import java.io.IOException;
import java.io.InputStream;
import java.net.URL;
import java.nio.file.Files;
import java.nio.file.OpenOption;
import java.nio.file.Path;
import java.nio.file.StandardOpenOption;
import java.nio.file.attribute.FileAttribute;
import java.nio.file.attribute.PosixFilePermissions;

public class URLCache {
    private Path temp;
    private static final FileAttribute[] defaultFileAttributes =
            { PosixFilePermissions.asFileAttribute(PosixFilePermissions.fromString("rwx------")) };
    private static final OpenOption[] defaultOpenOptions =
            { StandardOpenOption.TRUNCATE_EXISTING };
    private final int PRIVKEY(H){L} privkey;
    private final int PUBKEY(H){L} pubkey;

    public URLCache() {
	privkey = (PRIVKEY(H){L})42;
	pubkey = (PUBKEY(H){L})42;

        try {
            temp =
                    Files.createTempDirectory("browsercache",
                                              defaultFileAttributes);
        }
        catch (IOException e) {
            e.printStackTrace();
        }
    }

    public InputStream get(URL url) {
        int code = url.toString().hashCode();
        Path file = temp.resolve(Integer.toHexString(code));

        try {
            byte[] bytes = Files.readAllBytes(file);
            return new ByteArrayInputStream(CryptoLibrary.decrypt(privkey, bytes));
        }
        catch (IOException ioe) {
            return null;
        }
    }

    public void put(URL url, InputStream value) {
	byte[] bytes = CryptoLibrary.encrypt(pubkey, inputStreamToBytes(value));

	int code = url.toString().hashCode();
	Path file = temp.resolve(Integer.toHexString(code));
	Files.createFile(file, defaultFileAttributes);
	Files.write(file, bytes, defaultOpenOptions);
    }

    private static byte[] inputStreamToBytes(InputStream is) {
	ByteArrayOutputStream{L} byteArrayOutputStream =
				     new ByteArrayOutputStream();
	try {
	    int read;
	    read = is.read();
	    while (read != -1) {
		//	byteArrayOutputStream.write(read);
		try {
		read = is.read();
		} catch (Exception e) {}
	    }
	} catch (Exception e) {}

	byte[] retVal = byteArrayOutputStream.toByteArray();
	return retVal;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
include("macros.jl")

Ad = [
   [-1.0/6.0  3.0/6.0 -3.0/6.0 1.0/6.0];
   [ 3.0/6.0 -6.0/6.0  0.0/6.0 4.0/6.0];
   [-3.0/6.0  3.0/6.0  3.0/6.0 1.0/6.0];
   [ 1.0/6.0  0.0/6.0  0.0/6.0 0.0/6.0]
]

dAd = [
   [ 0.0 -0.5  1.0 -0.5];
   [ 0.0  1.5 -2.0  0.0];
   [ 0.0 -1.5  1.0  0.5];
   [ 0.0  0.5  0.0  0.0]
]

d2Ad = [
   [ 0.0 0.0 -1.0  1.0];
   [ 0.0 0.0  3.0 -2.0];
   [ 0.0 0.0 -3.0  1.0];
   [ 0.0 0.0  1.0  0.0]
]

function eval_UC_spline(smin, smax, orders, C, S)

    d = size(S,2)
    N = size(S,1)

    vals = zeros(N)
    eval_UC_spline!(smin, smax, orders, C, S, vals)
    return vals

end

function eval_UC_spline_G(a, b, orders, C, S)

    d = size(S,2)
    N = size(S,1)

    vals = zeros(N)
    grad = zeros(N,d)

    eval_UC_spline_G!(a, b, orders, C, S, vals, grad)

    return (vals, grad)

end

function eval_UC_multi_spline(a, b, orders, C, S)

    d = size(S,2)
    N = size(S,1)
    K = size(C,1) # number of splines to evaluate
    vals = zeros(K,N)
    eval_UC_multi_spline!(a, b, orders, C, S, vals)
    return vals

end

# problem with this approach: the functions don't get cached.

# fun = (create_function(1,"natural"))
# fun.args[2].args[2]
# fun.args[2].args[2]
#
# for d = 1:4
#     eval(create_function(d,"natural"))
# end
#
# for d = 1:4
#     eval(create_function_with_gradient(d,"natural"))
# end

# with this approach functions don't get cached either (yet)
# but it's cleaner

@generated function eval_UC_spline!(a, b, orders, C, S, V)
    d = C.parameters[2]
    # the speed penalty of extrapolating points when iterating over point
    # seems very small so this is the default
    fun = (create_function(d,"natural"))
    return fun.args[2].args[2]
end

@generated function eval_UC_spline_G!(a, b, orders, C, S, V, dV)
    d = C.parameters[2]
    fun = (create_function_with_gradient(d,"natural"))
    return fun.args[2].args[2]
end

@generated function eval_UC_multi_spline!(a, b, orders, C, S, V)
    d = C.parameters[2]-1 # first dimension of C indexes the splines
    # the speed penalty of extrapolating points when iterating over point
    # seems very small so this is the default
    fun = (create_function_multi_spline(d,"natural"))
    return fun.args[2].args[2]
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module TCASSimulatorImpl

export
    addObserver,

    simulate,
    TCASSimulator

using AbstractSimulatorImpl
using CommonInterfaces

using Base.Test

using Encounter
using PilotResponse
using DynamicModel
using WorldModel
using Sensor
using CollisionAvoidanceSystem
using RLESUtils, Observers

import CommonInterfaces.addObserver

type SimulationParameters

    em::AbstractEncounterModel
    pr::Vector{AbstractPilotResponse}
    dm::Vector{AbstractDynamicModel}
    wm::AbstractWorldModel
    sr::Vector{AbstractSensor}
    cas::Vector{AbstractCollisionAvoidanceSystem}

    number_of_aircraft::Int

    SimulationParameters() = new()
end

type TCASSimulator <: AbstractSimulator

    parameters::SimulationParameters

    time_step::Float64

    observer::Observer

    function TCASSimulator()

        obj = new()

        obj.parameters = SimulationParameters()
        obj.time_step = 1.

        obj.observer = Observer()

        return obj
    end
end

addObserver(sim::TCASSimulator, f::Function) = add_observer(sim.observer, f)
addObserver(sim::TCASSimulator, tag::AbstractString, f::Function) = add_observer(sim.observer, tag, f)

import Base.convert

convert(::Type{SimpleADMInitialState}, state::Union{CorrAEMInitialState, LLAEMInitialState}) = SimpleADMInitialState(state.t, state.x, state.y, state.h, state.v, state.psi, state.h_d)

convert(::Type{ASWMState}, state::SimpleADMOutputState) = ASWMState(state.t, state.x, state.y, state.h, state.vx, state.vy, state.vh)

convert(::Type{SimplePRCommand}, command::Union{CorrAEMCommand, LLAEMCommand}) = SimplePRCommand(command.t, command.v_d, command.h_d, command.psi_d, 0.0)

convert(::Type{SimpleADMCommand}, command::SimplePRCommand) = SimpleADMCommand(command.t, command.v_d, command.h_d, command.psi_d)

function convert(::Type{SimpleTCASSensorInput}, states::Vector{ASWMState})

    states_ = Array(SimpleTCASSRState, length(states))

    for i = 1:length(states)
        states_[i] = SimpleTCASSRState(states[i].t, states[i].x, states[i].y, states[i].h, states[i].vx, states[i].vy, states[i].vh)
    end

    return SimpleTCASSensorInput(states_)
end

convert(::Type{SimpleTCASInput}, output::SimpleTCASSensorOutput) = SimpleTCASInput(output.t, output.r, output.r_d, output.a, output.a_d, output.h, output.h_d)

convert(::Type{SimplePRResolutionAdvisory}, RA::SimpleTCASResolutionAdvisory) = SimplePRResolutionAdvisory(RA.h_d)
convert(::Type{SimplePRResolutionAdvisory}, RA::Void) = nothing

function simulate(sim::AbstractSimulator; bTCAS = false, sample_number = 0)

    aem = sim.parameters.em
    pr = sim.parameters.pr
    adm = sim.parameters.dm
    as = sim.parameters.wm
    sr = sim.parameters.sr
    cas = sim.parameters.cas

    number_of_aircraft = sim.parameters.number_of_aircraft
    @test number_of_aircraft == 2

    if sample_number == 0
        generateEncounter(aem)
    else
        generateEncounter(aem, sample_number = sample_number)
    end

    initial_coords = zeros(number_of_aircraft, 3)

    for i = 1:number_of_aircraft
        initial = Encounter.getInitialState(aem, i)
        initial_coords[i, :] = [initial.x, initial.y, initial.h]
    end

    # filter out the encounter within TCAS thresholds
#    for i = 1:number_of_aircraft
#        dmod, zthr = TCAS_thresholds(initial_coords[i, 3])
#
#        for j = (i + 1):number_of_aircraft
#            if norm(initial_coords[i, 1] - initial_coords[j, 1], initial_coords[i, 2] - initial_coords[j, 2]) <= dmod || abs(initial_coords[i, 3] - initial_coords[j, 3]) <= zthr
#                return -1
#            end
#        end
#    end

    for i = 1:number_of_aircraft
        initial = Encounter.getInitialState(aem, i)
        state = DynamicModel.initialize(adm[i], convert(SimpleADMInitialState, initial))
        WorldModel.initialize(as, i, convert(ASWMState, state))

        Sensor.initialize(sr[i])
        CollisionAvoidanceSystem.initialize(cas[i])
        PilotResponse.initialize(pr[i])
    end

    command = nothing

    while true
        for i = 1:number_of_aircraft
            command = Encounter.update(aem, i)

            if command == nothing
                break
            end

            if bTCAS
                states = WorldModel.getAll(as)
                output = Sensor.update(sr[i], convert(SimpleTCASSensorInput, states))
                RA = CollisionAvoidanceSystem.update(cas[i], convert(SimpleTCASInput, output))

                if RA != nothing
                    @notify_observer(sim.observer, "RA", [i, states[i].t, RA.h_d])
                end
            else
                RA = nothing
            end

            response = PilotResponse.update(pr[i], convert(SimplePRCommand, command), convert(SimplePRResolutionAdvisory, RA))
            state = DynamicModel.update(adm[i], convert(SimpleADMCommand, response))
            WorldModel.update(as, i, convert(ASWMState, state))
        end

        if command == nothing
            break
        end

        WorldModel.updateAll(as)
    end

    return 0
end

function TCAS_thresholds(h)

    if h < 1000
        dmod = 0
        zthr = 0
    elseif h < 2350
        dmod = 0.2 * 6076.12
        zthr = 600
    elseif h < 5000
        dmod = 0.35 * 6076.12
        zthr = 600
    elseif h < 10000
        dmod = 0.55 * 6076.12
        zthr = 600
    elseif h < 20000
        dmod = 0.80 * 6076.12
        zthr = 600
    elseif h < 42000
        dmod = 1.10 * 6076.12
        zthr = 700
    else
        dmod = 1.10 * 6076.12
        zthr = 800
    end

    return dmod, zthr
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
using StackOverflow

"""
    makequestionsarray(json_obj::Dict)

Returns a Vector of `Question`s extracted from a json object from stackexchange API.
"""
function makequestionsarray(json_obj::Dict)
    questionholder = Vector()
    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                link = ""
                view_count = 0
                creation_date = 0
                is_answered = false
                owner = Dict()
                last_activity_date = 0
                score = 0
                accepted_answer_id = 0
                question_id = 0
                title = 0
                answer_count = 0
                tags = []
                title = ""

                for (key, value) in item
                    if occursin("link", key)
                        link = value
                    elseif occursin("view_count", key)
                        view_count = value
                    elseif occursin("creation_date", key)
                        creation_date = value
                    elseif occursin("is_answered", key)
                        is_answered = value
                    elseif occursin("owner", key)
                        owner = value
                    elseif occursin("last_activity_date", key)
                        last_activity_date = value
                    elseif occursin("score", key)
                        score = value
                    elseif occursin("accepted_answer_id", key)
                        accepted_answer_id = value
                    elseif occursin("question_id", key)
                        question_id = value
                    elseif occursin("tags", key)
                        tags = value
                    elseif occursin("title", key)
                        title = value
                    elseif occursin("answer_count", key)
                        answer_count = value
                    end
                end

                question = Question(link, view_count, creation_date, is_answered, owner,
                        last_activity_date, score, accepted_answer_id, question_id, tags, title,
                        answer_count)
                push!(questionholder, question)

            end
        end
    end
    questionholder
end

"""
    getrecentquestionsfortag(;tag::String = "Julia", site::String = "stackoverflow", order::String = "desc",
                                    sort::String = "activity", todate::String = "", fromdate::String = "")

Will return the 30 StackOverflow questions given a certain tag and sort parameter.

Possible order values: "desc" or "asc"
Possible sort values: "activity", "votes", "creation", "hot", "week", or "month"
"""
function getrecentquestionsfortag(;tag::String = "Julia", site::String = "stackoverflow", order::String = "desc",
                                    sort::String = "activity", todate::String = "", fromdate::String = "")

    questionholder = Vector()

    if fromdate != "" && todate != ""
        datequery = "fromdate=$(fromdate)&todate=$(todate)&"  #?fromdate=1573776000&todate=1574467200&order
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions?$(datequery)order=$(order)&sort=$(sort)&tagged=$(tag)&site=$(site)")
    else
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions?order=$(order)&sort=$(sort)&tagged=$(tag)&site=$(site)")
    end

    json = convert_HTTP_Response_To_JSON(r)
    questionholder = makequestionsarray(json)
end

"""
    getquestionfromanswer(ans::Answer; sort::String = "activity", order::String = "desc")

Returns a `Vector` of `Question`s whose one of the answers is `ans` of type `Answer`.
"""
function getquestionfromanswer(ans::Answer; sort::String = "activity", order::String = "desc")
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/answers/$(ans.answer_id)/questions?order=$(order)&sort=$(sort)&site=stackoverflow")
    json = convert_HTTP_Response_To_JSON(r)
    questionholder = makequestionsarray(json)
end


text:
# Minimise energy functional to compute solutions
# Alex Tam, 12/10/2020

"Minimise energy functional to solve model"
function optimise_network(s_old, af, xl, mm, parN, parA, parM, Lxx, Lxy, Lyx, Lyy)
    dof = build_dof(s_old); # Construct DOF vector for initial guess
    @time result = optimize(x -> energy_functional(x, s_old, af, xl, mm, parN, parA, parM, Lxx, Lxy, Lyx, Lyy), dof, LBFGS(m = 20, alphaguess = LineSearches.InitialHagerZhang(α0 = 1.0)), Optim.Options(x_tol = parN.xTol, f_tol = parN.fTol, g_tol = parN.gTol); autodiff = :forward);
    if Optim.converged(result) == 0 # Check convergence
        @printf("Warning: Optimisation did not converge.\n")
    end
    return Optim.minimizer(result) # Return DOF vector that minimises energy
end


text:
using Codecs

function main()
  str_size = 10_000_000
  tries = 100

  str = repeat("a", str_size)
  str2 = ""

  print("encode: ")
  t = time()
  s = 0
  for i in range(0, tries)
    str2 = ASCIIString(encode(Base64, str))
    s += length(str2)
  end
  print(s, ", ", time() - t, "\n")

  print("decode: ")
  t = time()
  s = 0
  for i in range(0, tries)
    s += length(ASCIIString(decode(Base64, str2)))
  end
  print(s, ", ", time() - t, "\n")
end

main()

max_stars_repo_name: japgolly/kostya-benchmarks

max_issues_repo_name: japgolly/kostya-benchmarks

max_forks_repo_name: japgolly/kostya-benchmarks


text:
module CGI

import Base: EnvHash, AsyncStream, Pipe

# I'm exporting the internal helper functions in case anyone needs to use
# them in their own app for whatever reason.

export HttpFormPart, HttpFile, HttpInputCollection,
    CgiResponse, CgiApplication,
    header, CgiParseQuotedParameters, CgiParseHttpParameters,
    CgiPercentDecode, CgiParseSemicolonFields, CgiParseMultipartFormParts,
    CgiStartApplication

###
#   > HttpFormPart (type)
#
#       Contains data for a form part sent as multipart/form-data.
###

type HttpFormPart
    headers::Dict{UTF8String, Dict{UTF8String, UTF8String}}
    data::Array{Uint8}

    HttpFormPart() = new (Dict{UTF8String, Dict{UTF8String, UTF8String}}(), Uint8[])
end

###
#   > HttpFile (type)
#
#       Contains information about a file uploaded as multipart/form-data.
#
#       > name (field)
#           The name of the file as reported by the client's machine/browser.
#       > mime (field)
#           The MIME type of the file as reported by the client's machine/browser.
#       > data (field)
#           Binary file data.
###

type HttpFile
    name::UTF8String
    mime::UTF8String
    data::Array{Uint8}

    HttpFile() = new ("", "", Uint8[])
end

###
#   > HttpInputCollection (type)
#
#       Contains information about the current request/session.
#
#       > method (field)
#           The HTTP request method used to request the page (GET, POST, etc).
#       > get (field)
#           An array of string values specified in the URL query string (e.g. ?var=value).
#       > post (field)
#           An array of string values sent as form data via POST.
#       > files (field)
#           A string-indexed array of HttpFile objects, containing files uploaded in the request.
###

type HttpInputCollection
    method::UTF8String
    get::Dict{UTF8String, UTF8String} # A decision was made here to support only string keys and values
    post::Dict{UTF8String, UTF8String}
    files::Dict{UTF8String, HttpFile}

    HttpInputCollection() = new ("", Dict{UTF8String, UTF8String}(), Dict{UTF8String, UTF8String}(), Dict{UTF8String, HttpFile}())
end

###
#   > CgiResponse (type)
#
#       Data relevant to the HTTP response to be sent to the client.
#
#       > headers (field)
#           An array of HTTP headers (e.g. Content-Type).
###

type CgiResponse
    headers::Dict{UTF8String, UTF8String}

    CgiResponse() = new (Dict{UTF8String, UTF8String}())
end

###
#   > CgiApplication (type)
#
#       A container for an instance of a CGI application/session.
#
#       > input (field)
#           See HttpInputCollection.
#       > server (field)
#           An array of values defined by the environment variables.
#       > response (field)
#           The HTTP response to be sent (see CgiResponse).
#       > stdoutOriginal (field) [internal]
#           A copy of the handle to the original STDOUT stream.
#       > outRead (field) [internal]
#           The read end of the STDOUT pipe created by this module to capture (buffer) output.
#       > outWrite (field) [internal]
#           The write end of the STDOUT pipe created by this module to capture (buffer) output.
###

type CgiApplication
    input::HttpInputCollection
    server::Dict{UTF8String, UTF8String}
    response::CgiResponse

    ## Used internally

    stdoutOriginal::AsyncStream
    outRead::Pipe
    outWrite::Pipe

    ##

    CgiApplication() = new (HttpInputCollection(), Dict{UTF8String, UTF8String}(), CgiResponse())
end

###
#   > header (function)
#
#       Sets a HTTP header for the given CGI application.
###

function header(app::CgiApplication, field::UTF8String, value::UTF8String)
    app.response.headers[field] = value

    return app
end

function header(app::CgiApplication, field::UTF8String, value::Int)
    app.response.headers[field] = string(value)

    return app
end

function header(app::CgiApplication, field::String, value::String)
    app.response.headers[utf8(field)] = utf8(value)

    return app
end

function header(app::CgiApplication, field::String, value::Int)
    app.response.headers[utf8(field)] = string(value)

    return app
end

###
#   > CgiParseQuotedParameters (function)
#
#       An internal helper function to parse strings in "name=value" format.
###

function CgiParseQuotedParameters(data::UTF8String)
    tokens = split(data, "=", 2)

    if (length(tokens) == 2)
        return (tokens[1], tokens[2])
    end

    return null
end

###
#   > CgiParseHttpParameters (function)
#
#       An internal helper function to parse HTTP parameter strings.
###

function CgiParseHttpParameters(dataString::UTF8String)
    tokenArray::Dict{UTF8String, UTF8String} = Dict{UTF8String, UTF8String}()

    if (length(dataString) > 0)
        tokens = split(dataString, "&")

        for token in tokens
            tokenTokens = split(token, "=", 2)

            if (length(tokenTokens) == 2)
                tokenArray[tokenTokens[1]] = tokenTokens[2]
            else
                tokenArray[tokenTokens[1]] = ""
            end
        end
    end

    return tokenArray
end

###
#   > CgiPercentDecode (function)
#
#       An internal helper function to decode percent-encoded strings
#       into UTF-8 strings.
###

function CgiPercentDecode(dataString::UTF8String)
    finalString::UTF8String = ""

    segments = split(dataString, '%', 0, false)

    # If there's more than one segment, there must be a % somewhere that needs decoding.

    if (length(segments) > 1)
        for segment in segments
            # At least two characters are required after the % marker.

            if (length(segment) > 1)
                # Get the two characters and convert them into their hexadecimal/byte value.

                try
                    hex = hex2bytes(ascii(segment[1:2]))

                    finalString = string(finalString, utf8(hex))
                catch error
                    # This will happen if the characters following the % are not valid hexadecimal characters.
                end

                # Add the rest of the segment to the byte array.

                try
                    finalString = string(finalString, segment[3:end])
                catch error
                    # This will be a bounds error due to the segment being only 2 characters in length.
                    Base.showerror(STDOUT, error)
                end
            else
                # Add the segment to the byte array.

                finalString = string(finalString, segment)
            end
        end # for
    else
        finalString = string(finalString, segments[1])
    end # if

    return finalString
end

###
#   > CgiPercentDecodeUint8 (function)
#
#       An internal helper function to decode percent-encoded strings
#       into binary data (Uint8).
###

function CgiPercentDecodeUint8(dataString::UTF8String)
    data = Uint8[]

    segments = split(dataString, '%', 0, false)

    # If there's more than one segment, there must be a % somewhere that needs decoding.

    if (length(segments) > 1)
        for segment in segments
            # At least two characters are required after the % marker.

            if (length(segment) > 1)
                # Get the two characters and convert them into their hexadecimal/byte value.

                try
                    hex = hex2bytes(ascii(segment[1:2]))

                    push!(data, hex[1])
                catch error
                    # Doesn't matter.
                end

                # Add the rest of the segment to the byte array.

                try
                    append!(data, bytestring(segment[3:end]).data)
                catch error
                    # Doesn't matter.
                end
            else
                # Add the segment to the byte array.

                append!(data, bytestring(segment).data)
            end
        end # for
    else
        append!(data, bytestring(segments[1]).data)
    end # if

    return data
end

###
#   > CgiParseSemicolonFields (function)
#
#       An internal helper function to parse semicolon-separated data.
###

function CgiParseSemicolonFields(dataString::UTF8String)
    dataString = dataString * ";"

    data = Dict{UTF8String, UTF8String}()

    prevCharacter::Char = 0x00
    inSingleQuotes::Bool = false
    inDoubleQuotes::Bool = false
    ignore::Bool = false
    workingString::UTF8String = ""

    dataStringLength::Int = length(dataString)

    dataStringLengthLoop::Int = dataStringLength + 1

    charIndex::Int = 1

    while charIndex < dataStringLengthLoop
        character = dataString[charIndex]

        if (!inSingleQuotes && character == '"' && prevCharacter != '\\')
            inDoubleQuotes = !inDoubleQuotes
            ignore = true
        end

        if (!inDoubleQuotes && character == '\'' && prevCharacter != '\\')
            inSingleQuotes = !inSingleQuotes
            ignore = true
        end

        if (charIndex == dataStringLength || (character == ';' && !(inSingleQuotes || inDoubleQuotes)))
            workingString = strip(workingString)

            if (length(workingString) > 0)
                decoded = CgiParseQuotedParameters(workingString)

                if (decoded != null)
                    (key, value) = decoded

                    data[key] = value
                else
                    data[workingString] = workingString
                end

                workingString = ""
            end
        elseif (!ignore)
            workingString = workingString * string(character)
        end

        prevCharacter = character

        charIndex  = charIndex + 1

        ignore = false
    end

    return data
end

###
#   > CgiParseMultipartFormParts (function)
#
#       An internal helper function to process multipart form data.
###

function CgiParseMultipartFormParts(data::Array{Uint8}, boundary::UTF8String, boundaryLength::Int = length(boundary))
    formParts::Array{HttpFormPart} = HttpFormPart[]

    part::HttpFormPart = HttpFormPart()

    ### Go through each byte of data, parsing it into form parts, headers and data.

    # The loop is perhaps slightly ambitious, as I wanted to be able to parse all the data
    # in one pass - rather than one pass for boundaries, another for headers, etc.

    # According to the spec, the boundary chosen by the client must be a unique string
    # i.e. there should be no conflicts with the data within - so it should be safe to just do a basic string search.

    headerRaw::Array{Uint8} = Uint8[]

    captureAsData::Bool = false

    crOn::Bool = false
    hadLineEnding::Bool = false
    foundBoundary::Bool = false
    foundFinalBoundary::Bool = false

    bytes::Int = length(data)

    byteIndexOffset::Int = 0
    testIndex::Int = 1
    byteTestIndex::Int = 0

    byte::Uint8 = 0x00

    # Skip over the first boundary and CRLF

    byteIndex::Int = boundaryLength + 5

    while (!foundFinalBoundary && byteIndex <= bytes)
        byte = data[byteIndex]

        # Test for boundary.

        if (
            (byte == 0x0d && data[byteIndex + 1] == 0x0a && data[byteIndex + 2] == '-' && data[byteIndex + 3] == '-')
            || (byte == '-' && data[byteIndex + 1] == '-')
        )
            foundBoundary = true

            if (byte == 0x0d)
                byteIndexOffset = byteIndex + 3
            else
                byteIndexOffset = byteIndex + 1
            end

            byteTestIndex = byteIndexOffset

            testIndex = 1;

            while (testIndex < boundaryLength)
                byteTestIndex = byteIndexOffset + testIndex

                if (byteTestIndex > bytes || data[byteTestIndex] != boundary[testIndex])
                    foundBoundary = false
                    break
                end

                testIndex = testIndex + 1
            end

            if (foundBoundary)
                if (data[byteTestIndex + 2] == '-')
                    foundFinalBoundary = true
                    byteIndex = byteTestIndex + 5
                else
                    byteIndex = byteTestIndex + 3
                end

            end
        else
            foundBoundary = false
        end

        ## Otherwise, process data

        if (foundBoundary)
            captureAsData = false
            crOn = false
            hadLineEnding = false
            foundBoundary = false

            push!(formParts, part)

            part = HttpFormPart()
        else
            if (captureAsData)
                push!(part.data, byte)
            else
                ## Check for CR

                if (byte == 0x0d)
                    crOn = true
                else
                    ## Check for LF and previous CR

                    if (byte == 0x0a && crOn)
                        ## Check for CRLFCRLF

                        if (hadLineEnding)
                            ## End of headers

                            captureAsData = true

                            hadLineEnding = false
                        else
                            ## End of single-line header

                            header::UTF8String = utf8(headerRaw)

                            if (length(header) > 0)
                                headerParts = split(header, ": ", 2)

                                valueDecoded = CgiParseSemicolonFields(utf8(headerParts[2]));

                                if (length(valueDecoded) > 0)
                                    part.headers[headerParts[1]] = valueDecoded
                                end
                            end

                            headerRaw = Uint8[]

                            hadLineEnding = true
                        end
                    else
                        if (hadLineEnding)
                            hadLineEnding = false
                        end

                        push!(headerRaw, byte)
                    end

                    crOn = false
                end
            end
        end

        byteIndex = byteIndex + 1
    end

    # ...

    return formParts
end

###
#   > CgiStartApplication (function)
#
#       Configure and populate the given CGI application instance.
###

function CgiStartApplication(app::CgiApplication, environment::EnvHash, inputStream::AsyncStream)
    # Save environment variables.
    # Dev note: This loop may not be necessary to ensure that the strings are all UTF-8.

    for (variableName::UTF8String, variableValue::UTF8String) in environment
        app.server[variableName] = variableValue
    end

    # Save the request method

    app.input.method = get!(app.server, "REQUEST_METHOD", "GET")

    ## Parse query string

    app.input.get = CgiParseHttpParameters(app.server["QUERY_STRING"])

    ### Handle POST input (if applicable)

    if (
        app.input.method == "POST"
        && (postDataLength::Int32 = parseint(Int32, get!(app.server, "CONTENT_LENGTH", "0"))) > 0
    )
        # Read data as bytes from the socket.

        postDataBytes::Array{Uint8} = readbytes(inputStream, postDataLength)

        ### Process x-www-form-urlencoded POST data (if applicable)

        if (get!(app.server, "CONTENT_TYPE", "") == "application/x-www-form-urlencoded")
            # Convert POST data to UTF-8 string.

            postDataString::UTF8String = utf8(postDataBytes)

            # Split post data into tokens based on "&" and "=" separators - and then
            # iterate over each token, decoding any percent-encoded data.

            app.input.post = CgiParseHttpParameters(postDataString)

            if (length(app.input.post) > 0)
                for (name, dataString) in app.input.post
                    app.input.post[name] = CgiPercentDecode(dataString)
                end
            end

        ### Process multipart/form-data POST data (e.g. file upload/s) (if applicable)

        elseif (searchindex(get!(app.server, "CONTENT_TYPE", ""), "multipart/form-data") != 0)
            boundary::UTF8String = app.server["CONTENT_TYPE"][(searchindex(app.server["CONTENT_TYPE"], "boundary=") + 9):end]

            boundaryLength::Int = length(boundary)

            if (boundaryLength > 0)
                formParts::Array{HttpFormPart} = CgiParseMultipartFormParts(postDataBytes, boundary, boundaryLength)

                ### Process form parts

                if (length(formParts) > 0)
                    for part::HttpFormPart in formParts
                        hasFile::Bool = false
                        file::HttpFile = HttpFile()
                        fileFieldName::UTF8String = ""

                        for (field::UTF8String, values::Dict{UTF8String, UTF8String}) in part.headers
                            if (
                                field == "Content-Disposition"
                                && getkey(values, "form-data", null) != null
                            )

                                # Check to see whether this part is a file upload
                                # Otherwise, treat as basic POST data

                                if (getkey(values, "filename", null) != null)
                                    if (length(values["filename"]) > 0)
                                        fileFieldName = values["name"]
                                        file.name = values["filename"]
                                        hasFile = true
                                    end
                                elseif (getkey(values, "name", null) != null)
                                    app.input.post[values["name"]] = utf8(part.data)
                                end

                                ## ...
                            elseif (field == "Content-Type")
                                (file.mime, mime) = first(values)
                            end
                        end # for

                        if (hasFile)
                            file.data = part.data

                            app.input.files[fileFieldName] = file

                            fileFieldName = ""
                            file = HttpFile()
                            hasFile = false
                        end # if
                    end # for
                end # if
            end  # if

        ### ...

        end # if/elseif
    end # End POST

    # ...

    return app
end

###
#   > newApplication (function)
#
#       Create an instance of a CGI application.
###

function newApplication(environment::EnvHash = ENV, inputStream::AsyncStream = STDIN)
    # Create instance.

    cgi = CgiApplication()

    # Redirect STDOUT.

    cgi.stdoutOriginal = STDOUT

    (cgi.outRead, cgi.outWrite) = redirect_stdout()

    # Redirect STDERR to STDOUT.

    redirect_stderr(cgi.outWrite)

    ### Standard shutdown routine: write headers and flush output.

    atexit(function ()
        # A content type is required by the CGI spec; if none is set, we'll default to text/html.

        if (getkey(cgi.response.headers, "Content-Type", null) == null)
            cgi.response.headers["Content-Type"] = "text/html"
        end

        # Get buffered output.

        close(cgi.outWrite)

        start_reading(cgi.outRead)
        yield()

        bytesAvailable = nb_available(cgi.outRead)

        if bytesAvailable > 0
            data = readavailable(cgi.outRead)
        else
            data = ""
        end

        close(cgi.outRead)

        # Restore original STDOUT.

        redirect_stdout(cgi.stdoutOriginal)

        # Add Content-Length header.

        header(cgi, "Content-Length", bytesAvailable)

        # Output HTTP headers.

        for (header, value) in cgi.response.headers
            print(header * ": " * value * "\r\n")
        end

        print("\r\n")

        # Write output

        write(STDOUT, data)
    end)

    # Boot the default instance

    CgiStartApplication(cgi, environment, inputStream)

    # Return the CGI application instance.

    return cgi
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
module TEST_SAMSARA
using Samsara, Test, Conception
using Printf

""" global state_space - for testing """
state_space = MuExS()
the_active_state = Conception.the_active_event_of(state_space)

""" Setting up the test system: [Ω | a b c d e] , where | represents the gate.
Activating [c] causes the gate to open.
(e.g. going to SAT [Ω] from SAT [a] have to go via [c] (to open the gate ||) )
"""
function create_Ω_case()
    Conception.__purge_everything!!()

    global state_space = MuExS()
    first_link = LinkedCardinalNode(:Ω, in_MuEx=state_space)
    linear_list = Samsara.linked_list_factory(5, in_MuExS=state_space)
    # unlocking the gate by visiting [c]
    the_key = linear_list[3]
    Conception.set_connected_concept(the_key, SAT(:gate_open))
    # Connect the gate:
    the_gate = LinkedGate(first_link, conditional=SAT(:gate_open))
    Samsara._set_node_to_E!(first_link, the_gate)
    Samsara.gate_set_connected_node!(the_gate, linear_list[1])
    Samsara._set_node_to_W!(linear_list[1], the_gate)
    # Print
    println("Creating Ω test case: ", state_space)
    activate!(linear_list[1])  #activate element [a]    ((( eller er dette [e]? )))
    return state_space, the_gate, the_key
end

""" Convenience-function for testing the Ω case """
function traverse_MuExS(the_muex::Conception.MuExS; direction =:East)
    the_item = Conception.the_active_event_of(the_muex)
    number_of_nodes = 0
    while !isnothing(the_item)
        number_of_nodes+=1
        last_item = the_item
        direction == :West && activate!(Samsara.west_of(the_item))
        direction == :East && activate!(Samsara.east_of(the_item))
        if last_item == (the_item = Conception.the_active_event_of(the_muex))
            break
        end
    end
    return number_of_nodes
end

function go_west()
    current_state = Conception.the_active_event_of(state_space)
    activate!(Samsara.west_of(current_state))
end
function go_east()
    current_state = Conception.the_active_event_of(state_space)
    activate!(Samsara.east_of(current_state))
end

# # # # # # # # # TESTS # # # # # # # # # #
@testset "system initiation" begin
    DemoSys = Samsara.DemoSys

    case = DemoSys()
    @test isa(case, Samsara.AbstractSystem)

    @test isa(system_state(case), Tuple)

    case = DemoSys((1), 1)
    @test system_state(case) == (1)
    @test dimensionality(case) == 1
    " Observable state is defined by argument: In DemoSys, this is of type Tuple "

    case = DemoSys((1,1), 3)
    @test system_state(case) == (1, 1)
    @test dimensionality(case) == 2
    """
    Create a struct of type DemoSys, with observable_param = (1,1) and latent_param = 3:
        - system_state(case) = (observable_param)
        - dimentionality(case) = length(observable_param)
    """

    step_system_mechanics!(case) # update
    @test system_state(case) == (4, 4)
    step_system_mechanics!(case) # update
    @test system_state(case) == (6, 6)
    """
    DemoSys have defined system mechanics in step_system_mechanics!(DemoSys)
        --> other systems can be created by inheriting AbstractSystem and by
        multiple dispatch of function step_system_mechanics!
    """

    function Samsara.step_system_mechanics!(arg::DemoSys)
        arg._observable_parameters = (0)
    end
    Samsara.step_system_mechanics!(case) # update
    @test dimensionality(case) == 1
    " It is possible to redefine system mechanics, although this is discouraged.. "

    case = DemoSys( (10.0, 10.0), (0.0, 0.0) )
    function Samsara.step_system_mechanics!(arg::DemoSys)
        arg._observable_parameters = arg._observable_parameters .+ arg._latent_variables
    end
    Samsara.step_system_mechanics!(case) # update
    @test system_state(case) == (10.0, 10.0)
    case._latent_variables = (1.0, -5.0)
    Samsara.step_system_mechanics!(case) # update
    @test system_state(case) == (11.0, 5.0)
    " Døme på korleis pådrag settes ved latent_variables: Ingen krav om at det er skalar heller! "

    """
        Vi kan lage simuleringer ved å definere system mechanics i funksjonen
        Samsara.step_system_mechanics!(sys) -- som vidare åpner for multiple dispatch ved å
        lage andre subtyper SystemX <: AbstractSystem, og bare definere
        Sams.step_system_mechanics!(SysX)

        Funksjonane kan vidare undersøke om event-SAT er aktiv:  BAM! gjennom:
        HAL.still_active(SAT)
    """
end#testset

@testset "Tests involving LinkedLists with LinkedGate" begin
    the_muex, the_gate, the_key = create_Ω_case()
    first_link = the_muex._elements[1]
    @test isa(first_link, SAT)

    activate!(first_link)
    Samsara.open_gate!(the_gate, false)
    number_of_nodes = traverse_MuExS(the_muex)
    @test number_of_nodes == 1
    """ Gate:CLOSED -- iteration through the goes through all elements BEFORE gate. """

    activate!(first_link)
    activate!(SAT(:gate_open))
    number_of_nodes = traverse_MuExS(the_muex)
    @test number_of_nodes == 6
    """ Gate:OPEN -- iteration through the goes through all elements gate. """

    activate!(first_link)
    deactivate!(SAT(:gate_open))
    number_of_nodes = traverse_MuExS(the_muex)
    @test number_of_nodes == 1

    activate!(the_key)
    @test is_active(the_key._node)
    @test is_active(the_key._node._connectedConcept)
    activate!(first_link)
    number_of_nodes = traverse_MuExS(the_muex)
    @test number_of_nodes == 6

    # We are at last link: gate closed.
    #   Traverse list east again, passing SAT(c) => linked with SAT(:gate_open) => opened gate!
    deactivate!(SAT(:gate_open)) # CLOSE GATE
    number_of_nodes = traverse_MuExS(the_muex, direction=:West)
    @test number_of_nodes == 6   # .. still, we get 6 : Explanation:
    @test is_active(SAT(:gate_open))
    #               #       # activating the_key, along the way, opens the gate before we reach it.
    """ Traversing back activates the key, opening the gate again before we reach it. """
end
""" Ω test works as planned: See create_Ω_case() in line 10
Ω-case:     [Ω | a b c d e] , where | represents the gate.
The gate is initially locked, but can be opened by visiting [c]
Getting to [Ω] from [a] requires that we visit [c] first.
"""

@testset "First attempt to integrate Desire with experience/HAL: iSAT with linked actions" begin
    using Samsara: the_node, east_of, west_of
    using Conception: the_active_event_of

    global state_space
    the_muex, the_gate, the_key = create_Ω_case()
    action_set = MuExS()
    #   HER STÅR EG FAST: TA DET OPP I MORGON!
    aE = iSAT(:East, linked_function=()->go_east(), inMuExS=action_set)
    aW = iSAT(:West, linked_function=()->go_west(), inMuExS=action_set)
    @test length(action_set) == 2

    pre_state = the_active_event_of(state_space)
    go_east()
    @test the_active_event_of(state_space) == the_node(east_of(pre_state))
    go_west()
    @test the_active_event_of(state_space) == pre_state
    " verifying the linked functions for the SAT in the action_set.. "

    pre_state = the_active_event_of(state_space)
    activate!(aE)
    @test the_active_event_of(state_space) == the_node(east_of(pre_state))
    activate!(aW)
    @test the_active_event_of(state_space) == pre_state         # back again..
    pre_state = Conception.the_active_event_of(state_space)
    @test Conception.the_active_event_of(state_space) == pre_state  #gate closed..
    """ We can go to the east by activating iSAT :aE (with linked function go_east(),
    working on global state space) and we can go west by activating iSAT :aW (go_west())
    When the gate is closed, we stay put..
    """

    while(!is_active(SAT(:Ω)))
        activate!(rand(action_set))
    end
    @assert is_active(SAT(:Ω))
    " Do random-walk untill we reach Ω - gate is opened and we have gotten through! "

    # TODO
        #- SÅ kan vi teste: resett heile env, ved å deaktivere alt og lukker porten; for så
        #HENT UT DESIRE! DET ER SYKT SPENNENDES!

end
    # PLAN
    # - kjør random walk. Test IV for SAT(a) mtp. Ω , når gate er åpen, når gate er stengd.
    #   -> Du kan tenke: SAT(:door_open) som satA, og satB, satA som events på veg mot satOmega:
    #       (blir nesten som Q-learning! IV-learning med state-SAT som actions!)

end#module TEST_SAMSARA

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
@doc doc"""
    Sphere{N} <: Manifold

The unit sphere manifold $\mathbb S^n$ represented by $n+1$-Tuples, i.e. in by
vectors in $\mathbb R^{n+1}$ of unit length

# Constructor

    Sphere(n)

generates the $\mathbb S^{n}\subset \mathbb R^{n+1}$
"""
struct Sphere{N} <: Manifold end
Sphere(n::Int) = Sphere{n}()

function representation_size(::Sphere{N}) where N
    return (N+1,)
end

@doc doc"""
    manifold_dimension(S::Sphere)

Return the dimension of the manifold $\mathbb S^n$, i.e. $n$.
"""
manifold_dimension(S::Sphere{N}) where {N} = N

project_point!(S::Sphere, x) = (x ./= norm(x))

project_tangent!(S::Sphere, w, x, v) = (w .= v .- dot(x, v) .* x)

@doc doc"""
    inner(S::Sphere, x, w, v)

compute the inner product of the two tangent vectors `w,v` from the tangent
plane at `x` on the sphere `S=`$\mathbb S^n$ using the restriction of the
metric from the embedding, i.e. $ (v,w)_x = v^\mathrm{T}w $.
"""
@inline inner(S::Sphere, x, w, v) = dot(w, v)

norm(S::Sphere, x, v) = norm(v)

distance(S::Sphere, x, y) = acos(clamp(dot(x, y), -1, 1))

function exp!(S::Sphere, y, x, v)
    θ = norm(S, x, v)
    y .= cos(θ) .* x .+ usinc(θ) .* v
    return y
end

function log!(S::Sphere, v, x, y)
    cosθ = dot(x, y)
    if cosθ ≈ -1
        fill!(v, 0)
        if x[1] ≈ 1
            v[2] = 1
        else
            v[1] = 1
        end
        copyto!(v, v .- dot(x, v) .* x)
        v .*= π / norm(v)
    else
        cosθ = cosθ > 1 ? one(cosθ) : cosθ
        θ = acos(cosθ)
        v .= (y .- cosθ .* x) ./ usinc(θ)
    end
    project_tangent!(S, v, x, v)
    return v
end

injectivity_radius(S::Sphere, args...) = π

function zero_tangent_vector!(S::Sphere, v, x)
    fill!(v, 0)
    return v
end

function vector_transport_to!(M::Sphere, vto, x, v, y, ::ParallelTransport)
    v_xy = log(M, x, y)
    vl = norm(M, x, v_xy)
    vto .= v
    if vl > 0
        factor = 2*dot(v, y)/(norm(x + y)^2)
        vto .-= factor.*(x .+ y)
    end
    return vto
end

function flat!(M::Sphere, v::FVector{CotangentSpaceType}, x, w::FVector{TangentSpaceType})
    copyto!(v.data, w.data)
    return v
end

function sharp!(M::Sphere, v::FVector{TangentSpaceType}, x, w::FVector{CotangentSpaceType})
    copyto!(v.data, w.data)
    return v
end

"""
    uniform_distribution(S::Sphere, x)

Uniform distribution on given sphere. Generated points will be of similar
type to `x`.
"""
function uniform_distribution(S::Sphere, x)
    d = Distributions.MvNormal(zero(x), 1.0)
    return ProjectedPointDistribution(S, d, project_point!, x)
end

"""
    check_manifold_point(S,x; kwargs...)

checks, whether `x` is a valid point on the [`Sphere`](@ref) `S`, i.e. is a vector
of length [`manifold_dimension`](@ref)`(S)+1` (approximately) of unit length.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_manifold_point(S::Sphere{N},x; kwargs...) where {N}
    if length(x) != N+1
        return DomainError(size(x),"The point $(x) does not lie on $S, since its size is not $(N+1).")
    end
    if !isapprox(norm(x), 1.; kwargs...)
        return DomainError(norm(x), "The point $(x) does not lie on the sphere $(S) since its norm is not 1.")
    end
    return nothing
end

"""
    check_tangent_vector(S,x,v; kwargs... )

checks whether `v` is a tangent vector to `x` on the [`Sphere`](@ref) `S`, i.e.
after [`check_manifold_point`](@ref)`(S,x)`, `v` has to be of same dimension as `x`
and orthogonal to `x`.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_tangent_vector(S::Sphere{N},x,v; kwargs...) where N
    perr = check_manifold_point(S,x)
    perr === nothing || return perr
    if length(v) != N+1
        return DomainError(size(v),
            "The vector $(v) is not a tangent to a point on $S since its size does not match $(N+1).")
    end
    if !isapprox( abs(dot(x,v)), 0.; kwargs...)
        return DomainError(abs(dot(x,v)),
            "The vector $(v) is not a tangent vector to $(x) on $(S), since it is not orthogonal in the embedding."
        )
    end
    return nothing
end

"""
    normal_tvector_distribution(S::Sphere, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(S::Sphere, x, σ)
    d = Distributions.MvNormal(zero(x), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(S), x, d, project_vector!, x)
end


text:
macro addctx(t, func)
    n = length(func.args[2].args)
    insert!(func.args[2].args, n, :(pop!(ps.closer.cc)))
    pushfirst!(func.args[2].args, :(push!(ps.closer.cc, $(t))))
    return esc(func)
end

macro enterctx(ps, t)
    :(push!($(esc(ps)).closer.cc, $t))
end

macro exitctx(ps, ret)
    quote
        pop!($(esc(ps)).closer.cc)
        $(esc(ret))
    end
end

function accept_rparen(ps)
    if ps.nt.kind == Tokens.RPAREN
        return PUNCTUATION(next(ps))
    else
        push!(ps.errors, Error((ps.ws.startbyte:ps.ws.endbyte) .+ 1 , "Expected )."))
        return ErrorToken(PUNCTUATION(Tokens.RPAREN, 0, 0))
    end
end
accept_rparen(ps::ParseState, args) = push!(args, accept_rparen(ps))

function accept_rsquare(ps)
    if ps.nt.kind == Tokens.RSQUARE
        return PUNCTUATION(next(ps))
    else
        push!(ps.errors, Error((ps.t.startbyte:ps.t.endbyte) .+ 1 , "Expected ]."))
        return ErrorToken(PUNCTUATION(Tokens.RSQUARE, 0, 0))
    end
end
accept_rsquare(ps::ParseState, args) = push!(args, accept_rsquare(ps))

function accept_rbrace(ps)
    if ps.nt.kind == Tokens.RBRACE
        return PUNCTUATION(next(ps))
    else
        push!(ps.errors, Error((ps.t.startbyte:ps.t.endbyte) .+ 1 , "Expected }."))
        return ErrorToken(PUNCTUATION(Tokens.RBRACE, 0, 0))
    end
end
accept_rbrace(ps::ParseState, args) = push!(args, accept_rbrace(ps))

function accept_end(ps::ParseState)
    if ps.nt.kind == Tokens.END
        return KEYWORD(next(ps))
    else
        push!(ps.errors, Error((ps.t.startbyte:ps.t.endbyte) .+ 1 , "Expected end."))
        return ErrorToken(KEYWORD(Tokens.END, 0, 0))
    end
end
accept_end(ps::ParseState, args) = push!(args, accept_end(ps))

function accept_comma(ps)
    if ps.nt.kind == Tokens.COMMA
        return PUNCTUATION(next(ps))
    else
        return PUNCTUATION(Tokens.RPAREN, 0, 0)
    end
end
accept_comma(ps::ParseState, args) = push!(args, accept_comma(ps))

function recover_endmarker(ps)
    if ps.nt.kind == Tokens.ENDMARKER
        if !isempty(ps.closer.cc)
            closert = last(ps.closer.cc)
            if closert == :block
                return ErrorToken(KEYWORD(Tokens.END, 0, 0))
            elseif closert == :paren
                return ErrorToken(PUNCTUATION(Tokens.RPAREN, 0, 0))
            elseif closert == :square
                return ErrorToken(PUNCTUATION(Tokens.RSQUARE, 0, 0))
            elseif closert == :brace
                return ErrorToken(PUNCTUATION(Tokens.RBRACE, 0, 0))
            end
        end
    end
end


text:
using Documenter, MLJScientificTypes, ScientificTypes

makedocs(
    modules = [MLJScientificTypes, ScientificTypes],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        ),
    sitename = "MLJScientificTypes.jl",
    authors = "Anthony Blaom, Thibaut Lienart, and contributors.",
    pages = [
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/alan-turing-institute/MLJScientificTypes.jl",
    push_preview = true
)


text:
copy_from=".";
expr_target = r"jl_[A-Za-z0-9]{6}";
expr_removebackslash = r"^(\\|/){1,2}";

notice = Markdown.parse("""
1. Please apply the following first:
```
using Weave, GeneralTools
```

2. You can alternatively specify the weave option in the front matter. For example, add the following to the front matter (header) of the your julia-markdown file:
```
weave_options:
  fig_path: ./
  hugo:
    out_path: ./_index.md
  md2pdf:
    out_path : ../pdf
  md2html:
    out_path : ./_index.html
```
Notice that in this way the path is relative to the path of julia markdown, not the current directory. That is, for that in the front matter we have `out_path: ./_index.md` and doing `weave("$copy_from/your_file.md")`, the result is the same as `weave(...; out_path="$copy_from/_index.md")`

""");

default_outname = "_index.md";
"""
`cp2content(dst::AbstractString)` copy the generated hugo markdown file to your site's `content`. `cp2content` looks in `dst` for `$default_outname` or `index.md` and images.
# Example
```
lazyhugo();
cp2content("D:/GoogleDrive/Sites/learning-notes/content/en/my-project/MagTIP");
```

# keyword arguments
- `copy_from="$copy_from"` is defined by default; that is, the markdown file generated by `lazyhugo()` and all images in the directory `copy_from` will be copied and send to `dst`.
- `rename_to` allows you to rename the output markdown file. For example, you'd like to `rename_to="my_project.md"`. Default is `rename_to=false`.

# Notice
$notice

"""
function cp2content(dst::AbstractString; force=true, copy_from=copy_from, rename_to=false)

    allfiles = readdir(copy_from, join=true);

    bnames = splitext.(allfiles); # base names
    isfig = map(x -> occursin(r"\.(png|jpe?g|svg|eps|gif|apng)",x[end]), bnames);
    isindex = occursin.(r"_?index\.md", basename.(allfiles));

    figures = allfiles[isfig];
    theindex = allfiles[isindex];

    if length(theindex) != 1
        openit(joinpath(pwd(),copy_from));
        error("No or more than one `index.md` file(s). `weave` first and check the opened folder!");
    end

    src = [figures;theindex];	# make sure the "index.md" be the last element of src
    dsts = joinpath.(dst,basename.(src)); # hence "index.md" be the last element of dsts

    if isa(rename_to, String)
        dsts[end] = renamefile(dsts[end], rename_to);
    end

    cp.(src, dsts; force=force);

    println("Files copied to $dst except:");
    println.(allfiles[.!(isfig .| isindex)]);
end

default_outpath = "_index.md";
default_figpath = "";
"""
`lazyhugo([filename]; out_path="$default_outpath")` is equivalent to `weave("\$filename"; doctype="hugo", out_path="$default_outpath", fig_path="$default_figpath")`.

Simply `lazyhugo()` automatically search for a julia markdown file (`.jmd`) of any name first, then the markdown file not named as "_index.md" or "index.md".
Once a unique juliamarkdown file is found, it quit searching immediately; but if multiple or no (julia) markdown files were found, an error will be raised.

If `filename` is suffixed with `.md` (a markdown file instead of julia markdown), then weave option `informat="markdown"` will be added. If the first argument is a folder name, then `lazyhugo` search julia markdown then markdown files inside the folder.

`lazyhugo` also reformat `{{< figure src="/using_gadfly_14_1.png"  >}}`. The `/using_gadfly_14_1.png` is captured with the first slash being deleted, leading to `![](using_gadfly_14_1.png)`

# Notice
$notice

"""
function lazyhugo(; out_path=default_outpath)
  opt = [];
  allfiles = readdir();
  id = occursin.(r"\.jmd",allfiles); # check if there is .jmd file
  if all(.!id) # if there is not any .jmd exist, then find something with extension .md.
      id = occursin.(r"(?<!index)\.md",allfiles);
  end

  filenames = allfiles[id];
  if length(filenames) != 1
      error("No or multiple markdown files exists. Please explicitily specify the target file you want to weave. For example, `lazyhugo(Doc.md)` or `lazyhugo(Doc.jmd)`. Noted that your markdown file to be converted cannot be named as `index.md`.");
  end

  filename = filenames[1];
  println("`$filename` is going to be weaved:");
  _lazyhugo(filename, out_path, opt);
end

function lazyhugo(filename::AbstractString; out_path=default_outpath)
  opt = [];
  if !isfile(filename) # which means the first arg is expected to be a folder containing `.jmd` or `_index.md`
      currentdir = pwd();
      cd(filename);

      try
          lazyhugo(; out_path=out_path);
      catch e
          cd(currentdir);
          throw(e);
      end
      cd(currentdir);
      return
  end

  if splitext(filename)[end] == ".md"
      opt = [:informat => "markdown"];
  end
  _lazyhugo(filename, out_path, opt);
end

function _lazyhugo(filename::AbstractString, out_path, opt)
  weave("$filename"; doctype="hugo", out_path=out_path, fig_path=default_figpath, opt...);
  # Reformat figure referencing and remove extra first slash/backslash if there is any, because the slash/backslash makes referencing fail.
  txt = readlines(out_path);
  lines2replace = occursin.(r"(?<=figure src=)", txt);
  if !any(lines2replace)
    return
  end
  matches = match.(r"(?<=\").*(?=\")",txt[lines2replace]);
  newfigref = String[];
  for mt in matches
    str0 = mt.match;
    str1 = replace(str0, expr_removebackslash => "");
    newref = "![]($str1)"
    push!(newfigref, newref);
  end

  txt[lines2replace] .= newfigref;

  open(out_path, "w") do io
    for line in txt
      println(io, line);
    end
  end
end

"""
`getdoc(filepath::String; first_head=0)`
read the function descriptions and return a `Markdown.parse`d object.
Now it supports matlab's `.m` function only.
TODO: add julia support.

# Example (matlab)
- `y = getdoc("..\\..\\src\\statind.m")` get all the comments (starts with `% ` eachline) `before function statind(...)`
- `y = getdoc("..\\..\\src\\statind.m", first_head=3)` get all the comments before function section and rearrange all the heading levels based on forcing the first heading level to be 3. That is, take the first heading to be `#` for example, `#` becomes `###` and `##` becomes `####`
"""
function getdoc(filepath::String; first_head=0)

  s = open(filepath) do file
    readlines(file);
    end
  if splitext(filepath)[end] == ".m" # if it is the matlab file
    doc = _getrawcomment(s);
    x = join(lang_matlab(doc,first_head = first_head));
    y = Markdown.parse(x); # or @eval @md_str $x is the same.
  else
    y = Markdown.parse("<unsupported language>");
  end

  return y
end

function _getrawcomment(s)
  funbegin = findfirst(occursin.(r"^(function|classdef)\s",s));
  # find the line where function begin
  if funbegin == 1
    return "No documentation."
  end

  if isnothing(funbegin)
    docend = length(s);
  else
    docend = funbegin - 1;
  end

  doc = s[1:docend];
end

function lang_matlab(doc; first_head=0)
  # for Weave.jl, both default or pandoc2html output, mulitiple empty lines == one and results in a linebreak, and linebreak (but no empty line in between) is ignored. So I add linebreaks very generously here.
  br = "\n"; # linebreak
  # delete "%" in the beginning of the line.
  doc2 = replace.(doc, r"^%+\s?" => "");

  if first_head != 0 # reformat markdown head levels
    lenhashtags0 = match(r"^#+",doc2[1]).match |> length;
    adjust_diff = first_head - lenhashtags0;

    ishead = occursin.(r"^#+",doc2[1:end]);
    targetid = findall(ishead);
    for i in targetid
      lenhashtags_i = match.(r"^#+",doc2[i]).match |> length;
      new_head = "#"^(lenhashtags_i+adjust_diff);
      doc2[i] = replace(doc2[i], r"^#+" => new_head);
    end
  end

  for (i,val) in enumerate(doc2)
    doc2[i] = val*br; # add linebreak at the end of every line.
  end

  # codefence = findall(occursin.(r"^(```)",doc2));
  # codestart = codefence[1:2:end];
  # codeend = codefence[2:2:end];

  # numinserted = 0;

  # for id in codefence
  # 	id = id + numinserted;
  # 	# if iseven(numinserted)
  # 	# 	insertat = id; # insert before
  # 	# else
  # 	# 	insertat = id+1;  # insert after
  # 	# end
  # 	insert!(doc2, id, br); # add linebreak before
  # 	insert!(doc2, id+2, br); # add linebreak after
  # 	numinserted = numinserted + 2;
  # end
  return doc2
end

"""
`defolder(dir::AbstractString)`
Delete all folders left everytime we use the `weave(...;out_path=dir)`. The target folders are captured with regular expression pattern `r"jl_[A-Za-z0-9]{6}"`.

`defolder()` search targets at current directory.

"""
function defolder(dir::AbstractString;force=false)
  _defolder(dir, force);
end

function defolder(;force=false)
  dir = pwd();
  _defolder(dir, force);
end

function _defolder(dir, force::Bool)

  dir = abspath(dir); # reformat the path to prevent something like "dir/to/target\\jl_x3as5e"
  flist = readdir(dir, join=true);
  dirtrue = isdir.(flist);
  fmtmatched = .!isnothing.(match.(expr_target,flist));
  targets = flist[dirtrue .& fmtmatched];
  if isempty(targets)
      println("No target folder found. Do nothing.");
      return nothing
  end

  try
    rm.(targets; recursive=force);
  catch e
    str = sprint(showerror, e);
    openit(dir);
    @warn ("Some of the target folders are not empty [$str]. Set `force=true` to delete all of them without prompt.");
    println("Delete them anyway? (y/n)");
    if readline(stdin) == "y"
      _defolder(dir, true);
    end
  end

end

max_stars_repo_name: CGRG-lab/HypertextTools

max_issues_repo_name: CGRG-lab/HypertextTools

max_forks_repo_name: CGRG-lab/HypertextTools


text:
using Documenter

makedocs(
    format = Documenter.HTML(analytics = "UA-367975-10", mathengine = Documenter.MathJax()),
    sitename = "MaximinOPF",
    authors = "Sam",
    pages = [
        "Home" => "index.md",
        "Manual" => [
        	"Getting Started" => "gettingstarted.md",
        	"API Document" => "API.md",
        	"Examples" => "example.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/yim0331/yim0331.github.io.git",
)


text:
## Things that can't be loaded twice.g
# Some instances of data structures

## Symbol correspondence/translation between Julia and Symata

# This should be separated better. Some symbols are translated
# only on analyzing the ast, others, on both input and output.
# Fix. We delete the keys we don't want (unicode) from the reverse translation.

#### Input,Output translation table.

# These translations are made when encountered in the AST and on printing.
# A reversed dict is constructed for output.
# Excepted are the unicode characters. They are merged into the input dict.
# But, they are held in a separate table for output, so that unicode
# output can be disabled.

# Idea, but not implemented:
# We could use  a.b for contexts, modules.  In Mma
# a . b == a.b  == Dot[a,b], matrix and vector multiplication.
# a`b is symbol b in context a.  a .* b is illegal in Mma, but a .* b is legal syntax in Julia
# So: in Symata  a.b for contexts,  and a .* b for matrix multiplication

# We can use $ for (Symata) Function, but the translation to compound expression is broken
# We can use -> for (Symata) Function. This seems to work

## Symbols not yet in use
## .\
const JTOMSYM  =
 Dict(
     :(=) => :Set,
     :(:=) => :SetDelayed,
     :(^=) => :UpSet,
     :+ => :Plus,
     :- => :Minus,
     :* => :Times,
     :⋅ => :Dot,
     :^ => :Power,
     :(*=) => :TimesBy,
     :(+=) => :AddTo,
     # $ is good for Function below, but we might want it for interpolation instead
     # No, it does not interfere with interpolation, but it still wraps args in Comound expression. dont know why
     :(->) => :Function,
     :$ => :Function, #  $(1+~2) -->  Function(Plus(1,Slot(2))), but we don't do anything yet with this
     :~ => :Slot,     # we need to translate lone ~ , ie Slot into Slot(1)
     :(=>) => :Rule, # Mma uses ->. We use => because it has a higher precedence than ->, which is what we want.
     #      :(->) => :RuleDelayed,     # Mma uses :>. Julia parser does not allow this
     :(.>) => :RuleDelayed,      # Mma uses :>. Julia parser does not allow this.  .> has better precedence and parsing than ->
     :(./) => :ReplaceAll,       # Mma has /. for this !! But, /. is not legal Julia syntax
     :(.//) => :ReplaceRepeated,
     :(%) => :Map,
     :(.%) => :Apply,
     #      :(:) => :Span, # in v0.7, colon has become a call. We discriminate meaning by context
     :(&=) => :UpSetDelayed,  # This is available, at least. Not sure we want to take it.
     :... => :... ,  # We still need to decide what to do with this. Maybe evaluate the args and Apply Sequence
     :.. => :..,
     :! => :Not,
     :&& => :And,
     :|| => :Or,
     :| => :Alternatives,
     :(...) => :Repeated
 )

const JTOMSYM_ONEWAY  =
    Dict(
        :vcat => :List,
        :vect => :List,
        :ref => :Part,
        :cell1d => :List,   # curly brackets, but deprecated by julia
        :comparison => :Comparison,
    )

# FIXME  NB Alternatives: This is nary, but not Flat. So we need to parse
# a | b | c as Alternatives(a,b,c) rather than nested Alternatives as we do now.
# (a | b) | c is nested, and is not flattend.
# Other operators, such as Times are nary but Flat, so the flattening happens during the
# evaluations sequence.

# These symbols are unconditionally translated on input
# This is almost exactly the Dict unicode_translation. Maybe we can refactor.
const INSYMTRANS = Dict( :⇒  => :(=>),
                         :π => :Pi,
                         :≥ => :(>=),
                         :≤ => :(<=),
                         :γ => :EulerGamma,
                         :Γ => :Gamma,
                         :𝕖 => :E,
                         :𝕚 => :I,
                         :∞ => :Infinity,
                         :≠  =>  :!=,
                         :∈  => :Element
                         )

# Input translation. Output translation if unicode is enabled
const unicode_translation = Dict{Symbol,Symbol}(:π => :Pi,
                                                :γ => :EulerGamma,
                                                :∞ =>   :Infinity,
                                                :Γ =>   :Gamma,
                                                :𝕚 =>   :I,
                                                :≥  =>  :>=,
                                                :≤  =>  :<=,
                                                :≠  =>  :!=,
                                                :𝕖  =>  :E,
                                                :⇒  =>  :(=>),
                                                :→ => :Function,
                                                :∈ => :Element
)

const comparison_translation = Dict(
                                    :(==) => :Equal,
                                    :(!=) => :Unequal,
                                    :(>) => :Greater,
                                    :(<) => :Less,
                                    :(>=) => :GreaterEqual,
                                    :(<=) => :LessEqual,
                                    :(===) => :SameQ,
                                    :(!==) => :UnsameQ
                                    )

const inverse_comparison_translation = Dict{Symbol,Symbol}()
for x in comparison_translation
    inverse_comparison_translation[x[2]] = x[1]
end

# Output. Reverse dict for printing if unicode printing is enabled
const unicode_output = Dict{Symbol,Symbol}()
for (k,v) in unicode_translation unicode_output[v] = k end

# Reverse dict of ascii symbols for printing
const MTOJSYM = Dict{Symbol,Symbol}()
for (k,v) in JTOMSYM  MTOJSYM[v] = k end

# Finally, add the unicode symbols to the input translation table.
merge!(JTOMSYM, unicode_translation)
merge!(JTOMSYM, JTOMSYM_ONEWAY)

# This is only used for output. (Can't we detect it with the others on input ?)
MTOJSYM[:Span] = :(:)

merge!(MTOJSYM, inverse_comparison_translation)

function jtomsym(x::Symbol)
    check_autoload(x)  # trigger autoload of Symata language code
    if haskey(JTOMSYM,x)
        return JTOMSYM[x]
    end
    return x
end

jtomsym(x) = extomx(x)

function mtojsym(x::Symbol)
    # For output, reverse the rewrites done by preprocessing the input string befor parsing
    if haskey(REVERSE_PREPROCESS_SYMBOL_TRANSLATION,x)  # .>  --> :>   for output
         x = REVERSE_PREPROCESS_SYMBOL_TRANSLATION[x]   # =>  --> ->
     end
    if haskey(MTOJSYM,x)
        return MTOJSYM[x]
    end
    return x
end

# Non-Symbol Heads translate to themselves
mtojsym(x) = x

# These will not be translated back to Julia expression symbol on printing
# TODO. choose tuple or block, or both ?
# One unusual occurance of :block is filtered before we get here

# A compound expression is denoted by tuple () or by begin end.
# We could possible used tuple for something else.
for (k,v) = ( (:tuple,:CompoundExpression), (:block,:CompoundExpression), (:dict, :List), (:parameters, :TestExpression) )
    JTOMSYM[k] = v
end

#### Operator type: infix, prefix, postfix

const OPTYPE  = Dict{Symbol,Symbol}()

for op in (:(=), :(:=), :(=>), :Rule , :RuleDelayed, :Power, :(.>),
           :Set, :SetDelayed, :UpSet, :(*=), :(+=), :→, :(->), :Function,
           :TimesBy, :AddTo  # , :Element
           ) # need :Set here
    OPTYPE[op] = :binary
end

for op in keys(inverse_comparison_translation)
    OPTYPE[op] = :binary
end

for op in (:Plus, :Times, :Dot, :Span, :And, :Or, :Alternatives)
    OPTYPE[op] = :infix
end

# Nonsymbolic Heads, Integer, etc. assume they are prefix ops
getoptype(x) = :prefix

# function getoptype(x::Mxpr)
#     println("getting opt type for ", mhead(x), " args ", margs(x))
#     res = getoptype(mhead(x))
#     println("Resutt is $res")
#     res
# end

function getoptype(x::Symbol)
    if haskey(OPTYPE,x)
        return OPTYPE[x]
    end
    return :prefix
end

#### Comparison symbols

# This is only used in AST_tranlation.jl to make a Comparison Mxpr.
const COMPARISONSYMBOLS = Dict{Symbol,Bool}()

for op in ( :(==), :(<), :(>), :(>=), :(<=), :(===),
           :(!=), :(!==), :(<:))
    COMPARISONSYMBOLS[op] = true
end

function is_comparison_symbol(x::Symbol)
    if haskey(COMPARISONSYMBOLS,x)
        return true
    end
    return false
end

is_comparison_symbol(x) = false


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674

# %%
using DifferentialEquations
using Parameters
using StaticArrays
using Plots

# %%
function f(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    v = SVector(k1*A*B^2, kf*E*S, kr*ES, kcat*ES)
    M*v
end

# %%
M = [
    -1  0  0  0
    -2  0  0  0
    +1 -1 +1  0
     0 -1 +1 +1
     0 +1 -1 -1
     0  0  0 +1
]

p = (
    M = SMatrix{6, 4, Float64}(M),
    k1 = 1.0,
    kf = 1.0,
    kr = 1.0,
    kcat = 1.0,
)

u0 = SVector(
    #= A  =# 1.5,
    #= B  =# 1.2,
    #= S  =# 0.0,
    #= E  =# 1.0,
    #= ES =# 0.0,
    #= P  =# 0.0,
)

tspan = (0.0, 10.0)

# %%
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

# %%
function g(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    dA = -k1*A*B^2
    dB = -2k1*A*B^2
    dS = k1*A*B^2 - kf*E*S + kr*ES
    dE  = -kf*E*S + kr*ES + kcat*ES
    dES =  kf*E*S - kr*ES - kcat*ES
    dP = kcat*ES
    SVector(dA, dB, dS, dE, dES, dP)
end

# %%
prob = ODEProblem(g, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

# %%

max_stars_repo_path: 0019/Petersen matrix.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/Petersen matrix.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/Petersen matrix.jl

max_forks_repo_name: genkuroki/public


text:
##############################################################################################
module L96
##############################################################################################
# imports and exports
using Random, Distributions
using LinearAlgebra, SparseArrays
export dx_dt, jacobian, l96s_tay_2_step!, ρ, α

##############################################################################################
##############################################################################################
# Type union declarations for multiple dispatch
# and type aliases

# vectors and ensemble members of sample
VecA = Union{Vector{Float64}, SubArray{Float64, 1}}

# dictionary for model parameters
ParamDict = Union{Dict{String, Array{Float64}}, Dict{String, Vector{Float64}}}

##############################################################################################
##############################################################################################
# auxiliary function to return modular indices for the lorenz model

function mod_indx!(indx::Int64, dim::Int64)
    indx = mod(indx, dim)
    if indx==0
        indx = dim
    end
    return indx
end

##############################################################################################
# time derivative

function dx_dt(x::VecA, t::Float64, dx_params::ParamDict)
	"""Time derivative for Lorenz-96 model, x is a single model state of size state_dim."""

    # unpack the (only) derivative parameter for l96
    F = dx_params["F"][1]::Float64
    x_dim = length(x)
    dx = Vector{Float64}(undef, x_dim)

    for j in 1:x_dim
        # index j minus 2, modulo the system dimension
        j_m_2 = mod_indx!(j - 2, x_dim)

        # index j minus 1, modulo the system dimension
        j_m_1 = mod_indx!(j - 1, x_dim)

        # index j plus 1, modulo the system dimension
        j_p_1 = mod_indx!(j + 1, x_dim)

        dx[j] = (x[j_p_1] - x[j_m_2])*x[j_m_1] - x[j] + F
    end
    return dx
end

##############################################################################################
# vectorized time derivative for demonstration purposes only,
# note this is much less efficient in speed and in memory than the above version

function dx_dt_alt(x::VecA, t::Float64, dx_params::ParamDict)

    # unpack the (only) derivative parameter for l96
    F = dx_params["F"][1]::Float64
    x_dim = length(x)
    dx = Vector{Float64}(undef, x_dim)

    # shift minus and plus indices, array compute the derivative over arbitrary ensemble sizes
    x_m_2 = [x[end-1:end]; x[1:end-2]]
    x_m_1 = [x[end:end]; x[1:end-1]]
    x_p_1 = [x[2:end]; x[1:1]]

    # compute the vectorized derivative
    dx .= (x_p_1-x_m_2) .* x_m_1 - x .+ F
end

##############################################################################################
# linearized time derivative

function jacobian(x::Vector{Float64}, t::Float64, dx_params::ParamDict)
    """"Computes the Jacobian of Lorenz-96 about the state x.

    Note that this is designed to load entries in a zeros array and return a sparse array
    to make a compromise between memory and computational resources."""

    x_dim = length(x)
    dxF = zeros(x_dim, x_dim)

    # looping columns j of the jacobian, loading the standard matrix
    for j in 1:x_dim

        # index j minus 2, modulo the system dimension
        j_m_2 = mod_indx!(j - 2, x_dim)

        # index j minus 1, modulo the system dimension
        j_m_1 = mod_indx!(j - 1, x_dim)

        # index j plus 1, modulo the system dimension
        j_p_1 = mod_indx!(j + 1, x_dim)

        # index j plus 2, modulo the system dimension
        j_p_2 = mod_indx!(j + 2, x_dim)

        # load the jacobian entries in corresponding rows
        dxF[j_p_2, j] = -x[j_p_1]
        dxF[j_p_1, j] = x[j_p_2] - x[j_m_1]
        dxF[j, j] = -1.0
        dxF[j_m_1, j] = x[j_m_2]
    end
    return sparse(dxF)
end

##############################################################################################
# Auxiliary functions for the 2nd order Taylor-Stratonovich expansion below. These need
# to be computed once, only as a function of the order of truncation of the Fourier
# series, p for full timeseries.

function ρ(p::Int64)
    1.0/12.0 - 0.5 * π^(-2.0) * sum(1.0 ./ Vector{Float64}(1:p).^2.0)
end

function α(p::Int64)
    (π^2.0) / 180.0 - 0.5 * π^(-2.0) * sum(1.0 ./ Vector{Float64}(1:p).^4.0)
end

##############################################################################################
# 2nd order strong taylor SDE step
# This method is derived in
# Grudzien, C. et al.: On the numerical integration of the Lorenz-96 model,
# with scalar additive noise, for benchmark twin experiments,
# Geosci. Model Dev., 13, 1903–1924, https://doi.org/10.5194/gmd-13-1903-2020, 2020.
# This depends on rho and alpha as above
# NOTE: this Julia version still pending validation as in the above manuscript

function l96s_tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any})
    """One step of integration rule for l96 second order taylor rule

    The ρ and α are to be computed by the auxiliary functions, depending only on p,
    and supplied for all steps. This is the general formulation which includes,
    eg. dependence on the truncation of terms in the auxilliary function C with
    respect to the parameter p.  In general, truncation at p=1 is all that is
    necessary for order 2.0 convergence, and in this case C below is identically
    equal to zero.  This auxilliary function can be removed (and is removed) in other
    implementations for simplicity."""

    # Infer model and parameters
    sys_dim = length(x)
    dx_params = kwargs["dx_params"]::ParamDict
    h = kwargs["h"]::Float64
    diffusion = kwargs["diffusion"]::Float64
    p = kwargs["p"]::Int64
    ρ = kwargs["ρ"]::Float64
    α = kwargs["α"]::Float64

    # Compute the deterministic dxdt and the jacobian equations
    dx = dx_dt(x, 0.0, dx_params)
    Jac_x = jacobian(x, 0.0, dx_params)

    ## random variables
    # Vectors ξ, μ, ϕ are sys_dim X 1 vectors of iid standard normal variables,
    # ζ and η are sys_dim X p matrices of iid standard normal variables.
    # Functional relationships describe each variable W_j as the transformation of
    # ξ_j to be of variace given by the length of the time step h. Functions of random
    # Fourier coefficients a_i, b_i are given in terms μ / η and ϕ / ζ respectively.

    # draw standard normal samples
    rndm = rand(Normal(), sys_dim, 2*p + 3)
    ξ = rndm[:, 1]

    μ = rndm[:, 2]
    ϕ = rndm[:, 3]

    ζ = rndm[:, 4: p+3]
    η = rndm[:, p+4: end]

    ### define the auxiliary functions of random fourier coefficients, a and b

    # denominators for the a series
    denoms = repeat((1.0 ./ Vector{Float64}(1:p)), 1, sys_dim)

    # vector of sums defining a terms
    a = -2.0 * sqrt(h * ρ) * μ - sqrt(2.0*h) * sum(ζ' .* denoms, dims=1)' / π

    # denominators for the b series
    denoms = repeat((1.0 ./ Vector{Float64}(1:p).^2.0), 1, sys_dim)

    # vector of sums defining b terms
    b = sqrt(h * α) * ϕ + sqrt(h / (2.0 * π^2.0) ) * sum(η' .* denoms, dims=1)'

    # vector of first order Stratonovich integrals
    J_pdelta = (h / 2.0) * (sqrt(h) * ξ + a)

    ### auxiliary functions for higher order stratonovich integrals ###
    # C function is optional for higher precision but does not change order of convergence
    function C(l, j)
        if p == 1
            return 0.0
        end
        c = zeros(p, p)
        # define the coefficient as a sum of matrix entries where r and k do not agree
        indx = Set(1:p)

        for r in 1:p
            # vals are all values not equal to r
            vals = setdiff(indx, Set(r))
            for k in vals
                # and for column r, define all row entries below, with zeros on diagonal
                c[k, r] = (r / (r^2 - k^2)) * ((1.0 / k) * ζ[l, r] * ζ[j, k] + (1.0 / r) *
                                               η[l, r] * η[j, k])
            end
        end

        # return the sum of all values scaled by -1/(2pi^2)
        -0.5 * π^(-2.0) * sum(c)
    end

    function Ψ(l, j)
        # Ψ - generic function of the indicies l and j, define Ψ plus and Ψ minus index-wise
        h^2.0 * ξ[l] * ξ[j] / 3.0 + h * a[l] * a[j] / 2.0 +
        h^(1.5) * (ξ[l] * a[j] + ξ[j] * a[l]) / 4.0 -
        h^(1.5) * (ξ[l] * b[j] + ξ[j] * b[l]) / (2.0 * π) - h^2.0 * (C(l,j) + C(j,l))
    end

    # define the approximations of the second order Stratonovich integral
    Ψ_plus = Vector{Float64}(undef, sys_dim)
    Ψ_minus = Vector{Float64}(undef, sys_dim)
    for i in 1:sys_dim
        Ψ_plus[i] = Ψ(mod_indx!((i-1), sys_dim), mod_indx!((i+1), sys_dim))
        Ψ_minus[i] = Ψ(mod_indx!((i-2), sys_dim), mod_indx!((i-1), sys_dim))
    end

    # the final vectorized step forward is given as
    x .= collect(Iterators.flatten(
            x + dx * h + h^2.0 * 0.5 * Jac_x * dx +  # deterministic taylor step
            diffusion * sqrt(h) * ξ +                # stochastic euler step
            diffusion * Jac_x * J_pdelta +           # stochastic first order taylor step
            diffusion^2.0 * (Ψ_plus - Ψ_minus)       # stochastic second order taylor step
           ))
end

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
module PTools

include("srvr_task.jl")
include("j_shm.jl")
include("queued_tasks.jl")
include("pfork.jl")
include("worker_set.jl")

global _debug = false
function _set_debug(d)
    global _debug
    _debug = d
end

#export STask
#export Jshm

end


text:
# Unit tests for implementations/nbodysimulator/nbodysimulator_simulator.jl

@testset "nbodysimulator_simulator.jl" begin
    # TODO
end


text:
"""
    finite-difference velocity objective
"""
struct VelocityObjective <: Objective
    Q
    n
    h

    idx_angle
end

function velocity_objective(Q, n; h = 0.0, idx_angle = (1:0))
    VelocityObjective(Q, n, h, idx_angle)
end

function angle_diff(a1, a2)
    mod(a2 - a1 + pi, 2 * pi) - pi
end

function objective(Z, obj::VelocityObjective, idx, T)
    n = obj.n

    J = 0.0
    ϵ = 1.0e-5
    for t = 1:T-1
        q⁻ = view(Z, idx.x[t][n .+ (1:n)])
        q⁺ = view(Z, idx.x[t + 1][n .+ (1:n)])
        h = obj.h == 0.0 ? 1.0 : obj.h
        v = (q⁺ - q⁻) ./ h
        v[obj.idx_angle] = angle_diff.(view(q⁻, obj.idx_angle), view(q⁺, obj.idx_angle)) ./ h

        J += v' * obj.Q[t] * v
    end

    return J
end

function objective_gradient!(∇J, Z, obj::VelocityObjective, idx, T)
    # tmp(y) = objective(y, obj, idx, T)
    # ∇J .+= ForwardDiff.gradient(tmp, Z)
    n = obj.n

    J = 0.0
    ϵ = 1.0e-5
    for t = 1:T-1
        q⁻ = view(Z, idx.x[t][n .+ (1:n)])
        q⁺ = view(Z, idx.x[t + 1][n .+ (1:n)])
        h = obj.h == 0.0 ? 1.0 : obj.h

        v = (q⁺ - q⁻) ./ h
        v[obj.idx_angle] = angle_diff.(view(q⁻, obj.idx_angle), view(q⁺, obj.idx_angle)) ./ h

        dJdv = 2.0 * obj.Q[t] * v
        ∇J[idx.x[t][n .+ (1:n)]] += -1.0 ./ h * dJdv
        ∇J[idx.x[t + 1][n .+ (1:n)]] += 1.0 ./ h * dJdv
    end

    return nothing
end

# function objective(Z, obj::VelocityObjective, idx, T)
#     n = obj.n
#
#     J = 0.0
#
#     for t = 1:T-1
#         q⁻ = view(Z, idx.x[t][n .+ (1:n)])
#         q⁺ = view(Z, idx.x[t + 1][n .+ (1:n)])
#         h = obj.h == 0.0 ? view(Z, idx.u[t][end]) : obj.h
#         v = (q⁺ - q⁻) ./ h
#         v[obj.idx_angle] = angle_diff.(view(q⁻, obj.idx_angle), view(q⁺, obj.idx_angle)) ./ h
#         J += v' * obj.Q[t] * v
#     end
#
#     return J
# end
#
# function objective_gradient!(∇J, Z, obj::VelocityObjective, idx, T)
#     # _tmp(y) = objective(y, obj, idx, T)
#     # ∇J .= ForwardDiff.gradient(_tmp, Z)
#     n = obj.n
#
#     for t = 1:T-1
#         q⁻ = view(Z, idx.x[t][n .+ (1:n)])
#         q⁺ = view(Z, idx.x[t + 1][n .+ (1:n)])
#         h = obj.h == 0.0 ? view(Z, idx.u[t][end]) : obj.h
#
#         # v = (q⁺ - q⁻) ./ h
#         # # v[obj.idx_angle] = angle_diff.(view(q⁻, obj.idx_angle),
#         # #     view(q⁺, obj.idx_angle)) ./ h
#         #
#         # dJdv = 2.0 * obj.Q[t] * v
#
#         function tmp(x⁻, x⁺, g)
#             w = (x⁺ - x⁻) ./ g
#             w[obj.idx_angle] = angle_diff.(view(x⁻, obj.idx_angle), view(x⁺, obj.idx_angle)) ./ g
#             w' * obj.Q[t] * w
#         end
#
#         tmp1(y) = tmp(y, q⁺, h)
#         tmp2(y) = tmp(q⁻, y, h)
#         tmp3(y) = tmp(q⁻, q⁺, y)
#
#         ∇J[idx.x[t][n .+ (1:n)]] += ForwardDiff.gradient(tmp1, q⁻)
#         ∇J[idx.x[t + 1][n .+ (1:n)]] += ForwardDiff.gradient(tmp2, q⁺)
#         if obj.h == 0.0
#             ∇J[idx.u[t][end]] += ForwardDiff.gradient(tmp3, h)[1]
#         end
#     end
#
#     return nothing
# end

# using ForwardDiff
# b1 = -pi / 2.0
# b2 = pi / 2.0 + pi
#
# angle_diff(b1, b2)
#
# tmp1(y) = angle_diff(y, b2)
# tmp2(y) = angle_diff(b1, y)
# for i = 1:1000
#     b1 = rand(1)[1]
#     b2 = rand(1)[1]
#     h1 = rand(1)[1]
#     tmp1(y) = angle_diff(y, b2) / h1
#     tmp2(y) = angle_diff(b1, y) / h1
#     @assert ForwardDiff.derivative(tmp1, b1) == -1.0 / h1
#     @assert ForwardDiff.derivative(tmp2, b2) == 1.0 / h1
# end
#
# angle_diff.(ones(2), pi * ones(2))