Use correct way to find edgerxninds and corerxninds

Tricky thing here is that RMS rearrange the reactions while constructing Phase object in order to vectorize kinetics, so the old way of getting rxninds is not correct.
Preallocate vectors to help speeding up.

Author: hwpang

src/EdgeAnalysis.jl

@@ -234,13 +234,27 @@ Precalculate important indices and maps for use in edge analysis
 function getkeyselectioninds(coreedgedomain::AbstractDomain,coreedgeinters,domain,inters)
     corespcsinds = 1:length(domain.phase.species)
     edgespcsinds = length(domain.phase.species)+1:length(coreedgedomain.phase.species)
-    corerxninds = 1:length(domain.phase.reactions)
-    edgerxninds = length(domain.phase.reactions)+1:length(coreedgedomain.phase.reactions)
-    reactantindices = coreedgedomain.rxnarray[1:3,:]
-    productindices = coreedgedomain.rxnarray[4:6,:]
+    corerxninds = zeros(Int64,length(domain.phase.reactions))
+    corerxncount = 1
+    edgerxninds = zeros(Int64,length(coreedgedomain.phase.reactions))
+    edgerxncount = 1
+    coretoedgerxnmap = Dict{Int64,Int64}()
+    for (j,rxn) in enumerate(coreedgedomain.phase.reactions)
+        coreind = findfirst(x->rxn.reactants==x.reactants && rxn.products==x.products && rxn.kinetics==x.kinetics,domain.phase.reactions)
+        if coreind === nothing
+            edgerxninds[edgerxncount] = j
+            edgerxncount += 1
+        else
+            coretoedgerxnmap[coreind] = j
+            corerxninds[corerxncount] = j
+            corerxncount += 1
+        end
+    end
+
+    @views @inbounds reactantindices = coreedgedomain.rxnarray[1:3,:]
+    @views @inbounds productindices = coreedgedomain.rxnarray[4:6,:]
     coretoedgespcmap = Dict([i=>findfirst(isequal(spc),coreedgedomain.phase.species) for (i,spc) in enumerate(domain.phase.species)])
-    coretoedgerxnmap = Dict([i=>findfirst(isequal(rxn),coreedgedomain.phase.reactions) for (i,rxn) in enumerate(domain.phase.reactions)])
-    for j = 3:length(domain.indexes)
+    @simd for j = 3:length(domain.indexes)
         coretoedgespcmap[domain.indexes[j]] = coreedgedomain.indexes[j]
     end
     return corespcsinds,corerxninds,edgespcsinds,edgerxninds,reactantindices,productindices,coretoedgespcmap,coretoedgerxnmap