text:
description="""
Schoenauer vector triade benchmark
See e.g. https://blogs.fau.de/hager/archives/tag/benchmarking

Run multithreaded test e.g. with
JULIA_NUM_THREADS=4 julia --multithread-threads

Run multiporcess test e.g. with
julia -p 4 --multiprocess-distributed

"""

using ArgParse
using Distributed

using Printf

Distributed.@everywhere using SharedArrays
Distributed.@everywhere using LoopVectorization
################################################################
# Helper methods

# GFlops for vector triad
GFlops(N)=N*2.0/1.0e9

# Partition the range 1:N for a number of tasks
function partition(N,ntasks)
    loop_begin=zeros(Int64,ntasks)
    loop_end=zeros(Int64,ntasks)
    for itask=1:ntasks
        ltask=Int(floor(N/ntasks))
        loop_begin[itask]=(itask-1)*ltask+1
        if itask==ntasks # adjust last task length
            ltask=N-(ltask*(ntasks-1))
        end
        loop_end[itask]=loop_begin[itask]+ltask-1
    end
    return (loop_begin,loop_end)
end

# Create arrays for test
function make_arrays(N)
    a = Array{Float64,1}(undef,N)
    b = Array{Float64,1}(undef,N)
    c = Array{Float64,1}(undef,N)
    d = Array{Float64,1}(undef,N)

    Threads.@threads for i=1:N
        a[i]=i
        b[i]=N-i
        c[i]=i
        d[i]=-i
    end
    return(a,b,c,d)
end

# Create shared arrays for test
function make_shared_arrays(N)
    (_a,_b,_c,_d)=make_arrays(N)
    a = SharedArray(_a)
    b = SharedArray(_b)
    c = SharedArray(_c)
    d = SharedArray(_d)
    return(a,b,c,d)
end

# Create an array of vector lengths of size nrun
# starting with N0 and increasing in geometrica progression.
# ppomag denotes the number of
# elements per order of magnitude.
# Try to keep full powers of 10 at once.
function vsizes(;N0=1000,ppomag=8,nrun=41)
    vsz=[N0]
    N=N0
    N0*=10
    for irun=1:nrun-1
        N=N*10^(1.0/ppomag)
        if (irun%ppomag==0)
            N=N0
            N0*=10
        end
        push!(vsz,Int(ceil(N)))
    end
    return vsz
end

##################################################################################
# Tests

# Scalar operation
# In all cases, the triad is run nrepeat times, in order
# to use the same measurement method as in the C code.
function vtriad_scalar(N,nrepeat)
    (a,b,c,d)=make_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
            @inbounds @fastmath  for i=1:N
                d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Same with avx
function vtriad_scalar_avx(N,nrepeat)
    (a,b,c,d)=make_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
            @avx for i=1:N
                d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# fork-join style operation using  Threads.@threads
function vtriad_multithread_threads(N,nrepeat)
    (a,b,c,d)=make_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
             Threads.@threads for i=1:N
                 @inbounds @fastmath  d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Kernel for spawn bases operation
Distributed.@everywhere function _kernel_avx(a,b,c,d,n0,n1)
    @avx for i=n0:n1
        d[i]=a[i]+b[i]*c[i]
    end
    return 0
end

function _kernel(a,b,c,d,n0,n1)
    @inbounds @fastmath  for i=n0:n1
        d[i]=a[i]+b[i]*c[i]
    end
    return 0
end

# Operation using  Threads.@spawn
function vtriad_multithread_spawn(N,nrepeat)
    (a,b,c,d)=make_arrays(N)
    ntasks=Threads.nthreads()
    loop_begin,loop_end=partition(N,ntasks)
    t=@elapsed begin
        for j=1:nrepeat
            mapreduce(task->fetch(task),+,[Threads.@spawn _kernel(a,b,c,d,loop_begin[i],loop_end[i]) for i=1:ntasks])
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

function threadcall_fetch(f,args...)
    t=Threads.@spawn f(args...)
    fetch(t)
end

# Operation using  Threads.@spawn (slower than multithead_spawn...)
function vtriad_multithread_spawn_asyncmap(N,nrepeat)
    (a,b,c,d)=make_arrays(N)
    ntasks=Threads.nthreads()
    loop_begin,loop_end=partition(N,ntasks)
    function spwn(a,b,c,d,n0,n1)
        t=Threads.@spawn _kernel(a,b,c,d,n0,n1)
        fetch(t)
    end
    t=@elapsed begin
        for j=1:nrepeat
            res=asyncmap(i->threadcall_fetch(_kernel,a,b,c,d,loop_begin[i],loop_end[i]),1:ntasks)
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Run scalar benchmark using SharedArrays
function vtriad_scalar_shared(N,nrepeat)
    (a,b,c,d)=make_shared_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
            @inbounds @fastmath  for i=1:N
                d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Run scalar benchmark using SharedArrays
function vtriad_scalar_shared_avx(N,nrepeat)
    (a,b,c,d)=make_shared_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
            @avx for i=1:N
                d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Multiprocessing using Distributed.@spawn and SharedArrays
function vtriad_multiprocess_spawn(N,nrepeat)
    (a,b,c,d)=make_shared_arrays(N)

    ntasks=Distributed.nprocs()
    loop_begin,loop_end=partition(N,ntasks)

    t=@elapsed begin
        for j=1:nrepeat
            mapreduce(task->fetch(task),+,[Distributed.@spawn _kernel_avx(a,b,c,d,loop_begin[i],loop_end[i]) for i=1:ntasks])
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

# Multiprocessing using Distributed.@distributed and SharedArrays
function vtriad_multiprocess_distributed(N,nrepeat)
    (a,b,c,d)=make_shared_arrays(N)
    t=@elapsed begin
        for j=1:nrepeat
            Distributed.@sync Distributed.@distributed for i=1:N
                @avx  d[i]=a[i]+b[i]*c[i]
            end
        end
    end
    GC.gc()
    return [N,GFlops(N*nrepeat)/t]
end

##################################################################################
# Test driver

function run_vtriad(vtriad;N0=1000,ppomag=4,nrun=41,flopcount=5.0e8)
    # Create vector sizes
    vsz=vsizes(N0=N0,ppomag=ppomag,nrun=nrun)

    # Create result array
    result=zeros(2,length(vsz))

    # Run test with increasing sizes, adapting nrepeat as to keep the
    # number of operations constant and time a large enough batch
    for i=1:length(vsz)
        result[:,i].=vtriad(vsz[i],Int(ceil(flopcount/vsz[i])))
    end
    return result
end

function main(ARGS)

    settings = ArgParseSettings()
    settings.preformatted_description=true
    settings.description=description
    @add_arg_table settings begin

        "--nrun"   help= "number of runs"
        default=41

        "--flopcount"  help= "number of flops per run"
        default=5.0e8

        "--N0"  help= "smallest array size"
        default=1000

        "--ppomag"  help= "points per order of magnitude"
        default=8

        "--scalar" help= "scalar vtriad"
        action=:store_true

        "--scalar-avx" help= "scalar vtriad+avx"
        action=:store_true

        "--scalar-shared" help= "scalar vtriad with shared arrays"
        action=:store_true

        "--scalar-shared-avx" help= "scalar vtriad with shared arrays+avx"
        action=:store_true

        "--multithread-threads" help= "Threads.@threads"
        action=:store_true

        "--multithread-spawn"  help= "Threads.@spawn"
        action=:store_true

        "--multiprocess-spawn"  help= "Distributed.@spawn with SharedArrays"
        action=:store_true

        "--multiprocess-distributed" help= "Distributed.@distrubuted with SharedArrays"
        action=:store_true

    end
    parsed_args=parse_args(ARGS, settings)

    flopcount=parsed_args["flopcount"]
    ppomag=parsed_args["ppomag"]
    N0=parsed_args["N0"]
    nrun=parsed_args["nrun"]

    # Run test
    if parsed_args["multiprocess-spawn"]
        @printf("# multiprocess-spawn nprocs=%d\n",Distributed.nprocs())
        result=run_vtriad(vtriad_multiprocess_spawn;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["multiprocess-distributed"]
        @printf("# multiprocess-distributed nprocs=%d\n",Distributed.nprocs())
        result=run_vtriad(vtriad_multiprocess_distributed;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["scalar"]
        @printf("# scalar\n")
        result=run_vtriad(vtriad_scalar;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["scalar-avx"]
        @printf("# scalar-avx\n")
        result=run_vtriad(vtriad_scalar_avx;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["scalar-shared"]
        @printf("# scalar-shared\n")
        result=run_vtriad(vtriad_scalar_shared;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["scalar-shared-avx"]
        @printf("# scalar-shared-avx\n")
        result=run_vtriad(vtriad_scalar_shared_avx;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["multithread-spawn"]
        @printf("# multithread-spawn nthreads=%d\n",Threads.nthreads())
        result=run_vtriad(vtriad_multithread_spawn;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    elseif parsed_args["multithread-threads"]
        @printf("# multithread-threads nthreads=%d\n",Threads.nthreads())
        result=run_vtriad(vtriad_multithread_threads;N0=N0,ppomag=ppomag,nrun=nrun,flopcount=flopcount)
    else
        return
    end

    # Print result
    @printf("# Array size  GFlops/s\n")
    for i=1:size(result,2)
        @printf("% 10d %6.3f\n",result[1,i], result[2,i])
    end
end

main(ARGS)

max_stars_repo_name: j-fu/julia-tests

max_issues_repo_name: j-fu/julia-tests

max_forks_repo_name: j-fu/julia-tests


text:
#==============================================================================#
# AWSSDK.jl
#
# This file is generated by AWSAPI.jl from service decriptions at:
# https://github.com/aws/aws-sdk-js/tree/master/apis
#==============================================================================#

module AWSSDK

include("MigrationHub.jl")
include("ACM.jl")
include("ACMPCA.jl")
include("AlexaForBusiness.jl")
include("APIGateway.jl")
include("ApplicationAutoScaling.jl")
include("AppStream.jl")
include("AppSync.jl")
include("Athena.jl")
include("AutoScaling.jl")
include("AutoScalingPlans.jl")
include("Batch.jl")
include("Budgets.jl")
include("CostExplorer.jl")
include("Cloud9.jl")
include("CloudDirectory.jl")
include("CloudFormation.jl")
include("CloudFront.jl")
include("CloudHSM.jl")
include("CloudHSMV2.jl")
include("CloudSearch.jl")
include("CloudSearchDomain.jl")
include("CloudTrail.jl")
include("CodeBuild.jl")
include("CodeCommit.jl")
include("CodeDeploy.jl")
include("CodePipeline.jl")
include("CodeStar.jl")
include("CognitoIdentity.jl")
include("CognitoIdentityServiceProvider.jl")
include("CognitoSync.jl")
include("Comprehend.jl")
include("ConfigService.jl")
include("Connect.jl")
include("CUR.jl")
include("DataPipeline.jl")
include("DAX.jl")
include("DeviceFarm.jl")
include("DirectConnect.jl")
include("Discovery.jl")
include("DLM.jl")
include("DMS.jl")
include("DirectoryService.jl")
include("DynamoDB.jl")
include("EC2.jl")
include("ECR.jl")
include("ECS.jl")
include("EKS.jl")
include("ElastiCache.jl")
include("ElasticBeanstalk.jl")
include("EFS.jl")
include("ELB.jl")
include("ELBv2.jl")
include("EMR.jl")
include("ElasticTranscoder.jl")
include("SES.jl")
include("MarketplaceEntitlementService.jl")
include("ES.jl")
include("CloudWatchEvents.jl")
include("Firehose.jl")
include("FMS.jl")
include("GameLift.jl")
include("Glacier.jl")
include("Glue.jl")
include("Greengrass.jl")
include("GuardDuty.jl")
include("Health.jl")
include("IAM.jl")
include("ImportExport.jl")
include("Inspector.jl")
include("Iot.jl")
include("IotData.jl")
include("IoTJobsDataPlane.jl")
include("IoT1ClickDevicesService.jl")
include("IoT1ClickProjects.jl")
include("IoTAnalytics.jl")
include("Kinesis.jl")
include("KinesisVideoArchivedMedia.jl")
include("KinesisVideoMedia.jl")
include("KinesisAnalytics.jl")
include("KinesisVideo.jl")
include("KMS.jl")
include("Lambda.jl")
include("LexModelBuildingService.jl")
include("Lightsail.jl")
include("CloudWatchLogs.jl")
include("MachineLearning.jl")
include("Macie.jl")
include("MarketplaceCommerceAnalytics.jl")
include("MediaConvert.jl")
include("MediaLive.jl")
include("MediaPackage.jl")
include("MediaStore.jl")
include("MediaStoreData.jl")
include("MediaTailor.jl")
include("MarketplaceMetering.jl")
include("Mobile.jl")
include("MobileAnalytics.jl")
include("CloudWatch.jl")
include("MQ.jl")
include("MTurk.jl")
include("Neptune.jl")
include("OpsWorks.jl")
include("OpsWorksCM.jl")
include("Organizations.jl")
include("PI.jl")
include("Pinpoint.jl")
include("Polly.jl")
include("Pricing.jl")
include("RDS.jl")
include("Redshift.jl")
include("Rekognition.jl")
include("ResourceGroups.jl")
include("ResourceGroupsTaggingAPI.jl")
include("Route53.jl")
include("Route53Domains.jl")
include("LexRuntime.jl")
include("SageMakerRuntime.jl")
include("S3.jl")
include("SageMaker.jl")
include("SimpleDB.jl")
include("SecretsManager.jl")
include("ServerlessApplicationRepository.jl")
include("ServiceCatalog.jl")
include("ServiceDiscovery.jl")
include("Shield.jl")
include("SMS.jl")
include("Snowball.jl")
include("SNS.jl")
include("SQS.jl")
include("SSM.jl")
include("StepFunctions.jl")
include("StorageGateway.jl")
include("DynamoDBStreams.jl")
include("STS.jl")
include("Support.jl")
include("SWF.jl")
include("TranscribeService.jl")
include("Translate.jl")
include("WAF.jl")
include("WAFRegional.jl")
include("WorkDocs.jl")
include("WorkMail.jl")
include("WorkSpaces.jl")
include("XRay.jl")

end # module AWSSDK

#==============================================================================#
# End of file
#==============================================================================#

ext: jl


text:
using LESFilter
using InplaceRealFFT
using Base.Test

#=
Test function: f = sin(3x)+cos(3y)+sin(4z)
GaussianFilter(Δ,f) = exp(-2Δ^2/3) * (exp(7Δ^2/24) * (cos(3y) + sin(3x)) + sin(4z))
=#

lx = 1.0
ly = 1.0
lz = 1.0

nx = 32
ny = 32
nz = 32

x = reshape(linspace(0,lx*2π*(1-1/nx),nx),(nx,1,1))
y = reshape(linspace(0,ly*2π*(1-1/ny),ny),(1,ny,1))
z = reshape(linspace(0,lz*2π*(1-1/nz),nz),(1,1,nz))

field = @. sin(3*x) + cos(3*y) + sin(4*z)

correct(Δ::Real) = @. exp(-2Δ^2/3) * (exp(7Δ^2/24) * (cos(3y) + sin(3x)) + sin(4z))
dx = x[2]
for df in 1:0.5:10
  Δ = df*dx
  @test lesfilter(field,"G",Δ,(lx,ly,lz)) ≈ correct(Δ)
end

field2 = PaddedArray(field)

for df in 1:0.5:10
  Δ = df*dx
  @test lesfilter!(field2,"G",Δ,(lx,ly,lz)) ≈ correct(Δ)
  copy!(real(field2),field)
end

rfft!(field2)
@inferred LESFilter.loopgaussian!(complex(field2),rfftfreq(nx,lx),fftfreq(ny,ly),fftfreq(nz,lz),8*dx)


text:
using Primes: isprime

PrimesGen() = Iterators.filter(isprime, Iterators.countfrom(Int64(2)))

print("Sum of first 100,000 primes:  ")
println(Iterators.sum(Iterators.take(PrimesGen(), 100000)))
print("First 20 primes:  ( ")
foreach((p->print(p," ")), Iterators.take(PrimesGen(), 20))
println(")")
print("Primes between 100 and 150:  ( ")
for p in Iterators.filter((p->p>=100), PrimesGen()) p > 150 && break; print(p, " ") end
println(")")
let cnt = 0
    for p in PrimesGen()
        p > 8000 && break; if p > 7700 cnt += 1 end
    end; println("Number of primes between 7700 and 8000:  ", cnt)
end
println("The 10,000th prime:  ", Iterators.first(Iterators.drop(PrimesGen(), 9999)))
println()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DataStructures

queue = Queue(String)
@show enqueue!(queue, "foo")
@show enqueue!(queue, "bar")
@show dequeue!(queue) # -> foo
@show dequeue!(queue) # -> bar


text:
using Distributions, LinearAlgebra, Manifolds, Random, Statistics, StatsBase, Test

@testset "DiracDelta" begin
    M = Euclidean(3)
    @testset "constructor" begin
        p = randn(3)
        d = DiracDelta(M, p)
        @test d isa DiracDelta
        @test d isa Manifolds.MPointDistribution{typeof(M)}
        @test d.point == p
        @test d.manifold === M
    end

    @testset "eltype" begin
        p = randn(3)
        d = DiracDelta(M, p)
        @test eltype(d) === typeof(p)
    end

    @testset "rand!" begin
        p = randn(3)
        q = Vector{Float64}(undef, 3)
        d = DiracDelta(M, p)
        @test rand!(d, q) === q
        @test q == p
    end

    @testset "rand" begin
        p = randn(3)
        d = DiracDelta(M, p)
        @test rand(d) == p

        ps = rand(d, 2)
        @test typeof(ps) == Vector{typeof(p)}
        @test length(ps) == 2
        @test ps[1] == p
        @test ps[2] == p
    end

    @testset "insupport" begin
        p = randn(3)
        d = DiracDelta(M, p)
        @test insupport(d, p)
        @test !insupport(d, randn(3))
    end

    @testset "$f" for f in (mean, median, mode)
        p = randn(3)
        d = DiracDelta(M, p)
        @test f(d) == p
    end

    @testset "$f" for f in (var, std)
        p = randn(3)
        d = DiracDelta(M, p)
        @inferred f(d)
        @test isreal(f(d))
        @test iszero(f(d))

        M2 = Euclidean(3; field=ℂ)
        p2 = randn(ComplexF64, 3)
        d2 = DiracDelta(M2, p2)
        @inferred f(d2)
        @test isreal(f(d2))
        @test iszero(f(d2))
    end

    @testset "convolve" begin
        G = SpecialOrthogonal(3)
        p1, p2 = ntuple(_ -> group_exp(G, hat(G, Matrix(Diagonal(ones(3))), randn(3))), 2)
        d1 = DiracDelta(G, p1)
        d2 = DiracDelta(G, p2)
        d12 = convolve(d1, d2)
        @test d12 isa typeof(d1)
        @test isapprox(G, d12.point, compose(G, p1, p2))
    end

    @testset "inversion" begin
        p = randn(3)
        d = DiracDelta(M, p)
        @test inversion(d) === d
    end
end


text:
const units = Dict(  "k" => 1.0e3,
               "meg" => 1.0e6,
               "g" => 1.0e9,
               "t" => 1.0e12,
               "m" => 1.0e-3,
               "u" => 1.0e-6,
               "μ" => 1.0e-6, #965
               "µ" => 1.0e-6, #181
               "n" => 1.0e-9,
               "p" => 1.0e-12,
               "f" => 1.0e-15,
               )

const stringunits = Dict(  "k" => "*1e3",
               "meg" => "*1e6",
               "g" => "*1e9",
               "t" => "*1e12",
               "m" => "*1e-3",
               "u" => "*1e-6",
               "μ" => "*1e-6", #965
               "µ" => "*1e-6", #181
               "n" => "*1e-9",
               "p" => "*1e-12",
               "f" => "*1e-15",
               )


text:
abstract type AbstractNode end

#mutable struct Node{T} <: AbstractNode
#    parent::Union{Node{T},Nothing}
#    children::Union{Vector{Node{T}}, Nothing}
#    data::T
#
#    Node{T}() where T = new{T}(nothing, nothing, nothing)
#    Node{T}(data::T) where T = new{T}(nothing, nothing, data)
#end

struct BoxTreeOptions <: TreeOptions
    nmin
    maxlevel
end

function BoxTreeOptions(; nmin=1, maxlevel=64)
    return BoxTreeOptions(nmin, maxlevel)
end

mutable struct BoxTreeNode{T} <: AbstractNode
    parent::Union{BoxTreeNode{T},Nothing}
    children::Union{Vector{BoxTreeNode{T}}, Nothing}
    level::Integer
    boundingbox::BoundingBox
    data::T
end

function BoxTreeNode{T}(bbox::BoundingBox, data::T) where T
    return BoxTreeNode{T}(nothing, nothing, 0, bbox, data)
end

function BoxTreeNode{T}(bbox::BoundingBox, level, data::T) where T
    return BoxTreeNode{T}(nothing, nothing, level, bbox, data)
end

function add_child!(parent::BoxTreeNode, bbox::BoundingBox, data::T) where T
    childnode = BoxTreeNode(nothing, nothing, parent.level + 1, bbox, data)
    childnode.parent = parent
    if parent.children === nothing
        parent.children = [childnode]
    else
        push!(parent.children, childnode)
    end
end

function create_tree(points::Vector{SVector{D,T}}; treeoptions=BoxTreeOptions()) where {D,T}
    root = BoxTreeNode( nothing,
                        nothing,
                        0,
                        BoundingBox(points),
                        Vector(1:length(points)))

    fill_tree!(root, points, nmin=treeoptions.nmin, maxlevel = treeoptions.maxlevel)

    return root
end

function fill_tree!(node::BoxTreeNode, points::Vector{SVector{D,T}};
                    nmin = 1, maxlevel=-log2(eps(eltype(points)))) where {D,T}

    if length(node.data) <= nmin || node.level >= maxlevel - 1
        return
    end

    nchildren = D == 2 ? 4 : 8 # Number of (possible) children
    sorted_points = zeros(eltype(node.data), length(node.data)+1, nchildren)

    for pindex in node.data
        id = whichchildbox(node.boundingbox, points[pindex])
        sorted_points[1, id] +=  1
        sorted_points[1 + sorted_points[1, id], id] = pindex
    end

    for i = 1:nchildren
        if sorted_points[1, i] > 0
            childpointindices = sorted_points[2:(sorted_points[1, i]+1), i]
            add_child!(node, getchildbox(node.boundingbox, i), childpointindices)
            fill_tree!(node.children[end], points, nmin=nmin, maxlevel=maxlevel)
        end
    end
end

max_stars_repo_name: JoshuaTetzner/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

# Abaqus test file:
# *NODE
#  1, 0.0, 0.0, 0.0
#  2, 3.0, 0.0, 0.0
#  3, 0.0, 0.0, -3.0
#  4, 0.0, -4.0, 0.0
# *ELEMENT, TYPE=B33, ELSET=ELSET1
#  1, 1,2
# *ELEMENT, TYPE=B33, ELSET=ELSET2
#  2, 1,3
# *ELEMENT, TYPE=B33, ELSET=ELSET3
#  3, 1,4
# *BEAM GENERAL SECTION, ELSET=ELSET1, DENSITY=100.0, SECTION=GENERAL
#  2.5E-2, 1.0E-4, 0.0E-4, 1.5E-4, 3.0E-5
#  0.0, 1.0, 0.0
#  100.0E6, 80.0E6
# *BEAM GENERAL SECTION, ELSET=ELSET2, DENSITY=150.0, SECTION=GENERAL
#  1.5E-2, 1.2E-4, 0.0E-4, 2.0E-4, 4.0E-5
#  1.0, 0.0, 0.0
#  120.0E6, 90.0E6
# *BEAM GENERAL SECTION, ELSET=ELSET3, DENSITY=200.0, SECTION=GENERAL
#  1.8E-2, 1.3E-4, 0.0E-4, 1.4E-4, 5.0E-5
#  1.0, 0.0, 0.0
#  150.0E6, 70.0E6
# *BOUNDARY
#  2, 1,6
#  3, 1,6
#  4, 1,6
# *STEP
# *MATRIX GENERATE, STIFFNESS, MASS, LOAD
# *MATRIX OUTPUT, STIFFNESS, MASS, LOAD
# *DLOAD
#  1, P1, 1.0E3
#  2, P2, 2.0E3
#  3, PX, 1.0E3
#  3, PY, 2.0E3
#  3, PZ, 3.0E3
# *END STEP
# *STEP
# *STATIC
# *CLOAD
#  1,1, 1.0E3
#  1,2, 2.0E3
#  1,3, 3.0E3
#  1,4, 4.0E3
#  1,5, 5.0E3
#  1,6, 6.0E3
# *DLOAD
#  1, P1, 1.0E3
#  2, P2, 2.0E3
#  3, PX, 1.0E3
#  3, PY, 2.0E3
#  3, PZ, 3.0E3
# *NODE PRINT
#  U, COORD
# *END STEP

using FEMBeam
using FEMBase.Test

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [3.0, 0.0, 0.0],
    3 => [0.0, 0.0, -3.0],
    4 => [0.0, -4.0, 0.0])

beam1 = Element(Seg2, [1, 2])
update!(beam1, "geometry", X)
update!(beam1, "youngs modulus", 100.0e6)
update!(beam1, "shear modulus", 80.0e6)
update!(beam1, "cross-section area", 2.5e-2)
update!(beam1, "torsional moment of inertia 1", 1.0e-4)
update!(beam1, "torsional moment of inertia 2", 1.5e-4)
update!(beam1, "polar moment of inertia", 3.0e-5)
update!(beam1, "distributed load 1", 1.0e3)

beam2 = Element(Seg2, [1, 3])
update!(beam2, "geometry", X)
update!(beam2, "youngs modulus", 120.0e6)
update!(beam2, "shear modulus", 90.0e6)
update!(beam2, "cross-section area", 1.5e-2)
update!(beam2, "torsional moment of inertia 1", 1.2e-4)
update!(beam2, "torsional moment of inertia 2", 2.0e-4)
update!(beam2, "polar moment of inertia", 4.0e-5)
update!(beam2, "distributed load 2", 2.0e3)

beam3 = Element(Seg2, [1, 4])
update!(beam3, "geometry", X)
update!(beam3, "youngs modulus", 150.0e6)
update!(beam3, "shear modulus", 70.0e6)
update!(beam3, "cross-section area", 1.8e-2)
update!(beam3, "torsional moment of inertia 1", 1.3e-4)
update!(beam3, "torsional moment of inertia 2", 1.4e-4)
update!(beam3, "polar moment of inertia", 5.0e-5)
update!(beam3, "distributed load x", 1.0e3)
update!(beam3, "distributed load y", 2.0e3)
update!(beam3, "distributed load z", 3.0e3)

load1 = Element(Poi1, [1])
update!(load1, "point force 1",  1.0e3)
update!(load1, "point force 2",  2.0e3)
update!(load1, "point force 3",  3.0e3)
update!(load1, "point moment 1", 4.0e3)
update!(load1, "point moment 2", 5.0e3)
update!(load1, "point moment 3", 6.0e3)

update!(beam1, "normal", [0.0, 1.0, 0.0])
update!(beam2, "normal", [1.0, 0.0, 0.0])
update!(beam3, "normal", [1.0, 0.0, 0.0])

problem = Problem(Beam, "example 1", 6)
add_elements!(problem, [beam1, beam2, beam3, load1])
time = 0.0
assemble!(problem, time)
K = full(problem.assembly.K)[1:6,1:6]
f = full(problem.assembly.f)[1:6]
u = K\f
println("u = $u")
println("f = $f")


text:
module ReverseDiffCompat

using ..ReverseDiff: ReverseDiff, @grad, value, track, TrackedReal, TrackedVector,
    TrackedMatrix
using Requires, LinearAlgebra

using ..Bijectors: Log, SimplexBijector, maphcat, simplex_link_jacobian,
    simplex_invlink_jacobian, simplex_logabsdetjac_gradient, ADBijector,
    ReverseDiffAD, Inverse
import ..Bijectors: _eps, logabsdetjac, _logabsdetjac_scale, _simplex_bijector,
    _simplex_inv_bijector, replace_diag, jacobian, getpd, lower,
    _inv_link_chol_lkj, _link_chol_lkj

using Compat: eachcol
using Distributions: LocationScale

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:ReverseDiffAD}, Inverse{<:ADBijector{<:ReverseDiffAD}}},
    x::Real
)
    return ReverseDiff.gradient(x -> b(x[1]), [x])[1]
end
function jacobian(
    b::Union{<:ADBijector{<:ReverseDiffAD}, Inverse{<:ADBijector{<:ReverseDiffAD}}},
    x::AbstractVector{<:Real}
)
    return ReverseDiff.jacobian(b, x)
end

_eps(::Type{<:TrackedReal{T}}) where {T} = _eps(T)
function Base.minimum(d::LocationScale{<:TrackedReal})
    m = minimum(d.ρ)
    if isfinite(m)
        return d.μ + d.σ * m
    else
        return m
    end
end
function Base.maximum(d::LocationScale{<:TrackedReal})
    m = maximum(d.ρ)
    if isfinite(m)
        return d.μ + d.σ * m
    else
        return m
    end
end

logabsdetjac(b::Log{1}, x::Union{TrackedVector, TrackedMatrix}) = track(logabsdetjac, b, x)
@grad function logabsdetjac(b::Log{1}, x::AbstractVector)
    return -sum(log, value(x)), Δ -> (nothing, -Δ ./ value(x))
end
@grad function logabsdetjac(b::Log{1}, x::AbstractMatrix)
    return -vec(sum(log, value(x); dims = 1)), Δ -> (nothing, .- Δ' ./ value(x))
end
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Val{0})
    return track(_logabsdetjac_scale, a, value(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::Real, v::Val{0})
    return _logabsdetjac_scale(value(a), value(x), Val(0)), Δ -> (inv(value(a)) .* Δ, nothing, nothing)
end
# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Val{0})
    return track(_logabsdetjac_scale, a, value(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, v::Val{0})
    da = value(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, value(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Val{0})
    return track(_logabsdetjac_scale, a, value(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, v::Val{0})
    da = value(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, value(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    return track(_logabsdetjac_scale, a, value(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, v::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    da = value(a)
    J = inv.(da)
    return _logabsdetjac_scale(da, value(x), Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    return track(_logabsdetjac_scale, a, value(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, v::Val{1})
    da = value(a)
    Jᵀ = repeat(inv.(da), 1, size(x, 2))
    return _logabsdetjac_scale(da, value(x), Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end
function _simplex_bijector(X::Union{TrackedVector, TrackedMatrix}, b::SimplexBijector{1})
    return track(_simplex_bijector, X, b)
end
@grad function _simplex_bijector(Y::AbstractVector, b::SimplexBijector{1})
    Yd = value(Y)
    return _simplex_bijector(Yd, b), Δ -> (simplex_link_jacobian(Yd)' * Δ, nothing)
end
@grad function _simplex_bijector(Y::AbstractMatrix, b::SimplexBijector{1})
    Yd = value(Y)
    return _simplex_bijector(Yd, b), Δ -> begin
        maphcat(eachcol(Yd), eachcol(Δ)) do c1, c2
            simplex_link_jacobian(c1)' * c2
        end, nothing
    end
end

function _simplex_inv_bijector(X::Union{TrackedVector, TrackedMatrix}, b::SimplexBijector{1})
    return track(_simplex_inv_bijector, X, b)
end
@grad function _simplex_inv_bijector(Y::AbstractVector, b::SimplexBijector{1})
    Yd = value(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> (simplex_invlink_jacobian(Yd)' * Δ, nothing)
end
@grad function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector{1})
    Yd = value(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> begin
        maphcat(eachcol(Yd), eachcol(Δ)) do c1, c2
            simplex_invlink_jacobian(c1)' * c2
        end, nothing
    end
end

replace_diag(::typeof(log), X::TrackedMatrix) = track(replace_diag, log, X)
@grad function replace_diag(::typeof(log), X)
    Xd = value(X)
    f(i, j) = i == j ? log(Xd[i, j]) : Xd[i, j]
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = i == j ? ∇[i, j]/Xd[i, j] : ∇[i, j]
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

replace_diag(::typeof(exp), X::TrackedMatrix) = track(replace_diag, exp, X)
@grad function replace_diag(::typeof(exp), X)
    Xd = value(X)
    f(i, j) = ifelse(i == j, exp(Xd[i, j]), Xd[i, j])
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = ifelse(i == j, ∇[i, j]*exp(Xd[i, j]), ∇[i, j])
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

logabsdetjac(b::SimplexBijector{1}, x::Union{TrackedVector, TrackedMatrix}) = track(logabsdetjac, b, x)
@grad function logabsdetjac(b::SimplexBijector{1}, x::AbstractVector)
    xd = value(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, simplex_logabsdetjac_gradient(xd) * Δ)
    end
end
@grad function logabsdetjac(b::SimplexBijector{1}, x::AbstractMatrix)
    xd = value(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, maphcat(eachcol(xd), Δ) do c, g
            simplex_logabsdetjac_gradient(c) * g
        end)
    end
end

getpd(X::TrackedMatrix) = track(getpd, X)
@grad function getpd(X::AbstractMatrix)
    Xd = value(X)
    return LowerTriangular(Xd) * LowerTriangular(Xd)', Δ -> begin
        Xl = LowerTriangular(Xd)
        return (LowerTriangular(Δ' * Xl + Δ * Xl),)
    end
end
lower(A::TrackedMatrix) = track(lower, A)
@grad function lower(A::AbstractMatrix)
    Ad = value(A)
    return lower(Ad), Δ -> (lower(Δ),)
end

function find_alpha(wt_y::T, wt_u_hat::T, b::T) where {T<:TrackedReal}
    return track(find_alpha, wt_y, wt_u_hat, b)
end
@grad function find_alpha(wt_y::TrackedReal, wt_u_hat::TrackedReal, b::TrackedReal)
    α = find_alpha(data(wt_y), data(wt_u_hat), data(b))

    ∂wt_y = inv(1 + wt_u_hat * sech(α + b)^2)
    ∂wt_u_hat = - tanh(α + b) * ∂wt_y
    ∂b = ∂wt_y - 1
    find_alpha_pullback(Δ::Real) = (Δ * ∂wt_y, Δ * ∂wt_u_hat, Δ * ∂b)

    return α, find_alpha_pullback
end

end


text:
export Descent

mutable struct Descent
    eta::Float64
end

Descent() = Descent(0.1)

function apply!(o::Descent, x, ∆)
     ∆ .*= o.eta
end


text:
using Base: HasEltype, HasLength, EltypeUnknown

mutable struct LinkedNode{T}
    val::T
    next::LinkedNode{T}
    prev::LinkedNode{T}

    # All new created nodes are self referential only
    function LinkedNode{T}(k::T) where {T}
        x = new{T}(k)
        x.next = x
        x.prev = x
        return x
    end
end

mutable struct CyclicOrderedSet{T} <: AbstractSet{T}
    first::Union{LinkedNode{T}, Nothing}
    length::Int

    CyclicOrderedSet{T}() where {T} = new{T}(nothing, 0)
end

CyclicOrderedSet(itr) = _CyclicOrderedSet(itr, IteratorEltype(itr))

_CyclicOrderedSet(itr, ::HasEltype) = CyclicOrderedSet{eltype(itr)}(itr)
function _CyclicOrderedSet(itr, ::EltypeUnknown)
    T = Base.@default_eltype(itr)
    (isconcretetype(T) || T === Union{}) || return grow_to!(CyclicOrderedSet{T}(), itr)
    return CyclicOrderedSet{T}(itr)
end
CyclicOrderedSet{T}(itr) where {T} = union!(CyclicOrderedSet{T}(), itr)

function Base.iterate(s::CyclicOrderedSet, state = s.first)
    if state === nothing
        return nothing
    else
        return state.val, state.next == s.first ? nothing : state.next
    end
end
Base.IteratorSize(::Type{<:CyclicOrderedSet}) = HasLength()
Base.IteratorEltype(::Type{<:CyclicOrderedSet}) = HasEltype()
Base.length(s::CyclicOrderedSet) = s.length
Base.isempty(s::CyclicOrderedSet) = length(s) == 0

Base.empty(s::CyclicOrderedSet{T}, ::Type{U}=T) where {T,U} = CyclicOrderedSet{U}()

function Base.show(io::IO, s::CyclicOrderedSet{T}) where {T}
    print(io, "CyclicOrderedSet{", T, "}(")
    if s.first !== nothing
        f = s.first
        show(io, f.val)
        n = f.next
        while n !== f
            print(io, ", ")
            show(io, n.val)
            n = n.next
        end
    end
    print(io, ")")
end

function _findnode(s::CyclicOrderedSet, x)
    isempty(s) && return nothing
    n = s.first
    while !isequal(n.val, x)
        n = n.next
        n == s.first && return nothing
    end
    return n
end

# adding items
function _push!(s::CyclicOrderedSet{T}, n::LinkedNode{T}) where {T}
    if isempty(s)
        s.first = n
    else
        list = s.first
        n.next = list
        n.prev = list.prev
        list.prev.next = n
        list.prev = n
    end
    s.length += 1
    return s
end
function Base.push!(s::CyclicOrderedSet{T}, x) where T
    n = _findnode(s, x)
    if n === nothing
        _push!(s, LinkedNode{T}(x))
    else
        _delete!(s, n)
        _push!(s, n)
    end
    return s
end
Base.pushfirst!(s::CyclicOrderedSet, x) = rotate!(push!(s, x))

# removing items
function _delete!(s::CyclicOrderedSet{T}, n::LinkedNode{T}) where {T}
    n.next.prev = n.prev
    n.prev.next = n.next
    s.length -= 1
    if n == s.first
        s.first = s.length == 0 ? nothing : n.next
    end
    return s
end
function Base.pop!(s::CyclicOrderedSet)
    isempty(s) && throw(ArgumentError("collection must be non-empty"))
    n = s.first.prev
    _delete!(s, n)
    return n.val
end
function Base.pop!(s::CyclicOrderedSet, x)
    n = _findnode(s, x)
    if n === nothing
        throw(KeyError(x))
    else
        _delete!(s, n)
        return x
    end
end
function Base.pop!(s::CyclicOrderedSet, x, default)
    n = _findnode(s, x)
    if n === nothing
        return default
    else
        _delete!(s, n)
        return x
    end
end
function Base.delete!(s::CyclicOrderedSet, x)
    n = _findnode(s, x)
    if n !== nothing
        _delete!(s, n)
    end
    return s
end
function Base.empty!(s::CyclicOrderedSet)
    s.first = nothing
    s.length = 0
    return s
end

# Rotate one step forward, so last element is now first
function rotate!(s::CyclicOrderedSet)
    if length(s) > 1
        s.first = s.first.prev
    end
    return s
end

function _move_to_front!(s::CyclicOrderedSet{T}, n::LinkedNode{T}) where {T}
    if s.first !== n
        n.next.prev = n.prev
        n.prev.next = n.next
        n.next = s.first
        n.prev = s.first.prev
        s.first.prev.next = n
        s.first.prev = n
        s.first = n
    end
    return s
end


text:
module TestDeprecated

using MLJBase
using ..Models
using Test
using Statistics
using MLJScientificTypes
using StableRNGs

rng = StableRNG(698790187)

@load KNNRegressor

NN = 7
X = MLJBase.table(rand(rng,NN, 3));
y = 2X.x1 - X.x2 + 0.05*rand(rng,NN);
Xs = source(X); ys = source(y)

broadcast_mode(v) = mode.(v)

t = MLJBase.table
m = MLJBase.matrix
f = FeatureSelector()
h = OneHotEncoder()
k = KNNRegressor()
u = UnivariateStandardizer()
c = ConstantClassifier()

# build a linear network for training:
N = MLJBase.linear_learning_network(Xs, ys, nothing, u, nothing, f, k)
fM = machine(f, Xs)
Xt = transform(fM, Xs)

# build the same network by hand:
uM = machine(u, ys)
yt = transform(uM, ys)
kM = machine(k, Xt, yt)
zhat = predict(kM, Xt)

k.K = 3; f.features = [:x3,]
fit!(N)
hand_built = N();

## SIMPLE SUPERVISED PIPELINE WITH TARGET TRANSFORM

# test a simple pipeline prediction agrees with prediction of
# hand-built learning network built earlier:
p = @pipeline(Pipe(sel=FeatureSelector(),
                   knn=KNNRegressor(),
                   target=UnivariateStandardizer()))

p.knn.K = 3; p.sel.features = [:x3,]
mach = machine(p, X, y)
fit!(mach)
@test MLJBase.tree(mach.fitresult.predict).arg1.model.K == 3
MLJBase.tree(mach.fitresult.predict).arg1.arg1.model.features == [:x3, ]
@test predict(mach, X) ≈ hand_built

# test a simple probabilistic classifier pipeline:
X = MLJBase.table(rand(rng,7,3))
y = categorical(collect("ffmmfmf"))
Xs = source(X)
ys = source(y)
p = @pipeline(Pipe21(hot=OneHotEncoder(),
                    cnst=ConstantClassifier()),
              prediction_type=:probabilistic)
mach = machine(p, X, y)
fit!(mach)
@test p isa ProbabilisticComposite
pdf(predict(mach, X)[1], 'f') ≈ 4/7

# test a simple deterministic classifier pipeline:
X = MLJBase.table(rand(rng,7,3))
y = categorical(collect("ffmmfmf"))
Xs = source(X)
ys = source(y)
p = @pipeline(Piper3(hot=OneHotEncoder(), cnst=ConstantClassifier(),
                     broadcast_mode))
mach = machine(p, X, y)
fit!(mach)
@test predict(mach, X) == fill('f', 7)

# test pipelines with weights:
w = map(y) do η
    η == 'm' ? 100 : 1
end
mach = machine(p, X, y, w)
fit!(mach)
@test predict(mach, X) == fill('m', 7)

# test a pipeline with static transformation of target:
NN = 100
X = (x1=rand(rng,NN), x2=rand(rng,NN), x3=categorical(rand(rng,"abc", NN)));
y = 1000*abs.(2X.x1 - X.x2 + 0.05*rand(rng,NN))
# by hand:
Xs =source(X); ys = source(y)
hot = OneHotEncoder()
hot_=machine(hot, Xs)
W = transform(hot_, Xs)
sel = FeatureSelector(features=[:x1,:x3__a])
sel_ = machine(sel, W)
Wsmall = transform(sel_, W)
z = log(ys)
knn = KNNRegressor(K=4)
knn_ = machine(knn, Wsmall, z)
zhat = predict(knn_, Wsmall)
yhat = exp(zhat)
fit!(yhat)
pred1 = yhat()
# with pipeline:
p = @pipeline Pipe4(hot=OneHotEncoder(),
                    sel=FeatureSelector(),
                    knn=KNNRegressor(),
                    target=v->log.(v),
                    inverse=v->exp.(v))
p.sel.features = [:x1, :x3__a]
p.knn.K = 4
p_ = machine(p, X, y)
fit!(p_)
pred2 = predict(p_, X)
@test pred1 ≈ pred2

# and another:
age = [23, 45, 34, 25, 67]
X = (age = age,
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
height = [67.0, 81.5, 55.6, 90.0, 61.1]
p = @pipeline Pipe9(X -> coerce(X, :age=>Continuous),
                    hot = OneHotEncoder(),
                    knn = KNNRegressor(K=3),
                    target = UnivariateStandardizer())
fit!(machine(p, X, height))

# STATIC TRANSFORMERS IN PIPELINES

struct MyTransformer <: Static
    ftr::Symbol
end

MLJBase.transform(transf::MyTransformer, verbosity, X) =
    selectcols(X, transf.ftr)

p99 = @pipeline Pipe99(X -> coerce(X, :age=>Continuous),
                       hot = OneHotEncoder(),
                       selector = MyTransformer(:age))

mach  = machine(p99, X) |> fit!

@test transform(mach, X) == float.(X.age)

end
true


text:
export DataVariable, datavar, update!

import Base: show

mutable struct DataVariable{D, S} <: AbstractVariable
    name            :: Symbol
    collection_type :: AbstractVariableCollectionType
    messageout      :: S
    nconnected      :: Int
end

Base.show(io::IO, datavar::DataVariable) = print(io, "DataVariable(", indexed_name(datavar), ")")

function datavar(name::Symbol, ::Type{D}, collection_type::AbstractVariableCollectionType = VariableIndividual(); subject::S = RecentSubject(Union{Message{Missing}, Message{D}})) where { S, D }
    return DataVariable{D, S}(name, collection_type, subject, 0)
end

function datavar(name::Symbol, ::Type{D}, dims::Tuple; subject::S = RecentSubject(Union{Message{Missing}, Message{D}})) where { S, D }
    return datavar(name, D, dims...; subject = subject)
end

function datavar(name::Symbol, ::Type{D}, length::Int; subject::S = RecentSubject(Union{Message{Missing}, Message{D}})) where { S, D }
    vars = Vector{DataVariable{D, S}}(undef, length)
    @inbounds for i in 1:length
        vars[i] = datavar(name, D, VariableVector(i); subject = similar(subject))
    end
    return vars
end

function datavar(name::Symbol, ::Type{D}, dims::Vararg{Int}; subject::S = RecentSubject(Union{Message{Missing}, Message{D}})) where { S, D }
    vars = Array{DataVariable{D, S}}(undef, dims)
    @inbounds for i in CartesianIndices(axes(vars))
        vars[i] = datavar(name, D, VariableArray(i); subject = similar(subject))
    end
    return vars
end

Base.eltype(::Type{ <: DataVariable{D} }) where D = D
Base.eltype(::DataVariable{D})            where D = D

degree(datavar::DataVariable)          = nconnected(datavar)
name(datavar::DataVariable)            = datavar.name
collection_type(datavar::DataVariable) = datavar.collection_type
isconnected(datavar::DataVariable)     = datavar.nconnected !== 0
nconnected(datavar::DataVariable)      = datavar.nconnected

getlastindex(::DataVariable) = 1

messageout(datavar::DataVariable, ::Int) = datavar.messageout
messagein(datavar::DataVariable, ::Int)  = error("It is not possible to get a reference for inbound message for datavar")

update!(datavar::DataVariable, ::Missing)           = next!(messageout(datavar, 1), Message(missing, false, false))
update!(datavar::DataVariable, data::Number)        = update!(eltype(datavar), typeof(data), datavar, data)
update!(datavar::DataVariable, data::AbstractArray) = update!(eltype(datavar), typeof(data), datavar, data)

update!(::Type{ PointMass{D} }, ::Type{D}, datavar, data)   where { D }      = next!(messageout(datavar, 1), Message(PointMass(data), false, false))
update!(::Type{ PointMass{D1} }, ::Type{D2}, datavar, data) where { D1, D2 } = error("'$(name(datavar)) = datavar($D1, ...)' accepts data of type $D1, but $D2 has been supplied. Check 'update!($(name(datavar)), data::$D2)' and explicitly convert data to type $D1.")

resend!(datavar::DataVariable) = update!(datavar, Rocket.getrecent(messageout(datavar, 1)))

function update!(datavars::AbstractVector{ <: DataVariable }, data::AbstractVector)
    @assert size(datavars) === size(data) "Invalid update! call: size of datavar array and data should match"
    foreach(zip(datavars, data)) do (var, d)
        update!(var, d)
    end
end

finish!(datavar::DataVariable) = complete!(messageout(datavar, 1))

get_pipeline_stages(::DataVariable) = EmptyPipelineStage()

_getmarginal(datavar::DataVariable)              = datavar.messageout |> map(Marginal, as_marginal)
_setmarginal!(datavar::DataVariable, observable) = error("It is not possible to set a marginal stream for `DataVariable`")
_makemarginal(datavar::DataVariable)             = error("It is not possible to make marginal stream for `DataVariable`")

# Extension for _getmarginal
function Rocket.getrecent(proxy::ProxyObservable{ <: Marginal, S, M }) where { S <: Rocket.RecentSubjectInstance, D, M <: Rocket.MapProxy{D, typeof(as_marginal)} }
    return as_marginal(Rocket.getrecent(proxy.proxied_source))
end

function setmessagein!(datavar::DataVariable, ::Int, messagein)
    datavar.nconnected += 1
    return nothing
end


text:
using DifferentialEquations, ReservoirComputing, Surrogates, ModelingToolkit, Plots, MLJLinearModels, SparseArrays, ProgressBars

function ab(C::Array{Float64,1},V)
	# eq (13) from original paper
	return (C[1]*exp(C[2]*(V+C[3]))+C[4]*(V+C[5]))/(exp(C[6]*(V+C[3]))+C[7])
end

function Casmoother(Ca; ep=1.5e-10)
	return ep*0.5*(1.0 + tanh(1.0-(Ca/ep))) + Ca*0.5*(1.0 + tanh((Ca/ep)-1.0))
end

function BR!(dx,x,p,t)

	# spatially local currents
	IK = (exp(0.08*(x[1]+53.0)) + exp(0.04*(x[1]+53.0)))
	IK = 4.0*(exp(0.04*(x[1]+85.0)) - 1.0)/IK
	IK = IK+0.2*(x[1]+23.0)/(1.0-exp(-0.04*(x[1]+23.0)))
	IK = 0.35*IK
	Ix = x[3]*0.8*(exp(0.04*(x[1]+77.0))-1.0)/exp(0.04*(x[1]+35.0))
	INa= (4.0*x[4]*x[4]*x[4]*x[5]*x[6] + 0.003)*(x[1]-50.0)
	Is = 0.09*x[7]*x[8]*(x[1]+82.3+13.0287*log(x[2]))#Casmoother(x[2])))

	# these from Beeler & Reuter table:
	ax = ab([ 0.0005, 0.083, 50.0, 0.0, 0.0, 0.057, 1.0],x[1])
	bx = ab([ 0.0013,-0.06 , 20.0, 0.0, 0.0,-0.04 , 1.0],x[1])
	am = ab([ 0.0   , 0.0  , 47.0,-1.0,47.0,-0.1  ,-1.0],x[1])
	bm = ab([40.0   ,-0.056, 72.0, 0.0, 0.0, 0.0  , 0.0],x[1])
	ah = ab([ 0.126 ,-0.25 , 77.0, 0.0, 0.0, 0.0  , 0.0],x[1])
	bh = ab([ 1.7   , 0.0  , 22.5, 0.0, 0.0,-0.082, 1.0],x[1])
	aj = ab([ 0.055 ,-0.25 , 78.0, 0.0, 0.0,-0.2  , 1.0],x[1])
	bj = ab([ 0.3   , 0.0  , 32.0, 0.0, 0.0,-0.1  , 1.0],x[1])
	ad = ab([ 0.095 ,-0.01 , -5.0, 0.0, 0.0,-0.072, 1.0],x[1])
	bd = ab([ 0.07  ,-0.017, 44.0, 0.0, 0.0, 0.05 , 1.0],x[1])
	af = ab([ 0.012 ,-0.008, 28.0, 0.0, 0.0, 0.15 , 1.0],x[1])
	bf = ab([ 0.0065,-0.02 , 30.0, 0.0, 0.0,-0.2  , 1.0],x[1])

	# BR dynamics
	dx[1] = -(IK + Ix + INa + Is)/p[1]
	dx[2] = -10^-7 * Is + 0.07*(10^-7 - x[2])
	dx[3] = ax*(1.0 - x[3]) - bx*x[3]
	dx[4] = am*(1.0 - x[4]) - bm*x[4]
	dx[5] = ah*(1.0 - x[5]) - bh*x[5]
	dx[6] = aj*(1.0 - x[6]) - bj*x[6]
	dx[7] = ad*(1.0 - x[7]) - bd*x[7]
	dx[8] = af*(1.0 - x[8]) - bf*x[8]

	return nothing

end

# define original parameters and initial condition
# initial condition
u0 = [ -60.0,10^-7,0.01,0.01,0.99,0.99,0.01,0.99]
# parameters
p0 = (1.0)
tspan = (0.0, 1e3)
tt = 10.0.^collect(range(-3.0, +3.0; length=129))
modelODESize = length(u0)

# define ODEProblem, optimize it, and solve for original ODE solution at (u0, p0)
prob = ODEProblem(BR!, u0, tspan, p0)
#sol = solve(prob, Rosenbrock23(); abstol=1e-6, reltol=1e-6, saveat=tt)
sys = modelingtoolkitize(prob)
f   = eval(ModelingToolkit.generate_function(sys)[2])
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
prob = ODEProblem(ODEFunction(f, jac=jac), u0, tspan, p0)
sol = solve(prob, Rosenbrock23(); abstol=1e-10, reltol=1e-8, saveat=tt)

reservoirSize = 300
iterative = reservoirSize > 300 # this is kind of machine-specific?

W   = ReservoirComputing.init_reservoir_givendeg(reservoirSize, 0.8, 6)
Win = ReservoirComputing.init_dense_input_layer(reservoirSize, modelODESize, 0.1)
if iterative
	W = sparse(W)
end

function CTESN!(du, u, p, t)
	du .= tanh.(p[1]*u .+ p[2]*p[3](t))
	return nothing
end
esnprob = ODEProblem(CTESN!, randn(Float64, reservoirSize), tspan, (W, Win, sol))
r = solve(esnprob; saveat = sol.t)

@assert isapprox(sol.t, r.t)

# compute first Wout
Wout = zeros(Float64, modelODESize, reservoirSize)
function fitData!(Wout, xData, rData; beta = 0.01, iterative = iterative)
	fitModel = RidgeRegression(beta)
	#=
	Note:
		rData = hcat(r(sol.t)...)
		xData = hcat(sol(sol.t)...)
	But we avoid re-interpolating them every time by requiring the r solution use `saveat=sol.t`
	=#
	for n=1:modelODESize
		tmp = transpose(fit(fitModel, transpose(rData[:,:]), xData[n,:]; solver=Analytical(iterative=iterative)))
		Wout[n,:] .= tmp[1:reservoirSize]
	end
	return nothing
end
fitData!(Wout, sol[:,:], r[:,:]; beta=0.01)

# define lower and upper bounds for each variable
u_lower = minimum(sol, dims=2)
u_upper = maximum(sol, dims=2)
#u_lower = u0 .+ 1e-1.*abs.(u0)
#u_upper = u0 .- 1e-1.*abs.(u0)

# select set of Sobol-sampled state vectors
u	= sample(1000, u_lower, u_upper, SobolSample()) # this generates tuples; convert to arrays in prob_func
Wout 	= [zeros(Float64, modelODESize, reservoirSize) for _ in 1:length(u)]

# evaluate the ESN at each p, returning WOut(u)
function fitSamples!(Wout, u)

	function prob_func(prob, i, repeat)
		return remake(prob, u0=[u[i][j] for j=1:8])
	end
	ensprob = EnsembleProblem(prob, prob_func = prob_func)
	sim = solve(ensprob, Rosenbrock23(), EnsembleThreads(); trajectories=length(u), abstol=1e-10, reltol=1e-8, saveat=sol.t)
	for n in ProgressBar(1:length(u))
		#print("Fitting u[$(n)]...   \t")
		fitData!(Wout[n], sim[n], r; beta = 1e-6)
		#print("Done.\n")
	end
	return nothing
end
fitSamples!(Wout, u)

# form interpolant over sampled parameter vectors
WoutInterpolant = RadialBasis(u, Wout, u_lower, u_upper)

# define unseen p̂ ensure it is not in p
û= u0
issampled = false;
for n in 1:length(u), m in 1:length(û)
	if isapprox(û[m], u[n][m])
		global issampled = true
		break
	end
end
Ŵout = reshape(WoutInterpolant(û), size(Wout[1]))
x = Ŵout*r
if issampled
	û= u0 .+ 1e-3 .* randn(size(u0)) .* u0
	# form Wout(p̂) and approximate x(t;p̂) from WoutInterpolant at unseen p̂
	Ŵout = reshape(WoutInterpolant(û), size(Wout[1]))
	x = Ŵout*r
	#xInterpolant = RadialBasis(collect(transpose(sol.t)), x, [tspan[begin]], [tspan[end]])

	sol_hat = solve(remake(prob, u0=û), Rosenbrock23(); abstol=1e-6, reltol=1e-6, saveat=sol.t)

end

pp = []
for n=1:length(u0)
	plt = plot(sol, vars=(0,n))
	if issampled; plot!(plt, sol_hat, vars=(0,n), linestyle=:dash, linewidth=2); end
	plot!(plt, sol.t, x[n,:], linestyle=:dot, linewidth=3)
	plot!(plt, xscale=:log)
	plot!(legend=false)
	#if n==2; plot!(plt, yscale=:log); end
	push!(pp, plt)
end
plot(pp[1], pp[2], pp[3], pp[4], pp[5], pp[6], pp[7], pp[8], size=(1000,500), layout=(2,4))
savefig("./BR_CTESN.svg")

max_stars_repo_name: cmarcotte/CTESN

max_issues_repo_name: cmarcotte/CTESN

max_forks_repo_name: cmarcotte/CTESN


text:
numbers = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1 ]

#https://docs.julialang.org/en/latest/manual/faq/?highlight=splat#...-combines-many-arguments-into-one-argument-in-function-definitions-1
function createphonenumber(numbers)
    return "($(numbers[1:3]...)) $(numbers[4:6]...)-$(numbers[7:10]...)"
end

res = createphonenumber(numbers)

println(res)

max_stars_repo_path: kata-julia/Create Phone Number.jl

max_stars_repo_name: crimsonskyrem/codewars

max_issues_repo_path: kata-julia/Create Phone Number.jl

max_issues_repo_name: crimsonskyrem/codewars

max_forks_repo_path: kata-julia/Create Phone Number.jl

max_forks_repo_name: crimsonskyrem/codewars


text:
using Documenter
using DecisionMakingUtils

DocMeta.setdocmeta!(DecisionMakingUtils, :DocTestSetup, :(using DecisionMakingUtils); recursive=true)
makedocs(
    sitename = "DecisionMakingUtils",
    format = Documenter.HTML(),
    modules = [DecisionMakingUtils]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/DecisionMakingAI/DecisionMakingUtils.jl.git",
)


text:
function create_Bins!(x_min::Float64, x_max::Float64, n_pop::Int64, Bins::Array{Tuple{Float64,Float64}, 1}, type::String)
    #### Fill the arrays Bins with tuples of 2 Floats. Each tuple contains the lower and upper limit of the bin ####
    #### type can be either "log" of "lin", to create logarithmic (m, a) or linear (e) bins                     ####
    for k = 1:n_pop
        if (type == "log")
            x_down = 10^(log10(x_min) + (log10(x_max) - log10(x_min)) * (k - 1) / n_pop)
            x_up = 10^(log10(x_min) + (log10(x_max) - log10(x_min)) * k / n_pop)
        elseif (type == "lin")
            x_down = x_min + (x_max - x_min) * (k - 1) / n_pop
            x_up = x_min + (x_max - x_min) * k / n_pop
        else
            return("type sould be either lin or log")
        end

        Bins[k] = (x_down, x_up)
    end
end

function create_KLN(m_min::Float64, m_max::Float64, γₘ::Float64, n_pop_m::Int64,
                        a_min::Float64, a_max::Float64, γₐ::Float64, n_pop_a::Int64,
                        e_min::Float64, e_max::Float64, γₑ::Float64, n_pop_e::Int64;
                        G::Float64 = 1.0, M::Float64 = 1.0)
    #### Returns a tuple of 6 Arrays. They are boundaries of the bins of stellar population used by the distribution, ####
    #### and the K, L, N Arrays corresponding to these bins                                                           ####
    # getting the constants we will need
    n_pop = n_pop_m * n_pop_a * n_pop_e # Total number of populations

    # Preallocating the arrays we will need
    Bins_m = Array{Tuple{Float64,Float64}, 1}(undef, n_pop_m) # Will hold the bottom and top of the bins in mass
    Bins_a = Array{Tuple{Float64,Float64}, 1}(undef, n_pop_a) # Will hold the bottom and top of the bins in semi-major axis
    Bins_e = Array{Tuple{Float64,Float64}, 1}(undef, n_pop_e) # Will hold the bottom and top of the bins in eccentricity
    K = Array{Tuple{Float64,Float64,Float64}, 1}(undef, n_pop) # Will hold the m, a, e of each stellar populations
    L = zeros(n_pop) # Will hold the angular momentum of each stellar populations
    N = zeros(n_pop) # Will hold the number of particles in each stellar population

    # Let's compute the boundaries of the bins
    create_Bins!(m_min, m_max, n_pop_m, Bins_m, "log")
    create_Bins!(a_min, a_max, n_pop_a, Bins_a, "log")
    create_Bins!(e_min, e_max, n_pop_e, Bins_e, "lin")

    # Let's fill the bins with stars
    for k_m = 1:n_pop_m
        m_down, m_up = Bins_m[k_m]
        mₖ = 10^((log10(m_down) + log10(m_up)) / 2)
        vₘₖ = m_up - m_down # Volume of the bin in mass
        if (vₘₖ == 0.0)
            vₘₖ = 1.0
        end

        for k_a = 1:n_pop_a
            a_down, a_up = Bins_a[k_a]
            aₖ = 10^((log10(a_down) + log10(a_up)) / 2)
            vₐₖ = a_up - a_down # Volume of the bin in sma
            if (vₐₖ == 0.0)
                vₐₖ = 1.0
            end

            for k_e = 1:n_pop_e
                e_down, e_up = Bins_e[k_e]
                eₖ = (e_down + e_up) / 2
                vₑₖ = e_up - e_down # volume of the bin in eccentricity
                if (vₑₖ == 0.0)
                    vₑₖ = 1.0
                end

                index = k_e + n_pop_e * (k_a - 1) + n_pop_e * n_pop_a * (k_m - 1)

                K[index] = (mₖ, aₖ, eₖ)
                L[index] = mₖ * sqrt(G * M * aₖ * (1 - eₖ^2))
                N[index] = vₘₖ * vₐₖ * vₑₖ * mₖ^(γₘ) * aₖ^(γₐ) * eₖ^(γₑ) # population = volume * density
            end
        end
    end
    N /= sum(N) # We normalize all distributions to have 1 star in total

    # Return
    return((Bins_m, Bins_a, Bins_e, K, L, N))
end

function get_EL_for_dist(Etot_norm::Float64, s::Float64,
                        L::Array{Float64, 1}, N::Array{Float64, 1}, m_min::Float64, a_min::Float64;
                        G::Float64 = 1.0)
    #### Return a tuple of 2 Floats. They are the the unnormalized total energy and momentum that can be used by a distribution ####
    #### NOTE : the cluster provided normalized energy and momentum, but the distribution works with unnormalized ones          ####
    # getting the constants we will need
    n_pop = length(L)

    # E
    J_typ = G * m_min^2 / a_min # typical value of J (circular orbits and twice the same population) -- To be improved
    Etot = Etot_norm * J_typ # Etot_norm = Etot / (J N^2) but a distribution only has 1 star

    # L
    Ltot_max_dist = 0.0
    for k = 1:n_pop
        Ltot_max_dist += L[k] * N[k]
    end
    Ltot = s * Ltot_max_dist # Ltot_norm = s = Ltot / Ltot_max

    # Return
    return((Etot, Ltot))
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
using Bravais
using Test

import StaticArrays
using LinearAlgebra
using FastClosures

debug = false

let lattice = BravaisLattice{3}([4,6,8])
    @test length(lattice) == 4*6*8
    @test size(lattice) == (4*6*8,)

    # NOTE: assumes hypercubic by default, but BravaisLattice does not
    # define any concept of nearest neighbors
    @test realspace(lattice, [1,4,3]) == [1,4,3]
end

let lattice = TriangularLattice([4,6])
    @test length(lattice) == 4*6
    @test realspace(lattice, [1,0]) == [1.0, 0]
    @test realspace(lattice, [0,1]) == [0.5, sqrt(3)/2]
    @test realspace(lattice, 7) == [1.0, 0]
end

@test length(BravaisLattice(StaticArrays.@SVector[4,6])) == 4*6
@test length(HypercubicLattice{4}([4,6,8,10])) == 4*6*8*10
@test length(TriangularLattice([4,6])) == 4*6
@test length(HoneycombLattice([3,5])) == 3*5*2
@test length(KagomeLattice([8,6])) == 8*6*3

@test_throws BoundsError checkbounds(SquareLattice([4,6]), 0)
@test_throws BoundsError checkbounds(SquareLattice([4,6]), 4*6+1)

lattices = Any[
    "Basic BravaisLattice",
    @inferred(BravaisLattice{2}([4,6])),
    @inferred(BravaisLattice(StaticArrays.@SVector[4,6])),

    "1D and effectively 1D",
    @inferred(ChainLattice([8])),
    @inferred(CubicLattice([1, 1, 8])),
    @inferred(CubicLattice([1, 8, 1])),

    "1D Twisted",
    @inferred(ChainLattice([6], Diagonal([6]), [0//1])),
    @inferred(ChainLattice([6], Diagonal([6]), [1//2])),
    @inferred(ChainLattice([6], Diagonal([6]), [1//5])),

    "2D",
    @inferred(SquareLattice([4,6])),
    @inferred(SquareLattice([6,4])),

    "2D twisted",
    @inferred(SquareLattice([4,6], Diagonal([4,6]), [1//7, 3//13])),
    @inferred(SquareLattice([6,4], Diagonal([6,4]), [1//7, 3//13])),

    "2D twisted cylinder",
    @inferred(SquareLattice([6,4], Diagonal([6,0]), [1//7, 0//1])),

    "2D twisted helix",
    @inferred(SquareLattice([6,4], [6 0; 1 4], [1//7, 3//13])),

    "3D like",
    @inferred(CubicLattice([4,2,6])),
    @inferred(HypercubicLattice{4}([4,2,6,1])),
    @inferred(HypercubicLattice(StaticArrays.@MVector[4,2,6,1])),

    "3D open twisted helix",
    @inferred(CubicLattice([4,6,8], [4 0 0; 0 6 0; 1 0 8], [1//7, 0//1, 3//13])),

    "Other lattices",
    @inferred(TriangularLattice([4,6])),
    @inferred(HoneycombLattice([3,5])),
    @inferred(KagomeLattice([8,6])),
]

for lattice in lattices
    if isa(lattice, AbstractString)
        debug && println(lattice)
        continue
    end

    len = @inferred length(lattice)
    @test size(lattice) == (len,)
    @test eachindex(lattice) == 1:len

    @test_throws BoundsError lattice[0]
    @test_throws BoundsError lattice[len+1]

    @static if VERSION >= v"1.2.0-DEV"
        # Relies on https://github.com/JuliaLang/julia/pull/27516
        @test @inferred(Nothing, findfirst(isequal(maxcoords(lattice)), lattice)) === nothing
    else
        @test findfirst(isequal(maxcoords(lattice)), lattice) === nothing
    end
    @inferred(in(@inferred(getindex(lattice, 1)),lattice))

    last = false
    for (i, site) in enumerate(lattice)
        checkbounds(lattice, i)
        @test lattice[i] == site
        @test site in lattice
        @test findfirst(isequal(site), lattice) == i
        @test last == false
        if i == length(lattice)
            last = true
        end

        @inferred wraparound(lattice, site) # FIXME: actually test something here

        @inferred realspace(lattice, site) # FIXME: actually test something here
        realspace(lattice, site, ones(Int, ndimensions(lattice))) # FIXME: actually test something here
    end
    @test last == true

    @inferred dimensions(lattice)
    d = @inferred ndimensions(lattice)
    η = @inferred twist(lattice)
    M = @inferred repeater(lattice)

    @test ishelical(lattice) == !isdiag(repeater(lattice))

    @test isbravais(lattice) == (lattice === bravais(lattice))

    # Check that translating multiple times, in total by an $\mathbf(A)_i$
    # vector, brings us back to where we started (and picking up a phase due
    # to $η_i$).
    sz = length(lattice)
    for i in 1:d
        if M[i,i] == 0
            @test_throws ArgumentError LatticeTranslationCache(lattice, i)
            continue
        end

        # Initialize our array of each index, along with a beginning
        # phase of zero
        indices = Tuple{Int, Rational{Int}}[(z, 0//1) for z in 1:sz]

        # Translate as many times as we need to in each dimension
        for j in 1:d
            if M[i,j] != 0
                ltrc = LatticeTranslationCache(lattice, j)
                for z in 1:M[i,j]
                    new_indices = Tuple{Int, Rational{Int}}[]
                    sizehint!(new_indices, sz)
                    for (s, old_η) in indices
                        newidx, wrap_η = @inferred translateη(lattice, s, j)
                        @test (newidx, wrap_η) == @inferred translateη(ltrc, s)
                        push!(new_indices, (newidx, wrap_η + old_η))
                    end
                    indices = new_indices
                end
            end
        end

        # assert that everything is as expected (i.e., overall we got the
        # identity element while picking up some phase)
        for z in 1:sz
            @test indices[z] == (z, η[i])
        end
    end

    if !isbravais(lattice)
        continue
    end

    @test transpose(primvecs(lattice)) * recivecs(lattice) ≈ Matrix(2π * I, ndimensions(lattice), ndimensions(lattice))

    # Check the momenta across the boundary conditions
    n_k_idx = @inferred nmomenta(lattice)
    @test 1:n_k_idx == @inferred eachmomentumindex(lattice)
    for site in lattice
        for k_idx in 1:n_k_idx
            for i in 1:d
                if M[i,i] != 0
                    site2 = site + vec(M[i,:])
                    exp1 = exp(im * (kdotr(lattice, k_idx, site) + 2π * η[i]))
                    exp2 = exp(im * @inferred(kdotr(lattice, k_idx, site2)))
                    idx, η_wrap = @inferred wraparoundη(lattice, site2)
                    exp3 = exp(im * (kdotr(lattice, k_idx, idx) + 2π * η_wrap))
                    @test exp1 ≈ exp2
                    @test exp1 ≈ exp3
                    for charge in (1, 3)
                        k_total = @inferred momentum(lattice, k_idx, charge)
                        exp1 = exp(im * (kdotr(lattice, k_total, site) + 2π * η[i] * charge))
                        exp2 = exp(im * kdotr(lattice, k_total, site2))
                        exp3 = exp(im * (kdotr(lattice, k_total, idx) + 2π * η_wrap * charge))
                        @test exp1 ≈ exp2
                        @test exp1 ≈ exp3
                    end
                    k = @inferred momentumspace(lattice, k_idx)
                    exp1 = exp(im * (dot(k, realspace(lattice, site)) + 2π * η[i]))
                    exp2 = exp(im * dot(k, realspace(lattice, site2)))
                    site3, wrap = @inferred wraparound_site(lattice, site2)
                    exp3 = exp(im * (dot(k, realspace(lattice, site3)) + 2π * dot(wrap, twist(lattice))))
                    @test exp1 ≈ exp2
                    @test exp1 ≈ exp3
                end
            end
        end
    end
end

let lattice = KagomeLattice([2,3])
    # try a valid wraparound
    @inferred wraparound(lattice, [4,0,0])
    # invalid basis index
    @test_throws ArgumentError wraparound(lattice, [4,0,3])
    # attempt to "translate" in direction of the basis index
    @test_throws ArgumentError translateη(lattice, 1, 3)
end

# Invalid wraparound for OBC
let lattice = ChainLattice([8], Diagonal([0]))
    @inferred wraparound(lattice, [0])
    @inferred wraparound(lattice, [7])
    @test_throws ArgumentError wraparound(lattice, [19])
    @test_throws ArgumentError wraparound(lattice, [-1])
end

function test_neighbors(lattice, pairs, neigh=Val{1})
    mypairs = Set(pairs)
    @closure neighborsη(lattice, neigh) do i, j, η
        pop!(mypairs, (i, j))
        # FIXME: check that η is correct
    end
    @test isempty(mypairs)
end

pairs = [
    (1, 2),
    (2, 3),
    (3, 4),
    (4, 5),
    (5, 6),
    (6, 1),
]
test_neighbors(ChainLattice([6], Diagonal([6]), [1//5]), pairs)
test_neighbors(SquareLattice([6, 1]), pairs)
test_neighbors(SquareLattice([1, 6]), pairs)
test_neighbors(CubicLattice([6, 1, 1]), pairs)
test_neighbors(HypercubicLattice{4}([1, 1, 6, 1]), pairs)

pairs = [
    (1, 3),
    (2, 4),
    (3, 5),
    (4, 6),
    (5, 7),
    (6, 8),
    (7, 1),
    (8, 2),
    (1, 2),
    (3, 4),
    (5, 6),
    (7, 8),
    ]
# XXX:
#test_neighbors(SquareLattice([4, 2]), pairs)
#test_neighbors(HypercubicLattice{4}([4, 1, 2, 1]), pairs)
# make sure open boundary conditions in the 2 direction is equivalent
# (unless i decide to remove this feature that doesn't double count these
# links)
test_neighbors(SquareLattice([4, 2], Diagonal([4, 0])), pairs)

# 4x3 Open
pairs = [
    (1, 4),
    (2, 5),
    (3, 6),
    (4, 7),
    (5, 8),
    (6, 9),
    (7, 10),
    (8, 11),
    (9, 12),
    (10, 1),
    (11, 2),
    (12, 3),
    (1, 2),
    (2, 3),
    (4, 5),
    (5, 6),
    (7, 8),
    (8, 9),
    (10, 11),
    (11, 12),
]
test_neighbors(SquareLattice([4, 3], Diagonal([4, 0])), pairs)

# 4x4
pairs = [
    (1, 2),
    (2, 3),
    (3, 4),
    (4, 1),
    (5, 6),
    (6, 7),
    (7, 8),
    (8, 5),
    (9, 10),
    (10, 11),
    (11, 12),
    (12, 9),
    (13, 14),
    (14, 15),
    (15, 16),
    (16, 13),
    (1, 5),
    (2, 6),
    (3, 7),
    (4, 8),
    (5, 9),
    (6, 10),
    (7, 11),
    (8, 12),
    (9, 13),
    (10, 14),
    (11, 15),
    (12, 16),
    (13, 1),
    (14, 2),
    (15, 3),
    (16, 4),
]
test_neighbors(SquareLattice([4, 4]), pairs)
test_neighbors(HypercubicLattice{4}([1, 4, 1, 4]), pairs)
test_neighbors(HypercubicLattice{4}([4, 1, 1, 4]), pairs)

#=
TEST(HypercubicLattice, NearestNeighbors_2x4) (
    std::set<std::pair<size_t, size_t> > pairs = (
        (0, 1),
        (1, 2),
        (2, 3),
        (3, 0),
        (4, 5),
        (5, 6),
        (6, 7),
        (7, 4),
        (0, 4),
        (1, 5),
        (2, 6),
        (3, 7)
    )
    test_neighbors(HypercubicLattice((2, 4)), pairs)
    test_neighbors(HypercubicLattice((1, 2, 4, 1)), pairs)
)

template <class LatticeType>
static void test_next_nearest_neighbors(const LatticeType &lattice, const std::set<std::pair<size_t, size_t> > &pairs)
(
    auto mypairs = pairs
    lattice.next_nearest_neighbors([&mypairs](size_t i, size_t j, boost::rational<int> eta) (
        const size_t erased = mypairs.erase((i, j))
        if (!erased)
            std::cerr << i << ' ' << j << '\n'
        ASSERT_EQ(erased, 1)

        // fixme: check that eta is correct
        (void) eta
#if 0
        std::cerr << eta << std::endl
#endif
    ))
    ASSERT_EQ(mypairs.size(), 0)
)

TEST(HypercubicLattice, NextNearestNeighbors_OBC) (
    std::set<std::pair<size_t, size_t> > pairs = (
        (0, 2),
        (1, 3),
        (2, 4),
        (3, 5)
    )
    Eigen::Matrix<int, 1, 1> repeater
    repeater(0, 0) = 0
    test_next_nearest_neighbors(HypercubicLattice((6), (0), repeater), pairs)
)

TEST(HypercubicLattice, NextNearestNeighbors_PBC) (
    std::set<std::pair<size_t, size_t> > pairs = (
        (0, 2),
        (1, 3),
        (2, 4),
        (3, 5),
        (4, 0),
        (5, 1)
    )
    test_next_nearest_neighbors(HypercubicLattice((6)), pairs)
)
=#

function test_neighbor_sublattices(lattice, allowed, neigh=Val{1})
    @inferred isbipartite(lattice)
    @inferred istripartite(lattice)

    @closure neighbors(lattice, neigh) do i, j, wrap
        ind1 = @inferred sublattice_index(lattice, i)
        ind2 = sublattice_index(lattice, j)
        @test ind1 in allowed
        @test ind2 in allowed
        @test ind1 != ind2
    end
end

@test isbipartite(SquareLattice([4,4]))
@test !isbipartite(SquareLattice([3,3]))
@test isbipartite(SquareLattice([3,3], Diagonal([0,0])))
@test !isbipartite(SquareLattice([4,3]))
@test !isbipartite(SquareLattice([4,1])) # if two dimensions are specified, we don't consider this bipartite.
@test !istripartite(SquareLattice([4,4]))

let lattice = SquareLattice([4, 3], Diagonal([4,0]))
    test_neighbor_sublattices(lattice, [0,1])
    @test sublattice_index(lattice, 1) == 0
    @test sublattice_index(lattice, 2) == 1
    @test sublattice_index(lattice, 3) == 0
    @test sublattice_index(lattice, 4) == 1
    @test sublattice_index(lattice, 5) == 0
    @test sublattice_index(lattice, 6) == 1
    @test sublattice_index(lattice, 7) == 0
    @test sublattice_index(lattice, 8) == 1
    @test sublattice_index(lattice, 9) == 0
    @test sublattice_index(lattice, 10) == 1
    @test sublattice_index(lattice, 11) == 0
    @test sublattice_index(lattice, 12) == 1
end

let lattice = SquareLattice([4, 3])
    @test_throws ArgumentError sublattice_index(lattice, 1)
end

let lattice = SquareLattice([4, 2])
    test_neighbor_sublattices(lattice, [0,1])
    @test sublattice_index(lattice, 1) == 0
    @test sublattice_index(lattice, 2) == 1
    @test sublattice_index(lattice, 3) == 1
    @test sublattice_index(lattice, 4) == 0
    @test sublattice_index(lattice, 5) == 0
    @test sublattice_index(lattice, 6) == 1
    @test sublattice_index(lattice, 7) == 1
    @test sublattice_index(lattice, 8) == 0
end

let lattice = ChainLattice([6])
    test_neighbor_sublattices(lattice, [0,1])
    @test sublattice_index(lattice, 1) == 0
    @test sublattice_index(lattice, 2) == 1
    @test sublattice_index(lattice, 3) == 0
    @test sublattice_index(lattice, 4) == 1
    @test sublattice_index(lattice, 5) == 0
    @test sublattice_index(lattice, 6) == 1
end

let lattice = SquareLattice([6,1])
    @test_throws ArgumentError sublattice_index(lattice, 1)
end

let lattice = TriangularLattice([3,3])
    @test !isbipartite(lattice)
    @test istripartite(lattice)
end
let lattice = TriangularLattice([4,3])
    @test !isbipartite(lattice)
    @test !istripartite(lattice)
    @test_throws ArgumentError sublattice_index(lattice, 1)
end
let lattice = TriangularLattice([4,3], Diagonal([0,3]))
    @test !isbipartite(lattice)
    @test istripartite(lattice)
end
let lattice = TriangularLattice([4,3], Diagonal([0,0]))
    @test !isbipartite(lattice)
    @test istripartite(lattice)
end

let lattice = TriangularLattice([3,3])
    test_neighbor_sublattices(lattice, [0,1,2])
    @test sublattice_index(lattice, 1) == 0
    @test sublattice_index(lattice, 2) == 1
    @test sublattice_index(lattice, 3) == 2
    @test sublattice_index(lattice, 4) == 2
    @test sublattice_index(lattice, 5) == 0
    @test sublattice_index(lattice, 6) == 1
    @test sublattice_index(lattice, 7) == 1
    @test sublattice_index(lattice, 8) == 2
    @test sublattice_index(lattice, 9) == 0
end

let lattice = HoneycombLattice([5,5])
    @test isbipartite(lattice)
    @test !istripartite(lattice)
    test_neighbor_sublattices(lattice, [0,1])
    for (i, site) in enumerate(lattice)
        @test sublattice_index(lattice, i) == sublattice_index(lattice, site) == site[end] == ((i ⊻ 1) & 1)
    end
end

let lattice = KagomeLattice([4,3])
    @test !isbipartite(lattice)
    @test istripartite(lattice)
    test_neighbor_sublattices(lattice, [0,1,2])
    for (i, site) in enumerate(lattice)
        @test sublattice_index(lattice, i) == sublattice_index(lattice, site) == site[end]
    end
end

function test_neighbor_distances(lattice, neigh=Val{1}, expected_squared=1)
    @closure neighbors(lattice, neigh) do i, j, wrap
        diff = realspace(lattice, i) - realspace(lattice, j) - realspace(bravais(lattice), wrap .* dimensions(lattice))
        dist_squared = dot(diff, diff)
        @test expected_squared ≈ dist_squared atol=1e-8
    end
end

test_neighbor_distances(ChainLattice([8]))
test_neighbor_distances(ChainLattice([8], Diagonal([0])))
test_neighbor_distances(ChainLattice([8]), Val{1}, 1)
test_neighbor_distances(ChainLattice([8], Diagonal([0])), Val{1}, 1)
test_neighbor_distances(ChainLattice([8]), Val{2}, 4)
test_neighbor_distances(ChainLattice([8], Diagonal([0])), Val{2}, 4)
test_neighbor_distances(ChainLattice([8]), Val{5}, 25)
test_neighbor_distances(ChainLattice([8], Diagonal([0])), Val{5}, 25)
test_neighbor_distances(SquareLattice([4,6]))
test_neighbor_distances(SquareLattice([4,6], Diagonal([0,0])))
test_neighbor_distances(SquareLattice([4,6]), Val{2}, 2)
test_neighbor_distances(SquareLattice([4,6], Diagonal([0,0])), Val{2}, 2)
test_neighbor_distances(TriangularLattice([4,6]))
test_neighbor_distances(TriangularLattice([4,6], Diagonal([0,0])))
test_neighbor_distances(TriangularLattice([4,6]), Val{2}, 3)
test_neighbor_distances(TriangularLattice([4,6], Diagonal([0,0])), Val{2}, 3)
test_neighbor_distances(TriangularLattice([4,6]), Val{3}, 4)
test_neighbor_distances(TriangularLattice([4,6], Diagonal([0,0])), Val{3}, 4)
test_neighbor_distances(HoneycombLattice([4,6]))
test_neighbor_distances(HoneycombLattice([4,6], Diagonal([0,0])))
test_neighbor_distances(KagomeLattice([4,6]))
test_neighbor_distances(KagomeLattice([4,6], Diagonal([0,0])))


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Imports
import Base: string, print, show

# Export types
export Raster, PixelMap

# Export methods
export nlen, glen, blen, slen, map, bsize, bdim

"""
    Raster

A rasterisation of the bounding box of a geometry defined by the mesh, with
methods capable of transforming functions between nodal, and rasterised bases.

See PixelMap.
"""
@compat abstract type Raster end

@enum RasterBases Pixel CubicPixel GaussBlob BesselBlob HanningBlob RampBlob SplineBlob

"""
    PixelMap(mesh, dims; gscale=2)

Map from nodal coefficients to a rasterised grid of pixels of specified
dimension using linear interpolation via an intermediate grid containing an
integer multiple, gscale, of pixels.

# Arguments
* `mesh::Mesh`: a TOAST mesh
* `dims::NTuple{ndim,Integer}`: the dimensions of the raster basis
* `gscale=2::Int`: the scale factor used for the intermediate mapping raster
"""
mutable struct PixelMap <: Raster

  ptr::Cxx.CxxCore.CppPtr
  mesh::Mesh

  function PixelMap(mesh::Mesh, bdim::NTuple{N,Integer}; gscale::Integer=2) where {N}
    rasterptr = _rasterast_new(mesh, Pixel, bdim, gscale)
    pixelmap = new(rasterptr, mesh)
    finalizer(_raster_delete, pixelmap)
    return pixelmap
  end

end

# Create and initialise a new basis mapper
function _rasterast_new(mesh::Mesh,
                        basis::RasterBases,
                        dims::NTuple{N,Integer},
                        gscale::Integer;
                        blobarg::Float64 = 1.0,
                        blobrad::Float64 = 1.0) where {N}

  bdimarr = [Cint(i) for i in dims]
  gdimarr = [Cint(i*gscale) for i in dims]
  basiscode = Cint(basis)
  dim = Cint(dimensions(mesh))

  @assert length(bdimarr)==dim

  meshptr = mesh.ptr

  bd = @cxxnew IVector(Cint(2), pointer(bdimarr), icxx"""return DEEP_COPY;""")
  gd = @cxxnew IVector(Cint(2), pointer(gdimarr), icxx"""return DEEP_COPY;""")

  ptr = icxx"""

    RDenseMatrix *bb = 0;
    Raster *raster = 0;

    switch($basiscode) {
      case 0:
        raster = new Raster_Pixel (*$(bd), *$(gd), $meshptr, bb);
        break;
      case 1:
        raster = new Raster_CubicPixel (*$(bd), *$(gd), $meshptr, bb);
        break;
      case 2:
        raster = new Raster_GaussBlob (*$(bd),*$(gd), $meshptr, $blobarg, $blobrad, bb);
        break;
      case 3:
        raster = new Raster_BesselBlob (*$(bd), *$(gd), $meshptr, $blobarg, $blobrad, bb);
        break;
      case 4:
        raster = new Raster_HanningBlob (*$(bd), *$(gd), $meshptr, $blobrad, bb);
        break;
      case 5:
        raster = new Raster_RampBlob (*$(bd), *$(gd), $meshptr, $blobrad, bb);
        break;
      case 6:
        raster = new Raster_SplineBlob (*$(bd), *$(gd), $meshptr, $blobrad, bb);
        break;
    }

    return raster;
  """

  @assert ptr != C_NULL

  return ptr;

end

# Delete a raster map
_raster_delete(Raster::Raster) = finalize(Raster.ptr)

"""
    nlen(raster)

Return the number of nodal coefficients in the mesh basis associated with the
raster.
"""
nlen(Raster::Raster) = nodecount(Raster.mesh)

"""
    slen(raster)

Return the number of coefficients in the raster solution basis, which does
not include raster points outside of the support of the underlying mesh.
"""
slen(Raster::Raster) = @cxx Raster.ptr->SLen()

"""
    blen(raster)

Return the number of coefficients in the raster.
"""
blen(Raster::Raster) = @cxx Raster.ptr->BLen()

"""
    glen(raster)

Return the number of coefficients in the raster intermediate basis, which is
used to map between bases.
"""
glen(Raster::Raster) = @cxx Raster.ptr->GLen()

"""
    bdim(raster)

Return the number of elements in each dimension of the raster basis.
"""
function bdim(raster::Raster)

  ndim = dimensions(raster.mesh)
  dims = Vector{Int}(undef, ndim)
  dimp = pointer(dims)

  icxx"""
    for(int i=0; i<$(ndim); i++)
      $(dimp)[i] = $(raster.ptr)->BDim()[i];
  """

  return dims

end

"""
    bsize(raster)

The `size` (area/volume) of each element of the raster basis.
"""
function bsize(raster::Raster)
  bb = boundingbox(raster.mesh)
  prod((bb[2,:] .- bb[1,:])./ bdim(raster) )
end


text:
@testset "Simulation Sequence Correct Execution Order" begin
    problems = SimulationProblems(
        DAUC = OperationsProblem(MockOperationProblem; horizon = 48),
        HAUC = OperationsProblem(MockOperationProblem; horizon = 24),
        ED = OperationsProblem(MockOperationProblem; horizon = 12),
        AGC = OperationsProblem(MockOperationProblem; horizon = 6),
    )

    # Test RH sequences
    feedforward_chronologies = Dict(
        ("UC" => "HAUC") => Synchronize(periods = 24),
        ("HAUC" => "ED") => RecedingHorizon(),
        ("ED" => "AGC") => RecedingHorizon(),
    )
    ini_cond_chronology = InterProblemChronology()
    intervals = Dict(
        "DAUC" => (Hour(24), Consecutive()),
        "HAUC" => (Hour(1), RecedingHorizon()),
        "ED" => (Minute(5), RecedingHorizon()),
        "AGC" => (Minute(1), Consecutive()),
    )
    test_sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = feedforward_chronologies,
        intervals = intervals,
        ini_cond_chronology = ini_cond_chronology,
    )

    @test length(findall(x -> x == 4, test_sequence.execution_order)) == 24 * 60
    @test length(findall(x -> x == 3, test_sequence.execution_order)) == 24 * 12
    @test length(findall(x -> x == 2, test_sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, test_sequence.execution_order)) == 1

    for name in PSI.get_problem_names(problems)
        @test problems[name].internal.simulation_info.sequence_uuid == test_sequence.uuid
    end

    # Test single stage sequence
    test_sequence = SimulationSequence(
        problems = SimulationProblems(
            DAUC = OperationsProblem(MockOperationProblem; horizon = 48),
        ),
        intervals = Dict("DAUC" => (Hour(24), Consecutive())),
        ini_cond_chronology = InterProblemChronology(),
    )

    @test isa(test_sequence.ini_cond_chronology, IntraProblemChronology)
    @test test_sequence.execution_order == [1]

    # Test synchronized sequence
    problems = SimulationProblems(
        MD = OperationsProblem(MockOperationProblem; horizon = 3),
        DAUC = OperationsProblem(MockOperationProblem; horizon = 48),
        HAUC = OperationsProblem(MockOperationProblem; horizon = 24),
        ED = OperationsProblem(MockOperationProblem; horizon = 12),
        AGC = OperationsProblem(MockOperationProblem; horizon = 6),
    )

    feedforward_chronologies = Dict(
        ("MD" => "UC") => Synchronize(periods = 2),
        ("UC" => "HAUC") => Synchronize(periods = 24),
        ("HAUC" => "ED") => Synchronize(periods = 1),
        ("ED" => "AGC") => Synchronize(periods = 1),
    )
    ini_cond_chronology = InterProblemChronology()
    intervals = Dict(
        "MD" => (Hour(48), Consecutive()),
        "DAUC" => (Hour(24), Consecutive()),
        "HAUC" => (Hour(4), Consecutive()),
        "ED" => (Minute(5), Consecutive()),
        "AGC" => (Minute(1), Consecutive()),
    )
    test_sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = feedforward_chronologies,
        intervals = intervals,
        ini_cond_chronology = ini_cond_chronology,
    )

    @test length(findall(x -> x == 5, test_sequence.execution_order)) == 24 * 2 * 60
    @test length(findall(x -> x == 4, test_sequence.execution_order)) == 24 * 2 * 12
    @test length(findall(x -> x == 3, test_sequence.execution_order)) == 24 * 2 / 4
    @test length(findall(x -> x == 2, test_sequence.execution_order)) == 2
    @test length(findall(x -> x == 1, test_sequence.execution_order)) == 1

    for name in PSI.get_problem_names(problems)
        @test problems[name].internal.simulation_info.sequence_uuid == test_sequence.uuid
    end
end

@testset "Simulation Sequence invalid sequences" begin
    @test_throws ArgumentError SimulationSequence(
        problems = mock_uc_ed_simulation_problems(48, 12),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
        ini_cond_chronology = InterProblemChronology(),
    )

    @test_throws IS.ConflictingInputsError SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(2), Consecutive()),
            "ED" => (Hour(3), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
        ini_cond_chronology = InterProblemChronology(),
    )
end

@testset "Test print methods of sequence ascii art" begin
    sequence_2 = SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 2)),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    sequence_4 = SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 4)),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    sequence_3 = SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 3)),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    sequence_5 = SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 2)),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => RangeFF(
                variable_source_problem_ub = PSI.ON,
                variable_source_problem_lb = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    sequence_13 = SimulationSequence(
        problems = mock_uc_ed_simulation_problems(24, 12),
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon(periods = 13)),
        intervals = Dict(
            "UC" => (Hour(1), RecedingHorizon()),
            "ED" => (Minute(5), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    list = [sequence_2, sequence_3, sequence_4, sequence_5, sequence_13]
    _test_plain_print_methods(list)
    stage_1 = FakeStagesStruct(Dict(1 => 1, 2 => 2, 3 => 3, 4 => 4, 5 => 5, 6 => 6)) # testing 5 stages
    stage_3 = FakeStagesStruct(Dict(1 => 1, 2 => 100)) #testing 3 digits
    stage_4 = FakeStagesStruct(Dict(1 => 1, 2 => 1000)) #testing 4 digits
    stage_12 = FakeStagesStruct(Dict(1 => 1, 2 => 12, 3 => 5, 4 => 6))
    list = [stage_1, stage_3, stage_4, stage_12]
    _test_plain_print_methods(list)
end


text:
"""
$(TYPEDEF)

A `SinCosRevolute` joint type allows rotation about a fixed axis.

In contrast to the [`Revolute`](@ref) joint type, the configuration vector for the `SinCosRevolute` joint type
consists of the sine and cosine of the angle of rotation about the specified axis (in that order).
The velocity vector for the `SinCosRevolute` joint type is the same as for
the `Revolute` joint type, i.e., the time derivative of the angle about the axis.
"""
struct SinCosRevolute{T} <: JointType{T}
    axis::SVector{3, T}
    rotation_from_z_aligned::RotMatrix3{T}

    function SinCosRevolute{T}(axis::AbstractVector) where {T}
        a = normalize(axis)
        new{T}(a, rotation_between(SVector(zero(T), zero(T), one(T)), SVector{3, T}(a)))
    end
end

"""
$(SIGNATURES)

Construct a new `SinCosRevolute` joint type, allowing rotation about `axis`
(expressed in the frame before the joint).
"""
SinCosRevolute(axis::AbstractVector{T}) where {T} = SinCosRevolute{T}(axis)

Base.show(io::IO, jt::SinCosRevolute) = print(io, "SinCosRevolute joint with axis $(jt.axis)")

function Random.rand(::Type{SinCosRevolute{T}}) where {T}
    axis = normalize(randn(SVector{3, T}))
    SinCosRevolute(axis)
end

RigidBodyDynamics.flip_direction(jt::SinCosRevolute) = SinCosRevolute(-jt.axis)

num_positions(::Type{<:SinCosRevolute}) = 2
num_velocities(::Type{<:SinCosRevolute}) = 1
has_fixed_subspaces(jt::SinCosRevolute) = true
isfloating(::Type{<:SinCosRevolute}) = false

@propagate_inbounds function set_configuration!(q::AbstractVector, joint::Joint{<:Any, <:SinCosRevolute}, θ::Number)
    @boundscheck check_num_positions(joint, q)
    s, c = sincos(θ)
    @inbounds q[1] = s
    @inbounds q[2] = c
    q
end

@propagate_inbounds function set_velocity!(v::AbstractVector, joint::Joint{<:Any, <:SinCosRevolute}, θ̇::Number)
    @boundscheck check_num_velocities(joint, v)
    @inbounds v[1] = θ̇
    v
end

@propagate_inbounds function rand_configuration!(q::AbstractVector, jt::SinCosRevolute)
    q .= normalize(@SVector(randn(2)))
    nothing
end

@propagate_inbounds function zero_configuration!(q::AbstractVector, jt::SinCosRevolute)
    T = eltype(q)
    q[1] = zero(T)
    q[2] = one(T)
    nothing
end

@propagate_inbounds function joint_transform(jt::SinCosRevolute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D, q::AbstractVector)
    # from https://github.com/FugroRoames/Rotations.jl/blob/8d77152a76950665302b5a730b420973bb397f41/src/angleaxis_types.jl#L50
    T = eltype(q)

    axis = jt.axis
    s, c = q[1], q[2]
    c1 = one(T) - c

    c1x2 = c1 * axis[1]^2
    c1y2 = c1 * axis[2]^2
    c1z2 = c1 * axis[3]^2

    c1xy = c1 * axis[1] * axis[2]
    c1xz = c1 * axis[1] * axis[3]
    c1yz = c1 * axis[2] * axis[3]

    sx = s * axis[1]
    sy = s * axis[2]
    sz = s * axis[3]

    # Note that the RotMatrix constructor argument order makes this look transposed:
    rot = RotMatrix(
        one(T) - c1y2 - c1z2, c1xy + sz, c1xz - sy,
        c1xy - sz, one(T) - c1x2 - c1z2, c1yz + sx,
        c1xz + sy, c1yz - sx, one(T) - c1x2 - c1y2)

    Transform3D(frame_after, frame_before, rot)
end

@propagate_inbounds function joint_twist(jt::SinCosRevolute, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    angular = jt.axis * v[1]
    Twist(frame_after, frame_before, frame_after, angular, zero(angular))
end

@inline function bias_acceleration(jt::SinCosRevolute{T}, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector{X}, v::AbstractVector{X}) where {T, X}
    S = promote_type(T, X)
    zero(SpatialAcceleration{S}, frame_after, frame_before, frame_after)
end

@propagate_inbounds function joint_spatial_acceleration(jt::SinCosRevolute{T}, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector{X}, v::AbstractVector{X}, vd::AbstractVector{XD}) where {T, X, XD}
    S = promote_type(T, X, XD)
    angular = convert(SVector{3, S}, jt.axis * vd[1])
    SpatialAcceleration(frame_after, frame_before, frame_after, angular, zero(angular))
end

@inline function motion_subspace(jt::SinCosRevolute{T}, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector{X}) where {T, X}
    S = promote_type(T, X)
    angular = SMatrix{3, 1, S}(jt.axis)
    linear = zero(SMatrix{3, 1, S})
    GeometricJacobian(frame_after, frame_before, frame_after, angular, linear)
end

@inline function constraint_wrench_subspace(jt::SinCosRevolute{T}, joint_transform::Transform3D{X}) where {T, X}
    S = promote_type(T, X)
    R = convert(RotMatrix3{S}, jt.rotation_from_z_aligned)
    Rcols12 = R[:, SVector(1, 2)]
    angular = hcat(Rcols12, zero(SMatrix{3, 3, S}))
    linear = hcat(zero(SMatrix{3, 2, S}), R)
    WrenchMatrix(joint_transform.from, angular, linear)
end

@propagate_inbounds function joint_torque!(τ::AbstractVector, jt::SinCosRevolute, q::AbstractVector, joint_wrench::Wrench)
    τ[1] = dot(angular(joint_wrench), jt.axis)
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity!(v::AbstractVector, jt::SinCosRevolute, q::AbstractVector, q̇::AbstractVector)
    # q̇ = [c; -s] * v
    # [c -s] * [c; -s] = c^2 + s^2 = 1
    # v = [c -s] * q̇
    s, c = q[1], q[2]
    v[1] = c * q̇[1] - s * q̇[2]
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity_jacobian(jt::SinCosRevolute, q::AbstractVector)
    s, c = q[1], q[2]
    @SMatrix [c -s]
end

@propagate_inbounds function configuration_derivative_to_velocity_adjoint!(fq, jt::SinCosRevolute, q::AbstractVector, fv)
    qnorm = sqrt(q[1]^2 + q[2]^2)
    qnormalized = @SVector([q[1], q[2]]) / qnorm
    fq .= (configuration_derivative_to_velocity_jacobian(jt, qnormalized)' * fv) ./ qnorm
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative!(q̇::AbstractVector, jt::SinCosRevolute, q::AbstractVector, v::AbstractVector)
    s, c = q[1], q[2]
    q̇[1] = c * v[1]
    q̇[2] = -s * v[1]
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative_jacobian(jt::SinCosRevolute, q::AbstractVector)
    s, c = q[1], q[2]
    @SMatrix [c; -s]
end

# uses exponential coordinates centered around q0 (i.e, the relative joint angle)
@propagate_inbounds function local_coordinates!(ϕ::AbstractVector, ϕd::AbstractVector,
        jt::SinCosRevolute, q0::AbstractVector, q::AbstractVector, v::AbstractVector)
    # RΔ = R₀ᵀ * R
    s0, c0 = q0[1], q0[2]
    s, c = q[1], q[2]
    sΔ = c0 * s - s0 * c
    cΔ = c0 * c + s0 * s
    θ = atan(sΔ, cΔ)
    ϕ[1] = θ
    ϕd[1] = v[1]
    nothing
end

@propagate_inbounds function global_coordinates!(q::AbstractVector, jt::SinCosRevolute, q0::AbstractVector, ϕ::AbstractVector)
    # R = R₀ * RΔ:
    s0, c0 = q0[1], q0[2]
    θ = ϕ[1]
    sΔ, cΔ = sincos(θ)
    s = s0 * cΔ + c0 * sΔ
    c = c0 * cΔ - s0 * sΔ
    q[1] = s
    q[2] = c
    nothing
end

@propagate_inbounds function normalize_configuration!(q::AbstractVector, jt::SinCosRevolute)
    q ./= sqrt(q[1]^2 + q[2]^2)
    q
end

@propagate_inbounds function is_configuration_normalized(jt::SinCosRevolute, q::AbstractVector, rtol, atol)
    # TODO: base off of norm squared
    qnorm = sqrt(q[1]^2 + q[2]^2)
    isapprox(qnorm, one(eltype(q)); rtol = rtol, atol = atol)
end


text:
"""
    prediction( sys::StateSpace, N::Integer )

Form the prediction matrix ``Γ`` that computes the states of the linear system `sys` over the horizon of length `N`
by using all the inputs across the entire horizon.

# Extended help

This matrix consists of the blocks ``A^i B`` where ``i`` is the diagonal number, forming the matrix
```math
⎡ B   •   •  •  ⋯⎤
  ⎢ AB  B   •  •  ⋯⎥
  ⎢ A²B AB  B  •  ⋯⎥
  ⎢ A³B A²B AB B  ⋯⎥
  ⎣ ⋮   ⋮   ⋮  ⋮ ⋱ ⎦
```

It is used to compute the component of the state vector for a horizon that is caused by the system inputs
during the horizon using the equation ``x = Γu``.

The full state vector for a horizon can be computed as ``x = Γu + Φx₀``, where ``Φ`` is the initial condition
propagation matrix from [`initialpropagation`](@ref).
"""
function prediction( sys::StateSpace, N::Integer )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    nu = ninputs( sys )
    nx = nstates( sys )
    mt = eltype( sys.A )

    mainDiag  = kron( I(N), sys.B )

    otherDiag = kron( I(N), to_matrix( mt, I, nx ) )

    D = zeros( mt, N, N )
    D[diagind( D, -1 )] .= 1

    otherDiag = otherDiag + kron( D, -sys.A )

    Γ = otherDiag\mainDiag

    return BlockArray{ mt }( Γ, [nx for i = 1:N], [nu for i = 1:N] )
end

"""
    initialpropagation( sys::StateSpace, N::Integer )

Form the matrix ``Φ`` that propagates the initial condition ``x₀`` in the linear system `sys` across the horizon `N`.

# Extended help

This matrix consists of a column of the blocks ``A^i`` where ``i`` is the row number, forming the matrix
```math
⎡ A  ⎤
  ⎢ A² ⎥
  ⎢ A³ ⎥
  ⎢ A⁴ ⎥
  ⎣ ⋮  ⎦
```

It is used to compute the component of the state vector for a horizon that is caused by the intial state of the system using
the equation ``x = Φx₀``.

The full state vector for a horizon can be computed as ``x = Γu + Φx₀``, where ``Γ`` is the prediction matrix
from [`prediction`](@ref).
"""
function initialpropagation( sys::StateSpace, N::Integer )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    blockGen = (sys.A^i for i=1:N)

    # The mortar function requires the passed array to have ndims=2, so reshape to
    # have an array with 2 dimensions even though there is only one column
    blocks = reshape( collect( blockGen ),  (N, 1) )

    return mortar( blocks )
end

function inputprediction( sys::StateSpace, K::AbstractArray, N::Integer )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    nu = ninputs( sys )
    nx = nstates( sys )
    mt = eltype( sys.A )

    # Create the prediction matrix for the controlled system
    Γ = prediction( sys, N )
    K̄ = blockkron( 1.0*I(N), K )

    zerorow = zeros( Float64, (nu, nu*N) )

    Γ = -K̄*Γ
    Γ = vcat( zerorow, Γ )

    eye = blockkron( 1.0*I(N), Matrix( 1.0*I(nu) ) )
    eye = vcat( eye, zerorow )

    Γ = Γ + eye
    Γ = Γ[1:N*nu, 1:N*nu]

    return BlockArray{ mt }( Γ, [nu for i = 1:N], [nu for i = 1:N] )
end

function inputinitialpropagation( sys::StateSpace, K::AbstractArray, N::Integer )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    blockGen = (-K*sys.A^i for i=0:(N-1))

    # The mortar function requires the passed array to have ndims=2, so reshape to
    # have an array with 2 dimensions even though there is only one column
    blocks = reshape( collect( blockGen ),  (N, 1) )

    return mortar( blocks )
end

"""
    hessian( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )

Form the complete Hessian for the constrained time-invariant LQR problem `clqr` using the fully
condensed form.
"""
function hessian( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sys = getsystem( clqr )
    N   = clqr.N
    nx  = nstates( sys )
    nu  = ninputs( sys )

    Γ = prediction( sys, N )

    Q̅ = blockkron( 1.0*I(N), clqr.Qₖ )
    S̅ = blockkron( 1.0*I(N), clqr.S  )
    R̅ = blockkron( 1.0*I(N), clqr.R  )

    Q̅[Block( N, N )] = clqr.P
    S̅[Block( N, N )] = zeros( eltype( clqr.S ), size( clqr.S ) )

    # Pass the result through the Hermitian type to cleanup numerical errors that can occur
    # making it non-symmetric (they are on the order of 10e-17)
    return Hermitian( Γ'*Q̅*Γ + Γ'*S̅ + S̅'*Γ + R̅ )
end

"""
    linearcoefficients( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )

Form the vector of coefficients for the linear term of the fully condensed quadratic QP for the
CLQR problem defined by `clqr`.
"""
function linearcoefficients( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sys = getsystem( clqr )
    N   = clqr.N
    nx  = nstates( sys )
    nu  = ninputs( sys )

    Γ = prediction( sys, N )
    Φ = initialpropagation( sys, N )

    Q̅ = blockkron( 1.0*I(N), clqr.Qₖ )
    S̅ = blockkron( 1.0*I(N), clqr.S  )

    Q̅[Block( N, N )] = clqr.P
    S̅[Block( N, N )] = zeros( eltype( clqr.S ), size( clqr.S ) )

    return Γ'*Q̅*Φ + S̅'*Φ
end

function inequalityconstraints( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    # For the initial computations we want the non-pre-stabilized system
    sys = getsystem( clqr, prestabilized = false )

    # Find the number of constraints and the system size
    N   = clqr.N
    nC  = size( clqr.E, 1 );    # Number of constraints
    n, m  = size( sys.B );

    # Get the condensed initial propagation matrix
    Φ = initialpropagation( sys, N )

    # Get the prediction matrix for the system
    Γ = prediction( sys, N )

    # Create the component matrices
    Ē = blockkron( 1.0*I(N), clqr.E )
    F̄ = blockkron( 1.0*I(N), clqr.F )

    # Compute the coefficient matrix of the linear inequality constraints
    G = Ē*Γ + F̄;

    # Compute the initial state matrix for the RHS of the constraints
    L = -Ē*Φ;

    if isprestabilized( clqr )
        # Get the pre-stabilized matrices
        sysₖ = getsystem( clqr )
        Γₖ   = inputprediction( sysₖ, clqr.K, N )
        Φₖ   = inputinitialpropagation( sysₖ, clqr.K, N )

        # Modify the initial state coefficients to be in the new input space
        L = L - G*Φₖ

        # Modify the coefficient matrix to be in the new input space
        G = G * Γₖ
    end

    # Compute the constant vector for the RHS
    g = blockkron( ones( Float64, (N) ), clqr.g )

    return G, L, g
end

function equalityconstraints( ::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    # The fully condensed form does not have any equality constraints from the problem
    return nothing
end