text:
###################################################
# TESTES COM ARQUIVOS DA CUTEst
# Este script lê o diretório indicado, executa
# métodos para cada arquivo .SIF encontrado e
# tabela saídas em "resultados.txt".
#
# Autor: Leonardo D. Secchin
# Data : mai 2021
#
# Uso:
#   testesSIF(sifpath="caminho diretório SIFs")
# ou
#   testesSIF()   (lê diretório padrão "sif")
###################################################

using CUTEst, Printf

# carrega método do gradiente com busca linear por Armijo + interpolação quadrática + salvaguardas
include("gradiente_interp.jl")

# carrega método do gradiente espectral projetado (SPG)
# include("spg.jl")  #***** IMPLEMENTE SPG E INCLUA O ARQUIVO AQUI  *****

function testesSIF(; sifpath=pwd()*"/sif")

    sifdir = readdir(sifpath, join=true)

    # máximo de iterações para os métodos
    maxit = 10000

    # abre arquivo de saída em modo inserção
    arq = open("resultados.txt", "a")

    try

    # escreve cabeçalho no arquivo de saída
    write(arq, "\n\n        Prob |      n | método |    it |         f |     grad | st |\n"*
                   "====================================================================\n")
    flush(arq)

    compilar = true

    for i in 1:length(sifdir)

        # Verifica se o arquivo sifdir[i] tem extensão .SIF
        if uppercase(sifdir[i][max(end-3,1):end]) == ".SIF"

            # Na primeiro SIF encontrado, executa métodos com 1 iteração
            # somente para Julia compilar os códigos. Nenhuma saida é tabelada.
            # Assim, o tempo do primeiro problema não é afetado por tempos de
            # compilação.
            if compilar
                nlp = CUTEstModel(sifdir[i]);
                gradiente_interp(nlp, maxiter=1, saidas=false);
#                 spg(nlp, maxiter=1, saidas=false); #***** IMPLEMENTE SPG E DESCOMENTE ESTA LINHA  *****
                compilar = false
                finalize(nlp)
            end

            probname = basename(sifdir[i])[1:end-4]

            println("==============\nCarregando o problema "*probname*"...")

            # carrega arquivo SIF para estrutura NLPModels
            nlp = CUTEstModel(sifdir[i]);

            # pula problema caso o sentido de otimização não seja "minimização"
            if nlp.meta.minimize == true

                # escreve dados do problema no arquivo de saída
                write(arq, "-------------|--------|--------|-------|-----------|----------|----|\n")
                write(arq, @sprintf("%12s | %6d", probname, nlp.meta.nvar))
                flush(arq)

                # variável para impressão de resultados com vários métodos
                primeiro_res = true

                # conta variáveis que possuem limitante inferior ou superior
                boundedvars = count((nlp.meta.lvar .> -Inf) .| (nlp.meta.uvar .< Inf))

                ###############################################
                # APLICA MÉTODO DO GRADIENTE
                # se o problema for irrestrito.
                ###############################################
                if (nlp.meta.ncon == 0) & (boundedvars == 0)

                    println("Resolvendo "*probname*" pelo método do gradiente...")

                    # aplica gradiente a partir do ponto inicial fornecido no problema
                    ~, gra_f, gra_g, gra_it, gra_st = gradiente_interp(nlp, x0=nlp.meta.x0, maxiter=maxit, saidas=false);

                    # escreve resultado no arquivo de saida
                    write(arq, @sprintf(" |   grad | %5d | %9.2e | %8.2e | %2d |\n",
                            gra_it, gra_f, gra_g, gra_st))
                    flush(arq)

                    primeiro_res = false

                end

                ###############################################
                # APLICA MÉTODO DO GRAD. ESPECTRAL PROJ. (SPG)
                # se o problema for irrestrito ou só tiver
                # restrições de caixa.
                ###############################################
                if nlp.meta.ncon == 0

                    #**********************************************************
                    #***** IMPLEMENTE SPG E DESCOMENTE AS LINHAS A SEGUIR *****
                    #**********************************************************

#                     println("Resolvendo "*probname*" pelo método do gradiente espectral projetado...")
#
#                     # aplica SPG a partir do ponto inicial fornecido no problema
#                     ~, spg_f, spg_g, spg_it, spg_st = spg(nlp, x0=nlp.meta.x0, maxiter=maxit, saidas=false);
#
#                     # escreve resultado no arquivo de saida
#                     # se já há linha de método anterior, tabula adequadamente
#                     if !primeiro_res
#                         write(arq, "             |       ")
#                     end
#
#                     write(arq, @sprintf(" |    spg | %5d | %9.2e | %8.2e | %2d |\n",
#                             spg_it, spg_f, spg_g, spg_st))
#                     flush(arq)
#
#                     primeiro_res = false

                end

                # OBS:
                # Você pode inserir outros métodos a serem executados
                # em sequência. Adapte o bloco acima.

            end

            # fecha SIF atual
            finalize(nlp);

        end
    end

    close(arq);

    # tratamento de erros
    catch err
        close(arq);

        if err isa InterruptException
            println("\nCTRL+C pressionado. Tentando finalizar nlp...");
            finalize(nlp);
        else
            println("============================================\n"*
                    "  Um erro ocorreu :(\n"*
                    "  Isso pode ter sido causado por\n"*
                    "  1. um CTRL+C acionado no meio do processo.\n"*
                    "     Assim, o problema SIF não foi finalizado, e a estrutura nlp foi perdida.\n"*
                    "     Re-inclua testesSIF.jl e tente novamente. Caso não dê certo, finalize o\n"*
                    "     Julia e recarregue os scripts;\n"*
                    "  2. uso de espaços no caminho do diretório sif. Aparentemente o pacote\n"*
                    "     CUTEst.jl não lida bem com espaços. Certifique-se de que diretórios\n"*
                    "     e arquivos não possuam espaços.\n"*
                    "============================================");
            println("Erro:\n");
            println(err);
            println("\nObs: você pode descomentar o bloco try-catch para visualizar a saída completa do erro.")
        end
    end
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
immutable CacheStack{T}
  size::Int
  stack::Vector{T}
end

typealias CacheStackT{T} Type{CacheStack{T}}

(::CacheStackT{T}){T}(size = 1) = CacheStack{T}(size, T[])

function touch!{T}(c::CacheStack{T}, x::T)
  isloaded(x) && return x
  if isempty(c.stack)
    push!(c.stack, x)
  elseif c.stack[1] === x
    return x
  else
    cache = store!(c.stack[1])
    c.stack[1] = x
    return load!(x, cache)
  end
  return load!(x)
end

@enum CacheState Stored Loaded Modified

type CacheVector{T,A} <: AVector{T}
  cache::CacheStack{CacheVector{T,A}}
  data::A
  view::Vector{T}
  state::CacheState
end

CacheVector(c::CacheStack, xs::AVector) =
  CacheVector(c, xs, eltype(xs)[], Stored)

CacheVector{T,A}(c::CacheStack{CacheVector{T,A}}) = CacheVector(c, A())

isloaded(v::CacheVector) = v.state ≠ Stored

function load!(v::CacheVector)
  isloaded(v) && return v
  # info("loading data")
  v.view = collect(v.data)
  v.state = Loaded
  return v
end

function load!(v::CacheVector, cache)
  if length(cache) == length(v.data)
    copy!(cache, v.data)
    v.view = cache
    v.state = Loaded
  else
    return load!(v)
  end
  return load!(v)
end

function store!(v::CacheVector)
  isloaded(v) || return v
  cache = v.view
  if v.state == Modified
    # info("storing data")
    v.data[1:end] = slice(v.view, 1:endof(v.data))
    if length(cache) > length(v.data)
      append!(v.data, slice(cache, endof(v.data)+1:endof(cache)))
    end
    # TODO: catch shortened arrays as well
  end
  v.view = []
  v.state = Stored
  return cache
end

touch!(v::CacheVector) = touch!(v.cache, v)

for f in :[Base.getindex, Base.size].args
  @eval @inline function $f(xs::CacheVector, args...)
    touch!(xs)
    $f(xs.view, args...)
  end
end

for f in :[Base.setindex!, Base.push!].args
  @eval @inline function $f(xs::CacheVector, args...)
    touch!(xs)
    xs.state = Modified
    $f(xs.view, args...)
  end
end

Base.similar(xs::CacheVector, args...) =
  CacheVector(xs.cache, similar(xs.data, args...))

ext: jl


text:
# Event (experimental functionality)

function wgetch_events()
  # extern NCURSES_EXPORT(int) wgetch_events (WINDOW *, _nc_eventlist *);
end

function wgetnstr_events()
  # extern NCURSES_EXPORT(int) wgetnstr_events (WINDOW *,char *,int,_nc_eventlist *)
end

# * Function prototypes.  This is the complete X/Open Curses list of required
# * functions.  Those marked `generated' will have sources generated from the
# * macro definitions later in this file, in order to satisfy XPG4.2
# * requirements.

###  addch

# extern NCURSES_EXPORT(int) addch (const chtype);			/* generated */
# extern NCURSES_EXPORT(int) addchnstr (const chtype *, int);		/* generated */
# extern NCURSES_EXPORT(int) addchstr (const chtype *);			/* generated */
# extern NCURSES_EXPORT(int) mvaddch (int, int, const chtype);		/* generated */
# extern NCURSES_EXPORT(int) mvaddchnstr (int, int, const chtype *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvaddchstr (int, int, const chtype *);	/* generated */
# extern NCURSES_EXPORT(int) mvwaddch (WINDOW *, int, int, const chtype);	/* generated */
function mvwaddch( win, y::Int, x::Int, c )
    ccall((:mvwaddch, ncurses), Void,
        ( Ptr{Void}, Int, Int, Int ), win, y, x, c )
end
# extern NCURSES_EXPORT(int) mvwaddchnstr (WINDOW *, int, int, const chtype *, int);/* generated
# extern NCURSES_EXPORT(int) mvwaddchstr (WINDOW *, int, int, const chtype *);	/* generated */
# extern NCURSES_EXPORT(int) waddch (WINDOW *, const chtype);		/* implemented */
# extern NCURSES_EXPORT(int) waddchnstr (WINDOW *,const chtype *,int);	/* implemented */
# extern NCURSES_EXPORT(int) waddchstr (WINDOW *,const chtype *);		/* generated */

###  attr

# extern NCURSES_EXPORT(int) attroff (NCURSES_ATTR_T);			/* generated */
# extern NCURSES_EXPORT(int) attron (NCURSES_ATTR_T);			/* generated */
# extern NCURSES_EXPORT(int) attrset (NCURSES_ATTR_T);			/* generated */
# extern NCURSES_EXPORT(int) attr_get (attr_t *, NCURSES_PAIRS_T *, void *);
# extern NCURSES_EXPORT(int) attr_off (attr_t, void *);			/* generated */
# extern NCURSES_EXPORT(int) attr_on (attr_t, void *);			/* generated */
# extern NCURSES_EXPORT(int) attr_set (attr_t, NCURSES_PAIRS_T, void *);
# extern NCURSES_EXPORT(int) wattron (WINDOW *, int);			/* generated */
function wattron( win, attrs )
    ccall((:wattron, ncurses), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end
# extern NCURSES_EXPORT(int) wattroff (WINDOW *, int);			/* generated */
function wattroff( win, attrs )
    ccall((:wattroff, ncurses), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end
# extern NCURSES_EXPORT(int) wattrset (WINDOW *, int);			/* generated */
function wattrset( win, attrs )
    ccall((:wattrset, ncurses), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end
# extern NCURSES_EXPORT(int) wattr_get (WINDOW *, attr_t *, NCURSES_PAIRS_T *, void *);
# extern NCURSES_EXPORT(int) wattr_on (WINDOW *, attr_t, void *);		/* implemented */
# extern NCURSES_EXPORT(int) wattr_off (WINDOW *, attr_t, void *);	/* implemented */
# extern NCURSES_EXPORT(int) wattr_set (WINDOW *, attr_t, NCURSES_PAIRS_T, void *);	/*

##### win

# extern NCURSES_EXPORT(WINDOW *) newwin (int,int,int,int);	       	/* implemented */
function newwin( lines::Int, cols::Int, origy::Int, origx::Int )
    ccall((:newwin, ncurses), Ptr{Void}, ( Int, Int, Int, Int ),
        lines, cols, origy, origx )
end
# extern NCURSES_EXPORT(bool) isendwin (void);				/* implemented */
function isendwin()
    ccall((:isendwin, ncurses), Bool, () )
end
# extern NCURSES_EXPORT(int) endwin (void);				/* implemented */
function endwin()
    ccall((:endwin, ncurses), Int, () )
end
# extern NCURSES_EXPORT(WINDOW *) subwin (WINDOW *, int, int, int, int);	/* implemented */
function subwin(win::Ptr{Void}, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall((:subwin, ncurse ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end
# extern NCURSES_EXPORT(WINDOW *) derwin (WINDOW *,int,int,int,int);	/* implemented */
function derwin( win, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall((:derwin, ncurses ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end
# extern NCURSES_EXPORT(int) delwin (WINDOW *);				/* implemented */
function delwin( win::Ptr{Void} )
    ccall((:delwin, ncurses ), Void, ( Ptr{Void}, ), win )
end
# extern NCURSES_EXPORT(int) copywin (const WINDOW*,WINDOW*,int,int,int,int,int,int,int);	/* implemented */
# extern NCURSES_EXPORT(WINDOW *) dupwin (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(WINDOW *) getwin (FILE *);			/* implemented */
# extern NCURSES_EXPORT(int) untouchwin (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(int) touchwin (WINDOW *);				/* generated */
function touchwin( win::Ptr{Void} )
    ccall((:touchwin, ncurses ), Void, ( Ptr{Void}, ), win )
end
# extern NCURSES_EXPORT(int) redrawwin (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(int) putwin (WINDOW *, FILE *);			/* implemented */
# extern NCURSES_EXPORT(int) mvderwin (WINDOW *, int, int);		/* implemented */
# extern NCURSES_EXPORT(int) mvwin (WINDOW *,int,int);			/* implemented */
function mvwin( win, y::Int, x::Int )
    ccall((:mvwin, ncurses), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end

#### str

# extern NCURSES_EXPORT(int) addnstr (const char *, int);			/* generated */
# extern NCURSES_EXPORT(int) addstr (const char *);			/* generated */
# extern NCURSES_EXPORT(int) mvwaddnstr (WINDOW *, int, int, const char *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvwaddstr (WINDOW *, int, int, const char *);	/* generated */
# extern NCURSES_EXPORT(int) mvwgetnstr (WINDOW *, int, int, char *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvwgetstr (WINDOW *, int, int, char *);	/* generated */
# extern NCURSES_EXPORT(int) mvwin chnstr (WINDOW *, int, int, chtype *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvwinchstr (WINDOW *, int, int, chtype *);		/* generated */
# extern NCURSES_EXPORT(int) mvwinnstr (WINDOW *, int, int, char *, int);		/* generated */
# extern NCURSES_EXPORT(int) mvwinsnstr (WINDOW *, int, int, const char *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvwinsstr (WINDOW *, int, int, const char *);	/* generated */
# extern NCURSES_EXPORT(int) mvwinstr (WINDOW *, int, int, char *);		/* generated */
# extern NCURSES_EXPORT(int) mvinsnstr (int, int, const char *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvinsstr (int, int, const char *);		/* generated */
# extern NCURSES_EXPORT(int) mvinstr (int, int, char *);			/* generated */
# extern NCURSES_EXPORT(int) mvinchnstr (int, int, chtype *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvinchstr (int, int, chtype *);		/* generated */
# extern NCURSES_EXPORT(int) mvinnstr (int, int, char *, int);		/* generated */
# extern NCURSES_EXPORT(int) mvgetnstr (int, int, char *, int);		/* generated */
# extern NCURSES_EXPORT(int) mvgetstr (int, int, char *);			/* generated */
# extern NCURSES_EXPORT(int) mvaddnstr (int, int, const char *, int);	/* generated */
# extern NCURSES_EXPORT(int) mvaddstr (int, int, const char *);		/* generated */
# extern NCURSES_EXPORT(int) innstr (char *, int);			/* generated */
# extern NCURSES_EXPORT(int) inchnstr (chtype *, int);			/* generated */
# extern NCURSES_EXPORT(int) inchstr (chtype *);				/* generated */
# extern NCURSES_EXPORT(int) getnstr (char *, int);			/* generated */
# extern NCURSES_EXPORT(int) getstr (char *);				/* generated */
# extern NCURSES_EXPORT(int) insnstr (const char *, int);			/* generated */
# extern NCURSES_EXPORT(int) insstr (const char *);			/* generated */
# extern NCURSES_EXPORT(int) instr (char *);				/* generated */
# extern NCURSES_EXPORT(int) wgetnstr (WINDOW *,char *,int);		/* implemented */
# extern NCURSES_EXPORT(int) wgetstr (WINDOW *, char *);			/* generated */
# extern NCURSES_EXPORT(int) winchnstr (WINDOW *, chtype *, int);		/* implemented */
# extern NCURSES_EXPORT(int) winchstr (WINDOW *, chtype *);		/* generated */
# extern NCURSES_EXPORT(int) winnstr (WINDOW *, char *, int);		/* implemented */
# extern NCURSES_EXPORT(int) winsnstr (WINDOW *, const char *,int);	/* implemented */
# extern NCURSES_EXPORT(int) winsstr (WINDOW *, const char *);		/* generated */
# extern NCURSES_EXPORT(int) winstr (WINDOW *, char *);			/* generated */
# extern NCURSES_EXPORT(int) waddnstr (WINDOW *,const char *,int);	/* implemented */
# extern NCURSES_EXPORT(int) waddstr (WINDOW *,const char *);		/* generated */

### print/scan

# extern NCURSES_EXPORT(int) wprintw (WINDOW *, const char *,...)		/* implemented */
function wprintw( win::Ptr{Void} , str::String )
    ccall((:printw, ncurses), Void, (Ptr{Void}, Ptr{Uint8}), win, str )
end
# extern NCURSES_EXPORT(int) vwprintw (WINDOW *, const char *,va_list);	/* implemented */
# extern NCURSES_EXPORT(int) vw_printw (WINDOW *, const char *,va_list);	/* generated */
# extern NCURSES_EXPORT(int) printw (const char *,...)			/* implemented */å
function printw( str::String )
    ccall((:printw, ncurses, ncurses), Void, (Ptr{Uint8},), str )
end
# extern NCURSES_EXPORT(int) mvwprintw (WINDOW*,int,int, const char *,...)	/* implemented */
function mvwprintw( win::Ptr{Void}, row::Int, height::Int, fmt::String, str::String )
    ccall((:mvwprintw, ncurses), Void,
        ( Ptr{Void}, Int, Int, Ptr{Uint8}, Ptr{Uint8}),
        win, row, height, fmt, str )
end
# extern NCURSES_EXPORT(int) mvprintw (int,int, const char *,...)		/* implemented */
function mvprintw( row::Int, height::Int, str::String )
    ccall((:printw, ncurses), Void, (Int, Int, Ptr{Uint8}), row, height, str )
end
# extern NCURSES_EXPORT(int) wscanw (WINDOW *, NCURSES_CONST char *,...)	/* implemented */
# extern NCURSES_EXPORT(int) vwscanw (WINDOW *, NCURSES_CONST char *,va_list);	/* implemented */
# extern NCURSES_EXPORT(int) vw_scanw (WINDOW *, NCURSES_CONST char *,va_list);	/* generated */
# extern NCURSES_EXPORT(int) scanw (NCURSES_CONST char *,...)		/* implemented */
# extern NCURSES_EXPORT(int) mvwscanw (WINDOW *,int,int, NCURSES_CONST char *,...)
# extern NCURSES_EXPORT(int) mvscanw (int,int, NCURSES_CONST char *,...)	/* implemented */

#### Varia

# extern NCURSES_EXPORT(WINDOW *) initscr (void);				/* implemented */
function initscr()
    ccall((:initscr, ncurses ), Ptr{Void}, () )
end
# extern NCURSES_EXPORT(int) baudrate (void);				/* implemented */
function baudrate()
    ccall((:baudrate, ncurses), Int, () )
end
# extern NCURSES_EXPORT(int) beep  (void);				/* implemented */
function beep()
    ccall((:beep, ncurses ), Void, () )
end
# extern NCURSES_EXPORT(int) bkgd (chtype);				/* generated */
# extern NCURSES_EXPORT(void) bkgdset (chtype);				/* generated */
# extern NCURSES_EXPORT(int) border (chtype,chtype,chtype,chtype,chtype,chtype,chtype,chtype);	/* generated */
# extern NCURSES_EXPORT(int) box (WINDOW *, chtype, chtype);		/* generated */
function box( win, vchr, hchr )
    ccall((:box, ncurses ), Void, (Ptr{Void}, Char, Char), win, vchr, hchr )
end
# extern NCURSES_EXPORT(bool) can_change_color (void);			/* implemented */
# extern NCURSES_EXPORT(int) cbreak (void);				/* implemented */
function cbreak()
    ccall((:cbreak, ncurses), Void, ( ) )
end
# extern NCURSES_EXPORT(int) chgat (int, attr_t, NCURSES_PAIRS_T, const void *);	/* generated */
# extern NCURSES_EXPORT(int) clear (void);				/* generated */
# extern NCURSES_EXPORT(int) clearok (WINDOW *,bool);			/* implemented */
function clearok( win, bf )
    ccall((:clearok, ncurses), Int, ( Ptr{Void}, Bool ), win, bf )
end
# extern NCURSES_EXPORT(int) clrtobot (void);				/* generated */
# extern NCURSES_EXPORT(int) clrtoeol (void);				/* generated */
# extern NCURSES_EXPORT(int) color_content (NCURSES_COLOR_T,NCURSES_COLOR_T*,NCURSES_COLOR_T*,NCURSES_COLOR_T*);	/* implemented */
# extern NCURSES_EXPORT(int) color_set (NCURSES_PAIRS_T,void*);			/* generated */
# extern NCURSES_EXPORT(int) COLOR_PAIR (int);				/* generated */
# extern NCURSES_EXPORT(int) curs_set (int);				/* implemented */
function curs_set( vis )
    ccall((:curs_set, ncurses), Int, ( Int, ), vis )
end
# extern NCURSES_EXPORT(int) def_prog_mode (void);			/* implemented */
# extern NCURSES_EXPORT(int) def_shell_mode (void);			/* implemented */
# extern NCURSES_EXPORT(int) delay_output (int);				/* implemented */
# extern NCURSES_EXPORT(int) delch (void);				/* generated */
# extern NCURSES_EXPORT(void) delscreen (SCREEN *);			/* implemented */
# extern NCURSES_EXPORT(int) deleteln (void);				/* generated */
# extern NCURSES_EXPORT(int) doupdate (void);				/* implemented *
# extern NCURSES_EXPORT(int) echo (void);					/* implemented */
function echo()
    ccall((:echo, ncurses), Void, ( ) )
end
# extern NCURSES_EXPORT(int) echochar (const chtype);			/* generated */
# extern NCURSES_EXPORT(int) erase (void);				/* generated */
function erase()
    ccall((:erase, ncurses ), Void, () )
end
# extern NCURSES_EXPORT(char) erasechar (void);				/* implemented */
# extern NCURSES_EXPORT(void) filter (void);				/* implemented */
# extern NCURSES_EXPORT(int) flash (void);				/* implemented */
function flash()
    ccall((:flash, ncurses), Void, () )
end
# extern NCURSES_EXPORT(int) flushinp (void);				/* implemented */
# extern NCURSES_EXPORT(chtype) getbkgd (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(int) getch (void);				/* generated */
# extern NCURSES_EXPORT(int) halfdelay (int);				/* implemented */
# extern NCURSES_EXPORT(bool) has_colors (void);				/* implemented */
function has_colors()
    ccall((:has_colors, ncurses), Bool, () )
end
# extern NCURSES_EXPORT(bool) has_ic (void);				/* implemented */
# extern NCURSES_EXPORT(bool) has_il (void);				/* implemented */
# extern NCURSES_EXPORT(int) hline (chtype, int);				/* generated */
# extern NCURSES_EXPORT(void) idcok (WINDOW *, bool);			/* implemented */
# extern NCURSES_EXPORT(int) idlok (WINDOW *, bool);			/* implemented */
# extern NCURSES_EXPORT(void) immedok (WINDOW *, bool);			/* implemented */
function immedok( win, bf )
    ccall((:immedok, ncurses), Int, ( Ptr{Void}, Bool ), win, bf )
end
# extern NCURSES_EXPORT(chtype) inch (void);				/* generated */
# extern NCURSES_EXPORT(int) init_color (NCURSES_COLOR_T,NCURSES_COLOR_T,NCURSES_COLOR_T,NCURSES_COLOR_T);	/* implemented */
function init_color( color, r,g,b )
    ccall((:init_color, ncurses), Int, ( Int16, Int16, Int16, Int16 ),
        color, r,g,b )
end
# extern NCURSES_EXPORT(int) init_pair (NCURSES_PAIRS_T,NCURSES_COLOR_T,NCURSES_COLOR_T);
function init_pair( pair, f, b )
    ccall((:init_pair, ncurses), Int, ( Int16, Int16, Int16 ), pair, f, b )
end
# extern NCURSES_EXPORT(int) insch (chtype);				/* generated */
# extern NCURSES_EXPORT(int) insdelln (int);				/* generated */
# extern NCURSES_EXPORT(int) insertln (void);				/* generated */
# extern NCURSES_EXPORT(int) intrflush (WINDOW *,bool);			/* implemented */
# extern NCURSES_EXPORT(bool) is_linetouched (WINDOW *,int);		/* implemented */
# extern NCURSES_EXPORT(bool) is_wintouched (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(NCURSES_CONST char *) keyname (int);		/* implemented */
# extern NCURSES_EXPORT(int) keypad (WINDOW *,bool);			/* implemented */
function keypad( win, bf )
    ccall((:keypad, ncurses ), Int, (Ptr{Void}, Bool), win, bf )
end
# extern NCURSES_EXPORT(char) killchar (void);				/* implemented */
# extern NCURSES_EXPORT(int) leaveok (WINDOW *,bool);			/* implemented */
# extern NCURSES_EXPORT(char *) longname (void);				/* implemented */
# extern NCURSES_EXPORT(int) meta (WINDOW *,bool);			/* implemented */
# extern NCURSES_EXPORT(int) move (int, int);				/* generated */
# extern NCURSES_EXPORT(int) mvchgat (int, int, int, attr_t, NCURSES_PAIRS_T, const void *);	/* generated */
# extern NCURSES_EXPORT(int) mvcur (int,int,int,int);			/* implemented */
# extern NCURSES_EXPORT(int) mvdelch (int, int);				/* generated */
# extern NCURSES_EXPORT(int) mvgetch (int, int);				/* generated */
# extern NCURSES_EXPORT(int) mvhline (int, int, chtype, int);		/* generated */
# extern NCURSES_EXPORT(chtype) mvinch (int, int);			/* generated */
# extern NCURSES_EXPORT(int) mvinsch (int, int, chtype);			/* generated */
# extern NCURSES_EXPORT(int) mvvline (int, int, chtype, int);		/* generated */
# extern NCURSES_EXPORT(int) mvwchgat (WINDOW *, int, int, int, attr_t, NCURSES_PAIRS_T, const void *);/* generated */
# extern NCURSES_EXPORT(int) mvwdelch (WINDOW *, int, int);		/* generated */
# extern NCURSES_EXPORT(int) mvwgetch (WINDOW *, int, int);		/* generated */
# extern NCURSES_EXPORT(int) mvwhline (WINDOW *, int, int, chtype, int);	/* generated */
# extern NCURSES_EXPORT(chtype) mvwinch (WINDOW *, int, int);			/* generated */
# extern NCURSES_EXPORT(int) mvwinsch (WINDOW *, int, int, chtype);		/* generated */
# extern NCURSES_EXPORT(int) mvwvline (WINDOW *,int, int, chtype, int);	/* generated */
# extern NCURSES_EXPORT(int) napms (int);					/* implemented */
function napms( ms::Int )
    ccall((:napms, ncurses), Int, (Int, ), ms )
end
# extern NCURSES_EXPORT(WINDOW *) newpad (int,int);		       	/* implemented */
# extern NCURSES_EXPORT(SCREEN *) newterm (NCURSES_CONST char *,FILE *,FILE *);	/* implemented */
# extern NCURSES_EXPORT(int) nl (void);					/* implemented */
# extern NCURSES_EXPORT(int) nocbreak (void);				/* implemented */
function nocbreak()
    ccall((:nocbreak, ncurses), Void, ( ) )
end
# extern NCURSES_EXPORT(int) nodelay (WINDOW *,bool);			/* implemented */
function nodelay( win, bf )
    ccall((:nodelay, ncurses ), Int, (Ptr{Void}, Bool), win, bf )
end
# extern NCURSES_EXPORT(int) noecho (void);				/* implemented */
function noecho()
    ccall((:noecho, ncurses), Void, ( ) )
end
# extern NCURSES_EXPORT(int) nonl (void);					/* implemented */
# extern NCURSES_EXPORT(void) noqiflush (void);				/* implemented */
# extern NCURSES_EXPORT(int) noraw (void);				/* implemented */
function noraw()
    ccall((:noraw, ncurses ), Int, () )
end
# extern NCURSES_EXPORT(int) notimeout (WINDOW *,bool);			/* implemented */
function notimeout( win, bf )
    ccall((:notimeout, ncurses ), Int, (Ptr{Void}, Bool), win, bf )
end
# extern NCURSES_EXPORT(int) overlay (const WINDOW*,WINDOW *);		/* implemented */
# extern NCURSES_EXPORT(int) overwrite (const WINDOW*,WINDOW *);		/* implemented */
# extern NCURSES_EXPORT(int) pair_content (NCURSES_PAIRS_T,NCURSES_COLOR_T*,NCURSES_COLOR_T*);		/* implemented */
# extern NCURSES_EXPORT(int) PAIR_NUMBER (int);				/* generated */
# extern NCURSES_EXPORT(int) pechochar (WINDOW *, const chtype);		/* implemented */
# extern NCURSES_EXPORT(int) pnoutrefresh (WINDOW*,int,int,int,int,int,int);/* implemented */
# extern NCURSES_EXPORT(int) prefresh (WINDOW *,int,int,int,int,int,int);	/* implemented */
# extern NCURSES_EXPORT(void) qiflush (void);				/* implemented */
# extern NCURSES_EXPORT(int) raw (void);					/* implemented */
function raw()
    ccall((:raw, ncurses ), Int, () )
end
# extern NCURSES_EXPORT(int) refresh (void);				/* generated */
function refresh()
    ccall((:refresh, ncurses ), Void, ( ) )
end
# extern NCURSES_EXPORT(int) resetty (void);				/* implemented */
# extern NCURSES_EXPORT(int) reset_prog_mode (void);			/* implemented */
# extern NCURSES_EXPORT(int) reset_shell_mode (void);			/* implemented */
# extern NCURSES_EXPORT(int) ripoffline (int, int (*)(WINDOW *, int));	/* implemented */
# extern NCURSES_EXPORT(int) savetty (void);				/* implemented */
# extern NCURSES_EXPORT(int) scr_dump (const char *);			/* implemented */
# extern NCURSES_EXPORT(int) scr_init (const char *);			/* implemented */
# extern NCURSES_EXPORT(int) scrl (int);					/* generated */
# extern NCURSES_EXPORT(int) scroll (WINDOW *);				/* generated */
# extern NCURSES_EXPORT(int) scrollok (WINDOW *,bool);			/* implemented */
# extern NCURSES_EXPORT(int) scr_restore (const char *);			/* implemented */
# extern NCURSES_EXPORT(int) scr_set (const char *);			/* implemented */
# extern NCURSES_EXPORT(int) setscrreg (int,int);				/* generated */
# extern NCURSES_EXPORT(SCREEN *) set_term (SCREEN *);			/* implemented */
# extern NCURSES_EXPORT(int) slk_attroff (const chtype);			/* implemented */
# extern NCURSES_EXPORT(int) slk_attr_off (const attr_t, void *);		/* generated:WIDEC */
# extern NCURSES_EXPORT(int) slk_attron (const chtype);			/* implemented */
# extern NCURSES_EXPORT(int) slk_attr_on (attr_t,void*);			/* generated:WIDEC */
# extern NCURSES_EXPORT(int) slk_attrset (const chtype);			/* implemented */
# extern NCURSES_EXPORT(attr_t) slk_attr (void);				/* implemented */
# extern NCURSES_EXPORT(int) slk_attr_set (const attr_t,NCURSES_PAIRS_T,void*);	/* implemented */
# extern NCURSES_EXPORT(int) slk_clear (void);				/* implemented */
# extern NCURSES_EXPORT(int) slk_color (NCURSES_PAIRS_T);				/* implemented */
# extern NCURSES_EXPORT(int) slk_init (int);				/* implemented */
# extern NCURSES_EXPORT(char *) slk_label (int);				/* implemented */
# extern NCURSES_EXPORT(int) slk_noutrefresh (void);			/* implemented */
# extern NCURSES_EXPORT(int) slk_refresh (void);				/* implemented */
# extern NCURSES_EXPORT(int) slk_restore (void);				/* implemented */
# extern NCURSES_EXPORT(int) slk_set (int,const char *,int);		/* implemented */
# extern NCURSES_EXPORT(int) slk_touch (void);	      	       		/* implemented */
# extern NCURSES_EXPORT(int) standout (void);				/* generated */
# extern NCURSES_EXPORT(int) standend (void);				/* generated */
# extern NCURSES_EXPORT(int) start_color (void);				/* implemented */
function start_color()
    ccall((:start_color, ncurses), Void, () )
end
# extern NCURSES_EXPORT(WINDOW *) subpad (WINDOW *, int, int, int, int);	/* implemented */
# extern NCURSES_EXPORT(int) syncok (WINDOW *, bool);			/* implemented */
# extern NCURSES_EXPORT(chtype) termattrs (void);				/* implemented */
# extern NCURSES_EXPORT(char *) termname (void);				/* implemented */
# extern NCURSES_EXPORT(void) timeout (int);				/* generated */
function timeout( delay::Int )
    ccall((:timeout, ncurses ), Void, (Int,), delay )
end
# extern NCURSES_EXPORT(int) touchline (WINDOW *, int, int);		/* generated */
# extern NCURSES_EXPORT(int) typeahead (int);				/* implemented */
# extern NCURSES_EXPORT(int) ungetch (int);				/* implemented */
# extern NCURSES_EXPORT(void) use_env (bool);				/* implemented */
# extern NCURSES_EXPORT(void) use_tioctl (bool);				/* implemented */
# extern NCURSES_EXPORT(int) vidattr (chtype);				/* implemented */
# extern NCURSES_EXPORT(int) vidputs (chtype, NCURSES_OUTC);		/* implemented */
# extern NCURSES_EXPORT(int) vline (chtype, int);				/* generated */
# extern NCURSES_EXPORT(int) wbkgd (WINDOW *, chtype);			/* implemented */
function wbkgd( win, ch )
    ccall((:wbkgd, ncurses ), Void, ( Ptr{Void}, Uint32 ), win, ch )
end
# extern NCURSES_EXPORT(void) wbkgdset (WINDOW *,chtype);			/* implemented */
function wbkgdset( win, ch )
    ccall((:wbkgdset, ncurses ), Void, ( Ptr{Void}, Uint32 ), win, ch )
end
# extern NCURSES_EXPORT(int) wborder (WINDOW*,chtype,chtype,chtype,chtype,chtype,chtype,chtype,chtype);	/* implemented */
# extern NCURSES_EXPORT(int) wchgat (WINDOW *, int, attr_t, NCURSES_PAIRS_T, const void *);/* implemented */
# extern NCURSES_EXPORT(int) wclear (WINDOW *);				/* implemented */
function wclear( win::Ptr{Void} )
    ccall((:wclear, ncurses ), Void, (Ptr{Void},), win )
end
# extern NCURSES_EXPORT(int) wclrtobot (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(int) wclrtoeol (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(int) wcolor_set (WINDOW*,NCURSES_PAIRS_T,void*);		/* implemented */
# extern NCURSES_EXPORT(void) wcursyncup (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(int) wdelch (WINDOW *);				/* implemented */
# extern NCURSES_EXPORT(int) wdeleteln (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(int) wechochar (WINDOW *, const chtype);		/* implemented */
# extern NCURSES_EXPORT(int) werase (WINDOW *);				/* implemented */
function werase( win::Ptr{Void} )
    ccall((:werase, ncurses ), Void, (Ptr{Void},), win )
end
# extern NCURSES_EXPORT(int) wgetch (WINDOW *);				/* implemented */
function wgetch( win::Ptr{Void } )
    ccall((:wgetch, ncurses ), Uint32, (Ptr{Void},), win )
end
# extern NCURSES_EXPORT(int) whline (WINDOW *, chtype, int);		/* implemented */
# extern NCURSES_EXPORT(chtype) winch (WINDOW *);				/* implemented */
# extern NCURSES_EXPORT(int) winsch (WINDOW *, chtype);			/* implemented */
# extern NCURSES_EXPORT(int) winsdelln (WINDOW *,int);			/* implemented */
# extern NCURSES_EXPORT(int) winsertln (WINDOW *);			/* generated */
function wmove( win, y::Int, x::Int )
  # extern NCURSES_EXPORT(int) wmove (WINDOW *,int,int);			/* implemented */
  ccall((:wmove, ncurses), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end
# extern NCURSES_EXPORT(int) wnoutrefresh (WINDOW *);			/* implemented *
# extern NCURSES_EXPORT(int) wredrawln (WINDOW *,int,int);		/* implemented */
# extern NCURSES_EXPORT(int) wrefresh (WINDOW *);				/* implemented */
function wrefresh( win::Ptr{Void} )
    ccall((:wrefresh, ncurses ), Void, ( Ptr{Void}, ), win )
end
# extern NCURSES_EXPORT(int) wscrl (WINDOW *,int);			/* implemented */
# extern NCURSES_EXPORT(int) wsetscrreg (WINDOW *,int,int);		/* implemented */
# extern NCURSES_EXPORT(int) wstandout (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(int) wstandend (WINDOW *);			/* generated */
# extern NCURSES_EXPORT(void) wsyncdown (WINDOW *);			/* implemented */
# extern NCURSES_EXPORT(void) wsyncup (WINDOW *);				/* implemented */
# extern NCURSES_EXPORT(void) wtimeout (WINDOW *,int);			/* implemented */
function wtimeout( win, delay::Int )
    ccall((:wtimeout, ncurses ), Void, (Ptr{Void}, Int), win, delay )
end
# extern NCURSES_EXPORT(int) wtouchln (WINDOW *,int,int,int);		/* implemented */
# extern NCURSES_EXPORT(int) wvline (WINDOW *,chtype,int);		/* implemented */

####  not standard but convenient TermWin

function getwinmaxyx( win )
    maxy = ccall((:getmaxy, ncurses), Int, ( Ptr{Void}, ), win )
    maxx = ccall((:getmaxx, ncurses), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function getwinbegyx( win )
    maxy = ccall((:getbegy, ncurses), Int, ( Ptr{Void}, ), win )
    maxx = ccall((:getbegx, ncurses), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function is_term_resized( lines::Int, cols::Int )
    ccall((:is_term_resized, ncurses ), Bool, (Int, Int), lines, cols )
end

function wresize( win, lines::Int, cols::Int )
    ccall((:wresize), Int, (Ptr{Void}, Int, Int), win, lines, cols )
end

function getcuryx( win )
    cury = ccall((:getcury, ncurses), Cint, ( Ptr{Void}, ), win )
    curx = ccall((:getcurx, ncurses), Cint, ( Ptr{Void}, ), win )
    ( cury, curx )
end

function has_mouse()
    ccall((:has_mouse, ncurses), Bool, () )
end

function mousemask( mask )
    oldmm = Array( Uint64, 1 )
    resultmm = ccall((:mousemask, ncurses), Uint64, (Uint64, Ptr{Uint64}), mask, oldmm )
    ( resultmm, oldmm[1] )
end

function mouseinterval( n::Int )
    ccall((:mouseinterval, ncurses), Int, (Int, ) , n)
end

#hack!
const mouseByteString = bytestring( Array( Uint8, 64 ) )
function getmouse()
    #=
    type Mouse_Event_t
        id::Int8 # short
        x::Int32 # int
        y::Int32 # int
        z::Int32 # int
        bstate::Uint64 # unsigned long
    end
    =#
    # the 5th byte is x, 9th byte is y
    # 19th byte is x08 if scroll up, 20th byte is x08 if scroll down
    # 17th byte is x02 if button 1 pressed
    # 17th byte is x01 if button 1 released
    # 17-18th is 0xfffd if mousewheel is pressed down
    ccall((:getmouse, ncurses), Int, (Ptr{Uint8}, ), mouseByteString )
    bs = mouseByteString
    x = uint8(bs[5])
    y = uint8(bs[9])
    state=:unknown
    println(bs, x, y, state)
    # if bs[17] & 0x02 != 0
    #     state = :button1_pressed
    # elseif bs[19] & 0x08 != 0
    #     state = :scroll_up
    # elseif bs[20] & 0x08 != 0
    #     state = :scroll_down
    # end
    ( state, x, y, bs )
end

export initscr, endwin, isendwin, newwin, subwin, delwin, derwin, werase, erase, refresh, mvwaddch, wclear, mvwprintw, touchwin, wmove, box, wrefresh, wgetch, nocbreak, keypad, nodelay, noecho, cbreak, echo, raw, noraw, timeout, notimeout, mvwin, getwinbegyx, beep, wtimeout, flash, getwinmaxyx, is_term_resized, wresize, getcuryx, start_color, init_pair, init_color, has_colors, wattroff, wattron, wattrset, wbkgdset, wbkgd, curs_set, has_mouse, mousemask, mouseinterval, getmouse, mouseByteString, baudrate, clearok, immedok, napms


text:
#=
---------------------------------------------------------------------
This function plots the solution
Written by F.X. Giraldo on April 19 2019
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function plot_solution(q,coord_CG,periodicity_CG,DG_to_CG,space_method,time,Nx,Ny)

    #Allocate arrays
    Npoin_CG=Nx*Ny
    Npoin=length(q)
    q_sol=zeros(DFloat,Npoin_CG)
    lhowm=zeros(Int64,Npoin_CG)

    #Compute gridpoint solution
    for i=1:Npoin
        ip_CG=periodicity_CG[DG_to_CG[i]]
        q_sol[ip_CG] += q[i]
        lhowm[ip_CG] += 1
    end
    for i=1:Npoin_CG
        q_sol[i]=q_sol[i]/lhowm[i]
    end

    #Plot Solution
    x=reshape(coord_CG[1,:],(Nx,Ny))
    y=reshape(coord_CG[2,:],(Nx,Ny))
    qq=reshape(q_sol,(Nx,Ny))
    xx=x[:,1]
    yy=y[1,:]

    p1 = contourf(xx,yy,qq')
    timec=round(time/(2*π); digits=2)
    title!(p1,"$space_method solution at time = $timec")
    xlabel!(p1,"x")
    ylabel!(p1,"y")
    figure_1=plot(p1)
    display(figure_1)

end


text:
#Generic Branch Models
abstract type AbstractBranchFormulation <: AbstractDeviceFormulation end

#Abstract Line Models

abstract type AbstractLineFormulation <: AbstractBranchFormulation end

struct StaticLine <: AbstractLineFormulation end
struct StaticLineUnbounded <: AbstractLineFormulation end

struct FlowMonitoredLine <: AbstractLineFormulation end

#Abstract Transformer Models

abstract type AbstractTransformerFormulation <: AbstractBranchFormulation end

struct StaticTransformer <: AbstractTransformerFormulation end
struct StaticTransformerUnbounded <: AbstractTransformerFormulation end

# Not implemented yet
struct TapControl <: AbstractTransformerFormulation end
struct PhaseControl <: AbstractTransformerFormulation end

#################################### Branch Variables ##################################################
# Because of the way we integrate with PowerModels, most of the time PowerSimulations will create variables
# for the branch flows either in AC or DC.
function flow_variables!(canonical::Canonical,
                        system_formulation::Type{S},
                        devices::IS.FlattenIteratorWrapper{B}) where {B<:PSY.ACBranch,
                                                             S<:PM.AbstractPowerModel}
    return

end

function flow_variables!(canonical::Canonical,
                        system_formulation::Type{S},
                        devices::IS.FlattenIteratorWrapper{B}) where {B<:PSY.ACBranch,
                                                                      S<:StandardPTDFModel}

    var_name = Symbol("Fp_$(B)")

    add_variable(canonical,
                devices,
                var_name,
                false)
                #ub_value = d -> PSY.get_rate(d), # Add flow bounds in rate constraints
                #lb_value = d -> -1.0*PSY.get_rate(d))

    return

end

#################################### Flow Variable Bounds ##################################################
function branch_rate_bounds!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractDCPModel}

    range_data = [(PSY.get_name(h), (min = -1*PSY.get_rate(h), max = PSY.get_rate(h))) for h in devices]

    set_variable_bounds(canonical,
                        range_data,
                        Symbol("Fp_$(B)"))

    return

end

function branch_rate_bounds!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractActivePowerModel}

    range_data = [(PSY.get_name(h), (min = -1*PSY.get_rate(h), max = PSY.get_rate(h))) for h in devices]

    set_variable_bounds(canonical,
                        range_data,
                        Symbol("FpFT_$(B)"))

    set_variable_bounds(canonical,
                        range_data,
                        Symbol("FpTF_$(B)"))

    return

end

function branch_rate_bounds!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(h), (min = -1*PSY.get_rate(h), max = PSY.get_rate(h))) for h in devices]

    set_variable_bounds(canonical,
                        range_data,
                        Symbol("FpFT_$(B)"))

    set_variable_bounds(canonical,
                        range_data,
                        Symbol("FpTF_$(B)"))

    return

end

#################################### Rate Limits Constraints ##################################################
function branch_rate_constraint!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractDCPModel}

    range_data = [(PSY.get_name(h), (min = -1*PSY.get_rate(h), max = PSY.get_rate(h))) for h in devices]

    device_range(canonical,
                range_data,
                Symbol("RateLimit_$(B)"),
                Symbol("Fp_$(B)"))

    return

end

function branch_rate_constraint!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractActivePowerModel}

    range_data = [(PSY.get_name(h), (min = -1*PSY.get_rate(h), max = PSY.get_rate(h))) for h in devices]

    device_range(canonical,
                range_data,
                Symbol("RateLimitFT_$(B)"),
                Symbol("FpFT_$(B)"))

    device_range(canonical,
                range_data,
                Symbol("RateLimitTF_$(B)"),
                Symbol("FpTF_$(B)"))

    return

end

function branch_rate_constraint!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{B},
                                device_formulation::Type{D},
                                system_formulation::Type{S}) where {B<:PSY.ACBranch,
                                                                    D<:AbstractBranchFormulation,
                                                                    S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]

    rating_constraint!(canonical,
                        range_data,
                        Symbol("RateLimitFT_$(B)"),
                        (Symbol("FpFT_$(B)"), Symbol("FqFT_$(B)")))

    rating_constraint!(canonical,
                        range_data,
                        Symbol("RateLimitTF_$(B)"),
                        (Symbol("FpTF_$(B)"), Symbol("FqTF_$(B)")))

    return

end

#################################### Flow Limits Constraints ##################################################

function branch_flow_constraint!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{PSY.MonitoredLine},
                                device_formulation::Type{FlowMonitoredLine},
                                system_formulation::Union{Type{PM.DCPPowerModel}, Type{StandardPTDFModel}})

    flow_range_data = [(PSY.get_name(h), PSY.get_flowlimits(h)) for h in devices]

    var_name = Symbol("Fp_$(B)")

    device_range(canonical,
                range_data,
                Symbol("FlowLimit_$(B)"),
                var_name)

    return

end

function branch_flow_constraint!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{PSY.MonitoredLine},
                                device_formulation::Type{FlowMonitoredLine},
                                system_formulation::Type{S}) where {S<:PM.AbstractPowerModel}

    FTflow_range_data = [(PSY.get_name(h), PSY.get_flowlimits(h)) for h in devices]
    TFflow_range_data = [(PSY.get_name(h), (min = PSY.get_flowlimits(h).max, max = PSY.get_flowlimits(h).min)) for h in devices]

    device_range(canonical,
                FTflow_range_data,
                Symbol("FlowLimitFT_$(B)"),
                Symbol("FpFT_$(B)"))

    device_range(canonical,
                TFflow_range_data,
                Symbol("FlowLimitTF_$(B)"),
                Symbol("FpTF_$(B)"))

    return

end


text:
using StippleUI
using Test

@testset "StippleUI.jl" begin
    # Write your own tests here.
end


text:
#=
Define a discrete state space for the pedestrian MDP problem
=#

function n_states(mdp::PedMDP)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    n_ped = n_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    return n_ego*(n_ped + 1) # do not forget absent state
end

function POMDPs.states(mdp::PedMDP)
    state_space = PedMDPState[]
    ego_states = get_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    ped_states = get_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    for ego in ego_states
        for ped in ped_states
            crash =  is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(ped, mdp.ped_type, PED_ID))
            push!(state_space, PedMDPState(crash, ego, ped))
        end
        # add absent states
        push!(state_space, PedMDPState(false, ego, get_off_the_grid(mdp)))
    end
    return state_space
end

function POMDPs.stateindex(mdp::PedMDP, s::PedMDPState)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    n_ped = n_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    # step 1: find ego_index
    ego_i = ego_stateindex(mdp.env, s.ego, mdp.pos_res, mdp.vel_res)
    # step 2: find ped_index
    if s.ped.posG == mdp.off_grid
        ped_i = n_ped + 1
    else
        ped_i = ped_stateindex(mdp.env, s.ped, mdp.pos_res, mdp.vel_ped_res)
    end
    # step 3: find global index
    si = LinearIndices((n_ped+1, n_ego))[ped_i, ego_i]
end

function POMDPs.initialstate(mdp::PedMDP)
    ego = initial_ego_state(mdp)
    init_ped_dist = initial_ped_state_distribution(mdp)
    init_ped_states = init_ped_dist.vals
    init_states = Vector{PedMDPState}(length(init_ped_states))
    for i=1:length(init_ped_states)
        ped = init_ped_states[i]
        crash = is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(ped, mdp.ped_type, PED_ID))
        init_states[i] = PedMDPState(crash, ego, ped)
    end
    return SparseCat(init_states, init_ped_dist.probs)
end

function initial_ego_state(mdp::PedMDP)
    lanes = get_ego_route(mdp.env)
    posF = Frenet(mdp.env.roadway[lanes[1]], mdp.ego_start)
    v0 = 0.
    return VehicleState(posF, mdp.env.roadway, v0)
end

function initial_ped_state_distribution(mdp::PedMDP)
    init_ped_states = get_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    push!(init_ped_states, get_off_the_grid(mdp))
    # uniform (maybe add more weights to the states when pedestrians are not there?)
    probs = ones(length(init_ped_states))
    normalize!(probs, 1)
    return SparseCat(init_ped_states, probs)
end

function initial_ped_state(mdp::PedMDP, rng::AbstractRNG)
    init_dist = initial_ped_state_distribution(mdp)
    return rand(rng, init_dist)
end

function pedestrian_starting_states(mdp::PedMDP)
    # list of pedestrian starting states
    n_headings = 2
    lanes = get_ped_lanes(mdp.env)
    v_space = get_ped_vspace(mdp.env, mdp.vel_ped_res)
    ped_states = Vector{VehicleState}(length(lanes)*length(v_space)*n_headings)
    i = 1
    for lane in lanes
        for v in v_space
            ped_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], 0., 0., 0.), mdp.env.roadway, v)
            i += 1
            ped_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], get_end(mdp.env.roadway[lane]), 0., float(pi)), mdp.env.roadway, v)
            i += 1
        end
    end
    return ped_states
end


text:
import Base.exp
"""
computes log(x+y) given a = log(x) and b = log(y)
"""
log_add(a::Union{Float64, Array{Float64}},b::Union{Float64, Array{Float64}}) = max(a,b) + log(1 + exp(-abs(a-b)))

"""
Inplace version of log_add, stores log(x+y) into a and expects a = log(x), b = log(y) as inputs.
"""
function log_add!(a::Vector, i::Int64, b::Float64)
	a[i] = max(a[i],b) + safe_log(1 + exp(-abs(a[i]-b)))
end

function log_add2!(a::Vector, i::Int64, b::Float64)
	if a[i] < b
		a[i] = b + log(1+exp(a[i]-b))
	else
		a[i] = a[i] + log(1+exp(b-a[i]))
	end
end

function log_add!(a::Matrix, i::Int64, j::Int64, b::Float64)
	a[i,j] = max(a[i,j],b) + safe_log(1 + exp(-abs(a[i,j]-b)))
end

function log_add2!(a::Matrix, i::Int64, j::Int64, b::Float64)
	if a[i,j] < b
		a[i,j] = b + log(1+exp(a[i,j]-b))
	else
		a[i,j] = a[i,j] + log(1+exp(b-a[i,j]))
	end
end
"""
computes the sum log(a_1 + ... + a_n) given the array of
[log(a_1), ..., log(a_n)]
"""
function log_sum(a::Array{Float64, 1})
	r = a[1]
	for i in 2:length(a)
		r = log_add(r, a[i])
	end
	r
end

"""
computes log(x-y) given a = log(x) and b = log(y)
"""
function log_subtract(a::Float64,b::Float64)
	if a < b
		return -1000.0
	end
	a + safe_log(1 - exp(b-a))
end

function log_subtract(a::Array{Float64, 1},b::Array{Float64, 1})
	r = zeros(length(a))
	for i in 1:length(a)
		r[i] = log_subtract(a[i], b[i])
	end
	r
end

function log_subtract!(a::Matrix, col::Int64, b::Vector)
	for i in 1:length(b)
		a[i, col] = log_subtract(a[i,col], b)
	end
end

function safe_log(a)
	if a < 0
		return -1000.0
	end
	log(a + 1e-100)
end

"""
Computes m'*diagm(d)*m efficiently.
"""
function matrix_diagonal_product(m::Matrix, d::Vector)
	B = broadcast(.*, m, sqrt(d))
	B'*B
end

function chol_inv(m::Matrix{Float64})
	try
		L_inv = inv(chol(Hermitian(m))')
		return L_inv'*L_inv, prod(diag(L_inv))^2
	catch
		return inv(m), 1./det(m)
	end
end

function inverse_and_determinant(m::Matrix)
	(l,u,p) = lu(m)
	l = LowerTriangular(l)
	u = UpperTriangular(u)
	n = size(m,1)
	e_vec = zeros(n)
	inv_m = zeros(n,n)
	S = 0
	for i in 1:n
		e_vec[i] = 1.
		if i > 1 e_vec[i-1] = 0. end
		if p[i] != i S += 1 end
		y = l \ e_vec
		inv_m[p[i],:] = u\y
	end
	(inv_m, 1./(det(l)*det(u)*(-1)^S))
end

"""
Computes the inverse inv(A+B) in terms of inv(A) and B.
Based on the algorithm described in:
Miller, K. S. (1981). On the Inverse of the Sum of Matrices. Mathematics Magazine, 54(2), 67–72.

The above is simplified to the case where H is diagonal and hence each decomposition
is a multiple of the unit vector in the according dimension.
"""
function matrix_sum_inverse(a_inverse::Matrix, b::Matrix, det_a::Float64)
	#assume that b is a diagonal matrix (H in the paper)
	num_ranks = size(b,1)
	#C_1 is G
	C_k_inv = a_inverse
	result_det = det_a
	for r in 1:num_ranks
		#extract rank 1 decomposition
		# tmp =
		v_k = 1./(1+C_k_inv[r,r]*b[r,r])
		result_det /= v_k
		#prepare for next step
		C_k_inv -= v_k*C_k_inv.*b[:,r]'*C_k_inv
	end
	#return last C_k_inv = inv(G+H) and determinant of the inverse
	(C_k_inv, 1./result_det)
end

function matrix_sum_inverse2(a_inverse::Array{Float64,2}, decomposition::Array{Array{Float64,2},1})
	num_ranks = length(decomposition)
	C_k_inv = a_inverse
	result_det = prod(diag(a_inverse))
	for r in 1:num_ranks
		tmp = C_k_inv*decomposition[r]
		v_k = 1./(1+trace(tmp))
		result_det *= v_k
		C_k_inv -= v_k*tmp*C_k_inv
	end
	(C_k_inv, result_det)
end

function get_topic_distributions(m::GDTM)
	topic_probabilities = Vector{Matrix{Float64}}(m.K)
	for k in 1:m.K
		means = m.KnmKmmInv*m.μ[k]
		topic_probabilities[k] = zeros(m.T, m.V)
		Threads.@threads for w in 1:m.V
			topic_probabilities[k][:,w] = exp(means[:,w] - m.zeta[k,:])
		end
	end
	topic_probabilities
end

function get_topic_distributions(path::String, save_images::Bool=false)
	no_zeta = !isfile(joinpath(path, "zeta.mat"))
	if !no_zeta
		zeta_file = open(joinpath(path, "zeta.mat"), "r")
		zeta = read_array(zeta_file)
		close(zeta_file)
	end
	tp_raw = read_topic_distribution(joinpath(path, "topics.dat"))

	times = map(x->parse(Float64, x), readlines(joinpath(path, "times.dat")))
	inducing_times = map(x->parse(Float64, x), readlines(joinpath(path, "inducing_times.dat")))
	krn_file = open(joinpath(path, "kernel.krn"), "r")
	krn = deserialize(krn_file)
	close(krn_file)
	if length(times) < 50
		test_times = collect(linspace(times[1],times[end], 250))
		times = test_times
		no_zeta=true
	end
	Knm = Kernels.computeCrossCovarianceMatrix(krn, times, inducing_times)
	Kmm = Kernels.computeCovarianceMatrix(krn, inducing_times)
	KnmKmmInv = Knm*pinv(Kmm)

	for k in 1:length(tp_raw)
		tp_raw[k] = KnmKmmInv*tp_raw[k]
	end

	for k in 1:length(tp_raw)
		(T,V) = size(tp_raw[k])
		if no_zeta
			tp_raw[k] = exp(tp_raw[k])
			tp_raw[k] ./= vec(sum(tp_raw[k],2))
		else
			for w in 1:V
				tp_raw[k][:,w] = exp(tp_raw[k][:,w] - zeta[k,:])
			end
		end
	end
	tp_raw
end

function write_topic_words(m::GDTM, corpus::CorpusUtils.Corpus, target_file::String)
	topic_probabilities = get_topic_distributions(m)
	#extract top words for topics
	f = open(target_file, "w")
	for k in 1:m.K
		@printf(f, "topic %d\n", k)
		for t in 1:m.T
			cur_topic = vec(topic_probabilities[k][t,:])
			sorted_idx = sortperm(cur_topic,rev=true)
			@printf(f, "%s:\t", Dates.format(m.timestamps[t], "yyyy-mm-dd"))
			[@printf(f, "%s ", corpus.lexicon.words[word].value) for word in sorted_idx[1:10]]
			@printf(f, "\n")
		end
		@printf(f, "\n")
	end
	close(f)
end

function write_stats(m::GDTM; base_dir_suffix::String="", topics_only::Bool=false)
	base_dir = "experiments/$base_dir_suffix/$(split(string(typeof(m.krn)), ".")[2])$(m.krn.params)/alpha$(m.alpha)/datavar$(m.s_x)_prior_mean$(m.m_0)_prior_var$(m.s_0)/$(m.K)_topics"
	mkpath(base_dir)

	write_topic_words(m, m.corpus, joinpath(base_dir, "topics.txt"))
	if topics_only
		return
	end
	elbo_file = open(joinpath(base_dir,"elbo"), "w")
	lr_file = open(joinpath(base_dir,"learning_rates"), "w")
	for i in 1:length(m.likelihood_counts)
		@printf(elbo_file, "%f\t%f\n", m.likelihood_counts[i], m.likelihoods[i])
		@printf(lr_file, "%f\t%f\n", m.likelihood_counts[i], m.learning_rates[i])
	end
	close(elbo_file)
	close(lr_file)
	test_file = open(joinpath(base_dir,"test"), "w")
	for i in 1:length(m.test_counts)
		@printf(test_file, "%f\t%f\n", m.test_counts[i], m.test_perplexities[i])
	end
	close(test_file)
	times_file = open(joinpath(base_dir, "times.dat"), "w")
	times_to_store = map(Dates.datetime2julian,m.timestamps)
	for t in 1:m.T
		@printf(times_file, "%f\n", times_to_store[t])
	end
	close(times_file)
	inducing_times_file = open(joinpath(base_dir, "inducing_times.dat"), "w")
	for i in m.inducing_points
		@printf(inducing_times_file, "%f\n", i)
	end
	close(inducing_times_file)
	zeta_file = open(joinpath(base_dir, "zeta.mat"), "w")
	write_array(m.zeta, zeta_file)
	close(zeta_file)
	kernel_file = open(joinpath(base_dir, "kernel.krn"), "w")
	serialize(kernel_file, m.krn)
	close(kernel_file)
	write_topic_distribution(m, base_dir)
	write_document_distribution(m, base_dir)
	base_dir
end

function write_topic_distribution(m::GDTM, base_dir::String)
	#write the topics
	f = open(joinpath(base_dir, "topics.dat"), "w")
	#write num topics
	write(f, m.K)
	for k in 1:m.K
		write_array(m.μ[k], f)
	end
	close(f)
end

function read_topic_distribution(td_file::String)
	f = open(td_file, "r")
	K = read(f, Int64)
	tp = Vector{Array}(K)
	for k in 1:K
		tp[k] = read_array(f)
	end
	close(f)
	tp
end

function write_document_distribution(m::GDTM, base_dir::String)
	f = open(joinpath(base_dir, "documents.dat"), "w")
	# get distributions over topics for all documents
	(_, _, _, _, _, document_proportions) = e_step(m, collect(1:m.D))
	# write the array
	write_array(document_proportions, f)
	close(f)
end

function read_document_distribution(dd_file::String)
	f = open(dd_file, "r")
	dd = read_array(f)
	close(f)
	dd
end

function write_array(a, f)
	s = size(a)
	write(f, length(s))
	write(f, collect(s))
	write(f, a)
end

function read_array(f)
	num_dims = read(f, Int64)
	dims = tuple(read(f, Int64, num_dims)...)
	read(f, Float64, dims)
end

function tab_split(file::String)
	f = open(file,"r")
	c = Float64[]
	v = Float64[]
	for l in readlines(f)
		(cur_c, cur_v) = split(chomp(l), '\t')
		push!(c, parse(Float64, cur_c))
		push!(v, parse(Float64, cur_v))
	end
	close(f)
	(c,v)
end


text:
using FileIO
using JLD2
using Statistics
using Printf
using StatsPlots

include("usefulFunctions.jl")

WriteLatexString = false;

# Copy-pasted from "Analytic Continued Fractions for Regression: A Memetic Algorithm Approach"
# The only modification done was to fix the dataset names that were abreviated in the paper
Moscato2021Table3 = "ESL 0.274 0.319 0.272 0.268 fri_c3_1000_10 0.077 0.066 0.062 0.068 SWD 0.390 0.405 0.408 0.432 fri_c0_1000_5 0.042 0.074 0.071 0.106 LEV 0.425 0.424 0.422 0.353 fri_c3_100_5 0.250 0.182 0.093 0.125 ERA 2.514 2.585 2.567 2.450 fri_c1_1000_5 0.049 0.047 0.047 0.050 USCrime 3.93e2 2.57e2 3.78e2 2.20e2 fri_c3_250_5 0.126 0.110 0.108 0.086 FacultySalaries 4.035 8.071 4.111 1.277 fri_c4_250_10 0.148 0.173 0.172 0.133 vineyard 6.010 8.223 7.825 4.222 fri_c4_500_50 0.079 0.141 0.114 0.169 auto_price 5.89e6 3.89e6 4.03e6 6.01e6 fri_c3_500_5 0.103 0.085 0.071 0.078 autoPrice 4.17e6 5.29e6 2.87e6 4.83e6 fri_c3_1000_50 0.068 0.082 0.070 0.120 cloud 0.110 0.208 0.144 0.095 fri_c1_1000_25 0.057 0.070 0.067 0.076 elusage 1.35e2 1.99e2 1.37e2 65.797 fri_c0_100_10 0.149 0.319 0.307 0.226 machine_cpu 3.80e3 2.23e3 2.69e3 2.10e3 fri_c2_1000_50 0.063 0.069 0.076 0.118 analcatdata_vehicle 4.14e4 2.41e4 4.20e4 1.55e4 fri_c4_1000_10 0.050 0.074 0.060 0.063 vinnie 2.287 2.860 2.659 1.934 fri_c0_100_5 0.152 0.193 0.164 0.202 pm10 0.640 0.431 0.399 0.621 fri_c2_500_50 0.049 0.101 0.109 0.130 analcatdata_neavote 1.180 0.818 0.917 0.401 fri_c2_500_10 0.064 0.088 0.097 0.081 analcatdata_election2000 4.33e7 3.40e8 7.72e8 5.09e5 fri_c3_1000_5 0.059 0.049 0.048 0.046 pollution 1.87e3 2.19e3 1.67e3 1.42e3 fri_c1_500_5 0.068 0.077 0.074 0.075 no2 0.272 0.227 0.210 0.295 fri_c0_500_25 0.047 0.132 0.127 0.144 analcatdata_apnea2 1.12e6 9.42e5 7.86e5 6.09e5 fri_c2_100_10 0.750 0.321 0.233 0.320 analcatdata_apnea1 8.16e5 9.98e5 5.28e5 6.96e5 fri_c0_250_10 0.047 0.183 0.170 0.259 cpu 1.75e2 2.36e3 8.83e2 1.64e2 fri_c1_500_50 0.091 0.111 0.133 0.116 fri_c0_250_5 0.065 0.162 0.165 0.176 fri_c1_500_10 0.062 0.086 0.069 0.075 fri_c3_500_25 0.064 0.090 0.097 0.081 fri_c2_500_25 0.077 0.097 0.110 0.121 fri_c1_500_25 0.117 0.108 0.107 0.111 fri_c4_250_25 0.157 0.206 0.178 0.110 fri_c1_1000_50 0.066 0.073 0.077 0.104 fri_c3_500_50 0.101 0.132 0.131 0.156 fri_c4_500_25 0.127 0.111 0.095 0.111 fri_c3_500_10 0.062 0.079 0.066 0.075 fri_c3_1000_25 0.056 0.069 0.059 0.064 fri_c1_250_10 0.114 0.123 0.105 0.130 fri_c4_1000_100 0.047 0.087 0.072 0.227 fri_c1_250_50 0.098 0.171 0.154 0.140 fri_c2_1000_25 0.051 0.064 0.070 0.077 fri_c0_500_5 0.048 0.099 0.106 0.124 fri_c0_1000_50 0.044 0.111 0.113 0.141 fri_c0_500_50 0.053 0.131 0.170 0.211 fri_c1_100_10 0.485 0.245 0.154 0.233 fri_c0_100_25 0.259 0.419 0.383 0.365 fri_c4_1000_25 0.058 0.076 0.066 0.081 fri_c0_250_25 0.064 0.209 0.209 0.224 fri_c1_1000_10 0.047 0.051 0.049 0.059 fri_c0_500_10 0.048 0.153 0.138 0.202 fri_c2_100_5 0.509 0.277 0.255 0.259 fri_c1_100_5 0.490 0.208 0.176 0.241 fri_c0_1000_10 0.047 0.096 0.093 0.155 fri_c2_250_10 0.128 0.123 0.116 0.140 fri_c2_250_5 0.104 0.113 0.094 0.086 fri_c3_250_25 0.116 0.194 0.184 0.111 fri_c2_500_5 0.079 0.068 0.073 0.054 sleuth_ex1714 1.42e6 1.57e6 2.31e6 6.83e5 fri_c0_1000_25 0.039 0.092 0.093 0.074 rabe_266 7.107 7.316 3.043 2.643 fri_c2_1000_5 0.046 0.048 0.045 0.044 sleuth_case2002 75.772 56.240 72.363 41.742 fri_c1_250_5 0.084 0.103 0.093 0.093 rmftsa_ladata 3.012 3.511 3.212 2.764 fri_c3_250_10 0.086 0.150 0.120 0.084 visualizing_environmental 9.624 9.798 9.540 4.700 fri_c0_250_50 0.071 0.238 0.252 0.254 sleuth_ex1605 101.824 98.440 91.986 83.998 fri_c4_500_10 0.065 0.084 0.072 0.072 visualizing_galaxy 3.13e2 2.68e2 2.24e2 4.34e2 fri_c2_250_25 0.179 0.173 0.156 0.149 chatfield_4 2.82e2 3.84e2 2.88e2 1.89e2 fri_c2_1000_10 0.069 0.062 0.058 0.073 sleuth_case1202 3.29e3 3.42e3 3.20e3 1.39e3 fri_c4_1000_50 0.051 0.088 0.068 0.092 chscase_geyser1 38.888 39.902 42.887 31.053";
Moscato2021Table3 = split(Moscato2021Table3, ' ', keepempty=false);
@assert(rem(length(Moscato2021Table3),5)==0);
datasetNamesMoscato2021 = Moscato2021Table3[1:5:end];
datasetNamesMoscato2021 = replace.(datasetNamesMoscato2021, "_"=>" ");
resultsMoscato2021 = parse.(Float64, [Moscato2021Table3[2:5:end] Moscato2021Table3[3:5:end] Moscato2021Table3[4:5:end] Moscato2021Table3[5:5:end]])
function readResultsMoscato2021(datasetName)
    index = findfirst(datasetNamesMoscato2021.==datasetName);
    (index==nothing) && return nothing;
    return resultsMoscato2021[index, :];
end;

resultsDoME = Array{Float64,5}[];
datasets = datasetNames();
finishedExperiments = falses(length(datasets));
for numDatasetName in 1:length(datasets)
    datasetName = datasets[numDatasetName];
    (finished, resultsThisDataset) = experimentsDoMEFinished(datasetName);
    resultsMoscato2021 = readResultsMoscato2021(datasetNameFromFile(datasetName));
    if (finished) && (resultsMoscato2021!=nothing)
        finishedExperiments[numDatasetName] = true;
        push!(resultsDoME, resultsThisDataset);
    else
        println(datasetNameFromFile(datasetName), " ", finished ? "" : "DoME experiments not finished", " ", resultsMoscato2021!=nothing ? "" : "no results found in [Moscato2021]");
    end;
end;
datasets = datasets[finishedExperiments];

function findBestResultsParameters(results, training)
    @assert(!any(isnan.(view(results,:,:,:,:,[1,3,4,5]))));
    meanTrainingResults = mean(results[:,:,:,:,1], dims=4); meanTrainingResults = meanTrainingResults[:,:,:,1];
    stdTrainingResults = std(results[:,:,:,:,1], dims=4); stdTrainingResults = stdTrainingResults[:,:,:,1];
    medianTrainingResults = median(results[:,:,:,:,1], dims=4); medianTrainingResults = medianTrainingResults[:,:,:,1];

    meanTestResults = mean(results[:,:,:,:,2], dims=4); meanTestResults = meanTestResults[:,:,:,1];
    meanTestResults[isnan.(meanTestResults)] .= Inf;
    stdTestResults = std(results[:,:,:,:,2], dims=4); stdTestResults = stdTestResults[:,:,:,1];
    medianTestResults = median(results[:,:,:,:,2], dims=4); medianTestResults = medianTestResults[:,:,:,1];
    medianTestResults[isnan.(medianTestResults)] .= Inf;

    meanHeight = mean(results[:,:,:,:,4], dims=4); meanHeight = meanHeight[:,:,:,1];
    meanNumNodes = mean(results[:,:,:,:,5], dims=4); meanNumNodes = meanNumNodes[:,:,:,1];
    meanExecutionTime = mean(results[:,:,:,:,3], dims=4); meanExecutionTime = meanExecutionTime[:,:,:,1];
    medianHeight = median(results[:,:,:,:,4], dims=4); meanHeight = meanHeight[:,:,:,1];
    medianNumNodes = median(results[:,:,:,:,5], dims=4); meanNumNodes = meanNumNodes[:,:,:,1];
    medianExecutionTime = median(results[:,:,:,:,3], dims=4); meanExecutionTime = meanExecutionTime[:,:,:,1];
    if training
        (bestResult, index) = findmin(medianTrainingResults);
    else
        (bestResult, index) = findmin(medianTestResults);
    end;
    return (medianTrainingResults[index], stdTrainingResults[index], medianTestResults[index], stdTestResults[index], index[1], MinimumReductionsMSE[index[2]], MaxNumNodes[index[3]], medianHeight[index], medianNumNodes[index], medianExecutionTime[index])
end;

findBestTrainingResultsParameters(results) = findBestResultsParameters(results, true);
findBestTestResultsParameters(results) = findBestResultsParameters(results, false);

function fixText(text)
    if !WriteLatexString
        text = replace(text, "\\multirow{4}{*}{" => "");
        text = replace(text, "} & " => "\t");
        text = replace(text, "\\textbf{" => "");
        text = replace(text, " & " => "\t");
        text = replace(text, "& " => "\t\t");
        text = replace(text, "\$" => "");
        text = replace(text, "\\tiny" => "");
        text = replace(text, "\\normalsize" => "");
        text = replace(text, "{" => "(");
        text = replace(text, "}" => ")");
        text = replace(text, ") \\\\" => "");
        text = replace(text, " \\\\" => "");
        text = replace(text, "\\%" => "%");
        text = replace(text, " & " => " ");
    else
        text = replace(text, "e+0" => "e");
    end;
    return text;
end;

println("");
println("Table 5");

for numDataset in 1:length(datasets)
    datasetName = datasetNameFromFile(datasets[numDataset]);
    (medianTrainingResult, stdTrainingResults, medianTestResult, stdTestResults, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, bestHeight, bestNumNodes, bestTime) = findBestTestResultsParameters(resultsDoME[numDataset]);
    testResults = resultsDoME[numDataset][bestStrategy, findfirst(MinimumReductionsMSE.==bestMinimumReductionsMSE), findfirst(MaxNumNodes.==bestMaxNumNodes),:,2];
    stability = 100*median(abs.(testResults .- median(testResults))) / median(testResults);
    trainingResults = resultsDoME[numDataset][bestStrategy, findfirst(MinimumReductionsMSE.==bestMinimumReductionsMSE), findfirst(MaxNumNodes.==bestMaxNumNodes),:,1];
    stability = 100*median(abs.(trainingResults .- median(trainingResults))) / median(trainingResults);
    text = @sprintf("%s & %.3g & %.3g & \$10^{-%d}\$ & %.3g & %.3g & %.2f \\%% \\\\",
        datasetName, medianTestResult, bestStrategy, -log10(bestMinimumReductionsMSE), bestMaxNumNodes, bestTime, stability);
    text = fixText(text);
    println(text);
end;

println("");
println("Figures 5 and 6");

strategyBestMSEResults = zeros(4,length(datasets));
strategyBestTimes = zeros(4,length(datasets));
for numDataset in 1:length(datasets)
    datasetName = datasetNameFromFile(datasets[numDataset]);
    for strategy in 1:4
        newResults = copy(resultsDoME[numDataset]);
        newResults[setdiff(1:4,[strategy]), :, :, :, 2] .= Inf;
        (medianTrainingResult, stdTrainingResults, medianTestResult, stdTestResults, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, bestHeight, bestNumNodes, bestTime) = findBestTestResultsParameters(newResults);
        strategyBestMSEResults[strategy,numDataset] = medianTestResult;
        strategyBestTimes[strategy,numDataset] = bestTime;
    end;
end;

println("Read results from ", size(strategyBestMSEResults,2), " datasets");

strategyBestMSEResults ./= minimum(strategyBestMSEResults, dims=1);
boxplotStrategyMSE = boxplot(strategyBestMSEResults[1,:], label = "Exhaustive");
boxplot!(boxplotStrategyMSE, strategyBestMSEResults[2,:], label = "Exhaustive with constant optimization");
boxplot!(boxplotStrategyMSE, strategyBestMSEResults[3,:], label = "Selective with constant optimization");
boxplot!(boxplotStrategyMSE, strategyBestMSEResults[4,:], label = "Selective");
xaxis!(boxplotStrategyMSE, "Strategy", (0.5, 4.5))
yaxis!(boxplotStrategyMSE, "", (1, 1.75))
display(boxplotStrategyMSE);
savefig(boxplotStrategyMSE,string("results/boxplotStrategyMSE.pdf"));

strategyBestTimes ./= minimum(strategyBestTimes, dims=1);
boxplotStrategyTime = boxplot(strategyBestTimes[1,:], label = "Exhaustive");
boxplot!(boxplotStrategyTime, strategyBestTimes[2,:], label = "Exhaustive with constant optimization");
boxplot!(boxplotStrategyTime, strategyBestTimes[3,:], label = "Selective with constant optimization");
boxplot!(boxplotStrategyTime, strategyBestTimes[4,:], label = "Selective");
xaxis!(boxplotStrategyTime, "Strategy", (0.5, 4.5))
yaxis!(boxplotStrategyTime, "", (1, 100))
display(boxplotStrategyTime);
savefig(boxplotStrategyTime,string("results/boxplotStrategyTime.pdf"));

println("");
println("Table 6");

resultsComparison = zeros(length(datasets),5);
for numDataset in 1:length(datasets)
    datasetName = datasetNameFromFile(datasets[numDataset]);
    (medianTrainingResult, stdTrainingResults, medianTestResult, stdTestResults, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, bestHeight, bestNumNodes, bestTime) = findBestTestResultsParameters(resultsDoME[numDataset]);
    resultsComparisonThisDataset = readResultsMoscato2021(datasetName);
    push!(resultsComparisonThisDataset, medianTestResult);
    resultsComparison[numDataset,:] .= resultsComparisonThisDataset;
    function writeResult(numResult)
        result = resultsComparison[numDataset,numResult];
        if result==minimum(resultsComparisonThisDataset)
            return @sprintf(" & \\textbf{%.3g}", result);
        else
            return @sprintf(" & %.3g", result);
        end;
    end;
    text = @sprintf("%s%s%s%s%s%s \\\\",
        datasetName, writeResult(1), writeResult(2), writeResult(3), writeResult(4), writeResult(5));
    text = fixText(text);
    println(text);
end;

println("");
println("Figure 7");

# https://mirkobunse.github.io/CriticalDifferenceDiagrams.jl/dev/
using CriticalDifferenceDiagrams, DataFrames, PGFPlots

df = DataFrame(dataset_name = repeat(datasets, inner=5), SR_method = repeat(["eplex-1m",  "grad-b",  "xg-b",  "CFR",  "DoME"], size(resultsComparison,1)), MSE = resultsComparison'[:]);

criticalDifferencePlot = CriticalDifferenceDiagrams.plot(
    df,
    :SR_method,
    :dataset_name,
    :MSE;
    maximize_outcome=false,
    reverse_x=true,
    title="Critical Difference Diagram"
)
PGFPlots.save("results/criticalDifference.pdf", criticalDifferencePlot)

println("");
println("Figure 8");

resultsComparison ./= minimum(resultsComparison, dims=2);
boxplotMSEMethods = boxplot(resultsComparison[:,1], label = "eplex-1m");
boxplot!(boxplotMSEMethods, resultsComparison[:,2], label = "grad-b");
boxplot!(boxplotMSEMethods, resultsComparison[:,3], label = "xg-b");
boxplot!(boxplotMSEMethods, resultsComparison[:,4], label = "CFR");
boxplot!(boxplotMSEMethods, resultsComparison[:,5], label = "DoME");
yaxis!(boxplotMSEMethods, "", (1, 3))
display(boxplotMSEMethods);
savefig(boxplotMSEMethods,string("results/boxplotMSEMethods.pdf"));

println("");
println("Table 7");

for numDataset in Int64.(1:round(length(datasets)/2))
    datasetName = datasetNameFromFile(datasets[numDataset]);
    (medianTrainingResult, stdTrainingResults, medianTestResult, stdTestResults, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, bestHeight, bestNumNodes, bestTime) = findBestTestResultsParameters(resultsDoME[numDataset]);
    (inputs,) = loadDataset(datasets[numDataset]); numVariables = size(inputs,2);

    text = @sprintf("%s & %g & %g & %g",
        datasetName, bestNumNodes, bestHeight, numVariables);
    otherDatasetIndex = Int64(numDataset + round(length(datasets)/2));
    if otherDatasetIndex<=length(datasets)
        datasetName = datasetNameFromFile(datasets[otherDatasetIndex]);
        (medianTrainingResult, stdTrainingResults, medianTestResult, stdTestResults, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, bestHeight, bestNumNodes, bestTime) = findBestTestResultsParameters(resultsDoME[otherDatasetIndex]);
        (inputs,) = loadDataset(datasets[otherDatasetIndex]); numVariables = size(inputs,2);
        text = @sprintf("%s & %s & %g & %g & %g",
            text, datasetName, bestNumNodes, bestHeight, numVariables);
    else
        text = @sprintf("%s & & & & ", text);
    end;
    text = @sprintf("%s \\\\", text);
    text = fixText(text);
    println(text);
end;

# Finally, let's take some of the best expressions for some datasets
println("");
println("Table 8");

include("DoME.jl");

datasetNamesFinalExpressions = ["1029_LEV", "1089_USCrime", "192_vineyard", "210_cloud", "228_elusage", "519_vinnie", "522_pm10", "523_analcatdata_neavote", "678_visualizing_environmental", "687_sleuth_ex1605", "712_chscase_geyser1"];

# Redefine this function to write constants with up to three decimals
string(node::Constant) = node.semantic>=0 ? @sprintf("%.3f", node.semantic) : @sprintf("(%.3f)", node.semantic);

for datasetNameFinalExpression in datasetNamesFinalExpressions
    numDataset = findfirst(datasets .== datasetNameFinalExpression);
    (_, _, _, _, bestStrategy, bestMinimumReductionsMSE, bestMaxNumNodes, _, _, _) = findBestTestResultsParameters(resultsDoME[numDataset]);
    if bestMaxNumNodes<=80
        (inputs, targets) = loadDataset(datasets[numDataset]);

        (trainingMSE, validationMSE, testMSE, bestTree) = dome(inputs, targets;
            minimumReductionMSE = bestMinimumReductionsMSE ,
            maximumNodes = bestMaxNumNodes ,
            strategy = [Strategy1, Strategy2, Strategy4, Strategy4][bestStrategy]
        );

        if WriteLatexString
            text = latexString(bestTree);
            # Take the first and last parenthesis
            if (text[1:6]=="\\left(") && (text[end-6:end]=="\\right)")
                text = text[7:end-7];
            end;
            println(datasetNameFromFile(datasets[numDataset]), " & \${", text, "}\$ \\\\");
            println("\\midrule");
        else
            text = string(bestTree);
            println(datasetNameFromFile(datasets[numDataset]), ": ", string(bestTree));
        end;

    end;
end;


text:
# hann.jl
# Compare Matlab and Sound versions
# Note that Matlab has "hann" and "hanning" and they differ!

using Sound: hann
using MATLAB
using Test: @testset, @test

ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"

@testset "hann" begin
    for n in [8, 9] # test odd and even
        eval_string("h1 = hanning($n)")
        eval_string("h2 = hann($n+2)")
        @mget h1
        @mget h2
        @test [0; h1; 0] ≈ h2
        @test h1 ≈ hann(n)
    end
end


text:
export list_guild_emojis,
    get_guild_emoji,
    create_guild_emoji,
    modify_guild_emoji,
    delete_guild_emoji

"""
    list_guild_emojis(c::Client, guild::Integer) -> Vector{Emoji}

Get the [`Emoji`](@ref)s in a [`Guild`](@ref).
"""
list_guild_emojis(c::Client, guild::Integer) = Response{Vector{Emoji}}(c, :GET, "/guilds/$guild/emojis")

"""
    get_guild_emoji(c::Client, guild::Integer, emoji::Integer) -> Emoji

Get an [`Emoji`](@ref) in a [`Guild`](@ref).
"""
get_guild_emoji(c::Client, guild::Integer, emoji::Integer) = Response{Emoji}(c, :GET, "/guilds/$guild/emojis/$emoji")

"""
    create_guild_emoji(c::Client, guild::Integer; kwargs...) -> Emoji

Create an [`Emoji`](@ref) in a [`Guild`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/emoji#create-guild-emoji).
"""
create_guild_emoji(c::Client, guild::Integer; kwargs...) = Response{Emoji}(c, :POST, "/guilds/$guild/emojis"; body = kwargs)

"""
    modify_guild_emoji(c::Client, guild::Integer, emoji::Integer; kwargs...) -> Emoji

Edit an [`Emoji`](@ref) in a [`Guild`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/emoji#modify-guild-emoji).
"""
modify_guild_emoji(c::Client, guild::Integer, emoji::Integer; kwargs...) = Response{Emoji}(c, :PATCH, "/guilds/$guild/emojis/$emoji"; body = kwargs)

"""
    delete_guild_emoji(c::Client, guild::Integer, emoji::Integer)

Delete an [`Emoji`](@ref) from a [`Guild`](@ref).
"""
delete_guild_emoji(c::Client, guild::Integer, emoji::Integer) = Response(c, :DELETE, "/guilds/$guild/emojis/$emoji")


text:
## Test NarmaxRegressors constructor
@testset "Test NarmaxRegressors constructor" begin
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 4, 1:2)
    yterms = [1]
    uterms = [1]
    vterms = [1, 2]
    reg = NarmaxLasso.NarmaxRegressors(l, yterms, uterms, vterms)
    @test reg.n_vterms == 2
    @test reg.n_yterms == 1
    @test reg.n_uterms == 1
    for i = 1:(reg.noise_components_starts-1)
        @test all(reg.basis[i].exponents[3:4] .== 0)
    end
    for i = (reg.noise_components_starts):length(reg.basis)
        @test any(reg.basis[i].exponents[3:4] .> 0)
    end
end

## Test generate lag matrix
@testset "Test generate lag matrix" begin
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 4, 1:2)
    yterms = [1]
    uterms = [1]
    vterms = [1, 2]
    reg = NarmaxLasso.NarmaxRegressors(l, yterms, uterms, vterms)
    y = [1.0, 2.0, 3.0, 4.0]
    u = [10.0, 20.0, 30.0, 40.0]
    v = [100.0, 200.0, 300.0, 400.0]
    M = NarmaxLasso.lag_matrix(reg, y, u, v)
    @test M ≈ [1.0  10.0  100.0  100.0;
               1.0  10.0  100.0  100.0;
               2.0  20.0  200.0  100.0;
               3.0  30.0  300.0  200.0]
end

## Test generate regressor matrix
@testset "Test generate regressor matrix" begin
    l = NarmaxLasso.generate_all(NarmaxLasso.Monomial, 3, 1)
    yterms = [1]
    uterms = [1]
    vterms = [1]
    reg = NarmaxLasso.NarmaxRegressors(l, yterms, uterms, vterms)
    y = [1.0, 2.0, 3.0, 4.0]
    u = [10.0, 20.0, 30.0, 40.0]
    v = [100.0, 200.0, 300.0, 400.0]
    M = NarmaxLasso.lag_matrix(reg, y, u, v)
    @test M ≈ [1.0  10.0  100.0;
               1.0  10.0  100.0;
               2.0  20.0  200.0;
               3.0  30.0  300.0]
end

## Test generate all monomial terms
@testset "Test generate all monomial terms" begin
    reg = NarmaxLasso.generate_all(NarmaxLasso.NarmaxRegressors,
                                NarmaxLasso.Monomial, 2, 2, 2, 1:2)
    @test length(reg.basis) == 27
    reg = NarmaxLasso.generate_all(NarmaxLasso.NarmaxRegressors,
                                NarmaxLasso.Monomial, 2, 2, 2, 1:3)
    @test length(reg.basis) == 83
end


text:
# nextpow2(W) = vshl(one(W), vsub_fast(8sizeof(W), leading_zeros(vsub_fast(W, one(W)))))

# @inline _pick_vector(::StaticInt{W}, ::Type{T}) where {W,T} = Vec{W,T}
# @inline pick_vector(::Type{T}) where {T} = _pick_vector(pick_vector_width(T), T)
# @inline function pick_vector(::Val{N}, ::Type{T}) where {N, T}
#     _pick_vector(smin(nextpow2(StaticInt{N}()), pick_vector_width(T)), T)
# end
# pick_vector(N::Int, ::Type{T}) where {T} = pick_vector(Val(N), T)

@inline MM(::Union{Val{W},StaticInt{W}}) where {W} = MM{W}(0)
@inline MM(::Union{Val{W},StaticInt{W}}, i) where {W} = MM{W}(i)
@inline MM(::Union{Val{W},StaticInt{W}}, i::AbstractSIMDVector{W}) where {W} = i
@inline MM(::StaticInt{W}, i, ::StaticInt{X}) where {W,X} = MM{W,X}(i)
@inline gep(ptr::Ptr, i::MM) = gep(ptr, data(i))

@inline Base.one(::Type{MM{W,X,I}}) where {W,X,I} = one(I)
@inline staticm1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vsub_fast(data(i), one(I)))
@inline staticp1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vadd_nsw(data(i), one(I)))
@inline vadd_fast(i::MM{W,X}, j::IntegerTypesHW) where {W,X} = MM{W,X}(vadd_fast(data(i), j))
@inline vadd_fast(i::IntegerTypesHW, j::MM{W,X}) where {W,X} = MM{W,X}(vadd_fast(i, data(j)))
@inline vadd_fast(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd_fast(data(i), StaticInt{j}()))
@inline vadd_fast(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd_fast(StaticInt{i}(), data(j)))
@inline vadd_fast(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
@inline vadd_fast(::StaticInt{0}, j::MM{W,X}) where {W,X} = j
@inline vsub_fast(i::MM{W,X}, j::IntegerTypesHW) where {W,X} = MM{W,X}(vsub_fast(data(i), j))
@inline vsub_fast(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub_fast(data(i), StaticInt{j}()))
@inline vsub_fast(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i

@inline vadd_nsw(i::MM{W,X}, j::IntegerTypesHW) where {W,X} = MM{W,X}(vadd_nsw(data(i), j))
@inline vadd_nsw(i::IntegerTypesHW, j::MM{W,X}) where {W,X} = MM{W,X}(vadd_nsw(i, data(j)))
@inline vadd_nsw(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd_nsw(data(i), StaticInt{j}()))
@inline vadd_nsw(i::MM, ::Zero) = i
@inline vadd_nsw(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd_nsw(StaticInt{i}(), data(j)))
@inline vadd_nsw(::Zero, j::MM{W,X}) where {W,X} = j
@inline vsub_nsw(i::MM{W,X}, j::IntegerTypesHW) where {W,X} = MM{W,X}(vsub_nsw(data(i), j))
@inline vsub_nsw(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub_nsw(data(i), StaticInt{j}()))
@inline vsub(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
@inline vsub_fast(i::MM) = i * StaticInt{-1}()
@inline vsub_nsw(i::MM) = i * StaticInt{-1}()
@inline vsub(i::MM) = i * StaticInt{-1}()
@inline vadd(i::MM, j::IntegerTypesHW) = vadd_fast(i, j)
@inline vadd(j::IntegerTypesHW, i::MM) = vadd_fast(j, i)
@inline vsub(i::MM, j::IntegerTypesHW) = vsub_fast(i, j)
@inline vsub(j::IntegerTypesHW, i::MM) = vsub_fast(j, i)
@inline vadd(i::MM, ::StaticInt{j}) where {j} = vadd_fast(i, StaticInt{j}())
@inline vadd(::StaticInt{j}, i::MM) where {j} = vadd_fast(StaticInt{j}(), i)
@inline vsub(i::MM, ::StaticInt{j}) where {j} = vsub_fast(i, StaticInt{j}())
@inline vsub(::StaticInt{j}, i::MM) where {j} = vsub_fast(StaticInt{j}(), i)
@inline vadd(i::MM, ::Zero) = i
@inline vadd(::Zero, i::MM) = i
@inline vsub(::Zero, i::MM) = StaticInt{-1}()*i

@inline vmul_nsw(::StaticInt{M}, i::MM{W,X}) where {M,W,X} = MM{W}(vmul_nsw(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())
@inline vmul_nsw(i::MM{W,X}, ::StaticInt{M}) where {M,W,X} = MM{W}(vmul_nsw(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())

@inline vmul_fast(::StaticInt{M}, i::MM{W,X}) where {M,W,X} = MM{W}(vmul_fast(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())
@inline vmul_fast(i::MM{W,X}, ::StaticInt{M}) where {M,W,X} = MM{W}(vmul_fast(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())
@inline vmul(a, ::StaticInt{N}) where {N} = vmul_fast(a, StaticInt{N}())
@inline vmul(::StaticInt{N}, a) where {N} = vmul_fast(StaticInt{N}(), a)
@inline vmul(::StaticInt{N}, ::StaticInt{M}) where {N,M} = StaticInt{N}() * StaticInt{M}()

@inline vsub(a, ::StaticInt{N}) where {N} = vsub_fast(a, StaticInt{N}())
@inline vadd(a, ::StaticInt{N}) where {N} = vadd_fast(a, StaticInt{N}())
@inline vsub(::StaticInt{N}, a) where {N} = vsub_fast(StaticInt{N}(), a)
@inline vadd(::StaticInt{N}, a) where {N} = vadd_fast(StaticInt{N}(), a)
@inline vsub(::StaticInt{M}, ::StaticInt{N}) where {M,N} = StaticInt{M}() - StaticInt{N}()
@inline vadd(::StaticInt{M}, ::StaticInt{N}) where {M,N} = StaticInt{M}() + StaticInt{N}()

@inline vrem(i::MM{W,X,I}, ::Type{I}) where {W,X,I<:IntegerTypesHW} = i
@inline vrem(i::MM{W,X}, ::Type{I}) where {W,X,I<:IntegerTypesHW} = MM{W,X}(data(i) % I)
@inline veq(::AbstractIrrational, ::MM{W,<:Integer}) where {W} = zero(Mask{W})
@inline veq(x::AbstractIrrational, i::MM{W}) where {W} = x == Vec(i)
@inline veq(::MM{W,<:Integer}, ::AbstractIrrational) where {W} = zero(Mask{W})
@inline veq(i::MM{W}, x::AbstractIrrational) where {W} = Vec(i) == x

@inline vsub_nsw(i::MM, ::Zero) = i
@inline vsub_nsw(i::NativeTypes, j::MM{W,X}) where {W,X} = MM(StaticInt{W}(), vsub_nsw(i, data(j)), -StaticInt{X}())
@inline vsub_fast(i::NativeTypes, j::MM{W,X}) where {W,X} = MM(StaticInt{W}(), vsub_fast(i, data(j)), -StaticInt{X}())
@inline vsub_nsw(i::Union{FloatingTypes,IntegerTypesHW}, j::MM{W,X}) where {W,X} = MM(StaticInt{W}(), vsub_nsw(i, data(j)), -StaticInt{X}())
@inline vsub_fast(i::Union{FloatingTypes,IntegerTypesHW}, j::MM{W,X}) where {W,X} = MM(StaticInt{W}(), vsub_fast(i, data(j)), -StaticInt{X}())

@inline function Base.in(m::MM{W,X,<:Integer}, r::AbstractUnitRange) where {W,X}
    vm = Vec(m)
    (vm ≥ first(r)) & (vm ≤ last(r))
end

for op ∈ (:(+),:(-))
  @eval begin
    @inline Base.$op(vu::VecUnroll{N,1,T,T}, i::MM) where {N,T<:NativeTypes} = VecUnroll(fmap($op, data(vu), i))
    @inline Base.$op(i::MM, vu::VecUnroll{N,1,T,T}) where {N,T<:NativeTypes} = VecUnroll(fmap($op, i, data(vu)))
    @inline Base.$op(vu::VecUnroll{N,1,T,T}, i::VecUnroll{N,W,I,MM{W,X,I}}) where {N,W,T<:NativeTypes,I,X} = VecUnroll(fmap($op, data(vu), data(i)))
    @inline Base.$op(i::VecUnroll{N,W,I,MM{W,X,I}}, vu::VecUnroll{N,1,T,T}) where {N,W,T<:NativeTypes,I,X} = VecUnroll(fmap($op, data(i), data(vu)))
  end
end
# @inline Base.:(+)(vu::VecUnroll{N,W,T,T}, i::MM) = VecUnroll(fmap(+, data(vu), i))


text:
# issue 462

using Gadfly
using Base.Dates

a = [unix2datetime(100)]
b = [10]

plot(x=a, y=b, Geom.point)


text:
"""
  genGaussian(mean_in, cov_in, class_label, n)

Draw from a gaussian distribution
# Arguments
- `mean_in` : means
- `cov_in` : covariances
- `class_label` : class_label
- `n` : number of samples to draw
"""
function genGaussian(mean_in::Array, cov_in::Array, class_label::Int64, n::Int64)
  nv = Distributions.MvNormal(mean_in, cov_in)
  X = rand(nv, n)
  y = ones(Int64, n) .* class_label
  return nv, transpose(X), y
end

"""
  genZafarData(n = 10000; d = pi/4)

Generate synthetic data from Zafar et al., 2017 Fairness Constraints: Mechanisms for Fair Classification.
Code translated to julia from
https://github.com/mbilalzafar/fair-classification/blob/master/disparate_impact/synthetic_data_demo/generate_synthetic_data.py

# Arguments
- `n=10000` : number of samples
- `d=pi/4` : discrimination factor

# Returns
- `X` : DataFrame containing features and protected attribute z {"A", "B"} where
        z="B" is the protected group.
- `y` : Binary Target variable {-1, 1}
"""
function genZafarData(n = 10000; d = pi/4)
    mu1, sigma1 = [2., 2.],  [[5., 1.] [1., 5.]]
    mu2, sigma2 = [-2.,-2.], [[10., 1.]  [1., 3.]]
    npos = Int64(floor(n/2))
    nv1, X1, y1 = genGaussian(mu1, sigma1,  1, npos)     # positive class
    nv2, X2, y2 = genGaussian(mu2, sigma2, -1, n - npos) # negative class
    perm = shuffle(1:n)
    X = vcat(X1, X2)[perm,:]
    y = vcat(y1, y2)[perm]
    rotation_mult = [[cos(d), -sin(d)] [sin(d), cos(d)]]
    X_aux = X * rotation_mult
    """ Generate the sensitive feature here """
    x_control = Array{String, 1}(undef, n) # this array holds the sensitive feature value
    for i in 1:n
        x = X_aux[i,:]
        # probability for each cluster that the point belongs to it
        p1 = Distributions.pdf(nv1, x)
        p2 = Distributions.pdf(nv2, x)
        # normalize the probabilities from 0 to 1
        s = p1+p2
        p1 = p1/s # p2 = p2/s
        if rand(1)[1] < p1
            x_control[i] = "A" # majority class
        else
            x_control[i] = "B" # protected class
        end
    end
    X = DataFrame(X)
    X.z = categorical(x_control)
    coerce!(X, :z => Multiclass)
    y = categorical(y)
    return X, y
end

"""
Generate synthetic data from Zafar et al., 2017 Fairness Beyond Disparate Treatment & Disparate Impact
# Arguments
- `n=10000` : number of samples
# Returns
- `X` : DataFrame containing features and protected attribute z
- `y` : Binary Target variable
"""
function genZafarData2(n = 10000)

    y = rand(Distributions.Bernoulli(0.5), n)
    z = rand(Distributions.Bernoulli(0.5), n)
    # Depending on z and y, we draw from one of the 4 distributions in M.
    M = [
      [Distributions.MvNormal([1., 1.],    [[3., 1.] [1., 3.]]), Distributions.MvNormal([2., 2.], [[3., 1.] [1., 3.]])],
      [Distributions.MvNormal([-2., -2.],  [[3., 1.] [1., 3.]]), Distributions.MvNormal([2., 2.], [[3., 1.] [1., 3.]])]
      ]

    # Iterate over z and y in parallel, drawing from the appropritate Distribution.
    for i in 1:n
      if i == 1
        X = rand(M[z[i]+1][y[i]+1], 1)
      else
        X = hcat(X, rand(M[z[i]+1][y[i]+1], 1))
      end
    end

    perm = shuffle(1:n)
    X = DataFrame(transpose(X)[perm,:])
    y = y[perm]
    X.z = [x == 0 ? "A" : "B" for x in z[perm]]
    coerce!(X, :z => Multiclass)
    y = categorical(y)
    return X, y
end

"""
  genSubgroupData(n=10000, setting="B00")

Generate synthetic data from Loh et al., 2019 : Subgroup identification for precision medicine: A comparative review of 13 methods
# Arguments
- `n=10000` : number of samples
- `setting="B00"` : Simulation data setting: one of "B00", ..., "B02", "B1", ... , "B8"
For "B00", ..., "B02" there is no "bias" in the data, i.e. group membership has no effect on y.
whereas for  "B1", ... , "B8", there is a direct effect of group membership z on y, usually mediated by one or more features.
# Returns
- `X` : DataFrame containing features and protected attribute z
- `y` : Binary target variable
"""
function genSubgroupData(n = 10000; setting = "B00")
    xdists = [
      Distributions.Normal(0.,1.), # x1 is N(0,1)
      Distributions.MvNormal([0., 0.], [[1, 0.5] [0.5, 1]]), # x2 and x3 are N(0,1) with cor 0.5
      Distributions.Exponential(1.), # x4 is exp(1)
      Distributions.Bernoulli(0.5), # x5 is Ber(0.5)
      Distributions.Categorical(repeat([.1], 10)), # x6 is Multinom. w. equal
      Distributions.MvNormal([0., 0., 0., 0.], [[1, 0.5, 0.5, 0.5] [0.5, 1, 0.5, 0.5] [0.5, 0.5, 1, 0.5] [0.5, 0.5, 0.5, 1]]
      ) # x7 to 10 are N(0,1) with cor 0.5
    ]
    z = [x == 0 ? "A" : "B" for x in rand(Distributions.Bernoulli(0.5), n)]
    X = []
    for d in xdists
      if length(X) > 0
        xn = rand(d, n)
        if size(xn)[1] != n
          xn = transpose(xn)
        end
        X = hcat(X, xn)
      else
        X = rand(d, n)
      end
    end
    y = logit_fun(X, z, setting)
    X = DataFrame(X)
    X.z = z
    coerce!(X, :z => Multiclass)
    y = categorical(y)
    return X, y
end

# Logistic Link function
function log_link(x)
  return exp(x)/(1+exp(x))
end

"""
  logit_fun(X, z, setting)

Compute y from X and z according to a setting provided in Loh et al., 2019: Subgroup identification for precision medicine: A comparative review of 13 methods
# Arguments
- `X` : matrix of features
- `z` : vector of group assignments
- `setting` : Simulation data setting: one of "B00", ..., "B02", "B1", ... , "B8"
"""
function logit_fun(X, z, setting)
  if setting == "B00"
    logit = repeat([0.], length(z))
  elseif setting == "B01"
    logit = 0.5*(X[:,1] + X[:,2])
  elseif setting == "B02"
    logit = 0.5*(X[:,1] + X[:,2]) .^ 2
  elseif setting == "B1"
    logit = 0.5*(X[:,1] + X[:,2] - X[:,5]) + 2. * z .* map(x -> ifelse(isodd(Int64(x)), 1., 0.), X[:,6])
  elseif setting == "B2"
    logit = 0.5*X[:,2] + 2. * z .* map(x -> ifelse(x > 0, 1., 0.), X[:,1])
  elseif setting == "B3"
    logit = 0.3*(X[:,1] + X[:,2]) + 2. * z .* map(x -> ifelse(x > 0, 1., 0.), X[:,1])
  elseif setting == "B4"
    logit = 0.2*(X[:,2] + X[:,3] .- 2.) + 2. * z .* X[:,4]
  elseif setting == "B5"
    logit = 0.2*(X[:,1] + X[:,2] .- 3) + 2. * z .* map(x -> ifelse(x < 0, 1., 0.), X[:,1]) .* map(x -> ifelse(isodd(Int64(x)), 1., 0.), X[:,6])
  elseif setting == "B6"
    logit = 0.5*(X[:,2] .- 1.) + 2. * z .* map(x -> ifelse(abs(x) < 0.8, 1., 0.), X[:,1])
  elseif setting == "B7"
    logit = 0.2*(X[:,2] + X[:,2] .^ 2 .- 6.)  + 2. * z .* map(x -> ifelse(x < 0, 1., 0.), X[:,1])
  elseif setting == "B8"
    logit = 0.5 * X[:,2] + 2. * z .* X[:,5]
  end
  yprob = log_link.(logit)
  # Compute a random y with probability yprob
  u = rand(length(yprob))
  return yprob .> u
end

"""
  genBiasedSampleData(n=10000, sampling_bias=0.8)

Generate synthetic data: Biased sample

# Arguments
- `n=10000` : number of samples
- `sampling_bias=0.8` : Percentage of data belonging to majority group.

The idea behind this simulation is that algorithms might fit the process in the
majority group while disregarding the process in the minority group.

Two different processes for d1 and d2:
d1: logit(y) = 0.5*(    X1 + X2 + 0.3*X4) + 2*(I(X3 > 0))
d2: logit(y) = 0.5*(0.3*X1 + X2 +     X4) + 2*(I(X3 > 0.2))

# Returns
- `X` : DataFrame containing features and protected attribute z
- `y` : Binary Target variable
"""
function genBiasedSampleData(n = 10000; sampling_bias = 0.8)

  d1 = Distributions.MvNormal([0., 0., 0., 0.], [[1, 0.3, 0.3, 0.3] [0.3, 1, 0.3, 0.3] [0.3, 0.3, 1, 0.3] [0.3, 0.3, 0.3, 1]])
  d2 = Distributions.MvNormal([0.2, 0.2, 0.2, 0.2], [[1, 0.5, 0.5, 0.5] [0.5, 1, 0.5, 0.5] [0.5, 0.5, 1, 0.5] [0.5, 0.5, 0.5, 1]])

  s1 = Int64(floor(n * sampling_bias))

  X1 = transpose(rand(d1, s1))
  X2 = transpose(rand(d2, n-s1))

  perm = shuffle(1:n)
  yprob = vcat(
    log_link.(0.5 * (       X1[:,1] + X1[:,2] + 0.3 * X1[:,4]) + map(x -> ifelse(x > 0, 2., 0.), X1[:,3])),
    log_link.(0.5 * (0.3 .* X2[:,2] + X2[:,2] +       X2[:,4]) + map(x -> ifelse(x > .2, 2., 0.), X2[:,3]))
  )[perm]
  X = vcat(X1, X2)[perm,:]
  z = vcat(repeat(["A"], s1), repeat(["B"], n-s1))[perm]

  X = DataFrame(X)
  X.z = z
  coerce!(X, :z => Multiclass)
  u = rand(n)
  y = categorical(ifelse.(yprob .> u, 1, 0))
  return X, y
end