function [success,rsuccess]=checksol(csmall,gamma,m,ell,K,xi,A,SG)
% [success,rsuccess]=checksol(csmall,gamma,m,ell,K,xi,A,SG)
% takes the approximate solution csmall (in symmetry reduces variables)
% and computes the radii polynomial based on the approximate inverse A
% of the finite dimensional Jacobian.
% The physical parameter values are gamma, m and ell 
% (see the description of the Ohta-Kawasaki model in the paper).
% The integer K determines the size of the finite dimensional projection
% The symmetry group SG can be either 230 or 229.
% The default symmetry group is 230.
%
% Although the function works for the exponentially weighted norm
% sum(|c_k| xi^|k|) with |k|=sqrt(k_1^2+k_2^+k_3^2) on sequence space,
% the case xi=1 is the only one considered in the paper
%
% If m is an interval all computations will be rigorously verified
% using interval arithmetic in Intlab 
% (otherwise it is based on floating point arithmetic)
%
% If the radii polynomial is negative for some positive r, 
% then success=1 and rsuccess is the succesfull radius
% If the method is unsuccesfull then success=0 and rsuccess=NaN.

success=0;
rsuccess=NaN;

if nargin == 7
    SG = 230;
end;

% turn all data into intervals if m is an interval
if exist('intval','file') && isintval(m)
    disp('Including interval arithmetic')
    csmall=intval(csmall);
    gamma=intval(gamma);
    ell=intval(ell);
    xi=intval(xi);
    A=intval(A);
else
    disp('No intervals yet')
end

[kx,ky,kz]=ndgrid(-K:K);
% choose fundamental domain
ksym=(kx>=ky & ky>=kz & kz>=0);

% determine sum(alpha_g(k) , g in G)
symmetrizedcount=symmetrize(ksym,SG);
% determine the number of elements in the stabilizer
stabilizercount=abs(symmetrizedcount);
stabilizercount(stabilizercount==0)=1; % prevent division by zero

% determine number of elements in orbits (orbit-stabilizer theorem)
grouporder=symmetrizedcount(K+1,K+1,K+1);
orbitcount=grouporder./stabilizercount;

% find k that are not in Z_{triv} 
% within the ball of radius K 
% and exclude k=0
k2=kx.^2+ky.^2+kz.^2;
ksmall=(symmetrizedcount~=0 & k2<=K^2 & k2>0);
% find symmetry reduced indices 
ksymsmall=(ksym & ksmall);

c=ILzeros([2*K+1,2*K+1,2*K+1],csmall(1));

c(ksymsmall)=csmall;
% determine tensor representation 
c=symmetrize(c,SG);
c=c(:)./stabilizercount(:);
% reshape necessary to work with intlab
c=reshape(c,2*K+1,2*K+1,2*K+1);

lk=k2/ell^2;
if K>ell*max(gamma,sqrt(gamma))
    CK=(K/ell)^2/(gamma^(-2)*(K/ell)^4-(K/ell)^2+1);
else
    disp('K is too small; needs to be larger than ell*gamma');
    success = 0;
    return;
end

%% Additional counters and indices for larger tensors used for convolution products

fprintf('Collecting counters and indices. \n');
% D is for double the size
[Dkx,Dky,Dkz]=ndgrid(-2*K:2*K);
Dksym=(Dkx>=Dky & Dky>=Dkz & Dkz>=0);
Dk2=Dkx.^2+Dky.^2+Dkz.^2;

Dsymmetrizedcount=symmetrize(Dksym,SG);
Dstabilizercount=abs(Dsymmetrizedcount);
Dstabilizercount(Dstabilizercount==0)=1; % prevent division by zero
Dtildealpha=Dsymmetrizedcount./Dstabilizercount;
Dorbitcount=Dstabilizercount(2*K+1,2*K+1,2*K+1)./Dstabilizercount;

% DD iincludes k=0
% needed for products like c*c, where (c*c)_0 is nontrivial
DDksymball=(Dksym & Dsymmetrizedcount~=0 & Dk2<=(2*K)^2);
DDksymindexball=(1:length(find(DDksymball)))';
DDksymindex=zeros([4*K+1,4*K+1,4*K+1]);
DDksymindex(DDksymball)=DDksymindexball;
DDksymindex=abs(symmetrize(DDksymindex,SG))./Dstabilizercount;

% E is for triple the size
[Ekx,Eky,Ekz]=ndgrid(-3*K:3*K);
Eksym=(Ekx>=Eky & Eky>=Ekz & Ekz>=0);
Ek2=Ekx.^2+Eky.^2+Ekz.^2;
Elk=Ek2/ell^2;
Elambda=gamma^(-2)*Elk.^2-Elk+1;

Esymmetrizedcount=symmetrize(Eksym,SG);
Estabilizercount=abs(Esymmetrizedcount);
Estabilizercount(Estabilizercount==0)=1; % prevent division by zero
Etildealpha=Esymmetrizedcount./Estabilizercount;
Eorbitcount=Estabilizercount(3*K+1,3*K+1,3*K+1)./Estabilizercount;

Ekball=(Esymmetrizedcount~=0 & Ek2<=(3*K)^2 & Ek2>0);
Eksymball=(Eksym & Ekball);
Eksymindexball=(1:length(find(Eksymball)))';
Eksymindex=zeros([6*K+1,6*K+1,6*K+1]);
Eksymindex(Eksymball)=Eksymindexball;
Eksymindex=abs(symmetrize(Eksymindex,SG))./Estabilizercount;

% Eksymsmall keeps track of which indices are in the smaller ball
Eksymsmall=(Eksymball & Ek2<=K^2);
Eksymtail=(Eksymball & ~(Ek2<=K^2));
Eksymindextail=Eksymindex(Eksymtail);
Eksymindexsmall=Eksymindex(Eksymsmall);

% EE is for five times the size
[EEkx,EEky,EEkz]=ndgrid(-5*K:5*K);
EEksym=(EEkx>=EEky & EEky>=EEkz & EEkz>=0);
EEk2=EEkx.^2+EEky.^2+EEkz.^2;
EElk=EEk2/ell^2;
EElambda=gamma^(-2)*EElk.^2-EElk+1;

EEsymmetrizedcount=symmetrize(EEksym,SG);
EEstabilizercount=abs(EEsymmetrizedcount);
EEstabilizercount(EEstabilizercount==0)=1; % prevent division by zero
EEorbitcount=EEstabilizercount(5*K+1,5*K+1,5*K+1)./EEstabilizercount;

EEkball=(EEsymmetrizedcount~=0 & EEk2<=(5*K)^2 & EEk2>0);
EEksymball=(EEksym & EEkball);
EEksymindexball=(1:length(find(EEksymball)))';
EEksymindex=zeros([10*K+1,10*K+1,10*K+1]);
EEksymindex(EEksymball)=EEksymindexball;
EEksymindex=abs(symmetrize(EEksymindex,SG))./EEstabilizercount;

% Eksymsmall keeps track of which indices are in the smaller ball
EEksymsmall=(EEksymball & EEk2<=K^2);
EEksymtail=(EEksymball & ~(EEk2<=K^2));
EEksymindextail=EEksymindex(EEksymtail);
EEksymindexsmall=EEksymindex(EEksymsmall);

%% Norm used in the weights in the 1-norm

if exist('intval','file') && isintval(m)
    ktwo=sqrt(intval(k2));
    Ektwo=sqrt(intval(Ek2));
    EEktwo=sqrt(intval(EEk2));
    Dktwo=sqrt(intval(Dk2));
else
    ktwo=sqrt(k2);
    Ektwo=sqrt(Ek2);
    EEktwo=sqrt(EEk2);
    Dktwo=sqrt(Dk2);
end

% choose to use the 2-norm in the weights (irrelevant for xi=1)
kn=ktwo;  
Ekn=Ektwo;
EEkn=EEktwo;
Dkn=Dktwo;

% in case one wants to use the 1-norm in the weights:
% kn=abs(kx)+abs(ky)+abs(kz);  
% Ekn=abs(Ekx)+abs(Eky)+abs(Ekz);
% Dkn=abs(Dkx)+abs(Dky)+abs(Dkz);

%% Start of bounds %%


fprintf('Building F, J, Dcc and Eccc as intervals. \n');
[F,J,Eccc,Dcc]=FDF(csmall,gamma,m,ell,K,SG);
% we explicitely use the fact that we real Fourier coefficients here
Eccc=real(Eccc);

% abbreviations:
Ekta=Eksymtail;
Ekss=Eksymsmall;
EEkta=EEksymtail;
EEkss=EEksymsmall;
ks=ksymsmall;

display(['Verifying with gamma = ',num2str(ILsup(gamma)),...
    ' ; m = ', num2str(ILsup(m)), ' ; ell = ', num2str(ILsup(ell)),...
    ' and K = ',num2str(K)]);

%% Y bound %%

fprintf('Starting bounds. \n');
res=A*F;
Y1=abs(res).*xi.^(kn(ks)).*orbitcount(ks);
Y2=Elk(Ekta).*abs(Eccc(Ekta))./Elambda(Ekta).*xi.^(Ekn(Ekta)).*Eorbitcount(Ekta);
Y=sum(Y1)+sum(Y2);

%% Q1 bound %%

JM=length(find(ksymsmall));
IAJ=abs(eye(JM)-A*J);
v=(xi.^(kn(ks)).*orbitcount(ks))';
% norm of matrix IAJ with respect to symmetrized norm
Q1=max((v*IAJ)./v);

%% Q2 bound %%

%recall EEkta=EEksymtail and EEkss=EEksymsmall;
EEcombitail=1./EElambda(EEkta).*xi.^(EEkn(EEkta)).*EEorbitcount(EEkta);
EEcombismall=xi.^(EEkn(EEkss)).*EEorbitcount(EEkss);

normcc=sum(abs(Dcc(DDksymball)).*xi.^(Dkn(DDksymball)).*Dorbitcount(DDksymball));

% we use explicitly here that we have (x,y,z) reflection symmetries
% (exploiting that ksym lies in the first octant)
% to reduce the size of EEcc somewhat
EE11cc=ILzeros([11*K+1,11*K+1,11*K+1],m);
Mz=(11*K+1)^2;
My=11*K+1;

ES3c=3*K+1;
S3K2=ES3c+(-2*K:2*K);
EE11cc(S3K2,S3K2,S3K2)=Dcc;
ESx0=ES3c+EEkx(EEksymball);
ESy0=ES3c+EEky(EEksymball)-1;
ESz0=ES3c+EEkz(EEksymball)-1;

EEJM=length(find(EEksymball));
zerosEJM=ILzeros([EEJM,1],m);
EElkEEksymball=EElk(EEksymball);
Eksymsmallindex=find(Eksymsmall);
EJM=length(find(Eksymball));
maxcolumnnorm=0;

fprintf('  Work completed in first loop (of two):     ');
step=0;
for j=find(Eksymball)'
    step=step+1;
    if mod(step,100)==0
        % give an indication of how far along we are
        tt=floor(100*step/EJM);
        if tt < 10
            fprintf('\b\b%i%%',tt);
        else
            fprintf('\b\b\b%i%%',tt);
        end;
    end
    column=zerosEJM;
    for jj=find(Eksymindex==Eksymindex(j))'
        % uses shift property of the convolution 
        ESx=ESx0-Ekx(jj);
        ESy=ESy0-Eky(jj);
        ESz=ESz0-Ekz(jj);
        % ESx+My*ESy+Mz*ESz turns tensor indexing into linear indexing
        % so that the entire kj-th column of values corresponding to 
        % symmetry reduced indices can be extracted from the tensor EEcc
        column=column+Etildealpha(jj)*EE11cc(ESx+My*ESy+Mz*ESz);
    end
    
    column=column.*EElkEEksymball;
    columntail=column(EEksymindextail);
    columntail=abs(columntail).*EEcombitail;
    if ismember(j,Eksymsmallindex)
        % top left (the part absorbed into Q1) has no contribution in Q2
        columnnorm=sum(columntail)/(xi^(Ekn(j))*Eorbitcount(j));
        maxcolumnnorm=max(maxcolumnnorm,columnnorm);
    else
        columnsmall=column(EEksymindexsmall);
        columnsmall=abs(A*columnsmall).*EEcombismall;
        columnnorm=(sum(columnsmall)+sum(columntail))/(xi^(Ekn(j))*Eorbitcount(j));
        maxcolumnnorm=max(maxcolumnnorm,columnnorm);
    end
end
fprintf('\n');

Q2=3*max(maxcolumnnorm,CK*normcc);

%% Q3,1 bound %%

%recall Ekta=Eksymtail and Ekss=Eksymsmall;
Ecombitail=1./Elambda(Ekta).*xi.^(Ekn(Ekta)).*Eorbitcount(Ekta);
Ecombismall=xi.^(Ekn(Ekss)).*Eorbitcount(Ekss);

normc0=sum(abs(c(ks)).*xi.^(kn(ks)).*orbitcount(ks)); 
normc=normc0+abs(m);   % plug in c_0 explicitely

E8c=ILzeros([8*K+1,8*K+1,8*K+1],m);
My=8*K+1;
Mz=(8*K+1)^2;

ES2c=2*K+1;
S2K=ES2c+(-K:K);
E8c(S2K,S2K,S2K)=c;
ESx0=ES2c+Ekx(Eksymball);
ESy0=ES2c+Eky(Eksymball)-1;
ESz0=ES2c+Ekz(Eksymball)-1;

EJM=length(find(Eksymball));
zerosEJM=ILzeros([EJM,1],m);
ElkEksymball=Elk(Eksymball);
DJM=length(find(DDksymball));
maxcolumnnorm=0;

fprintf('  Work completed in second loop:     ');
step=0;
for j=find(DDksymball)'
        step=step+1;
    if mod(step,100)==0
        % give an indication of how far along we are
        tt=floor(100*step/DJM);
        if tt < 10
            fprintf('\b\b%i%%',tt);
        else
            fprintf('\b\b\b%i%%',tt);
        end;
    end
    column=zerosEJM;
    for jj=find(DDksymindex==DDksymindex(j))'
        % uses shift property of the convolution 
        DSx=ESx0-Dkx(jj);
        DSy=ESy0-Dky(jj);
        DSz=ESz0-Dkz(jj);
        column=column+Dtildealpha(jj)*E8c(DSx+My*DSy+Mz*DSz);
    end
    column=column.*ElkEksymball;
    columnsmall=column(Eksymindexsmall);
    columntail=column(Eksymindextail);
    columnsmall=abs(A*columnsmall).*Ecombismall;
    columntail=abs(columntail).*Ecombitail;
    columnnorm=(sum(columnsmall)+sum(columntail))/(xi^(Dkn(j))*Dorbitcount(j));
    maxcolumnnorm=max(maxcolumnnorm,columnnorm);
end
fprintf('\n');
Q31=6*max(maxcolumnnorm,CK*normc);

%% Q3,2 bound %%

Alk=abs(A)*diag(lk(ksymsmall));
% recall v=(xi.^(kn(ks)).*orbitcount(ks))';
normAlk=max((v*Alk)./v);
Q32=3*max([normAlk,CK]);

%% Compute radii polynomial and decision on successfulness %%

D=(Q1+Q2-1)^2-4*Y*Q31;
E=Q1+Q2-1;

if E<0 && D>0
    coeff=[Q32,Q31,E,Y];
    if exist('intval','file') && isintval(m)
        Y=ILsup(Y);
        E=ILsup(E);
        Q31=ILsup(Q31);
        Q32=ILsup(Q32);
        D=E^2-4*Y*Q31;
    end
    r0(1)= (-E-sqrt(D))/(2*Q31);
    r0(2)= (-E+sqrt(D))/(2*Q31);
    [r1,r2]=findroots([Q32,Q31,E,Y],r0);
    if r1>0 && r2>r1
        rr=(r1+r2)/2;
        if exist('intval','file') && isintval(m)
            rr=intval(rr);
        end
        result=sum(coeff.*rr.^(3:-1:0));
        if result<0
            success=1;
            rsuccess=ILsup(rr);
            display('Success!');
            display(['r = ',num2str(rsuccess),' ; the norm of the solution is ',num2str(-ILsup(-normc0))]);
            display(['Y = ',num2str(Y),' ; Q1+Q2-1 = ',num2str(E),' ; Q31 = ',num2str(Q31),' ; Q32 = ',num2str(Q32)]);
            display(['Only ',int2str(JM),' modes needed out of ',int2str((2*K+1)^3), ' possible.']);
            
        end
    end
    if ~success
        display('sorry, polynomial not negative')
        display(['Y = ',num2str(Y),' Q1+Q2-1 = ',num2str(E),' Q31 = ',num2str(Q31),' Q32 = ',num2str(Q32)]);
    end
else
    if ~(E<0)
        display('sorry, slope is positive')
        display(['Y = ',num2str(Y),' Q1+Q2-1 = ',num2str(E),' Q31 = ',num2str(Q31),' Q32 = ',num2str(Q32)]);
    else
        display('sorry, discriminant is negative')
        display(['Y = ',num2str(Y),' Q1+Q2-1 = ',num2str(E),' Q31 = ',num2str(Q31),' Q32 = ',num2str(Q32)]);
    end
end
    
end


