function [X,xs] = viss(nls,vars,appsol,tol)

[dflsys,dflsol] = ParaDfl(nls,vars,appsol,tol);

sys = myrep4inline(dflsys);
f = inline(sys);
[X,xs] = verifynlss(f,dflsol);

end

%%%%%%%%Parameterized  Deflation%%%%%%%%%

%function ParaDfl(pol,ord,mroots,tol)
function [dflsys,dflsol] = ParaDfl(pol,ord,mroots,tol)

%%%%%%%%Initialization%%%%%%%%

n = length(ord);%size of input
vars = myseq('x(i)','i',1:n);%variables
eqs = subs(pol,ord,vars);%equations
oprtzeros = [mroots(1:n)];
JacZero = subs(jacobian(eqs,vars),vars,mroots);%Jacobian matrix
[br,leftnull,rightnull] = AppNull(JacZero,tol);%corank
disp('ParaDfl:');
temp = (svd(JacZero));
disp(['  Breadth and smallest singular value of Jacobian: ',num2str(br),', ',num2str(temp(end))]);
% if br<n   
%     disp('ParaDfl:');
%     disp(['  Singular values around tolerance:',num2str(temp(end-br)),',',num2str(temp(end+1-br))]);
% else
%     disp('ParaDfl:');
%     disp(['  Singular values around tolerance:',num2str(temp(end+1-br))]);
% end
[rowind,colind] = CoRank(leftnull,rightnull);%rank-deficient information
numdfl = 0;%number of deflation
dfl = {eqs};%deflation system
nullvec = {};
perpos = [];
dflsz = n;%size of deflation system

%%%%%%%%Deflation%%%%%%%%

%%%%%%%%Special Case: Breadth One%%%%%%%%

if br==1
    appsol = mroots(:);
    for i = 1 : 4
        appsol = MultipleRootRefinerBreadthOne(eqs,vars,appsol,tol);%refine approximation to high accuracy
    end
    oprtzeros = appsol;
    JacZero = subs(jacobian(eqs,vars),vars,appsol);
    tmpneweqs = jacobian(eqs,vars(colind(1)));
    neweqs = tmpneweqs(:);
    compl = ComplSet(colind,n);
    while br==1
        numdfl = numdfl+1;
        %solve a least squares problem for new parameters
        tmpA = JacZero(1:end,compl(1:n-1));
        tmpb = -subs(neweqs,vars,oprtzeros);
        azeros = pinv(tmpA)*tmpb;
        oprtzeros = [oprtzeros;azeros];%add new zeros
        vars = [vars;myseq('x(i)','i',dflsz+1:dflsz+n-1)];%add new variables
        newvec = [];
        k = 1;
        for i = 1 : n %construct parameterized nullvector
            if i == colind(1)
                if numdfl==1 
                    newvec(i) = 1;
                    k = k+1;
                else
                    newvec(i) = 0;
                    k = k+1; 
                end
            else
                  newvec = symsplnv(newvec,i,myseq('x[i]','i',dflsz+i-k+1));
            end
        end
        newvec = newvec(:);
        nullvec = {nullvec{:},newvec};
        neweqs = neweqs+jacobian(eqs,vars(compl(1:n-1)))*newvec(compl(1:n-1));
        dflsz = dflsz+n-1;%add n-1 new parameters every deflation
        tmpep = expand(neweqs(1:n));
        dfl{length(dfl)+1} = tmpep(:).';
        neweqs = zeros(n,1);
        for i = 1 : numdfl%criterion for next order derivative
            neweqs = neweqs + i/(numdfl+1)*jacobian(dfl{numdfl-i+2},vars(1:n))*nullvec{i};
        end
        tempMatrix = [subs(neweqs,vars,oprtzeros),JacZero(1:end,compl(1:n-1))];
        [br,leftnull,rightnull]= AppNull(tempMatrix,tol);
        if br>0
            tmp = [subs(neweqs,vars,oprtzeros),JacZero(1:end,compl(1:n-1))];
            tmpsvd = svd(tmp);
%             disp('ParaDfl:');
%             disp(['  Singular values around tolerance: ',num2str(tmpsvd(end-br)),' , ',num2str(tmpsvd(end-br+1))]); 
        end
    end
    tmpsvd = svd([subs(neweqs,vars,oprtzeros),JacZero(1:end,compl(1:n-1))]);
    disp('ParaDfl:');
    disp(['  Number of deflation and smallest singular value of deflated system: ',num2str(numdfl),', ',num2str(tmpsvd(end))]);
    newper = 0;
    for i = 0 : numdfl-1
        newx = myseq('x[i]','i',dflsz+1+i);
        newper = newper+newx*vars(colind(1))^i;
    end
    for i = 1 : numdfl
        dfl{i}(rowind(1)) = dfl{i}(rowind(1))+newper;
        newper = 1/i*diff(newper,vars(colind(1)));
    end
    oprtzeros = [oprtzeros;zeros(numdfl,1)];
    dflsys = [dfl{:}];
    dflsol = oprtzeros;
    
%%%%%%%%General Case%%%%%%%%
    
else
    while br > 0 
        for i = 1 : br 
            tmpx = myseq('x[j]','j',dflsz+i);
            dfl{1}(rowind(i)) = dfl{1}(rowind(i))+(1/factorial(numdfl))*tmpx*vars(colind(i))^numdfl; %add perturbations to original system	
        end
		numdfl = numdfl+1;
        vars = [vars;myseq('x[i]','i',dflsz+1:dflsz+br)];%add new smoothing parameters
        perpos = [perpos,dflsz+1:dflsz+br];
        oprtzeros = [oprtzeros,zeros(1,br)];
        ranvec = ones(br,1);
        newvec = zeros(dflsz,1);
		k = 1;
		for i = 1 : n %construct parameterized nullvector
			if  k<=br && i==colind(k)
                newvec(i) = ranvec(k); 
                k = k+1; 
            else
                newvec = symsplnv(newvec,i,myseq('x[i]','i',dflsz+br+i+1-k));
            end
        end
		for i = n+1 : dflsz 
            newvec = symsplnv(newvec,i,myseq('x[i]','i',dflsz+i));
        end
        vars = [vars;myseq('x[i]','i',dflsz+br+1:2*dflsz)];
        nullvec = {nullvec{:},newvec};
		%solve a least squares problem for new parameters
		compl = ComplSet(colind,dflsz);
        tmpA = JacZero(1:end,compl(1:dflsz-br));
        tmpb = -JacZero(1:end,colind(1:br))*ranvec;
        azeros = pinv(tmpA)*tmpb;
        tmpazeros = azeros.';
        oprtzeros = [oprtzeros,tmpazeros(1:dflsz-br)];
		for i = 1 : numdfl 
            k = length(nullvec{i});
			if i < numdfl 
                tmpJacobian = jacobian(dfl{i},vars(1:k))*nullvec{i};
                tmpJacobian = tmpJacobian.';
                dfl{i+1} =[dfl{i},tmpJacobian(1:k)];
            else
                tmpJacobian = jacobian(dfl{i},vars(1:k))*nullvec{i};
                tmpJacobian = tmpJacobian.';
                dfl = {dfl{:},[dfl{i},tmpJacobian(1:k)]};
			end
		end
		dflsz = 2*dflsz;
        JacZero = subs(jacobian(dfl{end},vars),vars,oprtzeros);
		[br,leftnull,rightnull] = AppNull(JacZero,tol);
		if br>0 
			[rowind,colind] = CoRank(leftnull(end+1-n:end,1:end),rightnull(1:end,1:n));
% 			if br<n
%                 tmpsvd = svd(JacZero);
%                 disp('ParaDfl:');
%                 disp(['  Singular values around tolerance:',num2str(tmpsvd(end-br)),' , ',num2str(tmpsvd(end))]);
%             else
%                 tmpsvd = svd(JacZero);
%                 disp('ParaDfl:');
%                 disp(['  Singular values around tolerance:',num2str(tmpsvd(end-br+1))]);
% 			end
        end
        disp('ParaDfl:');
        disp(['  Number of deflation and breadth: ',num2str(numdfl),', ',num2str(br)]);
    end
    tmpsvd = svd(JacZero);
    disp('ParaDfl:');
    disp(['  Number of deflation and smallest singular value of deflated system: ',num2str(numdfl),', ',num2str(tmpsvd(end))]);
    dflsys = dfl{end};
    dflsol = oprtzeros(1:dflsz);
end
    
end

%%%%%%%%Refine Mutiple Roots Breadth One%%%%%%%%
 
% MultipleRootRefinerBreadthOne=proc(lsys,ord,approot,tol)
function  appsol = MultipleRootRefinerBreadthOne(lsys,ord,approot,tol)

n = length(ord);
appsol = approot;
Jac = jacobian(lsys,ord);
Jaclhs = subs(Jac,ord,appsol);
% disp('MultipleRootRefinerBreadthOne:');
% disp('  singular values of jacobian');
%disp(svd(Jaclhs));
tsvdJaclhs=svd(Jaclhs);
A=Jaclhs'*Jaclhs+tsvdJaclhs(n)*sparse(eye(n));
b=(-Jaclhs')*(subs(lsys,ord,appsol).');
appy = pinv(A)*b; %最小二乘法的暂时代替方法
appsol = appsol + appy;
% disp('MultipleRootRefinerBreadthOne:');
% disp('  roots after prenewton');
% disp(appsol);
Jaclhs = subs(Jac,ord,appsol);
[U,tsvdJaclhs,V]=svd(Jaclhs);
flag = tsvdJaclhs(n);
Vt = V';
temp = Vt(end,1:end)';%共轭转置
transtemp = temp.';%为了下步行优先转换向量，转置temp
nullvector = transtemp(:);
maindir = maindirection(nullvector);
maindir = maindir(2);
tmp = maindirection(nullvector);
tmp_1 = tmp(1);
nullvector = nullvector/tmp_1;
neword = exchange(ord,ord(1),ord(maindir));
newL = exchange(Jac*nullvector,ord(1),ord(maindir));
nlsys = exchange(lsys,ord(1),ord(maindir));
Jac = jacobian(nlsys,ord);
Jac = Jac(1:end,2:end);
Jaclhs = subs(Jac,neword,appsol);
if maindir == 1 
    A = nullvector(2:end).';
else
    A = nullvector([2:maindir-1,1,maindir+1:end]).';
end
L = newL;
multiplicity=2;
while flag < tol
    P = diff(newL,ord(1));
    Q = L*A;
    for i = 1 : n-1
        P = P + diff(Q(1:end,i),ord(i+1));
    end
    Jacrhs = subs(P,neword,appsol);
    tmpflag = svd([Jacrhs,Jaclhs]);
    flag = tmpflag(n);
%     disp('MultipleRootRefinerBreadthOne:');
%     disp(['  The smallest singular value: ',num2str(flag)]);
    if flag > tol 
        break;
    else
        [U,tsvdJaclhs,V] = svd([Jacrhs,Jaclhs]);
        Vt = V';
        newA = (Vt(end,1:end))';
        newA = newA(:);
        newA = newA(2:end)/newA(1);
        newL = (P + Jac*newA)/multiplicity;
        A = [newA.';A];
        L = [L,newL];
        multiplicity = multiplicity + 1;
%         disp('MultipleRootRefinerBreadthOne:');
%         disp(['  The multiplicity',',',num2str(multiplicity)]);
    end
end
% disp('MultipleRootRefinerBreadthOne:');
% disp(['  The multiplicity',',',num2str(multiplicity)]);
Jacrhs = subs(P,neword,appsol)/multiplicity;
A = [Jacrhs,Jaclhs];
b = -subs(newL,neword,appsol);
tmplsq = pinv(A)*b;
appsol = appsol + tmplsq(1)/multiplicity*nullvector;

end

%%%%%%%%交换变元%%%%%%%%

function new = exchange(p,para1,para2)
p = p(:);%将输入换成纵向的，原因;如果不换，在2011等版本会出问题，而2008不会
tmp = subs(p,para1,sym('t'));
tmp = subs(tmp,para2,para1);
new = subs(tmp,sym('t'),para2);
end

%%%%%%%%Approximate Nullspace%%%%%%%%

function [br,leftnull,rightnull] = AppNull(A,tol)
n = size(A,2);
[leftnull,SS,rightnull] = svd(A);
rightnull = rightnull';
SS = SS*ones(size(SS,2),1);
br = n;
for i = n :-1: 1
    if SS(i) > tol 
        br = n-i; 
        break;%count how many singular values are smaller than tol
    end
end
if br > 0
    leftnull = leftnull(1:end,end+1-br:end);
    rightnull = rightnull(end+1-br:end,1:end);
    return;
else
    leftnull = [];
    rightnull = [];
    return;
end
end

%%%%%%%%Complement%%%%%%%%

function C = ComplSet(A,n)
C = [];
k = 1;
for i = 1 : n
    if k<=length(A) && i==A(k) 
        k = k+1;
    else
        C = [C,i];
    end
end
end

%%%%%%%%Deficiency Information%%%%%%%%

function [rowind,colind] = CoRank(A,B)
n = size(A,1);
br = size(A,2);
rowind = zeros(1,br);
colind = zeros(1,br);
leftnull = A;
rightnull = B;
for i = br :-1: 1 
    lflag = 0;
    rflag = 0;
    for j = 1 : n
        if lflag < abs(leftnull(j,i)) 
            lflag = abs(leftnull(j,i)); 
            rowind(i) = j; 
        end
        if rflag < abs(rightnull(i,j)) 
           rflag = abs(rightnull(i,j)); 
           colind(i) = j;
        end
    end
    if i > 1
        for k = 1 : i-1
            for j = 1 : n
                leftnull(j,k) = leftnull(j,k)-leftnull(j,i)*leftnull(rowind(i),k)/leftnull(rowind(i),i);
                rightnull(k,j) = rightnull(k,j)-rightnull(i,j)*rightnull(k,colind(i))/rightnull(i,colind(i));
            end
        end
    end
end
colind = sort(colind);
rowind = sort(rowind);
end


%%%%%%%%Main Direction%%%%%%%%

function maindir = maindirection(nullvector)
temp = 0;
maindir = 1;
for i = 1 : length(nullvector)
    if abs(nullvector(i))>temp
        temp = abs(nullvector(i));
        maindir = i;
    end
end
maindir = [nullvector(maindir),maindir];

end  

%%%%%%%%用于实现数值向量中替换符号变元%%%%%%%%

function vec = symsplnv(numvec,i,symvar)
tmp = sym(numvec);
tmp(i) = symvar;
vec = tmp;
end


%%%%%%%%%将dflsys转换为verifynlss所需格式%%%%%%%%

function str = myrep4inline(e)
e = vpa(e);
newstrarr = myrep(char(e));
l = length(newstrarr);
for i = 1 : l
    if 1 == i
        str = [char(newstrarr(i))];
    else
        str = [str,';',char(newstrarr(i))];
    end
end
str = ['[',str,']'];
end

function newstr = myrep(oldstr)
pattern = 'x[\d]*';
oldarr = regexpi(oldstr,pattern,'match');
newarr = sprep(oldarr);
newstr = subs(oldstr,oldarr,newarr);

end

function new = sprep(old)
new = strrep(old,'x','x(');
new = strcat(new,')');
end

%%%%%%%%%%%%%%用于生成sym变元%%%%%%%%%%%%%

function vars = myseq(template,ind,vec)
% e.g
% Input : 
% template = 'x(i)'
% ind = 'i'
% vec = 1 : n

% Output :
% vars = [x1,x2,x3,...,xn]
if (~strcmp(class(template),'char') || ~strcmp(class(ind),'char'))
    error('Error,wrong type of template or ind');
end
workTemplate = template;
old = {'(',')','[',']'};
new = {''};
vars = [];
for k = 1 : length(old)
    workTemplate = strrep(workTemplate,old{k},new{1});
end
for k = 1 : length(vec)
    varv = strrep(workTemplate,ind,num2str(vec(k)));
    vars = [vars;sym(varv)];
end
end


