function [u, status, region, xk, xest] = mpcblock_pwa_evaluate(old_x, r, v, ym, ext_u, old_u, no_mv, nu, nx, nv, uoff, ...
    openloopflag, F, G, H, K, i1, i2, nr, np, xunfindex, refunfindex, uunfindex, mdunfindex,...
    thmin, thmax, CustomEstimation, RMYscale, myoff, xoff, A, Bu, Bv, Cm, Dvm, M, L, ...
    switch_in, switch_inport, no_switch, enable_value)
% Internal use (by "PWA function evaluation" MATLAB Function block)
%
% Note: r,v,old_u are scaled/offset-free, the output signal u is instead
% with offset and scaling.

% Author:  A. Bemporad
% Copyright 1994-2013 The MathWorks, Inc.

%#codegen

% Rules on returned MV signal
% 1. In normal operations, return with offset
% 2. When open-loop flag is on, return old_u + uoff (required by design from scratch or findop or linearize)

ZERO=zeros(1,1,'like',r);
ONE=ones(1,1,'like',r);

%% Initialize default outputs
u=zeros(nu,1,'like',r);
xk=zeros(nx,1,'like',r);
xest=zeros(nx,1,'like',r);
status=ONE;
region=ZERO;
v_aux = ones(nv,1,'like',r);
for ct=1:length(v)
    v_aux(ct) = v(ct);
end

%% if open loop, return u = old_u + uoff
if openloopflag
    u=old_u+uoff;
    return
end

% External MV override.
% NOTE: old_u and ext_mv input signals are dimensionless and offset-free.
if no_mv==ONE
    delmv = zeros(nu,1,'like',r);
else
    delmv = ext_u - old_u;
    old_u = ext_u;
end

% Obtain x[k|k]
if CustomEstimation
    % Remove offset from input state, which is x(k|k)
    xest=ym-xoff;
    y_innov = ZERO; % dummy for mex
else
    % Default state estimation.
    % Scale measured output and remove offset.
    ym=RMYscale.*ym-myoff;
    % Correct x(k|k-1) for possible external mv override.
    % NOTE:  Offset was removed from x[k|k-1] at k=0.
    xk = old_x + Bu*delmv;
    % Measurement upate to x(k|k)
    ym_est = Cm*xk + Dvm*v_aux;
    y_innov = ym - ym_est;
    xest = xk + M*y_innov;
end

evaluation_disabled = switch_inport==1 && no_switch~=1 && switch_in~=enable_value;
if ~evaluation_disabled
    
    % Prepare vector of parameters
    th=zeros(np,1,'like',r);
    i=int32(1);
    j=int32(0);
    if xunfindex(1)~=ZERO % In case it's empty, it gets stored as '0'
        j=int32(numel(xunfindex));
        for h=int32(0):j-int32(1)
            th(i+h)=xest(xunfindex(h+int32(1)));
        end
        i=j+int32(1);
    end
    if refunfindex(1)~=ZERO
        k=int32(numel(refunfindex));
        j=j+k;
        for h=int32(0):k-int32(1)
            th(i+h)=r(refunfindex(h+int32(1)));
        end
        i=j+int32(1);
    end
    if uunfindex(1)~=ZERO
        k=int32(numel(uunfindex));
        j=j+k;
        for h=int32(0):k-int32(1)
            th(i+h)=old_u(uunfindex(h+int32(1)));
        end
        i=j+int32(1);
    end
    if mdunfindex(1)~=ZERO
        k=int32(numel(mdunfindex));
        for h=int32(0):k-int32(1)
            th(i+h)=v(mdunfindex(h+int32(1)));
        end
    end
    
    for i=ONE:np
        if th(i)<thmin(i) || th(i)>thmax(i)
            status=ZERO;
            region=ZERO;
        end
    end
    
    du=zeros(nu,1,'like',r);  % Keep du=0 in case status==ZERO
    
    if abs(status-ONE)<eps(ONE)
        found=false;
        region=ZERO;
        regionmin=ZERO;
        %betamin=Inf(1,1,'like',r);
        betamin=realmax(class(r));
        while ~found && region<nr
            region=region+ONE;
            h=i1(int32(region));
            beta=-ONE;  % Compute beta=max(pwafun(i).H*th-pwafun(i).K)
            while h<=i2(int32(region))
                betah=H(h,:)*th-K(h);
                if betah>beta
                    beta=betah;
                end
                h=h+1;
            end
            if beta<=0 % th belongs to current region 'region'
                found=true;
            else
                if betamin>beta
                    betamin=beta;
                    regionmin=region;
                end
            end
        end
        
        if ~found
            region=regionmin;
            status=-ONE;
        end
        
        k=(int32(region)-int32(1))*int32(nu);
        for h=int32(1):int32(nu)
            du(h)=F(h+k,:)*th+G(h+k);
        end
        
    end
    
    % Compute control action including offset
    for i=int32(1):int32(nu)
        u(i)=old_u(i)+du(i)+uoff(i);
    end
else
    u=old_u+uoff;
end

if CustomEstimation
    xk=old_x;
else
    % update x[k+1|k], assuming that above u and v will be applied.
    xk = A*xk + Bu*(u-uoff) + Bv*v_aux + L*y_innov;
end

% return xest in original value
xest = xest + xoff;

