diff --git a/@memZono/@hybZono/getLeaves.m b/@memZono/@hybZono/getLeaves.m
deleted file mode 100644
index d067d36..0000000
--- a/@memZono/@hybZono/getLeaves.m
+++ /dev/null
@@ -1,54 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Method:
-%       Return the binary factors \xib \in {-1,1}^nGb that result in 
-%       non-empty constrained zonotopes 
-%   Syntax:
-%       [leaves] = getLeaves(Z,optSolver)
-%   Inputs:
-%       Z - 1D, 2D, or 3D hybrid zonotope in HCG-Rep (hybZono object)
-%       optSolver - solver options needed for mixed-integer linear propgrams
-%   Outputs:
-%       leaves - nGb x nLeaves matrix, each column denoting the \xi vector
-%                corresponding to one of the nLeaves non-empty constrained
-%                zonotopes
-%   Notes:
-%       If \xib denotes the i^th column of leaves, then the corresponding
-%       non-empty constrained zonotope is of the form:
-%       Z = { (c + Gb \xib) + Gc \xic | ||\xic||_inf <= 1, Ac \xi = b - Ab \xib }
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-function [leaves] = getLeaves(obj,optSolver)
-
-% Full tree when there are no constrains (nLeaves = 2^nGb)
-if isempty(obj.b)
-	combos = ff2n(obj.nGb);
-	combos(combos==0) = -1;
-	leaves = combos';
-	return
-end
-
-% Continue if set has equality constraints
-if nargin < 2 || isempty(optSolver)
-    optSolver = solverOptions;
-end
-
-if ~strcmp(optSolver.milpSolver,'gurobi')
-    error('Gurobi is currently required for this function.')
-end
-
-% Get up to 2^nGb solutions
-optSolver.nSolutions = 2^obj.nGb;
-optSolver.MIPFocus = 0; % Find feasible solutions
-% Problem data for mixed-integer linear program (MILP)
-Aeq = [obj.Ac 2*obj.Ab];
-beq = obj.b+obj.Ab*ones(obj.nGb,1);
-lb = -ones(obj.nGc+obj.nGb,1);
-ub =  ones(obj.nGc+obj.nGb,1);
-lb((obj.nGc+1):end) = 0; % Set lower bound on binary factors to zero
-vType(1:obj.nGc) = 'C';
-vType(obj.nGc+1:(obj.nGc+obj.nGb)) = 'B';
-[x,~,~] = solveMILP([],[],[],Aeq,beq,lb,ub,vType,optSolver);
-
-% convert binaries back to {-1 1}
-leaves = (x((obj.nGc+1):end,:)-0.5)/0.5;
-
-end
\ No newline at end of file
diff --git a/@memZono/@hybZono/hybZono.m b/@memZono/@hybZono/hybZono.m
deleted file mode 100644
index 2604616..0000000
--- a/@memZono/@hybZono/hybZono.m
+++ /dev/null
@@ -1,151 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Class:
-%       Hybrid zonotope of the form:
-%       Z = { c + Gc \xic +  Gb \xib | ||\xic||_inf <= 1, \xib \in {-1,1}^nGb, Ac \xic + Ab \xib = b }
-%   Syntax:
-%       Z = hybZono(Gc,Gb,c,Ac,Ab,b)
-%       Z = hybZono(z)
-%       Z = hybZono(H_collection)
-%       Z = hybZono({V,M})
-%   Inputs:
-%       Gc - n x nGc matrix to define hybrid zonotope in R^n with nGc continuous generators
-%       Gb - n x nGb matrix to define hybrid zonotope in R^n with nGb binary generators
-%       c  - n x 1 vector to define center
-%       Ac - nC x nGc matrix to define nC equality constraints (Ac \xic + Ab \xib = b)
-%       Ab - nC x nGb matrix to define nC equality constraints (Ac \xic + Ab \xib = b)
-%       b  - nC x 1 vector to define nC equality constraints (Ac \xic + Ab \xib = b)
-%       z - zono or conZono object to be recast as hybrid zonotope
-%       H_collection - N x 1 cell array with i^th cell [H_i f_i] defining H-rep set H_i x <= f_i
-%       V - n x nV matrix with each column defining a vertex in R^n
-%       M - nV x N matrix indicating which of the N sets the nV vertices belong
-%   Outputs:
-%       Z - hybrid zonotope as hybZono object
-%   Notes:
-%       Inherits methods from abstractZono class
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-classdef hybZono < abstractZono
-
-    properties
-        Gc      % Continuous generator matrix (n x nGc)
-        Gb      % Binary generator matrix (n x nGb)        
-        c       % Center (n x 1)
-        Ac      % Continuous constraint matrix (nC x nGc)
-        Ab      % Binary constraint matrix (nC x nGb)        
-        b       % Constraint vector (nC x 1)
-    end
-    properties (Dependent) % These properties get automatically updated when used
-        n       % Dimension
-        nGc     % Number of continuous generators
-        nGb     % Number of binary generators
-        nC      % Number of constraints
-    end
-    properties (Hidden) % Unused properties from superclass
-        G       % Generator matrix (n x nG)
-        A       % Constraint matrix (nC x nG)
-        nG      % Number of generators
-    end
-
-    methods
-        % Constructor
-        function obj = hybZono(varargin)
-            if nargin == 6
-                obj.Gc = varargin{1};
-                obj.Gb = varargin{2};
-                obj.c  = varargin{3};
-                obj.Ac = varargin{4};
-                obj.Ab = varargin{5};
-                obj.b  = varargin{6};
-                if sum(size(obj.Gc)) == 0
-                    obj.Gc = zeros(size(obj.c,1),0);
-                end
-                if sum(size(obj.Gb)) == 0
-                    obj.Gb = zeros(size(obj.c,1),0);
-                end
-                if sum(size(obj.Ac)) == 0
-                    obj.Ac = zeros(0,size(obj.Gc,2));
-                end
-                if sum(size(obj.Ab)) == 0
-                    obj.Ab = zeros(0,size(obj.Gb,2));
-                end
-            elseif nargin == 1
-                switch class(varargin{1})
-                    case 'zono'
-                        obj.Gc = varargin{1}.G;
-                        obj.Gb = zeros(size(obj.Gc,1),0);
-                        obj.c  = varargin{1}.c;
-                        obj.Ac = zeros(0,size(obj.Gc,2));
-                        obj.Ab = [];
-                        obj.b  = [];
-                    case 'conZono'
-                        obj.Gc = varargin{1}.G;
-                        obj.Gb = zeros(size(obj.Gc,1),0);
-                        obj.c  = varargin{1}.c;
-                        obj.Ac = varargin{1}.A;
-                        obj.Ab = zeros(size(obj.Ac,1),0);
-                        obj.b  = varargin{1}.b;
-                    case 'double'
-                        obj.Gc = zeros(length(varargin{1}),0);
-                        obj.Gb = zeros(length(varargin{1}),0);
-                        obj.c  = varargin{1};
-                        obj.Ac = [];
-                        obj.Ab = [];
-                        obj.b  = [];
-                    case 'cell'
-                        if size(varargin{1},2) == 1
-                            obj = hPoly2hybZono(varargin{1});
-                        else
-                            obj = vPoly2hybZono(varargin{1});
-                        end
-                end
-            else
-                error('Incorrect number of inputs.')
-            end
-            % Dimension compatibility checking
-            try [obj.c obj.Gc obj.Gb];
-            catch error('Center (c) and generator matrices (Gc and Gb) must have the same number of rows.')
-            end
-            try [obj.Gc; obj.Ac];
-            catch error('Continuous generator matrix (Gc) and constraint matrix (Ac) must have the same number of columns.')
-            end
-            try [obj.Gb; obj.Ab];
-            catch error('Binary generator matrix (Gb) and constraint matrix (Ab) must have the same number of columns.')
-            end
-            try [obj.Ac obj.Ab obj.b];
-            catch error('Constraint matrices (Ac and Ab) and vector (b) must have the same number of rows.')
-            end
-        end
-        
-        % Dependent property getters
-        % Dimension of the space
-        function dimension = get.n(obj)
-            dimension = size(obj.c,1);
-        end
-        % Number of continuous generators
-        function nContinuousGenerators = get.nGc(obj)
-            nContinuousGenerators = size(obj.Gc,2);
-        end
-        % Number of binary generators
-        function nBinaryGenerators = get.nGb(obj)
-            nBinaryGenerators = size(obj.Gb,2);
-        end
-        % Number of constraints
-        function nConstraints = get.nC(obj)
-            nConstraints = size(obj.Ac,1);
-        end
-
-        % Hybrid zonotope-specific methods
-        [v,f] = plotAsConZono(obj,optSolver);   % Plot as the union of constrained zonotopes
-        [v,f] = plotHybZono1D(obj,optSolver);   % Plot in 1 dimension
-        [v,f] = plotHybZono2D(obj,optSolver);   % Plot in 2 dimensions
-        [v,f] = plotHybZono3D(obj,optSolver);   % Plot in 3 dimensions
-        [leaves] = getLeaves(obj,optSolver)     % Identify non-empty combinations of binary factors
-        
-         % Export each leaf as a cell array of conZonos 
-        function out = hyb2leafArray(obj)
-            leaves = obj.getLeaves();
-            for leaf = 1: length(leaves)
-                out{leaf,1} = conZono(obj.Gc,obj.c+obj.Gb*leaves(:,leaf),obj.Ac,obj.b-obj.Ab*leaves(:,leaf));
-            end
-        end
-    end
-end
\ No newline at end of file
diff --git a/@memZono/@hybZono/plotAsConZono.m b/@memZono/@hybZono/plotAsConZono.m
deleted file mode 100644
index 1be39f8..0000000
--- a/@memZono/@hybZono/plotAsConZono.m
+++ /dev/null
@@ -1,64 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Method:
-%       Return vertices and faces for a hybrid zonotope in 1D, 2D, or 3D by
-%       treating hybrid zonotope at the union of constrained zonotopes
-%   Syntax:
-%       [v,f] = plotAsConZono(Z,optSolver)
-%   Inputs:
-%       Z - 1D, 2D, or 3D hybrid zonotope in HCG-Rep (hybZono object)
-%       optSolver - solver options needed for linear and mixed-integer linear propgrams
-%   Outputs:
-%       v - nV x 3 matrix, each row denoting the x (first column), y (second column), 
-%                          and z (third column) positions of the nV vertices
-%       f - nF x nMax matrix, each row denoting the vertices (up to nMax) contained
-%                          in the nF faces (padded with NaN if face
-%                          contains less than nMax vertices)
-%   Notes:
-%       Not intended to be called directly by user.
-%       Use [v,f] = plot(obj,varargin) instead (method of abstractZono)
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-function [v,f] = plotAsConZono(obj,optSolver)
-
-% Determine number of non-empty constrained zonotopes
-[leaves] = getLeaves(obj,optSolver);
-if isempty(leaves)
-    warning('zonoLAB:EmptyZonotope','Hybrid zonotope is empty and cannot be plotted.')
-    v = []; f = [];
-    return
-end
-nLeaves = size(leaves,2);
-
-optPlot = plotOptions('Display','off','SolverOpts',optSolver);
-v = [];
-f = [];
-nVerts = zeros(nLeaves,1);
-waitbarHandle = waitbar(0,['Plotting hybrid zonotope with ',num2str(nLeaves),' leaves.']);
-%  If \xib denotes the i^th column of leaves, then the corresponding
-%  non-empty constrained zonotope is of the form:
-%  Z = { (c + Gb \xib) + Gc \xic | ||\xic||_inf <= 1, Ac \xi = b - Ab \xib }
-emptyLeaves = false;
-for i = 1:nLeaves
-    Zi = conZono(obj.Gc,obj.c+obj.Gb*leaves(:,i),obj.Ac,obj.b-obj.Ab*leaves(:,i));
-    [vi,fi] = plot(Zi,optPlot);
-    if size(vi,1) == 0
-        emptyLeaves = true;
-    end
-    nVerts(i) = size(vi,1);
-    v = [v;vi];
-    if size(fi,2) > size(f,2)
-        f = [f nan(size(f,1),size(fi,2)-size(f,2))]; 
-    end
-    if size(fi,2) < size(f,2)
-        fi = [fi nan(size(fi,1),size(f,2)-size(fi,2))];
-    end
-    f = [f;fi+sum(nVerts(1:i-1))];
-    waitbar(i/nLeaves,waitbarHandle)
-end
-close(waitbarHandle)
-
-% Check for unplotted conZono leaves
-if emptyLeaves
-    warning('zonoLAB:Tolerance','Some leaves of the hybrid zonotope did not plot and may be caused by constraints that are nearly redundant (close to the plotting/optimization tolerances). Check the validity of other leaves.')
-end
-
-end
\ No newline at end of file
diff --git a/@memZono/@hybZono/plotHybZono1D.m b/@memZono/@hybZono/plotHybZono1D.m
deleted file mode 100644
index d52e0ed..0000000
--- a/@memZono/@hybZono/plotHybZono1D.m
+++ /dev/null
@@ -1,24 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Method:
-%       Return vertices and faces for a hybrid zonotope in 1D
-%   Syntax:
-%       [v,f] = plotHybZono1D(Z,optSolver)
-%   Inputs:
-%       Z - 1D hybrid zonotope in HCG-Rep (hybZono object)
-%       optSolver - solver options needed for linear and mixed-integer linear propgrams
-%   Outputs:
-%       v - nV x 2 matrix, each row denoting the x (first column) and y (second column) positions
-%                          of the nV vertices
-%       f - nF x 2 matrix, each row denoting the two vertices contained
-%                          in the nF faces
-%   Notes:
-%       Not intended to be called directly by user.
-%       Use [v,f] = plot(obj,varargin) instead (method of abstractZono)
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-function [v,f] = plotHybZono1D(obj,optSolver)
-
-% Identify vertices and faces by treating hybrid zonotope at the union of
-% constrained zonotopes
-[v,f] = plotAsConZono(obj,optSolver);
-
-end
\ No newline at end of file
diff --git a/@memZono/@hybZono/plotHybZono2D.m b/@memZono/@hybZono/plotHybZono2D.m
deleted file mode 100644
index 575772e..0000000
--- a/@memZono/@hybZono/plotHybZono2D.m
+++ /dev/null
@@ -1,25 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Method:
-%       Return vertices and faces for a hybrid zonotope in 2D
-%   Syntax:
-%       [v,f] = plotHybZono2D(Z,optSolver)
-%   Inputs:
-%       Z - 2D hybrid zonotope in HCG-Rep (hybZono object)
-%       optSolver - solver options needed for linear and mixed-integer linear propgrams
-%   Outputs:
-%       v - nV x 2 matrix, each row denoting the x (first column) and y (second column) positions
-%                          of the nV vertices
-%       f - nF x nMax matrix, each row denoting the vertices (up to nMax) contained
-%                          in the nF faces (padded with NaN if face
-%                          contains less than nMax vertices)
-%   Notes:
-%       Not intended to be called directly by user.
-%       Use [v,f] = plot(obj,varargin) instead (method of abstractZono)
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-function [v,f] = plotHybZono2D(obj,optSolver)
-
-% Identify vertices and faces by treating hybrid zonotope at the union of
-% constrained zonotopes
-[v,f] = plotAsConZono(obj,optSolver);
-
-end
\ No newline at end of file
diff --git a/@memZono/@hybZono/plotHybZono3D.m b/@memZono/@hybZono/plotHybZono3D.m
deleted file mode 100644
index 663cc7e..0000000
--- a/@memZono/@hybZono/plotHybZono3D.m
+++ /dev/null
@@ -1,25 +0,0 @@
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-%   Method:
-%       Return vertices and faces for a hybrid zonotope in 3D
-%   Syntax:
-%       [v,f] = plotHybZono3D(Z,optSolver)
-%   Inputs:
-%       Z - 3D hybrid zonotope in HCG-Rep (hybZono object)
-%       optSolver - solver options needed for linear and mixed-integer linear propgrams
-%   Outputs:
-%       v - nV x 3 matrix, each row denoting the x (first column), y (second column), 
-%                          and z (third column) positions of the nV vertices
-%       f - nF x nMax matrix, each row denoting the vertices (up to nMax) contained
-%                          in the nF faces (padded with NaN if face
-%                          contains less than nMax vertices)
-%   Notes:
-%       Not intended to be called directly by user.
-%       Use [v,f] = plot(obj,varargin) instead (method of abstractZono)
-% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
-function [v,f] = plotHybZono3D(obj,optSolver)
-
-% Identify vertices and faces by treating hybrid zonotope at the union of
-% constrained zonotopes
-[v,f] = plotAsConZono(obj,optSolver);
-
-end
\ No newline at end of file
diff --git a/@memZono/memZono.m b/@memZono/memZono.m
index 66f666b..1daeabe 100644
--- a/@memZono/memZono.m
+++ b/@memZono/memZono.m
@@ -205,11 +205,17 @@
         % Plotting -------------------------
         varargout = plot(obj, dims, varargin);
         % Projection & Export ------------------
-        out = projection(obj,dims); % projection according to dims
+        out = projection(obj,dims,options); % projection according to dims
         function out = Z(obj,dims), out = projection(obj,dims).Z_; end
+        function out = cartProdProj(obj1,obj2,outDims)
+            out = projection(cartProd(obj1,obj2),outDims,"removeExtraFactors",true);
+        end
+        % function out = addCons(obj1,obj2) %<== adds cons from 2nd object w/out changing dims of 1st one
+        %     out = cartProdProj(obj1,obj2,obj1.dimKeys);
+        % end
         % Additional methods (implimentations of abstractZono methods in memZono)
         [NN,Y] = reluNN(X,Ws,bs,a);
-        [s,x_out] = supportFunc(obj,dims,d_in)
+        [s,x_out] = supportFunc(obj,dims,d_in);
     end
 
     %% Indexing  ----------------------------
@@ -231,6 +237,9 @@
         % &, and() - overide memoryIntersection w/ checks
         function out = and(obj1,obj2,sharedDimLabels)
             if nargin == 2, error('Must Supply labels for shared dimensions'); end
+            if ~isa(obj2,'memZono')
+                obj2 = memZono(obj2,sharedDimLabels); %<== attempt to call memZono with labels
+            end
             out = memoryIntersection(obj1,obj2,sharedDimLabels);
         end
         % Overide for memoryUnion w/ checks ... NOT IMPLIMENTED
@@ -254,7 +263,7 @@
         end
         % all (extended and)
         function out = all(in,ptn)
-            out = in{1}; if isemtpy(ptn), ptn = 'all'; end
+            out = in{1}; if isempty(ptn), ptn = 'all'; end
             for i = 2:numel(in), out = and(out,in{i},sprintf('_%s_%d',ptn,i)); end %<========= not efficient
         end
         % any (extended or) ... NOT IMPLIMENTED
diff --git a/@memZono/memoryIntersection.m b/@memZono/memoryIntersection.m
index 1af8d59..29e892a 100644
--- a/@memZono/memoryIntersection.m
+++ b/@memZono/memoryIntersection.m
@@ -101,128 +101,4 @@
 
     %% Define memZono
     obj = memZono(G_,c_,A_,b_,vset_,keys_);
-
-
-
-
-
-
-
-
-
-
-
-    % % Input Conditioning
-    % if ~isa(obj1,'memZono') || ~isa(obj2,'memZono')
-    %     error('Inputs must both be memZono objects')
-    % end
-    
-    % %% Keys
-    % % shared factors
-    % [k1,ks,k2] = memZono.getUniqueKeys(obj1.factorKeys,obj2.factorKeys);
-    % [idxk1,idxks1,idxks2,idxk2] = memZono.getKeyIndices(obj1.factorKeys,obj2.factorKeys);
-    % % shared dims
-    % [d1,ds,d2] = memZono.getUniqueKeys(obj1.dimKeys,obj2.dimKeys);
-    % [idxd1,idxds1,idxds2,idxd2] = memZono.getKeyIndices(obj1.dimKeys,obj2.dimKeys);
-    % % shared cons
-    % [c1,cs,c2] = memZono.getUniqueKeys(obj1.conKeys,obj2.conKeys);
-    % [idxc1,idxcs1,idxcs2,idxc2] = memZono.getKeyIndices(obj1.conKeys,obj2.conKeys);
-
-    % %% Factor-based Memory Cartisian Product
-    % G_ = [
-    %     obj1.G_(idxd1,idxk1), obj1.G_(idxd1,idxks1), zeros(length(d1),length(k2));
-    %     zeros(length(d2),length(k1)), obj2.G_(idxd2,idxks2), obj2.G_(idxd2,idxk2)
-    %     ];
-    % c_ = [
-    %     obj1.c_(idxd1,:);
-    %     obj2.c_(idxd2,:)
-    %     ];
-    % A_ = [
-    %     obj1.A_(idxc1,idxk1), obj1.A_(idxc1,idxks1), zeros(length(c1),length(k2));
-    %     zeros(length(c2),length(k1)), obj2.A_(idxc2,idxks2), obj2.A_(idxc2,idxk2);
-    %     ];
-    % b_ = [
-    %     obj1.b_(idxc1,:);
-    %     obj2.b_(idxc2,:);
-    %     ];
-
-    % % hybrid Zono
-    % if obj1.vset_(idxks1) ~= obj2.vset_(idxks2)
-    %     error('c/d factors not lining up');
-    % end
-    % vset_ = [obj1.vset_(idxk1),obj1.vset_(idxks1),obj2.vset_(idxk2)];
-
-    % % Labeling
-    % keys_.factors = [k1,ks,k2];
-    % keys_.dims = [d1,d2];
-    % keys_.cons = [c1,c2];
-
-    % %% Shared Dimensions
-    % if ~isempty(ds)
-    %     % Intersection matrix definition
-    %     R = zeros(length(ds),obj1.n);
-    %     for k = 1:length(ds)
-    %         j = idxds1(k); 
-    %         R(k,j) = 1; 
-    %     end
-        
-    %     % Interesecting Matrices
-    %     G_ = [G_;
-    %         obj1.G_(idxds1,idxk1), obj1.G_(idxds1,idxks1), zeros(length(ds),length(k2));
-    %     ];
-    %     c_ = [c_;
-    %         obj1.c_(idxds1,:);
-    %     ];
-    %     A_ = [A_; 
-    %         R*obj1.G_(:,idxk1), R*obj1.G_(:,idxks1)-obj2.G_(idxds2,idxks2), -obj2.G_(idxds2,idxk2) 
-    %     ];
-    %     b_ = [b_;
-    %         obj2.c_(idxds2,:) - R*obj1.c_;
-    %     ];
-    %     % Labels
-    %     keys_.dims = [keys_.dims, ds];
-
-    %     % Constraint Labels
-    %     if ~isa(sharedDimLabels,'cell')
-    %         if ~(isstring(sharedDimLabels)||ischar(sharedDimLabels))
-    %             error('Intersection operation will add additional constraints but labels for new constraints are not given.'); 
-    %         end
-    %         cds{length(ds)} = [];
-    %         for k = 1:length(ds)
-    %             % TODO: Warning if these new conKey labels already exist in either obj1 or obj2
-    %             cds{k} = sprintf('%s_%s_%d',sharedDimLabels,ds{k},k);
-    %         end
-    %     elseif length(sharedDimLabels) ~= length(ds)
-    %         error('Intersection operation will add additional constraints but labels for new constraints are not given (or not enough given).');
-    %     else
-    %         cds = sharedDimLabels;
-    %     end
-    %     if ~isempty([obj1.conKeys,obj2.conKeys])
-    %         if any(ismember(cds,[obj1.conKeys,obj2.conKeys]))
-    %             error('Intersection labels already exists... provide new names')
-    %         end
-    %     end
-    %     keys_.cons = [keys_.cons, cds];
-    % end
-
-    % %% Shared Constraints
-    % if ~isempty(cs)
-    %     if all([isnumeric(obj1.A_),isnumeric(obj2.A_),isnumeric(obj1.b_),isnumeric(obj2.b_)])
-    %         if all(obj1.A_(idxcs1,idxks1) ~= obj2.A_(idxcs2,idxks2),'all') || all(obj1.b_(idxcs1) ~= obj2.b_(idxcs2),'all')
-    %                 error('Shared Constraints are not identical')
-    %         end
-    %     end
-    %     A_ = [A_;
-    %         zeros(length(cs),length(k1)), obj1.A_(idxcs1,idxks1), zeros(length(cs),length(k2))
-    %     ];
-    %     b_ = [b_;
-    %         obj1.b_(idxcs1,:)
-    %     ];
-    %     % Labeling
-    %     keys_.cons = [keys_.cons,cs];
-    % end
-
-    % %% Define memZono
-    % obj = memZono(G_,c_,A_,b_,vset_,keys_);
-    
 end
diff --git a/@memZono/projection.m b/@memZono/projection.m
index 03468da..98e4920 100644
--- a/@memZono/projection.m
+++ b/@memZono/projection.m
@@ -1,5 +1,10 @@
 % Projection is defined for internal use - subsref (indexing) is simpilar syntax
-function out = projection(obj,dims)
+function out = projection(obj,dims,options)
+    arguments
+        obj
+        dims
+        options.removeExtraFactors = false;
+    end
     if ~iscell(dims) % if not already in cell form
         if strcmp(dims,':'), dims = obj.dimKeys;
         else, dims = obj.keysStartsWith(dims).dims;
@@ -7,5 +12,14 @@
     end
     [~,idx] = ismember(dims,obj.dimKeys);
     keys_out = obj.keys_; keys_out.dims = dims;
-    out = memZono(obj.G_(idx,:),obj.c_(idx,:),obj.A_,obj.b_,obj.vset_,keys_out);
+    out = memZono(obj.G_(idx,:),obj.c_(idx,:),obj.A_,obj.b_,obj.vset_,keys_out);   
+
+    if options.removeExtraFactors
+        idx = any(vertcat(any(obj.G_),any(obj.A_))); %<= indices where anything is nonzero
+        if ~all(idx)
+            keys_out.factors = keys_out.factors(idx);
+            out = memZono(out.G_(:,idx),out.c_,out.A_(:,idx),out.b_,out.vset_(idx),keys_out);
+        end
+    end
+
 end
\ No newline at end of file