text:
# Uppermost types in the schema.  Most of these extend the BaseNodeType.
# Since there are only three types in the schema which extend the BaseNodeType,
# we simply reproduce the fields of the BaseNodeType directly within each
# struct.  This allows each field to be documented by DocStringExtensions
# which is not possible using Mixers.@pour.
#
# FIXME: Find a way to allow documented fields to be inserted into a struct.

@enumerated_struct(ChannelType, ("TRIGGERED", "CONTINUOUS", "HEALTH", "GEOPHYSICAL",
            "WEATHER", "FLAG", "SYNTHESIZED", "INPUT", "EXPERIMENTAL",
            "MAINTENANCE", "BEAM"))

"""
    Channel

A channel is a time series recording of a component of some observable,
often colocated with other channels at the same location of a station.

Equivalent to SEED blockette 52 and parent element for the related the response blockettes.

!!! note
    The presence of a `sample_rate_ratio` without a `sample_rate` field is
    not allowed in the standard, but it permitted by StationXML.jl.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Channel
    # BaseNodeType fields
    description::M{String} = missing
    "Persistent identifier for this `Channel`"
    identifier::Vector{Identifier} = Identifier[]
    comment::Vector{Comment} = Comment[]
    "A description of time series data availablility for this `Channel`"
    data_availability::M{DataAvailability} = missing
    "Channel's code (e.g., `\"BHE\"`)."
    code::String
    "Start date of operation of this channel of the containing station."
    start_date::M{DateTime} = missing
    "End date of operation of this channel.  May be `missing`."
    end_date::M{DateTime} = missing
    "A data source identifier in URI form"
    source_id::M{String} = missing
    "Information on whether the data in this channel are restricted, open,
     or otherwise.  (See [`RestrictedStatus`](@ref StationXML.RestrictedStatus).)"
    restricted_status::M{RestrictedStatus} = missing
    "A code used for display or association, alternate to the SEED-compliant code."
    alternate_code::M{String} = missing
    "A previously used code if different from the current code."
    historical_code::M{String} = missing
    # end of BaseNodeType fields
    "URI of any type of external report, such as data quality reports."
    external_reference::Vector{ExternalReference} = ExternalReference[]
    "Latitude coordinate of this channel's sensor."
    latitude::Latitude
    "Longitude coordinate of this channel's sensor."
    longitude::Longitude
    "Elevation of the sensor."
    elevation::Distance
    "The local depth or overburden of the instrument's location. For downhole
     instruments, the depth of the instrument under the surface ground level.
     For underground vaults, the distance from the instrument to the local ground level above."
    depth::Distance
    "Azimuth of the sensor in degrees from north, clockwise."
    azimuth::M{Azimuth} = missing
    "Dip of the instrument in degrees, down from horizontal."
    dip::M{Dip} = missing
    "Elevation of the water surface in meters for underwater sites, where 0 is sea level."
    water_level::M{Float} = missing
    "The type of data this channel collects. Corresponds to
     channel flags in SEED blockette 52. The SEED volume producer could
     use the first letter of an Output value as the SEED channel flag."
    type::Vector{ChannelType} = ChannelType[]
    # SampleRateGroup fields
    "Sampling rate of the channel in samples/s.  If `sample_rate_ratio` is
     also included, then `sample_rate` is more definitive."
    sample_rate::M{SampleRate} = missing
    "Sampling rate of the channel in number of samples per number of seconds."
    sample_rate_ratio::M{SampleRateRatio} = missing
    # end of SampleRateGroup fields
    "The storage format of the recorded data (e.g. SEED).
     **Removed in StationXML v1.1.**"
    storage_format::M{String} = missing
    "A tolerance value, measured in seconds per sample, used as a threshold for time
     error detection in data from the channel."
    clock_drift::M{ClockDrift} = missing
    calibration_units::M{Units} = missing
    sensor::M{Equipment} = missing
    pre_amplifier::M{Equipment} = missing
    data_logger::M{Equipment} = missing
    equipment::Vector{Equipment} = Equipment[]
    "The transfer function describing how ground motion is translated into
    digital counts."
    response::M{Response} = missing
    "Code defining the location of the channel (may be empty)."
    location_code::String
end

"""
    attribute_fields(T::Type)

Return the fields of type `T` which are stored as attributes in the
StationXML specification
"""
attribute_fields(::Type{Channel}) = (BASE_NODE_ATTRIBUTES..., :location_code)

"""
    removed_fields(T::Type)

Return the fields of type `T` which are no longer present in the latest StationXML
version.
"""
removed_fields(::Type{Channel}) = (:storage_format,)

"""
    Station

This type represents a Station epoch. It is common to only have a
single station epoch with the station's creation and termination dates as the epoch
start and end dates.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Station
    # BaseNodeType fields
    description::M{String} = missing
    "Persistent identifier for this `Station`"
    identifier::Vector{Identifier} = Identifier[]
    comment::Vector{Comment} = Comment[]
    "A description of time series data availablility for this `Station`"
    data_availability::M{DataAvailability} = missing
    "Station's code, e.g. `\"ANMO\"`."
    code::String
    "Start date of operation of the station."
    start_date::M{DateTime} = missing
    "End date of operation of the station."
    end_date::M{DateTime} = missing
    "A data source identifier in URI form"
    source_id::M{String} = missing
    "Information on whether the data from this station are restricted, open,
     or otherwise.  (See [`RestrictedStatus`](@ref StationXML.RestrictedStatus).)"
    restricted_status::M{RestrictedStatus} = missing
    "A code used for display or association, alternate to the SEED-compliant code."
    alternate_code::M{String} = missing
    "A previously used code if different from the current code."
    historical_code::M{String} = missing
    # end of BaseNodeType fields
    "Latitude coordinate of station in degrees."
    latitude::Latitude
    "Longitude coordinate of station in degrees."
    longitude::Longitude
    "Elevation of station above the local reference level in m."
    elevation::Distance
    "These fields describe the location of the station
     using geopolitical entities (country, city, etc.)."
    site::Site
    "Elevation of the water surface in meters for underwater sites, where 0 is sea level."
    water_level::M{Float} = missing
    "Type of vault, e.g. WWSSN, tunnel, transportable array, etc."
    vault::M{String} = missing
    "Type of rock and/or geologic formation."
    geology::M{String} = missing
    "Equipment used by all channels at a station."
    equipment::Vector{Equipment} = Equipment[]
    "An operating agency and associated contact persons. If
     there multiple operators, each one should be encapsulated within an
     Operator tag. Since the Contact element is a generic type that
     represents any contact person, it also has its own optional Agency
     element."
    operator::Vector{Operator} = Operator[]
    "Date and time (UTC) when the station was first installed.
     **Note that this field is mandatory in StationXML v1.0, but optional
     in StationXML v1.1.**"
    creation_date::M{DateTime} = missing
    # FIXME: Test whether blank values need to be handled as a M{DateTime} field
    "Date and time (UTC) when the station was terminated or
     will be terminated. A blank value should be assumed to mean that the
     station is still active."
    termination_date::M{DateTime} = missing
    "Total number of channels recorded at this station."
    total_number_channels::M{Int} = missing
    "Number of channels recorded at this station and
	 selected by the query that produced this document."
    selected_number_channels::M{Int} = missing
    "URI of any type of external report, such as IRIS data
     reports or dataless SEED volumes."
    external_reference::Vector{ExternalReference} = ExternalReference[]
    "Set of [`Channel`](@ref StationXML.Channel)s contained within this
     `Network` and returned as part of a request."
    channel::Vector{Channel} = Channel[]
end

attribute_fields(::Type{Station}) = (BASE_NODE_ATTRIBUTES...,)

"""
    Network

This type represents the Network layer, all station metadata is
contained within this element. The official name of the network or other descriptive
information can be included in the Description element. The Network can contain 0 or
more Stations.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Network
    # BaseNodeType fields
    description::M{String} = missing
    "Persistent identifier for this `Network`"
    identifier::Vector{Identifier} = Identifier[]
    comment::Vector{Comment} = Comment[]
    "A description of time series data availablility for this `Network`"
    data_availability::M{DataAvailability} = missing
    "Network's code (e.g., `\"IU\"`)."
    code::String
    start_date::M{DateTime} = missing
    end_date::M{DateTime} = missing
    "A data source identifier in URI form"
    source_id::M{String} = missing
    "Information on whether the data in this network are restricted, open,
     or otherwise.  (See [`RestrictedStatus`](@ref StationXML.RestrictedStatus).)"
    restricted_status::M{RestrictedStatus} = missing
    "A code used for display or association, alternate to the SEED-compliant code."
    alternate_code::M{String} = missing
    "A previously used code if different from the current code."
    historical_code::M{String} = missing
    # end of BaseNodeType fields
    "An operating agency and associated contact persons."
    operator::Vector{Operator} = Operator[]
    "The total number of stations contained in this
	 network, including inactive or terminated stations."
    total_number_stations::M{Int} = missing
    "The total number of stations in this network that
	 were selected by the query that produced this document, even if the
	 stations do not appear in the document. (This might happen if the
	 user wants a document that contains only information at
	 the Network level.)"
    selected_number_stations::M{Int} = missing
    "List of [`Station`](@ref StationXML.Station)s within a network included
     in this query."
    station::Vector{Station} = Station[]

    function Network(description, identifier, comment, data_availability, code,
            start_date, end_date,
            source_id, restricted_status, alternate_code, historical_code, operator,
            total_number_stations, selected_number_stations, station)
        if total_number_stations !== missing
            total_number_stations >= 0 ||
                throw(ArgumentError("total_number_stations must be 0 or more"))
        end
        if selected_number_stations !== missing
            selected_number_stations >= 0 ||
                throw(ArgumentError("selected_number_stations must be 0 or more"))
        end
        new(description, identifier, comment, data_availability, code, start_date, end_date,
            source_id, restricted_status, alternate_code, historical_code, operator,
            total_number_stations, selected_number_stations, station)
    end
end

attribute_fields(::Type{Network}) = (BASE_NODE_ATTRIBUTES...,)

"""
    FDSNStationXML

Top-level type for Station XML. Required field are `source` (network ID
of the institution sending the message), `created` (creation time of the document)
and one or more `Network`s containers.

This corresponds to the schema's RootType.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct FDSNStationXML
    "Network ID of the institution sending the message."
    source::String
    "Name of the institution sending this message."
    sender::M{String} = missing
    "Name of the software module that generated this document.
     [This differs from the StationXML schema because `module` is a reserved
     word in Julia.]"
    module_name::M{String} = missing
    "This is the address of the query that generated the document,
     or, if applicable, the address of the software that generated this document.
     Note that this string is not enforced to be a valid URI."
    module_uri::M{String} = missing
    "Date of creation of this StationXML record."
    created::DateTime
    "Set of [`Network`](@ref StationXML.Network)s contained within this
     document."
    network::Vector{Network} = Network[]
    "Version of the StationXML schema used by these data."
    schema_version::String
end

attribute_fields(::Type{FDSNStationXML}) = (:schema_version,)

"Types which should be compared using Base.==="
const COMPARABLE_TYPES = Union{Float64, String, DateTime}

# Hacky way to define == for all the types we create
for name in names(StationXML, all=true)
    name_str = String(name)
    if occursin(r"^[A-Z]", name_str) && !occursin(r"^BASE", name_str)
        name in (:M, :COMPARABLE_TYPES, :StationXML, :Numerator, :Denominator) && continue
        getfield(StationXML, name) isa DataType || continue
        @eval Base.:(==)(a::$name, b::$name)::Bool = a === b || local_equals(a, b)
    end
end

"""Local function to compare all types by each of their fields, apart from the types from
Base we use."""
local_equals(a::T, b::T) where {T<:COMPARABLE_TYPES} = a === b
local_equals(::Missing, ::Missing) = true
local_equals(::COMPARABLE_TYPES, ::Missing) = false
local_equals(::Missing, ::COMPARABLE_TYPES) = false
local_equals(a::T, b::T) where {T} =
    all(local_equals(va, vb) for (va, vb) in zip(all_values(a), all_values(b)))
local_equals(a::T1, b::T2) where {T1,T2} = false
local_equals(a::AbstractArray, b::AbstractArray)::Bool =
    size(a) == size(b) && all(local_equals(aa, bb) for (aa, bb) in zip(a,b))


text:
using sample
import other
include("ss")
include
include("another.jl")
include("subdir/yetanother.jl")

ext: jl

max_stars_repo_name: JuliaTagBot/deoplete-julia

max_issues_repo_name: JuliaTagBot/deoplete-julia

max_forks_repo_name: JuliaTagBot/deoplete-julia


text:
export ChanceConstraint

"""
Description:

    Chance constraint on the marginal q of the connected variable x, as
    ϵ ⩽ 1 - ∫_G q(x) dx, where G indicates a region. In other words,
    the probability mass of q is not allowed to overflow the region G
    by more than ϵ. Implementation according to (van de Laar et al.
    "Chance-Constrained Active Inference", MIT Neural Computation, 2021).

Interfaces:

    1. out

Construction:

    ChanceConstraint(out; G=(min,max), epsilon=epsilon, id=:my_node)
"""
mutable struct ChanceConstraint <: SoftFactor
    id::Symbol
    interfaces::Array{Interface,1}
    i::Dict{Symbol, Interface}

    G::Tuple # Range for G (1-D only)
    epsilon::Float64 # Threshold
    atol::Union{Float64, Nothing} # Tolerance for corrected belief

    function ChanceConstraint(out; G::Tuple, epsilon::Float64, atol=nothing, id=ForneyLab.generateId(ChanceConstraint))
        @ensureVariables(out)
        self = new(id, Vector{Interface}(undef, 1), Dict{Symbol,Interface}(), G, epsilon, atol)
        ForneyLab.addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

slug(::Type{ChanceConstraint}) = "P"

# A breaker message is required if interface is partnered with an expectation constraint
requiresBreaker(interface::Interface, partner_interface::Interface, partner_node::ChanceConstraint) = true

breakerParameters(interface::Interface, partner_interface::Interface, partner_node::ChanceConstraint) = (Message{GaussianMeanVariance, Univariate}, ()) # Univariate only

# Constraints do not contribute to average energy
averageEnergy(::Type{ChanceConstraint}, marg_out::ProbabilityDistribution) = 0.0


text:
module FractionalTimeDG

using ArgCheck
using OffsetArrays
import GaussQuadrature
import SpecialFunctions
import LinearAlgebra
import SparseArrays

export Store, coef_G, coef_K, coef_H0, coef_H1
export coef_H_uniform!, coef_Hn!
export FODEdG!, evaluate_pcwise_poly!, evaluate_pcwise_poly, reconstruction
export legendre_polys!
export SpatialGrid, spatial_L2_projection, FPDEDG, evaluate_pcwise_error,
       evaluate_fem1d_soln!

Γ(x) = SpecialFunctions.gamma(x)

struct Store{T<: AbstractFloat}
    # shared data for low-level routines
    # initialise with store = setup(α, rmax, Mmax)
    # or setup(α, rmax, Mmax, ppImax)
    α::T          # fraction diffusion exponent
    rmax::Integer # (max degree of piecewise polynomials) + 1
    Mmax::Integer # maximum number of Gauss points
    unitlegendre::Vector{Matrix{T}} # w(x) = 1 on [0,1]
    legendre::Vector{Matrix{T}}     # w(x) = 1 on [-1,1]
    jacobi1::Vector{Matrix{T}}      # w(x) = (1-x)^(α-1)
    jacobi2::Vector{Matrix{T}}      # w(x) = (1+x)^(α-1)
    jacobi3::Vector{Matrix{T}}      # w(x) = (1-x)^α
    jacobi4::Vector{Matrix{T}}      # w(x) = (1+x)^α
    jacobi5::Vector{Matrix{T}}      # w(x) = (1-x)(1+x)^(α-1)
    ppImax::Integer # max points per interval (second dimension of Ψ)
    Ψ::Matrix{T}
    dΨ::Matrix{T}
    A::Vector{T}
    B::Vector{T}
    C::Matrix{T}
end

include("parts/store.jl")
include("parts/legendre_utils.jl")
include("parts/coef_uniform.jl")
include("parts/coef.jl")
include("parts/scalar_fode.jl")
include("parts/FEM1D.jl")

end # module FractionalTimeDG


text:
module quality

using LinearAlgebra

export rankdata, qvalues , fractional_ranks

function rankdata( values_in::Array{Float64,1} , method::AbstractString="ordinal" )
    if method == "ordinal"
        return( sortperm(values_in) )
    else
        return( values_in )
    end
end

function fractional_ranks( values_in::Array{Float64,1} , method::AbstractString="ordinal" )
    if method == "ordinal"
        return( (rankdata(values_in,method).-0.5)./length(values_in) )
    else
        return( values_in )
    end
end

function qvalues( p_values_in::Array{Float64,1} ,
                  pi0::Float64 = 1.0 ,
		  method::AbstractString="ordinal" )
    """
    See my method with the same name in https://github.com/richardtjornhammar/impetuous/blob/master/src/impetuous/quantification.py
    """
    p_s = p_values_in
    qs_ = []
    m = length(p_s)
    ps = p_s
    frp_  = fractional_ranks( ps,method )
    ifrp_ = [];for (p,f) in zip(ps,frp_) append!(ifrp_, p<=f ? f : p ) end
    #ifrp_ = [ ( (p<=f)*f + p*(p>f) ) for p,f in zip(ps,frp_) ]
    for ip in 1:length(ps)
        p_ = ps[ ip ] ; f_ = frp_[ip]
        q_ = pi0 * p_ / ifrp_[ip]
        append!(qs_, (q_,p_) )
    end
    return(qs_)
end

end

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
function parsesigma(confidence_rate_symbol::Symbol = :sigma1)
    if confidence_rate_symbol == :sigma1
        n = 1
    elseif confidence_rate_symbol == :sigma2
        n = 2
    elseif confidence_rate_symbol == :sigma3
        n = 3
    else
        m = match( r"^sigma(\d+)$", string(confidence_rate_symbol) )
        if m == nothing
            n = 0
        end
        dstr = m.captured[1]
        n = try
            parseint(dstr)
        catch ex
            0
        end
    end
    if n <= 0
        throw(ErrorException("invalid parameter confidence_rate should be :sigmaN, where N is a positive integer"))
    end
    return n
end


text:
"""
    Youbot <: Environment

    Youbot robot designed by Kuka (URDF and mesh files from https://github.com/mas-group/youbot_description)
"""
struct Youbot end

using Dojo

# path = "DojoEnvironments/src/youbot/deps/youbot_arm_only.urdf"
path = "DojoEnvironments/src/youbot/deps/youbot_arm_only_fixed_gripper.urdf"
# path = "DojoEnvironments/src/youbot/deps/youbot_base_only.urdf"
# path = "DojoEnvironments/src/youbot/deps/youbot.urdf"
# path = "DojoEnvironments/src/youbot/deps/youbot_dual_arm.urdf"
# path = "DojoEnvironments/src/youbot/deps/youbot_with_cam3d.urdf"
mech = Mechanism(path; gravity = -9.81, timestep = 0.01)

joint1 = get_joint(mech,:arm_joint_1)
joint2 = get_joint(mech,:arm_joint_2)
joint3 = get_joint(mech,:arm_joint_3)
joint4 = get_joint(mech,:arm_joint_4)
joint5 = get_joint(mech,:arm_joint_5)

set_minimal_coordinates!(mech,joint1,[0])
set_minimal_coordinates!(mech,joint2,[0])
set_minimal_coordinates!(mech,joint3,[0])
set_minimal_coordinates!(mech,joint4,[0])
set_minimal_coordinates!(mech,joint5,[0])

storage = simulate!(mech, 5.0, record = true)
visualize(mech, storage)


text:
using Documenter, DrugInteractions

makedocs(sitename = "DrugInteractions Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
"""
# Description

A type that contains the results of translating an expression into a form
that can be used to evaluate that expression as an anonymous function in any
of the supported forms.

# Examples

```
julia> Expression(
    :x,
    :(a + b),
    :(x = a + b),
    (:a, :b),
     t -> t[1] + t[2],
    (a, b) -> a + b,
    (a, b) -> a + b,
)
```
"""
struct Expression
    alias::Symbol
    body::Any
    raw_form::Any
    input_columns::Tuple
    tuple_form::Any
    broadcast_form::Any
    vector_form::Any
end

# TODO: Base.show
function Base.show(io::IO, x::Expression)
    @printf(
        """
        == Expression object ==

        alias: %s
        body: %s
        raw_form: %s
        input_columns: %s
        tuple_form: %s
        broadcast_form: %s
        vector_form: %s""",
        x.alias,
        x.body,
        x.raw_form,
        x.input_columns,
        x.tuple_form,
        x.broadcast_form,
        x.vector_form,
    )
    nothing
end

function Base.print(io::IO, x::Expression)
    print(io, String(Symbol(x.raw_form)))
end

"""
# Description

...

# Arguments

1. ...

# Return Values

1. ...

# Examples

```
julia>

```
"""
# TODO: Document this
function expression_expr(@nospecialize(e::Any))
    # TODO: Support passes as an argument.

    alias, body = get_alias(e)

    column_names = find_column_names(body)

    quote
        Expression(
            $(QuoteNode(alias)),
            $(QuoteNode(body)),
            $(QuoteNode(e)),
            $(column_names_tuple_expr(column_names)),
            $(tuple_form(body, column_names)),
            $(broadcast_form(body, column_names)),
            $(vector_form(body, column_names)),
        )
    end
end

"""
# Description

...

# Arguments

1. ...

# Return Values

1. ...

# Examples

```
julia>

```
"""
# TODO: Document this
macro expression(e::Any)
    expression_expr(e)
end


text:
#! /usr/bin/env julia

using Documenter

using GeometryBasics

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

makedocs(format=Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true"),
         sitename="GeometryBasics.jl",
         pages=[
                "index.md",
                "primitives.md",
                "rectangles.md",
                "polygons.md",
                "meshes.md",
                "decomposition.md",
                "metadata.md",
                "api.md"
               ],
         modules=[GeometryBasics])

deploydocs(repo="github.com/JuliaGeometry/GeometryBasics.jl.git", push_preview=true)


text:
# CUR decomposition sketch
# LinearTimeCUR - Drineas, et  al. "Fast monte carlo algorithms for matrices iii: computing a compressed approximate matrix decomposition", SIAM Journal of Computing, 2006.
#
# Danny Perry (dperry@cs.utah.edu)
# April 2015

using Debug

# CUR decomposition (rank-k)
# A - the matrix to approximate
# fro_norm - frobenius norm of A
# num_cols - number of colums to sample for C (and U)
# num_rows - number of rows to sample for R (and U)
# k - rank k to approximate
#
# returns C,U,R matrices
@debug function LinearTimeCUR(A,fro_norm, num_cols,num_rows,k)
	n = size(A,1)
	d = size(A,2)

	c = num_cols
	r = num_rows

	C = zeros(n,c)
	U = zeros(c,r)
	R = zeros(r,d)

	for i=1:d
		pc = norm(A[:,i])^2 / fro_norm^2
		scale = 1/sqrt(c*pc)
		for j=1:c
			if rand() <= pc
				C[:,j] = scale * A[:,i]
			end
		end
	end
	CU,CS,CV = svd(C'*C)
	k = min(k, rank(diagm(CS)))
	for i=1:n
		pc = norm(A[i,:])^2 / fro_norm^2
		scale = 1/sqrt(r*pc)
		for j=1:r
			if rand() <= pc
				R[j,:] = scale * A[i,:]
				U[:,j] = scale * C[i,:]
			end
		end
	end
	U = CV[:,1:k] * pinv(diagm(CS[1:k])) * CU[:,1:k]' * U

	C,U,R
end

# This uses a decomposition to choose samples based on the leverage scores
# Uk,Vk - the first k columns of the matrix U and V, such that  A = USV'
@debug function LeverageCUR(A,fro_norm,Uk,Vk, num_cols,num_rows,k)
	n = size(A,1)
	d = size(A,2)

	c = num_cols
	r = num_rows

	C = zeros(n,c)
	U = zeros(c,r)
	R = zeros(r,d)

	col_lev = sum(Vk.^2,2) ./ k
	row_lev = sum(Uk.^2,2) ./ k

	for i=1:d
		pc = col_lev[i]
		scale = 1/sqrt(c*pc)
		for j=1:c
			if rand() <= pc
				C[:,j] = scale * A[:,i]
			end
		end
	end
	CU,CS,CV = svd(C'*C)
	k = min(k, rank(diagm(CS)))
	for i=1:n
		pc = row_lev[i]
		scale = 1/sqrt(r*pc)
		for j=1:r
			if rand() <= pc
				R[j,:] = scale * A[i,:]
				U[:,j] = scale * C[i,:]
			end
		end
	end
	U = CV[:,1:k] * pinv(diagm(CS[1:k])) * CU[:,1:k]' * U

	C,U,R
end


text:
using LinearAlgebra
using Memoize

"""
    pauli_vector(x, y, z)

Assemble the "Pauli vector" matrix.
"""
pauli_vector(x, y, z) = ComplexF64[z x - im * y; x + im * y -z]

"""
    binary_digits(x::Integer, M::Integer)

Return `M` bits of the binary representation of the integer `x` as array (least significant bit at first index).
"""
function binary_digits(M::Integer, x::Integer)
    binary_digits!(BitArray(undef, M), x)
end

"""
    binary_digits!(m::BitArray{1}, x::Integer)

Fill `m` with the binary representation of the integer `x` (least significant bit at first index).
"""
function binary_digits!(m::BitArray{1}, x::Integer=0)
    for i in 1:length(m)
        m[i] = x & 1
        x >>= 1
    end
    return m
end

"""
    binary_to_int(m::BitArray{1}, s=0)

converts BitArray `m` into Integer
"""
function binary_to_int(m::BitArray{1}, x=0)
    v = 1
    for i in Base.view(m, :)
        x += v*i
        v <<= 1
    end
    x
end

"""
    quaternary_digits(x::Integer, M::Integer)

Return `M` base-4 digits of the quaternary representation of the integer `x` as array (least significant digit at first index).
"""
function quaternary_digits(M::Integer, x::Integer)
    quaternary_digits!(Vector{Int}(undef, M), x)
end

"""
    quaternary_digits!(m::Vector{Int}, x::Integer)

Fill `m` with the quaternary representation of the integer `x` (least significant digit at first index).
"""
function quaternary_digits!(m::Vector{Int}, x::Integer)
    for i in 1:length(m)
        m[i] = x & 3
        x >>= 2
    end
    return m
end

"""
    intlog2(x::Integer)

Compute integer base-2 logarithm. Rounds to floor if x is not a power of 2.
"""
function intlog2(x::Integer)
    x == 0 && error("Logarithm of 0 is undefined")
    x < 0 && error("Logarithm of negative number is undefined")
    ret::Int = 0
    while x > 1
        x >>= 1
        ret += 1
    end
    ret
end

"""
    sliced_index(idx, targetwires::Tuple, N::Int)

Construct sliced index via target map (utility function).
"""
function sliced_index(idx, targetwires::Tuple, N::Int)
    islice = Vector{Any}(fill(Colon(), N))
    for k in 1:length(targetwires)
        islice[targetwires[k]] = idx[k] + 1
    end
    return islice
end

"""
    gramm_schmidt(a::Array{Float64})

Compute orthogonalized real vectors.
"""
function gramm_schmidt!(a::Array{Float64})
    num_vectors = size(a,2)

    a[:, end] = a[:, end] ./ norm(a[:, end])

    for i in num_vectors-1:-1:1
        for j in i+1:num_vectors
            a[:, i] = a[:, i] - (dot(a[:, i], a[:, j]) / dot(a[:, j], a[:, j])) * a[:, j]
        end
        a[:, i] = a[:, i] ./ norm(a[:, i])
    end

    return a
end

"""
    gramm_schmidt(a::Array{ComplexF64})

Compute orthogonalized complex vectors.
"""
function gramm_schmidt!(a::Array{ComplexF64})
    num_vectors = size(a,2)

    a[:, end] = a[:, end] ./ norm(a[:, end])

    for i in num_vectors-1:-1:1
        for j in num_vectors:-1:i+1
            a[:, i] = a[:, i] - conj(dot(a[:, i], a[:, j]) / dot(a[:, j], a[:, j])) * a[:, j]
        end
        a[:, i] = a[:, i] ./ norm(a[:, i])
    end
    return a
end

"""
    rdm(N, iwire, ψ, χ, d=2)

Compute the reduced density matrix ``tr_B[|ψ⟩⟨χ|]``, where the trace runs over
the subsystem complementary to the qubits specified by `iwire`.
"""
@views function rdm(N::Integer, iwire::NTuple{M,<:Integer}, ψ::Vector{T}, χ::Vector{T}, d::Int=2) where {M,T}
    M ≥ 1 || error("Need at least one wire to act on.")
    M ≤ N || error("Number of gate wires cannot be larger than total number of wires.")
    length(unique(iwire)) == M || error("Wire indices must be unique.")
    minimum(iwire) ≥ 1 || error("Wire index cannot be smaller than 1.")
    maximum(iwire) ≤ N || error("Wire index cannot be larger than total number of wires.")
    length(ψ) == d^N || error("Input vector 'ψ' has wrong length.")
    length(χ) == d^N || error("Input vector 'χ' has wrong length.")

    # convert to array
    iwire = collect(iwire)
    # complementary wires
    iwcompl = setdiff(1:N, iwire)
    @assert length(iwire) + length(iwcompl) == N

    ρ = zeros(eltype(ψ), d^M, d^M)
    # # Note: following the ordering convention of `kron` here, i.e.,
    # # last qubit corresponds to fastest varying index
    strides = [d^(j-1) for j in 1:N]
    wstrides = strides[iwire]
    cstrides = strides[iwcompl]

    if d == 2
        iw = BitArray{1}(undef, M)
        jw = BitArray{1}(undef, M)
        kw = BitArray{1}(undef, N - M)
        # TODO: optimize memory access pattern
        for k in 0:2^(N-M)-1
            binary_digits!(kw, k)
            koffset = dot(kw, cstrides)
            for i in 1:2^M
                binary_digits!(iw, i - 1)
                for j in 1:2^M
                    binary_digits!(jw, j - 1)
                    rowind = koffset + dot(iw, wstrides) + 1
                    colind = koffset + dot(jw, wstrides) + 1
                    @inbounds ρ[i, j] += ψ[rowind] * conj(χ[colind])
                end
            end
        end
    elseif d == 4
        iw = Vector{Int}(undef, M)
        jw = Vector{Int}(undef, M)
        kw = Vector{Int}(undef, N - M)
        # TODO: optimize memory access pattern
        for k in 0:4^(N-M)-1
            quaternary_digits!(kw, k)
            koffset = dot(kw, cstrides)
            for i in 1:4^M
                quaternary_digits!(iw, i - 1)
                for j in 1:4^M
                    quaternary_digits!(jw, j - 1)
                    rowind = koffset + dot(iw, wstrides) + 1
                    colind = koffset + dot(jw, wstrides) + 1
                    @inbounds ρ[i, j] += ψ[rowind] * conj(χ[colind])
                end
            end
        end
    else
        error("d = $d not supported yet.")
    end

    return ρ
end


text:
horizontal_position = 0
vertical_position = 0

commands = readlines("../inputs/day02.txt")

for command in commands
    direction = split(command)[1]
    magnitude = parse(Int, split(command)[2])

    if direction == "up"
        global vertical_position += magnitude
    elseif direction == "down"
        global vertical_position -= magnitude
    elseif direction == "forward"
        global horizontal_position += magnitude
    end
end

depth = -vertical_position

println("Day 02 Part 1 Solution: ", horizontal_position * depth)

horizontal_position = 0
vertical_position = 0
aim = 0

for command in commands
    direction = split(command)[1]
    magnitude = parse(Int, split(command)[2])

    if direction == "up"
        global aim += magnitude
    elseif direction == "down"
        global aim -= magnitude
    elseif direction == "forward"
        global horizontal_position += magnitude
        global vertical_position += aim * magnitude
    end
end

depth = -vertical_position

println("Day 02 Part 2 Solution: ", horizontal_position * depth)


text:
#-----------------------------------------------------------------------
# JuMPeR  --  JuMP Extension for Robust Optimization
# http://github.com/IainNZ/JuMPeR.jl
#-----------------------------------------------------------------------
# Copyright (c) 2016: Iain Dunning
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#-----------------------------------------------------------------------
# example/cap_budget.jl
# Apply the method from "Multistage Robust Mixed Integer Optimization
# with Adaptive Partitions" by Bertsimas and Dunning to the two-stage
# capital budgeting problem first described in Hanasusanto et. al.'s
# paper "Two-stage robust integer programming".
#
#   max  ∑ᵢ profitᵢ(ξ) [xᵢ + Θ yᵢ(ξ)]
#    st  ∑ᵢ costᵢ(ξ) [xᵢ + yᵢ(ξ)] ≤ B
#        xᵢ, yᵢ(ξ) ∈ {0,1}
#
# Requires a mixed-integer linear optimization problem solver.
#-----------------------------------------------------------------------

using JuMP, JuMPeR, LinearAlgebra, Random, Gurobi

"""
    TreeScenario

Stores the values of "active" uncertain parameters, as well as the
associated tree structure described in the paper.
"""
mutable struct TreeScenario
    ξ::Vector{Float64}
    parent
    children::Vector
end
is_leaf(t::TreeScenario) = isempty(t.children)

"""
    generate_data(seed, N)

Generates a set of data for an instance of size `N` using the distributions
described in the paper by Hanasusanto, et. al., where `seed` is the random
seed to use to help with replication.
"""
function generate_data(seed, N)
    Random.seed!(seed)
    c0 = rand(N)*10  # Nominal costs
    r0 = c0 ./ 5     # Nominal profit
    B  = sum(c0)/2   # Budget
    θ = 0.8          # Profit discount factor for late decisions
    Φ = rand(N,4)    # Factor matrices ...
    Ψ = rand(N,4)    # ... for the uncertainty set
    for i in 1:N
        Φ[i,:] ./= sum(Φ[i,:])
        Ψ[i,:] ./= sum(Ψ[i,:])
    end
    return θ, B, Φ, Ψ, c0, r0
end

"""
    solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenarios)

Solve the two-stage problem of size `N` with one partition of the uncertainty
set for every leaf scenario in the `scenario_tree`. At optimality, grows the
tree given the new scenarios obtained.
"""
function solve_partitioned_problem(N::Int, θ::Float64, B::Float64,
                                   Φ::Matrix{Float64}, Ψ::Matrix{Float64},
                                   c0::Vector{Float64}, r0::Vector{Float64},
                                   scenario_tree::Vector{TreeScenario})
    # scenario_tree is a vector containing every scenario in the tree
    # We will have one partition for every leaf scenario in the tree
    leaf_scenarios = filter(is_leaf, scenario_tree)
    P = length(leaf_scenarios)  # Number of partitions
    # Initialize the RO model
    rm = RobustModel(solver = GurobiSolver())
    # Decision variables: whether to pursue a project or not
    # First stage, here-and-now decision
    @variable(rm, x[1:N], Bin)
    # Second stage, wait-and-see decision
    # One set of variables per partition
    @variable(rm, y[1:P,1:N], Bin)
    # Define the uncertain parameters, which belong to a hypercube
    @uncertain(rm, -1 <= ξ[1:4] <= 1)
    # From the uncertain parameters we obtain the cost and profit vectors
    cost   = (1 + 0.5*Φ*ξ) .* c0
    profit = (1 + 0.5*Ψ*ξ) .* r0
    # The objective function will be the minimum of the objective function
    # across all the partitions. Put a default upper bound, just so we don't
    # start off unbounded if we are using a cutting plane method.
    @variable(rm, obj <= 10)
    @objective(rm, Max, obj)
    # For each partition...
    profit_con_refs = []
    budget_con_refs = []
    for p in 1:P
        # Define the partition uncertainty set. We define multiple hyperplanes
        # by walking up the scenario tree from this leaf.
        us = JuMPeR.BasicUncertaintySet()
        current_scenario = leaf_scenarios[p]
        parent_scenario = current_scenario.parent
        # We keep going until we hit the root of the tree, which is a scenario
        # that has no parent
        while parent_scenario != nothing
            for sibling_scenario in parent_scenario.children
                if current_scenario == sibling_scenario
                    continue  # Don't partition against ourself!
                end
                ξ_sub = sibling_scenario.ξ - current_scenario.ξ
                ξ_add = sibling_scenario.ξ + current_scenario.ξ
                @constraint(us, dot(ξ_sub, ξ) <= dot(ξ_sub,ξ_add)/2)
            end
            # Move up the scenario tree
            current_scenario  = parent_scenario
            parent_scenario = current_scenario.parent
        end
        # Constrain objective function for this partition
        profit_con_ref =
            @constraint(rm, obj <= sum{profit[i] * (x[i] + θ*y[p,i]), i in 1:N},
                                uncset=us)
        push!(profit_con_refs, profit_con_ref)
        # Constrain spending for this partition
        budget_con_ref =
            @constraint(rm, sum{cost[i] * (x[i] + y[p,i]), i in 1:N} <= B,
                                uncset=us)
        push!(budget_con_refs, budget_con_ref)
        # Can only choose each project once
        for i in 1:N
            @constraint(rm, x[i] + y[p,i] <= 1)
        end
    end
    # Solve, will use reformulation. We pass disable_cuts=true because
    # we want to be able to use solvers like CBC that don't support lazy
    # constraint callbacks.
    solve(rm, disable_cuts=true, active_scenarios=true)
    # Extend the scenario tree
    for p in 1:P
        # Extract the active uncertain parameter values
        profit_scen = getscenario(profit_con_refs[p])
        profit_scen_ξ = [uncvalue(profit_scen, ξ[i]) for i in 1:4]
        # Create a new child in the tree under this leaf
        profit_child = TreeScenario(profit_scen_ξ, leaf_scenarios[p], [])
        # Same for budget
        budget_scen = getscenario(budget_con_refs[p])
        budget_scen_ξ = [uncvalue(budget_scen, ξ[i]) for i in 1:4]
        budget_child = TreeScenario(budget_scen_ξ, leaf_scenarios[p], [])
        # Add to the tree
        push!(leaf_scenarios[p].children, profit_child)
        push!(leaf_scenarios[p].children, budget_child)
        push!(scenario_tree, profit_child)
        push!(scenario_tree, budget_child)
    end
    # Return the objective function value and the first-stage solution
    getobjectivevalue(rm), getvalue(x)
end

"""
    solve_problem(seed, N)

Solve the problem for the parameters corresponding to the given random seed
and instance size `N`. Prints results are each of five iterations.
"""
function solve_problem(seed, N)
    θ, B, Φ, Ψ, c0, r0 = generate_data(seed, N)
    # Start with no partitions (i.e., one scenario)
    scenario_tree = [ TreeScenario(zeros(4),nothing,[]) ]
    # Iteration 1
    z_iter1, x_iter1 = solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenario_tree)
    # Iteration 2
    z_iter2, x_iter2 = solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenario_tree)
    # Iteration 3
    z_iter3, x_iter3 = solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenario_tree)
    # Iteration 4
    z_iter4, x_iter4 = solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenario_tree)
    # Iteration 5
    z_iter5, x_iter5 = solve_partitioned_problem(N, θ, B, Φ, Ψ, c0, r0, scenario_tree)
    @show z_iter1
    println(round.(x_iter1))
    @show z_iter2
    println(round.(x_iter2))
    @show z_iter3
    println(round.(x_iter3))
    @show z_iter4
    println(round.(x_iter4))
    @show z_iter5
    println(round.(x_iter5))
    @show length(scenario_tree)
end

# By default, solve instance of size 5
solve_problem(106, 5)


text:
module Preferences

if VERSION < v"1.6.0-DEV"
    error("Preferences.jl can only be used on Julia v1.6+!")
end

using TOML
using Base: UUID, TOMLCache

export load_preference, @load_preference,
       has_preference, @has_preference,
       set_preferences!, @set_preferences!,
       delete_preferences!, @delete_preferences!

include("utils.jl")

"""
    load_preference(uuid_or_module, key, default = nothing)

Load a particular preference from the `Preferences.toml` file, shallowly merging keys
as it walks the hierarchy of load paths, loading preferences from all environments that
list the given UUID as a direct dependency.

Most users should use the `@load_preference` convenience macro which auto-determines the
calling `Module`.
"""
function load_preference(uuid::UUID, key::String, default = nothing)
    # Re-use definition in `base/loading.jl` so as to not repeat code.
    d = Base.get_preferences(uuid)
    if currently_compiling()
        Base.record_compiletime_preference(uuid, key)
    end
    # Drop any nested `__clear__` keys:
    function drop_clears(data::Dict)
        delete!(data, "__clear__")
        for (k, v) in data
            if isa(v, Dict)
                drop_clears(v)
            end
        end
        return data
    end
    drop_clears(x) = x

    return drop_clears(get(d, key, default))
end
function load_preference(m::Module, key::String, default = nothing)
    return load_preference(get_uuid(m), key, default)
end

"""
    @load_preference(key)

Convenience macro to call `load_preference()` for the current package.
"""
macro load_preference(key, default = nothing)
    return quote
        load_preference($(esc(get_uuid(__module__))), $(esc(key)), $(esc(default)))
    end
end

"""
    has_preference(uuid_or_module, key, default = nothing)

Return `true` if the particular preference is found, and `false` otherwise.

See the `has_preference` docstring for more details.
"""
function has_preference(uuid::UUID, key::String, default = nothing)
    value = load_preference(uuid, key, nothing)
    return !(value isa Nothing)
end
function has_preference(m::Module, key::String, default = nothing)
    return has_preference(get_uuid(m), key, default)
end

"""
    @has_preference(key)

Convenience macro to call `has_preference()` for the current package.
"""
macro has_preference(key, default = nothing)
    return quote
        has_preference($(esc(get_uuid(__module__))), $(esc(key)), $(esc(default)))
    end
end

"""
    process_sentinel_values!(prefs::Dict)

Recursively search for preference values that end in `nothing` or `missing` leaves,
which we handle specially, see the `set_preferences!()` docstring for more detail.
"""
function process_sentinel_values!(prefs::Dict)
    # Need to widen `prefs` so that when we try to assign to `__clear__` below,
    # we don't error due to a too-narrow type on `prefs`
    prefs = Base._typeddict(prefs, Dict{String,Vector{String}}())

    clear_keys = get(prefs, "__clear__", String[])
    for k in collect(keys(prefs))
        if prefs[k] === nothing
            # If this should add `k` to the `__clear__` list, do so, then remove `k`
            push!(clear_keys, k)
            delete!(prefs, k)
        else
            # `k` is not nothing, so drop it from `clear_keys`
            filter!(x -> x != k, clear_keys)
            if prefs[k] === missing
                # If this should clear out the mapping for `k`, do so
                delete!(prefs, k)
            elseif isa(prefs[k], Dict)
                # Recurse for nested dictionaries
                prefs[k] = process_sentinel_values!(prefs[k])
            end
        end
    end
    # Store the updated list of clear_keys
    if !isempty(clear_keys)
        prefs["__clear__"] = collect(Set(clear_keys))
    else
        delete!(prefs, "__clear__")
    end
    return prefs
end

# See the `set_preferences!()` docstring below for more details
function set_preferences!(target_toml::String, pkg_name::String, pairs::Pair{String,<:Any}...; force::Bool = false)
    # Load the old preferences in first, as we'll merge ours into whatever currently exists
    d = Dict{String,Any}()
    if isfile(target_toml)
        d = Base.parsed_toml(target_toml)
    end
    prefs = d
    if endswith(target_toml, "Project.toml")
        if !haskey(prefs, "preferences")
            prefs["preferences"] = Dict{String,Any}()
        end
        # If this is a `(Julia)Project.toml` file, we squirrel everything away under the
        # "preferences" key, while for a `Preferences.toml` file it sits at top-level.
        prefs = prefs["preferences"]
    end
    # Index into our package name
    if !haskey(prefs, pkg_name)
        prefs[pkg_name] = Dict{String,Any}()
    end
    # Set each preference, erroring unless `force` is set to `true`
    for (k, v) in pairs
        if !force && haskey(prefs[pkg_name], k) && (v === missing || prefs[pkg_name][k] != v)
            throw(ArgumentError("Cannot set preference '$(k)' to '$(v)' for $(pkg_name) in $(target_toml): preference already set to '$(prefs[pkg_name][k])'!"))
        end
        prefs[pkg_name][k] = v

        # Recursively scan for `nothing` and `missing` values that we need to handle specially
        prefs[pkg_name] = process_sentinel_values!(prefs[pkg_name])
    end
    open(target_toml, "w") do io
        TOML.print(io, d, sorted=true)
    end
    return nothing
end

"""
    set_preferences!(uuid_or_module, prefs::Pair{String,Any}...; export_prefs=false, force=false)

Sets a series of preferences for the given UUID/Module, identified by the pairs passed in
as `prefs`.  Preferences are loaded from `Project.toml` and `LocalPreferences.toml` files
on the load path, merging values together into a cohesive view, with preferences taking
precedence in `LOAD_PATH` order, just as package resolution does.  Preferences stored in
`Project.toml` files are considered "exported", as they are easily shared across package
installs, whereas the `LocalPreferences.toml` file is meant to represent local
preferences that are not typically shared.  `LocalPreferences.toml` settings override
`Project.toml` settings where appropriate.

After running `set_preferences!(uuid, "key" => value)`, a future invocation of
`load_preference(uuid, "key")` will generally result in `value`, with the exception of
the merging performed by `load_preference()` due to inheritance of preferences from
elements higher up in the `load_path()`.  To control this inheritance, there are two
special values that can be passed to `set_preferences!()`: `nothing` and `missing`.

* Passing `missing` as the value causes all mappings of the associated key to be removed
  from the current level of `LocalPreferences.toml` settings, allowing preferences set
  higher in the chain of preferences to pass through.  Use this value when you want to
  clear your settings but still inherit any higher settings for this key.

* Passing `nothing` as the value causes all mappings of the associated key to be removed
  from the current level of `LocalPreferences.toml` settings and blocks preferences set
  higher in the chain of preferences from passing through.  Internally, this adds the
  preference key to a `__clear__` list in the `LocalPreferences.toml` file, that will
  prevent any preferences from leaking through from higher environments.

Note that the behavior of `missing` and `nothing` is both similar (they both clear the
current settings) and diametrically opposed (one allows inheritance of preferences, the
other does not).  They can also be composed with a normal `set_preferences!()` call:

```julia
@set_preferences!("compiler_options" => nothing)
@set_preferences!("compiler_options" => Dict("CXXFLAGS" => "-g", LDFLAGS => "-ljulia"))
```

The above snippet first clears the `"compiler_options"` key of any inheriting influence,
then sets a preference option, which guarantees that future loading of that preference
will be exactly what was saved here.  If we wanted to re-enable inheritance from higher
up in the chain, we could do the same but passing `missing` first.

The `export_prefs` option determines whether the preferences being set should be stored
within `LocalPreferences.toml` or `Project.toml`.
"""
function set_preferences!(u::UUID, prefs::Pair{String,<:Any}...; export_prefs=false, kwargs...)
    # Find the first `Project.toml` that has this UUID as a direct dependency
    project_toml, pkg_name = find_first_project_with_uuid(u)
    if project_toml === nothing && pkg_name === nothing
        # If we couldn't find one, we're going to use `active_project()`
        project_toml = Base.active_project()

        # And we're going to need to add this UUID as an "extras" dependency:
        # We're going to assume you want to name this this dependency in the
        # same way as it's been loaded:
        pkg_uuid_matches = filter(d -> d.uuid == u, keys(Base.loaded_modules))
        if isempty(pkg_uuid_matches)
            error("Cannot set preferences of an unknown package that is not loaded!")
        end
        pkg_name = first(pkg_uuid_matches).name

        # Read in the project, add the deps, write it back out!
        project = Base.parsed_toml(project_toml)
        if !haskey(project, "extras")
            project["extras"] = Dict{String,Any}()
        end
        project["extras"][pkg_name] = string(u)
        open(project_toml, "w") do io
            TOML.print(io, project; sorted=true)
        end
    end

    # Finally, save the preferences out to either `Project.toml` or `Preferences.toml`
    # keyed under that `pkg_name`:
    target_toml = project_toml
    if !export_prefs
        target_toml = joinpath(dirname(project_toml), "LocalPreferences.toml")
    end
    return set_preferences!(target_toml, pkg_name, prefs...; kwargs...)
end
function set_preferences!(m::Module, prefs::Pair{String,<:Any}...; kwargs...)
    return set_preferences!(get_uuid(m), prefs...; kwargs...)
end

"""
    @set_preferences!(prefs...)

Convenience macro to call `set_preferences!()` for the current package.  Defaults to
setting `force=true`, since a package should have full control over itself, but not
so for setting the preferences in other packages, pending private dependencies.
"""
macro set_preferences!(prefs...)
    return quote
        set_preferences!($(esc(get_uuid(__module__))), $(esc(prefs...)), force=true)
    end
end

"""
    delete_preferences!(uuid_or_module, prefs::String...; block_inheritance::Bool = false, export_prefs=false, force=false)

Deletes a series of preferences for the given UUID/Module, identified by the
keys passed in as `prefs`.

See the docstring for `set_preferences!`for more details.
"""
function delete_preferences!(u::UUID, pref_keys::String...; block_inheritance::Bool = false, kwargs...)
    if block_inheritance
        return set_preferences!(u::UUID, [k => nothing for k in pref_keys]...; kwargs...)
    else
        return set_preferences!(u::UUID, [k => missing for k in pref_keys]...; kwargs...)
    end
end
function delete_preferences!(m::Module, pref_keys::String...; kwargs...)
    return delete_preferences!(get_uuid(m), pref_keys...; kwargs...)
end

"""
    @delete_preferences!(prefs...)

Convenience macro to call `delete_preferences!()` for the current package.  Defaults to
setting `force=true`, since a package should have full control over itself, but not
so for deleting the preferences in other packages, pending private dependencies.
"""
macro delete_preferences!(prefs...)
    return quote
        delete_preferences!($(esc(get_uuid(__module__))), $(esc(prefs...)), force=true)
    end
end

end # module Preferences


text:
function PlayGame(startingStr, endStep)
    startingNums = map(x->parse(Int, x),split(startingStr,","));

    memory = Dict{Int128, Int128}();
    currentStep::Int128 = 1;

    for num in startingNums
        memory[num] = currentStep;
        currentStep += 1;
    end

    lastNum = Int128(startingNums[end]);
    seenBefore::Int128 = -1;

    while currentStep <= endStep
        curNum::Int128 = 0;

        if seenBefore > 0
            curNum = currentStep - seenBefore - 1;
            seenBefore = -1;
        end
        #println("$(currentStep) $(curNum)");

        if haskey(memory, curNum)
            seenBefore = memory[curNum];
        end

        memory[curNum] = currentStep;
        currentStep += 1;
        lastNum = curNum;
    end

    println(lastNum);
    return lastNum;
end

PlayGame("0,3,6", 10);
PlayGame("1,3,2", 2020);

answer = PlayGame("9,6,0,10,18,2,1", 30000000)


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(Base.Broadcast.materialize),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(RecipesPipeline._scaled_adapted_grid),Tuple{Array{Function,1},Base.RefValue{Symbol},Base.RefValue{Symbol},Float64,Float64}}})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol,Any},AbstractArray{T,2} where T})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol,Any},Array{Function,1},Number,Number})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol,Any},RecipesPipeline.GroupBy,Any})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol,Any},Type{SliceIt},Any,Any,Any})
    Base.precompile(Tuple{typeof(RecipesPipeline._apply_type_recipe),Any,AbstractArray,Any})
    Base.precompile(Tuple{typeof(RecipesPipeline._apply_type_recipe),Any,Surface,Any})
    Base.precompile(Tuple{typeof(RecipesPipeline._extract_group_attributes),Array{String,1},Array{Float64,1}})
    Base.precompile(Tuple{typeof(RecipesPipeline._map_funcs),Function,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})
    Base.precompile(Tuple{typeof(RecipesPipeline._process_ribbon),Tuple{LinRange{Float64},LinRange{Float64}},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._process_seriesrecipe),Any,Any})
    Base.precompile(Tuple{typeof(RecipesPipeline._process_seriesrecipes!),Any,Any})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{AbstractArray{Float64,1},1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Array{Float64,1},1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Array{T,1} where T,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Float64,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Float64,2},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Function,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Int64,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Real,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{String,1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Array{Union{Missing, Int64},1},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Int64,Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),StepRange{Int64,Int64},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),Surface{Array{Int64,2}},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline._series_data_vector),UnitRange{Int64},Dict{Symbol,Any}})
    Base.precompile(Tuple{typeof(RecipesPipeline.filter_data),Array{Float64,1},Array{Int64,1}})
    Base.precompile(Tuple{typeof(RecipesPipeline.userrecipe_signature_string),Any,Vararg{Any,N} where N})
    Base.precompile(Tuple{typeof(RecipesPipeline.userrecipe_signature_string),Any})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{Array{Dates.DateTime,1},UnitRange{Int64},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{Array{Float64,1},Array{Float64,1},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{Array{String,1},Array{String,1},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Surface{Array{Int64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{DataType,RecipesPipeline.Formatted{Array{Int64,1}},UnitRange{Int64},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{RecipesPipeline.Formatted{Array{Int64,1}},UnitRange{Int64},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{RecipesPipeline.GroupBy,Array{Float64,1}}}}}}})
    Base.precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(string),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(typeof),Tuple{Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Surface{Array{Float64,2}}}}}}}})
    Base.precompile(Tuple{typeof(recipe_pipeline!),Any,Any,Any})
    Base.precompile(Tuple{typeof(unzip),Array{Tuple{Array{Float64,1},Array{Float64,1}},1}})
    Base.precompile(Tuple{typeof(unzip),Array{Tuple{Int64,Int64},1}})
    Base.precompile(Tuple{typeof(unzip),Array{Tuple{Int64,Real},1}})
    let fbody = try __lookup_kwbody__(which(RecipesPipeline._extract_group_attributes, (Array{String,1},Array{Float64,1},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Function,typeof(RecipesPipeline._extract_group_attributes),Array{String,1},Array{Float64,1},))
        end
    end
end


text:
function _add_relations_from_subfield(mL::NfToNfMor; use_aut = true, redo = false, bound::Int = -1)
  L = codomain(mL)
  K = domain(mL)
  OK = lll(maximal_order(L))
  c = create_ctx(OK, use_aut = use_aut, redo = redo, bound = bound)
  U = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(OK)
  set_attribute!(OK, :UnitGrpCtx => U)
  set_attribute!(OK, :ClassGrpCtx => c)

  lp = Set{NfOrdIdl}()
  for p in c.FB.ideals
    push!(lp, intersect_prime(mL, p))
  end
  lp1 = collect(lp)

  vals_subfield = val_from_subfield(c.FB, mL, lp1)
  cache_mL = Dict{nf_elem, nf_elem}()
  @vprint :ClassGroup 1 "Computing S-units of totally real subfield\n"
  S, mS = Hecke.sunit_group_fac_elem(lp1)
  @vprint :ClassGroup 1 "Embedding S-units of totally real subfield\n"
  for i = 1:ngens(S)
    @vprint :ClassGroup 1 "Embedding S-units $i/$(ngens(S))\n"
    sup = Dict{NfOrdIdl, fmpz}((mS.idl[i], v) for (i, v) in mS.valuations[i])
    u = compact_presentation(mS(S[i]), 2, decom = sup)
    if iszero(mS.valuations[i])
      if is_torsion_unit(u)[1]
        continue
      end
      add_unit!(U, u)
    else
      img_u = FacElem(Dict{nf_elem, fmpz}((_embed(mL, cache_mL, x), v) for (x,v) = u.fac if !iszero(v)))
      valofnewelement = mul(mS.valuations[i], vals_subfield)
      Hecke.class_group_add_relation(c, u, valofnewelement, add_orbit = false, always = false)
    end
  end
  c.finished = false
  d = root(abs(discriminant(OK)), 2)
  @show c.expect = class_group_expected(c, 100)
  class_group_via_lll(c)
  return nothing
end

function _embed(mL, D, a)
  b = get(D, a, nothing)
  if b === nothing
    c = mL(a)
    D[a] = c
    return c
  else
     return b
  end
end

function val_from_subfield(FB, mk, s)
  S = FB.ideals
  ZK = order(S[1])

  z = zero_matrix(SMat, FlintZZ, 0, length(S))

  zk = order(s[1])

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    v = Tuple{Int, fmpz}[]
    P = s[l]
    genofsl = elem_in_nf(mk(P.gen_two.elem_in_nf))
    pmin = minimum(P, copy = false)
    # Compute the number of ideals
    # Assume that L/K and L/Q are normal
    rele = divexact(ramification_index((FB.fb[pmin].lp)[1][2]), ramification_index(P))
    relf = divexact(degree((FB.fb[pmin].lp)[1][2]), degree(P))
    @assert mod(reldeg, rele * relf) == 0
    numb_ideal = divexact(reldeg, rele * relf)
    found = 0
    for k in 1:length(S)
      Q = S[k]
      if minimum(Q, copy = false) == pmin
        if genofsl in Q
          found += 1
          @assert mod(ramification_index(Q), ramification_index(s[l])) == 0
          push!(v, (k, divexact(ramification_index(Q), ramification_index(s[l]))))
        end
      end
    end
    sort!(v, by = x -> x[1])
    push!(z, sparse_row(FlintZZ, v))
  end
  return z
end

function class_group_cm(OK::NfOrd; redo = false, use_aut = true, bound::Int = Int(ceil((log(abs(discriminant(OK)))^2)*0.3)))
  K = nf(OK)
  O = lll(OK)
  fl, conj = is_cm_field(nf(O))
  L, mL = fixed_field1(K, [conj])
  _add_relations_from_subfield(mL, use_aut = use_aut, redo = redo, bound = bound)
  c, U, b = _class_unit_group(O, use_aut = true)
  return class_group(c, OK)
end

function create_ctx(OK::NfOrd; bound::Int = -1, method::Int = 3, large::Int = 1000, redo::Bool = false, use_aut::Bool = false)
  if !redo
    c = get_attribute(OK, :ClassGrpCtx)
    if c !== nothing
      return c::ClassGrpCtx{SMat{fmpz}}
    end
  end

  if bound == -1
    bound = Int(ceil(log(abs(discriminant(OK)))^2*0.3))
    (bound == 0) && (bound = 1)
  end

  c = class_group_init(OK, bound, complete = false, use_aut = use_aut)::ClassGrpCtx{SMat{fmpz}}
  @assert order(c) === OK
  return c
end


text:
using Test, Distributed

addprocs(4)
@everywhere begin
    using SimulatedAnnealing, ParallelDataTransfer
    # current, result, options = State(), Result(), Options(Ns=5)
end
# println(result)

current, result, options = SimulatedAnnealing.init(4)
options.print_status = true
@time sa!(current, result, options)

@testset verbose = true "All" begin

# Test collect_results!
# @testset verbose = true "collect_results" begin
#     current.x = [2.5,2.5,2.5,2.5]; current.f = -options.f(current.x)
#     result.fopt = current.f; result.xopt = current.x
#     @everywhere [2] current.x, current.f = [2,2,2,2], -options.f([2,2,2,2])
#     @everywhere [3] current.x, current.f = [3,3,3,3], -options.f([3,3,3,3])
#     @everywhere [4] current.x, current.f = [4,4,4,4], -options.f([4,4,4,4])
#     @everywhere [5] current.x, current.f = [5,5,5,5], -options.f([5,5,5,5])

#     SimulatedAnnealing.collect_results!(current, result)
#     @test all(current.x .== [2,2,2,2])
#

end


text:
# Sample script for stream tracing from a given starting point.
#
# Hongyang Zhou, hyzhou@umich.edu

using PyPlot, FieldTracer, Vlasiator
using Vlasiator: RE # Earth radius, [m]

file = "bulk.0000501.vlsv"
nameρ = "rho"
nameV = "rho_v"

meta = load(file)

pcolormesh(meta, nameρ)

v = readvariable(meta, nameV)
vx = reshape(v[1,:], meta.ncells[1], meta.ncells[2])
vy = reshape(v[2,:], meta.ncells[1], meta.ncells[2])
# tracing starting point
xstart, ystart = 12RE, 0RE
# regular Cartesian mesh
x = range(meta.coordmin[1], meta.coordmax[1], length=meta.ncells[1])
y = range(meta.coordmin[2], meta.coordmax[2], length=meta.ncells[2])

# RK4 scheme by default
x1, y1 = trace2d(vx, vy, xstart, ystart, x, y;
   ds=0.5, maxstep=3000, gridType="ndgrid")
x1 ./= RE
y1 ./= RE

plot(x1, y1)

streamplot(meta, nameV, comp="xy", color="w", density=1.0)


text:
using Base: @kwdef
import Base: show, getproperty

include("convenience_functions.jl")
# aliases
Maybe(T)=Union{T,Nothing}
const F64orFunc = Union{Float64,Function}

### PRIMITIVE VARIABLES ###
@kwdef struct PrimitiveVariables{Nvars}
    names::Union{Symbol,NTuple{Nvars,Symbol}}
    interpolations::Union{Interpolation,NTuple{Nvars,Interpolation}}
    function PrimitiveVariables{Nvars}(names,ips) where {Nvars}
        (typeof(names) == Symbol) && (@assert Nvars == 1)
        (typeof(names) == Tuple) && (@assert Nvars == length(names) == length(ips))
        return new(names,ips)
    end
end

refshape(::Type{Cell{1,N,2}}) where{N} = RefCube
refshape(::Type{Cell{2,N,4}}) where{N} = RefCube
refshape(::Type{Cell{3,N,6}}) where{N} = RefCube
refshape(::Type{Cell{2,N,3}}) where{N} = RefTetrahedron
refshape(::Type{Cell{3,N,4}}) where{N} = RefTetrahedron

function PrimitiveVariables{NV}(names,ip_order,el_geom::Type{Cell{dim,N,M}}) where {NV,dim,N,M}
    NV == 1 && return PrimitiveVariables{NV}(names, (Lagrange{dim,refshape(el_geom),ip_order[1]}(),))
    NV == 2 && return PrimitiveVariables{NV}(names, (Lagrange{dim,refshape(el_geom),ip_order[1]}(),
                                          Lagrange{dim,refshape(el_geom),ip_order[2]}()) )
    @error "More than two primitive variables (:u and :p or :D) not implemented yet"
end

# function Variables(names::Symbol,ips::Interpolation)
#     Variables{1}(names,ips)
# end
# function Variables(names::NTuple{N,Symbol},ips::NTuple{N,Interpolation}) where {N}
#     Variables{N}(names,ips)
# end

### BOUNDARY CONDITIONS ###

@kwdef struct BoundaryConditions
    dirichlet::Dict
    neumann::Maybe(Dict)
end
BoundaryConditions(::Nothing,::Dict) = @error "A minimal set of Dirichlet boundary conditions must be provided in order to prevent rigid motion of the material"

### BODY FORCES ###
struct BodyForces{T}
    components::T
    function BodyForces(c::T) where {T}
        if T <: AbstractVector
            dim = length(c)
            if T <: Vector
                c = Vec{dim}(c)
                new{typeof(c)}(c)
            elseif T <: Tensor
                return new{T}(c)
            end
        elseif T <: Function
            return new{T}(c)
        else
            @error "wrong field type. Must be an AbstractVector"
        end
    end
end

# other types :

include("material_behaviors_types.jl")
include("rheology_type.jl")
include("clock_type.jl")
include("solvers_types.jl")
include("material_state_types.jl")
include("output_writer_types.jl")
include("model_type.jl")

## TODO : Check whether a material parameters type is needed
### MATERIAL PARAMETERS ###

# struct MaterialProperties{T}
#     mp :: Vector{Rheology}
#     function MaterialProperties(mp)
#         T = eltype(mp)
#         new{T}(mp)
#     end
# end
# MaterialProperties{T}(mp) where T = MaterialProperties(mp)
# function Base.getproperty(MP::MaterialProperties,sym::Symbol)
#     if sym == :mp
#         Core.getfield(MP::MaterialProperties,sym::Symbol)
#     else
#         Core.getproperty(Core.getfield(MP,:mp),sym)
#     end
# end
#
# Base.setproperty!(MP::MaterialProperties,prop::Symbol,val) = Core.setproperty!(Core.getfield(MP,:mp),prop,val)
#
# function Base.propertynames(MP::MaterialProperties)
#     properties = Set{Symbol}()
#     for fields in fieldMP.mp
#         push!(properties, propertynames())
#     end
#
# end