def nmf_kfold(data, k, n_runs=10, **nmf_params):
    """
    Runs K-fold ensemble topic modeling (Belford et al. 2017)
    """
    # TODO
    nmf = NMF(k)
    W_list = []
    kf = KFold(n_splits=n_runs, shuffle=True)
    # TODO: randomly divide data into n_runs folds
    for train_index, test_index in kf.split(data.T):
        W = nmf.fit_transform(data[:,train_index])
        W_list.append(W)
    W_stacked = np.hstack(W_list)
    nmf_w = nmf.fit_transform(W_stacked)
    nmf_h = nmf.components_
    H_new = data.T.dot(nmf_w).T
    nmf2 = NMF(k, init='custom')
    nmf_w = nmf2.fit_transform(data, W=nmf_w, H=H_new)
    H_new = nmf2.components_
    #W_new = W_new/W_new.sum(0)
    return nmf_w, H_new