function treesketchpad
global m n
global CONFIGPROBS
global doplots domask
doplots = 1;
domask = 0;
gatherstats(6)
return
m = 2;
n = 3;
%T = randtree(n)
% make a chain tree
T = zeros(n); T(2:end,1:n-1) = eye(n-1)
psi = randpsi(m)
probconfigs(T,psi)
CONFIGPROBS
[b,A,C] = is_markov(CONFIGPROBS,m,n)
return
n = 6;
m = 2;
T = randtree(n)
psi = randpsi(m)
probconfigs(T,psi)
CONFIGPROBS
return

function gatherstats(n0)
global m n
global CONFIGPROBS
global VALS BEES ma
global Ainf Binf n
global MA
global doplots domask
alpha = .2;
%VALS = []; BEES = [];
n = n0;
%n = 4;
m = 2;
if ~domask
  [Ainf,Binf] = lipab(m,n);
end
%t = round(randinseg(1,n));
%t = 1;
%dimx = repmat([m],[1 t]);
bb = .001:.01:2*n;
F = zeros(m^n,1);
T = randtree(n)
% make a chain tree
T = zeros(n); T(2:end,1:n-1) = eye(n-1)
%xtind = round(randinseg(1,m^t));
rr = [];
yy = [];
ma = 0;
for i=1:10000000
  %for b=bb
  %b = randinseg(1,n*2);
  %psi = ones(m);
  psi = rand(m);
  psi = psi+.1;
  %psi = psi * .5;
  %psi(1,2) = b;
  %psi(2,1) = b;
  %psi = randpsi(m);
  %a = psi(1,1);
  %b = psi(1,2);
  %c = psi(2,2);
  %vals = [];
  T = randtree(n);
  t = round(randinseg(1,n));
  dimx = repmat([m],[1 t]);
  xtind = round(randinseg(1,m^n));
  Xt = ind2coord(xtind,dimx)
  
  probconfigs(T,psi);
  if 0
    [ism,A,C] = is_markov(CONFIGPROBS,m,n);
    if ~ism
      'not markov'
      keyboard
    end
    Pmark = markovfillprob(A,C,n);
    H = getH(A,t);
    [Ef1,K1] = condexp(Pmark,F,Xt         ,m,n);
    [Ef0,K0] = condexp(Pmark,F,Xt(1:end-1),m,n);
    Kmark = K1 - K0;
  end
  
  [Ef1,K1] = condexp(CONFIGPROBS,F,Xt         ,m,n);
  [Ef0,K0] = condexp(CONFIGPROBS,F,Xt(1:end-1),m,n);
  K = K1 - K0;

  if ~domask
    [phi,fv] = maxphi(K);
  else
    mask = optmask(K,m,n);
    phi = fillgmask(mask,m);
    fv = abs(phi*K);
  end
  %if fv0 + 1e-5 < fv1
  %  keyboard
  %end
  %fv = fv0;
  
  if 0
    VALS(end+1) = fv; BEES(end+1) = b;
    %a = -log(1-fv/n)/b;
    %ma = max(a,ma);
    %aa = -log(1-VALS/n)./BEES;
    aa = VALS/n./BEES;
    aa = max(aa);
    if aa > ma
      %ma = max(aa);
      load TREESTATS
      ma = max(aa,MA(n));
      MA(n) = ma;
      save TREESTATS MA
    end
    %fprintf('                          ma = %1.9f exp(-a)=%1.5f \n',ma,exp(-ma));
    fprintf('n=%d                          ma = %1.9f \n',n,ma);
    %if fv > H
    %    'bound dont hold'
    %    keyboard
    %end
  end
  
  %r = max(b,1/b)/2;
  Rmin = min(psi(:));
  Rmax = max(psi(:));
  r = Rmax/Rmin;
  %r = max([a/b,b/a,b/c,c/b]);
  %r = max([a/b,b/a]);
  %r = max([b/c,c/b]);
  %B = n*(1-exp(-alpha*r));
  if fv > r
    'bound fails'
    keyboard
  end
  
  % seems not too crude
  [fv1,xhat,K1] = boundYxhat(CONFIGPROBS,phi,Xt,m,n);
  if fv > r + 1e-9
    'too crude'
    keyboard
end
  % this is def. too crude
  %fv2 = n * sum(abs(K1));
  %fv2 = Psi(K1,m,n);
  %fv2 = .5*Psi(K,m,n);
  %if fv2 > r + 1e-9
   % 'too crude'
   % keyboard
  %end
  
  h = neostrict1(CONFIGPROBS,m,n);
  if h > 1-exp(1-r)
    'your gen thing dont work'
    keyboard
end

    VALS(end+1) = h; BEES(end+1) = r;


  if doplots
    hold off
    plot(BEES,VALS,'.b');
    %r = min(BEES,1./BEES);
    %r = min(b,1/b); rr(end+1) = r;
    %y = sum(r.^(0:n)); yy(end+1) = y;
    hold on
    %plot(rr,yy,'.r');
    x = 1:.01:max(BEES);
    f=1-exp(1-x);
    %f=n*(-exp(-ma*x)+exp(-ma))+.5;
    %f=n*(1-exp(-ma*x)); % USE THIS ONE!
    %f=n*(1-exp(-alpha*x)); % less precise, but universal
    %f=n*(ma*x); % even less precise
    plot(x,f,'.r')
    drawnow
  end
    %vals(end+1) = fv;
  %end
end
return


function myobj(psi)


return

function probconfigs(T,psi)
global n m % state space
global CONFIGPROBS Ztree
CONFIGPROBS = zeros(m^n,1);
dims = repmat([m],[1 n]);

[nodesi,nodesj]=find(T);

for xind = 1:m^n
  x = ind2coord(xind,dims);
  p = 1;
  for e=1:length(nodesi)
    i = nodesi(e);
    j = nodesj(e);
    p = p * psi(x(i),x(j));
  end
  CONFIGPROBS(xind) = p;
end
Ztree = sum(CONFIGPROBS);
CONFIGPROBS = CONFIGPROBS/Ztree;
return

function psi = randpsi(m)
psi = rand(m);
%psi = psi - diag(diag(psi));
psi = (psi + psi')/2;
return

function T = randtree(n)
T = zeros(n);
%T(1,2) = 1;
for i=2:n
  % attach each node to a random one in the tree
  jj = randperm(i-1); j = jj(1);
  T(i,j) = 1;
end
%used = 1:2;
%free = 3:n;
return

% nevermind this garbage
function tryex(T,psi,Xt)
global n m % state space
global CONFIGPROBS Ztree
probconfigs(T,psi);
F = zeros(m^n,1);
Rmin = min(psi(:));
Rmax = max(psi(:));
r = Rmax/Rmin;

Ppsi = CONFIGPROBS;
[Ef1,K1] = condexp(Ppsi,F,Xt         ,m,n);
[Ef0,K0] = condexp(Ppsi,F,Xt(1:end-1),m,n);
K = K1 - K0;
[phi,fv0] = maxphi(K);

K00 = K;

t = length(Xt);
dims = repmat([m],[1 n  ]);
dim0 = repmat([m],[1 n-t+1]);

tt0 = 1:t;
tt1 = 1:t-1;
xinds = getinds(m,n,Xt(tt1),tt1);

px0 = getPX(Ppsi,m,n,Xt(tt0),tt0);
px1 = getPX(Ppsi,m,n,Xt(tt1),tt1);

ysum = 0;
K = 0*F;
%for xi = xinds(:)'
for xi = 1:m^n
  %x0 = ind2coord(xi,dim0);
  x = ind2coord(xi,dims);
  %xt0 = [Xt(1:t)   x0(2:end)]
  xt0 = x; xt0(tt0) = Xt(tt0);
  %xt1 = [Xt(1:t-1) x0(1:end)]
  xt1 = x; xt0(tt1) = Xt(tt1);
  xi0 = coord2ind(xt0,dims);
  xi1 = coord2ind(xt1,dims);
  %f0 = F(xi0); f0 = F(xi0);
  p0 = Ppsi(xi0)/(px0+realmin);
  p1 = Ppsi(xi1)/(px1+realmin);
  K(xi0) = K(xi0) + p0;
  %K(xi1) = K(xi1) - p1;  
end

if totdif(K00,K) > 1e-6
  'Ks dont agree'
  keyboard
end


if 0
  % Ztree doesn't matter
  Ppsi = CONFIGPROBS * Ztree;
  [Ef1,K1] = condexp(Ppsi,F,Xt         ,m,n);
  [Ef0,K0] = condexp(Ppsi,F,Xt(1:end-1),m,n);
  K = K1 - K0;
  [phi,fv1] = maxphi(K);
  
  if abs(fv0-fv1) > 1e-6
    'Ztree matters'
    keyboard
  end
end
return


function inds = getinds(m,n,X,ii)
dims = repmat([m],[1 n]);
v = repmat({1:m},[1 n]);
for j=1:length(ii)
  v{ii(j)} = X(j);
end
xx = list_all_ind_g(v);
inds = coord2ind(xx,dims);
return