function ds = diffquot(x,r1,j,delta,m1,m2)
% This function calculates a difference quotient approximation of
% the derivative of the map f: R^(3^m1) -> R^(3^m2), which restricts a
% level m1 Y-circuit of resistors to a level m2 Y-circuit, with
% respect to a given edge
%
% support functions: restrict.m
%
% INPUT: 
% x - matrix of size 3^m1 x 3 holding resistances on a level m1
% Y network
% j - index of the edge to perturb
% r1 - restriction of x to level m2
% delta - delta value to use in difference quotient
% m1 - starting level of Y-network
% m2 - level of Y-network to restrict to

% 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
r2 = restrict(x+deltamatrix',m1,m2) ;
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