function ALL=testcostanza(mla)
% Finds parameters mm,ll where the ternary Hermite scheme of Costanza has maybe high regularity
% Input: 
%   Not needed
%
% Output: Text Output
%
% Written by: tommsch, 2018

if(nargin==0); 
    mla=zeros(0,3);
end
ALL={};

FAST=1;

%Initialize

k=100; %running idx;
MAX=0;
val=0;
it=0;
while(true)
    k=k+10;
    it=it+1;
    mla=sortrows(mla,3);
    %mmdiff=max(abs(mla(end,1))/k,abs(mla(end,1)-mla(end-1,1)));
    %mm=mla(end,1)+randn*mmdiff;
    try
        lldiff=max(abs(mla(end,2))/k,abs(mla(end,2)-mla(end-1,2)));
        ll=mla(end,2)+randn*lldiff;
    catch
        ll=randn/100+.276683836054697;
    end
    
    mm=-2.8125;
    %ll=randn/100+.278039805539948;
    
    if(val>MAX)
        fprintf('Maximum at \tmm= %.15g; \tll=%.15g; \talpha=%.15g;\n',mla(end,1),mla(end,2),mla(end,3));
        MAX=val;
        it=0;
    end
    
    %SCALE=50;
    %f = fit([mla(:,1), mla(:,2)],mla(:,3),'linearinterp'); plot( f )
    %%plot( f, [mla(:,1), mla(:,2)], mla(:,3))
    plot(mla(:,2),mla(:,3),'.')
    %stem3(mla(:,1),mla(:,2),mla(:,3),'o')
    
    %x0=mla(end,1); %xdiff=SCALE*abs(mla(end,1)-mla(end-1,1));
    %y0=mla(end,2); %ydiff=SCALE*abs(mla(end,2)-mla(end-1,2));
    %z0=mla(end,3); zdiff=SCALE*abs(mla(end,3)-mla(end-1,3));
    %axis([x0-mmdiff x0+mmdiff y0-lldiff y0+lldiff z0-zdiff z0+zdiff]);
    drawnow
    
    M=costanza_matrices(mm,ll,FAST);
    val=alpha(M);
    
%     try
%     if(val>MAX);
%         val=tjsr(costanza_matrices(mm,ll,FAST),'invariantsubspace','none','delta',.999999);
%         val=-log(val(2))/log(3);
%         if(val>MAX)
%             MAX=val;
%         end
%     end
%     catch
%         continue;
%     end
    mla(end+1,:)=[mm ll val]; %#ok<AGROW>
    %if(it==100)
    %    ALL{end+1}=mla; %#ok<AGROW>
    %    fprintf('\n\n==================================\n==RESTART========================\n============================\n\n');
    %    break;
    %end

    
end
    testcostanza %restart if algorithm stucks


end

function val=alpha(TV0);
    [c,~,val]=findsmp(TV0,'maxsmpdepth',20,'verbose',0);
    fprintf('%i ',c{1}');
    fprintf(' - ');
    val=val.jsrbound(1);
    val=-log(val)/log(3);
end


function dummy; end %#ok<DEFNU> %Generates an error, if the 'end' of a function is missing. 
