text:
@enum(LintCodes,IncorrectCallNargs,IncorrectIterSpec,NothingEquality,NothingNotEq)

function _typeof(x, state)
    if typof(x) in (CSTParser.Abstract, CSTParser.Primitive, CSTParser.Struct, CSTParser.Mutable)
        return getsymbolserver(state.server)["Core"].vals["DataType"]
    elseif typof(x) in (CSTParser.ModuleH, CSTParser.BareModule)
        return getsymbolserver(state.server)["Core"].vals["Module"]
    elseif typof(x) in (CSTParser.ModuleH, CSTParser.BareModule)
        return getsymbolserver["Core"].vals["Module"]
    elseif CSTParser.defines_function(x)
        return return getsymbolserver(state.server)["Core"].vals["Function"]
    end
end

function checks(x, server)
    # check_call_args(x)
    check_loop_iter(x, server)
    check_nothing_equality(x, server)
end

function check_call_args(x::EXPR)
    if typof(x) === Call
        if typof(x.args[1]) === IDENTIFIER && hasref(x.args[1])
            func = refof(x.args[1])
        else
            func = nothing
        end
        if func isa Binding
        elseif func isa SymbolServer.FunctionStore
            nargs = _get_call_nargs(x)
            ok = false
            for m in func.methods
                if nargs == length(m.args)
                    ok = true
                    break
                elseif nargs > length(m.args) && length(m.args) > 0 && endswith(last(m.args)[2], "...")
                    ok = true
                    break
                end
            end
            if !ok
                setref!(x, IncorrectCallNargs)
            end
        end
    end
end

function check_loop_iter(x::EXPR, server)
    if typof(x) === CSTParser.For
        if length(x.args) > 1 && typof(x.args[2]) === CSTParser.BinaryOpCall && refof(x.args[2]) === nothing
            rng = x.args[2].args[3]
            if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                setref!(x.args[2], IncorrectIterSpec)
            elseif typof(rng) === CSTParser.Call && refof(rng.args[1]) === getsymbolserver(server)["Base"].vals["length"]
                setref!(x.args[2], IncorrectIterSpec)
            end
        end
    elseif typof(x) === CSTParser.Generator
        for i = 3:length(x.args)
            if typof(x.args[i]) === CSTParser.BinaryOpCall && refof(x.args[i]) === nothing
                rng = x.args[i].args[3]
                if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                    setref!(x.args[i], IncorrectIterSpec)
                elseif typof(rng) === CSTParser.Call && valof(rng.args[1]) == "length" && refof(rng.args[1]) === getsymbolserver(server)["Base"].vals["length"]
                    setref!(x.args[i], IncorrectIterSpec)
                end
            end
        end
    end
end

function check_nothing_equality(x::EXPR, server)
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x.args[2]) === CSTParser.Tokens.EQEQ && valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbolserver(server)["Core"].vals["nothing"]
            setref!(x.args[2], NothingEquality)
        elseif kindof(x.args[2]) === CSTParser.Tokens.NOT_EQ && valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbolserver(server)["Core"].vals["nothing"]
            setref!(x.args[2], NothingNotEq)
        end
    end
end

function _get_call_nargs(x::EXPR)
    cnt = 0
    for i = 2:length(x.args)
        if typof(x.args[i]) === PUNCTUATION
        elseif typof(x.args[i]) === CSTParser.Parameters
            for j = 1:length(x.args[i].args)
                if typof(x.args[i].args[j]) !== PUNCTUATION
                    cnt += 1
                end
            end
        else
            cnt += 1
        end
    end
    return cnt
end

function _get_top_binding(x::EXPR, name::String)
    if scopeof(x) isa Scope
        return _get_top_binding(scopeof(x))
    elseif parentof(x) isa EXPR
        return _get_top_binding(parentof(x), name)
    else
        return nothing
    end
end

function _get_top_binding(s::Scope, name::String)
    if haskey(s.names, name)
        return s.names[name]
    elseif parentof(s) isa Scope
        return _get_top_binding(parentof(s), name)
    else
        return nothing
    end
end

function _get_global_scope(s::Scope)
    if !s.ismodule && parentof(s) isa Scope && parentof(s) != s
        return _get_global_scope(parentof(s))
    else
        return s
    end
end


text:
using UnPack
export @unpack_val, @unpack_units

"""
    This unpacks the numerical value of `GeoUnit` in a structure, without the units.
    All requested variables must be GeoUnits.

    This is a modification of the `@unpack` macro as implemented in the `UnPack.jl` package, which can be used to retrieve the full variables.

# Example
```julia
julia> struct Density{T}
        ρ::GeoUnit{T}
        α::GeoUnit{T}
       end
julia> r = Density(GeoUnit(100kg/m^3),GeoUnit(4e-5/K));
julia> @unpack_val ρ,α = r
julia> α
4.0e-5
julia> typeof(α)
Float64
```
"""
macro unpack_val(args)
    args.head!=:(=) && error("Expression needs to be of form `a, b = c`")
    items, suitecase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()

    # This extracts the value, but not the units
    kd = [:( $key = $UnPack.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()).val
            ) for key in items]

    kdblock      = Expr(:block, kd...)

    expr = quote
        local $suitecase_instance = $suitecase # handles if suitecase is not a variable but an expression
        $kdblock
        $suitecase_instance # return RHS of `=` as standard in Julia
    end
    esc(expr)
end

"""
    This unpacks the numerical value with units of `GeoUnit` parameters in a structure
    All requested variables must be GeoUnits.

    This is a modification of the `@unpack` macro as implemented in the `UnPack.jl` package, which can be used to retrieve the full variables.

# Example
```julia
julia> struct Density{T}
        ρ::GeoUnit{T}
        α::GeoUnit{T}
       end
julia> r = Density(GeoUnit(100kg/m^3),GeoUnit(4e-5/K));
julia> @unpack_units ρ,α = r
julia> α
4.0e-5 K⁻¹·⁰
julia> typeof(α)
Quantity{Float64, 𝚯⁻¹·⁰, Unitful.FreeUnits{(K⁻¹·⁰,), 𝚯⁻¹·⁰, nothing}}
```
"""
macro unpack_units(args)
    args.head!=:(=) && error("Expression needs to be of form `a, b = c`")
    items, suitecase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()

    # This extracts the value with units
    kd = [:( $key =   $UnPack.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()).val
                   .* $UnPack.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()).unit
            ) for key in items]

    kdblock      = Expr(:block, kd...)

    expr = quote
        local $suitecase_instance = $suitecase # handles if suitecase is not a variable but an expression
        $kdblock
        $suitecase_instance # return RHS of `=` as standard in Julia
    end
    esc(expr)
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e5013406-485e-11eb-1923-bf36ce0e5db9
using LightGraphs

# ╔═╡ 010aa6a0-485f-11eb-0fef-3963df8e593e
using GraphRecipes

# ╔═╡ 2c0483c6-485f-11eb-19a0-07f4e49962ff
using PGFPlotsX

# ╔═╡ 32dae7ee-485f-11eb-184a-85e8e38c594c
using Plots

# ╔═╡ f785d38e-485e-11eb-3937-13c79e021ccd
import Pkg

# ╔═╡ 5d70213c-485e-11eb-3141-27d44a424b79
md" ## _Types_
In this lesson we will learn what types are and how it is possible to define functions that work on types. We will learn which are the differences between abstract and concrete types, how to define immutable and mutable types and how to create a type constructor. We will give a brief introduction to multiple dispatch and see how types have a role in it.

We can think of types as containers for data only. Moreover, it is possible to define a type hierarchy so that functions that work for parent type work also for the children (if they are written properly). A parent type can only be an AbstractType (like Number), while a child can be both an abstract or concrete type.
"

# ╔═╡ 99ddb7ea-485e-11eb-153a-0511e9047049
begin
g = DiGraph(13)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
add_edge!(g, 3, 4)
add_edge!(g, 3, 7)
# Instead of adding an edge 5 -> 3, have an edge that goes from 5 to
# the ghost node, then to 3.
add_edge!(g, 4, 5)
add_edge!(g, 4, 6)
add_edge!(g, 7, 8)
add_edge!(g, 7, 9)
add_edge!(g, 8, 10)
add_edge!(g, 8, 11)
add_edge!(g, 9, 12)
add_edge!(g, 9, 13)

# Add an edge by itself that will force GraphRecipes to have a decent
# window size.
add_edge!(g, 7, 7)
add_edge!(g, 8, 8)
add_edge!(g, 9, 9)
add_edge!(g, 10, 10)
add_edge!(g, 11, 11)
add_edge!(g, 12, 12)
add_edge!(g, 13, 13)

names = [:Any, :Number, :Real, :abstractfloat, :float64, :float32, :Integer, :signed, :unsigned, :Int64, :Int32, :UInt64, :UInt32]
x = [0,0,0,0,-1,1,5,4,6,3.25,4.75,6,8]
y = [5,4,3,2,1,1,3,2,2,1,1,1,1]
graphplot(g, x=x, y=y, curvature_scalar=0.0, nodesize=1.0,
          node_weights=[1.0,1.0,1.0,1.0,1.0,1.0,1.0 ,1.0,1.0,1.0,1.0,1.0,1.0],
          names=names, nodecolor=:lightyellow, color=:black,
          nodeshape=:rect, self_edge_size=0.0)
end

# ╔═╡ 58d29492-485f-11eb-3feb-83d5046d9925
md" ## Implementation

To declare a Type we use either the type or struct keyword.

To declare an abstract type we use:"

# ╔═╡ 89001274-4860-11eb-115f-8ff6a4faa465
begin
	abstract type Person
	end

	abstract type Musician <: Person
	end
end

# ╔═╡ 91937b60-4860-11eb-12d2-592e8eeb7768
md"You may find it surprising, but apparently musicians are people, so Musician is a sub-type of Person. There are many kind of musicians, for example rock-stars and classic musicians, so we define two new concrete types (in particular this kind of type is called a composite type):"

# ╔═╡ 9bd7a6aa-4860-11eb-1c4b-c3be7d059e26
begin
	mutable struct Rockstar <: Musician
		name::String
		instrument::String
		bandName::String
		headbandColor::String
		instrumentsPlayed::Int
	end

	struct ClassicMusician <: Musician
		name::String
		instrument::String
	end
end

# ╔═╡ a3a8b202-4860-11eb-390e-95f6823e90aa
md"Notably rock-stars love to change the colour of their headband, so we have made Rockstar a mutable struct, which is a concrete type whose elements value can be modified. On the contrary, classic musicians are known for their everlasting love for their instrument, which will never change, so we have made ClassicMusician an immutable concrete type.

We can define another sub-type of Person, Physicist, as I am a physicist and I was getting envious of rock-stars:"

# ╔═╡ c7a5a778-4860-11eb-187e-bdf191565e2b
begin
	begin
		mutable struct Physicist <: Person
			name::String
			sleepHours::Float64
			favouriteLanguage::String
		end

		aure = Physicist("Aurelio", 6, "Julia")
	end
end

# ╔═╡ d3528442-4860-11eb-389b-c300063289aa
aure.name

# ╔═╡ fe554d64-4860-11eb-0f32-a178f0f21ff9
aure.sleepHours

# ╔═╡ 04a67940-4861-11eb-323d-49e8563f4a29
aure.favouriteLanguage

# ╔═╡ 111c30b6-4861-11eb-2281-0773165d8ecb
md"Luckily my exam session is over now and I finally have a little bit more time to sleep, so I’ll adjust my sleeping schedule to sleep eight hours:"

# ╔═╡ 1a4db196-4861-11eb-02f7-670caaeb9562
aure.sleepHours = 8

# ╔═╡ 21b56708-4861-11eb-1253-8f00dba8057c
aure.sleepHours

# ╔═╡ 2845cbd0-4861-11eb-1809-5b08f407c467
md"Incidentally I am also a ClassicMusician and I play violin, so I can create a new structure:"

# ╔═╡ 31e88ce0-4861-11eb-2c03-5976cdc914c7
aure_musician = ClassicMusician("Aurelio", "Violin")

# ╔═╡ 3ba8006c-4861-11eb-2da3-19bd73a95bbf
aure_musician.instrument = "Cello"

# ╔═╡ 46b194e6-4861-11eb-0c5d-39996e11e48e
md"As you can see, I love violin and I just can’t change my instrument, as ClassicMusician is an immutable struct.

I am not a rock-star, but my friend Ricky is one, so we’ll define:"

# ╔═╡ 5995feda-4861-11eb-12f5-bb97dd2330c8
ricky = Rockstar("Riccardo", "Voice", "Black Lotus", "red", 2)

# ╔═╡ 662bb806-4861-11eb-2850-39fb55ddbb6b
ricky.headbandColor

# ╔═╡ 68ff6be8-4861-11eb-10f1-1dda03a33073
md" ## Functions and Types: Multiple Dispatch
It is possible to write functions that operate on both abstract and concrete types. For example, every person is likely to have a name, so we can define the following function:"

# ╔═╡ bd6882e2-4862-11eb-3a7f-794c357d326c
function introduceMe(person::Person)
    println("Hello, my name is $(person.name).")
end

# ╔═╡ d7f9560e-4862-11eb-03ec-99bb75e49bf6
md"While only musicians play instruments, so we can define the following function:"

# ╔═╡ 40e20a1c-4863-11eb-2354-e9da49a7459a
function introduceMe(person::Musician)
    println("Hello, my name is $(person.name) and I play $(person.instrument).")
end

# ╔═╡ 623b7180-4863-11eb-32fa-4b749af28529
md" and for rockstar we can write:"

# ╔═╡ 77c6bf6e-4863-11eb-3354-55ba248b5c50
function introduceMe(person::Rockstar)
	if person.instrument == "Voice"
		println("Hello, my name is $(person.name) and I sing.")
	else
		println("Hello, my name is $(person.name) and I play $(person.instrument).")
	end

	println("My band name is $(person.bandName) and my favourite headband colour is $(person.headbandColor)!")
end

# ╔═╡ d012c286-4862-11eb-35d3-77692b72a1be
introduceMe(aure)

# ╔═╡ 4a3d5c60-4863-11eb-399f-617db3848bfe
introduceMe(aure_musician)

# ╔═╡ 7b0a5d54-4863-11eb-3db4-d3a04521f8c3
md"The ::SomeType notation indicates to Julia that person has to be of the aforementioned type or a sub-type. Only the most strict type requirement is considered (which is the lowest type in the type tree), for example ricky is a Person, but “more importantly” he is a Rockstar (Rockstar is placed lower in the type tree), thus introduceMe(person::Rockstar) is called. In other words, the function with the closest type signature will be called.

This is an example of multiple dispatch, which means that we have written a single function with different methods depending on the type of the variable. We will come back again to multiple dispatch in this lesson, as it is one of the most important features of Julia and is considered a more advanced topic, together with type annotations. As an anticipation ::Rockstar is a type annotation, the compiler will check if person is a Rockstar (or a sub-type of it) and if that is true it will execute the function.

 ## Type constructor

When a type is applied like a function it is called a constructor. When we created the previous types, two constructors were generated automatically (these are called default constructors). One accepts any arguments and calls convert to convert them to the types of the fields, and the other accepts arguments that match the field types exactly (String and String in the case of ClassicMusician). The reason both of these are generated is that this makes it easier to add new definitions without inadvertently replacing a default constructor.

Sometimes it is more convenient to create custom constructor, so that it is possible to assign default values to certain variables, or perform some initial computations."

# ╔═╡ 972fd872-4863-11eb-3b18-2b46a7bf85fb
mutable struct MyData
	x::Float64
	x2::Float64
	y::Float64
	z::Float64
	function MyData(x::Float64, y::Float64)
		x2=x^2
		z = sin(x2+y)
		new(x, x2, y, z)
	end
end

# ╔═╡ a9c2d0ac-4863-11eb-1640-e93121a06e98
MyData(2.0, 3.0)

# ╔═╡ b0a8f75a-4863-11eb-01e2-db6da2c6e3ee
md"Sometimes it may be useful to use other types for x, x2 and y, so it is possible to use parametric types (i.e. types that accept type information at construction time):"

# ╔═╡ ba37472e-4863-11eb-2f79-535974cc8e37
mutable struct MyData2{T<:Real}
	x::T
	x2::T
	y::T
	z::Float64
	function MyData2{T}(x::T, y::T) where {T<:Real}
		x2=x^2
		z = sin(x2+y)
		new(x, x2, y, z)
	end
end

# ╔═╡ c2660f98-4863-11eb-1133-1bfd03b8ca06
MyData2{Float64}(2.0,3.0)

# ╔═╡ c8f4ec08-4863-11eb-2301-f774e405a45d
MyData2{Int}(2,3)

# ╔═╡ db68e9ca-4863-11eb-2e46-975a727a08cf
md" ## Example
Mutable types are particularly useful when it comes to storing data that needs to be shared between some functions inside a module. It is not uncommon to define custom types in a module to store all the data which needs to be shared between functions and which is not constant."

# ╔═╡ f5949efc-4863-11eb-22bf-3f5547ce8624
circle1 = Circle{Float64}(5.0)

# ╔═╡ f5abb4d4-4863-11eb-1545-61a62964aa72
computePerimeter(circle1)

# ╔═╡ f5bf7d48-4863-11eb-1dc7-055b9a8c1bfe
circle1.perimeter

# ╔═╡ f5d0a514-4863-11eb-35e4-d7dabb0ae0b4
computeArea(circle1)

# ╔═╡ f5e298d2-4863-11eb-2178-53d4d230ed42
circle1.area

# ╔═╡ f6119e52-4863-11eb-1e46-a90051f25108
printCircleEquation(2, 3, circle1)

# ╔═╡ 29b3b0bc-4864-11eb-086f-e59856887cb7
md"This is a simple module which implements a Circle type which contains the radius, perimeter and area of the circle. There are three functions which respectively compute the perimeter and area of the circle and store them inside theCircle structure. The third function prints the equation of a circle with a given centre and the radius stored inside a Circle structure.

Notice that we could have simply computed the perimeter and area inside the type constructor, but I have chosen not to do so for educative purposes.

## Conclusions

This lesson has been a little bit more conceptually difficult than the previous ones, but you don’t need to remember everything right now! We will use types in the future lessons, so you will naturally get accustomed to how they works over time.

We have learnt how to define abstract and concrete types, and how to define mutable and immutable structures. We have then learnt how it is possible to define functions that work on custom types and we have introduced multiple dispatch. Furthermore, we have seen how to define an inner constructor, to aid the user create an instance of a composite type. Lastly, we have seen an example of a module which uses a custom type (Circle) to perform and store some specific computations."

# ╔═╡ f58a0762-4863-11eb-3f8a-fb10c45ee27f
#%%

using .TestModuleTypes1

# ╔═╡ e910d4a2-4863-11eb-3bd0-eb35783028b9
module TestModuleTypes1

export Circle, computePerimeter, computeArea, printCircleEquation

mutable struct Circle{T<:Real}
	radius::T
	perimeter::Float64
	area::Float64

	function Circle{T}(radius::T) where T<:Real
		# we initialize perimeter and area to -1.0, which is not a possible value
		new(radius, -1.0, -1.0)
	end
end

@doc raw"""
	computePerimeter(circle::Circle)

Compute the perimeter of `circle` and store the value.
"""
function computePerimeter(circle::Circle)
	circle.perimeter = 2*π*circle.radius
	return circle.perimeter
end

@doc raw"""
	computeArea(circle::Circle)

Compute the area of `circle` and store the value.
"""
function computeArea(circle::Circle)
	circle.area = π*circle.radius^2
	return circle.area
end

@doc raw"""
	printCircleEquation(xc::Real, yc::Real, circle::Circle )

Print the equation of a cricle with center at (xc, yc) and radius given by circle.
"""
function printCircleEquation(xc::Real, yc::Real, circle::Circle )
	println("(x - $xc)^2 + (y - $yc)^2 = $(circle.radius^2)")
	return
end
end # end module

# ╔═╡ Cell order:
# ╠═e5013406-485e-11eb-1923-bf36ce0e5db9
# ╠═f785d38e-485e-11eb-3937-13c79e021ccd
# ╠═010aa6a0-485f-11eb-0fef-3963df8e593e
# ╠═2c0483c6-485f-11eb-19a0-07f4e49962ff
# ╠═32dae7ee-485f-11eb-184a-85e8e38c594c
# ╠═5d70213c-485e-11eb-3141-27d44a424b79
# ╟─99ddb7ea-485e-11eb-153a-0511e9047049
# ╠═58d29492-485f-11eb-3feb-83d5046d9925
# ╠═89001274-4860-11eb-115f-8ff6a4faa465
# ╠═91937b60-4860-11eb-12d2-592e8eeb7768
# ╠═9bd7a6aa-4860-11eb-1c4b-c3be7d059e26
# ╠═a3a8b202-4860-11eb-390e-95f6823e90aa
# ╠═c7a5a778-4860-11eb-187e-bdf191565e2b
# ╠═d3528442-4860-11eb-389b-c300063289aa
# ╠═fe554d64-4860-11eb-0f32-a178f0f21ff9
# ╠═04a67940-4861-11eb-323d-49e8563f4a29
# ╠═111c30b6-4861-11eb-2281-0773165d8ecb
# ╠═1a4db196-4861-11eb-02f7-670caaeb9562
# ╠═21b56708-4861-11eb-1253-8f00dba8057c
# ╠═2845cbd0-4861-11eb-1809-5b08f407c467
# ╠═31e88ce0-4861-11eb-2c03-5976cdc914c7
# ╠═3ba8006c-4861-11eb-2da3-19bd73a95bbf
# ╠═46b194e6-4861-11eb-0c5d-39996e11e48e
# ╠═5995feda-4861-11eb-12f5-bb97dd2330c8
# ╠═662bb806-4861-11eb-2850-39fb55ddbb6b
# ╠═68ff6be8-4861-11eb-10f1-1dda03a33073
# ╠═bd6882e2-4862-11eb-3a7f-794c357d326c
# ╠═d012c286-4862-11eb-35d3-77692b72a1be
# ╠═d7f9560e-4862-11eb-03ec-99bb75e49bf6
# ╠═40e20a1c-4863-11eb-2354-e9da49a7459a
# ╠═4a3d5c60-4863-11eb-399f-617db3848bfe
# ╠═623b7180-4863-11eb-32fa-4b749af28529
# ╠═77c6bf6e-4863-11eb-3354-55ba248b5c50
# ╠═7b0a5d54-4863-11eb-3db4-d3a04521f8c3
# ╠═972fd872-4863-11eb-3b18-2b46a7bf85fb
# ╠═a9c2d0ac-4863-11eb-1640-e93121a06e98
# ╠═b0a8f75a-4863-11eb-01e2-db6da2c6e3ee
# ╠═ba37472e-4863-11eb-2f79-535974cc8e37
# ╠═c2660f98-4863-11eb-1133-1bfd03b8ca06
# ╠═c8f4ec08-4863-11eb-2301-f774e405a45d
# ╠═db68e9ca-4863-11eb-2e46-975a727a08cf
# ╠═e910d4a2-4863-11eb-3bd0-eb35783028b9
# ╠═f58a0762-4863-11eb-3f8a-fb10c45ee27f
# ╠═f5949efc-4863-11eb-22bf-3f5547ce8624
# ╠═f5abb4d4-4863-11eb-1545-61a62964aa72
# ╠═f5bf7d48-4863-11eb-1dc7-055b9a8c1bfe
# ╠═f5d0a514-4863-11eb-35e4-d7dabb0ae0b4
# ╠═f5e298d2-4863-11eb-2178-53d4d230ed42
# ╠═f6119e52-4863-11eb-1e46-a90051f25108
# ╠═29b3b0bc-4864-11eb-086f-e59856887cb7


text:
using MyPoly
p1=MyPoint(0,0)
p2=MyPoint(1,0)
p3=MyPoint(1,1)
p4=MyPoint(0,1)
p5=MyPoint(0,0)
poly1=MyPolyVar([p1,p2,p3,p4,p5],10)
a=area(poly1)
println(a)


text:
"""
    lastClick(X::AbstractArray; lookback = 60)

Generate last click attribution based on a lookback window. Default is 60.
"""
function lastclick(X::AbstractArray{T}) where T
    return X
end


text:
export acceptance, describe

# Compute acceptance or rejection rate of MCMCChain
function acceptance(c::MCMCChain; lags::Ranges=1:size(c.samples, 1), reject::Bool=false)
  rlen = length(lags)
  @assert lags[end] <= size(c.samples, 1) "Range of acceptance rate not within post-burnin range of MCMC chain"

  if reject
    return (rlen-sum(c.diagnostics["accept"][lags]))*100/rlen
  else
    return sum(c.diagnostics["accept"][lags])*100/rlen
  end
end

# TODO 1: Compute MCMC quantiles based on
# Flegal J.M, Galin L.J, Neath R.C. Markov Chain Monte Carlo Estimation of Quantiles. arXiv, 2013
# TODO 2: Include these MCMC estimates of quantiles in describe()

# describe() provides summary statistics for MCMCChain objects
describe(c::MCMCChain) = describe(STDOUT, c)

function describe(io::IO, c::MCMCChain)
  nsamples, npars = size(c.samples)
  if isa(c.task.sampler, ARS)
    indx = find(c.diagnostics["accept"])
  else
    indx = 1:nsamples
  end
  for i in 1:npars

    col = c.samples[indx, i]
    println(io, "Parameter $i")

    if sum(isnan(col)) != 0
      println(io, "Monte Carlo chain for parameter $i contains NaNs, which are not supported.")
      return
    end

    qs = quantile(col, [0, 1])
    varimse = mcvar_imse(col)
    variid = var(col)/length(indx)
    mcerror = sqrt(varimse)
    ss = length(indx)*variid./varimse
    act = varimse./variid

    statNames = ["Min", "Mean", "Max", "MC Error", "ESS", "AC Time"]
    statVals = [qs[1], mean(col), qs[2], mcerror, ss, act]
    for i = 1:6
        println(io, string(rpad(statNames[i], 10, " "), " ", string(statVals[i])))
    end

    println(io)
  end
end


text:
using DataFrames: DataFrame
using Iex

TEST_SLEEP_TIME = 15.0

profile_functions_test = [
    :company,
    :insider_roster,
    :insider_summary,
    :logo,
    :peers,
    :balance_sheet,
    :earnings,
    :income,
    :news,
]

params = (
    token="Tsk_37cebd73a857428f9e27b4cd3aa1bc4a",
    url="https://sandbox.iexapis.com/",
    version="stable",
)
iex = Connection(;url=params.url, version=params.version,
    token=params.token);

@testset "StockProfiles" begin
    for f in profile_functions_test
        @eval begin
            testname = string($f)
            @testset "Testing $testname" begin
                data = $f(iex, "AAPL")
                if testname == "company"
                    @test typeof(data) === Dict{String,Any}
                    @test length(data) === 20
                elseif testname in ["insider_roster", "insider_summary", "news"]
                    @test typeof(data) === Array{Any,1}
                    @test length(data) === 10
                elseif testname == "logo"
                    @test typeof(data) == Dict{String,Any}
                    @test length(data) == 1
                elseif testname == "peers"
                    @test typeof(data) === Array{Any,1}
                    @test length(data) === 4
                elseif testname in ["balance_sheet", "earnings", "income"]
                    @test typeof(data) == Dict{String,Any}
                    @test length(data) == 2
                elseif testname == "news"
                    @test typeof(data) == Array{Any,1}
                    @test length(data) == 10
                end
            end
        end
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
# note: this example only works on pure julia implementations of
# preprocessing elements and learners. Scikitlearn and other
# PyCall functions are not thread-safe and cannot be used inside
# the threads loop.

using AutoMLPipeline
using AutoMLPipeline.Utils
using Base.Threads

begin
  profbdata = getprofb()
  X = profbdata[:,2:end]
  Y = profbdata[:,1] |> Vector;
  head(x)=first(x,5)
  head(profbdata)
end

#### Column selector
catf = CatFeatureSelector();
numf = NumFeatureSelector()
ohe = OneHotEncoder()
#### Learners
rf = RandomForest();
ada = Adaboost()
dt=PrunedTree()

accuracy(X,Y)=score(:accuracy,X,Y)

acc=[]
learners=[rf,ada,dt]
@threads for i in 1:30
  @threads for lr in learners
    println(lr.name)
    pipe=@pipeline ((catf |> ohe) +(numf )) |> lr
    m=crossvalidate(pipe,X,Y,accuracy,10,true)
    push!(acc,(m...,name=lr.name))
    println(m)
  end
end
acc


text:
@Generated 1:32 immutable Row{Names, Types <: Tuple}
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Row parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
        error(str)
    end
    N = length(Names)
    if length(Types.parameters) != N || reduce((a,b) -> a | !isa(b, DataType), false, Types.parameters)
        str = "Row parameter 2 (Types) is expected to be a Tuple{} of $N DataTypes, got $Types"
        error(str)
    end

    exprs = [:( $(Names[i])::$(Types.parameters[i]) ) for i = 1:N]
    return Expr(:block, exprs...)
end

Base.call{Names,T <: Tuple}(::Type{Row{Names,T}}, data::T) = Row{Names,T}(data...) # Constructed from a tuple (signature must match...)

@generated function Base.call{Names}(::Type{Row{Names}}, data...)
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Row parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
        return :(error($str))
    end

    if length(Names) == length(data)
        return quote
            $(Expr(:meta,:inline))
            Row{Names, $(Expr(:curly, :Tuple, data...))}(data...)
        end
    elseif length(data) == 1 && data[1] <: Tuple && length(Names) == length(data[1].parameters)  # Constructed from a tuple (with correct length)
        return :(Row{Names, $(data[1])}(data[1]...))
    else
        return :(error("Can't construct Row with $(length(Names)) columns with input $data"))
    end
end

Base.call{Names, Types, Types_new}(::Type{Row{Names,Types_new}}, r::Row{Names,Types}) = convert(Row{Names,Types_new}, r)
@generated function Base.convert{Names, Names_new, Types, Types_new <: Tuple}(::Type{Row{Names_new,Types_new}}, r::Row{Names,Types})
    if !isa(Names_new, Tuple) || eltype(Names_new) != Symbol || length(Names_new) != length(Names) || length(Names_new) != length(unique(Names_new))
        str = "Cannot convert $(length(Names)) columns to new names $(Names_new)."
        return :(error($str))
    end
    if length(Types_new.parameters) != length(Types.parameters)
        str = "Cannot convert $(length(Types.pareters)) columns to $(length(Types_new.parameters)) new types."
        return :(error($str))
    end
    exprs = [:(convert(Types_new.parameters[$j], r.($j))) for j = 1:length(Names)]
    return Expr(:call, Row{Names_new,Types_new}, exprs...)
end

@inline colnames{Names}(::Row{Names}) = Names
@inline colnames{Names, Types <: Tuple}(::Type{Row{Names,Types}}) = Names
@inline colnames{Names}(::Type{Row{Names}}) = Names
@inline colnames{R <: Row}(::Type{R}) = colnames(super(R))
@inline eltypes{Names, Types <: Tuple}(::Row{Names,Types}) = Types
@inline eltypes{Names, Types <: Tuple}(::Type{Row{Names,Types}}) = Types
@inline eltypes{R <: Row}(::Type{R}) = eltypes(super(R))

@inline nrow(::Row) = 1
@inline ncol{Names}(::Row{Names}) = length(Names)

Base.start(r::Row) = false
Base.next(r::Row, state) = (r, true)
Base.done(r::Row, state) = state

getindex(r::Row) = r
getindex(r::Row, i) = i == 1 ? r : error("Cannot index Row at $i")

# reordering
@generated function permutecols{Names1,Names2,Types}(r::Row{Names1,Types}, ::Type{Val{Names2}})
    if Names1 == Names2
        return :(r)
    else
        if !(isa(Names2, Tuple)) || eltype(Names2) != Symbol || length(Names2) != length(Names1) || length(Names2) != length(unique(Names2))
            str = "New column names $Names2 do not match existing names $Names1"
            return :(error($str))
        end

        order = permutator(Names1, Names2)

        exprs = [:(r.($(order[j]))) for j = 1:length(Names1)]
        return Expr(:call, Row{Names1}, exprs...)
    end
end

#getindex(r::Row) = r
#@generated function getindex{Names}(r::Row{Names})
#    exprs = [:(r.($j)) for j = 1:length(Names)]
#    return Expr(:tuple, exprs...)
#end

# Horizontally concatenate cells and rows into rows
@generated Base.hcat{Name}(c::Cell{Name}) = :(Row{$((Name,))}(c.(1)))
Base.hcat(r::Row) = r

@generated function Base.hcat(r1::Union{Cell,Row}, r2::Union{Cell,Row})
    names1 = (r1 <: Cell ? (colname(r1),) : colnames(r1))
    names2 = (r2 <: Cell ? (colname(r2),) : colnames(r2))

    if length(intersect(names1, names2)) != 0
        str = "Column names are not distinct. Got $names1 and $names2"
        return :(error($str))
    end

    newnames = (names1..., names2...)
    exprs = vcat([:(r1.($j)) for j = 1:length(names1)], [:(r2.($j)) for j = 1:length(names2)])

    return Expr(:call, Row{newnames}, exprs...)
end

Base.hcat(r1::Union{Cell,Row}, r2::Union{Cell,Row}, rs::Union{Cell,Row}...) = hcat(hcat(r1, r2), rs...)

# copy
@generated function Base.copy{Names}(r::Row{Names})
    exprs = [:(copy(r.($j))) for j = 1:length(Names)]
    return Expr(:call, Row{Names}, exprs...)
end


text:
"""
`Maker.GenericTarget`

The type created by `task()`. Fields expected to be accessed publicly
include:

- `name::UTF8String`
- `dependencies::Vector{UTF8String}`
- `description::UTF8String`
"""
type GenericTarget <: AbstractTarget
    name::UTF8String
    dependencies::Vector{UTF8String}
    description::UTF8String
    action::Function
    timestamp::DateTime
    funhash::UInt64
    isstale::Bool
end

function isstale(t::GenericTarget)
    if t.isstale
        return true
    end
    ds = dependencies(t)
    isempty(ds) ? true : timestamp(t) < maximum([timestamp(d) for d in ds])
end

"""
```julia
task(action::Function, name::AbstractString, dependencies=[])
task(name::AbstractString, dependencies=[])
```
Define and register a `task` target for Maker.jl.

- `action` is the function that operates when the target is
  executed.

- `name` refers to the name of the task.

- `dependencies` refers to the name (`AbstractString`) or names
  (`Vector{AbstractString}`) of targets that need to be satisfied
  for this target before running the `action`. These are also referred
  to as prerequisites.

Targets are registered globally.

`task` targets are generic targets used to define actions and dependencies. If a
task does not have dependencies, it always runs. If a task has dependencies,
the task will run if any of the dependencies has a timestamp newer than the
last run of the task. If a file target depends on a task target without
dependencies, it will always update.

If the `action` or `dependencies` of a target are redefined, the
target will be marked as stale, and the action will be updated
at the next target check.

See also `make`, `file`, and `variable`. `task` registers a `GenericTarget` type.
"""
task(action::Function, name::AbstractString, dependencies=UTF8String[]) =
    target(GenericTarget, name, action, dependencies)
task(name::AbstractString, dependencies=UTF8String[]) =
    target(GenericTarget, name, ()->nothing, dependencies)


text:
using DataFrames
using LinearAlgebra
using Statistics: std, var, mean, cov
using HighFrequencyCovariance
using Random

brownian_corr_matrix = Hermitian([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
assets = [:BARC, :HSBC, :VODL, :RYAL]
twister = MersenneTwister(1)

ts1, true_covar, micro_noise, update_rates = generate_random_path(4, 5000; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister))
ts2, true_covar, micro_noise, update_rates = generate_random_path(4, 50000; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister))

# Testing constructors of SortedDataFrames
ts1 = SortedDataFrame(ts1.df, ts1.time, ts1.grouping, ts1.value,  ts1.groupingrows, ts1.time_period_per_unit)
# Subsetting to time.
ts1 = subset_to_time(ts1, 1670) # Will not delete anything
ts1 = subset_to_time(ts1, 1400) # Will delete some stuff.
# Subsetting to tick.
ts2 = subset_to_tick(ts2, 49000)

# Getting ticks per asset
tpa = ticks_per_asset(ts1)
all(values(tpa) .> 1)

iscloser(a,b) = (a.Correlation_error + a.Volatility_error < b.Correlation_error + b.Volatility_error)

# Preav Convergence
preav_estimate1 = preaveraged_covariance(ts1, assets)
preav_estimate2 = preaveraged_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(preav_estimate2, true_covar), calculate_mean_abs_distance(preav_estimate1, true_covar))
is_psd_matrix(preav_estimate1)
is_psd_matrix(preav_estimate2)
valid_correlation_matrix(preav_estimate1)
valid_correlation_matrix(preav_estimate2)

# simple Convergence
simple_estimate1 = simple_covariance(ts1, assets)
simple_estimate2 = simple_covariance(ts2, assets)
# The below fails. But no reason to  believe that the simple method would be consistant anyway.
# iscloser(calculate_mean_abs_distance(simple_estimate2, true_covar), calculate_mean_abs_distance(simple_estimate1, true_covar))
valid_correlation_matrix(simple_estimate1)
valid_correlation_matrix(simple_estimate2)

# bnhls Convergence
bnhls_estimate1 = bnhls_covariance(ts1, assets; regularisation = :eigenvalue_clean) # Need to change this regularsation because such little data means regu
bnhls_estimate2 = bnhls_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(bnhls_estimate2, true_covar), calculate_mean_abs_distance(bnhls_estimate1, true_covar))
valid_correlation_matrix(bnhls_estimate1)
valid_correlation_matrix(bnhls_estimate2)

# Preav Convergence
spectral_estimate1 = spectral_covariance(ts1, assets; num_blocks = 1) # not many observations so need to reduce the number of blocks here
spectral_estimate2 = spectral_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(spectral_estimate2, true_covar), calculate_mean_abs_distance(spectral_estimate1, true_covar))
valid_correlation_matrix(spectral_estimate1)
valid_correlation_matrix(spectral_estimate2)

# two scales Convergence
two_scales_estimate1 = two_scales_covariance(ts1, assets)
two_scales_estimate2 = two_scales_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(two_scales_estimate2, true_covar), calculate_mean_abs_distance(two_scales_estimate1, true_covar))
valid_correlation_matrix(two_scales_estimate1)
valid_correlation_matrix(two_scales_estimate2)

#############################
 # Serialisation and deserialisation

true_df = to_dataframe(true_covar)
reconstituted_df = dataframe_to_covariancematrix(true_df)
calculate_mean_abs_distance(true_covar, reconstituted_df).Correlation_error .< 10*eps()
calculate_mean_abs_distance(true_covar, reconstituted_df).Volatility_error .< 10*eps()

# Other regularisation algos:
two_scales_estimate_iden = two_scales_covariance(ts2, assets; regularisation = :identity_regularisation)
two_scales_estimate_nearest_corr = two_scales_covariance(ts2, assets; regularisation = :nearest_correlation_matrix)
two_scales_estimate_nearest_psd = two_scales_covariance(ts2, assets; regularisation = :nearest_psd_matrix)
two_scales_estimate_eigen = two_scales_covariance(ts2, assets; regularisation = :eigenvalue_clean)

# Running regularistation on a CovarianceMatrix's correlation matrix.
valid_correlation_matrix(two_scales_estimate_nearest_corr)

psd_mat = two_scales_estimate_nearest_corr

reg1 = identity_regularisation(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg1).Correlation_error .> 10*eps()
calculate_mean_abs_distance(psd_mat, reg1).Volatility_error .> 10*eps()

reg1_corr = identity_regularisation(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg1_corr).Correlation_error .> 10*eps()
calculate_mean_abs_distance(psd_mat, reg1_corr).Volatility_error .< 10*eps()

reg2 = nearest_correlation_matrix(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg2).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg2).Volatility_error .< 10*eps()

reg3 = nearest_psd_matrix(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg3).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg3).Volatility_error .< 1000*eps()

reg3_corr = nearest_psd_matrix(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg3_corr).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg3_corr).Volatility_error .< 10*eps()

reg4 = eigenvalue_clean(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg4).Correlation_error .> 2*eps()
calculate_mean_abs_distance(psd_mat, reg4).Volatility_error .> 10*eps()

reg4_cov = eigenvalue_clean(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg4_cov).Correlation_error .> 2*eps()
calculate_mean_abs_distance(psd_mat, reg4_cov).Volatility_error .< 10*eps()

# Testing blocking and regularisation.
blocking_dd = put_assets_into_blocks_by_trading_frequency(ts2, 1.3, :spectral_covariance)
block_estimate = blockwise_estimation(ts2, blocking_dd)
block_estimate = identity_regularisation(block_estimate, ts2)
iscloser(calculate_mean_abs_distance(block_estimate, true_covar), calculate_mean_abs_distance(spectral_estimate2, true_covar))

# Testing getting noise
noise = estimate_microstructure_noise(ts1, assets)
noise_from_2s = two_scales_volatility(ts1, assets)[2]
C = merge(-, noise, noise_from_2s)
maximum(collect(values(C))) < 1E-13

# Testing print functions
print(ts1)
show(ts1)
print(true_covar)
show(true_covar)
print(preav_estimate2, 3, 6)
show(preav_estimate2, 3, 6)
true # Otherwise it tends to throw from not evaluating to Bool in testing.


text:
# ---
# title: 327. Count of Range Sum
# id: problem327
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Divide and Conquer, Sort, Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/count-of-range-sum/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return the number of range sums that lie in
# `[lower, upper]` inclusive.
# Range sum `S(i, j)` is defined as the sum of the elements in `nums` between
# indices `i` and `j` (`i` ≤ `j`), inclusive.
#
# **Note:**
# A naive algorithm of _O_ ( _n_ 2) is trivial. You MUST do better than that.
#
# **Example:**
#
#
#
#     Input: _nums_ = [-2,5,-1], _lower_ = -2, _upper_ = 2,
#     Output: 3
#     Explanation: The three ranges are : [0,0], [2,2], [0,2] and their respective sums are: -2, -1, 2.
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module FSAtesting

using FixedSizeArrays
using FactCheck, Base.Test
using Compat

import FixedSizeArrays: similar_type

immutable Normal{N, T} <: FixedVector{N, T}
    values::NTuple{N, T}
end
immutable D3{N1, N2, N3, T} <: FixedArray{T, 3, Tuple{N1, N2, N3}}
    values::NTuple{N1, NTuple{N2, NTuple{N3, T}}}
end
immutable RGB{T} <: FixedVectorNoTuple{3, T}
    r::T
    g::T
    b::T
end

# subtyping:
immutable TestType{N,T} <: FixedVector{N,T}
    values::NTuple{N,T}
end

# Custom FSA with non-parameterized size and eltype
immutable Coord2D <: FixedVectorNoTuple{2,Float64}
    x::Float64
    y::Float64
end

typealias Vec1d Vec{1, Float64}
typealias Vec2d Vec{2, Float64}
typealias Vec3d Vec{3, Float64}
typealias Vec4d Vec{4, Float64}
typealias Vec3f Vec{3, Float32}
typealias Mat2d Mat{2,2, Float64}
typealias Mat3d Mat{3,3, Float64}
typealias Mat4d Mat{4,4, Float64}

# Compatibility hacks for 0.5 APL-style array slicing
if VERSION < v"0.5.0-dev+1195"
    # Remove all dimensions of size 1
    compatsqueeze(A) = squeeze(A,(find(collect(size(A)).==1)...))
else
    compatsqueeze(A) = A
end
function test()
facts("FixedSizeArrays") do

include("typeinf.jl")

context("fsa macro") do
    a = 1
    a1 = @fsa([a,2,3])
    a2 = @fsa([a 2 3])
    a3 = @fsa([a;2;3])
    a4 = @fsa([a 2;3 4])
    a5 = @fsa([a 2 3;4 5 6])
    a6 = @fsa([a 2;3 4;5 6])

    @fact a1 --> Vec(a,2,3)
    @fact a2 --> Mat((a,),(2,),(3,))
    @fact a3 --> Mat(((a,2,3),))
    @fact a4 --> Mat(((a,3),(2,4)))
    @fact a5 --> Mat(((a,4),(2,5),(3,6)))
    @fact a6 --> Mat(((a,3,5),(2,4,6)))
end

context("core") do
    context("ndims") do
        @fact ndims(D3) --> 3
        @fact ndims(Mat) --> 2
        @fact ndims(Vec) --> 1
        @fact ndims(Vec(1,2,3)) --> 1

        @fact ndims(D3{3,3,3}) --> 3
        @fact ndims(Mat{3,3}) --> 2
        @fact ndims(Vec{3}) --> 1

        @fact ndims(D3{3,3,3,Int}) --> 3
        @fact ndims(Mat{3,3,Int}) --> 2
        @fact ndims(Vec{3,Int}) --> 1
    end
    context("size_or") do
        @fact size_or(Mat, nothing) --> nothing
        @fact size_or(Mat{4}, nothing) --> nothing
        @fact size_or(Mat{4,4}, nothing) --> (4,4)
        @fact size_or(Mat{4,4, Float32}, nothing) --> (4,4)

        @fact size_or(Vec, nothing) --> nothing
        @fact size_or(Vec{4}, nothing) --> (4,)
        @fact size_or(Vec{4,Float32}, nothing) --> (4,)
        @fact size_or(FixedArray, nothing) --> nothing

    end
    context("eltype_or") do
        @fact eltype_or(Mat, nothing) --> nothing
        @fact eltype_or(Mat{4}, nothing) --> nothing
        @fact eltype_or(Mat{4,4}, nothing) --> nothing
        @fact eltype_or(Mat{4,4, Float32}, nothing) --> Float32

        @fact eltype_or(Vec, nothing) --> nothing
        @fact eltype_or(Vec{4}, nothing) --> nothing
        @fact eltype_or(Vec{4,Float32}, nothing) --> Float32

        @fact eltype_or(FixedArray, nothing) --> nothing

    end
    context("ndims_or") do
        @fact ndims_or(Mat, nothing) --> 2
        @fact ndims_or(Mat{4}, nothing) --> 2
        @fact ndims_or(Mat{4,4}, nothing) --> 2
        @fact ndims_or(Mat{4,4, Float32}, nothing) --> 2

        @fact ndims_or(Vec, nothing) --> 1
        @fact ndims_or(Vec{4}, nothing) --> 1
        @fact ndims_or(Vec{4, Float64}, nothing) --> 1

        @fact ndims_or(FixedArray, nothing) --> nothing
    end

    context("similar_type") do
        @fact similar_type(Vec{3,Int}, Float32) --> Vec{3, Float32}
        @fact similar_type(Vec{3}, Float32) --> Vec{3, Float32}
        @fact similar_type(Vec, Float32, (3,)) --> Vec{3, Float32}
        @fact similar_type(Vec, Float32, (1,2)) --> Mat{1,2, Float32}

        @fact similar_type(RGB, Float32) --> RGB{Float32}
        @fact similar_type(RGB{Float32}, Int) --> RGB{Int}
        @fact similar_type(RGB{Float32}, Int, (3,)) --> RGB{Int}
        @fact similar_type(RGB{Float32}, Int, (2,2)) --> Mat{2,2,Int}

        @fact similar_type(Mat{3,3,Int}, Float32) --> Mat{3,3,Float32}
        @fact similar_type(Mat, Float32, (3,3))   --> Mat{3,3,Float32}
        @fact similar_type(Mat{2,2,Int}, (3,3))   --> Mat{3,3,Int}

        @fact similar_type(Coord2D, Float64, (2,)) --> Coord2D
        @fact similar_type(Coord2D, Int, (2,))     --> Vec{2,Int}
        @fact similar_type(Coord2D, Float64, (3,)) --> Vec{3,Float64}
    end

    context("construct_similar") do
        @fact construct_similar(Vec{3,Int}, (1.0f0,2)) --> exactly(Vec{2,Float32}(1,2))
        @fact construct_similar(Vec{2}, (1,2,3))       --> exactly(Vec{3,Int}(1,2,3))
        @fact construct_similar(Vec, (1.0,2))          --> exactly(Vec{2,Float64}(1,2))

        @fact construct_similar(RGB, (1,2,3))                --> exactly(RGB{Int}(1,2,3))
        @fact construct_similar(RGB{Float32}, (1.0,2.0,3.0)) --> exactly(RGB{Float64}(1.0,2.0,3.0))

        @fact construct_similar(Mat{3,3,Int}, ((1.0f0,2),(1.0,2))) --> exactly(Mat{2,2,Float64}((1,2),(1,2)))
        @fact construct_similar(Mat, ((1,2),))                     --> exactly(Mat{2,1,Int}(((1,2),)))
    end

    context("nan") do
        for (p, r) in (
                (Point{2, Float32}(NaN, 1), true),
                (Point{2, Float64}(1, NaN), true),
                (Vec{11, Float64}(NaN), true),
                (Point{2, Float32}(1, 1), false),
                (RGB{Float32}(NaN), true),
            )
            @fact isnan(p) --> r
        end
    end

end

context("Array of FixedArrays") do

    N = 100
    a = Point{3, Float32}[Point{3, Float32}(0.7132) for i=1:N]
    b = RGB{Float32}[RGB{Float32}(52.293) for i=1:N]

    c = Point{3, Float64}[Point{3, Float64}(typemin(Float64)), a..., Point{3, Float64}(typemax(Float64))]
    d = RGB{Float64}[RGB(typemin(Float64)), b..., RGB(typemax(Float64))]

    context("reduce") do
        sa = sum(a)
        ma = mean(a)
        sb = sum(b)
        mb = mean(b)
        for i=1:3
            @fact sa[i]  --> roughly(Float32(0.7132*N))
            @fact ma[i]  --> roughly(Float32(0.7132*N)/ N)

            @fact sb[i]  --> roughly(Float32(52.293*N))
            @fact mb[i]  --> roughly(Float32(52.293*N)/ N)
        end

        @fact maximum(c) --> Point{3, Float32}(typemax(Float64))
        @fact minimum(c) --> Point{3, Float32}(typemin(Float64))

        @fact maximum(d) --> RGB(typemax(Float64))
        @fact minimum(d) --> RGB(typemin(Float64))

        @fact extrema(c) --> (minimum(c), maximum(c))
    end

    context("array ops") do
        for op in (.+, .-,.*, ./, .\, +, -)
            @fact typeof(op(a, 1f0)) --> typeof(a)
            @fact typeof(op(1f0, a)) --> typeof(a)
        end

        af = a + 1f0
        bf = b + 1f0
        aff = a + Point{3, Float32}(1)
        bff = b + RGB{Float32}(1)
        afd = a .+ 1f0
        bfd = b .+ 1f0
        @inferred(b .* 1f0)
        for i=1:N
            @fact a[1] + 1f0 --> af[i]
            @fact b[1] + 1f0 --> bf[i]
            @fact a[1] + 1f0 --> aff[i]
            @fact b[1] + 1f0 --> bff[i]
            @fact a[1] + 1f0 --> afd[i]
            @fact b[1] + 1f0 --> bfd[i]
        end
    end
    context("Show") do
        m1 = rand(Mat4d, 2)
        m2 = rand(RGB{Float32}, 2)
        m3 = rand(Vec3f, 2)
        println(m1)
        println(m2)
        println(m3)
        showcompact(Point(1,2,3))
    end
end

context("Constructor FixedVectorNoTuple") do
    for T=[Float32, Float64, Int, UInt, UInt32, UInt8]
        context("$T") do
            r = rand(T)
            x = RGB{Int}[RGB(1) for i=1:10]
            @fact RGB{Float32}(["0.222", "9.8822", "29.999"]) --> RGB{Float32}(0.222, 9.8822, 29.999)
            @fact RGB(["0.222", "9.8822", "29.999"]) --> RGB{Float64}(0.222, 9.8822, 29.999)
            @fact typeof(map(RGB{Float32}, x))  --> Vector{RGB{Float32}}
            @fact RGB{T}(r)                     --> RGB(r,r,r)
            @fact RGB{T}(Vec(r,r,r))            --> RGB(r,r,r)
            @fact RGB{T}([r,r,r])               --> RGB(r,r,r)
            @fact length(RGB{T}([r,r,r]))       --> 3
            @fact length(RGB{T})                --> 3
            @fact eltype(RGB{T}([r,r,r]))       --> T
            @fact eltype(RGB{T})                --> T
            @fact typeof(RGB(r,r,r))            --> RGB{T}
            @fact typeof(RGB{T}(1))             --> RGB{T}
            @fact typeof(RGB{T}(1,2,3))         --> RGB{T}
            @fact ndims(RGB{T}(1,2,3))          --> 1

            @fact typeof(RGB{T}(1f0))           --> RGB{T}
            @fact typeof(RGB{T}(1f0,2f0,3f0))   --> RGB{T}
            @fact typeof(RGB{T}(1f0, 2, 3.0))   --> RGB{T}
            @fact typeof(RGB(1f0, 2, 3.0))      --> RGB{Float64}
            @fact typeof(RGB{Int}(1f0, 2, 3.0)) --> RGB{Int}
            @fact_throws DimensionMismatch RGB((1,2,3), (2,3,4))
        end
    end
end

# A little brutal, but hey.... Better redudantant tests, than not enough tests
context("Constructor ") do
    context("Rand") do
        #Win32 seems to fail for rand(Vec4d)
        @fact typeof(rand(Vec4d)) --> Vec4d
        @fact typeof(rand(Mat4d)) --> Mat4d

        @fact typeof(rand(Mat{4,2, Int})) --> Mat{4,2, Int}
        @fact typeof(rand(Vec{7, Int})) --> Vec{7, Int}
        @fact typeof(rand(Vec{7, Int}, 1:7)) --> Vec{7, Int}
        @fact typeof(rand(Mat4d, -20f0:0.192f0:230f0)) --> Mat4d
        @fact typeof(rand(Mat{4,21,Float32}, -20f0:0.192f0:230f0)) --> Mat{4,21,Float32}

        x = rand(D3{4,4,4, Float32})
        @fact typeof(x) --> D3{4,4,4, Float32}
        @fact eltype(x) --> Float32
        @fact size(x) --> (4,4,4)
        @fact typeof(rand(Vec4d, 5,5)) --> Matrix{Vec4d}
    end
    context("Randn") do
        @fact typeof(randn(Base.Random.GLOBAL_RNG, Vec4d)) --> Vec4d
        @fact typeof(randn(Vec4d)) --> Vec4d
        @fact typeof(randn(Mat4d)) --> Mat4d
        @fact typeof(randn(Mat{4,2, Complex{Float64}})) --> Mat{4,2, Complex{Float64}}
        @fact typeof(randn(Vec{7, Complex{Float64}})) --> Vec{7, Complex{Float64}}
    end
    context("Zero") do
        @fact typeof(zero(Vec4d)) --> Vec4d
        @fact typeof(zero(Mat4d)) --> Mat4d

        @fact typeof(zero(Mat{4,2, Int})) --> Mat{4,2, Int}
        @fact typeof(zero(Vec{7, Int})) --> Vec{7, Int}
        @fact zero(Vec((1,2))) --> Vec((0,0))
        @fact zero(Vec((1.0,2.0))) --> Vec((0.0,0.0))
    end

    context("eye") do
        @fact typeof(eye(Mat4d)) --> Mat4d
        @fact typeof(eye(Mat{4,2, Int})) --> Mat{4,2, Int}
    end
    context("one") do
        x = one(Mat{4,2, Int})
        @fact typeof(one(Mat4d)) --> Mat4d
        @fact typeof(x) --> Mat{4,2, Int}
        @fact all(x-> x==1, x) --> true
    end

    context("unit") do
        u4 = unit(Vec4d, 1)
        u7 = unit(Vec{7, Int}, 7)
        @fact typeof(u4) --> Vec4d
        @fact typeof(u7) --> Vec{7, Int}
        @fact u4[1] --> 1.0
        @fact u4[2:end] --> (0.,0.,0.)

        @fact u7[end] --> 1
        @fact u7[1:end-1] --> (0,0,0,0,0,0)
    end
    for N=(1,10)
        context("construction, conversion, $N") do
            for VT=[Point, Vec], VT2=[Normal, Vec], ET=[Float32, Int, UInt], ET2=[Float64, UInt, Float32]
                rand_range  = ET(1):ET(10)
                rand_range2 = ET2(1):ET2(10)
                rn = rand(rand_range, N)
                v0 = VT(rn)
                # parse constructor:
                @fact VT{N, ET}(map(string, rn)) --> v0
                # multi constructor
                v1 = VT{N, ET}(rn...)
                @fact v1 --> v0
                @fact typeof(v1) --> VT{N, ET}
                @fact length(v1) --> N
                @fact eltype(v1) --> ET
                @fact ndims(v1) --> 1

                @fact length(typeof(v1)) --> N
                @fact eltype(typeof(v1)) --> ET

                for i=1:N
                    @fact v1[i] --> rn[i]
                end
                # from other FSA without parameters
                v2 = VT2(v1)

                @fact typeof(v2) --> VT2{N, ET}
                @fact length(v2) --> N
                @fact eltype(v2) --> ET
                for i=1:N
                    @fact v2[i] --> v1[i]
                end
                # from other FSA with parameters
                for i=1:N
                    v3 = VT2{i, ET2}(v1)
                    @fact typeof(v3) --> VT2{i, ET2}
                    @fact length(v3) --> i
                    @fact eltype(v3) --> ET2
                    for i=1:i
                        @fact v3[i] --> ET2(v2[i])
                    end
                end
                # from single
                r  = rand(rand_range)
                r2 = rand(rand_range2)
                v1 = VT{N, ET}(r)
                v2 = VT{N, ET2}(r)
                v3 = VT{N, ET}(r2)
                v4 = VT{N, ET2}(r2)

                for i=1:N
                    @fact v1[i] --> r
                    @fact v2[i] --> ET2(r)
                    @fact v3[i] --> r2
                    @fact v4[i] --> ET2(r2)
                end
                x = VT{N, ET}[VT{N, ET}(1) for i=1:10]
                x1 = VT2{N, ET}[VT{N, ET}(1) for i=1:10]
                x2 = map(VT2, x)
                x3 = map(VT, x2)
                @fact typeof(x)  --> Vector{VT{N, ET}}
                @fact typeof(x1) --> Vector{VT2{N, ET}}
                @fact typeof(x2) --> Vector{VT2{N, ET}}
                @fact typeof(x3) --> Vector{VT{N, ET}}
                @fact x3         --> x

                # Construction with only N, issue #56
                @fact VT{N}(ET(1)) --> Vec{N, ET}(1)
                @fact VT{N}(ntuple(x->ET(1), N)...) --> Vec{N, ET}(1)
            end
        end
    end
end

context("Constructors") do
    context("FixedVector: unary, from FixedVector") do
        @fact typeof(Vec3f(1,1,1))     --> Vec{3, Float32}
        @fact typeof(Vec3f(1,1f0,1))   --> Vec{3, Float32}
        @fact typeof(Vec3f(1f0,1,1.0)) --> Vec{3, Float32}
        @fact eltype(Vec3f(1f0,1,1.0)) --> Float32

        @fact typeof(Vec3f(1))      --> Vec{3, Float32}
        @fact typeof(Vec3f(0))      --> Vec{3, Float32}
        @fact Vec3f(1.0)             --> Vec(1f0,1f0,1f0)
        @fact Vec3f(1.0f0)             --> Vec(1f0,1f0,1f0)
        @fact Vec3f(1.0f0)             --> Vec3f(1)
        @fact Vec(1.0, 1.0, 1.0)     --> Vec3d(1)
        @fact Vec2d(Vec3d(1))         --> Vec(1.0, 1.0)
        @fact Vec(Vec3d(1), 1.0)     --> Vec4d(1)
        @fact Vec(Vec3d(1), 1)         --> Vec4d(1)
        @fact Vec3d(Vec3f(1.0))     --> Vec3d(1.0)
    end
end

context("map") do
    context("Vec and AbstractVector") do
        # Unary, binary & ternary map with specified output type
        @fact map(-, Vec{3,Float64}, Vec(1,2,3)) --> exactly(Vec{3,Float64}(-1,-2,-3))
        @fact map(+, Vec{3,Float64}, [1,2,3], Vec(1,2,3)) --> exactly(Vec{3,Float64}(2,4,6))
        @fact map(+, Vec{3,Float64}, [1,2,3], Vec(1,2,3), 1:3) --> exactly(Vec{3,Float64}(3,6,9))

        # Unary and binary map with deduced output types
        @fact map(-, Vec(1,2,3)) --> exactly(Vec{3,Int}(-1,-2,-3))
        @fact map(+, Vec(1,2,3), [1,2,3]) --> exactly(Vec{3,Int}(2,4,6))
        @fact map(+, [1,2,3], Vec(1,2,3)) --> exactly(Vec{3,Int}(2,4,6))
        @fact map(+, Vec(1,2,3), Vec(1,2,3)) --> exactly(Vec{3,Int}(2,4,6))
        # Some other `AbstractArray`s
        @fact map(+, Vec(1,2,3), 1:3) --> exactly(Vec{3,Int}(2,4,6))
        @fact map(+, 1:3, Vec(1,2,3)) --> exactly(Vec{3,Int}(2,4,6))

        # Binary map with mixed types
        @fact map(>, Vec(0.0,2.0), Vec(1,1)) --> exactly(Vec{2,Bool}(false,true))
        @fact map(+, Vec(0.0,0.0), Vec(1,1)) --> exactly(Vec{2,Float64}(1.0,1.0))
    end

    context("FixedVectorNoTuple") do
        # RGB with specified output
        @fact map(-, RGB{Float64}, RGB(1.0, 2.0, 3.0)) --> exactly(RGB{Float64}(-1.0, -2.0, -3.0))
        @fact map(-, RGB{Float64}, [1.0, 2.0, 3.0]) --> exactly(RGB{Float64}(-1.0, -2.0, -3.0))

        # RGB and AbstractVector interop
        @fact map(+, RGB(1.0, 2.0, 3.0), RGB(1.0, 2.0, 3.0)) --> exactly(RGB{Float64}(2.0, 4.0, 6.0))
        @fact map(+, RGB(1.0, 2.0, 3.0), [1.0, 2.0, 3.0]) --> exactly(RGB{Float64}(2.0, 4.0, 6.0))
        @fact map(+, [1.0, 2.0, 3.0], RGB(1.0, 2.0, 3.0)) --> exactly(RGB{Float64}(2.0, 4.0, 6.0))
        @fact map(+, RGB{Int}(1, 2, 3), RGB(1.0, 2.0, 3.0)) --> exactly(RGB{Float64}(2.0, 4.0, 6.0))
    end

    context("Mat and AbstractMatrix") do
        @fact map(+, Mat{2,2,Int}(((1,2),(3,4))), Mat{2,2,Int}(((1,2),(3,4)))) --> exactly(Mat{2,2,Int}(((2,4),(6,8))))
        @fact map(+, Mat{2,2,Int}(((1,2),(3,4))), [1 3; 2 4]) --> exactly(Mat{2,2,Int}(((2,4),(6,8))))
        @fact map(+, [1 3; 2 4], Mat{2,2,Int}(((1,2),(3,4)))) --> exactly(Mat{2,2,Int}(((2,4),(6,8))))
    end

    context("Size checking") do
        @fact_throws DimensionMismatch map(+, Vec(1,2,3), Vec(1,1))
        @fact_throws DimensionMismatch map(+, Vec(1,1), Vec(1,2,3))
        @fact_throws DimensionMismatch map(+, Vec(1,2,3), [1,1])
        @fact_throws DimensionMismatch map(+, [1,1], Vec(1,2,3))
        @fact_throws DimensionMismatch map(+, Vec(1,2,3), 1:2)
        @fact_throws DimensionMismatch map(+, 1:2, Vec(1,2,3))
        @fact_throws DimensionMismatch map(+, Vec(1,2,3), [1 2 3])
    end

    context("Broadcast of scalars") do
        # Arguably not the right thing for map(), but neither do we have a full
        # broadcast implementation yet...
        @fact map(+, Vec{3,Float64}, Vec(1,2,3), 1.0) --> exactly(Vec{3,Float64}(2,3,4))
        @fact map(+, Vec{3,Float64}, 1.0, Vec(1,2,3)) --> exactly(Vec{3,Float64}(2,3,4))
        @fact map(+, 1.0, Vec(1,2,3)) --> exactly(Vec{3,Float64}(2,3,4))
        @fact map(+, Vec(1,2,3), 1.0) --> exactly(Vec{3,Float64}(2,3,4))
    end
end

v2 = Vec(6.0,5.0,4.0)
v1 = Vec(1.0,2.0,3.0)
vi = Vec(1,2,3)
v2 = Vec(6.0,5.0,4.0)
v1c = Vec(6.0+3.0im,5.0-2im,4.0+0.0im)
v2c = v1 + v2*im
v2c = Vec(1.0 + 6.0im, 2.0 + 5.0im, 3.0 + 4.0im)

context("Complex Ops") do
    context("dot product") do
        @fact dot(v1c,v2c) --> dot([6.0+3.0im,5.0-2im,4.0+0.0im], [1.0,2.0,3.0] + [6.0,5.0,4.0]*im)
        @fact Vector(transpose(v1c)*v2c) --> [6.0+3.0im 5.0-2im 4.0+0.0im]*([1.0,2.0,3.0] + [6.0,5.0,4.0]*im)
        @fact Matrix(v2c*transpose(v1c)) --> ([1.0,2.0,3.0] + [6.0,5.0,4.0]*im)*[6.0+3.0im 5.0-2im 4.0+0.0im]
    end
end

context("Destructure") do
    rgb_ref = Int[1 2 3 4;
               2 4 6 8;
               3 6 9 12]
    rgb_ref_set = Int[1 10 10 4;
                   2 10 10 8;
                   3 10 10 12]
    # Test destructure
    rgb = RGB{Int}[RGB(i,2*i,3*i) for i=1:4]
    @fact destructure(rgb) --> rgb_ref
    destructure(rgb)[:,2:end-1] = 10
    @fact destructure(rgb) --> rgb_ref_set

    # Explicitly test DestructuredArray.  This wrapper type isn't used by
    # destructure() for plain old dense arrays, since a reinterpret is faster.
    rgb = RGB{Int}[RGB(i,2*i,3*i) for i=1:4]
    @fact FixedSizeArrays.DestructuredArray(rgb) --> rgb_ref
    destructure(rgb)[:,2:end-1] = 10
    @fact FixedSizeArrays.DestructuredArray(rgb) --> rgb_ref_set

    # destructure() with 2D FSA
    A = [@fsa([i 2*i; 3*i 4*i]) for i=1:2]
    @fact destructure(A) --> cat(3, [1 2; 3 4], [2 4; 6 8])
end

context("Indexing") do
    context("FixedVector") do
        @fact setindex(v1, 88.9, 1) --> Vec(88.9,2.0,3.0)
        @fact v1[1] --> 1.0
        @fact v1[2] --> 2.0
        @fact v1[3] --> 3.0
        @fact v1[1:3] --> (1.0, 2.0, 3.0)
        @fact v1[1:2] --> (1.0, 2.0)
        @fact v1[1:1] --> (1.0,)
        @fact v1[(1,2)] --> (1.0,2.0)
        @fact v1[(2,1)] --> (2.0,1.0)
        @fact_throws BoundsError v1[-1]
        @fact_throws BoundsError v1[0]
        @fact_throws BoundsError v1[4]
        @fact row(v1, 1) --> (1.0,)
    end
    m = Mat{4,4,Int}(
        (1,2,3,4),
        (5,6,7,8),
        (9,10,11,12),
        (13,14,15,16)
    )
    context("FixedMatrix") do
        @fact setindex(m, 42.0, 2,2) --> Mat{4,4,Int}(
            (1,2,3,4),
            (5,42.0,7,8),
            (9,10,11,12),
            (13,14,15,16)
        )
        @fact m[1] --> 1
        @fact m[2] --> 2
        @fact m[10] --> 10
        @fact m[2,2] --> 6
        @fact m[3,4] --> 15
        @fact m[1:4, 1] --> (1,5,9,13)
        @fact m[1, 1:4] --> (1,2,3,4)
        @fact_throws BoundsError m[-1]
        @fact_throws BoundsError m[0]
        @fact_throws BoundsError m[17]
        @fact_throws BoundsError m[5,1]
        @fact_throws BoundsError m[-1,1]
        @fact_throws BoundsError m[0,0]

        @fact row(m, 1) --> (1,5,9,13)

    end

    context("fslice") do
        context("getindex") do
            rgb = RGB{Int}[RGB(i,2*i,3*i) for i=1:10]

            # Plain indexing
            @fact @fslice(rgb[1,2]) --> rgb[2].r
            @fact @fslice(rgb[2,5]) --> rgb[5].g
            @fact @fslice(rgb[3,8]) --> rgb[8].b

            # Slicing along fixed dims
            @fact @fslice(rgb[:,1]) --> rgb[1]
            @fact @fslice(rgb[:,end]) --> rgb[end]

            # Slicing across fixed dims
            @fact compatsqueeze(@fslice(rgb[1,:]))  --> [c.r for c in rgb]
            @fact compatsqueeze(@fslice(rgb[2,:]))  --> [c.g for c in rgb]
            @fact compatsqueeze(@fslice(rgb[3,:]))  --> [c.b for c in rgb]
            # Slicing across fixed dims with field names
            @fact compatsqueeze(@fslice(rgb[:r,:])) --> [c.r for c in rgb]
            @fact compatsqueeze(@fslice(rgb[:g,:])) --> [c.g for c in rgb]
            @fact compatsqueeze(@fslice(rgb[:b,:])) --> [c.b for c in rgb]

            # Slicing FSAs with two fixed dimensions
            N = 3
            A = Mat{2,2,Int}[@fsa([i 2*i; 3*j 4*j]) for i=1:N, j=1:N]
            for i=1:N,j=1:N
                @fact compatsqueeze(@fslice(A[:,:,i,j])) --> A[i,j]
            end
            @fact compatsqueeze(@fslice(A[1,1,:,1])) --> [A[i,1][1,1] for i=1:N]
            @fact compatsqueeze(@fslice(A[end,end,end,:])) --> [A[end,j][end,end] for j=1:N]
            @fact compatsqueeze(@fslice(A[1,[1,end],1,1])) --> [A[1,1][1,1], A[1,1][1,end]]
        end

        context("setindex") do
            rgb = RGB{Int}[RGB(i,2*i,3*i) for i=1:10]

            @fslice rgb[:r,:] = -1
            @fslice rgb[:g,:] .+= 1
            @fslice rgb[3,:] = -3
            @fact rgb --> RGB{Int}[RGB(-1,2*i+1,-3) for i=1:10]
        end
    end
end

context("Ops") do
    context("revers") do
        @fact reverse(Vec(1,2,3,4,5)) --> Vec(5,4,3,2,1)
    end
    context("Negation") do
        @fact @inferred(-v1) --> Vec(-1.0,-2.0,-3.0)
        @fact isa(-v1, Vec3d) --> true
    end

    context("Addition") do
        @fact @inferred(v1+v2) --> Vec3d(7.0,7.0,7.0)
        @fact @inferred(RGB(1,2,3) + RGB(2,2,2)) --> exactly(RGB{Int}(3,4,5))
        @fact @inferred(Coord2D(1,2) + Coord2D(3,4)) --> exactly(Coord2D(4,6))
    end
    context("Subtraction") do
        @fact @inferred(v2-v1) --> Vec3d(5.0,3.0,1.0)
        @fact @inferred(RGB(1,2,3) - RGB(2,2,2)) --> exactly(RGB{Int}(-1,0,1))
        @fact @inferred(Coord2D(1,2) - Coord2D(3,4)) --> exactly(Coord2D(-2,-2))
    end
    context("Multiplication") do
        @fact @inferred(v1.*v2) --> Vec3d(6.0,10.0,12.0)
    end
    context("Mixed Type Multiplication") do
        @fact @inferred(vi.*v2) --> Vec3d(6.0,10.0,12.0)
    end
    context("Division") do
        @fact @inferred(v1 ./ v1) --> Vec3d(1.0,1.0,1.0)
    end

    context("Relational") do
        @fact Vec(1,3) .< Vec(2,2) --> exactly(Vec{2,Bool}(true,false))
        @fact RGB(1,2,3) .< RGB(2,2,2) --> exactly(RGB{Bool}(true,false,false))
        @fact Coord2D(1,3) .< Coord2D(2,2) --> exactly(Vec{2,Bool}(true,false))
        end

    context("Scalar") do
        @fact @inferred(1.0 + v1) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(1.0 .+ v1) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(v1 + 1.0) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(v1 .+ 1.0) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(1 + v1) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(1 .+ v1) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(v1 + 1) --> Vec3d(2.0,3.0,4.0)
        @fact @inferred(v1 .+ 1) --> Vec3d(2.0,3.0,4.0)

        @fact @inferred(v1 - 1.0) --> Vec3d(0.0,1.0,2.0)
        @fact @inferred(v1 .- 1.0) --> Vec3d(0.0,1.0,2.0)
        @fact @inferred(1.0 - v1) --> Vec3d(0.0,-1.0,-2.0)
        @fact @inferred(1.0 .- v1) --> Vec3d(0.0,-1.0,-2.0)
        @fact @inferred(v1 - 1) --> Vec3d(0.0,1.0,2.0)
        @fact @inferred(v1 .- 1) --> Vec3d(0.0,1.0,2.0)
        @fact @inferred(1 - v1) --> Vec3d(0.0,-1.0,-2.0)
        @fact @inferred(1 .- v1) --> Vec3d(0.0,-1.0,-2.0)

        @fact @inferred(2.0 * v1) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(2.0 .* v1) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(v1 * 2.0) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(v1 .* 2.0) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(2 * v1) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(2 .* v1) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(v1 * 2) --> Vec3d(2.0,4.0,6.0)
        @fact @inferred(v1 .* 2) --> Vec3d(2.0,4.0,6.0)

        @fact @inferred(v1 / 2.0) --> Vec3d(0.5,1.0,1.5)
        @fact @inferred(v1 ./ 2.0) --> Vec3d(0.5,1.0,1.5)
        @fact @inferred(v1 / 2) --> Vec3d(0.5,1.0,1.5)
        @fact @inferred(v1 ./ 2) --> Vec3d(0.5,1.0,1.5)

        @fact @inferred(12.0 ./ v1) --> Vec3d(12.0,6.0,4.0)
        @fact @inferred(12 ./ v1) --> Vec3d(12.0,6.0,4.0)

        @fact @inferred((v1 .^ 2)) --> Vec3d(1.0,4.0,9.0)
        @fact @inferred((v1 .^ 2.0)) --> Vec3d(1.0,4.0,9.0)
        @fact @inferred((2.0 .^ v1)) --> Vec3d(2.0,4.0,8.0)
        @fact @inferred((2 .^ v1)) --> Vec3d(2.0,4.0,8.0)

                a = Vec(3.2f0)
                @fact @inferred(a+0.2) --> Vec1d(3.2f0+0.2)
                @fact @inferred(0.2+a) --> Vec1d(3.2f0+0.2)
                @fact @inferred(a*0.2) --> Vec1d(3.2f0*0.2)
                @fact @inferred(0.2*a) --> Vec1d(3.2f0*0.2)
                @fact @inferred(a+0.2f0) --> Vec{1,Float32}(3.4f0)
                @fact @inferred(0.2f0+a) --> Vec{1,Float32}(3.4f0)
                @fact @inferred(a*0.2f0) --> Vec{1,Float32}(3.2f0*0.2f0)
                @fact @inferred(0.2f0*a) --> Vec{1,Float32}(3.2f0*0.2f0)
    end
    context("vector norm+cross product") do

        @fact norm(Vec3d(1.0,2.0,2.0)) --> 3.0
        @fact norm(Vec3d(1.0,2.0,2.0),2) --> 3.0
        @fact norm(Vec3d(1.0,2.0,2.0),Inf) --> 2.0
        @fact norm(Vec3d(1.0,2.0,2.0),1) --> 5.0

        # cross product
        @fact cross(v1,v2) --> Vec3d(-7.0,14.0,-7.0)
        @fact isa(cross(v1,v2), Vec3d)  --> true

        @fact cross(vi,v2) --> Vec3d(-7.0,14.0,-7.0)
        @fact isa(cross(vi,v2),Vec3d)  --> true

        a,b = Vec2d(0,1), Vec2d(1,0)
        @fact cross(a,b) --> -1.0
        @fact isa(cross(a,b), Float64) --> true
    end

    context("hypot") do
        a = Vec{2,Int}(1,2)
        b = Vec{2,Float64}(1.,2.)
        @fact hypot(a) --> 2.23606797749979
        @fact hypot(b) --> 2.23606797749979
        @fact hypot(a) == hypot(b) --> true
    end
    context("normalize") do
        a = Vec(3,4)
        b = Vec(3.,4.)
        @fact normalize(a) --> Vec(0.6,0.8)
        @fact normalize(b) --> Vec(0.6,0.8)
    end

    context("reduce") do
        a = rand(Vec{7, Float32})
        x = reduce(+, a)
        y = 0f0
        for elem in a
            y += elem
        end
        @fact y --> x

        a = rand(Mat{7, 9, Cuint})
        x2 = reduce(+, a)
        y2 = Cuint(0)
        for elem in a
            y2 += elem
        end
        @fact y2 --> x2
    end
end

context("Promotion") do
    @fact promote_type(Vec{2,Float64}, Int) --> Vec{2,Float64}
end

# type conversion
context("Conversion 2") do
    @fact isa(convert(Vec3f,v1), Vec3f)  --> true

    @fact isa(convert(Vector{Float64}, v1), Vector{Float64})  --> true
    @fact convert(Vector{Float64}, v1) --> [1.0,2.0,3.0]
end

for T in [UInt, Int, Float32, Float64]
    context("Conversion to Vec{N,$T}") do
        X = map(T, (1,2,3,4,5))

        context("single value conversion") do
            x = X[1]
            for N in 1:4
                @fact convert(Vec{N,T}, x) --> Vec{N,T}(repeated(x,N)...)
            end
        end

        context("conversion from vararg, tuple & array") do
            for N in 1:4
                tup = X[1:N]
                arr = [tup...]
                @fact convert(Vec{N,T}, tup...) --> Vec{N,T}(tup...) "Vec{$N,$T} from vararg"
                @fact convert(Vec{N,T}, tup)    --> Vec{N,T}(tup...) "Vec{$N,$T} from tuple"
                @fact convert(Vec{N,T}, arr)    --> Vec{N,T}(tup...) "Vec{$N,$T} from array"
                @fact convert(Vec, tup...) --> Vec{N,T}(tup...) "Vec from vararg"
                @fact convert(Vec, tup)    --> Vec{N,T}(tup...) "Vec from tuple"
                @fact convert(Vec, arr)    --> Vec{N,T}(tup...) "Vec from array"
            end
        end

        context("conversion from too many args should fail") do
            for N in 1:4
                tup = X[1:N+1]
                arr = [tup...]
                @fact_throws convert(Vec{N,T}, tup...)
                @fact_throws convert(Vec{N,T}, tup)
                @fact_throws convert(Vec{N,T}, arr)
            end
        end

        context("conversion from too few args should fail") do
            for N in 3:5
                tup = X[1:N-1]
                arr = [tup...]
                @fact_throws convert(Vec{N,T}, tup...)
                @fact_throws convert(Vec{N,T}, tup)
                @fact_throws convert(Vec{N,T}, arr)
            end
        end
    end
end

# matrix operations

#typealias Mat1d Matrix1x1{Float64}

zeromat = Mat2d((0.0,0.0),(0.0,0.0))

context("Matrix") do
    @fact map(Float64, zeromat) --> zeromat
    @fact length(Mat2d) --> 4
    @fact length(zeromat) --> 4

    @fact size(Mat2d) --> (2,2)
    @fact size(zeromat) --> (2,2)

    @fact zero(Mat2d) --> zeromat

    for i=1:4, j=1:4
        x1 = rand(i,j)
        @fact @inferred(ctranspose(Mat(x1))) --> Mat(x1')
    end

    v = Vec(1.0,2.0,3.0,4.0)
    r = row(v)
    c = column(v)

    #@fact r' --> c
    #@fact c' --> r

    a = c*r

    b = Mat(
        (1.0,2.0,3.0,4.0),
        (2.0,4.0,6.0,8.0),
        (3.0,6.0,9.0,12.0),
        (4.0,8.0,12.0,16.0)
    )

    x = Mat(
        (1,1,1,),
        (2,2,2,),
        (3,3,3,),
    )
    @fact transpose(x) --> Mat(
        (1,2,3),
        (1,2,3),
        (1,2,3),
    )
    @fact transpose(b) --> b

    @fact length(b) --> 16

    @fact a-->b
    mat30 = Mat(((30.0,),))
    @fact r*c --> mat30

    #@fact row(r, 1) --> v
    #@fact column(c,1) --> v
    #@fact row(r+c',1) --> 2*v
    @fact sum(r) --> sum(v)
    @fact prod(c) --> prod(v)

    @fact eye(Mat3d) --> Mat((1.0,0.0,0.0),
                                (0.0,1.0,0.0),
                                (0.0,0.0,1.0))
    #@fact v*eye(Mat4d)*v --> 30.0
    @fact -r --> -1.0*r
    #@fact diag(diagm(v)) --> v

    # type conversion
    #@fact isa(convert(Matrix1x4{Float32},r),Matrix1x4{Float32})
    jm = rand(4,4)
    im = Mat(jm)
    for i=1:4*2
        @fact jm[i] --> im[i]
    end
    #im = Matrix4x4(jm)
    @fact isa(im, Mat4d)  --> true

    jm2 = convert(Array{Float64,2}, im)
    @fact isa(jm2, Array{Float64,2})  --> true
    @fact jm --> jm2

    #Single valued constructor
    @fact Mat4d(0.0) --> zero(Mat4d)

    a = Vec4d(0)
    b = Vec4d(0,0,0,0)
    @fact a --> b

    v = rand(4)
    m = rand(4,4)
    vfs = Vec(v)
    mfs = Mat(m)
    @fact typeof(vfs) --> Vec4d
    @fact typeof(mfs) --> Mat4d

    # issue #65, wrong
    a = Mat((1,2), (3,4))
    @fact Mat(a) --> a
    b = Mat([1,2,3,4])
    @fact b --> Mat((1,2,3,4))
    @fact b --> Mat([1,2,3,4]'')
end
context("Matrix Math") do
    for i=1:5, j=1:5
        v = rand(j)
        m = rand(i,j)
        m2 = rand(i,j)
        mc = rand(i,j) + im*rand(i,j)
        vfs = Vec(v)
        mfs = Mat(m)
        m2fs = Mat(m2)
        mfsc = Mat(mc)

        vi = randperm(j)
        mi = reshape(randperm(i*j), i, j)
        mi2 = reshape(randperm(i*j), i, j)
        vifs = Vec(vi)
        mifs = Mat(mi)
        mi2fs = Mat(mi2)

        context("Matrix{$i, $j} * Vector{$j}") do
            vm = m * v
            @fact isapprox(@inferred(mfs * vfs), vm)  --> true
            @fact isapprox(@inferred(Matrix(mfs) * vfs), vm)  --> true
            @fact isapprox(@inferred(mfs * Vector(vfs)), vm)  --> true
        end
        context("Matrix{$i, $j} * Matrix{$j, $i}") do
            mm = m * m2'
            @fact isapprox(@inferred(mfs * m2fs'), mm)  --> true
            @fact isapprox(@inferred(Matrix(mfs) * m2fs'), mm)  --> true
            @fact isapprox(@inferred(mfs * Matrix(m2fs')), mm)  --> true
        end
        context("Matrix{$i, $j}*(2I)") do
            mm = m*(2)
            @fact isapprox(@inferred(m*(2I)), mm)  --> true
        end

        # test different element types
        context("Matrix{$i, $j, T} * Vector{$j, U}") do
            vmi = mi * v
            @fact isapprox(@inferred(mifs * vfs), vmi)  --> true
            vmi = m * vi
            @fact isapprox(@inferred(mfs * vifs), vmi)  --> true
            # Custom vector types
            @fact @inferred(eye(Mat{3,3,Float64}) * RGB{Int}(1,2,3)) --> exactly(RGB{Float64}(1,2,3))
        end
        context("Matrix{$i, $j, T} * Matrix{$j, $i, U}") do
            mmi = mi * m2'
            @fact isapprox(@inferred(mifs * m2fs'), mmi)  --> true
            mmi = m * mi2'
            @fact isapprox(@inferred(mfs * mi2fs'), mmi)  --> true
        end

        if i == j
            context("(2*I + I*M)\\v") do
                mm = (2*I+I*m) \ v
                @fact isapprox(@inferred((2*I+I*mfs) \ vfs), mm)  --> true
            end
            context("det(M)") do
                mm = det(m)
                fmm = det(mfs)
                @fact isapprox(fmm, mm)  --> true
            end
            context("trace(M)") do
                mm = trace(m)
                fmm = trace(mfs)
                @fact isapprox(fmm, mm)  --> true
            end
            context("inv(M)") do
                mm = inv(m)
                fmm = inv(mfs)
                @fact isapprox(fmm, mm)  --> true
            end
            context("expm(M)") do
                mm = expm(m)
                fmm = expm(mfs)
                @fact isapprox(fmm, mm)  --> true

                mm = expm(mc)
                fmm = expm(mfsc)
                @fact isapprox(fmm, mm)  --> true
            end
            context("lyap(M,M2*M2')") do
                mm = lyap(m, m2*m2')
                fmm = lyap(mfs, m2fs*m2fs')
                @fact isapprox(fmm, mm) --> true
            end
            context("chol(M2*M2')") do
                mm = full(chol(m2*m2'))
                mm2 = full(chol(map(Mat, m2*m2'))) # Matrix of Mat's
                fmm = chol(m2fs*m2fs')
                @fact isapprox(fmm, mm) --> true
                @fact isapprox(mm, map(first, mm2)) --> true

            end

        else
            context("Matrix{$i, $j} * Matrix{$i, $j}") do
                @fact_throws DimensionMismatch mfs * mfs
            end
        end
        context("transpose M") do
            mm = m'
            fmm = mfs'
            @fact isapprox(fmm, mm)  --> true
        end

        context("ctranspose M") do
            mm = mc'
            fmm = mfsc'
            @fact isapprox(fmm, mm)  --> true
        end
    end
    context("expm(M::Mat{3,3,Float64})") do
        # in practice the precision is eps(), if m has not a triple eigenvalue
        for i in 1:30
            m = (rand(0:1,3,3).*randn(3,3) .+ rand(-3:3,3,3)) # some entries are natural numbers to have higher chance of multiple eigenvalues to trigger all branches
            @fact norm(Matrix(expm(Mat(m))) -  expm(m))/norm(expm(m)) <= 1E-9 --> true
            m = m + m' # symmetric
            @fact norm(Matrix(expm(Mat(m))) -  expm(m))/norm(expm(m)) <= 1E-9 --> true
            m = 1. *rand(-1:1,3,3) # eigenvalues equal with high probability to test worse case
            @fact norm(Matrix(expm(Mat(m))) -  expm(m))/norm(expm(m)) <= 1E-9 --> true
            m = m + m'
            @fact norm(Matrix(expm(Mat(m))) -  expm(m))/norm(expm(m)) <= 1E-9 --> true
        end
    end
    context("expm(M::Mat{3,3, BigFloat})") do
        @fact norm(Matrix(expm(Mat(big([0.0 0.0 1.0; -1.0 1.0 0.0; -1.0 0.0 2.0])))) - big([-0.0 0.0  1.0; -0.5  1.0  -0.5; -1.0  0.0  2.0])*e,1) <  10eps(big(1.)) --> true
    end

    context("Matrix * FixedVectorNoTuple") do
        rgb = rand(3)
        m = rand(3,3)
        rgbfs = RGB(rgb)
        mfs = Mat(m)
        @fact isapprox(mfs * rgbfs, m * rgb) --> true
    end

    context("Outer product  Vec{N} * Mat{1,M}") do
        v1 = Vec(1,2)
        v2 = Vec(1,2,3)
        @fact v1*v2' --> Vector(v1)*Vector(v2)'
    end

    context("Large matrix multiply") do
        M = rand(9,9)
        v = rand(9)
        @fact isapprox(Mat(M)*Vec(v), M*v) --> true
    end
end

ac = rand(3)
aci = randperm(3)
bc = rand(3)

a = rand(4)
b = rand(4)
c = rand(4,4)

d = cross(ac, bc)
di = cross(aci, bc)
d2 = a+b
f = c*a
g = c*b
h = c*f
i = dot(f, a)
j = dot(a, g)
k = abs(f)
l = abs(-f)

acfs = Vec(ac)
acifs = Vec(aci)
bcfs = Vec(bc)

afs = Vec(a)
bfs = Vec(b)
cfs = Mat(c)

dfs = cross(acfs, bcfs)
difs = cross(acifs, bcfs)
d2fs = afs+bfs
ffs = cfs*afs
gfs = cfs*bfs
hfs = cfs*ffs
ifs = dot(ffs, afs)
jfs = dot(afs, gfs)
kfs = abs(ffs)
lfs = abs(-ffs)

context("Meta") do
    sym, expr = FixedSizeArrays.gen_functor(:+, 2)
    @fact typeof(sym) --> Symbol
    @fact typeof(expr) --> Expr
end

context("Vector Math") do
    context("all") do
        @fact isapprox(acfs, ac)  --> true
        @fact isapprox(bcfs, bc)  --> true

        @fact isapprox(afs, a) --> true
        @fact isapprox(bfs, b) --> true
        @fact isapprox(cfs, c) --> true

        @fact isapprox(dfs, d) --> true
        @fact isapprox(difs, di) --> true
        @fact isapprox(d2fs, d2) --> true
        @fact isapprox(ffs, f) --> true
        @fact isapprox(gfs, g) --> true
        @fact isapprox(hfs, h) --> true
        @fact isapprox(ifs, i) --> true
        @fact isapprox(jfs, j) --> true
        @fact isapprox(kfs, k) --> true
        @fact isapprox(lfs, l) --> true
        @fact isapprox(lfs, lfs) --> true
    end
end

context("Equality") do
    @fact Vec{3, Int}(1) --> Vec{3, Float64}(1)
    @fact Vec{2, Int}(1) --> not(Vec{3, Float64}(1))
    @fact Vec(1,2,3) --> Vec(1.0,2.0,3.0)
    @fact Vec(1,2,3) --> not(Vec(1.0,4.0,3.0))
    @fact Vec(1,2,3) --> [1,2,3]
    @fact Mat((1,2),(3,4)) --> Mat((1,2),(3,4))
    @fact one(Mat{4,1, Float32}) --> one(Vec{4, Float32})
    @fact isapprox(Vec(1.0,0.0), Vec(1.0,1e-14)) --> true
end
#=
#don't have this yet
let
    a = rand(16)
    b = Mat4d(a)
    @fact b --> reshape(a, (4,4))
    @fact reshape(a, (4,4)) --> b
    @fact b --> not(reshape(a, (2,8)))
end
=#

const unaryOps = (
    -, ~, conj, abs,
    sin, cos, tan, sinh, cosh, tanh,
    asin, acos, atan, asinh, acosh, atanh,
    sec, csc, cot, asec, acsc, acot,
    sech, csch, coth, asech, acsch, acoth,
    sinc, cosc, cosd, cotd, cscd, secd,
    sind, tand, acosd, acotd, acscd, asecd,
    asind, atand, rad2deg, deg2rad,
    log, log2, log10, log1p, exponent, exp,
    exp2, expm1, cbrt, sqrt, erf,
    erfc, erfcx, erfi, dawson,

    trunc, round, ceil, floor,
    significand, lgamma,
    gamma, lfact, frexp, modf, airy, airyai,
    airyprime, airyaiprime, airybi, airybiprime,
    besselj0, besselj1, bessely0, bessely1,
    eta, zeta, digamma, real, imag
)

# vec-vec and vec-scalar
const binaryOps = (
    .+, .-, .*, ./, .\, /,
    .==, .!=, .<, .<=, .>, .>=, +, -,
    min, max,

    atan2, besselj, bessely, hankelh1, hankelh2,
    besseli, besselk, beta, lbeta
)

context("mapping operators") do
    context("binary: ") do
        test1 = (Vec(1,2,typemax(Int)), Mat((typemin(Int),2,5), (2,3,5), (-2,3,6)), Vec{4, Float32}(0.777))
        test2 = (Vec(1,0,typemax(Int)), Mat((typemin(Int),77,1), (2,typemax(Int),5), (-2,3,6)), Vec{4, Float32}(-23.2929))
        for op in binaryOps
            for i=1:length(test1)
                v1 = test1[i]
                v2 = test2[i]
                context("$op with $v1 and $v2") do
                    try # really bad tests, but better than nothing...
                        if applicable(op, v1[1], v2[1]) && typeof(op(v1[1], v2[1])) == eltype(v1)
                            r = op(v1, v2)
                            for j=1:length(v1)
                                @fact r[j] --> op(v1[j], v2[j])
                            end
                        end
                    end
                end
            end
        end
    end
    context("unary: ") do
        test = (Vec(1,2,typemax(Int)), Mat((typemin(Int),2,5), (2,3,5), (-2,3,6)), Vec{4, Float32}(0.777))
        for op in unaryOps
            for t in test
                context("$op with $t") do
                    try
                        if applicable(op, t[1]) && typeof(op(t[1])) == eltype(t)
                            v = op(t)
                            for i=1:length(v)
                                @fact v[i] --> op(t[i])
                            end
                        end
                    end
                end
            end
        end
    end
end

context("typed round/floor/ceil/trunc") do
    v = Vec((0.8,1.2,-0.3))
    @fact floor(Int, v) --> Vec((0,1,-1))
    @fact ceil( Int, v) --> Vec((1,2,0))
    @fact trunc(Int, v) --> Vec((0,1,0))
    @fact round(Int, v) --> Vec((1,1,0))
end

context("shift, push...") do
    v = Vec(1,2,3)
    p = Point(1,2.)
    @fact @inferred(shift(v)       ) --> Vec(2,3)
    @fact @inferred(shift(p)       ) --> Point(2.)

    @fact @inferred(unshift(v, 42) ) --> Vec(42, 1,2,3)
    @fact @inferred(unshift(v, 42.)) --> Vec(42., 1,2,3)
    @fact @inferred(unshift(p, 42.)) --> Point(42, 1, 2.)

    @fact @inferred(push(v, 42)    ) --> Vec(1,2,3, 42)
    @fact @inferred(push(v, 42.)   ) --> Vec(1,2,3, 42.)
    @fact @inferred(push(p, 42.)   ) --> Point(1,2,42.)

    @fact @inferred(pop(v)      ) --> Vec(1,2)
    @fact @inferred(pop(p)      ) --> Point(1.)

    @fact @inferred(deleteat(v,1)  ) --> Vec(2,3)
    @fact @inferred(deleteat(v,2)  ) --> Vec(1,3)
    @fact @inferred(deleteat(v,3)  ) --> Vec(1,2)
    @fact @inferred(deleteat(p,2)  ) --> Point(1.)

    @fact_throws BoundsError deleteat(v,5)
    @fact_throws BoundsError deleteat(v,-9)

    @fact @inferred(insert(v, 1, 42) ) --> Vec(42,1,2,3)
    @fact @inferred(insert(v, 2, 42) ) --> Vec(1,42,2,3)
    @fact @inferred(insert(v, 3, 42) ) --> Vec(1,2,42,3)
    @fact @inferred(insert(v, 4, 42) ) --> Vec(1,2,3,42)
    @fact @inferred(insert(p, 3, 42.)) --> Point(1,2,42.)

    @fact_throws BoundsError insert(v, 5, 42)
    @fact_throws BoundsError insert(v, 0, 42)
end

context("Base.Test") do
    a = rand(2)
    @fact Base.Test.@test_approx_eq(a, Vec(a)) --> nothing
end

facts("show for subtype") do

    Base.show(io::IO, x::TestType) = print(io, "$(Tuple(x))")  # show for new type

    x = TestType(1, 2)
    @fact string(x) --> "(1,2)"
end

end

end

test()

end

FactCheck.exitstatus()


text:
"""
Functionality for handling the defautl conent of the various Genie files (migrations, models, controllers, etc).
"""
module FileTemplates

using Genie.Inflector

"""
    newtask(module_name::String) :: String

Default content for a new Genie Toolbox task.
"""
function newtask(module_name::String) :: String
  """
  module $module_name

  using Genie, Genie.Toolbox

  \"\"\"
  Description of the task here
  \"\"\"
  function runtask()
    # Build something great
  end

  end
  """
end

"""
    newcontroller(controller_name::String) :: String

Default content for a new Genie controller.
"""
function newcontroller(controller_name::String) :: String
  """
  module $(controller_name)Controller
    # Build something great
  end
  """
end

"""
    newtest(plural_name::String, singular_name::String) :: String

Default content for a new test file.
"""
function newtest(plural_name::String, singular_name::String) :: String
  """
  using Genie, App.$(plural_name)

  ### Your tests here
  @test 1 == 1
  """
end

"""
    appmodule(path::String)

Generates a custom app module when a new app is bootstrapped.
"""
function appmodule(path::String)
  path = replace(path, "-"=>"_") |> strip
  appname = split(path, "/", keepempty = false)[end] |> String |> Inflector.from_underscores

  content = """
  module $appname

  using Genie, Genie.Router, Genie.Renderer, Genie.AppServer

  function main()
    Base.eval(Main, :(const UserApp = $appname))

    include(joinpath("..", "genie.jl"))

    Base.eval(Main, :(const Genie = $appname.Genie))
    Base.eval(Main, :(using Genie))
  end; main()

  end
  """

  (appname, content)
end

end


text:
#TODO: Some kind of generator/iterator/resumable function would be prefered here...
function read_return(episode, i, γ, non_termination_penalty)
    s, a, r = episode[i]
    g = r
    w = 1
    for j = (i+1):length(episode)-1
        w = w * γ
        _, _, r1 = episode[j]
        g += w * r1
    end
    s_fin, _, _ = episode[end]
    if s_fin != TERMINAL_STATE_CODE
        g -= non_termination_penalty
    end
    return s, a, g
end

function update_state_action_values_repository!(R, episode, γ::Real, non_termination_penalty)
    for i = 1:length(episode)-1
        s, a, g = read_return(episode, i, γ, non_termination_penalty)
        ∑r, n = R[s, a]
        R[s, a] = (∑r + g, n + 1)
    end
end

function Q_from_state_action_values_repository!(Q, R)
    for s = 1:size(Q, 1)
        for a = 1:size(Q, 2)
            ∑r, n = R[s, a]
            if n != 0
                Q[s, a] = ∑r / n
            else
                Q[s, a] = NaN
            end
        end
    end
end

"""
    mk_evaluate_policy!(V, Q, mdp, policy, γ; tol, maxiter)

Evaluates (deterministic or stochastic) decision policy `policy` for the given Markov
decision process `mdp` and given discount factor `γ`.

By saying that a policy has been "evaluated", we imply that state values `V` and
state-action values `Q` are computed for this policy. This function evaluates policies in-place,
approximately, in an iterrative manner, until maximal difference between state values in consequtive
iterations is less than the specified `tol`, or until a maximal number of iterations `maxiter`
is performed.

The function returns a 2-tuple containing: the number of iterations performed, and a boolean
indicator of convergence (which is `true` if the vector of state values is evaluated with the
required precission, and `false` otherwise).
"""
function mk_evaluate_policy!(Q::AbstractMatrix{<:Real}, policy_simulator, γ::Real; maxiter = 100, non_termination_penalty = 1000)
    R = fill((0.0, 0), size(Q))
    for _ = 1:maxiter
        s0 = rand(1:size(Q, 1))
        a0 = rand(1:size(Q, 2))
        episode = policy_simulator(s0, a0)
        update_state_action_values_repository!(R, episode, γ, non_termination_penalty)
    end # for: iterations
    Q_from_state_action_values_repository!(Q, R)
end

function mk_update_Q!(Q, ω, policy_simulator, γ; non_termination_penalty = 1000)
    s0 = rand(1:size(Q, 1))
    a0 = rand(1:size(Q, 2))
    episode = policy_simulator(s0, a0)
    Δ = -Inf
    for i = 1:length(episode)-1
        s, a, g = read_return(episode, i, γ, non_termination_penalty)
        new_value = ω * Q[s, a] + (1 - ω) * g
        Δ = max(Δ, abs(new_value - Q[s, a]))
        Q[s, a] = new_value
    end
    return Δ
end


text:
"""
Boolean functions of the type B^(n+1) -> B
with one non-essential variable
"""
struct BooleanFunction{T} <: AbstactBoolean
    m :: T
    dom :: BCube{Int}
    rng :: BCube{Int}
end

"""
Boolean map of type B^n -> B^m
"""
struct BooleanMap{T} <: AbstactBoolean
    m :: T
    dom :: BCube{Int}
    rng :: BCube{Int}
end

const BoolFun = BooleanFunction
const BoolMap = BooleanMap

domain(bf :: AbstactBoolean) = bf.dom
codomain(bf :: AbstactBoolean) = bf.rng

# standard constructor from domain (and range)
BoolFun(bc :: BCube) = BoolFun(zeros(Int, length(bc)), bc, BCube(1))

function BoolMap(dom :: BCube, rng :: BCube)
    bv = BoolVec(0, ndims(rng))
    m = [bv for i in 1 : length(dom)]
    return BoolMap(m, dom, rng)
end

# construct from arities - number
BoolFun(arity :: Int) = BoolFun(BCube(arity))
BoolMap(m :: Int, n :: Int) = BoolMap(BCube(m), BCube(n))

# construct from values
function BoolFun(vals :: Union{NTuple, Vector})
    arity = Int(log2(length(vals)))
    dom = BCube(arity)
    cod = BCube(1)
    return BoolFun(vals, dom, cod)
end

# apply to BoolVec
(f :: AbstactBoolean)(v :: BoolVec) = f.m[value(v) + 1]
# apply to Iterable
function (f :: AbstactBoolean)(it :: Union{NTuple, Vector})
    x = domain(f)(it)
    return f(x)
end
# apply to bunch of numbers
(fs :: AbstactBoolean)(v...) = fs(tuple(v...))

setindex!(f :: AbstactBoolean, v, k :: BoolVec) = (f.m[value(k) + 1] = domain(f)(v))
setindex!(f :: AbstactBoolean, v, k) = (f[domain(f)(k)] = v)
setindex!(f :: AbstactBoolean, v, k...) = (f[tuple(k...)] = v)

function direct(fs :: Union{NTuple{N, T}, Vector{T}}) where T <: BoolFun where N
    #dim_dom = map(x -> ndims(domain(x)), +, fs)
    dom = domain(fs[1])
    dim_dom = ndims(dom)
    dim_rng = length(fs)
    fprod = BoolMap(dim_dom, dim_rng)
    for x in dom
        y = Tuple(f(x) for f in fs)
        fprod[x] = y
    end
    return fprod
end

mat_repr(f :: BoolFun) = f.m
function mat_repr(fs :: BoolMap)
    rows = length(domain(fs))
    cols = ndims(codomain(fs))
    res = zeros(Int, rows, cols)
    row = zeros(Int, cols)
    for (i, x) in enumerate(domain(fs))
        fx = fs(x)
        fill_bits!(row, value(fx))
        res[i, :] .= row
    end
    return res
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
using Gen
import Random

#########
# model #
#########

@gen function datum(x::Float64, @ad(inlier_std), @ad(outlier_std), @ad(slope), @ad(intercept))
    is_outlier = @addr(bernoulli(0.5), :z)
    std = is_outlier ? inlier_std : outlier_std
    y = @addr(normal(x * slope + intercept, std), :y)
    return y
end

data = plate(datum)

@gen function model(xs::Vector{Float64})
    n = length(xs)
    inlier_std = exp(@addr(normal(0, 2), :inlier_std))
    outlier_std = exp(@addr(normal(0, 2), :outlier_std))
    slope = @addr(normal(0, 2), :slope)
    intercept = @addr(normal(0, 2), :intercept)
    if all([@change(addr) == NoChange() for addr in [:slope, :intercept, :inlier_std, :outlier_std]])
        change = NoChange()
    else
        change = nothing
    end
    ys = @addr(data(xs, fill(inlier_std, n), fill(outlier_std, n), fill(slope, n), fill(intercept, n)), :data, change)
    return ys
end

#######################
# inference operators #
#######################

@gen function is_outlier_proposal(prev, i::Int)
    prev_z = get_assignment(prev)[:data => i => :z]
    @addr(bernoulli(prev_z ? 0.0 : 1.0), :data => i => :z)
end

@gen function observer(ys::Vector{Float64})
    for (i, y) in enumerate(ys)
        @addr(dirac(y), :data => i => :y)
    end
end

Gen.load_generated_functions()

#####################
# generate data set #
#####################

Random.seed!(1)

prob_outlier = 0.5
true_inlier_noise = 0.5
true_outlier_noise = 5.0
true_slope = -1
true_intercept = 2
xs = collect(range(-5, stop=5, length=200))
ys = Float64[]
for (i, x) in enumerate(xs)
    if rand() < prob_outlier
        y = true_slope * x + true_intercept + randn() * true_inlier_noise
    else
        y = true_slope * x + true_intercept + randn() * true_outlier_noise
    end
    push!(ys, y)
end

##################
# run experiment #
##################

selection = DynamicAddressSet()
push_leaf_node!(selection, :slope)
push_leaf_node!(selection, :intercept)
push_leaf_node!(selection, :inlier_std)
push_leaf_node!(selection, :outlier_std)

function do_inference(n)
    observations = get_assignment(simulate(observer, (ys,)))

    # initial trace
    (trace, _) = generate(model, (xs,), observations)

    for i=1:n
        trace = mala(model, selection, trace, 0.0001)
        trace = hmc(model, selection, trace)

        # step on the outliers
        for j=1:length(xs)
            trace = mh(model, is_outlier_proposal, (j,), trace)
        end

        score = get_call_record(trace).score

        # print
        assignment = get_assignment(trace)
        slope = assignment[:slope]
        intercept = assignment[:intercept]
        inlier_std = assignment[:inlier_std]
        outlier_std = assignment[:outlier_std]
        println("score: $score, slope: $slope, intercept: $intercept, inlier_std: $inlier_std, outlier_std: $outlier_std")
    end
end

@time do_inference(1000)
@time do_inference(100)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


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

using FEMBasis
using FEMBasis: create_basis
using Test

X = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
basis = [:(1.0 - u - v), :(1.0u), :(1.0v)]
dbasis = [-1.0 -1.0; 1.0 0.0 ; 0.0 1.0]'

code = create_basis(:TestTriangle1, "test triangle 1", X, basis, dbasis)
eval(code)
N = zeros(length(TestTriangle1, 2))
X = get_reference_element_coordinates(TestTriangle1)
for i=1:length(TestTriangle1)
    eval_basis!(TestTriangle1, N, X[i])
    N_expected = zeros(length(TestTriangle1))
    N_expected[i] = 1.0
    @test isapprox(N, N_expected)
end
dN = [Vec{zeros(size(TestTriangle1)...)
eval_dbasis!(TestTriangle1, dN, X[1])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle1, dN, X[2])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle1, dN, X[3])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])

code = create_basis(:TestTriangle2, "test triangle 2", X, basis)
eval(code)
N = zeros(1,size(TestTriangle2, 2))
X = get_reference_element_coordinates(TestTriangle2)
for i=1:length(TestTriangle2)
    eval_basis!(TestTriangle2, N, X[i])
    N_expected = zeros(1, length(TestTriangle2))
    N_expected[i] = 1.0
    @test isapprox(N, N_expected)
end
dN = zeros(size(TestTriangle2)...)
eval_dbasis!(TestTriangle2, dN, X[1])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle2, dN, X[2])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle2, dN, X[3])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])

p = :(1 + u + v)
code = create_basis(:TestTriangle3, "test triangle 3", X, p)
eval(code)
N = zeros(1, size(TestTriangle3, 2))
X = get_reference_element_coordinates(TestTriangle3)
for i=1:length(TestTriangle3)
    eval_basis!(TestTriangle3, N, X[i])
    N_expected = zeros(1, length(TestTriangle3))
    N_expected[i] = 1.0
    @test isapprox(N, N_expected)
end
dN = zeros(size(TestTriangle3)...)
eval_dbasis!(TestTriangle3, dN, X[1])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle3, dN, X[2])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
eval_dbasis!(TestTriangle3, dN, X[3])
@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])


text:
using CoordinateTransformations

@testset "Indexing" begin
    x = GeoArray(rand(5, 5, 5))
    @test length(x[1]) == 1
    @test length(x[1, 2]) == 5
    @test size(x[1:3, 1:3]) == (3, 3, 5)
    @test size(x[1:3, 1:3, 1:3]) == (3, 3, 3)

    @inferred size(x)
    @inferred size(x[1, 1, 1])

    x = GeoArray(rand(10, 10, 2))
    xs = x[1:2:end, 1:2:end]
    @test size(xs) == (5, 5, 2)
    @test bbox(xs) == bbox(x)
end

@testset "Concrete" begin
    ga = GeoArray(rand(5, 5))
    @test isconcretetype(typeof(ga))
    @test isconcretetype(typeof(ga.A))
    @test isconcretetype(typeof(ga.f))
    @test isconcretetype(typeof(ga.crs))
end

@testset "Reading rasters" begin
    ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
    @test bbox(ga) == (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
    @inferred bbox(ga)
    @test bboxes(ga)[1] == (min_x=440720.0, max_x=440780.0, min_y=3.75126e6, max_y=3.75132e6)
    @test bboxes(ga)[end] == (min_x=446660.0, max_x=446720.0, min_y=3.74532e6, max_y=3.74538e6)
    @inferred bboxes(ga)
end

@testset "Coords" begin
    straight = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.0; 0.0 -1.0], [375000.03, 380000.03]), "")
    rot = GeoArray(rand(5, 5, 1), AffineMap([1.0 0.5; 0.1 1.0], [0.0, 0.0]), "")
    @inferred coords(straight, :x, GeoArrays.Vertex())
    @test coords(straight, :x, GeoArrays.Vertex()) == collect(375000.03:1:375005.03)
    @test coords(straight, :y, GeoArrays.Vertex()) == collect(380000.03:-1:379995.03)
    @test_throws ErrorException coords(straight, :z)
    @test_throws ErrorException coords(rot, :x)
    @test_throws ErrorException coords(rot, :y)
end

@testset "GeoArray constructors" begin
    x, y = range(4, stop=8.0, length=10), range(0, stop=1, length=9)
    ga2 = GeoArray(rand(10, 9), x, y)
    @inferred GeoArray(rand(10, 9), x, y)
    ga2 = GeoArray(rand(10, 9), x, y, "")
    ga3 = GeoArray(rand(10, 9, 8), x, y)
    ga3 = GeoArray(rand(10, 9, 8), x, y, "")
    for i in 1:length(x), j in 1:length(y)
        @test GeoArrays.coords(ga2, [i, j]) ≈ [x[i], y[j]]
    end
    for i in 1:length(x), j in 1:length(y)
        @test GeoArrays.coords(ga3, [i, j]) ≈ [x[i], y[j]]
    end
    x, y = range(4, stop=8.0, length=11), range(0, stop=1, length=9)
    @test_throws ErrorException GeoArray(rand(10, 9), x, y)
end

@testset "Conversions" begin
    ga = GeoArray(rand(1:32000, 5, 5))
    @inferred GeoArrays.write!(joinpath(testdatadir, "test_conversion.tif"), ga)
    @inferred GeoArrays.write(joinpath(testdatadir, "test_conversion.tif"), ga)

    ga = GeoArray(rand(Bool, 5, 5))
    @test_throws ErrorException GeoArrays.write!(joinpath(testdatadir, "test_conversion.tif"), ga)
end

@testset "Similar" begin
    ga = GeoArrays.read(joinpath(testdatadir, remotefiles[end-1]))
    gg = ga[250:end-2, 250:end-250]
    @test gg[1, 1] == ga[250, 250]
    @test coords(gg, [1, 1]) == coords(ga, [250, 250])
    GeoArrays.write!("test.tif", gg)
end

@testset "Indexing" begin
    ga = GeoArray(rand(10, 10))
    @inferred ga[Float32(1.0), Float32(2.0)]
    i, j = 2, 5
    x, y = coords(ga, (i, j))
    ii, jj = indices(ga, (x, y))
    @inferred indices(ga, (x, y))
    @test ii == i
    @test jj == j
end