diff --git a/.gitignore b/.gitignore index 1427f59..6ca2367 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ inst/doc docs __pycache__ +smoke_test* diff --git a/DESCRIPTION b/DESCRIPTION index b0c0044..90d5c55 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,8 +24,7 @@ Depends: R (>= 3.5.0) Imports: Rcpp, generics, - matrixStats, - pracma + matrixStats Suggests: keras, tensorflow, reticulate, luz, torch, @@ -37,7 +36,7 @@ LinkingTo: RcppArmadillo VignetteBuilder: knitr Roxygen: list(markdown = TRUE) +Config/roxygen2/version: 8.0.0 Config/testthat/edition: 3 URL: https://ibidat.github.io/nn2poly/, https://github.com/IBiDat/nn2poly BugReports: https://github.com/IBiDat/nn2poly/issues -Config/roxygen2/version: 8.0.0 diff --git a/R/RcppExports.R b/R/RcppExports.R index 65354af..8bc3178 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,14 +1,30 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -select_allowed_partitions <- function(equivalent_label, q_previous_layer, labels, partitions) { - .Call(`_nn2poly_select_allowed_partitions`, equivalent_label, q_previous_layer, labels, partitions) +obtain_final_poly_order <- function(max_order, taylor_orders) { + .Call(`_nn2poly_obtain_final_poly_order`, max_order, taylor_orders) +} + +obtain_taylor_vector <- function(taylor_orders, af_string_list) { + .Call(`_nn2poly_obtain_taylor_vector`, taylor_orders, af_string_list) +} + +obtain_derivatives_list <- function(taylor_orders, af_string_list) { + .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list) +} + +obtain_partitions_with_labels <- function(p, q_max) { + .Call(`_nn2poly_obtain_partitions_with_labels`, p, q_max) } alg_non_linear <- function(coeffs_input, labels_input, labels_output, taylor_orders, current_layer, g, partitions_labels, partitions) { .Call(`_nn2poly_alg_non_linear`, coeffs_input, labels_input, labels_output, taylor_orders, current_layer, g, partitions_labels, partitions) } +nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders) { + .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders) +} + combinations_with_repetition <- function(n, k) { .Call(`_nn2poly_combinations_with_repetition`, n, k) } @@ -17,7 +33,3 @@ generate_partitions <- function(p, q_max) { .Call(`_nn2poly_generate_partitions`, p, q_max) } -generate_partitions_full <- function(p, q_max) { - .Call(`_nn2poly_generate_partitions_full`, p, q_max) -} - diff --git a/R/nn2poly.R b/R/nn2poly.R index d0c9d6a..b92a54e 100644 --- a/R/nn2poly.R +++ b/R/nn2poly.R @@ -47,13 +47,6 @@ NULL #' #' @param ... Ignored. #' -#' @param all_partitions Optional argument containing the needed multipartitions -#' as list of lists of lists. If set to \code{NULL}, nn2poly will compute said -#' multipartitions. This step can be computationally expensive when the chosen -#' polynomial order or the dimension are too high. In such cases, it is -#' encouraged that the multipartitions are stored and reused when possible. -#' Default set to \code{NULL}. -#' #' @return Returns an object of class `nn2poly`. #' #' If \code{keep_layers = FALSE} (default case), it returns a list with two @@ -121,49 +114,18 @@ nn2poly <- function(object, max_order = 2, keep_layers = FALSE, taylor_orders = 8, - ..., - all_partitions = NULL - ) { + ...) { UseMethod("nn2poly") } #' @export -nn2poly.list <- function(object, ...) { - result_raw <- nn2poly_algorithm(object, names(object), ...) - - # Check if object is a single polynomial or a list of polynomials. - # This is controlled by keep_layers. - bool_layers = is.null(result_raw$labels) - - if(bool_layers){ - n_items <- length(result_raw) - n_layers <- ceiling(n_items/2) - result <- list() - for (i in 1:n_layers){ - layer_name <- paste0("layer_", i) - result[[layer_name]] <- list() - result[[layer_name]][["input"]] <- result_raw[[(2*i-1)]] - # Transpose to have polynomials as columns - result[[layer_name]][["input"]][["values"]] <- - t(result[[layer_name]][["input"]][["values"]]) - if (2*i <= n_items){ - result[[layer_name]][["output"]] <- result_raw[[(2*i)]] - # Transpose to have polynomials as columns - result[[layer_name]][["output"]][["values"]] <- - t(result[[layer_name]][["output"]][["values"]]) - } else { - # If there is a linear output, i.e. single polynomial in final layer - # and odd number of items, then we repeat the input as the output, as - # the activation functions takes no effect on the polynomial - result[[layer_name]][["output"]] <- result[[layer_name]][["input"]] - } - } - - } else { - result_raw$values <- t(result_raw$values) - result = result_raw - } - +nn2poly.list <- function(object, + max_order = 2, + keep_layers = FALSE, + taylor_orders = 8, + ...) { + result <- nn2poly_algorithm(object, names(object), + max_order, keep_layers, taylor_orders) class(result) <- "nn2poly" result } diff --git a/R/nn2poly_algorithm.R b/R/nn2poly_algorithm.R deleted file mode 100644 index 12de066..0000000 --- a/R/nn2poly_algorithm.R +++ /dev/null @@ -1,399 +0,0 @@ -#' NN2Poly algorithm full algorithm. -#' -#' Internal function that performs the full NN2Poly algorithm given a list of -#' weights and a list of activation functions (af), -#' used inside nn2poly S3 method. -#' -#' @inheritParams nn2poly -#' -#' @param weights_list \code{list} of length L (number of hidden layers + 1) -#' containing the weights matrix for each layer. See nn2poly object argument -#' documentation. -#' -#' @param af_string_list \code{list} of length L containing \code{character} -#' strings with the names of the activation function used at each layer as the -#' names of the list expected by nn2poly. -#' -#' @return A list as expected in \code{nn2poly} output but without the nn2poly -#' class. -#' -#' @noRd -nn2poly_algorithm <- function(weights_list, - af_string_list, - max_order = 2, - keep_layers = FALSE, - taylor_orders = 8, - ..., - all_partitions = NULL - ) { - - if (!check_weights_dimensions(weights_list)) { - stop("The list of weights has incorrect dimensions. - Please, check the right dimmensions in the documentation.", - call. = FALSE - ) - } - - # Obtain number of variables (dimension p) - p <- dim(weights_list[[1]])[1] - 1 - - # Obtain number of layers L (L-1 hidden + 1 output, input is denoted by 0) - L <- length(af_string_list) - - # Initialize current layer in algorithm: - current_layer <- 1 - - # Check if the last layer is linear - last_linear <- (af_string_list[[L]] == "linear") - - # The list with the results of coefficients at each layer, - # depending on the last layer being linear or not - if (last_linear) { - results <- vector(mode = "list", length = 2 * L - 1) - } else { - results <- vector(mode = "list", length = 2 * L) - } - - # Create a default taylor_orders if it is not given by the user - taylor_orders <- obtain_taylor_vector( - taylor_orders = taylor_orders, - af_string_list = af_string_list - ) - - # Obtain all the derivatives up to the desired Taylor degree at each layer - af_derivatives_list <- obtain_derivatives_list( - af_string_list = af_string_list, - taylor_orders = taylor_orders - ) - - # Obtain the maximum degree of the final polynomial: - q_max <- obtain_final_poly_order(max_order, taylor_orders) - - # Check if partitions have not been given as an input - if (is.null(all_partitions)) { - all_partitions <- obtain_partitions_with_labels(p, q_max) - } - - ################# current_layer = 1, linear ################# - - # Starting point for the algorithm: Set weights as coefficients - # of an order 1 polynomial. - - # The labels for each coefficient vector at the same layer and linear or - # same layer and non linear case will have be the same, so they can be - # stored only once as the element `labels` of the coeffs_list_output, - # so we define length h+1 - coeffs_list_output <- vector(mode = "list", length = 0) - - # generate and store the labels (as a list of integer vectors) - # In this case the integer vectors are of length 1. - labels_output <- vector(mode = "list", length = p+1) - for (i in 0:p) { - labels_output[[i+1]] <- i - } - coeffs_list_output$labels <- labels_output - - # For each neuron in the first hidden layer, when computing the activation - # potentials (u_j), each column of the weight matrix represents the - # coefficients of an order 1 polynomial for that neuron potential. - # The first element will be the bias, and the rest the coefficient - # associated with each variable from x_1 to x_p. - coeffs_list_output$values <- t(weights_list[[1]]) - - # Store the results - results[[1]] <- coeffs_list_output - - # Stop if last layer and the last layer is linear - if (current_layer == L && last_linear) { - return(coeffs_list_output) - } - - ################# Loop over all layers ################# - - # Note that the loop will iterate the current layer from 1 to L - # and compute the linear and then non linear situation. - # The linear case at layer 1 has been computed outside so we skip it - - for (current_layer in 1:L) { - - ########## Linear case ########## - if (current_layer != 1){ - - # Treat the previous coefficients output as input - coeffs_list_input <- coeffs_list_output - - # Note that the polynomial in this case does not increase its order - # from the one in the non linear previous layer, so the labels - # will be the same and are already stored in $labels output. - # Only the matrix of $values will change its number of rows - - - ####### New version alg linear START ---------------------------------- - # apply the linear algorithm - values <- coeffs_list_input$values - coeffs_list_output$values <- t(weights_list[[current_layer]]) %*% - rbind(c(1,rep(0,ncol(values)-1)),values) - - ####### New version alg linear END ---------------------------------- - - # Save results from this layer: - results[[2 * (current_layer) - 1]] <- coeffs_list_output - - # Stop if last layer and the last layer is linear - if (current_layer == L && last_linear) { - if (keep_layers) { - return(results) - } else { - return(coeffs_list_output) - } - } - - } - - - ########## Non linear case ########## - # The output dimension remains the same as in the previous linear case, - # because the same number of neurons is considered - - # Treat the previous coeff output as input - coeffs_list_input <- coeffs_list_output - - # In the non linear case the polynomial order increases (unless max_order is - # reached), so the new labels need to be computed. However, the previous - # ones can be reused. - # The new labels will be for monomials of orders between the total order - # of the previous polynomial and the total order of the new polynomial: - if (current_layer == 1) { - previous_total_order <- 1 - } else { - # Get the total order used in the previous iteration - previous_total_order <- new_total_order - } - - # Compute the new total order with the product of taylor_orders. - # If a max_order value is used, its taken as the minimum between both. - if (is.null(max_order)) { - new_total_order <- previous_total_order * taylor_orders[current_layer] - } else { - new_total_order <- min(previous_total_order * taylor_orders[current_layer], - max_order) - } - - # If the order has increased, create new needed labels. - # If not, max_order has been reached and no new labels are needed. - if (previous_total_order != new_total_order){ - # Loop over each of the new orders up to the maximum one - for (order in (previous_total_order+1):new_total_order){ - combinations_indexes <- combinations_with_repetition(p, order) - - # Number of different combinations - n_combinations <- nrow(combinations_indexes) - - # Intialize list with new labels for the given order - new_labels <- vector(mode = "list", length = n_combinations) - - for (i in 1:n_combinations){ - new_labels[[i]] <- combinations_indexes[i, ] - } - - # update labels with the new ones - coeffs_list_output$labels <- c(coeffs_list_output$labels, new_labels) - - } - } - - # Parallel lapply - # The output index is already computed in the linear case - # but not for l=1 #REVISETHISLATER - values <- alg_non_linear( - coeffs_list_input$values, - labels_input = coeffs_list_input$labels, - labels_output = coeffs_list_output$labels, - taylor_orders = taylor_orders, - current_layer = current_layer, - g = af_derivatives_list[[current_layer]], - partitions_labels = all_partitions$labels, - partitions = all_partitions$partitions - ) - - # Join the coefficients with the labels as first list element: - coeffs_list_output$values <- values - - # Save results from this layer: - results[[2 * (current_layer)]] <- coeffs_list_output - - # Check if this is the last layer and if the last layer is not linear - if (current_layer == L && !last_linear) { - if (keep_layers) { - return(results) - } else { - return(coeffs_list_output) - } - } - - } -} - -#' This functions will generate the partitions obtained with Knuth's algorithm, -#' compute their labels and store both things in a list of length 2. -#' -#' -#' @param p Number of variables. -#' -#' @param q_max Maximum degree of the final polynomial. -#' -#' @return List of length 2 where the first element is a list with the labels -#' and the second element is a list with the partitions. -#' -#' @noRd -obtain_partitions_with_labels <- function(p, q_max) { - # This function will return a list with 2 elements: - # - # * The partitions obtained with Knuth's algorithm - # * The actual coefficient's "labels" for which the partitions are obtained - # - - if (missing(p) && missing(q_max)) { - stop("Missing both arguments p and q_max.", call. = FALSE) - } else if (missing(p)) { - stop("Missing argument p.", call. = FALSE) - } else if (missing(q_max)) { - stop("Missing argument q_max.", call. = FALSE) - } - - - partitions <- generate_partitions(as.integer(p), as.integer(q_max)) - - labels <- vector(mode = "list", length = length(partitions)) - # Obtain labels: - for (i in 1:length(partitions)){ - # Here it is used that the first partition of the multiset is always - # the multiset itself. This could be generalized in case we change the - # generation order. #REVISETHISLATER - labels[[i]] <- partitions[[i]][[1]][[1]] - } - - return(list("labels" = labels, "partitions" = partitions)) -} - -#' Computes the maximum polynomial order allowed by max_order and -#' Taylor orders at each layer -#' -#' Internal function used in nn2poly_algorithm. -#' Computes the final polynomial order by checking if max_order is an integer -#' and that the product of the Taylor order values allows for it. If it is -#' not reached, a warning is provided. -#' -#' @inheritParams nn2poly -#' @return An integer with the final polynomial order -#' -#' @noRd -obtain_final_poly_order <- function(max_order, taylor_orders){ - if(is.numeric(max_order) & (max_order %% 1)==0){ - # This condition allows for integers and also - # integers written as numeric - max_order <- as.integer(max_order) - poly_order <- min(prod(taylor_orders),max_order) - - # Warning if max_order has not been reached. (Very rare situation) - if (poly_order < max_order){ - warning("Argument `max_order` has not been reached due to chosen taylor_orders") - } - } else { - stop("Argument `max_order` is not an integer value", call. = FALSE) - } - - return(poly_order) -} - -#' Obtain a vector containing the Taylor order to be applied at each layer -#' -#' Internal function used in nn2poly_algorithm. -#' It allows the user to specify a single numeric value, in which case that -#' value will be employed at all no linear layers and 1 at linear layers. The -#' user can also provide their own vector, where in that case the dimensions -#' are checked so they match the number of layers provided by af_string_list. -#' -#' @inheritParams nn2poly -#' -#' @return An integer vector -#' -#' @noRd -obtain_taylor_vector <- function(taylor_orders, af_string_list){ - # Get the number of layers - L <- length(af_string_list) - - # First check if taylor_orders are integers or numeric with no decimal part: - if(!(is.numeric(taylor_orders) & all((taylor_orders %% 1)==0))){ - stop("Argument `taylor_orders` is non numeric", call. = FALSE) - } else if(length(taylor_orders)==1) { - # Single value case, set 1 in linear, taylor_orders in other AF - taylor_orders <- ifelse(af_string_list=="linear", 1, taylor_orders) - } else { - # Vector provided by user, check if dimensions match - if(!(length(taylor_orders)==length(af_string_list))){ - stop("Argument `taylor_orders` length does not match provided number of layers", - call. = FALSE) - } - } - return(taylor_orders) -} - -#' Obtain needed derivatives up to the chosen order (q Taylor) -#' -#' This function is internally used in nn2poly_algorithm to obtain the -#' derivatives of the given activation function at 0 up to the desired order $q$ -#' for each layer. -#' -#' @inheritParams nn2poly -#' @inheritParams nn2poly_algorithm -#' -#' @return list of vectors with the derivatives -#' -#' @noRd -obtain_derivatives_list <- function(af_string_list, taylor_orders) { - - n <- length(af_string_list) - af_function_list <- string_to_function(af_string_list) - af_derivatives_list <- vector(mode = "list", length = n) - - for (i in 1:n) { - # Obtain the vector with the derivatives of the activation function up to the given degree: - # centered at 0 - # and use rev to reverse and match our notation. - af_derivatives_list[[i]] <- rev(pracma::taylor(af_function_list[[i]], 0, taylor_orders[i])) - - # here we have a problem: if the last term of the taylor expansion is 0, - # the previous method deletes that entry and then the dimensions willm not macth later - # therefore, we add 0's if needed: - diff_len <- (taylor_orders[i] + 1) - length(af_derivatives_list[[i]]) - if (diff_len > 0) { - af_derivatives_list[[i]] <- c(af_derivatives_list[[i]], rep(0, diff_len)) - } - } - - return(af_derivatives_list) -} - -#' Checks that the weights dimensions are correct. -#' -#' This means that each matrix has the same number of rows as the number -#' of columns in the previous matrix + 1. This is because the number of -#' output neurons in the previous layer is the same as the number of inputs -#' to the current layer + the bias. -#' -#' @inheritParams nn2poly_algorithm -#' -#' @return `TRUE` if the dimensiones are correct, `FALSE` if not. -#' -#' @noRd -check_weights_dimensions <- function(weights_list) { - for (matrix_index in 2:length(weights_list)) { - nrows_current <- nrow(weights_list[[matrix_index]]) - ncols_prev <- ncol(weights_list[[matrix_index - 1]]) - - if (nrows_current != ncols_prev + 1) - return(FALSE) - } - return(TRUE) -} diff --git a/R/plot_taylor_and_activation_potentials.R b/R/plot_taylor_and_activation_potentials.R index af345a0..0de340e 100644 --- a/R/plot_taylor_and_activation_potentials.R +++ b/R/plot_taylor_and_activation_potentials.R @@ -64,10 +64,10 @@ plot_taylor_and_activation_potentials.list <- function(object, # } # Create a default taylor_orders if it is not given by the user (as in nn2poly) - taylor_orders <- obtain_taylor_vector( - taylor_orders = taylor_orders, - af_string_list = af_string_list - ) + taylor_orders <- obtain_taylor_vector(taylor_orders, af_string_list) + + # Obtain all the derivatives up to the desired Taylor degree at each layer + af_derivatives_list <- obtain_derivatives_list(taylor_orders, af_string_list) # The number of plots that we want to obtain is the number of hidden layers # (L-1) plus the output layer (L in total). @@ -82,7 +82,6 @@ plot_taylor_and_activation_potentials.list <- function(object, # get the AF as R functions: af_function_list <- string_to_function(af_string_list) - # We have to store the output of each layer to use it as input in the next one # and use it to compute the activation potentials. # Therefore, we initialize the variable "output" with data so the loop starts @@ -136,8 +135,9 @@ plot_taylor_and_activation_potentials.list <- function(object, # compute the true function yf <- fun(x) # compute the Taylor approximation - pol <- pracma::taylor(fun, 0, min(taylor_orders[k],max_order)) - yp <- pracma::polyval(pol, x) + q <- min(taylor_orders[[k]], max_order) + pol <- af_derivatives_list[[k]][1:(q + 1)] + yp <- as.vector(outer(x, 0:(length(pol) - 1), "^") %*% pol) # compute the error as the absolute value of the difference error <- abs(yf - yp) diff --git a/man/nn2poly.Rd b/man/nn2poly.Rd index d5d7cb1..776c105 100644 --- a/man/nn2poly.Rd +++ b/man/nn2poly.Rd @@ -4,14 +4,7 @@ \alias{nn2poly} \title{Obtain polynomial representation} \usage{ -nn2poly( - object, - max_order = 2, - keep_layers = FALSE, - taylor_orders = 8, - ..., - all_partitions = NULL -) +nn2poly(object, max_order = 2, keep_layers = FALSE, taylor_orders = 8, ...) } \arguments{ \item{object}{An object for which the computation of the NN2Poly algorithm is @@ -51,13 +44,6 @@ value is used, that value is set for each non linear layer and 1 for linear at each layer activation function. Default set to \code{8}.} \item{...}{Ignored.} - -\item{all_partitions}{Optional argument containing the needed multipartitions -as list of lists of lists. If set to \code{NULL}, nn2poly will compute said -multipartitions. This step can be computationally expensive when the chosen -polynomial order or the dimension are too high. In such cases, it is -encouraged that the multipartitions are stored and reused when possible. -Default set to \code{NULL}.} } \value{ Returns an object of class \code{nn2poly}. diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..22ebc63 --- /dev/null +++ b/src/Makevars @@ -0,0 +1 @@ +PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..22ebc63 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1 @@ +PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1360b43..cd0f901 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,6 +1,7 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#include "nn2poly_types.h" #include #include @@ -11,40 +12,89 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif -// select_allowed_partitions -std::vector> select_allowed_partitions(IntegerVector equivalent_label, int q_previous_layer, ListOf labels, List partitions); -RcppExport SEXP _nn2poly_select_allowed_partitions(SEXP equivalent_labelSEXP, SEXP q_previous_layerSEXP, SEXP labelsSEXP, SEXP partitionsSEXP) { +// obtain_final_poly_order +int obtain_final_poly_order(int max_order, const Term& taylor_orders); +RcppExport SEXP _nn2poly_obtain_final_poly_order(SEXP max_orderSEXP, SEXP taylor_ordersSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerVector >::type equivalent_label(equivalent_labelSEXP); - Rcpp::traits::input_parameter< int >::type q_previous_layer(q_previous_layerSEXP); - Rcpp::traits::input_parameter< ListOf >::type labels(labelsSEXP); - Rcpp::traits::input_parameter< List >::type partitions(partitionsSEXP); - rcpp_result_gen = Rcpp::wrap(select_allowed_partitions(equivalent_label, q_previous_layer, labels, partitions)); + Rcpp::traits::input_parameter< int >::type max_order(max_orderSEXP); + Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP); + rcpp_result_gen = Rcpp::wrap(obtain_final_poly_order(max_order, taylor_orders)); + return rcpp_result_gen; +END_RCPP +} +// obtain_taylor_vector +Term obtain_taylor_vector(const Term& taylor_orders, const Functions& af_string_list); +RcppExport SEXP _nn2poly_obtain_taylor_vector(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP); + Rcpp::traits::input_parameter< const Functions& >::type af_string_list(af_string_listSEXP); + rcpp_result_gen = Rcpp::wrap(obtain_taylor_vector(taylor_orders, af_string_list)); + return rcpp_result_gen; +END_RCPP +} +// obtain_derivatives_list +CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list); +RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP); + Rcpp::traits::input_parameter< const Functions& >::type af_string_list(af_string_listSEXP); + rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list)); + return rcpp_result_gen; +END_RCPP +} +// obtain_partitions_with_labels +PartitionsList obtain_partitions_with_labels(int p, int q_max); +RcppExport SEXP _nn2poly_obtain_partitions_with_labels(SEXP pSEXP, SEXP q_maxSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type p(pSEXP); + Rcpp::traits::input_parameter< int >::type q_max(q_maxSEXP); + rcpp_result_gen = Rcpp::wrap(obtain_partitions_with_labels(p, q_max)); return rcpp_result_gen; END_RCPP } // alg_non_linear -arma::mat alg_non_linear(arma::mat coeffs_input, ListOf labels_input, ListOf labels_output, IntegerVector taylor_orders, int current_layer, arma::vec g, ListOf partitions_labels, List partitions); +Weights alg_non_linear(const Weights& coeffs_input, const Terms& labels_input, const Terms& labels_output, const Term& taylor_orders, int current_layer, const Coeffs& g, const Terms& partitions_labels, const Partitions& partitions); RcppExport SEXP _nn2poly_alg_non_linear(SEXP coeffs_inputSEXP, SEXP labels_inputSEXP, SEXP labels_outputSEXP, SEXP taylor_ordersSEXP, SEXP current_layerSEXP, SEXP gSEXP, SEXP partitions_labelsSEXP, SEXP partitionsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type coeffs_input(coeffs_inputSEXP); - Rcpp::traits::input_parameter< ListOf >::type labels_input(labels_inputSEXP); - Rcpp::traits::input_parameter< ListOf >::type labels_output(labels_outputSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type taylor_orders(taylor_ordersSEXP); + Rcpp::traits::input_parameter< const Weights& >::type coeffs_input(coeffs_inputSEXP); + Rcpp::traits::input_parameter< const Terms& >::type labels_input(labels_inputSEXP); + Rcpp::traits::input_parameter< const Terms& >::type labels_output(labels_outputSEXP); + Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP); Rcpp::traits::input_parameter< int >::type current_layer(current_layerSEXP); - Rcpp::traits::input_parameter< arma::vec >::type g(gSEXP); - Rcpp::traits::input_parameter< ListOf >::type partitions_labels(partitions_labelsSEXP); - Rcpp::traits::input_parameter< List >::type partitions(partitionsSEXP); + Rcpp::traits::input_parameter< const Coeffs& >::type g(gSEXP); + Rcpp::traits::input_parameter< const Terms& >::type partitions_labels(partitions_labelsSEXP); + Rcpp::traits::input_parameter< const Partitions& >::type partitions(partitionsSEXP); rcpp_result_gen = Rcpp::wrap(alg_non_linear(coeffs_input, labels_input, labels_output, taylor_orders, current_layer, g, partitions_labels, partitions)); return rcpp_result_gen; END_RCPP } +// nn2poly_algorithm +List nn2poly_algorithm(const Layers& layers, const Functions& af_list, double max_order, bool keep_layers, const Term& taylor_orders); +RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Layers& >::type layers(layersSEXP); + Rcpp::traits::input_parameter< const Functions& >::type af_list(af_listSEXP); + Rcpp::traits::input_parameter< double >::type max_order(max_orderSEXP); + Rcpp::traits::input_parameter< bool >::type keep_layers(keep_layersSEXP); + Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP); + rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders)); + return rcpp_result_gen; +END_RCPP +} // combinations_with_repetition -IntegerMatrix combinations_with_repetition(int n, int k); +Terms combinations_with_repetition(int n, int k); RcppExport SEXP _nn2poly_combinations_with_repetition(SEXP nSEXP, SEXP kSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -56,7 +106,7 @@ BEGIN_RCPP END_RCPP } // generate_partitions -List generate_partitions(int p, int q_max); +Partitions generate_partitions(int p, int q_max); RcppExport SEXP _nn2poly_generate_partitions(SEXP pSEXP, SEXP q_maxSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -67,25 +117,16 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// generate_partitions_full -List generate_partitions_full(int p, int q_max); -RcppExport SEXP _nn2poly_generate_partitions_full(SEXP pSEXP, SEXP q_maxSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< int >::type p(pSEXP); - Rcpp::traits::input_parameter< int >::type q_max(q_maxSEXP); - rcpp_result_gen = Rcpp::wrap(generate_partitions_full(p, q_max)); - return rcpp_result_gen; -END_RCPP -} static const R_CallMethodDef CallEntries[] = { - {"_nn2poly_select_allowed_partitions", (DL_FUNC) &_nn2poly_select_allowed_partitions, 4}, + {"_nn2poly_obtain_final_poly_order", (DL_FUNC) &_nn2poly_obtain_final_poly_order, 2}, + {"_nn2poly_obtain_taylor_vector", (DL_FUNC) &_nn2poly_obtain_taylor_vector, 2}, + {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 2}, + {"_nn2poly_obtain_partitions_with_labels", (DL_FUNC) &_nn2poly_obtain_partitions_with_labels, 2}, {"_nn2poly_alg_non_linear", (DL_FUNC) &_nn2poly_alg_non_linear, 8}, + {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 5}, {"_nn2poly_combinations_with_repetition", (DL_FUNC) &_nn2poly_combinations_with_repetition, 2}, {"_nn2poly_generate_partitions", (DL_FUNC) &_nn2poly_generate_partitions, 2}, - {"_nn2poly_generate_partitions_full", (DL_FUNC) &_nn2poly_generate_partitions_full, 2}, {NULL, NULL, 0} }; diff --git a/src/algorithms.cpp b/src/algorithms.cpp deleted file mode 100644 index c128300..0000000 --- a/src/algorithms.cpp +++ /dev/null @@ -1,191 +0,0 @@ -// [[Rcpp::depends(RcppArmadillo)]] -#include -#include "utils.h" -using namespace Rcpp; - - -// [[Rcpp::export]] -std::vector> select_allowed_partitions( - IntegerVector equivalent_label, int q_previous_layer, - ListOf labels, List partitions) -{ - //REVISETHISLATER This function could be omitted if we already include it when - // generating the partitions. - - // Obtain chosen label position from the partitions labels list: - int pos; for (pos = 0; pos < labels.size(); pos++) { - if (labels[pos].size() != equivalent_label.size()) - continue; - if (is_true(all(labels[pos] == equivalent_label))) - break; - } - - // Extract the list with the partitions for the chosen label: - ListOf> all_partitions_for_this_label = - as>>(partitions[pos]); - - // number of partitions - int n_partitions = all_partitions_for_this_label.size(); - - // initialize list to store the output. - std::vector> output; - - for (int i = 0; i < n_partitions; i++) { - // Select one single partition - ListOf partition = - as>(all_partitions_for_this_label[i]); - - // Check that the given partition has all elements allowed by q_previous_layer. - // With sapply we sum the lengths of the vectors associated with each partition and see if - // they are less or equal than q_previous_layer. - if (is_true(all(sapply(partition, Rf_length) <= q_previous_layer))) - output.push_back(clone(partition)); - } - - return output; -} - - - -// [[Rcpp::export]] -arma::mat alg_non_linear(arma::mat coeffs_input, - ListOf labels_input, - ListOf labels_output, - IntegerVector taylor_orders, - int current_layer, arma::vec g, - ListOf partitions_labels, List partitions) -{ - // Extract the needed parameters and values: - int q_layer = taylor_orders[current_layer - 1]; - int q_previous_layer = 1; - if (current_layer != 1) - q_previous_layer = taylor_orders[current_layer - 2]; - - // Obtain total number of terms in the polynomial from labels - int n_poly_terms = labels_output.size(); - - // Obtain number of neurons - int h_l = coeffs_input.n_rows; - - // We define the vector that will contain all the output coefficients - arma::mat coeffs_output(h_l,n_poly_terms); - - ////////// Intercept ////////// - - for (int n = 0; n <= q_layer; n++) { - coeffs_output.col(0) = coeffs_output.col(0) + g[n] * arma::pow(coeffs_input.col(0), n); - // we have to use g[n] to obtain g^(n)/n!, - // because the function taylor already includes the term 1/n! - } - - ////////// Rest of the coefficients ////////// - - // As we already have all the coefficient labels, we can loop over them - // Note that the intercept has to be skipped so start at 1 - for (int coeff_index = 1; coeff_index < n_poly_terms; coeff_index++) { - IntegerVector label = labels_output[coeff_index]; - - // Find the equivalence between label and a the ones needed for the - // reduced partitions list - std::multiset mset(label.begin(), label.end()); - IntegerVector comp = unique(label).sort(); - IntegerVector mult(comp.size()); - for (int i = 0; i < comp.size(); i++) - mult[i] = mset.count(comp[i]); - comp = comp[order(mult, true)]; //decreasing - IntegerVector seq = Range(1, comp.size()); - // IntegerVector equivalent_label = - // concat(seq, label)[match(label, concat(comp, label)) - 1]; - IntegerVector equivalent_label = match(label, comp); - equivalent_label.sort(); - - // Obtain all allowed partitions of the equivalent term - auto allowed_partitions = select_allowed_partitions( - equivalent_label, q_previous_layer, partitions_labels, partitions); - - // Number of partitions - int n_allowed_partitions = allowed_partitions.size(); - - // Replace again all the partitions to match the original indexes - for (int p_index = 0; p_index < n_allowed_partitions; p_index++) { - ListOf aux = allowed_partitions[p_index]; - for (int i = 0; i < aux.size(); i++) { - IntegerVector auxv = aux[i]; - aux[i] = concat(comp, auxv)[match(auxv, concat(seq, auxv)) - 1]; - aux[i].sort(); - } - } - - // Now, use the correctly renamed partitions - for (int n = 1; n <= q_layer; n++) { - - arma::vec summatory(h_l); - - for (int p_index = 0; p_index < n_allowed_partitions; p_index++) { - // Extract the chosen partition (a list) from the allowed partitions - ListOf partition = allowed_partitions[p_index]; - - // We now need to check that each partition does not exceed n elements - // so we have the condition m_0 + ... + m_C = n satisfied. - // We also need the difference between the n_terms_in_partition - // with respect to n, so we can add that difference as the exponent - // of the intercept term. Then we compute this diff: - int difference = n - partition.size(); - - // If this diff is <0, we skip the partition - // This is due to the second restriction to the allowed partitions, that - // depends on n - if (difference < 0) continue; - - // We need to obtain the m_index values to compute the multinomial - // coefficient - - // This is simply counting how many times each unique term appears, - // obtaining the factorials and then doing the product. The terms that - // do not appear dont need to be counted as they will be 0, their - // factorial 1 and at the end will, not affect the total product. - - // This can be done as follows. //REVISETHISLATER when this used - // string vectors, it was easier to count all with table(). - // Now with vectors of different lengths in the vector this no - // longer works and a not so efficient workaround is used - ListOf unique_in_partition = Function("unique")(partition); - NumericVector m(unique_in_partition.size() + 1); - for (int i = 0; i < unique_in_partition.size(); i++) - m[i + 1] = sum(as(Function("%in%")( - partition, List::create(unique_in_partition[i])))); - m[0] = difference; - - // Compute the multinomial coefficient - NumericVector fm = factorial(m); - double multinomial_coef = std::tgamma(n + 1) / prod(fm); - - - // Now we need to use the labels to get the needed coefficients: - LogicalVector needed = Function("%in%")(labels_input, partition); - arma::mat coeffs_input_needed = coeffs_input.cols(find(as(needed) == 1)); - for (unsigned int i = 0; i < coeffs_input_needed.n_cols; i++) - coeffs_input_needed.col(i) = arma::pow(coeffs_input_needed.col(i), m[i + 1]); - - - // Finally compute the product of coefficients according to multinomial - // theorem and add it to the summatory - // For the product, it is sufficient to call prod(coeffs_input_needed) - // without including the exponent m, as this vector will contain - // each coefficient as many times as its exponent would indicate. - // REVISETHISLATER esto debería poder hacerse sin bucle con row product - summatory += multinomial_coef * - arma::prod(coeffs_input_needed,1) % arma::pow(coeffs_input.col(0), difference); - // Note that coeffs_input[0] is the intercept - - - } - // After the summatory over the partitions has been computed, we need to - // get its result and multiply by the correspondent derivative value, and - // add to the already stored values, here we are computing the summatory - // over n. - coeffs_output.col(coeff_index) = coeffs_output.col(coeff_index) + g[n] * summatory; - } - } - return coeffs_output; -} diff --git a/src/combinatorics.cpp b/src/combinatorics.cpp deleted file mode 100644 index 360ef81..0000000 --- a/src/combinatorics.cpp +++ /dev/null @@ -1,88 +0,0 @@ -#include -#include "multiset.h" -using namespace Rcpp; - -// [[Rcpp::export]] -IntegerMatrix combinations_with_repetition(int n, int k) { - int n_row = Rf_choose(n + k - 1, k); - IntegerMatrix out(n_row, k); - IntegerVector pos(k, 1); - - for (int row = 0; row < n_row; row++) { - for (int i = k - 1; i >= 0; i--) { - if (pos[i] > n) { - pos[i - 1]++; - for (int j = i; j < k; j++) - pos[j] = pos[j - 1]; - } - } - out(row, _) = pos; - pos[k - 1]++; - } - - return out; -} - -// [[Rcpp::export]] -List generate_partitions(int p, int q_max) { - // Initialize the output - std::list>> output; - - // Iterate over the degrees - for (int i = 1; i <= q_max; i++) { - // Generate all the combinations for the coeffs of the given order - IntegerMatrix comb = combinations_with_repetition(p, i); - - // Iterate over those combinations obtaining its partitions when needed - for (int j = 0; j < comb.nrow(); j++) { - // Count the occurences of each variable in the combination - std::multiset mset(comb(j, _).begin(), comb(j, _).end()); - - // Check if occurrences are sorted in descending order, which - // eliminates the equivalent situations (i.e., 112 is equivalent to 223, - // to 332, to 221 and so on) - int k = 1; - for (size_t last = mset.count(k++); k <= p; k++) { - size_t cur = mset.count(k); - if (last < cur) break; - last = cur; - } - if (k <= p) continue; - - // Add the list of partitions for that combination to the output - std::list> tmp; - for (auto partition: multiset_partitions(mset)) - tmp.push_back(partition); - output.push_back(tmp); - } - } - - return wrap(output); -} - - -// [[Rcpp::export]] -List generate_partitions_full(int p, int q_max) { - // Initialize the output - std::list>> output; - - // Iterate over the degrees - for (int i = 1; i <= q_max; i++) { - // Generate all the combinations for the coeffs of the given order - IntegerMatrix comb = combinations_with_repetition(p, i); - - // Iterate over those combinations obtaining its partitions - for (int j = 0; j < comb.nrow(); j++) { - // Count the occurences of each variable in the combination - std::multiset mset(comb(j, _).begin(), comb(j, _).end()); - - // Add the list of partitions for that combination to the output - std::list> tmp; - for (auto partition: multiset_partitions(mset)) - tmp.push_back(partition); - output.push_back(tmp); - } - } - - return wrap(output); -} diff --git a/src/multiset.h b/src/multiset.h index 2046fa6..b4b3eb7 100644 --- a/src/multiset.h +++ b/src/multiset.h @@ -7,16 +7,13 @@ #include #include -template -using Partition = std::vector>; - template class MultisetPartitions { std::vector comp; std::vector mult; public: - using value_type = Partition; + using value_type = std::vector>; using pointer = value_type*; using reference = value_type&; @@ -45,9 +42,12 @@ class MultisetPartitions { explicit operator bool() const { return !done; } value_type operator*() { return get(); } - pointer operator->() const { return &get(); } - iterator& operator++() { done = next(); return *this; } + iterator& operator++() { + done = next(); + Rcpp::checkUserInterrupt(); + return *this; + } iterator operator++(int) { auto it = *this; ++*this; return it; } friend bool operator==(const iterator& lhs, const iterator& rhs) { @@ -142,14 +142,13 @@ class MultisetPartitions { value_type get() { value_type partition; + partition.reserve(l + 1); for (int i = 0; i <= l; i++) { std::vector part; - for (int j = f[i]; j < f[i + 1]; j++) { + for (int j = f[i]; j < f[i + 1]; j++) for (int k = 0; k < v[j]; k++) part.push_back(obj->comp[c[j]]); - Rcpp::checkUserInterrupt(); - } - partition.push_back(part); + partition.push_back(std::move(part)); } return partition; } diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp new file mode 100644 index 0000000..e5082fe --- /dev/null +++ b/src/nn2poly.cpp @@ -0,0 +1,381 @@ +#include "partitions.h" +#include "taylor.h" + +// Forward declaration from partitions.cpp +Partitions generate_partitions(int p, int q_max); +Terms combinations_with_repetition(int n, int k); + +// [[Rcpp::export]] +int obtain_final_poly_order(int max_order, const Term& taylor_orders) { + long long product = 1; + for (int order : taylor_orders) + product *= order; + int poly_order = std::min(static_cast(product), static_cast(max_order)); + + // Warning if max_order has not been reached. (Very rare situation) + if (poly_order < max_order) + warning("Argument `max_order` has not been reached due to chosen taylor_Worders"); + + return poly_order; +} + +// [[Rcpp::export]] +Term obtain_taylor_vector(const Term& taylor_orders, + const Functions& af_string_list) { + if (taylor_orders.size() == 1) { + // Single value case, set 1 in linear, taylor_orders in other AF + Term out(af_string_list.size()); + const int base = taylor_orders[0]; + for (size_t i = 0; i < af_string_list.size(); i++) + out[i] = (af_string_list[i] == "linear") ? 1 : base; + return out; + } + + if (taylor_orders.size() != af_string_list.size()) + stop("Argument `taylor_orders` length does not match provided number of layers"); + + return taylor_orders; +} + +// [[Rcpp::export]] +CoeffsList obtain_derivatives_list(const Term& taylor_orders, + const Functions& af_string_list) { + if (taylor_orders.size() != af_string_list.size()) + stop("Argument `taylor_orders` length does not match provided number of layers"); + + CoeffsList out(af_string_list.size()); + for (size_t i = 0; i < af_string_list.size(); i++) { + if (taylor_orders[i] < 0) + stop("Argument `taylor_orders` must be non-negative"); + // Obtain the vector with the derivatives of the activation function up to + // the given degree centered at 0 + out[i] = coeffs_taylor(af_string_list[i], taylor_orders[i]); + } + + return out; +} + +// [[Rcpp::export]] +PartitionsList obtain_partitions_with_labels(int p, int q_max) { + // This function will return a list with 2 elements: + // * The partitions obtained with Knuth's algorithm + // * The actual coefficient's "labels" for which the partitions are obtained + + PartitionsList out; + out.partitions = generate_partitions(p, q_max); + out.labels.reserve(out.partitions.size()); + + for (size_t i = 0; i < out.partitions.size(); i++) { + if (out.partitions[i].empty() || out.partitions[i][0].empty()) + stop("Internal error while extracting partition labels."); + // Here it is used that the first partition of the multiset is always + // the multiset itself. This could be generalized in case we change the + // generation order. #REVISETHISLATER + out.labels.push_back(out.partitions[i][0][0]); + } + + return out; +} + +// [[Rcpp::export]] +Weights alg_non_linear(const Weights& coeffs_input, + const Terms& labels_input, + const Terms& labels_output, + const Term& taylor_orders, + int current_layer, + const Coeffs& g, + const Terms& partitions_labels, + const Partitions& partitions) { + // Extract the needed parameters and values: + const int q_layer = taylor_orders[current_layer - 1]; + int q_previous_layer = 1; + if (current_layer != 1) + q_previous_layer = taylor_orders[current_layer - 2]; + + // Obtain total number of terms in the polynomial from labels + const int n_poly_terms = static_cast(labels_output.size()); + + // Obtain number of neurons + const int h_l = static_cast(coeffs_input.n_rows); + + // We define the vector that will contain all the output coefficients + Weights coeffs_output(h_l, n_poly_terms, arma::fill::zeros); + + ////////// Intercept ////////// + + for (int n = 0; n <= q_layer; n++) { + coeffs_output.col(0) = coeffs_output.col(0) + g[n] * arma::pow(coeffs_input.col(0), n); + // we have to use g[n] to obtain g^(n)/n!, + // because the function taylor already includes the term 1/n! + } + + ////////// Rest of the coefficients ////////// + + // As we already have all the coefficient labels, we can loop over them + // Note that the intercept has to be skipped so start at 1 + for (int coeff_index = 1; coeff_index < n_poly_terms; coeff_index++) { + const Term& label = labels_output[coeff_index]; + + // Find the equivalence between label and a the ones needed for the + // reduced partitions list + const TermEquivalence equivalence = summarize_label_equivalence(label); + + // Obtain all allowed partitions of the equivalent term + Partition allowed_terms = filter_allowed_terms( + equivalence.signature, + q_previous_layer, + partitions_labels, + partitions + ); + + // Replace again all the partitions to match the original indexes + for (Terms& terms : allowed_terms) + terms = rename_terms(terms, equivalence.canonical_order); + + // Now, use the correctly renamed partitions + for (int n = 1; n <= q_layer; n++) { + arma::vec summatory(h_l, arma::fill::zeros); + + for (const Terms& terms : allowed_terms) { + // We now need to check that each partition does not exceed n elements + // so we have the condition m_0 + ... + m_C = n satisfied. + // We also need the difference between the n_terms_in_partition + // with respect to n, so we can add that difference as the exponent + // of the intercept term. Then we compute this diff: + const int difference = n - static_cast(terms.size()); + + // If this diff is <0, we skip the partition + // This is due to the second restriction to the allowed partitions, that + // depends on n + if (difference < 0) + continue; + + // We need to obtain the m_index values to compute the multinomial + // coefficient + + // This is simply counting how many times each unique term appears, + // obtaining the factorials and then doing the product. The terms that + // do not appear dont need to be counted as they will be 0, their + // factorial 1 and at the end will, not affect the total product. + TermSummary term_summary = summarize_terms(terms); + Term mult(term_summary.unique_terms.size() + 1, 0); + for (size_t i = 0; i < term_summary.unique_terms.size(); i++) { + auto it = term_summary.counts.find(term_summary.unique_terms[i]); + if (it == term_summary.counts.end()) + stop("Internal error while counting partition terms."); + mult[i + 1] = it->second; + } + mult[0] = difference; + + // Compute the multinomial coefficient + double multinomial_coef = std::tgamma(static_cast(n) + 1.0); + for (int m : mult) + multinomial_coef /= std::tgamma(static_cast(m) + 1.0); + + // Now we need to use the labels to get the needed coefficients: + const std::vector idx = in_terms_positions(labels_input, term_summary); + Weights coeffs_input_needed(h_l, idx.size()); + if (!idx.empty()) { + coeffs_input_needed = coeffs_input.cols(to_arma_indices(idx)); + for (unsigned int i = 0; i < coeffs_input_needed.n_cols; i++) + coeffs_input_needed.col(i) = arma::pow(coeffs_input_needed.col(i), mult[i + 1]); + } + + // Finally compute the product of coefficients according to multinomial + // theorem and add it to the summatory + // For the product, it is sufficient to call prod(coeffs_input_needed) + // without including the exponent m, as this vector will contain + // each coefficient as many times as its exponent would indicate. + // REVISETHISLATER esto debería poder hacerse sin bucle con row product + summatory += multinomial_coef * + arma::prod(coeffs_input_needed, 1) % arma::pow(coeffs_input.col(0), difference); + // Note that coeffs_input[0] is the intercept + } + // After the summatory over the partitions has been computed, we need to + // get its result and multiply by the correspondent derivative value, and + // add to the already stored values, here we are computing the summatory + // over n. + coeffs_output.col(coeff_index) = coeffs_output.col(coeff_index) + g[n] * summatory; + } + } + + return coeffs_output; +} + +inline void check_weights_dimensions(const Layers& layers) { + for (size_t i = 1; i < layers.size(); i++) { + if (layers[i].n_rows != layers[i - 1].n_cols + 1) + stop("The list of weights has incorrect dimensions. Please, check the right dimmensions in the documentation."); + } +} + +// [[Rcpp::export]] +List nn2poly_algorithm(const Layers& layers, + const Functions& af_list, + double max_order, + bool keep_layers, + const Term& taylor_orders) { + if (layers.empty()) + stop("Argument `layers` is empty"); + if (af_list.empty()) + stop("Activation functions are missing in `af_list`"); + if (layers.size() != af_list.size()) + stop("`layers` and `af_list` must have the same length"); + + check_weights_dimensions(layers); + + // Obtain number of variables (dimension p) + const int p = static_cast(layers[0].n_rows) - 1; + + // Obtain number of layers L (L-1 hidden + 1 output, input is denoted by 0) + const int L = static_cast(af_list.size()); + + // Check if the last layer is linear + const bool last_linear = (af_list.back() == "linear"); + + // The list with the results of coefficients at each layer, + // depending on the last layer being linear or not + WeightsLists results(static_cast(last_linear ? 2 * L - 1 : 2 * L)); + + // Create a default taylor_orders if it is not given by the user + const Term taylor = obtain_taylor_vector(taylor_orders, af_list); + + // Obtain all the derivatives up to the desired Taylor degree at each layer + const CoeffsList af_derivatives_list = obtain_derivatives_list(taylor, af_list); + + // Obtain the maximum degree of the final polynomial: + const int q_max = obtain_final_poly_order(max_order, taylor); + + // Obtain the partitions and their labels for the given p and q_max + const PartitionsList partition_data = obtain_partitions_with_labels(p, q_max); + + ////////////////// current_layer = 1, linear /////////////////// + + // Starting point for the algorithm: Set weights as coefficients + // of an order 1 polynomial. + + // The labels for each coefficient vector at the same layer and linear or + // same layer and non linear case will have be the same, so they can be + // stored only once as the element `labels` of the coeffs_list_output, + // so we define length h+1 + // For each neuron in the first hidden layer, when computing the activation + // potentials (u_j), each column of the weight matrix represents the + // coefficients of an order 1 polynomial for that neuron potential. + // The first element will be the bias, and the rest the coefficient + // associated with each variable from x_1 to x_p. + Weights coeffs_list_output = arma::trans(layers[0]); + + // generate and store the labels (as a list of integer vectors) + // In this case the integer vectors are of length 1. + Terms labels_output; + labels_output.reserve(p + 1); + for (int i = 0; i <= p; i++) + labels_output.push_back(Term{ i }); + + // Store the results + results[0] = WeightsList{labels_output, coeffs_list_output}; + + // Stop if last layer and the last layer is linear + int new_total_order = 1; + if (L == 1 && last_linear) + goto out; + + ////////////////// Loop over all layers /////////////////// + + // Note that the loop will iterate the current layer from 1 to L + // and compute the linear and then non linear situation. + // The linear case at layer 1 has been computed outside so we skip it + + for (int current_layer = 1; current_layer <= L; current_layer++) { + + ///////////////// Linear case ////////////////// + if (current_layer != 1) { + + // Treat the previous coefficients output as input + const Weights coeffs_list_input = coeffs_list_output; + + // Note that the polynomial in this case does not increase its order + // from the one in the non linear previous layer, so the labels + // will be the same and are already stored in $labels output. + // Only the matrix of $values will change its number of rows + + /////// New version alg linear START ---------------------------------- + // apply the linear algorithm + arma::rowvec first_row(coeffs_list_input.n_cols, arma::fill::zeros); + first_row[0] = 1.0; + Weights stacked(coeffs_list_input.n_rows + 1, coeffs_list_input.n_cols); + stacked.row(0) = first_row; + stacked.rows(1, coeffs_list_input.n_rows) = coeffs_list_input; + coeffs_list_output = arma::trans(layers[current_layer - 1]) * stacked; + + /////// New version alg linear END ---------------------------------- + + // Save results from this layer: + results[2 * current_layer - 2] = WeightsList{labels_output, coeffs_list_output}; + + // Stop if last layer and the last layer is linear + if (current_layer == L && last_linear) + goto out; + } + + ////////// Non linear case ////////// + // The output dimension remains the same as in the previous linear case, + // because the same number of neurons is considered + + // Treat the previous coeff output as input + const Weights coeffs_list_input = coeffs_list_output; + const Terms labels_input = labels_output; + + // In the non linear case the polynomial order increases (unless max_order is + // reached), so the new labels need to be computed. However, the previous + // ones can be reused. + // The new labels will be for monomials of orders between the total order + // of the previous polynomial and the total order of the new polynomial: + int previous_total_order = 1; + if (current_layer != 1) + previous_total_order = new_total_order; + + // Compute the new total order with the product of taylor_orders. + // If a max_order value is used, its taken as the minimum between both. + new_total_order = std::min(previous_total_order * taylor[current_layer - 1], + static_cast(max_order)); + + // If the order has increased, create new needed labels. + // If not, max_order has been reached and no new labels are needed. + if (previous_total_order != new_total_order) { + // Loop over each of the new orders up to the maximum one + for (int order = previous_total_order + 1; order <= new_total_order; order++) { + Terms comb = combinations_with_repetition(p, order); + + // update labels with the new ones + for (const Term& term : comb) + labels_output.push_back(term); + } + } + + // The output index is already computed in the linear case + // but not for l=1 #REVISETHISLATER + coeffs_list_output = alg_non_linear( + coeffs_list_input, + labels_input, + labels_output, + taylor, + current_layer, + af_derivatives_list[current_layer - 1], + partition_data.labels, + partition_data.partitions + ); + + // Save results from this layer: + results[2 * current_layer - 1] = WeightsList{labels_output, coeffs_list_output}; + + // Check if this is the last layer and if the last layer is not linear + if (current_layer == L && !last_linear) + goto out; + } + +out: + if (keep_layers) + return wrap(results); + return wrap(results.back()); +} diff --git a/src/nn2poly_types.h b/src/nn2poly_types.h new file mode 100644 index 0000000..d10d82a --- /dev/null +++ b/src/nn2poly_types.h @@ -0,0 +1,64 @@ +#ifndef nn2poly__types_h +#define nn2poly__types_h + +#include +#include + +using Term = std::vector; +using Terms = std::vector; +using Partition = std::vector; +using Partitions = std::vector; + +struct PartitionsList { + Terms labels; + Partitions partitions; +}; + +struct TermHash { + std::size_t operator()(Term const& key) const noexcept { + std::size_t seed = key.size(); + for(auto x : key) { + x = ((x >> 16) ^ x) * 0x45d9f3b; + x = ((x >> 16) ^ x) * 0x45d9f3b; + x = (x >> 16) ^ x; + seed ^= x + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } +}; + +struct TermSummary { + Terms unique_terms; + std::unordered_map counts; +}; + +struct TermEquivalence { + Term signature; + Term canonical_order; +}; + +using Coeffs = std::vector; +using CoeffsList = std::vector; + +// [[Rcpp::depends(RcppArmadillo)]] +#include +using namespace Rcpp; + +using Weights = arma::mat; +using Layers = std::vector; +using Functions = std::vector; + +struct WeightsList { + Terms labels; + Weights values; +}; + +using WeightsLists = std::vector; + +namespace Rcpp { +List wrap(const PartitionsList& data); +List wrap(const WeightsList& data); +List wrap(const WeightsLists& data); +} + +#endif \ No newline at end of file diff --git a/src/partitions.cpp b/src/partitions.cpp new file mode 100644 index 0000000..8d0cc2a --- /dev/null +++ b/src/partitions.cpp @@ -0,0 +1,69 @@ +#include "nn2poly_types.h" +#include "multiset.h" + +// [[Rcpp::export]] +Terms combinations_with_repetition(int n, int k) { + if (n <= 0 || k < 0) + stop("Arguments `n` and `k` must satisfy n > 0 and k >= 0"); + + Terms out; + if (k == 0) { + out.push_back(Term{}); + return out; + } + + Term pos(k, 1); + while (true) { + out.push_back(pos); + + int pivot = k - 1; + while (pivot >= 0 && pos[pivot] == n) + pivot--; + + if (pivot < 0) + break; + + pos[pivot]++; + for (int j = pivot + 1; j < k; j++) + pos[j] = pos[pivot]; + } + + return out; +} + +// [[Rcpp::export]] +Partitions generate_partitions(int p, int q_max) { + // Initialize the output + Partitions output; + + // Iterate over the degrees + for (int i = 1; i <= q_max; i++) { + // Generate all the combinations for the coeffs of the given order + Terms comb = combinations_with_repetition(p, i); + + // Iterate over those combinations obtaining its partitions when needed + for (size_t j = 0; j < comb.size(); j++) { + // Count the occurences of each variable in the combination + std::multiset mset(comb[j].begin(), comb[j].end()); + + // Check if occurrences are sorted in descending order, which + // eliminates the equivalent situations (i.e., 112 is equivalent to 223, + // to 332, to 221 and so on) + int k = 1; + for (size_t last = mset.count(k++); k <= p; k++) { + size_t cur = mset.count(k); + if (last < cur) break; + last = cur; + } + if (k <= p) continue; + + // Add the list of partitions for that combination to the output + Partition tmp; + for (auto terms: multiset_partitions(mset)) + tmp.push_back(terms); + output.push_back(tmp); + } + } + + return output; +} diff --git a/src/partitions.h b/src/partitions.h new file mode 100644 index 0000000..5fef007 --- /dev/null +++ b/src/partitions.h @@ -0,0 +1,126 @@ +#ifndef nn2poly__partitions_h +#define nn2poly__partitions_h + +#include "nn2poly_types.h" + +inline TermEquivalence summarize_label_equivalence(const Term& label) { + std::unordered_map counts; + for (int value : label) + counts[value]++; + + std::vector> frequency; + frequency.reserve(counts.size()); + for (const auto& entry : counts) + frequency.emplace_back(entry.first, entry.second); + + std::sort(frequency.begin(), frequency.end(), [](const auto& left, const auto& right) { + if (left.second != right.second) + return left.second > right.second; + return left.first < right.first; + }); + + Term canonical_order; + canonical_order.reserve(frequency.size()); + std::unordered_map rank_by_value; + for (size_t i = 0; i < frequency.size(); i++) { + canonical_order.push_back(frequency[i].first); + rank_by_value[frequency[i].first] = static_cast(i) + 1; + } + + Term signature; + signature.reserve(label.size()); + for (int value : label) + signature.push_back(rank_by_value[value]); + std::sort(signature.begin(), signature.end()); + + return {signature, canonical_order}; +} + +inline Terms rename_terms(const Terms& terms, const Term& canonical_order) { + Terms renamed; + renamed.reserve(terms.size()); + + for (const Term& term : terms) { + Term mapped; + mapped.reserve(term.size()); + for (int value : term) { + const int index = value - 1; + if (index < 0 || index >= static_cast(canonical_order.size())) { + stop("Internal error while renaming partition terms."); + } + mapped.push_back(canonical_order[index]); + } + std::sort(mapped.begin(), mapped.end()); + renamed.push_back(std::move(mapped)); + } + + return renamed; +} + +inline Partition filter_allowed_terms(const Term& equivalent_label, + int q_previous_layer, + const Terms& labels, + const Partitions& partitions) { + //REVISETHISLATER This function could be omitted if we already include it when + // generating the partitions. + + // Obtain chosen label position from the partitions labels list: + auto label_it = std::find(labels.begin(), labels.end(), equivalent_label); + if (label_it == labels.end()) + stop("Internal error while locating the equivalent partition label."); + const size_t label_index = static_cast(std::distance(labels.begin(), label_it)); + + Partition output; + for (const Terms& terms : partitions[label_index]) { + bool allowed = true; + // Check that the given partition has all elements allowed by q_previous_layer. + for (const Term& term : terms) { + if (static_cast(term.size()) > q_previous_layer) { + allowed = false; + break; + } + } + if (allowed) + output.push_back(terms); + } + + return output; +} + +inline arma::uvec to_arma_indices(const std::vector& positions) { + arma::uvec indices(positions.size()); + for (size_t i = 0; i < positions.size(); i++) + indices[i] = positions[i]; + return indices; +} + +inline TermSummary summarize_terms(const Terms& terms) { + TermSummary out; + out.unique_terms.reserve(terms.size()); + out.counts.reserve(terms.size()); + + for (size_t i = 0; i < terms.size(); i++) { + const Term& term = terms[i]; + Term key = term; + auto it = out.counts.find(key); + if (it == out.counts.end()) { + out.counts.emplace(std::move(key), 1); + out.unique_terms.push_back(term); + } else it->second++; + } + + return out; +} + +inline std::vector in_terms_positions(const Terms& labels_input, + const TermSummary& term_summary) { + std::vector needed; + needed.reserve(labels_input.size()); + for (size_t i = 0; i < labels_input.size(); i++) { + if (term_summary.counts.find(labels_input[i]) != term_summary.counts.end()) + needed.push_back(i); + } + return needed; +} + +#endif diff --git a/src/taylor.h b/src/taylor.h new file mode 100644 index 0000000..16845f9 --- /dev/null +++ b/src/taylor.h @@ -0,0 +1,60 @@ +#ifndef nn2poly__taylor_h +#define nn2poly__taylor_h + +#include "nn2poly_types.h" + +Coeffs coeffs_taylor_sigmoid(int order) { + Coeffs a(order + 1, 0.0); + a[0] = 0.5; + for (int n = 0; n < order; n++) { + double conv = 0.0; + for (int k = 0; k <= n; k++) + conv += a[k] * a[n - k]; + a[n + 1] = (a[n] - conv) / static_cast(n + 1); + } + return a; +} + +Coeffs coeffs_taylor_tanh(int order) { + Coeffs a(order + 1, 0.0); + for (int n = 0; n < order; n++) { + double conv = 0.0; + for (int k = 0; k <= n; k++) + conv += a[k] * a[n - k]; + const double rhs = (n == 0 ? 1.0 : 0.0) - conv; + a[n + 1] = rhs / static_cast(n + 1); + } + return a; +} + +Coeffs coeffs_taylor_linear(int order) { + Coeffs a(order + 1, 0.0); + if (order >= 1) + a[1] = 1.0; + return a; +} + +Coeffs coeffs_taylor_softplus(int order) { + Coeffs a(order + 1, 0.0); + a[0] = std::log(2.0); + if (order == 0) + return a; + Coeffs sig = coeffs_taylor_sigmoid(order - 1); + for (int n = 1; n <= order; n++) + a[n] = sig[n - 1] / static_cast(n); + return a; +} + +Coeffs coeffs_taylor(const std::string& af, int order) { + if (af == "sigmoid") + return coeffs_taylor_sigmoid(order); + if (af == "tanh") + return coeffs_taylor_tanh(order); + if (af == "softplus") + return coeffs_taylor_softplus(order); + if (af == "linear") + return coeffs_taylor_linear(order); + throw std::invalid_argument("Function '" + af + "' not supported"); +} + +#endif diff --git a/src/utils.h b/src/utils.h deleted file mode 100644 index 232b581..0000000 --- a/src/utils.h +++ /dev/null @@ -1,38 +0,0 @@ -#ifndef nn2poly__utils_h -#define nn2poly__utils_h - -#include -using namespace Rcpp; - -template -IntegerVector order(Vector& x, bool decreasing=false) { - IntegerVector out(x.size()); - std::iota(out.begin(), out.end(), 0); - std::sort(out.begin(), out.end(), [&](size_t a, size_t b) { - if (decreasing) - return x[a] > x[b]; - return x[a] < x[b]; - }); - return out; -} - -template -Vector concat(Vector& x, Vector&y) { - int xn = x.size(); - int yn = y.size(); - Vector out(xn + yn); - for (int i = 0; i < xn; i++) - out[i] = x[i]; - for (int i = 0; i < yn; i++) - out[xn + i] = y[i]; - return out; -} - -template -typename traits::storage_type::type prod(Vector& x) { - using type = typename traits::storage_type::type; - type init = 1; - return std::accumulate(x.begin(), x.end(), init, std::multiplies()); -} - -#endif diff --git a/src/wrap.cpp b/src/wrap.cpp new file mode 100644 index 0000000..a20a329 --- /dev/null +++ b/src/wrap.cpp @@ -0,0 +1,53 @@ +#include "nn2poly_types.h" + +namespace Rcpp { + +List wrap(const PartitionsList& data) { + return List::create( + _["labels"] = data.labels, + _["partitions"] = data.partitions + ); +} + +List wrap(const WeightsList& data) { + return List::create( + _["labels"] = data.labels, + // Transpose to have polynomials as columns + _["values"] = arma::trans(data.values) + ); +} + +List wrap(const WeightsLists& data) { + const int n_items = static_cast(data.size()); + const int n_layers = static_cast(std::ceil(static_cast(n_items) / 2.0)); + + List out(n_layers); + CharacterVector out_names(n_layers); + + for (int i = 0; i < n_layers; i++) { + const int input_index = 2 * i; + const int output_index = 2 * i + 1; + + List layer(2); + CharacterVector layer_names = CharacterVector::create("input", "output"); + + layer[0] = wrap(data[input_index]); + if (output_index < n_items) { + layer[1] = wrap(data[output_index]); + } else { + // If there is a linear output, i.e. single polynomial in final layer + // and odd number of items, then we repeat the input as the output, as + // the activation functions takes no effect on the polynomial + layer[1] = layer[0]; + } + + layer.attr("names") = layer_names; + out[i] = layer; + out_names[i] = "layer_" + std::to_string(i + 1); + } + + out.attr("names") = out_names; + return out; +} + +} diff --git a/tests/testthat/_snaps/combinatorics.md b/tests/testthat/_snaps/combinatorics.md index f0e58c2..97ffbd5 100644 --- a/tests/testthat/_snaps/combinatorics.md +++ b/tests/testthat/_snaps/combinatorics.md @@ -3,70 +3,180 @@ Code combinations_with_repetition(5, 3) Output - [,1] [,2] [,3] - [1,] 1 1 1 - [2,] 1 1 2 - [3,] 1 1 3 - [4,] 1 1 4 - [5,] 1 1 5 - [6,] 1 2 2 - [7,] 1 2 3 - [8,] 1 2 4 - [9,] 1 2 5 - [10,] 1 3 3 - [11,] 1 3 4 - [12,] 1 3 5 - [13,] 1 4 4 - [14,] 1 4 5 - [15,] 1 5 5 - [16,] 2 2 2 - [17,] 2 2 3 - [18,] 2 2 4 - [19,] 2 2 5 - [20,] 2 3 3 - [21,] 2 3 4 - [22,] 2 3 5 - [23,] 2 4 4 - [24,] 2 4 5 - [25,] 2 5 5 - [26,] 3 3 3 - [27,] 3 3 4 - [28,] 3 3 5 - [29,] 3 4 4 - [30,] 3 4 5 - [31,] 3 5 5 - [32,] 4 4 4 - [33,] 4 4 5 - [34,] 4 5 5 - [35,] 5 5 5 + [[1]] + [1] 1 1 1 + + [[2]] + [1] 1 1 2 + + [[3]] + [1] 1 1 3 + + [[4]] + [1] 1 1 4 + + [[5]] + [1] 1 1 5 + + [[6]] + [1] 1 2 2 + + [[7]] + [1] 1 2 3 + + [[8]] + [1] 1 2 4 + + [[9]] + [1] 1 2 5 + + [[10]] + [1] 1 3 3 + + [[11]] + [1] 1 3 4 + + [[12]] + [1] 1 3 5 + + [[13]] + [1] 1 4 4 + + [[14]] + [1] 1 4 5 + + [[15]] + [1] 1 5 5 + + [[16]] + [1] 2 2 2 + + [[17]] + [1] 2 2 3 + + [[18]] + [1] 2 2 4 + + [[19]] + [1] 2 2 5 + + [[20]] + [1] 2 3 3 + + [[21]] + [1] 2 3 4 + + [[22]] + [1] 2 3 5 + + [[23]] + [1] 2 4 4 + + [[24]] + [1] 2 4 5 + + [[25]] + [1] 2 5 5 + + [[26]] + [1] 3 3 3 + + [[27]] + [1] 3 3 4 + + [[28]] + [1] 3 3 5 + + [[29]] + [1] 3 4 4 + + [[30]] + [1] 3 4 5 + + [[31]] + [1] 3 5 5 + + [[32]] + [1] 4 4 4 + + [[33]] + [1] 4 4 5 + + [[34]] + [1] 4 5 5 + + [[35]] + [1] 5 5 5 + --- Code combinations_with_repetition(3, 5) Output - [,1] [,2] [,3] [,4] [,5] - [1,] 1 1 1 1 1 - [2,] 1 1 1 1 2 - [3,] 1 1 1 1 3 - [4,] 1 1 1 2 2 - [5,] 1 1 1 2 3 - [6,] 1 1 1 3 3 - [7,] 1 1 2 2 2 - [8,] 1 1 2 2 3 - [9,] 1 1 2 3 3 - [10,] 1 1 3 3 3 - [11,] 1 2 2 2 2 - [12,] 1 2 2 2 3 - [13,] 1 2 2 3 3 - [14,] 1 2 3 3 3 - [15,] 1 3 3 3 3 - [16,] 2 2 2 2 2 - [17,] 2 2 2 2 3 - [18,] 2 2 2 3 3 - [19,] 2 2 3 3 3 - [20,] 2 3 3 3 3 - [21,] 3 3 3 3 3 + [[1]] + [1] 1 1 1 1 1 + + [[2]] + [1] 1 1 1 1 2 + + [[3]] + [1] 1 1 1 1 3 + + [[4]] + [1] 1 1 1 2 2 + + [[5]] + [1] 1 1 1 2 3 + + [[6]] + [1] 1 1 1 3 3 + + [[7]] + [1] 1 1 2 2 2 + + [[8]] + [1] 1 1 2 2 3 + + [[9]] + [1] 1 1 2 3 3 + + [[10]] + [1] 1 1 3 3 3 + + [[11]] + [1] 1 2 2 2 2 + + [[12]] + [1] 1 2 2 2 3 + + [[13]] + [1] 1 2 2 3 3 + + [[14]] + [1] 1 2 3 3 3 + + [[15]] + [1] 1 3 3 3 3 + + [[16]] + [1] 2 2 2 2 2 + + [[17]] + [1] 2 2 2 2 3 + + [[18]] + [1] 2 2 2 3 3 + + [[19]] + [1] 2 2 3 3 3 + + [[20]] + [1] 2 3 3 3 3 + + [[21]] + [1] 3 3 3 3 3 + # multiset partitions are correctly generated diff --git a/tests/testthat/_snaps/nn2poly_methods_plot/top-null.svg b/tests/testthat/_snaps/nn2poly_methods_plot/top-null.svg deleted file mode 100644 index a14879c..0000000 --- a/tests/testthat/_snaps/nn2poly_methods_plot/top-null.svg +++ /dev/null @@ -1,65 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -2.5 -5.0 -7.5 - - - - - - - - - - - -0 -1 -2 -1,2 -1,1 -2,2 -Variables or interactions -Coefficient (absolute) values -Sign - -+ -top NULL - - diff --git a/tests/testthat/helper-check_constraints.R b/tests/testthat/helper-check_constraints.R index dbce5a8..b2d1092 100644 --- a/tests/testthat/helper-check_constraints.R +++ b/tests/testthat/helper-check_constraints.R @@ -1,3 +1,13 @@ +Norm <- function (x, p = 2) { + stopifnot(is.numeric(x) || is.complex(x), is.numeric(p), length(p) == 1) + if (p > -Inf && p < Inf) + sum(abs(x)^p)^(1/p) + else if (p == Inf) + max(abs(x)) + else if (p == -Inf) + min(abs(x)) + else return(NULL) +} check_weight_constraints <- function(weights, maxnorm) { #### Compute the norm for the full matrix #### @@ -19,17 +29,17 @@ check_weight_constraints <- function(weights, maxnorm) { # Loop over each column for (j in 1:n_cols) { if (maxnorm[[1]] == "l1_norm") { - norm <- pracma::Norm(weights[[i]][, j], 1) + norm <- Norm(weights[[i]][, j], 1) weights_norms[[i]][j] <- norm } else if (maxnorm[[1]] == "l2_norm") { - norm <- pracma::Norm(weights[[i]][, j], 2) + norm <- Norm(weights[[i]][, j], 2) weights_norms[[i]][j] <- norm } else if (maxnorm[[1]] == "unit") { - norm <- pracma::Norm(weights[[i]][, j], 2) + norm <- Norm(weights[[i]][, j], 2) weights_norms[[i]][j] <- norm } else { print("Imprecise norm. Computing the l1 norm...") - norm <- pracma::Norm(weights[[i]][, j], 1) + norm <- Norm(weights[[i]][, j], 1) weights_norms[[i]][j] <- norm } diff --git a/tests/testthat/helper-testing_data.R b/tests/testthat/helper-testing_data.R index c78444d..c88cc18 100644 --- a/tests/testthat/helper-testing_data.R +++ b/tests/testthat/helper-testing_data.R @@ -8,7 +8,7 @@ testing_helper_1 <- function(){ testing_data_1$weights_list[[1]] <- matrix(1,3,2) testing_data_1$weights_list[[2]] <- matrix(1,3,2) testing_data_1$weights_list[[3]] <- matrix(1,3,1) - testing_data_1$af_string_list <- list("softplus", "softplus", "linear") + testing_data_1$af_string_list <- c("softplus", "softplus", "linear") testing_data_1$taylor_orders <- c(2, 2, 1) return(testing_data_1) } diff --git a/tests/testthat/test-nn2poly.R b/tests/testthat/test-nn2poly.R index fc811bf..b15fa20 100644 --- a/tests/testthat/test-nn2poly.R +++ b/tests/testthat/test-nn2poly.R @@ -50,7 +50,7 @@ test_that("nn2poly with list input against precomputed example with label <- result$labels[[4]] coeff <- result$values[4,1] expect_equal(label,c(1,1)) - expect_equal(round(coeff,3),-4.429) + expect_equal(round(coeff,2),-4.55) }) diff --git a/tests/testthat/test-nn2poly_algorithm.R b/tests/testthat/test-nn2poly_algorithm.R index 295f686..f4f9f04 100644 --- a/tests/testthat/test-nn2poly_algorithm.R +++ b/tests/testthat/test-nn2poly_algorithm.R @@ -10,23 +10,24 @@ test_that("nn2poly_algorithm against precomputed example", { taylor_orders <- testing_data$taylor_orders result <- nn2poly_algorithm( - weights_list = weights_list, - af_string_list = af_string_list, + weights_list, + af_string_list, max_order = 3, keep_layers = TRUE, taylor_orders = taylor_orders ) - n_terms <- length(result[[length(result)]]$labels) - order <- length(result[[length(result)]]$labels[[n_terms]]) + n_layers <- length(result) + n_terms <- length(result[[n_layers]]$output$labels) + order <- length(result[[n_layers]]$output$labels[[n_terms]]) expect_equal(order, 3) - # Desired coefficient in output polynomial at layer 2 (element 2*2=4 in list), + # Desired coefficient in output polynomial at layer 2, # neuron 1, coefficient "1,1" - label <- result[[4]]$labels[[4]] - coeff <- result[[4]]$values[1,4] - expect_equal(label,c(1,1)) - expect_equal(coeff,0.63351833) + label <- result$layer_2$output$labels[[4]] + coeff <- result$layer_2$output$values[4, 1] + expect_equal(label, c(1, 1)) + expect_equal(coeff, 0.63351833, tolerance = 1e-6) }) @@ -135,15 +136,6 @@ test_that("Final order warns if max_order is not reached", { expect_equal(output, 6L) }) - -test_that("Final order stops with an error if max_order is not an integer", { - max_order <- 4.3 - # Note that Taylor orders will always be a vector in the desired form - # as it is built inside nn2poly_algorithm - taylor_orders <- c(2,1,3) - expect_error(obtain_final_poly_order(max_order, taylor_orders)) -}) - test_that("Taylor vector obtained with single value and mutiple linear layers", { af_string_list <- c("softplus", "linear", "softplus", "linear") taylor_orders <- 5L @@ -163,10 +155,3 @@ test_that("Taylor vector gets error because of dimension missmatch", { taylor_orders <- c(5,1) expect_error(obtain_taylor_vector(taylor_orders, af_string_list)) }) - - -test_that("Taylor vector gets error becauseof non numeric value", { - af_string_list <- c("softplus", "softplus", "linear") - taylor_orders <- c(5.4, 2.3, 1) - expect_error(obtain_taylor_vector(taylor_orders, af_string_list)) -})