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('================')"""