diff --git a/README.md b/README.md
index ea1aa8ea8..60f9bb307 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,7 @@ Authors:
 ------------
 
 ####  Install Rcpp package  
- 
+
 * Install package-dependencies: `Rcpp`, `RcppEigen`, `BH`.  
 
 1. Then use devtools package to install `volesti` Rcpp package. In folder `/root/R-prog` Run:
@@ -116,7 +116,6 @@ The function prints the volume.
 `
 ./vol -file <filename> -sample
 `
-  
 - The default settings are: `100` uniformly distributed points from the uniform distribution using billiard walk with walk length `1`.  
 - You can use the following flags: i) `-walk_length <walk_length>` to set the walk length of the random walk, ii) `-N <number_of_points>` to set the number of points to sample, iii) `-boltz` to sample from the Boltzmann distribution, iv) `-rdhr` to sample with random directions Hit and Run, `-cdhr` to sample with coordinate directions Hit and Run, `-hmc` to sample with Hamiltonian Monte Carlo with reflections, `-temperature <variance_of_boltzmann_distribution>` to set the variance of the Boltzmann distribution.  
 
@@ -130,10 +129,16 @@ The function prints the sampled points.
 ./vol -file <filename> -sdp
 `
 - The default settings are: `20` iterations are performed, with HMC sampling with walk length equal to `1`.  
+
 - You can use the following flags: i) `-N <number_of_iterations>` to set the number of iterations, ii) `-walk_length  <walk_length>` to set the walk length of the random walk iii) `-rdhr` to sample with random directions Hit and Run, `-hmc` to sample with Hamiltonian Monte Carlo with reflections.  
 
 - Example:  
-`./generate -n 10 -m 16`  
-`./vol -file sdp_prob_10_16.txt -sdp -N 30 -hmc -walk_length 3`  
-The function prints the values of the objective function of each iteration.
+  `./generate -n 10 -m 16`  
+  `./vol -file sdp_prob_10_16.txt -sdp -N 30 -hmc -walk_length 3`  
+  The function prints the values of the objective function of each iteration.
+
+  
+
+## Volume estimation of correlation matrices (matlab code)
 
+In folder `matlab_code` we include the matlab code to estimate the volume of correlation matrices. To use the main function `volume` use the script `volume_example`. 
diff --git a/matlab_code/utils/get_billiard_L.m b/matlab_code/utils/get_billiard_L.m
new file mode 100644
index 000000000..06db55693
--- /dev/null
+++ b/matlab_code/utils/get_billiard_L.m
@@ -0,0 +1,59 @@
+function L = get_billiard_L(m, N)
+
+    n = m^2 / 2 - m/2;
+    
+    [upper, lower] = initialize_sampler(m);
+    x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+    
+    bc = ones(2 * n, 1);
+    L = 0;
+    
+    for j = 1:N
+            v = get_direction(n);
+            
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+                
+            B(lower) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            eigenvalues = eig(B, -A);
+            l_max = 1 / max(eigenvalues);
+            l_min = 1 / min(eigenvalues);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            l_max_temp = 1 / max(lambdas);
+            l_min_temp = 1 / min(lambdas);
+                
+            % decide which boundary the line hits first
+            if (l_max_temp < l_max)
+                l_max = l_max_temp;
+            end
+            if (l_min_temp > l_min)
+                l_min = l_min_temp;
+            end
+            
+            if (L < (l_max + l_min))
+                L = l_max + l_min;
+            end
+            
+            % pick a uniformly distributed point from the segment
+            %lambda = l_min + rand * (l_max - l_min);
+            
+            % update the current point of the random walk
+            %x = x + lambda * v;
+            
+    end
+    
+    L = 4 * L;
+
+end
diff --git a/matlab_code/utils/get_direction.m b/matlab_code/utils/get_direction.m
new file mode 100644
index 000000000..332eb065d
--- /dev/null
+++ b/matlab_code/utils/get_direction.m
@@ -0,0 +1,6 @@
+function [v] = get_direction(n)
+
+    v = randn(n,1);
+    v = v / norm(v);
+
+end
\ No newline at end of file
diff --git a/matlab_code/utils/get_gradient.m b/matlab_code/utils/get_gradient.m
new file mode 100644
index 000000000..7d69bd990
--- /dev/null
+++ b/matlab_code/utils/get_gradient.m
@@ -0,0 +1,7 @@
+function [s] = get_gradient(q)
+    
+    s = nchoosek(q, 2);
+    s = s(:,1) .* s(:,2);
+    s = s / norm(s);
+    
+end
\ No newline at end of file
diff --git a/matlab_code/utils/initialize_sampler.m b/matlab_code/utils/initialize_sampler.m
new file mode 100644
index 000000000..37bc4450e
--- /dev/null
+++ b/matlab_code/utils/initialize_sampler.m
@@ -0,0 +1,7 @@
+function [upper, lower] = initialize_sampler(m)
+    
+    At = zeros(m, m);
+    upper  = (1:size(At,1)).' < (1:size(At,2));
+    lower = (1:size(At,1)).' > (1:size(At,2));
+    
+end
\ No newline at end of file
diff --git a/matlab_code/utils/randsphere.m b/matlab_code/utils/randsphere.m
new file mode 100644
index 000000000..191c0b7cb
--- /dev/null
+++ b/matlab_code/utils/randsphere.m
@@ -0,0 +1,18 @@
+function X = randsphere(m,n,r)
+ 
+% This function returns an m by n array, X, in which 
+% each of the m rows has the n Cartesian coordinates 
+% of a random point uniformly-distributed over the 
+% interior of an n-dimensional hypersphere with 
+% radius r and center at the origin.  The function 
+% 'randn' is initially used to generate m sets of n 
+% random variables with independent multivariate 
+% normal distribution, with mean 0 and variance 1.
+% Then the incomplete gamma function, 'gammainc', 
+% is used to map these points radially to fit in the 
+% hypersphere of finite radius r with a uniform % spatial distribution.
+% Roger Stafford - 12/23/05
+ 
+X = randn(m,n);
+s2 = sum(X.^2,2);
+X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);
\ No newline at end of file
diff --git a/matlab_code/volume/billiard_walk_intersection.m b/matlab_code/volume/billiard_walk_intersection.m
new file mode 100644
index 000000000..19a66bc17
--- /dev/null
+++ b/matlab_code/volume/billiard_walk_intersection.m
@@ -0,0 +1,127 @@
+function [points] = billiard_walk_intersection(m, J, L, R, N)
+
+    n = sum(sum(J>0));
+    points = zeros(n, N);
+    
+    [upper, ~] = initialize_sampler_vol(m);
+    %x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+    
+    bc = ones(2 * n, 1);
+    pos = 1:(2*n);
+    pos = mod(pos, n);
+    pos(pos==0) = n;
+    
+    for i = 1:N
+        
+        for jj = 1:3
+        
+        T = -log(rand) * L;
+        v = get_direction(n);
+            
+        rho = 0;
+        
+        while (true)
+            
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+                
+            B(J) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            [Q, eigenvalues] = eig(B, -A);
+            [max_eig, pos_max_eig] = max(diag(eigenvalues));
+            l_max = 1 / max_eig;
+            
+            x = A(J);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            [l_max_temp, pos_max] = max(lambdas);
+            l_max_temp = 1 / l_max_temp;
+            
+            vrc = v'*x;
+            v2 = v'*v;
+            rc2 = x'*x;
+
+            disc_sqrt = sqrt(vrc^2 - v2 * (rc2 - R^2));
+            
+            lR_max = max((-vrc + disc_sqrt)/v2, (-vrc - disc_sqrt)/v2);
+            
+            [l_max, lmax_ind] = min([l_max l_max_temp lR_max]);
+                
+            % decide which boundary the line hits first
+            if (lmax_ind == 2)
+                %l_max = l_max_temp;
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+            
+                % update the current point of the random walk
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T * B;
+                    break;
+                end
+                A = A + lambda * B;
+                %x = x + lambda * v;
+                    
+                %reflevt the ray
+                %v = v - (2*AA(pos_max, :)*v)*AA(pos_max, :)'
+                p = pos(pos_max);
+                v(p) = -v(p);
+            elseif (lmax_ind == 1)
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                s = get_gradient(Q(:, pos_max_eig));
+                % update the current point of the random walk
+                %x = x + lambda * v;
+                A = A + lambda * B;
+                %reflect the ray
+                v = v - (2*(v'*s))*s;
+            else
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                A = A + lambda * B;
+                x = x + lambda * v;
+                
+                x = x / norm(x);
+                v = v - (2 * (v' * x)) * x;
+            end
+            rho = rho + 1;
+            T = T - lambda;
+        end
+        end
+        %A(lower) = x;
+        %q = triu(A',1);
+        %A(upper) = q(upper);
+        x = A(J);
+        points(:, i) = x;
+        
+        %is_in = false;
+        %if (norm(x) <= R2)
+        %    is_in = true;
+        %end
+        
+        %conv = update_window(ratio_parameters, is_in);
+        
+    end
+    
+    %ratio = ratio_parameters.count_in / ratio_parameters.tot_count;
+    
+end
\ No newline at end of file
diff --git a/matlab_code/volume/check_convergence_for_ball.m b/matlab_code/volume/check_convergence_for_ball.m
new file mode 100644
index 000000000..52868c9f6
--- /dev/null
+++ b/matlab_code/volume/check_convergence_for_ball.m
@@ -0,0 +1,46 @@
+function [conv, ratio, too_few] = check_convergence_for_ball(nu, X, r0, lb, ub, lastball)
+
+    conv = false;
+    too_few = false;
+    
+    N = size(X, 2);
+    
+    mm = N / nu;
+    countsIn = 0;
+    ratios = [];
+    
+    for i=1:N
+        
+        x = X(:,i);
+        
+        if (norm(x) <= r0)
+            countsIn = countsIn + 1;
+        end
+        
+        if (mod(i, mm) == 0)
+            ratios = [ratios countsIn/mm];
+            countsIn = 0;
+        end
+    end
+    
+    ratio = mean(ratios);
+    t_a = tinv(0.90, nu);
+    rs = std(ratios);
+    T = rs * (t_a / sqrt(nu));
+    
+    if (ratio > lb + T)
+        if (lastball)
+            conv = true;
+            return
+        end
+        
+        if (ratio < ub + T)
+            conv = true;
+            return
+        end
+        return
+    end
+    too_few = true;
+
+end
+
diff --git a/matlab_code/volume/check_convergence_in_spectra.m b/matlab_code/volume/check_convergence_in_spectra.m
new file mode 100644
index 000000000..774a345c4
--- /dev/null
+++ b/matlab_code/volume/check_convergence_in_spectra.m
@@ -0,0 +1,70 @@
+function [conv, ratio, too_few] = check_convergence_in_spectra(m, J, X, nu, lb, ub)
+
+    too_few = false;
+    conv = false;
+    
+    %n = sum(sum(J>0));
+    
+    [upper, ~] = initialize_sampler_vol(m);
+    A = eye(m);
+    
+    N = size(X, 2);
+    mm = N / nu;
+    countsIn = 0;
+    ratios = [];
+    
+    for i=1:N
+        
+        A(J) = X(:,i);
+        q = triu(A',1);
+        A(upper) = q(upper);
+        
+        if (min(eig(A)) > 0)
+            countsIn = countsIn + 1;
+        end
+        
+        if (mod(i, mm) == 0)
+            ratios = [ratios countsIn/mm];
+            countsIn = 0;
+            if (length(ratios) > 1)
+                
+                ratio = mean(ratios);
+                t_a = tinv(0.995, nu);
+                rs = std(ratios);
+                T = rs * (t_a / sqrt(length(ratios)));
+                
+            
+                %[h, ~] = ttest(ratios, lb, 'Alpha', 0.005, 'Tail', 'right');
+                if (ratio + T < lb)
+                    too_few = true;
+                    conv = false;
+                    return
+                end
+            
+                %[h, ~] = ttest(ratios, ub, 'Alpha', 0.005, 'Tail', 'left');
+                if (ratio - T > ub)
+                    conv = false;
+                    return
+                end
+            end
+        end
+    end
+    
+    ratio = mean(ratios);
+    t_a = tinv(0.95,nu);
+    rs = std(ratios);
+    T = rs * (t_a / sqrt(nu));
+    
+    if (ratio > lb + T)
+        if (ratio < ub - T)
+            conv = true;
+            return
+        end
+        conv = false;
+        return
+    end
+    too_few = true;
+    conv = false;
+
+end
+
diff --git a/matlab_code/volume/estimate_ratio.m b/matlab_code/volume/estimate_ratio.m
new file mode 100644
index 000000000..9370724a4
--- /dev/null
+++ b/matlab_code/volume/estimate_ratio.m
@@ -0,0 +1,129 @@
+function [ratio, points, numpoints] = estimate_ratio(m, J, L, R1, R2, p, ratio, error, W, N_nu)
+
+    n = sum(sum(J>0));
+    points = p;
+    
+    [upper, ~] = initialize_sampler(m);
+    %x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+    
+    bc = ones(2 * n, 1);
+    pos = 1:(2*n);
+    pos = mod(pos, n);
+    pos(pos==0) = n;
+    
+    conv = false;
+    ratio_parameters = initialize_parameters(W, N_nu, ratio);
+    
+    while (~conv)
+        
+        for jj=1:3
+        T = -log(rand) * L;
+        v = get_direction(n);
+            
+        rho = 0;
+        
+        while (true)
+            
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+                
+            B(J) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            [Q, eigenvalues] = eig(B, -A);
+            [max_eig, pos_max_eig] = max(diag(eigenvalues));
+            l_max = 1 / max_eig;
+            
+            x = A(J);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            [l_max_temp, pos_max] = max(lambdas);
+            l_max_temp = 1 / l_max_temp;
+            
+            vrc = v'*x;%.dot(r);
+            v2 = v'*v;%.dot(v);
+            rc2 = x'*x;%.dot(r);
+
+            disc_sqrt = sqrt(vrc^2 - v2 * (rc2 - R1^2));
+            
+            lR_max = max((-vrc + disc_sqrt)/v2, (-vrc - disc_sqrt)/v2);
+            
+            [l_max, lmax_ind] = min([l_max l_max_temp lR_max]);
+                
+            % decide which boundary the line hits first
+            if (lmax_ind == 2)
+                %l_max = l_max_temp;
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+            
+                % update the current point of the random walk
+                if (T <= lambda)
+                    %x = x + T * v;
+                    A = A + T * B;
+                    break;
+                end
+                A = A + lambda * B;
+                %x = x + lambda * v;
+                    
+                %reflevt the ray
+                %v = v - (2*AA(pos_max, :)*v)*AA(pos_max, :)'
+                p = pos(pos_max);
+                v(p) = -v(p);
+            elseif (lmax_ind == 1)
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                s = get_gradient(Q(:, pos_max_eig));
+                % update the current point of the random walk
+                %x = x + lambda * v;
+                A = A + lambda * B;
+                %reflect the ray
+                v = v - (2*(v'*s))*s;
+            else
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T * B;
+                    break;
+                end
+                A = A + lambda * B;
+                x = x + lambda * v;
+                
+                x = x / norm(x);
+                v = v - (2 * (v' * x)) * x;
+            end
+            rho = rho + 1;
+            T = T - lambda;
+        end
+        end
+        %A(lower) = x;
+        %q = triu(A',1);
+        %A(upper) = q(upper);
+        x = A(J);
+        
+        is_in = false;
+        if (norm(x) <= R2)
+            is_in = true;
+        end
+        
+        [conv, ratio, ratio_parameters] = update_window(ratio_parameters, is_in, error);
+        
+    end
+    numpoints = ratio_parameters.tot_count;
+    %ratio = ratio_parameters.count_in / ratio_parameters.tot_count;
+    
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/estimate_ratio_first.m b/matlab_code/volume/estimate_ratio_first.m
new file mode 100644
index 000000000..c8e6deb8b
--- /dev/null
+++ b/matlab_code/volume/estimate_ratio_first.m
@@ -0,0 +1,30 @@
+function ratio = estimate_ratio_first(m, J, R, ratio, error, W, N_nu)
+
+    n = sum(sum(J>0));
+    
+    [upper, ~] = initialize_sampler(m);
+    A = eye(m);
+
+    conv = false;
+    
+    ratio_parameters = initialize_parameters(W, N_nu, ratio);
+    
+    while (~conv)
+        
+        x = randsphere(1, n, R)';
+        A(J) = x;
+        q = triu(A',1);
+        A(upper) = q(upper);
+        
+        is_in = false;
+        if (min(eig(A)) > 0)
+            is_in = true;
+        end
+        
+        [conv, ratio, ratio_parameters] = update_window(ratio_parameters, is_in, error);
+       
+    end
+    
+    %ratio = ratio_parameters.count_in / ratio_parameters.tot_count;
+    
+end
\ No newline at end of file
diff --git a/matlab_code/volume/estimate_ratio_spectra_ball.m b/matlab_code/volume/estimate_ratio_spectra_ball.m
new file mode 100644
index 000000000..0f8407c8e
--- /dev/null
+++ b/matlab_code/volume/estimate_ratio_spectra_ball.m
@@ -0,0 +1,111 @@
+function [ratio, points, numpoints] = estimate_ratio_spectra_ball(m, J, L, R, p, ratio, error, W, N_nu)
+
+    n = sum(sum(J>0));
+    points = p;
+    
+    [upper, ~] = initialize_sampler(m);
+    %x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+    
+    bc = ones(2 * n, 1);
+    pos = 1:(2*n);
+    pos = mod(pos, n);
+    pos(pos==0) = n;
+    
+    conv = false;
+    ratio_parameters = initialize_parameters(W, N_nu, ratio);
+    
+    while (~conv)
+        
+        for jj=1:3
+        T = -log(rand) * L;
+        v = get_direction(n);
+            
+        rho = 0;
+        
+        while (true)
+            
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+                
+            B(J) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            [Q, eigenvalues] = eig(B, -A);
+            [max_eig, pos_max_eig] = max(diag(eigenvalues));
+            l_max = 1 / max_eig;
+            
+            x = A(J);
+            
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            [l_max_temp, pos_max] = max(lambdas);
+            l_max_temp = 1 / l_max_temp;
+            
+            
+            
+            [l_max, lmax_ind] = min([l_max l_max_temp]);
+                
+            % decide which boundary the line hits first
+            if (lmax_ind == 2)
+                %l_max = l_max_temp;
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+            
+                % update the current point of the random walk
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T * B;
+                    break;
+                end
+                A = A + lambda * B;
+                %x = x + lambda * v;
+                    
+                %reflevt the ray
+                %v = v - (2*AA(pos_max, :)*v)*AA(pos_max, :)'
+                p = pos(pos_max);
+                v(p) = -v(p);
+            elseif (lmax_ind == 1)
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                s = get_gradient(Q(:, pos_max_eig));
+                % update the current point of the random walk
+                %x = x + lambda * v;
+                A = A + lambda * B;
+                %reflect the ray
+                v = v - (2*(v'*s))*s;
+            end
+            rho = rho + 1;
+            T = T - lambda;
+        end
+        end
+        %A(lower) = x;
+        %q = triu(A',1);
+        %A(upper) = q(upper);
+        x = A(J);
+        
+        is_in = false;
+        if (norm(x) <= R)
+            is_in = true;
+        end
+        
+        [conv, ratio, ratio_parameters] = update_window(ratio_parameters, is_in, error);
+        
+    end
+    numpoints = ratio_parameters.tot_count;
+    %ratio = ratio_parameters.count_in / ratio_parameters.tot_count;
+    
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/exact_volume.m b/matlab_code/volume/exact_volume.m
new file mode 100644
index 000000000..8d6ff6223
--- /dev/null
+++ b/matlab_code/volume/exact_volume.m
@@ -0,0 +1,12 @@
+function log_vol = exact_volume(m)
+    
+    a = (m+1)/2;
+    log_vol = -m * log_gamma_function((m+1)/2);
+
+    log_vol = log_vol + ((m*(m-1))/4)*log(pi);
+    
+    for i=1:m
+        log_vol = log_vol + log_gamma_function(a - ((i-1)/2));
+    end
+end
+
diff --git a/matlab_code/volume/get_L.m b/matlab_code/volume/get_L.m
new file mode 100644
index 000000000..a8ace57a7
--- /dev/null
+++ b/matlab_code/volume/get_L.m
@@ -0,0 +1,10 @@
+function L = get_L(m, J, radius)
+
+    n = m^2 / 2 - m/2;
+    L = 5 * sqrt(n) * radius;
+    
+    points = billiard_walk_intersection(m, J, L, 10^10, 10 * (ceil(log(n))*10));
+    
+    L = max([L, 2 * max(sqrt(sum(points.^2, 1)))]);
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/get_inscribed_radius.m b/matlab_code/volume/get_inscribed_radius.m
new file mode 100644
index 000000000..2b1f31989
--- /dev/null
+++ b/matlab_code/volume/get_inscribed_radius.m
@@ -0,0 +1,38 @@
+function radius = get_inscribed_radius(m, J)
+
+    n = sum(sum(J));
+    
+    [upper, ~] = initialize_sampler_vol(m);
+    A = eye(m);
+    B = zeros(m);
+    
+    radius = Inf;
+    
+    for i=1:n
+       
+        v = zeros(n,1);
+        v(i) = 1;
+        
+        B(J) = v;
+        q = triu(B',1);
+        B(upper) = q(upper);
+        
+        [~, eigenvalues] = eig(B, -A);
+        max_eig = max(diag(eigenvalues));
+        l_pos = 1 / max_eig;
+        
+        [~, eigenvalues] = eig(-B, -A);
+        max_eig = max(diag(eigenvalues));
+        l_neg = 1 / max_eig;
+        
+        ell = min([l_pos, l_neg]);
+        
+        if (ell < radius)
+            radius = ell;
+        end
+        
+    end
+    
+    radius = radius / sqrt(n);
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/get_next_ball.m b/matlab_code/volume/get_next_ball.m
new file mode 100644
index 000000000..92fd72a94
--- /dev/null
+++ b/matlab_code/volume/get_next_ball.m
@@ -0,0 +1,33 @@
+function [rad, ratio] = get_next_ball(X, nu, rmax, r0, lb, ub)
+
+    rad0 = r0;
+    rad_m = rmax;
+    tolerance = 0.00000000001;
+   
+    while (true)
+        
+        rad_med = 0.5 * (r0 + rmax);
+        
+        [conv, ratio, too_few] = check_convergence_for_ball(nu, X, rad_med, lb, ub, false);
+    
+        if (conv)
+            rad = rad_med;
+            return
+        end
+        
+        if (too_few)
+            r0 = rad_med;
+        else
+            rmax = rad_med;
+        end
+        
+        if (rmax - r0 < tolerance)
+            r0 = rad0;
+            rmax = rad_m;
+        end
+        
+    end
+
+end
+
+
diff --git a/matlab_code/volume/get_sequence_of_balls.m b/matlab_code/volume/get_sequence_of_balls.m
new file mode 100644
index 000000000..de3eb5e00
--- /dev/null
+++ b/matlab_code/volume/get_sequence_of_balls.m
@@ -0,0 +1,49 @@
+function [Rs, ratios] = get_sequence_of_balls(m, J, L, radius, N, nu, lb, ub)
+    
+    Rs = [];
+    ratios = [];
+    n = sum(sum(J>0));
+
+    [r0, ratio0] = get_smallest_ball(m, J, radius, 1200, nu, lb, ub); %ready
+    %r0
+    %ratio0
+    
+    %n = sum(sum(J>0));
+    Ntot = N * nu;
+    
+    rmax = Inf;
+    rad = Inf;
+    while (true)
+        
+        if (isinf(rmax))
+            [points] = billiard_walk_intersection(m, J, L, 4*sqrt(n), Ntot);
+        else
+            L2 = min(L, rmax);
+            [points] = billiard_walk_intersection(m, J, L2, rmax, Ntot);
+        end
+        
+        [conv, ratio] = check_convergence_for_ball(nu, points, r0, lb, ub, true); %ready
+        
+        if (conv)
+            ratios = [ratios ratio];
+            Rs = [Rs r0];
+            ratios = [ratios ratio0];
+            return
+        end
+        
+        if (isinf(rmax))
+            rmax = max(sqrt(sum(points.^2,1)));
+        else
+            rmax = rad;
+        end
+        
+        [rad, ratio] = get_next_ball(points, nu, rmax, r0, lb, ub); %ready
+        Rs = [Rs rad];
+        ratios = [ratios ratio];
+        rmax = rad;
+    end
+    
+    %Rs = [Rs r0];
+    %ratios = [ratios ratio0];
+    
+end
\ No newline at end of file
diff --git a/matlab_code/volume/get_smallest_ball.m b/matlab_code/volume/get_smallest_ball.m
new file mode 100644
index 000000000..21ee38bad
--- /dev/null
+++ b/matlab_code/volume/get_smallest_ball.m
@@ -0,0 +1,57 @@
+function [r, ratio] = get_smallest_ball(m, J, radius, N, nu, lb, ub)
+
+    n = sum(sum(J>0));
+    tolerance = 0.00000000001;
+    
+    sqrt_n = sqrt(n);
+    
+    rad1 = radius;
+    rmax = 2 * sqrt_n * radius;
+    conv = false;
+    
+    while (~conv)
+      
+        X = randsphere(N, n, rmax)';
+        %size(X)
+        
+        [conv, ratio, too_few] = check_convergence_in_spectra(m, J, X, nu, lb, ub);
+
+        if (conv)
+            r = rmax;
+            return
+        end
+        if (too_few)
+            break
+        end
+        rad1 = rmax;
+        rmax = rmax + 2 * sqrt_n * radius;
+    end
+    
+    rad0 = rad1;
+    rad_m = rmax;
+    while (true)
+        
+        rad_med = 0.5 * (rad1 + rmax);
+        
+        X = randsphere(N, n, rad_med)';
+        [conv, ratio, too_few] = check_convergence_in_spectra(m, J, X, nu, lb, ub);
+        
+        if (conv)
+            r = rad_med;
+            return
+        end
+        
+        if (too_few)
+            rmax = rad_med;
+        else
+            rad1 = rad_med;
+        end
+        
+        if (rmax - rad1 < tolerance)
+            rad1 = rad0;
+            rmax = rad_m;
+        end
+        
+    end
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/initialize_parameters.m b/matlab_code/volume/initialize_parameters.m
new file mode 100644
index 000000000..b2f3de5b8
--- /dev/null
+++ b/matlab_code/volume/initialize_parameters.m
@@ -0,0 +1,17 @@
+function parameters = initialize_parameters(W, N, ratio)
+
+    parameters = struct;
+    
+    parameters.min_val = -Inf;
+    parameters.max_val = Inf;
+    parameters.min_index = W;
+    parameters.max_index = W;
+    parameters.W = W;
+    parameters.index =1;
+    parameters.tot_count = N;
+    parameters.count_in = floor(N * ratio);
+    %parameters.tot_count = 0;
+    %parameters.count_in = 0;
+    parameters.last_W = zeros(1, W);
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/initialize_sampler_vol.m b/matlab_code/volume/initialize_sampler_vol.m
new file mode 100644
index 000000000..00232898c
--- /dev/null
+++ b/matlab_code/volume/initialize_sampler_vol.m
@@ -0,0 +1,7 @@
+function [upper, lower] = initialize_sampler_vol(m)
+    
+    At = zeros(m, m);
+    upper  = (1:size(At,1)).' < (1:size(At,2));
+    lower = (1:size(At,1)).' > (1:size(At,2));
+    
+end
\ No newline at end of file
diff --git a/matlab_code/volume/log_gamma_function.m b/matlab_code/volume/log_gamma_function.m
new file mode 100644
index 000000000..6a7d52fd3
--- /dev/null
+++ b/matlab_code/volume/log_gamma_function.m
@@ -0,0 +1,9 @@
+function y = log_gamma_function(x) 
+
+    if (x <= 100)
+        y = log(gamma(x));
+        return
+    end
+    
+    y = log(x - 1) + log_gamma_function(x - 1);
+end
\ No newline at end of file
diff --git a/matlab_code/volume/randsphere.m b/matlab_code/volume/randsphere.m
new file mode 100644
index 000000000..191c0b7cb
--- /dev/null
+++ b/matlab_code/volume/randsphere.m
@@ -0,0 +1,18 @@
+function X = randsphere(m,n,r)
+ 
+% This function returns an m by n array, X, in which 
+% each of the m rows has the n Cartesian coordinates 
+% of a random point uniformly-distributed over the 
+% interior of an n-dimensional hypersphere with 
+% radius r and center at the origin.  The function 
+% 'randn' is initially used to generate m sets of n 
+% random variables with independent multivariate 
+% normal distribution, with mean 0 and variance 1.
+% Then the incomplete gamma function, 'gammainc', 
+% is used to map these points radially to fit in the 
+% hypersphere of finite radius r with a uniform % spatial distribution.
+% Roger Stafford - 12/23/05
+ 
+X = randn(m,n);
+s2 = sum(X.^2,2);
+X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);
\ No newline at end of file
diff --git a/matlab_code/volume/run_experiments.m b/matlab_code/volume/run_experiments.m
new file mode 100644
index 000000000..61044e32c
--- /dev/null
+++ b/matlab_code/volume/run_experiments.m
@@ -0,0 +1,26 @@
+mat_sizes = [8 10 12];
+errors = zeros(3, 10);
+runtimes = zeros(3, 10);
+log_vols = zeros(3, 10);
+total_points = zeros(3,10);
+
+for j = 1:length(mat_sizes)
+    
+    m = mat_sizes(j);
+    [upper, lower] = initialize_sampler(m);
+    J = lower;
+
+    for i=1:10
+        disp([m i])
+        tic
+        [log_vol, Rs, ratios_out, tot_points] = volume(m, 0.1, 1000*j, J);
+        tim = toc;
+        log_vols(j, i) = log_vol;
+        runtimes(j, i) = tim;
+        total_points(j, i) = tot_points;
+        errors(j, i) = abs(exp(log_vol) - exp(exact_volume(m)))/ exp(exact_volume(m))
+    end
+    
+end
+save('elliptope_volume_results.mat', 'log_vols', 'runtimes', 'errors')
+
diff --git a/matlab_code/volume/update_window.m b/matlab_code/volume/update_window.m
new file mode 100644
index 000000000..681806f90
--- /dev/null
+++ b/matlab_code/volume/update_window.m
@@ -0,0 +1,60 @@
+function [conv, val, ratio_parameters] = update_window(ratio_parameters, is_in, error)
+
+    conv = false;
+    
+    if (is_in)
+        ratio_parameters.count_in = ratio_parameters.count_in + 1;
+    end
+    
+    ratio_parameters.tot_count = ratio_parameters.tot_count + 1;
+    val = ratio_parameters.count_in / ratio_parameters.tot_count;
+    ratio_parameters.last_W(ratio_parameters.index) = val;
+
+    if (val <= ratio_parameters.min_val)
+        ratio_parameters.min_val = val;
+        ratio_parameters.min_index = ratio_parameters.index;
+    elseif (ratio_parameters.min_index == ratio_parameters.index)
+        [ratio_parameters.min_val, ratio_parameters.min_index] = min(ratio_parameters.last_W);
+        %ratio_parameters.minmaxIt = std::min_element(ratio_parameters.last_W.begin(), ratio_parameters.last_W.end());
+        %ratio_parameters.min_val = (*ratio_parameters.minmaxIt);
+        %ratio_parameters.min_index = std::distance(ratio_parameters.last_W.begin(), ratio_parameters.minmaxIt);
+    end
+
+    if (val >= ratio_parameters.max_val)
+        ratio_parameters.max_val = val;
+        ratio_parameters.max_index = ratio_parameters.index;
+    elseif (ratio_parameters.max_index == ratio_parameters.index)
+        [ratio_parameters.max_val, ratio_parameters.max_index] = max(ratio_parameters.last_W);
+        %ratio_parameters.minmaxIt = std::max_element(ratio_parameters.last_W.begin(), ratio_parameters.last_W.end());
+        %ratio_parameters.max_val = (*ratio_parameters.minmaxIt);
+        %ratio_parameters.max_index = std::distance(ratio_parameters.last_W.begin(), ratio_parameters.minmaxIt);
+    end
+    
+    %ratio_parameters.last_W
+    %ratio_parameters.min_val
+    %ratio_parameters.max_val
+    
+    %ratio_parameters.last_W
+    %ratio_parameters.max_val
+    %ratio_parameters.max_index
+    
+    %ratio_parameters.min_val
+    %ratio_parameters.min_index
+    
+    %error
+    if ( (ratio_parameters.max_val - ratio_parameters.min_val) / ratio_parameters.max_val <= error/2 )
+        %(ratio_parameters.max_val - ratio_parameters.min_val) / ratio_parameters.max_val
+        %val
+        %ratio_parameters.tot_count
+        conv = true;
+        return
+    end
+       
+
+    ratio_parameters.index = mod(ratio_parameters.index, ratio_parameters.W) + 1;
+    %if (ratio_parameters.index == ratio_parameters.W + 1) 
+    %    ratio_parameters.index = 1;
+    %end
+    
+
+end
\ No newline at end of file
diff --git a/matlab_code/volume/volume.m b/matlab_code/volume/volume.m
new file mode 100644
index 000000000..f4f24a227
--- /dev/null
+++ b/matlab_code/volume/volume.m
@@ -0,0 +1,53 @@
+function [log_vol, Rs, ratios_out, total_points] = volume(m, e, W, J)
+    
+    n = sum(sum(J>0));
+    nu = 10;
+    N = 300;
+    N_nu = nu * N;
+    lb = 0.10;
+    ub = 0.15;
+    ratios_out = [];
+
+    radius = get_inscribed_radius(m, J);
+
+    L = get_billiard_L(m, 1000)/2;
+    
+    [Rs, ratios] = get_sequence_of_balls(m, J, L, radius, N, nu, lb, ub);
+    
+    log_vol = log(pi) * (n/2) + n * log(Rs(end)) - log_gamma_function(n / 2 + 1);
+    
+    mm = length(Rs) + 1;
+    er0 = e / (2.0 * sqrt(mm));
+    er1 = (e * sqrt(4.0 * mm - 1)) / (2.0 * sqrt(mm));
+    
+    ratio = estimate_ratio_first(m, J, Rs(end), ratios(end), er0, W, 1200); 
+    ratio_small = ratio;
+    log_vol = log_vol + log(ratio);
+    
+    er1 = er1 / sqrt(mm - 1.0);
+    total_points = 0;
+    
+    if (ratios(1) < 0.9999)
+        [ratio, points, numpoints] = estimate_ratio_spectra_ball(m, J, L, Rs(1), [], ratios(1), er1, W, N_nu);
+        total_points = total_points + numpoints;
+        ratios_out = [ratios_out ratio];
+        log_vol = log_vol - log(ratio);
+    end
+    
+    
+    for i=1:(mm-2)
+        
+        L2 = min(Rs(i), L);
+        [ratio, points, numpoints] = estimate_ratio(m, J, L2, Rs(i), Rs(i+1), points, ratios(i+1), er1, W, N_nu); %ready
+        %ratio
+        total_points = total_points + numpoints;
+        ratios_out = [ratios_out ratio];
+        log_vol = log_vol - log(ratio);
+        
+    end
+    ratios_out = [ratios_out ratio_small];
+
+end
+
+
+
diff --git a/matlab_code/volume_example.m b/matlab_code/volume_example.m
new file mode 100644
index 000000000..b8e5cfdc1
--- /dev/null
+++ b/matlab_code/volume_example.m
@@ -0,0 +1,19 @@
+addpath('./utils')
+addpath('./volume')
+
+%size of correlation matrices
+m=8;
+%initialize billiard walk sampler
+[upper, lower] = initialize_sampler(m);
+J = lower;
+%estimate the logarithm of the volume
+[log_vol, Rs, ratios_out, tot_points] = volume(m, 0.1, 1500, J);
+
+esti_vol = exp(log_vol);
+disp(strcat('estimated volume = ', num2str(esti_vol)));
+
+exact_vol = exp(exact_volume(m));
+disp(strcat('exact volume = ', num2str(exact_vol)));
+
+rel_error = abs(exp(log_vol) - exp(exact_volume(m)))/ exp(exact_volume(m));
+disp(strcat('relative error = ', num2str(rel_error)));