function mccheck
global n m A C P
while 1
  A = pnormdim(rand(m,m,n-1),1)
  C = pnormdim(rand(m,1,1  ),1);
  P = markovfillprob(A,C,n);
  
  [hh,gHt,gH] = gethhn(P,m,n);
  mHt = [];
  for t=1:n
    %mHt = max(mH,getH(A,t));
    mHt(t) = getH(A,t);
  end
  %good = (gH <= mH)
  good = all(gHt <= mHt);
  if ~good
    'big nono'
    keyboard
  end

  
  if 1  
    theta = getTH(A);
    for i=1:n
      for j=i+1:n
        markovaux(i,j);
        if 0
          etaij = hh(i,j);
          %thaij = prod(theta(i:j-1));
          thaij = prod(theta(i:end));
          
          %good = (abs(etaij-thaij)<1e-9)
          good = etaij <= thaij
          if ~good
            'nono'
            keyboard
          end
        end
      end
    end
  end
  
end
return

function markovaux(i,j)
global n m A C P
%load DEBUG
dims = repmat([m],[1 i]);
for yi = 1:m^i
  y = ind2coord(yi,dims);
  h0 = etaijX(P,m,n,i,j,y);
  h1 = markovauxy1(i,j,y);
  %compsumz(i,j,y)
  good = ((h0-h1) < 1e-9);
  if ~good
    'gotcha'
    keyboard
  end
  hhat = markovauxyhat(i,j,y);
  good = h1 <= hhat
  if ~good
    'gotcha2'
    keyboard
  end  
end
return

% works
function S = sumz0(i,j,x,y)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
sz1 = 0;
for zi = 1:m^(j-i-1)
  z = ind2coord(zi,dimz);
  Pzx = getPX(P,m,n,[y z x],[1:n])/getPX(P,m,n,[y],[1:i]);
  sz1 = sz1 + Pzx;
end
sz2 = 0;
for w=1:m
  for zi = 1:m^(j-i-1)
    z = ind2coord(zi,dimz);
    Pwzx = getPX(P,m,n,[y(1:i-1) w z x],[1:n])/getPX(P,m,n,[y(1:i-1)],[1:i-1]);
    sz2 = sz2 + Pwzx;
  end
end
S = sz1 - sz2;
return

function S = sumz01(i,j,x,y)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
sz1 = 0;
for zi = 1:m^(j-i-1)
  z = ind2coord(zi,dimz);
  %Pzx0 = getPX(P,m,n,[y z x],[1:n])/getPX(P,m,n,[y],[1:i]);
  u = [z x];
  Pzx = A(u(1),y(i),i);
  for k=i+2:n
    Pzx = Pzx * A(u(k-i),u(k-i-1),k-1);
  end
  sz1 = sz1 + Pzx;
end
sz2 = 0;
for w=1:m
  for zi = 1:m^(j-i-1)
    z = ind2coord(zi,dimz);
    %Pwzx0 = getPX(P,m,n,[y(1:i-1) w z x],[1:n])/getPX(P,m,n,[y(1:i-1)],[1:i-1]);
    u = [z x];
    if i > 1
      Pwzx = A(w,y(i-1),i-1);
    else
      Pwzx = C(w);
    end
    Pwzx = Pwzx * A(u(1),w,i);
    for k=i+2:n
      Pwzx = Pwzx * A(u(k-i),u(k-i-1),k-1);
    end
    sz2 = sz2 + Pwzx;
  end
end
S = sz1 - sz2;
return


% works
function S = sumz2(i,j,x,y)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
S = 0;
for zi = 1:m^(j-i-1)
  z = ind2coord(zi,dimz);
  u = [z x];
  %comzx = 1;
  %for k=i+2:n
  %  comzx = comzx * A(u(k-i),u(k-i-1),k-1);
  %end
  comzx = marku(u,i);
  
  left = A(u(1),y(i),i); 

  right = 0;
  for w=1:m
    if i > 1
      pw = A(w,y(i-1),i-1);
    else
      pw = C(w);
    end
    pw = pw * A(u(1),w,i);
    right = right + pw;
  end
  S = S + comzx * (left-right);
end
return


function S = sumz3(i,j,x,y)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
S = 0;
for zi = 1:m^(j-i-1)
  z = ind2coord(zi,dimz);
  comzx = marku(z,i);
  comzx = comzx * A(x(1),z(end),j-1);
  left = A(z(1),y(i),i); 
  right = 0;
  for w=1:m
    if i > 1
      pw = A(w,y(i-1),i-1);
    else
      pw = C(w);
    end
    pw = pw * A(z(1),w,i);
    right = right + pw;
  end
  S = S + comzx * (left-right);
end
return


function S = sumzhat(i,j,x,y,w)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
%w = getwhat(y,i);
S = 0;
for zi = 1:m^(j-i-1)
  z = ind2coord(zi,dimz);
  comzx = marku(z,i);
  comzx = comzx * A(x(1),z(end),j-1);
  left = A(z(1),y(i),i); 
  right= A(z(1),w   ,i);
  S = S + comzx * (left-right);
end
return


function w = getwhat(y,i)
global n m A C P
Pi = A(:,:,i);
y = y(i);
Py = repmat(Pi(:,y),[1 m]);
dd = sum(abs(Pi-Py));
[x,w] = max(dd);
return



function p = marku(u,i)
global n m A C P
p = 1;
for k=2:length(u)
  p = p * A(u(k),u(k-1),k+i-1);
end
return



% works
function S = markovauxy0(i,j,y)
global n m A C P
S = 0;
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
for xi = 1:m^(n-j+1)
  x = ind2coord(xi,dimx);
  if 0
    sz1 = 0;
    for zi = 1:m^(j-i-1)
      z = ind2coord(zi,dimz);
      %zxi = coord2ind([z x],dims);
      %PZX = getcondP(P,m,n,[i+1:n],1:i,y);
      %Pzx = PZX(zxi);
      Pzx = getPX(P,m,n,[y z x],[1:n])/getPX(P,m,n,[y],[1:i]);
      sz1 = sz1 + Pzx;
    end
    
    sz2 = 0;
    for w=1:m
      for zi = 1:m^(j-i-1)
        z = ind2coord(zi,dimz);
        %wzxi = coord2ind([w z x],dims1);
        %PWZX = getcondP(P,m,n,[i:n],1:i-1,y(1:i-1));
        %Pwzx = PWZX(wzxi);
        Pwzx = getPX(P,m,n,[y(1:i-1) w z x],[1:n])/getPX(P,m,n,[y(1:i-1)],[1:i-1]);
        %sw = sw + Pwzx;
        sz2 = sz2 + Pwzx;
      end
    end
    %S = S + Pzx - sw;
    S = S + abs(sz1 - sz2);
  end
  S = S + abs(sumz2(i,j,x,y));
end
S = S/2;
return


function compsumz(i,j,y)
global n m A C P
S = 0;
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
for xi = 1:m^(n-j+1)
  x = ind2coord(xi,dimx);
  sz0 = sumz0(i,j,x,y);
  sz01 = sumz01(i,j,x,y);
  if abs(sz0-sz01)>1e-9
      'dang'
  end
end
S = S/2;
return


function S = markovauxy1(i,j,y)
global n m A C P
S = 0;
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
for xi = 1:m^(n-j+1)
  x = ind2coord(xi,dimx);
  if j-i > 1
    px = marku(x,j-1);
    S = S + px*abs(sumz3(i,j,x,y));
  else
    S = S + abs(sumz2(i,j,x,y));
  end
end
S = S/2;
return

function S = markovauxyhat(i,j,y)
global n m A C P
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
dims = repmat([m],[1 n-i]);
dims1 = repmat([m],[1 n-i+1]);
Sw = [];
for w=1:m
  S = 0;
  for xi = 1:m^(n-j+1)
    x = ind2coord(xi,dimx);
    if j-i > 1
      px = marku(x,j-1);
      S = S + px*abs(sumzhat(i,j,x,y,w));
    else
      S = S + abs(sumz2(i,j,x,y));
    end
  end
  S = S/2;
  Sw(w) = S;
end
S = max(Sw);
return



function S = markovauxy(i,j,y)
global n m A C P
S = 0;
dimx = repmat([m],[1 n-j+1]);
dimz = repmat([m],[1 j-i-1]);
for xi = 1:m^(n-j+1)
  x = ind2coord(xi,dimx);
  pxx = 1;
  for l=j+1:n
    pxx = pxx * A(x(l-j+1),x(l-j),l-1);
  end
  sz = 0;
  if j > i+1
    for zi = 1:m^(j-i-1)
      z = ind2coord(zi,dimz)
      pxz = A(x(1),z(end),j-1);
      pzz = 1;
      for k=i+2:j-1
        pzz = pzz * A(z(k-i),z(k-i-1),k-1);
      end
      pzy = A(z(1),y(i),i);
      sw = 0;
      for w=1:m
        if i > 1
          pwy = A(w,y(i-1),i-1);
        else
          pwy = C(w);
        end
        sw = pwy * A(z(1),w,i);
      end
      sz = sz + pxx * pxz * pzz * (pzy - sw);
    end
  else
    pxy = A(x(1),y(i),i);
    sw = 0;
    for w=1:m
      if i > 1
        pwy = A(w,y(i-1),i-1);
      else
        pwy = C(w);
      end
      sw = pwy * A(x(1),w,i);
    end
    sz = sz + pxx * (pxy - sw);
  end
  S = S + abs(sz);
end
S = S/2;
return

function etaijconv
global n mhid mobs Po Ph Poh

while 1
  Ph   = pnormdim(rand(mhid^n,     1),1);
  Poh0 = pnormdim(rand(mobs,mhid),1);
  %Poh1 = pnormdim(rand(mobs,mhid),1);
  %Poh2 = pnormdim(rand(mobs,mhid),1);
  
  for i=1:n
    dims = repmat([mobs],[1 i]);
    for j=i+1:n
      rho1 = pnormdim(rand(mobs,1),1);
      rho2 = pnormdim(rand(mobs,1),1);
      Poh1 = Poh0; Poh1(:,h) = rho1;
      Poh2 = Poh0; Poh2(:,h) = rho2;
      a = rand;
      b = 1-a;
      Pohc = a*Poh1 + b*Poh2;
      Po1 = latentiifillprob(Ph,Poh1,mhid,mobs,n);
      Po2 = latentiifillprob(Ph,Poh2,mhid,mobs,n);
      Poc = latentiifillprob(Ph,Pohc,mhid,mobs,n);
      
      for xi = 1:mobs^i
        X = ind2coord(xi,dims);
        h1 = etaijX(Po1,mobs,n,i,j,X);
        h2 = etaijX(Po2,mobs,n,i,j,X);    
        hc = etaijX(Poc,mobs,n,i,j,X);
        g = (hc < a*h1 + b*h2);
        if ~g
          'not convex'
          keyboard
        end  
      end
    end
  end
end
return