% this FULLY VINDICATES the alf-bet inductive step
function [A,gam,fval] = alfbetcheck(c1,c2,h)

%n=10;
m=2;

global nstate lenseq
nstate = m;
%lenseq = n;

T = 1;

adim = m^2*T;
gdim = m^2;
xdim = adim + gdim;


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
LB(adim+1:end) = LB(adim+1:end) - Inf;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + Inf;  % no box constraints on phi


% normalization constraints
Aeq = zeros(m*T,xdim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);


% stricture constraints
SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(size(SIGNS,1)*m*(m-1)/2*T,m^2*T);
rind = 0;

for t=1:T
  for i1=1:m
    for i2=i1+1:m
      for s = 1:size(SIGNS,1)
        rind = rind+1;
        ss = SIGNS(s,:);
        for j=1:m
          ind1 = coord2ind([j i1 t],[m m T]);
          ind2 = coord2ind([j i2 t],[m m T]);
          Aina(rind,ind1) =  ss(j);
          Aina(rind,ind2) = -ss(j);
        end
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h*2;

Aing = zeros(0,gdim);
rind=0;
% gamma sum constraint
for s=[-1 1]
  for i=1:2
    rind = rind+1;
    cind1 = coord2ind([i 1],[m m]);
    cind2 = coord2ind([i 2],[m m]);
    Aing(rind,cind1) = s;
    Aing(rind,cind2) = -s;
    Bing(rind,1) = c1;
  end
  
  rind = rind+1;
  cind1 = coord2ind([1 1],[m m]);
  cind2 = coord2ind([2 2],[m m]);
  Aing(rind,cind1) = s;
  Aing(rind,cind2) = -s;
  Bing(rind,1) = c1+c2;

  rind = rind+1;
  cind1 = coord2ind([1 2],[m m]);
  cind2 = coord2ind([2 1],[m m]);
  Aing(rind,cind1) = s;
  Aing(rind,cind2) = -s;
  Bing(rind,1) = c1+c2;
 
end

Ain = assembleblock({Aina,Aing});
Bin = [Bina;Bing];

A0 = pnormdim(ones(m,m,T),1);
A0 = A0+randn(size(A0))*1e-2;
A0 = pnormdim(A0,1);

G0 = randn(m,m);

x0 = [A0(:); G0(:)];

[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));

y = myobjm(x);

global A FN

A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);
return


function y = myobjm(x)
global nstate lenseq

m = nstate;
%n = lenseq;
T = 1;

adim = m^2*T;
gdim = m^2;
xdim = adim + gdim;

A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);

alf0 = gam(1,1);
alf1 = gam(2,1);
bet0 = gam(1,2);
bet1 = gam(2,2);

a = A(1,1);
b = A(1,2);

sA = alf0-alf1;
sB = bet1-bet0;
sC = alf1-bet1;

s = a*sA + b*sB + sC;

y = -s^2;
return