function pwa = hys2pwa(S, verbose, constr, plot)
% HYS2PWA Transformation of HYSDEL model into PWA model
%===============================================================================
%
% Title:        hys2pwa
%                                                             
% Project:      Transformation of HYSDEL model into PWA model
%
% Input:        S: structure containing MLD model generated by HYSDEL compiler 
%               verbose: =0 silent
%                         1 verbose important information
%                         2 verbose detailed information
%               constr: In hysdel, bounds need to be specified for real states 
%                       and inputs.
%                       =0 don't use these bounds as additional constraints
%                       =1 use these bounds as additional constraints. Check them 
%                          a priori (and not a posteriori as constraints in the 
%                          must section) thus speeding up the transformation. Do 
%                          not add them to Hi{i}, Ki{i} in order to have an
%                          efficient representation of the polyhedra.
%                       =2 as 1, but add the bounds to Hi{i}, Ki{i}
%               plot: =0: don't plot anything (=default)
%                     =1: plot resulting polyhedral partition over real state-
%                         input space with bounds as specified in HYSDEL
%                     =2: plot resulting polyhedral partition as in 1 and 
%                         additionally the PWA output functions
%                     Plot 1 and 2 work only for dimensions 1 and 2.
%
%               Remark: The plotting feature is not contained in the publicly 
%               distributed version. However, the files are available from the 
%               authors on request (see below).
%
% Output:       P: structure containing the PWA model, where P{k} is the k-th 
%               PWA submodel with binary state xb = P{k}.xb and binary input 
%               ub = P{k}.ub.
%
%               Let's consider now the PWA dynamics for a given xb and ub:
%
%               The real state/input space is partitioned into a polyhedral 
%               partition with N polyhedra 
%               Hi{i}*[xr; ur] <= Ki{i}, i=1,...,N.
%
%               The PWA dynamics are defined on the polyhedral partition by
%               x' = fx(i) * xr + fu(i) * ur + f0(i)
%               y  = gx(k) * xr + gu(i) * ur + g0(i)
%
%               Additionally, the z and d-variables (as defined in HYSDEL) 
%               are given by
%               z  = Zx(i) * xr + Zu(i) * ur + Z0(i)
%               d  = delta(i)
%
%               P{k}.delta_AD denotes the delta variables, that are defined in the
%               AD section, i.e. that is part of a hyperplane arrangement (in 
%               contrast to delta variables defined in the logic section which
%               are not part of an hyperplane arrangement as they depend via a
%               boolean expression on other delta variables). 
%               If P{k}.delta(i) > 0, then the i-th delta variable has been defined
%               in the AD section and it has computational order 'P{k}.delta(i)'.
%
%               The delta variables of the i-th region are given by P{k}.delta{i}.
%
%               According to Ziegler, we define the markings in the following way:
%               Example: delta = [-1 -1 +1;
%                                 -1 +1 +1] 
%                        implies the existence of 2 hyperplanes generating 3 polyhedra: 
%                               1. A x <= B;
%                               2. A(1,:) <= B(1),  A(2,:) >= B(2);
%                               3. A(1,:) >= B(1),  A(2,:) >= B(2);
%               This is according to the definition of the paper, but contrary to 
%               the way the deltas are internally processed in hys2pwa! In the 
%               very end of hys2pwa, delta is multiplied with (-1) to get the 
%               proper definition.
%
%               In general, the polyhedral partition is not bounded. Bounds
%               resulting from the bounds defined in HYSDEL are given by 
%
% Authors:      Tobias Geyer <geyer@control.ee.ethz.ch>, Fabio D. Torrisi
%
% Papers:       'Efficient Mode Enumeration of Compositional Hybrid Systems',
%               in A. Pnueli & O. Maler, eds. 'Hybrid Systems: Computation 
%               and Control', Vol. 2623 of Lecture Notes in Computer Science,
%               Springer Verlag, pp. 216-232.
%       
%               see also the extended version:
%               'Efficient Mode Enumeration of Compositional Hybrid Systems',
%               Technical Report, AUT03-01, ETH Zurich, 2003
                                                                       
% History:      date        ver.    subject                                       
%               2002.09.20  1.0.0   initial version 
%               2002.11.04  1.0.0   release of beta version
%               2002.11.20  1.0.1   plot>0 does not change Hi and Ki any longer
%               2002.11.29  1.0.2   added progress bar
%               2002.12.09  1.0.3   bug fixed occuring when converting linear systems
%               2003.01.17  1.1.0   features now algebraic loops
%               2003.01.19  1.1.1   code rewritten and documented
%               2003.02.05  1.1.2   improved plotting features
%               2003.03.14  1.1.3   restructured plotANDanalyse
%               2003.05.12  1.2     improved data management (making hys2pwa 2-8 times faster)
%               2003.05.16  1.3     management of hyperplane arrangement computation
%               2003.05.20  1.3.1   option constraining state-input space to bounds specified in hysdel
%               2003.09.23  1.3.2   handles now also degenerated hyperplanes (like 0*[xr; ur] <= b) 
%                                   this new feature is not yet really debugged!!!
%
% Requires:     remAL.m 
%               findFAS.m
%               hyparr.m
%               syminfo.m (part of the HYSDEL package)
%               mb_polyreduce.m 
%               polyinnerball.m
%               lpsolve.m (calls an LP solver like linprog)
%
%               plotANDanalyse requires:
%               PWAwelldefined
%
%
% Contact:      Tobias Geyer
%               Automatic Control Laboratory
%               ETH Zentrum, CH-8092 Zurich, Switzerland
%
%               geyer@aut.ee.ethz.ch
%
%               Comments and bug reports are highly appreciated
%
%===============================================================================



global mptOptions

if ~isstruct(mptOptions)
    mpt_error;
end

clear global pwa
% call mpt_hyparr2 with no input arguments => it will clear it's persistent
% variable, otherwise we get errors on subsequent runs with different hysdel
% models
mpt_hyparr;

% choose LP solver 
lpsolver = mptOptions.lpsolver;
% lpsolver=0: uses E04MBF.m
% lpsolver=1: uses linprog.m (MATLAB)
% lpsolver=2: uses CPLEX


% complement the inputs
if nargin < 1, pwa = '1.3.1'; return; end
if nargin == 1, verbose = 0; constr = 1; plot = 0; end;
if nargin == 2, constr = 1; plot = 0; end;
if nargin == 3, plot = 0; end;

plot = 0;

% check inputs
if (plot > 0) & (S.nxr + S.nur > 2), error('Can plot only over 1- or 2-dimensional real state-input space'); end;
if ~isstruct(S), error('No valid MLD system provided'); end
if ~isfield(S,'MLDstructver'), error('No valid MLD system provided'); end

if verbose >= 0,
    if constr==0, disp('derive PWA representation for x, u \in\RR^n'); end;
    if constr, disp('restrict x-u space to bounds given in the x-u definition in HSYDEL'); end;
    if constr==1, disp('do not add these bounds to the polyhedral representation (efficient representation)'); end;
    if constr==2, disp('add these bounds to the polyhedral representation'); end;
else
    verbose=0;
end

% From now on we assume that the input is corrrect.
if (S.MLDstructver < 2), error('Upgrade to HYSDEL 2.0.x (http://control.ethz.ch/~hybrid/hysdel)'); end

if ~isfield(S, 'B5'), S.B5 = zeros(S.nx,1); end;
if ~isfield(S, 'D5'), S.D5 = zeros(S.ny,1); end;

t_start = clock;


% get the xr-ur domain that is specified by the bounds in HYSDEL,
% where dom.H * [xr; ur] <= dom.K
dom.H = zeros(2*(S.nxr+S.nur), S.nxr+S.nur);
dom.K = zeros(2*(S.nxr+S.nur), 1);
for i=1:length(S.symtable)
    if (S.symtable{i}.kind == 'x') & (S.symtable{i}.type == 'r') 
        r = S.symtable{i}.index;
        dom.H(2*r-1, r) =  1;
        dom.H(2*r  , r) = -1;
        dom.K(2*r-1)    =  S.symtable{i}.max;
        dom.K(2*r)      = -S.symtable{i}.min;
    end
    if (S.symtable{i}.kind == 'u') & (S.symtable{i}.type == 'r') %& ~isfield(S.symtable{i}, 'aux_input')
        r = S.symtable{i}.index + S.nxr;
        dom.H(2*r-1, r) =  1;
        dom.H(2*r  , r) = -1;
        dom.K(2*r-1)    =  S.symtable{i}.max;
        dom.K(2*r)      = -S.symtable{i}.min;
    end
end
dom.constr = constr;


% remove algebraic loops by adding auxiliary inputs and equality constraints
S = remAL(S);
% now, S contains in this order 
%  * nxr real states
%  * nxb binary states
%  * nur real inputs
%  * nub binary inputs,
%  where n?or are original real ones,
%        n?ar are aux. real ones, 
%        n?ob are original binary ones and
%        n?ab are aux. binary ones.
%  in general: nxr = nxor + nxar etc.


% derive map for z- and d-variables to speed up the data management and
% derive computational order
[map, order] = buildMap(S, verbose);


% number of binary combinations for the states and the input
% remark: binary aux. states and inputs are included
nxb_comb = 2^S.nxb;
nub_comb = 2^S.nub;

% initialize wait bar
progr.h = waitbar(0, 'Progress', 'Resize', 'off', 'Name', 'Hysdel to PWA');
progr.t0 = clock;
progr.d_perc = 1/(nxb_comb * nub_comb);
try
    
    % transform into pwa model
    pwa = {};
    k = 1;
    for i = 0:nxb_comb-1,
        xb = (dec2bin(i,S.nxb)-'0')';
        for j = 0:nub_comb-1,
            ub = (dec2bin(j,S.nub)-'0')';
            
            if verbose
                disp('binary state-input combination')
                disp([xb' ub'])
            end;
            
            % fix and check feasibility of the binary state-input combination (xb, ub)
            [feas, S_XUb_fixed] = fix_XUb(S, xb, ub, 1);
            if feas,
                progr.perc = (k - 1) * progr.d_perc;
                
                % if (xb, ub) is feasible, find the PWA representation
                pwa{k} = intern_hys2pwa(S_XUb_fixed, map, order, dom, progr, verbose, lpsolver);
                
                % add binary state-input combination (without the aux. inputs).
                pwa{k}.xb = xb;
                pwa{k}.ub = ub(1:S.nub);
                
                % add bounds of state-input space
                pwa{k}.dom.A = dom.H; 
                pwa{k}.dom.B = dom.K;
                
                % adapt defintion of markings to Ziegler
                for s=1:length(pwa{k}.delta)
                    pwa{k}.delta{s} = (-1)*pwa{k}.delta{s};
                end;
                
                k = k+1;
            end
        end
    end
catch
    drawnow
    close(progr.h);
    drawnow
    error(lasterr);
end
drawnow
close(progr.h);
drawnow
fprintf('Conversion from HYSDEL to PWA form finished (%1.2f sec)\n', etime(clock,t_start))

% plot and analyse the pwa partition
plotANDanalyse(S, pwa, dom.H, dom.K, plot, lpsolver);

if length(pwa) == 1
    pwa = pwa{1};
end;

return




% --------------------------------------------------------------------

function [map, order] = buildMap(S, verbose);

% Build a map that, given the index of a z- or d-variable, holds
% * the number of the entry in the symtable, 
% * the rownumbers of the ineq-constraints that define this variable, and
% * the HYSDEL item in which the variable is defined.
% The purpose of this is to avoid the (extensive) use of the functions 
% 'findineq' and 'syminfo' and to thus speed up the code by a factor of 2.
%
% map is a structure with the elements z and d, 
% where each structure is a cell array:
% e.g. map.z{i}, where i is the index i
%
% then given the index i, map.z{i} contains the following:
% map.z{i}.symtableEntry: number of the entry in 'S.symtable' 
% map.z{i}.definedInRow:  rows of 'S.rowinfo.ineq' that define the variable 
% map.z{i}.itemType:      {Linear, DA, Cont_must, AL_must}
%
% and accordingly for map.d{i}
%
% additionally, map.rows holds for Logic_must, Cont_must and AL_must the 
% rows of the ineq. constraints which are of the respective type of
% constraint
%
% futhermore, the computation order for the z- and d-variables is derived

% get the maximal indices as well as the maximal computation order
maxZindex = 0;
maxDindex = 0;
maxord = 0;
for i=1:length(S.symtable)
    if S.symtable{i}.kind == 'z'
        maxZindex = max(maxZindex, S.symtable{i}.index); 
    elseif S.symtable{i}.kind == 'd'
        maxDindex = max(maxDindex, S.symtable{i}.index); 
    end;
    maxord = max(maxord,S.symtable{i}.computable_order); 
end

% build map.z and map.d
map.z = cell(maxZindex, 1);
map.d = cell(maxDindex, 1);
for i=1:length(S.symtable)
    if S.symtable{i}.kind == 'z'
        
        % index
        index = S.symtable{i}.index;
        if ~isempty(map.z{index}), error('corrupted symbol table'); end;
        map.z{index}.symtableEntry = i;
        
        name = S.symtable{i}.name;
        
        % definedInRow, itemType
        map.z{index}.definedInRow = [];
        map.z{index}.itemType = [];
        for j=1:length(S.rowinfo.ineq)
            if strcmp(name, S.rowinfo.ineq{j}.defines)
                map.z{index}.definedInRow(end+1) = j;
                if isempty(map.z{index}.itemType)
                    map.z{index}.itemType = S.rowinfo.ineq{j}.item_type;
                elseif ~strcmp(map.z{index}.itemType, S.rowinfo.ineq{j}.item_type)
                    % apparently, the z-variable is defined in two different 
                    % sections. This is not allowed. However, if the 2nd
                    % section is the newly (in remAL) introduced AL_must
                    % section, then everything is ok. AL_must constraints
                    % have been introduced in remAL to break the loops
                    if ~strcmp(S.rowinfo.ineq{j}.item_type, 'AL_must')
                        error('z-variable defined in two different sections');
                    end;
                end;
            end;
        end;        
        
    elseif S.symtable{i}.kind == 'd'
        
        % index
        index = S.symtable{i}.index;
        if ~isempty(map.d{index}), error('corrupted symbol table'); end;
        map.d{index}.symtableEntry = i;
        
        name = S.symtable{i}.name;
        
        % definedInRow, itemType
        map.d{index}.definedInRow = [];
        map.d{index}.itemType = [];
        for j=1:length(S.rowinfo.ineq)
            if strcmp(name, S.rowinfo.ineq{j}.defines)
                map.d{index}.definedInRow(end+1) = j;
                if isempty(map.d{index}.itemType)
                    map.d{index}.itemType = S.rowinfo.ineq{j}.item_type;
                elseif ~strcmp(map.d{index}.itemType, S.rowinfo.ineq{j}.item_type)
                    error('d-variable defined in two different sections');
                end;
            end;
        end;        
        
    end;
end;

% build map.rows (with Logic_must, Cont_must and AL_must)
map.rows.Logic_must = [];
map.rows.Cont_must = [];
map.rows.AL_must = [];
for i=1:length(S.rowinfo.ineq)
    if strcmp(S.rowinfo.ineq{i}.item_type, 'Logic_must')
        map.rows.Logic_must(end+1) = i;
    elseif strcmp(S.rowinfo.ineq{i}.item_type, 'Cont_must')
        map.rows.Cont_must(end+1) = i;
    elseif strcmp(S.rowinfo.ineq{i}.item_type, 'AL_must')
        map.rows.AL_must(end+1) = i;
    end;
end;

% build 'order.d' and 'order.z'
% they hold the indices of the d- and z-variables sorted by the 
% computational order
% the indices are defined in S.symtable
order.d = cell(maxord,1);    
order.z = cell(maxord,1);
for i=1:length(S.symtable)
   if S.symtable{i}.kind == 'd'
      order.d{S.symtable{i}.computable_order}(end+1) = S.symtable{i}.index; 
   end
   if S.symtable{i}.kind == 'z'
      order.z{S.symtable{i}.computable_order}(end+1) = S.symtable{i}.index; 
   end
end

% show some information
if verbose
    
    % how many d-variables are defined in the AD section?
    d_ad = 0;
    for i=1:length(map.d)
        if strcmp(map.d{i}.itemType, 'AD'), d_ad = d_ad + 1; end;
    end;
    
    disp('model complexity:')
    fprintf('   # real/binary states: %1i/%1i\n', S.nxr, S.nxb);
    fprintf('   # real/binary inputs: %1i/%1i\n', S.nur, S.nub);
    fprintf('   # d-variables: %1i (%1i induce hyperplane arrangement)\n', S.nd, d_ad);
    
    disp('computational order:')
    for i=1:maxord
        % d:
        fprintf('  %2i: d:', i);
        if length(order.d{i}) > 0
            d_ad = 0;
            for j=1:length(order.d{i})
                if strcmp(map.d{order.d{i}(j)}.itemType, 'AD'), d_ad = d_ad + 1; end;
            end;
            fprintf(' %1i (%1i): ', length(order.d{i}), d_ad);
            for j=1:length(order.d{i})
                fprintf(' %1i', order.d{i}(j));
            end;
        end;
        fprintf('\n');
        
        % z:
        fprintf('      z:');
        for j=1:length(order.z{i})
            fprintf(' %1i', order.z{i}(j));
        end;        
        fprintf('\n');
    end;
end;

disp(' ')

return




% --------------------------------------------------------------------

function [feas, S] = fix_XUb(S, xb, ub, severity)
% fix the binary state-input combination (xb, ub):
%   * update the system matrices:
%     map (xb, ub) into B5, D5 and removing the binary parts of A, B1, C
%     and D1. nxb and nub are set to 0. nx and nu are reduced accordingly.
%   * update inequality constraints:
%     map (xb, ub) into E5 and remove the binary parts of E1 and E4.
% --> the binary states and inputs disappear completely
% 
% check if the system is feasible:
% if severity = 1: check, if the (logic) constraints that depend only 
%                  on (xb, ub) are satisfied for the given combination.
% if severity = 2: same as severity 1 plus feasibility test checking whether 
%                  the relaxed MLD contraints are satisfied.
% if severity = 3: same as severity 2 but performs a mixed integer feasibility 
%                  test instead of the feasibility test.

% update the state-update matrices
Ab  = S.A(:,S.nxr+1:end);
B1b = S.B1(:,S.nur+1:end);    
S.B5 = S.B5 + Ab*xb + B1b*ub;

% and delete the columns of A and B1 corresponding to xb and ub
S.A (:,S.nxr+1:end) = [];
S.B1(:,S.nur+1:end) = [];


% update the output matrices
Cb   = S.C(:,S.nxr+1:end);
D1b  = S.D1(:,S.nur+1:end);
S.D5 = S.D5 + Cb*xb + D1b*ub;

% and delete the columns of C and D1 corresponding to xb and ub
S.C (:,S.nxr+1:end) = [];
S.D1(:,S.nur+1:end) = [];


% update the inequality constraints
E1b  = S.E1(:,S.nur+1:end);    
E4b  = S.E4(:,S.nxr+1:end);
S.E5 = S.E5 + E1b*ub + E4b*xb;

% and delete the columns of E1 and E4 corresponding to xb and ub
S.E1(:,S.nur+1:end) = [];
S.E4(:,S.nxr+1:end) = [];


% update the system dimensions
S.nub = 0; S.nu = S.nur;
S.nxb = 0; S.nx = S.nxr;


% check feasibilty:
% we look for all the lines $rows$ of [E1,E2,E3,E4] where all the elements are 0,
% if for any row \in $rows$, E5(row) < 0 then 0 <= E5(row) -> infeasible
Z1 = all(S.E1 == 0,2);
Z2 = all(S.E2 == 0,2);
Z3 = all(S.E3 == 0,2);
Z4 = all(S.E4 == 0,2);
Z_index = find(all([Z1,Z2,Z3,Z4],2));
if ~isempty(Z_index),
   if any(S.E5(Z_index) < 0),
      feas = 0;
      return
   end
end

if (severity > 1),
    if ~is_feasible([S.E2,S.E3,-S.E1,-S.E4],[S.E5]),
        feas = 0;
        return
    end
end

if (severity > 2),
    warning(sprintf('Severity level %d not yet implemented',severity));
end

feas = 1;

return



% --------------------------------------------------------------------

function pwa = intern_hys2pwa(S, map, order, dom, progr, verbose, lpsolver)
% inputs: 
% S: MLD system, where the binary states and inputs have been fixed. Thus, 
%    the state-input space comprises only real states and real (original and
%    auxiliary) inputs. Also the aux. binary inputs have been fixed. The 
%    constraint matrices together with the state-update and output matrices
%    are updated as we move on to higher computational orders.
% map: map of d- and z-variables to avoid extensive searches in the rowinfo.ineq
% order: computational order of d- and z-variables
% dom: domain, constraints on the xr-ur space
% progr: .perc is the percentage covered by other 
%        .d_perc is the percentage that this conversion will cover
%        .t0 is the time when the conversion started
% verbose: verbosity of the operation
% lpsolver: index of lp solver to be used

global pwa

% initialize the problem
level=1;

% initialize guardlines 
GL.A = [];
GL.B = [];
    
% initialize z-matrices
S.Zx = zeros(S.nz,S.nx);
S.Zu = zeros(S.nz,S.nu);
S.Z0 = zeros(S.nz,1);

% initialize delta marker
S.delta = zeros(S.nd,1);    % this stores the value assigned to delta during 
                            % the computations
S.delta_AD = zeros(S.nd,1); % this is >0 if the delta is coming from a AD 
                            % converter, the numeric value is the 
                            % computational order of the delta variable.

% initialize problem counter
S.counter = 0;

% initialize stack for regions
pwa.Hi = [];
pwa.Ki = [];
pwa.fx = [];
pwa.fu = [];
pwa.f0 = [];
pwa.gx = [];
pwa.gu = [];
pwa.g0 = [];
pwa.Zx = [];
pwa.Zu = [];
pwa.Z0 = [];
pwa.delta = {};

if verbose, fprintf('level = '); end
recursive(S, map, order, level, GL, dom, progr, verbose, lpsolver)
if verbose, fprintf('\n'); end

return




% --------------------------------------------------------------------

function recursive(S, map, order, level, GL, dom, progr, verbose, lpsolver)
% inputs:
%    S: system, where the binary states and inputs have been fixed
%    map: mapping to speed up data management
%    order: computational order of the z- and d-variables
%    level: current computational order
%    GL: the polyhedron we are in defined by the guardlines or constraints 
%        GL.A(i,:)*[xr; ur] <= GL.B(i), for i = 1...rows(GL.A) 
%    progr: .perc is the percentage covered by other 
%           .d_perc is the percentage that this conversion will cover
%           .t0 is the time when the conversion started
%    verbose: verbosity of the operation
% 
% overview:
%   1.) Given the level, represent the z's with computational order 'level'
%       by real states and real inputs. 
%       remark: The z's at a level l do not depend on the deltas at level l.
%               Thus we can first modify the z's and then branch on the 
%               deltas.
%   remark: At a given level, the deltas with indices order_delta(delta)
%           can be computed. Some are uniquely defined by the logic section
%           depending on formerly fixed deltas. The other deltas are not
%           defined. The hyperplane arrangement has to be computed and we 
%           have to branch on the different markings.
%   2.) Fix the deltas defined by the logic section.
%   remark: rem_delta_j are the indices of the deltas, that have are defined 
%           in the the logic section. 
%   3.) Find the hyperplane arrangement and the markings.
%   4.) Go through all the markings. Fix one at a time, update the matrices
%       and move to the next computational order.


global pwa

% progress information (adapt to your needs)
if verbose==1 
    if level <= 3,
        fprintf('%1i, ', level)
    end;
    if level == 3
        fprintf(' #reg=%1i\n', length(pwa.Ki));
    end;
end;
if verbose==2
    fprintf('%1i, ', level)
end;

% represent z by xr and ur:
% z(req_z_i) = Kx * xr + Ku * ur + H 
req_z_i = order.z{level};
if ~isempty(req_z_i),
   S = Zxtract(S, req_z_i, map); 
end; 

% fix the deltas defined in the logic section
req_delta_j = order.d{level};       
if ~isempty(req_delta_j),
   [S, rem_delta_j] = LOGICxtract(S, req_delta_j, map);
else 
   rem_delta_j = [];
end;
% rem_delta_j are the indices of the deltas, that have are defined in the
% the logic section. 

if ~isempty(rem_delta_j),
   
   % get the hyperplane arrangement 'Deltas' and the corresponding guardlines
   % 'GLnew'
   [Deltas, GLnew] = ADxtract(S, rem_delta_j, GL, dom, map, lpsolver, verbose);
   % GL.A, B: all the guardlines up to now
   % GLnew.A, B: new guardlines 
   
   % put the new guardlines together with the current level on the stack
   % (we need this information later for the merging)
   % do we really need it for mering? These guardlines do not seem to be updated
   % when resolving the equality constraints resulting from the removal of loops
   % --> useless information?!?
   pwa.GL{level}.A = GLnew.A;
   pwa.GL{level}.B = GLnew.B;
   
   % number of markings (feasible regions within GL)
   n_expl = size(Deltas,2);
   
   if n_expl > 0
       progr_rec.t0 = progr.t0;
       progr_rec.d_perc = progr.d_perc/n_expl;
       progr_rec.h = progr.h;
   end;
   
   % branch on the regions of the hyperplane arrangement
   for choice = 1:n_expl,
      delta_values = Deltas(:,choice);    % corresponding to req_delta_j
      
      % we are now in the region defined by 'delta_values'
      
      % set hyperplane arrangement according to our choice
      GLmod.A = diag(delta_values) * GLnew.A;
      GLmod.B = diag(delta_values) * GLnew.B;
      
      % remove redundant guardlines
      %GLmod = MinGL(GLmod, Deltas, choice);
      
      % update list of guardlines 
      % (= add the guardlines of the current level to the guardlines of the 
      % previous level)
      GLmod.A = [GL.A; GLmod.A]; 
      GLmod.B = [GL.B; GLmod.B];
      
      % fix deltas
      Smod = S;
      Smod.delta(rem_delta_j) = delta_values;
      Smod.delta_AD(rem_delta_j) = level;
      
      % update ineq. constraint
      Smod.E5 = Smod.E5 - Smod.E2(:,rem_delta_j) * max(0,delta_values);
      Smod.E2(:,req_delta_j) = 0; 
      
      % update state-update function
      Smod.B5 = Smod.B5 + Smod.B2(:,rem_delta_j) * max(0,delta_values);
      Smod.B2(:,req_delta_j) = 0; 
      
      % update output function
      Smod.D5 = Smod.D5 + Smod.D2(:,rem_delta_j) * max(0,delta_values);
      Smod.D2(:,req_delta_j) = 0; 
      
      % move to next computational order
      maxord = length(order.d);
      if level < maxord
         progr_rec.perc = progr.perc + progr_rec.d_perc * (choice - 1);
         recursive(Smod, map, order, level+1, GLmod, dom, progr_rec, verbose, lpsolver)
      else
         error('corrupted computational order');
      end;
   end;
   
else
   
   maxord = length(order.d);
   if level < maxord
      recursive(S, map, order, level+1, GL, dom, progr, verbose, lpsolver)
   else
      % we are done with this branch. Add the results to 'pwa'.
      push(S, GL, map, progr, dom, verbose, lpsolver);
   end;
   
end;

return



% --------------------------------------------------------------------

function S = Zxtract(S, req_z_i, map) 
% inputs:
%    S: internal MLD structure 
%    req_z_i: indizes of z-variables
%
% purpose:
%    represent the z's with indices 'req_z_i' by xr and ur:
%    z(req_z_i) = Kx * xr + Ku * ur + H 
%    and substitute it into the ineq. constraint, state-update and output
%    matrices
%
% comments: 
%    z-variables are either defined via the Linear or the AD section

for i = req_z_i,
    
    % get the rows of the ineq. constraints that define the variable
    % with index i
    row_all = map.z{i}.definedInRow;
        
    % find Kx, Ku and H, 
    % where z(i) = Kx * xr + Ku * ur + H
    
    if strcmp( map.z{i}.itemType, 'Linear' )
        
        % filter out the rows with subindex==2
        row_subindex2 = [];
        for j=1:length(row_all)
            if S.rowinfo.ineq{row_all(j)}.subindex == 2
                row_subindex2(end+1) = row_all(j);
            end;
        end;
        
        Kx = S.E4(row_subindex2,:);
        Ku = S.E1(row_subindex2,:);
        H  = S.E5(row_subindex2,:);
        
        % E2(row,:) has to be zero 
        if sum(abs(S.E2(row_subindex2,:))) > 0, 
            error('item-type=Linear: the E2-row has to be 0');
        end;
        
    elseif strcmp( map.z{i}.itemType, 'DA' )
        
        % we have 4 inequalities, 2 ineq. are the same. Find these two.
        
        % the sign-vector of z(i)
        sign = S.E3(row_all,i);
        
        % correct Kx, Ku and H by the sign-vector (element-wise multiplication)
        Kx = diag(sign) * S.E4(row_all,:);
        Ku = diag(sign) * S.E1(row_all,:);
        H  = diag(sign) * S.E5(row_all,:);
        
        % find the 2 ineq., that are the same
        one_two = 0; three_four = 0;
        if (all(toleq(Kx(1,:),Kx(2,:)))) & (all(toleq(Ku(1,:),Ku(2,:)))) & (toleq(H(1),H(2))),
            one_two = 1;
            Kx_n=Kx(1,:); Ku_n=Ku(1,:); H_n=H(1);
        end
        if (all(toleq(Kx(3,:),Kx(4,:)))) & (all(toleq(Ku(3,:),Ku(4,:)))) & (toleq(H(3),H(4))),
            three_four = 1;
            Kx_n=Kx(3,:); Ku_n=Ku(3,:); H_n=H(3);
        end
        if ~xor(one_two,three_four),
            error('item-type=DA: choice ambigeous or no choice fits');
        end       
        Kx = Kx_n;
        Ku = Ku_n;
        H = H_n;
        
    else
        
        error('Zxtract: unknown case.');
        
    end;  
    
    
    % replace z(i) by z(i) = Kx * x + Ku * u + H
    
    % update ineq. constraint matrices
    S.E1 = S.E1 - S.E3(:,i) * Ku;
    S.E4 = S.E4 - S.E3(:,i) * Kx;
    S.E5 = S.E5 - S.E3(:,i) * H;
    S.E3(:,i) = 0;
    
    % update state-update matrices
    S.A  = S.A  + S.B3(:,i) * Kx;
    S.B1 = S.B1 + S.B3(:,i) * Ku;
    S.B5 = S.B5 + S.B3(:,i) * H;
    S.B3(:,i) = 0;
    
    % update output matrices
    S.C  = S.C  + S.D3(:,i) * Kx;
    S.D1 = S.D1 + S.D3(:,i) * Ku;
    S.D5 = S.D5 + S.D3(:,i) * H;  
    S.D3(:,i) = 0;
    
    % update Z matrices
    S.Zx(i,:) = Kx;
    S.Zu(i,:) = Ku;
    S.Z0(i,:) = H;  
    
end;

return



% --------------------------------------------------------------------

function [S, rem_d_i] = LOGICxtract (S, req_d_i, map) 
% inputs:
%    S:  internal MLD structure 
%    req_d_i: indices of delta variables
%
% outputs:
%    S: internal MLD structure with the 'updated' matrices
%    rem_d_i: indices of deltas that we were not able to fix 
%             (because they are not defined by the logic section)
%
% purpose:
%    Determine all deltas(i) with index i \in req_d_i and which are 
%    defined in the logic section. Substitute them into the ineq. constraint,
%    state-update and output matrices.

rem_d_i = [];
for i = req_d_i,
    % is the d-variable with index i defined in the Logic section?
    if ~strcmp(map.d{i}.itemType, 'Logic'),
        % the delta is not defined in a logic statement. Add it to the deltas,
        % that we can not fix in this subfunction (remaining deltas)
        rem_d_i(end+1) = i;
    else    
        % get the rows of the ineq. constraints that define the variable
        % with index i
        row = map.d{i}.definedInRow;
        
        one_feasible = 0;
        zero_feasible = 0;
        if all(S.E2(row,i) <= S.E5(row,1))
            one_feasible = 1;
        end
        if all(0 <= S.E5(row,1))
            zero_feasible = 1;
        end
        if ~xor(one_feasible, zero_feasible)   
            % 0 and 1 are either both feasible or both infeasible
            error(sprintf('item-type=Logic: Delta %d can not be fixed.',i))
        else 
            % either 0 or 1 is feasible (but not both)
                                    
            % update ineq. constraint
            S.E5 = S.E5 - S.E2(:,i) * one_feasible;
            S.E2(:,i) = 0;
            
            % update state-update function
            S.B5 = S.B5 + S.B2(:,i) * one_feasible;
            S.B2(:,i) = 0;
            
            % update output function
            S.D5 = S.D5 + S.D2(:,i) * one_feasible;
            S.D2(:,i) = 0;
            
            % fix delta variable
            if one_feasible == 0, S.delta(i) = -1;
            else                  S.delta(i) = +1;
            end;
        end
    end
end

return



% --------------------------------------------------------------------

function [Deltas, GLnew] = ADxtract(S, req_d_i, GL, dom, map, lpsolver, verbose) 
% purpose:
%   In the AD section thresholds or guardlines on real variables are specified 
%   that set the corresponding delta variables to zero if the threshold is not
%   fulfilled or to one, respectively. The thresholds given by 'req_d_i'
%   depend only on xr and ur and they are defined within the current polyhedron
%   given by GL. First, we extract these guardlines or thresholds. They are 
%   denoted by GLnew. Using the new guardlines, we next compute the hyperplane
%   arrangement 'Deltas' within the polyhedron given by GL. 
%
% inputs:
%    S: internal MLD structure
%    req_d_i: indizes of delta variables
%    GL: guardlines of the polyhedron, where the j-th guardline is given by 
%        GL.A(j,:) * [xr; ur] <= GL.B(j)
%
% outputs:
%    Deltas: hyperplane arrangement (or markings of the regions generated by 
%            the hyperplanes), where
%            Deltas(:,j) is the marking of the j-th polyhedron of the 
%            arrangement. This marking is a {0,1} vector. His indices are given 
%            by 'req_d_i'.
%            Example: delta = [1  1 -1;
%                              1 -1 -1] 
%            implies the existence of 2 hyperplanes generating 3 polyhedra: 
%               1. A x <= B;
%               2. A(1,:) <= B(1),  A(2,:) >= B(2);
%               3. A(1,:) >= B(1),  A(2,:) >= B(2);
%   GLnew: guardlines corresponding to the delta variables with indices 
%          'req_d_i', where the j-th guardline is given by 
%          GLnew.A(j,:) * [xr; ur] <= GLnew.B(j)
%
% remark: Note, that the hyperplane arrangement corresponds to the 
%         d-variables with the indices 'req_d_i', and that it is only 
%         computed for the polyhedron {w: GL.A*w <= GL.B}, where w=[xr; ur].


% extract guardlines
GLnew.A = [];
GLnew.B = [];
for i = req_d_i,
    
    % get the rows of the ineq. constraints that define the d-variable
    % with index i
    row_all = map.d{i}.definedInRow; 
    
    % filter out the rows with subindex==2
    row_subindex2 = [];
    for j=1:length(row_all)
        if S.rowinfo.ineq{row_all(j)}.subindex == 2
            row_subindex2(end+1) = row_all(j);
        end;
    end;
    row = row_subindex2;
    
    if isempty(row), error('item-type=AD: variable definition not found'); end
    
    % the row with subindex=2 is: a*[xr; ur] - b <= M*( 1-d(i) )
    
    % E3(row,:) has to be zero
    if sum(abs(S.E3(row,:))) > 0, error('item-type=AD: the E3-row has to be 0');  end;
    
    % extract the guardline a * [xr; ur] -b <= 0, 
    % where a = [-E4r(row) -E1r(row)]
    %       b = E5(row)-M, M=E2(row)
    % remark: the guardline does not depend on z, as all the z are functions 
    %         of x and u
    a = [-S.E4(row,1:S.nx), -S.E1(row,1:S.nu)]; 
    b = S.E5(row)-S.E2(row,i);
    
    GLnew.A = [GLnew.A; a];
    GLnew.B = [GLnew.B; b];
end;

% check dimensions of new constraints
if (size(GLnew.A,2) ~= S.nx+S.nu) | (size(GLnew.A,1) ~= length(GLnew.B)) 
    error('dimensions of new guardlines are wrong');
end;

% check which hyperplanes are well-defined (have at least one non-zero element)
wellDef=[]; illDef=[];
for i=1:length(req_d_i)
    if all(GLnew.A(i,:)==0)
        % all elements in a are zero, i.e. 0*[xr; ur] <= b 
        var_name = S.symtable{map.d{req_d_i(i)}.symtableEntry}.name;
        fprintf('hyperplane corresponding to %s is badly defined.\n', var_name);
        illDef(end+1) = i;
    else
        wellDef(end+1) = i;
    end;
end;
% the hyperplanes with indices wellDef are well-defined

% now, we split the problem into two subproblems:
% 1.) get the hyperplane arrangement for well-defined hyperplanes
% 2.) determine the deltas for the ill-defined hyperplanes
% 3.) merge everthing into Deltas

% 1.) get hyperplane arrangement (for the well-defined hyperplanes)
if wellDef
    GLnewWell.A = GLnew.A(wellDef,:);
    GLnewWell.B = GLnew.B(wellDef);
    DeltasWell = mpt_hyparr(GLnewWell, GL, dom, lpsolver, verbose);
else
    DeltasWell = [];
end;

% 2.) get the deltas of the ill-defined hyperplanes
if illDef
    % determine, if the delta is true (+1) or false (-1)
    GLnewIll.A = GLnew.A(illDef,:);
    GLnewIll.B = GLnew.B(illDef);
    
    % the number of regions
    nR = max(1, length(DeltasWell));
    
    % as the ill-defined hyperplane is given by 0*[xr; ur] <= b, we only
    % need to check if b is positive or negative.
    % b>=0: delta=+1 (true)
    % b<0:  delta=-1 (false)
    for i=1:length(GLnewIll.B)
        value = (GLnewIll.B(i)>=0) *2-1;
        DeltasIll(i,:) = value*ones(1,nR); 
    end;
end;

% 3.) merge both solutions
if ~isempty(wellDef) & isempty(illDef)
    Deltas = DeltasWell;
elseif isempty(wellDef) & ~isempty(illDef)
    Deltas = DeltasIll;
elseif ~isempty(wellDef) & ~isempty(illDef)
    Deltas = NaN*ones(length(req_d_i), nR);
    Deltas(wellDef,:) = DeltasWell;
    Deltas(iLLDef,:)  = DeltasIll;
else
    error('ADxtract: it seems that there is no hyperplane')
end;

return



% --------------------------------------------------------------------

function push(S, GL, map, progr, dom, verbose, lpsolver)
% purpose:
%   S contains exactely one polyhedron with one pwa dynamics. This is one region
%   of the partition. Add this one to the global variable 'pwa'.
%
% inputs:
%    S: internal MLD structure
%    GL: guardlines of the polyhedron, where the i-th guardline is given by 
%        GL.A(i,:) * [xr; ur] <= GL.B(i)
%    progr: .perc is the percentage covered by other 
%           .d_perc is the percentage that this conversion will cover
%           .t0 is the time when the conversion started
%    verbose: verbosity of the operation
%
% overview:
%   1.) add constraints of the MUST section
%   2.) add equality constraints resulting from the opening of algebraic loops
%   3.) remove redundant constraints (or guardlines)
%   4.) add to 'pwa'
%

global pwa


% check:

% the matrices E2 and E3 have to be completely zero
if any(any(S.E2)) > 0, error('push: E2 has to be zero'); end;
if any(any(S.E3)) > 0, error('push: E3 has to be zero'); end;

% the matrices B2 and B3 have to be completely zero
if any(any(S.B2)) > 0, error('push: B2 has to be zero'); end;
if any(any(S.B3)) > 0, error('push: B3 has to be zero'); end;

% the matrices D2 and D3 have to be completely zero
if any(any(S.D2)) > 0, error('push: D2 has to be zero'); end;
if any(any(S.D3)) > 0, error('push: D3 has to be zero'); end;


% have all d-variables been fixed? (are -1 or 1 and not as initially 0?)
if find(S.delta==0), error('push: d-variable(s) not fixed'); end;


% check logic constraints
% (the rows of the logic constraints are given by map.rows.Logic_must)
if any(S.E5(map.rows.Logic_must)<0),
    %disp('push: region does not satisfy the logic must constraints');
    % logic constraints should be checked in the very beginning when fixing the
    % binary state-input combination to not slow down the code.
    return
end;

% add constraints of the MUST section (if specified)
rows = map.rows.Cont_must;
if ~isempty(rows)
    % We are looking for constraints on xr and ur: Kx*xr + Ku*ur <= H.
    % At this point, there should be no constraints on z and delta anymore.
    % This has been checked already a few lines above...
    Kx = -S.E4(rows,:);
    Ku = -S.E1(rows,:);
    H  =  S.E5(rows);
    
    % add the constraints
    GL.A = [GL.A; [Kx Ku]];
    GL.B = [GL.B; H];
end;

% add equality constraints resulting from the opening of algebraic loops
[S, GL, ALsolved] = resolveEqConstr(S, GL, map, verbose);

% if set, add bounds to the polyhedral representation
if dom.constr == 2
    GL.A = [GL.A; dom.H];
    GL.B = [GL.B; dom.K];
end;

if ALsolved
    % algebraic loop has been solved
    
    % remove redundant constraints
    % this yields Hi*[xr; ur] <= Ki
    if ~isempty(GL.A) & ~isempty(GL.B)
        %[Hi, Ki, isemptypoly, keptrows] = mb_polyreduce(GL.A, GL.B, 1e-6, lpsolver, 1);
        Poly = polytope(GL.A, GL.B);
        [Hi, Ki] = double(Poly);
        isemptypoly = ~isfulldim(Poly);
    else
        % GL.A and GL.B are empty, i.e. we have only one region encompassing
        % the whole xr-ur-space
        Hi = [];
        Ki = [];
        isemptypoly = 0;
    end;
else
    % algebraic loops has not been solved

    isemptypoly = 1;
end;

if verbose==2
    if isemptypoly,
        fprintf('-,\n        ');
    else
        fprintf('*,\n        ');
    end  
end

t = etime(progr.t0,clock);
% tt total time
% t : tt = perc + d_perc : 1; => tt = t / (perc + d_perc)
% time remaining 
tr = t - t / (progr.perc + progr.d_perc);
moretime = ['Progress (Time remaining: ', num2str(tr,'%.0f'), ' s)'];
waitbar(progr.perc + progr.d_perc,progr.h,moretime)

% add region with the pwa dynamics to the solution variable 'pwa'
if ~isemptypoly,
    pwa.delta{end+1} = S.delta;
    pwa.Hi{end+1} = Hi;
    pwa.Ki{end+1} = Ki;
    pwa.fx{end+1} = S.A;
    pwa.fu{end+1} = S.B1;
    pwa.f0{end+1} = S.B5;
    pwa.gx{end+1} = S.C;
    pwa.gu{end+1} = S.D1;
    pwa.g0{end+1} = S.D5;
    pwa.Zx{end+1} = S.Zx;
    pwa.Zu{end+1} = S.Zu;
    pwa.Z0{end+1} = S.Z0;
    pwa.delta_AD = S.delta_AD; % this is clearly the same for all the regions
end;

return



% --------------------------------------------------------------------

function [S, GL, ALsolved] = resolveEqConstr(S, GL, map, verbose);

% The algebraic loops were removed by adding auxiliary inputs and adding
% additional equality constraints with the label 'AL_must'. Now, we use these 
% equality constraints to determine the auxiliary inputs, i.e. express the 
% auxiliary inputs as functions of the xr and ur. As a result, the auxiliary 
% inputs are removed and the input space is reduced to its original size.

% u is structured in the following way: u = [ur; uar; ub; uab], 
% where uor = original real inputs      S.nuor
%       uar = aux. real inputs          S.nuar
%       uob = original binary inputs    S.nuob has to be zero
%       uab = aux. binary inputs        S.nuab has to be zero
% and S.nur = S.nuor + S.nuar

if S.nub > 0, error('resolveEqConstr: the number of binary inputs has to be 0'); end;

rows = map.rows.AL_must;
if ~isempty(rows)
        
    % take only every second row, 
    % as the equality constraints are represented as pairs of ineq. constraints
    rows = rows(1:2:end);
    
    % derive the equality constraints 
    % [Ko Ka]*[uor; uar] + L*xr = H
    % from 0 <= S.E1*u + S.E4*x + S.E5
    Ko = -S.E1(rows, 1:S.nuor);
    Ka = -S.E1(rows, S.nuor+1:S.nur);
    L  = -S.E4(rows, :);   % columns corresponding to binary states have been removed
    H  =  S.E5(rows);
    
    if det(Ka) ~= 0
        
        % the determinant of Ka is not zero
        % --> we can derive an expression for uar
    
        % solve for uar:
        % uar = inv(Ka) * (H - L*xr - Ko*uor)
        %     = Fu * uor + Fx * xr + F0
        Fu = -inv(Ka) * Ko;
        Fx = -inv(Ka) * L;
        F0 =  inv(Ka) * H;
        
        % decompose GL.A and GL.B defining the polyhedron:
        % [Hx Hu_or Hu_ar] * [xr; uor; uar] <= K
        Hx    = GL.A(:, 1:S.nxr);
        Hu_or = GL.A(:, S.nxr+1:S.nxr+S.nuor);
        Hu_ar = GL.A(:, S.nxr+S.nuor+1:S.nxr+S.nur);
        if size([Hx Hu_or Hu_ar],2) ~= size(GL.A,2)
            error('resolveEqConstr: error when decomposing GL.A')
        end;
        
        % update GL.A and GL.B
        % using uar = Fu * uor + Fx * xr + F0
        % yields Hx*xr + Hu_or*uor + Hu_ar*(Fx*xr+Fu*uor+F0) <= K
        GL.A = [Hx + Hu_ar*Fx  Hu_or + Hu_ar*Fu];
        GL.B = GL.B - Hu_ar*F0;
        
        % update x' = S.A * xr + S.B1 * ur + S.B5
        B1_r  = S.B1(:,1:S.nuor);
        B1_ar = S.B1(:,S.nuor+1:end);
        
        S.A  = S.A  + B1_ar*Fx;
        S.B1 = B1_r + B1_ar*Fu;
        S.B5 = S.B5 + B1_ar*F0;
        
        % update y = S.C * xr + S.D1 * ur + S.D5
        D1_r  = S.D1(:,1:S.nuor);
        D1_ar = S.D1(:,S.nuor+1:end);
        
        S.C  = S.C  + D1_ar*Fx;
        S.D1 = D1_r + D1_ar*Fu;
        S.D5 = S.D5 + D1_ar*F0;
        
        % update z  = S.Zx * xr + S.Zu * ur + S.Z0
        Zu_r  = S.Zu(:,1:S.nuor);
        Zu_ar = S.Zu(:,S.nuor+1:end);
        
        S.Zx = S.Zx + Zu_ar*Fx;
        S.Zu = Zu_r + Zu_ar*Fu;
        S.Z0 = S.Z0 + Zu_ar*F0;    
    
    else
        
        % the determinant of Ka is zero
        % --> we can not derive an expression for uar as there are either
        %     no or an inf. number of solutions
        
        ALsolved = 0;
        %warning('det(Ka) = 0, couldn''t resolve algebraic loop for a polyhedron')
        return;
        
    end;
    
end;

ALsolved = 1;
return



% --------------------------------------------------------------------

function C = toleq(A,B)
% returns 1 if A-B is below a given tolerance
tol = 1e-8;  % 1e-10
C = abs(A-B) <= tol;
return




% --------------------------------------------------------------------

function plotANDanalyse(S, pwa, Hdom, Kdom, plot, lpsolver);


% plot over the domain Hdom * [xr; ur] <= Kdom
if plot > 0
    
     for k=1:length(pwa)
        % the k-th pwa partition
        Hi = pwa{k}.Hi;
        Ki = pwa{k}.Ki;
        
        % and the k-th output functions
        gx = pwa{k}.gx;
        gu = pwa{k}.gu;
        g0 = pwa{k}.g0;
        
        % add plotting boundaries to the regions
        Hi_p = {};  Ki_p = {};
        for i = 1:length(Hi);
            [Hi_p{i}, Ki_p{i}] = mb_polyreduce([Hi{i};Hdom], [Ki{i};Kdom], 1e-6, lpsolver);
        end
        
        % plot
        if S.nxr + S.nuor == 1
            if plot==1, tg_plotRegions1dim(Hi_p, Ki_p);
            else        tg_plotRegions1dim(Hi_p, Ki_p, gx, gu, g0);
            end;
        elseif S.nxr + S.nuor == 2
            if plot==1, tg_plotRegions(Hi_p, Ki_p);
            else        tg_plotRegions(Hi_p, Ki_p, gx, gu, g0);
            end;
        else
            error('can not plot over this dimension');
        end;
        %title('partition of the x-u space')
    end;
end;



% if the original problem contains algebraic loops:
if isfield(S, 'algLoop')
    disp('Checking well-definedness of the PWA system:')
    PWAwelldefined(pwa, Hdom, Kdom, lpsolver, plot);
end;

return


% --------------------------------------------------------------------

function out=is_feasible(A,B)
% returns true is \exists x s.t. Ax \le B

nv=size(A,2);
nc=size(A,1);

x = mpt_solveLP([zeros(1,nv) 1]', [A -ones(nc,1); zeros(1,nv) -1], [B; 0]);
out=~([zeros(1,nv),1]*x(:));
