text:
BUILTINS = []

function builtin (name_list)
    function add(fn)
        for name in name_list
            BUILTINS.append((name, fn, names[0]))
        end
        return fn
    end
    return add
end

function scheme_stringp(x)
    return (typeof(x) == String) and x.startswith('"')
end

scheme_stringp = builtin("string?")(scheme_stringp)

function scheme_symbolp(x):
    return (typeof(x) == String) and ~scheme_stringp(x)
end

scheme_symbolp = builtin("boolean?")(scheme_symbolp)

function validate_type(val, predicate, k, name):
    """Returns VAL.  Raises a SchemeError if not PREDICATE(VAL)
    using "argument K of NAME" to describe the offending value."""
    if ~predicate(val):
        type = type(val)
        if scheme_symbolp(val):
            type = Symbol
        end
        raise SchemeError("argument $k of $name has wrong type $type_name")
    return val

function scheme_booleanp(x)
    return x == True || x == False
end

scheme_booleanp = builtin("boolean?")(scheme_booleanp)

function is_true_primitive(val)
    return val != False
end

function is_false_primitive(val)
    return val == False
end

function scheme_not(x)
    return ~is_true_primitive(x)
end

scheme_not = builtin("not")(scheme_not)

function scheme_pairp(x)
    return typeof(x) == Pair #TODO see if this can be overloaded
end
scheme_pairp = builtin("pair?")(scheme_pairp)

function scheme_nullp(x)
    return type(x) == Void
end
scheme_nullp = builtin("null?")(scheme_nullp)

function scheme_promisep(x)
    return type(x) == Promise
end
scheme_promisep = builtin("promise?")(scheme_promisep)

function scheme_valid_cdrp(x)
    return scheme_pairp(x) || scheme_nullp(x) || scheme_promisep(x)
end
scheme_valid_cdrp = builtin("pair?")(scheme_valid_cdrp)

function scheme_equalp(x, y)
    if scheme_pairp(x) && scheme_pairp(y):
        return scheme_equalp(x.first, y.first) && scheme_equalp(x.rest, y.rest)
    elseif scheme_numberp(x) && scheme_numberp(y)
        return x == y
    else
        return type(x) == type(y) && x == y
    end
end

scheme_equalp = builtin("equal?")(scheme_equalp)

function scheme_eqvp(x)
    if scheme_pairp(x) && scheme_pairp(y)
        return x == y
    elseif scheme_symbolp(x) && scheme_symbolp(y)
        return x == y
    else
        return x === y
    end
end
scheme_eqvp = builtin("eqv?")(scheme_eqvp)

"eq?")
function scheme_eqp(x, y)
    if scheme_symbolp(x) && scheme_symbolp(y)
        return x == y
    else
        return x === y
    end
end
scheme_eqp = builtin("eq?")(scheme_eqp)

# Streams

function scheme_force(x)
    validate_type(x, scheme_promisep, 0, 'promise')
    return x.evaluate() #TODO what is evaluate; it's not eval
end
scheme_force = builtin("force?")(scheme_force)

function scheme_cdr_stream(x)
    lambda_x = scheme_pairp(x) and scheme_promisep(x.rest)
    validate_type(x, lambda_x, 0, 'cdr-stream')
    return scheme_force(x.rest)
end
scheme_cdr_stream = builtin("cdr_stream?")(scheme_cdr_stream)

function scheme_listp(x)
    """Return whether x is a well-formed list. Assumes no cycles."""
    while x is not Void
        if not isinstance(x, Pair)
            return False
        end
        x = x.rest
    end
    return True
end
scheme_listp = builtin("list?")(scheme_listp)

function scheme_length(x)
    validate_type(x, scheme_listp, 0, 'length')
    if x is Void
        return 0
    end
    return len(x)
end
scheme_length = builtin("length?")(scheme_length)

function scheme_cons(x, y)
    return Pair(x, y)
end
scheme_cons = builtin("cons?")(scheme_cons)

function scheme_car(x)
    validate_type(x, scheme_pairp, 0, 'car')
    return x.first
end
scheme_car = builtin("car?")(scheme_car)

function scheme_cdr(x)
    validate_type(x, scheme_pairp, 0, 'cdr')
    return x.rest
end
scheme_cdr = builtin("cdr?")(scheme_cdr)

function scheme_set_car(x, y)
    validate_type(x, scheme_pairp, 0, 'set-car!')
    x.first = y
end
scheme_set_car = builtin("set_car?")(scheme_set_car)

function scheme_set_cdr(x, y)
    validate_type(x, scheme_pairp, 0, 'set-cdr!')
    validate_type(y, scheme_valid_cdrp, 1, 'set-cdr!')
    x.rest = y
end
scheme_set-cdr = builtin("set-cdr?")(scheme_set-cdr)

"list")
function scheme_list(*vals)
    result = nil
    for e in reversed(vals)
        result = Pair(e, result)
    return result
end
scheme_p = builtin("?")(scheme_p)

"append")
function scheme_append(*vals)
    if len(vals) == 0:
        return nil
    result = vals[-1]
    for i in range(len(vals)-2, -1, -1)
        v = vals[i]
        if v is not nil:
            validate_type(v, scheme_pairp, i, 'append')
            r = p = Pair(v.first, result)
            v = v.rest
            while scheme_pairp(v)
                p.rest = Pair(v.first, result)
                p = p.rest
                v = v.rest
            result = r
    return result
end
scheme_p = builtin("?")(scheme_p)

"string?")

"number?")
function scheme_numberp(x)
    return isinstance(x, numbers.Real) and not scheme_booleanp(x)
scheme_p = builtin("?")(scheme_p)
end

"integer?")
function scheme_integerp(x)
    return scheme_numberp(x) and (isinstance(x, numbers.Integral) or int(x) == x)
scheme_p = builtin("?")(scheme_p)
end

function _check_nums(*vals)
    """Check that all arguments in VALS are numbers."""
    for i, v in enumerate(vals)
        if not scheme_numberp(v)
            msg = "operand {0} ({1}) is not a number"
            raise SchemeError(msg.format(i, v))
end

function _arith(fn, init, vals)
    """Perform the FN operation on the number values of VALS, with INIT as
    the value when VALS is empty. Returns the result as a Scheme value."""
    _check_nums(*vals)
    s = init
    for val in vals:
        s = fn(s, val)
    s = _ensure_int(s)
    return s
end

function _ensure_int(x)
    if int(x) == x:
        x = int(x)
    return x
end

"+")
function scheme_add(*vals)
    return _arith(operator.add, 0, vals)
end
scheme_p = builtin("?")(scheme_p)

"-")
function scheme_sub(val0, *vals)
    _check_nums(val0, *vals) # fixes off-by-one error
    if len(vals) == 0:
        return _ensure_int(-val0)
    return _arith(operator.sub, val0, vals)
end
scheme_p = builtin("?")(scheme_p)

"*")
function scheme_mul(*vals)
    return _arith(operator.mul, 1, vals)
end
scheme_p = builtin("?")(scheme_p)

"/")
function scheme_div(val0, *vals)
    _check_nums(val0, *vals) # fixes off-by-one error
    try:
        if len(vals) == 0:
            return _ensure_int(operator.truediv(1, val0))
        return _arith(operator.truediv, val0, vals)
    except ZeroDivisionError as err:
        raise SchemeError(err)
end
scheme_p = builtin("?")(scheme_p)

"expt")
function scheme_expt(val0, val1)
    _check_nums(val0, val1)
    return pow(val0, val1)
end
scheme_p = builtin("?")(scheme_p)

"abs")
function scheme_abs(val0)
    return abs(val0)
end
scheme_p = builtin("?")(scheme_p)

"quotient")
function scheme_quo(val0, val1)
    _check_nums(val0, val1)
    try:
        return -(-val0 // val1) if (val0 < 0) ^ (val1 < 0) else val0 // val1
    except ZeroDivisionError as err:
        raise SchemeError(err)

end
scheme_p = builtin("?")(scheme_p)

"modulo")
function scheme_modulo(val0, val1)
    _check_nums(val0, val1)
    try:
        return val0 % val1
    except ZeroDivisionError as err:
        raise SchemeError(err)

end
scheme_p = builtin("?")(scheme_p)

"remainder")
function scheme_remainder(val0, val1)
    _check_nums(val0, val1)
    try:
        result = val0 % val1
    except ZeroDivisionError as err:
        raise SchemeError(err)
    while result < 0 and val0 > 0 or result > 0 and val0 < 0:
        result -= val1
    return result
end
scheme_p = builtin("?")(scheme_p)

function number_fn(module, name, fallback=None)
    """A Scheme built-in procedure that calls the numeric Python function named
    MODULE.FN."""
    py_fn = getattr(module, name) if fallback is None else getattr(module, name, fallback)
    function scheme_fn(*vals)
        _check_nums(*vals)
        return py_fn(*vals)
    return scheme_fn
end

for _name in ["acos", "acosh", "asin", "asinh", "atan", "atan2", "atanh",
              "ceil", "copysign", "cos", "cosh", "degrees", "floor", "log",
              "log10", "log1p", "radians", "sin", "sinh", "sqrt",
              "tan", "tanh", "trunc"]:
    builtin(_name)(number_fn(math, _name))
builtin("log2")(number_fn(math, "log2", lambda x: math.log(x, 2)))  # Python 2 compatibility

# Add number functions in the math module as built-in procedures in Scheme

function _numcomp(op, x, y)
    _check_nums(x, y)
    return op(x, y)
end

"=")
function scheme_eq(x, y)
    return _numcomp(operator.eq, x, y)
end
scheme_p = builtin("?")(scheme_p)

"<")
function scheme_lt(x, y)
    return _numcomp(operator.lt, x, y)
end
scheme_p = builtin("?")(scheme_p)

">")
function scheme_gt(x, y)
    return _numcomp(operator.gt, x, y)
end
scheme_p = builtin("?")(scheme_p)

"<=")
function scheme_le(x, y)
    return _numcomp(operator.le, x, y)
end
scheme_p = builtin("?")(scheme_p)

">=")
function scheme_ge(x, y)
    return _numcomp(operator.ge, x, y)
end
scheme_p = builtin("?")(scheme_p)

"even?")
function scheme_evenp(x)
    _check_nums(x)
    return x % 2 == 0
end
scheme_p = builtin("?")(scheme_p)

"odd?")
function scheme_oddp(x)
    _check_nums(x)
    return x % 2 == 1
end
scheme_p = builtin("?")(scheme_p)

"zero?")
function scheme_zerop(x)
    _check_nums(x)
    return x == 0
end
scheme_p = builtin("?")(scheme_p)

##
## Other operations
##

"atom?")
function scheme_atomp(x)
    return (scheme_booleanp(x) or scheme_numberp(x) or scheme_symbolp(x) or
            scheme_nullp(x) or scheme_stringp(x))

end
scheme_p = builtin("?")(scheme_p)

"display")
function scheme_display(*vals)
    vals = [repl_str(val[1:-1] if scheme_stringp(val) else val) for val in vals]
    print(*vals, end="")
end
scheme_p = builtin("?")(scheme_p)

"print")
function scheme_print(*vals)
    vals = [repl_str(val) for val in vals]
    print(*vals)
end
scheme_p = builtin("?")(scheme_p)

"displayln")
function scheme_displayln(*vals)
    scheme_display(*vals)
    scheme_newline()
end
scheme_p = builtin("?")(scheme_p)

"newline")
function scheme_newline()
    print()
    sys.stdout.flush()
end
scheme_p = builtin("?")(scheme_p)

"error")
function scheme_error(msg=None)
    msg = "" if msg is None else repl_str(msg)
    raise SchemeError(msg)
end
scheme_p = builtin("?")(scheme_p)

"exit")
function scheme_exit()
    raise EOFError
    return TkCanvas(1000, 1000, init_hook=_title)
end

scheme_p = builtin("?")(scheme_p)


text:
module NFP
using GLM, StatsBase, DataFrames, MultivariateStats, CSV
using Combinatorics, Parameters, DataFramesMeta, Lazy, Plots
pyplot()

struct Loss
    mVar::Array{Int64}
    mScore::Array{Float64}
    mNames::Array{String}
end

function Loss(m::Array,mScore::Array,iBest::Int64,vNames::Vector)
    mBest = m[:,1:iBest]
    mBscore =  mScore[mBest]
    mBnames = vNames[mBest]
    return Loss(mBest,mBscore,mBnames)
end

struct Score
    vLoss::Array{Loss}
end

struct UnivariateSelection
    vVar::Vector{Int64}
    vNames::Vector{String}
    s::Score
end

struct Results
    U::UnivariateSelection
    Best_comb
    Factor_in
end

struct Results_f
    U::Array{String}
    Best_comb
    Factor_in
    m3dFore
end

function UnivariateSelection(mX::Array,vY::Vector,vNames,H::Vector{Int64},iStart::Int64,iBest::Int64,fLoss::Array{Function})
    F = length(fLoss)
    mFore,mLossScore = get_score(mX,vY,H,iStart,fLoss)

    mLossScore_b = [Loss(get_variable_idx(mLossScore[:,:,f]::Array,false),mLossScore[:,:,f],iBest,vNames) for f = 1:F]
    s = Score(mLossScore_b)

    vBest = get_best(mLossScore_b,iBest,fLoss)
    vBest_names = vNames[vBest]

    return UnivariateSelection(vBest,vBest_names,s)
end

function UnivariateSelection(mX::Array,vY::Vector,vNames,H::Vector{Int64},iStart::Int64,iBest::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    F = length(fLoss)
    mFore,mLossScore = get_score(mX,vY,H,iStart,l,fLoss)

    mLossScore_b = [Loss(get_variable_idx(mLossScore[:,:,f]::Array,false),mLossScore[:,:,f],iBest,vNames) for f = 1:F]
    s = Score(mLossScore_b)

    vBest = get_best(mLossScore_b,iBest,fLoss)
    vBest_names = vNames[vBest]

    return UnivariateSelection(vBest,vBest_names,s)
end

function get_best(mLossScore_b::Array{Loss},iBest::Int64,fLoss::Array{Function})
    F = length(fLoss)
    vBest = unique(sort(vcat([unique(mLossScore_b[f].mVar) for f = 1:F]...)))
end

function get_score(mX::Array{Float64,2},vY::Vector,H::Vector{Int64},iStart::Int64,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    F = length(fLoss)
    # Out of sample
    mFore = zeros(T,P,K)
    mLossScore  = zeros(P,K,F)
    @simd for j = 1:K
        println(j)
        mFore_s   = zeros(T,P)
        vLossScore  = zeros(P,F)
        count = 0
        for h = H
            count = count+1
            mFore_s[:,count] = out_of_sample_forecast(mX[:,j:j],vY,iStart,h)
            [vLossScore[count,f] = fLoss[f](vY[iStart+h:end],mFore_s[iStart+h:end,count]) for f = 1:F]
        end
        mFore[:,:,j] = mFore_s
        mLossScore[:,j,:] = vLossScore
    end
    return mFore,mLossScore
end

function get_score(mX::Array{Float64,2},vY::Vector,H::Vector{Int64},iStart::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    F = length(fLoss)
    # Out of sample
    mFore = zeros(T,P,K)
    mLossScore  = zeros(P,K,F)
    @simd for j = 1:K
        println(j)
        mFore_s   = zeros(T,P)
        vLossScore  = zeros(P,F)
        count = 0
        for h = H
            count = count+1
            mFore_s[:,count] = out_of_sample_forecast(mX[:,j:j],vY,iStart,h,l)
            [vLossScore[count,f] = fLoss[f](vY[iStart+h:end],mFore_s[iStart+h:end,count]) for f = 1:F]
        end
        mFore[:,:,j] = mFore_s
        mLossScore[:,j,:] = vLossScore
    end
    return mFore,mLossScore
end

function out_of_sample_forecast(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64)
    T,K = size(mX)::Tuple{Int64,Int64}
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:K)
        vY_n = view(vY,(1:j))
        t = size(vY_n,1)
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,1:K])
        vYY_n = vY_n[h+1:t]
        md = fit(LinearModel,mXX_n,vYY_n)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function out_of_sample_forecast(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,l::GLM.ProbitLink)
    T,K = size(mX)::Tuple{Int64,Int64}
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:K)
        vY_n = view(vY,(1:j))
        t = size(vY_n,1)
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,1:K])
        vYY_n = vY_n[h+1:t]
        md = GLM.fit(GLM.GeneralizedLinearModel,mXX_n,vYY_n,Binomial())
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function get_variable_idx(mLoss_s::Array,bRev::Bool)
    mLoss = zeros(Int64, size(mLoss_s))
    for i = 1:size(mLoss_s,1)
        mLoss[i,:] = sortperm(mLoss_s[i,:],rev = bRev)
    end
    return mLoss
end

function par_get_best_comb(mX::Array,vY::Vector,H::Vector{Int64},U::UnivariateSelection,iStart::Int64,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    N = size(U.vVar,1)
    F = length(fLoss)
    counter = 0
    for comb = 1:size(U.vVar,1)
        println("Starting models with $comb variables")
        comb_index = get_comb(U.vVar,comb)
        n = size(comb_index,1)::Int64
        mLossScore = SharedArray{Float64}(P,n,F)
        mLossScore_pc = SharedArray{Float64}(P,n,F)
        @time @sync @parallel for j = 1:n
            j%10 == 0 && println("#Combination: $j")
            mFore,mFore_pc = zeros(Float64,T,P),zeros(Float64,T,P)
            vLossScore     = zeros(Float64,P,F)
            vLossScore_pc  = zeros(Float64,P,F)
            count = 0
            for h = H
                count += 1
                mFore[:,count]    = out_of_sample_forecast(mX[:,comb_index[j,:]],vY,iStart,h)
                [vLossScore[count,f] = fLoss[f](mFore[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
                mFore_pc[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,comb_index[j,:],1:K)
                [vLossScore_pc[count,f] = fLoss[f](mFore_pc[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
            end
            mLossScore[:,j,:] = vLossScore
            mLossScore_pc[:,j,:] = vLossScore_pc
        end

        for f = 1:F
            dfLossScore = get_df_score([mLossScore[:,:,f] mLossScore[:,:,f]], n)
            filename = joinpath(Pkg.dir("NFP"),"test","$(string(fLoss[f]))_best_$(size(U.vVar,1))_comb_$comb.csv")
            writetable(filename,dfLossScore, header = true)
        end
    end
end

function par_get_best_comb(mX::Array,vY::Vector,H::Vector{Int64},U::UnivariateSelection,iStart::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    N = size(U.vVar,1)
    F = length(fLoss)
    counter = 0
    for comb = 1:size(U.vVar,1)
        println("Starting models with $comb variables")
        comb_index = get_comb(U.vVar,comb)
        n = size(comb_index,1)::Int64
        mLossScore = SharedArray{Float64}(P,n,F)
        mLossScore_pc = SharedArray{Float64}(P,n,F)
        @time @sync @parallel for j = 1:n
            j%10 == 0 && println("#Combination: $j")
            mFore,mFore_pc = zeros(Float64,T,P),zeros(Float64,T,P)
            vLossScore     = zeros(Float64,P,F)
            vLossScore_pc  = zeros(Float64,P,F)
            count = 0
            for h = H
                count += 1
                mFore[:,count]    = out_of_sample_forecast(mX[:,comb_index[j,:]],vY,iStart,h,l)
                [vLossScore[count,f] = fLoss[f](mFore[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
                mFore_pc[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,comb_index[j,:],1:K,l)
                [vLossScore_pc[count,f] = fLoss[f](mFore_pc[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
            end
            mLossScore[:,j,:] = vLossScore
            mLossScore_pc[:,j,:] = vLossScore_pc
        end

        for f = 1:F
            dfLossScore = get_df_score([mLossScore[:,:,f] mLossScore[:,:,f]], n)
            filename = joinpath(Pkg.dir("NFP"),"test","$(string(fLoss[f]))_best_$(size(U.vVar,1))_comb_$comb.csv")
            writetable(filename,dfLossScore, header = true)
        end
    end
end

function get_df_score(m::DenseArray,n::Int64)
    mm = convert(Array,m)::Array
    dfM =  DataFrame(mm)
    names!(dfM,repmat(map(Symbol,1:n),2),makeunique = true)
    return dfM
end

function get_comb(vVar::Vector{Int64},comb::Int64)
    N = size(vVar,1)
    aComb_idx = collect(combinations(1:N, comb))
    vComb_idx = hcat(aComb_idx...)'::Array{Int64}
    return vVar[vComb_idx]::Array{Int64}
end

function out_of_sample_forecast_pca(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,vVar_idx::Vector{Int64},rPC_idx::UnitRange{Int64})
    T,k = size(mX)::Tuple{Int64,Int64}
    K = size(vVar_idx,1)
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:k)
        vY_n = view(vY,(1:j))
        mPC = pca(mX_n[:,rPC_idx])
        t = size(vY_n,1)
        vYY_n = vY_n[h+1:t]
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,vVar_idx], mPC[1:t-h,1:1])
        md = GLM.lm(mXX_n,vYY_n)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K), view(mPC,t:t,1:1))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function out_of_sample_forecast_pca(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,vVar_idx::Vector{Int64},rPC_idx::UnitRange{Int64},l::GLM.ProbitLink)
    T,k = size(mX)::Tuple{Int64,Int64}
    K = size(vVar_idx,1)
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:k)
        vY_n = view(vY,(1:j))
        mPC = pca(mX_n[:,rPC_idx])
        t = size(vY_n,1)
        vYY_n = vY_n[h+1:t]
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,vVar_idx], mPC[1:t-h,1:1])
        md = GLM.fit(GLM.GeneralizedLinearModel,mXX_n,vYY_n,Binomial(),l)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K), view(mPC,t:t,1:1))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function pca(data::Array{Float64}; n::Int64 = 1)
    # Need MultivariateStats
    mX = data.-mean(data,1)
    M = fit(PCA, mX, method = :cov)
    return projection(M)[:,1:n]
end

function get_df_score(m::DenseArray,n::Int64)
    mm = convert(Array,m)::Array
    dfM =  DataFrame(mm)
    names!(dfM,repmat(map(Symbol,1:n),2),makeunique = true)
    return dfM
end

function get_independent(dfData::DataFrame, vVar::Vector{Symbol})
    mX = @>> dfData[.~[(x in vVar) for x in names(dfData)]] convert(Array{Float64,2})
    vNames = @> dfData[.~[(x in vVar) for x in names(dfData)]] names()
    return mX, vNames
end

function get_dependent(vY::Vector,fThreshold::Float64)
    T = size(vY,1)
    vY_prob = zeros(Int64,T)
    for i = 1:T
        vY[i] < fThreshold ? (vY_prob[i] = 1) : (vY_prob[i] = 0)
    end
    return vY_prob
end

function get_dependent(vY::Vector,vThreshold::Vector)
    T = size(vY,1)
    vY_prob = zeros(Int64,T)
    for i = 1:T
        vY[i] < vThreshold[i] ? (vY_prob[i] = 1) : (vY_prob[i] = 0)
    end
    return vY_prob
end

function load_score(sCrit::String,ncomb_load::Int64,H::Vector,vNames::Vector{String})
    ncomb = size(vNames,1)
    vCrit = H
    comb_identifier = "Models" # preallocate comb_identifier for concatenation
    for comb = 1:ncomb_load
        comb_index = collect(combinations(1:size(vNames,1), comb))
        comb_index = transpose(hcat(comb_index...))
        dfCrit = readtable(joinpath(Pkg.dir("NFP"),"test","$(sCrit)_best_$(ncomb)_comb_$comb.csv"), header =  true)
        comb_identifier = @>> map(string,names(dfCrit)).*"_comb$comb" vcat(comb_identifier)
        vCrit = @>> convert(Array,dfCrit) hcat(vCrit)
        println(size(vCrit))
    end
    best = [[findmin(vCrit[j,2:end])[2]+1 findmin(vCrit[j,2:end])[1]] for j in 1:size(vCrit,1)] # +1 exclude first column and reset indexes
    best = vcat(best...)
    best_comb = @>> comb_identifier[convert(Array{Int64},best[:,1])] map(x -> [Symbol(split(x,"_c")[1]) split(x,"_comb")[2]])
    best_comb = vcat(best_comb...)
    factor_in = [sum(collect(string(best_comb[i,1])) .== '_') for i in 1:size(best_comb,1)]
    return best_comb, factor_in
end

function variable_selection(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,fLoss::Array{Function})
    mX,vNames = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64, 1},dfData[iSymbol])
    U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,fLoss)
    const ncomb = size(U.vVar,1)
    print("In the second step you are using $ncomb different variables \n")
    # sPath folder is the folder where the Results are saved and has to be created
    par_get_best_comb(mX,vY,H,U,iStart,fLoss)
    return mX,vY,U,vNames
end

function variable_selection(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    mX,vNames = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64, 1},dfData[iSymbol])
    U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    const ncomb = size(U.vVar,1)
    print("In the second step you are using $ncomb different variables \n")
    # sPath folder is the folder where the Results are saved and has to be created
    par_get_best_comb(mX,vY,H,U,iStart,l,fLoss)
    return mX,vY,U,vNames
end

function check_binary(vY::Array)
    T = size(vY,1)
    if sum(isequal.(vY,0) .| isequal.(vY,1)) != T
        error("The dependent variable must be binary")
    end
end

function sforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,fLoss::Array{Function})
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1))
    rm_var = size(vSymbol,1)
    mX,vY,U,vNames = variable_selection(dfData,vSymbol,iSymbol,H,iStart,iBest,fLoss)
    #line_st = Dict(best_comb_vMae => (:solid,"green",2), best_comb_vRmse => (:dash, "red",2))
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,U.vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(U.vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2))
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && plot!(l_plot,dta,vY[iStart+h:end], linecolor = "black",
            yaxis = true, legend = false, subplot = count, yticks = true, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    CSV.write(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), DataFrame(names = r.U.vNames), header = false)
    return l_plot,r
end

function fforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,fLoss::Array{Function})
    vNames = readtable(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), header = false)
    vNames = string.(vNames[1])
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1));
    rm_var = size(vSymbol,1)
    mX, vNames_full = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64,1}, dfData[iSymbol])
    # U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2))
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && plot!(l_plot,dta,vY[iStart+h:end], linecolor = "black",
            yaxis = true, legend = false, subplot = count, yticks = true, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    return l_plot,r
end

function sforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1))
    rm_var = size(vSymbol,1)
    mX,vY,U,vNames = variable_selection(dfData,vSymbol,iSymbol,H,iStart,iBest,l,fLoss)
    check_binary(vY::Array)
    #line_st = Dict(best_comb_vMae => (:solid,"green",2), best_comb_vRmse => (:dash, "red",2))
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,U.vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(U.vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2),l::GLM.ProbitLink)
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h,l::GLM.ProbitLink)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && bar!(l_plot,dta,vY[iStart+h:end], fill = (0,"lightgray"), linecolor = "lightgrey",
            yaxis = false, legend = false, subplot = count, yticks = false, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    CSV.write(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), DataFrame(names = r.U.vNames), header = false)
    return l_plot,r
end

function fforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    vNames = readtable(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), header = false)
    vNames = string.(vNames[1])
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1))
    rm_var = size(vSymbol,1)
    mX, vNames_full = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64,1}, dfData[iSymbol])
    check_binary(vY::Array)
    # U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2),l::GLM.ProbitLink)
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h,l::GLM.ProbitLink)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && bar!(l_plot,dta,vY[iStart+h:end], fill = (0,"lightgray"), linecolor = "lightgrey",
            yaxis = false, legend = false, subplot = count, yticks = false, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    return l_plot,r
end

using GLM, StatsBase, DataFrames, MultivariateStats
using Combinatorics, Parameters, DataFramesMeta, Lazy, Plots

struct Loss
    mVar::Array{Int64}
    mScore::Array{Float64}
    mNames::Array{String}
end

function Loss(m::Array,mScore::Array,iBest::Int64,vNames::Vector)
    mBest = m[:,1:iBest]
    mBscore =  mScore[mBest]
    mBnames = vNames[mBest]
    return Loss(mBest,mBscore,mBnames)
end

struct Score
    vLoss::Array{Loss}
end

struct UnivariateSelection
    vVar::Vector{Int64}
    vNames::Vector{String}
    s::Score
end

struct Results
    U::UnivariateSelection
    Best_comb
    Factor_in
end

function UnivariateSelection(mX::Array,vY::Vector,vNames,H::Vector{Int64},iStart::Int64,iBest::Int64,fLoss::Array{Function})
    F = length(fLoss)
    mFore,mLossScore = get_score(mX,vY,H,iStart,fLoss)

    mLossScore_b = [Loss(get_variable_idx(mLossScore[:,:,f]::Array,false),mLossScore[:,:,f],iBest,vNames) for f = 1:F]
    s = Score(mLossScore_b)

    vBest = get_best(mLossScore_b,iBest,fLoss)
    vBest_names = vNames[vBest]

    return UnivariateSelection(vBest,vBest_names,s)
end

function UnivariateSelection(mX::Array,vY::Vector,vNames,H::Vector{Int64},iStart::Int64,iBest::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    F = length(fLoss)
    mFore,mLossScore = get_score(mX,vY,H,iStart,l,fLoss)

    mLossScore_b = [Loss(get_variable_idx(mLossScore[:,:,f]::Array,false),mLossScore[:,:,f],iBest,vNames) for f = 1:F]
    s = Score(mLossScore_b)

    vBest = get_best(mLossScore_b,iBest,fLoss)
    vBest_names = vNames[vBest]

    return UnivariateSelection(vBest,vBest_names,s)
end

function get_best(mLossScore_b::Array{Loss},iBest::Int64,fLoss::Array{Function})
    F = length(fLoss)
    vBest = unique(sort(vcat([unique(mLossScore_b[f].mVar) for f = 1:F]...)))
end

function get_score(mX::Array{Float64,2},vY::Vector,H::Vector{Int64},iStart::Int64,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    F = length(fLoss)
    # Out of sample
    mFore = zeros(T,P,K)
    mLossScore  = zeros(P,K,F)
    @simd for j = 1:K
        println(j)
        mFore_s   = zeros(T,P)
        vLossScore  = zeros(P,F)
        count = 0
        for h = H
            count = count+1
            mFore_s[:,count] = out_of_sample_forecast(mX[:,j:j],vY,iStart,h)
            [vLossScore[count,f] = fLoss[f](vY[iStart+h:end],mFore_s[iStart+h:end,count]) for f = 1:F]
        end
        mFore[:,:,j] = mFore_s
        mLossScore[:,j,:] = vLossScore
    end
    return mFore,mLossScore
end

function get_score(mX::Array{Float64,2},vY::Vector,H::Vector{Int64},iStart::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    F = length(fLoss)
    # Out of sample
    mFore = zeros(T,P,K)
    mLossScore  = zeros(P,K,F)
    @simd for j = 1:K
        println(j)
        mFore_s   = zeros(T,P)
        vLossScore  = zeros(P,F)
        count = 0
        for h = H
            count = count+1
            mFore_s[:,count] = out_of_sample_forecast(mX[:,j:j],vY,iStart,h,l)
            [vLossScore[count,f] = fLoss[f](vY[iStart+h:end],mFore_s[iStart+h:end,count]) for f = 1:F]
        end
        mFore[:,:,j] = mFore_s
        mLossScore[:,j,:] = vLossScore
    end
    return mFore,mLossScore
end

function out_of_sample_forecast(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64)
    T,K = size(mX)::Tuple{Int64,Int64}
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:K)
        vY_n = view(vY,(1:j))
        t = size(vY_n,1)
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,1:K])
        vYY_n = vY_n[h+1:t]
        md = fit(LinearModel,mXX_n,vYY_n)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function out_of_sample_forecast(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,l::GLM.ProbitLink)
    T,K = size(mX)::Tuple{Int64,Int64}
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:K)
        vY_n = view(vY,(1:j))
        t = size(vY_n,1)
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,1:K])
        vYY_n = vY_n[h+1:t]
        md = GLM.fit(GLM.GeneralizedLinearModel,mXX_n,vYY_n,Binomial())
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function get_variable_idx(mLoss_s::Array,bRev::Bool)
    mLoss = zeros(Int64, size(mLoss_s))
    for i = 1:size(mLoss_s,1)
        mLoss[i,:] = sortperm(mLoss_s[i,:],rev = bRev)
    end
    return mLoss
end

function par_get_best_comb(mX::Array,vY::Vector,H::Vector{Int64},U::UnivariateSelection,iStart::Int64,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    N = size(U.vVar,1)
    F = length(fLoss)
    counter = 0
    for comb = 1:size(U.vVar,1)
        println("Starting models with $comb variables")
        comb_index = get_comb(U.vVar,comb)
        n = size(comb_index,1)::Int64
        mLossScore = SharedArray{Float64}(P,n,F)
        mLossScore_pc = SharedArray{Float64}(P,n,F)
        @time @sync @parallel for j = 1:n
            j%10 == 0 && println("#Combination: $j")
            mFore,mFore_pc = zeros(Float64,T,P),zeros(Float64,T,P)
            vLossScore     = zeros(Float64,P,F)
            vLossScore_pc  = zeros(Float64,P,F)
            count = 0
            for h = H
                count += 1
                mFore[:,count]    = out_of_sample_forecast(mX[:,comb_index[j,:]],vY,iStart,h)
                [vLossScore[count,f] = fLoss[f](mFore[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
                mFore_pc[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,comb_index[j,:],1:K)
                [vLossScore_pc[count,f] = fLoss[f](mFore_pc[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
            end
            mLossScore[:,j,:] = vLossScore
            mLossScore_pc[:,j,:] = vLossScore_pc
        end

        for f = 1:F
            dfLossScore = get_df_score([mLossScore[:,:,f] mLossScore[:,:,f]], n)
            filename = joinpath(Pkg.dir("NFP"),"test","$(string(fLoss[f]))_best_$(size(U.vVar,1))_comb_$comb.csv")
            writetable(filename,dfLossScore, header = true)
        end
    end
end

function par_get_best_comb(mX::Array,vY::Vector,H::Vector{Int64},U::UnivariateSelection,iStart::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    T,K = size(mX)::Tuple{Int64,Int64}
    P = size(H,1)
    N = size(U.vVar,1)
    F = length(fLoss)
    counter = 0
    for comb = 1:size(U.vVar,1)
        println("Starting models with $comb variables")
        comb_index = get_comb(U.vVar,comb)
        n = size(comb_index,1)::Int64
        mLossScore = SharedArray{Float64}(P,n,F)
        mLossScore_pc = SharedArray{Float64}(P,n,F)
        @time @sync @parallel for j = 1:n
            j%10 == 0 && println("#Combination: $j")
            mFore,mFore_pc = zeros(Float64,T,P),zeros(Float64,T,P)
            vLossScore     = zeros(Float64,P,F)
            vLossScore_pc  = zeros(Float64,P,F)
            count = 0
            for h = H
                count += 1
                mFore[:,count]    = out_of_sample_forecast(mX[:,comb_index[j,:]],vY,iStart,h,l)
                [vLossScore[count,f] = fLoss[f](mFore[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
                mFore_pc[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,comb_index[j,:],1:K,l)
                [vLossScore_pc[count,f] = fLoss[f](mFore_pc[iStart+h:end,count],vY[iStart+h:end]) for f = 1:F]
            end
            mLossScore[:,j,:] = vLossScore
            mLossScore_pc[:,j,:] = vLossScore_pc
        end

        for f = 1:F
            dfLossScore = get_df_score([mLossScore[:,:,f] mLossScore[:,:,f]], n)
            filename = joinpath(Pkg.dir("NFP"),"test","$(string(fLoss[f]))_best_$(size(U.vVar,1))_comb_$comb.csv")
            writetable(filename,dfLossScore, header = true)
        end
    end
end

function get_df_score(m::DenseArray,n::Int64)
    mm = convert(Array,m)::Array
    dfM =  DataFrame(mm)
    names!(dfM,repmat(map(Symbol,1:n),2),makeunique = true)
    return dfM
end

function get_comb(vVar::Vector{Int64},comb::Int64)
    N = size(vVar,1)
    aComb_idx = collect(combinations(1:N, comb))
    vComb_idx = hcat(aComb_idx...)'::Array{Int64}
    return vVar[vComb_idx]::Array{Int64}
end

function out_of_sample_forecast_pca(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,vVar_idx::Vector{Int64},rPC_idx::UnitRange{Int64})
    T,k = size(mX)::Tuple{Int64,Int64}
    K = size(vVar_idx,1)
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:k)
        vY_n = view(vY,(1:j))
        mPC = pca(mX_n[:,rPC_idx])
        t = size(vY_n,1)
        vYY_n = vY_n[h+1:t]
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,vVar_idx], mPC[1:t-h,1:1])
        md = GLM.lm(mXX_n,vYY_n)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K), view(mPC,t:t,1:1))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function out_of_sample_forecast_pca(mX::Array{Float64,2},vY::Vector,iStart::Int64,h::Int64,vVar_idx::Vector{Int64},rPC_idx::UnitRange{Int64},l::GLM.ProbitLink)
    T,k = size(mX)::Tuple{Int64,Int64}
    K = size(vVar_idx,1)
    vFore = zeros(T)
    for j = iStart:T-h
        mX_n = view(mX,1:j,1:k)
        vY_n = view(vY,(1:j))
        mPC = pca(mX_n[:,rPC_idx])
        t = size(vY_n,1)
        vYY_n = vY_n[h+1:t]
        mXX_n = hcat(ones(t-h,1), mX_n[1:t-h,vVar_idx], mPC[1:t-h,1:1])
        md = GLM.fit(GLM.GeneralizedLinearModel,mXX_n,vYY_n,Binomial(),l)
        iXX_n = hcat(ones(1,1), view(mX_n,t:t,1:K), view(mPC,t:t,1:1))
        vFore[j+h] = first(GLM.predict(md,iXX_n))
    end
    return vFore
end

function pca(data::Array{Float64}; n::Int64 = 1)
    # Need MultivariateStats
    mX = data.-mean(data,1)
    M = fit(PCA, mX, method = :cov)
    return projection(M)[:,1:n]
end

function get_df_score(m::DenseArray,n::Int64)
    mm = convert(Array,m)::Array
    dfM =  DataFrame(mm)
    names!(dfM,repmat(map(Symbol,1:n),2),makeunique = true)
    return dfM
end

function get_independent(dfData::DataFrame, vVar::Vector{Symbol})
    mX = @>> dfData[.~[(x in vVar) for x in names(dfData)]] convert(Array{Float64,2})
    vNames = @> dfData[.~[(x in vVar) for x in names(dfData)]] names()
    return mX, vNames
end

function get_dependent(vY::Vector,fThreshold::Float64)
    T = size(vY,1)
    vY_prob = zeros(Int64,T)
    for i = 1:T
        vY < fThreshold ? (vY_prob[i] = 1) : (vY_prob[i] = 0)
    end
    return vY_prob
end

function get_dependent(vY::Vector,vThreshold::Vector)
    T = size(vY,1)
    vY_prob = zeros(Int64,T)
    for i = 1:T
        vY[i] < vThreshold[i] ? (vY_prob[i] = 1) : (vY_prob[i] = 0)
    end
    return vY_prob
end

function load_score(sCrit::String,ncomb_load::Int64,H::Vector,vNames::Vector{String})
    ncomb = size(vNames,1)
    vCrit = H
    comb_identifier = "Models" # preallocate comb_identifier for concatenation
    for comb = 1:ncomb_load
        comb_index = collect(combinations(1:size(vNames,1), comb))
        comb_index = transpose(hcat(comb_index...))
        dfCrit = readtable(joinpath(Pkg.dir("NFP"),"test","$(sCrit)_best_$(ncomb)_comb_$comb.csv"), header =  true)
        comb_identifier = @>> map(string,names(dfCrit)).*"_comb$comb" vcat(comb_identifier)
        vCrit = @>> convert(Array,dfCrit) hcat(vCrit)
        println(size(vCrit))
    end
    best = [[findmin(vCrit[j,2:end])[2]+1 findmin(vCrit[j,2:end])[1]] for j in 1:size(vCrit,1)] # +1 exclude first column and reset indexes
    best = vcat(best...)
    best_comb = @>> comb_identifier[convert(Array{Int64},best[:,1])] map(x -> [Symbol(split(x,"_c")[1]) split(x,"_comb")[2]])
    best_comb = vcat(best_comb...)
    factor_in = [sum(collect(string(best_comb[i,1])) .== '_') for i in 1:size(best_comb,1)]
    return best_comb, factor_in
end

function variable_selection(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,fLoss::Array{Function})
    mX,vNames = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64, 1},dfData[iSymbol])
    U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,fLoss)
    const ncomb = size(U.vVar,1)
    print("In the second step you are using $ncomb different variables \n")
    # sPath folder is the folder where the Results are saved and has to be created
    par_get_best_comb(mX,vY,H,U,iStart,fLoss)
    return mX,vY,U,vNames
end

function variable_selection(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    mX,vNames = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64, 1},dfData[iSymbol])
    U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    const ncomb = size(U.vVar,1)
    print("In the second step you are using $ncomb different variables \n")
    # sPath folder is the folder where the Results are saved and has to be created
    par_get_best_comb(mX,vY,H,U,iStart,l,fLoss)
    return mX,vY,U,vNames
end

function check_binary(vY::Array)
    T = size(vY,1)
    if sum(isequal.(vY,0) .| isequal.(vY,1)) != T
        error("The dependent variable must be binary")
    end
end

function sforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,fLoss::Array{Function})
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1))
    rm_var = size(vSymbol,1)
    mX,vY,U,vNames = variable_selection(dfData,vSymbol,iSymbol,H,iStart,iBest,fLoss)
    #line_st = Dict(best_comb_vMae => (:solid,"green",2), best_comb_vRmse => (:dash, "red",2))
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,U.vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(U.vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2))
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && plot!(l_plot,dta,vY[iStart+h:end], linecolor = "black",
            yaxis = true, legend = false, subplot = count, yticks = true, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    CSV.write(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), DataFrame(names = r.U.vNames), header = false)
    return l_plot,r
end

function fforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,fLoss::Array{Function})
    vNames = readtable(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), header = false)
    vNames = string.(vNames[1])
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(iH,1));
    rm_var = size(vSymbol,1)
    mX, vNames_full = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64,1}, dfData[iSymbol])
    # U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    m3dFore =  zeros(size(vY,1),size(H,1),F)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2))
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && plot!(l_plot,dta,vY[iStart+h:end], linecolor = "black",
            yaxis = true, legend = false, subplot = count, yticks = true, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results_f(vNames,Best_comb,Factor_in,m3dFore)
    return l_plot,r
end

function sforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(length(H),1))
    rm_var = size(vSymbol,1)
    mX,vY,U,vNames = variable_selection(dfData,vSymbol,iSymbol,H,iStart,iBest,l,fLoss)
    check_binary(vY::Array)
    #line_st = Dict(best_comb_vMae => (:solid,"green",2), best_comb_vRmse => (:dash, "red",2))
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,U.vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(U.vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = U.vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2),l::GLM.ProbitLink)
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h,l::GLM.ProbitLink)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && bar!(l_plot,dta,vY[iStart+h:end], fill = (0,"lightgray"), linecolor = "lightgrey",
            yaxis = false, legend = false, subplot = count, yticks = false, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        m3dFore[:,:,f] = mFore
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results(U,Best_comb,Factor_in)
    CSV.write(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), DataFrame(names = r.U.vNames), header = false)
    return l_plot,r
end

function fforecast(dfData::DataFrame,vSymbol::Array{Symbol,1},iSymbol::Symbol,H::Vector,iStart::Int64,iBest::Int64,ncomb_load::Int64,l::GLM.ProbitLink,fLoss::Array{Function})
    vNames = readtable(joinpath(Pkg.dir("NFP"),"test","best_univariate_names.csv"), header = false)
    vNames = string.(vNames[1])
    F = length(fLoss)
    iH = length(H)
    l_plot = plot(layout = grid(iH,1))
    rm_var = size(vSymbol,1)
    mX, vNames_full = get_independent(dfData,vSymbol)
    vY = convert(Array{Float64,1}, dfData[iSymbol])
    check_binary(vY::Array)
    # U = UnivariateSelection(mX,vY,vNames,H,iStart,iBest,l,fLoss)
    counter = 0
    Factor_in = zeros(Int64, iH,F)
    Best_comb = Array{Array}(F,1)
    m3dFore =  zeros(size(vY,1),size(H,1),F)
    for f = 1:F
        counter += 1
        best_comb,factor_in = load_score(string(fLoss[f]),ncomb_load,H,vNames)
        model_index = Array{Vector}(size(best_comb,1))
        variables = Array{Vector}(size(best_comb,1))
        mFore    = zeros(size(vY,1),size(H,1))
        count = 0
        for i = 1:size(best_comb,1)
            comb_index = collect(combinations(1:size(vNames,1), parse(Int64,best_comb[i,2])))
            comb_index = transpose(hcat(comb_index...))
            if @> sum(collect(string(best_comb[i,1])) .== '_') == 0
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),"x")[2]),:]]
            else
                model_index[i] = comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]
                variables[i] = vNames[comb_index[parse(Int64,split(string(best_comb[i,1]),r"x|_")[2]),:]]
            end
            model_index[i] = map(x-> find(x .== map(string,names(dfData)))[1],variables[i])-rm_var
            println("Starting models with $(variables[i]) variables and factor $(Bool(factor_in[i]))")
            h = H[i]
            count += 1
            if factor_in[i] == 1 # test if best use factors
                @inbounds mFore[:,count] = out_of_sample_forecast_pca(mX,vY,iStart,h,model_index[i],1:size(mX,2),l::GLM.ProbitLink)
            else
                @inbounds mFore[:,count] = out_of_sample_forecast(mX[:,model_index[i]],vY,iStart,h,l::GLM.ProbitLink)
            end
            dta = Dates.format(dfData[:date][iStart+h:end], "yyyy-mm-dd")
            f == 1 && bar!(l_plot,dta,vY[iStart+h:end], fill = (0,"lightgray"), linecolor = "lightgrey",
            yaxis = false, legend = false, subplot = count, yticks = false, ylabel = "percentage %")
            plot!(l_plot, dta, mFore[iStart+h:end,count],
            title = "Horizon $h", subplot = count, legend = false, ylabel = string(iSymbol))
        end
        m3dFore[:,:,f] = mFore
        Best_comb[f] = best_comb
        Factor_in[:,f] = factor_in
    end
    r = Results_f(vNames,Best_comb,Factor_in,m3dFore)
    return l_plot,r
end

export sforecast, fforecast

end # module ends

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2021, DLR Institute of System Dynamics and Control

module ModiaPlot_GLMakie

# Constants
const headingSize = 10

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.3.3"
const Date = "2021-11-23"

println("Importing ModiaPlot_GLMakie Version $Version ($Date) - this takes some time due to GLMakie import")

import ModiaResult
import Colors
import Measurements
import MonteCarloMeasurements
using  Unitful

using  GLMakie
include("$(ModiaResult.path)/src/plot.jl")

const showFigureStringInDiagram = true
const callDisplayFunction = true
const reusePossible = true

const Makie_Point2f = isdefined(GLMakie, :Point2f) ? Point2f : Point2f0
include("$(ModiaResult.path)/src/makie.jl")

function showFigure(figureNumber::Int)::Nothing
    #println("... in showFigure")
    if haskey(figures, figureNumber)
        matrixFigure = figures[figureNumber]
        fig = matrixFigure.fig
        display(fig)
    else
        @warn "showFigure: figure $figureNumber is not defined."
    end
    return nothing
end

function saveFigure(figureNumber::Int, fileName; kwargs...)::Nothing
    #println("... in saveFigure")
    if haskey(figures, figureNumber)
        fig = figures[figureNumber].fig
        fullFileName = joinpath(pwd(), fileName)
        println("... save plot in file: \"$fullFileName\"")
        save(fileName, fig; kwargs...)
        display(fig)
    else
        @warn "saveFigure: figure $figureNumber is not defined."
    end
    return nothing
end

function closeFigure(figureNumber::Int)::Nothing
    #println("... in closeFigure")
    delete!(figures,figureNumber)
    if length(figures) > 0
        dictElement = first(figures)
        display(dictElement[2].fig)
    else
        fig = Figure()
        display(fig)
    end
    return nothing
end

"""
    closeAllFigures()

Close all figures.
"""
function closeAllFigures()::Nothing
    #println("... in closeAllFigures")
    if length(figures) > 0
        empty!(figures)
        fig = Figure()
        display(fig)
    end
    return nothing
end

end


text:
# m6.4s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
scale!(df, [:kcal_per_g, :perc_fat, :perc_lactose])

stan6_4 = "
data{
    int <lower=1> N;              // Sample size
    vector[N] K;
    vector[N] L;
}
parameters{
    real a;
    real bL;
    real<lower=0> sigma;
}
model{
    vector[N] mu;
    sigma ~ exponential( 1 );
    a ~ normal( 0 , 0.2 );
    bL ~ normal( 0 , 0.5 );
    mu = a + bL * L;
    K ~ normal( mu , sigma );
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m6_4s = SampleModel("m6.4s", stan6_4);
m6_4_data = Dict("N" => size(df, 1), "L" => df[:, :perc_lactose_s],
    "K" => df[!, :kcal_per_g_s]);
rc6_4s = stan_sample(m6_4s, data=m6_4_data);
if success(rc6_4s)
    part6_4s = read_samples(m6_4s, :particles)
    part6_4s |> display
end

# End of m6.4s.jl


text:
function _logsumexp(a, b)
    m = max(a, b)
    return m + log(exp(a - m) + exp(b - m))
end


text:
default_array(::Type{S}, d) where {S} = Array{S}(undef, d)

struct StructArrayInitializer{F, G}
    unwrap::F
    default_array::G
end
StructArrayInitializer(unwrap = t->false) = StructArrayInitializer(unwrap, default_array)

const default_initializer = StructArrayInitializer()

function (s::StructArrayInitializer)(S, d)
    ai = ArrayInitializer(s.unwrap, s.default_array)
    buildfromschema(typ -> ai(typ, d), S)
end

struct ArrayInitializer{F, G}
    unwrap::F
    default_array::G
end
ArrayInitializer(unwrap = t->false) = ArrayInitializer(unwrap, default_array)

(s::ArrayInitializer)(S, d) = s.unwrap(S) ? buildfromschema(typ -> s(typ, d), S) : s.default_array(S, d)

_reshape(v, itr) = _reshape(v, itr, Base.IteratorSize(itr))
_reshape(v, itr, ::Base.HasShape) = reshapestructarray(v, axes(itr))
_reshape(v, itr, ::Union{Base.HasLength, Base.SizeUnknown}) = v

# temporary workaround before it gets easier to support reshape with offset axis
reshapestructarray(v::AbstractArray, d) = reshape(v, d)
reshapestructarray(v::StructArray{T}, d) where {T} =
    StructArray{T}(map(x -> reshapestructarray(x, d), fieldarrays(v)))

"""
`collect_structarray(itr, fr=iterate(itr); initializer = default_initializer)`

Collects `itr` into a `StructArray`. The user can optionally pass a `initializer`, that is to say
a function `(S, d) -> v` that associates to a type and a size an array of eltype `S`
and size `d`. By default `initializer` returns a `StructArray` of `Array` but custom array types
may be used. `fr` represents the moment in the iteration of `itr` from which to start collecting.
"""
collect_structarray(itr; initializer = default_initializer) =
    collect_structarray(itr, iterate(itr); initializer = initializer)

collect_structarray(itr, fr; initializer = default_initializer) =
    collect_structarray(itr, fr, Base.IteratorSize(itr); initializer = initializer)

collect_structarray(itr, ::Nothing; initializer = default_initializer) =
    collect_empty_structarray(itr; initializer = initializer)

function collect_empty_structarray(itr::T; initializer = default_initializer) where {T}
    S = Core.Compiler.return_type(first, Tuple{T})
    res = initializer(S, (0,))
    _reshape(res, itr)
end

function collect_structarray(itr, elem, sz::Union{Base.HasShape, Base.HasLength};
                             initializer = default_initializer)
    el, i = elem
    S = typeof(el)
    dest = initializer(S, (length(itr),))
    dest[1] = el
    v = collect_to_structarray!(dest, itr, 2, i)
    _reshape(v, itr, sz)
end

function collect_to_structarray!(dest::AbstractArray, itr, offs, st)
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = offs
    while true
        elem = iterate(itr, st)
        elem === nothing && break
        el, st = elem
        if iscompatible(el, dest)
            @inbounds dest[i] = el
            i += 1
        else
            new = widenstructarray(dest, i, el)
            @inbounds new[i] = el
            return collect_to_structarray!(new, itr, i+1, st)
        end
    end
    return dest
end

function collect_structarray(itr, elem, ::Base.SizeUnknown; initializer = default_initializer)
    el, st = elem
    dest = initializer(typeof(el), (1,))
    dest[1] = el
    grow_to_structarray!(dest, itr, iterate(itr, st))
end

function grow_to_structarray!(dest::AbstractArray, itr, elem = iterate(itr))
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = length(dest)+1
    while elem !== nothing
        el, st = elem
        if iscompatible(el, dest)
            push!(dest, el)
            elem = iterate(itr, st)
            i += 1
        else
            new = widenstructarray(dest, i, el)
            push!(new, el)
            return grow_to_structarray!(new, itr, iterate(itr, st))
        end
    end
    return dest
end

widenstructarray(dest::AbstractArray{S}, i, el::T) where {S, T} = widenstructarray(dest, i, _promote_typejoin(S, T))

function widenstructarray(dest::StructArray, i, ::Type{T}) where {T}
    sch = hasfields(T) ? staticschema(T) : nothing
    sch !== nothing && fieldnames(sch) == propertynames(dest) || return widenarray(dest, i, T)
    types = ntuple(x -> fieldtype(sch, x), fieldcount(sch))
    cols = Tuple(fieldarrays(dest))
    newcols = map((a, b) -> widenstructarray(a, i, b), cols, types)
    return StructArray{T}(newcols)
end

widenstructarray(dest::AbstractArray, i, ::Type{T}) where {T} = widenarray(dest, i, T)

widenarray(dest::AbstractArray{T}, i, ::Type{T}) where {T} = dest
function widenarray(dest::AbstractArray, i, ::Type{T}) where T
    new = similar(dest, T, length(dest))
    copyto!(new, 1, dest, 1, i-1)
    new
end