text:
import FileIO.load
using AndorSIF
using CSV
using Dates
using DelimitedFiles
using JSON: json
import LightXML
using HDF5

"""
    list_spectra(; exact="", inexact="", date::Tuple{Int64,Int64,Int64}, group=false)

Return a `DataFrame` with matching spectra.
First `grab(datafolder)` the spectra to genereate a `.spectralist` file inside your current working
directory. Note that you have to have the `.spectralist` file in your current working directory
for `list_spectra` to work.

Filter the dataframe:
* `exact="name"`: Match the exact name of the spectrum
* `inexact="part"`: Check if the name contains part of the string
* `date=(2011, 01, 31)`: Filter for specific date
* `group=true`: Groups spectra by name

The result with some useful information is stored in a DataFrame. No spectra are
loaded hereby to save time. Load the actual spectra via `load_spectra(id)`
where `id` is `list_spectra().id` See the DataFrames package for more info.
"""
function list_spectra(; exact=""::AbstractString,
                        inexact=""::AbstractString,
                        date=(0,0,0)::Tuple{Int64,Int64,Int64},
                        group=false)

  dir = "./"

  spectrafile = joinpath(dir, ".spectralist")
  if isfile(spectrafile)
    dat = readdlm(spectrafile; comments=false)
  else
    error("The file $spectrafile does not exist.")
  end

  df = CSV.File(spectrafile) |> DataFrame

  # Filter the dataframe
  if exact != ""
      df = df[lowercase.(df.name) .== lowercase(exact), :]
  end

  if inexact != ""
      df = df[occursin.(lowercase(inexact), lowercase.(df.name)), :]
  end

  if !all(iszero.(date))
      df = df[Dates.yearmonthday.(df.date) .== [date], :]
  end

  sort!(df)

  if group
      df = combine(df -> DataFrame(
            N = length(df.id),
            sizes = [unique(df.sizes)],
            dates = [unique(floor.(df.date, Dates.Day(1)))],
            id = [df.id]
        ),
        groupby(df, :name)
      )
  end

  return df

end

"""
    grab(dir="./"; getall=false, singledata="",delete_dirs=false)

Make a data file that contains information where to find spectra connected to
an ID. All spectra are internally referenced by this id. `dir` is the main data
directory (`dir/2011-01-31/Name/...`).

By default this function writes only newly found spectra to the .spectralist file.
If you want to rewrite the whole file pass `getall=true` as a keyword argument
to the function.

By default grab checks if there are as many .sif files as .xml files in a directory,
if there arent this directory is not grabbed. Use the flag delete_dirs=true if you want
to delete these directorys. Directorys containing the keyword "DARK" wont be removed.

Returns the number of added spectra and the number of spectra in total.

Optional: specify a single data folder to grab with singledata="path/to/datafolder"
"""
function grab(dir="./"; getall=false, singledata= "",delete_dirs=false)

    missing_sifs = []

    if delete_dirs ==false

        push!(missing_sifs,".sif Files are missing in the following directories:")

    elseif delete_dirs == true

        push!(missing_sifs,""".sif Files are missing in the following directories and therefore deleted, except "DARK" dirs.""")
    end

    #grab single data folder in $singledata?
    if singledata !== ""

        idexisting = Int64[]
        idlist = Int64[]
        namelist = String[]
        dirlist = String[]
        datelist = DateTime[]
        sizelist = []
        numlist = Int64[]

        for (root, dirs, files) in walkdir(singledata)
            for dir in dirs
                if dir == "raw"
                    # returns a boolean if there are as many .sif files as .xml files in dir
                    if check_sif_files(joinpath(root, dir)) == true

                        # search for xml files first
                        format = :xml
                        mlist = searchdir(joinpath(root, dir), ".xml")

                        # if still empty throw an error
                        isempty(mlist) && error("No meta files found in $(root)/$(dir)")

                        mdict = get_metadata(joinpath(root, dir, mlist[1]))

                        # Get ID
                        id = Int64(Dates.value(DateTime(mdict["timestamp"])))
                        if !any(idexisting .== id)
                            push!(idlist, id)
                            push!(namelist, splitdir(splitdir(root)[1])[2])
                            push!(dirlist, joinpath(abspath(root), dir))
                            push!(datelist, DateTime(mdict["timestamp"]))
                            if format == :xml
                                sif_files = searchdir(joinpath(root, dir), ".sif")
                                filestats = stat(joinpath(root, dir, sif_files[1]))
                                specsize = filestats.size
                                push!(sizelist, specsize)
                            elseif format == :txt
                                push!(sizelist, Int64[N_PIXEL/mdict["x_binning"], N_PIXEL/mdict["y_binning"], length(mlist)])
                            end

                            push!(numlist, parse(Int64, splitdir(splitdir(dirlist[end])[1])[2]))
                        end
                    elseif check_sif_files(joinpath(root, dir)) == false
                        push!(missing_sifs,root)
                    end
                end
            end
        end

        # # create new SPECTRALIST file
        # df = DataFrame(id=idlist, name=namelist, path=dirlist, date=datelist, sizes=sizelist, number=numlist)
        # CSV.write(".spectralist", df; append=false)

    # grab all spectra files
    else

        # Check if the .spectralist file exists.
        # If it exist read its contents.
        if !isfile(".spectralist") || getall
            idexisting = Int64[]
        else
            df = CSV.read(".spectralist", DataFrame)
            idexisting = convert.(Int64, df.id)
        end

        idlist = Int64[]
        namelist = String[]
        dirlist = String[]
        datelist = DateTime[]
        sizelist = []
        numlist = Int64[]

        for (root, dirs, files) in walkdir(dir)
            for dir in dirs
                if dir == "raw"
                    if check_sif_files(joinpath(root, dir)) == true
                        # search for xml files first
                        format = :xml
                        mlist = searchdir(joinpath(root, dir), ".xml")
                        # if the list is empty search for txt files
                        if isempty(mlist)
                            mlist = searchdir(joinpath(root, dir), ".txt")
                            format = :txt
                        end

                        # if still empty throw an error
                        isempty(mlist) && error("No meta files found in $(root)/$(dir)")

                        mdict = get_metadata(joinpath(root, dir, mlist[1]))

                        # Get ID
                        id = Int64(Dates.value(DateTime(mdict["timestamp"])))
                        if !any(idexisting .== id)
                            push!(idlist, id)
                            push!(namelist, splitdir(splitdir(root)[1])[2])
                            push!(dirlist, joinpath(abspath(root), dir))
                            push!(datelist, DateTime(mdict["timestamp"]))
                            if format == :xml
                                sif_files = searchdir(joinpath(root, dir), ".sif")
                                filestats = stat(joinpath(root, dir, sif_files[1]))
                                specsize = filestats.size
                                push!(sizelist, specsize)
                            elseif format == :txt
                                push!(sizelist, Int64[N_PIXEL/mdict["x_binning"], N_PIXEL/mdict["y_binning"], length(mlist)])
                            end
                            push!(numlist, parse(Int64, splitdir(splitdir(dirlist[end])[1])[2]))
                        end
                    elseif check_sif_files(joinpath(root, dir)) == false
                        push!(missing_sifs,root)
                    end
                end
            end
        end
    end

    df = DataFrame(id=idlist, name=namelist, path=dirlist, date=datelist, sizes=sizelist, number=numlist)

    # open(".spectralist", "a") do f
    #     # writedlm(f, [idlist namelist dirlist datelist sizelist numlist])
    #     CSV.write(f, df)
    # end

    if !isfile(".spectralist") || getall || singledata != ""
        CSV.write(".spectralist", df; append=false)
    else
        CSV.write(".spectralist", df; append=true)
    end

    # println("Collected $(length(idlist)) spectra. ($(length(idexisting) + length(idlist)) overall)")
    if length(missing_sifs) > 1
        println.(missing_sifs);
        if delete_dirs == true
            function isdark!(dirlist)
                filter!(x->occursin("DARK",x) == false,dirlist)
            end
            isdark!(missing_sifs)
            rm.(missing_sifs,recursive=true)
        end
    end

    length(idlist), length(idexisting) + length(idlist)

end

"""
    load_spectra(id::Int64)

Load the spectra specified by `id`.
"""
function load_spectra(id::Int64; format=:sif)
  # Load grabbed data
  dir = getdir(id)

  # Load the spectrum files
  s = read_as_3D(dir; format=format)

  sfspectrum = SFSpectrum[SFSpectrum(id, s)]
end

"""
    load_spectra(id::Array{Int64})

Load the spectra specified by `id`.
"""
function load_spectra(id::AbstractArray{Int64}; format=:sif)
  sfspectra = Array{SFSpectrum,1}(undef, length(id))
  for i in 1:length(id)
    sfspectrum = load_spectra(id[i], format=format)
    sfspectra[i] = sfspectrum[1]
  end
  return sfspectra
end

"""
Get an attribute of a spectrum or an array of spectra. To get available
attribute list these via `get_metadata(spectrum)`.
"""
function get_attribute(s::Array{SFSpectrum}, attr::AbstractString)
    values = []
    for i in s
        val = get_attribute(i, attr)
        push!(values, val)
    end
    return values
end

function get_attribute(id::Int64, attr::AbstractString)

    """
    Flattens an array of dicts into a dict of arrays.
    """
    function flattendictarray(dictarray)
    typeof(dictarray) <: Dict && (return dictarray)
    d = Dict()
        for k in keys(dictarray[1])
            d[k] = [x[k] for x in dictarray]
        end
    d
    end

    dict = get_metadata(id)
    # check if we deal with an old txt metadatafile or a new xml metadata file
    # (the old one should have an x_binning key somewhere)
    haskey(dict, "x_binning") ? (format = :txt) : (format = :xml)
    if format == :txt
        # metafile is in txt format. just return the value for the key
        attr in keys(dict) ? val = dict[attr] : val = nothing
        return val
    elseif format == :xml
        # format is xml lookup the dictionary to get to the right entry
        translate = Dict(
            "name" => splitpath(getdir(id))[end-2],
            "comment" => dict["comment"],
            "grating" => "n/a", #dict["spectra pro 300"]["grating"],
            "spectrometer_wavelength" => dict["spectra pro 300"]["wavelength set"],
            "mirror_position" => "n/a", #dict["spectra pro 300"]["mirror position"],
            "pump_dl_position" => dict["ls-180"]["pos"],
            "ccd_exposure_time" => dict["ixon"]["exposure_time"],
            "ccd_readout_speed" => "n/a",
            "x_binning" => dict["ixon"]["horizontal_bin"],
            "y_binning" => dict["ixon"]["vertical_bin"],
            "ccd_temperature" => dict["ixon"]["temperature"],
            "timestamp" => dict["timestamp"],
            "twin1_wavelength" => dict["twins"]["twin1 wavelength"],
            "twin2_wavelength" => dict["twins"]["twin2 wavelength"],
            "twin1_shutter" => dict["twins"]["twin1 shutter"],
            "twin2_shutter" => dict["twins"]["twin2 shutter"],
            "pump_ir_wavelength" => dict["ekspla laser"]["ekspla wavelength"],
            "micos_position_x" => dict["smc stages"]["xpos"],
            "micos_position_y" => dict["smc stages"]["ypos"],
            "micos_position_z" => dict["smc stages"]["zpos"],
            "probe_ir_power" => dict["ir power meters"]["probe ir power"] |> flattendictarray,
            "pump_ir_power" => dict["ir power meters"]["pump ir power"] |> flattendictarray,
            "vis_power" => dict["pm100usb"]["sample power"],
        )
        return translate[attr]
    end
end

get_attribute(s::SFSpectrum, attr::AbstractString) = get_attribute(s.id, attr)

"""
Check if a spectrum has the attribute `attr`.
"""
function is_attribute(id::Int64, attr::AbstractString)
    dict = get_metadata(id)
    attr in keys(dict)
end

function is_attribute(s::Array{SFSpectrum}, attr::AbstractString)
    values = []
    for i in s
        val = is_attribute(i, attr)
        push!(values, val)
    end
    return values
end

is_attribute(s::SFSpectrum, attr::AbstractString) = is_attribute(s.id, attr)

"""
Read tiff files and put them in a 3D matrix
"""
function read_as_3D(path::AbstractString; format=:sif)
    filelist = String[]
    format == :sif && (filelist = searchdir(path, ".sif"))
    # if there are no sif files check for tiff files instead
    isempty(filelist) && (format = :tiff)
    format == :tiff && (filelist = searchdir(path, ".tiff"))
    if isempty(filelist)
        error("Could not find a file with format sif or tiff in $directory.")
    end

    if format == :sif
        # transpose the image that's loaded from the sif file
        I = load(joinpath(path, filelist[1]))
        C = Array{Float64,3}(undef, size(I,1), size(I,2), length(filelist))
        C[:,:,1] .= I[:,:,1]
    elseif format == :tiff
        I = load(joinpath(path, filelist[1]))
        C = Array{eltype(I),3}(undef, size(I,1), size(I,2), length(filelist))
        C[:,:,1] .= I
    end
    @inbounds for i = 2:length(filelist)
        if format == :sif
            C[:,:,i] .= load(joinpath(path, filelist[i]))[:,:,1]
        elseif format == :tiff
            C[:,:,i] .= load(joinpath(path, filelist[i]))
        end
    end
    if format == :sif
        F = C
    elseif format == :tiff
        F = (C .|> Float64) .* typemax(UInt16) .|> round
    end
    return F
end

"""
Get the directory of a spectrum with a given id. If the ID does not exist return an empty string.
"""
function getdir(id::Int64)
    df = CSV.read(".spectralist", DataFrame)
    idx = findall((in)(id), df.id)
    isempty(idx) && error("Could not find spectrum with id $id.")
    dir = df.path[idx[1]]
end

"""
Reads an xml metadata file
"""
function read_xml(path::String)

    function loopelements(element, dict)
        # println("Looping")
        for e in LightXML.child_elements(element)
            # println("Current Element Type: $(LightXML.name(e))")
            if haskey(lvtypes, LightXML.name(e))
                # Go deeper
                if any(LightXML.name(e) .== ["Cluster", "Array"])
                    # Module Data
                    clustername = LightXML.find_element(e, "Name") |> LightXML.content
                    # Get rid of the Module Data in the Name of the Cluster,
                    # strip whitespace, make lowercase
                    clustername_clean =
                        replace(clustername, "Module Data" => "") |> strip |> lowercase
                    # println("Found Cluster: $clustername")
                    dict[clustername_clean] = Dict()
                    loopelements(e, dict[clustername_clean])
                else
                    n, v = parseelement(e)
                    # println("Found Element $n with value $v")
                    global dict[n] = v
                end
            elseif any(LightXML.name(e) .== ["NumElts", "Name", "Dimsize"])
                # Ignore entry
                # println("Ignoring")
                continue
            else
                @warn "Could not find type for $(LightXML.name(e))"
            end
        end
    end

    function parseelement(element)
        typ = lvtypes[LightXML.name(element)]
        n = LightXML.find_element(element, "Name") |> LightXML.content |> strip |> lowercase
        v_str = LightXML.find_element(element, "Val") |> LightXML.content
        if typ == String
            global v = v_str
        elseif typ == Bool
            if v_str == "1"
                v = true
            elseif v_str == "0"
                v = false
            end
        else
            v = parse(typ, v_str)
        end
        n, v
    end

    lvtypes = Dict(
        "Cluster" => Dict,
        "Array"   => Array,
        "Refnum"  => String,
        "DBL"     => Float32,
        "I32"     => Int32,
        "I64"     => Int64,
        "U8"      => UInt8,
        "String"  => String,
        "Boolean" => Bool,
    )

    dict = Dict()

    xdoc = LightXML.parse_file(path)
    xroot = LightXML.root(xdoc)

    loopelements(xroot, dict)

    LightXML.free(xdoc)

    dict
end

"""
Get all available metadata for a spectrum. Takes a SFSpectrum struct or an id as
an Int64.
"""
get_metadata(s::SFSpectrum) = get_metadata(s.id)

function get_metadata(id::Int64)
    dir = getdir(id)
    mdict = get_metadata(dir)
end

function get_metadata(path::AbstractString)
    path
    path == "" && (return Dict())

    mdict = Dict{String, Any}()
    mfiles = Array{String}[]

    if splitext(path)[2] != ".txt" && splitext(path)[2] != ".xml"
        # The path is a folder
        @assert isdir(path)
        # Search for xml files by default
        mfiles = searchdir(path, ".xml")
        format = :xml
        # If no xml files were found search for txt files
        if length(mfiles) == 0
            mfiles = searchdir(path, ".txt")
            format = :txt
        end
        paths = Array{String}(undef, size(mfiles, 1))
        [paths[i] = joinpath(path, mfiles[i]) for i = 1:length(mfiles)]
    else
        # Single file
        @assert isfile(path)
        paths = [path]
        mfiles = [path]
        # Get the file extenstion of the input file
        format_str = splitext(path)[2]
        format = replace(format_str, "." => "") |> lowercase |> Symbol
    end

    if format == :xml
        # get the first metadatafile to extract the keys and initialize arrays
        # with the values of the right type
        dict_first = read_xml(paths[1])
        dictkeys = keys(dict_first) #the values of these keys can be dictionaries
        if length(paths) > 1
            for dictkey in dictkeys, k in keys(dict_first[dictkey])
                # So we loop through these dictionaries if dictkey is indeed a
                # dictonary. Otherwise we just put it in on the first layer of
                # the dictionary
                if typeof(dict_first[dictkey]) <: Dict
                    val = dict_first[dictkey][k]
                    valtype = typeof(val)
                    !haskey(mdict, dictkey) &&
                        (mdict[dictkey] = Dict{String, Array{Any}}())
                    # Put the values
                    mdict[dictkey][k] = valtype[val]
                else
                    val = dict_first[dictkey]
                    valtype = typeof(val)
                    # Put the values
                    mdict[dictkey] = valtype[val]
                    # continue # otherwise we loop through all the keys
                end
            end
            for i = 2:length(paths), dictkey in dictkeys, k in keys(dict_first[dictkey])
                dict = read_xml(paths[i])
                if typeof(dict[dictkey]) <: Dict
                    dictkey, k
                    val = dict[dictkey][k]
                    push!(mdict[dictkey][k], val)
                else
                    val = dict[dictkey]
                    push!(mdict[dictkey], val)
                    # continue # otherwise we loop through all the keys
                end
            end
        elseif length(paths) == 1
            mdict = dict_first
        end

        # next we read the metadata from the .sif file
        # First search for the sif files in the folder
        if length(paths) == 1
            # We deal with a single file which has a .xml extension that we need
            # to get rid of and change it to .sif.
            #file  = searchdir(paths, "sif")[1]

            file_dir,filename = splitdir(paths[1])

            sif_files = [joinpath(file_dir, searchdir(file_dir, ".sif")[1])]

        else
            # We deal with a path and perhaps multiple files
            file_dir,filename = splitdir(paths[1])
            sif_files = joinpath.(file_dir, searchdir(file_dir, ".sif"))

        end
        if length(sif_files) == 1
            # load the first metadata
            ixon_meta = load(sif_files[1]).ixon
            mdict["ixon"] = ixon_meta
        elseif length(sif_files) > 1
            # load the first metadata
            ixon_meta = load(sif_files[1]).ixon
            mdict["ixon"] = Dict{String, Any}()
            for key in keys(ixon_meta)
                val = ixon_meta[key]
                typ = typeof(val)
                # make an array of the proper type for value val
                mdict["ixon"][key] = typ[val]
            end
            # and loop through the remaining spectra in the folder
            for i = 2:length(sif_files), key in keys(ixon_meta)
                dict = load(joinpath(path, sif_files[i])).ixon
                push!(mdict["ixon"][key], dict[key])
            end
        else
            @warn "No corresponding sif file found in $path"
        end

    elseif format == :txt
        values = Array{Any,2}
        data = readdlm(paths[1], '\t'; comments=false)
        keyz = Array{String}
        keyz = data[:,1]
        value = data[:,2]
        if length(mfiles) == 1
          for (i, key) in enumerate(keyz)
              mdict[key] = value[i]
          end
        elseif length(mfiles) > 1
          valuemat = Array{Any,2}(undef, length(keyz), length(mfiles))
          valuemat[:,1] = data[:,2]
          for i = 2:length(mfiles)
              data = readdlm(paths[i], '\t'; comments=false)
              valuemat[:,i] = data[:,2]
          end
          for (i, key) in enumerate(keyz)
              mdict[key] = valuemat[i,:]
          end
        end
    end

    # We need the comment only once
    # length(mfiles) > 1 && (mdict["comment"] = mdict["comment"][1])
    !haskey(mdict, "comment") && (mdict["comment"] = "")

    return mdict
  end

function searchdir(directory::AbstractString, key::AbstractString)
    filter!(x->occursin(key, x), readdir(directory))
end

function failed_sifs(dir = "./")

    failed_sifs = []

    for (root,dirs,files) in walkdir(dir)
        for dir in dirs
            if dir == "raw"
                filenames= searchdir(joinpath(root, dir),".sif")
                    for file in filenames
                        ind_counter = last(findfirst("data",file)) + 1

                        if typeof(try parse(Int,file[ind_counter]) catch end) == Int64 && try parse(Int,file[ind_counter]) catch end > 0

                            push!(failed_sifs,joinpath(root,dir,file))
                        end
                    end
            end
         end
    end

    if length(failed_sifs) > 0

        return failed_sifs

    elseif length(failed_sifs) == 0

       return "Hurray!! Labview managed to save all .sif files at her first attempt."

    end
end

"""
`savejson(path::String, spectra::Array{SFSpectrum{T,1},1})`
Saves the array `spectra` to `path` in JSON format.

This is primarily to read the data set by the Matlab fit program.
"""
function savejson(path::String, spectra::Array{SFSpectrum{T,1},1}) where T
    splitext(path)[end] != ".json" && (path = path * ".json")

    for i = 1:length(spectra)
        length(size(spectra[i])) > 1 && error("Please make the Spectra 1D")
    end

    ids = Int64[]
    wavelengths = Array{Float64}[]
    wavenumbers = Array{Float64}[]
    signals = Array{Float64}[]
    names = String[]

    for i = 1:length(spectra)
        λ = get_ir_wavelength(spectra[i])
        ν = 1 ./ λ .* 1e7

        push!(ids, spectra[i].id)
        push!(wavelengths, λ)
        push!(wavenumbers, ν)
        push!(signals, spectra[i].s)
        push!(names, get_attribute(spectra[i], "name"))
    end

    dict = Dict(
        "ids"         => ids,
        "wavelengths" => wavelengths,
        "wavenumbers" => wavenumbers,
        "signals"     => signals,
        "names"       => names
    )

    jsonstring = json(dict)

    open(path, "w+") do io
        write(io, jsonstring)
    end

end

"""
Check a dir if .sif files are missing. If delete_dirs = true directories with missing .sif Files are removed
    check_sif_files(dirlist; delete_dirs = false)
"""

function check_sif_files(dir)
    xml_files = searchdir(dir,".xml")
    sif_files = searchdir(dir,".sif")

    if length(xml_files) == length(sif_files)
        return true
    else
        return false
    end
end

    #function eq2one(array)
    #        all(x->x==1,array)
    #end
      # compare if each xml file also has a matching sif file
        #compare_xml_sif =    occursin.(_xml_files,_sif_files)

"""
Convert a date (yyyy,mm,dd) of type Tuple{Int64,Int64,Int64} to a yyyymmdd string.
"""
function date2dashboard(date::Tuple{Int64,Int64,Int64})
    if ndigits(date[1]) == 4
        yyyy = string(date[1])
    else
        error("$string(date[1]) is not a year. Format must be yyyy.")
    end

    if ndigits(date[2]) == 2
        mm = string(date[2])
    elseif  ndigits(date[2]) == 1
        mm = "0"*string(date[2])
    else
        error("$string(date[2]) is not a month. Format must be either mm or m.")
    end

    if ndigits(date[3]) == 2
        dd = string(date[3])
    elseif  ndigits(date[3]) == 1
        dd = "0"*string(date[3])
    else
        error("$string(date[3]) is not a day. Format must be either dd or d.")
    end

    return yyyy*mm*dd
end

"""
Change the - or + to ⁻ or ⁺
"""

function convert2unicode(scantype)
    if occursin("-",scantype)
        return replace(scantype,"-"=> "⁻")
    elseif occursin("+",scantype)
        return replace(scantype,"+" => "⁺")
    else
        return scantype
    end
end

""" Save the data in a HDF5 File. \n

Register a new Sample name or choose a sample name from the SFGDashboard Sample Dropdown, e.g. "CaAra" or "CaAra d4". If you are not sure visit --> 132.252.80.59:8050/\n
Example:\n
    \tsample = "CaAra" \n
    \tsample_prep = "20200504"\n
    \tfoldername =  "CaAra_20200504_1_DL-Scan_r--pumped_ppp"\n
    save_data(sample,32,"ppp","DL d-")\n
    save_data(  sample::String, surface_density_value::Int, polarisation_comb::String, scan::String;\n
                date=date, sample_prep= sample_prep, foldername= foldername, raw_spectra = raw,\n
                sigmatrix = sig03, refmatrix = ref03, probe_wavenumbers = ν, delay_time = dltime_sorted,\n
                pump_wavenumbers = nothing, mode_name = mode_name, \n
                sig_bleaches = [mean(sig03[:,pixel[i]], dims=2)[:,1] for i in 1:length(mode_name)],\n
                ref_bleaches = [mean(ref03[:,pixel[i]], dims=2)[:,1] for i in 1:length(mode_name)],\n
                add_comment= ""
    )
"""
function save_data( sample::String, surface_density_value::Int, polarisation_comb::String, scan::String;
                    date=Main.date, sample_prep= Main.sample_prep, foldername= Main.foldername, raw_spectra = Main.raw,
                    sigmatrix = Main.sig03, refmatrix = Main.ref03, probe_wavenumbers = Main.ν, delay_time = Main.dltime_sorted,
                    pump_wavenumbers = nothing, mode_name = Main.mode_name,
                    sig_bleaches = [mean(Main.sig03[:,Main.pixel[i]], dims=2)[:,1] for i in 1:length(Main.mode_name)],
                    ref_bleaches = [mean(Main.ref03[:,Main.pixel[i]], dims=2)[:,1] for i in 1:length(Main.mode_name)],
                    add_comment= ""
        )

    # Check if date has the right type
    if typeof(date) == String
        directory = date * "/" * foldername
        dashboard_date = date
    elseif  typeof(date) !== String
        dashboard_date = date2dashboard(date)
        directory = dashboard_date * "/" * foldername
    elseif length(date) !== 8
        error("$date has not the type  yyyymmdd")
    end

    # sample surface density
    surface_density = "$surface_density_value mNm⁻¹"

    # scan type (delay or wavenumber scan)
    if occursin("dl",lowercase(scan)) == true
        scan_type = "delay_scan"
        pump_resonance= convert2unicode(split(scan)[2])*"pumped"
    elseif occursin("wl",lowercase(scan)) == true || occursin("wn",lowercase(scan)) == true
        scan_type = "wavenumber_scan"
    else
        error("""Scan type cannot be $scan !!\n
        Example for delay scan with d- pumped scan ="dl d-"\n
        Example for wavelength scan scan = "wl".
        """)
    end

    # generate filename
    filename = sample *"_"* "$surface_density_value"*"mNm-1"*"_"* polarisation_comb *"_"*split(scan)[1]*"_"*split(scan)[2] *"_"* sample_prep * dashboard_date*".h5"

    # calculate pump wavenumber (Ekspla) for delay scan
    ekspla_wavelength = get_metadata(raw_spectra[1])["ekspla laser"]["ekspla wavelength"][1]
    ekspla_wavenumber = round(10^7 / ekspla_wavelength, digits=2)

    h5open(filename, "w") do fid
        g0 = g_create(fid, sample)
        g1 = g_create(g0, surface_density)
        g2 = g_create(g1, polarisation_comb)
        g3 = g_create(g2, scan_type)

        if scan_type == "delay_scan"

            if first(size(sigmatrix)) !== first(size(delay_time))
                error("Dimensions of sigmatrix and delay_time must match!!\nYou got $(first(size(sigmatrix)))-elements in sigmatrix and $(first(size(delay_time)))-elements in delay_time! ")
            elseif last(size(sigmatrix)) !== first(size(probe_wavenumbers))
                error("Dimensions of sigmatrix and probe_wavenumbers must match!!\nYou got $(last(size(sigmatrix)))-elements in sigmatrix and $(first(size(probe_wavenumbers)))-elements in probe_wavenumbers!")
            elseif first(size(sig_bleaches[1])) !== first(size(delay_time))
                error("Dimensions of sig_bleaches[i] and delay_time must match!!\nYou got $(first(size(sig_bleaches[1])))-elements in sig_bleaches[i] and $(first(size(delay_time)))-elements in delay_time! ")
            else

                g4 = g_create(g3, pump_resonance)
                g5 = g_create(g4, dashboard_date)

                g6 = g_create(g5, "Data")
                g6["sig_matrix"] = sigmatrix
                g6["ref_matrix"] = refmatrix
                g6["pump_wavenumber"] = ekspla_wavenumber
                g6["wavenumber"] = probe_wavenumbers
                g6["dltime"] = delay_time

                for i in 1:length(mode_name)
                    g6["sig_mean_$(mode_name[i])"] = sig_bleaches[i]
                end

                for i in 1:length(mode_name)
                    g6["ref_mean_$(mode_name[i])"] = ref_bleaches[i]
                end

                g6["comment"] = get_metadata(raw_spectra[1])["comment"][1]*add_comment
                g6["folder_name"] = directory
            end

        else scan_type == "wavenumber_scan"

            if pump_wavenumbers === nothing
                pump_wavenumbers = [10^7 ./ get_metadata(raw_spectra[i])["ekspla laser"]["ekspla wavelength"][1] for i in 1:length(raw_spectra)]
            end

            if first(size(sigmatrix)) !== first(size(pump_wavenumbers))
                error("Dimensions of sigmatrix and pump_wavenumbers must match!!\nYou got $(first(size(sigmatrix)))-elements in sigmatrix and $(first(size(pump_wavenumbers)))-elements in pump_wavenumbers! ")
            elseif last(size(sigmatrix)) !== first(size(probe_wavenumbers))
                error("Dimensions of sigmatrix and probe_wavenumbers must match!!\nYou got $(last(size(sigmatrix)))-elements in sigmatrix and $(first(size(probe_wavenumbers)))-elements in probe_wavenumbers!")
            elseif first(size(sig_bleaches[1])) !== first(size(pump_wavenumbers))
                error("Dimensions of sig_bleaches[i] and pump_wavenumbers must match!!\nYou got $(first(size(sig_bleaches[1])))-elements in sig_bleaches[i] and $(first(size(pump_wavenumbers)))-elements in pump_wavenumbers! ")
            else

                g4 = g_create(g3, dashboard_date)

                g5 = g_create(g4, "Data")
                g5["sig_matrix"] = sigmatrix
                g5["ref_matrix"] = refmatrix
                g5["pump_wavenumber"] = pump_wavenumbers # Ekspla pump wavenumber for wavenumber scan
                g5["wavenumber"] = probe_wavenumbers
                g5["dltime"] = delay_time

                for i in 1:length(mode_name)
                    g5["sig_mean_$(mode_name[i])"] = sig_bleaches[i]
                end

                for i in 1:length(mode_name)
                    g5["ref_mean_$(mode_name[i])"] = ref_bleaches[i]
                end

                g5["comment"] = get_metadata(raw_spectra[1])["comment"][1]*add_comment
                g5["folder_name"] = directory
            end
        end
    end
end

"""
Save Spectra in a Matlab file.
"""
# function save_mat(filename::AbstractString, s::SFSpectrum)
#     save_mat(filename, SFSpectrum[s])
# end
#
# function save_mat(filename::AbstractString, s::Array{SFSpectrum})
#     # Check if filename has a proper extension
#     if splitext(filename)[end] != ".mat"
#         filename *= ".mat"
#     end
#
#     f = matopen(filename, "w")
#
#     length(s)
#
#     for i = 1:length(s)
#
#         # Put Stuff in dict
#         name = ""
#         try
#             name = get_attribute(s[i], "name")
#         catch
#             name = string(s[i].id)
#         end
#
#         data = Dict()
#         try
#             data = Dict(
#                 "name" => get_attribute(s[i], "name"),
#                 "signal" => s[i].s,
#                 "wavelength" => get_ir_wavelength(s[i]),
#                 "wavenumber" => get_ir_wavenumber(s[i]),
#             )
#         catch
#             data = Dict(
#                 "name" => name,
#                 "signal" => s[i].s
#             )
#         end
#         data
#
#         # Write to file
#         write(f, "data$i", data)
#     end
#     close(f)
# end


text:
module TestTearing

println("TestTearing: Tests tearing algorithm of the symbolic handling.")

using Modia

# Desired:
#   using ModiaMath: plot
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

# Tearing1
# No tearing is performed since coefficients of the block equations are not 1 or -1.
@model Tearing1 begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(start=1.0)
@equations begin
     2.0*x1 + 3.0*x2 = x3
    -1.4*x1 + 5.0*x2 = 2.0*x3
    der(x3) = -x3
    end
end
result = simulate(Tearing1, 1.0; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3"))

# Tearing1B
#
@model Tearing1B begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(start=1.0)
@equations begin
     2.0*x1 + x2 = x3
    -1.4*x1 + 5.0*x2 = 2.0*x3
    der(x3) = -x3
    end
end
result = simulate(Tearing1B, 1.0; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3"))

# Tearing2
# A correct solution is:
#    x2 := (x3 - 1.0) - 2.0 * (sin)(x1)
#    -1.4*sin(x1) + 5.0*cos(x2) = 2.0*x3
#     der(x3) = -x3
#
@model Tearing2 begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(start=1.0)
@equations begin
     2.0*sin(x1) + x2 = x3-1.0
    -1.4*sin(x1) + 5.0*cos(x2) = 2.0*x3
    der(x3) = -x3
    end
end
result = simulate(Tearing2, 0.3; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3"))

# Tearing3
# Variables that are explicitely solved due to tearing: x1, x2
#
@model Tearing3 begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(size=())
    x4 = Float(start=1.0)
@equations begin
     2.0*sin(x1) + x2 = x3-1.0 - 0.01*sin(x4)
    -1.4*sin(x1) + 5.0*cos(x2) + (0.02*x3)^3 = 2.0*x4
     sin(0.1*x4) = abs(x3) - x1
    der(x4) = -x4
    end
end
result = simulate(Tearing3, 0.3; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3", "x4"))

# Tearing4
# Variables that are explicitely solved due to tearing: x1, x2
#
@model Tearing4 begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(size=())
    x4 = Float(size=())
    x5 = Float(start=1.0)
@equations begin
     2.0*sin(x1) + x2 = x3-1.0 - 0.01*sin(x4)
    -1.4*sin(x1) + 5.0*cos(x2) + (0.02*x3)^3 = 2.0*x4
     sin(0.1*x4) = abs(x3) - x1
     -2*x4 + 3.0*x1 -3*x2 = 0.0
    der(x5) = -x5
    end
end
result = simulate(Tearing4, 0.3; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3", "x4", "x5"))

@model TearingCombined begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(start=1.0)

    x11 = Float(size=())
    x12 = Float(size=())
    x13 = Float(start=1.0)

    x21 = Float(size=())
    x22 = Float(size=())
    x23 = Float(start=1.0)

    x31 = Float(size=())
    x32 = Float(size=())
    x33 = Float(start=1.0)
    x34 = Float(start=1.0)
    x35 = Float(size=())
@equations begin
     2.0*x1 + 3.0*x2 = x3
    -1.4*x1 + 5.0*x2 = 2.0*x3
    der(x3) = -x3

     2.0*x11 + x12 = x13
    -1.4*x11 + 5.0*x12 = 2.0*x13
    der(x13) = -x13

     2.0*sin(x21) + x22 = x23-1.0
    -1.4*sin(x21) + 5.0*cos(x22) = 2.0*x23
    der(x23) = -x23

     2.0*sin(x31) + x32 = x33-1.0
    -1.4*sin(x35) + 5.0*cos(x34) = 2.0*x33
    x34 = 1.0*x32
    2*x35 = 2*x31
    der(x33) = -x33
    end
end
result = simulate(TearingCombined, 1.0; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3"))

# Tearing5
# RemoveAuxiliary should remove equation "2*x1 = 4*x2"
# Tearing should then eliminate x1 = 2*x2
#
# The log shows, that all this works. However, the generated code
# has a state vector of [x1, x2, x3, x4]. Since tearing can compute x1 from x2,
# it would be possible that the state vector has only elements [x2,x3,x4].
# Otherwise, tearing has basically no effect, because the system of equations
# was not reduced in the generated code.
@model Tearing5 begin
    x1 = Float(size=())
    x2 = Float(size=())
    x3 = Float(size=())
    x4 = Float(start=1.0)
@equations begin
     x1 = 2*x2
   2*x1 = 4*x2
     x1 + x2 = x3 + x4
     3.1*x2 + 1.2*x3 + 1.3*x4 = 0
    der(x4) = -x4
    end
end
result = simulate(Tearing5, 0.3; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("x1", "x2", "x3", "x4"), figure=5)

# Tearing6
@model Tearing6 begin
    C1 = 1e-3
    C2 = 2e-3
    u1 = Float(size=(), start=1.0)
    u2 = Float(size=(), state=false)
    i1 = Float(size=())
    v1 = Float(size=())
    v0 = Float(size=())
@equations begin
    C1*der(u1) = i1
    C2*der(u2) = -i1
    u1 = v1 - v0
    u2 = v1 - v0
    v0 = 0
    end
end
result = simulate(Tearing6, 1.0; logTranslation=true, logSimulation=true, tearing=true, removeSingularities=true)
plot(result, ("u1", "u2", "i1", "v1", "v0"), figure=6)

end

ext: jl


text:
module Functional

using Purple

function fn(x::Int)
    y = [ 0.0, 0.0, 0.0]
    return map(e -> x + e^2, y)
end

src = lift(fn,
           Tuple{Int};
           clean = true,
           jit = false,
           opt = false)
display(src)

src = lift(fn,
           Tuple{Int};
           jit = false,
           opt = false)
display(src)

end # module


text:
mutable struct Network
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
end

function Network(ns::Array{Int64, 1})
	nl = length(ns)
	Network(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[rand(ns[i])/10 for i in 2:nl])
end

#network with sparse (hacked) connectivy
mutable struct Network_sparseconnect
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	w_connect::Array{Array{Float64, 2}, 1} #connectivity matrix with only ones or zeros
	b::Array{Array{Float64, 1}, 1}
end
function Network_sparseconnect(ns::Array{Int64, 1})
	nl = length(ns)
	Network_sparseconnect(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[ones(ns[i+1], ns[i]) for i in 1:nl - 1],
			[rand(ns[i])/10 for i in 2:nl])
end

#init for sigmoid nonlinearity for reproduction of Lillicrap paper
function Network_LC(ns::Array{Int64, 1})
	nl = length(ns)
	Network(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[(rand(ns[i+1], ns[i])-0.5)/2. for i in 1:nl - 1],#2
			[(rand(ns[i])-0.5)/10. for i in 2:nl]) #-0.5
end

#init for Oja's rule (no biases)
function Network_Oja(ns::Array{Int64, 1})
	nl = length(ns)
	Network(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[zeros(ns[i]) for i in 2:nl])
end

#  intialize network (wi, bi) with random projections with:
#- learned weights (w) and offsets (b) by backprop for the first weights
#- random weights and offsets in output layer
function Network(ns::Array{Int64, 1},w::Array{Array{Float64, 2}, 1},b::Array{Array{Float64, 1}, 1})
	nl = length(ns)
	wi = [w[i] for i in 1:nl - 2]
	push!(wi,randn(ns[nl], ns[nl-1])/(10*sqrt(ns[nl-1])))
	bi = [b[i] for i in 1:nl - 2]
	push!(bi,rand(ns[nl])/10)
	Network(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			wi,
			bi)
end

#Network initialized with principal components (stored in w and b)
#ONLY ONE HIDDEN LAYER!!!
function Network(ns::Array{Int64, 1},w::Array{Float64, 2},b::Array{Float64, 1})
	nl = length(ns)
	wi = [w]
	push!(wi,randn(ns[nl], ns[nl-1])/(10*sqrt(ns[nl-1])))
	bi = [b]
	#push!(bi,zeros(ns[nl]))
	push!(bi,rand(ns[nl])/10)
	Network(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			wi,
			bi)
end

#Network with random feedback weights and offsets for "backprop"
mutable struct Network_RF
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	w_rf::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
end

function Network_RF(ns::Array{Int64, 1})
	nl = length(ns)
	Network_RF(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1], #feedback
			[rand(ns[i])/10 for i in 2:nl])
end

mutable struct Network_RF_sparseconnect
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	w_rf::Array{Array{Float64, 2}, 1}
	w_connect::Array{Array{Float64, 2}, 1} #connectivity matrix with only ones or zeros
	b::Array{Array{Float64, 1}, 1}
end
function Network_RF_sparseconnect(ns::Array{Int64, 1})
	nl = length(ns)
	Network_RF_sparseconnect(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1], #feedback
			[ones(ns[i+1], ns[i]) for i in 1:nl - 1], #connectivity
			[rand(ns[i])/10 for i in 2:nl])
end

#network with random feedback for reproduction of Lillicrap paper
function Network_RF_LC(ns::Array{Int64, 1})
	nl = length(ns)
	Network_RF(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[(rand(ns[i+1], ns[i])-0.5)/2 for i in 1:nl - 1],
			[(rand(ns[i+1], ns[i])-0.5)/1 for i in 1:nl - 1], #feedback
			[(rand(ns[i])-0.5)/2 for i in 2:nl])
end

#Network with weights respecting Dale's law
mutable struct Network_Dale
	nl::Int64
	ntypes::Array{Array{Bool, 1}, 1} #determines if neuron is excit. (true) or inhib. (false)
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
end

#function for initializating a Network respecting Dale's law
#type_ratio: n_inhibit/n_all for all layers
function Network_Dale(ns::Array{Int64, 1}, inhib_ratio::Array{Float64, 1})
	nl = length(ns)
	#draw inhibitory neurons for each layer
	ntypes = Array{Array{Bool, 1}, 1}()
	wi = Array{Array{Float64, 2}, 1}()
	for i in 1:nl-1
		inhibitory_neurons = [inhib_ratio[i] < rand(1)[1] for j in 1:ns[i]]
		wi_temp = randn(ns[i+1], ns[i])/(10*sqrt(ns[i]))
		for j in 1:ns[i]
			if inhibitory_neurons[j]
				wi_temp[:,j] = abs(wi_temp[:,j]) #excitatory -> all weights positive
			else
				wi_temp[:,j] = - 1. * abs(wi_temp[:,j]) #inhibitory -> all weights negative
			end
		end
		push!(ntypes,inhibitory_neurons)
		push!(wi,wi_temp)
	end
	Network_Dale(nl - 1,
			ntypes,
		  [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			wi,
			[rand(ns[i])/10 for i in 2:nl])
end

#Network with random feedback weights and offsets for "backprop" AND DALE'S LAW!
mutable struct Network_Dale_RF
	nl::Int64
	ntypes::Array{Array{Bool, 1}, 1} #determines if neuron is excit. (true) or inhib. (false)
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	w_rf::Array{Array{Float64, 2}, 1} #feedback
	b::Array{Array{Float64, 1}, 1}
end

function Network_Dale_RF(ns::Array{Int64, 1}, inhib_ratio::Array{Float64, 1})
	nl = length(ns)
	#draw inhibitory neurons for each layer
	ntypes = Array{Array{Bool, 1}, 1}()
	wi = Array{Array{Float64, 2}, 1}()
	for i in 1:nl-1
		inhibitory_neurons = [inhib_ratio[i] < rand(1)[1] for j in 1:ns[i]]
		wi_temp = randn(ns[i+1], ns[i])/(10*sqrt(ns[i]))
		for j in 1:ns[i]
			if inhibitory_neurons[j]
				wi_temp[:,j] = abs(wi_temp[:,j]) #excitatory -> all weights positive
			else
				wi_temp[:,j] = - 1. *abs(wi_temp[:,j]) #inhibitory -> all weights negative
			end
		end
		push!(ntypes,inhibitory_neurons)
		push!(wi,wi_temp)
	end
	Network_Dale_RF(nl - 1,
			ntypes,
		  [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			wi,
			wi, #feedback
			[rand(ns[i])/10 for i in 2:nl])
end

#Network for sparse coding with "activity memory" A,B
mutable struct Network_sparse
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
	A::Array{Array{Float64, 2}, 1}
	B::Array{Array{Float64, 2}, 1}
end

#init for sparse coding network (no biases)
function Network_sparse(ns::Array{Int64, 1})
	nl = length(ns)
	Network_sparse(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i], ns[i]) for i in 2:nl],
			[zeros(ns[i], ns[i+1]) for i in 1:nl - 1])
end

mutable struct NetworkWTA
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	ax_tr::Array{Array{Float64, 1}, 1} # low pass filtered internal variable
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
end

function NetworkWTA(ns::Array{Int64, 1})
	nl = length(ns)
	NetworkWTA(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[rand(ns[i])/10 for i in 2:nl])
end

mutable struct config
	task::String
	hidden_sizes::Array{Int64, 1}
  n_inits::Int64
	iterations::Int64
	learningrates::Array{Float64, 1}
	weight_decay::Float64
	nonlinearity::String
end

mutable struct config_iter
	task::String
	hidden_sizes::Array{Int64, 1}
  n_inits::Int64
	iterations::Array{Int64, 1}
	learningrates::Array{Float64, 1}
	weight_decay::Float64
	nonlinearity::String
end

mutable struct config_iterations
	task::String
	hidden_sizes::Array{Int64, 1}
  n_inits::Int64
	iterations::Array{Int64, 1}
	n_sampling::Int64
	learningrates::Array{Float64, 1}
	weight_decay::Float64
	nonlinearity::String
end

mutable struct config_sparse
	task::String
	hidden_sizes::Array{Int64, 1}
	n_inits::Int64
	iterations::Int64
	learningrates::Array{Float64, 1}
	lambda::Array{Float64, 1}
	weight_decay::Float64
	nonlinearity::String
end

#############################################################
#Data Import

using Knet
include(Knet.dir("data", "cifar.jl"))
include(Knet.dir("data", "mnist.jl"))
function import_data(dataset::String)
	if dataset == "CIFAR10"
		smallimgs, labels, smallimgstest, labelstest = cifar10()
	    smallimgs = reshape(smallimgs, 32^2 * 3, 50000)
	    smallimgstest = reshape(smallimgstest, 32^2 * 3, 10000)
	elseif dataset == "MNIST"
		smallimgs, labels, smallimgstest, labelstest = mnist()
	    smallimgs = reshape(smallimgs, 28^2, 60000)
	    smallimgstest = reshape(smallimgstest, 28^2, 10000)
	end
	smallimgs, labels .- 1, smallimgstest, labelstest .- 1
end

#############################################################
# data is dxn array of n d-dimensional random vectors
using LinearAlgebra, Distributions
function subtractmean!(data)
        m = mean(data, dims=2)
        d, n = size(data)
        for i in 0:n-1
                BLAS.axpy!(-1., m, 1:d, data, i*d + 1: (i+1) * d)
        end
end
function subtractmean(data)
        m = deepcopy(data)
        subtractmean!(m)
        m
end

# data is dxn array of n d-dimensional random vectors with mean 0
function whiten(data; method = :ZCA)
        f = svd(data)
        eval(Symbol(method, "Whitening"))(f) * sqrt(size(data, 2) - 1)
end

ZCAWhitening(f::LinearAlgebra.SVD) = f.U * f.Vt
PCAWhitening(f::LinearAlgebra.SVD) = f.U

#####################################################
#Helpers

function getsmallimg()
    global patternindex = rand(1:size(smallimgs)[2])
    smallimgs[:, patternindex]
end

function getsample()
    global patternindex = rand(1:n_samples)
    smallimgs[:, patternindex]
end

function getlabel(x)
    [labels[patternindex] == i for i in 0:9]
end

function gethiddenrep()
    global patternindex = rand(1:60000)
    hidden_rep_train[hidden_size_index][:, patternindex]
end

function gethiddensample()
    global patternindex = rand(1:n_samples)
    hidden_rep_train[hidden_size_index][:, patternindex]
end

function geterrors(net, imgs, labels; nonlinearity = nonlinearity, getwrongindices = false)
	error = 0
	noftest = size(imgs)[2]
	if getwrongindices
		wrongindices = []
		for i in 1:noftest
			net.x[1] = imgs[:, i]
			forwardprop!(net, nonlinearity = nonlinearity)
			if findmax(net.x[end])[2] != labels[i] + 1
				error += 1
				push!(wrongindices,i)
			end
		end
		return error/noftest, wrongindices
	else
		for i in 1:noftest
			net.x[1] = imgs[:, i]
			forwardprop!(net, nonlinearity = nonlinearity)
			error += findmax(net.x[end])[2] != labels[i] + 1
		end
		error/noftest
	end
end

# reduce full random connectivity to patches with random connectivity
# by multiplying with mask of linear dimension "patch_size"
# mode:
# - spherical: random vectors for patches lie on the surface of a hypersphere
# - randn: conventional variance conserving initialisation with gaussian random number
function _rescale_weights!(net, i, input_dim, patch_size, mode; colorimage = false)
	if mode == "randn"
		net.w[1][i,:] .*= input_dim./(patch_size) #rescale weights for same input strength as full connectivity
		!colorimage && (net.w[1][i,:] .* 3)
	elseif mode == "spherical"
		net.w[1][i,:] .*= 1. / norm(net.w[1][i,:]) # normalize random patches: since sampled from randn they lie on the surface of a hypersphere now! norm grows as sqrt(eff. in-dim)
		net.w[1][i,:] .*= 3.#sqrt(patch_size) # heuristic rescaling for best performance
	end
end
function set_connectivity!(net,patch_size; mode = "randn", colorimage = false)
	if colorimage
		  input_size = Int(length(net.x[1])/3)
	else
		  input_size = Int(length(net.x[1]))
	end
	input_dim = Int(sqrt(input_size))
	nr_hidden_neurons = length(net.x[2])
	mask = zeros(input_dim,input_dim) #initialize mask
	mask[1:patch_size,1:patch_size] .= 1.
	if colorimage
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			for j in 1:3
			  net.w[1][i,(j-1)*input_size+1:j*input_size] = reshape(reshape(net.w[1][i,(j-1)*input_size+1:j*input_size],input_dim,input_dim).*circshift(mask,shifts),input_size)
			  _rescale_weights!(net, i, input_dim, patch_size, mode; colorimage = colorimage)
			end
		end
	else
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			net.w[1][i,:] = reshape(reshape(net.w[1][i,:],input_dim,input_dim).*circshift(mask,shifts),input_size)
			_rescale_weights!(net, i, input_dim, patch_size, mode; colorimage = colorimage)
		end
	end
end
function set_connectivity_gabor!(net, patchsize; colorimage = false,
      λrange = [patchsize/4,2*patchsize], ψupperbound = 2π,
      σrange = [patchsize/8,patchsize], γrange = 0, amplitude = 1)
    if colorimage
          input_size = Int(length(net.x[1])/3)
    else
          input_size = Int(length(net.x[1]))
    end
	input_dim = Int(sqrt(input_size))
    hiddensize = length(net.x[2])
    if colorimage
        for i in 1:hiddensize
            mask = zeros(input_dim,input_dim)
            mask[1:patchsize,1:patchsize] =
                randMNISTgaborfilterBayesOpt(patchsize, λrange, ψupperbound, σrange, γrange, amplitude)
            shifts = rand(0:input_dim-patchsize,2)
            for j in 1:3
              net.w[1][i,(j-1)*input_size+1:j*input_size] = circshift(mask,shifts)[:]
            end
        end
    else
        for i in 1:hiddensize
            mask = zeros(input_dim,input_dim)
            mask[1:patchsize,1:patchsize] =
                randMNISTgaborfilterBayesOpt(patchsize, λrange, ψupperbound, σrange, γrange, amplitude)#randMNISTgaborfilter(patchsize)#randgaborfilter(patchsize)
            shifts = rand(0:input_dim-patchsize,2)
            net.w[1][i,:] = circshift(mask,shifts)[:]
        end
    end
end
function set_connectivity_matrix!(net,patch_size; mode = "randn", colorimage = false) # for Network_sparse & Network_RF_sparseconnect
	if colorimage
		input_size = Int(length(net.x[1])/3)
	else
		input_size = length(net.x[1])
	end
	input_dim = Int(sqrt(input_size))
	nr_hidden_neurons = length(net.x[2])
	mask = zeros(input_dim,input_dim) #initialize mask
	mask[1:patch_size,1:patch_size] .= 1.
	if colorimage
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			for j in 1:3
				net.w[1][i,(j-1)*input_size+1:j*input_size] = reshape(reshape(net.w[1][i,(j-1)*input_size+1:j*input_size],input_dim,input_dim).*circshift(mask,shifts),input_size)
				net.w_connect[1][i,(j-1)*input_size+1:j*input_size] = clamp.(ceil.(abs.(net.w[1][i,(j-1)*input_size+1:j*input_size])),0,1)
			end
			_rescale_weights!(net, i, input_dim, patch_size, mode; colorimage = colorimage)
		end
	else
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			net.w[1][i,:] = reshape(reshape(net.w[1][i,:],input_dim,input_dim).*circshift(mask,shifts),input_size)
			net.w_connect[1][i,:] = clamp.(ceil.(abs.(net.w[1][i,:])),0,1)
			_rescale_weights!(net, i, input_dim, patch_size, mode; colorimage = colorimage)
		end
	end
end

"""
    gaborfilter
"""
function gaborfilter(λ, θ, ψ, σ, γ, N)
    w = zeros(N, N)
    for x in 1:N
        for y in 1:N
            xᴵ, yᴵ = [cos(θ) sin(θ); -sin(θ) cos(θ)] * ([x; y] .- (1 + (N-1)/2))
            w[x, y] = exp(-(xᴵ^2 + γ^2 * yᴵ^2)/(2σ^2)) * cos(2π * xᴵ/λ + ψ)
        end
    end
    w
end
"""
    randgaborfilter
"""
function randgaborfilter(N)
    w = gaborfilter(2N*rand(), 2π*rand(), 2π*rand(), N*rand(), N*rand(), N)
    w ./ sum(abs.(w))
end
function randMNISTgaborfilter(N)
    w = gaborfilter((2N-N/4)*rand()+N/4, 2π*rand(), 2π*rand(), (N-N/8)*rand()+N/8, 1, N)
    w ./ sum(abs.(w))
end
function randMNISTgaborfilterBayesOpt(N, λrange, ψupperbound, σrange, γrange, amplitude)
    w = gaborfilter(λrange[1] + (λrange[2] - λrange[1]) * rand(),
                    2π*rand(),
                    ψupperbound * rand(),
                    σrange[1] + (σrange[2] - σrange[1]) * rand(),
                    1 - γrange * rand(), N)
    w .* amplitude / sum(abs.(w))
end

# reduce full random connectivity to sparse random connectivity
# by connecting with Bernoulli-sampling
function set_connectivity_sparse!(net,in_degree)
	input_size = length(net.x[1])
	if in_degree != 1.
		nr_hidden_neurons = length(net.x[2])
		for i in 1:nr_hidden_neurons
			no_connection = rand(Bernoulli(in_degree),input_size)
			net.w[1][i,:] .*= no_connection*3/sqrt(in_degree) #rescale weights for same input strength as full connectivity
		end
	end
end

#re_initialize random weights differently for RP
function re_initialize!(net; w_bool = false,b_bool = true)
	if w_bool
		net.w[1] = (rand(size(net.w[1])[1], size(net.w[1])[2])-0.5)/sqrt(size(net.w[1])[2]) #(10*sqrt(size(net.w[1])[2]))
	elseif b_bool
		net.b[1] = zeros(size(net.w[1])[1]) #reference = rand(size(net.w[1])[1])/10
	end
end

#####################################################

#forwardprop for different nonlinearities
function forwardprop!(net; nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		BLAS.gemv!('N', 1., net.w[i], net.x[i], 0., net.ax[i])
		BLAS.axpy!(1., net.b[i], net.ax[i])
		nonlinearity[i](net.ax[i], net.x[i+1])
	end
end

#forwardprop for random feedback weights (functionality = same)
function forwardprop!(net::Network_RF; nonlinearity = Array{Function, 1})
	for i in 1:net.nl
		BLAS.gemv!('N', 1., net.w[i], net.x[i], 0., net.ax[i])
		BLAS.axpy!(1., net.b[i], net.ax[i])
		nonlinearity[i](net.ax[i], net.x[i+1])
	end
end

Z(x) = map(f -> (f(pi*x[1]) + f(pi*(x[1]+x[2])) + 2)/4, [cos; sin])

function heaviside!(x::Array{Float64, 1})
	for i in 1:length(x)
		x[i] = x[i] > 0
	end
end
function heaviside!(inp, outp, threshold = zeros(inp))
	for i in 1:length(inp)
		outp[i] = Int(inp[i] >= threshold[i])
	end
end
function heaviside!(inp, outp; threshold = 0.)
	for i in 1:length(inp)
		outp[i] = Int(inp[i] >= threshold)
	end
end

function relu!(inp, outp)
	for i in 1:length(inp)
		outp[i] = max(0, inp[i])
	end
end

function relu_diff!(ax, error)
	for i in 1:length(ax)
		error[i] *= ax[i] > 0
	end
end

function lin!(inp, outp)
	for i in 1:length(inp)
		outp[i] = inp[i]
	end
end

function lin_diff!(ax, error)
end

function tanh!(x::Array{Float64, 1})
	for i in 1:length(x)
		x[i] = tanh(x[i])
	end
end

function tanh_diff!(ax, error)
	for i in 1:length(ax)
		error[i] *= sech(ax[i])^2
	end
end

function sigm!(inp,outp)
	for i in 1:length(inp)
		outp[i] = 1. / (1. + exp(-inp[i]))
	end
end

function sigm!(x::Array{Float64, 1})
	for i in 1:length(x)
		x[i] = 1. / (1. + exp(-x[i]))
	end
end

function sigm_diff!(ax, error)
	for i in 1:length(ax)
		error[i] *= 1. / (1 + exp(-ax[i])) * (1. - 1. / (1 + exp(-ax[i])))
	end
end

#for different nonlinearities
function backprop!(net::Network, learningrate, weight_decay;
				   nonlinearity_diff = Array{Function, 1})
	for i in net.nl:-1:2
		nonlinearity_diff[i](net.ax[i], net.e[i])
		BLAS.gemv!('T', 1., net.w[i], net.e[i], 0., net.e[i-1])
	end
	nonlinearity_diff[1](net.ax[1], net.e[1])
	for i in 1:net.nl
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
		if weight_decay != 0.
			scale!(net.w[i],1. - weight_decay)
			scale!(net.b[i],1. - weight_decay)
		end
	end
end

#for different nonlinearities
function backprop!(net::Network_sparseconnect, learningrate, weight_decay;
				   nonlinearity_diff = Array{Function, 1})
	for i in net.nl:-1:2
		nonlinearity_diff[i](net.ax[i], net.e[i])
		BLAS.gemv!('T', 1., net.w[i], net.e[i], 0., net.e[i-1])
	end
	nonlinearity_diff[1](net.ax[1], net.e[1])
	for i in 1:net.nl
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		net.w[i] .*= net.w_connect[i]
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
		if weight_decay != 0.
			scale!(net.w[i],1. - weight_decay)
			scale!(net.b[i],1. - weight_decay)
		end
	end
end

#backprop with random feedback WITH derivative of activation function
#nonlin_diff: Use or don't use diff for calculating errors
function backprop!(net::Network_RF, learningrate, weight_decay;
				   nonlinearity_diff = Array{Function, 1},
					 nonlin_diff = Bool(1))
	if nonlin_diff
		for i in net.nl:-1:2
			nonlinearity_diff[i](net.ax[i], net.e[i])
			BLAS.gemv!('T', 1., net.w_rf[i], net.e[i], 0., net.e[i-1])
		end
		nonlinearity_diff[1](net.ax[1], net.e[1])
	else
		for i in net.nl:-1:2
			#nonlinearity_diff(net.ax[i], net.e[i])
			BLAS.gemv!('T', 1., net.w_rf[i], net.e[i], 0., net.e[i-1])
		end
		#nonlinearity_diff(net.ax[1], net.e[1])
	end
	for i in 1:net.nl
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
		if weight_decay != 0.
			scale!(net.w[i],1. - weight_decay)
			scale!(net.b[i],1. - weight_decay)
		end
	end
end
function backprop!(net::Network_RF_sparseconnect, learningrate, weight_decay;
				   nonlinearity_diff = Array{Function, 1},
					 nonlin_diff = Bool(1))
	if nonlin_diff
		for i in net.nl:-1:2
			nonlinearity_diff[i](net.ax[i], net.e[i])
			BLAS.gemv!('T', 1., net.w_rf[i], net.e[i], 0., net.e[i-1])
		end
		nonlinearity_diff[1](net.ax[1], net.e[1])
	else
		for i in net.nl:-1:2
			#nonlinearity_diff(net.ax[i], net.e[i])
			BLAS.gemv!('T', 1., net.w_rf[i], net.e[i], 0., net.e[i-1])
		end
		#nonlinearity_diff(net.ax[1], net.e[1])
	end
	for i in 1:net.nl
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		net.w[i] .*= net.w_connect[i]
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
		if weight_decay != 0.
			scale!(net.w[i],1. - weight_decay)
			scale!(net.b[i],1. - weight_decay)
		end
	end
end

#backprop for random projections and different nonlinearities
function backprop_rp!(net::Network, learningrate;
				   nonlinearity_diff = Array{Function, 1})
	for i in net.nl:-1:net.nl
		nonlinearity_diff[i](net.ax[i], net.e[i])
		BLAS.gemv!('T', 1., net.w[i], net.e[i], 0., net.e[i-1])
	end
	nonlinearity_diff[1](net.ax[1], net.e[1])
	for i in net.nl:net.nl
		BLAS.ger!(learningrate[i], net.e[i], net.x[i], net.w[i])
		BLAS.axpy!(learningrate[i], net.e[i], net.b[i])
	end
end

function getlrates(layerdims)
	1 ./ sqrt(layerdims[1:end-1] + 1)
end

using ProgressMeter

#target must be one-hot
function _seterror_mse!(net, target)
	net.e[end] = target - net.x[end]
end

function _loss_mse(net, target)
	return norm(net.e[end])
end

function _softmax(input)
	input = deepcopy(input)
	exps = exp.(input .- maximum(input))
  return exps / sum(exps)
end

#target must be one-hot
function _seterror_crossentropysoftmax!(net, target)
	probs = _softmax(net.x[end])
	net.e[end] = target - probs
end

function _loss_crossentropy(net, target)
	probs = _softmax(net.x[end])
	return -target'*log.(probs)
end

#learn function for different nonlinearities
function learn!(net::Network,
				inputfunction::Function,
				targetfunction::Function,
				iterations::Int64,
				learningrates::Array{Float64, 1},
				weight_decay::Float64;
				nonlinearity = Array{Function, 1},
				nonlinearity_diff = Array{Function, 1},
				nonlin_diff = true,
				lossfunction = _loss_mse,
				set_error_lossderivative! = _seterror_mse!)
	losses = Float64[]
	loss = 0.
	@showprogress for i in 1:iterations
		net.x[1] = inputfunction()
		forwardprop!(net, nonlinearity = nonlinearity)
		target = targetfunction(net.x[1])
		set_error_lossderivative!(net, target)#net.e[end] = target - net.x[end]
		loss += lossfunction(net, target)#loss += norm(net.e[end])
		if i % div(iterations, 100) == 0
			push!(losses, loss/div(iterations, 100))
			loss = 0.
		end
		backprop!(net, learningrates, weight_decay, nonlinearity_diff = nonlinearity_diff)
	end
	losses
end

function learn!(net::Network_sparseconnect,
				inputfunction::Function,
				targetfunction::Function,
				iterations::Int64,
				learningrates::Array{Float64, 1},
				weight_decay::Float64;
				nonlinearity = Array{Function, 1},
				nonlinearity_diff = Array{Function, 1},
				nonlin_diff = true,
				lossfunction = _loss_mse,
				set_error_lossderivative! = _seterror_mse!)
	losses = Float64[]
	loss = 0.
	@showprogress for i in 1:iterations
		net.x[1] = inputfunction()
		forwardprop!(net, nonlinearity = nonlinearity)
		target = targetfunction(net.x[1])
		set_error_lossderivative!(net, target)#net.e[end] = target - net.x[end]
		loss += lossfunction(net, target)#loss += norm(net.e[end])
		if i % div(iterations, 100) == 0
			push!(losses, loss/div(iterations, 100))
			loss = 0.
		end
		backprop!(net, learningrates, weight_decay, nonlinearity_diff = nonlinearity_diff)
	end
	losses
end

#learn for random projections and different nonlinearities
function learn_rp!(net::Network,
				inputfunction::Function,
				targetfunction::Function,
				iterations::Int64,
				learningrates::Array{Float64, 1},
				weight_decay::Float64;
				nonlinearity = Array{Function, 1},
				nonlinearity_diff = Array{Function, 1},
				nonlin_diff = Bool(1),
				lossfunction = _loss_mse,
				set_error_lossderivative! = _seterror_mse!)
	losses = Float64[]
	loss = 0.
	@showprogress for i in 1:iterations
		net.x[1] = inputfunction()
		forwardprop!(net, nonlinearity = nonlinearity)
		target = targetfunction(net.x[1])
		set_error_lossderivative!(net, target)#net.e[end] = target - net.x[end]
		loss += lossfunction(net, target)#loss += norm(net.e[end]))
		if i % div(iterations, 100) == 0
			push!(losses, loss/div(iterations, 100))
			loss = 0.
		end
		backprop_rp!(net, learningrates, nonlinearity_diff = nonlinearity_diff)
	end
	losses
end

#learn for random feedback weights
function learn!(net::Network_RF,
				inputfunction::Function,
				targetfunction::Function,
				iterations::Int64,
				learningrates::Array{Float64, 1},
				weight_decay::Float64;
				nonlinearity = Array{Function, 1},
				nonlinearity_diff = Array{Function, 1},
				nonlin_diff = Bool(1),
				lossfunction = _loss_mse,
				set_error_lossderivative! = _seterror_mse!)
	losses = Float64[]
	loss = 0.
	@showprogress for i in 1:iterations
		net.x[1] = inputfunction()
		forwardprop!(net, nonlinearity = nonlinearity)
		target = targetfunction(net.x[1])
		set_error_lossderivative!(net, target)#net.e[end] = target - net.x[end]
		loss += lossfunction(net, target)#loss += norm(net.e[end]))
		if i % div(iterations, 10) == 0
			push!(losses, loss/div(iterations, 10))
			loss = 0.
		end
		backprop!(net, learningrates, weight_decay, nonlinearity_diff = nonlinearity_diff,
		 					nonlin_diff = nonlin_diff)
	end
	losses
end

function learn!(net::Network_RF_sparseconnect,
				inputfunction::Function,
				targetfunction::Function,
				iterations::Int64,
				learningrates::Array{Float64, 1},
				weight_decay::Float64;
				nonlinearity = Array{Function, 1},
				nonlinearity_diff = Array{Function, 1},
				nonlin_diff = Bool(1),
				lossfunction = _loss_mse,
				set_error_lossderivative! = _seterror_mse!)
	losses = Float64[]
	loss = 0.
	@showprogress for i in 1:iterations
		net.x[1] = inputfunction()
		forwardprop!(net, nonlinearity = nonlinearity)
		target = targetfunction(net.x[1])
		set_error_lossderivative!(net, target)#net.e[end] = target - net.x[end]
		loss += lossfunction(net, target)#loss += norm(net.e[end]))
		if i % div(iterations, 10) == 0
			push!(losses, loss/div(iterations, 10))
			loss = 0.
		end
		backprop!(net, learningrates, weight_decay, nonlinearity_diff = nonlinearity_diff,
		 					nonlin_diff = nonlin_diff)
	end
	losses
end


text:
@doc doc"""
XY model with energy $E = -\sum_{ij} J_{ij} \cos(\theta_i - \theta_j)$,
where $\theta_i = 2\pi \sigma_i$ and $\sigma_i \in [0, 1)$.
"""
mutable struct XY <: Model
    lat :: Lattice
    spins :: Matrix{Float64}
    rng :: Random.MersenneTwister

    function XY(lat::Lattice, rng::Random.AbstractRNG)
        spins = rand(rng, 1, numsites(lat))
        return new(lat, spins, rng)
    end
end

XY(lat::Lattice) = XY(lat, Random.seed!(Random.MersenneTwister(0)))
XY(lat::Lattice, seed) = XY(lat, Random.seed!(Random.MersenneTwister(0), seed...))

@doc doc"""
   XY(param)

Generates `XY` using `param["Lattice"]`,  and `param["Seed"]` (if defined).
Each spin $\sigma_i$ will be initialized randomly and independently.
"""
function XY(param::Parameter)
    lat = generatelattice(param)
    if "Seed" in keys(param)
        return XY(lat, param["Seed"])
    else
        return XY(lat)
    end
end

include("update.jl")
include("estimator.jl")
include("postproc.jl")


text:
@testset "Model matrix" begin

    using StatsBase: StatisticalModel

    using SparseArrays, DataFrames, Tables

    sparsetype = SparseMatrixCSC{Float64,Int}

    d = DataFrame(y = 1:4,
                  x1 = 5:8,
                  x2 = 9:12,
                  x3 = 13:16,
                  x4 = 17:20)
    d.x1p = categorical(d.x1)

    d_orig = deepcopy(d)

    x1 = [5.:8;]
    x2 = [9.:12;]
    x3 = [13.:16;]
    x4 = [17.:20;]
    f = @formula(y ~ 1 + x1 + x2)
    mf = ModelFrame(f, d)
    @test coefnames(mf) == ["(Intercept)","x1","x2"]
    @test response(mf) == [1:4;]
    mm = ModelMatrix(mf)
    smm = ModelMatrix{sparsetype}(mf)
    @test mm.m[:,1] == ones(4)
    @test mm.m[:,2:3] == [x1 x2]
    @test mm.m == smm.m

    @test isa(mm.m, Matrix{Float64})
    @test isa(smm.m, sparsetype)

    f_implint = @formula(y ~ x1 + x2)
    @test ModelMatrix(ModelFrame(f_implint, d)).m == mm.m

    @test ModelMatrix(ModelFrame(f_implint, d, model=Nothing)).m == hcat(x1, x2)

    #test_group("expanding a nominal array into a design matrix of indicators for each dummy variable")

    mf = ModelFrame(@formula(y ~ 1 + x1p), d)
    mm = ModelMatrix(mf)

    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test coefnames(mf)[2:end] == ["x1p: 6", "x1p: 7", "x1p: 8"]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    #test_group("Creating a model matrix using full formulas: y => x1 + x2, etc")

    d = deepcopy(d_orig)
    f = @formula(y ~ 1 + x1 & x2)
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1.*x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    f = @formula(y ~ 1 + x1 * x2)
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1 x2 x1.*x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    d.x1 = categorical(x1)
    x1e = [[0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
    f = @formula(y ~ 1 + x1 * x2)
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1e x2 [0, 10, 0, 0] [0, 0, 11, 0] [0, 0, 0, 12]]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    # additional tests from Tom
    y = [1., 2, 3, 4]
    mf = ModelFrame(@formula(y ~ 1 + x2), d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m
    @test response(mf) == y''

    d = deepcopy(d_orig)
    d.x1 = CategoricalArray{Union{Missing, Float64}}(d.x1)

    f = @formula(y ~ 1 + x2 + x3 + x3*x2)
    mm = ModelMatrix(ModelFrame(f, d))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x3*x2 + x2 + x3), d))
    @test mm.m == [ones(4) x3 x2 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x1 + x2 + x3 + x4), d))
    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test mm.m[:,5] == x2
    @test mm.m[:,6] == x3
    @test mm.m[:,7] == x4

    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x2 + x3 + x4), d))
    @test mm.m == [ones(4) x2 x3 x4]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x2 + x2), d))
    @test mm.m == [ones(4) x2]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x2*x3 + x2&x3), d))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x2*x3*x4), d))
    @test mm.m == [ones(4) x2 x3 x4 x2.*x3 x2.*x4 x3.*x4 x2.*x3.*x4]
    mm = ModelMatrix(ModelFrame(@formula(y ~ 1 + x2&x3 + x2*x3), d))
    @test mm.m == [ones(4) x2 x3 x2.*x3]

    f = @formula(y ~ 1 + x2 & x3 & x4)
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2.*x3.*x4]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    f = @formula(y ~ 1 + x1 & x2 & x3)
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2:end] == Matrix(Diagonal(x2.*x3))
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    ## Distributive property of :& over :+
    d = deepcopy(d_orig)
    f = @formula(y ~ 1 + (x1+x2) & (x3+x4))
    mf = ModelFrame(f, d)
    mm = ModelMatrix(mf)
    @test mm.m == hcat(ones(4), x1.*x3, x1.*x4, x2.*x3, x2.*x4)
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    ## Condensing nested :+ calls
    f = @formula(y ~ 1 + x1 + (x2 + (x3 + x4)))
    @test ModelMatrix(ModelFrame(f, d)).m == hcat(ones(4), x1, x2, x3, x4)

    ## Extra levels in categorical column
    mf_full = ModelFrame(@formula(y ~ 1 + x1p), d)
    mm_full = ModelMatrix(mf_full)
    @test size(mm_full) == (4,4)

    mf_sub = ModelFrame(@formula(y ~ 1 + x1p), d[2:4, :])
    mm_sub = ModelMatrix(mf_sub)
    ## should have only three rows, and only three columns (intercept plus two
    ## levels of factor) - it also catches when some levels do not occur in the data
    @test size(mm_sub) == (3,3)

    ## Missing data
    d.x1m = [5, 6, missing, 7]
    mf = ModelFrame(@formula(y ~ 1 + x1m), d)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2] == d[completecases(d), :x1m]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    ## Same variable on left and right side
    mf = ModelFrame(@formula(x1 ~ x1), d)
    mm = ModelMatrix(mf)
    mm.m == response(mf)

    ## Promote non-redundant categorical terms to full rank
    @testset "non-redundant categorical terms" begin
        d = DataFrame(x = repeat([:a, :b], outer = 4),
                      y = repeat([:c, :d], inner = 2, outer = 2),
                      z = repeat([:e, :f], inner = 4))
        cs = Dict([Symbol(name) => EffectsCoding() for name in names(d)])
        d.n = 1.:8

        ## No intercept
        mf = ModelFrame(@formula(n ~ 0 + x), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == [1 0
                       0 1
                       1 0
                       0 1
                       1 0
                       0 1
                       1 0
                       0 1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["x: a", "x: b"]

        ## promotion blocked when we block default model=StatisticalModel
        mf = ModelFrame(@formula(n ~ 0 + x), d, model=Nothing, contrasts=cs)
        mm = ModelMatrix(mf)
        @test all(mm.m .== ifelse.(d.x .== :a, -1, 1))
        @test coefnames(mf) == ["x: b"]

        ## No first-order term for interaction
        mf = ModelFrame(@formula(n ~ 1 + x + x&y), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m[:, 2:end] == [-1 -1  0
                                 1  0 -1
                                 -1  1  0
                                 1  0  1
                                 -1 -1  0
                                 1  0 -1
                                 -1  1  0
                                 1  0  1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["(Intercept)", "x: b", "x: a & y: d", "x: b & y: d"]

        ## When both terms of interaction are non-redundant:
        mf = ModelFrame(@formula(n ~ 0 + x&y), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == [1 0 0 0
                       0 1 0 0
                       0 0 1 0
                       0 0 0 1
                       1 0 0 0
                       0 1 0 0
                       0 0 1 0
                       0 0 0 1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                                "x: a & y: d", "x: b & y: d"]

        # only a three-way interaction: every term is promoted.
        mf = ModelFrame(@formula(n ~ 0 + x&y&z), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == Matrix(1.0I, 8, 8)
        @test mm.m == ModelMatrix{sparsetype}(mf).m

        # two two-way interactions, with no lower-order term. both are promoted in
        # first (both x and y), but only the old term (x) in the second (because
        # dropping x gives z which isn't found elsewhere, but dropping z gives x
        # which is found (implicitly) in the promoted interaction x&y).
        mf = ModelFrame(@formula(n ~ 0 + x&y + x&z), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == [1 0 0 0 -1  0
                       0 1 0 0  0 -1
                       0 0 1 0 -1  0
                       0 0 0 1  0 -1
                       1 0 0 0  1  0
                       0 1 0 0  0  1
                       0 0 1 0  1  0
                       0 0 0 1  0  1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                                "x: a & y: d", "x: b & y: d",
                                "x: a & z: f", "x: b & z: f"]

        # ...and adding a three-way interaction, only the shared term (x) is promoted.
        # this is because dropping x gives y&z which isn't present, but dropping y or z
        # gives x&z or x&z respectively, which are both present.
        mf = ModelFrame(@formula(n ~ 0 + x&y + x&z + x&y&z), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == [1 0 0 0 -1  0  1  0
                       0 1 0 0  0 -1  0  1
                       0 0 1 0 -1  0 -1  0
                       0 0 0 1  0 -1  0 -1
                       1 0 0 0  1  0 -1  0
                       0 1 0 0  0  1  0 -1
                       0 0 1 0  1  0  1  0
                       0 0 0 1  0  1  0  1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                                "x: a & y: d", "x: b & y: d",
                                "x: a & z: f", "x: b & z: f",
                                "x: a & y: d & z: f", "x: b & y: d & z: f"]

        # two two-way interactions, with common lower-order term. the common term x is
        # promoted in both (along with lower-order term), because in every case, when
        # x is dropped, the remaining terms (1, y, and z) aren't present elsewhere.
        mf = ModelFrame(@formula(n ~ 0 + x + x&y + x&z), d, contrasts=cs)
        mm = ModelMatrix(mf)
        @test mm.m == [1 0 -1  0 -1  0
                       0 1  0 -1  0 -1
                       1 0  1  0 -1  0
                       0 1  0  1  0 -1
                       1 0 -1  0  1  0
                       0 1  0 -1  0  1
                       1 0  1  0  1  0
                       0 1  0  1  0  1]
        @test mm.m == ModelMatrix{sparsetype}(mf).m
        @test coefnames(mf) == ["x: a", "x: b",
                                "x: a & y: d", "x: b & y: d",
                                "x: a & z: f", "x: b & z: f"]

        ## FAILS: When both terms are non-redundant and intercept is PRESENT
        ## (not fully redundant). Ideally, would drop last column. Might make sense
        ## to warn about this, and suggest recoding x and y into a single variable.
        mf = ModelFrame(@formula(n ~ 1 + x&y), d[1:4, :], contrasts=cs)
        @test_broken ModelMatrix(mf).m == [1 1 0 0
                                           1 0 1 0
                                           1 0 0 1
                                           1 0 0 0]
        @test_broken coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                                       "x: a & y: d", "x: b & y: d"]

        ## note that R also does not detect this automatically. it's left to glm et al.
        ## to detect numerically when the model matrix is rank deficient, which is hard
        ## to do correctly.
        # > d = data.frame(x = factor(c(1, 2, 1, 2)), y = factor(c(3, 3, 4, 4)))
        # > model.matrix(~ 1 + x:y, d)
        #   (Intercept) x1:y3 x2:y3 x1:y4 x2:y4
        # 1           1     1     0     0     0
        # 2           1     0     1     0     0
        # 3           1     0     0     1     0
        # 4           1     0     0     0     1

    end

    @testset "arbitrary functions in formulae" begin
        d = deepcopy(d_orig)
        mf = ModelFrame(@formula(y ~ log(x1)), d, model=Nothing)
        @test coefnames(mf) == ["log(x1)"]
        mm = ModelMatrix(mf)
        @test all(mm.m .== log.(x1))

        # | is not special in base formula:
        d = DataFrame(x = [1,2,3], y = [4,5,6])
        mf = ModelFrame(@formula(y ~ 1 + (1 | x)), d)
        @test coefnames(mf) == ["(Intercept)", "1 | x"]

        mf = ModelFrame(@formula(y ~ 0 + (1 | x)), d)
        @test all(ModelMatrix(mf).m .== float.(1 .| d.x))
        @test coefnames(mf) == ["1 | x"]
    end

    # Ensure X is not a view on df column
    d = DataFrame(x = [1.0,2.0,3.0], y = [4.0,5.0,6.0])
    mf = ModelFrame(@formula(y ~ 0 + x), d)
    X = ModelMatrix(mf).m
    X[1] = 0.0
    @test mf.data[:x][1] === 1.0

    # Ensure string columns are supported
    d1 = DataFrame(A = 1:4, B = categorical(["M", "F", "F", "M"]))
    d2 = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])
    d3 = DataFrame(Any[1:4, ["M", "F", "F", "M"]], [:A, :B])

    M1 = ModelMatrix(ModelFrame(@formula(A ~ B), d1))
    M2 = ModelMatrix(ModelFrame(@formula(A ~ B), d2))
    M3 = ModelMatrix(ModelFrame(@formula(A ~ B), d3))

    @test (M1.m, M1.assign) == (M2.m, M2.assign) == (M3.m, M3.assign)

    @testset "row-wise model matrix construction" begin
        d = DataFrame(r = rand(8),
                      w = rand(8),
                      x = repeat([:a, :b], outer = 4),
                      y = repeat([:c, :d], inner = 2, outer = 2),
                      z = repeat([:e, :f], inner = 4))

        f = apply_schema(@formula(r ~ 1 + w*x*y*z), schema(d))
        modelmatrix(f, d)
        @test reduce(vcat, last.(modelcols.(Ref(f), Tables.rowtable(d)))') == modelmatrix(f,d)
    end

    @testset "modelmatrix and response set schema if needed" begin
        d = DataFrame(r = rand(8),
                      w = rand(8),
                      x = repeat([:a, :b], outer = 4),
                      y = repeat([:c, :d], inner = 2, outer = 2),
                      z = repeat([:e, :f], inner = 4))

        f = @formula(r ~ 1 + w*x*y*z)

        mm1 = modelmatrix(f, d)
        mm2 = modelmatrix(apply_schema(f, schema(d)), d)
        @test mm1 == mm2

        r1 = response(f, d)
        r2 = response(apply_schema(f, schema(d)), d)
        @test r1 == r2
    end

    @testset "setcontrasts!" begin
        @testset "#95" begin
            tbl = (Y = randn(8),
                   A = repeat(['N','Y'], outer=4),
                   B=repeat(['-','+'], inner=2, outer=2),
                   C=repeat(['L','H'], inner=4))

            contrasts = Dict(:A=>HelmertCoding(), :B=>HelmertCoding(), :C=>HelmertCoding())

            mf = ModelFrame(@formula(Y ~ 1 + A*B*C), tbl)
            mf_helm = ModelFrame(@formula(Y ~ 1 + A*B*C), tbl, contrasts = contrasts)

            @test size(modelmatrix(mf)) == size(modelmatrix(mf_helm))

            mf_helm2 = setcontrasts!(ModelFrame(@formula(Y ~ 1 + A*B*C), tbl), contrasts)

            @test size(modelmatrix(mf)) == size(modelmatrix(mf_helm2))
            @test modelmatrix(mf_helm) == modelmatrix(mf_helm2)

        end
    end

    @testset "#136" begin
        t = (x = rand(100), y = randn(100));
        f = @formula(y ~ x)
        @test_throws ArgumentError modelcols(f, t)
    end

    @testset "#185 - interactions of scalar terms for row tables" begin
        t = (a = rand(10), b = rand(10), c = rand(10))
        f = apply_schema(@formula(0 ~ a&b&c), schema(t))
        @test vec(modelcols(f.rhs, t)) == modelcols.(Ref(f.rhs), Tables.rowtable(t))
    end

end


text:
# Estimation/Simulation options, part of each model
@with_kw mutable struct estimOptions
	algo = NelderMead()
	innerAlgo = NLopt.LN_BOBYQA # used only for MultiStart
	lb = nothing
	ub = nothing
	autodiff = GalacticOptim.AutoFiniteDiff(;fdtype = Val(:central),fdhtype = Val(:hcentral))
	nDraws::Int64 = 500
	useHalton::Bool = false
	showEvery::Int64 = 5
	maxIter::Int64 = 100_000
	f_tol::Float64 = 0.0
	x_tol::Float64 = 0.0
	g_tol::Float64 = 1e-8
end

# Define abstract type, where each model is subtype of this
abstract type Model end

#################################################
# Model definitions
# Uses Parameter.jl pkg which allows to set default values
@with_kw mutable struct Weitzman <: Model
	β::Array{Float64,1}
	cs::Float64
	cd::Float64
	nd::Int64
	nA0::Int64
	dChars::Distribution
	outDummy::Bool
	dV::Normal			= Normal()
	dE::Normal			= Normal()
	cfun::String 		= "linear"
	seed::Int64 = 23909
	smo::Float64
	options::estimOptions = estimOptions()
end

@with_kw mutable struct FullInfo <: Model
	β::Array{Float64,1}
	dChars::Distribution
	outDummy::Bool
	dV::Normal
	dE::Normal
	seed::Int64 = 23909
	smo::Float64
	options::estimOptions = estimOptions()
end

@with_kw mutable struct RandomSearch <: Model
	β::Array{Float64,1}
	cs::Float64
	dChars::Distribution
	outDummy::Bool
	dV::Normal
	dE::Normal
	cfun::String
	seed::Int64 = 23909
	smo::Float64
	options::estimOptions = estimOptions()
end

@with_kw mutable struct SearchDiscovery <: Model
	β::Array{Float64,1}
	cs::Float64
	cd::Float64
	dze::Float64
	nd::Int64
	nA0::Int64 = 0
	dChars::Distribution
	outDummy::Bool
	dV::Normal
	dE::Normal
	cfun::String
	zdfun::String
	seed::Int64 = 23909
	options::estimOptions = estimOptions()
end


text:
module CreateTableChirps

import Migration: create_table, column, column_id, add_index, drop_table

function up()
  create_table(:chirps) do
    [
      column_id()
      column(:content, :text)
      column(:created_at, :datetime)
    ]
  end

  add_index(:chirps, :created_at)
end

function down()
  drop_table(:chirps)
end

end

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
using FactCheck
using ValidatedNumerics

facts("Operations on boxes") do
    A = IntervalBox(1..2, 3..4)
    B = IntervalBox(0..2, 3..6)

    @fact 2*A --> IntervalBox(2..4, 6..8)
    @fact A + B --> IntervalBox(1..4, 6..10)
    @fact dot(A, B) --> @interval(9, 28)

    @fact A ⊆ B --> true

    X = IntervalBox(1..2, 3..4)
    Y = IntervalBox(3..4, 3..4)

    @fact isempty(X ∩ Y) --> true

    v = [@interval(i, i+1) for i in 1:10]
    V = IntervalBox(v...)
    @fact length(V) --> 10

    Y = IntervalBox(1..2)  # single interval
    @fact isa(Y, IntervalBox) --> true
    @fact length(Y) --> 1
    @fact Y --> IntervalBox( (Interval(1., 2.),) )

end

facts("@intervalbox tests") do
    @intervalbox f(x, y) = (x + y, x - y)

    X = IntervalBox(1..1, 2..2)
    @fact f(X) --> IntervalBox(3..3, -1 .. -1)

    @intervalbox g(x, y) = x - y
    @fact isa(g(X), IntervalBox) --> true

    @fact emptyinterval(X) --> IntervalBox(∅, ∅)

end

facts("setdiff for IntervalBox") do
    X = IntervalBox(2..4, 3..5)
    Y = IntervalBox(3..5, 4..6)
    @fact setdiff(X, Y) --> [ IntervalBox(3..4, 3..4),
                              IntervalBox(2..3, 3..5) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, 4..5)
    @fact setdiff(X, Y) --> [ IntervalBox(2..5, 3..4),
                              IntervalBox(2..5, 5..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(4..6, 4..5)
    @fact setdiff(X, Y) --> [ IntervalBox(4..5, 3..4),
                              IntervalBox(4..5, 5..6),
                              IntervalBox(2..4, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(3..4, 4..5)
    @fact setdiff(X, Y) --> [ IntervalBox(3..4, 3..4),
                              IntervalBox(3..4, 5..6),
                              IntervalBox(2..3, 3..6),
                              IntervalBox(4..5, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(2..4, 10..20)
    @fact setdiff(X, Y) --> typeof(X)[X]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, -10..10)
    @fact setdiff(X, Y) --> typeof(X)[]

    X = IntervalBox(1..4, 3..6, 7..10)
    Y = IntervalBox(2..3, 4..5, 8..9)
    @fact setdiff(X, Y) --> [ IntervalBox(2..3, 4..5, 7..8),
                              IntervalBox(2..3, 4..5, 9..10),
                              IntervalBox(2..3, 3..4, 7..10),
                              IntervalBox(2..3, 5..6, 7..10),
                              IntervalBox(1..2, 3..6, 7..10),
                              IntervalBox(3..4, 3..6, 7..10) ]

end