classdef testPolytopeOperation < TestCase
    properties
        simpleVrep
        simpleHrep
    end
    
    methods
        %% Setup functions
        function self = testPolytopeOperation(name)
            self = self@TestCase(name);
        end
        
        function setUp(self)
            % simple example: cube
            self.simpleVrep = obtainVrep('cube');
            self.simpleHrep.A = [1 0 0; -1  0  0; ...
                          0 1 0;  0 -1  0; ...
                          0 0 1;  0  0 -1];
            self.simpleHrep.h = [1 0 1 0 1 0]';
        end
        
        function tearDown(self)
            
        end
        
        %% test move functions
        function testMoveVrep(self)
            Vrep = moveVrep(self.simpleVrep, [1 1 1]);
            assertElementsAlmostEqual(Vrep.V, self.simpleVrep.V + 1, ...
                'relative', eps, eps, 'moving Vrep failed');
        end
        
        function testMoveHrep(self)
            v = [1 2 3]';
            Hrep = moveHrep(self.simpleHrep, v);
            assertElementsAlmostEqual(Hrep.h, self.simpleHrep.h + self.simpleHrep.A*v, ...
                'relative', eps, eps, 'moving Hrep failed');
            assertEqual(Hrep.A, self.simpleHrep.A, 'moving Hrep failed');
        end
        
        function testMovePolytope(self)
            v = [1 2 3]';
            Pt = movePolytope(self.simpleHrep, v);
            assertElementsAlmostEqual(Pt.h, self.simpleHrep.h + self.simpleHrep.A*v, ...
                'relative', eps, eps, 'moving Hrep failed');
            assertEqual(Pt.A, self.simpleHrep.A, 'moving polytope failed');
        end
        
        %% test normalizeHrep
        function testNormalizeHrep(self)
            Hrep = self.simpleHrep;
            Hrep.A = 2*Hrep.A;
            Hrep = normalizeHrep(Hrep);
            assertElementsAlmostEqual(Hrep.h, 0.5*self.simpleHrep.h, ...
                'relative', eps, eps, 'normalizing Hrep failed');
            assertEqual(Hrep.A, self.simpleHrep.A, 'normalizing Hrep failed');
        end
        
        %% test stretch functions
        function testStretchVrep(self)
            scale = [1 2 3];
            Vrep = stretchVrep(self.simpleVrep, scale);
            minMax = [min(Vrep.V(:,1)) max(Vrep.V(:,1)); ...
                      min(Vrep.V(:,2)) max(Vrep.V(:,2)); ...
                      min(Vrep.V(:,3)) max(Vrep.V(:,3))];
            assertElementsAlmostEqual(minMax, ...
                [0 1;0 2; 0 3], 'relative', eps, eps, 'stretching Vrep failed');
        end
        
        function testStretchVrepFixpoint(self)
            scale = [1 2 3];
            x0 = [1 1 1]*0.5;
            Vrep = stretchVrep(self.simpleVrep, scale, x0);
            minMax = [min(Vrep.V(:,1)) max(Vrep.V(:,1)); ...
                      min(Vrep.V(:,2)) max(Vrep.V(:,2)); ...
                      min(Vrep.V(:,3)) max(Vrep.V(:,3))];
            assertElementsAlmostEqual(minMax, ...
                [0 1;-0.5 1.5;-1 2], 'relative', eps, eps, 'stretching Vrep with fisxpoint failed');
        end
        
        function testStretchHrep(self)
            scale = [1 2 3];
            Hrep = stretchHrep(self.simpleHrep, scale);
            Hrep = normalizeHrep(Hrep);
            assertElementsAlmostEqual(Hrep.h, [1 0 2 0 3 0]', ...
                'relative', eps, eps, 'stretching Hrep failed');
            assertEqual(Hrep.A, self.simpleHrep.A, 'stretching Hrep failed');
        end
        
        function testStretchHrepFixpoint(self)
            scale = [1 2 3];
            x0 = [1 1 1]*0.5;
            Hrep = stretchHrep(self.simpleHrep, scale, x0);
            Hrep = normalizeHrep(Hrep);
            assertElementsAlmostEqual(Hrep.h, [1 0 1.5 0.5 2 1]', ...
                'relative', eps, eps, 'stretching Hrep with fisxpoint failed');
            assertEqual(Hrep.A, self.simpleHrep.A, 'stretching Hrep with fisxpoint failed');
        end
        
        function testStretchPolytope(self)
            scale = [1 2 3]';
            Pt = stretchPolytope(self.simpleHrep, scale);
            Pt = normalizeHrep(Pt);
            assertElementsAlmostEqual(Pt.h, [1 0 2 0 3 0]', ...
                'relative', eps, eps, 'stretching Polytope failed');
            assertEqual(Pt.A, self.simpleHrep.A, 'stretching Polytope failed');
        end
        
        %% rotatePolytope3d (Vrep, Hrep only)
        
        function testRotateHrep2d(self)
            % we try rotation from cube to diamond shape. Note that the
            % size of the object is not considered.
            HrepOriginal = obtainHrep('square', [], 2);
            HrepResult = obtainHrep('circle', 4, 2);
            
            Hrep = rotateHrep2d(HrepOriginal, pi/4, [0.5 0.5]);
            assertEqual(...
                eliminateRows(HrepResult.A, Hrep.A, 2*eps), [], ...
                'mismatch of normal vectors after rotation');
        end
        
        function testRotateHrep3d(self)
            % we try rotation from cube to diamond shape. Note that the
            % size of the object is not considered.
            HrepOriginal = obtainHrep('cube');
            HrepResult = obtainHrep('prism', 4);
            
            Hrep = rotateHrep3d(HrepOriginal, 'z', pi/4, [0 0 0]);
            assertEqual(...
                eliminateRows(HrepResult.A, Hrep.A, 2*eps), [], ...
                'mismatch of normal vectors after rotation');
        end
        
        function testRotateVrep2d(self)
            % we try rotation from cube to diamond shape.
            VrepOriginal = obtainVrep('square', [], 2);
            VrepResult = obtainVrep('circle', 4, 2);
            VrepResult.V = VrepResult.V ./ sqrt(2);
            VrepResult = moveVrep(VrepResult, [0.5, 0.5]);
            
            Vrep = rotateVrep2d(VrepOriginal, pi/4, [0.5 0.5]);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        function testRotateVrep3d(self)
            % we try rotation from cube to diamond shape.
            VrepOriginal = obtainVrep('cube');
            VrepResult = obtainVrep('prism', 4);
            VrepResult = stretchVrep(VrepResult, [1/sqrt(2) 1/sqrt(2) 1], [0 0 0]);
            VrepResult = moveVrep(VrepResult, [0.5, 0.5, 0]);
            
            Vrep = rotateVrep3d(VrepOriginal, 'z', pi/4, [0.5 0.5 0]);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        %% rotatePolytope (repetition from above)
        function testRotatePolyHrep2d(self)
            % we try rotation from cube to diamond shape. Note that the
            % size of the object is not considered.
            HrepOriginal = obtainHrep('square', [], 2);
            HrepResult = obtainHrep('circle', 4, 2);
            
            Hrep = rotatePolytope(HrepOriginal, '+', pi/4, 'center', [0.5 0.5]);
            assertEqual(...
                eliminateRows(HrepResult.A, Hrep.A, 2*eps), [], ...
                'mismatch of normal vectors after rotation');
        end
        
        function testRotatePolyHrep3d(self)
            % we try rotation from cube to diamond shape. Note that the
            % size of the object is not considered.
            HrepOriginal = obtainHrep('cube');
            HrepResult = obtainHrep('prism', 4);
            
            Hrep = rotatePolytope(HrepOriginal, 'z', pi/4);
            assertEqual(...
                eliminateRows(HrepResult.A, Hrep.A, 2*eps), [], ...
                'mismatch of normal vectors after rotation');
        end
        
        function testRotatePolyVrep2d(self)
            % we try rotation from cube to diamond shape.
            VrepOriginal = obtainVrep('square', [], 2);
            VrepResult = obtainVrep('circle', 4, 2);
            VrepResult.V = VrepResult.V ./ sqrt(2);
            VrepResult = moveVrep(VrepResult, [0.5, 0.5]);
            
            Vrep = rotatePolytope(VrepOriginal,  'center', [0.5 0.5], '+', pi/4);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        function testRotatePolyVrep3d(self)
            % we try rotation from cube to diamond shape.
            VrepOriginal = obtainVrep('cube');
            VrepResult = obtainVrep('prism', 4);
            VrepResult = stretchVrep(VrepResult, [1/sqrt(2) 1/sqrt(2) 1], [0 0 0]);
            VrepResult = moveVrep(VrepResult, [0.5, 0.5, 0]);
            
            Vrep = rotatePolytope(VrepOriginal, 'center', [0.5 0.5 0], 'z', pi/4);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        %% rotatePolytope
        % poly without fixed position, see above
        
        % poly with subsequent rotation
        function testRotatePolytopeMulti(self)
            % we rotate a cube into itself
            VrepOriginal = obtainVrep('cube');
            VrepResult = VrepOriginal;
                        
            Vrep = rotatePolytope(VrepOriginal, 'x', pi/2, 'y', -1*pi/2, 'z', pi);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        %% rotate to Vector
        function testRotateHrep3dVector(self)
            % we try rotation from cube to diamond shape. Note that the
            % size of the object is not considered.
            HrepOriginal = obtainHrep('cube');
            HrepResult = obtainHrep('prism', 4);
            
            HrepResult = rotateHrep3d(HrepResult, 'vector', [1 1 0]);
            Hrep = rotateHrep3d(HrepOriginal, 'vector', [1 0 0]);
            assertEqual(...
                eliminateRows(HrepResult.A, Hrep.A, 2*eps), [], ...
                'mismatch of normal vectors after rotation');
        end
        
        function testRotateVrep3dVector(self)
            % we compare manual rotation by automatic rotation
            VrepOriginal = obtainVrep('bipyramid', 7);
            VrepResult = rotatePolytope(VrepOriginal, 'y', -pi/4);
            
            Vrep = rotateVrep3d(VrepOriginal, 'vector', [1 0 1]);
            assertEqual(...
                eliminateRows(VrepResult.V, Vrep.V, 2*eps), [], ...
                'mismatch after rotation');
        end
        
        %% Scaling polytopes
        % The main focus of scaling is to ensure:
        %   (1) the origin is contained in the inner of the polytope
        %   (2) a reasonable sizing
        % The first condition is rather strict, while the second depends on
        %   requirements.
        function testScaleHrep(self)
            hrepLarge = moveHrep(self.simpleHrep, [100 0 0]);
            hrepLarge = stretchHrep(hrepLarge, 1e4);
            hrepLarge = normalizeHrep(hrepLarge);
            % scaling
            [hrepScaled scale x0]= scaleHrep(hrepLarge);
            % was origin outside
            assertFalse(all(hrepLarge.A * [0 0 0]' < hrepLarge.h), ...
                'Test polytope should have its origin outside');
            % is origin contained?
            assertTrue(all(hrepScaled.A * [0 0 0]' < hrepScaled.h), ...
                'Scaled polytope does not have its origin inside');
            % we know the origingal mean width (1.5), so the scaled should 
            %   be between 1 and 2
            meanWidthScaled = computeMeanWidthHrep(hrepScaled);
            assertTrue(1 < meanWidthScaled && meanWidthScaled > 2, ...
                'volume of scaled unit cube exceeds [1 2]');
            % check the scaling result by volume
            meanWidthLarge = computeMeanWidthHrep(hrepLarge);
            assertElementsAlmostEqual(meanWidthLarge, scale*meanWidthScaled, ...
                'relative', 10*eps, eps, 'Scaling for scaleHrep does not match (mean width)')
            % check x0 and scaling by manual reverse transformation
            hrepScaledBack = stretchHrep(hrepScaled, scale);
            hrepScaledBack = moveHrep(hrepScaledBack, x0);
            hrepScaledBack = normalizeHrep(hrepScaledBack);
            % rough check
            assertElementsAlmostEqual(...
                sortrows(hrepLarge.A), sortrows(hrepScaledBack.A), ...
                'relative', eps, eps, 'Manual reverse scaling does not match (1)');
            assertElementsAlmostEqual(...
                sortrows(hrepLarge.h), sortrows(hrepScaledBack.h), ...
                'relative', eps, 3e3*eps, 'Manual reverse scaling does not match (2)');
        end
        
        function testScaleVrep(self)
            vrepLarge = moveVrep(self.simpleVrep, [100 0 0]);
            vrepLarge = stretchVrep(vrepLarge, 1e-4);
            % scaling
            [vrepScaled scale x0]= scaleVrep(vrepLarge, [], 1);
            % is origin contained.. ..according to the used polytope the
            %   following test should suffice
            assertTrue(all(sum([vrepScaled.V;vrepScaled.V] > 0, 1)), ...
                'Scaled polytope does not have its origin inside');
            % we know the origingal mean width (1.5), so the scaled should 
            %   be between 1 and 2
            meanWidthScaled = computeMeanWidthVrep(vrepScaled);
            assertTrue(1 < meanWidthScaled && meanWidthScaled > 2, ...
                'volume of scaled unit cube exceeds [1 2]');
            % check the scaling result by volume
            meanWidthLarge = computeMeanWidthVrep(vrepLarge);
            assertElementsAlmostEqual(meanWidthLarge, scale*meanWidthScaled, ...
                'relative', 10*eps, eps, 'Scaling for scaleVrep does not match (mean width)')
            % check x0 and scaling by manual reverse transformation
            vrepScaledBack = stretchVrep(vrepScaled, scale);
            vrepScaledBack = moveVrep(vrepScaledBack, x0);
            % rough check
            assertElementsAlmostEqual(...
                sortrows(vrepLarge.V), sortrows(vrepScaledBack.V), ...
                'relative', eps, eps, 'Manual reverse scaling does not match (1)');
        end
        
        %% Dual polytope
        function testDualFromHrep(self)
            hrep = obtainHrep('prism', 4);
            hrep = rotateHrep3d(hrep, 'z', pi/4);
            hrep = moveHrep(hrep, [0 0 -0.5]);
            hrep = stretchHrep(hrep, [sqrt(2) sqrt(2) 2]);
            dualVrepResult = obtainVrep('bipyramid', 4);
            
            dualVrep = computeDualPolytope(hrep);
            dualVrepInfo = identifyPolytope(dualVrep);
            
            assertEqual(dualVrepInfo.rep, 'Vrep', ...
                'representation should be identified to ''Vrep''');
            
            assertEqual(...
                eliminateRows(dualVrepResult.V, dualVrep.V, 2*eps), [], ...
                'mismatch after dual operation');
        end
        
        function testDualFromVrep(self)
            vrep = obtainVrep('prism', 4);
            vrep = rotateVrep3d(vrep, 'z', pi/4);
            vrep = moveVrep(vrep, [0 0 -0.5]);
            vrep = stretchVrep(vrep, [sqrt(2) sqrt(2) 2]);
            dualHrepResult = obtainHrep('bipyramid', 4);
            
            dualHrep = computeDualPolytope(vrep);
            dualHrepInfo = identifyPolytope(dualHrep);
            
            assertEqual(dualHrepInfo.rep, 'Hrep', ...
                'representation should be identified to ''Hrep''');
            
            assertEqual(...
                eliminateRows(dualHrepResult.A, dualHrep.A, 2*eps), [], ...
                'mismatch after dual operation');
        end
        
        function testDualBackAndForth(self)
            vrepOriginal = obtainVrep('prism', 7);
            dualHrep = computeDualPolytope(vrepOriginal);
            vrep = computeDualPolytope(dualHrep);
            
            assertElementsAlmostEqual(vrepOriginal.V, vrep.V, ...
                'relative', eps, eps, ...
                'mismatch after dual operation');
        end
                    
        %% Minkowski addition (Vrep, Hrep)
        % zerotol and reduce options not tested here
        %   zerotol is called in the section below with
        %   addMinkowskiPolytope
        function testAddMinkowskiVrep(self)
           vrepOne = obtainVrep('line');
           vrepOne = rotateVrep3d(vrepOne, 'y', -pi/2);
           vrepTwo = obtainVrep('polygon', 5);
           vrepResult = obtainVrep('prism', 5);
           
           vrep = addMinkowskiVrep(vrepOne, vrepTwo);
           
           assertEqual(...
               eliminateRows(vrep.V, vrepResult.V, 2*eps), [], ...
               'mismatch after Minkowski addition');
            
        end
        
        function testAddMinkowskiHrep(self)
           hrepOne = obtainHrep('line');
           hrepOne = rotateHrep3d(hrepOne, 'y', -pi/2);
           hrepTwo = obtainHrep('polygon', 5);
           hrepResult = obtainHrep('prism', 5);

           hrep = addMinkowskiHrep(hrepOne, hrepTwo);
           
           assertEqual(...
               eliminateRows(hrep.A, hrepResult.A, 4*eps), [], ...
               'mismatch after Minkowski addition');
            
        end
        
        %% Minkowski addidion (poly)
        
        % multiple Vrep
        function testAddMinkowskiPolyMultipleVrep(self)
           vrepBaseline = obtainVrep('line', [], 2);
           vrepBaseline = moveVrep(vrepBaseline, [-0.5 0]);
           vrepOne = vrepBaseline;
           vrepTwo = rotateVrep2d(vrepBaseline, pi/2);
           % these both above result in a square
           vrepThree = rotateVrep2d(vrepBaseline, pi/4);
           vrepFour = rotateVrep2d(vrepBaseline, -pi/4);
           % with these two, we have a octagon
           vrepResult = obtainVrep('polygon', 8, 2);
           vrepResult.V = 1/sqrt(2-sqrt(2))*vrepResult.V;
           vrepResult = rotateVrep2d(vrepResult, pi/8);
           
           vrep = addMinkowskiPolytope('zerotol', 1e-8, ...
               vrepOne, vrepTwo, vrepThree, vrepFour);
           
           assertEqual(...
               eliminateRows(vrep.V, vrepResult.V, 3*eps), [], ...
               'mismatch after Minkowski addition');
        end
        
        % mixed Vrep/Hrep, first Hrep
        function testAddMinkowskiPolyMixVrepHrep(self)
           vrepBaseline = obtainVrep('line', [], 2);
           vrepBaseline = moveVrep(vrepBaseline, [-0.5 0]);
           hrepOne = convertVrepToHrep(vrepBaseline);
           vrepTwo = rotateVrep2d(vrepBaseline, pi/2);
           % these both above result in a square
           vrepThree = rotateVrep2d(vrepBaseline, pi/4);
           hrepFour = convertVrepToHrep(...
               rotateVrep2d(vrepBaseline, -pi/4));
           % with these two, we have a octagon
           hrepResult = obtainHrep('polygon', 8, 2);
           hrepResult = stretchHrep(hrepResult, [1 1]*1/sqrt(2-sqrt(2)));
           hrepResult = rotateHrep2d(hrepResult, pi/8);
           hrepResult = normalizeHrep(hrepResult);
           
           hrep = addMinkowskiPolytope('zerotol', 1e-8, ...
               hrepOne, vrepTwo, vrepThree, hrepFour);
           
           assertEqual(...
               eliminateRows(hrep.A, hrepResult.A, 4*eps), [], ...
               'mismatch after Minkowski addition');
        end
        
        %% Projection and Support set
        function testProjectionVrep(self)
            vrep = projectVrep3d(self.simpleVrep, [0 0 1]);
            vrepSquare = obtainVrep('square', [], 2);
            assertEqual(...
                eliminateRows(vrep.V, vrepSquare.V, 2*eps), [], ...
                'mismatch after projection');
        end
        
         function testSupportSetVrep(self)
            vrep = computeSupportSetVrep(self.simpleVrep, [0 0 1]);
            vrepSquare = obtainVrep('square', [], 2);
            vrep.V = sortrows(vrep.V);
            vrep = moveVrep(vrep, -1*vrep.V(1,:));
            assertEqual(...
                eliminateRows(vrep.V,  vrepSquare.V, 2*eps), [], ...
                'mismatch for support set');
            
            vrep = computeSupportSetVrep(self.simpleVrep, [0 1 1]);
            assertElementsAlmostEqual(vrep.V, [-0.5, 0.5]', ...
                'mismatch for support set');
            
            vrep = computeSupportSetVrep(self.simpleVrep, [1 1 1]);
            assertTrue(isempty(vrep.V), 'mismatch for support set');
         end
        
         %% reduceVrepDimension
         % is contained in computeSupportSet*. It is not in private to
         %   support triangulation (which might be in polytope-operation)
         
         %% reduceHrepDimension & liftHrepDimension
         % both functions are checked together
         
         function testReduceHrepNarrowFullDim(self)
             % this polytope should result in a full-dimensional polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); ...
                       1/sqrt(2) 1/sqrt(2) 0; ...
                       1/sqrt(2) -1/sqrt(2) 0];
             hrep.h = [1 1 1 1 1 1 1 1]'*9e-10;
             
             hrepRed = reduceHrepDimension(hrep, 1e-10);
             
             assertEqual(hrep.A, hrepRed.A, 'The result should exaclty be the input');
             assertEqual(hrep.h, hrepRed.h, 'The result should exaclty be the input');
         end
         
         function testReduceHrepEmpty(self)
             % this polytope should result in an empty polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); ...
                       1/sqrt(2) 1/sqrt(2) 0; ...
                       1/sqrt(2) -1/sqrt(2) 0];
             hrep.h = -1*[1 1 1 1 1 1 1 1]'*3e-10;
             
             hrepRed = reduceHrepDimension(hrep, 1e-10);
             
             assertTrue(isempty(hrepRed), 'This polytope should be empty.');
         end
         
         function testReduceHrepZeroDim(self)
             % this polytope should result in an empty polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); ...
                       1/sqrt(2) 1/sqrt(2) 0; ...
                       1/sqrt(2) -1/sqrt(2) 0];
             hrep.h = -1*[1 1 1 1 1 1 1 1]'*1e-11;
             hrep = rotateHrep3d(hrep, 'vector', [1 2 3]');
             hrep = movePolytope(hrep, [4 5 10]);
             
             [hrepRed rot x0] = reduceHrepDimension(hrep, 1e-10);
             
             assertTrue(isempty(hrepRed.A), 'This polytope should be empty.');
             assertTrue(isempty(hrepRed.h), 'This polytope should be empty.');
%              x0 - [4 5 10]'
             assertElementsAlmostEqual(x0, [4 5 10]', 'absolute', 5e-11, eps);
         end
         
         function testReduceHrepOneDimThreeBounds(self)
             % this polytope should result in an empty polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); ...
                       1/sqrt(2) 1/sqrt(2) 0; ...
                       1/sqrt(2) -1/sqrt(2) 0];
             hrep.h = [1 1 1e-2    1e-11 1 10    0 0]';
             % >> this should have a length of 10.01
             hrep = rotateHrep3d(hrep, 'vector', [1 2 3]');
             hrep = movePolytope(hrep, [4 5 10]);
             
             [hrepRed rot x0] = reduceHrepDimension(hrep, 1e-10);
             
             assertEqual(size(hrepRed.A, 2), 1, 'Result should be 1-dimensional');
             assertElementsAlmostEqual(...
                 hrepRed.h(1)+hrepRed.h(2), ...
                 10.01, 'absolute', 8*eps, eps);
         end
         
         function testReduceHrepOneDimTwoBounds(self)
             % this polytope should result in an empty polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); 1 1 0; 1 -1 0];
             hrep.h = [0 0 3 0 0 5 0 0]';
             % >> this should have a length of 8
             hrep = rotateHrep3d(hrep, 'vector', [1 2 3]');
             hrep = movePolytope(hrep, [4 5 10]);
             
             [hrepRed rot x0] = reduceHrepDimension(hrep, 1e-10);
             hrepRec = liftHrepDimension(hrep, rot, x0);
             
             assertEqual(size(hrepRed.A, 2), 1, 'Result should be 1-dimensional');
             assertElementsAlmostEqual(...
                 hrepRed.h(1)+hrepRed.h(2), ...
                 8, 'absolute', 8*eps, eps);
         end
         
         function testReduceHrepTwoDim(self)
             % this polytope should result in an empty polytope
             %   according to zerotol
             hrep.A = [eye(3); -1*eye(3); ...
                       1/sqrt(2) 1/sqrt(2) 0; ...
                       1/sqrt(2) -1/sqrt(2) 0];
             hrep.h = [1 2 0 4 5 0 0.5 0.25]';
             % we keep the orioginal which only lives in the xy-plane so
             %   that we can reduce dimensions by hand
             hrepRef.A = hrep.A(:, 1:2);
             hrepRef.h = hrep.h(any(abs(hrepRef.A) > 1e-10, 2));
             hrepRef.A = hrepRef.A(any(abs(hrepRef.A) > 1e-10, 2), :);
             
             % but the example shall be rotated and moved
             hrep = rotateHrep3d(hrep, 'vector', [1 2 3]');
             hrep = movePolytope(hrep, [4 5 10]);
             
             [hrepRed rot x0] = reduceHrepDimension(hrep, 1e-10);
             hrepRec = liftHrepDimension(hrep, rot, x0);
             
%              computeVolumeHrep(hrepRed) - computeVolumeHrep(hrepRef)
             assertEqual(size(hrepRed.A, 2), 2, 'Result should be 2-dimensional');
             assertElementsAlmostEqual(...
                 computeVolumeHrep(hrepRed), ...
                 computeVolumeHrep(hrepRef), 'absolute', 90*eps, eps);
         end
         
         
    end
end
