text:
%%*************************************************************************
%% sqlp: main solver
%%
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

  function [obj,X,y,Z,info,runhist] = sqlpmain(blk,At,C,b,par,parbarrier,X0,y0,Z0);

   global spdensity smallblkdim  printlevel msg
   global solve_ok  use_LU  exist_analytic_term numpertdiagschur
   global schurfun  schurfun_par
%%
   randstate = rand('state');  randnstate = randn('state');
   rand('state',0);   randn('state',0);
%%
   vers          = par.vers;
   predcorr      = par.predcorr;
   gam           = par.gam;
   expon         = par.expon;
   gaptol        = par.gaptol;
   inftol        = par.inftol;
   steptol       = par.steptol;
   maxit         = par.maxit;
   printlevel    = par.printlevel;
   stoplevel     = par.stoplevel;
   scale_data    = par.scale_data;
   spdensity     = par.spdensity;
   rmdepconstr   = par.rmdepconstr;
   smallblkdim   = par.smallblkdim;
   schurfun      = par.schurfun;
   schurfun_par  = par.schurfun_par;
   ublksize      = par.ublksize;
%%
   tstart = clock;
   X = X0; y = y0; Z = Z0;
   for p = 1:size(blk,1)
      if strcmp(blk{p,1},'u'); Z{p} = zeros(blk{p,2},1); end
   end
%%
%%-----------------------------------------
%% convert unrestricted blk to linear blk.
%%-----------------------------------------
%%
   convertlen = 0;
   [blk,At,C,X,Z,u2lblk,ublkidx] = sqlpu2lblk(blk,At,C,X,Z,par,convertlen);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if (u2lblk(p) == 1)
         n = 2*blk{p,2};
         blk{p,1} = 'l';  blk{p,2} = n;
         parbarrier{p} = zeros(1,n);
         At{p} = [At{p}; -At{p}];
         tau = max(1,norm(C{p}));
         C{p} = [C{p}; -C{p}];
         msg = 'convert ublk to lblk';
         if (printlevel); fprintf(' *** %s',msg); end
         b2 = 1 + abs(b');
         normCtmp = 1+norm(C{p});
         normAtmp = 1+sqrt(sum(At{p}.*At{p}));
         if (n > 1000)
            const = sqrt(n);
         else
	    const = n;
         end
         if (par.startpoint == 1)
            X{p} = const* max([1,b2./normAtmp]) *ones(n,1);
            Z{p} = const* max([1,normAtmp/sqrt(n),normCtmp/sqrt(n)]) *ones(n,1);
            X{p} = X{p}.*(1+1e-10*rand(n,1));
            Z{p} = Z{p}.*(1+1e-10*rand(n,1));
	 else
            const = max(abs(X{p})) + 100;
            X{p} = [X{p}+const; const*ones(n/2,1)];
            %%old: const = 100; Z{p} = [const*ones(n/2,1); const*ones(n/2,1)];
            Z{p} = [abs(Z0{p}); abs(Z0{p})] + 1e-4;
         end
      end
   end
%%-----------------------------------------
%% check whether {A1,...,Am} is
%% linearly independent.
%%-----------------------------------------
%%
   m0 = length(b);
   [At,b,y,indeprows,par.depconstr,feasible,par.AAt] = ...
    checkdepconstr(blk,At,b,y,rmdepconstr);
   if (~feasible)
      obj = []; X = cell(size(blk,1),1); y = []; Z = cell(size(blk,1),1);
      runhist = [];
      msg = 'SQLP is not feasible';
      if (printlevel); fprintf('\n %s \n',msg); end
      return;
   end
   par.normAAt = norm(par.AAt,'fro');
%%
%%-----------------------------------------
%% scale SQLP data. Note: must be done only
%% after checkdepconstr
%%-----------------------------------------
%%
   normA2 = 1+ops(At,'norm');
   normb2 = 1+norm(b);
   normC2 = 1+ops(C,'norm');
   normX0 = 1+ops(X0,'norm');
   normZ0 = 1+ops(Z0,'norm');
   if (scale_data)
      [At,C,b,normA,normC,normb,X,y,Z] = scaling(blk,At,C,b,X,y,Z);
   else
      normA = 1; normC = 1; normb = 1;
   end
%%
%%-----------------------------------------
%% find the combined list of non-zero
%% elements of Aj, j = 1:k, for each k.
%% IMPORTANT NOTE: Ak, C are permuted.
%%-----------------------------------------
%%
   par.numcolAt = length(b);
   [At,C,X,Z,par.permA,par.permZ] = sortA(blk,At,C,b,X,Z);
   [par.isspA,par.nzlistA,par.nzlistAsum,par.isspAy,par.nzlistAy] = nzlist(blk,At,par);
%%
%%-----------------------------------------
%% create an artifical non-negative block
%% for a purely log-barrier problem
%%-----------------------------------------
%%
   numblkold = size(blk,1);
   nn = 0;
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      idx = find(parbarrier{p}==0);
      if ~isempty(idx);
         if strcmp(pblk{1},'l')
            nn = nn + length(idx);
         elseif strcmp(pblk{1},'s') | strcmp(pblk{1},'q')
            nn = nn + sum(pblk{2}(idx));
         end
      end
   end
   if (nn==0)
      analytic_prob = 1;
      numblk = size(blk,1)+1;
      blk{numblk,1} = 'l'; blk{numblk,2} = 1;
      At{numblk,1} = sparse(1,length(b));
      C{numblk,1} = 1;
      X{numblk,1} = 1e3;
      Z{numblk,1} = 1e3;
      parbarrier{numblk,1} = 0;
      u2lblk(numblk,1) = 0;
      nn = nn + 1;
   else
      analytic_prob = 0;
   end
%%
   exist_analytic_term = 0;
   for p = 1:size(blk,1);
      idx = find(parbarrier{p} > 0);
      if ~isempty(idx);
         exist_analytic_term = 1;
      end
   end
%%-----------------------------------------
%% initialization
%%-----------------------------------------
%%
   EE = ops(blk,'identity');
   normE2 = ops(EE,'norm'); Zpertold = 1;
   for p = 1:size(blk,1)
      normCC(p) = 1+ops(C(p),'norm');
      normEE(p) = 1+ops(EE(p),'norm');
   end
   [Xchol,indef(1)] = blkcholfun(blk,X);
   [Zchol,indef(2)] = blkcholfun(blk,Z);
   if any(indef)
      msg = 'stop: X or Z not positive definite';
      if (printlevel); fprintf('\n  %s\n',msg); end
      info.termcode = -3;
      info.msg1 = msg;
      obj = []; X = cell(size(blk,1),1); y = []; Z = cell(size(blk,1),1);
      runhist = [];
      return;
   end
   AX = AXfun(blk,At,par.permA,X);
   rp = b-AX;
   ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
   ZpATynorm = ops(ZpATy,'norm');
   Rd = ops(C,'-',ZpATy);
   objadd0 = 0;
   if (scale_data)
      for p = 1:size(blk,1)
         pblk = blk(p,:);
         objadd0 = objadd0 + sum(parbarrier{p}.*pblk{2})*log(normA{p});
      end
   end
   objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
   obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
   gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
   relgap = gap/(1+sum(abs(obj)));
   prim_infeas = norm(rp)/normb2;
   dual_infeas = ops(Rd,'norm')/normC2;
   infeas = max(prim_infeas,dual_infeas);
   if (scale_data)
      infeas_org(1) = prim_infeas*normb;
      infeas_org(2) = dual_infeas*normC;
   else
      infeas_org = [0,0];
   end
   trXZ = blktrace(blk,X,Z,parbarrier);
   if (nn > 0); mu  = trXZ/nn; else; mu = gap/ops(X,'getM'); end
   normX = ops(X,'norm');
%%
   termcode = 0; restart = 0;
   pstep = 1; dstep = 1; pred_convg_rate = 1; corr_convg_rate = 1;
   prim_infeas_min  = prim_infeas;
   dual_infeas_min  = dual_infeas;
   prim_infeas_best = prim_infeas;
   dual_infeas_best = dual_infeas;
   infeas_best = infeas;
   relgap_best = relgap;
   homRd = inf; homrp = inf; dy = zeros(length(b),1);
   msg = []; msg2 = []; msg3 = [];
   runhist.pobj    = obj(1);
   runhist.dobj    = obj(2);
   runhist.gap     = gap;
   runhist.relgap  = relgap;
   runhist.pinfeas = prim_infeas;
   runhist.dinfeas = dual_infeas;
   runhist.infeas  = infeas;
   runhist.step    = 0;
   runhist.normX   = normX;
   runhist.cputime = etime(clock,tstart);
   ttime.preproc   = runhist.cputime;
   ttime.pred = 0; ttime.pred_pstep = 0; ttime.pred_dstep = 0;
   ttime.corr = 0; ttime.corr_pstep = 0; ttime.corr_dstep = 0;
   ttime.pchol = 0; ttime.dchol = 0; ttime.misc = 0;
%%
%%-----------------------------------------
%% display parameters and initial info
%%-----------------------------------------
%%
   if (printlevel >= 2)
      fprintf('\n********************************************');
      fprintf('***********************\n');
      fprintf('   SDPT3: Infeasible path-following algorithms');
      fprintf('\n********************************************');
      fprintf('***********************\n');
      [hh,mm,ss] = mytime(ttime.preproc);
      if (printlevel>=3)
         fprintf(' version  predcorr  gam  expon  scale_data\n');
         if (vers == 1); fprintf('   HKM '); elseif (vers == 2); fprintf('    NT '); end
         fprintf('     %1.0f      %4.3f',predcorr,gam);
         fprintf('   %1.0f        %1.0f    %1.0f\n',expon,scale_data);
         fprintf('\nit pstep dstep pinfeas dinfeas  gap')
         fprintf('      mean(obj)   cputime\n');
         fprintf('------------------------------------------------');
         fprintf('-------------------\n');
         fprintf('%2.0f|%4.3f|%4.3f|%2.1e|%2.1e|',0,0,0,prim_infeas,dual_infeas);
         fprintf('%2.1e|%- 7.6e| %s:%s:%s|',gap,mean(obj),hh,mm,ss);
      end
   end
%%
%%---------------------------------------------------------------
%% start main loop
%%---------------------------------------------------------------
%%
   param.termcode    = termcode;
   param.iter        = 0;
   param.obj         = obj;
   param.relgap      = relgap;
   param.prim_infeas = prim_infeas;   param.dual_infeas = dual_infeas;
   param.homRd       = homRd;         param.homrp       = homrp;
   param.AX          = AX;            param.ZpATynorm   = ZpATynorm;
   param.normA       = normA;
   param.normb       = normb;         param.normC       = normC;
   param.normX0      = normX0;        param.normZ0      = normZ0;
   param.m0          = m0;            param.indeprows   = indeprows;
   param.prim_infeas_bad = 0;
   param.dual_infeas_bad = 0;
   param.prim_infeas_min = prim_infeas;
   param.dual_infeas_min = dual_infeas;
   param.gaptol      = gaptol;
   param.inftol      = inftol;
   param.maxit       = maxit;
   param.scale_data  = scale_data;
   param.printlevel  = printlevel;
   param.ublksize    = ublksize;
   Xbest = X; ybest = y; Zbest = Z;
%%
   for iter = 1:maxit;
      tstart  = clock;
      timeold = tstart;
      update_iter = 0; breakyes = 0;
      pred_slow = 0; corr_slow = 0; step_short = 0;
      par.parbarrier = parbarrier;
      par.iter    = iter;
      par.obj     = obj;
      par.relgap  = relgap;
      par.pinfeas = prim_infeas;
      par.dinfeas = dual_infeas;
      par.rp      = rp;
      par.y       = y;
      par.dy      = dy;
      par.normX   = normX;
      par.ZpATynorm = ZpATynorm;
      %%if (printlevel > 2); fprintf(' %2.1e',par.normX); end
      if (iter == 1 | restart); Cpert = min(1,normC2/ops(EE,'norm')); end
      if (runhist.dinfeas(1) > 1e-3) & (~exist_analytic_term) ...
         & (relgap > 1e-4)
         if (par.normX > 5e3 & iter < 20)
            Cpert = Cpert*0.5;
         elseif (par.normX > 5e2 & iter < 20);
            Cpert = Cpert*0.3;
         else;
            Cpert = Cpert*0.1;
         end
         Rd = ops(Rd,'+',EE,Cpert);
         %%if (printlevel > 2); fprintf('|%2.1e',Cpert); end
      end
%%---------------------------------------------------------------
%% predictor step.
%%---------------------------------------------------------------
%%
      if (predcorr)
         sigma = 0;
      else
         sigma = 1-0.9*min(pstep,dstep);
         if (iter == 1); sigma = 0.5; end;
      end
      sigmu = cell(size(blk,1),1);
      for p = 1:size(blk,1)
         sigmu{p} = max(sigma*mu, parbarrier{p}');
      end
      invXchol = cell(size(blk,1),1);
      invZchol = ops(Zchol,'inv');
      if (vers == 1);
         [par,dX,dy,dZ,coeff,L,hRd] = ...
          HKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol);
      elseif (vers == 2);
         [par,dX,dy,dZ,coeff,L,hRd] = ...
          NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);
      end
      if (solve_ok <= 0)
         msg = 'stop: difficulty in computing predictor directions';
         if (printlevel); fprintf('\n  %s',msg); end
         runhist.pinfeas(iter+1) = runhist.pinfeas(iter);
         runhist.dinfeas(iter+1) = runhist.dinfeas(iter);
         runhist.relgap(iter+1)  = runhist.relgap(iter);
         runhist.cputime(iter+1) = etime(clock,tstart);
         termcode = -4;
         break; %% do not ues breakyes = 1
      end
      timenew = clock;
      ttime.pred = ttime.pred + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%% step-lengths for predictor step
%%-----------------------------------------
%%
      if (gam == 0)
         gamused = 0.9 + 0.09*min(pstep,dstep);
      else
         gamused = gam;
      end
      [Xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);
      pstep = min(1,gamused*full(Xstep));
      timenew = clock;
      ttime.pred_pstep = ttime.pred_pstep + etime(timenew,timeold); timeold = timenew;
      Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
      dstep = min(1,gamused*full(Zstep));
      trXZnew = trXZ + pstep*blktrace(blk,dX,Z,parbarrier) ...
                 + dstep*blktrace(blk,X,dZ,parbarrier) ...
                 + pstep*dstep*blktrace(blk,dX,dZ,parbarrier);
      if (nn > 0); mupred  = trXZnew/nn; else; mupred = 1e-16; end
      mupredhist(iter) = mupred;
      timenew = clock;
      ttime.pred_dstep = ttime.pred_dstep + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for predictor step.
%%-----------------------------------------
%%
      if (min(pstep,dstep) < steptol) & (stoplevel) & (iter > 10)
         msg = 'stop: steps in predictor too short';
         if (printlevel)
            fprintf('\n  %s',msg);
            fprintf(': pstep = %3.2e,  dstep = %3.2e\n',pstep,dstep);
         end
         runhist.cputime(iter+1) = etime(clock,tstart);
         termcode = -2;
         breakyes = 1;
      end
      if (~predcorr)
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            pred_slow = all(mupredhist(idx)./mupredhist(idx-1) > 0.4);
            idx = [max(2,iter-5) : iter];
            pred_convg_rate = mean(mupredhist(idx)./mupredhist(idx-1));
            pred_slow = pred_slow + (mupred/mu > 5*pred_convg_rate);
         end
         if (max(mu,infeas) < 1e-6) & (pred_slow) & (stoplevel)
            msg = 'stop: lack of progress in predictor';
            if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mupred/mu = %3.2f, pred_convg_rate = %3.2f.',...
               mupred/mu,pred_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%---------------------------------------------------------------
%% corrector step.
%%---------------------------------------------------------------
%%
      if (predcorr) & (~breakyes)
         step_pred = min(pstep,dstep);
         if (mu > 1e-6)
            if (step_pred < 1/sqrt(3));
               expon_used = 1;
            else
               expon_used = max(expon,3*step_pred^2);
            end
         else
            expon_used = max(1,min(expon,3*step_pred^2));
         end
         if (nn==0)
             sigma = 0.2;
         elseif (mupred < 0)
             sigma = 0.8;
         else
            sigma = min(1, (mupred/mu)^expon_used);
         end
         sigmu = cell(size(blk,1),1);
         for p = 1:size(blk,1)
            sigmu{p} = max(sigma*mu, parbarrier{p}');
         end
         if (vers == 1)
            [dX,dy,dZ] = HKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         elseif (vers == 2)
            [dX,dy,dZ] = NTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         end
         if (solve_ok <= 0)
            msg = 'stop: difficulty in computing corrector directions';
            if (printlevel); fprintf('\n  %s',msg); end
            runhist.pinfeas(iter+1) = runhist.pinfeas(iter);
            runhist.dinfeas(iter+1) = runhist.dinfeas(iter);
            runhist.relgap(iter+1)  = runhist.relgap(iter);
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -4;
            break; %% do not ues breakyes = 1
         end
         timenew = clock;
         ttime.corr = ttime.corr + etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------
%% step-lengths for corrector step
%%-----------------------------------
%%
         if (gam == 0)
            gamused = 0.9 + 0.09*min(pstep,dstep);
         else
            gamused = gam;
         end
         Xstep = steplength(blk,X,dX,Xchol,invXchol);
         pstep = min(1,gamused*full(Xstep));
         timenew = clock;
         ttime.corr_pstep = ttime.corr_pstep+etime(timenew,timeold); timeold = timenew;
         Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
         dstep = min(1,gamused*full(Zstep));
         trXZnew = trXZ + pstep*blktrace(blk,dX,Z,parbarrier) ...
                    + dstep*blktrace(blk,X,dZ,parbarrier)...
                    + pstep*dstep*blktrace(blk,dX,dZ,parbarrier);
         if (nn > 0); mucorr  = trXZnew/nn; else; mucorr = 1e-16; end
         timenew = clock;
         ttime.corr_dstep = ttime.corr_dstep+etime(timenew,timeold); timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for corrector step
%%-----------------------------------------
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            corr_slow = all(runhist.gap(idx)./runhist.gap(idx-1) > 0.8);
            idx = [max(2,iter-5) : iter];
            corr_convg_rate = mean(runhist.gap(idx)./runhist.gap(idx-1));
            corr_slow = corr_slow + (mucorr/mu > max(min(1,5*corr_convg_rate),0.8));
         end
	 if (max(relgap,infeas) < 1e-6) & (iter > 20) ...
            & (corr_slow > 1) & (stoplevel)
            msg = 'stop: lack of progress in corrector';
   	    if (printlevel)
               fprintf('\n  %s',msg);
               fprintf(': mucorr/mu = %3.2f, corr_convg_rate = %3.2f',...
               mucorr/mu,corr_convg_rate);
            end
            runhist.cputime(iter+1) = etime(clock,tstart);
            termcode = -2;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%---------------------------------------------------------------
%% udpate iterate
%%---------------------------------------------------------------
      indef = [1,1];
      if (update_iter)
         for t = 1:5
            [Xchol,indef(1)] = blkcholfun(blk,ops(X,'+',dX,pstep));
            timenew = clock;
            ttime.pchol = ttime.pchol + etime(timenew,timeold); timeold = timenew;
            if (indef(1)); pstep = 0.8*pstep; else; break; end
         end
	 if (t > 1); pstep = gamused*pstep; end
	 for t = 1:5
            [Zchol,indef(2)] = blkcholfun(blk,ops(Z,'+',dZ,dstep));
            timenew = clock;
            ttime.dchol = ttime.dchol + etime(timenew,timeold); timeold = timenew;
            if (indef(2)); dstep = 0.8*dstep; else; break; end
         end
	 if (t > 1); dstep = gamused*dstep; end
         %%-------------------------------------------
         AXtmp = AX + pstep*AXfun(blk,At,par.permA,dX);
         prim_infeasnew = norm(b-AXtmp)/normb2;
         if (relgap < 5*infeas); alpha = 1e2; else; alpha = 1e3; end
         if any(indef)
            if indef(1); msg = 'stop: X not positive definite'; end
            if indef(2); msg = 'stop: Z not positive definite'; end
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = -3;
            breakyes = 1;
         elseif (prim_infeasnew > max([1e-8,relgap,20*prim_infeas]) & iter > 10) ...
            | (prim_infeasnew > max([1e-7,1e3*prim_infeas,0.1*relgap]) & relgap < 1e-2) ...
            | (prim_infeasnew > alpha*max([1e-9,param.prim_infeas_min]) ...
               & (prim_infeasnew > max([3*prim_infeas,0.1*relgap])) ...
               & (iter > 25) & (dual_infeas < 1e-6) & (relgap < 0.1)) ...
            | ((prim_infeasnew > 1e3*prim_infeas & prim_infeasnew > 1e-12) ...
               & (max(relgap,dual_infeas) < 1e-8))
            if (stoplevel)
               msg = 'stop: primal infeas has deteriorated too much';
               if (printlevel); fprintf('\n  %s, %2.1e',msg,prim_infeasnew); end
               termcode = -7;
               breakyes = 1;
            end
         elseif (trXZnew > 1.05*runhist.gap(iter)) & (~exist_analytic_term) ...
	    & ((infeas < 1e-5) & (relgap < 1e-4) & (iter > 20) ...
	       | (max(infeas,relgap) < 1e-7) & (iter > 10))
            if (stoplevel)
               msg = 'stop: progress in duality gap has deteriorated';
               if (printlevel); fprintf('\n  %s, %2.1e',msg,trXZnew); end
               termcode = -8;
               breakyes = 1;
            end
         else
            X = ops(X,'+',dX,pstep);
            y = y + dstep*dy;
            Z = ops(Z,'+',dZ,dstep);
         end
      end
%%---------------------------------------------------------------
%% adjust linear blk arising from unrestricted blk
%%---------------------------------------------------------------
      if (~breakyes)
         for p = 1:size(blk,1)
            if (u2lblk(p) == 1)
               len = blk{p,2}/2;
               xtmp = min(X{p}([1:len]),X{p}(len+[1:len]));
               alpha = 0.8;
               X{p}([1:len])     = X{p}([1:len]) - alpha*xtmp;
               X{p}(len+[1:len]) = X{p}(len+[1:len]) - alpha*xtmp;
               if (mu < 1e-4) %% old: (mu < 1e-7)
                  Z{p} = 0.5*mu./max(1,X{p}); %% good to keep this step
               else
                  ztmp = min(1,max(Z{p}([1:len]),Z{p}(len+[1:len])));
                  if (dual_infeas > 1e-4 & dstep < 0.2)
                     beta = 0.3;
                  else
                     beta = 0.0;
                  end
                  %% important to set beta = 0 at later stage.
                  Z{p}([1:len])     = Z{p}([1:len]) + beta*ztmp;
                  Z{p}(len+[1:len]) = Z{p}(len+[1:len]) + beta*ztmp;
               end
            end
         end
      end
%%--------------------------------------------------
%% perturb Z: do this step before checking for break
%%--------------------------------------------------
      if (~breakyes) & (~exist_analytic_term)
         trXZtmp = blktrace(blk,X,Z);
         trXE  = blktrace(blk,X,EE);
         Zpert = max(1e-12,0.2*min(relgap,prim_infeas)).*normC2./normE2;
         Zpert = min(Zpert,0.1*trXZtmp./trXE);
         Zpert = min([1,Zpert,1.5*Zpertold]);
         if (infeas < 0.1)
            Z = ops(Z,'+',EE,Zpert);
            [Zchol,indef(2)] = blkcholfun(blk,Z);
            if any(indef(2))
               msg = 'stop: Z not positive definite';
               if (printlevel); fprintf('\n  %s',msg); end
               termcode = -3;
               breakyes = 1;
            end
            %%if (printlevel > 2); fprintf(' %2.1e',Zpert); end
         end
         Zpertold = Zpert;
      end
%%---------------------------------------------------------------
%% compute rp, Rd, infeasibities, etc
%%---------------------------------------------------------------
%%
      AX  = AXfun(blk,At,par.permA,X);
      rp  = b-AX;
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
      ZpATynorm = ops(ZpATy,'norm');
      Rd  = ops(C,'-',ZpATy);
      objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
      obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
      gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
      relgap = gap/(1+sum(abs(obj)));
      prim_infeas = norm(rp)/normb2;
      dual_infeas = ops(Rd,'norm')/normC2;
      infeas = max(prim_infeas,dual_infeas);
      if (scale_data)
         infeas_org(1) = prim_infeas*normb;
         infeas_org(2) = dual_infeas*normC;
      end
      homRd = inf; homrp = inf;
      if (ops(parbarrier,'norm') == 0)
         if (obj(2) > 0); homRd = ZpATynorm/(obj(2)); end
         if (obj(1) < 0); homrp = norm(AX)/(-obj(1))/(normC); end
      end
      trXZ = blktrace(blk,X,Z,parbarrier);
      if (nn > 0); mu = trXZ/nn; else; mu = gap/ops(X,'getM'); end
      normX = ops(X,'norm');
%%
      runhist.pobj(iter+1)  = obj(1);
      runhist.dobj(iter+1)  = obj(2);
      runhist.gap(iter+1)   = gap;
      runhist.relgap(iter+1)  = relgap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas;
      runhist.step(iter+1)    = min(pstep,dstep);
      runhist.normX(iter+1)   = normX;
      runhist.cputime(iter+1) = etime(clock,tstart);
      timenew = clock;
      ttime.misc = ttime.misc + etime(timenew,timeold); timeold = timenew;
      [hh,mm,ss] = mytime(sum(runhist.cputime));
      if (printlevel>=3)
         fprintf('\n%2.0f|%4.3f|%4.3f',iter,pstep,dstep);
         fprintf('|%2.1e|%2.1e|%2.1e|',prim_infeas,dual_infeas,gap);
         fprintf('%- 7.6e| %s:%s:%s|',mean(obj),hh,mm,ss);
      end
%%--------------------------------------------------
%% check convergence
%%--------------------------------------------------
      param.use_LU      = use_LU;
      param.stoplevel   = stoplevel;
      param.termcode    = termcode;
      param.iter        = iter;
      param.obj         = obj;
      param.gap         = gap;
      param.relgap      = relgap;
      param.prim_infeas = prim_infeas;
      param.dual_infeas = dual_infeas;
      param.mu        = mu;
      param.homRd     = homRd;
      param.homrp     = homrp;
      param.AX        = AX;
      param.ZpATynorm = ZpATynorm;
      param.normX     = ops(X,'norm');
      param.normZ     = ops(Z,'norm');
      param.numpertdiagschur = numpertdiagschur;
      if (~breakyes)
         [param,breakyes,restart,msg2] = sqlpcheckconvg(param,runhist);
      end
      if (restart)
         [X,y,Z] = infeaspt(blk,At,C,b,2,1e5);
         rp  = b-AXfun(blk,At,par.permA,X);
         ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
         Rd  = ops(C,'-',ZpATy);
         trXZ = blktrace(blk,X,Z,parbarrier);
         mu   = trXZ/nn;
         gap  =  (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
         prim_infeas = norm(rp)/normb2;
         dual_infeas = ops(Rd,'norm')/normC2;
         infeas = max(prim_infeas,dual_infeas);
         [Xchol,indef(1)] = blkcholfun(blk,X);
         [Zchol,indef(2)] = blkcholfun(blk,Z);
         stoplevel = 3;
      end
%%--------------------------------------------------
%% check for break
%%--------------------------------------------------
      if ((prim_infeas < 1.5*prim_infeas_best) ...
         | (max(relgap,infeas) < 0.8*max(relgap_best,infeas_best))) ...
         & (max(relgap,dual_infeas) < 0.8*max(relgap_best,dual_infeas_best))
         Xbest = X; ybest = y; Zbest = Z;
         prim_infeas_best = prim_infeas;
         dual_infeas_best = dual_infeas;
         relgap_best = relgap; infeas_best = infeas;
         update_best(iter+1) = 1;
         %%fprintf('#')
      else
         update_best(iter+1) = 0;
      end
      if (max(relgap_best,infeas_best) < 1e-4 ...
          & norm(update_best(max(1,iter-1):iter+1)) == 0)
         msg = 'lack of progress in infeas';
         if (printlevel); fprintf('\n  %s',msg); end
         termcode = -9;
         breakyes = 1;
      end
      if (breakyes); break; end
   end
%%---------------------------------------------------------------
%% end of main loop
%%---------------------------------------------------------------
%%
   use_bestiter = 1;
   if (use_bestiter) & (param.termcode <= 0)
      X = Xbest; y = ybest; Z = Zbest;
      Xchol = blkcholfun(blk,X);
      Zchol = blkcholfun(blk,Z);
      AX = AXfun(blk,At,par.permA,X);
      rp = b-AX;
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
      Rd = ops(C,'-',ZpATy);
      objadd = blkbarrier(blk,X,Z,Xchol,Zchol,parbarrier) + objadd0;
      obj = (normb*normC)*[blktrace(blk,C,X), b'*y] + objadd;
      gap = (normb*normC)*blktrace(blk,X,Z) - diff(objadd);
      relgap = gap/(1+sum(abs(obj)));
      prim_infeas = norm(rp)/normb2;
      dual_infeas = ops(Rd,'norm')/normC2;
      infeas = max(prim_infeas,dual_infeas);
      runhist.pobj(iter+1)  = obj(1);
      runhist.dobj(iter+1)  = obj(2);
      runhist.gap(iter+1)   = gap;
      runhist.relgap(iter+1)  = relgap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas;
   end
%%---------------------------------------------------------------
%% unscale and produce infeasibility certificates if appropriate
%%---------------------------------------------------------------
   if (iter >= 1)
      [X,y,Z,termcode,resid,reldist,msg3] = ...
      sqlpmisc(blk,At,C,b,X,y,Z,par.permZ,param);
   end
%%---------------------------------------------------------------
%% recover unrestricted blk from linear blk
%%---------------------------------------------------------------
%%
   for p = 1:size(blk,1)
      if (u2lblk(p) == 1)
         n = blk{p,2}/2;
         X{p} = X{p}(1:n)-X{p}(n+[1:n]);
         Z{p} = Z{p}(1:n);
      end
   end
   for p = 1:size(ublkidx,1)
      if ~isempty(ublkidx{p,2})
         n0 = ublkidx{p,1}; idxB = setdiff([1:n0]',ublkidx{p,2});
         tmp = zeros(n0,1); tmp(idxB) = X{p}; X{p} = tmp;
         tmp = zeros(n0,1); tmp(idxB) = Z{p}; Z{p} = tmp;
      end
   end
   if (analytic_prob)
      X = X(1:numblkold); Z = Z(1:numblkold);
   end
%%---------------------------------------------------------------
%% print summary
%%---------------------------------------------------------------
%%
   maxC = 1+ops(ops(C,'abs'),'max');
   maxb = 1+max(abs(b));
   if (scale_data)
      dimacs = [infeas_org(1)*normb2/maxb; 0; infeas_org(2)*normC2/maxC; 0];
   else
      dimacs = [prim_infeas*normb2/maxb; 0; dual_infeas*normC2/maxC; 0];
   end
   dimacs = [dimacs; [-diff(obj); gap]/(1+sum(abs(obj)))];
   info.dimacs   = dimacs;
   info.termcode = termcode;
   info.iter     = iter;
   info.obj      = obj;
   info.gap      = gap;
   info.relgap   = relgap;
   info.pinfeas  = prim_infeas;
   info.dinfeas  = dual_infeas;
   info.cputime  = sum(runhist.cputime);
   info.time     = ttime;
   info.resid    = resid;
   info.reldist  = reldist;
   info.normX    = ops(X,'norm');
   info.normy    = norm(y);
   info.normZ    = ops(Z,'norm');
   info.normb    = normb2; info.maxb = maxb;
   info.normC    = normC2; info.maxC = maxC;
   info.normA    = normA2;
   info.msg1     = msg;
   info.msg2     = msg2;
   info.msg3     = msg3;
   sqlpsummary(info,ttime,infeas_org,printlevel);
   rand('state',randstate);
   randn('state',randnstate);
%%*****************************************************************************

repo: REST


text:
init_fig;

% c
axes('position', [.05, .57, .95, .37]);
hold on;
plot(y, 'color', col{8}/255);
alpha(.7);
plot(true_c, 'color', col{3}/255, 'linewidth', 1.5);
plot(c_oasis, '-.', 'color', col{5}/255);
if plot_cvx && exist('c_cvx', 'var')
    plot(c_cvx, '-.', 'color', col{7}/255);
end
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', []);
set(gca, 'ytick', 0:2);
ylabel('Fluor.');
box off;
if plot_cvx
    legend('Data', 'Truth', 'OASIS', 'CVX', 'location', 'northeast', 'orientation', 'horizental');
else
    legend('Data', 'Truth', 'OASIS', 'location', 'northeast', 'orientation', 'horizental');
end
% s
axes('position', [.05, .18, .95, .37]);
hold on;
plot(true_s, 'color', col{3}/255, 'linewidth', 1.5);
plot(s_oasis, '-.', 'color', col{5}/255);
if plot_cvx && exist('s_cvx', 'var')
    plot(s_cvx, '-.', 'color', col{7}/255);
end
axis tight;
xlim([0, 2000]);
set(gca, 'xtick', [0, 25, 50, 75]*30);
set(gca, 'xticklabel', get(gca, 'xtick')/30);
set(gca, 'ytick', [0,1]);
xlabel('Time [s]');
ylabel('Activity.');

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [psi] = ftH2O2psi(ftH2O)
% Convert pressure from feet of water column at 4 degrees to pounds per
% square inch
% Chad Greene 2012
psi = ftH2O*0.433528;

author: Sable

repo: mcbench-benchmarks


text:
function [y,fs,wmode,fidx]=v_readaif(filename,mode,nmax,nskip)
%V_READAIF  Read a .AIF format sound file [Y,FS,WMODE,FIDX]=(FILENAME,MODE,NMAX,NSKIP)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .AIF extension) or alternatively
%                 can be the FIDX output from a previous call to READAIF
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%   File I/O: 'f'    Do not close file on exit
%             'd'    Look in data directory: v_voicebox('dir_data')
%
%	NMAX     maximum number of samples to read (or -1 for unlimited [default])
%	NSKIP    number of samples to skip from start of file
%               (or -1 to continue from previous read when FIDX is given instead of FILENAME [default])
%
% Output Parameters:
%
%	Y        data matrix of dimension (samples,channels)
%	FS       sample frequency in Hz
%	WMODE    mode string needed to recreate the data file
%	FIDX     Information row vector containing the element listed below.
%
%           (1)  file id
%				(2)  current position in file
%           (3)  dataoff	byte offset in file to start of data
%           (4)  nsamp	number of samples
%           (5)  nchan	number of channels
%           (6)  nbyte	bytes per data value
%           (7)  bits	number of bits of precision
%           (8)  code	Data format: 1=PCM, 2=ADPCM, 6=A-law, 7=Mu-law (always 1 for AIFF)
%           (9)  fs	sample frequency
%           (10) offset
%           (11) block size
%
%   If no output parameters are specified, header information will be printed.
%
%   For stereo data, y(:,1) is the left channel and y(:,2) the right

%  future enhancements: (1) handle LIST and CAT files as well
%                       (2) deal with offset and blocksize properly
%                       (3) handle other chunk types sensibly

%	   Copyright (C) Mike Brookes 1998
%      Version: $Id: v_readaif.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<1 error('Usage: [y,fs,wmode,fidx]=READAIFF(filename,mode,nmax,nskip)'); end

if nargin<2 mode='p';
else mode = [mode(:).' 'p'];
end
k=find((mode>='p') & (mode<='s'));
mno=all(mode~='o');                      % scale to input limits not output limits
sc=mode(k(1));

info=zeros(1,11);
if ischar(filename)
    if any(mode=='d')
        filename=fullfile(v_voicebox('dir_data'),filename);
    end
    fid=fopen(filename,'rb','b');
    if fid == -1
        fn=[filename,'.aif'];
        fid=fopen(fn,'rb','b');
        if fid ~= -1 filename=fn; end
    end
    if fid == -1 error(sprintf('Can''t open %s for input',filename)); end
    info(1)=fid;
else
    info=filename;
    fid=info(1);
end

if ~info(3)
    fseek(fid,0,-1);						% read file header
    header=fread(fid,4,'uchar');
    if header' ~= 'FORM' fclose(fid); error(sprintf('File %s does not begin with a FORM group ID',filename)); end
    filen=fread(fid,1,'ulong')-4;
    header=fread(fid,4,'uchar');
    if header' ~= 'AIFF' fclose(fid); error(sprintf('File %s does not begin with a AIFF type ID',filename)); end

    fmt=0;
    data=0;
    while filen>=8						% loop to read all chunks
        header=fread(fid,4,'char');
        len=fread(fid,1,'ulong');
        lenx=len+rem(len,2);    % round up to an even number
        filen=filen-lenx-8;
        %fprintf(1,'%s chunk, %d bytes, %d bytes remaining\n',char(header'),len,filen);
        if header' == 'COMM'					% ******* found COMM chunk *********
            fmt=1;
            info(5)=fread(fid,1,'ushort');			%number of channels
            info(4)=fread(fid,1,'ulong');			%number of samples
            info(7)=fread(fid,1,'ushort');				% bits per sample
            info(6)=ceil(info(7)/8);            % bytes per sample
            info(8)=1;
            ieeex=fread(fid,5,'ushort');               % read a 10-byte extended ieee format number
            fsign=ieeex(1)>32767;
            ieeex(1)=ieeex(1)-32768*fsign;
            if ~ieeex info(9)=0;
            elseif ieeex(1)==32767 info(9)=nan;
            else
                info(9)=(1-2*fsign)*sum(pow2(ieeex(2:5)',ieeex(1)-(16398:16:16446)));
            end
            fseek(fid,lenx-18,0);				% skip to end of header
        elseif header' == 'SSND'				% ******* found DATA chunk *********
            info(10)=fread(fid,1,'ulong');			%number of channels
            info(11)=fread(fid,1,'ulong');			%number of channels
            info(3)=ftell(fid);
            data=1;
            fseek(fid,lenx-8,0);				% skip to end of chunk
        else							% ******* found unwanted chunk *********
            fseek(fid,lenx,0);
        end
    end
    if filen~=0
        fprintf(2,'READAIF warning: Inconsistent file length - %d extra bytes\n',filen);
    end
end
fs=info(9);
if ~fmt | ~data fclose(fid); error(sprintf('File %s does not contain COMM and SSND chunks',filename)); end

if nargin<4 nskip=info(2);
elseif nskip<0 nskip=info(2);
end

ksamples=info(4)-nskip;
if nargin>2
    if nmax>=0
        ksamples=min(nmax,ksamples);
    end
elseif ~nargout
    ksamples=min(5,ksamples);
end
if ksamples>0
    info(2)=nskip+ksamples;
    pk=pow2(0.5,8*info(6))*(1+(mno/2-all(mode~='n'))/pow2(0.5,info(7)));  % use modes o and n to determine effective peak
    fseek(fid,info(3)+info(6)*info(5)*nskip,-1);
    nsamples=info(5)*ksamples;
    if info(6)<3
        if info(6)<2
            y=fread(fid,nsamples,'schar');
        else
            y=fread(fid,nsamples,'short');
        end
    else
        if info(6)<4
            y=fread(fid,3*nsamples,'uchar');
            y=reshape(y,3,nsamples);
            y=[1 256 65536]*y-pow2(fix(pow2(y(3,:),-7)),24);
        else
            y=fread(fid,nsamples,'long');
        end
    end
    if sc ~= 'r'
        if sc=='s' sf=1/max(max(abs(y)),1);
        elseif sc=='p' sf=1/pk;
        else   sf=2.03033976/pk;
        end
        y=sf*y;
    else
        y=y*pow2(1,info(7)-8*info(6)); % shift to get the data into the LSB end
    end

    if info(5)>1 y = reshape(y,info(5),ksamples).'; end
else
    y=[];
end

if all(mode~='f') fclose(fid); end

if nargout>2
    if ~mno wmode=[wmode 'o']; end
    if any(mode=='n') wmode=[wmode 'n']; end
    wmode=[sc num2str(info(7))];
    fidx=info;
elseif ~nargout
    fprintf(1,'\n%d Hz sample rate\n%d channel x %d samples x %d bit = %.3g seconds\n',info([9 5 4]), info(7),info(4)/info(9));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function m = collatz_max ( seed )

%*****************************************************************************80
%
%% COLLATZ_MAX returns the largest element of a Collatz sequence.
%
%  Discussion:
%
%    The Collatz sequence is defined recursively as follows:
%
%      Let T be the current entry of the sequence, and U the next:
%
%      If T = 1, the sequence terminates (or U = 1, your choice);
%      Else if T is even, U = T / 2;
%      Else (if T is odd, and greater than 1) U = 3 * T + 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, the starting point for the sequence.
%    SEED must be positive.
%
%    Output, integer M, the maximum element of the sequence.
%
  m = seed;

  if ( seed <= 0 )
    return
  end

  t = seed;

  while ( t ~= 1 )

    if ( 2 * round ( t / 2 ) == t )
      t = t / 2;
    else
      t = 3 * t + 1;
    end

    m = max ( m, t );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function W = threshold_proportional(W, p)
%THRESHOLD_PROPORTIONAL     Proportional thresholding
%
%   W_thr = threshold_proportional(W, p);
%
%   This function "thresholds" the connectivity matrix by preserving a
%   proportion p (0<p<1) of the strongest weights. All other weights, and
%   all weights on the main diagonal (self-self connections) are set to 0.
%
%   Inputs: W,      weighted or binary connectivity matrix
%           p,      proportion of weights to preserve
%                       range:  p=1 (all weights preserved) to
%                               p=0 (no weights preserved)
%
%   Output: W_thr,  thresholded connectivity matrix
%
%
%   Mika Rubinov, U Cambridge,
%   Roan LaPlante, Martinos Center, MGH
%   Zitong Zhang, Penn Engineering

%   Modification history:
%   2010: Original (MR)
%   2012: Bug fix for symmetric matrices (RLP)
%   2015: Improved symmetricity test (ZZ)

n=size(W,1);                                %number of nodes
W(1:n+1:end)=0;                             %clear diagonal

if max(max(abs(W-W.'))) < 1e-10;            %if symmetric matrix
    W=triu(W);                              %ensure symmetry is preserved
    ud=2;                                   %halve number of removed links
else
    ud=1;
end

ind=find(W);                                %find all links
E=sortrows([ind W(ind)], -2);               %sort by magnitude
en=round((n^2-n)*p/ud);                     %number of links to be preserved

W(E(en+1:end,1))=0;                         %apply threshold

if ud==2                                    %if symmetric matrix
    W=W+W.';                                %reconstruct symmetry
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dep = cfg_vout_fileparts(job)

% Define virtual outputs for cfg_run_fileparts. The path names can either
% be assigned to a cfg_files input or to an evaluated cfg_entry. File names
% and extensions can only be assigned to an evaluated cfg_entry.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_vout_fileparts.m 4899 2012-09-05 13:44:17Z volkmar $

rev = '$Rev: 4899 $'; %#ok
dep(1)            = cfg_dep;
dep(1).sname      = sprintf('Directories');
dep(1).src_output = substruct('.','p');
dep(1).tgt_spec   = cfg_findspec({{'class','cfg_files', 'strtype','e'}});
dep(2)            = cfg_dep;
dep(2).sname      = sprintf('Directories (unique)');
dep(2).src_output = substruct('.','up');
dep(2).tgt_spec   = cfg_findspec({{'class','cfg_files', 'strtype','e'}});
dep(3)            = cfg_dep;
dep(3).sname      = sprintf('Filenames');
dep(3).src_output = substruct('.','n');
dep(3).tgt_spec   = cfg_findspec({{'strtype','e'}});
dep(4)            = cfg_dep;
dep(4).sname      = sprintf('Extensions');
dep(4).src_output = substruct('.','e');
dep(4).tgt_spec   = cfg_findspec({{'strtype','e'}});

author: spm


text:
%SETFEATDOM Reset feature domains of dataset
%
%   A = SETFEATDOM(A,FEATDOM,K)
%
% FEATDOM has to be a cell-array, defining for each feature
% the desired domain. The following formats are supported:
% []         : empty, all numbers are allowed
% R size(N,2): Defines N ranges on the real axis,
%              lowerbounds are R(:,1), upperbounds R(:,2).
%              N SHOULD BE AT LEAST 2. E.g:
%              R = [0 1; 0 1] defines a single region between 0 and 1.
%              R = [-1 1; 5 10; 11 11] defines two regions and
%              a single value (11).
% J size(1,N): Set of N integer values
% S size(N,F): String array of N strings. The feature values in the
%              DATA field of A will be coded by integers from 1 to N.
%              A feature like this may be assigned to the dataset
%              by A(:,j) = S, in which S is a string array, containing
%              a string for each object.
% By this routine and also by all changes in the dataset, TESTFEATDOM
% will be called to check consistency between data and feature domain
% descriptions.
%
% If given, K has to be a set of indices with length equal to
% length(FEATDOM) pointing to the feature domains to be reset.
% Checking is done as well. So reset data first if needed.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end


text:
function [lList,rList]=genAllLRBinTreesAlt(n)
%%GENALLLRBINTREESALT Generate representations of all binary trees having n
%                     nodes. l(0+1) gives the root node index. l(i+1) and
%                     r(i) then give the left and right descendents of node
%                     i. If there is a zero, then there is no descendent.
%                     The structure of the (left-right) relationships is
%                     what is unique, not the specific nodes. That is, the
%                     identity of the root node does not really matter. For
%                     example, node 1 will never have a left neighbor and
%                     node 4 will never have a right neighbor. However,
%                     independent of node labels, all binary trees are
%                     formed. The format of the output differs from the
%                     function genAllLRBinTrees.
%
%INPUTS: n The number of nodes that will be in the trees. n>=1.
%
%OUTPUTS: lList The (n+1)XnumTrees lists of left links, where lList(i+1,k)
%               is the left link of node i in tree k with lList(0+1,k)
%               being the root.
%         rList The nXnumTrees lists of right links.
%
%This function implements Algorithm L of Section 7.2.1.6 of [1]. The number
%of binary trees is CatalanNumber(n).
%
%EXAMPLE:
%To recreate the list of links in Table 3 of Section 7.2.1.6 of [1], one
%can use
% [lList,rList]=genAllLRBinTreesAlt(4)
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numBinTrees=CatalanNumber(n);

%Allocate space
lList=zeros(n+1,numBinTrees);
rList=zeros(n,numBinTrees);
kList=zeros(n,numBinTrees);

k=zeros(n,1);
o=zeros(n+1,1);
l=zeros(n+1,1);
r=zeros(n,1);

%Step L1
l((1:n)+1)=0;
r(1:(n-1))=2:n;
k(1:n)=0:(n-1);
o((1:n)+1)=-1;
l(0+1)=1;
o(0+1)=1;
r(n)=0;

for curTree=1:numBinTrees
    %Step L2
    lList(:,curTree)=l;
    rList(:,curTree)=r;
    kList(:,curTree)=k;
    j=n;
    p=0;

    while(1)
        %Step L3
        if(o(j+1)>0)
            m=l(j+1);

            if(m~=0)
               %Step L5
                if(j==0)
                    return;
                end
                l(j+1)=r(m);
                r(m)=j;
                k(j)=m;
                x=k(m);
                if(x==0)
                    l(p+1)=m;
                else
                    r(x)=m;
                end
                break;
            else
                %p=j;
                o(j+1)=-o(j+1);
                j=j-1;
                continue;
            end
        else%o(j+1)<0
            m=k(j);

            if(m~=0)
                %Step L4
                r(m)=l(j+1);
                l(j+1)=m;
                x=k(m);
                k(j)=x;
                if(x==0)
                   l(p+1)=j;
                else
                   r(x)=j;
                end
                break;
            else
                p=j;
                o(j+1)=-o(j+1);
                j=j-1;
                continue;
            end
        end
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [cpulse, verbose] = tapas_physio_get_oxy_pulses_filtered(c, t, ...
            dt120, verbose)
% Determines peaks of pulse oximeter data after Gaussian and high pass
% filtering, thresholding and assuming maximum heart rate (minimum peak
% distance)
%
%   [cpulse, verbose] = tapas_physio_get_oxy_pulses_filtered(c, t, ...
%            dt120, verbose);
%
% IN
%   c               [nSamples, 1] raw pulse oximeter samples
%   t               [nSamples, 1] time vector corresponding to samples (in seconds)
%   dt120           number of samples corresponding to a heart rate of 120 beats
%                   per minutes, i.e. number of samples acquiredi 0.5 seconds
%   verbose         Substructure of PhysIO, holding verbose.level and
%                   verbose.fig_handles with plotted figure handles
%                   debugging plots for thresholding are only provided, if verbose.level >=2
%
% OUT
%   cpulse          time points (seconds) of detected cardiac pulses
%   (heartbeat events, e.g. R-peaks)
%   verbose         Substructure of PhysIO, augmentedy by the additional
%                   figure handles created during this function
%
% EXAMPLE
%   tapas_physio_get_oxy_pulses_filtered
%
%   See also

% Author: Lars Kasper
% Created: 2014-08-03
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

dt = t(2) - t(1);
c = c-mean(c); c = c./max(c); % normalize time series

% smooth noisy pulse oximetry data to detect peaks
w = tapas_physio_gausswin(2*floor(dt120/2)+1, 1);  % Odd number of samples
sc = tapas_physio_conv(c, w, 'symmetric');
sc = sc-mean(sc); sc = sc./max(sc); % normalize time series

% Highpass filter to remove drifts
cutoff = 1/dt; %1 seconds/per sampling units
forder = 2;
[b,a] = butter(forder,2/cutoff, 'high');
sc =filter(b,a, sc);
sc = sc./max(sc);

[tmp, cpulse] = tapas_physio_findpeaks(sc, 'minpeakheight',...
    thresh_cardiac.min, 'minpeakdistance', dt120);

if verbose.level >=2 % visualise influence of smoothing on peak detection
    verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params();
    set(gcf, 'Name', 'Preproc: PPU-OXY: Tresholding Maxima for Heart Beat Detection');
    [tmp, cpulse2] = tapas_physio_findpeaks(c,'minpeakheight',thresh_cardiac.min,'minpeakdistance', dt120);
    plot(t, c, 'k');
    hold all;
    plot(t, sc, 'r', 'LineWidth', 2);
    hold all
    hold all;stem(t(cpulse2),c(cpulse2), 'k--');
    hold all;stem(t(cpulse),sc(cpulse), 'm', 'LineWidth', 2);
    plot(t, repmat(thresh_cardiac.min, length(t),1),'g-');
    legend('Raw PPU time series', 'Smoothed PPU Time Series', ...
        'Detected Heartbeats in Raw Time Series', ...
        'Detected Heartbeats in Smoothed Time Series', ...
        'Threshold (Min) for Heartbeat Detection');
end

cpulse = t(cpulse);

author: translationalneuromodeling

repo: tapas


text:
classdef TestDecomposeEssentialMat
    %TestDecomposeEssentialMat

    properties (Constant)
        fields = {'R1', 'R2', 't'};
    end

    methods (Static)
        function test_1
            E = [0 0 0; 0 0 -1; 0 1 0];
            S = cv.decomposeEssentialMat(E);
            validateattributes(S, {'struct'}, {'scalar'});
            assert(all(ismember(TestDecomposeEssentialMat.fields, fieldnames(S))));
            validateattributes(S.R1, {class(E)}, {'size',[3 3]});
            validateattributes(S.R2, {class(E)}, {'size',[3 3]});
            validateattributes(S.t, {class(E)}, {'vector', 'numel',3});
        end

        function test_error_argnum
            try
                cv.decomposeEssentialMat();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Releases

% SUN'2012
%    - 397 scene categories full image set
%    - 17000 annotated images

% SUN'2009
%    - 12000 annotated images
%
% 397-SUN

%addpath('/Users/torralba/atb/MatlabTools/LabelMeToolbox')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Building dataset of annotated images
HOMEANNOTATIONS = 'http://labelme.csail.mit.edu/Annotations';
HOMEIMAGES = 'http://labelme.csail.mit.edu/Images';
baseFolderDestination = '/csail/vision-torralba/datasets/';

% Build SUN index
disp('Reading folder list')
folder = folderlist(HOMEANNOTATIONS, 'users/antonio/static_sun_database');

disp('create index')
D = LMdatabase(HOMEANNOTATIONS, folder);

% Remove MT images
disp('remove Mechanical turk')
D = LMquery(D, 'object.polygon.username', '-mt_');

% Remove duplicate images
% get list of duplicate images
disp('remove duplicate files')
load dupFiles
dup = []; i = 0;
for n = 1:size(Sfiles,1)
    for m = 1:6
        f = Sfiles{n,m};
        if ~isempty(f)
            i = i + 1;
            k=max(strfind(f, '/'));
            f = f(k+1:end);
            dup{i} = f;
        end
    end
end
j_dup = [];
for n = 1:length(D)
    k = strmatch(D(n).annotation.filename, dup);
    if ~isempty(k)
        j_dup = [j_dup n];
    end
end

D = D(setdiff(1:length(D), j_dup));

% Remove empty images
%counts = LMcountobject(D);
%D = D(counts>0);

% Remove unreadable images
%readable = zeros(length(D),1);
%for i = 1:length(D)
%    disp(i)
%    try
%        img = LMimread(D, i, HOMEIMAGES);
%        D(i).annotation.imagesize.nrows = size(img,1);
%        D(i).annotation.imagesize.ncols = size(img,2);
%        readable(i) = 1;
%    catch
%        readable(i) = 0;
%    end
%    disp(sum(1-readable(1:i)))
%end

% Remove partially labeled images
disp('remove partially labeled images')
relativearea = LMlabeledarea(D);
good = find(relativearea>.8);
D = D(good);

% Add cropped tag
D = addcroplabel(D);

% Tags
%D = LMaddtags(D, 'tagsSUN.txt');

% Create list of images and paths
clear images folders
for i = 1:length(D)
    images{i} = D(i).annotation.filename;
    tmp = D(i).annotation.folder;
    tmp = strrep(tmp, 'users/antonio/static_sun_database/', '');
    folders{i} = tmp;
end

%destImages = '/Users/torralba/atb/Databases/SUN2012/Images';
%destAnnotations = '/Users/torralba/atb/Databases/SUN2012/Annotations/';
disp('Install SUN')
destImages = fullfile(baseFolderDestination, 'SUN2012/Images/');
destAnnotations = fullfile(baseFolderDestination, 'SUN2012/Annotations/');
LMinstall(folders, images, destImages, destAnnotations, 'http://labelme.csail.mit.edu/Images/users/antonio/static_sun_database/', 'http://labelme.csail.mit.edu/Annotations/users/antonio/static_sun_database/');

% make sure folder names are coherent
ADMINrenamefolder(destAnnotations);

Dtest = LMdatabase(destAnnotations, destImages);

% Make labelMe format:
% pascalfolder = '/Users/torralba/atb/Databases';
disp('create PASCAL release')
pascalfolder = baseFolderDestination;
databasename = 'SUN2012pascalformat';

labelme2pascal(Dtest, databasename, destImages, pascalfolder)

author: CSAILVision

repo: LabelMeToolbox


text:
function OUT = lt(varargin)

funname = 'lt';
narginchk(2,2)
OUT = builtincaller(funname,varargin{:});

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function h = setpapersize(s)

h = 1;

set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [0 0 s(1) s(2)]);

author: andrewssobral


text:
% OtQueryEmpty Check whether the record set is empty.
%
% Usage:
%    mt = OtQueryEmpty(hquery)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-13 Created (Johan Meijdam, VORtech)

function mt = OtQueryEmpty(hquery)
   mt = Ot('query.empty', hquery);
end


text:
%%*****************************************************************************
%% sqlp: solve an semidefinite-quadratic-linear program
%%       by infeasible path-following method.
%%
%%  [obj,X,y,Z,info,runhist] = sqlp(blk,At,C,b,OPTIONS,X0,y0,Z0);
%%
%%  Input: blk: a cell array describing the block diagonal structure of SQL data.
%%          At: a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%         b,C: data for the SQL instance.
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%     OPTIONS: a structure that specifies parameters required in sqlp.m,
%%              (if it is not given, the default in sqlparameters.m is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%----------------------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  scale_data ...
%%      (all have default values set in sqlparameters.m).
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

function [obj,X,y,Z,info,runhist] = sqlp(blk,At,C,b,OPTIONS,X0,y0,Z0)

if (nargin < 5); OPTIONS = []; end

isemptyAtb = 0;
if isempty(At) && isempty(b);
    %% Add redundant constraint: <-I,X> <= 0
    b = 0;
    At = ops(ops(blk,'identity'),'*',-1);
    numblk = size(blk,1);
    blk{numblk+1,1} = 'l'; blk{numblk+1,2} = 1;
    At{numblk+1,1} = 1; C{numblk+1,1} = 0;
    isemptyAtb = 1;
end
%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%

matlabversion = sscanf(version,'%f');
if strcmp(computer,'PCWIN64') || strcmp(computer,'GLNXA64')
    par.computer = 64;
else
    par.computer = 32;
end
par.matlabversion = matlabversion(1);
par.vers        = 0;
par.predcorr    = 1;
par.gam         = 0;
par.expon       = 1;
par.gaptol      = 1e-8;
par.inftol      = 1e-8;
par.steptol     = 1e-6;
par.maxit       = 100;
par.printlevel  = 3;
par.stoplevel   = 1;
par.scale_data  = 0;
par.spdensity   = 0.4;
par.rmdepconstr = 0;
par.smallblkdim = 50;
par.schurfun     = cell(size(blk,1),1);
par.schurfun_par = cell(size(blk,1),1);
%%
parbarrier  = cell(size(blk,1),1);
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'s') || strcmp(pblk{1},'q')
        parbarrier{p} = zeros(1,length(pblk{2}));
    elseif strcmp(pblk{1},'l') || strcmp(pblk{1},'u' )
        parbarrier{p} = zeros(1,sum(pblk{2}));
    end
end
parbarrier_0 = parbarrier;
%%
if nargin > 4,
    if isfield(OPTIONS,'vers');        par.vers     = OPTIONS.vers; end
    if isfield(OPTIONS,'predcorr');    par.predcorr = OPTIONS.predcorr; end
    if isfield(OPTIONS,'gam');         par.gam      = OPTIONS.gam; end
    if isfield(OPTIONS,'expon');       par.expon    = OPTIONS.expon; end
    if isfield(OPTIONS,'gaptol');      par.gaptol   = OPTIONS.gaptol; end
    if isfield(OPTIONS,'inftol');      par.inftol   = OPTIONS.inftol; end
    if isfield(OPTIONS,'steptol');     par.steptol  = OPTIONS.steptol; end
    if isfield(OPTIONS,'maxit');       par.maxit    = OPTIONS.maxit; end
    if isfield(OPTIONS,'printlevel');  par.printlevel  = OPTIONS.printlevel; end
    if isfield(OPTIONS,'stoplevel');   par.stoplevel   = OPTIONS.stoplevel; end
    if isfield(OPTIONS,'scale_data');  par.scale_data  = OPTIONS.scale_data; end
    if isfield(OPTIONS,'spdensity');   par.spdensity   = OPTIONS.spdensity; end
    if isfield(OPTIONS,'rmdepconstr'); par.rmdepconstr = OPTIONS.rmdepconstr; end
    if isfield(OPTIONS,'smallblkdim'); par.smallblkdim = OPTIONS.smallblkdim; end
    if isfield(OPTIONS,'parbarrier');
        parbarrier = OPTIONS.parbarrier;
        if isempty(parbarrier); parbarrier = parbarrier_0; end
        if ~iscell(parbarrier);
            tmp = parbarrier; clear parbarrier; parbarrier{1} = tmp;
        end
        if (length(parbarrier) < size(blk,1))
            len = length(parbarrier);
            parbarrier(len+1:size(blk,1)) = parbarrier_0(len+1:size(blk,1));
        end
    end
    if isfield(OPTIONS,'schurfun');
        par.schurfun = OPTIONS.schurfun;
        if ~isempty(par.schurfun); par.scale_data = 0; end
    end
    if isfield(OPTIONS,'schurfun_par'); par.schurfun_par = OPTIONS.schurfun_par; end
    if isempty(par.schurfun);     par.schurfun = cell(size(blk,1),1); end
    if isempty(par.schurfun_par); par.schurfun_par = cell(size(blk,1),1); end
end
if (size(blk,2) > 2); par.smallblkdim = 0; end
%%
%%-----------------------------------------
%% convert matrices to cell arrays.
%%-----------------------------------------
%%
if ~iscell(At); At = {At}; end;
if ~iscell(C);  C = {C}; end;
if all(size(At) == [size(blk,1), length(b)]);
    convertyes = zeros(size(blk,1),1);
    for p = 1:size(blk,1)
        if strcmp(blk{p,1},'s') && all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
        end
    end
    if any(convertyes)
        if (par.printlevel);
            fprintf('\n sqlp: converting At into required format');
        end
        At = svec(blk,At,ones(size(blk,1),1));
    end
end
%%
%%-----------------------------------------
%% validate SQLP data.
%%-----------------------------------------
%%
% tstart = cputime;
[blk,At,C,b,blkdim,numblk,parbarrier] = validate(blk,At,C,b,par,parbarrier);
[blk,At,C,b,iscmp] = convertcmpsdp(blk,At,C,b);
if (iscmp) && (par.printlevel>=2);
    fprintf('\n SQLP has complex data');
end
if (nargin <= 5) || (isempty(X0) || isempty(y0) || isempty(Z0));
    par.startpoint = 1;
    [X0,y0,Z0] = infeaspt(blk,At,C,b);
else
    par.startpoint = 2;
    if ~iscell(X0);  X0 = {X0}; end;
    if ~iscell(Z0);  Z0 = {Z0}; end;
    y0 = real(y0);
    if (length(y0) ~= length(b));
        error('sqlp: length of b and y0 not compatible');
    end
    [X0,Z0] = validate_startpoint(blk,X0,Z0,par.spdensity,iscmp);
end
if (par.printlevel>=2)
    fprintf('\n num. of constraints = %2.0d',length(b));
    if blkdim(1);
        fprintf('\n dim. of sdp    var  = %2.0d,',blkdim(1));
        fprintf('   num. of sdp  blk  = %2.0d',numblk(1));
    end
    if blkdim(2);
        fprintf('\n dim. of socp   var  = %2.0d,',blkdim(2));
        fprintf('   num. of socp blk  = %2.0d',numblk(2));
    end
    if blkdim(3); fprintf('\n dim. of linear var  = %2.0d',blkdim(3)); end
    if blkdim(4); fprintf('\n dim. of free   var  = %2.0d',blkdim(4)); end
end
%%
%%-----------------------------------------
%% detect unrestricted blocks in linear blocks
%%-----------------------------------------
%%
user_supplied_schurfun = 0;
for p = 1:size(blk,1)
    if ~isempty(par.schurfun{p}); user_supplied_schurfun = 1; end
end
if (user_supplied_schurfun == 0)
    [blk2,At2,C2,ublkinfo,parbarrier2,X02,Z02] = ...
        detect_ublk(blk,At,C,parbarrier,X0,Z0,par.printlevel);
else
    blk2 = blk; At2 = At; C2 = C;
    parbarrier2 = parbarrier; X02 = X0; Z02 = Z0;
    ublkinfo = cell(size(blk2,1),1);
end
ublksize = blkdim(4);
for p = 1:size(ublkinfo,1)
    ublksize = ublksize + length(ublkinfo{p});
end
%%
%%-----------------------------------------
%% detect diagonal blocks in semidefinite blocks
%%-----------------------------------------
%%
if (user_supplied_schurfun==0)
    [blk3,At3,C3,diagblkinfo,diagblkchange,parbarrier3,X03,Z03] = ...
        detect_lblk(blk2,At2,C2,b,parbarrier2,X02,Z02,par.printlevel);
else
    blk3 = blk2; At3 = At2; C3 = C2;
    parbarrier3 = parbarrier2; X03 = X02; Z03 = Z02;
    diagblkchange = 0;
    diagblkinfo = cell(size(blk3,1),1);
end
%%
%%-----------------------------------------
%% main solver
%%-----------------------------------------
%%
%    exist_analytic_term = 0;
%    for p = 1:size(blk3,1);
%       idx = find(parbarrier3{p} > 0);
%       if ~isempty(idx); exist_analytic_term = 1; end
%    end
%
if (par.vers == 0);
    if blkdim(1); par.vers = 1; else par.vers = 2; end
end
par.blkdim = blkdim;
par.ublksize = ublksize;
[obj,X3,y,Z3,info,runhist] = ...
    sqlpmain(blk3,At3,C3,b,par,parbarrier3,X03,y0,Z03);
%%
%%-----------------------------------------
%% recover semidefinite blocks from linear blocks
%%-----------------------------------------
%%
if any(diagblkchange)
    X2 = cell(size(blk2,1),1); Z2 = cell(size(blk2,1),1);
    count = 0;
    for p = 1:size(blk2,1)
        pblk = blk2(p,:);
        n = sum(pblk{2});
        blkno      = diagblkinfo{p,1};
        idxdiag    = diagblkinfo{p,2};
        idxnondiag = diagblkinfo{p,3};
        if ~isempty(idxdiag)
            len = length(idxdiag);
            Xtmp = [idxdiag,idxdiag,X3{end}(count+1:count+len); n, n, 0];
            Ztmp = [idxdiag,idxdiag,Z3{end}(count+1:count+len); n, n, 0];
            if ~isempty(idxnondiag)
                [ii,jj,vv] = find(X3{blkno});
                Xtmp = [Xtmp; idxnondiag(ii),idxnondiag(jj),vv]; %#ok
                [ii,jj,vv] = find(Z3{blkno});
                Ztmp = [Ztmp; idxnondiag(ii),idxnondiag(jj),vv]; %#ok
            end
            X2{p} = spconvert(Xtmp);
            Z2{p} = spconvert(Ztmp);
            count = count + len;
        else
            X2(p) = X3(blkno); Z2(p) = Z3(blkno);
        end
    end
else
    X2 = X3; Z2 = Z3;
end
%%
%%-----------------------------------------
%% recover linear block from unrestricted block
%%-----------------------------------------
%%
numblk = size(blk,1);
numblknew = numblk;
X = cell(numblk,1); Z = cell(numblk,1);
for p = 1:numblk
    n = blk{p,2};
    if isempty(ublkinfo{p,1})
        X{p} = X2{p}; Z{p} = Z2{p};
    else
        Xtmp = zeros(n,1); Ztmp = zeros(n,1);
        Xtmp(ublkinfo{p,1}) = max(0,X2{p});
        Xtmp(ublkinfo{p,2}) = max(0,-X2{p});
        Ztmp(ublkinfo{p,1}) = max(0,Z2{p});
        Ztmp(ublkinfo{p,2}) = max(0,-Z2{p});
        if ~isempty(ublkinfo{p,3})
            numblknew = numblknew + 1;
            Xtmp(ublkinfo{p,3}) = X2{numblknew};
            Ztmp(ublkinfo{p,3}) = Z2{numblknew};
        end
        X{p} = Xtmp; Z{p} = Ztmp;
    end
end
%%
%%-----------------------------------------
%% recover complex solution
%%-----------------------------------------
%%
if (iscmp)
    for p = 1:numblk
        pblk = blk(p,:);
        n = sum(pblk{2})/2;
        if strcmp(pblk{1},'s');
            X{p} = X{p}(1:n,1:n) + sqrt(-1)*X{p}(n+1:2*n,1:n);
            Z{p} = Z{p}(1:n,1:n) + sqrt(-1)*Z{p}(n+1:2*n,1:n);
            X{p} = 0.5*(X{p}+X{p}');
            Z{p} = 0.5*(Z{p}+Z{p}');
        end
    end
end
if (isemptyAtb)
    X = X(1:end-1); Z = Z(1:end-1);
end
%%*****************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [estSig, cost] = ILRMAISS(mixSig, nSrc, sampFreq, nBases, fftSize, shiftSize, windowType, nIter, ilrmaType, refMic, applyNormalize, applyWhitening, drawConv)
% Blind source separation with independent low-rank matrix analysis (ILRMA) based on iterative source steering update rule
%
% Coded by D. Kitamura (d-kitamura@ieee.org)
%
% Copyright 2021 Daichi Kitamura
%
% These programs are distributed only for academic research at
% universities and research institutions.
% It is not allowed to use or modify these programs for commercial or
% industrial purpose without our permission.
% When you use or modify these programs and write research articles,
% cite the following references:
%
% # Original paper (The algorithm was called "Rank-1 MNMF" in this paper)
% D. Kitamura, N. Ono, H. Sawada, H. Kameoka, H. Saruwatari, "Determined
% blind source separation unifying independent vector analysis and
% nonnegative matrix factorization," IEEE/ACM Trans. ASLP, vol. 24,
% no. 9, pp. 1626-1641, September 2016.
%
% # Book chapter (The algorithm was renamed as "ILRMA")
% D. Kitamura, N. Ono, H. Sawada, H. Kameoka, H. Saruwatari, "Determined
% blind source separation with independent low-rank matrix analysis,"
% Audio Source Separation. Signals and Communication Technology.,
% S. Makino, Ed. Springer, Cham, pp. 125-155, March 2018.
%
% # Original AuxIVA-ISS paper
% S. Robin and N. Ono,
% "Fast and stable blind source separation with rank-1 updates,"
% Proc. ICASSP, pp.236–240, 2020.
%
% See also:
% http://d-kitamura.net
% http://d-kitamura.net/demo-ILRMA_en.html
%
% [syntax]
%   [estSig,cost] = ILRMAISS(mixSig,nSrc,sampFreq,nBases,fftSize,shiftSize,windowType,nIter,ilrmaType,refMic,applyNormalize,applyWhitening,drawConv)
%
% [inputs]
%         mixSig: observed mixture (sigLen x nCh)
%           nSrc: number of sources in the mixture (scalar)
%       sampFreq: sampling frequency [Hz] of mixSig (scalar)
%         nBases: number of bases in NMF model (scalar, # of bases for "each" source when ilrmaType=1, and # of bases for "all" sources when ilrmaType=2, default: 4)
%        fftSize: window length [points] in STFT (scalar, default: next higher power of 2 that exceeds 0.256*sampFreq)
%      shiftSize: shift length [points] in STFT (scalar, default: fftSize/2)
%     windowType: window function used in STFT (name of window function, default: 'hamming')
%          nIter: number of iterations in the parameter update in ILRMA (scalar, default: 100)
%      ilrmaType: without or with partitioning function (1: ILRMA without partitioning function (ILRMA type1), 2: ILRMA with partitioning function (ILRMA type2), default: 1)
%         refMic: reference microphone for applying back projection (default: 1)
% applyNormalize: normalize parameters in each iteration to avoid numerical divergence (0: do not apply, 1: average-power-based normalization, 2: back projection (only for ILRMA type1), normalization may collapse monotonic decrease of the cost function, default: 0)
% applyWhitening: apply whitening to the observed multichannel spectrograms or not (true or false, default: true)
%       drawConv: plot cost function values in each iteration or not (true or false, default: false)
%
% [outputs]
%         estSig: estimated signals (sigLen x nCh x nSrc)
%           cost: convergence behavior of cost function in ILRMA (nIter+1 x 1)
%

% Arguments check and set default values
arguments
    mixSig (:,:) double
    nSrc (1,1) double {mustBeInteger(nSrc)}
    sampFreq (1,1) double
    nBases (1,1) double {mustBeInteger(nBases)} = 4
    fftSize (1,1) double {mustBeInteger(fftSize)} = 2^nextpow2(0.256*sampFreq)
    shiftSize (1,1) double {mustBeInteger(shiftSize)} = fftSize/2
    windowType char {mustBeMember(windowType,{'hamming','hann','rectangular','blackman','sine'})} = 'hamming'
    nIter (1,1) double {mustBeInteger(nIter)} = 100
    ilrmaType (1,1) double {mustBeInteger(ilrmaType)} = 1
    refMic (1,1) double {mustBeInteger(refMic)} = 1
    applyNormalize (1,1) double {mustBeInteger(applyNormalize)} = 0
    applyWhitening (1,1) logical = true
    drawConv (1,1) logical = false
end

% Error check
[sigLen, nCh] = size(mixSig); % sigLen: signal length, nCh: number of channels
if sigLen < nCh; error("The size of mixSig might be wrong.\n"); end
if nCh < nSrc || nSrc < 2; error("The number of channels must be equal to or grater than the number of sources in the mixture.\n"); end
if sampFreq <= 0; error("The sampling frequency (sampFreq) must be a positive value.\n"); end
if nBases < 1; error("The number of bases (nBases) must be a positive integer value.\n"); end
if fftSize < 1; error("The FFT length in STFT (fftSize) must be a positive integer value.\n"); end
if shiftSize < 1; error("The shift length in STFT (shiftSize) must be a positive integer value.\n"); end
if nIter < 1; error("The number of iterations (nIter) must be a positive integer value.\n"); end
if ilrmaType ~= 1 && ilrmaType ~= 2; error("The ILRMA type (ilrmaType) must be set to 1 or 2.\n"); end
if refMic < 1 || refMic > nCh; error("The reference microphone must be an integer between 1 and nCh.\n"); end
if applyNormalize ~= 0 && applyNormalize ~= 1 && applyNormalize ~= 2; error("The normalization type (applyNormalize) must be set to 0, 1, or 2.\n"); end
if applyNormalize == 2 && ilrmaType == 2; error("The back-projection-based normalization only supports ILRMA type 1.\n"); end

% Apply multichannel short-time Fourier transform (STFT)
[mixSpecgram, windowInStft] = STFT(mixSig, fftSize, shiftSize, windowType);

% Apply whitening (decorrelate X so that the correlation matrix becomes an identity matrix) based on principal component analysis
if applyWhitening
    inputMixSpecgram = whitening(mixSpecgram, nSrc); % apply whitening, where dimension is reduced from nCh to nSrc when nSrc < nCh
else
    inputMixSpecgram = mixSpecgram(:,:,1:nSrc); % when nSrc < nCh, only mixSpecgram(:,:,1:nSrc) is input to ILRMA so that the number of microphones equals to the number of sources (determined condition)
end

% Apply ILRMA
if ilrmaType == 1
    [estSpecgram, cost] = local_ILRMA1ISS(inputMixSpecgram, nIter, nBases, applyNormalize, drawConv, mixSpecgram(:,:,refMic));
else
    [estSpecgram, cost] = local_ILRMA2ISS(inputMixSpecgram, nIter, nBases, applyNormalize, drawConv);
end

% Apply back projection (fix the scale ambiguity using the reference microphone channel)
scaleFixedSepSpecgram = backProjection(estSpecgram, mixSpecgram(:,:,refMic)); % scale-fixed estimated signal

% Inverse STFT for each source
estSig = ISTFT(scaleFixedSepSpecgram, shiftSize, windowInStft, sigLen);
end

%% Local function for ILRMA type1 (without pertitioning function) based on ISS
function [Y, cost] = local_ILRMA1ISS(X, nIter, L, applyNormalize, drawConv, refMixSpecgram)
% [inputs]
%              X: observed multichannel spectrograms (I x J x M)
%          nIter: number of iterations of the parameter updates
%              L: number of bases in NMF model for each source
% applyNormalize: normalize parameters in each iteration to avoid numerical divergence (0: do not apply, 1: average-power-based normalization, 2: back projection, normalization may collapse monotonic decrease of the cost function)
%       drawConv: plot cost function values in each iteration or not (true or false)
% refMixSpecgram: observed reference spectrogram before apply whitening (I x J)
%
% [outputs]
%              Y: estimated spectrograms of sources (I x J x N)
%           cost: convergence behavior of cost function in ILRMA (nIter+1 x 1)
%
% [scalars]
%              I: number of frequency bins,
%              J: number of time frames
%              M: number of channels (microphones)
%              N: number of sources (equals to M)
%              L: number of bases in NMF model for each source
%
% [matrices]
%              X: observed multichannel spectrograms (I x J x M)
%             pX: permuted observed multichannel spectrograms (M x J x I)
%             pY: permuted separated multisource spectrograms (N x J x I)
%              W: frequency-wise demixing matrix (N x M x I)
%              Y: estimated multisource spectrograms (I x J x N)
%              P: estimated multisource power spectrograms (I x J x N)
%              T: sourcewise basis matrix in NMF (I x L x N)
%              V: sourcewise activation matrix in NMF (L x J x N)
%              R: sourcewise low-rank model spectrogram constructed by T and V (I x J x N)
%              E: identity matrix (N x N)
%              U: model-spectrogram-weighted sample covariance matrix of the mixture (M x M)
%

% Initialization
[I,J,M] = size(X); % I:frequency bins, J: time frames, M: channels
pX = permute(X, [3,2,1]); % permuted X whose dimensions are M x J x I
N = M; % N: number of sources, which equals to M in ILRMA
W = zeros(N,M,I); % frequency-wise demixing matrix
Y = zeros(I,J,N); % estimated spectrograms of sources (Y(i,:,n) =  W(n,:,i)*pX(:,:,i))
for i = 1:I
    W(:,:,i) = eye(N); % initial demixing matrices are set to identity matrices
    Y(i,:,:) = (W(:,:,i)*pX(:,:,i)).'; % initial estimated spectrograms
end
P = max(abs(Y).^2, eps); % power spectrogram of Y
T = max(rand( I, L, N ), eps); % sourcewise basis matrix in NMF
V = max(rand( L, J, N ), eps); % sourcewise activation matrix in NMF
R = zeros(I,J,N); % sourcewise low-rank model spectrogram constructed by T and V (R(:,:,n) = T(:,:,n)*V(:,:,n))
for n = 1:N
    R(:,:,n) = T(:,:,n)*V(:,:,n); % initial source model defined by T and V
end
v = zeros(N,1); % v vector for iterative source steering
cost = zeros(nIter+1, 1);

% Calculate initial cost function value
if drawConv
    cost(1,1) = local_calcCostFunction( P, R, W, I, J );
end

% Optimize parameters in ILRMA (W, T, and V)
fprintf('Iteration:    ');
for iIter = 1:nIter
    fprintf('\b\b\b\b%4d', iIter);

    %%%%% Update parameters %%%%%
    for n = 1:N
        %%%%% Update rule of T %%%%%
        T(:,:,n) = T(:,:,n) .* sqrt((P(:,:,n)./(R(:,:,n).^2))*V(:,:,n).' ./ ( (1./R(:,:,n))*V(:,:,n).' ));
        T(:,:,n) = max(T(:,:,n), eps);
        R(:,:,n) = T(:,:,n)*V(:,:,n);
        %%%%% Update rule of V %%%%%
        V(:,:,n) = V(:,:,n) .* sqrt(T(:,:,n).'*(P(:,:,n)./(R(:,:,n).^2)) ./ ( T(:,:,n).'*(1./R(:,:,n)) ));
        V(:,:,n) = max(V(:,:,n), eps);
        R(:,:,n) = T(:,:,n)*V(:,:,n);
        %%%%% Update rule of Y %%%%%
        YY = Y .* conj(Y(:,:,n)); % I x J x N, using implicit expansion (IxJxN .* IxJx1)
        for i = 1:I
            for nn = 1:N % calculate v vector BEGIN
                d = sum((1./(R(i,:,nn))).*real(YY(i,:,n)), 2) / J; % scalar
                if nn ~= n
                    u = sum((1./(R(i,:,nn))).*YY(i,:,nn), 2) / J; % scalar
                    v(nn,1) = u / d;
                else
                    v(nn,1) = 1 - 1/sqrt(d);
                end
            end % calculate v vector END
            Y(i,:,:) = Y(i,:,:) - permute(v.*Y(i,:,n), [3,2,1]); % update Y, usnig implicit expansion (Nx1 .* 1xJ = NxJ -> 1xJxN)
        end
    end
    P = max(abs(Y).^2, eps); % power spectrogram of Y

    %%%%% Normalization %%%%%
    if applyNormalize == 1 % average-power-based normalization
        lambda = sqrt(sum(sum(P,1),2)/(I*J)); % 1 x 1 x N
        Y = Y./lambda; % I x J x N (use implicit expansion)
        lambdaPow = lambda.^2; % 1 x 1 x N
        P = P./lambdaPow; % I x J x N (use implicit expansion)
        R = R./lambdaPow; % I x J x N (use implicit expansion)
        T = T./lambdaPow; % I x L x N (use implicit expansion)
    elseif applyNormalize == 2 % back projection
        lambda = local_backProjection(Y, refMixSpecgram, I, N); % N x 1 x I
        pLambda = permute(lambda, [3,2,1]); % I x 1 x N
        Y = Y.*pLambda; % I x J x N (use implicit expansion)
        lambdaPow = abs(pLambda).^2; % I x 1 x N
        P = P.*lambdaPow; % I x J x N (use implicit expansion)
        R = R.*lambdaPow; % I x J x N (use implicit expansion)
        T = T.*lambdaPow; % I x L x N (use implicit expansion)
    end

    %%%%% Calculate cost function value %%%%%
    if drawConv
        pY = permute(Y, [3,2,1]); % IxJxN -> NxJxI
        for i = 1:I
            W(:,:,i) = pY(:,:,i)*pX(:,:,i)' / (pX(:,:,i)*pX(:,:,i)'); % derived by "Y(:,:,i) = W(:,:,i)*X(:,:,i)"
        end
        cost(iIter+1,1) = local_calcCostFunction( P, R, W, I, J );
    end
end

% Draw convergence behavior
if drawConv
    figure; plot((0:nIter), cost);
    set(gca, 'FontName', 'Times', 'FontSize', 16);
    xlabel('Number of iterations', 'FontName', 'Arial', 'FontSize', 16);
    ylabel('Value of cost function', 'FontName', 'Arial', 'FontSize', 16);
end
fprintf(' ILRMA1-ISS done.\n');
end

%% Local function for ILRMA type2 (with pertitioning function) based on ISS
function [Y, cost] = local_ILRMA2ISS(X, nIter, K, applyNormalize, drawConv)
% [inputs]
%              X: observed multichannel spectrograms (I x J x M)
%          nIter: number of iterations of the parameter updates
%              K: number of bases in NMF model shared for all the source
% applyNormalize: normalize parameters in each iteration to avoid numerical divergence (0: do not apply, 1: average-power-based normalization, normalization may collapse monotonic decrease of the cost function)
%       drawConv: plot cost function values in each iteration or not (true or false)
%
% [outputs]
%              Y: estimated spectrograms of sources (I x J x N)
%           cost: convergence behavior of cost function in ILRMA (nIter+1 x 1)
%
% [scalars]
%              I: number of frequency bins,
%              J: number of time frames
%              M: number of channels (microphones)
%              N: number of sources (equals to M)
%              K: number of bases in NMF model shared for all the source
%
% [matrices]
%              X: observed multichannel spectrograms (I x J x M)
%             pX: permuted observed multichannel spectrograms (M x J x I)
%             pY: permuted separated multisource spectrograms (N x J x I)
%              W: frequency-wise demixing matrix (N x M x I)
%              Y: estimated multisource spectrograms (I x J x N)
%              P: estimated multisource power spectrograms (I x J x N)
%              T: basis matrix in NMF (I x K)
%              V: activation matrix in NMF (K x J)
%              Z: partitioning function matrix in NMF that clusters K bases into N sources (N x K)
%              R: sourcewise low-rank model spectrogram constructed by T and V (I x J x N)
%              E: identity matrix (N x N)
%              U: model-spectrogram-weighted sample covariance matrix of the mixture (M x M)
%

% Initialization
[I,J,M] = size(X); % I:frequency bins, J: time frames, M: channels
pX = permute(X, [3,2,1]); % permuted X whose dimensions are M x J x I
N = M; % N: number of sources, which equals to M in ILRMA
W = zeros(N,M,I); % frequency-wise demixing matrix
Y = zeros(I,J,N); % estimated spectrograms of sources (Y(i,:,n) =  W(n,:,i)*pX(:,:,i))
for i = 1:I
    W(:,:,i) = eye(N); % initial demixing matrices are set to identity matrices
    Y(i,:,:) = (W(:,:,i)*pX(:,:,i)).'; % initial estimated spectrograms are set to the observed mixture spectrograms
end
P = max(abs(Y).^2, eps); % power spectrogram of Y
T = max(rand( I, K ), eps); % basis matrix in NMF shared for all the sources
V = max(rand( K, J ), eps); % activation matrix in NMF shared for all the sources
Z = max(rand( N, K ), eps); % partitioning function matrix in NMF that clusters K bases into N sources
Z  = Z./sum(Z,1); % ensure sum_n z_{nk} = 1 (use implicit expansion)
tmpT = zeros(size(T)); % temporal variable used in the update rule of T
tmpV = zeros(size(V)); % temporal variable used in the update rule of V
tmpZ = zeros(size(Z)); % temporal variable used in the update rule of Z
R = zeros(I,J,N); % sourcewise low-rank model spectrogram constructed by T, V, and Z (R(:,:,n) = T(:,:,n)*V(:,:,n))
for n = 1:N
    R(:,:,n) = (Z(n,:).*T)*V; % initial source model defined by T, V, and Z (use implicit expansion)
end
v = zeros(N,1); % v vector for iterative source steering
cost = zeros(nIter+1, 1);

% Calculate initial cost function value
if drawConv
    cost(1,1) = local_calcCostFunction( P, R, W, I, J );
end

% Optimize parameters in ILRMA (W, T, and V)
fprintf('Iteration:    ');
for iIter = 1:nIter
    fprintf('\b\b\b\b%4d', iIter);

    %%%%% Update parameters %%%%%
    %%%%% Update rule of Z %%%%%
    for n = 1:N
        tmpZ(n,:) = (( sum((T.'*(P(:,:,n)./(R(:,:,n).^2))).*V, 2) )./( sum((T.'*(1./R(:,:,n))).*V, 2) )).';
    end
    Z = Z .* sqrt(tmpZ);
    Z  = max(Z./sum(Z,1), eps); % ensure sum_n z_{nk} = 1 (use implicit expansion)
    for n = 1:N
        R(:,:,n) = (Z(n,:).*T)*V; % initial source model defined by NMF (use implicit expansion)
    end
    %%%%% Update rule of T %%%%%
    for i = 1:I
        Pi = squeeze(P(i,:,:)); % J x N
        Ri = squeeze(R(i,:,:)); % J x N
        tmpT(i,:) = (( sum((V*(Pi./(Ri.^2))).*(Z.'), 2) )./( sum((V*(1./Ri)).*(Z.'), 2) )).';
    end
    T = max(T.*sqrt(tmpT), eps);
    for n = 1:N
        R(:,:,n) = (Z(n,:).*T)*V; % initial source model defined by NMF (use implicit expansion)
    end
    %%%%% Update rule of V %%%%%
    for j = 1:J
        Pj = squeeze(P(:,j,:)); % I x N
        Rj = squeeze(R(:,j,:)); % I x N
        tmpV(:,j) = ( sum((T.'*(Pj./(Rj.^2))).*(Z.'), 2) )./( sum((T.'*(1./Rj)).*(Z.'), 2) );
    end
    V = max(V.*sqrt(tmpV), eps);
    for n = 1:N
        R(:,:,n) = (Z(n,:).*T)*V; % initial source model defined by NMF (use implicit expansion)
    end
    %%%%% Update rule of Y %%%%%
    for n=1:N
        YY = Y .* conj(Y(:,:,n)); % I x J x N, using implicit expansion (IxJxN .* IxJx1)
        for i = 1:I
            for nn = 1:N % calculate v vector BEGIN
                d = sum((1./(R(i,:,nn))).*real(YY(i,:,n)), 2) / J; % scalar
                if nn ~= n
                    u = sum((1./(R(i,:,nn))).*YY(i,:,nn), 2) / J; % scalar
                    v(nn,1) = u / d;
                else
                    v(nn,1) = 1 - 1/sqrt(d);
                end
            end % calculate v vector END
            Y(i,:,:) = Y(i,:,:) - permute(v.*Y(i,:,n), [3,2,1]); % update Y, usnig implicit expansion (Nx1 .* 1xJ = NxJ -> 1xJxN)
        end
    end
    P = max(abs(Y).^2, eps); % power spectrogram of Y

    %%%%% Normalization %%%%%
    if applyNormalize == 1
        lambda = sqrt(sum(sum(P,1),2)/(I*J)); % 1 x 1 x N
        Y = Y./lambda; % I x J x N (use implicit expansion)
        P = P./lambda.^2; % I x J x N (use implicit expansion)
        R = R./lambda.^2; % I x J x N (use implicit expansion)
        Zlambda = Z./(squeeze(lambda).^2); % N x K
        ZlambdaSum = sum(Zlambda,1); % 1 x K
        T = T.*ZlambdaSum; % I x K (use implicit expansion)
        Z = Zlambda./ZlambdaSum; % N x K (use implicit expansion)
    end

    %%%%% Calculate cost function value %%%%%
    if drawConv
        pY = permute(Y, [3,2,1]); % IxJxN -> NxJxI
        for i = 1:I
            W(:,:,i) = pY(:,:,i)*pX(:,:,i)' / (pX(:,:,i)*pX(:,:,i)'); % derived by "Y(:,:,i) = W(:,:,i)*X(:,:,i)"
        end
        cost(iIter+1,1) = local_calcCostFunction( P, R, W, I, J );
    end
end

% Draw convergence behavior
if drawConv
    figure; plot((0:nIter), cost);
    set(gca, 'FontName', 'Times', 'FontSize', 16);
    xlabel('Number of iterations', 'FontName', 'Arial', 'FontSize', 16);
    ylabel('Value of cost function', 'FontName', 'Arial', 'FontSize', 16);
end
fprintf(' ILRMA2-ISS done.\n');
end

%% Local function for calculating cost function value in ILRMA
function [ cost ] = local_calcCostFunction(P, R, W, I, J)
logDetAbsW = zeros(I,1);
for i = 1:I
    logDetAbsW(i,1) = log(max(abs(det(W(:,:,i))), eps));
end
cost = sum(sum(sum(P./R+log(R),3),2),1) - 2*J*sum(logDetAbsW, 1);
end

%% Local function for applying back projection and returns frequency-wise coefficients
function [ D ] = local_backProjection(Y, X, I, N)
D = zeros(I,N);
for i = 1:I
    Yi = squeeze(Y(i,:,:)).'; % N x J
    D(i,:) = X(i,:,1)*Yi'/(Yi*Yi'); % 1 x N
end
D(isnan(D) | isinf(D)) = 0; % replace NaN and Inf to 0
D = permute(D, [2,3,1]); % N x 1 x I
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: d-kitamura

repo: ILRMA

save_path: github-repos/MATLAB/d-kitamura-ILRMA


text:
%%generate data config
param.split_data_num=1;
param.file_path_cam1=fullfile(fileparts(pwd),'prid_2011','multi_shot','cam_a');
param.file_path_cam2=fullfile(fileparts(pwd),'prid_2011','multi_shot','cam_b');
param.save_traindata_filename='train_data_Prid';
param.save_testdata_filename='test_data_Prid';
%% generate ten split train/test data
for split_num_start=1:param.split_data_num
    if ~exist(strcat(param.save_traindata_filename,num2str(split_num_start)),'file')
        mkdir(strcat(param.save_traindata_filename,num2str(split_num_start)));
    end
    if ~exist(strcat(param.save_testdata_filename,num2str(split_num_start)),'file')
        mkdir(strcat(param.save_testdata_filename,num2str(split_num_start)));
    end
    subdir_cam1=dir(param.file_path_cam1);
    subdir_cam2=dir(param.file_path_cam2);
    %% select images >27 frame
    person_file_cam1=[];
    person_file_cam2=[];
    for i=3:202
        image_path_a=fullfile(param.file_path_cam1,subdir_cam1(i).name,'/');
        image_list_a=dir(image_path_a);
        len=length(image_list_a);
        if len<29
            continue;
        else
            person_file_cam1=[person_file_cam1;i-2];
        end
    end
    for i=3:202
        image_path_b=fullfile(param.file_path_cam2,subdir_cam2(i).name,'/');
        image_list_b=dir(image_path_b);
        len=length(image_list_b);
        if len<29
            continue;
        else
            person_file_cam2=[person_file_cam2;i-2];
        end
    end
    %%
    person_num=intersect(person_file_cam1,person_file_cam2);
    select_list_len=size(person_num,1);
    train_data_select=randperm(select_list_len,select_list_len/2);
    train_data_index=person_num(train_data_select);
    test_data_index=setdiff(person_num,train_data_index);
    %generate train data
    train_data_cam1=[];
    label_train_cam1=[];
    train_image_name_cam1={};
    train_data_cam2=[];
    label_train_cam2=[];
    train_image_name_cam2={};
    for i=1:(length(train_data_index))
        fprintf('process train data:%d/%d\n',i,(length(train_data_index)));
        index=train_data_index(i);
        image_path_cam1=fullfile(param.file_path_cam1,subdir_cam1(index+2).name,'/');
        image_list_cam1=dir(image_path_cam1);
        for j=3:length(image_list_cam1);
            image_name=strcat(image_path_cam1,image_list_cam1(j).name);
            image_data=imread(image_name);
            train_data_cam1=[train_data_cam1;reshape(image_data,1,64*128*3)];
            label_train_cam1=[label_train_cam1;i];
            train_image_name_cam1=[train_image_name_cam1,image_name];
        end
        image_path_cam2=fullfile(param.file_path_cam2,subdir_cam2(index+2).name,'/');
        image_list_cam2=dir(image_path_cam2);
        for j=3:length(image_list_cam2);
            image_name=strcat(image_path_cam2,image_list_cam2(j).name);
            image_data=imread(image_name);
            train_data_cam2=[train_data_cam2;reshape(image_data,1,64*128*3)];
            label_train_cam2=[label_train_cam2;i];
            train_image_name_cam2=[train_image_name_cam2,image_name];
        end
    end
    save(strcat(param.save_traindata_filename,num2str(split_num_start),'/train_data.mat'),'train_data_cam1','train_data_cam2','label_train_cam1','label_train_cam2','train_image_name_cam1','train_image_name_cam2');
    %% generate test data
    test_data_cam1=[];
    label_test_cam1=[];
    test_image_name_cam1={};
    test_data_cam2=[];
    label_test_cam2=[];
    test_image_name_cam2={};
    for i=1:(length(test_data_index))
        fprintf('process test data:%d/%d\n',i,(length(test_data_index)));
        index=test_data_index(i);
        image_path_cam1=fullfile(param.file_path_cam1,subdir_cam1(index+2).name,'/');
        image_list_cam1=dir(image_path_cam1);
        for j=3:length(image_list_cam1);
            image_name=strcat(image_path_cam1,image_list_cam1(j).name);
            image_data=imread(image_name);
            test_data_cam1=[test_data_cam1;reshape(image_data,1,64*128*3)];
            label_test_cam1=[label_test_cam1;i+89];
            test_image_name_cam1=[test_image_name_cam1,image_name];
        end
        image_path_cam2=fullfile(param.file_path_cam2,subdir_cam2(index+2).name,'/');
        image_list_cam2=dir(image_path_cam2);
        for j=3:length(image_list_cam2);
            image_name=strcat(image_path_cam2,image_list_cam2(j).name);
            image_data=imread(image_name);
            test_data_cam2=[test_data_cam2;reshape(image_data,1,64*128*3)];
            label_test_cam2=[label_test_cam2;i+89];
            test_image_name_cam2=[test_image_name_cam2,image_name];
        end
    end
    save(strcat(param.save_testdata_filename,num2str(split_num_start),'/test_data.mat'),'test_data_cam1','test_data_cam2','label_test_cam1','label_test_cam2','test_image_name_cam1','test_image_name_cam2');
end

author: liuyuisanai

repo: Quality-Aware-Network


text:
function s=get_name(a)

s=[get_name(a.algorithm)];  %% print name of algorithm

if iscell(a.child)
  s=[s '(' get_name(a.child{1}) ',...)'];
else
  s=[s '(' get_name(a.child) ')'];
end

eval_name                   %% print values of hyperparameters

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [varargout]=tubeplot(x,y,z,varargin)

% TUBEPLOT - plots a tube r along the space curve x,y,z.
%
% tubeplot(x,y,z) plots the basic tube with radius 1
% tubeplot(x,y,z,r) plots the basic tube with variable radius r (either a vector or a value)
% tubeplot(x,y,z,r,v) plots the basic tube with coloring dependent on the values in the vector v
% tubeplot(x,y,z,r,v,s) plots the tube with s tangential subdivisions
% (default is 6)
%
% [X,Y,Z]=tubeplot(x,y,z) returns [Nx3] matrices suitable for mesh or surf
%
% Note that the tube may pinch at points where the normal and binormal
% misbehaves. It is suitable for general space curves, not ones that
% contain straight sections. Normally the tube is calculated using the
% Frenet frame, making the tube minimally twisted except at inflexion points.
%
% To deal with this problem there is an alternative frame:
% tubeplot(x,y,z,r,v,s,vec) calculates the tube by setting the normal to
% the cross product of the tangent and the vector vec. If it is chosen so
% that it is always far from the tangent vector the frame will not twist unduly
%
% Example:
%
%  t=0:(2*pi/100):(2*pi);
%  x=cos(t*2).*(2+sin(t*3)*.3);
%  y=sin(t*2).*(2+sin(t*3)*.3);
%  z=cos(t*3)*.3;
%  tubeplot(x,y,z,0.14*sin(t*5)+.29,t,10)
%
% Written by Anders Sandberg, asa@nada.kth.se, 2005

  subdivs = 6;

  N=size(x,1);
  if (N==1)
    x=x';
    y=y';
    z=z';
    N=size(x,1);
  end

  if (nargin == 3)
    r=x*0+1;
  else
    r=varargin{1};
    if (size(r,1)==1 & size(r,2)==1)
      r=r*ones(N,1);
    end
  end
  if (nargin > 5)
    subdivs=varargin{3}+1;
  end
  if (nargin > 6)
    vec=varargin{4};
    [t,n,b]=frame(x,y,z,vec);
  else
    [t,n,b]=frenet(x,y,z);
  end

  X=zeros(N,subdivs);
  Y=zeros(N,subdivs);
  Z=zeros(N,subdivs);

  theta=0:(2*pi/(subdivs-1)):(2*pi);

  for i=1:N
    X(i,:)=x(i) + r(i)*(n(i,1)*cos(theta) + b(i,1)*sin(theta));
    Y(i,:)=y(i) + r(i)*(n(i,2)*cos(theta) + b(i,2)*sin(theta));
    Z(i,:)=z(i) + r(i)*(n(i,3)*cos(theta) + b(i,3)*sin(theta));
  end

  if (nargout==0)
    if (nargin > 4)
      V=varargin{2};
      if (size(V,1)==1)
	V=V';
      end
      V=V*ones(1,subdivs);
      surf(X,Y,Z,V);
    else
      surf(X,Y,Z);
    end
  else
    varargout(1) = {X};
    varargout(2) = {Y};
    varargout(3) = {Z};
  end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft