classdef testPolytopeConversion < TestCase
    properties
        simpleVrep
        simpleHrep
        simpleVolume
        flatHrep
        flatVrep
        flatVolume
        redundantHrep
        redundantVrep
        redundantVolume
        unboundedHrep
        unboundedInfSize
        unboundedVrep
        backAndForthVrep
    end
    
    methods
        %% Setup functions
        function self = testPolytopeConversion(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]';
            self.simpleVolume = 1;
            
            % flat example: line in 3D
            self.flatVrep.V = [-2 0 0; 1  0 0];
            self.flatHrep.A = [1 0 0; -1 0 0;0 1 0;0 -1 0;0 0 1;0 0 -1];
            self.flatHrep.h = [1 2 0 0 0 0]';
            self.flatVolume = @(flatness)(3*flatness^2);
            
            % redundant H-rep: square in 2D
            self.redundantHrep.A = [1 0;-1 0;0 1;0 -1; 1 1];
            self.redundantHrep.h = [1 1 1 1 4]';
            self.redundantVrep.V = [-1 -1;-1 1;1 -1;1 1];
            self.redundantVolume = 1;
            
            % unbounded H-rep: first quadrant in 2D
            self.unboundedHrep.A = [-1, 0; 0, -1];
            self.unboundedHrep.h = [0 0]';
            self.unboundedInfSize = 100;
            self.unboundedVrep.V = [0 0; 0 1; 1 0; 1 1] * self.unboundedInfSize;
            
            % back and forth sphere
            self.backAndForthVrep = obtainVrep('bipyramid', 15);
            
        end
        
        function tearDown(self)
            
        end
        
        %% test identifyPolytope function
        function testIdentifyVrep(self)
            info = identifyPolytope(self.simpleVrep);
            assertEqual(info.rep, 'Vrep', ...
                'representation should be identified to ''Vrep''');
            assertEqual(info.dim, 3, ...
                'dimension should be identified to 3');
            assertEqual(info.m, 8, ...
                'number of extreme points should be identified to 8');
        end
        
        function testIdentifyHrep(self)
            info = identifyPolytope(self.simpleHrep);
            assertEqual(info.rep, 'Hrep', ...
                'representation should be identified to ''Hrep''');
            assertEqual(info.dim, 3, ...
                'dimension should be identified to 3');
            assertEqual(info.m, 6, ...
                'number of facets should be identified to 6');
        end
        
        function testIdentifyError(self)
            assertExceptionThrown(@()(identifyPolytope(@(a)(a))), ...
                'elk:polytope:wrongType', 'Anonymous function as input should throw an error')
        end
        
        function testIdentifyHrepErrors(self)
            Hrep = self.simpleHrep;
            Hrep.h = Hrep.h';
            assertExceptionThrown(@()(identifyPolytope(Hrep)), ...
                'elk:polytope:wrongType', 'Transposed h-vector should throw an error.')
            Hrep.h = Hrep.h';
            Hrep.h(end) = [];
            assertExceptionThrown(@()(identifyPolytope(Hrep)), ...
                'elk:polytope:wrongType', 'Dimension mismatch should throw an error.')
        end
        
        %% Self convertions
        function testConvertVrepToVrep(self)
            PtOut = convertPolytope(self.simpleVrep, 'Vrep');
            info = identifyPolytope(PtOut);
            assertEqual(info.rep, 'Vrep', ...
                'representation should be converted to ''Vrep''');
        end
        
        function testConvertHrepToHrep(self)
            PtOut = convertPolytope(self.simpleHrep, 'Hrep');
            info = identifyPolytope(PtOut);
            assertEqual(info.rep, 'Hrep', ...
                'representation should be converted to ''Hrep''');
            assertEqual(...
                eliminateRows(PtOut.A, self.simpleHrep.A, eps), [], ...
                'mismatch of normal vectors during conversion');
            assertEqual(...
                eliminateRows(PtOut.h, self.simpleHrep.h, eps), [], ...
                'mismatch of facet distances during conversion');
        end
        
        %% Hrep versus Vrep
        function testConvertHrepToVrep(self)
            warning('off','elk:polytope:method');
            for useCddlib = 0:1
                PtOut = convertPolytope(self.simpleHrep, 'Vrep', ...
                    'useCddlib', useCddlib);
                info = identifyPolytope(PtOut);
                assertEqual(info.rep, 'Vrep', ...
                    'representation should be converted to ''Vrep''');
                assertEqual(...
                    eliminateRows(PtOut.V, self.simpleVrep.V, eps), [], ...
                    'mismatch of vertices during conversion')
            end
            warning('on','elk:polytope:method');
        end
        
        function testConvertVrepToHrep(self)
            for useCddlib = 0:1
                PtOut = convertPolytope(self.simpleVrep, 'Hrep', ...
                    'useCddlib', useCddlib);
                info = identifyPolytope(PtOut);
                assertEqual(info.rep, 'Hrep', ...
                    'representation should be converted to ''Hrep''');
                assertEqual(...
                    eliminateRows(PtOut.A, self.simpleHrep.A, 3*eps), [], ...
                    'mismatch of normal vectors during conversion');
                assertEqual(...
                    eliminateRows(PtOut.h, self.simpleHrep.h, 2*eps), [], ...
                    'mismatch of facet distances during conversion');
            end
        end
        
        %% Redundant representation
        function testConvertHrepToVrepRedundant(self)
            warning('off','elk:polytope:method');
            for useCddlib = 0:1
                PtOut = convertPolytope(self.redundantHrep, 'Vrep', ...
                    'useCddlib', useCddlib);
                PtOut.V = sortrows(round(1e15*PtOut.V)*1e-15);
                info = identifyPolytope(PtOut);
                assertEqual(info.rep, 'Vrep', ...
                    'representation should be converted to ''Vrep''');
                assertElementsAlmostEqual(PtOut.V, self.redundantVrep.V, ...
                    'mismatch of vertices during conversion')
            end
            warning('on','elk:polytope:method');
        end
        
        %% Degenerated Vrep
        function testConvertVrepToHrepFlat(self)
            for useCddlib = 0:1
                PtOut = convertPolytope(self.flatVrep, 'Hrep', ...
                    'useCddlib', useCddlib);
                info = identifyPolytope(PtOut);
                assertEqual(info.rep, 'Hrep', ...
                    'representation should be converted to ''Hrep''');
                assertElementsAlmostEqual(PtOut.A, self.flatHrep.A, ...
                    'mismatch of normal vectors during conversion');
                assertElementsAlmostEqual(PtOut.h, self.flatHrep.h, ...
                    'mismatch of facet distances during conversion');
            end
        end
        
        %% Unbounded Hrep
        function testConvertHrepToVrepUnbounded(self)
            for useCddlib = 0:1
                PtOut = convertPolytope(self.unboundedHrep, 'Vrep', ...
                    'infinityBoxSize', self.unboundedInfSize, ...
                    'useCddlib', useCddlib);
                PtOut.V = sortrows(PtOut.V);

                info = identifyPolytope(PtOut);
                assertEqual(info.rep, 'Vrep', ...
                    'representation should be converted to ''Vrep''');
                assertElementsAlmostEqual(PtOut.V, self.unboundedVrep.V, ...
                    'mismatch of vertices during conversion')
            end
        end
        
        %% testing obtainHrep (contained conversion)
        function testObtainHrepSimple(self)
            Hrep = obtainHrep('cube');
        end
        
        %% flatness option
        function testConvertHrepToPolyFlatnessConvinience(self)
            Hrep = convertPolytope(self.flatVrep, 'Hrep', 'flatness', 1e-3);
            Vrep = convertPolytope(Hrep, 'Vrep');
            info = identifyPolytope(Vrep);
            assertEqual(info.rep, 'Vrep', ...
                'representation should be converted to ''Vrep''');
            assertEqual(info.m, 8, ...
                'flatness does not seem to be applied');
            [K, V] = convhulln(Vrep.V);
            assertElementsAlmostEqual(V, self.flatVolume(1e-3), ...
                'relative', eps, eps, ...
                'volume mismatch: flatness might be applied with a wrong value');
        end
        
        function testConvertHrepToPolyFlatnessManual(self)
            Hrep = convertVrepToHrep(self.flatVrep, 1e-10, 1e-1);
            Vrep = convertPolytope(Hrep, 'Vrep');
            info = identifyPolytope(Vrep);
            assertEqual(info.rep, 'Vrep', ...
                'representation should be converted to ''Vrep''');
            assertEqual(info.m, 8, ...
                'flatness does not seem to be applied');
            [K, V] = convhulln(Vrep.V);
            assertElementsAlmostEqual(V, self.flatVolume(1e-1), ...
                'relative', eps, eps, ...
                'volume mismatch: flatness might be applied with a wrong value');
        end
        
        %% Back and forth
        function testConvertBackAndForth(self)
            mOriginal = size(self.backAndForthVrep.V, 1);
            
            % round one
            stepOneHrep = convertPolytope(self.backAndForthVrep, 'Hrep');
            stepTwoVrep = convertPolytope(stepOneHrep, 'Vrep');
            
            assertEqual(mOriginal, size(stepTwoVrep.V, 1), ...
                'vertices were eliminated in the first round');
            assertEqual(eliminateRows(stepTwoVrep.V, self.backAndForthVrep.V, 500*eps), ...
                [], 'sphere back and forth conversion failed in first round');
            
            % round two
            stepThreeHrep = convertPolytope(stepTwoVrep, 'Hrep');
            stepFourVrep = convertPolytope(stepThreeHrep, 'Vrep');
            
            assertEqual(mOriginal, size(stepFourVrep.V, 1), ...
                'vertices were eliminated in the second round');
            assertEqual(eliminateRows(stepFourVrep.V, self.backAndForthVrep.V, 1e3*eps), ...
                [], 'sphere back and forth conversion failed in second round');
            
            % round three
            stepFiveHrep = convertPolytope(stepFourVrep, 'Hrep');
            stepSixVrep = convertPolytope(stepFiveHrep, 'Vrep');
            
            assertEqual(mOriginal, size(stepSixVrep.V, 1), ...
                'vertices were eliminated in the third round');
            assertEqual(eliminateRows(stepSixVrep.V, self.backAndForthVrep.V, 2e3*eps), ...
                [], 'sphere back and forth conversion failed in third round');
            
        end
        
        %% Debug examples
        function testConvertVrepDualProblem(self)
            vrepOriginal.V = [1 0;0 0;0 1];
            hrep = convertVrepToHrep(vrepOriginal);
            vrep = convertHrepToVrep(hrep);
            
            assertEqual(size(vrep.V, 1), 3, 'vertices were eliminated');
            
            assertEqual(eliminateRows(vrep.V, vrepOriginal.V, eps), ...
                [], 'conversion of vrep does not map back to itself');
        end
        
        % the following case is an example for which cddlib throws an error
        %   so that convertHrepToVrep throws a warning and switches to an
        %   alternative method (convhulln via dual polytope)
        % the second interesting part is that cddlib reduces the number of
        %   vertices in subsequent reductions. See details in code
        function testCddlibErrors(self)
            % load polytope
            load hrepCddlibError_1.mat
            
            % note that the reduction of the h-representation indeed has
            %   less facets.
            hrepReduced = reduceHrep(hrep);
            assertEqual(size(hrepReduced.A, 1), 16);
            
            % for the following conversion, a warning is thrown, because
            %   cddlib fails! We supress this temporarly.
            % This case failed in REV 340 because the backup alternative
            %   was not triggered properly, resulting in an error.
            % For failure zerotolCddlib=1e-10 is required.
%             s = warning('off', 'elk:polytope:conversionMethod');
            vrep_c1 = convertHrepToVrep(hrep);
%             warning(s);
            
            % conversion back to hrep gives the original reduced hrep
            hrep_c2 = convertVrepToHrep(vrep_c1);
            assertEqual(eliminateRows(hrep_c2.A, ...
                hrepReduced.A, 100*eps), [], ['Original hrep is different ' ...
                'from converted hrep']);
            
            % the next interesting part is that hrep_c2 is reduced in the
            %   number of vertices if converted to vrep! This is an error 
            %   in cddlib! Not using cddlib converts properly.
            % However, this only happens for zerotolCddlib < 1e-10, while
            %   rather large values are recommended by the cddlib 
            %   documentation.
            vrep_c3 = convertHrepToVrep(hrep_c2);
            vrep_c3_noCddlib = convertHrepToVrep(hrep_c2, elkZerotol, 0, 0);
            assertEqual(size(vrep_c1.V, 1), size(vrep_c3_noCddlib.V, 1), ...
                'Vrep does not match original vrep after double conversion');
            assertEqual(size(vrep_c1.V, 1), size(vrep_c3.V, 1), ...
                ['You''have selected a rather small zerotolCddlib, check ' ...
                 'cddlib manual for reference, 1e-9 should succeed here, ' ...
                 'even though it might fail in other cases.']);
        end
    end
end
