docs: def parse(self, handler):
          """
          Main method to parse the Sonata files and call the appropriate methods
          in the handler
          """
  
          ########################################################################
          #  load the main configuration scripts    
      
          main_config_filename = os.path.abspath(self.parameters['filename'])
          
          config = load_json(main_config_filename)
          
          self.init_substitutes = {'.':'%s/'%os.path.dirname(main_config_filename),
                                   '../':'%s/'%os.path.dirname(os.path.dirname(main_config_filename))}
          
          self.substitutes = {'${configdir}':'%s'%os.path.dirname(main_config_filename)}
          
          if 'network' in config:
              self.network_config = load_json(self.subs(config['network']))
          else:
              self.network_config = config
              
          if 'simulation' in config:
              self.simulation_config = load_json(self.subs(config['simulation']))
          else:
              self.simulation_config = None
              
          for m in self.network_config['manifest']:
              path = self.subs(self.network_config['manifest'][m])
              self.substitutes[m] = path
          
          if 'id' in self.parameters:
              id = self.parameters['id']
          else:
              id = 'SonataNetwork'
              
          if id[0].isdigit():  # id like 9_cells is not a valid id for NeuroML
              id='NML2_%s'%id
          
          ########################################################################
          #  Feed the handler the info on the network   
          
          self.handler = handler
      
          notes = "Network read in from Sonata: %s"%main_config_filename
          handler.handle_document_start(id, notes)
          
          handler.handle_network(id, notes)
          
          self.node_types = {}
          
          ########################################################################
          #  Get info from nodes files    
          
          for n in self.network_config['networks']['nodes']:
              nodes_file = self.subs(n['nodes_file'])
              node_types_file = self.subs(n['node_types_file'])
              
              print_v("\nLoading nodes from %s and %s"%(nodes_file,node_types_file))
  
              h5file=tables.open_file(nodes_file,mode='r')
  
              print_v("Opened HDF5 file: %s"%(h5file.filename))
              self.parse_group(h5file.root.nodes)
              h5file.close()
              self.node_types[self.current_sonata_pop] = load_csv_props(node_types_file)
              self.current_sonata_pop = None
              
          
          
          ########################################################################
          #  Get info from edges files    
          
          self.edges_info = {}
          self.conn_info = {}
          
          if 'edges' in self.network_config['networks']:
              for e in self.network_config['networks']['edges']:
                  edges_file = self.subs(e['edges_file'])
                  edge_types_file = self.subs(e['edge_types_file'])
  
                  print_v("\nLoading edges from %s and %s"%(edges_file,edge_types_file))
  
                  h5file=tables.open_file(edges_file,mode='r')
  
                  print_v("Opened HDF5 file: %s"%(h5file.filename))
                  self.parse_group(h5file.root.edges)
                  h5file.close()
                  self.edges_info[self.current_edge] = load_csv_props(edge_types_file)
                  self.current_edge = None
              
              
          ########################################################################
          #  Use extracted node/cell info to create populations
              
          for sonata_pop in self.cell_info:
              types_vs_pops = {}
              
              for type in self.cell_info[sonata_pop]['type_count']:
                  node_type_info = self.node_types[sonata_pop][type]
                  
                  model_name_type = node_type_info['model_name'] if 'model_name' in node_type_info \
                                    else (node_type_info['pop_name'] if 'pop_name' in node_type_info else node_type_info['model_type'])
                  
                  model_type = node_type_info['model_type']
                  model_template = node_type_info['model_template'] if 'model_template' in node_type_info else '- None -'
                  
                  nml_pop_id = '%s_%s_%s'%(sonata_pop,model_name_type,type) 
                      
                  print_v(" - Adding population: %s which has model info: %s"%(nml_pop_id, node_type_info))
                  
                  size = self.cell_info[sonata_pop]['type_count'][type]
                  
                  if model_type=='point_process' and model_template=='nrn:IntFire1':
                      raise Exception('Point process model not currently supported: %s\nTry expressing the I&F cell in NEST format with nest:iaf_psc_alpha'%model_template)
                      pop_comp = 'cell_%s'%nml_pop_id #model_template.replace(':','_')
                      self.pop_comp_info[pop_comp] = {}
                      self.pop_comp_info[pop_comp]['model_type'] = model_type
                      
                      dynamics_params_file = self.subs(self.network_config['components']['point_neuron_models_dir']) +'/'+node_type_info['dynamics_params']
                      self.pop_comp_info[pop_comp]['dynamics_params'] = load_json(dynamics_params_file)
                      
                  elif model_type=='point_process' and model_template=='nest:iaf_psc_alpha':
                      pop_comp = 'cell_%s'%nml_pop_id # = model_template.replace(':','_')
                      self.pop_comp_info[pop_comp] = {}
                      self.pop_comp_info[pop_comp]['model_type'] = model_type
                      self.pop_comp_info[pop_comp]['model_template'] = model_template
                      
                      dynamics_params_file = self.subs(self.network_config['components']['point_neuron_models_dir']) +'/'+node_type_info['dynamics_params']
                      self.pop_comp_info[pop_comp]['dynamics_params'] = load_json(dynamics_params_file)
                  else:
                      pop_comp = DUMMY_CELL
                      self.pop_comp_info[pop_comp] = {}
                      self.pop_comp_info[pop_comp]['model_type'] = pop_comp
                      
                  self.nml_pop_vs_comps[nml_pop_id] = pop_comp
                      
                  properties = {}
  
                  properties['type_id']=type
                  properties['sonata_population']=sonata_pop
                  properties['region']=sonata_pop
                  
                  for i in node_type_info:
                      properties[i]=node_type_info[i]
                      if i=='ei':
                          properties['type']=node_type_info[i].upper()
                          
                  color = '%s %s %s'%(self.myrandom.random(),self.myrandom.random(),self.myrandom.random())
                  
                  try:
                      import opencortex.utils.color as occ
                      interneuron = 'SOM' in nml_pop_id or 'PV' in nml_pop_id
                      if 'L23' in nml_pop_id:
                          color = occ.L23_INTERNEURON if interneuron else occ.L23_PRINCIPAL_CELL
                          pop.properties.append(neuroml.Property('region','L23'))
                      if 'L4' in nml_pop_id:
                          color = occ.L4_INTERNEURON if interneuron else occ.L4_PRINCIPAL_CELL
                          pop.properties.append(neuroml.Property('region','L4'))
                      if 'L5' in nml_pop_id:
                          color = occ.L5_INTERNEURON if interneuron else occ.L5_PRINCIPAL_CELL
                          pop.properties.append(neuroml.Property('region','L5'))
                      if 'L6' in nml_pop_id:
                          color = occ.L6_INTERNEURON if interneuron else occ.L6_PRINCIPAL_CELL
                          pop.properties.append(neuroml.Property('region','L6'))
                  except:
                      pass # Don't specify a particular color, use random, not a problem...
                      
                  properties['color']=color
                  if True or not 'locations' in self.cell_info[sonata_pop]['0']:
                      
                      properties={}  #############  temp for LEMS...
                      
                  if model_type != 'virtual':
                  
                      self.handler.handle_population(nml_pop_id, 
                                               pop_comp, 
                                               size,
                                               component_obj=None,
                                               properties=properties)
                                           
                  types_vs_pops[type] = nml_pop_id
              self.cell_info[sonata_pop]['pop_count'] = {}  
              self.cell_info[sonata_pop]['pop_map'] = {}   
              
              for i in self.cell_info[sonata_pop]['types']:
                  
                  pop = types_vs_pops[self.cell_info[sonata_pop]['types'][i]]
                  
                  if not pop in self.cell_info[sonata_pop]['pop_count']:
                      self.cell_info[sonata_pop]['pop_count'][pop] = 0
                      
                  index = self.cell_info[sonata_pop]['pop_count'][pop]
                  self.cell_info[sonata_pop]['pop_map'][i] = (pop, index)
                  
                  if not pop in self.nml_ids_vs_gids:
                      self.nml_ids_vs_gids[pop] = {}
                  self.nml_ids_vs_gids[pop][index] = (sonata_pop, i)
                  
                  if i in self.cell_info[sonata_pop]['0']['locations']: 
                      if not pop in self.nml_pops_having_locations: 
                          self.nml_pops_having_locations.append(pop)
                      pos = self.cell_info[sonata_pop]['0']['locations'][i]
                      #print('Adding pos %i: %s'%(i,pos))
                      self.handler.handle_location(index, 
                                                   pop, 
                                                   pop_comp, 
                                                   pos['x'] if 'x' in pos and pos['x'] is not None else 0, 
                                                   pos['y'] if 'y' in pos and pos['y'] is not None else 0, 
                                                   pos['z'] if 'z' in pos and pos['z'] is not None else 0)
                  
                  self.cell_info[sonata_pop]['pop_count'][pop]+=1
                  
                  
          ########################################################################
          #  Load simulation info into self.simulation_config
          
          if self.simulation_config:
              
              if self.simulation_config:
                  for m in self.simulation_config['manifest']:
                      path = self.subs(self.simulation_config['manifest'][m])
                      self.substitutes[m] = path
                      
              for s1 in ['output']:
                  for k in self.simulation_config[s1]:
                      self.simulation_config[s1][k] = self.subs(self.simulation_config[s1][k])
                      
              for s1 in ['inputs']:
                  for s2 in self.simulation_config[s1]:
                      for k in self.simulation_config[s1][s2]:
                          self.simulation_config[s1][s2][k] = self.subs(self.simulation_config[s1][s2][k])
                          
              if 'node_sets_file' in self.simulation_config:
                  node_sets = load_json(self.subs(self.simulation_config['node_sets_file']))
                  self.simulation_config['node_sets'] = node_sets 
                  
              if not 'node_sets' in self.simulation_config:
                  self.simulation_config['node_sets'] = {}
                  
              for sonata_pop in self.cell_info:
                  self.node_set_mappings[sonata_pop] = {}
                  for sindex in self.cell_info[sonata_pop]['pop_map']:
                      nml_pop = self.cell_info[sonata_pop]['pop_map'][sindex][0]
                      nml_index = self.cell_info[sonata_pop]['pop_map'][sindex][1]
  
                      # Add all in this sonata_pop to a 'node_set' named after the sonata_pop
                      if not nml_pop in self.node_set_mappings[sonata_pop]:
                          self.node_set_mappings[sonata_pop][nml_pop] = []
                      self.node_set_mappings[sonata_pop][nml_pop].append(nml_index)
                  
              #pp.pprint(self.simulation_config)
              #pp.pprint(self.pop_comp_info)
                  
              for node_set in self.simulation_config['node_sets']:
                  self.node_set_mappings[node_set] = {}
                  node_set_props = self.simulation_config['node_sets'][node_set]
                  #print_v('===========Checking which cells in pops match node_set: %s = %s'%(node_set,node_set_props))
                  
                  for sonata_pop in self.cell_info:
                      for sindex in self.cell_info[sonata_pop]['pop_map']:
                          #print('Does %s %s match %s?'%(sonata_pop, sindex, node_set_props))
                          
                          type = self.cell_info[sonata_pop]['types'][sindex]
                          type_info = self.node_types[sonata_pop][type]
                          nml_pop = self.cell_info[sonata_pop]['pop_map'][sindex][0]
                          nml_index = self.cell_info[sonata_pop]['pop_map'][sindex][1]
                          
                          if 'population' in node_set_props and node_set_props['population'] == sonata_pop:
                              if 'node_id' in node_set_props and sindex in node_set_props['node_id']:
                                  if not nml_pop in self.node_set_mappings[node_set]:
                                      self.node_set_mappings[node_set][nml_pop] = []
                                  self.node_set_mappings[node_set][nml_pop].append(nml_index)
                                  
                          matches = _matches_node_set_props(type_info, node_set_props)
                          #print_v('Node %i in %s (NML: %s[%i]) has type %s (%s); matches: %s'%(sindex, sonata_pop, nml_pop, nml_index, type, type_info, matches))
                          if matches:
                              if not nml_pop in self.node_set_mappings[node_set]:
                                  self.node_set_mappings[node_set][nml_pop] = []
                              self.node_set_mappings[node_set][nml_pop].append(nml_index)
              
              ##pp.pprint(self.node_set_mappings)
              
      
          ########################################################################
          #  Extract info from inputs in simulation_config
          #pp.pprint(self.simulation_config)
          
          for input in self.simulation_config['inputs']:
              info = self.simulation_config['inputs'][input]
              #print_v(" - Adding input: %s which has info: %s"%(input, info)) 
              
              self.input_comp_info[input] = {}
              self.input_comp_info[input][info['input_type']] = {}
              
              node_set = info['node_set']
              
              if info['input_type'] == 'current_clamp':
                  
                  comp = 'PG_%s'%input
                  self.input_comp_info[input][info['input_type']][comp] = {'amp':info['amp'],'delay':info['delay'],'duration':info['duration']}
                  
                  for nml_pop_id in self.node_set_mappings[node_set]:
                      
                      input_list_id = 'il_%s_%s'%(input,nml_pop_id)
                      indices = self.node_set_mappings[node_set][nml_pop_id]
                      
                      self.handler.handle_input_list(input_list_id, 
                                                     nml_pop_id, 
                                                     comp, 
                                                     len(indices))
                      count = 0
                      for index in indices:
                          self.handler.handle_single_input(input_list_id, 
                                                           count, 
                                                           cellId = index, 
                                                           segId = 0, 
                                                           fract = 0.5)
                          count+=1
                      
                  
              elif info['input_type'] == 'spikes':
                  node_info = self.cell_info[node_set]
                  
                  from pyneuroml.plot.PlotSpikes import read_sonata_spikes_hdf5_file
  
                  ids_times = read_sonata_spikes_hdf5_file(self.subs(info['input_file']))
                  
                  for id in ids_times:
                      times = ids_times[id]
                      if id in node_info['pop_map']:
                          nml_pop_id, cell_id = node_info['pop_map'][id] 
                          print_v("Cell %i in Sonata node set %s (cell %s in nml pop %s) has %i spikes"%(id, node_set, nml_pop_id, cell_id, len(times)))
  
                          component = '%s__%i'%(nml_pop_id,cell_id)
  
                          self.input_comp_info[input][info['input_type']][component] ={'id': cell_id, 'times': times}
  
                          '''
                          input_list_id = 'il_%s_%i'%(input,cell_id)
                          self.handler.handle_input_list(input_list_id, 
                                                         nml_pop_id, 
                                                         component, 
                                                         1)
  
                          self.handler.handle_single_input(input_list_id, 
                                                            0, 
                                                            cellId = cell_id, 
                                                            segId = 0, 
                                                            fract = 0.5)
                                                           '''
                      else:
                          print_v("Cell %i in Sonata node set %s NOT FOUND!"%(id, node_set))
                  
              else:
                  raise Exception("Sonata input type not yet supported: %s"%(info['input_type']))
  
  
          ########################################################################
          #  Use extracted edge info to create connections
          
          projections_created = []
          for conn in self.conn_info:
              
              pre_node = self.conn_info[conn]['pre_node']
              post_node = self.conn_info[conn]['post_node']
              
              for i in range(len(self.conn_info[conn]['pre_id'])):
                  pre_id = self.conn_info[conn]['pre_id'][i]
                  post_id = self.conn_info[conn]['post_id'][i]
                  nsyns = self.conn_info[conn]['nsyns'][i] if 'nsyns' in self.conn_info[conn] else 1
                  type = self.conn_info[conn]['edge_type_id'][i]
                  #print_v('   Conn with %i syns, type %s: %s(%s) -> %s(%s)'%(nsyns,type,pre_node,pre_id,post_node,post_id))
                  pre_pop,pre_i = self.cell_info[pre_node]['pop_map'][pre_id]
                  post_pop,post_i = self.cell_info[post_node]['pop_map'][post_id]
                  #print_v('   Mapped: Conn %s(%s) -> %s(%s)'%(pre_pop,pre_i,post_pop,post_i))
                  # print self.edges_info[conn][type]
                  #print self.cell_info[pre_node]
                  #print 11
                  #print self.node_types[pre_node]
                  #print 22
                  cell_type_pre = self.cell_info[pre_node]['types'][pre_id]
                  #print cell_type_pre
                  #print 444
                  pop_type_pre = self.node_types[pre_node][cell_type_pre]['model_type']
                  #print pop_type_pre
                  #print 333
                  
                  synapse = self.edges_info[conn][type]['dynamics_params'].split('.')[0]
                  self.syn_comp_info[synapse] = {}
                  #print self.edges_info[conn][type]
                  #pp.pprint(self.init_substitutes)
                  #pp.pprint(self.substitutes)
                  dynamics_params_file = self.subs(self.network_config['components']['synaptic_models_dir']) +'/'+self.edges_info[conn][type]['dynamics_params']
                  #print_v('Adding syn %s (at %s)'%(self.edges_info[conn][type]['dynamics_params'], dynamics_params_file))
                  
                  #TODO: don't load this file every connection!!!
                  self.syn_comp_info[synapse]['dynamics_params'] = load_json(dynamics_params_file)
                  proj_id = '%s_%s_%s'%(pre_pop,post_pop,synapse)
                  
                  sign = self.syn_comp_info[synapse]['dynamics_params']['sign'] if 'sign' in self.syn_comp_info[synapse]['dynamics_params'] else 1
                  weight = self.edges_info[conn][type]['syn_weight'] if 'syn_weight' in self.edges_info[conn][type] else 1.0
                  
                  syn_weight_edge_group_0 = self.conn_info[conn]['syn_weight_edge_group_0'][i] if 'syn_weight_edge_group_0' in self.conn_info[conn] else None
                  
                  # Assume this overrides value from csv file...
                  if syn_weight_edge_group_0:
                      weight = syn_weight_edge_group_0
                  
                  #print_v('Adding syn %s (at %s), weight: %s, sign: %s, nsyns: %s'%(self.edges_info[conn][type]['dynamics_params'], dynamics_params_file, weight, sign, nsyns))
                  
                  weight_scale = 0.001
                  if 'level_of_detail' in self.syn_comp_info[synapse]['dynamics_params']:
                      weight_scale = 1
                  
                  
                  weight=weight_scale * sign * weight * nsyns
                  
                  delay = self.edges_info[conn][type]['delay'] if 'delay' in self.edges_info[conn][type] else 0
                  
                  if not pop_type_pre == 'virtual':
                      if not proj_id in projections_created:
  
                          self.handler.handle_projection(proj_id, 
                                               pre_pop, 
                                               post_pop, 
                                               synapse)
  
                          projections_created.append(proj_id)
  
                      self.handler.handle_connection(proj_id, 
                                                   i, 
                                                   pre_pop, 
                                                   post_pop, 
                                                   synapse, \
                                                   pre_i, \
                                                   post_i, \
                                                   weight=weight, \
                                                   delay=delay)
                                                   
                  else:
                      component = '%s__%i'%(pre_pop,pre_i)
                      #print_v('   ---  Connecting %s to %s[%s]'%(component, post_pop, post_i))
                      
                      
                      #self.input_comp_info[input][info['input_type']][component] ={'id': cell_id, 'times': times}
  
  
                      input_list_id = 'il_%s_%s_%i_%i'%(component,post_pop,post_i,i)
                      
                      self.handler.handle_input_list(input_list_id, 
                                                     post_pop, 
                                                     component, 
                                                     1)
  
                      self.handler.handle_single_input(input_list_id, 
                                                        0, 
                                                        cellId = post_i, 
                                                        segId = 0, 
                                                        fract = 0.5,
                                                        weight=weight)
                          
          """
          print('~~~~~~~~~~~~~~~')
          print('node_types:')
          pp.pprint(self.node_types)
          print('~~~~~~~~~~~~~~~')
          print('cell_info:')
          pp.pprint(self.cell_info)
          print('================')"""