function mcd_testmain
NN = 9;
global nstate
nstate = 2;
global FN
Fdmax = 0;
fname = '/usr1/lkontor/MCDBAD.mat';
for t=1:1000000
  %d = rand*10;
  d = 1;
  FN = randFN(nstate,NN,d,'exact');
  %Fmax = max(FN{NN});
  %Fmin = min(FN{NN});
  %Fdel = Fmax - Fmin;
  %Fdmax = max(Fdel,Fdmax)
  %n = ceil(rand*(NN-2))+1;
  n = round(randinseg(2,NN-1));
  %n = 3;
  n = 1;
  a = rand; 
  %b = rand;
  b = randinseg(0,a);
  %x = round(rand)+1;
  x = [];
  %L = NN - n;
  %L = NN;
  %if (L>0)
    %x0 = round(rand(1,L));
    %x1 = 1-x0;
  %  x0 = round(rand(1,L))
  %  x1 = round(rand(1,L))
  %else
  %  x0 = [];
  %  x1 = [];
  %end
  
  for a=.001:.05:1
    for b=0.001:.05:1
  %    a = 1;
  %    b = 0;
  
      %[alf0, bet0] = alfbet(n,[0 x0],a,b);
      %[alf1, bet1] = alfbet(n,[1 x1],a,b);
      %[alf0, bet0] = alfbet(n,[x0],a,b);
      %[alf1, bet1] = alfbet(n,[x1],a,b);
      [alf, bet] = alfbet(n,[],a,b);
      
      nn=0:n-1;
      h = abs(a-b);
      H0 = sum(h.^nn);
      %H0 = Hn(3,a,b);
      %g00 = ( abs(alf0-bet0) <= d*H0+1e-10);
      %g11 = ( abs(alf1-bet1) <= d*H0+1e-10);
      g = ( abs(alf-bet) <= d*H0+1e-10);
      %if ~(g00 & g11)
      if ~g
        'its all over'
        good = robustsave(fname,10,FN,n,x,a,b);
        keyboard
      end
      
      if 0
        g00 = ( abs(alf0-bet0) <= d*4+1e-10);
        g01 = ( abs(alf0-bet1) <= d*4+1e-10);
        g10 = ( abs(alf1-bet0) <= d*4+1e-10);
        g11 = ( abs(alf1-bet1) <= d*4+1e-10);
        %if ~(g00 & g01 & g10 & g11)
        %  'meod muzar'
        %  keyboard
        %end
      end

    end
  end
  
  
  %fprintf('good t = %d  a=%1.3f b=%1.3f  FF=[%s] \n',t,a,b,num2str(FN{3}));
  fprintf('good t = %d n = %d \n',t,n);
  %fprintf('good t = %d n = %d maxE = %1.9f \n',t,n,maxE);
  %fprintf('good t = %d n = %d minG = %1.7f maxG = %1.7f \n',t,n,minG,maxG);
  %if (rand<0.005)
  %end
end
return



function h = Hn(n,a,b)
h = 3;
return


function h = Hn_old(n,a,b)
if (n==1)
  h = 1;
else
  H0 = Hn(n-1,a,b);
  h = H0 + (a-b)^(n-1);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function [alf, bet] = alfbet(n,x,a,b)
if (n==1)
  alf = F([0 x]);
  bet = F([1 x]);
else
  [alf0x, bet0x] = alfbet(n-1,[0 x],a,b);
  [alf1x, bet1x] = alfbet(n-1,[1 x],a,b);
  alf =    a  * alf0x -    b  * bet0x;
  bet = (a-1) * alf1x - (b-1) * bet1x;
end
return



function [alf, bet] = alfbet_old(n,x,a,b)
%if (n==2)
%  alf =  a   *F([1 1 x]) - b   *F([2 1 x]);
%  bet = (a-1)*F([1 2 x]) -(b-1)*F([2 2 x]);
if (n==1)
  alf = F([1 x]);
  bet = F([2 x]);
else
  [alf0x, bet0x] = alfbet(n-1,[1 x],a,b);
  [alf1x, bet1x] = alfbet(n-1,[2 x],a,b);
  alf =    a  * alf0x -    b  * bet0x;
  bet = (a-1) * alf1x - (b-1) * bet1x;
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


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


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function ind = coord2ind(coords,dims)
% dims is a ROW vector of dimensions
% coords is a matrix of ROW coordinates vectors
% inds is a list of indices
% this is a recent (16/02/2002) change

% compare with matlab's sub2ind routine

if (length(dims)<size(coords,2))
  dims = repmat(dims,[1,size(coords,2)]);
end

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

return


function ind = coord2ind_old(coords,dims)
% NOTE: coords and dims are both ROW vectors

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

% REALLY OLD VERSION, to illustrate the idea
%ind = 0;
%coords = coords - 1;
%f = 1;
%for ij = 1:length(dims)
%  ind = ind + coords(ij)*f;
%  f = f*dims(ij);
%end
%ind = ind+1;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------





%--------------------------------------------------------------
function FN = randFN(nstate,NN,d,estr)
%if ~exist('d','var')
%  d=1;
%end
if ~exist('estr','var')
  exact = 0;
else
  exact = strcmp(estr,'exact');
end
FN = {};
for j=1:NN
  FN{j} = randF(nstate,j,d,exact);
end
return

function F = randF(nstate,lenseq,d,exact)
% 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;
    if exact
      F(i) = round(F(i));
    end
    good = 1;
    for j = 1:i-1
      good = goodij(X,F,i,j,d,exact);
      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,exact)
x = X(i,:);
y = X(j,:);
if (length(find(x~=y))==1)
  if exact
    b = (abs(F(i)-F(j))==1);
  else
    b = (abs(F(i)-F(j))<=d);
  end
else
  b = 1;
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function PROGend
global PROGsuppress
%if isempty(PROGsuppress), PROGsuppress = 0; end;
if PROGsuppress
  return
end

fprintf('\n');

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function PROGstart(tot,str)
%global nohupmode
global PROGsuppress
if isempty(PROGsuppress), PROGsuppress = 0; end;
if PROGsuppress
  return
end


if ~exist('str'), str = ''; end
%if nohupmode, str = ['(nohup)',str]; end

if ~isempty(str)
  fprintf('%s \n',str);
end

global PROGtot
PROGtot = tot;
global PROGdots
PROGdots = 0;
global PROGstr
PROGstr = str;

%tic;  % for nohupmode
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function PROGupdate(iter)
global PROGtot
global PROGdots
global PROGstr
global nohupmode
%if isempty(nohupmode)
%  nohupmode = 0;
%end
global PROGsuppress
%if isempty(PROGsuppress), PROGsuppress = 0; end;
if PROGsuppress
  return
end


curdots = round(iter/PROGtot*100);
if curdots > PROGdots
  % need more dots!
  moredots = curdots - PROGdots;
  str = repmat('.',1,moredots);
  fprintf('%s',str);
  PROGdots = PROGdots + moredots;
  %if nohupmode
  %  PSTR = PROGstr;
  %  PTIME = toc;
  %  PDOTS = PROGdots;
  %  save PROGdata PDOTS PSTR PTIME
  %  tic;
  %end %if ~nohupmode
end

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
% list all n-tuples of the dimension vector v

function list = list_all_ind(v)
list = [];
if ~isempty(v)
  d1 = v(1);

  L1 = list_all_ind(v(2:end));
  [m,n] = size(L1);
  m = max(m,1);
  D1 = repmat(1:d1,[m,1]); D1 = D1(:);
  L2 = repmat(L1,[d1,1]);
  list = [D1,L2];

end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function r = randinseg(a,b)
r = rand*(b-a) + a;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function good = robustsave(fname,ntries,varargin)
% sample usage:
%>> a='a'                                         
%>> b=1
%>> good = robustsave('scrap.mat',1,a,b)    


%if mod(nargin,2)
%  error('robustsave requires an even number of arguments after fname,ntries');
%end

vnams = {};
for j=3:nargin
  vnam = inputname(j);
  if isempty(vnam)
    error('bad input name');
  end
  vnams{end+1} = vnam;
end
vvals = varargin;
varargin = [];

for itry=1:ntries
  good = 1;
  rehash
  try
    save(fname,'vnams','vvals');
  catch
    good = 0;
    fprintf('robustsave fname = %s, try = %d out of %d, error = %s \n',fname,itry,ntries,lasterr);
  end
  if good
    break
  else
    %pause(rand*10);
  end
end

if ~good
  fprintf('robustsave fname = %s, failed after %d tries, error = %s \n',fname,ntries,lasterr);
  %error('could not execute robustsave');
end
return



%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------