fix: keep each subsystem with its reaction when readYAMLmodel flattens them

src/io/readYAMLmodel.m

@@ -813,9 +813,12 @@
 % end
 
 if isfield(model,'subSystems')
+% A reaction without subsystem holds [] or '' here: give it an empty entry,
+% so every reaction keeps its place when the entries are flattened
+    model.subSystems(cellfun(@isempty, model.subSystems)) = {{''}};
 % If all entries are 1x1, then flatten
     if all(cellfun(@(x) numel(x) <= 1, model.subSystems))
-        model.subSystems = transpose([model.subSystems{:}]);
+        model.subSystems = cellfun(@(x) x{1}, model.subSystems, 'UniformOutput', false);
     end
 end
 

tests/function_tests/tIO.m

@@ -207,6 +207,77 @@ function writeReadYAMLRoundTrip(testCase)
             testCase.verifyEqual(numel(m2.rxns), numel(testCase.model.rxns));
         end
 
+        function readYAMLmodelKeepsSubsystemsAligned(testCase)
+            % Reactions without a subsystem sit between, before and after
+            % those with one: each subsystem stays with its own reaction,
+            % and the result is a column of chars that checkModelStruct
+            % accepts.
+            m = testCase.model;
+            n = numel(m.rxns);
+            m.subSystems = repmat({{''}}, n, 1);
+            withSubsystem = [2 5 n-3];
+            names = {'Alpha'; 'Beta'; 'Gamma'};
+            for i = 1:numel(withSubsystem)
+                m.subSystems{withSubsystem(i)} = names(i);
+            end
+            f = [tempname '.yml'];
+            testCase.addTeardown(@() delete(f));
+            evalc('writeYAMLmodel(m, f);');
+            evalc('m2 = readYAMLmodel(f);');
+            testCase.verifyTrue(iscolumn(m2.subSystems));
+            testCase.verifyEqual(numel(m2.subSystems), n);
+            testCase.verifyTrue(iscellstr(m2.subSystems));
+            expected = repmat({''}, n, 1);
+            expected(withSubsystem) = names;
+            testCase.verifyEqual(string(m2.subSystems), string(expected));
+            ws = warning('off', 'RAVEN:modelError');
+            testCase.addTeardown(@() warning(ws));
+            try
+                evalc('checkModelStruct(m2, ''throwErrors'', true);');
+            catch err
+                testCase.verifyFail(['checkModelStruct refused the model: ' err.message]);
+            end
+        end
+
+        function readYAMLmodelKeepsNestedSubsystemsAligned(testCase)
+            % Reactions with several subsystems sit next to reactions with
+            % one and with none: every entry is a cell of chars, each stays
+            % with its own reaction, and checkModelStruct accepts the model.
+            m = testCase.model;
+            n = numel(m.rxns);
+            m.subSystems = repmat({{''}}, n, 1);
+            nested = {3, {'Alpha'; 'Delta'}; ...
+                      6, {'Beta'}; ...
+                      n-2, {'Gamma'; 'Epsilon'; 'Zeta'}};
+            for i = 1:size(nested,1)
+                m.subSystems{nested{i,1}} = nested{i,2};
+            end
+            f = [tempname '.yml'];
+            testCase.addTeardown(@() delete(f));
+            evalc('writeYAMLmodel(m, f);');
+            evalc('m2 = readYAMLmodel(f);');
+            testCase.verifyTrue(iscolumn(m2.subSystems));
+            testCase.verifyEqual(numel(m2.subSystems), n);
+            testCase.verifyTrue(all(cellfun(@iscell, m2.subSystems)));
+            for i = 1:n
+                row = find([nested{:,1}] == i);
+                if isempty(row)
+                    testCase.verifyTrue(isempty([m2.subSystems{i}{:}]), ...
+                        sprintf('reaction %d has no subsystem', i));
+                else
+                    testCase.verifyEqual(string(m2.subSystems{i}), ...
+                        string(nested{row,2}), sprintf('reaction %d', i));
+                end
+            end
+            ws = warning('off', 'RAVEN:modelError');
+            testCase.addTeardown(@() warning(ws));
+            try
+                evalc('checkModelStruct(m2, ''throwErrors'', true);');
+            catch err
+                testCase.verifyFail(['checkModelStruct refused the model: ' err.message]);
+            end
+        end
+
         function writeYAMLmodelQuotesYAML11Booleans(testCase)
             % A YAML 1.1 loader such as PyYAML reads a plain y, n, yes, no,
             % on or off in any case as a boolean, so these values are