function phylosketch(m0,n0)
global m n
randomize
m = m0;
n = n0;

global pstr
pstr = mfilename
global lstr rstr

global maxr
maxr = 0;

global A C T P Pleaf leafs ll

while 1
  
  %T = randdtree(n);
  T=zeros(n); T(1,2:n)=ones(1,n-1);
  %T = bintree(2);
  %n = size(T,1);
  [nodesi,nodesj]=find(T);
  A=[];
  A0 = pnormdim(rand(m,m),1);
  tha = strict(A0);
  for e=1:length(nodesi)
    e1 = nodesi(e);
    e2 = nodesj(e);
    Ae = pnormdim(rand(m,m),1);
    Ae = reshape(Ae,1,1,m,m);
    A(e1,e2,:,:) = Ae;
    %A(e1,e2,:,:) = A0;
  end
  C = pnormdim(rand(m,1),1);
  P = dtreefillprob1(m,n,T,A,C);
  [hh0,Ht,H] = gethhn(P,m,n);
  
  
  leafs = leaves(T);
  [w,levels,ll] = treewidth(T);
  
  Pleaf = getmargX(P,m,n,leafs);
  [hh,Ht,H] = gethhn(Pleaf,m,length(leafs));
  
  ll = length(leafs);
  for i=1:ll
    for j=i+1:ll
      %checketaij(i,j)
      
      iL = leafs(i);
      jL = leafs(j);
      hub = dtreehijub(A,C,T,m,n,iL,jL);
      hij = hh(i,j);
      %thij = tha^(j-i);
      global thij
      %thij = (1-(1-tha)^w)^floor((j-i)/w);
      r = hij/hub;
      %r = hij/thij;
      %r = hij/h0ij;
      if r>2,keyboard,end
      maxr = max(r,maxr);
    end
  end


  fprintf('%s: PARAMS: n=%d m=%d \n',pstr,n,m);
  fprintf('RATIO_________________________ maxr = %1.9f \n',maxr);

end


return

function checketaij(i,j)
global A C T P Pleaf leafs ll m n
leafs = leafs(:)';
for yi=1:m^(i-1)
  y = i2c(yi,m,i-1);
  for w0 = 1:m
    for w1 = w0+1:m
      checketayw(i,j,y,w0,w1);
    end
  end
end
return


function checketayw(i,j,y,w0,w1)
global A C T P Pleaf leafs ll m n
P0 = getcondP(Pleaf,m,ll,[j:ll],[1:i],[y w0]);
P1 = getcondP(Pleaf,m,ll,[j:ll],[1:i],[y w1]);
h0 = .5*sum(abs(P0-P1));

Int = setdiff(1:n,leafs);
nI = length(Int);

xsum = 0;
for xi=1:m^(ll-j+1)
  x = i2c(xi,m,ll-j+1);
  zsum = 0;
  for zi=1:m^(j-i-1)
    z = i2c(zi,m,j-i-1);
    for shi = 1:m^nI
      sh = i2c(shi,m,nI);
      psh = mu(sh,Int);
      pco = muc([y z x],leafs([1:i-1,i+1:ll]),sh,Int);
      f0 = muc(w0,leafs(i),sh,Int)/rho([y w0],1:i);
      f1 = muc(w1,leafs(i),sh,Int)/rho([y w1],1:i);
      zsum = zsum + psh * pco * (f0 - f1);
    end
  end
  xsum = xsum + abs(zsum);
end
h1 = xsum/2;
good = (abs(h0-h1)<1e-9);
if ~good
    'bad1'
  keyboard
end

xsum = 0;
for xi=1:m^(ll-j+1)
  x = i2c(xi,m,ll-j+1);
  for shi = 1:m^nI
    sh = i2c(shi,m,nI);
    psh = mu(sh,Int);
    zsum = 0;
    for zi=1:m^(j-i-1)
      z = i2c(zi,m,j-i-1);
      pco = muc([y z x],leafs([1:i-1,i+1:ll]),sh,Int);
      f0 = muc(w0,leafs(i),sh,Int)/rho([y w0],1:i);
      f1 = muc(w1,leafs(i),sh,Int)/rho([y w1],1:i);
      zsum = zsum + pco * (f0 - f1);
    end
  xsum = xsum + psh*abs(zsum);
  end
end
h2 = xsum/2;
good = (h1+1e-9<h2)
global thij
good = good & (h2+1e-9< thij)
if ~good
    'bad2'
  keyboard
end


%%%%%%%%%%%%%%%%%%%%%%%

xsum = 0;
  for shi = 1:m^nI
    sh = i2c(shi,m,nI);
    psh = mu(sh,Int);
    zsum = 0;
    for zi=1:m^(j-i-1)
      z = i2c(zi,m,j-i-1);
      pco = muc([y z],leafs([1:i-1,i+1:j-1]),sh,Int);
      f0 = muc(w0,leafs(i),sh,Int)/rho([y w0],1:i);
      f1 = muc(w1,leafs(i),sh,Int)/rho([y w1],1:i);
      zsum = zsum + pco * (f0 - f1);
    end
  xsum = xsum + psh*abs(zsum);
  end
h3 = xsum/2;
good = (abs(h2-h3)<1e-9);
if ~good
    'bad3'
  keyboard
end

if 0 %%%%%%%%%%%%%% FALSE %%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%
  xsum = 0;
  for shi = 1:m^nI
    sh = i2c(shi,m,nI);
    psh = mu(sh,Int);
    zsum = 0;
    for zi=1:m^(j-i-1)
      z = i2c(zi,m,j-i-1);
      pco = muc([z],leafs([i+1:j-1]),sh,Int);
      f0 = muc(w0,leafs(i),sh,Int)/rho([w0],i);
      f1 = muc(w1,leafs(i),sh,Int)/rho([w1],i);
      zsum = zsum + pco * (f0 - f1);
    end
    xsum = xsum + psh*abs(zsum);
  end
  %%%%%%%%%%%%% BAD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  h4 = xsum/2;
  good = (abs(h2-h4)<1e-9);
  if ~good
    'bad4'
    keyboard
  end
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

return

function p = muc(x,Ix,s,Is)
global A C T P Pleaf leafs ll m n
Pjj = getcondP(P,m,n,Ix,Is,s);
xi = c2i(x,m,length(Ix));
p = Pjj(xi);
return

function p = mu(sh,I)
% hidden tree measure
% indexed by I\subset[n]
global A C T P Pleaf leafs ll m n
PI = getmargX(P,m,n,I);
si = c2i(sh,m,length(I));
p = PI(si);
return

function p = rho(x,I)
% observed tree measure
% indexed by I\subset leafs
global A C T P Pleaf leafs ll m n
PI = getmargX(Pleaf,m,ll,I);
xi = c2i(x,m,length(I));
p = PI(xi);
return



function i=c2i(c,m,n)
dims = repmat(m,[1 n]);
i = coord2ind(c,dims);
return

function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return

function [w,levels,ll] = treewidth(T)
[n,n] = size(T);
global LEVELS
LEVELS = repmat({[]},n,1);
getLevels(T,1,0);
w = 0;
ll = [];
for ell = 1:n
  w = max(w,length(LEVELS{ell}));
  if ~isempty(LEVELS{ell})
    ll(end+1)=length(LEVELS{ell});
  end
end
levels = LEVELS;
return


function getLevels(T,u,ell)
global LEVELS
ell = ell+1;
LEVELS{ell}(end+1) = u;
ukids = find(T(u,:));
for vi=1:length(ukids)
  v = ukids(vi);
  getLevels(T,v,ell);
end

return


function [dp,j0] = dpath(T,i,j)
% bounds hh(i,j) for a tree

% find j0, where
% j0 is the minimum descendant of i s.t. j0>=j


Ti = getSubt(T,i);

jj = Ti(Ti>=j);
if isempty(jj)
  dp = [];
  j0 = 0;
else
  j0 = min(jj);
  dp = parentseq(T,j0,i);
end
return

function pseq = parentseq(T,u,r)
v = u;
el = 0;
pseq = [u];
while v~=r
  el = el+1;
  rents = find(T(:,v));
  v = rents;
  pseq = [v pseq];
end
return



function sT = getSubt(T,v)
% a list of v's descendents in T
% includes v itself
global SUBT
SUBT = [];
Pi(T,v);
sT = SUBT;
return

function Pi(T,v)
global SUBT
SUBT(end+1) = v;
vkids = find(T(v,:));
for vi=1:length(vkids)
  v1 = vkids(vi);
  Pi(T,v1);
end