diff --git a/NAMESPACE b/NAMESPACE
index 927db13..c80ce2b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@
 export(admix1)
 export(calculatePopLogml)
 export(colon)
+export(computeAllFreqs2)
 export(computeIndLogml)
 export(computePersonalAllFreqs)
 export(computeRows)
diff --git a/R/admix1.R b/R/admix1.R
index 28ce731..c50d4fe 100644
--- a/R/admix1.R
+++ b/R/admix1.R
@@ -5,7 +5,7 @@
 admix1 <- function(tietue) {
     if (!is.list(tietue)) {
 #         c(filename, pathname) = uigetfile('*.mat', 'Load mixture result file');
-#         if (filename==0 & pathname==0), return; 
+#         if (filename==0 & pathname==0), return;
 #         else
 #             disp('---------------------------------------------------');
 #             disp(['Reading mixture result from: ',[pathname filename],'...']);
@@ -13,7 +13,7 @@ admix1 <- function(tietue) {
 #         pause(0.0001);
 #         h0 = findobj('Tag','filename1_text');
 #         set(h0,'String',filename); clear h0;
-        
+
 #         struct_array = load([pathname filename]);
 #         if isfield(struct_array,'c')  #Matlab versio
 #             c = struct_array.c;
@@ -31,7 +31,7 @@ admix1 <- function(tietue) {
 #             disp('Incorrect file format');
 #             return;
 #         end
-        
+
 #         if isfield(c, 'gene_lengths') && ...
 #                 (strcmp(c.mixtureType,'linear_mix') | ...
 #                 strcmp(c.mixtureType,'codon_mix')) # if the mixture is from a linkage model
@@ -40,7 +40,7 @@ admix1 <- function(tietue) {
 #             linkage_admix(c);
 #             return
 #         end
-        
+
 #         PARTITION = c.PARTITION; COUNTS = c.COUNTS; SUMCOUNTS = c.SUMCOUNTS;
 #         alleleCodes = c.alleleCodes; adjprior = c.adjprior; popnames = c.popnames;
 #         rowsFromInd = c.rowsFromInd; data = c.data; npops = c.npops; noalle = c.noalle;
@@ -148,7 +148,7 @@ admix1 <- function(tietue) {
 #     for iterationNum = 1:iterationCount
 #         disp(['Iter: ' num2str(iterationNum)]);
 #         allfreqs = simulateAllFreqs(noalle);   # Allele frequencies on this iteration.
-        
+
 #         for ind=to_investigate
 #             #disp(num2str(ind));
 #             omaFreqs = computePersonalAllFreqs(ind, data, allfreqs, rowsFromInd);
@@ -172,7 +172,7 @@ admix1 <- function(tietue) {
 #                     PARTITION(ind)=isoimman_indeksi;
 #                 end
 #                 logml = computeIndLogml(omaFreqs, osuusTaulu);
-                
+
 #                 for osuus = [0.5 0.25 0.05 0.01]
 #                     [osuusTaulu, logml] = etsiParas(osuus, osuusTaulu, omaFreqs, logml);
 #                 end
@@ -227,23 +227,23 @@ admix1 <- function(tietue) {
 #     for pop = admix_populaatiot'
 
 #         for level = 1:n_missing_levels(pop)
-            
+
 #             potential_inds_in_this_pop_and_level = ...
 #                 find(PARTITION==pop & missing_level_partition==level &...
 #                 likelihood>3);  # Potential admix individuals here.
-            
+
 #             if ~isempty(potential_inds_in_this_pop_and_level)
-            
+
 #                 #refData = simulateIndividuals(nrefIndsInPop,rowsFromInd,allfreqs);
 #                 refData = simulateIndividuals(nrefIndsInPop, rowsFromInd, allfreqs, ...
 #                     pop, missing_levels(pop,level));
-            
+
 #                 disp(['Analysing the reference individuals from pop ' num2str(pop) ' (level ' num2str(level) ').']);
 #                 refProportions = zeros(nrefIndsInPop,npops);
 #                 for iter = 1:iterationCountRef
 #                     #disp(['Iter: ' num2str(iter)]);
 #                     allfreqs = simulateAllFreqs(noalle);
-                
+
 #                     for ind = 1:nrefIndsInPop
 #                         omaFreqs = computePersonalAllFreqs(ind, refData, allfreqs, rowsFromInd);
 #                         osuusTaulu = zeros(1,npops);
@@ -312,7 +312,7 @@ admix1 <- function(tietue) {
 #         end
 #     end
 
-#     tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount); 
+#     tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount);
 
 #     viewPartition(proportionsIt, popnames);
 
@@ -335,7 +335,7 @@ admix1 <- function(tietue) {
 
 
 #         if (~isstruct(tietue))
-#             c.proportionsIt = proportionsIt; 
+#             c.proportionsIt = proportionsIt;
 #             c.pvalue = uskottavuus; # Added by Jing
 #             c.mixtureType = 'admix'; # Jing
 #             c.admixnpops = npops;
@@ -355,7 +355,7 @@ admix1 <- function(tietue) {
 
 # function [npops] = poistaLiianPienet(npops, rowsFromInd, alaraja)
 # % Muokkaa tulokset muotoon, jossa outlier yksilöt on
-# % poistettu. Tarkalleen ottaen poistaa ne populaatiot, 
+# % poistettu. Tarkalleen ottaen poistaa ne populaatiot,
 # % joissa on vähemmän kuin 'alaraja':n verran yksilöit?
 
 # global PARTITION;
@@ -400,30 +400,4 @@ admix1 <- function(tietue) {
 # global COUNTS; COUNTS = [];
 # global SUMCOUNTS; SUMCOUNTS = [];
 # global PARTITION; PARTITION = [];
-# global POP_LOGML; POP_LOGML = [];
-
-# %--------------------------------------------------------
-
-
-# function allFreqs = computeAllFreqs2(noalle)
-# % Lisää a priori jokaista alleelia
-# % joka populaation joka lokukseen j 1/noalle(j) verran.
-
-# global COUNTS;
-# global SUMCOUNTS;
-
-# max_noalle = size(COUNTS,1);
-# nloci = size(COUNTS,2);
-# npops = size(COUNTS,3);
-
-# sumCounts = SUMCOUNTS+ones(size(SUMCOUNTS));
-# sumCounts = reshape(sumCounts', [1, nloci, npops]);
-# sumCounts = repmat(sumCounts, [max_noalle, 1 1]);
-
-# prioriAlleelit = zeros(max_noalle,nloci);
-# for j=1:nloci
-#     prioriAlleelit(1:noalle(j),j) = 1/noalle(j);
-# end
-# prioriAlleelit = repmat(prioriAlleelit, [1,1,npops]);
-# counts = COUNTS + prioriAlleelit;
-# allFreqs = counts./sumCounts;
\ No newline at end of file
+# global POP_LOGML; POP_LOGML = [];
\ No newline at end of file
diff --git a/R/computeAllFreqs2.R b/R/computeAllFreqs2.R
new file mode 100644
index 0000000..e3f68f8
--- /dev/null
+++ b/R/computeAllFreqs2.R
@@ -0,0 +1,29 @@
+#' @title Compute all freqs - version 2
+#' @description Lisää a priori jokaista alleelia joka populaation joka lokukseen
+#' j 1/noalle(j) verran.
+#' @param noalle noalle
+#' @param COUNTS counts
+#' @param SUMCOUNTS sumcounts
+#' @export
+computeAllFreqs2 <- function (noalle, COUNTS = matrix(NA, 0, 0),
+    SUMCOUNTS = NA) {
+
+    max_noalle <- size(COUNTS, 1)
+    nloci <- size(COUNTS,2)
+    npops <- size(COUNTS,3)
+
+    sumCounts <- SUMCOUNTS + ones(size(SUMCOUNTS))
+    sumCounts <- reshape(t(sumCounts), c(1, nloci, npops))
+    sumCounts <- repmat(sumCounts, c(max_noalle, 1, 1))
+
+    prioriAlleelit <- zeros(max_noalle, nloci)
+    if (nloci > 0) {
+        for (j in 1:nloci) {
+            prioriAlleelit[1:noalle[j], j] <- 1 / noalle[j]
+        }
+    }
+    prioriAlleelit <- repmat(prioriAlleelit, c(1, 1, npops))
+    counts <- COUNTS + prioriAlleelit
+    allFreqs <- counts / sumCounts
+    return(allFreqs)
+}
\ No newline at end of file
diff --git a/R/repmat.R b/R/repmat.R
index be72325..5c5e92f 100644
--- a/R/repmat.R
+++ b/R/repmat.R
@@ -3,15 +3,16 @@
 #' @details This function was created to replicate the behavior of a homonymous
 #' function on Matlab
 #' @param mx matrix
-#' @param n either a scalar with the number of replications in both rows and columns or a 2-length vector with individual repetitions.
+#' @param n either a scalar with the number of replications in both rows and
+#' columns or a <= 3-length vector with individual repetitions.
 #' @return matrix replicated over `ncol(mx) * n` columns and `nrow(mx) * n` rows
 #' @note The Matlab implementation of this function accepts `n` with length > 2.
-#' 
+#'
 #' It should also be noted that a concatenated vector in R, e.g. `c(5, 2)`, becomes a column vector when coerced to matrix, even though it may look like a row vector at first glance. This is important to keep in mind when considering the expected output of this function. Vectors in R make sense to be seen as column vectors, given R's Statistics-oriented paradigm where variables are usually disposed as columns in a dataset.
 #' @export
 repmat <- function (mx, n) {
     # Validation
-    if (length(n) > 2) warning("Extra dimensions of n ignored")
+    if (length(n) > 3) warning("Extra dimensions of n ignored")
     if (length(n) == 1) n <- rep(n, 2)
     if (class(mx) != "matrix") mx <- as.matrix(mx)
 
@@ -23,6 +24,9 @@ repmat <- function (mx, n) {
         for (i in seq(n[1] - 1)) out <- rbind(out, mx_col)
     }
 
+    # Replicating 3rd dimension
+    if (!is.na(n[3]) & n[3] > 1) out <- array(out, c(dim(out), n[3]))
+
     # Output
-    return(unname(as.matrix(out)))
+    return(unname(as.array(out)))
 }
\ No newline at end of file
diff --git a/R/reshape.R b/R/reshape.R
new file mode 100644
index 0000000..824bbae
--- /dev/null
+++ b/R/reshape.R
@@ -0,0 +1,24 @@
+#' @title Reshape array
+#' @description Reshapes a matrix according to a certain number of dimensions
+#' @param A input matrix
+#' @param sz vector containing the dimensions of the output vector
+#' @details This function replicates the functionality of the `reshape()`
+#' function on Matlab. This function is basically a fancy wrapper for the
+#' `array()` function in R, but is useful because it saves the user translation
+#' time. Moreover, it introduces validation code that alter the behavior of
+#' `array()` and makes it more similar to `replicate()`.
+#' @note The Matlab function also accepts as input the dismemberment of sz as
+#' scalars.
+reshape <- function(A, sz) {
+    # Validation
+    if (prod(sz) != prod(dim(A))) {
+        stop("To RESHAPE the number of elements must not change.")
+    }
+    if (length(sz) == 1) {
+        stop("Size vector must have at least two elements.")
+    }
+
+    # Reshaping A
+    A <- array(A, sz)
+    return(A)
+}
\ No newline at end of file
diff --git a/man/computeAllFreqs2.Rd b/man/computeAllFreqs2.Rd
new file mode 100644
index 0000000..7bb8b04
--- /dev/null
+++ b/man/computeAllFreqs2.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/computeAllFreqs2.R
+\name{computeAllFreqs2}
+\alias{computeAllFreqs2}
+\title{Compute all freqs - version 2}
+\usage{
+computeAllFreqs2(noalle, COUNTS = matrix(NA, 0, 0), SUMCOUNTS = sum(COUNTS))
+}
+\arguments{
+\item{noalle}{noalle}
+
+\item{COUNTS}{counts}
+
+\item{SUMCOUNTS}{sumcounts}
+}
+\description{
+Lisää a priori jokaista alleelia joka populaation joka lokukseen
+j 1/noalle(j) verran.
+}
diff --git a/man/reshape.Rd b/man/reshape.Rd
new file mode 100644
index 0000000..39d74f9
--- /dev/null
+++ b/man/reshape.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/reshape.R
+\name{reshape}
+\alias{reshape}
+\title{Reshape array}
+\usage{
+reshape(A, sz)
+}
+\arguments{
+\item{A}{input matrix}
+
+\item{sz}{vector containing the dimensions of the output vector}
+}
+\description{
+Reshapes a matrix according to a certain number of dimensions
+}
+\details{
+This function replicates the functionality of the `reshape()`
+function on Matlab. This function is basically a fancy wrapper for the
+`array()` function in R, but is useful because it saves the user translation
+time. Moreover, it introduces validation code that alter the behavior of
+`array()` and makes it more similar to `replicate()`.
+}
+\note{
+The Matlab function also accepts as input the dismemberment of sz as
+scalars.
+}
diff --git a/man/simulateAllFreqs.Rd b/man/simulateAllFreqs.Rd
index 97bf0e6..9a0682b 100644
--- a/man/simulateAllFreqs.Rd
+++ b/man/simulateAllFreqs.Rd
@@ -4,7 +4,7 @@
 \alias{simulateAllFreqs}
 \title{Simulate All Frequencies}
 \usage{
-simulateAllFreqs(noalle, COUNTS = matrix())
+simulateAllFreqs(noalle, COUNTS = matrix(NA, 0, 0))
 }
 \arguments{
 \item{noalle}{noalle}
diff --git a/tests/testthat/test-admix1.R b/tests/testthat/test-admix1.R
index b5cd9e3..12ae460 100644
--- a/tests/testthat/test-admix1.R
+++ b/tests/testthat/test-admix1.R
@@ -114,7 +114,7 @@ test_that("computeIndLogml works like on Matlab", {
     expect_equivalent(computeIndLogml(1, 0), -Inf)
     expect_equivalent(computeIndLogml(0, 0), -Inf)
     expect_equivalent(computeIndLogml(-pi, -8), 3.2242, tol = .0001)
-    expect_equivalent(computeIndLogml(2:3, 2), 2.3026, tol = .0001)    
+    expect_equivalent(computeIndLogml(2:3, 2), 2.3026, tol = .0001)
     expect_equivalent(computeIndLogml(matrix(8:5, 2), 100), 14.316, tol = .001)
     expect_equivalent(
         object    = computeIndLogml(matrix(8:5, 2), matrix(c(1, 3), 1)),
@@ -223,4 +223,8 @@ test_that("simulateAllFreqs works as expected", {
         object = suppressWarnings(simulateAllFreqs(matrix(1:4, 2))),
         expected = empty_mt
     )
+})
+
+test_that("computeAllFreqs2 works as expected", {
+    expect_equivalent(computeAllFreqs2(10), matrix(NA, 0, 0))
 })
\ No newline at end of file
diff --git a/tests/testthat/test-convertedBaseFunctions.R b/tests/testthat/test-convertedBaseFunctions.R
index f4bf346..7f75050 100644
--- a/tests/testthat/test-convertedBaseFunctions.R
+++ b/tests/testthat/test-convertedBaseFunctions.R
@@ -28,6 +28,10 @@ test_that("repmat works properly", {
         object = repmat(mx2, c(4, 1)),
         expected = rbind(mx2, mx2, mx2, mx2)
     )
+    expect_equal(
+        object = repmat(mx2, c(1, 1, 2)),
+        expected = array(mx2, c(2, 2, 2))
+    )
 })
 
 test_that("zeros and ones work as expected", {
@@ -84,4 +88,15 @@ test_that("size works as on MATLAB", {
     expect_equal(size(mx, 2), 3)
     expect_equal(size(ra, 2), 3)
     expect_equal(size(ra, 3), 4)
+})
+
+test_that("reshape reshapes properly", {
+    mx <- matrix(1:4, 2)
+    ra <- array(1:12, c(2, 3, 2))
+    expect_equal(reshape(mx, c(1, 4)), matrix(1:4, 1))
+    expect_equal(reshape(mx, c(2, 2)), mx)
+    expect_equal(reshape(mx, c(1, 1, 4)), array(mx, c(1, 1, 4)))
+    expect_error(reshape(mx, c(1, 2, 3)))
+    expect_error(reshape(ra, c(1, 2, 3)))
+    expect_equal(reshape(ra, c(3, 2, 2)), array(ra, c(3, 2, 2)))
 })
\ No newline at end of file