analyze_polycyclic_reactions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "\n",
    "from rmgpy.chemkin import loadChemkinFile\n",
    "from rmgpy.data.kinetics.family import TemplateReaction\n",
    "from rmgpy.data.kinetics.library import LibraryReaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load chemkin files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "chem_path_1 = '/home/mjliu/Documents/Models/Hexylbenzene/e_run2/chem_annotated.inp'\n",
    "dict_path_1 = '/home/mjliu/Documents/Models/Hexylbenzene/e_run2/species_dictionary.txt'\n",
    "\n",
    "chem_path_2 = '/home/mjliu/Documents/Models/Hexylbenzene/e_run3/chem_annotated.inp'\n",
    "dict_path_2 = '/home/mjliu/Documents/Models/Hexylbenzene/e_run3/species_dictionary.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "spc_list_1, rxn_list_1 = loadChemkinFile(chem_path_1, dict_path_1)\n",
    "spc_list_2, rxn_list_2 = loadChemkinFile(chem_path_2, dict_path_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print len(spc_list_1), len(rxn_list_1)\n",
    "print len(spc_list_2), len(rxn_list_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "common = 0\n",
    "for spc1 in spc_list_1:\n",
    "    for spc2 in spc_list_2:\n",
    "        if spc1.isIsomorphic(spc2):\n",
    "            common += 1\n",
    "            break\n",
    "print common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "common = 0\n",
    "rxn_list_2_copy = list(rxn_list_2)\n",
    "for rxn1 in rxn_list_1:\n",
    "    i = 0\n",
    "    while i < len(rxn_list_2_copy):\n",
    "        if rxn1.isIsomorphic(rxn_list_2_copy[i]):\n",
    "            common += 1\n",
    "            del rxn_list_2_copy[i]\n",
    "            break\n",
    "        else:\n",
    "            i += 1\n",
    "print common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "common = 0\n",
    "rxn_list_1_copy = list(rxn_list_1)\n",
    "for rxn2 in rxn_list_2:\n",
    "    i = 0\n",
    "    while i < len(rxn_list_1_copy):\n",
    "        if rxn2.isIsomorphic(rxn_list_1_copy[i]):\n",
    "            common += 1\n",
    "            del rxn_list_1_copy[i]\n",
    "            break\n",
    "        else:\n",
    "            i += 1\n",
    "print common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "len(rxn_list_2_copy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "common"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyze reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def analyze_reaction_list(rxn_list):\n",
    "    families = {}\n",
    "    library_rxns = []\n",
    "    for rxn in rxn_list:\n",
    "        if isinstance(rxn, TemplateReaction):\n",
    "            if rxn.family not in families:\n",
    "                families[rxn.family] = [rxn]\n",
    "            else:\n",
    "                families[rxn.family].append(rxn)\n",
    "        elif isinstance(rxn, LibraryReaction):\n",
    "            library_rxns.append(rxn)\n",
    "            \n",
    "    family_counts = {}\n",
    "    for family, rxns in families.iteritems():\n",
    "        family_counts[family] = len(rxns)\n",
    "    \n",
    "    return families, family_counts, library_rxns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "families_1, family_counts_1, library_rxns_1 = analyze_reaction_list(rxn_list_1)\n",
    "print family_counts_1, len(library_rxns_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "families_2, family_counts_2, library_rxns_2 = analyze_reaction_list(rxn_list_2)\n",
    "print family_counts_2, len(library_rxns_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get percentages for kinetics families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "counts_1 = [(family, count) for family, count in family_counts_1.iteritems()]\n",
    "counts_1.sort(key=lambda x: x[1], reverse=True)\n",
    "for family, c in counts_1:\n",
    "    print '{0} {1}'.format(family, c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "counts_1 = [(family, count) for family, count in family_counts_2.iteritems()]\n",
    "counts_1.sort(key=lambda x: x[1], reverse=True)\n",
    "for family, c in counts_1:\n",
    "    print '{0} {1}'.format(family, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Count all duplicates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def count_duplicates(families):\n",
    "    output = {}\n",
    "    for family, reactions in families.iteritems():\n",
    "        sorted_reactions = []\n",
    "        for reaction in reactions:\n",
    "            if not sorted_reactions:\n",
    "                sorted_reactions.append([reaction])\n",
    "            else:\n",
    "                found = False\n",
    "                for rxn_list in sorted_reactions:\n",
    "                    for rxn in rxn_list:\n",
    "                        if reaction.isIsomorphic(rxn):\n",
    "                            rxn_list.append(reaction)\n",
    "                            found = True\n",
    "                            break\n",
    "                    if found:\n",
    "                        break\n",
    "                else:\n",
    "                    sorted_reactions.append([reaction])\n",
    "        output[family] = sorted_reactions\n",
    "\n",
    "    all_stats = {}\n",
    "    for family, sorted_reactions in output.iteritems():\n",
    "        stats = {}\n",
    "        stats['rxn_num'] = len(sorted_reactions)\n",
    "        stats['total_rxns'] = 0\n",
    "        for rxn_list in sorted_reactions:\n",
    "            stats['total_rxns'] += len(rxn_list)\n",
    "        stats['rxns'] = sorted_reactions\n",
    "        \n",
    "        all_stats[family] = stats\n",
    "    \n",
    "    return all_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stats_1 = count_duplicates(families_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_total = 0\n",
    "total_rxn = 0\n",
    "for family, stats in stats_1.iteritems():\n",
    "    print family\n",
    "    print stats['total_rxns'], stats['rxn_num']\n",
    "    total_total += stats['total_rxns']\n",
    "    total_rxn += stats['rxn_num']\n",
    "print total_total, total_rxn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stats_2 = count_duplicates(families_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_total = 0\n",
    "total_rxn = 0\n",
    "for family, stats in stats_2.iteritems():\n",
    "    print family\n",
    "    print stats['total_rxns'], stats['rxn_num']\n",
    "    total_total += stats['total_rxns']\n",
    "    total_rxn += stats['rxn_num']\n",
    "print total_total, total_rxn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Count duplicates in intra add families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def analyze_intra_add_rxns(families):\n",
    "    target = ['Intra_R_Add_Endocyclic', 'Intra_R_Add_Exocyclic', 'Intra_R_Add_Polycyclic']\n",
    "    output = {}\n",
    "    for family, reactions in families.iteritems():\n",
    "        if family in target:\n",
    "            sorted_reactions = []\n",
    "            for reaction in reactions:\n",
    "                if not sorted_reactions:\n",
    "                    sorted_reactions.append([reaction])\n",
    "                else:\n",
    "                    found = False\n",
    "                    for rxn_list in sorted_reactions:\n",
    "                        for rxn in rxn_list:\n",
    "                            if reaction.isIsomorphic(rxn):\n",
    "                                rxn_list.append(reaction)\n",
    "                                found = True\n",
    "                                break\n",
    "                        if found:\n",
    "                            break\n",
    "                    else:\n",
    "                        sorted_reactions.append([reaction])\n",
    "            output[family] = sorted_reactions\n",
    "    \n",
    "    cross_duplicates = []\n",
    "    \n",
    "    if 'Intra_R_Add_Endocyclic' in output and 'Intra_R_Add_Exocyclic' in output:\n",
    "        for rxn_list_1 in output['Intra_R_Add_Endocyclic']:\n",
    "            for rxn_list_2 in output['Intra_R_Add_Exocyclic']:\n",
    "                if rxn_list_1[0].isIsomorphic(rxn_list_2[0]):\n",
    "                    cross_duplicates.append(rxn_list_1 + rxn_list_2)\n",
    "                    break\n",
    "    \n",
    "    return output, cross_duplicates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output_1, cross_1 = analyze_intra_add_rxns(families_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_num = 0\n",
    "for rxn_list in cross_1:\n",
    "    display(rxn_list[0])\n",
    "    print len(rxn_list)\n",
    "    total_num += len(rxn_list)\n",
    "print len(cross_1), total_num\n",
    "\n",
    "# display(cross_1[22][0])\n",
    "# for rxn in cross_1[22]:\n",
    "#     print rxn.kinetics.comment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for family, sorted_reactions in output_1.iteritems():\n",
    "    print family, len(sorted_reactions)\n",
    "    total_num = 0\n",
    "    for rxn_list in sorted_reactions:\n",
    "        total_num += len(rxn_list)\n",
    "        display(rxn_list[0])\n",
    "        print len(rxn_list)\n",
    "    print 'Total reactions = {}'.format(total_num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output_2, cross_2 = analyze_intra_add_rxns(families_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for rxn_list in cross_2:\n",
    "    display(rxn_list[0])\n",
    "    print len(rxn_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for family, sorted_reactions in output_2.iteritems():\n",
    "    print family\n",
    "    for rxn_list in sorted_reactions:\n",
    "        display(rxn_list[0])\n",
    "        print len(rxn_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "display(output_2['Intra_R_Add_Polycyclic'][5][0])\n",
    "for rxn in output_2['Intra_R_Add_Polycyclic'][5]:\n",
    "    print rxn.kinetics.comment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_num = 0\n",
    "for rxn_list in output_2['Intra_R_Add_Polycyclic']:\n",
    "    total_num += len(rxn_list)\n",
    "print len(output_2['Intra_R_Add_Polycyclic']), total_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find reactions which are no longer generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "old_intra_add = [rxn_list[0] for family in output_1.itervalues() for rxn_list in family]\n",
    "new_intra_add = [rxn_list[0] for family in output_2.itervalues() for rxn_list in family]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diff_reactions = []\n",
    "for rxn in old_intra_add:\n",
    "    for rxn2 in new_intra_add:\n",
    "        if rxn.isIsomorphic(rxn2):\n",
    "            break\n",
    "    else:\n",
    "        diff_reactions.append(rxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for rxn in diff_reactions:\n",
    "    display(rxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diff_reactions = []\n",
    "for rxn in new_intra_add:\n",
    "    for rxn2 in old_intra_add:\n",
    "        if rxn.isIsomorphic(rxn2):\n",
    "            break\n",
    "    else:\n",
    "        diff_reactions.append(rxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for rxn in diff_reactions:\n",
    "    display(rxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:rmg_env]",
   "language": "python",
   "name": "conda-env-rmg_env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}