text:
@defcomp GDP begin
# GDP: Gross domestic product $M
# GRW: GDP growth rate %/year
    region            = Index()

    # Variables
    gdp               = Variable(index=[time, region], unit="\$M")
    cons_consumption  = Variable(index=[time, region], unit="\$million")
    cons_percap_consumption = Variable(index=[time, region], unit="\$/person")
    cons_percap_consumption_0 = Variable(index=[region], unit="\$/person")
    yagg_periodspan = Variable(index=[time], unit="year")

    # Parameters
    y_year_0          = Parameter(unit="year")
    y_year            = Parameter(index=[time], unit="year")
    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year") #From p.32 of Hope 2009
    gdp_0             = Parameter(index=[region], unit="\$M") #GDP in y_year_0
    save_savingsrate  = Parameter(unit="%", default=15.00) #pp33 PAGE09 documentation, "savings rate".
    pop0_initpopulation = Parameter(index=[region], unit="million person")
    pop_population    = Parameter(index=[time,region],unit="million person")

    # Saturation, used in impacts
    isat0_initialimpactfxnsaturation = Parameter(unit="unitless", default=33.333333333333336) #pp34 PAGE09 documentation
    isatg_impactfxnsaturation = Variable(unit="unitless")

    function init(p, v, d)

        v.isatg_impactfxnsaturation = p.isat0_initialimpactfxnsaturation * (1 - p.save_savingsrate/100)
        for rr in d.region
            v.cons_percap_consumption_0[rr] = (p.gdp_0[rr] / p.pop0_initpopulation[rr])*(1 - p.save_savingsrate / 100)
        end
    end

    function run_timestep(p, v, d, t)

        # Analysis period ranges - required for abatemnt costs and equity weighting, from Hope (2006)
        if is_first(t)
            ylo_periodstart = p.y_year_0
        else
            ylo_periodstart = (p.y_year[t] + p.y_year[t-1]) / 2
        end

        if t.t == length(p.y_year)
            yhi_periodend = p.y_year[t]
        else
            yhi_periodend = (p.y_year[t] + p.y_year[t+1]) / 2
        end

        v.yagg_periodspan[t] = yhi_periodend- ylo_periodstart

        for r in d.region
            #eq.28 in Hope 2002
            if is_first(t)
                v.gdp[t, r] = p.gdp_0[r] * (1 + (p.grw_gdpgrowthrate[t,r]/100))^(p.y_year[t] - p.y_year_0)
            else
                v.gdp[t, r] = v.gdp[t-1, r] * (1 + (p.grw_gdpgrowthrate[t,r]/100))^(p.y_year[t] - p.y_year[t-1])
            end
            v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
            v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]
        end
    end
end

ext: jl


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: managedblockchain
using AWS.Compat
using AWS.UUIDs
"""
    CreateMember()

Creates a member within a Managed Blockchain network.

# Required Parameters
- `ClientRequestToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the operation. An idempotent operation completes no more than one time. This identifier is required only if you make a service request directly using an HTTP client. It is generated automatically if you use an AWS SDK or the AWS CLI.
- `InvitationId`: The unique identifier of the invitation that is sent to the member to join the network.
- `MemberConfiguration`: Member configuration parameters.
- `networkId`: The unique identifier of the network in which the member is created.

"""

create_member(ClientRequestToken, InvitationId, MemberConfiguration, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/members", Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "InvitationId"=>InvitationId, "MemberConfiguration"=>MemberConfiguration); aws_config=aws_config)
create_member(ClientRequestToken, InvitationId, MemberConfiguration, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/members", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "InvitationId"=>InvitationId, "MemberConfiguration"=>MemberConfiguration), args)); aws_config=aws_config)

"""
    CreateNetwork()

Creates a new blockchain network using Amazon Managed Blockchain.

# Required Parameters
- `ClientRequestToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the operation. An idempotent operation completes no more than one time. This identifier is required only if you make a service request directly using an HTTP client. It is generated automatically if you use an AWS SDK or the AWS CLI.
- `Framework`: The blockchain framework that the network uses.
- `FrameworkVersion`: The version of the blockchain framework that the network uses.
- `MemberConfiguration`: Configuration properties for the first member within the network.
- `Name`: The name of the network.
- `VotingPolicy`:  The voting rules used by the network to determine if a proposal is approved.

# Optional Parameters
- `Description`: An optional description for the network.
- `FrameworkConfiguration`:  Configuration properties of the blockchain framework relevant to the network configuration.
"""

create_network(ClientRequestToken, Framework, FrameworkVersion, MemberConfiguration, Name, VotingPolicy; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks", Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "Framework"=>Framework, "FrameworkVersion"=>FrameworkVersion, "MemberConfiguration"=>MemberConfiguration, "Name"=>Name, "VotingPolicy"=>VotingPolicy); aws_config=aws_config)
create_network(ClientRequestToken, Framework, FrameworkVersion, MemberConfiguration, Name, VotingPolicy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "Framework"=>Framework, "FrameworkVersion"=>FrameworkVersion, "MemberConfiguration"=>MemberConfiguration, "Name"=>Name, "VotingPolicy"=>VotingPolicy), args)); aws_config=aws_config)

"""
    CreateNode()

Creates a peer node in a member.

# Required Parameters
- `ClientRequestToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the operation. An idempotent operation completes no more than one time. This identifier is required only if you make a service request directly using an HTTP client. It is generated automatically if you use an AWS SDK or the AWS CLI.
- `NodeConfiguration`: The properties of a node configuration.
- `memberId`: The unique identifier of the member that owns this node.
- `networkId`: The unique identifier of the network in which this node runs.

"""

create_node(ClientRequestToken, NodeConfiguration, memberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/members/$(memberId)/nodes", Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "NodeConfiguration"=>NodeConfiguration); aws_config=aws_config)
create_node(ClientRequestToken, NodeConfiguration, memberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/members/$(memberId)/nodes", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientRequestToken"=>ClientRequestToken, "NodeConfiguration"=>NodeConfiguration), args)); aws_config=aws_config)

"""
    CreateProposal()

Creates a proposal for a change to the network that other members of the network can vote on, for example, a proposal to add a new member to the network. Any member can create a proposal.

# Required Parameters
- `Actions`: The type of actions proposed, such as inviting a member or removing a member. The types of Actions in a proposal are mutually exclusive. For example, a proposal with Invitations actions cannot also contain Removals actions.
- `ClientRequestToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the operation. An idempotent operation completes no more than one time. This identifier is required only if you make a service request directly using an HTTP client. It is generated automatically if you use an AWS SDK or the AWS CLI.
- `MemberId`: The unique identifier of the member that is creating the proposal. This identifier is especially useful for identifying the member making the proposal when multiple members exist in a single AWS account.
- `networkId`:  The unique identifier of the network for which the proposal is made.

# Optional Parameters
- `Description`: A description for the proposal that is visible to voting members, for example, \"Proposal to add Example Corp. as member.\"
"""

create_proposal(Actions, ClientRequestToken, MemberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/proposals", Dict{String, Any}("Actions"=>Actions, "ClientRequestToken"=>ClientRequestToken, "MemberId"=>MemberId); aws_config=aws_config)
create_proposal(Actions, ClientRequestToken, MemberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/proposals", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Actions"=>Actions, "ClientRequestToken"=>ClientRequestToken, "MemberId"=>MemberId), args)); aws_config=aws_config)

"""
    DeleteMember()

Deletes a member. Deleting a member removes the member and all associated resources from the network. DeleteMember can only be called for a specified MemberId if the principal performing the action is associated with the AWS account that owns the member. In all other cases, the DeleteMember action is carried out as the result of an approved proposal to remove a member. If MemberId is the last member in a network specified by the last AWS account, the network is deleted also.

# Required Parameters
- `memberId`: The unique identifier of the member to remove.
- `networkId`: The unique identifier of the network from which the member is removed.

"""

delete_member(memberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/networks/$(networkId)/members/$(memberId)"; aws_config=aws_config)
delete_member(memberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/networks/$(networkId)/members/$(memberId)", args; aws_config=aws_config)

"""
    DeleteNode()

Deletes a peer node from a member that your AWS account owns. All data on the node is lost and cannot be recovered.

# Required Parameters
- `memberId`: The unique identifier of the member that owns this node.
- `networkId`: The unique identifier of the network that the node belongs to.
- `nodeId`: The unique identifier of the node.

"""

delete_node(memberId, networkId, nodeId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)"; aws_config=aws_config)
delete_node(memberId, networkId, nodeId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)", args; aws_config=aws_config)

"""
    GetMember()

Returns detailed information about a member.

# Required Parameters
- `memberId`: The unique identifier of the member.
- `networkId`: The unique identifier of the network to which the member belongs.

"""

get_member(memberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)"; aws_config=aws_config)
get_member(memberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)", args; aws_config=aws_config)

"""
    GetNetwork()

Returns detailed information about a network.

# Required Parameters
- `networkId`: The unique identifier of the network to get information about.

"""

get_network(networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)"; aws_config=aws_config)
get_network(networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)", args; aws_config=aws_config)

"""
    GetNode()

Returns detailed information about a peer node.

# Required Parameters
- `memberId`: The unique identifier of the member that owns the node.
- `networkId`: The unique identifier of the network to which the node belongs.
- `nodeId`: The unique identifier of the node.

"""

get_node(memberId, networkId, nodeId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)"; aws_config=aws_config)
get_node(memberId, networkId, nodeId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)", args; aws_config=aws_config)

"""
    GetProposal()

Returns detailed information about a proposal.

# Required Parameters
- `networkId`: The unique identifier of the network for which the proposal is made.
- `proposalId`: The unique identifier of the proposal.

"""

get_proposal(networkId, proposalId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals/$(proposalId)"; aws_config=aws_config)
get_proposal(networkId, proposalId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals/$(proposalId)", args; aws_config=aws_config)

"""
    ListInvitations()

Returns a listing of all invitations for the current AWS account.

# Optional Parameters
- `maxResults`: The maximum number of invitations to return.
- `nextToken`: The pagination token that indicates the next set of results to retrieve.
"""

list_invitations(; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/invitations"; aws_config=aws_config)
list_invitations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/invitations", args; aws_config=aws_config)

"""
    ListMembers()

Returns a listing of the members in a network and properties of their configurations.

# Required Parameters
- `networkId`: The unique identifier of the network for which to list members.

# Optional Parameters
- `isOwned`: An optional Boolean value. If provided, the request is limited either to members that the current AWS account owns (true) or that other AWS accounts own (false). If omitted, all members are listed.
- `maxResults`: The maximum number of members to return in the request.
- `name`: The optional name of the member to list.
- `nextToken`: The pagination token that indicates the next set of results to retrieve.
- `status`: An optional status specifier. If provided, only members currently in this status are listed.
"""

list_members(networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members"; aws_config=aws_config)
list_members(networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members", args; aws_config=aws_config)

"""
    ListNetworks()

Returns information about the networks in which the current AWS account has members.

# Optional Parameters
- `framework`: An optional framework specifier. If provided, only networks of this framework type are listed.
- `maxResults`: The maximum number of networks to list.
- `name`: The name of the network.
- `nextToken`: The pagination token that indicates the next set of results to retrieve.
- `status`: An optional status specifier. If provided, only networks currently in this status are listed.
"""

list_networks(; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks"; aws_config=aws_config)
list_networks(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks", args; aws_config=aws_config)

"""
    ListNodes()

Returns information about the nodes within a network.

# Required Parameters
- `memberId`: The unique identifier of the member who owns the nodes to list.
- `networkId`: The unique identifier of the network for which to list nodes.

# Optional Parameters
- `maxResults`: The maximum number of nodes to list.
- `nextToken`: The pagination token that indicates the next set of results to retrieve.
- `status`: An optional status specifier. If provided, only nodes currently in this status are listed.
"""

list_nodes(memberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)/nodes"; aws_config=aws_config)
list_nodes(memberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/members/$(memberId)/nodes", args; aws_config=aws_config)

"""
    ListProposalVotes()

Returns the listing of votes for a specified proposal, including the value of each vote and the unique identifier of the member that cast the vote.

# Required Parameters
- `networkId`:  The unique identifier of the network.
- `proposalId`:  The unique identifier of the proposal.

# Optional Parameters
- `maxResults`:  The maximum number of votes to return.
- `nextToken`:  The pagination token that indicates the next set of results to retrieve.
"""

list_proposal_votes(networkId, proposalId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals/$(proposalId)/votes"; aws_config=aws_config)
list_proposal_votes(networkId, proposalId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals/$(proposalId)/votes", args; aws_config=aws_config)

"""
    ListProposals()

Returns a listing of proposals for the network.

# Required Parameters
- `networkId`:  The unique identifier of the network.

# Optional Parameters
- `maxResults`:  The maximum number of proposals to return.
- `nextToken`:  The pagination token that indicates the next set of results to retrieve.
"""

list_proposals(networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals"; aws_config=aws_config)
list_proposals(networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("GET", "/networks/$(networkId)/proposals", args; aws_config=aws_config)

"""
    RejectInvitation()

Rejects an invitation to join a network. This action can be called by a principal in an AWS account that has received an invitation to create a member and join a network.

# Required Parameters
- `invitationId`: The unique identifier of the invitation to reject.

"""

reject_invitation(invitationId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/invitations/$(invitationId)"; aws_config=aws_config)
reject_invitation(invitationId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("DELETE", "/invitations/$(invitationId)", args; aws_config=aws_config)

"""
    UpdateMember()

Updates a member configuration with new parameters.

# Required Parameters
- `memberId`: The unique ID of the member.
- `networkId`: The unique ID of the Managed Blockchain network to which the member belongs.

# Optional Parameters
- `LogPublishingConfiguration`: Configuration properties for publishing to Amazon CloudWatch Logs.
"""

update_member(memberId, networkId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("PATCH", "/networks/$(networkId)/members/$(memberId)"; aws_config=aws_config)
update_member(memberId, networkId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("PATCH", "/networks/$(networkId)/members/$(memberId)", args; aws_config=aws_config)

"""
    UpdateNode()

Updates a node configuration with new parameters.

# Required Parameters
- `memberId`: The unique ID of the member that owns the node.
- `networkId`: The unique ID of the Managed Blockchain network to which the node belongs.
- `nodeId`: The unique ID of the node.

# Optional Parameters
- `LogPublishingConfiguration`: Configuration properties for publishing to Amazon CloudWatch Logs.
"""

update_node(memberId, networkId, nodeId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("PATCH", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)"; aws_config=aws_config)
update_node(memberId, networkId, nodeId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("PATCH", "/networks/$(networkId)/members/$(memberId)/nodes/$(nodeId)", args; aws_config=aws_config)

"""
    VoteOnProposal()

Casts a vote for a specified ProposalId on behalf of a member. The member to vote as, specified by VoterMemberId, must be in the same AWS account as the principal that calls the action.

# Required Parameters
- `Vote`:  The value of the vote.
- `VoterMemberId`: The unique identifier of the member casting the vote.
- `networkId`:  The unique identifier of the network.
- `proposalId`:  The unique identifier of the proposal.

"""

vote_on_proposal(Vote, VoterMemberId, networkId, proposalId; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/proposals/$(proposalId)/votes", Dict{String, Any}("Vote"=>Vote, "VoterMemberId"=>VoterMemberId); aws_config=aws_config)
vote_on_proposal(Vote, VoterMemberId, networkId, proposalId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = managedblockchain("POST", "/networks/$(networkId)/proposals/$(proposalId)/votes", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Vote"=>Vote, "VoterMemberId"=>VoterMemberId), args)); aws_config=aws_config)


text:
struct Periodic
    rs
end

function calcLeftFlux(bc::Periodic,
                      state::HydroSnapshot,
                      edge_velocities::Array{Float64,1},
                      eos)::Conserved
    return calcFlux(bc.rs,
                    state.cells[end],
                    state.cells[1],
                    edge_velocities[1],
                    eos)
end

function calcRightFlux(bc::Periodic,
                       state::HydroSnapshot,
                       edge_velocities::Array{Float64,1},
                       eos)::Conserved
    return calcFlux(bc.rs,
                    state.cells[end],
                    state.cells[1],
                    edge_velocities[end],
                    eos)
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
module CrypticServer

using Mux: Mux, @app, page, serve
using CrypticCrosswords: CrypticCrosswords, solve, derive!, DerivedArc, DerivedSolution, explain
using HTTP: queryparams
using Sockets: @ip_str
using JSON: json, JSON

export main

include("Templates.jl")
using .Templates

function split_query(app, req)
    req[:query_params] = queryparams(req[:query])
    app(req)
end

const STATIC_DIR = normpath(joinpath(@__DIR__, "..", "frontend", "dist"))

function maybe(f, dict, key, default=nothing)
    if key in keys(dict)
        f(dict[key])
    else
        default
    end
end

decode_clue(s::AbstractString) = replace(s, '+' => ' ')

function handle_solve(request, timeout=15.0)
    @show request
    params = request[:query_params]
    @show params
    clue = decode_clue(params["clue"])
    length = tryparse(Int, get(params, "length", ""))
    pattern = get(params, "pattern", "")
    @show clue length pattern
    start_time = time()
    timed_out = Ref(false)
    should_continue = let timed_out=timed_out, start_time=start_time, timeout=timeout
        function ()
            elapsed = time() - start_time
            if elapsed > timeout
                timed_out[] = true
                return false
            else
                return true
            end
        end
    end
    solutions, state = solve(clue,
                             length=length,
                             pattern=Regex(pattern),
                             should_continue=should_continue
                             )
    derivations = Iterators.flatten([derive!(state, s) for s in Iterators.take(solutions, 10)])
    Templates.Index(
        Templates.HomeLink() *
        Templates.ClueInput(clue, length, pattern) *
        (timed_out[] ? "
Timed out
" : Templates.Results(derivations))
    )
end

function handle_home(request)
    Templates.Index(
        Templates.HomeLink() *
        Templates.Intro() *
        Templates.ClueInput() *
        Templates.Examples()
    )
end

function main(; host=ip"127.0.0.1", port=8000)
    @app server = (
        Mux.stack(Mux.todict, Mux.basiccatch, Mux.splitquery, Mux.toresponse),
        page("/solve",
             split_query,
             req -> Base.invokelatest(handle_solve, req)),
        page("/",
             req -> Base.invokelatest(handle_home, req)),
        Mux.notfound())
    serve(server, host, port)
end

end


text:
module TheAlgorithms

# Usings/Imports (keep sorted)
using DataFrames
using DifferentialEquations
using GLM
using LinearAlgebra
using Plots
using Random

## Exports
# Please keep the folders/functions sorted

# Exports: knapsack
export ZeroOnePack!
export CompletePack!

# Exports: math
export abs_max
export abs_min
export abs_val
export area_circle
export area_ellipse
export area_heron_triangle
export area_parallelogram
export area_rectangle
export area_rhombus
export area_square
export area_trapezium
export area_triangle
export ceil_val, floor_val
export collatz_sequence
export euler_method
export is_armstrong
export line_length
export mean
export median
export mode
export SIR # TODO: make the name lowercase if possible
export surfarea_cube
export surfarea_sphere
export trapazoidal_area

# Exports: matrix
export determinant
export lu_decompose

# Exports: project-rosalind
export count_nucleotides
export dna2rna
export reverse_complement

# Exports: searches
export binarysearch
export search

# Exports: statistics
export OLSbeta
export PearsonCorrelation
# Exports: strings
export is_palindrome

# Exports: scheduling
export fcfs

## Includes
# Please keep the folders/files sorted (by dependencies then alphabetical order)

# Includes: knapsack
include("knapsack/knapsack.jl")

# Includes: math
include("math/abs.jl")
include("math/area_under_curve.jl")
include("math/armstrong_number.jl")
include("math/area.jl")
include("math/average_mean.jl")
include("math/average_median.jl")
include("math/average_mode.jl")
include("math/ceil_floor.jl")
include("math/collatz_sequence.jl")
include("math/euler_method.jl")
include("math/line_length.jl")
include("math/sir_model.jl")

# Includes: matrix
include("matrix/lu_decompose.jl") # used by determinant.jl
include("matrix/determinant.jl")

# Includes: project-rosalind
include("project-rosalind/count_nucleotide.jl")
include("project-rosalind/dna2rna.jl")
include("project-rosalind/reverse_complement.jl")

# Includes: search
include("searches/binary_search.jl")
include("searches/linear_search.jl")

# Includes: statistics
include("statistics/ordinary_least_squares.jl")

# Includes: strings
include("strings/is_palindrome.jl")

# Includes: scheduling
include("scheduling/fcfs.jl")

end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
using PkgBenchmark
using Dates

if (length(ARGS) == 1) && (ARGS[1] == "-b")
    path =  joinpath(@__DIR__, "results", "$(Dates.format(now(), dateformat"yyyymmddTHHMMSS"))-baseline.json")
else
    path =  joinpath(@__DIR__, "results", "$(Dates.format(now(), dateformat"yyyymmddTHHMMSS")).json")
end

benchmarkpkg(
    dirname(@__DIR__),
    BenchmarkConfig(
        env = Dict(
            "JULIA_NUM_THREADS" => "1",
            "OMP_NUM_THREADS" => "1",
        ),
    ),
    resultfile = path,
)


text:
using GenericLinearAlgebra
export adjust_for_errtype!, solve_linlsqr!

"""
    adjust_for_errtype!(A, b, objfun_vals, errtype)

Adjusts the matrix `A` and vector `b` to `errtype`. `A` is a matrix, e.g., the
Jacobian or system matrix in a least squares problem. `b` is a vector, e.g., the
residuals or the right-hand side in a lieast squares problem. `objfun_vals` is a
vector of objective function values, and `errtype` can be `:abserr` or
`:relerr`, i.e., absolute error or relative error.
"""
function adjust_for_errtype!(A, b, objfun_vals, errtype)
    # Assumes a Jacobian and residual-vector that is computed in absolute error.
    # Adjusts Jacobian and residual-vector from absolute, to relative error,
    # i.e., objective function 1/2 sum_i (Z(x_i)-f(x_i))^2 / f(x_i)^2 and
    if (errtype == :abserr)
        # Do nothing
    elseif (errtype == :relerr)
        objfun_vals[objfun_vals.==0] .= eps() * 100 # Avoid division by zero
        D = Diagonal(1 ./ objfun_vals)
        A[:] = D * A
        b[:] = D * b
    else
        error("Unknown errtype '", errtype, "'.")
    end
    return (A, b)
end

"""
    d = solve_linlsqr!(A, b, linlsqr, droptol)

Solves the linear least squares problem

    Ad=b.

The argument `linlsqr` determines how the linear least squares problem is
solved. It can be `:backslash`, `:real_backslash`, `:nrmeq`, `:real_nrmeq`,
`:svd`, or `:real_svd`. For the latter two options singular values below
`droptol` are disregarded. The `:real_X` options optimizes `d` in the space of
real vectors. The input matrix `A` is sometimes overwritten.
"""
function solve_linlsqr!(A, b, linlsqr, droptol)
    if (linlsqr == :backslash)
        d = A \ b
    elseif (linlsqr == :real_backslash)
        d = vcat(real(A), imag(A)) \ vcat(real(b), imag(b))
    elseif (linlsqr == :nrmeq)
        d = (A' * A) \ (A' * b)
    elseif (linlsqr == :real_nrmeq)
        Ar = real(A)
        Ai = imag(A)
        br = real(b)
        bi = imag(b)
        d = (Ar' * Ar + Ai' * Ai) \ (Ar' * br + Ai' * bi)
    elseif (linlsqr == :svd) || (linlsqr == :real_svd)
        if (linlsqr == :real_svd)
            A = vcat(real(A), imag(A))
            b = vcat(real(b), imag(b))
        end
        if (eltype(A) == BigFloat || eltype(A) == Complex{BigFloat})
            Sfact = svd!(A; full = false, alg = nothing)
        else
            Sfact = svd(A)
        end
        d = Sfact.S
        # Use pseudoinverse if droptol>0
        II = (d / d[1]) .< droptol
        dinv = 1 ./ d
        dinv[II] .= 0
        # No explicit construction, only multiplication
        # JJ0=Sfact.U*Diagonal(d)*Sfact.Vt
        d = Sfact.V * (dinv .* (Sfact.U' * b))
    end
    return d
end


text:
"""
Check if vertex has property
"""
function has_prop_vertex(value,g,prop)
    gvs = [ g[v,prop] for v in vertices(g) ]
    return value in gvs
end

"""
Ordering tuples of paths
1: tick
2: score
"""
function isless_paths(a,b)
    # return @match a,b begin
    #     a[2] < b[2]
    # end
    if a[2] == b[2]
        return a[1] > b[1]
    else
        return a[2] < b[2]
    end
end

"""
Join two subgraphs assuming they are both part of a global graph.
The id in the global graph is given by property :eid.
"""
function join_subgraphs(g1,g2)
    gt = deepcopy(g1)
    eqv = []
    for v in vertices(g2)

        gv = g2[v,:eid]

        lv = to_local_vertex(gt,gv)

        if lv == 0
            add_vertex!(gt,:eid,gv)
            push!(eqv,(v,nv(gt)))
        else
            push!(eqv,(lv,gv))
        end
    end

    for e in edges(g2)

        src_t = to_local_vertex(gt,g2[src(e),:eid])
        dst_t = to_local_vertex(gt,g2[dst(e),:eid])
        add_edge!(gt,src_t,dst_t)

        # add_edge!(gt,
        #     first([last(x) for x in eqv if first(x) == src(e) ]),
        #     first([last(x) for x in eqv if first(x) == dst(e) ]),
        # )
    end
    return gt
end

#function (lg, predictions_traffic, predictions_rul)

"""
Search for a path between nodes s and d in the local graph lg
It assumes property :eid of each vertex is global id of vertex
"""

# function query_paths(lg::SimpleWeightedGraph,s,d)
function query_paths(lg,s,d)
    ls = to_local_vertex(lg,s)
    ld = to_local_vertex(lg,d)
    paths = []
    scores = []
    result =   YenState{Float64,Int64}(scores,paths)

    if ls > 0 && ld > 0
            #slg = SimpleGraph(lg)
            #return yen_k_shortest_paths(slg,ls,ld, weights(slg),2,Inf)
            result = yen_k_shortest_paths(lg,ls,ld)
    end

    #gvs = [ lg[v,:eid] for v in vertices(lg)]
    #log_info("network contains: gvs: $gvs")
    # log_info("query_path:  g v: $(vertices(lg)), s: $(s) - ls: $(ls), d: $d - ld $ld result ==> $(result)")

    for path in result.paths
        # path = !isempty(path) && typeof(path) == Array{Array{Int64,1},1} ? first(path) : path
        #convert paths to global graph (eids)
        cpath = [ lg[v,:eid] for v in path]
        push!(paths,cpath)
        push!(scores,score_path(cpath))
    end

    result =   YenState{Float64,Int}(scores,paths)

    return result
end

"""
It gives a score to the given path, initially only based on length of path
"""
function score_path(path)
    return length(path)
end

"""
Local search receiving source and destination in a tuple
"""
function query_paths(lg,t)
    query_paths(lg,t...)
end

"""
obtains local id of a vertex given its global id in property :eid
"""
function to_local_vertex(lg,gv)
    lv = [ x for x=1:nv(lg) if lg[x,:eid] == gv]
    return isempty(lv) ? 0 : first(lv)
end

"""
obtains local id of a vertex given its global id in property in gid
"""
function to_local_vertex(lg,gv,gid::Symbol)
    lv = [ x for x=1:nv(lg) if lg[x,gid] == gv]
    return isempty(lv) ? 0 : first(lv)
end

"""
Creates a subgraph (MetaGraph) for the given
adjacency matrix (m) and vector of equivalences (eqv).
In eqv, every pair has the form: (lv,gv) where lv is the
local vertex id and gv is the global vertex id.
"""
function create_subgraph(m,eqv)
    gw = SimpleWeightedGraph(m)
    g = MetaGraph(m)
    [ set_prop!(g, r, c, :weight, SimpleWeightedGraphs.weights(gw)[r,c]) for r=1:size(SimpleWeightedGraphs.weights(gw),1),c=1:size(SimpleWeightedGraphs.weights(gw),2) if SimpleWeightedGraph.weights(gw)[r,c] >0]

    for eq in eqv
        set_props!(g,first(eq),Dict(:eid=>last(eq)))
    end
    set_indexing_prop!(g, :eid)
    return g
end

"""
Creates a subgraph (MetaGraph) for the given
edge list and vector of equivalences (eqv).
In eqv, every pair has the form: (lv,gv) where lv is the
local vertex id and gv is the global vertex id.
"""
function create_subgraph(egs,eqv,gid_prop)
    #log_info("Creating subgraph egs: $(egs) and eqv: $eqv")
    g = MetaGraph()
    set_indexing_prop!(g, gid_prop)

    #create vertices
    n_v = max([ src(e) > dst(e) ? src(e) : dst(e) for e in egs]...)

    for v=1:n_v
        gid = last(first([ x for x in eqv if first(x) == v]))
        add_vertex!(g,gid_prop,gid)
    end

    #create edges

    for e in egs
        add_edge!(g,src(e),dst(e))
        add_edge!(g,dst(e),src(e))
    end

    return g
end

function soft_remove_vertex(g::AbstractGraph,dpn_id::Int)

    new_g = deepcopy(g)
    nbs₀ = deepcopy(all_neighbors(new_g,dpn_id))

    for nb in nbs₀
        rem_edge!(new_g,dpn_id,nb)
        rem_edge!(new_g,nb,dpn_id)
    end
    return new_g#
end

function remove_vertex(g::AbstractGraph,dpn_id::Int)
    sm_g = sparse(g)
    sm_new_g = spzeros((nv(g)-1),(nv(g)-1))
    for i=1:nv(g)
        for j=1:nv(g)
            #log_info(" $i,$j value: $(sparse(ntw)[i,j])")
                x,y =   i < dpn_id && j < dpn_id ? (i,j) :
                        i < dpn_id && j > dpn_id ? (i,j-1) :
                        i > dpn_id && j < dpn_id ? (i-1,j) :
                        i > dpn_id && j > dpn_id ? (i-1,j-1) : (0,0)

                if x > 0 && y > 0
                    sm_new_g[x,y] = sm_g[i,j]
                    sm_new_g[y,x] = sm_g[j,i]
                end
        end
    end
    #[i >=dpn_id ? labels[i] = i+1 : labels[i] = i  for i in keys(labels)]
    return MetaGraph(sm_new_g)
end

function remove_vertices(g::AbstractGraph,dpn_ids::Array{Int})
    new_g = g
    for dpn_id in dpn_ids
        new_g = remove_vertex!(new_g,dpn_id)
    end
    return new_g
end

"""
Add edges between the local vertex lv and the list of gids (global ids passed) using property gid
"""
function add_edges_gids(g,lv,nbs_gids,gid)
    nb_lvs = []
    for v in vertices(g)
        if g[v,gid] in nbs_gids
            push!(nb_lvs,v)
        end
    end

    for nb in nb_lvs
        add_edge!(g,lv,nb)
    end

    return g
end

"""
get_graph

Get underlying graph, passing the adj matrix and separator.
"""
function get_graph(seed,size,topo;k=0,B=0,adj_m_csv=nothing,sep=';')
    Random.seed!(seed)
    ntw = @match topo begin
        GraphModel(0)=> load_graph_from_csv(adj_m_csv;sep=sep)#custom_topo
        GraphModel(2) => MetaGraph( LightGraphs.watts_strogatz(size,2,0))
        GraphModel(3) => MetaGraph(LightGraphs.complete_graph(size))
        GraphModel(4) => MetaGraph( Int.(Laplacians.grid2(Int(sqrt(size)))))
        GraphModel(5) => MetaGraph( SimpleGraph(Int.(hcat(vcat(zeros(1),ones(size-1)),vcat(transpose(ones(size-1)),zeros(size-1,size-1))))) )
        GraphModel(6) => MetaGraph(barabasi_albert(size,k,seed=seed))
        GraphModel(7) => MetaGraph(watts_strogatz(size,k,B,seed=seed))
    end
end
"""
It generates subgraph for the vector of nodes passed. This is similar to egonet but it keeps custom indexes given in id_prop parameter.

- g: graph
- It receives vector of controlled assets (nodes)
- id_prop: :eid (simNE) or :aid (agent)
"""
function get_subgraph(g,nodes,id_prop)
    # calculate local subgraph for the underlying network

    nbs = []

    for i=1:length(nodes)
        # get vertex for node id (assumes only one)
        v = first(filter(v->get_prop(g,v,id_prop) == nodes[i],1:nv(g)))

        #subgraph
        push!(nbs,[ get_prop(g,j,id_prop) for j in neighbors(g,v)])
        #push!(nbs,neighbors(g,nodes[i]))
        push!(nbs,[nodes[i]])
    end
    nnbs = vcat(nbs...)
    sub_g = deepcopy(g)
    vs = [ get_prop(sub_g,v,id_prop) for v in collect(vertices(sub_g))]
    to_del = [v for v ∈ vs if v ∉ nnbs]

    for d in to_del
        for v in collect(vertices(sub_g))
            if !has_prop(sub_g,v,id_prop) || get_prop(sub_g,v,id_prop) == d
                rem_vertex!(sub_g,v)
            end
        end
    end
    return sub_g
end

"""
Load the graph of the network to control
"""
function load_network_graph(graph::MetaGraph)
    ntw = deepcopy(graph)
    set_indexing_prop!(ntw,:eid)
    return ntw
end

"""
Load the graph of the control system
"""
function load_control_graph(graph::MetaGraph)
    ntw = deepcopy(graph)
    #indexing can't be done here because aid has not been assigned
    #set_indexing_prop!(ntw,:aid)
    return ntw
end

"""
Find a set of paths in the network such as a given coverage of node is ensured. e.g. if 0.95 is the coverage, it means that 95% of vertices of the graph are included in the resulting set of paths.
"""

function find_paths_by_seed(seed,g::G,coverage::Float64)where G<: AbstractGraph
    Random.seed!(seed)

    cc = closeness_centrality(g)
    cci = sort([ (i,cc[i]) for i=1:length(cc) ],by=last,rev=true)

    pending = cci
    cp =[]

        while length(pending) >= (1- coverage) * nv(g)
            pending_i = [ first(p) for p in pending]

            #node with the most closeness_centrality
            s = first(first(pending))

            #max distance to any other pending node
            ds = gdistances(g,s)
            sds = sort([(i,ds[i]) for i=1:length(ds) if i in pending_i ],by=last,rev=true)
            d = first(first(sds))

            #shortest path between these two nodes
            sp = first(yen_k_shortest_paths(g,s,d).paths)
            if length(sp) > 2 #only paths with more than 2 nodes
                push!(cp, sp)
            end

            #remove nodes in the shortest path from pending list
            pending_i = collect(setdiff([ first(p) for p in pending ],Set(sp)))
            # println(pending_i)
            pending = [ p for p in pending  if first(p) in pending_i ]
        end
    return cp
end

"""
Given a set of paths, return the start and end vertices
"""
function get_end_points(seed,g::G,coverage::Float64)where G<:AbstractGraph
    return [ (first(p),last(p)) for p in find_paths_by_seed(seed,g,coverage) if first(p) != last(p) ]
end

export get_graph


text:
using Distributed

if nworkers() == 1 && 1 in workers()
    addprocs(1)
end

@everywhere include("update.jl")

@everywhere function randomwalk_generator(rc, prob=0.5)
    x = 0
    put!(rc, x)
    while true
        x = update(x, prob)
        put!(rc, x)
    end
end

function main(nsteps, prob)
    rc = RemoteChannel()
    remote_do(randomwalk_generator, 2, rc, 0.5)
    for i in 1:nsteps
        x = take!(rc)
        println(i-1, " ", x)
    end
end
if abspath(PROGRAM_FILE) == @__FILE__
    nsteps = length(ARGS)>0 ? parse(Int, ARGS[1]) : 100
    prob = length(ARGS)>1 ? parse(Float64, ARGS[1]) : 0.5
    main(nsteps, prob)
end

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
using ProperOrthogonalDecomposition
using Test
using Random
using Statistics
using Pkg
using DelimitedFiles

# Define test matrix X with dimensions n×m, where n is number of data poitns and
# m is the number of snapshots
# W is the weight matrix representing the cell volume

Random.seed!(1)
X = rand(1000,10)
X .+= 10
Random.seed!(1)
W = rand(1:0.2:10,1000)

@testset "Standard POD" begin

    PODbase, Σ = POD(X)
    PODbaseEig, Σeig = PODeigen(X)
    PODbaseSvd, Σsvd = PODsvd(X)

    meanPODbaseEig, meanΣeig = PODeigen(copy(X), subtractmean = true)
    meanPODbaseSvd, meanΣsvd = PODsvd(copy(X), subtractmean = true)

    Σ₁ = 1050.362168606664
    Σ₂ = 8.301301177004774
    meanΣ₁ = 9.878480670625322
    meanΣ₂ = 8.30589586913458

    @testset "POD using eigenvalue decomposition" begin
        @test Σeig[1] ≈ Σ₁
        @test Σeig[end] ≈ Σ₂
        @test meanΣeig[1] ≈ meanΣ₁
        @test meanΣeig[end-1] ≈ meanΣ₂

    end

    @testset "POD using SVD" begin
        @test Σsvd[1] ≈ Σ₁
        @test Σsvd[end] ≈ Σ₂
        @test meanΣsvd[1] ≈ meanΣ₁
        @test meanΣsvd[end-1] ≈ meanΣ₂
    end

    @testset "Default method" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbase.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbase.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σsvd.-Σ)) ≈ 0 atol=1e-8
    end

    @testset "Method equality of SVD and Eig" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbaseEig.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbaseEig.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(meanPODbaseSvd.coefficients) .- abs.(meanPODbaseEig.coefficients))) ≈ 0 atol=1e-7
        @test maximum(abs.(abs.(meanPODbaseSvd.modes[:,1:end-1]) .- abs.(meanPODbaseEig.modes[:,1:end-1]))) ≈ 0 atol=1e-7
        @test maximum(abs.(Σeig.-Σsvd)) ≈ 0 atol=1e-8
    end

    @testset "Rebuild solution" begin
        @test PODbaseEig.modes*PODbaseEig.coefficients ≈ X
        @test PODbaseSvd.modes*PODbaseSvd.coefficients ≈ X
        @test meanPODbaseEig.modes*meanPODbaseEig.coefficients ≈ X .- mean(X,dims=2)
        @test meanPODbaseSvd.modes*meanPODbaseSvd.coefficients ≈ X .- mean(X,dims=2)
    end
end

@testset "Weighted POD" begin

    PODbase, Σ = POD(X, W)
    PODbaseEig, Σeig = PODeigen(X,W)
    PODbaseSvd, Σsvd = PODsvd(X,W)

    meanPODbaseEig, meanΣeig = PODeigen(copy(X), W, subtractmean = true)
    meanPODbaseSvd, meanΣsvd = PODsvd(copy(X), W, subtractmean = true)

    Σ₁ = 2445.2807019537136
    Σ₂ = 19.41742998859931
    meanΣ₁ = 23.11141372908026
    meanΣ₂ = 19.444770050580598

    @testset "POD using eigenvalue decomposition" begin
        @test Σeig[1] ≈ Σ₁
        @test Σeig[end] ≈ Σ₂
        @test meanΣeig[1] ≈ meanΣ₁
        @test meanΣeig[end-1] ≈ meanΣ₂

    end

    @testset "POD using SVD" begin
        @test Σsvd[1] ≈ Σ₁
        @test Σsvd[end] ≈ Σ₂
        @test meanΣsvd[1] ≈ meanΣ₁
        @test meanΣsvd[end-1] ≈ meanΣ₂
    end

    @testset "Default method with weights" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbase.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbase.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σsvd.-Σ)) ≈ 0 atol=1e-8
    end

    @testset "Method equality of SVD and Eig with weights" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbaseEig.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbaseEig.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σeig.-Σsvd)) ≈ 0 atol=1e-8
    end

    @testset "Rebuild solution with weights" begin
        @test PODbaseEig.modes*PODbaseEig.coefficients ≈ X
        @test PODbaseSvd.modes*PODbaseSvd.coefficients ≈ X
    end
end

@testset "Mode convergence" begin
    A = [   0.010197212000627108 0.007070856617510225 0.004345509205119795 0
            0.6179627074745591 0.5476731024048103 0.12030902309667643 0
            1.377162205777348 0.8636384160114929 0.3197534098872419 0   ]
    Amean =[0.6322816398505274 0.5534565187905752 0.13120270384723295 0
            1.3821479511006756 0.8593667661723656 0.3083628970324381 0
            1.3536726312063883 1.1047589549024917 0.9957138077649103 0  ]
    W₂ = ones(1000)

    pkgpath = abspath(joinpath(dirname(Base.find_package("ProperOrthogonalDecomposition")), ".."))
    testdataPath = joinpath(pkgpath, "test","testdata.csv")

    @testset "Convergence of number of included snapshots" begin

        @test modeConvergence(X,PODeigen,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence(X,PODsvd,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),PODeigen!,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),PODsvd!,[1:4,1:6,1:8,1:10],3) ≈ A

    end

    @testset "Convergence of number of included snapshots with one weights" begin

        @test modeConvergence(X,x->PODeigen(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence(X,x->PODsvd(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),x->PODeigen!(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),x->PODsvd!(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A

    end

end


text:
export Dirichlet,
    Neumann, GeneralBoundaryConditions, SurfaceDrivenWaterBoundaryConditions

abstract type AbstractBoundaryFunctions end

"""
    struct Dirichlet{Fs, Fb} <: AbstractBoundaryFunctions

A concrete type to hold the surface state and bottom state variable
values/functions, if Dirichlet boundary conditions are desired.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Dirichlet{Fs, Fb} <: AbstractBoundaryFunctions
    "Surface state boundary condition"
    surface_state::Fs = nothing
    "Bottom state boundary condition"
    bottom_state::Fb = nothing
end

"""
    struct Neumann{Fs, Fb} <: AbstractBoundaryFunctions

A concrete type to hold the surface and/or bottom diffusive flux
values/functions, if Neumann boundary conditions are desired.

Note that these are intended to be scalar values. In the boundary_state!
functions, they are multiplied by the `ẑ` vector (i.e. the normal vector `n̂`
to the domain at the upper boundary, and -`n̂` at the lower boundary. These
normal vectors point out of the domain.)

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Neumann{Fs, Fb} <: AbstractBoundaryFunctions
    "Surface flux boundary condition"
    surface_flux::Fs = nothing
    "Bottom flux boundary condition"
    bottom_flux::Fb = nothing
end

"""
   AbstractBoundaryConditions
"""
abstract type AbstractBoundaryConditions end

"""
    GeneralBoundaryConditions{D, N} <: AbstractBoundaryConditions

A concrete example of the abstract type `AbstractBoundaryConditions`; to be used
when the user wishes to supply specific functions of space and time for
boundary conditions, at either the top or bottom of the domain, and either
of type Neumann or Dirichlet.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct GeneralBoundaryConditions{D, N} <: AbstractBoundaryConditions
    "Place to store Dirichlet Boundary conditions"
    dirichlet_bc::D
    "Place to store Neumann Boundary conditions"
    neumann_bc::N
end

"""
    GeneralBoundaryConditions(
        dirichlet_bc::AbstractBoundaryFunctions = Dirichlet(),
        neumann_bc::AbstractBoundaryFunctions = Neumann(),
    )

Constructor for the GeneralBoundaryConditions object. The default
is `nothing` boundary functions for the top and bottom of the domain,
for both Dirichlet and Neumann types. The user must change these in order
to apply boundary conditions.
"""
function GeneralBoundaryConditions(
    dirichlet_bc::AbstractBoundaryFunctions = Dirichlet(),
    neumann_bc::AbstractBoundaryFunctions = Neumann(),
)
    args = (dirichlet_bc, neumann_bc)
    return GeneralBoundaryConditions{typeof.(args)...}(args...)
end

"""
    SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <: AbstractBoundaryConditions

A concrete example of AbstractBoundaryConditions; this is to be used when the user wishes to
apply physical fluxes of water at the top of the domain (according to precipitation, runoff, and
evaporation rates) and zero flux at the bottom of the domain.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <:
       AbstractBoundaryConditions where {FT, PD, RD}
    "Precipitation model"
    precip_model::PD
    "Runoff model"
    runoff_model::RD
end

"""
    SurfaceDrivenWaterBoundaryConditions(
        ::Type{FT};
        precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(),
        runoff_model::AbstractSurfaceRunoffModel{FT} = NoRunoff(),
    ) where {FT}

Constructor for the SurfaceDrivenWaterBoundaryConditions object. The default
is a constant precipitation rate on the subgrid scale, and no runoff.
"""
function SurfaceDrivenWaterBoundaryConditions(
    ::Type{FT};
    precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(
        (t) -> (0.0),
    ),
    runoff_model::AbstractSurfaceRunoffModel = NoRunoff(),
) where {FT}
    args = (precip_model, runoff_model)
    return SurfaceDrivenWaterBoundaryConditions{FT, typeof.(args)...}(args...)
end

max_forks_repo_name: jm-c/CLIMA


text:
using WeavePynb
using Mustache

mmd(fname) = mmd_to_html(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")
## uncomment to generate just .md files
#mmd(fname) = mmd_to_md(fname, BRAND_HREF="../toc.html", BRAND_NAME="Calculus with Julia")

fnames = ["polar_coordinates",
          "vectors",
          "vector_valued_functions",
          "scalar_functions",
          "scalar_functions_applications",
          "vector_fields"
]

function process_file(nm, twice=false)
    include("$nm.jl")
    mmd_to_md("$nm.mmd")
    markdownToHTML("$nm.md")
    twice && markdownToHTML("$nm.md")
end

process_files(twice=false) = [process_file(nm, twice) for nm in fnames]

"""
## TODO differential_vector_calcululs

### Add questions for scalar_function_applications
* Newton's method??
* optimization. Find least squares for perpendicular distance using the same 3 points...??

"""


text:
# Loggers to record properties throughout a simulation

export
    TemperatureLogger,
    log_property!,
    CoordinateLogger,
    EnergyLogger,
    StructureWriter

"""
    TemperatureLogger(n_steps)
    TemperatureLogger(T, n_steps)

Log the temperature throughout a simulation.
"""
struct TemperatureLogger{T} <: Logger
    n_steps::Int
    temperatures::Vector{T}
end

TemperatureLogger(T::Type, n_steps::Integer) = TemperatureLogger(n_steps, T[])

TemperatureLogger(n_steps::Integer) = TemperatureLogger(DefaultFloat, n_steps)

function Base.show(io::IO, tl::TemperatureLogger)
    print(io, "TemperatureLogger{", eltype(tl.temperatures), "} with n_steps ",
                tl.n_steps, ", ", length(tl.temperatures),
                " temperatures recorded")
end

"""
    log_property!(logger, simulation, step_n)

Log a property thoughout a simulation.
Custom loggers should implement this function.
"""
function log_property!(logger::TemperatureLogger, s::Simulation, step_n::Integer)
    if step_n % logger.n_steps == 0
        push!(logger.temperatures, temperature(s))
    end
end

"""
    CoordinateLogger(n_steps; dims=3)

Log the coordinates throughout a simulation.
"""
struct CoordinateLogger{T} <: Logger
    n_steps::Int
    coords::Vector{Vector{T}}
end

function CoordinateLogger(T, n_steps::Integer; dims::Integer=3)
    return CoordinateLogger(n_steps,
                            Array{SArray{Tuple{dims}, T, 1, dims}, 1}[])
end

function CoordinateLogger(n_steps::Integer; dims::Integer=3)
    return CoordinateLogger(DefaultFloat, n_steps, dims=dims)
end

function Base.show(io::IO, cl::CoordinateLogger)
    print(io, "CoordinateLogger{", eltype(eltype(cl.coords)), "} with n_steps ",
                cl.n_steps, ", ", length(cl.coords), " frames recorded for ",
                length(cl.coords) > 0 ? length(first(cl.coords)) : "?", " atoms")
end

function log_property!(logger::CoordinateLogger, s::Simulation, step_n::Integer)
    if step_n % logger.n_steps == 0
        push!(logger.coords, deepcopy(s.coords))
    end
end

"""
    EnergyLogger(n_steps)

Log the energy of the system throughout a simulation.
"""
struct EnergyLogger{T} <: Logger
    n_steps::Int
    energies::Vector{T}
end

EnergyLogger(T::Type, n_steps::Integer) = EnergyLogger(n_steps, T[])

function EnergyLogger(n_steps::Integer)
    return EnergyLogger(DefaultFloat, n_steps)
end

function Base.show(io::IO, el::EnergyLogger)
    print(io, "EnergyLogger{", eltype(el.energies), "} with n_steps ",
                el.n_steps, ", ", length(el.energies), " energies recorded")
end

function log_property!(logger::EnergyLogger, s::Simulation, step_n::Integer)
    if step_n % logger.n_steps == 0
        push!(logger.energies, energy(s))
    end
end

"""
    StructureWriter(n_steps, filepath)

Write 3D output structures to the PDB file format throughout a simulation.
"""
mutable struct StructureWriter <: Logger
    n_steps::Int
    filepath::String
    structure_n::Int
end

StructureWriter(n_steps::Integer, filepath::AbstractString) = StructureWriter(
        n_steps, filepath, 1)

function Base.show(io::IO, sw::StructureWriter)
    print(io, "StructureWriter with n_steps ", sw.n_steps, ", filepath \"",
                sw.filepath, "\", ", sw.structure_n - 1, " frames written")
end

function log_property!(logger::StructureWriter, s::Simulation, step_n::Integer)
    if step_n % logger.n_steps == 0
        append_model(logger, s)
        logger.structure_n += 1
    end
end

function append_model(logger::StructureWriter, s::Simulation)
    open(logger.filepath, "a") do output
        println(output, "MODEL     ", lpad(logger.structure_n, 4))
        for (i, coord) in enumerate(s.coords)
            at_rec = atomrecord(s.atoms[i], i, coord)
            println(output, pdbline(at_rec))
        end
        println(output, "ENDMDL")
    end
end

atomrecord(at::Atom   , i, coord) = AtomRecord(false, i, at.name, ' ', at.resname, "A", at.resnum,
                                                ' ', 10 .* coord, 1.0, 0.0, "  ", "  ")

atomrecord(at::AtomMin, i, coord) = AtomRecord(false, i, "??"   , ' ', "???"     , "A", 0        ,
                                                ' ', 10 .* coord, 1.0, 0.0, "  ", "  ")


text:
if isdefined(Base, :datatype_name)
    datatype_name(x) = Base.datatype_name(x)
else
    datatype_name(t) = t.name.name
end

const Str = String
import Dates
import Dates: CONVERSION_SPECIFIERS, TimeType, DatePart, tryparsenext, character_codes, genvar, CONVERSION_TRANSLATIONS, CONVERSION_DEFAULTS, _directives, DateLocale
include("date-tryparse-internal.jl")
const ISODateFormat = Dates.ISODateFormat
const ISODateTimeFormat = Dates.ISODateTimeFormat
const RFC1123Format = Dates.RFC1123Format


text:
include("Mortgage.jl")
mutable struct houseMortgage <: Mortgage

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
## ------------------------------------------------------------------------------------------------------------------
const REGEXS = Dict{Symbol, Regex}()

function _fill_regexs!()

    ## --------------------------------------------------------
    _blank = "\\t\\s"
    REGEXS[:TEX_COMMENT_REGEX] = Regex("^(?<uncomment>[^%]*)(?<comment>%.*\$)")

    ## --------------------------------------------------------
    # ASS SECTION
    _head = "%%"
    _blank = "\\t\\s"
    _fill = "-"
    texass_cmd_reg_str(reg_str_to_insert) = string(
        "^[$(_blank)]*",
        "$(_head)[$(_blank)]*",
            "TexAssistant:", "[$(_blank)]+",
            reg_str_to_insert,
            "[$(_blank)]*",
        "[$(_blank)$(_fill)]*\$"
    )
    # Ex: "   %% TexAssistant: section{Section Name} ----------- "
    reg_str = texass_cmd_reg_str(string("section{", "(?<label>.*)", "}"))
    REGEXS[:TEXASS_SECTION_REGEX] = Regex(reg_str)
    # Ex: "   %% TexAssistant: end{label} ----------- "
    reg_str = texass_cmd_reg_str("end{(?<label>.*)}")
    REGEXS[:TEXASS_END_REGEX] = Regex(reg_str)

    ## --------------------------------------------------------
    # INPUT_FILE_REGEX
    # Ex: \input{sub.tex}
    _blank = "\\t\\s"
    _back = "\\\\"
    REGEXS[:TEX_INPUT_CMD_REGEX] = Regex("^[$(_blank)]*$(_back)input{(?<label>.*?)}$(_back)*\$")

    ## --------------------------------------------------------
    # TEX SECTIONS
    # Ex: \section{Cool section}
    # Ex: \subsection{Cool section}
    # Ex: \subsubsection{Cool section}
    _blank = "\\t\\s"
    _back = "\\\\"
    REGEXS[:TEX_SECTION_CMD_REGEX] = Regex("^[$(_blank)]*$(_back)section{(?<label>.*?)}.*\$")
    REGEXS[:TEX_SUBSECTION_CMD_REGEX] = Regex("^[$(_blank)]*$(_back)subsection{(?<label>.*?)}.*\$")
    REGEXS[:TEX_SUBSUBSECTION_CMD_REGEX] = Regex("^[$(_blank)]*$(_back)subsubsection{(?<label>.*?)}.*\$")

    ## --------------------------------------------------------
    # TEX BEGIN
    # Ex: \begin{document}
    # Ex: \end{document}
    _blank = "\\t\\s"
    _brakets = "(?:\\[(?<arg>.*)\\])"
    _back = "\\\\"
    REGEXS[:TEX_END_REGEX] = Regex("^[$(_blank)]*$(_back)end{(?<label>.*)}$(_back)*\$")
    REGEXS[:TEX_BEGIN_REGEX] = Regex("^[$(_blank)]*$(_back)begin{(?<label>.*)}$(_back)*$(_brakets)?$(_back)*\$")

    ## --------------------------------------------------------
    # INDENT
    _blank = "\\t\\s"
    REGEXS[:INDENT] = Regex("^(?<indent>[$(_blank)]*)?.*")

end


text:
# JSONWrapper for Swagger models handles
# - null fields
# - field names that are Julia keywords
struct JSONWrapper{T<:SwaggerModel} <: AbstractDict{Symbol, Any}
    wrapped::T
    flds::Vector{Symbol}
end

JSONWrapper(o::T) where {T<:SwaggerModel} = JSONWrapper(o, filter(n->hasproperty(o,n), propertynames(T)))

getindex(w::JSONWrapper, s::Symbol) = getproperty(w.wrapped, s)
keys(w::JSONWrapper) = w.flds
length(w::JSONWrapper) = length(w.flds)

function iterate(w::JSONWrapper, state...)
    result = iterate(w.flds, state...)
    if result === nothing
        return result
    else
        name,nextstate = result
        val = getproperty(w.wrapped, name)
        return (name=>val, nextstate)
    end
end

lower(o::T) where {T<:SwaggerModel} = JSONWrapper(o)

to_json(o) = JSON.json(o)

from_json(::Type{Union{Nothing,T}}, json::Dict{String,Any}) where {T} = from_json(T, json)
from_json(::Type{T}, json::Dict{String,Any}) where {T} = from_json(T(), json)
from_json(::Type{T}, json::Dict{String,Any}) where {T <: Dict} = convert(T, json)
from_json(::Type{T}, j::Dict{String,Any}) where {T <: String} = to_json(j)
from_json(::Type{Any}, j::Dict{String,Any}) = j

function from_json(o::T, json::Dict{String,Any}) where {T <: SwaggerModel}
    jsonkeys = [Symbol(k) for k in keys(json)]
    for name in intersect(propertynames(T), jsonkeys)
        from_json(o, name, json[String(name)])
    end
    o
end

function from_json(o::T, name::Symbol, json::Dict{String,Any}) where {T <: SwaggerModel}
    ftype = property_type(T, name)
    fval = from_json(ftype, json)
    setfield!(o, field_name(T,name), convert(ftype, fval))
    o
end
from_json(o::T, name::Symbol, v) where {T} = (setfield!(o, field_name(T,name), convert(property_type(T, name), v)); o)
function from_json(o::T, name::Symbol, v::Vector) where {T}
    # in Julia we can not support JSON null unless the element type is explicitly set to support it
    ftype = property_type(T, name)
    vtype = isa(ftype, Union) ? ((ftype.a === Nothing) ? ftype.b : ftype.a) : isa(ftype, Vector) ? ftype : Union{}
    (Nothing <: eltype(vtype)) || filter!(x->x!==nothing, v)
    setfield!(o, field_name(T,name), convert(property_type(T, name), v))
    o
end
from_json(o::T, name::Symbol, v::Nothing) where {T} = o


text:
# a bunch of interpolation helpers
function interpolate_state(mdp::Union{CarMDP, PedMDP}, state::VehicleState)
    # interpolate longitudinal position and velocity
    if state.posG == mdp.off_grid
        return VehicleState[state], Float64[1.0]
    end
    lane = get_lane(mdp.env.roadway, state)
    grid = mdp._car_grid[lane.tag]
    real_state = SVector{2, Float64}(state.posF.s, state.v)
    idx, weights = interpolants(grid, real_state)
    n_pts = length(idx)
    states = Vector{VehicleState}(n_pts)
    probs = zeros(n_pts)
    for i=1:n_pts
        sg, vg = ind2x(grid, idx[i])
        states[i] = VehicleState(Frenet(lane, sg), mdp.env.roadway, vg)
        probs[i] = weights[i]
    end
    return states, probs
end

function conservative_interpolation(mdp::Union{CarMDP, PedMDP}, state::VehicleState)
    # interpolate position and velocity separately
    if state.posG == mdp.off_grid
        return VehicleState[state], Float64[1.0]
    end
    lane = get_lane(mdp.env.roadway, state)
    route = get_ego_route(mdp.env)
    l_space = get_discretized_lane(lane.tag, mdp.env.roadway, mdp.pos_res)
    pgrid = mdp._l_grid[lane.tag]
    pidx, pweights = interpolants(pgrid, [state.posF.s])
    junction = state.posF.s > l_space[end] && lane.tag != route[end] ? 1 : 0
    jweight = 0.
    next_lane = nothing
    if junction == 1
        if length(pweights) != 1
            println(pweights)
            println([ind2x(pgrid, p) for p in pidx])
        end
        @assert length(pweights) == 1
        @assert pweights[1] == 1.
        dist_to_end = abs(get_end(lane) - state.posF.s)
        dist_to_prev = abs(state.posF.s - ind2x(pgrid, pidx[1])[1])
        jweight = dist_to_end/(dist_to_end + dist_to_prev)
        pweights[1] =  dist_to_prev/(dist_to_end + dist_to_prev)
        for (i, l) in enumerate(route)
            if l == lane.tag
                next_lane = mdp.env.roadway[route[i+1]]
                break
            end
        end
    end
    # println([ind2x(pgrid, pi) for pi in pidx])
    vgrid = mdp._v_grid
    vidx, vweights = interpolants(vgrid, [state.v])
    n_pts = length(pidx)*length(vidx)+junction*length(vidx)
    states = Vector{VehicleState}(n_pts)
    probs = zeros(n_pts)
    i = 1
    for j = 1:length(pidx)
        for k = 1:length(vidx)
            vg = ind2x(vgrid, vidx[k])[1]
            sg = ind2x(pgrid, pidx[j])[1]
            states[i] = VehicleState(Frenet(lane, sg), mdp.env.roadway, vg)
            probs[i] = pweights[j]*vweights[k]
            i += 1
        end
    end
    if junction == 1
        for l = 1:length(vidx)
            vg = ind2x(vgrid, vidx[l])[1]
            sg = 0.
            states[i] = VehicleState(Frenet(next_lane, sg), mdp.env.roadway, vg)
            probs[i] = jweight*vweights[l]
            i += 1
        end
    end
    normalize!(probs, 1)
    return states, probs
end

function handle_junction(mdp, state)
    lane = get_lane(mdp.env.roadway, state)

end

# take into account heading as well
function interpolate_pedestrian(mdp::Union{PedMDP}, state::VehicleState)
    # interpolate longitudinal position and velocity
    if state.posG == mdp.off_grid
        return (state,), (1.0,)
    end
    lane = get_lane(mdp.env.ped_roadway, state)
    grid = mdp._ped_grid[lane.tag]
    real_state = SVector{3, Float64}(state.posF.s, state.v, state.posF.ϕ)
    idx, weights = interpolants(grid, real_state)
    n_pts = length(idx)
    states = Vector{VehicleState}(n_pts)
    probs = zeros(n_pts)
    for i=1:n_pts
        sg, vg, phig = ind2x(grid, idx[i])
        states[i] = VehicleState(Frenet(lane, sg, 0., phig), mdp.env.ped_roadway, vg)
        probs[i] = weights[i]
    end
    return states, probs
end

function init_ped_grid(env::UrbanEnv, pos_res::Float64, vel_res::Float64)
    d = Dict{LaneTag, RectangleGrid{3}}()
    phi_space = SVector{2, Float64}(0., float(pi))
    v_space = get_ped_vspace(env, vel_res)
    for lane in env.crosswalks
        l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
        grid = RectangleGrid(l_space, v_space,  phi_space)
        d[lane.tag] = grid
    end
    return d
end

function init_car_grid(env::UrbanEnv, pos_res::Float64, vel_res::Float64)
    d = Dict{LaneTag, RectangleGrid{2}}()
    v_space = get_car_vspace(env, vel_res)
    for i=1:length(env.roadway.segments)
        for lane in env.roadway.segments[i].lanes
            l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
            grid = RectangleGrid(l_space, v_space)
            d[lane.tag] = grid
        end
    end
    return d
end
function init_l_grid(env::UrbanEnv, pos_res::Float64)
    d = Dict{LaneTag, RectangleGrid{1}}()
    for i=1:length(env.roadway.segments)
        for lane in env.roadway.segments[i].lanes
            l_space = get_discretized_lane(lane.tag, env.roadway, pos_res)
            grid = RectangleGrid(l_space)
            d[lane.tag] = grid
        end
    end
    return d
end

function init_v_grid(env::UrbanEnv, vel_res::Float64)
    v_space = get_car_vspace(env, vel_res)
    return RectangleGrid(v_space)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
__precompile__()

module CFET_OrderConditions

using Giac
using Combinatorics

export commutator, LCCC
export gen_exp_tBt_derivatives_coeffs
export gen_C_derivatives_at_0
export gen_C_derivatives_at_0_in_terms_of_A
export gen_exp_tBt_derivatives_at_0_in_terms_of_A
export expand_commutators
export gen_CFET_order_conditions

export number_commutators_of_grade_equal_to
export number_commutators_of_grade_equal_or_less_than

export integer_vector_weighted, lyndon_words, gen_hall_basis_elements_of_grade
export gen_hall_basis_elements_of_grades_leq
export extend_by_rightmost_subwords
export gen_CFET_order_conditions_without_redundancies

const Commutator = Array{Int64,1}
LCCC{T} = Dict{Array{Array{Int64,1},1},T} # Linear Combination of Commutator Chains

function gen_exp_tBt_derivatives_coeffs(N::Integer)
    A = [Dict{Array{Int64,1},Int64}([]=>1),    #Id
         Dict{Array{Int64,1},Int64}([0]=>1)]
    for n = 2:N
        a0 = A[n]
        # a1 = a0*C:
        a1 = [vcat(key, 0)=>val for (key,val) in a0]
        # a1 = a1 + (d/dt)a0:
        for (p,c) in a0
            for j = 1:length(p)
                q = copy(p)
                q[j] += 1
                d = get(a1, q, 0)
                a1[q] = d+c
            end
        end
        push!(A, a1)
    end
    A
end

function gen_C_derivatives_at_0(N::Integer)
    A = Array{Dict{Array{Int64,1},Int64},1}[[Dict(vcat(zeros(Int64,k),1)=>1)] for k=1:N-1]
    for k = 1:N-1
        for l = 1:N-k-1
            a0 = A[k][l]
            a1 = Dict{Array{Int64,1},Int64}()
            for (p,c) in a0
                for j=1:length(p)
                    q = copy(p)
                    q[j] += 1
                    if q[end] != q[end-1]
                        d = get(a1, q, 0)
                        a1[q] = d+c
                    end
                end
            end
            push!(A[k], a1)
        end
    end
    C = [Dict{Array{Int64,1},Int64}([n]=>(n+1)) for n=0:N]
    for n = 2:N
        a1 = C[n+1]
        for k=1:n-1
            f = binomial(n,k+1)
            a0 = A[k][n-k-1+1]
            for (p,c) in a0
                d = get(a1, p, 0//1)
                a1[p] = d+f*c
            end
        end
        C[n+1] = a1
    end
    C
end

function gen_C_derivatives_at_0_in_terms_of_A{T}(N::Integer, b::Array{T,1})
    C_derivatives_at_0=gen_C_derivatives_at_0(N)
    C1 = [Dict{Array{Int64,1},T}() for n=0:N]
    for n = 0:N
        for (p,coeff) in C_derivatives_at_0[n+1]
            C1[n+1][p] = coeff * prod([b[p[j]+1]  for j =1:length(p)])
        end
    end
    C1
end

Id(T::Type; f=one(T)) = LCCC{T}(Array[Int64[]]=>f)
is_c_number(x) = length(x)==1 && length(first(keys(x))[1])==0

function mult{T}(x::LCCC{T}, y::Dict{Array{Int64,1},T})
    r = LCCC{T}()
    for (p,c) in x
        for (q, d) in y
            pq = copy(p)
            push!(pq, q)
            r[pq] = c*d
        end
    end
    r
end

function mult{T}(x::LCCC{T}, y::LCCC{T})
    if is_c_number(x)
        c = first(values(x))
        return LCCC{T}([p=>c*d for (p,d) in y])
    elseif is_c_number(y)
        c = first(values(y))
        return LCCC{T}([p=>c*d for (p,d) in x])
    else
        r = LCCC{T}()
        for (p,c) in x
            for (q, d) in y
                pq = vcat(p,q)
                r[pq] = c*d
            end
        end
        return r
    end
end

function add{T}(x::LCCC{T}, y::LCCC{T};
                fx=one(T), fy=one(T))
    r = copy(x)
    for (p,c) in y
        d = get(r, p, 0)
        r[p] = fy*c + fx*d
    end
    r
end

function gen_exp_tBt_derivatives_at_0_in_terms_of_A{T}(N::Integer, b::Array{T,1})
    C_derivatives_at_0_in_terms_of_A=gen_C_derivatives_at_0_in_terms_of_A(N-1, b)
    exp_tBt_derivatives_coeffs=gen_exp_tBt_derivatives_coeffs(N)
    r = [LCCC{T}(Array[Int64[]]=>one(T))] # Id
    for n = 1:N
        y = LCCC{T}()
        for (p,d) in exp_tBt_derivatives_coeffs[n+1]
            x = LCCC{T}(
            [Array[p1]=>d*d1 for (p1,d1) in C_derivatives_at_0_in_terms_of_A[p[1]+1]])
            for i = 2:length(p)
                x = mult(x, C_derivatives_at_0_in_terms_of_A[p[i]+1])
            end
            y = add(y, x)
        end
        push!(r,y)
    end
    r
end

function expand_commutators{T}(a::Array{Int64,1}, f::T; coeff_type=T)
    if length(a)>=3
        b = LCCC{coeff_type}(Array[[a[1]]]=>one(T))
        c = expand_commutators(a[2:end], f, coeff_type=coeff_type)
        return add(mult(b,c), mult(c,b), fy=-one(T))
    elseif length(a)==2
        return LCCC{coeff_type}(Array[[a[1]], [a[2]]]=>f, Array[[a[2]], [a[1]]]=>-f)
    else
        return LCCC{coeff_type}(Array[a]=>f)
    end
end

function expand_commutators{T}(a::LCCC{T})
    b = LCCC{T}()
    for (p,c) in a
        q = expand_commutators(p[1], one(T), coeff_type=T)
        for i = 2:length(p)
            q = mult(q, expand_commutators(p[i], one(T), coeff_type=T))
        end
        b = add(b, q, fy=c)
    end
    b
end

multinomial_coeff(q::Int, k::Array{Int,1}) = div(factorial(q), prod([factorial(i) for i in k]))

function gen_CFET_order_conditions(N::Integer,J::Integer)
    b = [giac[giac(string("b",j,"_",n)) for n=0:N] for j=1:J]
    C = gen_exp_tBt_derivatives_at_0_in_terms_of_A(N, b[1])
    C = [expand_commutators(c) for c in C]
    C = [LCCC{giac}([key=>factor(val) for (key,val) in c]) for c in C]
    CC = Array[C]
    for j=2:J
        push!(CC, [LCCC{giac}([key=>subst(val,b[1],b[j]) for (key,val) in c]) for c in C])
    end
    Y = Dict{Array{Int64,1},giac}[]
    for q=0:N-1
        y = LCCC{giac}()
        for k0 in Combinatorics.WithReplacementCombinations(1:J,q+1)
            k = zeros(Int,J)
            for j in k0
                k[j]+=1
            end
            coeff = multinomial_coeff(q+1, k)
            x = Id(giac, f=giac(coeff))
            for j = J:-1:1
                x = mult(x, CC[j][k[j]+1])
            end
            y = add(y, x)
        end
        for k0 in Combinatorics.WithReplacementCombinations(1:J+1,q)
            k = zeros(Int,J+1)
            for j in k0
                k[j]+=1
            end
            coeff = multinomial_coeff(q, k)
            x = LCCC{giac}(Array[[k[J+1]]]=>coeff)
            for j = J:-1:1
                x = mult(x, CC[j][k[j]+1])
            end
            y = add(y, x, fy=-one(giac))

        end
        push!(Y, Dict{Array{Int64,1},giac}([[k[1] for k in key]=>factor(val) for (key,val) in y]))
    end
    Y
end

function number_commutators_of_grade_equal_to(n::Integer)
    if n==1
        return 1
    end
    f = factor(n)
    p = collect(keys(f))
    k = collect(values(f))
    nu = 0
    for l = 0:(2^length(p)-1)
        lb = bin(l,length(p))
        mu = (-1)^count_ones(l)
        pot = prod([p[j]^(k[j]-(lb[j]=='0'?0:1)) for j=1:length(p)])
        nu += mu*2^pot
    end
    div(nu,n)
end

number_commutators_of_grade_equal_or_less_than(n::Integer) =
    sum([num_of_eqs(k) for k=1:n])

# The following functions are stolen from MuPAD-Combinat,
# see http://mupad-combinat.sourceforge.net

# mupad/COMBINAT/words.mu/evaluationTable:
function evaluation_table(word::Vector)
    t = Dict{Any, Int}()
    for l in word
        if haskey(t, l)
            t[l] += 1
        else
            t[l] = 1
        end
    end
    t
end

# mupad/COMBINAT/words.mu/standard
function standard(word::Vector)
    t = evaluation_table(word)
    kk = sort(collect(keys(t)))
    offset = 0
    tmp = 0
    for k in kk
        tmp = t[k]
        t[k] = offset
        offset += tmp
    end
    r = zeros(Int, length(word))
    i = 1
    for l in word
        t[l] += 1
        r[i] = t[l]
        i += 1
    end
    r
end

# mupad/COMBINAT/integerVectorsWeighted.mu/
# (There it is a local function inside function list
function recfun(n::Int, l1::Vector{Int})
    result = []
    l = copy(l1)
    w = pop!(l)
    if length(l)==0
        d = div(n, w)
        if mod(n, w)==0
            return Vector{Int}[[d]]
        else
            return Vector{Int}[]
        end
    end
    for d = div(n, w):-1:0
        result = vcat(result, [vcat(x,d) for x in recfun(n-d*w, l)])
    end
    result
end

# mupad/COMBINAT/permutations.mu/mult2
function mult2(x::Vector{Int}, y::Vector{Int})
    nx = length(x)
    ny = length(y)
    # complete with fixed points if the permutations are not
    # of the same length
    x = vcat(x, nx+1:ny)
    y = vcat(y, ny+1:nx)
    result = x
    for i=1:length(x)
        result[y[i]] = x[i]
    end
    result
end

# mupad/COMBINAT/integerVectorsWeighted.mu/list
function integer_vector_weighted(n::Int, l::Vector{Int})
    if length(l)==0
        if n==0
            return Vector{Int}[[]]
        else
            return Vector{Int}[]
        end
    end

    # Heuristic: we sort the weights increasingly to limit
    # the backtracking, and reorder back the result afterward
    perm = standard(l)
    sort!(l)

    [mult2(x, perm) for x in recfun(n,l)]
end

# For the following 2 functions, see:
#   Sawada, Joe.  "A fast algorithm to generate necklaces with fixed content",
#   Source Theoretical Computer Science archive Volume 301 , Issue 1-3 (May 2003)
# and sage/src/sage/combinat/necklace.py

function simple_fixed_content1(res::Vector{Vector{Int}}, a::Vector{Int},
            content::Vector{Int}, t::Int, p::Int, k::Int)
    n = length(a)
    if t>n
        if n==p
            push!(res, a)
        end
    else
        for j=a[t-p]:k-1
            if content[j+1]!=0
                a[t] = j
                content[j+1] -= 1
                if j==a[t-p]
                    simple_fixed_content1(res, copy(a), content, t+1, p, k)
                else
                    simple_fixed_content1(res, copy(a), content, t+1, t, k)
                end
                content[j+1] += 1
            end
        end
    end
end

function simple_fixed_content(content::Vector{Int})
    a = zeros(Int, sum(content))
    content[1] -= 1
    k = length(content)
    res = Vector{Int}[]
    simple_fixed_content1(res, a, content, 2, 1, k)
    res
end

#The following function is stolen from sage/combinat/necklace.py

function lyndon_words(e::Vector{Int})
    k=1
    while e[k]==0
        k += 1
    end
    Vector{Int}[[i+k for i in z] for z in simple_fixed_content(e[k:end])]
end

#The following is stolen from http://mathoverflow.net/questions/97703/list-of-hall-basis

function gen_hall_basis_elements_of_grade(g::Int)
    L = integer_vector_weighted(g, collect(1:g))
    res = Vector{Int}[]
    for l in L
        for w in lyndon_words(l)
            push!(res, w-1)
        end
    end
    res
end

gen_hall_basis_elements_of_grades_leq(p) = vcat([gen_hall_basis_elements_of_grade(k) for k=1:p]...)

function gen_CFET_order_conditions_without_redundancies(N::Integer,J::Integer)
    C0 = gen_CFET_order_conditions(N, J)
    Dict{Array{Int, 1}, giac}[Dict{Array{Int, 1}, giac}(
        [b=>C0[k][b] for b in gen_hall_basis_elements_of_grade(k)]) for k=1:length(C0)]
end

coeff_rhs(w::Array{Int64,1}) = 1//prod([sum(w[j:end]+1) for j=1:length(w)])

function extend_by_rightmost_subwords(W::Array{Array{Int64,1},1})
    WW=Dict{Array{Int64,1},Int}(Int64[]=>1)
    for w in W
        for l=1:length(w)
            WW[w[l:end]] = 1
        end
    end
    return collect(keys(WW))
end

function gen_CFET_order_conditions(W::Array{Array{Int64,1},1}, J::Integer)
    W1 = extend_by_rightmost_subwords(W)

    qmax = maximum([maximum(w) for w in W])
    b = giac[giac(string("b",j, "_", q)) for q=0:qmax, j=1:J]
    bj = giac[giac(string("bj_", q)) for q=0:qmax]

    m=length(W1)
    M=zeros(giac, m, m)
    for k=1:m
        for l=1:m
            r = length(W1[k])-length(W1[l])
            if r>=0 && W1[k][r+1:end]==W1[l]
                w = W1[k][1:r]
                M[k,l]  = factor(prod(bj[w+1]))/factorial(length(w))
            end
        end
    end

    c = map(x->subst(x, bj, b[:,1]), M[:,findfirst(W1, Int64[])])
    for j=2:J
        Mj = map(x->subst(x, bj, b[:,j]), M)
        c = Mj*c
    end

    c = [c[findfirst(W1, w)] for w in W]

    eqs = [ c[j]-coeff_rhs(W[j]) for j=1:length(W) ]
    vars = vcat(b...)
    eqs, vars
end

end # module


text:
# Method Declarations
# ===================

# Method declaration macro.
macro declare(names)
    @assert names.head == :tuple
    base_names = Set(Base.names(Base))
    ex = Expr(:block)
    for name in names.args
        if name ∈ base_names
            error("method name '$(name)' is already used in Base")
        end
        push!(ex.args, :(function $(name) end))
        push!(ex.args, :(export $(name)))
    end
    return esc(ex)
end

# Base functions used in (and exported from) BioJulia packages.
@declare (
    distance,
    seqname,
    hasseqname,
    sequence,
    hassequence,
    metadata,
    leftposition,
    hasleftposition,
    rightposition,
    hasrightposition,
    isoverlapping,
    metainfoval,
    metainfotag,
    header,
    isfilled
)


text:
#Communication avoiding LU

import Base: start, next, done, size, getindex
import DistributedArrays: DArray

#              |
#   A(ma x na) | B(ma x nb)
# -------------+------------
#   C(mb x na) | D(mb x nb)

immutable RowBlockedMatrix{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    A::S
    m::Int
    n::Int
    k::Int
end

size(A::RowBlockedMatrix) = ceil(Int, A.m/A.k), 1
size(A::RowBlockedMatrix, i::Int) = if i==1 ceil(Int, A.m/A.k) elseif i==2 1 else error() end

function getindex(A::RowBlockedMatrix, i, j=1)
    rowstart = (i-1)*A.k + 1
    rowend = min(i*A.k, A.m)
    view(A.A, rowstart:rowend, 1:A.n)
end

RowBlockedMatrix{T}(A::AbstractMatrix{T}, k) = RowBlockedMatrix{T,typeof(A)}(A, size(A, 1), size(A, 2), k)

blockrows(A, k) = RowBlockedMatrix(A, k)

#RowBlockedMatrix is iterable over blocked rows
start(A::RowBlockedMatrix) = 1
next(A::RowBlockedMatrix, iblock) = (A[iblock], iblock+1)
done(A::RowBlockedMatrix, iblock) = iblock == 1+ceil(Int, A.m/A.k)

function whoowns{T,N}(A::SubArray{T,N,DArray{T,2,Matrix{T}}})
    rowrange = A.indexes[1]
    colrange = A.indexes[2]

    #Find out which pid owns the most data
    myid = myndata = 0
    paridx = A.parent.indexes
    for i = 1:size(paridx, 1), j = 1:size(paridx, 2)
        localrowrange = rowrange ∩ paridx[i,j][1]
        localcolrange = colrange ∩ paridx[i,j][2]

        ndata = length(localrowrange)*length(localcolrange)
        if ndata > myndata
            myid = A.parent.pids[i,j]
            myndata = ndata
        end
    end
    #info("whoown darray: $myid ($(myndata/(length(rowrange)*length(colrange))))")
    return myid
end

#Tall and skinny LU
#k - row block size
#b - number of pivot rows
function tslu!(A, piv, k, b)
    #Step 1: Find set of good pivot rows
    #Serial LU on each block row
    t = time()
    #DEBUG info("Time $(round(time() - t, 3)): Step 1: broadcast")
    lus = []
    for blockrow in blockrows(A, k)
        push!(lus, @spawnat whoowns(blockrow) lufact!(Array(blockrow)))
    end
    #DEBUG info("Time $(round(time() - t, 3)): Step 1: collect")
    #lus = map(fetch, lus)
    #map(wait, lus) #Wait but don't copy data back to master
    for l in lus
        if isa(l,RemoteException)
            showerror(STDERR, l)
            error()
        end
    end
    #DEBUG info("Time $(round(time() - t, 3)): Step 1: collect done")
    while length(lus) > 1
        #DEBUG info("Time $(round(time() - t, 3)): Step 1: recursion with $(length(lus)) work to do")
        newlus = []
        nlu = length(lus)
        npairs, isodd = divrem(nlu, 2)
        for ipair in 1:npairs
            w = @spawnat lus[2ipair-1].where begin
                Upair = [fetch(lus[2ipair-1])[:U]; fetch(lus[2ipair])[:U]]
                lufact!(Upair)
            end
            if isa(w,RemoteException)
                showerror(STDERR, w)
                error()
            end
            push!(newlus, w)
            #push!(newlus, @spawnat(@show(whoowns(blockrow)), lufact(blockrow)))
        end
        if isodd==1
            #push!(newlus, @spawnat(@show(whoowns(blockrow)), lufact(blockrow)))
            push!(newlus, lus[end])
        end
        #DEBUG info("$(length(lus)) -> $(length(newlus))")
        lus = newlus
    end
    #DEBUG info("Time $(round(time() - t, 3)): Step 1: recursion done")

    #Step 2: Permute pivot rows into first b rows of the panel
    #DEBUG info("Time $(round(time() - t, 3)): Step 2: pivot")
    perm = view(LinAlg.ipiv2perm(fetch(lus[1])[:p], size(A, 2)), 1:min(b, size(A, 2)))
    permuterows!(A, perm)
    #DEBUG info("Time $(round(time() - t, 3)): Step 2: pivot done")

    #Step 3: Unpivoted LU on panel
    #DEBUG info("Time $(round(time() - t, 3)): Step 3: LU on panel")
    F = if b ≥ size(A, 2)
        lufact!(A, Val{false})
    else
        lufact!(view(A, :, 1:b), Val{false})
    end
    #DEBUG info("Time $(round(time() - t, 3)): Step 3: LU done")
end

function pidmap{T}(A::SubArray{T,2,DArray{T,2,Matrix{T}}})
    pids = Dict()

	#XXX use dispatch to make the subrowrange, subcolrange computations cleaner
    subrowrange = A.indexes[1]
    if isa(subrowrange, Colon)
         subrowrange = 1:size(A, 1)
    elseif isa(subrowrange, Integer)
         subrowrange = subrowrange:subrowrange
    end

    subcolrange = A.indexes[2]
    if isa(subcolrange, Colon)
        subcolrange = 1:size(A, 2)
    elseif isa(subcolrange, Integer)
        subcolrange = subcolrange:subcolrange
    end

    for (p, (localrowrange, localcolrange)) in zip(A.parent.pids, A.parent.indexes)
        #Compute intersection in global address space
        rows = subrowrange ∩ localrowrange
        cols = subcolrange ∩ localcolrange
        if length(rows) > 0 && length(cols) > 0
            #Compute blocks in process-local address space
            localrows = rows-localrowrange.start+1
            localcols = cols-localcolrange.start+1

            subrows = rows-subrowrange.start+1
            subcols = cols-subcolrange.start+1

            pids[p] = rows, cols, localrows, localcols, subrows, subcols

            #@assert A.parent[rows, cols] ==
            #fetch(@spawnat p localpart(A.parent)[localrows, localcols]) ==
            #A[subrows, subcols] #Not implemented
        end
    end
    pids
end

function Base.lufact!{T,piv}(
        A::SubArray{T,2,DArray{T,2,Matrix{T}}},
        ::Type{Val{piv}}
    )

    #Copy to local memory
    tmpA = Array(A)

    #Do computation locally
    lufact!(tmpA, Val{piv})

    #Redistribute results
    pids = pidmap(A)
    @sync for (p, (gr, gc, lr, lc, sr, sc)) in pids
        @spawnat p localpart(A.parent)[lr, lc] = view(tmpA, sr, sc)
    end
end

function Base.A_ldiv_B!{T}(
        F::LinAlg.LU{T, DArray{T, 2, Array{T,2}}},
        b::AbstractVector{T}
    )

    Ff = Array(F.factors)
    U = UpperTriangular(Ff)
    L = LinAlg.UnitLowerTriangular(Ff)
    bp = b[LinAlg.ipiv2perm(F.ipiv, length(b))]
    x = L \ bp
    return U \ x
end

#Intesection of index ranges
function Base.intersect(A::Tuple{UnitRange{Int64},UnitRange{Int64}},
                        B::Tuple{UnitRange{Int64},UnitRange{Int64}})
    a1, a2 = A
    b1, b2 = B
    return (a1∩b1, a2∩b2)
end

function DistributedArrays.localpart{T,N}(A::SubArray{T,N,DArray{T,2,Matrix{T}}})
    subindexes = A.indexes
    local globallocalindex, globalsubindex
    for (p, localindexes) in zip(A.parent.pids, A.parent.indexes)
        if p == myid()
            #Compute intersection of indexes in global address space
            globallocalindex = localindexes
            globalsubindex = (localindexes ∩ subindexes)
            break
        end
    end

    #Convert local index from global address space to process-local address space
    localindex1 = globalsubindex[1] - globallocalindex[1].start + 1
    localindex2 = globalsubindex[2] - globallocalindex[2].start + 1
    return view(localpart(A.parent), localindex1, localindex2)
end

#Swap rows of a matrix in place
function swaprows!(A::AbstractMatrix, r1, r2)
    nc = size(A, 2)
    for i=1:nc
        A[r1, i], A[r2, i] = A[r2,i], A[r1,i]
    end
end

#Swap rows of a DArray in place
# e.g. A = drandn(4000, 4000)
#      swaprows!(A, [(1,2), (1000,1001)])
# swaps the first two rows of A in place, and also rows 1000 and 1001
function swaprows!{T}(A::SubArray{T,2,DArray{T,2,Matrix{T}}}, r1, r2)
    pidmap1 = pidmap(view(A, r1:r1, :))
    pidmap2 = pidmap(view(A, r2:r2, :))

    @sync for (i,((p1, (gr1, gc1, lr1, lc1, sr1, sc1)),
         (p2, (gr2, gc2, lr2, lc2, sr2, sc2)))) in enumerate(zip(pidmap1, pidmap2))

        #Need to look up which parts in pidmap2 align
        # info("FROM proc $p1, chunk $gr1,$gc1 ($lr1, $lc1) $sr1, $sc1")
        # info("TO   proc $p2, chunk $gr2,$gc2 ($lr2, $lc2) $sr2, $sc2")
            if p1 == p2 #Do all the work locally
                # info("Local work")
                @assert lc2 == lc2
                @spawnat p1 swaprows!(localpart(A), sr1, sr2)
            else #copy from remote, swap, send
                # info("Local work $p2 -> $p1")
                @spawnat p1 begin
                     lA = localpart(A)
                     B = Array(view(A, sr1, sc1))

                     for (i1, i2) in enumerate(sr1), j in sr2
                         B[i1, j], lA[i2, j] = lA[i2, j], B[i1, j]
                     end

                     @spawnat p2 localpart(A)[sr2, sc2] = B
                end
        end
    end
    #DEBUG info("swaprows! done")
end

function permuterows!{T}(A::SubArray{T,2,DArray{T,2,Matrix{T}}}, perm)
    # info("permuterows! with permutation $perm")
    for (a, b) in enumerate(perm)
        if a!=b
            swaprows!(A, a, b)
        end
    end
end

function calu!(A, k=64, b=64)
    m, n = size(A)
    piv = collect(1:m)
    for i = 1:ceil(Int, n/k)
        #DEBUG info("calu!: tslu! on [$((1+(i-1)*k):m), $((1+(i-1)*k):n)]")
        tslu!(view(A, (1+(i-1)*k):m, (1+(i-1)*k):n), piv, k, b)
    end
    LinAlg.LU(A, piv, 0)
end


text:
module MultigridBase

export AffineMap, LinearProblem, compute_residual!, coarsen_problem!, jsmooth!
export cartesian_1d_crs, cartesian_2d_crs
export MultigridSolver, add_layer!, vcycle!

using LinearAlgebra, SparseArrays

include("./base.jl")
include("./cartesian.jl")
include("./solver.jl")

end


text:
using HypothesisTests, Base.Test

@test_throws DimensionMismatch HypothesisTests.check_same_length([1], [])
@test_throws ArgumentError HypothesisTests.check_alpha(0.0)

type TestTest <: HypothesisTests.HypothesisTest end
result = HypothesisTests.population_param_of_interest(TestTest())
@test result[1] == "not implemented yet"
@test isnan(result[2])


text:
"""
$(TYPEDSIGNATURES)

Generates the Liouville transformed set of ODEs, which is the original
ODE system with a new variable `trJ` appended, corresponding to the
-tr(Jacobian). This variable is used for properties like uncertainty
propagation from a given initial distribution density.

For example, if ``u'=p*u`` and `p` follows a probability distribution
``f(p)``, then the probability density of a future value with a given
choice of ``p`` is computed by setting the inital `trJ = f(p)`, and
the final value of `trJ` is the probability of ``u(t)``.

Example:

```julia
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t α β γ δ
@variables x(t) y(t)
D = Differential(t)

eqs = [D(x) ~ α*x - β*x*y,
       D(y) ~ -δ*y + γ*x*y]

sys = ODESystem(eqs)
sys2 = liouville_transform(sys)
@variables trJ

u0 = [x => 1.0,
      y => 1.0,
      trJ => 1.0]

prob = ODEProblem(sys2,u0,tspan,p)
sol = solve(prob,Tsit5())
```

Where `sol[3,:]` is the evolution of `trJ` over time.

Sources:

Probabilistic Robustness Analysis of F-16 Controller Performance: An
Optimal Transport Approach

Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya
https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf
"""
function liouville_transform(sys::AbstractODESystem)
      t = get_iv(sys)
      @variables trJ
      D = ModelingToolkit.Differential(t)
      neweq = D(trJ) ~ trJ*-tr(calculate_jacobian(sys))
      neweqs = [equations(sys);neweq]
      vars = [states(sys);trJ]
      ODESystem(neweqs,t,vars,parameters(sys))
end


text:
module Coriolis

export
    FPlane, ConstantCartesianCoriolis, BetaPlane, NonTraditionalBetaPlane,
    HydrostaticSphericalCoriolis, VectorInvariantEnergyConserving, VectorInvariantEnstrophyConserving,
    x_f_cross_U, y_f_cross_U, z_f_cross_U

using Printf
using Oceananigans.Grids
using Oceananigans.Operators

# Physical constants for constructors.
const Ω_Earth = 7.292115e-5 # [s⁻¹] https://en.wikipedia.org/wiki/Earth%27s_rotation#Angular_speed
const R_Earth = 6371.0e3    # Mean radius of the Earth [m] https://en.wikipedia.org/wiki/Earth

"""
    AbstractRotation

Abstract supertype for parameters related to background rotation rates.
"""
abstract type AbstractRotation end

include("no_rotation.jl")
include("f_plane.jl")
include("constant_cartesian_coriolis.jl")
include("beta_plane.jl")
include("non_traditional_beta_plane.jl")
include("hydrostatic_spherical_coriolis.jl")

end


text:
#!/usr/bin/env julia

import REPL: REPLCompletions

const snipmate_dir = "snippets"
const snipmate_filename = "_.snippets"
const ultisnips_dir = "UltiSnips"
const ultisnips_filename = "all_unicode.snippets"

mkpath(snipmate_dir)
mkpath(ultisnips_dir)

snipmate_file = open("$(snipmate_dir)/$(snipmate_filename)", "w")
ultisnips_file = open("$(ultisnips_dir)/$(ultisnips_filename)", "w")

for file in (snipmate_file, ultisnips_file)
    println(file,
            "# This file was generated by `$(basename(Base.source_path()))`")
    println(file, "# using Julia $VERSION")
    println(file)
end

for (keyword, unicode) in sort!(collect(Iterators.flatten(
        (REPLCompletions.latex_symbols,
         REPLCompletions.emoji_symbols))),
        by=x -> x[2])
    println(snipmate_file, "snippet ", keyword, " ", unicode)
    println(snipmate_file, "\t", unicode)

    println(ultisnips_file, "snippet ", keyword, " \"", unicode, "\" i")
    println(ultisnips_file, unicode)
    println(ultisnips_file, "endsnippet")
end

close(snipmate_file)
close(ultisnips_file)

max_stars_repo_name: danielwe/vim-unicode-snippets

max_issues_repo_name: danielwe/vim-unicode-snippets

max_forks_repo_name: danielwe/vim-unicode-snippets