text:
#========================================================================================#
#	Laboratory 5
#
# Exploratory data analysis, datasets and Dataframes.
#
# Author: Niall Palfreyman, 04/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we use Exploratory Data Analysis (EDA) and the Julia DataFrames
		package to explore one of the over 750 industry standard statistical RDatasets.

		Since we haven't used the RDatasets package before, we must first add it to our Julia
		system. To do this, type ']' at the Julia prompt - this starts Julia's Package Manager,
		which you can exit at any time by typing <Backspace>. Within the Package Manager, enter:

		add RDatasets

		This will download and compile around 200 packages that are relevant to the RDatasets
		package (this may take a couple of minutes). When it has finished, exit Package Manager,
		then enter:

		using RDatasets						# Takes a few seconds ...
		RDatasets.datasets()

		How many Rows are in the "affairs" dataset?
		""",
		"\"affairs\" is the first dataset in the list of datasets that datasets() prints out.",
		x -> x==601
	),
	Activity(
		"""
		OK, now we've got RDatasets up and running, we'll look specifically at the Iris dataset
		compiled by Edgar Anderson and Ronald Fisher in the early days of genetics research in
		1936. This dataset contains the length and width (in cm) of the sepals and petals of 50
		samples from each of several iris flower species. Enter this code to load the Iris dataset:

		iris = dataset("datasets","iris")

		What is the SepalLength of the first specimen in iris?
		""",
		"This is the specimen in Row 1",
		x -> x==5.1
	),
	Activity(
		"""
		The Iris dataset is an object iris of type DataFrame: a table of rows and columns. Each
		row represents a single specimen, and each column represents an attribute field of that
		specimen. To view the fields of the Iris dataset, enter names(iris). How many fields
		does iris have?
		""",
		"",
		x -> x==5
	),
	Activity(
		"""
		What is the size() of the dataset?
		""",
		"size(iris)",
		x -> x==size(Main.iris)
	),
	Activity(
		"""
		We can view the first five rows of iris using the usual array indexing with square
		brackets []:

		iris[1:5,:]

		Report back to me the last 5 rows of the Iris dataset:
		""",
		"iris[end-4:end,:]",
		x -> x==Main.iris[end-4:end,:]
	),
	Activity(
		"""
		We've seen that the first rows of iris concern the species setosa, and the last rows
		concern the species virginica. What other species does iris contain. To find out, we
		want to group the dataset by Species using the command:

		species_groups = groupby(iris,"Species")

		How many groups are contained in this grouping of the data?
		""",
		"size(species_groups)",
		x -> x==3
	),
	Activity(
		"""
		OK, so we've seen that the first of three Species groups contains "iris setosa" specimens,
		and the last Species group contains "iris virginica" specimens. To inspect all three
		Species groups, we will combine specimens from each group into a single row:

		combine(species_groups,nrow)

		What is the name of the second iris species?
		""",
		"",
		x -> occursin(lowercase(x),"versicolor")
	),
	Activity(
		"""
		Don't worry if you don't yet understand how to use the combine() function - we'll soon
		find out! For now, we'd like to understand the Iris dataset a little better. A good
		first step is to describe() the dataset - this provides an overview of each attribute
		column in the table. What is the median petal length of all 150 iris specimens?
		""",
		"describe(iris)",
		x -> x==4.35
	),
	Activity(
		"""
		We can also extend the describe() function using additional descriptors such as the
		25th percentile (:q25) or the standard deviation (:std). Which iris attribute has the
		greatest standard deviation?
		""",
		"describe(iris,:std)",
		x -> x=="PetalLength"
	),
	Activity(
		"""
		We can also calculate these values ourselves by loading and applying the functions
		in Julia's Statistics package:

		using Statistics
		std(iris[:,"PetalLength"])

		What is the median sepal length?
		""",
		"median(iris[:,\"SepalLength\"])",
		x -> x==5.8
	),
	Activity(
		"""
		Let's investigate the covariance attributes in the Iris database. Enter this code:

		attributes = names(iris)[1:end-1]				# Leave out the non-numeric Species
		for x in attributes, y in attributes
			if x > y
				println( "\$x \t \$y \t \$(cov(iris[:,x],iris[:,y]))")
			end
		end

		Look at the results: Unsurprisingly, PetalWidth covaries with PetalLength, but
		SepalLength covaries with both PetalLength and PetalWidth. With which of these
		does SepalLength covary least strongly?
		""",
		"",
		x -> x=="PetalWidth"
	),
	Activity(
		"""
		This is what we might expect - after all, there seems no particular reason why SepalLength
		should covary strongly PetalWidth. But wait! Covariance is a measure that depends on the
		absolute size of the variables! Of PetalLength and PetalWidth, which has the greater
		mean value?
		""",
		"mean(iris[:,\"PetalLength\"]",
		x -> x=="PetalLength"
	),
	Activity(
		"""
		Hm. This is tricky. SepalLength covaries with PetalLength more strongly than with
		PetalWidth, however it may be that the smaller covariance with PetalWidth arises simply
		from the fact that PetalWidth values are much smaller than PetalLength. We can solve this
		question using correlation instead of covariance, because correlation is a relative
		measure rather than an absolute one.

		Write code to compare the two correlation values for SepalLength-PetalLength and for
		SepalLength-PetalWidth. Then find out the percentage difference between these two
		values:
		""",
		"100 * (cor(SL,PL)-cor(SL,PW))/cor(SL,PL)",
		x -> abs(x-6.173) < 0.1
	),
	Activity(
		"""
		We see that the correlation of PetalWidth is definitely major - almost as large as the
		correlation PetalLength to SepalLength. What biological function of sepals might explain
		why SepalLength correlates so highly with PetalWidth?
		""",
		"Both length and width of petals increases the bud volume that sepals must enclose.",
		x -> true
	),
	Activity(
		"""
		We can use the columns of iris as ranges for generating random values, for example:

		rand(iris[:,"SepalLength"])

		Create a Vector containing a random sample of 21 PetalWidths:
		""",
		"rand(range,10)",
		x -> length(x)==21 && all(map(x) do pw in(pw,Main.iris[:,"PetalWidth"]) end)
	),
	Activity(
		"""
		Let's create our own dataframe of the Group I chemical elements. It is a bad idea to
		clutter our global namespace with lots of trivial names, so we'll create the dataframe
		inside a function using local names:

		function group1()
			symbol = ["Li","Na","K","Rb","Cs","Fr"]
			protons = [3,11,19,37,55,87]
			nucleons = [7,23,39,85,133,223]
			electronegativity = [1.0,0.9,0.8,0.8,0.7,0.7]
			DataFrame(; symbol, protons, nucleons, electronegativity)
		end

		You may first need to add the DataFrames package and load it with the keyword "using".
		Now enter the function group1(), use it to create a dataframe named grp1, then extract
		from grp1 the subtable consisting of the 3rd and 4th rows:
		""",
		"",
		x -> size(x)==(2,4) && x[1:2,"symbol"] == ["K","Rb"]
	),
	Activity(
		"""
		Finally, we'll learn to save and load dataframes to and from a file. First we'll do it
		simply, using comma-separated variable (CSV) files. Add and load the CSV package, then
		write the group1 dataframe to a CSV file:

		CSV.write( "group1.csv", group1())

		Now give me the result of reading the csv file:

		str = read("group1.csv",String);
		""",
		"",
		x -> x[43:44] == "Li"
	),
	Activity(
		"""
		Sometimes we want to write to other file formats, such as Excel files, and Julia offers
		us a range of packages such as XLSX for doing this. However, the advantage of csv files
		for small projects is that we can read them into text strings and with a text editor.

		We can also read our csv file as a DataFrame. Give me the following dataframe:

		CSV.read("group1.csv",DataFrame)
		""",
		"",
		x -> typeof(x)==Main.DataFrame && size(x)==(6,4)
	),
]


text:
using news4
using Test

# either this if i don't export the function with the export keyword
#news4.my_f(2,3)

# or directly if there is an export statement.
my_f(2,3)

# adding a test for 2x + 3y (if that was the case)
#@test my_f(2,3) == 13

@testset "news4.jl" begin
    # Write your own tests here.
    @test my_f(2,3) == 13
    @test my_f(1,2) == 8
    @test my_f(3,4) == 18

end

ext: jl


text:
module RocketSwitchMapToOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: switch_map_to()" begin

    println("Testing: operator switch_map_to()")

    struct DummyType end

    @testset begin
        @test_throws InvalidSubscribableTraitUsageError of(0) |> switch_map_to(DummyType())
    end

    run_testset([
        (
            source      = from([ 0, 0 ]) |> switch_map_to(from([ 1, 2 ])),
            values      = @ts([ 1, 2, 1, 2, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> switch_map_to(of(0.0)),
            values      = @ts([ 0.0, 0.0, 0.0, 0.0, 0.0, c ]),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map_to(faulted(Float64, "err")),
            values      = @ts(e("err")),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map_to(completed(Float64)),
            values      = @ts(c),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map_to(never(Float64)),
            values      = @ts(),
            source_type = Float64
        ),
        (
            source      = completed() |> switch_map_to(of(0)),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(String, "e") |> switch_map_to(of(0)),
            values      = @ts(e("e")),
            source_type = Int
        ),
        (
            source      = never() |> switch_map_to(of(0)),
            values      = @ts(),
            source_type = Int
        ),
        (
            source      = from([ 0, 0 ]) |> async(0) |> switch_map_to(from([ 1, 2 ])),
            values      = @ts([ 1, 2 ] ~ [ 1, 2 ] ~ c),
            source_type = Int
        ),
        (
            source      = from([ 0, 0 ]) |> async(0) |> switch_map_to(from([ 1, 2 ]) |> async(0)),
            values      = @ts([ 1 ] ~ [ 2 ] ~ c),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> switch_map_to(of(0)),
            values      = @ts([ 0, 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> switch_map(Int) |> switch_map_to(of(0)),
            values      = @ts([ 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> switch_map_to(of(0)),
            values      = @ts([ 0, 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> switch_map(Int) |> switch_map_to(of(0)),
            values      = @ts([ 0, e("err") ]),
            source_type = Int
        )
    ])

end

end


text:
#========================================================================================#
#	Laboratory 11
#
# The course project!!!
#
# Author: Niall Palfreyman, 23/03/2022
#========================================================================================#
[
	Activity(
		"""
		This is the last laboratory in this course, and it presents the assessed project for
		the course. In this laboratory, you will implement a program to demonstrate chaos in a
		gravitational system containing three bodies of equal mass in two dimensions. You will use
		Julia to analyse the execution of a complex simulation program, and adapt this program to
		use Runge-Kutta integration to (a) simulate 3-body motion and (b) evaluate its chaoticity.
		You will present your methods and results in two deliverables: Julia demonstration code and
		a poster (a single Powerpoint slide readable from the back of the room, containing title,
		author, text blocks, lines of code and diagrams).

		As an introduction to the topic, please view this video-clip now:

		https://www.youtube.com/watch?v=LwkvO3t1b30&t=113s
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Demonstrate simple logistic chaos.
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
	Activity(
		"""
		""",
		"",
		x -> true
	),
]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random

# test the basic floating point functions

@testset "flipsign" begin
    for elty in (Float32,Float64)
        x = convert(elty,-2.0)
        x = flipsign(x,-1.0)
        @test flipsign(x,big(-1.0)) == convert(elty,-2.0)
    end
end

@testset "maxintfloat" begin
    @test maxintfloat(Float16) === Float16(2048f0)
    for elty in (Float16,Float32,Float64)
        @test maxintfloat(rand(elty)) === maxintfloat(elty)
    end
    @test maxintfloat() === maxintfloat(Float64)
    @test maxintfloat(Float64, Int32) === 2147483647.0
    @test maxintfloat(Float32, Int32) === maxintfloat(Float32)
    @test maxintfloat(Float64, Int16) === 32767.0
    @test maxintfloat(Float64, Int64) === maxintfloat(Float64)
end

@testset "isinteger" begin
    for elty in (Float16, Float32, Float64)
        @test !isinteger(elty(1.2))
        @test isinteger(elty(12))
        @test isinteger(zero(elty))
        @test isinteger(-zero(elty))
        @test !isinteger(nextfloat(zero(elty)))
        @test !isinteger(prevfloat(zero(elty)))
        @test isinteger(maxintfloat(elty))
        @test isinteger(-maxintfloat(elty))
        @test !isinteger(elty(Inf))
        @test !isinteger(-elty(Inf))
        @test !isinteger(elty(NaN))
    end
end

@testset "round" begin
    for elty in (Float32, Float64)
        x = rand(elty)
        A = fill(x,(10,10))
        @test round.(A,RoundToZero) == fill(trunc(x),(10,10))
        @test round.(A,RoundUp) == fill(ceil(x),(10,10))
        @test round.(A,RoundDown) == fill(floor(x),(10,10))
        A = fill(x,(10,10,10))
        @test round.(A,RoundToZero) == fill(trunc(x),(10,10,10))
        @test round.(A,RoundUp) == fill(ceil(x),(10,10,10))
        @test round.(A,RoundDown) == fill(floor(x),(10,10,10))
        for elty2 in (Int32,Int64)
            A = fill(x,(10,))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,))
            A = fill(x,(10,10))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10))
            A = fill(x,(10,10,10))
            @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10,10))
            @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10,10))
            @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10,10))
            @test round.(elty2,A) == fill(round(elty2,x),(10,10,10))
        end
    end
end

@testset "Types" begin
    for x in (Int16(0), 1, 2f0, pi, 3//4, big(5//6), 7.8, big(9), big(ℯ))
        @test float(typeof(x)) == typeof(float(x))
        @test float(typeof(complex(x, x))) == typeof(float(complex(x, x)))
    end
end

@testset "significant digits" begin
    # (would be nice to have a smart vectorized
    # version of signif)
    @test round(123.456, sigdigits=1) ≈ 100.
    @test round(123.456, sigdigits=3) ≈ 123.
    @test round(123.456, sigdigits=5) ≈ 123.46
    @test round(123.456, sigdigits=8, base = 2) ≈ 123.5
    @test round(123.456, sigdigits=2, base = 4) ≈ 128.0
    @test round(0.0, sigdigits=1) === 0.0
    @test round(-0.0, sigdigits=1) === -0.0
    @test round(1.2, sigdigits=2) === 1.2
    @test round(1.0, sigdigits=6) === 1.0
    @test round(0.6, sigdigits=1) === 0.6
    @test round(7.262839104539736, sigdigits=2) === 7.3
    @test isinf(round(Inf, sigdigits=3))
    @test isnan(round(NaN, sigdigits=3))
    @test round(1.12312, sigdigits=1000) === 1.12312
    @test round(Float32(7.262839104539736), sigdigits=3) === Float32(7.26)
    @test round(Float32(7.262839104539736), sigdigits=4) === Float32(7.263)
    @test round(Float32(1.2), sigdigits=3) === Float32(1.2)
    @test round(Float32(1.2), sigdigits=5) === Float32(1.2)
    @test round(Float16(0.6), sigdigits=2) === Float16(0.6)
    @test round(Float16(1.1), sigdigits=70) === Float16(1.1)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
module PandorasBoxDreamDashController

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/dreamDashController" DreamDashController(x::Integer, y::Integer, allowSameDirectionDash::Bool=false, allowDreamDashRedirect::Bool=true, overrideDreamDashSpeed::Bool=false, neverSlowDown::Bool=false, useEntrySpeedAngle::Bool=false, bounceOnCollision::Bool=false, stickOnCollision::Bool=false, overrideColors::Bool=false, sameDirectionSpeedMultiplier::Number=1.0, dreamDashSpeed::Number=240.0, activeBackColor::String="Black", disabledBackColor::String="af2e2d", activeLineColor::String="White", disabledLineColor::String="6a8480", particleLayer0Colors::String="ffef11,ff00d0,08a310", particleLayer1Colors::String="5fcde4,7fb25e,e0564c", particleLayer2Colors::String="5b6ee1,CC3B3B,7daa64")

const placements = Ahorn.PlacementDict(
    "Dream Dash Controller (Pandora's Box)" => Ahorn.EntityPlacement(
        DreamDashController
    )
)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

Ahorn.editingOptions(entity::DreamDashController) = Dict{String, Any}(
    "activeBackColor" => colors,
    "disabledBackColor" => colors,
    "activeLineColor" => colors,
    "disabledLineColor" => colors
)

const iconMatrix = [
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 1 1 2 2 2 5 2 2 4 4 4 2 1 1 0 0 0 0 0
    0 0 0 0 1 1 2 2 2 2 2 2 4 4 2 4 4 2 1 1 0 0 0 0
    0 0 0 1 2 2 2 7 2 6 2 2 2 4 4 4 2 2 2 2 1 0 0 0
    0 0 1 7 2 2 2 2 6 6 6 2 2 2 4 2 2 2 6 2 2 1 0 0
    0 1 1 2 2 2 2 2 2 6 2 2 2 2 2 2 2 2 2 2 2 1 1 0
    0 1 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 7 2 2 1 0
    1 2 2 2 2 4 2 2 1 1 2 1 1 2 1 1 2 2 7 7 7 2 2 1
    1 2 2 5 2 2 2 2 1 1 1 1 1 1 1 1 2 7 7 2 7 7 2 1
    1 2 5 5 5 2 2 2 3 1 1 1 1 1 1 3 2 2 7 7 7 2 2 1
    1 2 2 5 2 2 2 1 1 1 1 3 3 1 1 1 1 2 2 7 2 2 2 1
    1 2 2 2 2 2 2 1 1 1 1 2 2 1 1 1 1 2 2 2 2 2 2 1
    1 2 8 2 2 2 2 3 3 1 1 1 1 1 1 3 3 2 2 2 2 2 2 1
    1 8 8 8 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 7 2 2 2 1
    1 2 8 2 2 2 2 2 1 1 3 1 1 3 1 1 2 2 2 2 2 2 2 1
    1 2 2 2 2 2 7 2 3 3 2 1 1 2 3 3 2 2 2 8 2 2 2 1
    1 2 2 2 2 2 2 2 2 2 2 3 3 2 2 2 2 2 2 2 2 2 2 1
    0 1 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 2 2 5 1 0
    0 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 2 2 1 1 0
    0 0 1 2 2 7 2 2 2 2 2 2 2 2 2 2 2 4 2 2 2 1 0 0
    0 0 0 1 2 2 2 2 2 7 2 2 2 2 6 2 2 2 2 2 1 0 0 0
    0 0 0 0 1 1 2 2 7 7 7 2 2 2 2 2 2 2 1 1 0 0 0 0
    0 0 0 0 0 1 1 2 2 7 2 2 2 2 2 2 2 1 1 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
]

const iconSurface = Ahorn.matrixToSurface(
    iconMatrix,
    [
        (0 / 255, 0 / 255, 0 / 255, 0.0),
        (255 / 255, 255 / 255, 255 / 255, 1.0), # White
        (0 / 255, 0 / 255, 0 / 255, 1.0), # Black
        (148 / 255, 165 / 255, 165 / 255, 1.0), # Gray
        (0 / 255, 161 / 255, 4 / 255, 1.0), # Green
        (243 / 255, 242 / 255, 5 / 255, 1.0), # Yellow
        (247 / 255, 0 / 255, 213 / 255, 1.0), # Magenta
        (84 / 255, 207 / 255, 231 / 255, 1.0), # Cyan
        (212 / 255, 79 / 255, 68 / 255, 1.0), # Red

    ]
)

sprite = "objects/pandorasBox/DreamDashController/DreamDashController0"

function Ahorn.selection(entity::DreamDashController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamDashController, room::Maple.Room) = Ahorn.drawImage(ctx, iconSurface, -12, -12)

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
@testset "PDA" begin
    p = @pdata [9, 9, 8, missing, 1, 1]
    pcopy = copy(p)
    @test levels(p) == [1, 8, 9]
    @test levels(setlevels(p, ["a", "b", "c"])) == ["a", "b", "c"]
    @test collect(skipmissing(setlevels(p, (@data ["a", "b", missing])))) == ["b", "a", "a"]
    @test collect(skipmissing(setlevels(p, (@data ["a", "b", "a"])))) == ["a", "a", "b", "a", "a"]
    @test levels(setlevels(p, (@data ["a", "b", "a"]))) == ["a", "b"]
    @test levels(setlevels(p, Dict([(1, 111)]))) == [111, 8, 9]
    @test levels(setlevels(p, Dict([(1, 111), (8, missing)]))) == [111, 9]
    @test levels(PooledDataArray(p, [9, 8, 1])) == [9, 8, 1]
    @test levels(PooledDataArray(p, [9, 8])) == [9, 8]
    @test collect(skipmissing(PooledDataArray(p, [9, 8]))) == [9, 9, 8]
    @test levels(PooledDataArray(p, levels(p)[[3,2,1]])) == [9,8,1]
    v = collect(1:6)
    @test isequal(p, reorder(p))
    # @test levels(reorder(p, v)) == [9,8,1]
    @test isequal(p, pcopy)

    @test levels(setlevels!(copy(p), [10,80,90])) == [10, 80, 90]
    @test levels(setlevels!(copy(p), [1,8,1])) == [1, 8]
    @test levels(setlevels!(copy(p), (@data [1, 8, missing]))) == [1, 8]
    @test levels(setlevels!(copy(p), [1,8,9, 10])) == [1, 8, 9, 10]
    @test levels(setlevels!(copy(p), Dict([(1, 111)]))) == [111, 8, 9]
    @test levels(setlevels!(copy(p), Dict([(1, 111), (8, missing)]))) == [111, 9]
    # issue #201
    @test levels(setlevels!(@pdata([1.0, 2.0]), [3,4])) == [3.0, 4.0]

    y = @pdata [1, missing, -2, 1, missing, 4, missing]
    @test isequal(unique(y), @pdata [1, missing, -2, 4])
    @test isequal(unique(reverse(y)), @data [missing, 4, 1, -2])
    @test isequal(unique(skipmissing(y)), @data [1, -2, 4])
    @test isequal(unique(reverse(collect(skipmissing(y)))), @data [4, 1, -2])

    z = @pdata ["frank", missing, "gertrude", "frank", missing, "herbert", missing]
    @test isequal(unique(z), @pdata ["frank", missing, "gertrude", "herbert"])
    @test isequal(unique(reverse(z)), @pdata [missing, "herbert", "frank", "gertrude"])
    @test isequal(unique(skipmissing(z)), @pdata ["frank", "gertrude", "herbert"])
    @test isequal(unique(reverse(collect(skipmissing(z)))), @pdata ["herbert", "frank", "gertrude"])

    # check case where only missing occurs in final position
    @test isequal(unique(@pdata [1, 2, 1, missing]), @pdata [1, 2, missing])

    pp = PooledDataArray(Any[])
    @test length(pp) == 0
    @test length(levels(pp)) == 0

    # test construction with unordered types
    pim = @pdata [1 + im, 2 + im, 3 + im, 2 + im, 1 + im]
    @test levels(pim) == [1 + im, 2 + im, 3 + im]

    # Test explicitly setting refs type
    testarray = [1, 1, 2, 2, 0, 0, 3, 3]
    testdata = @data [1, 1, 2, 2, 0, 0, 3, 3]
    for t in Any[testarray, testdata]
        for R in [UInt8, UInt16, UInt32, UInt64]
            @test eltype(PooledDataArray(t, R).refs) == R
            @test eltype(PooledDataArray(t, [1,2,3], R).refs) == R
            @test eltype(PooledDataArray(t, [1,2,3], t .== 0, R).refs) == R
        end
    end

    pcopy = copy(p)
    @test levels(append!(pcopy, @pdata [4, missing, 6, 5])) == [1, 8, 9, 4, 5, 6]

    x = PooledDataArray([9, 9, 8])
    y = PooledDataArray([1, 9, 3, 2, 2])
    @test append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = PooledDataArray([9, 9, 8])
    y = [1, 9, 3, 2, 2]
    @test append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    # convert methods
    for from in (@pdata(ones(5, 5)), @data(ones(5, 5)), ones(5, 5))
        for (to, totype) in ((PooledDataArray{Float32,UInt16,2}, PooledDataArray{Float32,UInt16,2}),
                             (PooledDataArray{Float32,UInt32,2}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64,UInt16,2}, PooledDataArray{Float64,UInt16,2}),
                             (PooledDataArray{Float64,UInt32,2}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray{Float32,UInt16}, PooledDataArray{Float32,UInt16,2}),
                             (PooledDataArray{Float32,UInt32}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64,UInt16}, PooledDataArray{Float64,UInt16,2}),
                             (PooledDataArray{Float64,UInt32}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray{Float32}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray, PooledDataArray{Float64,UInt32,2}))
            rettype = typeof(convert(to, from))
            if rettype != totype
                error("convert($to, ::$(typeof(from))) returned $rettype (expected $totype)")
            end
        end
    end

    da = convert(DataArray, @pdata(ones(5, 5)))
    @test isequal(da, @data(ones(5, 5)))
    @test typeof(da) == DataArray{Float64,2}
    da = convert(DataArray{Float32}, @pdata(ones(5, 5)))
    @test isequal(da, @data(ones(Float32, 5, 5)))
    @test typeof(da) == DataArray{Float32,2}
    da = convert(DataArray{Float32,2}, @pdata(ones(5, 5)))
    @test isequal(da, @data(ones(Float32, 5, 5)))
    @test typeof(da) == DataArray{Float32,2}

    # permute
    pda = @pdata([missing, "A", "B", "C", "A", "B"])
    @test isequal(Base.permute!!(copy(pda), [2, 5, 3, 6, 4, 1]), @pdata(["A", "A", "B", "B", "C", missing]))
    @test isequal(Base.ipermute!!(copy(pda), [6, 1, 3, 5, 2, 4]), @pdata(["A", "A", "B", "B", "C", missing]))

    a1 = 1:200
    a2 = 100:300
    pa1 = PooledDataArray(a1)
    pa2 = PooledDataArray(a2)
    ca1 = compact(pa1)
    ca2 = compact(pa2)
    r = vcat(ca1, ca2)
    @test r == vcat(a1, a2)
    @test isa(r, PooledDataArray{Int,DataArrays.DEFAULT_POOLED_REF_TYPE})
    @test isa(vcat(ca1, pa2), PooledDataArray{Int,DataArrays.DEFAULT_POOLED_REF_TYPE})

    a1 = Array{Int}(undef, 2,3,4,5)
    a2 = Array{Int}(undef, 3,3,4,5)
    a1[1:end] = length(a1):-1:1
    a2[1:end] = (1:length(a2)) + 10
    ca1 = compact(PooledDataArray(a1))
    ca2 = compact(PooledDataArray(a2))
    r = vcat(ca1, ca2)
    @test r == vcat(a1, a2)
    @test isa(r, PooledDataArray{Int,DataArrays.DEFAULT_POOLED_REF_TYPE,4})

    # Issue #265:
    #   Ensure that any levels that have are not expressed in a
    #   PooledDataArray are handled correctly and not left uninitialised when
    #   using unique()
    x = ["A", "B", "A"];
    masks = [[1], [2], [3], [1, 3]]
    for mask in masks
        y = PooledDataArray(x)
        y[mask] = missing
        @test isequal(sort(unique(y)), sort(DataArray(unique(y))))
    end
    z = PooledDataArray([1, 2], [1, 2, 3])
    @test sort(unique(z)) == DataArray([1, 2])

    # Issue #273
    #   Ensure that PooledDataArray works for arrays with non leaf type elements.
    #   Those two examples used to throw exceptions.
    x = AbstractString["a"]
    m = Bool[false]
    r = UInt32
    @test isa(PooledDataArray(x, m, r), DataArrays.PooledDataArray{AbstractString, UInt32, 1})
    x = Integer[1]
    @test isa(PooledDataArray(x, m, r), DataArrays.PooledDataArray{Integer, UInt32, 1})
end


text:
function optimize(f::Function,
				  lb, ub,
				  method::AbstractMetaheuristic;
				  options = Options())

    # Parameters
    it = 0
	it_unchanged, old_obj = 0, NaN
	timedout, t0 = false, time()

    # Initialize pop
	println("Population initialization...")
    initialize_population!(f, lb, ub, method, options)

    # Algorithm
	println("Starting optimization...")
    while it < options.iterations && it_unchanged < options.itmax_unchanged && !timedout
		# Iteration
        it += 1

		# Update population
		update_population!(f, lb, ub, method, options)

		# Check convergence
		it_unchanged = isconverged!(method.population[1,end], old_obj, options.reltol, it_unchanged)

		# Store the objective value for convergence checking
		old_obj = method.population[1,end]

		# Check timing
		timedout = istimedout(time(), t0, options.time_limit)

		# Verbose
		options.verbose ? show_verbose(it, time() - t0, method.population[1,end]) : nothing
    end

	# Results
	results = MetaheuristicResults(method,
	get_status(it, it_unchanged, timedout, options),
	method.population[1,1:length(lb)],
	method.population[1,end],
	it,
	options)

	# Print results to the REPL
	show_results(results)

    return results
end


text:
function test_deferrableLoad()

    T = 2
    load = [2.0, 3.0]

    pmin = zeros(T)
    pmax = ones(T)

    cycle_begin = [1]
    cycle_end = [T]

    emin = [1.2]
    emax = [1.5]

    l = DR.DER.DeferrableLoad(
        index=0, T=T, ncycles=1,
        cycle_begin=cycle_begin, cycle_end=cycle_end,
        cycle_energy_min=emin, cycle_energy_max=emax,
        pwr_min=pmin, pwr_max=pmax,
        binflag=false
    )

    # GLPK does not support -in-Interval constraints, so we bridge
    model = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

    lmin = [0.0, 0.0]
    lmax = [10.0, 10.0]
    price = [1.0, -0.01]

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [l]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.2 - 0.01

    # Check net load
    pnet1 = var2idx[(:house, h.index, :pnet, 1)]
    pnet2 = var2idx[(:house, h.index, :pnet, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), pnet1) ≈ 0.2
    @test MOI.get(model, MOI.VariablePrimal(), pnet2) ≈ 1.0

    return nothing

end

@testset "Def." begin
    test_deferrableLoad()
end


text:
# Run the notebook tests in separate modules

ENV["JUAFEM_TESTING"] = true

module TestStiffness
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/stiffness_example.ipynb")
        end
    end
end

module Cantilever
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/cantilever.ipynb")
        end
    end
end

module HyperElasticity
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/hyperelasticity.ipynb")
        end
    end
end

module HeatSquare
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/heat_square.ipynb")
        end
    end
end

module Helmholtz
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/helmholtz.ipynb")
        end
    end
end

module Cook
    using NBInclude
    mktempdir() do dir
        cd(dir) do
            nbinclude("../examples/cooks_membrane_mixed_up.ipynb")
        end
    end
end


text:
abstract type Space{T, D} end
struct InfNorm{T, D, δ} <: Space{T, D} end
struct Product{T, D1, D2, D3, S1 <: Space{T, D1}, S2 <: Space{T, D2}} <: Space{T, D3}
	s1::S1
	s2::S2
	Product(s1::S1, s2::S2) where {T, D1, D2, S1 <: Space{T, D1}, S2 <: Space{T, D2}} = new{T, D1, D2, D1+D2, S1, S2}(s1, s2)
end
abstract type Cone{T, D} <: Space{T, D} end
dim(::Cone{T,D}) where {T,D} = D
dim(::Type{T}) where {D, Tc, T<:Cone{Tc,D}} = D
struct Reals{T, D} <: Cone{T, D} end # R^D
struct Zeros{T, D} <: Cone{T, D} end # 0^D
struct POCone{T, D} <: Cone{T, D} end # { x | forall i, x_i >= 0}
struct NOCone{T, D} <: Cone{T, D} end # { x | forall i, x_i <= 0}
struct SOCone{T, D} <: Cone{T, D} end # { [t,x] | |x| <= t }
struct NSOCone{T, D} <: Cone{T, D} end # { [t,x] | |x| <= -t }

struct PTCone{T, Cs<:Tuple{Vararg{C where C<:Cone{T}}}, D} <: Cone{T, D}
	cones::Cs
	@generated PTCone{T}(cs::Cs) where {T, Cs <: Tuple{Vararg{C where C<:Cone{T}}}} = Expr(:new, PTCone{T, Cs, sum(dim.(Cs.parameters); init=0)}, :(cs))
end
cone_dim(::Cone{T,D}) where {T,D} = D
cone_dim(::Type{C}) where {T, D, C<:Cone{T,D}} = D

# projections
project(::Reals{T, D}, x::SVector{D, T}) where {D, T} = x
project(::Zeros{T, D}, x::SVector{D, T}) where {D, T} = zeros(SVector{D})
project(::InfNorm{T, D, δ}, x::SVector{D, T}) where {D, T, δ} = clamp.(x, -δ, δ)

function project(c::POCone{T, D}, x::SVector{D, T}) where {D, T}
	return max.(x, zero(T))
end

function project(c::NOCone{T, D}, x::SVector{D, T}) where {D, T}
	return min.(x, zero(T))
end

@generated function project(c::SOCone{T, D}, x::SVector{D, T}) where {T, D}
	onev = one(T)
	two = 2*onev
	zv = zeros(SVector{D,T})
	vect_inds = SVector{D-1}(2:D)
	return quote
		xnorm = norm(x[$vect_inds])
		r = x[1]
		if xnorm <= r
			return x
		elseif xnorm <= -r
			return $zv
		else
			scalefact = (xnorm + r)/($two)
			component_factor = (scalefact)/xnorm
			return vcat(SVector(scalefact), component_factor * x[$vect_inds])
		end
	end
end

@generated function project(c::NSOCone{T, D}, x::SVector{D, T}) where {T, D}
	vect_inds = SVector{D-1}(2:D)
	return quote
		po_equiv = project(SOCone{T,D}(), vcat(-x[1], x[$vect_inds]))
		return vcat(-po_equiv[1], po_equiv[$vect_inds])
	end
end

@generated function project(c::PTCone{T, Cs, D}, x::SVector{D, T}) where {T, Cs, D}
	if length(Cs.parameters) == 0 || D == 0 return :(SVector{0,T}()) end
	idxes = prepend!(accumulate(+, dim.(Cs.parameters); init=1), 1)
	ranges = map(x->UnitRange(x...), zip(idxes, Iterators.drop(idxes, 1) .- 1))
	projections = [:(project(c.cones[$i], x[SVector{$(dim(Cs.parameters[i])), Int64}($(ranges[i]))])) for i in 1:length(Cs.parameters)]
	return :(vcat($(projections...)))
end

@generated function project(c::Product{T, D1, D2, D3, C1, C2}, x::SVector{D3, T}) where {D1, D2, D3, C1, C2, T}
	r1 = SVector{D1, Int64}(1:D1)
	r2 = SVector{D2, Int64}(D1+1:D3)
	return :(vcat(project(c.s1, x[$r1]), project(c.s2, x[$r2])))
end

# polars
polar(::Reals{T,D}) where {T,D} = Zeros{T,D}()
polar(::Zeros{T,D}) where {T,D} = Reals{T,D}()
polar(::POCone{T,D}) where {T,D} = NOCone{T,D}()
polar(::NOCone{T,D}) where {T,D} = POCone{T,D}()
polar(::SOCone{T,D}) where {T,D} = NSOCone{T,D}()
polar(::NSOCone{T,D}) where {T,D} = SOCone{T,D}()

@generated function polar(c::PTCone{T, Cs, D}) where {T, Cs, D}
	tup = Expr(:tuple, (:(polar(c.cones[$i])) for i in 1:length(Cs.parameters))...)
	return :(PTCone{T}($tup))
end

# problem
# minimize 1/2 z^t P z + q^T z
# s.t. H z - g ∈ K, z ∈ D

struct Problem{T,N,M,K<:Cone{T,M},D<:Space{T,N}}
	k::K
	d::D
	H::SparseMatrixCSC{T, Int}
	P::SparseMatrixCSC{T, Int}
	q::MVector{N, T}
	g::MVector{M, T}
	c::T
	Problem(k::K, d::D, H::SparseMatrixCSC{T, Int}, P::SparseMatrixCSC{T, Int}, q::MVector{N, T}, g::MVector{M, T}, c::T) where {T,N,M,K<:Cone{T,M},D<:Space{T,N}} =
		new{T,N,M,K,D}(k,d,transpose(H),P,q,g,c)
end

abstract type Diagnostics{T} end
struct NoDiagnostics{T} <: Diagnostics{T} end
struct LogDiagnostics{T} <: Diagnostics{T} end
record_diagnostics(d::NoDiagnostics{T}, i, w, z, v, w_delta, z_delta) where T = nothing
function record_diagnostics(d::LogDiagnostics{T}, i, w, z, v, w_delta, z_delta) where T
	if i == 1
		println("iteration | w | w_delta | z | z_delta")
	end
	if i % 10000 == 0
		println("$i | $w | $w_delta | $z | $z_delta")
	end
end

abstract type Scaling{T, M, N} end
struct GeoMean{T, M, N} <: Scaling{T, M, N}
	colmax::MVector{N, T}
	colmin::MVector{N, T}
	GeoMean(::Problem{T,N,M}) where {T,M,N} = new{T, M, N}(MVector{N, T}(zeros(T, N)), MVector{N, T}(zeros(T, N)))
end
struct ArithMean{T, M, N} <: Scaling{T, M, N}
	colsum::MVector{N, T}
	ArithMean(::Problem{T,N,M}) where {T,M,N} = new{T, M, N}(MVector{N, T}(zeros(T, N)))
end
struct Equilibration{T, M, N} <: Scaling{T, M, N}
	colmax::MVector{N, T}
	Equilibration(::Problem{T,N,M}) where {T,M,N} = new{T, M, N}(MVector{N, T}(zeros(T, N)))
end

@enum SolverState INDEFINITE OPTIMAL PRIMAL_INFEASIBLE DUAL_INFEASIBLE TIMEOUT
mutable struct State{T,N,M,K,D,P <: Problem{T,N,M,K,D}, G <: Diagnostics{T}, SCs <: Tuple{Vararg{<:Scaling{T, M, N}}}}
	w_i1::MVector{M, T}
	z_i1::MVector{N, T}
	z_i2::MVector{N, T}
	v_i1::MVector{M, T}

	col_scale::MVector{N, T}
	row_scale::MVector{M, T}

	solver_state::SolverState
	primal::MVector{N, T}
	dual::MVector{M, T}

	diagnostics::G
	scaling::SCs
	function State(p::P; diag::G = NoDiagnostics{T}(), scaling::SCs = (GeoMean(p), Equilibration(p))) where {
			T,N,M,K,D,
			P<:Problem{T,N,M,K,D},
			G<:Diagnostics{T},
			SCs<:Tuple{Vararg{Scaling{T, M, N}}}}
	return new{T,N,M,K,D,P,G,SCs}(
		MVector{M,T}(zeros(T, M)), MVector{N,T}(zeros(T, N)), MVector{N,T}(zeros(T, N)), MVector{M,T}(zeros(T, M)),
		MVector{N,T}(ones(T, N)), MVector{M, T}(ones(T, M)),
		INDEFINITE, MVector{N, T}(zeros(T, N)), MVector{M, T}(zeros(T, M)), diag, scaling)
	end
end

function objective_value(p::P, s::State{T,N,M,K,D,P}) where {T,N,M,K,D,P<:Problem{T,N,M,K,D}}
	return transpose(s.primal) * p.P * s.primal * 0.5 + dot(s.primal, p.q)
end
function dual_objective_value(p::P, s::State{T,N,M,K,D,P}) where {T,N,M,K,D,P<:Problem{T,N,M,K,D}}
	return -transpose(s.primal) * p.P * s.primal * 0.5 - dot(s.dual, p.g)
end

function compute_α(p::P, γ::T) where {T,N,M,K,D,P<:Problem{T,N,M,K,D}}
	λ = max_singular_value(p.P)
	ν = max_singular_value(p.H)
	return (8.0-4.0/γ)/(sqrt(λ^2 + 16*ν^2) + λ)
end

function spmul!(res::MVector{N,T}, A::SparseMatrixCSC{T, Int}, B::SVector{M,T}, α::T, β::T) where {T<:Number,N,M}
	@assert N == A.m
	@assert M == A.n
	res .*= β
	for j in 1:A.n
		vv = B[j] * α
		for i in A.colptr[j]:A.colptr[j+1]-1
			@inbounds res[A.rowval[i]] += A.nzval[i] * vv
		end
	end
end
function spmul!(res::MVector{N,T}, At::Transpose{T, SparseMatrixCSC{T, Int}}, B::SVector{M,T}, α::T, β::T) where {T<:Number,N,M}
	A = At.parent
	@assert N == A.n
	@assert M == A.m
	res .*= β
	for j in 1:A.n
		acc = zero(T)
		for i in A.colptr[j]:A.colptr[j+1]-1
			@inbounds acc += A.nzval[i] * B[A.rowval[i]] * α
		end
		@inbounds res[j] += acc
	end
end

function pipg(p::P,s::State{T,N,M,K,D,P,G}, iters::Int, α::T, ϵ::T, z::SVector{N,T}, v::SVector{M,T}) where {T,N,M,K,D,P<:Problem{T,N,M,K,D}, G<:Diagnostics{T}}
	w = v
	w_delta = zero(T)
	z_delta = zero(T)
	w_prev_delta = zero(T)
	z_prev_delta = zero(T)
	niters = 0
	for i=1:iters
		w_prev = w
		# w_raw = v + α*(transpose(p.H) * z - p.g)
		s.w_i1 .= p.g
		spmul!(s.w_i1, transpose(p.H), z, α, -α)
		w_raw = v + s.w_i1
		w = project(polar(p.k), w_raw)
		z_prev = z
		# z_raw = z - α*(p.P * z + p.q + p.H * w)
		s.z_i1 .= p.q
		spmul!(s.z_i1, p.P, z, α, α)
		spmul!(s.z_i1, p.H, w, α, 1.0)
		z_raw = z - s.z_i1
		z = project(p.d, z_raw)
		# v = w + α * transpose(p.H) * (z - z_prev)
		spmul!(s.v_i1, transpose(p.H), z - z_prev, α, 0.0)
		v = w + s.v_i1

		w_prev_delta = w_delta
		z_prev_delta = z_delta
		w_delta = norm(w - w_prev)
		z_delta = norm(z - z_prev)
		record_diagnostics(s.diagnostics, i, w, z, v, w_delta, z_delta)
		niters = i
		if w_delta < ϵ && z_delta < ϵ
			break
		end
		if isnan(w_delta) || isnan(z_delta)
			# halt iteration with the previous values of w and z
			w = w_prev
			z = z_prev
			w_delta = w_prev_delta
			z_delta = z_prev_delta
			break
		end
	end
	if w_delta < ϵ && z_delta < ϵ
		s.solver_state = OPTIMAL
		s.primal .= z
		s.dual .= v
		return niters
	elseif w_delta > ϵ
		s.solver_state = PRIMAL_INFEASIBLE
		s.dual .= w ./ norm(w)
		return niters
	elseif z_delta > ϵ
		s.solver_state = DUAL_INFEASIBLE
		s.primal .= z ./ norm(z)
		return niters
	end
	s.solver_state = TIMEOUT
	s.primal .= z
	return niters
end

function max_singular_value(m::SparseMatrixCSC{T, Int}; iters=10, ϵ=1e-8) where {T}
	rows, cols = size(m)
	vect = rand(cols)
	temp = zeros(rows)
	for i=1:iters
		if norm(vect) <= ϵ return norm(vect) end
		mul!(temp, m, vect)
		mul!(vect, transpose(m), temp, 1/norm(vect), 0.0)
	end
	return sqrt(norm(vect))
end
#=
function test()
	P_ex = spzeros(2,2)
	q = MVector(-50.0, -120.0)
	H_ex = sparse([100.0 200; 10 30; 1 1; -1 0; 0 -1])
	g = MVector(10000.0, 1200, 110, -0.0, -0.0)
	prob = PIPG.Problem(PIPG.NOCone{Float64, 5}(), PIPG.Reals{Float64, 2}(), H_ex, P_ex, q, g, 0.0)
	state = PIPG.State(prob)

	λ = PIPG.max_singular_value(P_ex)
	ν = PIPG.max_singular_value(H_ex)
	γ = 0.9
	α = (8.0-4.0/γ)/(sqrt(λ^2 + 16*ν^2) + λ)
	result = @time PIPG.pipg(prob, state, 2000000, α, 0.0001, SVector(1.0, 1.0), SVector(0.0,0.0,0.0,0,0))
	println(result)
end
=#


text:
using CoupledSystems
using Documenter

makedocs(;
    modules=[CoupledSystems],
    pages=[
        "Home" => "index.md",
        "Getting Started" => "guide.md",
        "Theory" => "theory.md",
        "Library" => "library.md"
    ],
    sitename="CoupledSystems.jl",
    authors="Taylor McDonnell <taylor.golden.mcdonnell@gmail.com> and contributors",
)

deploydocs(
    repo = "github.com/byuflowlab/CoupledSystems.jl.git",
    devbranch = "main",
)


text:
function B_0_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_equation = Equation(0, 1, 0, cur_G_T)

  cur_equation
end


text:
using GRAPE
using Pkg
using Documenter
using Plots

gr()
ENV["GKSwstype"] = "100"
ENV["GRAPE_LINESEARCH_ANALYSIS_VERBOSE"] = "0"

include(joinpath(@__DIR__, "download_dumps.jl"))

# Generate examples
include("generate.jl")

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

PROJECT_TOML = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))
VERSION = PROJECT_TOML["version"]
NAME = PROJECT_TOML["name"]
AUTHORS = join(PROJECT_TOML["authors"], ", ") * " and contributors"
GITHUB = "https://github.com/JuliaQuantumControl/GRAPE.jl"

println("Starting makedocs")

makedocs(;
    modules=[GRAPE],
    authors=AUTHORS,
    repo="$GITHUB/blob/{commit}{path}#{line}",
    sitename="GRAPE.jl",
    format=Documenter.HTML(;
        prettyurls = true,
        canonical  = "https://juliaquantumcontrol.github.io/GRAPE.jl",
        assets     = String[],
        mathengine = KaTeX(),
        footer     = "[$NAME.jl]($GITHUB) v$VERSION docs powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl)."
    ),
    pages=[
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Examples" => [
            "List of Examples" => "examples/index.md",
            "Example 1 (TLS)" => "examples/simple_state_to_state.md",
            "Example 2 (PE)" => "examples/perfect_entanglers.md",
        ],
        "API" => "api.md",
    ]
)

println("Finished makedocs")

rm(joinpath(@__DIR__, "build", "examples", ".gitignore"))

deploydocs(; repo="github.com/JuliaQuantumControl/GRAPE.jl", devbranch="master")


text:
include("test_actorsworld.jl")
include("test_kara_01.jl")
include("test_xml.jl")
include("test_kara_gui.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1VolumeNodeAffinity <: SwaggerModel
    required::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSelector } # spec name: required

    function IoK8sApiCoreV1VolumeNodeAffinity(;required=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1VolumeNodeAffinity, Symbol("required"), required)
        setfield!(o, Symbol("required"), required)
        o
    end
end # type IoK8sApiCoreV1VolumeNodeAffinity

const _property_map_IoK8sApiCoreV1VolumeNodeAffinity = Dict{Symbol,Symbol}(Symbol("required")=>Symbol("required"))
const _property_types_IoK8sApiCoreV1VolumeNodeAffinity = Dict{Symbol,String}(Symbol("required")=>"IoK8sApiCoreV1NodeSelector")
Base.propertynames(::Type{ IoK8sApiCoreV1VolumeNodeAffinity }) = collect(keys(_property_map_IoK8sApiCoreV1VolumeNodeAffinity))
Swagger.property_type(::Type{ IoK8sApiCoreV1VolumeNodeAffinity }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1VolumeNodeAffinity[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1VolumeNodeAffinity }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1VolumeNodeAffinity[property_name]

function check_required(o::IoK8sApiCoreV1VolumeNodeAffinity)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1VolumeNodeAffinity }, name::Symbol, val)
end


text:
@testset "Test Objects" begin
    int_obj = m.IntegerObj(123)
    str_obj = m.StringObj("str")
    true_obj = m._TRUE
    false_obj = m._FALSE
    null_obj = m._NULL
    error_obj = m.ErrorObj("error")
    arr_obj = m.ArrayObj([m.IntegerObj(1), m.IntegerObj(2)])
    hash_obj = m.HashObj(
        Dict(m.IntegerObj(1) => m.IntegerObj(2), m.IntegerObj(3) => m.IntegerObj(4)),
    )
    function_obj = m.FunctionObj(
        [m.Identifier(m.Token(m.IDENT, "x"), "x")],
        m.BlockStatement(
            m.Token(m.LBRACE, "{"),
            [
                m.ExpressionStatement(
                    m.Token(m.IDENT, "x"),
                    m.InfixExpression(
                        m.Token(m.PLUS, "+"),
                        m.Identifier(m.Token(m.IDENT, "x"), "x"),
                        "+",
                        m.IntegerLiteral(m.Token(m.INT, "2"), 2),
                    ),
                ),
            ],
        ),
        m.Environment(),
    )
    macro_obj = m.MacroObj(
        [m.Identifier(m.Token(m.IDENT, "x"), "x")],
        m.BlockStatement(m.Token(m.LBRACE, "{"), []),
        m.Environment(),
    )
    len = m.get_builtin_by_name("len")
    return_value = m.ReturnValue(m._TRUE)
    quote_obj = m.QuoteObj(m.NullLiteral(m.Token(m.NULL, "null")))
    compiled_function_obj = m.CompiledFunctionObj(m.Instructions([]), 0, 0, true)
    closure_obj = m.ClosureObj(compiled_function_obj, [])

    @testset "Is Truthy" begin
        @test m.is_truthy(m.IntegerObj(1))
        @test m.is_truthy(m.IntegerObj(0))
        @test m.is_truthy(int_obj)
        @test m.is_truthy(m.StringObj(""))
        @test m.is_truthy(m.StringObj("0"))
        @test m.is_truthy(m.StringObj("1"))
        @test m.is_truthy(str_obj)
        @test m.is_truthy(true_obj)
        @test !m.is_truthy(false_obj)
        @test !m.is_truthy(null_obj)
        @test m.is_truthy(error_obj)
        @test m.is_truthy(arr_obj)
        @test m.is_truthy(hash_obj)
        @test m.is_truthy(function_obj)
        @test m.is_truthy(macro_obj)
        @test m.is_truthy(len)
        @test m.is_truthy(return_value)
        @test m.is_truthy(quote_obj)
        @test m.is_truthy(compiled_function_obj)
        @test m.is_truthy(closure_obj)
    end

    @testset "Get Type" begin
        @test m.type_of(int_obj) == "INTEGER"
        @test m.type_of(str_obj) == "STRING"
        @test m.type_of(true_obj) == "BOOLEAN"
        @test m.type_of(false_obj) == "BOOLEAN"
        @test m.type_of(null_obj) == "NULL"
        @test m.type_of(error_obj) == "ERROR"
        @test m.type_of(arr_obj) == "ARRAY"
        @test m.type_of(hash_obj) == "HASH"
        @test m.type_of(function_obj) == "FUNCTION"
        @test m.type_of(macro_obj) == "MACRO"
        @test m.type_of(len) == "BUILTIN"
        @test m.type_of(return_value) == "RETURN_VALUE"
        @test m.type_of(quote_obj) == "QUOTE"
        @test m.type_of(compiled_function_obj) == "COMPILED_FUNCTION"
        @test m.type_of(closure_obj) == "CLOSURE"
    end

    @testset "Stringify" begin
        @test string(int_obj) == "123"
        @test string(true_obj) == "true"
        @test string(false_obj) == "false"
        @test string(null_obj) == "null"
        @test string(str_obj) == "\"str\""
        @test string(error_obj) == "ERROR: error"
        @test string(arr_obj) == "[1, 2]"
        @test string(hash_obj) ∈ ["{1:2, 3:4}", "{3:4, 1:2}"]
        @test string(function_obj) == "fn(x) {\n(x + 2)\n}"
        @test string(macro_obj) == "macro(x) {\n\n}"
        @test string(len) == "builtin function"
        @test string(return_value) == "true"
        @test string(quote_obj) == "QUOTE(null)"
        @test string(compiled_function_obj) == "compiled function"
        @test string(closure_obj) == "closure"
    end

    @testset "Equality" begin
        @test int_obj == m.IntegerObj(123)
        @test int_obj != m.IntegerObj(122)
        @test str_obj == m.StringObj("str")
        @test str_obj != m.StringObj("srt")
        @test true_obj == m.BooleanObj(true)
        @test false_obj == m.BooleanObj(false)
        @test null_obj == m.NullObj()
        @test true_obj != false_obj
        @test false_obj != null_obj
        @test true_obj != null_obj
        @test arr_obj == m.ArrayObj([m.IntegerObj(1), m.IntegerObj(2)])
        @test arr_obj != m.ArrayObj([m.IntegerObj(2), m.IntegerObj(1)])
        @test hash_obj == m.HashObj(
            Dict(m.IntegerObj(1) => m.IntegerObj(2), m.IntegerObj(3) => m.IntegerObj(4)),
        )
        @test hash_obj == m.HashObj(
            Dict(m.IntegerObj(3) => m.IntegerObj(4), m.IntegerObj(1) => m.IntegerObj(2)),
        )
        @test hash_obj != m.HashObj(
            Dict(m.IntegerObj(3) => m.StringObj("4"), m.IntegerObj(1) => m.IntegerObj(2)),
        )
        @test hash_obj != m.HashObj(
            Dict(m.IntegerObj(3) => m.IntegerObj(4), m.IntegerObj(1) => m._TRUE),
        )
    end
end


text:
using StructArrays, ProgressMeter, KitBase

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "sod"
    D[:space] = "1d1f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "fix"
    D[:cfl] = 0.5
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 0

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 64
    D[:vMeshType] = "rectangle"
    D[:nug] = 0

    D[:knudsen] = 0.01
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 0.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

ks = SolverSet(D)

ctr = Array{ControlVolume1D1F}(undef, ks.pSpace.nx)
for i in eachindex(ctr)
    prim = [10.0 * rand(), randn(), 1 / rand()]
    w = prim_conserve(prim, ks.gas.γ)
    g = maxwellian(ks.vs.u, prim)

    ctr[i] = ControlVolume1D1F(ks.pSpace.x[i], ks.pSpace.dx[i], w, prim, g)
end
ctr = ctr |> StructArray

face = Array{Interface1D1F}(undef, ks.pSpace.nx + 1)
for i = 1:ks.pSpace.nx+1
    face[i] = Interface1D1F(deepcopy(ks.ib.wL), deepcopy(ks.ib.fL)) # necessary
end
face = face |> StructArray

plot_line(ks, ctr)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ks.ib.wL)

@showprogress for iter = 1:100#nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))
    KitBase.update!(
        ks,
        ctr,
        face,
        dt,
        res;
        coll = Symbol(ks.set.collision),
        bc = Symbol(ks.set.boundary),
    )

    t += dt
    #if t > ks.set.maxTime || maximum(res) < 5.e-7
    #    break
    #end
end

plot_line(ks, ctr)


text:
module btest

type BackTraceElement
    mod::Symbol
    func::Symbol
    file::Symbol
    line::Int
end

function getbacktrace()
    bt = backtrace()
    btout = BackTraceElement[]
    mod = current_module()
    for b in bt
        code = ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Cint), b, true)
        if !code[6]
            push!(btout, BackTraceElement(symbol(mod),code[1],code[2],code[3]))
        end
    end
    return btout
end

function test1()
    println("Calling test 1")
    test2("Calling test 2")
    a = pwd()
    return joinpath(a,"a")
end

function test2(a)
    println("$a")
    bts = getbacktrace()
    for bt in bts
        println(bt)
    end
end

test1()
end


text:
using BinDeps
using Compat

import BinDeps: bindir, downloadsdir

@BinDeps.setup
chromium = library_dependency("chromium")

# Chromium revision to download
const revision = 511134

const googlehost = "https://storage.googleapis.com/chromium-browser-snapshots/"

const DOWNLOADS = Dict(
    "Linux-x86_64"   => ("Linux_x64/$revision/chrome-linux.zip", "chrome-linux/chrome"),
    "Darwin-x86_64"  => ("Mac/$revision/chrome-mac.zip",         "chrome-mac/Chromium.app/Contents/MacOS/Chromium"),
    "Windows-x86_64" => ("Win_x64/$revision/chrome-win32.zip",   "chrome-win32/chrome.exe"),
    "Windows-i686"   => ("Win/$revision/chrome-win32.zip",       "chrome-win32/chrome.exe")
)

const SYSTEM = string(BinDeps.OSNAME, '-', Sys.ARCH)

if haskey(DOWNLOADS, SYSTEM)
    suburl, exename = DOWNLOADS[SYSTEM]
    durl = googlehost * suburl
else
    error("No precompiled binaries found for your system. Sorry...")
end

dwnlfile = joinpath(downloadsdir(chromium), "chromium.zip")
unzipdir = bindir(chromium)

run(@build_steps begin
  CreateDirectory(downloadsdir(chromium))
  FileDownloader(durl, dwnlfile)
  # CreateDirectory(bindir(chromium))
  FileUnpacker(dwnlfile, unzipdir, ".")
end)

### write deps.jl file
exepath = joinpath(dirname(@__FILE__),"usr/bin",exename)
open(joinpath(dirname(@__FILE__), "deps.jl"), "w") do io
    print(io, """
        const chromium = "$(escape_string(exepath))"
    """)
end


text:
#!/usr/bin/env julia

const N = 440_000_000
const cache = (0, ntuple(i -> i^i, 9)...)

function is_munchausen(number)
    n = number
    total = 0
    while n > 0
        (n, digit) = divrem(n, 10)
        @inbounds total += cache[digit+1]
        if total > number
            return false
        end
    end

    return total == number
end

function main()
    @sync for chunk ∈ Iterators.partition(0:N, N÷(4*Threads.nthreads()))
        Threads.@spawn for n ∈ chunk
            if is_munchausen(n)
                println(n)
            end
        end
    end
end

main()

max_stars_repo_name: uninhm/SpeedTests

max_issues_repo_name: uninhm/SpeedTests

max_forks_repo_name: uninhm/SpeedTests


text:
using JLD2, Plots, Printf, LinearAlgebra, Statistics
using LaTeXStrings
include(pwd() * "/analysis_scripts/" * "post_analysis.jl")

filename = "/Weno_20_zonal_averages.jld2"
file = jldopen( pwd() * filename)
data_dictionary = get_data(file)
fine_z  = file["grid"]["zC"][2:end-1]
fine_y  = file["grid"]["yC"][2:end-1]
fine_yF = file["grid"]["yF"][2:end-1]
fine_zF = file["grid"]["zF"][2:end-1]
fine_b = mean(data_dictionary["b"][2:end], dims = 1)[1]
fine_u = mean(data_dictionary["u"][2:end], dims = 1)[1]
fine_v = mean(data_dictionary["v"][2:end], dims = 1)[1]
fine_v = (fine_v[2:end,:] + fine_v[1:end-1,:]) ./ 2
data_dictionary = get_data(file)
close(file)

filename = "/Weno_CPU_20_zonal_averages.jld2"
file = jldopen( pwd() * filename)
data_dictionary = get_data(file)
coarse_z  = file["grid"]["zC"][2:end-1]
coarse_y  = file["grid"]["yC"][2:end-1]
coarse_yF = file["grid"]["yF"][2:end-1]
coarse_zF = file["grid"]["zF"][2:end-1]
coarse_b = mean(data_dictionary["b"][2:end], dims = 1)[1]
coarse_u = mean(data_dictionary["u"][2:end], dims = 1)[1]
coarse_v = mean(data_dictionary["v"][2:end], dims = 1)[1]
coarse_v = (coarse_v[2:end,:] + coarse_v[1:end-1,:]) ./ 2
close(file)

filename = "/Weno_CPU_Hack_20_zonal_averages.jld2"
file = jldopen( pwd() * filename)
data_dictionary = get_data(file)
hack_z  = file["grid"]["zC"][2:end-1]
hack_y  = file["grid"]["yC"][2:end-1]
hack_yF = file["grid"]["yF"][2:end-1]
hack_zF = file["grid"]["zF"][2:end-1]
hack_b = mean(data_dictionary["b"][50:end], dims = 1)[1]
hack_u = mean(data_dictionary["u"][50:end], dims = 1)[1]
hack_v = mean(data_dictionary["v"][50:end], dims = 1)[1]
hack_v = (hack_v[2:end,:] + hack_v[1:end-1,:]) ./ 2
close(file)

function plot_field(a, y, yF, z, zF; name = " ", clims = nothing)
    py = get_y(a, y, yF) # global scope
    pz = get_z(a, z, zF) # global scope
    if clims == nothing
        clims = (minimum(a), maximum(a))
    end
    further_label = "Zonal and Temporal Average of "
    p1 = contourf(py, pz, a',
    color = :thermometer, title = " " * name,
    xlabel = "South to North [m]", ylabel = "Depth [m]"
    , clims = clims, linewidth = 0)
    display(p1)
    return p1
end

coarse_plot = plot_field(coarse_b, coarse_y, coarse_yF, coarse_z, coarse_zF, name = "Coarse Buoyancy")
fine_plot = plot_field(fine_b, fine_y, fine_yF, fine_z, fine_zF, name = "Coarse Buoyancy")

function avgy(Φ, n)
    m = size(Φ)[1]
    scale = Int(floor(m/n))
    if ( abs(Int(floor(m/n)) - m/n) > eps(1.0))
        return error
    end
    ny, nz = size(Φ)
    Φ2 = zeros(n, nz)
    for i in 1:n
        Φ2[i,:] .= 0
            for j in 1:scale
                Φ2[i,:] .+= Φ[scale*(i-1) + j,:] / scale
            end
    end
    return Φ2
end

coarse_grained_b = avgy(fine_b, 12)
coarse_grained_plot = plot_field(coarse_grained_b, coarse_y, coarse_yF, coarse_z, coarse_zF, name = "Coarse Grained Buoyancy")

norm(coarse_grained_b - coarse_b) / norm(coarse_grained_b)
norm(coarse_grained_b - hack_b) / norm(coarse_grained_b)

residual_plot = plot_field(coarse_b-coarse_grained_b , coarse_y, coarse_yF, coarse_z, coarse_zF, name = "Slack Buoyancy")

## Zonal Velocity
coarse_plot = plot_field(coarse_v, coarse_y, coarse_yF, coarse_z, coarse_zF, name = "Coarse Meridional Velocity")
fine_plot = plot_field(fine_v, fine_y, fine_yF, fine_z, fine_zF, name = "Meridional Velocity")
coarse_grained_v = avgy(fine_v, 12)
clims = (minimum(coarse_v), maximum(coarse_v))
coarse_grained_plot = plot_field(coarse_grained_v - coarse_v, coarse_y, coarse_yF,
                                 coarse_z, coarse_zF,
                                 name = "Difference in Meridional Velocity")

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale