text:
using FileIO
using Test

# Stub readers---these might bork any existing readers, so don't
# run these tests while doing other things!
module TestLoadSave
import FileIO: File, @format_str
load(file::File{format"PBMText"})   = "PBMText"
load(file::File{format"PBMBinary"}) = "PBMBinary"
load(file::File{format"JLD"})       = "JLD"
load(file::File{format"GZIP"})       = "GZIP"
end
module TestLoadSave2
import FileIO: File, @format_str
fileio_load(file::File{format"HDF5"})      = "HDF5"
end

sym2loader = copy(FileIO.sym2loader)
sym2saver = copy(FileIO.sym2saver)

try
    empty!(FileIO.sym2loader)
    empty!(FileIO.sym2saver)
    file_dir = joinpath(dirname(@__FILE__), "files")
    @testset "Load" begin

        add_loader(format"PBMText", :TestLoadSave)
        add_loader(format"PBMBinary", :TestLoadSave)
        add_loader(format"HDF5", :TestLoadSave2)
        add_loader(format"JLD", :TestLoadSave)
        add_loader(format"GZIP", :TestLoadSave)

        @test load(joinpath(file_dir,"file1.pbm")) == "PBMText"
        @test load(joinpath(file_dir,"file2.pbm")) == "PBMBinary"

        # Regular HDF5 file with magic bytes starting at position 0
        @test load(joinpath(file_dir,"file1.h5")) == "HDF5"
        # This one is actually a JLD file saved with an .h5 extension,
        # and the JLD magic bytes edited to prevent it from being recognized
        # as JLD.
        # JLD files are also HDF5 files, so this should be recognized as
        # HDF5. However, what makes this more interesting is that the
        # magic bytes start at position 512.
        @test load(joinpath(file_dir,"file2.h5")) == "HDF5"
        # JLD file saved with .jld extension
        @test load(joinpath(file_dir,"file.jld")) == "JLD"
        # GZIP file saved with .gz extension
        @test load(joinpath(file_dir,"file.csv.gz")) == "GZIP"
        @test_throws Exception load("missing.fmt")
    end
finally
    merge!(FileIO.sym2loader, sym2loader)
    merge!(FileIO.sym2saver, sym2saver)
end

# A tiny but complete example
# DUMMY format is:
#  - n, a single Int64
#  - a vector of length n of UInt8s

add_format(format"DUMMY", b"DUMMY", ".dmy")

module Dummy

using FileIO

mutable struct DummyReader{IOtype}
    stream::IOtype
    ownstream::Bool
    bytesleft::Int64
end

function DummyReader(stream, ownstream)
    read(stream, 5) == magic(format"DUMMY") || error("wrong magic bytes")
    DummyReader(stream, ownstream, read(stream, Int64))
end

function Base.read(stream::DummyReader, n=stream.bytesleft)
    toread = min(n, stream.bytesleft)
    buf = read(stream.stream, toread)
    stream.bytesleft -= length(buf)
    buf
end

Base.eof(stream::DummyReader) = stream.bytesleft == 0 || eof(stream.stream)
Base.close(stream::DummyReader) = stream.ownstream && close(stream.stream)

mutable struct DummyWriter{IOtype}
    stream::IOtype
    ownstream::Bool
    headerpos::Int64
    byteswritten::Int
end

function DummyWriter(stream, ownstream)
    write(stream, magic(format"DUMMY"))  # Write the magic bytes
    # store the position where we'll need to write the length
    pos = position(stream)
    # write a dummy length value
    write(stream, 0xffffffffffffffff)
    DummyWriter(stream, ownstream, pos, 0)
end

function Base.write(stream::DummyWriter, data)
    udata = convert(Vector{UInt8}, data)
    n = write(stream.stream, udata)
    stream.byteswritten += n

    n
end

function Base.close(stream::DummyWriter)
    here = position(stream.stream)
    # go back and write the header
    seek(stream.stream, stream.headerpos)
    write(stream.stream, convert(Int64, stream.byteswritten))
    seek(stream.stream, here)
    stream.ownstream && close(stream.stream)

    nothing
end

loadstreaming(s::Stream{format"DUMMY"}) = DummyReader(s, false)
loadstreaming(file::File{format"DUMMY"}) = DummyReader(open(file), true)
savestreaming(s::Stream{format"DUMMY"}) = DummyWriter(s, false)
savestreaming(file::File{format"DUMMY"}) = DummyWriter(open(file, "w"), true)

# we could implement `load` and `save` in terms of their streaming versions
function load(file::File{format"DUMMY"}; extra=UInt8[])
    open(file) do s
        load(s; extra=extra)
    end
end

function metadata(file::File{format"DUMMY"})
    s = open(file)
    skipmagic(s)
    n = read(s, Int64)
    close(s)
    return n
end

function load(s::Stream{format"DUMMY"}; extra=UInt8[])
    skipmagic(s)
    n = read(s, Int64)
    out = Vector{UInt8}(undef, n)
    read!(s, out)
    # verify that the extradata is as expected. This is just to test that
    # the keyword arguments are handled properly in loading and saving
    extradata = read(s, length(extra))
    if extradata != extra
        throw(ErrorException("Got extra data $extradata instead of $extra"))
    end
    close(s)
    out
end

function save(file::File{format"DUMMY"}, data; extra=UInt8[])
    open(file, "w") do s
        save(s, data; extra=extra)
    end
end

function save(s::Stream{format"DUMMY"}, data; extra=UInt8[])
    write(s, magic(format"DUMMY"))  # Write the magic bytes
    write(s, convert(Int64, length(data)))
    udata = convert(Vector{UInt8}, data)
    write(s, udata)
    write(s, extra)
end

add_loader(format"DUMMY", :Dummy)
add_saver(format"DUMMY", :Dummy)

end # module Dummy

@testset "Save" begin
    a = [0x01,0x02,0x03]
    fn = string(tempname(), ".dmy")
    save(fn, a)
    @test metadata(fn) == 3

    # Test for absolute paths
    cd(dirname(fn)) do
        fnrel = basename(fn)
        f = query(fnrel)
        @test isabspath(filename(f))
        @test endswith(filename(f),fn) # TravisOSX prepends "/private"
        f = File(format"DUMMY", fnrel)
        @test !(isabspath(filename(f)))
        open(f) do s
            @test isabspath(filename(s))
            @test endswith(filename(s), fn)
        end
    end

    fn2 = string(tempname(), ".dmy")
    a |> save(fn2)

    # Test for absolute paths
    cd(dirname(fn2)) do
        fnrel = basename(fn2)
        f = query(fnrel)
        @test isabspath(filename(f))
        @test endswith(filename(f),fn2) # TravisOSX prepends "/private"
        f = File(format"DUMMY", fnrel)
        @test !(isabspath(filename(f)))
        open(f) do s
            @test isabspath(filename(s))
            @test endswith(filename(s),fn2)
        end
    end
    rm(fn2)

    # Test IO
    b = load(query(fn))
    @test a == b

    b = open(query(fn)) do s
        load(s)
    end
    @test a == b

    # low-level I/O test
    open(query(fn)) do s
        @test position(s) == 0
        skipmagic(s)
        @test position(s) == length(magic(format"DUMMY"))
        seek(s, 1)
        @test position(s) == 1
        seekstart(s)
        @test position(s) == 0
        seekend(s)
        @test eof(s)
        skip(s, -position(s)+1)
        @test position(s) == 1
        @test isreadonly(s)
        @test isopen(s)
        @test read(s,2) == b"UM"
    end
    rm(fn)

    # streaming I/O with filenames
    fn = string(tempname(), ".dmy")
    save(fn, a)
    loadstreaming(fn) do reader
        @test read(reader) == a
    end
    rm(fn)
    savestreaming(fn) do writer
        write(writer, a)
    end
    @test load(fn) == a
    rm(fn)

    # force format
    fn = string(tempname(), ".dmy")
    savestreaming(format"DUMMY", fn) do writer
        write(writer, a)
    end
    @test load(fn) == a
    rm(fn)

    # streaming I/O with streams
    save(fn, a)
    open(fn) do io
        loadstreaming(io) do reader
            @test read(reader) == a
        end
        @test isopen(io)
    end
    rm(fn)
    open(fn, "w") do io
        savestreaming(format"DUMMY", io) do writer
            write(writer, a)
        end
        @test isopen(io)
    end
    @test load(fn) == a
    rm(fn)

    # test keyword arguments

    a = [0x01,0x02,0x03]
    fn = string(tempname(), ".dmy")
    save(fn, a; extra=[0x42, 0x43])
    # the loader verifies that the extra data was written properly
    load(fn; extra=[0x42, 0x43])

    @test_throws Exception save("missing.fmt",5)
end

del_format(format"DUMMY")

# PPM/PBM can be either binary or text. Test that the defaults work,
# and that we can force a choice.
module AmbigExt
import FileIO: File, @format_str, Stream, stream, skipmagic

load(f::File{format"AmbigExt1"}) = open(f) do io
    skipmagic(io)
    read(stream(io), String)
end
load(f::File{format"AmbigExt2"}) = open(f) do io
    skipmagic(io)
    read(stream(io), String)
end

fileio_save(f::File{format"AmbigExt1"}, testdata) = open(f, "w") do io
    s = stream(io)
    print(s, "ambigext1")
    print(s, testdata)
end
fileio_save(f::File{format"AmbigExt2"}, testdata) = open(f, "w") do io
    s = stream(io)
    print(s, "ambigext2")
    print(s, testdata)
end
end

@testset "Ambiguous extension" begin
    add_format(format"AmbigExt1", "ambigext1", ".aext", [:AmbigExt])
    add_format(format"AmbigExt2", "ambigext2", ".aext", [:AmbigExt])
    A = "this is a test"
    fn = string(tempname(), ".aext")
    # Test the forced version first: we wouldn't want some method in Netpbm
    # coming to the rescue here, we want to rely on FileIO's logic.
    # `save(fn, A)` will load Netpbm, which could conceivably mask a failure
    # in the next line.
    save(format"AmbigExt2", fn, A)

    B = load(fn)
    @test B == A
    @test typeof(query(fn)) == File{format"AmbigExt2"}
    rm(fn)

    save(fn, A)
    B = load(fn)
    @test B == A
    @test typeof(query(fn)) == File{format"AmbigExt1"}

    rm(fn)
    del_format(format"AmbigExt1")
    del_format(format"AmbigExt2")
end

@testset "Absent file" begin
    @test_throws SystemError load("nonexistent.oops")
end


text:
using .NorahPhd
using Polyhedra
import CairoMakie as CMakie
include("./Example2.jl")
c1 = [1; 0]
c = [1; 1]
z1 = [1; 0]
z2 = [0; 1]
C = [c1 c]
Z = [z1 z2]
# v0 = [0; 0]
# w, u = D2(v0)
# P1(z1)
# D2(v0)
Xbar, Tbar, V, p = primal(Γ, P1, P2, D2, C, Z, ϵ = 0.01)
# print(p)
m = Polyhedra.Mesh(p)
CMakie.mesh(m, color = :blue)
CMakie.wireframe(m)

ext: jl


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.Orientations
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Interpolation
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics: total_energy, air_density
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using Distributions: Uniform
using LinearAlgebra
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Planet:
    R_d, day, grav, cp_d, planet_radius, Omega, kappa_d, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import CLIMAParameters
CLIMAParameters.Planet.Omega(::EarthParameterSet) = 0.0
CLIMAParameters.Planet.planet_radius(::EarthParameterSet) = 6.371e6 / 125.0
CLIMAParameters.Planet.MSLP(::EarthParameterSet) = 1e5

function init_isothermal_zonal_flow!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)
    param_set = parameter_set(bl)
    φ = latitude(bl.orientation, aux)
    z = altitude(bl.orientation, param_set, aux)

    _grav::FT = grav(param_set)
    _a::FT = planet_radius(param_set)
    _R_d::FT = R_d(param_set)
    _MSLP::FT = MSLP(param_set)

    u₀ = FT(20)
    T₀ = FT(300)

    shallow_atmos = false
    if shallow_atmos
        f1 = z
        f2 = FT(0)
        shear = FT(1)
    else
        f1 = z
        f2 = z / _a + z^2 / (2 * _a^2)
        shear = 1 + z / _a
    end

    u_sphere = SVector{3, FT}(u₀ * shear * cos(φ), 0, 0)
    u_init = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    prefac = u₀^2 / (_R_d * T₀)
    fac1 = prefac * f2 * cos(φ)^2
    fac2 = prefac * sin(φ)^2 / 2
    fac3 = _grav * f1 / (_R_d * T₀)
    exparg = fac1 - fac2 - fac3
    p = _MSLP * exp(exparg)

    ρ = air_density(param_set, T₀, p)

    e_pot = gravitational_potential(bl.orientation, aux)
    e_kin = u_init' * u_init / 2

    state.ρ = ρ
    state.ρu = ρ * u_init
    state.energy.ρe = ρ * total_energy(param_set, e_kin, e_pot, T₀)
end

function config_isothermal_zonal_flow(
    FT,
    poly_order,
    cutoff_order,
    resolution,
    ref_state,
)
    # Set up a reference state for linearization of equations

    domain_height = FT(10e3)

    # Set up the atmosphere model
    exp_name = "IsothermalZonalFlow"

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        moisture = DryModel(),
    )
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_isothermal_zonal_flow!,
        source = (Gravity(),),
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_isothermal_zonal_flow!;
        model = model,
        numerical_flux_first_order = RoeNumericalFlux(),
        Ncutoff = cutoff_order,
    )

    return config
end

function config_diagnostics(FT, driver_config)
    interval = "40000steps" # chosen to allow a single diagnostics collection

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(10), FT(10), FT(100)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 5                           # discontinuous Galerkin polynomial order
    cutoff_order = 4
    n_horz = 10                              # horizontal element number
    n_vert = 5                               # vertical element number
    timestart = FT(0)                        # start time (s)
    timeend = FT(3600)                       # end time (s)

    # set up the reference state for linearization of equations
    temp_profile_ref = IsothermalProfile(param_set, FT(300))
    ref_state = HydrostaticState(temp_profile_ref)

    # Set up driver configuration
    driver_config = config_isothermal_zonal_flow(
        FT,
        poly_order,
        cutoff_order,
        (n_horz, n_vert),
        ref_state,
    )

    # Set up experiment
    ode_solver_type = ClimateMachine.IMEXSolverType(
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
        split_explicit_implicit = false,
        discrete_splitting = true,
    )
    CFL = FT(0.4)
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        init_on_cpu = true,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
    )

    # save the initial condition for testing
    Q0 = copy(solver_config.Q)

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    filterorder = 64
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )

    relative_error = norm(solver_config.Q .- Q0) / norm(Q0)
    @info "Relative error = $relative_error"
    @test relative_error < 1e-7

end

main()


text:
δ(i, j) = Int64(==(i, j))

diagm(x::Vector) = [δ(i, j)*x[i] for i in 1:length(x), j in 1:length(x)]

I(n) = [δ(i, j) for i = 1:n, j = 1:n]

vtanh(x::Vector) = tanh.(x)

dtanh(x::Vector) = I(length(x)) - diagm([i^2 for i in x])

relu(z::Vector) = max.(0, z)

leasSquareCost(o, t) = sum((o - t).^2)/2
dLeastSquareCost(o, t) = o-t

softmaxLikelihood(o, t) = -log(o[argmax(t)])
dSoftmaxLikelihood(o, t) = o - t

function drelu(z::Vector)
  dvec = Int64.(>(0).(z))
  return diagm(dvec[:, 1])
end

function softmax(z::Vector)
  shift = maximum(z)

  z = exp.(z .- shift)
  total = sum(z)
  return z ./ total
end

function dsoftmax(z::Vector)
  n = length(z)

  w = softmax(z)
  e = ones(n, 1)

  return w * e' .* (I(n) - e * w')
end

mean(x) = sum(x)/length(x)

activations = Dict(:vtanh => (vtanh, dtanh), :relu => (relu, drelu), :softmax => (softmax, dsoftmax))

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition


text:
begin
    using CSV
    using Revise
    using BNPRegressionGGA2021
    using DataFrames
    using DataStructures
    using Distributions
    using Random
    using RCall
    using Statistics
    using StatsBase
    using LinearAlgebra
    const BNP = BNPRegressionGGA2021
end;

R"""
library(dplyr)
library(leaps)

# Simulate a dataset
N <- 50
K <- 5
df <-
    data.frame(
        x1 = rnorm(N),
        x2 = runif(N) <= 0.5,
        x3 = runif(N) <= 0.5,
        x4 = runif(N) <= 0.5,
        x5 = runif(N) <= 0.5
    ) |>
    dplyr::mutate(y = sn::rsn(xi = x1, omega = 1.5, alpha = 4))

# Fit the model
fitted_model <- leaps::regsubsets(y ~ ., data = df, nvmax = 5)
fit_summary <- summary(fitted_model, matrix.logical = TRUE)
kstar <- which.min(fit_summary$bic)
best_gamma <- fit_summary$which[kstar, ]
"""

# Helper functions -------------------------------------------------------------

function run_experiment_bnp(dy; N0, Nrep, Niter, id, filename)
    Random.seed!(1) # seed
    samples = [generate_sample(dy; N0, Nrep) for _ in 1:Niter] # data
    gammas = [get_gamma_bnp(samples[id]) for id in 1:Niter] # map gamma
    reduced_gammas = reduce_gammas(gammas) # reduced df
    reduced_gammas[!, :method] .= "bnp"
    reduced_gammas[!, :id] .= id
    reduced_gammas[!, :N0] .= N0
    CSV.write(filename, reduced_gammas) # csv file
end

function run_experiment_freq(dy; N0, Nrep, Niter, id, filename)
    Random.seed!(1) # seed
    samples = [generate_sample(dy; N0, Nrep) for _ in 1:Niter] # data
    gammas = [get_gamma_freq(samples[id]) for id in 1:Niter] # map gamma
    reduced_gammas = reduce_gammas(gammas) # reduced df
    reduced_gammas[!, :method] .= "freq"
    reduced_gammas[!, :id] .= id
    reduced_gammas[!, :N0] .= N0
    CSV.write(filename, reduced_gammas) # csv file
end

function generate_sample(dy; N0 = 50, Nrep = 10)
    # Simulate the data
    x0 = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    X0d = rand([0, 1], length(x0), 4)
    y0 = @. rand(dy(x0, X0d[:, 1]))

    # Generate the grid
    N1 = 0
    x1 = LinRange(-1, 1, N1)
    y1 = LinRange(-3, 3, 50) |> collect

    # Expand the grid
    grid = rcopy(R"""expand.grid(y1 = $y1, x1 = $x1)""")
    y1 = grid[!, :y1]
    x1 = grid[!, :x1]

    # Expand x0-x1 using splines
    R"""
    x0min <- min($x0)
    x0max <- max($x0)
    X0 <- splines::bs($x0, df = 6, Boundary.knots = c(x0min, x0max))
    """
    @rget X0

    # Standardise responses and predictors
    mX0, sX0 = mean_and_std(X0, 1)
    for col in 1:size(X0, 2)
        if sX0[col] > 0
            X0[:, col] = (X0[:, col] .- mX0[col]) ./ sX0[col]
        end
    end
    my0, sy0 = mean_and_std(y0, 1)
    y0 = (y0 .- my0) ./ sy0

    # Add a constant term to the design matrices
    X0d = ((X0d .- mean(X0d, dims = 1)) ./ std(X0d, dims = 1))
    X0 = [ones(size(X0, 1)) X0 X0d]
    X1 = zeros(0, size(X0, 2))

    # Generate the true (standardized) density over the grid
    f1 = pdf.(dy.(x1), my0 .+ y1 .* sy0) .* sy0

    # Set the mapping
    mapping = [[1], collect(2:7), [8], [9], [10], [11]]

    # Return the preprocessed results
    return y0, x0, y1, x1, X0, X1, f1, mapping
end;

# Get the MAP estimator of gamma using BNP
function get_gamma_bnp(sample)
    y0, x0, y1, x1, X0, X1, f1, mapping = sample
    smpl = BNP.DGSBPNormalDependent(; y0, X0, y1, X1, mapping)
    _, _, chaing = BNP.sample!(smpl; mcmcsize = 4000)
    gamma = begin
        chaing |>
        DataStructures.counter |>
        collect |>
        x -> sort(x, by = x -> x[2], rev = true) |>
        x -> getindex(x, 1) |>
        x -> getindex(x, 1)
    end
    return gamma
end

# Estimate gamma using a frequentist alternative
function get_gamma_freq(sample)
    y0, x0, y1, x1, X0, X1, f1, mapping = sample
    R"""
    library(dplyr)
    data = data.frame(
        y = $y0,
        x2 = $x0,
        x3 = $X0[, 8],
        x4 = $X0[, 9],
        x5 = $X0[, 10],
        x6 = $X0[, 11]
    )
    fitted_model <- leaps::regsubsets(y ~ ., data = data, nvmax = 5)
    fitsum <- summary(fitted_model, matrix.logical = TRUE)
    kstar <- which.min(fitsum$bic)
    best_gamma <- fitsum$which[kstar, ]
    """
    @rget best_gamma
    return best_gamma
end

# Reduce the MAP estimators after their estimation
function reduce_gammas(gammas)
    names = "g" .* string.(1:6)
    data = reduce(vcat, gammas')
    DataFrame(data, names)
end

# Experiment 1, N0 = 50 (bnp)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-bnp-50.csv"
    run_experiment_bnp(dy; N0 = 50, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 50 (freq)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-freq-50.csv"
    run_experiment_freq(dy; N0 = 50, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 100 (bnp)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-bnp-100.csv"
    run_experiment_bnp(dy; N0 = 100, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 100 (freq)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-freq-100.csv"
    run_experiment_freq(dy; N0 = 100, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 200 (bnp)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-bnp-200.csv"
    run_experiment_bnp(dy; N0 = 200, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 1, N0 = 200 (freq)
begin
    dy(x1, x2) = Normal(x1 + x2, 1)
    filename = "data/final/regression-gamma-normal-freq-200.csv"
    run_experiment_freq(dy; N0 = 200, Nrep = 5, Niter = 100, id = 1, filename)
end

# Experiment 2, N0 = 50 (bnp)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-bnp-50.csv"
    run_experiment_bnp(dy; N0 = 50, Nrep = 5, Niter = 100, id = 2, filename)
end

# Experiment 2, N0 = 50 (freq)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-freq-50.csv"
    run_experiment_freq(dy; N0 = 50, Nrep = 5, Niter = 100, id = 2, filename)
end

# Experiment 2, N0 = 100 (bnp)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-bnp-100.csv"
    run_experiment_bnp(dy; N0 = 100, Nrep = 5, Niter = 100, id = 2, filename)
end

# Experiment 2, N0 = 100 (freq)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-freq-100.csv"
    run_experiment_freq(dy; N0 = 100, Nrep = 5, Niter = 100, id = 2, filename)
end

# Experiment 2, N0 = 200 (bnp)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-bnp-200.csv"
    run_experiment_bnp(dy; N0 = 200, Nrep = 5, Niter = 100, id = 2, filename)
end

# Experiment 2, N0 = 200 (freq)
begin
    dy(x1, x2) = SkewNormal(x1 + x2, 1.5, 4)
    filename = "data/final/regression-gamma-skewnormal-freq-200.csv"
    run_experiment_freq(dy; N0 = 200, Nrep = 5, Niter = 100, id = 2, filename)
end

# Summary of the results
R"""
data <-
    list.files(
        path = "data/final",
        pattern = "regression-gamma*",
        full.names = TRUE
    ) |>
    purrr::map(readr::read_csv, show_col_types = FALSE) |>
    dplyr::bind_rows() |>
    dplyr::group_by(method, id, N0, g1, g2, g3, g4, g5, g6) |>
    dplyr::count(name = "frequency") |>
    dplyr::ungroup() |>
    dplyr::mutate(gamma = paste0(g1, g2, g3, g4, g5, g6)) |>
    dplyr::select(gamma, frequency, id, method, N = N0)
p <-
    data |>
    dplyr::mutate(
        method = method |>
            dplyr::recode(bnp = "BNP", freq = "BSS (Best Subset Selection)"),
        distribution = id |>
        dplyr::recode(`1` = "Normal", `2` = "Skew Normal"),
        N = factor(N, ordered = TRUE)
    ) |>
    ggplot2::ggplot(ggplot2::aes(y = gamma, x = N, fill = frequency)) +
    ggplot2::geom_tile(color = "grey90") +
    ggplot2::geom_text(ggplot2::aes(label = frequency)) +
    ggplot2::facet_grid(
        distribution ~ method,
        labeller = ggplot2::labeller(
            distribution = ggplot2::label_both,
            method = ggplot2::label_both
        )
    ) +
    ggplot2::theme_classic() +
    ggplot2::theme(legend.position = "top") +
    ggplot2::scale_fill_gradient(low = "white", high = "black") +
    ggplot2::labs(
        x = "N (sample size)",
        y = "gamma (hypothesis vector), as a string of 0s and 1s"
    )
ggplot2::ggsave("figures/final/heatmap-regression-zeta-1div4.png")
"""


text:
# tensor.jl
#
# Tensor provides a dense implementation of an AbstractTensor type without any
# symmetry assumptions, i.e. it describes tensors living in the full tensor
# product space of its index spaces.

#++++++++++++++
# Tensor type:
#++++++++++++++
# Type definition and constructors:
#-----------------------------------
immutable Tensor{S,T,N} <: AbstractTensor{S,ProductSpace,T,N}
    data::Array{T,N}
    space::ProductSpace{S,N}
    function Tensor(data::Array{T},space::ProductSpace{S,N})
        if length(data)!=dim(space)
            throw(DimensionMismatch("data not of right size"))
        end
        if promote_type(T,eltype(S))!=eltype(S)
            warn("For a tensor in $(space), the entries should not be of type $(T)")
        end
        return new(reshape(data,map(dim,space)),space)
    end
end

# Show method:
#-------------
function Base.show{S,T,N}(io::IO,t::Tensor{S,T,N})
    print(io," Tensor ∈ $T")
    print(io,"[")
    for n=1:N
        n==1 || print(io, " ⊗ ")
        show(io,space(t,n))
    end
    println(io,"]:")
    Base.showarray(io,t.data;header=false)
end

# Basic methods for characterising a tensor:
#--------------------------------------------
space(t::Tensor,ind::Integer)=t.space[ind]
space(t::Tensor)=t.space

# General constructors
#---------------------
# with data
tensor{T<:Real,N}(data::Array{T,N})=Tensor{CartesianSpace,T,N}(data,mapreduce(CartesianSpace,⊗,size(data)))
function tensor{T<:Complex}(data::Array{T,1})
    warning("for complex array, consider specifying Euclidean index spaces")
    Tensor{ComplexEuclideanSpace,T,1}(data,⊗(ComplexSpace(size(data,1))))
end
function tensor{T<:Complex}(data::Array{T,2})
    warning("for complex array, consider specifying Euclidean index spaces")
    Tensor{ComplexEuclideanSpace,T,2}(data,ComplexSpace(size(data,1))⊗ComplexSpace(size(data,2))')
end

tensor{S,T,N}(data::Array{T},P::ProductSpace{S,N})=Tensor{S,T,N}(data,P)
tensor(data::Array,V::ElementarySpace)=tensor(data,⊗(V))

# without data
tensor{T}(::Type{T},P::ProductSpace)=tensor(Array(T,dim(P)),P)
tensor{T}(::Type{T},V::IndexSpace)=tensor(T,⊗(V))
tensor(V::Union{ProductSpace,IndexSpace})=tensor(Float64,V)

Base.similar{S,T,N}(t::Tensor{S},::Type{T},P::ProductSpace{S,N}=space(t))=tensor(similar(t.data,T,dim(P)),P)
Base.similar{S,T}(t::Tensor{S},::Type{T},V::S)=similar(t,T,⊗(V))
Base.similar{S,N}(t::Tensor{S},P::ProductSpace{S,N}=space(t))=similar(t,eltype(t),P)
Base.similar{S}(t::Tensor{S},V::S)=similar(t,eltype(t),V)

Base.zero(t::Tensor)=tensor(zero(t.data),space(t))

Base.zeros{T}(::Type{T},P::ProductSpace)=tensor(zeros(T,dim(P)),P)
Base.zeros{T}(::Type{T},V::IndexSpace)=zeros(T,⊗(V))
Base.zeros(V::Union{ProductSpace,IndexSpace})=zeros(Float64,V)

Base.rand{T}(::Type{T},P::ProductSpace)=tensor(rand(T,dim(P)),P)
Base.rand{T}(::Type{T},V::IndexSpace)=rand(T,⊗(V))
Base.rand(V::Union{ProductSpace,IndexSpace})=rand(Float64,V)

Base.eye{S<:ElementarySpace,T}(::Type{T},::Type{ProductSpace},V::S)=tensor(eye(T,dim(V)),V⊗dual(V))
Base.eye{S<:ElementarySpace}(::Type{ProductSpace},V::S)=eye(Float64,ProductSpace,V)

Base.eye{S<:ElementarySpace,T}(::Type{T},P::ProductSpace{S,2})=(P[1]==dual(P[2]) ? eye(T,ProductSpace,P[1]) : throw(SpaceError("Cannot construct eye-tensor when second space is not the dual of the first space")))
Base.eye{S<:ElementarySpace}(P::ProductSpace{S,2})=eye(Float64,P)

# tensors from concatenation
function tensorcat{S}(catind, X::Tensor{S}...)
    catind = collect(catind)
    isempty(catind) && error("catind should not be empty")
    # length(unique(catdims)) != length(catdims) && error("every dimension should appear only once")

    nargs = length(X)
    numindX = map(numind, X)

    all(n->(n == numindX[1]), numindX) || throw(SpaceError("all tensors should have the same number of indices for concatenation"))

    numindC = numindX[1]
    ncatind = setdiff(1:numindC,catind)
    spaceCvec = Array(S, numindC)
    for n = 1:numindC
        spaceCvec[n] = space(X[1],n)
    end
    for i = 2:nargs
        for n in catind
            spaceCvec[n] = directsum(spaceCvec[n], space(X[i],n))
        end
        for n in ncatind
            spaceCvec[n] == space(X[i],n) || throw(SpaceError("space mismatch for index $n"))
        end
    end
    spaceC = ⊗(spaceCvec...)
    typeC = mapreduce(eltype, promote_type, X)
    dataC = zeros(typeC, map(dim,spaceC))

    offset = zeros(Int,numindC)
    for i=1:nargs
        currentdims=ntuple(n->dim(space(X[i],n)),numindC)
        currentrange=[offset[n]+(1:currentdims[n]) for n=1:numindC]
        dataC[currentrange...] = X[i].data
        for n in catind
            offset[n]+=currentdims[n]
        end
    end
    return tensor(dataC,spaceC)
end

# Copy and fill tensors:
#------------------------
function Base.copy!(tdest::Tensor,tsource::Tensor)
    # Copies data of tensor tsource to tensor tdest if compatible
    space(tdest)==space(tsource) || throw(SpaceError("tensor spaces don't match"))
    copy!(tdest.data,tsource.data)
    return tdest
end
Base.fill!{S,T}(tdest::Tensor{S,T},value::Number)=fill!(tdest.data,convert(T,value))

# Vectorization:
#----------------
Base.vec(t::Tensor)=vec(t.data)
# Convert the non-trivial degrees of freedom in a tensor to a vector to be passed to eigensolvers etc.

# Conversion and promotion:
#---------------------------
Base.full(t::Tensor)=t.data

Base.promote_rule{S,T1,T2,N}(::Type{Tensor{S,T1,N}},::Type{Tensor{S,T2,N}})=Tensor{S,promote_type(T1,T2),N}
Base.promote_rule{S,T1,T2,N1,N2}(::Type{Tensor{S,T1,N1}},::Type{Tensor{S,T2,N2}})=Tensor{S,promote_type(T1,T2)}
Base.promote_rule{S,T1,T2}(::Type{Tensor{S,T1}},::Type{Tensor{S,T2}})=Tensor{S,promote_type(T1,T2)}

Base.promote_rule{S,T1,T2,N}(::Type{AbstractTensor{S,ProductSpace,T1,N}},::Type{Tensor{S,T2,N}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2),N}
Base.promote_rule{S,T1,T2,N1,N2}(::Type{AbstractTensor{S,ProductSpace,T1,N1}},::Type{Tensor{S,T2,N2}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2)}
Base.promote_rule{S,T1,T2}(::Type{AbstractTensor{S,ProductSpace,T1}},::Type{Tensor{S,T2}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2)}

Base.convert{S,T,N}(::Type{Tensor{S,T,N}},t::Tensor{S,T,N})=t
Base.convert{S,T1,T2,N}(::Type{Tensor{S,T1,N}},t::Tensor{S,T2,N})=copy!(similar(t,T1),t)
Base.convert{S,T}(::Type{Tensor{S,T}},t::Tensor{S,T})=t
Base.convert{S,T1,T2}(::Type{Tensor{S,T1}},t::Tensor{S,T2})=copy!(similar(t,T1),t)

Base.float{S,T<:AbstractFloat}(t::Tensor{S,T})=t
Base.float(t::Tensor)=tensor(float(t.data),space(t))

Base.real{S,T<:Real}(t::Tensor{S,T})=t
Base.real(t::Tensor)=tensor(real(t.data),space(t))
Base.complex{S,T<:Complex}(t::Tensor{S,T})=t
Base.complex(t::Tensor)=tensor(complex(t.data),space(t))

for (f,T) in ((:float32,    Float32),
              (:float64,    Float64),
              (:complex64,  Complex64),
              (:complex128, Complex128))
    @eval (Base.$f){S}(t::Tensor{S,$T}) = t
    @eval (Base.$f)(t::Tensor) = tensor(($f)(t.data),space(t))
end

# Basic algebra:
#----------------
function Base.conj!(t1::Tensor,t2::Tensor)
    space(t1)==conj(space(t2)) || throw(SpaceError())
    copy!(t1.data,t2.data)
    conj!(t1.data)
    return t1
end

# transpose inverts order of indices, compatible with graphical notation
function Base.transpose!(tdest::Tensor,tsource::Tensor)
    space(tdest)==space(tsource).' || throw(SpaceError())
    N=numind(tsource)
    TensorOperations.tensorcopy!(tsource.data,1:N,tdest.data,reverse(1:N))
    return tdest
end
function Base.ctranspose!(tdest::Tensor,tsource::Tensor)
    space(tdest)==space(tsource)' || throw(SpaceError())
    N=numind(tsource)
    TensorOperations.tensorcopy!(tsource.data,1:N,tdest.data,reverse(1:N))
    conj!(tdest.data)
    return tdest
end

Base.scale!{S,T,N}(t1::Tensor{S,T,N},t2::Tensor{S,T,N},a::Number)=(space(t1)==space(t2) ? scale!(t1.data,t2.data,a) : throw(SpaceError());return t1)
Base.scale!{S,T,N}(t1::Tensor{S,T,N},a::Number,t2::Tensor{S,T,N})=(space(t1)==space(t2) ? scale!(t1.data,a,t2.data) : throw(SpaceError());return t1)

Base.LinAlg.axpy!(a::Number,x::Tensor,y::Tensor)=(space(x)==space(y) ? Base.LinAlg.axpy!(a,x.data,y.data) : throw(SpaceError()); return y)

-(t::Tensor)=tensor(-t.data,space(t))
+(t1::Tensor,t2::Tensor)= space(t1)==space(t2) ? tensor(t1.data+t2.data,space(t1)) : throw(SpaceError())
-(t1::Tensor,t2::Tensor)= space(t1)==space(t2) ? tensor(t1.data-t2.data,space(t1)) : throw(SpaceError())

# Scalar product and norm: only valid for EuclideanSpace
Base.dot{S<:EuclideanSpace}(t1::Tensor{S},t2::Tensor{S})= (space(t1)==space(t2) ? dot(vec(t1),vec(t2)) : throw(SpaceError()))
Base.vecnorm{S<:EuclideanSpace}(t::Tensor{S})=vecnorm(t.data) # frobenius norm

# Indexing
#----------
# # linear indexing using ProductBasisVector
# Base.getindex{S,T,N}(t::Tensor{S,T,N},b::ProductBasisVector{N,S})=getindex(t.data,Base.to_index(b))
# Base.setindex!{S,T,N}(t::Tensor{S,T,N},value,b::ProductBasisVector{N,S})=setindex!(t.data,value,Base.to_index(b))

@generated function Base.getindex{G,T,N}(t::Tensor{AbelianSpace{G},T,N},s::NTuple{N,G})
    quote
        @nexprs $N n->(r_n=to_range(s[n],space(t,n)))
        return @ncall $N slice t.data r
    end
end

Base.setindex!{G,T,N}(t::Tensor{AbelianSpace{G},T,N},v::Array,s::NTuple{N,G})=(length(v)==length(t[s]) ? copy!(t[s],v) : throw(DimensionMismatch()))
Base.setindex!{G,T,N}(t::Tensor{AbelianSpace{G},T,N},v::Number,s::NTuple{N,G})=fill!(t[s],v)

# Tensor Operations
#-------------------
TensorOperations.scalar(t::Tensor)=iscnumber(space(t)) ? t.data[1] : throw(SpaceError("Not a scalar"))

function TensorOperations.add!{CA}(α, A::Tensor, ::Type{Val{CA}}, β, C::Tensor, indCinA)
    # Implements C = β*C + α*permute(op(A))
    NA = numind(A)

    spaceA = CA == :C ? conj(space(A)) : space(A)

    for i = 1:NA
        spaceA[indCinA[i]] == space(C,i) || throw(SpaceError("incompatible index spaces of tensors"))
    end

    if NA == 0 #scalars
        C.data[1] = β*C.data[1] + α*A.data[1]
    elseif indCinA == collect(1:NA) #trivial permutation
        if β == 0
            scale!(copy!(C,A), α)
        else
            Base.LinAlg.axpy!(α, A, scale!(C, β))
        end
    else
        TensorOperations.add!(α, A.data, Val{CA}, β, C.data, indCinA)
    end

    return C
end

function TensorOperations.trace!{CA}(α, A::Tensor, ::Type{Val{CA}}, β, C::Tensor, indCinA, cindA1, cindA2)
    NA = numind(A)
    NC = numind(C)

    spaceA = CA == :C ? conj(space(A)) : space(A)

    for i = 1:NC
        spaceA[indCinA[i]] == space(C,i) || throw(SpaceError("space mismatch"))
    end

    for i = 1:div(NA-NC, 2)
        spaceA[cindA1[i]] == dual(spaceA[cindA2[i]]) || throw(SpaceError("space mismatch"))
    end

    TensorOperations.trace!(α, A.data, Val{CA}, β, C.data, indCinA, cindA1, cindA2)

    return C
end

function TensorOperations.contract!{CA,CB,ME}(α, A::Tensor, ::Type{Val{CA}}, B::Tensor, ::Type{Val{CB}}, β, C::Tensor, oindA, cindA, oindB, cindB, indCinoAB, ::Type{Val{ME}}=Val{:BLAS})
    # check size compatibility
    spaceA = CA == :C ? conj(space(A)) : space(A)
    spaceB = CB == :C ? conj(space(B)) : space(B)
    spaceC = space(C)

    cspaceA = spaceA[cindA]
    cspaceB = spaceB[cindB]

    ospaceA = spaceA[oindA]
    ospaceB = spaceB[oindB]

    ospaceAB = ospaceA ⊗ ospaceB

    for i = 1:length(cspaceA)
        cspaceA[i] == dual(cspaceB[i]) || throw(SpaceError("A and B have incompatible index spaces"))
    end

    for i in 1:length(indCinoAB)
        spaceC[i] == ospaceAB[indCinoAB[i]] || throw(SpaceError("C has incompatible index space"))
    end

    if length(indCinoAB) == 0
        Cscal = pointer_to_array(pointer(C.data, 1), ())
        TensorOperations.contract!(α, A.data, Val{CA}, B.data, Val{CB}, β, Cscal, oindA, cindA, oindB, cindB, indCinoAB, Val{ME})
    else
        TensorOperations.contract!(α, A.data, Val{CA}, B.data, Val{CB}, β, C.data, oindA, cindA, oindB, cindB, indCinoAB, Val{ME})
    end

    return C
end

function TensorOperations.similar_from_indices{T,CA}(::Type{T}, indices, A::Tensor, ::Type{Val{CA}}=Val{:N})
    spaceA = CA == :C ? conj(space(A)) : space(A)
    return similar(A, T, spaceA[indices])
end

function TensorOperations.similar_from_indices{T,CA,CB}(::Type{T}, indices, A::Tensor, B::Tensor,
                                                        ::Type{Val{CA}}=Val{:N}, ::Type{Val{CB}}=Val{:N})
    spaceA = CA == :C ? conj(space(A)) : space(A)
    spaceB = CB == :C ? conj(space(B)) : space(B)
    spaceAB = spaceA ⊗ spaceB
    return similar(A, T, spaceAB[indices])
end

# Methods below are only implemented for CartesianTensor or ComplexTensor:
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
typealias ComplexTensor{T,N} Tensor{ComplexSpace,T,N}
typealias CartesianTensor{T,N} Tensor{CartesianSpace,T,N}

# Index methods
#---------------
@eval function insertind{S}(t::Tensor{S},ind::Int,V::S)
    N=numind(t)
    0<=ind<=N || throw(IndexError("index out of range"))
    iscnumber(V) || throw(SpaceError("can only insert index with c-number index space"))
    spacet=space(t)
    newspace=spacet[1:ind] ⊗ V ⊗ spacet[ind+1:N]
    return tensor(t.data,newspace)
end
@eval function deleteind(t::Tensor,ind::Int)
    N=numind(t)
    1<=ind<=N || throw(IndexError("index out of range"))
    iscnumber(space(t,ind)) || throw(SpaceError("can only delete index with c-number index space"))
    spacet=space(t)
    newspace=spacet[1:ind-1] ⊗ spacet[ind+1:N]
    return tensor(t.data,newspace)
end

for (S,TT) in ((CartesianSpace,CartesianTensor),(ComplexSpace,ComplexTensor))
    @eval function fuseind(t::$TT,ind1::Int,ind2::Int,V::$S)
        N=numind(t)
        ind2==ind1+1 || throw(IndexError("only neighbouring indices can be fused"))
        1<=ind1<=N-1 || throw(IndexError("index out of range"))
        fuse(space(t,ind1),space(t,ind2),V) || throw(SpaceError("index spaces $(space(t,ind1)) and $(space(t,ind2)) cannot be fused to $V"))
        spacet=space(t)
        newspace=spacet[1:ind1-1]*V*spacet[ind2+1:N]
        return tensor(t.data,newspace)
    end
    @eval function splitind(t::$TT,ind::Int,V1::$S,V2::$S)
        1<=ind<=numind(t) || throw(IndexError("index out of range"))
        fuse(V1,V2,space(t,ind)) || throw(SpaceError("index space $(space(t,ind)) cannot be split into $V1 and $V2"))
        spacet=space(t)
        newspace=spacet[1:ind-1]*V1*V2*spacet[ind+1:N]
        return tensor(t.data,newspace)
    end
end

# Factorizations:
#-----------------
for (S,TT) in ((CartesianSpace,CartesianTensor),(ComplexSpace,ComplexTensor))
    @eval function svd!(t::$TT,n::Int,::NoTruncation)
        # Perform singular value decomposition corresponding to bipartion of the
        # tensor indices into the left indices 1:n and remaining right indices,
        # thereby destroying the original tensor.
        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)
        newdim=length(F[:S])
        newspace=$S(newdim)
        U=tensor(F[:U],leftspace*newspace')
        S=tensor(diagm(F[:S]),newspace*newspace')
        V=tensor(F[:Vt],newspace*rightspace)
        return U,S,V,abs(zero(eltype(t)))
    end

    @eval function svd!(t::$TT,n::Int,trunc::Union{TruncationDimension,TruncationSpace})
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        dim(trunc) >= min(leftdim,rightdim) && return svd!(t,n)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)
        sing=F[:S]

        # compute truncation dimension
        if eps(trunc)!=0
            sing=F[:S]
            normsing=vecnorm(sing)
            truncdim=0
            while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
                truncdim+=1
            end
            truncdim=min(truncdim,dim(trunc))
        else
            truncdim=dim(trunc)
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
        U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
        S=tensor(diagm(sing[1:truncdim]),newspace*newspace')
        V=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
        return U,S,V,truncerr
    end

    @eval function svd!(t::$TT,n::Int,trunc::TruncationError)
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)

        # find truncdim based on trunctolinfo
        sing=F[:S]
        normsing=vecnorm(sing)
        truncdim=0
        while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
            truncdim+=1
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/normsing
        U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
        S=tensor(diagm(sing[1:truncdim]),newspace*newspace')
        V=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
        return U,S,V,truncerr
    end

    @eval function leftorth!(t::$TT,n::Int,::NoTruncation)
        # Create orthogonal basis U for indices 1:n, and remainder C for right
        # indices, thereby destroying the original tensor.
        # Decomposition should be unique, such that it always returns the same
        # result for the same input tensor t. UC = QR is fastest but only unique
        # after correcting for phases.
        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        if leftdim>rightdim
            newdim=rightdim
            tau=zeros(eltype(data),(newdim,))
            Base.LinAlg.LAPACK.geqrf!(data,tau)

            C=zeros(eltype(t),(newdim,newdim))
            for j in 1:newdim
                for i in 1:j
                    @inbounds C[i,j]=data[i,j]
                end
            end
            Base.LinAlg.LAPACK.orgqr!(data,tau)
            U=data
            for i=1:newdim
                tau[i]=sign(C[i,i])
            end
            scale!(U,tau)
            scale!(conj!(tau),C)
        else
            newdim=leftdim
            C=data
            U=eye(eltype(data),newdim)
        end

        newspace=$S(newdim)
        return tensor(U,leftspace*newspace'), tensor(C,newspace*rightspace), abs(zero(eltype(t)))
    end

    @eval function leftorth!(t::$TT,n::Int,trunc::Union{TruncationDimension,TruncationSpace})
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        dim(trunc) >= min(leftdim,rightdim) && return leftorth!(t,n)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)
        sing=F[:S]

        # compute truncation dimension
        if eps(trunc)!=0
            sing=F[:S]
            normsing=vecnorm(sing)
            truncdim=0
            while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
                truncdim+=1
            end
            truncdim=min(truncdim,dim(trunc))
        else
            truncdim=dim(trunc)
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
        U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
        C=tensor(scale!(sing[1:truncdim],F[:Vt][1:truncdim,:]),newspace*rightspace)
        return U,C,truncerr
    end

    @eval function leftorth!(t::$TT,n::Int,trunc::TruncationError)
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)

        # compute truncation dimension
        sing=F[:S]
        normsing=vecnorm(sing)
        truncdim=0
        while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
            truncdim+=1
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
        U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
        C=tensor(scale!(sing[1:truncdim],F[:Vt][1:truncdim,:]),newspace*rightspace)

        return U,C,truncerr
    end

    @eval function rightorth!(t::$TT,n::Int,::NoTruncation)
        # Create orthogonal basis U for right indices, and remainder C for left
        # indices. Decomposition should be unique, such that it always returns the
        # same result for the same input tensor t. CU = LQ is typically fastest but only
        # unique after correcting for phases.
        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        if leftdim<rightdim
            newdim=leftdim
            tau=zeros(eltype(data),(newdim,))
            datat=transpose(data)
            Base.LinAlg.LAPACK.geqrf!(datat,tau)

            C=zeros(eltype(t),(newdim,newdim))
            for j in 1:newdim
                for i in 1:j
                    @inbounds C[j,i]=datat[i,j]
                end
            end
            Base.LinAlg.LAPACK.orgqr!(datat,tau)
            Base.transpose!(data,datat)
            U=data

            for i=1:newdim
                tau[i]=sign(C[i,i])
            end
            scale!(tau,U)
            scale!(C,conj!(tau))
        else
            newdim=rightdim
            C=data
            U=eye(eltype(data),newdim)
        end

        newspace=$S(newdim)
        return tensor(C,leftspace ⊗ dual(newspace)), tensor(U,newspace ⊗ rightspace), abs(zero(eltype(t)))
    end

    @eval function rightorth!(t::$TT,n::Int,trunc::Union{TruncationDimension,TruncationSpace})
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        dim(trunc) >= min(leftdim,rightdim) && return rightorth!(t,n)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)
        sing=F[:S]

        # compute truncation dimension
        if eps(trunc)!=0
            sing=F[:S]
            normsing=vecnorm(sing)
            truncdim=0
            while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
                truncdim+=1
            end
            truncdim=min(truncdim,dim(trunc))
        else
            truncdim=dim(trunc)
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
        C=tensor(scale!(F[:U][:,1:truncdim],sing[1:truncdim]),leftspace*newspace')
        U=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
        return C,U,truncerr
    end

    @eval function rightorth!(t::$TT,n::Int,trunc::TruncationError)
        # Truncate rank corresponding to bipartition into left indices 1:n
        # and remain right indices, based on singular value decomposition,
        # thereby destroying the original tensor.
        # Truncation parameters are given as keyword arguments: trunctol should
        # always be one of the possible arguments for specifying truncation, but
        # truncdim can be replaced with different parameters for other types of tensors.

        N=numind(t)
        spacet=space(t)
        leftspace=spacet[1:n]
        rightspace=spacet[n+1:N]
        leftdim=dim(leftspace)
        rightdim=dim(rightspace)
        data=reshape(t.data,(leftdim,rightdim))
        F=svdfact!(data)

        # compute truncation dimension
        sing=F[:S]
        normsing=vecnorm(sing)
        truncdim=0
        while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
            truncdim+=1
        end
        newspace=$S(truncdim)

        # truncate
        truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
        C=tensor(scale!(F[:U][:,1:truncdim],sing[1:truncdim]),leftspace*newspace')
        U=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
        return C,U,truncerr
    end
end

# Methods below are only implemented for CartesianMatrix or ComplexMatrix:
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
typealias ComplexMatrix{T} ComplexTensor{T,2}
typealias CartesianMatrix{T} CartesianTensor{T,2}

function Base.pinv(t::Union{ComplexMatrix,CartesianMatrix})
    # Compute pseudo-inverse
    spacet=space(t)
    data=copy(t.data)
    leftdim=dim(spacet[1])
    rightdim=dim(spacet[2])

    F=svdfact!(data)
    Sinv=F[:S]
    for k=1:length(Sinv)
        if Sinv[k]>eps(Sinv[1])*max(leftdim,rightdim)
            Sinv[k]=one(Sinv[k])/Sinv[k]
        end
    end
    data=F[:V]*scale(F[:S],F[:U]')
    return tensor(data,spacet')
end

function Base.eig(t::Union{ComplexMatrix,CartesianMatrix})
    # Compute eigenvalue decomposition.
    spacet=space(t)
    spacet[1] == spacet[2]' || throw(SpaceError("eigenvalue factorization only exists if left and right index space are dual"))
    data=copy(t.data)

    F=eigfact!(data)

    Lambda=tensor(diagm(F[:values]),spacet)
    V=tensor(F[:vectors],spacet)
    return Lambda, V
end

function Base.inv(t::Union{ComplexMatrix,CartesianMatrix})
    # Compute inverse.
    spacet=space(t)
    spacet[1] == spacet[2]' || throw(SpaceError("inverse only exists if left and right index space are dual"))

    return tensor(inv(t.data),spacet)
end


text:
module ImageProcessing

using Images
using Base.Cartesian

include("blockproc.jl")
include("filtering.jl")

export
  im2col,
  im2col!,
  col2im,
  col2im!,
  imfiltermtx,
  psf2otf,
  padforfilter

end