function mcd_nstate
% writeup sanity check
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
global nstate0
nstate0 = 3;
nstate = nstate0;
%lenseq = 2;
%global A
%global ec
global FN
global NN
NN=5;
maxR = 0;
global ksi

global AA
global C
sc = [0 1];
for T=1:10000
  FN = randFN(nstate0,NN);
  for lenseq=2:NN
    if 0
      [AA,C] = initAAC(sc);
      Y0 = Y_gam;
      Y1 = Y_gamc;
      g = (abs(Y0)<=abs(Y1)+1e-10);
    end
    
    nstate = nstate0;
    [AA,C] = initAAC(sc);
    h0 = SC(AA);
    Y0 = Y_gam;
    nstate = 2;
    
    YY1 = [];
    h0 = SC(AA);
    [x,y,z] = size(AA);
    AA = [];    
    C(1) = 0;
    C(2) = 1;
    aa = h0:.01:1;
    for aj = 1:length(aa)
      for t=1:z
        a = aa(aj);
        b = a-h0;
        AA(:,:,t) = [[a,b];[1-a,1-b]];
      end
      YY1(end+1) = Y_gam;
    end
    h1 = SC(AA);

    %[AA,C] = smartcollapse2(AA,C);
    %h1 = SC(AA);
    %Y1 = Y_gam;
    %gy = (abs(Y0)<=abs(Y1)+1e-10);
    gy = (abs(Y0)<=max(abs(YY1))+1e-10);
    gh = (h1 <= h0 +1e-10);
    if ~(gy & gh)
      'problem!!!!'
      keyboard
    end
    fprintf('T=%d n=%d \n',T,lenseq);

    
    % BAD APPROACH (FALSE)
    %nstate = nstate0;
    %[AA,C] = initAAC(sc);
    %sc0 = SC(AA);
    %Y0 = Y_gam;
    %nstate = 2;
    %[AA,C] = smartcollapse(AA,C);
    %sc1 = SC(AA);
    %Y1 = Y_gam;
    %g = (abs(Y1)<=abs(Y0)+1e-10);
    %fprintf('T=%d g=%d sc(A)=%1.3f  SC(A1)=%1.3f  \n',T,g,sc0,sc1);
    %if ~g
    %  'problem!!!!'
    %  keyboard
    %end
    fprintf('good T=%d \n',T);
  end
end
return


function A = initA(sc)
global nstate
while 1
  A = rand(nstate);
  A = pnormdim(A,1);
  if ( (sc(1)<=SC(A)) & (SC(A)<=sc(2)) )
    break
  end
end
return

function [AA,C] = initAAC(sc)
global NN
global nstate
C = pnormdim(rand(nstate,1),1);
%AA = initA(sc);
%return
AA = zeros(nstate,nstate,NN-1);
A0 = initA(sc);
for i=1:NN-1
  AA(:,:,i) = initA(sc);
  %AA(:,:,i) = A0;
end
%C(1) = 0;
%C(2) = 1;
return



function FN = randFN(nstate,NN,d)
if ~exist('d','var')
  d=1;
end
FN = {};
for j=1:NN
  FN{j} = randF(nstate,j,d);
end
return

function F = randF(nstate,lenseq,d)
% F is a random function $\X^n \to [0,1]$ that sat. ibd w/c=1
X = list_all_ind(repmat(nstate,1,lenseq));
K = 10;
%K = 1;
%K = lenseq;
%F(1) = 0;
PROGstart(size(X,1),'generating randF...');
for i=1:size(X,1)
%for i=2:size(X,1)
  good = 0;
  while ~good
    F(i) = K*rand;
    %F(i) = round(K*rand);
    good = 1;
    for j = 1:i-1
      good = goodij(X,F,i,j,d);
      if ~good, break, end;
    end
  end
  %fprintf('randF: done %d of %d \n',i,size(X,1));
  PROGupdate(i);
end
PROGend;
return

function b = goodij(X,F,i,j,d)
x = X(i,:);
y = X(j,:);
if (length(find(x~=y))==1)
  b = (abs(F(i)-F(j))<=d);
  %b = (abs(F(i)-F(j))==1);
else
  b = 1;
end
return



function Y=Y_gam(MC,FN)
global lenseq
n = lenseq;
gam = gamma(n,[],MC,FN);
%Y = gam(1)-sum(gam(2:end));
Y=sum(gam);
return


function gam = gamma(n,x,MC,FN)
global nstate
gam = zeros(nstate,1);
if (n==2)
  for k=1:nstate
    for x1=1:nstate
      gam(k) = gam(k) + p_0(x1,MC)*(p_n(k ,1 ,1,MC)*F([1 k x],FN) - p_n(k ,x1,1,MC)*F([x1 k x],FN));
    end
  end
else
  for k=1:nstate
    gamkx = gamma(n-1,[k x],MC,FN);
    for k1=1:nstate
      gam(k) = gam(k) + p_n(k,k1,n-1,MC) * gamkx(k1);
    end
  end
end  
return

function gam = gammac(n,x)
global nstate
global AA
gam = zeros(nstate,1);
if (n==2)
  for k=1:nstate
    for x1=1:nstate
      gam(k) = gam(k) + p_0(x1)*(p_n(k ,1 ,1)*F([1 k x]) - p_n(k ,x1,1)*F([x1 k x]));
    end
  end
else
  [h,k1,k2] = SCA(AA(:,:,n-1));
  for k=1:nstate
    gamkx = gammac(n-1,[k x]);
    tha = pnormdim(rand(nstate-1,1),1);
    p1 = p_n(k,k1,n-1);
    p2 = 0;
    sd = setdiff(1:nstate,k1);
    for j=1:length(sd)
      p2 = p2 + tha(j)*p_n(k,sd(j),n-1);
    end
    gam(k) = p1 * gamkx(k1) + p2 *  gamkx(k1);
  end
end  
return


function p=p_0(k,MC)
%global C
p=MC.C(k);
return

function p=p_n(k,k1,n,MC)
%global AA
p=MC.AA(k,k1,n);
%p=AA(k,k1);
return

function r = F(x,FN)
%global nstate0
%global lenseq
nstate = FN.nstate;
if ~x(end), x=x(1:end-1); end;
lenseq = length(x);
%global FN
%FF = FN{lenseq};
FF = FN.FFF{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
return


function [h,i0,j0] = SCA(A)
[n,n] = size(A);
h = 0;
for i=1:n
  for j=i+1:n
    x = A(:,i);
    y = A(:,j);
    hxy = max(abs(x-y));
    if (hxy>h)
      i0=i;j0=j;h=hxy;
    end
  end
end
return



%obsolete
function Y=Y_gamc
global lenseq
global AA
n = lenseq;
gam = gammac(n,[]);
Y=sum(gam);
return