Accessing and reformatting output

The purpose of this tutorial is to provide an explanation of how to access the raw output of CatMAP. The plotting functions included in CatMAP are not meant to provide publication-quality figures by default, but rather to provide a quick way of analyzing the output and determining whether or not the solution is correct. If you want to create more complex plots, or more beautiful plots, then you will likely want to reformat the raw data and plot it with software of your choice (MATLAB, etc.) or at least know how to access the matplotlib figure object if you are brave enough to directly edit the figure in python.

The simplest approach to post-processing is to convert the CatMAP output into a data table and use methods you are familiar with. We will use the CO oxidation example from tutorial 3 to provide some concrete context. If you have run Tutorial 3 you should have the file "CO_oxidation.log" in the tutorial 3 directory. Assuming you finished the tutorial, this file should have outputs for the coverage, rate, production_rate, and rate_control. We will look at how to read in each of these and convert them to more conventional tables.

First, lets take a look at how output is stored natively by CatMAP. Create a script with the following commands in the tutorial 3 directory:

model = ReactionModel(setup_file='CO_oxidation.log')```

which reads the model (along with outputs) into the script. The model, after its run, will have "map" attributes for each output variable. The map is a python list structured like:

```[
[[x_0,y_0],[out_0_0, out_0_1, ... , out_0_n]],
[[x_1,y_1],[out_1_0, out_1_1, ... , out_1_n]],
...,
[[x_m,y_m],[out_m_0, out_m_1, ... , out_m_n]]
]```
where x and y represent descriptor values, and out_i_j represents the output vector at each point. For example, the coverage map on the CO oxidation model. Add the following to the script:

```for pt, cvgs in model.coverage_map:
    print 'descriptors:', pt
    print 'coverages', cvgs```

If you run this you will get a bunch of numbers like:
```descriptors [1.9289202008928572, 2.232142857142857]
coverages [mpf('2.2800190488939763e-5'), mpf('0.020278026143493406'), mpf('2.2031908471388691e-7'), mpf('1.4713059814195941e-5'), mpf('0.18181715627634512')]
descriptors [0.1428571428571428, 2.553571428571429]
coverages [mpf('2.7893704760253176e-27'), mpf('0.050000000000000003'), mpf('2.8850535476689032e-25'), mpf('1.4466885932549904e-12'), mpf('0.94999999999855326')]
descriptors [3.0, 1.267857142857143]
coverages [mpf('0.049999987577682677'), mpf('4.4864510401468476e-25'), mpf('0.69752183650972651'), mpf('4.653706806446857e-10'), mpf('1.789519988373698e-17')]
descriptors [2.1830357142857144, 2.8750000000000004]
coverages [mpf('6.4052886383458481e-12'), mpf('0.024816133631270126'), mpf('3.4777261296001614e-13'), mpf('6.8031802914834851e-9'), mpf('0.32310023705530582')]
...```

You should note that the coverages are "mpf" objects, which indicates that they are multiple precision. You probably also don't know which intermediate each coverage corresponds to. Try the following:

```labels = model.output_labels['coverage']
for pt,cvg in model.coverage_map:
    print 'descriptors',pt
    print 'intermediates',labels
    print 'coverages', [float(c) for c in cvg]```

which gives the slightly more readable output:

```descriptors [1.9289202008928572, 2.232142857142857]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [2.2800190488939762e-05, 0.020278026143493406, 2.203190847138869e-07, 1.471305981419594e-05, 0.1818171562763451]
descriptors [0.1428571428571428, 2.553571428571429]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [2.7893704760253176e-27, 0.049999999999999996, 2.885053547668903e-25, 1.4466885932549903e-12, 0.9499999999985532]
descriptors [3.0, 1.267857142857143]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [0.049999987577682675, 4.486451040146847e-25, 0.6975218365097264, 4.653706806446857e-10, 1.7895199883736977e-17]
descriptors [2.1830357142857144, 2.8750000000000004]
intermediates ('CO_s', 'O_s', 'CO_t', 'O2_t', 'O_t')
coverages [6.405288638345848e-12, 0.024816133631270124, 3.477726129600161e-13, 6.803180291483484e-09, 0.3231002370553058]
...```

Based on this, you can probably see how to create a text table containing coverage outputs. All the other outputs follow the same basic format; however, there are a few tricky situations when looking at other outputs. For example, the "labels" for reaction-specific quantities (rates, rate constants, etc.) are actually lists which need to be flattened into strings. Even more difficult are "matrix" outputs like rate control, where the output is a list of lists rather than a single list. To make life easier I have created the following script which should create a tab-separated text table from any output (.log) file created by CatMAP. Just place this script into the output directory, and run it with the name of the output of interest as its first argument.

```from catmap.model import ReactionModel
import sys

output_variable = sys.argv[1]

model = ReactionModel(setup_file='CO_oxidation.log')

if output_variable == 'rate_control':
    dim = 2
else:
    dim = 1

labels = model.output_labels[output_variable]

def flatten_2d(output):
    "Helper function for flattening rate_control output"
    flat = []
    for x in output:
        flat+= x
    return flat

#flatten rate_control labels
if output_variable == 'rate_control':
    flat_labels = []
    for i in labels[0]:
        for j in labels[1]:
            flat_labels.append('d'+i+'/d'+j)
    labels = flat_labels

#flatten elementary-step specific labels
if output_variable in ['rate','rate_constant','forward_rate_constant','reverse_rate_constant']:
    str_labels = []
    for label in labels:
        states = ['+'.join(s) for s in label]
        if len(states) == 2:
            new_label = '<->'.join(states)
        else:
            new_label = states[0]+'<->'+states[1]+'->'+states[2]
        str_labels.append(new_label)
    labels = str_labels

table = '\t'.join(list(['descriptor-'+d for d in model.descriptor_names])+list(labels))+'\n'

for pt, output in getattr(model,output_variable+'_map'):
    if dim == 2:
        output = flatten_2d(output)
    table += '\t'.join([str(float(i)) for i in pt+output])+'\n'

f = open(output_variable+'_table.txt','w')
f.write(table)
f.close()```