function ds = diffquotSeq(x,col,j,delta,m1,m2)
% This (yes, crudely commented) function computes a difference quotient
% using the original calculated sequence of 1-step restrictions 
% to save computation (compare diffquot.m). We accept a column
% matrix holding each intermediate restriction down to level m1
% Y-network to a level m2 Y-network, stacked on top of one
% another. Since we only perturb one edge, at each step of the new
% restriction (with edge number j perturbed) 2/3 of the restriction
% is the same, so we save this computation.
%
% support functions: smartrestrict.m
%
% INPUT: 
% x - matrix of size [3^m1 3] holding resistances on a level m1
% Y network
% col - all steps of the restriction of x to level m2
% j - index of the edge to perturb
% delta - delta value to use in difference quotient
% m1 - starting level of Y-network
% m2 - level of Y-network to restrict to
%
% OUTPUT: 
% 
% ds - approximation of derivative matrix (size [3^m2 3])

ds = [] ;

% create a vector of indices which demarcate the stored steps of the
% restriction of x to level m2
n = m1 - m2 ;
for i = 1:n+1
  breakpoints(i) = 3^(m1-i+1)*(3^i-1)/(3-1) ;
end

% create delta matrix for difference quotient
% (for proper ordering of entries we use its transpose)
deltamatrix = zeros(3,3^m1) ;
deltamatrix(j) = delta ;
 
% calculate difference quotient: (f(x+d) - f(x)) / ||d||_\infty

% pick out the level m1 cell (I think this is just x as I write
% comments a month later ...)
r1 = col(breakpoints(n)+1:breakpoints(n+1),:) ;
% calculate restriction of perturbed network to level m2
r2 = smartrestrict(x+deltamatrix',col,j,m1,m2,breakpoints) ;
% calculate difference quotient
r = (delta)^-1 * (r2 - r1) ;
r = r' ;

% we return data as a vector due to later modifications in how we
% store data on SG
for k = 1:3^(m2+1)
  ds(k) = r(k) ;
end