docs: def indirectInitialMatrix(self, initialState):
          """
          Given some initial state, this iteratively determines new states.
          We repeatedly call the transition function on unvisited states in the frontier set.
          Each newly visited state is put in a dictionary called 'mapping' and the rates are stored in a dictionary.
          """
          mapping                 = {}   
          rates                   = OrderedDict()
          
          #Check whether the initial state is defined and of the correct type, and convert to a tuple or int. 
          convertedState           = self.checkInitialState(initialState)        
  
          if isinstance(convertedState,set):
              #If initialstates is a set, include all states in the set in the mapping.
              frontier = set( convertedState )
              for idx,state in enumerate(convertedState):  
                  mapping[state] = idx
                  if idx == 0: #Test the return type of the transition function (dict or numpy).
                      usesNumpy = self.checkTransitionType(initialState)                
          else:
              #Otherwise include only the single state.
              frontier                = set( [convertedState] )
              usesNumpy               = self.checkTransitionType(initialState)
              mapping[convertedState] = 0       
       
  
          while len(frontier) > 0:
              fromstate = frontier.pop()
              fromindex = mapping[fromstate]    
              
              if usesNumpy: #If numpy is used, convert to a dictionary with tuples and rates.
                  transitions = self.transition(np.array(fromstate))
                  transitions = self.convertToTransitionDict(transitions) 
              else:
                  transitions = self.transition(fromstate)
                  
              for tostate,rate in transitions.items():
                  if tostate not in mapping:
                      frontier.add(tostate)
                      mapping[tostate] = len(mapping)
                  toindex                     = mapping[tostate]
                  rates[(fromindex, toindex)] = rate
  
          #Inverse the keys and values in mapping to get a dictionary with indices and states.
          self.mapping = {value: key for key, value in list(mapping.items())}
          
          #Use the `rates` dictionary to fill a sparse dok matrix.
          D = dok_matrix((self.size,self.size))
          D.update(rates)
          return D.tocsr()