function [u, cost, useq, status, iAout, slack] = mpcblock_optimizer(...
    rseq, vseq, umin, umax, ymin, ymax, switch_in, x, old_u, iA, ...
    isQP, nmv, ny, degrees, Hinv, Kx, Ku1, Kut, Kr, Kv, Mlim, ...
    Mx, Mu1, Mv, utarget, p, uoff, voff, yoff, ...
    CustomSolver, CustomSolverCodeGen, UseSuboptimalSolution, ...
    QPAlgorithm, ASOptions, IPOptions, ADOptions, MIQPOptions, nxQP, openloopflag, ...
    no_umin, no_umax, no_ymin, no_ymax, no_cc, switch_inport, no_switch, ...
    enable_value, return_cost, H, return_sequence, Linv, Ac, ...
    ywt, uwt, duwt, rhoeps, no_ywt, no_uwt, no_duwt, no_rhoeps,...
    Wy, Wdu, Jm, SuJm, Su1, Sx, Hv, Wu, I1, ...
    isAdaptive, isLTV, A, Bu, Bv, C, Dv, Mrows, nCC, Ecc, Fcc, Scc, Gcc,...
    RYscale, RMVscale, m, isHyb, Mdis, Ndis, Vdis, numdis, maxdis)
% This function is unsupported and might change or be removed without
% notice in a future version.
%
%   MPCBLOCK_OPTIMIZER implements the fixed prediction and control horizon
%   feature used by "mpcmove"/"mpcmoveAdaptive"/"mpcmoveCodeGeneration"
%   commands as well as MPC block and Adaptive MPC block.  It only supports
%   "fixed-size data" for code generation.
%
%   Input definition:
%
%   From online signals
%       x: current controller state
%       old_u: last mv
%       rseq: reference preview sequence
%       vseq: md preview sequence
%       utarget: MV target preview sequence
%       umin: mv lower bound across prediction horizon
%       umax: mv upper bound across prediction horizon
%       ymin: ov lower bound across prediction horizon
%       ymax: ov upper bound across prediction horizon
%       ywt: OV weights
%       uwt: MV weights
%       duwt: MV Rate weights
%       rhoeps: ECR weight
%       A: prediction model A
%       Bu: prediction model B related to MV
%       Bv: prediction model B related to MD plus DXOffset
%       C: prediction model C
%       Dv: prediction model D related to MD plus DXoffset
%   
%   From constant parameters from @mpc
%       p: prediction horizon
%       degrees: # of decision variables
%       isQP: has constraint => no analytical solution ==> slack exists
%       isLTV: whether LTV plant model
%       H: Hessian
%       Wy: original OV weights
%       Wdu: original MV Rate weights
%       Wu: original MV weights
%       Jm: used to update H and Ac and calculate sequence
%       SuJm: used to update H
%       Su1: used to update H and calculate cost
%       Sx: used to update H and calculate cost
%       Hv: used to update H and calculate cost
%       I1: used to update H and calculate cost
%       Hinv: inverse Hessian
%       Linv: LT inverse Hessian
%       Kx: used to compute f 
%       Ku1: used to compute f
%       Kut: used to compute f
%       Kr: used to compute f
%       Kv: used to compute f
%       Ac: constraint coefficient matrix
%       Mlim: used to compute Bc
%       Mx: used to compute Bc
%       Mu1: used to compute Bc
%       Mv: used to compute Bc
%       Mrows: total contraints index
%       
%       nCC, Ecc, Fcc, Scc, Gcc are related to mixed constraints
%
%   Other utility constants
%       nmv: number of MV
%       ny: number of OV
%       nxQP: number of states in the plant and disturbance model
%       uoff: MV offset
%       voff: MD offset
%       yoff: OV offset
%       no_umin, no_umax, no_ymin, no_ymax, no_cc: whether online constraints
%       no_ywt, no_uwt, no_duwt, no_rhoeps: whether online weights
%       CustomSolver: whether use custom solver for simulation
%       CustomSolverCodeGen: whether use custom solver for code generation
%       UseSuboptimalSolution: whether use suboptimal solution
%       return_cost, return_sequence: return information
%
%       switch_in, switch_inport, no_switch, enable_value and openloopflag
%       are used to determine whether optimization is disabled.
%
%       RYscale, RMVscale: scaling factors for output and input constraints
%       m: free control moves (required for time-varying input constraints)

% Authors:  A. Bemporad, N.L. Ricker, R. Chen
% Copyright 1994-2018 The MathWorks, Inc.

%#codegen
ONE = ones(1,1,'like',rseq);
ZERO = zeros(1,1,'like',rseq);
ONE32 = ones('int32');
ZERO32 = zeros('int32');
SMALLHESSIAN = ONE*10*sqrt(eps);

% Rules regarding signal scaling;
% 1. Input signals "rseq"/"vseq"/"u_old" are scaled and offset-free
% 2. Input signals "x" are offset-free
% 3. Input signals "umin"/"umax"/"ymin"/"ymax" are in original units
% 4. Input signals of weights are dimensionless
% 5. Output signals (u, useq) are returned dimensionless.
% 6. Output "useq" is (p+1)-by-nmv with last row duplicated

% Rules regarding MV and USEQ values
% 1. In normal operations, return old_u + zopt
% 2. When optimization is disabled (either manually from block dialog or by the MMPC block), return old_u
% 3. When open-loop flag is on (no MPC object, design from scratch, findop or linearize), return old_u

%% Initialize default outputs
u=zeros(nmv,1,'like',rseq);
cost=ZERO;
useq=zeros(p+1,nmv,'like',rseq);
status=ONE;
iAout=false(size(iA));
slack=ZERO; 

TotalFreeMoves=degrees-int32(isQP);
if isHyb
    Mdisinv=linsolve(Mdis,eye(TotalFreeMoves),struct('LT',true)); % Mdis is unit lower triangular
else
    Mdisinv=ZERO;
end

%% if open loop, return u = old_u + uoff
optimization_disabled = switch_inport && ~no_switch && switch_in~=enable_value;
if openloopflag || optimization_disabled
    u=old_u+uoff;
    for i=1:nmv
        useq(:,i)=u(i);
    end
    return
end

%% initialize xQP
xQP=x(1:nxQP);

%% Impact of prediction model update (adaptive mode)
pny = p*ny;
pnu = p*nmv;
nCon = ONE32*size(Ac,1);   % Total number of QP constraints
if isAdaptive || isLTV
    % Update model coefficients used to define constraints
    if isAdaptive
        [SuJm, Sx, Su1, Hv] = mpc_constraintcoef(A(1:nxQP,1:nxQP,1), Bu(1:nxQP,:,1), Bv(1:nxQP,:,1), C(:,1:nxQP,1), Dv(:,:,1), Jm);
    else
        [SuJm, Sx, Su1, Hv] = mpc_constraintcoefLTV(A(1:nxQP,1:nxQP,:), Bu(1:nxQP,:,:), Bv(1:nxQP,:,:), C(:,1:nxQP,:), Dv, Jm);
    end
    % Update matrices affected directly by Su update.
    % Make sure Mrows has been defined.  It might have been [] or zero.
    if ~isempty(Mrows) && (Mrows(1)>ZERO)     
        for i = 1:length(Mrows)
            % Mrows indicates the indices of the general bounds that were
            % retained in the QP.  Use Mrows to update these constraints
            % one-by-one.
            ix = Mrows(i);
            % Only first 2*pny constraint rows depend on model coeffients.  
            if ix <= pny    % ymax bounds
                % Special handling of Ac updates because last column is for the
                % slack variable and its coefficients don't change.
                Ac(i,1:end-1) = -SuJm(ix,:);
                if isHyb
                    Ac(i,1:end-1)=Ac(i,1:end-1)*Mdisinv;
                end
                Mx(i,:)       = -  Sx(ix,:);
                Mu1(i,:)      = - Su1(ix,:);
                Mv(i,:)       = -  Hv(ix,:);
            elseif ix <= 2*pny      % ymin bounds
                ix = ix - pny;
                Ac(i,1:end-1)  = SuJm(ix,:);
                if isHyb
                    Ac(i,1:end-1)=Ac(i,1:end-1)*Mdisinv;
                end
                Mx(i,:)        =   Sx(ix,:);
                Mu1(i,:)       =  Su1(ix,:);
                Mv(i,:)        =   Hv(ix,:);
            else
                break   % All subsequent Mrow values must be > 2*pny
            end
        end   
    end
end

% update existing custom constraint when (1) adaptive/ltv mode or (2) online change
if nCC>ZERO
    if (isAdaptive || isLTV) || ~no_cc
        [MuCC, MvCC, Mu1CC, MxCC, MlimCC] = mpc_customconstraintcoef(SuJm, Sx, Su1, Hv, C(:,1:nxQP,:), Dv, Jm, nCC, Ecc, Fcc, Scc, Gcc, uoff, voff, yoff);
        % Update custom (mixed) constraint coefficients.  In general, 
        % all such constraints depend on the prediction model through the 
        % F*y term.
        nCCrows = nCC*(p+1);
        rows = (nCon-nCCrows+1):nCon;
        Ac(rows,1:end-1) = -MuCC;  
        if isHyb
            Ac(rows,1:end-1)=Ac(rows,1:end-1)*Mdisinv;
        end
        Mx(rows,:)       =  MxCC;
        Mu1(rows,:)      =  Mu1CC;
        Mv(rows,:)       =  MvCC;
        Mlim(rows,:)     =  MlimCC;
    end
end

%% Online weight adjustments 
if ~no_ywt || ~no_duwt || ~no_uwt || isAdaptive || isLTV
    if ~no_ywt
        % Online uniform or time-varying OV weights: 
        % ywt is a column vector of ny*(1..p) elements, the first ny
        % elements for time k+1, next ny elements for time k+2, ...        
        % When offline uniform, Wy is a row vector of ny elements
        % When offline time-varying, Wy is a row vector p*ny elements
        Wy = updateWeights(Wy, ywt, ny, p, ZERO);
    end
    if ~no_uwt
        % Online uniform or time-varying MV weights: 
        % uwt is a column vector of nmv*(1..p) elements, the first nmv
        % elements for time k+1, next nmv elements for time k+2, ...        
        % When offline uniform, Wu is a row vector of nmv elements
        % When offline time-varying, Wu is a row vector p*nmv elements
        Wu = updateWeights(Wu, uwt, nmv, p, ZERO);
    end
    if ~no_duwt
        % Online uniform or time-varying MVRate weights: 
        % duwt is a column vector of nmv*(1..p) elements, the first nmv
        % elements for time k+1, next nmv elements for time k+2, ...        
        % When offline uniform, Wdu is a row vector of nmv elements
        % When offline time-varying, Wdu is a row vector p*nmv elements
        Wdu = updateWeights(Wdu, duwt, nmv, p, ZERO);
    end
    % Recompute Hessian and cost function linear terms
    I2Jm = kron(tril(ones(p,'like',rseq)),eye(nmv,'like',rseq))*Jm;
    [H(1:TotalFreeMoves,1:TotalFreeMoves), Ku1, Kut, Kx, Kv, Kr] = mpc_calculatehessian(...
        Wy, Wu, Wdu, SuJm, I2Jm, Jm, I1, Su1, Sx, Hv, nmv, ny);
    
    if isHyb
        % Transform H in new coordinates
        H(1:TotalFreeMoves,1:TotalFreeMoves)=Mdisinv'*H(1:TotalFreeMoves,1:TotalFreeMoves)*Mdisinv;
    end
    
    % Compute L.  If Hessian is not positive definite, return with -2
    [Linv, BadH] = mpc_checkhessian(H);
    if BadH>ONE32
        u=old_u+uoff;
        for i=1:nmv
            useq(:,i)=u(i);
        end
        status = ONE*-2;
        return
    end
    % Compute Linv and Hinv
    Linv = linsolve(Linv,eye(degrees,'like',rseq),struct('LT',true)); 
    Hinv = Linv'*Linv;
end
% ECR weight is online
if ~no_rhoeps && isQP
    if rhoeps<SMALLHESSIAN
        rhoeps=SMALLHESSIAN; % Forces ECR weight to be > SMALLHESSIAN
    end
    H(degrees,degrees)=rhoeps;
    Hinv(degrees,degrees)=ONE/rhoeps;
    Linv(degrees,degrees)=ONE/sqrt(rhoeps);
end

%% Update constraint bounds, if necessary
if isQP % handle linear constraints, if they have been defined
    % Mlim contains the constant bounds defined in the object (including
    % custom constraint bounds).  See mpc_buildmat.  The following
    % defines the rhs in Ac*z >= Bc iff the bounds remain constant.  Note
    % that custom-constraint bounds depend on offsets.
    Bc = -(Mlim + Mx*xQP + Mu1*old_u + Mv*vseq);
    if ~(no_umin && no_umax && no_ymin && no_ymax)
        Bc0 = -(Mx*xQP + Mu1*old_u + Mv*vseq); % constnt part of Bc
        % The ymax, ymin, umax, and umin bounds can be updated.  Only those
        % with finite bounds in the original QP are eligible.  
        %
        % For each class of constraint and each I/O variable,
        % when the provided bound is a vector with ny (or nmv) entries,
        % the update is compared to the first finite value in the horizon.  The
        % difference, if any, is stored and used to increment all
        % subsequent finite bounds of the same class and variable.  This
        % maintains the constraint profile specified in the design.
        %
        % Otherwise, if the provided bound is has k entries per channel, 2<=k<=p,
        % such entries are used to redefine the corresponding bound, but
        % only for the bounds that were specified as finite in the MPC
        % object. In this case, when k<p the last (k-th) value is repeated over
        % the entire horizon.
        %
        % Mlim=[ymax(:);-ymin(:);umax(:);-umin(:);dumax;-dumin]; (see MPC_BUILDMAT.M)
        %
        % If no on-line bound has been supplied, the one in the object will
        % be used.  Initialize to indicate that increments are yet to be
        % determined.
        
        %% preprocess online OV bound types if present
        % OVMax
        ymax_incr_flag = false(1, ny);
        if int32(numel(ymax))==ny
            % profile-keeping
            ymax_incr = zeros(1, ny, 'like', Mlim);
            nymax = ONE32;              % step 1
        else
            % time-varying
            ymax_incr = [];         % not used
            nymax = int32(size(ymax,1));   % how many steps are provided
        end
        % OVMin
        ymin_incr_flag = false(1, ny);
        if int32(numel(ymin))==ny
            % profile-keeping
            ymin_incr = zeros(1, ny, 'like', Mlim);
            nymin = ONE32;              % step 1
        else
            % time-varying
            ymin_incr = [];         % not used
            nymin = int32(size(ymin,1));   % how many steps are provided
        end
        %% preprocess online MV bound types if present
        % for time-varying bounds, compute block moves column vector from m
        if numel(umax)>nmv || numel(umin)>nmv
            % time-varying
            if int32(numel(m))==ONE32
                moves = [ones(m(1)-1,1,'int32');p-m(1)+1];
            else
                moves = m(:);
            end
        else
            % profile-keeping
            moves=[]; % not used
        end
        % MVMax
        umax_incr_flag = false(1, nmv);
        if int32(numel(umax))==nmv
            % profile-keeping
            umax_incr = zeros(1, nmv, 'like', Mlim);
            numax = ONE32;
        else
            % time-varying
            umax_incr = [];         % not used
            numax = int32(size(umax,1));   % how many steps are provided
            umax = tighten(umax,numax,moves,1); % tighten bounds within block moves
        end
        % MVMin
        umin_incr_flag = false(1, nmv);
        if int32(numel(umin))==nmv
            % profile-keeping
            umin_incr = zeros(1, nmv, 'like', Mlim);
            numin = ONE32;
        else
            % time-varying
            umin_incr = [];         % not used
            numin = int32(size(umin,1));   % how many steps are provided
            umin = tighten(umin,numin,moves,0); % tighten bounds within block moves
        end
        %% Calculate new bounds
        if ~isempty(Mrows) && (Mrows(1)>ZERO32)
            % Make sure Mrows has been defined.  It might have been [] or zero.
            for i = 1:length(Mrows)
                % The bounds for the first 2*(pny+pnu) constraints can be
                % updated.  Get the index of the constraint in the full QP
                % (before removal of redundant and infinite constraints).
                ix = Mrows(i);      % constraint index in full QP
                % Need to scale bound inputs to compute delta bounds
                if ix <= pny
                    % ymax constraint
                    if ~no_ymax
                        % On-line ymax updating possible
                        ii = mod(ix-1, ny) + ONE32;     % OV index for this constraint
                        if nymax==ONE32
                            % profile-keeping
                            [DelBound, ymax_incr(ii), ymax_incr_flag(ii)] = BoundIncrement(-Mlim(i), ...
                                -(RYscale(ii)*ymax(ii) - yoff(ii)), ymax_incr(ii), ymax_incr_flag(ii));
                            Bc(i) = Bc(i) + DelBound;
                        else
                            % time-varying
                            % Mrows only contains indices of constraints that are kept
                            % The number ix-ii = multiple of ny
                            k = (ix-ii)/ny+1; % time index for constraint indexed by ix
                            %DelBound = Mlim(i)-(RYscale(ii)*ymax(min(nymax,k),ii) - yoff(ii));
                            Bc(i) = Bc0(i) - (RYscale(ii)*ymax(min(nymax,k),ii) - yoff(ii));
                        end
                    end
                elseif ix <= 2*pny
                    % ymin constraint
                    if ~no_ymin
                        % On-line ymin updating possible
                        ii = mod(ix-pny-1, ny) + ONE32;     % OV index for this constraint
                        if nymin==ONE32
                            % profile-keeping
                            [DelBound, ymin_incr(ii), ymin_incr_flag(ii)] = BoundIncrement(-Mlim(i), ...
                                (RYscale(ii)*ymin(ii) - yoff(ii)), ymin_incr(ii), ymin_incr_flag(ii));
                            Bc(i) = Bc(i) + DelBound;
                        else
                            % time-varying
                            k = (ix-pny-ii)/ny+1;
                            %DelBound = Mlim(i)+(RYscale(ii)*ymin(min(nymin,k),ii) - yoff(ii));
                            Bc(i) = Bc0(i) + (RYscale(ii)*ymin(min(nymin,k),ii) - yoff(ii));
                        end
                    end
                elseif ix <= 2*pny + pnu
                    % umax constraint
                    if ~no_umax
                        % On-line umax updating possible
                        ii = mod(ix-2*pny-1, nmv) + ONE32;     % MV index for this constraint
                        if numax==ONE32
                            % profile-keeping
                            [DelBound, umax_incr(ii), umax_incr_flag(ii)] = BoundIncrement(-Mlim(i), ...
                                -(RMVscale(ii)*umax(ii) - uoff(ii)), umax_incr(ii), umax_incr_flag(ii));
                            Bc(i) = Bc(i) + DelBound;
                        else
                            % time-varying
                            k = (ix-2*pny-ii)/nmv+1;
                            %DelBound = Mlim(i)-(RMVscale(ii)*umax(min(numax,k),ii) - uoff(ii));
                            Bc(i) = Bc0(i) - (RMVscale(ii)*umax(min(numax,k),ii) - uoff(ii));
                        end
                    end
                elseif ix <= 2*pny + 2*pnu
                    % umin constraint
                    if ~no_umin
                        % On-line umin updating possible
                        ii = mod(ix-2*pny-pnu-1, nmv) + ONE32;     % MV index for this constraint
                        if numin==ONE32
                            % profile-keeping
                            [DelBound, umin_incr(ii), umin_incr_flag(ii)] = BoundIncrement(-Mlim(i), ...
                                (RMVscale(ii)*umin(ii) - uoff(ii)), umin_incr(ii), umin_incr_flag(ii));
                            Bc(i) = Bc(i) + DelBound;
                        else
                            % time-varying
                            k = (ix-2*pny-pnu-ii)/nmv+1;
                            %DelBound = Mlim(i)+(RMVscale(ii)*umin(min(numin,k),ii) - uoff(ii));
                            Bc(i) = Bc0(i) + (RMVscale(ii)*umin(min(numin,k),ii) - uoff(ii));
                        end
                    end
                else
                    break   % remaining constraints have constant bounds
                end
            end
        end
    end
else
    Bc = zeros(size(Mlim),'like',rseq);
end

%% Option to disable the warm start
if ~ASOptions.UseWarmStart
    iA = false(size(iA));
end

%% Compute optimal control action zopt and update iA
[zopt, f, status, iAout] = mpc_solveQP(xQP,nCon,degrees,isQP,Kx,Kr,rseq,Ku1,...
    old_u,Kv,vseq,Kut,utarget,Linv,Hinv,Ac,Bc,iA,CustomSolver,CustomSolverCodeGen,...
    H,UseSuboptimalSolution,QPAlgorithm,ASOptions,IPOptions,ADOptions,MIQPOptions,isHyb,Mdis,Ndis,Vdis,uoff,...
    numdis,maxdis);
if isQP
    slack = zopt(end);
end
% zopt is analytical solution: when isQP is false or when isQP is true and max iterations is 0.  status is 1.
% zopt is 0 (freezing): when status<0 (infeasible QP) or when (status=0 && UseSuboptimalSolution is false).  isQP is true.

%% Compute control action including offset
for i=1:nmv
    u(i)=old_u(i)+zopt(i)+uoff(i);
end

%% Clip control action when sub-optimal solution returned or max iterations is 0
Suboptimal = ((int32(status) == ZERO32 || (int32(status) == -3*ONE32 && isHyb)) && UseSuboptimalSolution);
if (isQP || isHyb) && (Suboptimal || (QPAlgorithm==int32(1) && ASOptions.MaxIterations == ZERO32) || (QPAlgorithm==int32(3) && IPOptions.MaxIterations == ZERO32))
    u = mpc_clipping(u, old_u+uoff, uoff, Mlim, Mrows, p, ny, degrees);
end

%% Compute cost (scalar)
% If QP is feasible and cost is required, compute optimal cost, otherwise
% cost is 0 (not required or QP solver failed)
if return_cost && (int32(status)>ZERO32 || Suboptimal)
    % cost of U and Y
    if int32(length(rseq)) < int32(size(Hv,1))
        aux2=Sx(ny+1:end,:)*xQP + Su1(ny+1:end,:)*old_u + ...
            Hv(ny+1:end,:)*vseq - rseq;
        aux3=I1(nmv+1:end,:)*old_u-utarget;
    else
        aux2=Sx*xQP + Su1*old_u + Hv*vseq - rseq;
        aux3=I1*old_u-utarget;
    end
    if int32(size(Wu,1)) == ONE32
        if length(Wu)==nmv
            cost = aux3'*(kron(ones(p,1,'like',rseq),Wu').*aux3);
        else
            cost = aux3'*(Wu'.*aux3);
        end
    else
        cost = aux3'*Wu*aux3;
    end
    if int32(size(Wy,1)) == ONE32
        if length(Wy)==ny
            cost = cost + aux2'*(kron(ones(p,1,'like',rseq),Wy').*aux2);
        else
            cost = cost + aux2'*(Wy'.*aux2);
        end
    else
        cost = cost + aux2'*Wy*aux2;
    end
    % cost of dU and slack (if any)
    aux=zeros(degrees,1,'like',rseq);
    for i=1:degrees
        aux(i)=zopt(i);
    end
    f=H*aux+2*f;
    cost=cost+aux'*f;
end

%% Compute optimal control sequence (p-by-nu matrix)
% If required, compute optimal control sequence, otherwise 0 (default)
if return_sequence
    % U_optimal = I2*Fm*du_optimal + I1*(old_u + uoff)
    useq(1:p,:) = (reshape(kron(tril(ones(p,'like',rseq)),eye(nmv,'like',rseq))*Jm*zopt(1:(end-isQP),1) + ...
            I1*(old_u(:) + uoff(:)), nmv, p))';
    useq(p+1,:) = useq(p,:);
end

function [Del_Bound, Del_Save, Del_Save_Flag] = BoundIncrement(MlimVal, NewVal, Del_Save0, Del_Save_Flag0)
% If Del_Save_Flag0 is false, compute increment.  Else pass Del_Save0.
if ~Del_Save_Flag0
    Del_Bound = NewVal - MlimVal;
    Del_Save_Flag0 = true;
else
    Del_Bound = Del_Save0;
end
Del_Save = Del_Bound;
Del_Save_Flag = Del_Save_Flag0;

function W = updateWeights(W, signal, channels, p, ZERO)
% W is a row vector, signal is a column vector
if length(W)==channels
    % uniform weights
    for ct = 1:channels
        if signal(ct) < ZERO
            W(ct) = ZERO; %  weight must be non-negative
        else
            W(ct)= signal(ct)^2;
        end
    end
else
    % time-varying weights 
    rows = length(signal); % signal length <= W length
    for ct = 1:rows
        if signal(ct) < ZERO
            W(ct) = ZERO; %  weight must be non-negative
        else
            W(ct)= signal(ct)^2;
        end
    end
    if rows < p*channels
        W((rows+1):end) = repmat(W((rows-channels+1):rows),1,p-rows/channels);
    end
end

function lim = tighten(lim,nlim,moves,ismax)
% Tighten constraints within blocked move
% ismax=1: lim contains upper bounds
% ismax=0: lim contains lower bounds
ONE32 = ones('int32');
ZERO32 = zeros('int32');
cmoves=cumsum(moves);
h=ONE32; % free move index
kh=ZERO32; % step at which blocked move changes
k=ZERO32;
mh=ZERO32;
while k<nlim
    k=k+ONE32;
    if k>kh
        kh=cmoves(h);
        mh=moves(h);
        h=h+ONE32;
    end
    k2=min(k+mh-ONE32,nlim);
    if ismax
        % only lim(k,:) is actually consumed, as Mrows does not contains indices corresponding to moves already blocked
        for j=k+1:k2
            lim(k,:) = min(lim(k,:),lim(j,:));
        end
    else
        for j=k+1:k2
            lim(k,:)=max(lim(k,:),lim(j,:));
        end
    end
    k=k2;
end

