text:
"""
A structure that stores the discrete semantic vectors:
pS is the discrete semantic matrix;
f2i is a dictionary returning the indices for features;
i2f is a dictionary returning the features for indices.
"""
struct PS_Matrix_Struct
    pS::Union{Matrix,SparseMatrixCSC}
    f2i::Dict
    i2f::Dict
end

"""
A structure that stores Lexome semantic vectors:
L is Lexome semantic matrix;
f2i is a dictionary returning the indices for features;
i2f is a dictionary returning the features for indices.
"""
struct L_Matrix_Struct
    L::Matrix{Float64}
    f2i::Dict{AbstractString,Int64}
    i2f::Vector{AbstractString}
    ncol::Int64
end

"""
Make discrete semantic matrix.
"""
function make_pS_matrix end

"""
Make simulated semantic matrix.
"""
function make_S_matrix end

"""
Make simulated lexome matrix.
"""
function make_L_matrix end

"""
Make combined simulated S matrices, where combined features from both training datasets and validation datasets
"""
function make_combined_S_matrix end

"""
Make combined simulated Lexome matrix, where combined features from both training datasets and validation datasets
"""
function make_combined_L_matrix end

"""
    make_pS_matrix(utterances)

Create a discrete semantic matrix given a dataframe.

# Obligatory Arguments
- `data::DataFrame`: the dataset

# Optional Arguments
- `features_col::Symbol=:CommunicativeIntention`: the column name for target
- `sep_token::String="_"`: separator

# Examples
```julia
s_obj_train = JudiLing.make_pS_matrix(
    utterance,
    features_col=:CommunicativeIntention,
    sep_token="_")
```
"""
function make_pS_matrix(
    utterances;
    features_col = :CommunicativeIntention,
    sep_token = "_",
)

    # find out all possible features in this dataset
    features = unique(vcat(split.(utterances[:, features_col], sep_token)...))

    # using dict to store feature names
    f2i = Dict(v => i for (i, v) in enumerate(features))
    i2f = Dict(i => v for (i, v) in enumerate(features))

    # find out features for each utterance
    vs = unique.(split.(utterances[:, features_col], sep_token))

    # total number of feature in the entire dataset
    # to initialize a sparse matrix
    n_f = sum([length(v) for v in vs])

    # initialize sparse matrix components
    m = size(utterances, 1)
    n = length(i2f)
    I = zeros(Int64, n_f)
    J = zeros(Int64, n_f)
    V = ones(Int64, n_f)

    # calculate each cell in sparse matrix
    cnt = 0
    for (i, v) in enumerate(vs)
        for (j, f) in enumerate(v)
            cnt += 1
            I[cnt] = i
            J[cnt] = f2i[f]
        end
    end

    # create sparse matrix
    pS = sparse(I, J, V, m, n, *)

    PS_Matrix_Struct(pS, f2i, i2f)
end

"""
    make_pS_matrix(utterances, utterances_train)

Construct discrete semantic matrix for the validation datasets given by the
exemplar in the dataframe, and given the S matrix for the training datasets.

# Obligatory Arguments
- `utterances::DataFrame`: the dataset
- `utterances_train::PS_Matrix_Struct`: training PS object

# Optional Arguments
- `features_col::Symbol=:CommunicativeIntention`: the column name for target
- `sep_token::String="_"`: separator

# Examples
```julia
s_obj_val = JudiLing.make_pS_matrix(
    utterance_val,
    s_obj_train,
    features_col=:CommunicativeIntention,
    sep_token="_")
```
"""
function make_pS_matrix(
    utterances,
    utterances_train;
    features_col = :CommunicativeIntention,
    sep_token = "_",
)

    # find out all possible features in this dataset
    features = unique(vcat(split.(utterances[:, features_col], sep_token)...))

    # using dict to store feature names
    f2i = utterances_train.f2i
    i2f = utterances_train.i2f

    # find out features for each utterance
    vs = unique.(split.(utterances[:, features_col], sep_token))

    # total number of feature in the entire dataset
    # to initialize a sparse matrix
    n_f = sum([length(v) for v in vs])

    # initialize sparse matrix components
    m = size(utterances, 1)
    n = length(i2f)
    I = zeros(Int64, n_f)
    J = zeros(Int64, n_f)
    V = ones(Int64, n_f)

    # calculate each cell in sparse matrix
    cnt = 0
    for (i, v) in enumerate(vs)
        for (j, f) in enumerate(v)
            cnt += 1
            I[cnt] = i
            J[cnt] = f2i[f]
        end
    end

    # create sparse matrix
    pS = sparse(I, J, V, m, n, *)

    PS_Matrix_Struct(pS, f2i, i2f)
end

"""
    make_S_matrix(data::DataFrame, base::Vector, inflections::Vector)

Create simulated semantic matrix for the training datasets, given the input
data of a vector specified contex lexemes and a vector specified gramatic
lexemes. The semantic vector of a word form is constructed summing semantic
vectors of content and gramatic lexemes.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train = JudiLing.make_S_matrix(
    french,
    ["Lexeme"],
    ["Tense","Aspect","Person","Number","Gender","Class","Mood"],
    ncol=200)

# deep mode
S_train = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    isdeep=true,
    ...)

# non-deep mode
S_train = JudiLing.make_S_matrix(
    ...
    isdeep=false,
    ...)

# add additional Gaussian noise
S_train = JudiLing.make_S_matrix(
    ...
    add_noise=true,
    sd_noise=1,
    ...)

# further control of means and standard deviations
S_train = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    sd_base=4,
    sd_inflection=4,
    sd_noise=1,
    ...)
```
"""
function make_S_matrix(
    data::DataFrame,
    base::Vector,
    inflections::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_L_matrix(
        data,
        base,
        inflections,
        sd_base_mean = sd_base_mean,
        sd_inflection_mean = sd_inflection_mean,
        sd_base = sd_base,
        sd_inflection = sd_inflection,
        ncol = ncol,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data,
        base,
        inflections,
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_S_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector, inflections::Vector)

Create simulated semantic matrix for the validation datasets, given the input
data of a vector specified contex lexemes and a vector specified gramatic
lexemes. The semantic vector of a word form is constructed summing semantic
vectors of content and gramatic lexemes.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_S_matrix(
    french,
    french_val,
    ["Lexeme"],
    ["Tense","Aspect","Person","Number","Gender","Class","Mood"],
    ncol=200)

# deep mode
S_train, S_val = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    isdeep=true,
  ...)

# non-deep mode
S_train, S_val = JudiLing.make_S_matrix(
    ...
    isdeep=false,
    ...)

# add additional Gaussian noise
S_train, S_val = JudiLing.make_S_matrix(
    ...
    add_noise=true,
    sd_noise=1,
    ...)

# further control of means and standard deviations
S_train, S_val = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    sd_base=4,
    sd_inflection=4,
    sd_noise=1,
    ...)
```
"""
function make_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector,
    inflections::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_L_matrix(
        data_train,
        base,
        inflections,
        sd_base_mean = sd_base_mean,
        sd_inflection_mean = sd_inflection_mean,
        sd_base = sd_base,
        sd_inflection = sd_inflection,
        ncol = ncol,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data_train,
        data_val,
        base,
        inflections,
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_S_matrix(data::DataFrame, base::Vector)

Create simulated semantic matrix for the training datasets with only base
features, given the input data of a vector specified contex lexemes and a
vector specified gramatic lexemes. The semantic vector of a word form is
constructed summing semantic vectors of content and gramatic lexemes.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_base::Int64=4`: the sd of base features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train = JudiLing.make_S_matrix(
    french,
    ["Lexeme"],
    ncol=200)

# deep mode
S_train = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    isdeep=true,
    ...)

# non-deep mode
S_train = JudiLing.make_S_matrix(
    ...
    isdeep=false,
    ...)

# add additional Gaussian noise
S_train = JudiLing.make_S_matrix(
    ...
    add_noise=true,
    sd_noise=1,
    ...)

# further control of means and standard deviations
S_train = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    sd_base=4,
    sd_inflection=4,
    sd_noise=1,
    ...)
```
"""
function make_S_matrix(
    data::DataFrame,
    base::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_base = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_L_matrix(
        data,
        base,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        ncol = ncol,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data,
        base,
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_S_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector)

Create simulated semantic matrix for the validation datasets with only base
features, given the input data of a vector specified contex lexemes and a
vector specified gramatic lexemes. The semantic vector of a word form is
constructed summing semantic vectors of content and gramatic lexemes.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_base::Int64=4`: the sd of base features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_S_matrix(
    french,
    french_val,
    ["Lexeme"],
    ncol=200)

# deep mode
S_train, S_val = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    isdeep=true,
    ...)

# non-deep mode
S_train, S_val = JudiLing.make_S_matrix(
    ...
    isdeep=false,
    ...)

# add additional Gaussian noise
S_train, S_val = JudiLing.make_S_matrix(
    ...
    add_noise=true,
    sd_noise=1,
    ...)

# further control of means and standard deviations
S_train, S_val = JudiLing.make_S_matrix(
    ...
    sd_base_mean=1,
    sd_inflection_mean=1,
    sd_base=4,
    sd_inflection=4,
    sd_noise=1,
    ...)
```
"""
function make_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_base = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_L_matrix(
        data_train,
        base,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        ncol = ncol,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data_train,
        data_val,
        base,
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_S_matrix(data_train::DataFrame, data_val::DataFrame, pyndl_weights::Pyndl_Weight_Struct, n_features_columns::Vector)

Create semantic matrix for pyndl mode
"""
function make_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    pyndl_weights::Pyndl_Weight_Struct,
    n_features_columns::Vector,
)

    f2i = Dict(v => i for (i, v) in enumerate(pyndl_weights.outcomes))
    i2f = Dict(i => v for (i, v) in enumerate(pyndl_weights.outcomes))

    n_f = length(pyndl_weights.outcomes)

    St_train = zeros(Float64, n_f, size(data_train, 1))
    for i = 1:size(data_train, 1)
        for f in data_train[i, n_features_columns]
            St_train[f2i[f], i] = 1
        end
    end

    St_val = zeros(Float64, n_f, size(data_val, 1))
    for i = 1:size(data_val, 1)
        for f in data_val[i, n_features_columns]
            St_val[f2i[f], i] = 1
        end
    end

    St_train', St_val'
end

"""
    make_S_matrix(data_train::DataFrame, base::Vector, inflections::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix where lexome matrix is available.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes
- `L::L_Matrix_Struct`: the lexome matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S1 = JudiLing.make_S_matrix(
    latin,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
     L1,
     add_noise=true,
    sd_noise=1,
    normalized=false
    )
```
"""
function make_S_matrix(
    data_train::DataFrame,
    base::Vector,
    inflections::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    n_train = size(data_train, 1)
    n_base = length(base)
    n_infl = length(inflections)

    St_train = make_St(L, n_train, data_train, base, inflections)
    add_noise && add_St_noise!(St_train, sd_noise)
    normalized && normalize_St!(St_train, n_base, n_infl)

    Array(St_train')
end

"""
    make_S_matrix(data_train::DataFrame, data_val::Union{DataFrame, Nothing}, base::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix where lexome matrix is available.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `L::L_Matrix_Struct`: the lexome matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S1, S2 = JudiLing.make_S_matrix(
     latin,
    latin_val,
    ["Lexeme"],
    L1,
    add_noise=true,
    sd_noise=1,
    normalized=false
    )
```
"""
function make_S_matrix(
    data_train::DataFrame,
    data_val::Union{DataFrame,Nothing},
    base::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    JudiLing.make_S_matrix(
        data_train,
        data_val,
        base,
        [],
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )

end

"""
    make_S_matrix(data::DataFrame, base::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix where lexome matrix is available.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes
- `L::L_Matrix_Struct`: the lexome matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S1 = JudiLing.make_S_matrix(
    latin,
    ["Lexeme"],
    L1,
    add_noise=true,
    sd_noise=1,
    normalized=false
    )
```
"""
function make_S_matrix(
    data::DataFrame,
    base::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    JudiLing.make_S_matrix(
        data,
        base,
        [],
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )

end

"""
    make_S_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector, inflections::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix where lexome matrix is available.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes
- `L::L_Matrix_Struct`: the lexome matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S1, S2 = JudiLing.make_S_matrix(
    latin,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    L1,
    add_noise=true,
    sd_noise=1,
    normalized=false
    )
```
"""
function make_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector,
    inflections::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    n_train = size(data_train, 1)
    n_val = size(data_val, 1)
    n_base = length(base)
    n_infl = length(inflections)

    St_train = make_St(L, n_train, data_train, base, inflections)
    add_noise && add_St_noise!(St_train, sd_noise)
    normalized && normalize_St!(St_train, n_base, n_infl)

    St_val = make_St(L, n_val, data_val, base, inflections)
    add_noise && add_St_noise!(St_val, sd_noise)
    normalized && normalize_St!(St_val, n_base, n_infl)

    Array(St_train'), Array(St_val')
end

"""
    make_L_matrix(data::DataFrame, base::Vector, inflections::Vector)

Create Lexome Matrix with simulated semantic vectors.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized

# Examples
```julia
# basic usage
L = JudiLing.make_L_matrix(
    latin,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    ncol=200)
```
"""
function make_L_matrix(
    data::DataFrame,
    base::Vector,
    inflections::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
)

    # collect all features and f2i mappings
    base_f, base_f2i = process_features(data, base)
    infl_f, infl_f2i = process_features(data, inflections)

    # setup seed
    Random.seed!(seed)

    n_base_f = length(base_f)
    n_infl_f = length(infl_f)

    # pre-allocated L matrix
    L = Matrix{Float64}(undef, (n_base_f + n_infl_f, ncol))

    if isdeep # deep mode random means for each feature
        return L_Matrix_Struct(
            L,
            sd_base,
            sd_base_mean,
            sd_inflection,
            sd_inflection_mean,
            base_f,
            infl_f,
            base_f2i,
            infl_f2i,
            n_base_f,
            n_infl_f,
            ncol,
        )
    else # otherwise use mean=0 for all features
        return L_Matrix_Struct(
            L,
            sd_base,
            sd_inflection,
            base_f,
            infl_f,
            base_f2i,
            infl_f2i,
            n_base_f,
            n_infl_f,
            ncol,
        )
    end
end

"""
    make_L_matrix(data::DataFrame, base::Vector)

Create Lexome Matrix with simulated semantic vectors where there are only base features.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `base::Vector`: context lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_base::Int64=4`: the sd of base features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized

# Examples
```julia
# basic usage
L = JudiLing.make_L_matrix(
    latin,
    ["Lexeme"],
    ncol=200)
```
"""
function make_L_matrix(
    data::DataFrame,
    base::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_base = 4,
    seed = 314,
    isdeep = true,
)

    make_L_matrix(
        data,
        base,
        [];
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        seed = seed,
        isdeep = isdeep,
    )
end

"""
    make_combined_L_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector, inflections::Vector)

Create Lexome Matrix with simulated semantic vectors, where features are
combined from both training datasets and validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized

# Examples
```julia
# basic usage
L = JudiLing.make_combined_L_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    ncol=n_features)
```
"""
function make_combined_L_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector,
    inflections::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
)

    data_combined = copy(data_train)
    append!(data_combined, data_val)

    make_L_matrix(
        data_combined,
        base,
        inflections;
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        seed = seed,
        isdeep = isdeep,
    )

end

"""
    make_combined_L_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector)

Create Lexome Matrix with simulated semantic vectors, where features are
combined from both training datasets and validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized

# Examples
```julia
# basic usage
L = JudiLing.make_combined_L_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ncol=n_features)
```
"""
function make_combined_L_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
)

    make_combined_L_matrix(
        data_train,
        data_val,
        base,
        [];
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        seed = seed,
        isdeep = isdeep,
    )
end

"""
    make_combined_S_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector, inflections::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix for the training datasets and validation
datasets with existing Lexome matrix, where features are combined from both
training datasets and validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes
- `L::L_Matrix_Struct`: the Lexome Matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_combined_S_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    L)
```
"""
function make_combined_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector,
    inflections::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    make_S_matrix(
        data_train,
        data_val,
        base,
        inflections,
        L;
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_combined_S_matrix(data_train::DataFrame, data_val::Union{DataFrame, Nothing}, base::Vector, L::L_Matrix_Struct)

Create simulated semantic matrix for the training datasets and validation
datasets with existing Lexome matrix, where features are combined from both
training datasets and validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `L::L_Matrix_Struct`: the Lexome Matrix

# Optional Arguments
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_combined_S_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    L)
```
"""
function make_combined_S_matrix(
    data_train::DataFrame,
    data_val::Union{DataFrame,Nothing},
    base::Vector,
    L::L_Matrix_Struct;
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    make_S_matrix(
        data_train,
        data_val,
        base,
        L;
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_combined_S_matrix(  data_train::DataFrame, data_val::DataFrame, base::Vector, inflections::Vector)

Create simulated semantic matrix for the training datasets and validation
datasets, where features are combined from both training datasets and
validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes
- `inflections::Vector`: grammatic lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_combined_S_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    ncol=n_features)
```
"""
function make_combined_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector,
    inflections::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_combined_L_matrix(
        data_train,
        data_val,
        base,
        inflections;
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data_train,
        data_val,
        base,
        inflections,
        L;
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    make_combined_S_matrix(data_train::DataFrame, data_val::DataFrame, base::Vector)

Create simulated semantic matrix for the training datasets and validation
datasets, where features are combined from both training datasets and
validation datasets.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `base::Vector`: context lexemes

# Optional Arguments
- `ncol::Int64=200`: dimension of semantic vectors, usually the same as that of cue vectors
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `seed::Int64=314`: the random seed
- `isdeep::Bool=true`: if true, mean of each feature is also randomized
- `add_noise::Bool=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Bool=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd

# Examples
```julia
# basic usage
S_train, S_val = JudiLing.make_combined_S_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person","Number","Tense","Voice","Mood"],
    ncol=n_features)
```
"""
function make_combined_S_matrix(
    data_train::DataFrame,
    data_val::DataFrame,
    base::Vector;
    ncol = 200,
    sd_base_mean = 1,
    sd_inflection_mean = 1,
    sd_base = 4,
    sd_inflection = 4,
    seed = 314,
    isdeep = true,
    add_noise = true,
    sd_noise = 1,
    normalized = false,
)

    L = make_combined_L_matrix(
        data_train,
        data_val,
        base,
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_base = sd_base,
        seed = seed,
        isdeep = isdeep,
    )

    make_S_matrix(
        data_train,
        data_val,
        base,
        L,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized,
    )
end

"""
    L_Matrix_Struct(L, sd_base, sd_base_mean, sd_inflection, sd_inflection_mean, base_f, infl_f, base_f2i, infl_f2i, n_base_f, n_infl_f, ncol)

Construct L_Matrix_Struct with deep mode.
"""
function L_Matrix_Struct(
    L,
    sd_base,
    sd_base_mean,
    sd_inflection,
    sd_inflection_mean,
    base_f,
    infl_f,
    base_f2i,
    infl_f2i,
    n_base_f,
    n_infl_f,
    ncol,
)
    comp_f_M!(L, sd_base, sd_base_mean, n_base_f, ncol, 0)
    comp_f_M!(L, sd_inflection, sd_inflection_mean, n_infl_f, ncol, n_base_f)

    L_Matrix_Struct(
        L,
        merge_f2i(base_f2i, infl_f2i, n_base_f, n_infl_f),
        vcat(base_f, infl_f),
        ncol,
    )
end

"""
    L_Matrix_Struct(L, sd_base, sd_inflection, base_f, infl_f, base_f2i, infl_f2i, n_base_f, n_infl_f, ncol)

Construct L_Matrix_Struct without deep mode.
"""
function L_Matrix_Struct(
    L,
    sd_base,
    sd_inflection,
    base_f,
    infl_f,
    base_f2i,
    infl_f2i,
    n_base_f,
    n_infl_f,
    ncol,
)
    comp_f_M!(L, sd_base, n_base_f, ncol, 0)
    comp_f_M!(L, sd_inflection, n_infl_f, ncol, n_base_f)

    L_Matrix_Struct(
        L,
        merge_f2i(base_f2i, infl_f2i, n_base_f, n_infl_f),
        vcat(base_f, infl_f),
        ncol,
    )
end

"""
    process_features(data, feature_cols)

Collect all features given datasets and feature column names.
"""
function process_features(data, feature_cols)
    features =
        [f for fc in feature_cols for f in skipmissing(unique(data[:, fc]))]
    base_f2i = Dict(v => i for (i, v) in enumerate(features))

    features, base_f2i
end

"""
    comp_f_M!(L, sd, sd_mean, n_f, ncol, n_b)

Compose feature Matrix with deep mode.
"""
function comp_f_M!(L, sd, sd_mean, n_f, ncol, n_b)
    means = rand(Normal(0, sd_mean), n_f)
    for i = 1:n_f
        L[i+n_b, :] = rand(Normal(means[i], sd), ncol)
    end
end

"""
    comp_f_M!(L, sd, n_f, ncol, n_b)

Compose feature Matrix without deep mode.
"""
function comp_f_M!(L, sd, n_f, ncol, n_b)
    for i = 1:n_f
        L[i+n_b, :] = rand(Normal(0, sd), ncol)
    end
end

"""
    merge_f2i(base_f2i, infl_f2i, n_base_f, n_infl_f)

Merge base f2i dictionary and inflectional f2i dictionary.
"""
function merge_f2i(base_f2i, infl_f2i, n_base_f, n_infl_f)
    # add infl_f2i into base_f2i
    # infl_f2i need shift n_base_f
    f2i = copy(base_f2i)
    for (k, v) in infl_f2i
        f2i[k] = v + n_base_f
    end
    f2i
end

"""
    lexome_sum(L, features)

Sum up semantic vector, given lexome vector.
"""
function lexome_sum(L, features)
    ls = [L.L[L.f2i[f], :] for f in skipmissing(features)]
    if length(ls) > 0
        return sum(ls)
    end
    zeros(Float64, L.ncol)
end

"""
    make_St(L, n, data, base, inflections)

Make S transpose matrix with inflections.
"""
function make_St(L, n, data, base, inflections)
    St = Array{Float64,2}(undef, L.ncol, n)
    for i = 1:n
        ls_base = lexome_sum(L, data[i, base])
        ls_infl = lexome_sum(L, data[i, inflections])
        ls = ls_base + ls_infl
        St[:, i] = ls
    end
    St
end

"""
    make_St(L, n, data, base)

Make S transpose matrix without inflections.
"""
function make_St(L, n, data, base)
    St = Array{Float64,2}(undef, L.ncol, n)
    for i = 1:n
        ls_base = lexome_sum(L, data[i, base])
        St[:, i] = ls_base
    end
    St
end

"""
    add_St_noise!(St, sd_noise)

Add noise.
"""
function add_St_noise!(St, sd_noise)
    Noise = rand(Normal(0, sd_noise), size(St))
    broadcast!(+, St, St, Noise)
end

"""
    normalize_St!(St, n_base, n_infl)

Normalize S transpose with inflections.
"""
function normalize_St!(St, n_base, n_infl)
    broadcast!(/, St, St, n_base + n_infl)
end

"""
    normalize_St!(St, n_base)

Normalize S transpose without inflections.
"""
function normalize_St!(St, n_base)
    St = St ./ n_base
end


text:
using SymbolicUtils: Symbolic, Sym, FnType, Term, Add, Mul, Pow, symtype, operation, arguments
using SymbolicUtils
using IfElse: ifelse
using Test

@testset "@syms" begin
    let
        @syms a b::Float64 f(::Real) g(p, h(q::Real))::Int

        @test a isa Sym{Number}
        @test a.name === :a

        @test b isa Sym{Float64}
        @test b.name === :b

        @test f isa Sym{FnType{Tuple{Real}, Number}}
        @test f.name === :f

        @test g isa Sym{FnType{Tuple{Number, FnType{Tuple{Real}, Number}}, Int}}
        @test g.name === :g

        @test f(b) isa Term
        @test symtype(f(b)) === Number
        @test_throws ErrorException f(a)

        @test g(b, f) isa Term
        @test_throws ErrorException g(b, a)

        @test symtype(g(b, f)) === Int

        # issue #91
        @syms h(a,b,c)
        @test isequal(h(1,2,3), h(1,2,3))
    end
end

@testset "hashing" begin
    @syms a b f(x, y)
    @test hash(a) == hash(a)
    @test hash(a) != hash(b)
    @test hash(a+1) == hash(a+1)
    @test hash(sin(a+1)) == hash(sin(a+1))
    @test hash(f(1,a)) == hash(f(1, a))

    c = a
    g = f
    @syms a f(x, y)
    @test hash(a) == hash(c)
    @test hash(g(a, b)) == hash(f(a,b))
    @test hash(f(a, b)) == hash(f(c,b))
    @test hash(sin(a+1)) == hash(sin(c+1))

    ex = sin(a+1)
    h = hash(ex, UInt(0))
    @test ex.hash[] == h
    ex1 = sin(a+1)
    hash(asin(ex1), UInt(0))
    @test ex1.hash[] == h
end

struct Ctx1 end
struct Ctx2 end

@testset "metadata" begin
    @syms a b c
    for a = [a, sin(a), a+b, a*b, a^3]

        a′ = setmetadata(a, Ctx1, "meta_1")

        @test hasmetadata(a′, Ctx1)
        @test !hasmetadata(a′, Ctx2)

        a′ = setmetadata(a′, Ctx2, "meta_2")

        @test hasmetadata(a′, Ctx1)
        @test hasmetadata(a′, Ctx2)

        @test getmetadata(a′, Ctx1) == "meta_1"
        @test getmetadata(a′, Ctx2) == "meta_2"
    end

    # In substitute #283
    #
    @syms f(t) t
    f = setmetadata(f(t), Ctx1, "yes")
    hasmetadata(f, Ctx1) # true
    newf = substitute(f, Dict(a=>b)) # unrelated substitution
    @test hasmetadata(newf, Ctx1)
    @test getmetadata(newf, Ctx1) == "yes"

    @test isequal(substitute(1+sqrt(a), Dict(a => 2), fold=false),
                  1 + term(sqrt, 2, type=Number))
    @test substitute(1+sqrt(a), Dict(a => 2), fold=true) isa Float64
end

@testset "Base methods" begin
    @syms w::Complex z::Complex a::Real b::Real x

    @test isequal(w + z, Add(Number, 0, Dict(w=>1, z=>1)))
    @test isequal(z + a, Add(Number, 0, Dict(z=>1, a=>1)))
    @test isequal(a + b, Add(Real, 0, Dict(a=>1, b=>1)))
    @test isequal(a + x, Add(Number, 0, Dict(a=>1, x=>1)))
    @test isequal(a + z, Add(Number, 0, Dict(a=>1, z=>1)))

    foo(w, z, a, b) = 1.0
    SymbolicUtils.promote_symtype(::typeof(foo), args...) = Real
    @test SymbolicUtils._promote_symtype(foo, (w, z, a, b,)) === Real

    # promote_symtype of identity
    @test isequal(Term(identity, [w]), Term{Complex}(identity, [w]))
    @test isequal(+(w), w)
    @test isequal(+(a), a)

    @test isequal(rem2pi(a, RoundNearest), Term{Real}(rem2pi, [a, RoundNearest]))

    # bool
    for f in [(==), (!=), (<=), (>=), (<), (>)]
        @test isequal(f(a, 0), Term{Bool}(f, [a, 0]))
        @test isequal(f(0, a), Term{Bool}(f, [0, a]))
        @test isequal(f(a, a), Term{Bool}(f, [a, a]))
    end

    @test symtype(ifelse(true, 4, 5)) == Int
    @test symtype(ifelse(a < 0, b, w)) == Union{Real, Complex}
    @test SymbolicUtils.promote_symtype(ifelse, Bool, Int, Bool) == Union{Int, Bool}
    @test_throws MethodError w < 0
    @test isequal(w == 0, Term{Bool}(==, [w, 0]))

    @eqtest x // 5 == (1 // 5) * x
    @eqtest (1//2 * x) / 5 == (1 // 10) * x
    @eqtest x // Int16(5) == Rational{Int16}(1, 5) * x
    @eqtest 5 // x == 5 / x
    @eqtest x // a == x / a
end

@testset "array-like operations" begin
    abstract type SquareDummy end
    Base.:*(a::Symbolic{SquareDummy}, b) = b^2
    @syms s t a::SquareDummy A[1:2, 1:2]

    @test isequal(ndims(A), 2)
    @test_broken isequal(a.*[1 (s+t); t pi], [1 (s+t)^2; t^2 pi^2])
    @test isequal(s.*[1 (s+t); t pi], [s s*(s+t); s*t s*pi])
end

@testset "err test" begin
    @syms t()
    @test_throws ErrorException t(2)
end

@testset "substitute" begin
    @syms a b
    @test substitute(a, Dict(a=>1)) == 1
    @test isequal(substitute(sin(a+b), Dict(a=>1)), sin(b+1))
    @test substitute(a+b, Dict(a=>1, b=>3)) == 4
    @test substitute(exp(a), Dict(a=>2)) ≈ exp(2)
end

@testset "occursin" begin
    @syms a b c
    @test occursin(a, a + b)
    @test !occursin(sin(a), a + b + c)
    @test occursin(sin(a),  a * b + c + sin(a^2 * sin(a)))
end

@testset "printing" begin
    @syms a b c
    @test repr(a+b) == "a + b"
    @test repr(-a) == "-a"
    @test repr(-a + 3) == "3 - a"
    @test repr(-(a + b)) == "-a - b"
    @test repr((2a)^(-2a)) == "(2a)^(-2a)"
    @test repr(1/2a) == "1 / (2a)"
    @test repr(2/(2*a)) == "1 / a"
    @test repr(Term(*, [1, 1])) == "1"
    @test repr(Term(*, [2, 1])) == "2*1"
    @test repr((a + b) - (b + c)) == "a - c"
    @test repr(a + -1*(b + c)) == "a - b - c"
    @test repr(a + -1*b) == "a - b"
    @test repr(-1^a) == "-(1^a)"
    @test repr((-1)^a) == "(-1)^a"
end

@testset "similarterm" begin
    @syms a b c
    @test isequal(SymbolicUtils.similarterm((b + c), +, [a,  (b+c)]).dict, Dict(a=>1,b=>1,c=>1))
    @test isequal(SymbolicUtils.similarterm(b^2, ^, [b^2,  1//2]), b)
end

toterm(t) = Term{symtype(t)}(operation(t), arguments(t))

@testset "diffs" begin
    @syms a b c
    @test isequal(toterm(-1c), Term{Number}(*, [-1, c]))
    @test isequal(toterm(-1(a+b)), Term{Number}(+, [-1a, -b]))
    @test isequal(toterm((a + b) - (b + c)), Term{Number}(+, [a, -1c]))
end

@testset "hash" begin
    @syms a b
    @test hash(a + b, UInt(0)) === hash(a + b) === hash(a + b, UInt(0)) # test caching
    @test hash(a + b, UInt(2)) !== hash(a + b)
end

@testset "methoderror" begin
    @syms a::Any b::Any

    @test_throws MethodError a * b
    @test_throws MethodError a + b
end

@testset "canonical form" begin
    @syms a b c
    for x in [a, a*b, a+b, a-b, a^2, sin(a)]
        @test isequal(x * 1, x)
        @test x * 0 === 0
        @test isequal(x + 0, x)
        @test isequal(x + x, 2x)
        @test isequal(x + 2x, 3x)
        @test x - x === 0
        @test isequal(-x, -1x)
        @test isequal(x^1, x)
        @test isequal((x^-1)*inv(x^-1), 1)
    end
end

@testset "isequal" begin
    @syms a b c
    @test isequal(a + b, a + b + 0.01 - 0.01)
end

@testset "subtyping" begin
    T = FnType{Tuple{T,S,Int} where {T,S}, Real}
    s = Sym{T}(:t)
    @syms a b c::Int
    @test isequal(arguments(s(a, b, c)), [a, b, c])
end

@testset "div" begin
    @syms x::SafeReal y::Real
    @test (2x/2y).num isa Sym
    @test (2x/3y).num.coeff == 2
    @test (2x/3y).den.coeff == 3
    @test (2x/-3x).num.coeff == -2
    @test (2x/-3x).den.coeff == 3
    @test (2.5x/3x).num.coeff == 2.5
    @test (2.5x/3x).den.coeff == 3
    @test (x/3x).den.coeff == 3

    @syms x y
    @test (2x/2y).num isa Sym
    @test (2x/3y).num.coeff == 2
    @test (2x/3y).den.coeff == 3
    @test (2x/-3x) == -2//3
    @test (2.5x/3x).num == 2.5
    @test (2.5x/3x).den == 3
    @test (x/3x) == 1//3
end

@testset "LiteralReal" begin
    @syms x::LiteralReal y::LiteralReal z::LiteralReal
    @test repr(x+x) == "x + x"
    @test repr(x*x) == "x*x"
    @test repr(x*x + x*x) == "x*x + x*x"
    for ex in [sin(x), x+x, x*x, x\x, x/x]
        @test typeof(sin(x)) <: Term{LiteralReal}
    end
    @test repr(sin(x) + sin(x)) == "sin(x) + sin(x)"
end

ext: jl


text:
# Based on Lars Ruthotto's initial implementation.
function get_div_grad(n1 :: Int, n2 :: Int, n3 :: Int)

  # Divergence
  D1 = kron(Matrix(1.0I, n3, n3), kron(Matrix(1.0I, n2, n2), ddx(n1)))
  D2 = kron(Matrix(1.0I, n3, n3), kron(ddx(n2),   Matrix(1.0I, n1, n1)))
  D3 = kron(ddx(n3),   kron(Matrix(1.0I, n2, n2), Matrix(1.0I, n1, n1)))

  # DIV from faces to cell-centers
  Div = [D1 D2 D3]

  return Div * Div';
end

# 1D finite difference on staggered grid
function ddx(n :: Int)
  e = ones(n);
  return d = sparse([1:n; 1:n], [1:n; 2:n+1], [-e; e])
end