def Xstr(gas, fuel_dict, phi, oxid_dict):


	mech = gas#ct.Solution(mech_file)
	element = mech.element_names
	species = mech.species
	reaction = mech.reaction

	C_fuel = 0
	H_fuel = 0
	O_fuel = 0
	O_oxid = 0

	Xstr = ''
	for fuel_sp in fuel_dict.keys():
		val = fuel_dict[fuel_sp]
		Xstr += (fuel_sp + ':' + str(1.0*val*phi) + ', ')
		atoms = species(fuel_sp).composition
		if 'C' in atoms.keys():
			C_fuel = C_fuel + atoms['C'] * val
		if 'c' in atoms.keys():
			C_fuel = C_fuel + atoms['c'] * val
		if 'H' in atoms.keys():
			H_fuel = H_fuel + atoms['H'] * val
		if 'h' in atoms.keys():
			H_fuel = H_fuel + atoms['h'] * val
		if 'O' in atoms.keys():
			O_fuel = O_fuel + atoms['O'] * val
		if 'o' in atoms.keys():
			O_fuel = O_fuel + atoms['o'] * val


	for oxid_sp in oxid_dict.keys():
		val = oxid_dict[oxid_sp]
		atoms = species(oxid_sp).composition
		if 'O' in atoms.keys():
			O_oxid = O_oxid + atoms['O'] * val
		if 'o' in atoms.keys():
			O_oxid = O_oxid + atoms['o'] * val
		if 'C' in atoms.keys():
			O_oxid -= atoms['C'] * val * 2
		if 'c' in atoms.keys():
			O_oxid -= atoms['c'] * val * 2
		if 'H' in atoms.keys():
			O_oxid -= atoms['H'] * val * 0.5
		if 'h' in atoms.keys():
			O_oxid -= atoms['h'] * val * 0.5


	O_stoic = C_fuel * 2.0 + H_fuel / 2.0 - O_fuel
	n_oxid = O_stoic / O_oxid

	for oxid_sp in oxid_dict.keys():
		val = oxid_dict[oxid_sp]
		Xstr += (oxid_sp + ':' + str(n_oxid*val) + ', ')

	Xstr = Xstr[:-2]
	return Xstr def senkin(soln, atm, T0, X0, if_half=True, dir_raw=None, if_fine=False):

	print 'if_half = '+str(if_half)

	#if soln.n_species > 100:
	#    verbose=True
	#else:
	#    verbose=False

	dT_ign = 400

   # if_half=True

	verbose=True

	cpu0 = time.time()

	print '>'*30
	print 'senkin for ['+ X0 + '] at '+ str(atm)+'atm' + ' and '+str(T0)+'K'
	print '<'*30
	
	p = ct.one_atm * atm

	"""
	find the ignition process

	:param mech:    mechanism
	:param p:       pressure (Pa)
	:param T0:      temperature (K) of the inlet flow
	:param T:       initial temperature (K) of the reactor
	:param tau:     residence time (s) of the reactor

	:return reactor:

	"""

	if if_half:
		max_dT = 40
	else:
		max_dT = 100

	if if_fine:
		n_loop = 2
	else:
		n_loop = 1

	i_loop = 0


	# @@@@@@@@@@@@@@@@@@@@@@@@@@
	# loop until grid resolution found
	while i_loop < n_loop:
		i_loop += 1
		print '-'*5
		print 'i_loop = '+str(i_loop)

		soln.TPX = T0, p, X0
		reactor = ct.IdealGasConstPressureReactor(soln)
		network = ct.ReactorNet([reactor])

		if i_loop == 2:
			dt = 1.0*tau_ign/20
			max_dT /= 4.0
		else:
			dt = 0.1
		t = 0

		if verbose:
			print 'dt0 = '+str(dt)
		T_prev = T0

		ii_add = []
		#TT = []
		raw = None
		raw_all = None

		# @@@@@@@@@@@@@@@@@@@@@@@@@@
		# loop until ignited/equil reached
		while True:
			t_prev = t

			# @@@@@@@@@@@@@@@@@@@@@@@@@@
			# loop until find the proper time step
			while True:
				if i_loop == 2:
					dt = min(dt, 1.0*tau_ign/20)
				t += dt
				t = float(str(t))

				if verbose:
					print 'now, '+str(t)

				# find T at tried t ------------
				# corresponding data stored in raw_all[i_add]

				new_sim = True
				if raw_all is not None:
					if (t in raw_all['axis0']):
						new_sim = False
						if verbose:
							print 'calculated before, loaded !!!!!!'
						i_add = raw_all['axis0'].index(t)
						T = raw_all['temperature'][i_add]

				if new_sim:
					network.advance(t)
					if verbose:
						print 'advanced, T = '+str(soln.T)
					raw_all = soln2raw(t, 'time', soln, raw_all)
					i_add = len(raw_all['axis0'])-1
					T = soln.T

				# decide if tried t is proper -------
				
				dT = abs(T - T_prev)
				if dT > max_dT or (if_half and T-T0>dT_ign+4):
					# if dT is too large, this t is not okay, decrease time step and start over
					t = t_prev
					dt /= 2.0
					soln.TPX = T0, p, X0
					reactor = ct.IdealGasConstPressureReactor(soln)
					network = ct.ReactorNet([reactor])
				
				else:
					# if dT is too small, this t is okay, we increase time step for next t
					if (soln.T - T0 < dT_ign and dT < 10):
						dt *= 2.0
					if (soln.T - T0 >= dT_ign and dT < 50):
						dt *= 5.0
					break

			ii_add.append(i_add)
			T_prev = T

			if if_half:
				if soln.T - T0 > dT_ign:
					print 'ignited'
					tau_ign = t
					break
			else:
				if soln.T - T0 > 1500 and dT < 1:
					print 'equilibrium reached'
					break




	raw = slice_raw(raw_all, ii_add)
	print 'n_points = ' + str(len(raw['axis0'])) + '/' + str(len(raw_all['axis0']))
	print 'CPU time = '+str(time.time() - cpu0)


	if dir_raw is not None:

		path_raw = os.path.join(dir_raw,'raw.npz')
		raw = save_raw_npz(raw, path_raw)
		save_raw_csv(raw, soln, dir_raw)


	return raw#['axis0'][-1] def S_curve(soln_in, soln_out, atm, T0, X0, dir_raw=None):

    #raise

    print '>'*30
    print 'psr for ['+ X0 + '] at '+ str(atm)+'atm' + ' and '+str(T0)+'K'
    print '<'*30

    if soln_in.n_species > 100:
        verbose = True
    else:
        verbose = False


    p = ct.one_atm * atm

    dT_turn = 500.0
    dT_froze = 50.0
    dT_try = 10.0

    status = 0

    tau = 1
    #tt = []
    #TT = []
    raw = None
    T_ini = 2500
    T_burn = None
    tau_burn = tau
    X_ini = X0

    tau_r = 0.5

    #soln_in = ct.Solution(mech)
    #soln_out = ct.Solution(mech)

    while True:

        soln_i = psr_ss(soln_in, soln_out, p, T0, T_ini, X0, X_ini, tau)
        if soln_i is None:
            break

        T = soln_i.T
        X = soln_i.concentrations
        if verbose:
            print str(int(T_ini)) + ', ' + str(int(T)) + ', ' + str(tau)
        if T_burn is None:
            T_burn = T

        if abs(T_burn - T) > 50:
            # if extinction happens or dT too large
            
            if tau_r > 0.999:
                if verbose:            
                    print 'finished, tau_r = '+str(tau_r)
                break
            else:
                tau = tau_burn
                tau_r = tau_r + (1-tau_r)*0.5
                if verbose:            
                    print 'refined, tau_r = '+str(tau_r)
        else:
            T_burn = T
            X_ini = X
            tau_burn = tau

            raw = soln2raw(tau, 'residence_time', soln_i, raw)
            save_raw_npz(raw, os.path.join(dir_raw,'raw_temp.npz'))

        T_ini = T_burn
        tau *= tau_r

    raw = save_raw_npz(raw, os.path.join(dir_raw,'raw.npz'))
    save_raw_csv(raw, soln_in, dir_raw)

    return raw def skeletal(detailed_folder, sk_folder, species_kept, notes=None):

	if not os.path.exists(sk_folder):
		os.makedirs(sk_folder)

	if detailed_folder != sk_folder:
		shutil.copyfile(os.path.join(detailed_folder,'therm.dat'), os.path.join(sk_folder,'therm.dat'))
	trandat = os.path.join(detailed_folder,'tran.dat')
	try:
		ft = open(trandat,'r')
		ft.close()
		shutil.copyfile(trandat, os.path.join(sk_folder,'tran.dat'))
	except IOError:
		pass

	sk_inp = os.path.join(sk_folder,'chem.inp')

	mech = build_mech(detailed_folder,overwrite=False)
	rxn_all = mech['reaction']

	f = open(sk_inp,'w')
	if notes is not None:
		for note in notes:
			f.write(note+'\n')
	f.write('\n')

	f.write('ELEMENTS\n')
	for e in mech['element'].keys():
		f.write(e + ' ')
	f.write('\nEND\n\n')



	f.write('SPECIES\n')
	n = 0
	for s in species_kept:
		f.write(s + ' ')
		n += 1
		if n == 5:
			f.write('\n')
			n = 0
	if n != 0:
		f.write('\n')
	f.write('END\n\n')


	f.write('REACTIONS\n')
	rxn_kept = []
	for rxn in rxn_all:
		if all(member in species_kept for member in rxn_all[rxn]['member'].keys()):
			n_ln = 0
			for info in rxn_all[rxn]['info']:
				if n_ln > 0:
					f.write('    ')


				if '/' in info and \
					('LOW' not in info.upper()) and ('TROE' not in info.upper()) \
					and ('REV' not in info.upper()):

					# this line describes three-body collision * efficiency *
					# we should remove these not included in mech

					ss = info.split('/')
					info = ''
					for i in range(len(ss)):
						s = ss[i].strip()
						if s in species_kept:
							info += (ss[i] + '/' + ss[i+1] + '/')
					

				f.write(info.strip() + '\n')
				n_ln += 1
			if n_ln > 1:
				f.write('\n')
	f.write('END\n\n')
	f.close() def ck2cti(folder):

    ok = True

    options = dict()
    options['--input'] = os.path.join(folder,'chem.inp')
    options['--thermo'] = os.path.join(folder,'therm.dat')
    try:
        open(os.path.join(folder,'tran.dat'))
        options['--transport'] = folder + 'tran.dat'
    except IOError:
        pass
    options['--permissive'] = True
    main(options)

    return ok def load_raw(path_raw):
    raw_file = np.load(open(path_raw, 'rb'))
    raw = dict()
    for key in raw_file.keys():
        if key in ['mole_fraction','net_reaction_rate']:
            raw[key] = np.matrix(raw_file[key])
        else:
            raw[key] = raw_file[key]
    return raw def GPS_algo(soln, flux_graph, source, target, path_save=None, \
    K=1, alpha=0.1, beta=0.5, normal='max', iso=None, overwrite=False, \
    raw='unknown',notes=None, gamma=None):

    reaction = soln.reaction

    """
    """

    """ --------------------------------
    check if results already exist
    if so, load and return
    -------------------------------- """

    if path_save is not None:
        if overwrite is False:
            try:
                GPS_results = json.load(open(path_save, 'r'))
                return GPS_results
            except IOError:
                pass

    """ --------------------------------
    if not, then compute, save, and return
    -------------------------------- """

    print 'performing GPS'+str(path_save)

    arrow = ' --> '
    global_path = {}
    hub_GP = {}

    # ----------------------------------------------------------
    # score = max(in_deg, out_deg)

    in_deg = flux_graph.in_degree(weight='flux')
    out_deg = flux_graph.out_degree(weight='flux')

    if normal is 'source':
        norm_deg = out_deg[source]
    elif normal is 'max':
        try:
            norm_deg = max([max(out_deg.values()), max(in_deg.values())])
        except:
            return {'species':dict()}

    score = {node: max(in_deg[node], out_deg[node]) / norm_deg for node in in_deg.keys()}

    # ----------------------------------------------------------
    # hubs = species whose score is above crit_score

    hubs, iso_kept = find_hubs(score, iso, alpha, gamma)
    n_hub = len(hubs)

    species_kept = dict()
    for hub in hubs:
        species_kept[hub] = dict()        
        species_kept[hub]['by_alpha'] = True
        species_kept[hub]['by_K'] = []
        species_kept[hub]['by_beta'] = []

    if target is not None:


        for hub in hubs:
            if hub not in hub_GP.keys():
                hub_GP[hub] = {}
                hub_GP[hub]['score'] = score[hub]
                hub_GP[hub]['global_path'] = []

            # ----------------------------------------------------------
            # a global path(way), GP, is the pathway, p0, from source to hub, plus the pathway, p1, from hub to target
            # however, we should further ensure there's no cycle in GP

            pp0 = []
            if K > 1:
                dd0, pp0 = yen(flux_graph, source, hub, k=K, weight='1/flux')
                dd1, pp1 = yen(flux_graph, hub, target, k=K, weight='1/flux')
            else:
                try:
                    p0 = nx.shortest_path(flux_graph, source, hub, weight='1/flux')
                    p1 = nx.shortest_path(flux_graph, hub, target, weight='1/flux')
                    pp0 = [p0]
                    pp1 = [p1]
                except:# nx.NetworkXNoPath:
                    print 'no path found for: '+str(source)+' --> '+str(hub)+' --> '+str(target)
                    pass

            if bool(pp0):
                k = 0
                for p0 in pp0:

                    k += 1
                    for p1 in pp1:

                        # although what returned from shortest path is ensured to not contain any cycle
                        # cycle may appear when you concatenate p0 and p1

                        no_cycle = True
                        p1 = p1[1:]
                        for s in p1:
                            if s in p0:
                                no_cycle = False
                                break

                        if no_cycle:
                            p = p0 + p1
                            GP = arrow.join(p)
                            hub_GP[hub]['global_path'].append(GP)

                            if GP not in global_path.keys():
                                global_path[GP] = {}
                                global_path[GP]['member'] = p

        # ---------------------------------     
        st_list = []
        for GP in global_path.keys():
            for sp in global_path[GP]['member']:
                if sp not in species_kept.keys():
                    species_kept[sp] = {}
                    species_kept[sp]['by_alpha'] = False
                    species_kept[sp]['by_K'] = []
                    species_kept[sp]['by_beta'] = []

                species_kept[sp]['by_K'].append(GP)

            if beta == 0:
                continue

            for isp in range(len(global_path[GP]['member']) - 1):
                s = global_path[GP]['member'][isp]
                t = global_path[GP]['member'][isp + 1]
                st = s + arrow + t
                if st not in st_list:
                    #print '-'*10
                    #print st
                    st_list.append(st)
                    rxn_beta = flux_graph[s][t]['member']
                    sum_beta = sum(rxn_beta.values())
                    rxn_top = keys_sorted(rxn_beta)
                    total_beta = 0.0
                    for rxn in rxn_top:
                        #print reaction(int(rxn)).equation + ': ' + str(int(rxn_beta[rxn]/sum_beta*100)) +'%'
                        st_rxn = st + ', ' + reaction(int(rxn)).equation

                        for sp in (reaction(int(rxn)).products.keys() + reaction(int(rxn)).reactants.keys()):
                            if sp not in species_kept.keys():
                                species_kept[sp] = {}
                                species_kept[sp]['by_alpha'] = False
                                species_kept[sp]['by_K'] = []
                                species_kept[sp]['by_beta'] = []   

                            if st_rxn not in species_kept[sp]['by_beta']:
                                species_kept[sp]['by_beta'].append(st_rxn)

                        total_beta += rxn_beta[rxn]/sum_beta
                        if total_beta >= beta:
                            break


    GPS_results = dict()

    GPS_results['parameter'] = {}
    GPS_results['parameter']['K'] = K
    GPS_results['parameter']['alpha'] = alpha
    GPS_results['parameter']['beta'] = beta
    GPS_results['parameter']['gamma'] = gamma
    GPS_results['parameter']['source'] = source
    GPS_results['parameter']['target'] = target
    GPS_results['parameter']['raw'] = raw
    GPS_results['parameter']['notes'] = notes

    GPS_results['isomer'] = iso_kept
    GPS_results['global_path'] = global_path
    GPS_results['hubs'] = hub_GP
    GPS_results['species'] = species_kept

    GPS_results['summary'] = {}
    GPS_results['summary']['n_hub'] = n_hub
    GPS_results['summary']['n_species_kept'] = len(species_kept.keys())
    GPS_results['summary']['n_global_path'] = len(global_path.keys())

    if path_save is not None:
        json.dump(GPS_results, open(path_save, 'w'))
    return GPS_results def build_flux_graph(soln, raw, traced_element, path_save=None, overwrite=False, i0=0, i1='eq', constV=False):

	"""
	:param mechanism:        type = dict, keys include "species", "reaction", "element", etc
	:param raw:              type = dict, keys include "mole_fraction", "net_reaction_rate", etc
	:param traced_element:   type = str
	:param i0:               type = int, specifying the starting point of the considered interval of the raw data
	:param i1:               type = int or str, specifying the ending point of the considered interval of the raw data
	:return flux graph:      type = networkx object, will be also saved as a .json file,
	"""

	element = soln.element_names
	species = soln.species
	reaction = soln.reaction
	n_rxn = soln.n_reactions

	""" --------------------------------
	check if results already exist, if so, load
	-------------------------------- """

	if path_save is not None:
		if overwrite is False:
			try:
				data = json.load(open(path_save, 'r'))
				flux_graph = json_graph.node_link_graph(data)
				return flux_graph
			except IOError:
				pass

	""" --------------------------------
	if not, then compute, and save
	-------------------------------- """

	# ---------------------------------------------
	# check if traced_element is legal

	if traced_element not in element:
		raise('traced element ' + traced_element + ' is not listed in mechanism')

	# ---------------------------------------------
	# find the reaction rate during the considered interval
	# unit will be converted to mole/sec


	rr = np.reshape(raw['net_reaction_rate'][i0,:],[n_rxn,1]) 
	flux_graph = DiGraph()

	# -------------------------------------
	# adding edge from reactions
	# one edge may contribute from multiple reactions, the list of the contributors will be stored in edge['member']

	# note though in .cti id_rxn starts from 1, in soln.reaction, id_rxn starts from 0
	for id_rxn in range(n_rxn):


		# sp_mu is a dict, where key is species, val is net stoichiometric coefficient
		sp_mu = reaction(id_rxn).products
		for sp in reaction(id_rxn).reactants.keys():
			mu = reaction(id_rxn).reactants[sp]
			if sp in sp_mu.keys():
				sp_mu[sp] -= mu
			else:
				sp_mu[sp] = -mu

	   
		# -----------------------
		# produced is a dict, where key is sp, val is number of traced atoms
		# being transferred when this sp is produced
		produced = {}
		consumed = {}

		for sp in sp_mu.keys():
			atoms = species(sp).composition
			if traced_element in atoms.keys():
				n = int(sp_mu[sp] * atoms[traced_element] * np.sign(rr[id_rxn]))
				if n > 0:
					produced[sp] = abs(n)
				elif n < 0:
					consumed[sp] = abs(n)

		# -----------------------
		# consider this reaction only when traced element is transferred
		# note "if bool(consumed)" works the same way
		if bool(produced):
			n_sum = sum(produced.values())
			for target in produced.keys():
				for source in consumed.keys():

					n_i2j = 1.0 * produced[target] * consumed[source] / n_sum

					# note that the direction (source-->target) is already assured
					# therefore we use abs(RR) here
					dw = float(n_i2j * abs(rr[id_rxn]))

					try:
						flux_graph[source][target]['flux'] += dw
					except KeyError:
						# if this edge doesn't exist, create it
						flux_graph.add_edge(source, target)
						flux_graph[source][target]['flux'] = dw
						flux_graph[source][target]['member'] = {}

					flux_graph[source][target]['member'][str(id_rxn)] = dw
					flux_graph[source][target]['1/flux'] = 1.0 / flux_graph[source][target]['flux']

	# -------------------------------------
	# save the graph using json, which is fast, and human-readable

	data = json_graph.node_link_data(flux_graph)
	json.dump(data, open(path_save, 'w'))
	#print 'graph saved as',path_save

	return flux_graph def st2name(i_pnt, e, source, target):
    name_gps = 'point '+ str(i_pnt) + ', graph '+ e
    if source is not None:
        name_gps += ', '+source + ' to '+target
    return name_gps def	find_GPSA(dir_raw, GP_dict, soln, dnR, fuel_comp, n_break=0):

	GP_name = GP_dict['name']
	traced = GP_dict['traced']
	dir_save = os.path.join(dir_raw,'GPSA',traced)
	if not os.path.exists(dir_save):
		os.makedirs(dir_save)
	path_save = os.path.join(dir_save, shorten_GP_name(GP_name)+'.json')
	print('[GPSA] to be saved: '+str(path_save))
	
	# =====================================
	# if previously computed, return False

	GPSA_all = dict()
	if os.path.exists(path_save):
		GPSA_all = json.load(open(path_save,'r'))
		if GP_name in GPSA_all.keys():
			return False

	# =====================================
	# if not loaded, compute these results


	print 'computing GPSA for '+GP_name

	traced = GP_dict['traced']
	GP_member = GP_dict['member']
	dir_graph = os.path.join(dir_raw,'graph')	
	path_raw = os.path.join(dir_raw,'raw.npz')
	raw = load_raw(path_raw)
	rr_mat = raw['net_reaction_rate']

	GPSA = dict()
	GPSA['R_GP'] = []		# net radical production rate associated with a Global Pathway (GP)
	GPSA['Q_GP'] = []		# net heat release rate associated with a GP
	GPSA['D_GP'] = []		# dominancy of a GP

	GPSA['R_ij'] = dict()	# net radical production rate associated with a conversion step (from the i-th species to the j-th species)
	GPSA['Q_ij'] = dict()
	GPSA['a_iji'] = dict()
	#GPSA['perc_ij'] = dict()

	for i in range(len(GP_member)-1):
		edge = find_edge_name(GP_member, i)

		GPSA['R_ij'][edge] = dict()
		GPSA['R_ij'][edge]['member'] = []
		GPSA['R_ij'][edge]['net'] = []

		GPSA['Q_ij'][edge] = dict()
		GPSA['Q_ij'][edge]['member'] = []
		GPSA['Q_ij'][edge]['net'] = []

		GPSA['a_iji'][edge] = dict()
		GPSA['a_iji'][edge]['member'] = []
		GPSA['a_iji'][edge]['net'] = []

		#GPSA['perc_ij'][edge] = dict()
		#GPSA['perc_ij'][edge]['member'] = []
		#GPSA['perc_ij'][edge]['net'] = []

	source = GP_dict['member'][0]
	traced = GP_dict['traced']
	if source not in fuel_comp.keys():
		perc_from_source = 0.0
	else:
		total_atom = 0
		for k in fuel_comp.keys():
			sp = soln.species(k)
			atom = 0.0
			if traced in sp.composition.keys():
				atom += fuel_comp[k] * sp.composition[traced]
			if source == k:
				atom_source = atom
			total_atom += atom

		perc_from_source = 1.0 * atom_source / total_atom

		#print 'total '+traced+' atoms for '+str(fuel_comp)+' is '+str(total_atom)
		#print 'source '+str(source)+' has '+str(atom_source)+' atms'
		#print 'so perc_from_source = '+str(perc_from_source)

	#GPSA['perc_ij']['from_source'] = perc_from_source


	# for each point -----------

	i_pnt = 0
	while True:

			

		path_graph = os.path.join(dir_graph, traced+'_'+str(i_pnt) + '.json')
		if not os.path.exists(path_graph): 			
			if i_pnt>n_break:
				print 'break as cannot find: '+str(path_graph)
				break
			
			else:			
				# fill this with None ---------------
				if i_pnt%10 == 0:
					print '   fill None GPSA for '+str(path_graph)

				for i in range(len(GP_member)-1):
					s = GP_member[i]
					t = GP_member[i+1]
					edge = find_edge_name(GP_member, i)
					for k in ['member','net']:
						GPSA['R_ij'][edge][k].append(None)
						GPSA['Q_ij'][edge][k].append(None)
						GPSA['a_iji'][edge][k].append(None)
						#GPSA['perc_ij'][edge][k].append(None)

				GPSA['R_GP'].append(None)
				GPSA['D_GP'].append(None)

				i_pnt += 1
				continue

		soln = raw2soln(soln, raw, i_pnt)	
		if i_pnt%10 == 0:
			print '   finding GPSA for '+str(path_graph)

		# fill this with real value ---------------
		#norm_Rpro = 0.0
		#norm_Rcon = 0.0
		for id_rxn in range(soln.n_reactions):
			dR = dnR[id_rxn] * rr_mat[i_pnt, id_rxn]
			#norm_Rpro += max(0, dR)
			#norm_Rcon += max(0, -dR)


		flux_graph = json_graph.node_link_graph(json.load(open(path_graph, 'r')))
		out_deg = flux_graph.out_degree(weight='flux')
		norm_out_deg = sum([out_deg[m] for m in fuel_comp.keys() if m in out_deg.keys()])

		flux = []
		rxn_involved = []
		sum_OMEGA_R = 0
		sum_OMEGA_Q = 0
		perc_ij_list = []

		# for each edge (conversion step) -----------------


		for i in range(len(GP_member)-1):
			s = GP_member[i]
			t = GP_member[i+1]
			edge = find_edge_name(GP_member, i)
			GPSA['R_ij'][edge]['member'].append(dict())
			GPSA['Q_ij'][edge]['member'].append(dict())
			GPSA['a_iji'][edge]['member'].append(dict())
			#GPSA['perc_ij'][edge]['member'].append(dict())

			# ------------------------------

			try:
				st = flux_graph[s][t]
			except KeyError:
				st = None

			perc_ij = None
			if st is not None:
				flux.append(st['flux'])
				perc_ij = 1.0 * st['flux'] / out_deg[s]
				perc_ij_list.append(perc_ij)

				for id_rxn_s in st['member'].keys():
					id_rxn = int(id_rxn_s)
					rxn = soln.reaction_equation(id_rxn)
					rr = rr_mat[i_pnt, id_rxn]
					if rr < 0:
						sign_rxn = -id_rxn
					else:
						sign_rxn = id_rxn

					GPSA['a_iji'][edge]['member'][i_pnt][sign_rxn] = st['member'][id_rxn_s]
					#GPSA['perc_ij'][edge]['member'][i_pnt][sign_rxn] = 1.0 * st['member'][id_rxn_s]/out_deg[s]

					OMEGA_R = float(rr * dnR[id_rxn])
					GPSA['R_ij'][edge]['member'][i_pnt][sign_rxn] = OMEGA_R					

					OMEGA_Q = float(rr * soln.delta_enthalpy[id_rxn])
					GPSA['Q_ij'][edge]['member'][i_pnt][sign_rxn] = OMEGA_Q

					if id_rxn not in rxn_involved:
						sum_OMEGA_R += OMEGA_R
						sum_OMEGA_Q += OMEGA_Q
						rxn_involved.append(id_rxn)




			# ------------------------------

			try:
				ts = flux_graph[t][s]
			except KeyError:
				ts = None

			if ts is not None:
				for id_rxn_s in ts['member'].keys():
					id_rxn = int(id_rxn_s)
					rxn = soln.reaction_equation(id_rxn)
					rr = rr_mat[i_pnt, id_rxn]
					if rr < 0:
						sign_rxn = -id_rxn
					else:
						sign_rxn = id_rxn
					GPSA['a_iji'][edge]['member'][i_pnt][sign_rxn] = - ts['member'][id_rxn_s]

			# ------------------------------

			GPSA['R_ij'][edge]['net'].append(sum(GPSA['R_ij'][edge]['member'][i_pnt].values()))
			GPSA['Q_ij'][edge]['net'].append(sum(GPSA['Q_ij'][edge]['member'][i_pnt].values()))
			GPSA['a_iji'][edge]['net'].append(sum(GPSA['a_iji'][edge]['member'][i_pnt].values()))
			#GPSA['perc_ij'][edge]['net'].append(perc_ij)


		domi_perc = gmean(perc_ij_list)*perc_from_source
		if bool(flux):
			min_flux = min(flux)
			if norm_out_deg>0:
				domi_flux = 1.0 * min_flux/norm_out_deg
			else:
				domi_flux = float('nan')
		else:
			min_flux = 0.0
			domi_flux = 0.0

		R_GP = domi_perc * sum_OMEGA_R
		Q_GP = domi_perc * sum_OMEGA_Q
		GPSA['R_GP'].append(R_GP)
		GPSA['Q_GP'].append(Q_GP)
		GPSA['D_GP'].append(domi_perc)


		i_pnt += 1


	GPSA_all[GP_name] = GPSA
	json.dump(GPSA_all, open(path_save,'w'))
	return True
import sys
import os
import json
import time
import cantera as ct
import shutil
import copy
from PyQt4 import uic
from PyQt4.QtGui import * 
from PyQt4.QtCore import * 
from src.core.def_tools import *
from src.ct.def_ct_tools import Xstr
from src.ct.senkin import senkin
from src.ct.psr import S_curve
from src.ck.def_cheminp import skeletal
from src.ct.ck2cti_GPS import ck2cti
from dialog_GPS import dialog_GPS
from dialog_database import dialog_database
from dialog_mech import dialog_mech
from dialog_view_mech import dialog_view_mech
from find_tau_ign import find_tau_ign
from src.ct.def_ct_tools import load_raw
from src.core.def_GPS import GPS_algo
from src.core.def_build_graph import build_flux_graph
from networkx.readwrite import json_graph
from src.core.def_tools import st2name
from src.core.def_GPSA import find_GPSA
							if progress.stop:
								progress.close()
								return False

							fuel = parent.project['fuel'][fuel_name]
							oxid = parent.project['oxid'][oxid_name]

							dir_raw = cond2dir(dir_desk, fuel_name, oxid_name, phi, atm, T0, \
								reactor, parent.n_digit)
							path_raw = os.path.join(dir_raw,'raw.npz')

							if os.path.exists(path_raw):
								progress.set_info('<already exists, skipped> '+dir_raw)
							else:
								progress.set_info(dir_raw)
								find_raw(soln, soln_in,\
									dir_desk, fuel,	oxid, phi, atm, T0, reactor, parent.n_digit)

							v += dv_raw
							progress.set_value(bar, v)

				#print '@'*20
				if 'autoignition' in reactor:
					fld = os.path.join(dir_desk,'raw',
						'['+fuel_name.strip('[').strip(']')+'] + ['+oxid_name.strip('[').strip(']')+']')
					find_tau_ign(fld)
					print
					print 'find_tau_ign'
					print 'fld = '+str(fld)
					print

				#print dir_desk
				#print '@'*20
	#sys.exit()
	#fld = os.path.join(dir_desk,)




	return v















def run_graph(parent, progress, task):
	if task == 'GPS':
		obj = progress.w.label_GPS
	elif task == 'GPSA':
		obj = progress.w.label_GPSA

	obj.setText('building graph')

		
	dir_public = parent.project['dir_public']
	soln = parent.soln['detailed']

	list_train = []
	train_name = parent.train_name
	for db_name in parent.project['database'].keys():
		if parent.project['database'][db_name]['train']:
			list_train.append(db_name)

	traced_list = []
	for GPS_name in parent.project['GPS'].keys():
		if parent.project['GPS'][GPS_name]['checked']:
			GPS = parent.project['GPS'][GPS_name]
			es_name = GPS['es']
			es = parent.project['es'][es_name]
			for e in es['element'].keys():
				if es['element'][e]['traced']:
					if e not in traced_list:
						traced_list.append(e)

	if not bool(traced_list):
		msg = 'no traced element in selected GPS'
		QMessageBox.information(QWidget(),'',msg)
		return False

	v = 0.0
	dv_db = 100.0/len(list_train)
	for db_name in list_train:

		database = parent.project['database'][db_name]
		phi_list = database['phi']
		T0_list = database['T0']
		atm_list = database['atm']
		fuel_list = database['fuel']
		oxid_list = database['oxid']
		reactor = database['reactor']

		dv_raw = dv_db / (len(phi_list) * len(atm_list) * len(T0_list) *\
			len(fuel_list) * len(oxid_list))
		
		for fuel_name in fuel_list:
			for oxid_name in oxid_list:
				for phi in phi_list:
					for atm in atm_list:
						for T0 in T0_list:

							if progress.stop:
								progress.close()
								return False

							dir_de = os.path.join(dir_public,'detailed')
							dir_raw = cond2dir(dir_de, fuel_name, \
								oxid_name, phi, atm, T0, \
								reactor, parent.n_digit)							
							if 'DNS' in reactor:
								dir_raw = os.path.join(dir_raw, database['case'][0])

							raw = load_raw(os.path.join(dir_raw,'raw.npz'))