diff --git a/day3/Interp_ds107_sub012_t1r2.nii b/day3/Interp_ds107_sub012_t1r2.nii new file mode 100644 index 0000000..0ebb8be Binary files /dev/null and b/day3/Interp_ds107_sub012_t1r2.nii differ diff --git a/day3/slice_time_image.py b/day3/slice_time_image.py index ff73e3d..2423310 100644 --- a/day3/slice_time_image.py +++ b/day3/slice_time_image.py @@ -21,10 +21,12 @@ """ # Add any extra imports here +import os import sys - +import numpy as np +import matplotlib.pyplot as plt import nibabel as nib - +from scipy.interpolate import InterpolatedUnivariateSpline as IUS def slice_time_image(img, slice_times, TR): """Run slice-timing correction on nibabel image 'img'. @@ -47,8 +49,32 @@ def slice_time_image(img, slice_times, TR): A new copy of the input image with slice-time interpolation applied """ # Get the image data as an array; + data = img.get_data() # Make a new empty array "interp_data" to hold the interpolated data; + interp_data = np.zeros(data.shape) + # Do the interpolation; + time_offset = slice_times/slice_times.shape * TR + print "Time Offsets: %s" % time_offset + slice_0_times = TR * np.arange(data.shape[-1]) + + # For each z coordinate (slice index): + for z in range(data.shape[2]): + print "Starting Slice %s of %s" % (z+1, data.shape[2]) + slice_z_times = slice_0_times + time_offset[z] + # For each x coordinate: + for x in range(data.shape[0]): + # For each y coordinate: + for y in range(data.shape[1]): + # extract the time series at this x, y, z coordinate; + current_timecourse = data[x,y,z,:] + # make a linear interpolator object with the `slice_z_times` and the extracted time series; + current_interp = IUS(slice_z_times, current_timecourse, k=1) + # resample this interpolator at the slice 0 times; + current_timecourse_estimate = current_interp(slice_0_times) + # put this new resampled time series into the new data at the same position + interp_data[x,y,z,:] = current_timecourse_estimate + # This is how to make a new image with the interpolated data new_img = nib.Nifti1Image(interp_data, img.affine, img.header) return new_img @@ -71,7 +97,13 @@ def slice_time_file(fname, slice_times, TR): TR: double Repetition time in seconds. """ - # Hint: use os.path.split and os.path.join to make the new filename + img = nib.load(fname) + interp_img = slice_time_image(img, slice_times, TR) + + path, tail = os.path.split(fname) + new_tail = 'Interp_' + tail # I decided to use a different new filename + # I thought it would provide more information than just adding 'a' + new_fname = os.path.join(path,new_tail) nib.save(interp_img, new_fname) @@ -82,12 +114,32 @@ def main(): fname = sys.argv[1] # Assume the TR TR = 2.0 + # Assume the slices were acquired even slices first, inferior to # superior, then odd slices, inferior to superior, where the most inferior # slice is index 0 and 0 is an even number. What are the slice acquisition # times in seconds, where the first value is the acquisition time of slice # 0, the second is acquisition time of slice 1, etc? - slice_times = ? + + # have to load in image to know the number of slices there are + img = nib.load(fname) + num_slices = img.shape[2] + print "Number of slices: %s" % num_slices + + # use ceil in the next two lines in case of an odd number of slices + first_half = np.arange(0, np.ceil(num_slices/2.)) + second_half = np.arange(np.ceil(num_slices/2.), num_slices) + slice_times = np.array([]) + + for i in first_half: + slice_times = np.append(slice_times,first_half[i]) + if i < len(second_half): + # if there are an odd number of slices + # len(second_half)]" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXl8W9Wd9/85sixZlmRbtuXdiRPHdkggCQkkYSk4tGEp\nBCi004FhLfMrU1qmhZm2v87zzI/QoQvTB55p6QwtS1gLPPQptA3blKVmCRBDiLNjO04c27LlTbus\nzfL5/XF1ru/VvVos27KcnPfrlVfso6t7j+R7z/d8P9/v+R5CKQWHw+FwTk00C90BDofD4Swc3Ahw\nOBzOKQw3AhwOh3MKw40Ah8PhnMJwI8DhcDinMNwIcDgczilMUiNACNlBCBkmhByIa7+TEHKEEHKQ\nEHJ/rG0rIeRTQsj+2P9bJMdvIIQcIIR0E0J+OT8fhcPhcDgzJZUn8ASAS6UNscH9SgBrKKWnA/hf\nsZdGAVxBKV0D4GYAz0je9jCA2yilTQCaCCGyc3I4HA5nYUhqBCil7wNwxjV/C8DPKKWR2DGjsf87\nKKX22DGHARgIIfmEkGoAZkppe+y1pwFcPVcfgMPhcDiZk0lMoAnABYSQjwkhbYSQs1SOuRbAnpih\nqAUwIHnNFmvjcDgczgKjzfA9FkrpZkLI2QBeBLCcvUgIWQ3g5wC2zk0XORwOhzNfZGIEBgC8BACU\n0k8IIVOEkDJK6TghpC722o2U0uOx420A6iTvr4u1KSCE8EJGHA6HkwGUUpLJ+zKRg/4I4CIAIIQ0\nA9DFDEAJgFcB/JBS+pGkY0MAPISQTYQQAuDG2DlUoZTm1L977rlnwfvA+3Ry9Yv3ifdprv/NhlQp\nos8D+BBAMyGknxByK4AdAJbH0kafB3BT7PDvAGgEcA8hZG/sX3nstTsAPAagG8BRSukbs+o1h8Ph\ncOaEpHIQpfS6BC/dqHLsfQDuS3CePQDOmHHvOBwOhzOv8BXDKWhtbV3oLijgfUqfXOwX71N68D5l\nBzJbPWkuIYTQXOoPh8PhLAYIIaBZDAxzOPPO3qG9GPQOLnQ3OJyTHm4EODnJL3f/Eq90vbLQ3eBw\nTnq4EeDkJOFoGOFoeKG7weGc9HAjwMlJwtEwQpOhhe4Gh3PSw40AJyfhngCHkx24EeDkJNwIcDjZ\ngRsBTk4SjoYRinI5iMOZb7gR4OQk3BPgcLIDNwKcnIQbAQ4nO3AjwMlJeHYQh5MduBHg5CThaBjh\nKe4JcDjzDTcCnJyEewIcTnbgRoCTk/CYAIeTHbgR4OQk3AhwONkh1c5iOwghw7FdxKTtdxJCjhBC\nDhJC7pe0/4gQ0k0I+ZwQcrGkfQMh5EDstV/O/cfgnGzwdQIcTnZI5Qk8AeBSaQMhZAuAKwGsoZSe\nDuB/xdpXAfg6gFWx9/xXbE9hAHgYwG2U0iYATYQQ2Tk5nHi4J8DhZIekRoBS+j4AZ1zztwD8jFIa\niR0zGmu/CsDzlNIIpbQXwFEAmwgh1QDMlNL22HFPA7h6jvrPOUnhgWEOJztkEhNoAnABIeRjQkgb\nIeSsWHsNgAHJcQMAalXabbF2DkcVSikiUxHuCXA4WSDpRvNJ3mOhlG4mhJwN4EUAy+eqQ9u3bxd/\nbm1tPSn39OQkJzIVAQBuBDicBLS1taGtrW1OzpWJERgA8BIAUEo/IYRMEULKIczw6yXH1cWOtcV+\nlrbbEp1cagQ4pyZs8OeBYQ5HnfgJ8r333pvxuTKRg/4I4CIAIIQ0A9BRSscA/BnA3xJCdISQZRBk\no3ZKqR2AhxCyKRYovjF2Dg5HFWYEuCfA4cw/ST0BQsjzAC4EUEYI6Qfw/wHYAWBHLG00DOAmAKCU\nHiaEvAjgMIBJAHdQSmnsVHcAeBKAAcBrlNI35uGzcE4SRE+AB4Y5nHmHTI/TCw8hhOZSfzgLQ5+7\nDy2/boEx34ixH4wtdHc4nJyHEAJKKUl9pBK+YpiTc4SjYZh1Zi4HcThZgBsBTs4RjoZh0pl4YJjD\nyQLcCHByDmYEwtEwuDzI4cwv3Ahwco5wNIwCbQG0Gq24ZoDD4cwP3Ahwco5wNAxdng66PB2PC3A4\n8ww3ApycgxkBfZ6eGwEOZ57hRoCTc0g9Ab5WgMOZX7gR4OQcoieg5Z4AhzPfcCPAyTl4TIDDyR7c\nCHByDpkcxNcKcDjzCjcCnJyDB4Y5nOzBjQAn5+CBYQ4ne3AjwMk5eGCYw8ke3Ahwco5wNAx9np4H\nhjmcLMCNACfn4IFhDid7cCPAyTl4YJjDyR5JjQAhZAchZDi2ixhr204IGSCE7I39uyzWXkAIeZ4Q\nsp8QcpgQ8v9K3rOBEHKAENJNCPnl/H0czskADwxzONkjlSfwBIBL49oogAcppWfG/r0ea/9bAKCU\nrgGwAcDthJAlsdceBnAbpbQJQBMhJP6cHI4IDwxzONkjqRGglL4PwKnykto2ZkMAjISQPABGCPsP\newgh1QDMlNL22HFPA7g68y5zTnZET0CzcIHhJzue5HsZcE4JMo0J3EkI2UcIeZwQUgIAlNL/BuCB\nYAx6AfyCUuoCUAtgQPJeW6yNw1FloQPDgUgAt/7pVvjCvqxfm8PJNtoM3vMwgB/Hfv43AA8AuI0Q\ncgMAA4BqAKUA3ieEvD3Tk2/fvl38ubW1Fa2trRl0kbOYWWg5yBFwAABG/CMw681Zvz6Hk4q2tja0\ntbXNyblmbAQopSPsZ0LIYwB2xn49F8DLlNIogFFCyC4IsYEPANRJTlEHwRtQRWoEOKcmCx0YHg+M\nAwBGJ0bRWNqY9etzOKmInyDfe++9GZ9rxnJQTONnfAUAyxz6HMBFsWOMADYD+JxSaocQG9hECCEA\nbgTwx4x7zDnpWegUUeYJjPpHs35tDifbJPUECCHPA7gQQDkhpB/APQBaCSHrIGQJHQdwe+zw3wJ4\nPJZOqgGwg1J6MPbaHQCehCAXvUYpfWOuPwjn5EHqCUxEJrJ+fakcxOGc7CQ1ApTS61SadyQ4NgTg\nhgSv7QFwxox7xzklkRoBZ1AtOW1+GZ+YloM4nJOdTALDHE5SpugUpugUtJrMbq9cCAzna/IXRA4a\n9g0jSqMoLyyHLk+X9etzTj142QjOnPPSkZdw+87bUx+YgIUODDsCDqwoXYGRiezKQR8PfIwl/7EE\nax5eg+v/cH1Wr805deFGgDPn2Dw22P32jN8vCwxPZe4J/N1Lf5fRbH48MI6W8pasewI2jw1fbvoy\ner/Xi7eOvYVh33BWr885NeFGgDPnuIIueEKejN8/V3sM7+zciXZbe+oD43AEHFhZtjLrMYHxwDjK\nDGUw6Uy4euXV+N2B32X1+pxTE24EOLPio/6PFOUVnEHnnBmBTOWgcDQMb9iLDnuHrN0X9mH/8P6k\n73UEHFhZvjLr2UHjE4IRAIBb1t2CJzqeWDSlK/bZ9yEQCSx0NzgZwI0AJ2MopWh9qhVjE2Oy9rky\nArMJDLMMn45huRF4pesV3PHqHUnf6wg4RDkom4PweGAcZYWCEbhg6QXwhX3Ya9+btevPhu++8V28\ne+Ldhe4GJwO4EeBkjDPoRDgahjvklrXnghw0HhiHPk+v8AT63H3oHO9M+d66ojrkafKyWj9obGIM\n5YXlAAAN0eCaldfg9e7XU7wrN/CFffCH/QvdDU4GcCPAyRi7Twj+xg/4zoAT7qA741m0NDCcaQG5\n8YlxrK9ej0HvILwhr9je5+7D2MSYuCBMDUfAgVJDKSqMFVmVhFhMgFFeWA5v2JvkHbmDP+JfkIV9\nnNnDjQAnYxIagaATURpFYDIzjXiuPIFKUyVWW1fLYgD9nn4AQPd4t+r7mK5dmF8Ia6E1q8Hh8Ylp\nOQgAjDrjopld+8N++COLo68cOdwIcDImkRFwBV2q7emSKjAciUbwbq9Sfz48eljsEwuyrqtaJ5OE\n+tx9aC5rTigJjQfGUWooBQBYjdasponGewLGfOOiGVj9Ef+iMVgcOdwIcDKGDbjuoDwm4Aw4UWWq\nmrURSBQY7rB34NoXr1XITf++69/x7P5nAQj6upoR6Hf340vLvoSu8S7VazMpCED25aCJcTEmAAje\nSCqJhVIKZ2C6tEYgEliQLB1/mMtBixVuBDgZo+YJRKIRBCeDqDHXKIxDuqSSg5xBJ8YD4+hx9sja\nXUEXbB6hSjnLtDmz6kwxQ2giMgF/xI/zlpyX0BOQGoFsykHRqShcQRcsBovYZtSl9gQ+HfwUVzx/\nhfj7Ax89gH/f9e/z1k81olNRhKKhReO1cOQsaiPwz3/5Z7x17K2F7sYpi91nV8z4XUEXSgpKUKwv\nzsgTiE5FQUGRp8lLGBhmM9/dA7tl7a6gCwNeYRO78YAwqz6j8gwcHj2MSDSCfnc/6orqsLJ8ZUJP\nQJqrby3MnhzkCrpg1ptl9ZaM+caUs+vRiVFZkNsRcGDYn92Vxmzw53LQ4mRRG4F3jr+DA8MHUh/I\nmRfsPjuay5plg70z6ITFYEGRvigjIxCZioiF0xJ5Aq6gCxqiwW6b0giInsDE9OrbZSXLsH94P/rc\nfVhSvARNpU3oHu/GFJ1SnFshB81x/aAde3fgyOgRRXt8PAAQ5KBUA6s76JalsfrCvnmtvNo51olH\n9zwqa2N95HLQ4mTRGgFKKbrGuzDoHRTbXu16ddGssDwZsPvsaC5tlq0TcAacgidQkJknwKQgAAkD\nw86gExtrN6obAa9cDgKATbWb0G5rR7+nH/VF9TDrzbAYLBjwDCjOLZODEgSGD44cxDHnsRl/NgB4\n4eALeL/vfUV7fGYQkJ4c5Al5FEYgWfrrbNkztAePfPaIrE30BLgctChZtEZg0DsIf8SPId8QAEFG\nuOqFq/hGIPPM2MSYODtP6AkUWFCky8wTkBqBRIFhZ8CJLy37Eg6OHJQZCVfQhUHvIKbolBgYBoBN\ndZuw27Zb9AQAoLmsWVUScgQccjlIJSbwHx//B37z6W9m/NkA4fsZ8g4p2qULxRjxcpDdZ1dMctwh\nwRNg7b6wTxYonmt8YR86xzpl/WCeQCojMBGZyDhOxJk/khoBQsgOQshwbLcw1radEDJACNkb+3ep\n5LU1hJCPCCEHCSH7CSG6WPsGQsgBQkg3IeSXc9HxrvEuaDVa0RMY8g0hSqOLZnHNYuXbr30bT3U8\nhUg0AmfQiRWlKxQxASYHxa8kTod4TyCRHFRbVIum0ibsG94HQNjDwBf2oVhfjBH/iGxmval22gjU\nF9UDAFrKWtA5pgwOS1NEq83VMk+TYfPaFCuR04UZKrXrppKDtj2/TVEQzx10Y3JqUvye5tsT8If9\n8Ia9srgDG/xTyUG/+fQ32N62fd76xsmMVJ7AEwAujWujAB6klJ4Z+/cGABBCtACeAfBNSunpELal\nnIy952EAt1FKmwA0SQ1HpnSNd+HsmrNFT6DP3Qcg89x0TnrYfXZ8NPARRidGUWYog8VgkXsCASdK\n9CUZxwTUjIBagTpLgUUY3GPBYU/IA5POhCXFS9Dn7oMr6BIH89UVq9Hv7seBkQOiJ9BU2oRuh3LB\nmFQOqjHXwBV0KQa3Ac8AOuwdqtLj2MQYLnzywoSfzxlwYtCnYgQmlEYgXg5yBByi3MVghpZJQvMd\nE2DXkRpQf9gvrGlIEb8Y8AwsyE5xM+Gh3Q/ht5/+dqG7kVWSGgFK6fsA1P5qRKXtYgD7KaUHYu91\nUkqnYhvTmymlbArzNICrZ9FnAEDneCdaG1rFWVW/W1gJyo3A/DI2MYbdtt1iZlD8jH+2gWGpEdAQ\nDfI0eZicmpQd4wwKcYdNdZvQPijcViwrqbaoFodHD8sybbQaLdZXr8eng5+KRqDGXCOmuEqRGgEN\n0aChpEGh/9s8NvjCPnECImXYN4w9g3tUPxulFK6gS1UOksYwGGydADM23pBX0Wf2HcuMQMCpGvSe\nC9h1pFKaP+KH1WhNKQfZffac99R7nD2qk4OTmUxjAncSQvYRQh4nhJTE2poAUELIG4SQPYSQ78fa\nawFII3C2WNus6BrvwsbajQhHw/CH/aInIK0To8YHfR+oBgQ56TE2MYbOsU50jXeJRkAhBxVY5iQw\nDMSCw3FpokxyWlayTPy7MyNQZ67DPvs+xax6U+0mAEB9sSAHVZmqFANqOBqGzWuTDcaNlkb0OKbX\nI0xEJhCYDODc+nNVJSFPyAN/xI9INKJ4zRf2IUqj6nJQ3EIxQDBeWo1W/PyekEfRZzVPgILO22TI\nH/GjxlwjW2fhD/tRYaxIKQfZffaUz+dC4w/7xRXvpwqZbAL7MIAfx37+NwAPALgNQD6A8wGcBSAA\n4G1CyB4AMxKGt2/fLv7c2tqK1tZW1eM6xzvRUtaCalM1hnxDYk2YVDf/PW334Ourv45vbvjmTLo1\nZ/S6etFQ0rAg154tU3QKjoADZ9WchVe6XkGVqUqxHsAZcKLR0jgnngAAYXexuLiAMyDIQSadSUwE\nkHoC7xx/RzGr3lS3SXwPoDQC+4f344aXbsBp5adhZflKsb3R0ihblGbz2FBrrhUWodk78OWmL8uu\nwz6zM+hEhbFC3u+gE5XGSoxOjCI6FUWeJg/HncexzLIMY4ExheECIMosBAShaEhpBIJKI2DSmeAI\nOFBSUKI432zxhX1YX71e4QlUGCtwwnUi6XvtPvu89Gku8Uf8C7Kv9Uxpa2tDW1vbnJxrxkaAUiqm\n3xBCHgOwM/ZrP4D3KKWO2GuvAVgP4FkAdZJT1EHwBlSRGoFEhKNh9Lv70VjaiBpzDQa9g+hz98Fa\naE068FBK0WHvwCWNl6S8xnwwEZlA80PN8P7IC71WvyB9mA2uoAsmnQnn1Z+HHR07cPuG2xWDvVQO\nmm1gGFAPDjM5iBAibsEoegJFddg3vE+c+TNaG1rxD2f9g/h7vBH43hvfw81rb8bd59wNQqbVzuWW\n5bIBb8AzgNqiWqyrWoc/df5J0X/2XTgCDqURCDhhNVpBQTHiH0GpoRQtv27Bnm/uUU0RBSSSEARJ\nSM0TKNYXy4xAU1mTkCFkUZxu1vgjfqyvWo//c+j/TLeF/bAWpicH5Wny5r5Tc8hiqYEUP0G+9957\nMz7XjOWgmMbP+AoAljn0FwBnEEIMsSDxhQAOUUrtADyEkE1EeLpuBPDHjHsM4JjzGOqK6qDL04kZ\nHP2efqyuWJ1UcxzwDMARcMxrCl0yRvwjiExFRAljscHSGDfVbYIr6EKVqUo0ZsHJIIBpOWjOPAGt\nXpYGSimFO+hGSUEJSg2l8Ia9CEfD056AuRaOgEMhrZQXluOnX/yp+HtJQQkCk9N1dnpdvbhq5VUy\nAwAAjaVxnoDXhrqiOkVNIga7/9TuMfbd1JhrMOQbQtd4FyJTETy17ynV7CBgOjjMZBS1mEBtUS18\nYZ8QRAdFlalq3jKEfGEf1latRa+rV5S8/BHBCEjjF/GEJkNwBp1cDspBUqWIPg/gQwAthJB+Qsg3\nANwfS//cB2GgvwsQAsEAHgTwCYC9APZQStmOGHcAeAxAN4CjLKMoU7rGu9BS3gIAqDHVYMg7hD53\nH1ZbVycdeNhDu1AZCmzhUa+rd0GuP1uYEdhYuxGAMJsGIBvw2Sx9roxAvCfgDXthyDcgPy8fGqJB\neWE5Rv2jMk8AgOqAKoUQgipTFYb9w5iiU+LgHk8iOailvAUDngHFoCb1BOJhXlK1SZi4HBo9hLWV\na/Hs/mcx7BtW9QSYHOQJeVCkL1KVg2rNghFgUlCpoXTe7nFf2IcyQxlqzDXifewP+1FcUCyLX8Qz\n7B+GQWvI+cCwP3LqGYGkchCl9DqV5h1Jjv8dAMXu2JTSPQDOmHHvEtA51onm0mYAQpZHj7MH3pAX\nK0pXJF3J2WHvQENJw8IZgYnFbwSshVYsLV6KCmOFwghUGCsEvd5gybh2UKrAMFuRzKg0VmLYPyyL\nCQBQHVDjYZKQPk+PkoISFGgLFMcsswjB58mpSWg1Wgx4BrCidAW0Gi1WVwh7FZy35DzxeKkxjIf1\nPV+Tj0HvIAY8A7iy5Uq8ffxtfNj/oarhklYSbSptwoGRA6CUghAieEUhN2rMNTIjYCmwiJ7Icwee\nQ4e9AwatAf964b/KahNlgj/sh0lnEstxN5U1wR/xo8pUJRoste/R7rOjsbRRdW1GLsE9gUXCcddx\nLLcsByAs6Gm3taOuqC7lwNMx3IHWhtYF+yOzIOZiNQKj/lGUF5aDEIJnv/KsqLtLv3cmeZj1ZnhC\nHlV54MjokYSyQarAMDs/o9JUiRH/iGgEivRFMOlMKT0BYNoI9Hv6xdTReAq0BagwVogpyFKPYbll\nuZiQwGCegZonwPpebarGkHcIh0YPYbV1NW5ZewsK8wthyDco3sPkIGZkDVqDaGACkwHkkTyUGcoU\nngC7/j1t9yCP5OGh9odUs5JmCruGdMW1P+yHUWdMWvra7rNjafFSAFAtBTIXRKIRWSZXJjBP4FQq\nP7MojYD0oa0x16DD3oH64noU6YuSupsd9g5sadgi02vn4o+d7jlG/aOwFlrR6+6d9TUXAmlpg62N\nW8VBq0hfBHfQjSk6BU/Ig+KCYujydNBqtKq7i53z+DnY1b9L9Rqp5CAmqTAqjBUY9k17AgBQV1SX\nnidgFIyAdCWxGo2WRtHDZIFhACgtKFUM9lKPKB62yI0lMxwaOYRV1lX4+ulfx12b71K9Nptde8Ne\nFOmLZAFtd9CN4oJimHQmpScQdGJyahL97n5sb92OuqK6tGNhye5nf0QY8JdbluO48/h0W74xaa0j\ntq7ErDfP277Nbb1tuPVPt87qHP6wH1EaPaXqIC1KI9Dn7hPzvatN1YhMRbCkeElSHdoddGPYN4yz\na86WuerNv26eVaC4c6wT5+44N61jRydGsbF246L1BNTq2wDTcpAn5EFhfqEoOaj9PYKTQbhDbjzV\n8ZTqNeKNQElBiVgZFFCXg6SeAACcX38+msuaU34e0RNwJ/YEAHlcQOoJlBWWYXxiXHasJ+zB0uKl\n6jGBWN9rzDXodfei19WL5rJmFOmLcN9F96lem82uWUxAagQ8IQ+K9UojwDyBfnc/KowV0Gv1sBgs\nacugl/3uMkWZbga7RoWxQpQ3/RE/CvMLk64aZkbApDPNW1zAGXTOOiDODNqpJAktSiMgfWhrzDUA\nIFaHTGQE9g/vxxmVZ6CssEz8A0eiERx1HJ1V0bkjY0fSzvYZ8Y/g7JqzF4UROOE6gb//898jOhUV\n28YCyY1A/KYoakZg1D8Kk86E/3vk/6pKB+FoGDrNtBG49rRr8dzB58Tf4+WgCmOFLCYAAI9e+SjW\nVa1L+RnZgJrSEygVFoxNTk1i1D8qxkKksgvDE/JgaclS1QHXFRK+n2pzNXb17cIyy7KUqcJsi0k1\nI+AOuUX5S+YJxAb8HmcPGksbAUAWJ0jFMecxfD72uaKdUgpf2AdjvlFWXC9dOajKVAWzzjxvGUKu\noGtW8T5KKSYiE6gtquVGIJfxhX0ITgZFzZcF9JgnkOgG22vfi3WV62ApsIiaH7uJZ3Pj9Lp60364\nRidGsbZqLcYmxtLSRR/d86jioRryDuH3h36fUV9nwpGxI3h87+P48bs/FttYYDgeFhOIn6WrxWhG\nJ0bRaGnE5rrNePnIy4pzxXsCf7P6b/D2sbcxNjEGYFpSYah5AukiGgFPX0pP4MjYEdh9dpQXloue\nTqmhFI6g3Ah4Q140FDck9ASYHOSP+LHKuiplH9lm856QB2adObEcFFF6Aj2OHjRaYkbAYEl7YBvx\njyhiHQAQioag1WiRn5cPq9EqTp5mKgdl4gm8f+J97LPvS3qMO+ielVcfnAxCq9GivLCcG4Fcpt/d\nj/riejGfmxCCalM16ovqk8pBL3/+MrY2bkV+Xj70eXr4wj5xodFsbpxeVy8Ck4G0BvVR/yiqTdWo\nK6pL6T2Eo2F85/XvKOrQvH70dfx8189TXivTKpcMT8iDLyz5Ah7b+xjePvY2gOnAcDxsYdiRsSNo\nKm2StXtCHvS6ekVjNuofRYWxAresvQVPdDyhOFe8ESguKMYVzVfg+QPPA1DKQWqeQLpI5SAmL6qx\ntXEr9tr34pE9j8jSSMsMZTPyBFj6bKWxEgQEq62rU/aRza69IZWYQGyhmMwTyJ/ODupxThuBEn1J\nWpOd0GQI7pBb9f5kRgYQvneW8jwRmYBRZ0wqBw37h2flCTy7/1n87oAi8VCGO+RWfRaPOo6Kk4hk\nMGNWUlCSNSPQbmvPeG+KuWLRGQE11/2uzXfh7NqzExqB487jODhyEFc0C3uxlhQID4S05EA6qC0n\nZ9JOonNIq2CO+EdgNVrRUNIgk4SiU1FFkbR99n0IR8OKmveHRg6llJNCkyFseGTDrAJwnpAHTaVN\n+MXWX+DBjx8EkDomsHtgt2ylbpG+COMT49j6zFY8s+8ZANPfwZUtV2JX/y7FpujxRgAAbll3C57c\n9yQAKCSnSlOlIjCcLlI5KJknUGooxfPXPo+fffAzMSjM2hUxgVDimADrO5tJp+UJJJGDWBA+UUxA\nJgcZ0pODmHes5gmw9FBAWHw3HhjHFJ0Sq4imJQdlGBh2Bp0JtwRlsGcw/ln83hvfw6r/XIU/HP5D\n0vczWSuZEYh/TmfLQ+0PYVefepJEtlh0RkAtne/OTXeivLAcJp0J/ohfUUHx6X1P47rTrxMHF+Ya\ns5roMwmYfdD3gayNDciJznHps5finePvABAeMGuhFQ3FciPwUPtD+OGbP5S9b7dtNzREo7jxD48d\nhiPgSJoKO+AZEOv8ZAqbebY2tGL3wG5QSlMbAdtubKqTG4Hf7PkNjjqOikX72HdgyDeguawZh0YP\nyc4ViAQUOvmWhi0Y9g3j87HPVeUgu88OX9iHIn3RjD5jpakSg95BsaZPMs5fcj4evPhBXLDkArFN\nLSbgDXsFT0AtOygw3fety7cqSluoIcpBYQ/MeqUcVKRLEhOQykEF6clBo/5R6PP0CT0BY74RgJC1\nZcw3whlwihlDzGDFQymddWDYFXSlNAKsTEn8s2j32fE/L/if+Nar38KhkUNqbwUQ82jyjSjRqxuB\n17tfR/NDqRMOZkImk5e5ZtEZgWSzNg3RKDbimKJTeGrfU7hl3S1iG3OXmScwk4BZ/E3U6+rFitIV\nqueIRCN3clQzAAAgAElEQVT4aOAj7LXvhT8sGCeTzqTwBHbbdisGw9223bi48WJZtUZA8ATMOnPS\nYl3sAZ6NzOUJCYNOjbkGhnwDOsc74Y/4VW/YIn0RRidGcWDkADZUb5C1v3P8HVx/xvViHXwmBwFQ\nLb0QH/gFgDxNHrY1b8POzp2KFFGr0Yph/zBMOhM0ZGa3M8vNrzHXpFXT5s5Nd+Kuc6ZTORMGhmOe\ngNo+COz7e/aaZ7HMsiytPibKDnKH1FNEzTozgpNBdDu6RU+Aeb+pGPGPYG3VWvS5+xT9l8pBAMQM\nIeYJGHVGVU/AF/aBgIh9i5eDbnz5RtVtPKU4g04cdRxNOhNnxfTiB/Bh/zCuarkKG2o24IQ78XPD\njJmaJzDgGcCtf7oVfe4+WbLEbOFGIAPYPrGJiJeE3jr2Fkw6E86sOlNsYzOlYd9w2kEgSikGvYOy\nEgKuoAtTdAqNlkbVc+wf3o/gZBCHRg+JM2BCiGAEJGsFOuwdsvMCwO6B3bhxzY2y2Y8n5MHYxBjO\nrT83qSTEXPnZeAJs0AGEMsyvdb+GMkOZorYOIOj2H/R9gKbSJhh1xul2fTEuWHoBbjjjBtETGPGP\niMHldZVKI+AIOmSDPGNbyzbs7NqpeGh0eTpYCiwZP0hVpqqk91MyWHkGNliGJkOglIolFKQDYnAy\niOhUFIX5hTO6Bptdq8YEgioxAZ0JhBAxYYJ9L+mmiI74R9BU2gQCoigAyAZJBgsOT0QmUJhfqJiA\nMZgXAEAwAhJPIDoVxYuHXsRRx9Gk/WJ7JLDJz6/bf60wJq6gS1hMF7cOaMQ/Iq5wV9tDQvx84cQx\ngZtevgnf3fRdFBcUZxQvcAVdeGTPI4r2+BjXQrDojEAq/dasm04TdQQc+ObOb+LftvybbPBif+Rh\n/zBaylrSejicQSfC0bBssGZloRM9YLttu7GyfCUOjRwSFooZhcFP6gn4w370unrR7+4XZzmOgAN2\nnx1XtlyJY85jYvuR0SNYWb4SjZbGpEZA9ARmkfUkNQIbazfi1e5XVaUgAGJNG1ZTiPEPZ/0Dnrr6\nKdQV1U17AhPT34OaJ+AMOMVNXaRctOwidNg7cNx5XOEpVBgrZmUEkt1PycjPy4dBaxDvN+l3Fn9P\nsHiAmhFNhnTFcJG+CNZCK5xBJyLRiOgJGHVGmREABAPFVtUD6aeIsgGT7dAmJd4TsBZa0efug16r\nR54mL6EcJDMCerknYPPaEI6GZdtVquEMOnFG5RnoGu/CRGQC//SXf8KBkQOyY9whtyIo7w65oc/T\nw5BvEBcHJiKRJ0Apxft97+Puc+5W9f7S4fDoYfziw18o2rknkAGpMjmYJ0ApxS1/vAXXnHYNrlp5\nlewYqRyUrhEY9A5CQzSySD4zAiX6EtUHrN3WjlvX3YojY0cw7B8WZZBV1lU4OHIQwckgDo4cxGnl\np8lKE7Tb2rGhZgNMOhOqTFXi7OfQ6CGsrlitkJPUviMN0cwuJhBboQoInsB7J95LagTYcVLqi+vR\nUNKAuqI6WUyAfQ9rq9Zi//B+WQzHEXAoBnlAKN+wZdkWDPmGFJ5CpalyQTwBQC4JSb+z+MFCGg+Y\nCVI5yKwzI0+Th4aSBnw+9rloGOI9AUC4x1k8AEDaGS/MCNQX14v3Y/xeBYwKYwV6Xb1inEAqB0kz\ndD7o+wBnVAilw+I9AVbmIdlaHVY5dlPtJnSOd+LjgY8Fw+GTGw530I2GkgbZ5xzxj6DSJMR7MvUE\nvGEvCrQF0Gv1GRsBf9ivmqHEjcAMoZSmJQd5w17s6t+FrvEu/PxLynRKtqx+2D+MleUr05ohDXmH\nsLZyLXocPaL7n44ncEnjJSjSF2HP4B5RBikrLMMZFWegrbcNHfYOrKtaJytZLM2ykdZoOTQi1JqJ\nl5Pi6fP0oam0afYxAZ0ZALChZgOm6FRqI1CnHugsKShBJBqBN+QVS2cAwkBpMVhkhtUZVPcEAGBb\n8zbxfFJm4wlcsOQCnFuf3opvNcoKp9NEZZ5A3Mw704c9voooAFy49EK8e+JdMUXUoDWI5bRFI2CQ\nG4G05aCJmCdQJHgCA54B1P/velkWEMNaaMVx53FRImJykC/sQ9UDVTjmPAZKKZ7c9yRuXnczAIgG\ni8Hu+fgBXQqrHLvauhpd411o620DoCyr7Qq60FDcIPveh33Tk6+URiCBJyDdcrTUUIrxwHiiUyQ9\ntyvoku04F4lGEJwMygzrQpBzRuD6P1yfcPvHsYkxYXm6RJeMh3kCR0aP4Jz6cxTphkAsSMY8gfKW\ntGZIg95BrK5YjTxNnmjRRSOgknnhCrow4BnA6orVWG1djXdPvCtbaMUCnaIRkGxj+LHtY1FaaSlr\nEYPDh8cOY5V1VVqewLqqdXMmB5l0JpxecbrqQjFASBesNFbitPLTVF8nhKC2qBY2r01MEWXES0KO\ngHpMAAAub7pcLKImpdKYuSfw7Y3fxraWbRm9F5APCiyYztplnkBcQDtdpFIPO3drQyvaetvExWKE\nEBjzjWKAHABqzbVYXTG9DmGmclB9cT36Pf14tetVuIIujPpHlXKQUaiDJXoCMTmoe7wbrqALP3nv\nJ/h44GMQEHFSE79YrMfRg/LC8qRyENPN2YSorbcN59SdIxvQp+jUdGaW5L4f9g+LmV+ZegJSI6C2\nNiQdWKxEakCYnDdTiXCuyTkj8PLnL8tqxUhJFQ8AIJaOYNtPqmExWOAIOjDqH0VzWXNag+WQbwg1\nphpZHRlRDipQykGf2D7B+ur10Gq0WGVdhQ/7P5TtNHVF8xXY2bVTWMnMjIBTKE3wYf+H+MKSLwBI\n4gmkiAmsrVyruFn3DO7Bi4deTPlZAbkRAASpJ5EnUFJQgr67+pJm2NQV1eGY8xiCk0EU64vF9vjg\ncKKYACBUjO2/q1/x0FSZqhbMpZbJQSFvwpjAbOSgsYkxsSAfIBiBd0+8C1fQJX6XJp0Jdp9dHKR/\nc8Vv8Len/614ngJtASiouPlPIuJjAju7hI0DbV4b/BG/uhykk8tBneOd+OKyL+JPnX/Cve/ei1vW\n3SL+zeKzg3qcPTi3/tykchDLGGsua8aBkQP4bOgzXHPaNbIBnaWvlheWy57FEf9I+kYgwWKx8Ynp\nDX8yloNisRKpJJQLUhCQelOZHYSQYULIAUnbdkLIACFkb+zfpXHvWUII8RFC/knStoEQcoAQ0k0I\n+WWyawYng6qVJ4HUmUEAUKQTSkd0jXclLCJmKbCIbmylsTKtGdKgdxDV5mqxjgyQXA7abduNjTXC\nbH61dTVC0ZBsBrzKugpajRafDH6CtZVrsdyyHD3OHuwd2oslxUvEY1vKBU9g1D+KsYkxNJQ0oLyw\nHMHJoOrKS3fQDQqKZZZlij79/vDv8c9/+WfFOgo1vGGvOPMEgO2t2/Gdjd9JeLyaxyWl1lyLDnsH\nrEarbBCXegLBySCiNKqY6ae6zu0bbsf3z/1+0uvPF9KZodRwxlcYnY0cNOIfkRnkuqI6lBSU4Kjj\nKIoLpo2AI+AQB2mtRitLmSWEpOUNsOyt+qJ6dI534r0T7+Hc+nMx4BmQrRMApgPDrK0wvxD+iB9d\n413YVLsJ3z7723jz2Ju4cc2N4nsUnoCzB+fWnZvcE4h5UUuKl8AddGNN5RqsKF0Bu396QHcFXSgu\nKBYG8ND0AD4jOSjBYjGFHDShLgd12DsSGjPmCSw6IwDgCQCXxrVRAA9SSs+M/YvfJexBAK/GtT0M\n4DZKaROApnjDEU+iVYcnXCdSegJMDuoa70roCZQUlKBzvBMVxgoUFwj7s6YaGAe9g6gxq3sCanLQ\n/uH9OLNaSEtlbrlUTiGEYFvzNjSUNKC4oFg0Lm29bWhd2ioe11zWjM+GPsP6R9bjHzf9I/I0eWKa\nqVrOM1tRrTZj6bB3wO6z46/H/5r0swJKT6DGXCMG2DKh1lyLvfa9CklpZflK0dNhXsBM3WOr0Spb\nyZtNpN+zJ+RBkU7iCUgG3PhFbuli1BlBQRUL4VqXtsra2eCfTF9OFReQplMuKV6CTwc/xVk1Z+F0\n6+mweWyqctDk1OS0JxCLX3SOd6K5rBl3n3M3dly5Q/a3kXoClFL0OHpwTv05ST0BJgflafKwonQF\nWhtaFQM6S5eNN3TSwHCRvgiRqYg4IMc/H8wTYGmgLPYXbwTUPAFKKW546Qa8dOQl1c+waD0BSun7\nANTuGtWnlBByNYBjAA5L2qoBmCml7bGmpwFcney68aUEGL2uXiwrSb7Axqw3Yzwwjl5Xr7hQJh6L\nwQJHwIFKYyU0RAOTziQuNEnEkG9IZgT2Du2FUWeEpcCi+nBJ5ShWHiB+4/Eb196Ivzvj7wBMlyv+\na+9f0drQKh5TX1SP1oZWPHnVk7IgdyJJiK2oLjWUKmZ9HfYO/OOmfxRLMCSCUioLDM8FdUV12Du0\nV/EdsP12gcSZQbmMdGYojQlYCiyyweLI2JGMUlHZugKpVwYIklAeyRNn4ekYATXZUoov7IOGaGDU\nGcUaSduat6G2qBYDngFVOQiAIjuIbf9aXFAsBoQZ0hXDjoADFBSrrKuSBoalCwi/tupruOa0axRG\ngA2o8c+iNCtPuqXoMecxLPvlMpk8xjwBXZ4Oeq1eHLjHA+PymEBQaQT2DO3BodFDCT2NxewJJOJO\nQsg+QsjjhJASACCEmAD8AMD2uGNrAUgjvbZYmyp5JC+hHNTrFmbeySjSF2H/8H7UmGtUt7kDIN5Q\nbIaQTubEoHcQ1aZpOejH7/0YPzj3B+LCHOnDNUWn0D3eLcpRJQUlaClrUaS2nlVzFn685cdiH/I1\n+fhr719xwdLp0gR5mjy8/PWX8cXlX5S9N770BIPFTeIHIbvPjshUBD8874fY2bkzadmJUDSEPJKX\nsszxTKgtqkWPs0cmiQHC3ys6FYU35M04eLqQSCuJxqeIsnvKE/Lg1a5Xce2qa2d8fq1GC12eTuEJ\nbFm2BQ0lDaLXZNKZQEBUdydjxHus7xx/Bz948wfi78wLAAC9Vo/msmZc2XKluM7DF/bJkjJYjEiW\nHRTxC9u/JpBipbWDWIG7UkMpvGGvam0uQO5F3dN6D86qOUssF8Jm6yzIGi/lSAPDwLQk9NnQZ/CE\nPGKmETDtCQDylFpHwCGLCajJQU92PIn6ovrERiB2bunKaFfQhRL94jQCDwNYBmAdgCEAD8TatwP4\n35TSCSTwFNLh9IrTk3oC6RiBPUN7km4qwgaaikLhhk+llVJKMeQdEmIClkbsGdqDdls7vrnhm+L7\npTeezWNDSUGJbPZ25NtHxL0PEtFY2ogVpSsUA6UaDSUNqlvp9buFuEm8YWNZSFajFVuWbUnotgLy\nGe1cwWaW8XIQIUT0BqRu92IhYUzAUIoT7hOglOL3h36Pi5ZdlDCwngpjvlFhBGrMNei+s1v83aQz\nwagzJi2dEX9P7LPvw177XvF3qREAhHu2sbQRteaYJxCWewK6PB2K9cWy7KB+d7+YT6+GVA7qcQgF\n7jREg/LC8oSlI9i+1VIM+QYYtAbxuUsmB0k/EzMCHfYOlBnK8ErXK+Jr0hXR8UYgmRwUmgzhhYMv\n4O5z7k7qCSwtWZqTnsCMd52mlIriHSHkMQA7Y79uBHAtIeTfAZQAmCKEBAC8BKBOcoo6CN6AKqSN\n4KVPX4KtzobW1la0traKr6VrBFxBV8J4AAAYtAbka/JFTyDVQhpn0IkCbQEK8wtRUFQASil+eN4P\nxVlXcUExvGEvolNR5GnyBCmoXH79dHTuRktjwjTMeM6tPxd3vHaH+Pt7J95DjbkGfZ4+fGnZl1Cs\nL4Y3NN2nDnsH1lUKG62sr1qfdJl+fDxgLqg1C85fvBwEQNxuMdMMmoUkPibAJLQLll4AX9iHxz57\nDE/vf3pWgevC/ELVv4f0njLpTCnzzeMXNdq8tqQDJjMozBOwFloV16gwVsjkoMhUJOkEjHmXocmQ\nrNR1pbESw/5h1diOK+hCU1mTop0N6KwgZElBieJZHPYNy2JZbNVwh70Dd22+C4989ggeuuwhEEJk\n6yASGQHpuhDGzq6dWFu1FptqN+GFgy+ofm5/xI+lxUsxFpgbI9DW1oa2traM3hvPjI0AIaSaUjoU\n+/UrAA4AAKX0Askx9wDwUkr/K/a7hxCyCUA7gBsB/CrR+S++7WKUGkrxw/PlVTVdQRcmpyZTzhTZ\nQ5jsRiSEwGKwiG5iKjloyDskzuI1RINnvvKMLLdcQzQw68xwh9woNZQKmUmlM682eOfGO9OeCW+u\n2wybx4YTrhOoK6rDdX+4DsHJIPI1+fjGum8gT5MHs366Tx32DlzedLn4eY+MHkl47vkwApUmIf6i\nZuSqzcLG64vRE5DKA1I5yJBvwO+/9nucv+N85GnycNmKyzK+hlFnTBmfSccIxG8sM+AZkA1oI/4R\n0TuWwmICBq1Blh0ECMFhqRwEJH/2WF+9YS+OOY9hc91mAIIxSRQcThRUZ0bgNOtp4sI56bNo0BoQ\nioZkKclST+Chyx7C43sfx4GRA1hTuSahJzAeGBf3rFbzBHZ27cRXT/tq0uyjicgElhYvxTHX9MLI\nVJPVZMRPkO+9996MzgOkThF9HsCHAFoIIf2EkG8AuJ8Qsp8Qsg/AhQDUd8iWcweAxwB0AziqklEk\nYsg3qGYHnXCdkGmgiWAPYaob0VJgEWc9qeQglh7K+NrqryniDdJskM4xpSeQDuctOQ+nWdUXXMWT\np8nDl5u+jFe7X8Wbx95ErbkWH9/2MbYs24IzKoUl+tLgMJODgFjQUiW4xZDmu88VWo0W1aZqVamr\nxlQjlnRebJ5AohXDgJD59Oi2R3HPhfcgPy8/42uoyUHxpGUECuSTHZvXJvtdWtJDSrG+WJBEfUOK\na1gLraJhYPGLVAMbqx+0175X3FiH7QuhRqJYUXwxPWmxPFfQJXo20jGjylSFAyMH4Av70FDSIC7a\nBCDzBMoLy0WjJJ2csN3yWCVRSinaettw0bKLUGmSxymk+CN+NJQ0LD45iFJ6nUrzjlQnpZTeG/f7\nHgBnpNMhg9aguiw7HSkImDYCqQbhprImrChdASC1HMTSQ5MhjQt0ObqwtXFryr7Olm3N27CjYweK\n9cW4Zd0taCprwvPXPi/rkyPgQFW4Cn3uPqwsXwkAqplDUuY6M4ixpnKNrKgZg8UEgpPBjGdGCwUb\nWFlGVfxg/ZXTvjLraySSg6Sk6wkcHD0o/j7gGYAr6BKlkxH/CJYWL1W8j6347hrvUlyjpaxFJuEY\n840pJ2BmnVncm5uVGmHbhKqRqNKm1Ai4gi6xNDeb1E3RKcU+EVWmKrx97G2sr14PQggub74cP3n/\nJ/gfF/wPmSewvGS5WM5EagTyNHmi5FxWWIZeVy/C0TCay5pBCIEuTwd3yK3oL4sJKALDuW4EFgJD\nvgEBjzIw3OvqRUNxQ8r3lxSUwJhvlG0DqMbO63aKP8fPkOJhq4WTIZWUki1Um0suWXEJbvvzbdAQ\nDf7r8v9SvM4yVA6MHMDK8pXibDSV/DUfchAAvPZ3r6m2V5ursde+V5TpFhOskqg37IU35J3zgDow\nd3KQNIuNJTsYtAZRMhzxj+DsmrNV31trFoxAfMmW+7feL/u9SF8kTjYSYdab8cLBF3DZisvEVdBs\nm1A11PaYAKDcWyEm+7C9E4KTQYVnU2WqgjfsFb3iNZVrRGlU6gk0ljbi1e5XQSlVyJTM+ysrLBPW\n9TS0it4G65PCCLCYQA56AjlXNsKgNaimiKbrCViNVhy84+CMNhhJtfVev7tfJgepwR6w0GQINo8t\n5XqGuaBIX4RNdZtwcePFqlo6Ww/RbmuXPdyplr7PlxFIBAsML8aYACAMCsecx+bte5tLOYh5q2MT\nY2KVWnbvxweGpbBJVXxMIJ53b3k35baZZp0ZL3/+slgQEJgODKuRVA7yTxuBeDkoPijM3gNANAKV\nxkoEJ4PiDmksrsFqebEKotKV6tJ6Ue+eeFe2uDNRXMAf9qPCWAEKKsrduWIEcs4TKMwvVDcC7l6c\nv+T8tM6RjrGQwmYObBs86YB/ePQwXjz8It675b2k52APWI+zB0tLls5KA54Jv9j6i4SzxNICQfbZ\nbduNLy6bXmeQKgYiDXBmg2pTtVAiusCy6GICAPC9Td/D1me2Cts9zsP3dvc5d8sqgqpxSeMlKaU0\nqQc44BlAbVEt8jX5cAQcaERjUiNQa66FPk+f8r5eWqKUk+Ix6UwITgZx6YrpwgFJA8MJssbi5SBW\nQoPd32MTY4pANzMKzAgQQsT6XGyxGACxqq/axEQ6iWrrbcOPzv+Rap+kMKmpvLAcYxNjWFK8JGeM\nQO55AvkG1XUC6XoCmcAG8Cc6nsDmxzeLgR1/2I+v/f5ruP9L96cM2DJJqXMsceG6+YCVoVbtU8wT\niN8Ans2UEpXKmK+YQCKknsBik4MA4Lubv4u3bnwLN6+9eV6+t/OXnJ/SE60tqsV5S85LeoxUDrJ5\nbagrqpMZhlSeQLLqvTPBrDfjgqUXiIM2kDgwHJwMgoKqLvxU22UNEJ7FQe8gnjv4HM6qOUv2ngJt\nAe7ceKfMW2kpb8HBkYOI0ij0eUIKq7XQinA0jGPOY+JCMQYzAr2uXoSiIZn0m8wTYAXumCTEjUAC\nZisHZYLFYMGIfwT3vXcfRv2jOO46DgB44KMHsKZyDW5dd2vKc7AH7M9df8aayjXz0s+ZYimwoNvR\njRH/iEyn1Wq0KMwvVC1AB2RfDirSF2GKTmHAM7Ao5SBA2CDn0SsfTWuv4oWivLAc44FxRKIRwRMw\n14oD2hSdwnhgPOGCttqi2jmre7/CsgLXnS7POUkUGGZegFpWYHxMgA2oJQUl+MWHv8DayrX46qqv\nKt73q8t+JZN3mkubsde+F8Z8o3gdQggaLY34xPaJ4p4sM5RhfGIcbxx9A1satiiyj+KNAKVU3ILT\nWmjFqH80Z/YSAHJQDlJLEU13jUCmlBSUYK99L760/EvYXLcZbb1tWG5Zjpc/fxm/uvRXaS30shgs\neHzv48jT5GHPN/fMSz9nSqmhFH/p+QvOqjlLMTixh186G2NkWw4ihKDaVI0eZ8+ilIMWCyadCStK\nV+Czoc9g89hQa66FhmjgDDjhCDhQpC9KKPewkulzwb9e+K+KNqvRivHAOMYmxmSGKFkpEWuhVdT+\npXJQpakSdUV1+M0Vv0nr2W0ua8Zfjv1F4ekstyxH+2B7QjnozWNv4l++8C+y16pMVbJ9wQEgMBkQ\nt+BknkCu7CUA5KonECcHHXceT2uNQKaIdUkuvEfcsKPf3Y8+dx/OqT8n7XMMegfx4ldfzOoAmgyL\nwYIh35Bi20f2WqIMoWx7AoAgCZl0pqzFUk5VWpcK9zeTg9iAlkwKAoC1lWvx3zf897z1S5enw/c2\nfQ83vXyTTKZMJpnkafLw/XO/j+tfuh7haFgMWt+09ia0/z/tac+yW8pbsM++TxH0TuQJlBpK8dHA\nRzjmPIZLGi+RvabmCcSvPxibGMsZKQjIRSOQr5SDnux4Ehcvv3jerllhrMDLX38Z5y85XzQCr3S9\nIkthS8VlTZfhrZveEstH5wLs5lXb9jG+wJyU+agdlIoacw33ArJAa0Mr2k60iYFhFstKZQRYBc75\n5L6L7oM75MYvdk1vyJ6qlMj21u2YnJpEkb5InCSqFdxLRlNpk2yNAKOxtBH9nn7VmMCbx97EDWtu\nUExaVI2A5NzcCKRBvCcw6B3EM/ufwffPm79NQwghuHqlUN26qbQJk1OT+M9P/lOWwpaKkoKStLOX\nsgV7eNQ8gWQLxhbCE6g2VS/aeMBi4oKlF2BX3y70unpn5Alkg/y8fDy67VH8+pNfi22pKstqNVo8\nd81z+Kdz/inhMakw682oNlWregIAVGMCAHDzWnmZbCC1J2AttHIjkIr4FNH7P7gft667dd5nIQxC\nCFobWtE53olLVlyS+g05TF1RHbY0bFHNLEnmCcxH2YhU1JhrFmVm0GKjrLAMyyzL0O3oRq25VpQF\nE9UNyjbLLcsx4h8RM/TSKSpYW1Sr0OZnSkt5i6onACiNwHLLclzedLlYnkWKtVCIbbCyEoDcE1hm\nWYb3+95XXVC2UOScEZCmiHpCHjy176l59QLU+OKyL+LCpRfmzB8pU8oKy/DOze+oviatdx9PtlNE\nASG/PH6JP2d+aF3aisL8QpQUlOSUJwBAXJjF9ruQbgoznzSXNis8gSXFS6DVaMXicYyW8ha8cv0r\nUCM/Lx+WAgtGJ6bLQ0g9gUsaL8Eq6yr86O0f5cReAkAuGgFJiqjdZ4fVaM2aF8C4ed3N+MPf/CGr\n18w2yVZJL4QcdM1p1+CRbY9k9ZqnKq0Nrag118r2Hc4VIwAIMTo2iEo3ip9PmsuaFZ6AVqPF0uKl\nM5Ypq0xVGPIOib9LPQFCCB7d9ij0efqcmWTmXIpogbYAockQpuiUbAFINtFqtKqpkycTpYZS1U1p\npuiUYhvBbKDVaHMmq+pk58tNX5aVRmZyUDqbGWUDlku/onRF1jyBG9feiIsblcknT179ZMJ6SolY\nX70eu/p3iUkiUk8AEPYfeeumt1SrjS4EOecJEEKg1+oRnAzmVPDkZCNR0Tx/WKifksuLnjizQ6/V\ni1uYslXlueQJWI1WceGYdKP4+aTCWKGq8Z+/5PwZb7O6rXlbwh3LGA0lDWLV04Um54wAMJ0hxBZU\ncOaeREXkhv3DPEvnFMKYbxRXEOeKEagonJaDhn3Diy5WtLVxKz7s/1BckR/vCeQauWkEYmsFFkoO\nOhVItFjsE9sn2FC9YQF6xFkIWPnuPndfzhgBq9Eq1t3Plhw0lxTpi7C5bjPePPYmAPkG9rlIqp3F\ndhBChgkhByRt2wkhA4SQvbF/l8batxJCPo3tOvYpIWSL5D0bCCEHCCHdhJBfpupUYX4hApEAl4Pm\nkUSewG7bbtV1BZyTF0uBBXmavJx51qyFghzkD/tBKc2J+jozZVvzNuzskuxYNkfF9+aDVJ7AEwAu\njVxIkFkAACAASURBVGujAB6klJ4Z+8e2ihwFcAWldA2AmwE8I3nPwwBuo5Q2AWhihiMRLENIulEE\nZ+ZMTQEffaT+WqJy0u22dmys3TjPPePkEqWGUlgLrTPag2M+YdlBzAvIhfo6M2Vbyza82vUqolPR\nxe0JUErfB6CWR6j4q1BKOyilbKncYQAGQkg+IaQagJlS2h577WkAVye7Lisi5w7ymMBs6OwEtmwB\ngkHla2a9GRORCUSiEbEtHA1j3/A+RfldzsmNxWDJKcnFarRidGI0a0Hh+aChpAFmvRndju5F7wkk\n4k5CyD5CyOOEEDUf8loAeyilEQC1AAYkr9libQlhgWFXiMtBs2FgAAiFgN27la9piEaxt/L+4f1Y\nblme9bpBnIWl1FCaW0YgJgcN+xZfPEDKautqHB49vLg9gQQ8DGAZgHUAhgA8IH2RELIawM8B3J5p\np3hgeG6w2YT/29rUX2fpgYx2WzuPB5yCWApyyxOoMFZg1C/IQYstM0jKautqHBo5JO4lkKvMeLEY\npVTc+YEQ8hiAnZLf6wC8BOBGSunxWLMNgHTX97pYmyrbt29H78FePPnhk+g196J4EzcCmTIwAKxd\nKxiBe+5Rvh5fOmK3bTfOr8+tInic+afMUJZ2tdxswOSgxZgeKmV1xWq80vWK6jqB2dLW1oa2RLO7\nGTLjvzwhpJpSytZEfwXAgVh7CYBXAfyQUiqGIymlQ4QQDyFkE4B2ADcC+FWi82/fvh3dL3XjshWX\n4YGPHuBy0Cyw2YCvfx34yU+EuEBB3A598UXkdg/sxt2b785yLzkLzbfO/hYmpyYXuhsirH7QUedR\nrK9av9DdyZhV1lW4f9f9KNYXz7kc1NraitbWVvH3e++9N+NzpUoRfR7AhwBaCCH9hJBvALg/lga6\nD8CFAO6KHf4dAI0A7pGkj7Itgu4A8BiAbgBHJRlFqhRqhRRRLgfNjoEB4LTTgFWr1OMCNeYa2DyC\nUxaIBNDr6sXqitVZ7mVyXC4gR1bXn7RUGCtQY65Z6G7IqDBW4ODIwUUbGAaAleUrcdRxFO6QW+EJ\nBIPAxESCN2aZpJ4ApfQ6leYdCY69D8B9CV7bA0C5JjsBLCbA1wnMDpsNqKsDWlsFSejCC+WvN5c1\ni1vhHXUcxTLLspySBQDgq18Ftm8Hzucq1SmFtdCK/cP7cypWMVMK8wtRa65F51inwhP4j/8QJjg/\n//kCdU5CbiQGx2HQGuAP+xekmuXJhM0G1NYCX/gCsGuX8vXmsmZ0OQQj0DXehZayliz3MDUeDzA0\nlPo4zsmF1WhFYDKwqGMCgCAJRaYiCk9geBjo7V2YPsWTm0Yg34CxiTHotXq+52yGhEKA0wlUVAAr\nVgAnTiiPaSlrQedYJwCgc7wTzWXNWe5lakIhYGxsoXvByTZsg5vFLAcBQoYQAIUn4HROZ+8tNLlp\nBLQG2P25s/POYmRoCKiqAvLyBG/AZlNq642ljeh19WJyajJnPYFQCBgfX+hecLKN1WhFHslb9MUM\nWYwt3hNwOLgRSIoh3wC7z86DwrNgYEAY/AGgqAjQaAC3W35MgbYA1eZq9Lp6c9oTkBoBuz3xsZyT\nB2uhFVZj7pSyyJRV1lXQarTQ5elk7cwTyIWkh5z8hg1aA4a8Q7xkxCxgQWEG8wbiYZJQ13gXWspz\n2xPweIDGRqEmEufkpsJYsaiDwoxV1lX4xrpvKNodDiAczg2pMyeNQGF+YU5txLwYYUFhRl2d4B3E\n01zWjA/7P0R0KgprYW7sLCVFagSGhoS0Ovbg/PWvwLe/vXB948wfjaWNOK38tIXuxqwp0Bbgt9t+\nq2h3OoGyMvVnMtvkpBEw5BvgDDq5HDQLBgbS9wRe6X4FzWXNOVmtURoYZlLQ4KDw/8GDQI9yh0zO\nScC59efiha++sNDdmDccDuCMM3IjLpCbRkBrAADuCWTA668Lskm8J1Bbm9gT2D+8PyelIEDuCcQb\ngYGB3Flww+GkSyAg/L9iRW4YgdxaGRTDkC8YAe4JzJw77gDOPltdDuroUB7PgsHNpbkXFJ6aAiYn\nlUaArRuw2bgR4Cw+HA6gtDSxRJttctoT4IHhmePxAHv3CpvJxMtBajdcfXE9CrQFOekJhEKATgd4\nvYIxsNuB/PxpT4AbAc5ixOkELJbEEm22yU0jkM/loEygVEgD/dOfgDVrgBpJOZi6OvUbTkM0uHDp\nhVhXtW5G17roIuDIkVl2OAWhEGAwACUlwoNjtws6KpeDOIsZZgQSPZPZJiflIFZ7m8tBM2NiQpgp\nr1qllH4SeQIA8MYNSev5KYhGBU+jp0coUDdfhEKAXi8YgbExwQiceaYgB1EqPECmxbf9LOcUh8lB\nyZ7JbJKbngCXgzLC7QaKE3xlVqsgq6htNTlTjh8XzjMykvrY2cCMQFmZEBew24H16wVPwOUSJCLu\nCXAWG7nmCeSmEeByUEYkMwIaDVBdLQygdvvsSjEcOiT8Pzyc+TnSIZkRGBgAli8XMi1yYdUlh5Mu\nzBMoKREWjPl8C9uf3DQCWp4dlAkeT2IjAExnI1x/PfDII5lf59AhQauXGoH5GIiZESgvF7yOsTEh\n1jEyAvT3A0uXAlqt8CAtNrjhOnVhngAhuREczkkjUKAVtsDinsDMcLuFOkGJqK0FnntOWGnr92d+\nncOHgfPOm5aDHA5hVj7XSD2Bzz8XHpzCQsHQdXQIn6ewcHFKQps3A52dC90LzkLAPAEgNyShVDuL\n7SCEDBNCDkjathNCBiS7h10mee1HhJBuQsjnhJCLJe0bCCEHYq/9MlWnCCEoM5TBYrBk+rnQ3Q3c\nfYrtlJhMDgKEG+6RR4TZNFuwkgmHDgFbtkx7AseOCbXRvd7Mz6kGSxEtKxOuWVUltNfUAJ98Inye\nxWgEKAUOHBC8Gc6pB/MEgNwIDqfyBJ4AcGlcGwXwIKX0zNi/1wGAELIKwNcBrIq957/IdB2ChwHc\nRiltAtBECIk/p4LuO7th0mWe+vHCC8Brr2X89kVJKiNQWyvM2P/+7zMfOKNRYQbb2jrtCbDBbK4H\nNaknIDUC1dWCEVisnsDQkGCEeYnsUxOpJ/DTnwLbti1sf5IaAUrp+wCcKi+pFZm5CsDzlNIIpbQX\nwFEAmwgh1QDMlNL22HFPA7g6Vcdm4wUAwM6dgsU9lUgVE7j5ZuF7MZszHziPHxc2qlm+fNoT6OuT\n/z9XSI3AwIDcE2Arog2GxWcEWL0jbgROTaSewJIl0z8vFJnGBO4khOwjhDxOCGHCfQ0AqWMzAKBW\npd0Wa583hoaE2arDcWoF4FLFBMrKhLz+2cyeDx0CVq8WgrUOh5CmOd+eQHm58LvUEwCm5aDZSFsL\nATMCuVBGmJN9nM5pTyAXyGSx2MMAfhz7+d8APADgtrnq0Pbt28WfW1tb0draOuNzvPYacOmlwqzX\n7z91FhS53cLMIhWzNQKrVglZORaLMJvt6wNWrpxfTwCQewLA4pWDenqE9EDuCZyaOByzn/23tbWh\nra1tTvozYyNAKRWXCBFCHgOwM/arDUC95NA6CB6ALfaztD1hPFxqBDJl507gq18VNld3Ok8tI5BM\nDmLMZuD8/HMhHgAAlZWCJNTfL2QLxXsCzz0nFIG74YbMrpXMCOh0goewEEbgt78VArs/+5kgrc2U\nnh6hyN98GIGnnhK2FM30O88GDz8sTFYuv3yhe5J9pqaEhY6zNQLxE+R7770343PNWA6KafyMrwBg\nmUN/BvC3hBAdIWQZgCYA7ZRSOwAPIWRTLFB8I4A/ZtxjAO+/Dzz7rLrUEwoJKZCXXSZ80Q7HbK6U\ne9hswO9+p/5aqpgAYyYD55EjwCuvTP8+Ojo9GFdUCEagr08wAvGewCefCOUlMiWZEaitFfKsF8II\nvPUW8OGHQpZVb+/M39/TA2zcOD9G4A9/AD79dO7PO5e0twN79ix0L7KL2y1k5nm9wj2rzaGCPalS\nRJ8H8CGAFkJIPyHkGwDuJ4TsJ4TsA3AhgLsAgFJ6GMCLAA4DeB3AHZSKw/QdAB4D0A3gKKV0ZsVq\n4nj2WeBb3xJmEvHB364uYYAoKxN0t5MtONzWJsyk1EgVE2DMREf/7/8WZpcMl0uQMgDBE7DZBG17\n0yalJ+BwzG5VMTMCOp3gzTEjcNZZQvYX+yzZNgK9vcLf4OyzgQ8+mPn758sIUArs3p37Ex+f79SL\nh+zaBdx+u/D85lI8AEghB1FKr1Np3pHk+J8C+KlK+x4AZ8y4dwno7RWkhl/9SpB+brpp+rXOTqA5\nVhr/ZPQEensTf6Z05aCZZNTY7fLrSTMbKiuBzz4T/m9oEIzA1JRQogIQBjmPJ73rqMGMAADcd5+w\nvzAgyB0bNwo/z5cRCAaBggL113p7hc9bUyN4RjPB7RY+16pVcz8QnjghpO1KJz4+H2A0Cl5TrnAq\nGoGeHuHZ/P73Fz4bKJ6cXDGcit5eoKlJyHRxueSvdXUBLbHS+CejJ9Dbm/gzzUdMwG6XX09qBCoq\nBOlhyRLhnGazfFCcK08AAL773emfpcxHiujRo4LXofY9+3xCskFFhVCUb6ZF9Hp6hPRaVg9pLmlv\nn87aYlx2GfDuu3N7ndni9c7ceC52enqExatuNzcCs2ZqStCelywRZAk1I3AqeAJq8ZD5NgKUCj9L\n5aCODqA+lg5QXy+XhByO2VUalRqBRMxHiuh99wmDvZrUc+KEULOIEMEIsMGMUmGPhVTlOHp6BI+m\nuFg4NhKZu37v3g1cconceNlsghSRS5yKnsCxY8JeGP/yL/J9PnKBRWcEhocF3buwUN0IxMtBJ6Mn\nEA6rD3weT/oxgUzkoEBAkGKYTFJRIbSxtNQlS+TB4fFx4e8TCqV3rXjSNQJz6QkcPQr8/+19e3hV\n1bXvb5AQJCB5ESAhQUJ4IxBERduqaAWt5eHRetRTX+g5vdVz6q21tWr1yNHWR1tbvafHHj8tWK3S\n6/VaK59XRFsjtrRisWAQERB5hCTkIZDwMEAy7h9jDdfca6+13y9g/r4vX5K1195r7LXmHL85nvOV\nV2TV9uab4a+rKwiQ768ksHu3nB8tEK6WQJ8+qV+kvPOOrPzNz2xrk+OA3M8bbwQuvxx48MHUXTde\ndHWlhwR27kyuMWI6oeR/883AokCHenaQcyTw9tuRV0fmJAyyBEx30NFkCaxfH9pH5M03pU2DordX\nVtp+3+vQIblvhYXRr1NQIJ8byyq0pUXI5ciRUFcQIJYA4G8J9PbK+Ym4TBTZIIEf/hD41reAiy+W\nIJ4X5vgzv5vudhYtdVuVASCum1Qpw8OHxSqbNUvuO7PENfbtExJglvG0cqWc8/Ofp+a6iUAtAbVm\n33ordJwnipUrgSeeSP5zUo3eXqm0HzVKLMiCgmxLFIqcI4F/+RfZIzcIkUigo0MGU3m5/H+0WQJ3\n3QU86rTXO3QImDNH8tEVzc1CABUV4d9LXUGxBAA1tTKaG6WnR+7poEFyn01XEOCSgJ8loKlwVVW5\nSQKHD4dbKMyyNec3viEZSJs2hd9nLwmoJdDcLPGJaCSwaRMwerT8ncq4QEODuKmGDBFr7cAB+exh\nw8Tq2L5dUn2vvBK44QZxRXkXUJnCvn0yvtV19o//CKxbl/znNjYmT6oHDqRm4yUTzc0yh3K1Xinn\nSODTT0WhBSESCagrSBXh0WYJvPOOq0T+9jcZkM3N7uuffCLf3Y/cYo0HKGIhgbY2uYfl5XI9ryUw\nZIj89rMEOjpEyWktQSJIJwk89hhwySWhx5qbZXvOYcNktXbmmWKZmghyBzU1SZX6mjWR4wKmpZpK\nEli3DqhztonW8dHeLs9uxgwZW0uXSrMyIpknGzem5trxoKdHlGxlpcinO9SlIlC8c2fy9/OBB4Cf\n/jR5WUyY1l8u4pgiAXOCAbFZAuefn1waY6rQ2ChKb8MG+f5KBupmANzv7kdusRaKKWLJqmlpEYWo\nvmtvpWP//iKPPg/TEtBOiUOH5qYl8NFH0l5E/eWA2xdJMXNm+MreHH9FRUKk3d1CIKNHA1OnBscF\nurrkHlY59fOpJIGODtcC1vFhksCTT4qFMHGinDNuXOZI4Nln3bbu+/dLyuqQISKfuj9TQQKNjXKP\nk9lkqKUl9T2wLAnEiZ6eyGZqNBLQoDAQW4ron/6U/m0SY8E778hGI6efLjLV1wMnn+xPAkGWQCxB\nYUUsylNJQO+j1xIAJOtBj5mWgEkC6bQEEk0R3bpVig3Navv1610lCUQnASLx67e1yXOqrPR/j2Lj\nRklt1jqKwYNTRwLms9Hx0dYm15gxA3j9ddcKADJnCaxbJy7eVU4P4X37xC2i8RAdL0ELhXffBZ56\nKvz4j34UvhDSzVmSsf4//VTGfSphSSABJOMO8loCkQbEoUOiaCJdL1m8/XZsvtd33pHJOnMmsHy5\nrCYvvzzUHRTJEkjEHRSPJaAkYMYEgNAYREWFKJ5Dh1wSGDIk/ZaA1621YkV0627rVuDee0VJqTXg\ntQSmThXLrLdX/jdrBBTqEmpqku8/c6Z/VhEQvkgpK4vPh71iRbCryWxKpqTd3i7K9tRT5TmZfevH\njQve2ewPf3DHBnNo25B48Nln4u+/6SZ30dLVJfUkSgJqOQZZAq++KjEMb63DL34R3h5j507JXEuG\nWNNBAlu2WBKIG0FKU2sETjpJ/h80SCa7TlLvJCsulkEXlHmgO2GliwS6uyXL5L77op9rksCiRbJi\nnDQpPksgXSSgpONnCZjIz5fzm5rcmEC6LQG/73HHHZELpJjlXo4dKwrq6afluJcEBgyQ76CKyqwR\nUGhwuLlZLIEzzpDEhiNHwq/rXaTE6w76p38CHn7Y/d8kELM9sS5+lAQGDQJ+9SvgnHPc84Msgbfe\nAmbPFiIA5NnNm+fOsXjw3nuilL/9bXe8mpZAW5tYAiUlwSTQ3CzX//rX3cVET4/8bZIYs4y7SZOS\nCw53dFhLICcQpJTNGgFAlE5hoQws5vCbnZcnAy7o83S1mK4siVdeEf/vU09FHlhHjsiEOe00cQcd\nOSJkUFGRvphAIpZALN0PNS5gWgKZJoHduyMr17Y2UfADB4qCWbpUxo+2yTZhKkvTClVomqi6g048\nUXpX+a2y/SyBWEmgrU3u6X/+pzyHX/86lFAiWQIAsGCBBL2938tU7q2tomwnTgzdKIhZ5li8+Phj\nkbGkxJ1jXV2h7qDt24Hp04OtxaYmaQtz4YXA4sVyrL1d5DZJrL1dnmlVVWosgVTuQ2JJIAEEKW2/\nSaguob17hRS8rX0jxQWUBNJlCTz1FHDLLcDVVwMPPSQTeMKEcMvkgw9EgZSUiNKbPVuKfiorXXdQ\nT4+smnQnokzHBNQS8LqDvNC4QDYDw3v2RFYE5jiaMEGU47JlkhGkwVWF6TbZvBmoqQl9Xd1Bzc3u\nZjfTpkmWkBfJkMDatbJA+OpXpeDou991+xABwZaA9/soBg2SH11k9PbKOL3mGuCKK8I3Ckpkjqjy\nKyyUhU13t5DJiSeKXCYJBFkCSq5Tpriy6ILKJFrdaS7ZYLsurpS0LrssPEMsHnR1ibsy6DnkAnKO\nBEpLg1fmW7e6riCFkoCa415Eigsk6g665Zbo7QF27RIf7te+Bnz/+8B//7d0vmxtDfXzAxL80oZo\ngOSqz5olCnjXLpmgTU1yb044ITUxgf79o6eI+sUEYrUEspkiGs0S8AZ3584VkjZdQQrTEnj3XVFY\nJsrLRRmdcIJrodbVhZMAc2g1OxBfYHjNGvncu+6S5okPPSQEq8rTtAS8geEgmN/twQflPt57b2iW\nl/5OhgSIZJ7u3u0fGI5EAjqvhw93A78tLVJ4ZVoCjY1iBUS6p/fcEzkL6eBBmWsnneRaA6+/Lvc8\nHsvg3ntDs56GDMmtBn5e5BwJ1NQEDzi/bdmUBDQw50WkNNFo7qCGBuD990OPNTUBjzzibhEYhOee\nk3jAwIEi17p1QgrjxoX3oN+wwV8BFRSIYm9vl/frOZmOCQSliPrBawmUl8uk7O2V1XY08jSRCAl8\n9pn8RFKuWm+hmDNHfOFeVxAQqig1bmOivFxW6ebY8yOBlhaXwBUaGG5tFfdOJEWjJDB6tIy9BQtC\ni9XMuWGmiEYiAbVyXnlFOvIuWSLW9IgR4ZZAIi5T0w2iLiEzMNzWFtkSUN//0KFCAqpYW1rkObS0\nuAsZ0xIIigk8+WRwMBxwx2xFhdszS2MNfhlfjY3hVsKhQ0LQWvymi6FcxlFFAgcOiN/PRDRLIFLB\nWCR3UGur+CEfeij0uAYcTV+9H1avloZiitpaiVGMHBlOApF8hpWVci1VAkDms4MipYh64Y0J9O0r\ncq1ZI/735ctjlzGeFFFVoKqsIgUHvW7Fs88W14gfEaui7OgQi2bChNDXhwyRxYI59urqhBhMpe51\nBQHuc5wyRYLZkfremM/fbGDX2uo29vNaAtFIYOxYKY765jeBF15w6xeqq1NrCZgymZbA5s3yPUaO\nFHLwtjFpb5f5XVAgspmWQFWVWAObN8uxnTvlWJA7qKdH3heJzHTMDhsm5378sZDu3XcDCxeGk/Tv\nfgf8+7+HHtMiTyU1SwIJoKYm+EFpoYmJZC2B/Pzw66l/VCstTdTXizL3unS8aG0NTSVU+JFApBQy\nDQ6bSsDvO8XirzcRjQQOHpSf4uLIKaJemO4gXZkOGSIutH79JBc/VsRCAnl5QjRaIKT3JVZ3ECBK\n5gc/AL785fBz1TWwYoWkWublhb5eXi7j0iQBHYfmGPFmBgEi9403Ai+9JGmld93lH0s4eFDGiNdS\n0XhEV5dYGRr4jdUSmDVLqqYbGoAvfck9Pny4yK5xqJNOip8EurrkR++F1x1UXi7XGDFC6ib8VvDm\nnB42TL7r4cPu4sS00hobI8cE2trcflZB8COB2lrJylq7Nly+rVvleZnkoBaDntvREfkZ5AKi7Sy2\niIh2EVGDz2u3ElEvEZU6/59AREucXcfWE9HtxrnTiaiBiDYR0aORrhnJEohGAkGWQNCD7+qSgeO9\n3tKl8vCWLBFlbj78+npp1xvNEmhr8w8GeUnAL6vJhAaHTRLQ72xmdrS0+JNgEKKRwK5dMhmIYk8R\nBULdQboCGjpULKOFCyUIHitiIQEg9Lvs3i3/x0MCAHDbbeFKGpBFQk2NuPfMuI1Cn7F574nCXULv\nvefvbvrFLyStdNw44Mc/lnRKLz74QBSet/GYuoO8btKSEiHigoLgjXEAqYN45JFwYtctPZub5XOm\nTInfHbRli9w3LYzThYS6g1RebTmi32XrVuCUU+SYad3n5wvptbS4JGAG7aMFhnW+RiOBsrJwEtDU\nZz8S2LMntHNufb1k+em57e1HvyWwGMCF3oNEVA1gFoBtxuErAICZpwCYDuB/EJHTWgy/BHADM48B\nMIaIwj5TMXKkrND9/KORSCCRwHBnpwxCLwls2SIro379ZPWn1Y5NTfJQZ8+ObgloQMjv+5kk0Noq\nEzXIlVNRIaudHTtcJZWfL6spsyBKJ0asiEYCLS1ugzjN4+7ujt4Eq7RUVuUdHS5hVFYC//qvkqee\nCRIYPTqYBJjdfP9YMW4c8PLL4fEAwH3G3rFnkkB3t7hbvva1yNeZNy98ZQmELgBMqDvIDAoD8gxa\nW5NbgY4YIa6WTz8Fxo8PXpgtXgy8+GL4ca91qzEBtQT69pW5q80H9busWiV1Frt3h1v3GhcIsgQi\nBYZ1vkYiM7VeTRIYNUpe8yOXrVvl+etzPnRIijwvuSTUEjiqSYCZ3wbgx50/A3Cb51gzgAFElAdg\nAIBDkA3mKwCcyMyOKsXTAC4OumZ5uShFv7zkRNxB0WIC1dXhA8P8rBkzXBJ46y3xH1dVRbYEmGO3\nBKLlEFdWyj6/EyeG5nmb5KZtnuNJQ4vWQG73bnfwFha6GR7RshyI5J7qRAekM+p994k/fdMm/0Iq\nPyRCAnv2yP3UVsXMbn45IIpGawRixdixMsH9SKCoSL6nd+zV1Uk2ESCW5ZQp4daHF2VlMva9YyuI\nBNQd5LUEdGWfLAn89a/yvUpKgkngt78V1+mHH4Ye945rryWg8qkloN9FFerateHWvcYFTBLwswT8\n4kGxWgJ+7iAgmASUuAGJB4wZIz/HDAn4gYjmA2hk5pC8GWZ+DUAnhAy2AvgJM+8BMByA0SUfO51j\nvigtlYnlN+gScQdNnuwqcS+CLAHTqjDjAvX1/kVcfp8bZIpr5oXWCpirDT9UVsqE8CoB083V1iYD\nzeuvjoRoPXf27HGtEyKZxLFuizdiROjAHzxY7seAATLBtmzxv55pVjMnbglUVMj19u0TpX/99W46\nsJ8rKBrGjRMF5De+tH+Q97ULLhASWLZM6kWuuy62a02cGG4tvfuu6yIxEWQJ5OXJs0uGBKqrZUey\n6mr/fTsUGzfK/guXXRY6nrwk4I0JACKfaQkoCYwYIb+91r3XEhg3TjLrXnxRSLqkxE0x91Y4NzXJ\n+1NFAp2dkoV2/vkuCah+0Mwn4BgkASIqBHAngHvMw85rVwHoD6ACQA2A7xJRTdiHRMHjjy/EoUML\ncf/9C1Hvycvavz980xQdXGaxjonTTpMH8skn4a91dcmA85KAnyXw4Ycy2ObMCS3iAsIHXJArCBBi\nUH8rEN0SUDm8JGAGh+N1BQHR3UHebKN4SKC6OjyVVzFxYnhw+JVXJDNnwQL32JEj4k+Ohdi8JFBc\n7E5atbo01TEREjj3XIkZBOHf/k2a/ZkoK5PumddeK8r00ktju9akSaEk0N0t6Ybe+gQgOCYAyLNK\n1hL485/ld9Ci7LPPZBz/8Idy3vPPu68FWQImCVxzjRuQVkJbs0aOr1kTbt1XVcnnHjggn1deLmnY\nixdLhhOR6yr1klZzs4y9aNlBGhPYtk1W82bHV9PC2LZNxpEWBvb0AL/5jegHc7OgaMH5RFFfX4+F\nCxd+/pMM8uM8vxbASABrSfwCVQBWE9EMAF8A8Dtm7gHQRkR/hsQG/uScp6iCWAO+eOCBhXj7tRaa\n6wAAFuhJREFUbSlf/+IXQ18LsgS2bXMfvhd9+kiV5dKlUmlporNTHvKePbLyVFeHaVVUVMg1L7wQ\nuP9+Gdjd3aJ4e3tFoV16qZCEBsFaWyO7ZtQlpIPaTCX1QuXwswTUHZQuEjADhqWl4fc+CN5tJk2o\nkrvYcQiuWyer5McfFxLo6RHF390d+w5MZuHbnj3yzPxIYOLExEhg1ChZ7Qbhzjv9j59zjhQKtrTE\nfu8mTZIgumLNGnEv+L1fXSheSwBwlWSiqK4WpV1dHUwCmzfLvczPl2aHL7/sWjx+JGDWCQCSGWV+\nl9deE2KZN086jxYUhFsCv/+9xKp0rpquPoXGBUxi1L5C770X/J31PeXlohvGjHEXId5Ygxau1tbK\nM3jiCZkv554bmkySLktg5syZmDlz5uf//4fZDjdOxGUJMHMDMw9l5hpmroG4eU5h5l0ANgA4DwCI\naACAMwBsYOYWSGxgBglzXA3gpUCB+sTvDtqwIfLmzXPnCgl40dnpuipMheg1Q7/4RVmx/PM/y//9\n+gnhdHSIq2jjxtCGZUHxAIUZF4jWYXDYMPnOU6aEHk+3JWC6g/R6sVoCo0YFy+N1d3z0kdzbSy4R\nRbBhgxyP1RUEhFsCJSXhJKCklAgJJIPvfEeyfmKF9/6sWuUfiwDc1bOfJVBamrwloL+D3EHm/h0X\nXQS88YbblbepKTT47ucO8n6XP/5RFjsnnyyfvXVruCWwZk30se4XF1BLIBZ3UH6+yGPOS687SMdR\nXp7MzVtvlew3Ijl3925ZJB717iAiWgJgJYCxRLSDiBZ4TjHzGB4HUOCkk64CsIiZddO4mwA8CWAT\ngM3MvCzSdYuL4yOBgwcjp0fOmiVBLm974a4uKRIySWf/fvEvmgrwySel06QZFDVTN2trQ3ueR3IH\nAaEkEM0dVFAgCszbFygTloB5D0pLYyeBK66Qnbv8MGlSqDtI+yEBofGXZEhA3UHt7XKfy8uTcwdl\nEnp/NEPIr0pZUVTkumS8zyZZEjB3izPnx/btUiQFhBbAlZeL8q6vlx5Zl10W+vxMd5C3v5e+v7NT\nSKB/f1lItLWFjuvhw2Wux0IC3iBuU1PsJADINWIhAUBknjpV9Azg9jDTHlZHNQkw85XMXMnM/Zi5\nmpkXe14fxcyfOn93M/NVzDyZmScx88PGeaud46OZ+WbvdbwoKvJfeQSRABDZEhg4UFbz3mrVzk6X\nBPR6agWYCn/AgHDftAaH164Vn+jvf+8GH2O1BPbtExmi5ff7+dczERMw3UHxWAJ9+/pPdEAyhDZu\ndAPj27e7CicVJKCtLUxL4KyzsmcJxAtvhlAkEtCq4Y0bw5/NXXcB//APicsxZIgsQLwxgeXLpZbB\nrxfSnDkSB3n0Ubm+CdMd5GcJ6KJJ3Z51da6VrhjupJPESwLafmL8+Ogpoqqw/UjAtC7McXT77VJH\nYuqMwYNl4XHkSOyuwGwh5yqGgfjcQbpCjqZIZ88O7//R2SnKyrQ8grKMvKislLS0tWtlBTBzpuSC\nA7HFBDZuFFPfLKiJB8OHSywESIwEojWQ87qD5s8Xkz9Z6NaCGqg3LYHTT0+NJaBBUSWBs892WyLH\nWyOQDahLqKPDVV5B0AZ23oXClCnJxQT69BH3xtixoe6gHTvkXjY0hG/nOncu8Mwz0gHX2yLDLzDs\n/R5AKAl452FhoXxOvCTQ1uZmDh08GN6eQmFaAgsWSIZX0GeaJDBiRPjCYvBgeS6DB+d28zggR0kg\nyAfpRwJqekVT3NOmicJWMLtBKpN0guoNvKisBFaulPeWlUkgW0kgmjto0iRRRlddJYGkRDB1qvt9\nMuEOOu88WVGnAub+tqYlUFcndQT796cmJtDeLvf5rLNEeSVSI5ANaPB81Sr/VhUmystjq+ROBHfc\nIYuFgQPF7XTkiDyvsjKJsXn7IU2aJIuFu+8O/6xBg9x9P/yea0mJzAUlvHPPDd0ER1FVFX2sm9k5\ngLuw01oXP91y8KBYDJp9eMUVoX2i/ALDkSzKwYMlvpXrriAgR0nAzxLo6RFfff/+4ecXF0cnAVWa\nms558KCYmtrgzOsOioaKCtn6TlcuM2ZIpSMQ3R2kq/imJmkbkAjGj5cJuX9/ZtxBqYRZ6bl9u2sJ\n9OsnfuXVq5NzB2lMYP16UfrjxwsJeLuH5iqmTpW9j6+5BvjCFyKfq4uNoJTcVIBIFkudnXIfr79e\nVvyHDrlV5XreSy9JVo0XffoIEQQRcJ8+EhhW98/06dLZ1ItRo9zxEoTy8tBNnMz0cXVL7d0rcQuz\n51RpafCqXWNwzHIfou0RUF5uSSAp+JGA1gj4PaTi4tj86iUlrhtC4wH6/kQsgeZmlwSqq2W1tGtX\ndHdQKtC3r6xUGhoykx2USmjPl+5umVim7BoXiJcE9u+XhcL+/fJcy8okHXDkSHl94EApujoaSOC6\n66T69C9/Cfete6HjLB2WgAldQW/fLhZsU5OQeTyujpKS5K2w556L7pY87TRJBFGYLl7NUvrwQ7Hc\ndWvRaAHcggJZgHZ2+m8z6oW6gywJJAg/k83PFaRYtCj6igkI7eei8QAg3B0UqyWgnwm4TcPWro3u\nDkoV6urEJdXdHb/CLigQpRnkH82EJbBzp9xH091x8skyQeMhgZoamXB79ggBaFfKvXtD/bYrVhwd\nJJCfL/2PRo+Ofg/Ky2XspYuwFWot79ghAdMLLvBvuBcJxcXBCQOxorAwegxt6lRZoOlmRubCTmMT\nH38sFsuPfiTWgBkPCIK6GP06wnphSSBJBFkCQSRw6qmhfXWCYJKApofq9eJ1B/kVcdXViUsomjso\nVairkwIb7fYZD4iC+wf19IiVkC7fufZ82bHDjQcoamuldiIeElDrwfSN6+RTpV9dLRuAHA0kEA+G\nDBHlmkhyQTwoKhLFWVgo8/DWW0MrvGNBKiyBWJCXJ3GgFSvkf3NOqzvo44/FHTRmjKS0rl8fGwl0\ndPjvRe3F4MGis3K9jTSQoyTgVycQiQRihdcSSNYdNGdOqFKpq5OBF62Fb6pQVydFavG6ghRBLqHO\nTpms6VIsI0bIZFq/Pty/W1srEzQeEjjpJIn1NDS41ouXBEaMkJXhsUYC5eXpdwUBcl8bGlzSPuMM\n//0XIqGkJHlLIFbMnCnZgMziGtIAtrqDtD7n/vulDujhh4NTcRUmCfhtQGRClb+1BBKEX51AOkjA\nzx0UqyVQUCAZEqaiVKWcCVcQIGmA3d2Jk0BQmmg6XUGA3LPRoyUQ6LUEqqvFktq7N3YSIJIJvGyZ\nqxAHDRK3imkJAMceCZx0kps/n04UFQkJRAvKRkJxceYys5QEli6V8aGEZbqDamvFi/Dhh9IC4447\nIn+mZgitX29JIO2I1x0UK3Qbu/b2UEtASWffPvGRJ+pfHT9e3p8JVxAgckZq0RANQZZAvFtVJoJx\n44QEvEolL0+OffRR7CQAhJOAdvc0LQEg92sE4sW0adKuId1QEvCSdjzIlDsIkLhAUxPwve/JBvPq\nLjXdQZG69/qhrMztLhotJqA6wJJAghg4UFaoZt/5VJAAkQyONWtCYwLqDlJXUKLFHQUFskLIFAkA\nYn2kmgTSmRmkGDtWgnF+SmXUKFltxUsC27eHWjDLlrm+2+pqIYVcrxFIBLE22ksGxcUSq0nGEsik\nO0jjAoWFUrugKC6WhIQ9e+K3oMrKxLVUVeWfqm7iaLIE4u0imhFoTvHeve5NTAUJAMCZZ4rLZsCA\ncEvgvvtCqwQTQV1dZisEb745cddNJEsgne4gwPXR+imV2loJeMejcE47zd33QDF1qvt3XR3w0EOJ\nyWohc4Q5OUtg7tz4t6lMBrffLrE5cz6WlEgdysiR8ce8ysok5hdLJmJRkRDR0RAYzkkSAFzfXapJ\nYM4c6f9+0UWhMQFtaubdWD5ezJsX3qgunfCrqowV2XYHAf5KpbZW6jnisQSKisQdFxQkHTBAipws\nEoOOh2QsgcmTUyNLrPBT1iUl4mpMpAVKWZnEq6JlBgFCPFddlZl4TbLIWRLwbgt54EDqLIGdOyXC\nf/75ckx3hnr++fBNa+KF9sk/GnDiif5dFTPhDho/XgjAT2nX1kq2TzwkAEgwsKoq6mkWCUAtw2RI\nIBdQUiIWTaTOvUHQVX20oLDC7Cycy8jJmAAQ2iUTSJ0lkJcnDa5efdV1Bw0cKP5ks1fI8YDTT5eq\nVC8y4Q4qKZHKSz/XmU7QeEngv/5L9ru1SD2KisR9EkvmXC5Dx3UiJKBeiVhJ4GhBzpKAuYcukDoS\nAMQ3efhwaI/+XO/0lw5oGp0XmXAHAcH3XLM24iWB4/EZZgpFRUIA+TnrO4gNankmSgJ9+sRfKZ3r\nyFkSKCkJdQelkgRmz3a7jx7PmDZNLCDdFFuRCXdQJBQWSpZWvCRgkT5Mngz89KfZliJ56LhOhASq\nqqRFTbTMoKMN0XYWW0REu5zdwryv3UpEvURUahybQkR/IaJ1RPQ+ERU4x6cTUQMRbSKiR2MRLJ2W\nQFGRlL17e54fb8jPl60dtbxekQl3UDTU1loSyCUUFso+wkc78vNlm9h4awQAcSVfe23qZco2olkC\niwFc6D1IRNUAZgHYZhzLB/AMgG8w88kAzgGgmf6/BHADM48BMIaIwj7Ti3RaAgDw4INHf5ArFfBz\nCWXKHRQJ06eHtim2sEgVnngiM7UVRwuibS/5NgC/XTl/BuA2z7HZAN5n5gbnvbuZuZeIKgCcyMyr\nnPOeBhA1hyadloCFCz8S0J782cQjjwCXXppdGSwsjgfEHRMgovkAGpn5fc9LYwAwES0jotVE9D3n\n+HAAjcZ5O51jEZFuS8BCoHGBu+8WxdvTkxuWgIWFRWYQV6yfiAoB3AlxBX1+2PndF8CXAJwK4CCA\nPxDRagA+uwUHY+HChQBk+7ZPPpkJYCYASwLpQn4+8PjjsgvSokXS/M6SgIVFbqO+vh71fql9CYCY\nOfIJRCMBLGXmyUQ0GcAbALTOtAqysp8B0dZfYebrnPfdBeAzAL8B8CYzT3COXwngHGb+ps+1WOV5\n/33Zt7fBCUlPngw8+6x0zrRID954Q6qpt20TK+xYy4KwsDhWQURg5oSSpONyBzFzAzMPZeYaZq6B\nuHlOYeZdAF4DMJmI+jtB4nMAfMDMLQA6iWgGERGAqwG8FO1a6SoWswjGl78sVZE9PZnZD8HCwiL7\niJYiugTASgBjiWgHEXn3EvrcjGDmPZCA8bsA/g5gNTO/6rx8E4AnAWwCsJmZl0UTzAaGMw8iYOFC\nKYqxhVcWFscHorqDMgnTHcQseeJdXfJ74EDZ8OV4L/DKBGLdZ9nCwiI3kDF3UCZB5FoDzNJALtnm\nbhaxwRKAhcXxg5wlAcBNEz14UKyBvLxsS2RhYWFxbCGnSUAtARsPsLCwsEgPcpoE1BKwJGBhYWGR\nHuQ8CVhLwMLCwiJ9yGkS0N3FLAlYWFhYpAc5TQJqCXz0ke34aWFhYZEO5DQJaGB46VLgq1/NtjQW\nFhYWxx5ymgRKSoBdu4Dlyy0JWFhYWKQDOU0CpaXAa6/Jnp52gxELCwuL1COnSaCkRNoaz52bbUks\nLCwsjk3kNAmUOrsXWxKwsLCwSA9ymgSGDwfmz7d7CFhYWFikCznbRdTCwsLCIjYck11ELSwsLCzS\nj2ibyiwiol1E1ODz2q1E1EtEpZ7jI4hoHxHdahybTkQNRLSJiB5NnfgWFhYWFskgmiWwGMCF3oNE\nVA3ZbH6bz3t+BuAVz7FfAriBmccAGENEYZ+Zq0jVZs6phJUpduSiXFam2GBlygwikgAzvw1gt89L\nPwNwm/cgEV0MYAuA9caxCgAnMvMq59DTAC5OVOBMIxcfupUpduSiXFam2GBlygzijgkQ0XwAjcz8\nvuf4QAgxLPS8ZThkQ3rFTueYhYWFhUWWkR/PyURUCOBOiCvo88PO74UAfs7MB4jsNuUWFhYWRwOi\npogS0UgAS5l5MhFNBvAGgAPOy1WQlf0MAP8HQLVzvBhAL4C7AbwI4E1mnuB83pUAzmHmb/pcy+aH\nWlhYWCSARFNE47IEmLkBwOddfIjoEwDTmflTAGcbx+8B0MXMjzn/dxLRDACrAFwN4H8FfL61ICws\nLCwyiGgpoksArAQwloh2ENECzymxrtxvAvAkgE0ANjPzsrgltbCwsLBIOXKqYtjCwsLCIrPIiYph\nIrqQiDY4xWTfz5IM1UT0JhF9QETriOhm53gpEb1ORBuJaDkRFWdBtjwi+jsRLc0hmYqJ6AUi+pCI\n1hPRjGzLRUR3OM+vgYieI6J+mZbJr8AykgyOzJuc8T87gzL9xHl2a4noRSIqyqRMQXIZr4UVo2br\nXjnHv+Xcr3VE9FC2ZSKi04lolaMX3iWi0xKWiZmz+gMgD8BmACMB9AWwBsCELMgxDECd8/dAAB8B\nmADgxwBuc45/H8CDWZDtOwCeBfCy838uyPRrANc7f+cDKMqmXM742QKgn/P//wZwbaZlAnAWgGkA\nGoxjvjIAmOiM976O/JsB9MmQTLP0WgAezLRMQXI5x6sBLAPwCYDSHLhX5wJ4HUBf5//yHJCpHsAF\nzt9fgSTfJCRTLlgCp0PiBFuZ+TCA3wKYn2khmLmFmdc4f+8D8CGknmEeROHB+Z3RQjciqgJwESSm\nooHzbMtUBOAsZl4EAMx8hJn3ZlmuTgCHARQSUT6AQgBNmZaJ/Qssg2SYD2AJMx9m5q2QCXt6JmRi\n5teZudf59x1Ipl/GZAqSy4FfMWrW7hWAGwE84OgnMHNbDsjUDFl4AZKNuTNRmXKBBIYD2GH834gs\nF5M5abHTIJNjKDPvcl7aBSM7KkP4OYDvQVJuFdmWqQZAGxEtJqL3iOgJIhqQTblYMtQeBrAdovz3\nMPPr2ZTJQJAMlQgtpMzW2L8ewP9z/s6qTEHFqFmWawyAs4nor0RUT0Sn5oBMtwN4mIi2A/gJgDsS\nlSkXSCCnItNO5fP/BfA/mbnLfI3F3sqYvEQ0B0ArM/8drhUQgkzL5CAfwCkAHmPmUwDshwzKrMlF\nRLUAvg0xgSsBDCSiq7Ipkx9ikCGj8hHRDwAcYubnIpyWEZnILUa9xzwc4S2Zulf5AEqY+QzIguz5\nCOdmSqZfAbiZmUcAuAXAogjnRpQpF0hgJ9wiMzh/Nwacm1YQUV8IATzDzC85h3cR0TDn9QoArRkU\n6QsA5jn1GEsAnEdEz2RZJkCeTyMzv+v8/wKEFFqyKNepAFYycwczH4EUKZ6ZZZkUQc/LO/a1+DIj\nIKLrIK7GrxuHsylTLYTE1zpjvgrAaiIammW5GiHjCc6Y7yWiwVmW6XRm/p3z9wtwXT5xy5QLJPA3\nSGfRkURUAOByAC9nWggiIgi7rmfmR4yXXoYEGOH8fsn73nSBme9k5mpmrgFwBYA/MvPV2ZTJkasF\nwA4iGuscOh/ABwCWZlGuDQDOIKL+zrM8H9LIMJsyKYKe18sAriCiAiKqgbgdVvm8P+Ug6eT7PQDz\nmfkzj6xZkYmZG5h5KDPXOGO+EcApjista3JBntd5AOCM+QJmbs+yTJuJ6Bzn7/MAbHT+jl+mVEey\nE4x+fwWSjbMZwB1ZkuFLEL/7GgB/d34uBFAKaZWxEcByAMVZku8cuNlBWZcJwFQA7wJYC1klFWVb\nLkgw8QMADZAAbN9MywSx2JoAHILEuhZEkgHi/tgMIbELMiTT9ZDCzW3GWH8skzJ55OrWe+V5fQuc\n7KAs3KvPZXLG0TPOuFoNYGaWn98CiOX7jqOv/gJgWqIy2WIxCwsLi+MYueAOsrCwsLDIEiwJWFhY\nWBzHsCRgYWFhcRzDkoCFhYXFcQxLAhYWFhbHMSwJWFhYWBzHsCRgYWFhcRzDkoCFhYXFcYz/D2YX\nMOCLkpA/AAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Slice out time series for voxel (23, 19, 0)\n", + "voxel_0_timecourse = good_data[23,19,0,:]\n", + "# Slice out time series for voxel (23, 19, 1)\n", + "voxel_1_timecourse = good_data[23,19,1,:]\n", + "# Plot both these time series against volume number, on the same graph\n", + "plt.plot(voxel_0_timecourse,'b-')\n", + "plt.plot(voxel_1_timecourse,'g-')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The scanner collected slices for these data in an \"ascending interleaved\" order. That is, the scanner first collected z slice 0, then z slice 2, up to z slice 28. It then went back to collect z slice 1, 3, 5 up to z slice 29.\n", + "\n", + "That means the scanner started collecting slice 0 in each volume, at the beginning of the TR.\n", + "\n", + "The TR is 2.5 seconds." + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# - the time between scans\n", + "TR = 2.5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make a time vector, length 172, that corresponds to the start time in seconds of each volume. This also gives the slice 0 start times." + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Make time vector containing start times in second of each volume,\n", + "# relative to start of first volume.\n", + "# Call this `slice_0_times`\n", + "slice_0_times = TR*np.arange(172)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The scanner starts to collect z slice 1 exactly half way through the TR. Make a new vector that is the start time of acquisition of slice 1." + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Make time vector containing start times in seconds of z slice 1,\n", + "# relative to start of first volume.\n", + "# Call this `slice_1_times`\n", + "slice_1_times = TR*np.arange(172) + TR/2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now plot the first 10 values for the slice 0 times, against the first 10 values of the slice 0 time series.\n", + "\n", + "To the same plot for the first 10 values of the slice 1 times, against the first 10 values of the slice 1 time series.\n", + "\n", + "Use the `:+` line marker for the plots to get the actual position of the points, and dotted lines betweeen them." + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 82, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFPWd//HXh4HhVi4FBBE8IGpIVIxoItJm4xVj1LiJ\nYaPGxCS6nnHZuNHoMpg1m5iARxKJB2rQiD9N4rUCHmh7ER0FhBHEgDLKIPfhcDPMfH5/VM1MM8xB\nz3R39fF+Ph7zoOpb1VWfadvPfPtT3/qWuTsiIlI42kUdgIiIZJYSv4hIgVHiFxEpMEr8IiIFRolf\nRKTAKPGLiBSYZhO/md1vZqvMrCyhrcTMKsxsbvhzRtjeycymmtl8M1toZj9PeM0IMyszs8Vmdkf6\nfh0REWlJSz3+B4DTG7Q5MNHdjw5/poft3wVw9y8AI4BLzWxQuG0ScIm7HwYcZmYNjykiIhnSbOJ3\n99eADY1sskbaVgBdzawI6ArsBCrNrD/Q3d1Lw/2mAOe0PmQREWmL1tb4rzKzeWY22cx6ALj7c0Al\nwR+AcuC37r4RGABUJLx2edgmIiIRaE3inwQMAY4iSPITAMzsAqAz0D/c/p9mNiRFcYqISIq0T/YF\n7r66dtnM7gOeCVe/DDzh7tXAGjN7g6DW/zowMOEQAwl6/XswM00cJCLSCu7eWAm+UUn3+MOafa1z\ngdoRP4uAr4b7dAWOBxa5+0qCWv9IMzPgQuDJpo7v7vpxZ9y4cZHHkC0/ei/0Xui9aP4nWc32+M1s\nKjAa6GNmy4BxQMzMjiIY3bMUuDTc/W5gcjj0sx1wv7u/F267HHiQoBQ0zd1nJB2piIikRLOJ393H\nNNJ8fxP77gAuaGLbbGB40tGJiEjK6c7dLBWLxaIOIWvovain96Ke3ovWs9bUh9LFzDyb4hERyQVm\nhqfz4q6IiOQ2JX4RkQKjxC8iUmCU+EVECowSv4hIgVHiFxEpMEr8IiIFRolfRKTAKPGLiBQYJX4R\nkQKjxC8iUmCU+EVECowSv4hIgVHiF5GMi5fHow6hoCnxi0jGKfFHS4lfRDJq3sp5LFm/JOowClqz\nj14UEUmFbVXbeGv5W8TL41TuqOQvZX/h0F6HAhAbHCM2OBZtgBGKl8cz/vsr8YtIWtV4Dcfeeywz\nL5pZl+D26bgPJbESANZtXRddcFkgisSvUo+IpNzHGz/mw/UfAtDO2vH2j9+mX7d+e+y3bus6vjrl\nq1RVV2U6xKyw7LNlkZxXPX4RSbnnP3yeLh26cEivQwDo0qHLbttre7i9u/TmnR+/Q4eiDpkOMVLx\n8jjx8jgzlszgreVv1bVnquylh62LSJut3rKae2bfw40n3djqY7g7Fz15ETeOupFhfYalMLrsUeM1\nLFm/hKG9h9a1lcRL6sperZXSh62b2f1mtsrMyhLaSsyswszmhj+nJ2z7gpn9w8zeM7P5ZlYcto8w\nszIzW2xmd7TmFxOR7OLu1HbUenTqQa/OvWhLx83MuPq4q+u+JeSjxesWM/b5sW16n1Kh2R6/mY0C\nNgNT3H142DYO2OTuExvs2x6YDVzg7mVm1hP4zN1rzKwUuNLdS81sGnCnu89o5Hzq8YvkiIueuIhL\nR1zKVwZ9JS3Hf+OTN9in4z4M7zs8LcfPlC07twDQtbgrEPzBNKvvnKfi4m5Ke/zu/hqwobHzNNJ2\nKjDf3cvC124Ik35/oLu7l4b7TQHO2dsARSQ7VNdUs2LTirr18bHxHD/w+LSd79NNn7Jqy6q0HT9T\n/vvl/+aJRU/UrScmfSCSoaytvbh7lZldBLwDjHX3jcBhgJvZDGA/4FF3/y0wAKhIeO3ysE1EcsjT\nHzzNix+9yB/P/CMAQ3oOSev5vn3kt+uWq2uqWb9tPft13S+t50yV7bu206l9JwB+c8pvaN8uu8bR\ntCaaScDN4fIvgQnAJUAH4ETgWGAbMNPMZgOfJXPwkpKSuuVYLEYsFmtFiCLSVlXVVfyl7C98/4vf\nx8w453PncM7novmy/vonr/On2X9i6nlTIzl/MrZVbePou4/mnZ+8Q7fibmlJ+vF4nHg83urXtziq\nx8wGA8/U1vib2mZm5wNnuPvF4bYbge3Aw8DL7n542D4GGO3ulzVyPNX4RbKEuzP2+bGUxErYp+M+\nUYdDdU01Re2Kog6jSYm1+8odlRl9z1Ja42/iBP0TVs8Fakf8PA8MN7PO4YXe0cACd18JVJrZSAve\nlQuBJ5M9r4ik372z7+WpRU8BQTKZeNrErEj6QF3S37h9I6c9fBo7du2IOKJ6f1v4N6597tq69Wx5\nz5rS0qieqQQJvA+wChgHxICjAAeWApe6+6pw/+8B14fbnnX3n4ftI4AHgc7ANHe/uonzqccvkmE7\nq3dSXFQMwLsr36Vnp54c1OOgiKNqmrszd+Vcjul/TNSh1Nm0YxNVNVX06twrkvMn2+PXDVwiBaTh\n0MFFaxfxw6d+yKxLZkUXVBs9UvYIXzv4a+zfdf+MnveCv1/AL0b9gsP3Ozyj521M2ks9IpK74uVx\n5qyYw87qnQAM6z2M6d+bHnFUbVO+sZzqmuqMn3fsCWPrZhjNNUr8IgXmzrfurJsP38zYt9O+EUfU\nNjeMuoH+3YNLj2u3rmX7ru1pOc+H6z9k7HNj69aP7n90zs4xpMQvkudmLJnBT575CSXxEsa/Mp7B\nPQbz2ILH8vIpWJPnTObe2fem5dgD9xnIV4d8NS3HzrTsuqtARFJu/67707GoY91EYG2dECybXfeV\n63BSd53w0fce5aB9D+KEA0+gY/uOnDn0zJQdO0rq8YvkoU83fcrWqq0AHNP/GH7/9d9HHFFmmBnt\nLEhrb1a8yY+e/lGbjtercy+6FXdLRWhZRYlfJA/d/MrNvPHJG3u0F9IjDkf0H8HPvvyzpF6zY9cO\nJr09qW72zFMPOTXnJ4lrjBK/SJ5Yv2193fKkMydxyiGn7LFPISX+DkUd6ub1r66p5o4372jxpq/2\n7dpTUVlR920pXynxi+SBdVvXMfrB0XWPMGw4A2Sh275rO5U7KnebN6f24nZFZQVvL38bCO4OvuVf\nbqmbQjlfKfGL5Ch3rxuP37tLb+b8ZE7ODi9Mt67FXblp9E110z6Ubyzn5aUvA7BwzcLdHn9YCHTn\nrkiO+kPpH1i7dW1ej9JJl3P/37ns32V/7j7r7qhDSQlN2SCSxxJngNyycwvFRcXq5Seh9iHn7s7N\nr97MuNHjgMw95DxdlPhF8pS7E/tzjIfOfYhB+w6KOpycl4qHnGcLzdUjkqfMjAfPfpAD9zkw6lAk\nxynxi2Sx2Z/O5opnr6hbH9JziEbspEgul3baSqUekSy2rWobH6z7gKP6HRV1KJLFVOoRyXG3vnFr\n3bjyzh06K+lLyinxp0E+znoomTNywEgG7jMw6jAkjynxp4ESvyRj4/aN3PLqLXXzw4wePLpufnmR\ndFDiT7FdNbsieRqQ5K6uHbrSrbgb1a7PjWSG5uNPkdobQ2Ytm8ULH71Qd2t4rt8YIukxZ8Ucduza\nwQkHnkCHog5cc/w1UYckBUSJv402bt/Ibf+4jZJYCbHBMXZW7+SWV2/JmxtDJHUSH3S+ZsuatD0i\nUKQlSvytUFVdRVG7ItpZO7oXd2e/rvtR7dW0t/YUFxVrnLXswd25d/a9jBo0iqJ2RZx26GlRhyQF\nTDX+VvjWY99i1rJZQDCN65XHXbnbdK+xwTGe+eAZPlj7QVQhShZ6f+37rNy8MuowRJrv8ZvZ/cCZ\nwGp3Hx62lQA/AtaEu13v7jMSXjMIWAiMc/cJYdsI4EGgEzDN3XOqoPn+mvdZu3Utow4aBcDU86Y2\n+zi22OAYj5Q9wuadmzMVomShtVvX8viCx1m1ZRUAc1fO5d45wYPAde1HotRSqecB4PfAlIQ2Bya6\n+8QmXjMReLZB2yTgEncvNbNpZnZ64h+LbLdy80qWb1rOKILEvzfP4Py34f+W7rAkyy1au4h129bt\ndr1H134kGzRb6nH314ANjWxqtIhtZucAHxH0+Gvb+gPd3b00bJoCnNOqaDNk1eZVjHpgFDVeA8DJ\nQ07mgi9c0KpjaQqKwjJjyYy6x/udOOhEbjzpxogjEtlTa2v8V5nZPDObbGY9AMysG3AdUNJg3wFA\nRcL68rAtq8z+dDbbqrYB0LdbX+476z7aWdsugbg7ox8czYfrP0xFiJIDZiyZwbLKZY1uU2lHskVr\nRvVMAm4Ol38JTAAuIUj4t7n7VmvDsJaSkpK65VgsRiwWa+2hknLP7Hu48rgrGd53OEDdQ5rbwsx4\n5LxHGNA96/7OSYqs3bqWBasXMHrwaABuP/32JvdV4pdUicfjxOPxVr++xdk5zWww8Eztxd2mtpnZ\nq0DtROE9gBrgJuDvwMvufnj4mjHAaHe/rJHjZWx2zqc/eJrVW1bzo2N+lJHzSX4qW1XG4wsf5+aT\nb255Z5E0SfvsnGHNvta5QBmAu5/k7kPcfQhwO3CLu9/l7iuBSjMbGX4TuBB4MtnzpkLiULoj9juC\nEwaekJHzVlVXMfOjmRk5l6TftMXTqNxRCcDwvsOV9CXnNJv4zWwqMAsYZmbLzOyHwG/MbL6ZzQNG\nA9fuxXkuB+4DFgNLohjRU1FZwTce+UbdxdZDex3KkfsfmZFz76jewd2z76aquioj55P0Kl1eSkVl\nRcs7imSpvHsQS+Jt8b+b9TsuPupi+nTpA0CN17T5gq0UnnVb1/HW8rf4+mFfjzoUkUYV/INYXlr6\nUt1yny592Fm9s25dSV9aY/POzbz28WtRhyGSMnnV47/zrTv5+/t/J35xPHVBpdCcFXN4aN5D3Hb6\nbVGHIi2Yvng6x/Q/hr7d+kYdikiLku3x58UkbbVTIu+s3skrH79CSbwEyL7b4of2Hqo7enPE+2vf\np3/3/kr8kpfyqscPUBIv0W3xkrQN2zbw0tKXOO+I86IORSRpBV/jzwVbq7ayavOqqMOQBLtqdvH2\np29rig0pCHnX408c1ZOtbvtHUOO/9oS9GQkr6TJjyQyG9R7GkJ5Dog5FpE2S7fHnXeLPBe6uh7Vk\ngclzJvOFvl/gSwO+FHUoIm2ixC/SQO23wM+2f8aTi57k+0d9P+qQRFJKNf4cctfbd7FwzcKWd5Q2\niZfHgeBpaR+s+6Buum2RQqXEH6ED9zmQzu07Rx1GXlm/bX3dPDoAv3n9N3zy2SdA8ACdX/3Lr3Qj\nnxS8vBjHn6vOGnZW1CHkHHenxmsoalcEwOMLHueA7gfwlUFfAWB8fDyxwTF6du5JvDzO8srlPPDu\nAwzadxCQffd2iERBiT8LbKvaRucO6vk35s2KN2ln7ThuwHEAXD/zevp168dPj/8pAD0799ztUZh3\nnHFH3XJtgh+wzwDd2yGSQN95I1bjNRx777G7TRmdT2rr600p31jOe6vfq1u/Z/Y9/Pr1X9etr926\nlnVb19Wt33zyzVwz8pq69a8d/DW+2O+LqQtYpACoxx+xdtaO0h+V0rW4a9ShpMVzS57jkJ6HcOC+\nwTN6pi+ezsI1Cxn75bEAzFs5jxWbV/D5/T8PwHmHn7dbDf4bQ7+x2/GKi4qTjkGlHZHdaTinpNWY\nv45h4D4D+e2pvwWC5yJU7qjkiP2OiDgykfyhcfw5qnJHJc/+81nGDB8TdShtFi+P8/LSlzEzxr8y\nnnGjxwG6sCqSLgU5O2c+aGftmL1iNud//vycH264q2YXq7esZtI3JgHowqpIllHizxLdirvxu1N/\nF3UYKXHSQSdxSM9Dog5DRJqQ213LPJXr5a7iouK6ic9U2hHJPkr8WWbmRzO59P8ujTqMpFXXVHPW\n1LNYsWnFbu1K/CLZRxd3s8zmnZvZtGMT/bv3jzqUpL3z6TuM6D9CM4+KZJhG9YiIFBjNzpknln22\njI83fhx1GC2aWjaV++bcF3UYIpKEZhO/md1vZqvMrCyhrcTMKsxsbvhzeth+ipm9Y2bzw39PTnjN\nCDMrM7PFZnZHY+eS3T2x6AlmLZsVdRgtOnHQiYwaNCrqMEQkCc2WesxsFLAZmOLuw8O2ccAmd5/Y\nYN+jgJXuvtLMjgSec/eB4bZS4Ep3LzWzacCd7j6jkfOp1CMikqSUlnrc/TVgQ2PnaWTfd929dqax\nhUBnM+tgZv2B7u5eGm6bApyztwFK9tm0YxOX/d9l7Ni1I+pQRKQVWlvjv8rM5pnZZDPr0cj284DZ\n7l4FDAAqErYtD9tkL1z3wnVZ95SuLh26cOohp9KxfceoQxGRVmjNnbuTgJvD5V8CE4BLajeGZZ5f\nA6e0JqCSkpK65VgsRiwWa81h8sZZQ8+if7fsGtpZ1K6Ibx3+rajDEClY8XiceDze6te3OJzTzAYD\nz9TW+JvbZmYDgZnAxe7+j7CtP/CSux8ero8BRrv7ZY0cTzX+LPbguw/yuT6f4/iBx0cdiogkSPtw\nzjCR1zoXKAvbewDPAv9Vm/QB3H0FUGlmIy24s+dC4Mlkz1vosuFBLQd0P4D9uuwXdRgi0kYtjeqZ\nCowG+gCrgHFADDgKcGApcKm7rzKzG4GfA4sTDnGKu681sxHAg0BnYJq7X93E+dTjb0RVdRXH3HMM\nr//gdfbttG/U4YhIltGdu3mqxmsima553dZ1TJk3hZ8e/1NNxSCSpXTnbp6Kao7+aq+mc4fOSvoi\neUQ9/hyyYtMKHlvwGNccf03LO4tIwVCPP49179i99j9w2s/16HuPsnrL6rSfR0QyT0/gyiHdirtx\n9chGr4un3CeffcL2Xdszci4RySyVenJUVXUVHYo6RB2GiGQBlXoKwOMLHufq6anv+a/avIqZH81M\n+XFFJLuox5+DtlVtw8zo1L5TSo/7zqfv8Er5K4z98tiUHldE0kvj+EVECoxKPQXkrYq3WLxuccs7\ntuCFD1+guqY6BRGJSC5Q4s9hZavL+OSzT9p0jF01u3ho/kOs37Y+RVGJSLZTqUdEJMep1CN7paKy\ngorKipZ3FJG8o8SfB87/6/m8v+b9pF7z0tKXeHKRZscWKUQq9eSBBasXMKzPMNq3043YIoVIpZ4C\ndOT+R+510s+25/eKSOYp8ecJd+fdle82u8+GbRu4/NnL2bFrR4aiEpFspFJPntiycwvfmPoNpv3b\nNDp36Nzkfu6uufVF8ozu3JU9rNi0gt5delNcVBx1KCKSBqrxyx5ufeNWnvngmajDEJEsoR5/nvlg\n7Qc8UvYI408eT7w8TmxwLLLn9YpIZqjHX+AG7DOAEQeMAGDa4mlAdM/rFZHspB5/nlq4ZiFf/8vX\nWXrNUl3MFclzyfb4dcdPnomXx4mXxwH4+LOPGf/KeABig2PEBseiC0xEskazid/M7gfOBFa7+/Cw\nrQT4EbAm3O0Gd58ebrse+CFQDVzt7s+H7SOAB4FOwDR3vyblv4kAeyb4klhJZLGISHZqqfj7AHB6\ngzYHJrr70eFPbdI/AjgfOCJ8zV1WX2OYBFzi7ocBh5lZw2OKiEiGNJv43f01YEMjmxqrJZ0NTHX3\nKncvB5YAI82sP9Dd3UvD/aYA57Q+ZNlbKu2ISGNaO9zjKjObZ2aTzaxH2HYAkDjPbwUwoJH25WG7\npJkSv4g0pjUXdycBN4fLvwQmAJekKqCSkpK65VgsRiwWS9WhRXJePA76X0Li8TjxeLzVr29xOKeZ\nDQaeqb2429Q2M/s5gLv/Otw2AxgHfAy87O6Hh+1jgNHuflkjx9NwTpFmlJQEPyKJ0n4DV1izr3Uu\nUBYuPw1818yKzWwIcBhQ6u4rgUozGxle7L0Q0BNARFph507YuDHaGNrQ0ZQs0dJwzqnAaKCPmS0j\n6MHHzOwogtE9S4FLAdx9oZk9BiwEdgGXJ3TfLycYztmZYDjnjDT8LiJ56aWX4NVXg+X//d9g/fTT\ng5LPp59CaSncfnuwfcYMeP99uPbaYL20FD75BP71X4P1Dz+E9evhS18K1jdsCP6Y9O279/Go3JT7\nmk387j6mkeb7m9n/V8CvGmmfDexRKhKR5lVWBkm8tBQ6dgzaEks927fDmWfWr3/+8zAgYehEx47Q\ntWv9enk5LFlSn/iffx7KyuB//idYv/feYP3OO4P1J54I9v/Zz4L1V1+F+fNT+RtKFDRlg0gWqq6G\noqJg+dNP4YADguV01/irq6GqCjp1CtZXroTNm6GiIujpr1wJd98N48YF22Mx9f6zgaZsEMlxEyYE\nCfi664L12qQP6U+yRUX1f3AA+vUL/j300Ppz9+unC8y5TolfJAvs2gXtw/8bf/hD6NKl8f2yqXf9\ni1/ABRfA4YdHHYkkS/P1NqARC5Jp27fDF78ImzYF6z171tfzs1HtH59YbPdvI5I7lPgbUOKXTKmq\nCv7t1Cn43HXvHmk4e6028Z9yCuy7b6ShSCsp8SeoqYEtW6KOQgrBn/4E48fXr++3X3SxtNWqVXD1\n1cF1CckNqvET9LbicVi7Fv74x6C+aqYRC5JaO3dCcfi8+wsugA4doo0nVXr3htGjd78oLNlNwzkb\nGDdu956YSCrs2hXU8V95Bfr0iToayTd65m4rvPAC1P69qX2CwJo1wUU3kbbYsSP4t317eOON/E/6\nTz0FU6dGHYW0pOAT/44d8OCD9fOf1JZ2brkFnn02qqgkHzz8cP1YfIAePZreN18ceigccUTUUUhL\nVOppQk0NtCv4P4uSrG3boHPnYHn79uAbZDYPzZT8oFLPXnrqqWA0QlMSk/66demPR3KfO5xwAixb\nFqx36lS4Sd89mN9nxYqoI0mvXB3+XbCJf/HiYJbClmzcCCefXF+rFWlo69bgXzN4/XU48MBo48kG\nZnD88fk/zl+JP8f853/u3a3mPXrAO+8Ubs9NmvfUU3DllfXr3bpFF0u2Oe+8pqeeyAe1N+DlooKq\n8cfjQe/s619v3eurq4Ox2LU1XCkciXPQb95cn+B37Qo+F+oYNG3p0mAARb4Mk66972fePHjyyeyY\nqVSzczajrQn7j3+Ezz6Dm25KTTySOxIT/2mnweTJ8LnPBcM02xfU/0XJ69MHhufJ0zjc6xO8e/DH\nLBdnKi2oj+zIkW17/WWX1Y/zl8KSeI3nxRf1rS8Z3bvXPwEs1116afC7nHpqbueCvE/8//wnPPBA\n8Mi6tqq93R52v/1e8lPtV/qPPw5KFbXlHE3l0XoPPQQDBwYDJnLRTTdB/4Snjufq5yDvE3///sE8\nIqlUXR2MWHj22d0/BJJfEhP8gAG5+ZU+2xx0UHLP941aVVXQabzuumB4bsMRW0r8Wap79+DB1KlU\nVBRM89C7d2qPK9lLdfzUOOmkqCNITvv20KtXcBE/n+TlcM4tW+Css4LRF+mSmPRratJ3HonOnDnB\nBX3I3Z5dttq5E/7rv+ofPpNNampg4cJg2SwYrptvw3TzMvF37Qo33JCZ/1gvvAA//nH6zyOZ16sX\nDB0aLCvxp1aHDsGoqGwcBrt4Mdx4Y/3Ejfkor8bxu2f+SntVVTCPv2r9Irktl+fnSulcPWZ2v5mt\nMrOyRraNNbMaM+sVrncys6lmNt/MFprZzxP2HWFmZWa22MzuSOYXSsa3vw2zZ6fr6I3r0KE+6WfR\n31Bpg40bgz/mkhkLFsAf/hBtDA88UH8jViFo6e/bA8Ael0bN7EDgFODjhObvArj7F4ARwKVmNijc\nNgm4xN0PAw4zsxRfbg389rdw9NHpOPLe+d73Mv+HR1LvxRfh1lujjqJw9OwZjJqK0nnnBZPKFYoW\nSz1mNhh4xt2HJ7Q9DvwSeAoY4e7rzew04ArgXKAn8AYwEugMvOTuh4ev/S4Qc/fLGjlX0qWeqqrg\n61k2PPbtn/+EQw7JjlikbaIoG0pm/eEPwV3Yhx0WdSRtl/Zpmc3sbKDC3ecntrv7c0AlsAIoB37r\n7huBAUBFwq7Lw7aUuPVW+P3vU3W0thk6VEk/XyjpR2PSpKD0kwkDBhTuTZhJjU42sy7ADQRlnrrm\ncNsFBL37/kAv4DUzm5lsQCUJd8nEYjFiLQynuPba7Lsgs2RJcIffI48ogeSShx4KhgJftsd3UcmU\n/v3TOxpv7tz6cvC556bvPOkWj8eJt2FO6KRKPWY2HHgRCGcgZyBBD34kMA6Y5e4Ph6+bDEwHXgde\nTij1jAFGt6XU4w4bNgTD7bJRdXXwATv22KgjkWSsWgWVlfnx1V/2tGMHfPOb8OijwXWFfJLWUo+7\nl7l7X3cf4u5DCEo4x7j7KmAR8NUwiK7A8cAid18JVJrZSDMz4ELgyWTO29Drr8NPftKWI6RXUZGS\nfi7q21dJP1ts3hxMkZGKO2arq4N/O3aE557Lv6TfGi0N55wKzAKGmtkyM/tBg10Su+d3A8Xh0M9S\n4H53fy/cdjlwH7AYWOLuM9oS9KhRQRklF0yYADPa9NtKui1aFEzEJtmjY0fYf/+2H6esDM44o+3H\nyTc5dQPXkiVw6KEZDCgF5s4NLiKl4kMs6fHnPwfXYi66KOpIJNXcgxJev35RR5JeyZZ6cibxb9wY\nPDnr5Zez8zZvEUmP0lIoL4fvfGfv9n/++WAOoPPOS2tYWSXtwzmj0qMHvPFG7ib9ykq44go9tD2b\nZFGfR5rRtWtyI3323x8OOCB98eSDrE/8b71Vf4Enl4dGdusWPBegUMcNZ6Prr4eHH446CmnJkUe2\n/JzsefNg27Zg+aij4IQT0h9XLsvqxO8e3NDx0UdRR9J27doFX1Vz+Y9Xvrn++mD6bskdt98On34a\nLCcOY7/33uBCruydrH68hFnwyLt8M3NmcM2ikGqQ2WjffaOOQJK17771N2zOnFk/XXbUk7zlmqy8\nuFtaGjzo5JBDoo4oPd59NxinfOKJUUdSmKZPh4MPhmHDoo5EWmv7dhg0KBjpt88+UUcTvWQv7mZl\nj/+994IhkPma+I86KuoICtvKlbn13FepF4/Xl3jWrIGJE4PlxOcjS8uyssdfKGpq4LbbgrlhunaN\nOhqR3FJSEvxIHg/nzEdmQb2yqirqSApDOp/BLJJLsq7HP25cEI++ukmqjRkDP/gBnHpq1JFIKsTj\nyhG18vbO3Xy3bBnMmQNnnx11JPlr+/bgPopsm8ZbpK1U6slRmzbBJ5/Ur7dhqm1pQqdOSvoioMSf\nNY44Aq4ltVHlAAAJO0lEQVS6qn5diT91JkyADz6IOgqR7KHEn4WmTw/KEpIaAwfCfvtFHYVI9sjK\ncfyFqnaMcjwOr7wSlCZAF7rb6vzzo45AJLvo4m6W0hjltlu6FAYP1vxIkv90cTeP7NoVTEqlcf6t\nc8UVwdO1RGR3KvVkqVgsuLN3y5bgD0CHDlFHlHuefVa9fZHGqNQjIpLjVOrJU+Xl8MtfRh1F9nOH\n730vuCFORBqnxJ8jevbUNMJ7wyyY9G7AgKgjEcleKvWIiOQ4lXoKwF//Cn/6U9RRZJ94vP75zCLS\ntGYTv5ndb2arzGyPp1ma2VgzqzGzXgltXzCzf5jZe2Y238yKw/YRZlZmZovN7I7U/xqF5bjjdENX\nQ7t2wV13BXMeiUjzmi31mNkoYDMwxd2HJ7QfCNwLDANGuPt6M2sPzAYucPcyM+sJfObuNWZWClzp\n7qVmNg24091nNHI+lXqS5K4hiyKFLqWlHnd/DdjQyKaJwHUN2k4F5rt7WfjaDWHS7w90d/fScL8p\nwDl7G6A079//HV59NeoooqW+gkhykq7xm9nZQIW7z2+w6TDAzWyGmc02s5+F7QOAioT9lodtkgL/\n8R9w/PFRRxGdbdvg2GNV4hFJRlJ37ppZF+AG4JTE5vDfDsCJwLHANmCmmc0GPkvmHCUJE9TEYjFi\nKmY3a+jQ+uVdu6B9gd2L3bkz/P3v0L171JGIZE48HifehrnbWxzOaWaDgWfcfbiZDQdeBLaGmwcS\n9OBHAjHgDHe/OHzdjcB24GHgZXc/PGwfA4x298saOZdq/K20cyeccAI89xz06RN1NCKSSWkdzunu\nZe7e192HuPsQghLOMe6+CngOGG5mncMLvaOBBe6+Eqg0s5FmZsCFwJPJnFdaVlwczE1TKEm/pgbu\nuUcT2Im0RkvDOacCs4ChZrbMzH7QYJe67rm7byS46Ps2MBeY7e7Tw82XA/cBi4EljY3okbbr169+\nefPm6OLIhM2bg2ksioqijkQk9+jO3Ty0eDFcfDG8/rqGeooUgmRLPUr8eWrLFujaNeoo0iOffzeR\n1tCUDQLUJ8YdO2DNmmhjSaWVK4Phq9XVUUcikruU+PPc44/DhAlRR5E6/frBW2+pti/SFir15Dn3\nYASMEqVI/lKpR3ZjVp/0Fy8O6uO5aMMG+MUvND2DSCoo8ReQu++GN96IOorWO/JIjVISSQWVeiTr\nxeOahlqkOSr1yF6ZMSM3RsbU1MBTT0UdhUh+UeIvQNXVwWifXBjmuWBBMBWFiKSOSj2SdR55BCor\ngzH7AOPHw7hxwXIsprKPSEPJlnoKbBJfaWjrVnjssWCKh6j87W/ByKNzwsfz9O4NX/xicDG3VsJs\n3SLSRkr8BW7bNvjoo6CW3i5Dhb/p04OhpVdfHawffPDu9xmcdlpm4hApVCr1SNrNmhV8q7j99mD9\n44+DJ2Z9/vN793qN6hFpniZpk1abOxfmzWtd2Sfxoe+LFsHYsfUXZTdsgNWrYdiwlIUqIgk0nFNa\nrXt36NWrfr25J7vt3Fm/vG5d8AjI2r/ZBx8Md91Vv71nTyV9kWyixC91Dj0UvvnN+vXExL9xY3Ad\nAILhoIMGBReGIbgY++ab9T3+4mI46KCMhCwiraDEL436/e+Di761Tj4ZPvkkWC4qCur0XbrUb+/d\nO7PxiUjraVSP7CYeD36WL4eHHgrKNgATJ8LgwfX7dewYQXAikhK6uCtNKinR+HmRXKCLuyIi0iwl\nfmmSxs6L5CeVekREclxKSz1mdr+ZrTKzska2jTWzGjPr1aB9kJltNrOxCW0jzKzMzBab2R17G5yI\niKReS6WeB4DTGzaa2YHAKcDHjbxmItBwIt1JwCXufhhwmJntcUzZXby5u6cKjN6Lenov6um9aL1m\nE7+7vwZsaGTTROC6ho1mdg7wEbAwoa0/0N3dS8OmKcA5rQ24UOhDXU/vRT29F/X0XrRe0hd3zexs\noMLd5zdo70bwx6CkwUsGABUJ68vDNhERiUBSN3CZWRfgBoIyT11z+G8JcJu7bzXTI7FFRLJVi6N6\nzGww8Iy7Dzez4cCLQDhLCwMJevAjgceBA8P2HkANcBPwd+Bldz88PN4YYLS7X9bIuTSkR0SkFdL2\nBC53LwP61q6b2VJghLuvB05KaB8HbHL3u8L1SjMbCZQCFwJ3tjVwERFpnZaGc04FZgFDzWyZmf2g\nwS5720O/HLgPWAwscfcZSUcqIiIpkVU3cImISPplxZQNZna6mS0Kb/D6r6jjiZKZlZvZfDOba2al\nLb8ifzR2w6CZ9TKzF8zsn2b2vJn1iDLGTGnivSgxs4rwszG3UO6HMbMDzexlM1tgZu+Z2dVhe8F9\nNpp5L5L6bETe4zezIuAD4GsEF4rfBsa4+/uRBhaRBtdNCoqZjQI2A1PcfXjYdiuw1t1vDTsFPd39\n51HGmQlNvBe1184mRhpchplZP6Cfu78bDhufTXAv0A8osM9GM+/Fd0jis5ENPf7jCOr+5e5eBTwK\nnB1xTFEryIvcTdww+E3gz+HynymQm/+auXmy4D4b7r7S3d8NlzcD7xPcC1Rwn41m3gtI4rORDYl/\nALAsYb2Cwr7By4EXzewdM/tx1MFkgb7uvipcXkXCqLICdZWZzTOzyYVQ2mgoHF5+NPAWBf7ZSHgv\n3gyb9vqzkQ2JX1eXd/cVdz8aOAO4IvzKL0A4dWshf14mAUOAo4AVwIRow8mssLTxN+Aad9+UuK3Q\nPhvhe/FXgvdiM0l+NrIh8S+n/sYvwuWKJvbNe+6+Ivx3DfAEQSmskK0K65q18z6tjjieyLj7ag8R\nDI8umM+GmXUgSPoPufuTYXNBfjYS3ouHa9+LZD8b2ZD43yGYsXOwmRUD5wNPRxxTJMysi5l1D5e7\nAqcCe0yJXWCeBr4fLn8feLKZffNamNxqnUuBfDbCKWAmAwvd/faETQX32WjqvUj2sxH5qB4AMzsD\nuB0oAia7+/9GHFIkzGwIQS8fgruq/1JI70V4w+BooA9Bzfa/gaeAx4BBQDnwHXffGFWMmdLIezEO\niBF8lXdgKXBpQo07b5nZicCrwHzqyznXE8wEUFCfjSbeixuAMSTx2ciKxC8iIpmTDaUeERHJICV+\nEZECo8QvIlJglPhFRAqMEr+ISIFR4hcRKTBK/CIiBUaJX0SkwPx/IupEaiJjL+sAAAAASUVORK5C\nYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot first 10 values of slice 0 times against first 10 of slice 0 time series\n", + "# Plot first 10 values of slice 1 times against first 10 of slice 1 time series\n", + "# Use ':+' marker\n", + "plt.plot(slice_0_times[0:10],voxel_0_timecourse[0:10],'b:+')\n", + "plt.plot(slice_1_times[0:10],voxel_1_timecourse[0:10],'g:+')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Import `InterpolatedUnivariateSpline` from `scipy.interpolate`. Make a new linear (`k=1`) interpolation object for slice 1, with the slice 1 times and values." + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Import `InterpolatedUnivariateSpline` from `scipy.interpolate`\n", + "from scipy.interpolate import InterpolatedUnivariateSpline as IUS\n", + "# Make a new linear (`k=1`) interpolation object for slice 1, with slice 1 times and values.\n", + "interp = IUS(slice_1_times, voxel_1_timecourse, k=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Call the object you got with the slice 0 times, to get the estimated time series values for slice 1, if slice 1 had been collected at the same time as slice 0:" + ] + }, + { + "cell_type": "code", + "execution_count": 84, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Call interpolator with `slice_0_times` to get estimated values\n", + "voxel_1_timecourse_est = interp(slice_0_times)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Repeat the plot of the first 10 values of the time series. This time, on the same plot, add the estimated values for slice 1, if they had been collected at the same time as slice 0. Use a black `x` for the estimated points (marker `kx'`):" + ] + }, + { + "cell_type": "code", + "execution_count": 85, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 85, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPX18PHPIRBZNWwKAhqqQF1QFC1SQMb2cWnVqvXX\nWlpttdhq3VofilXREvUX2xrBrZW6IFar+NQuogJRREcEqsgiRAEFTYQghC0h7NnO88e9SSYh2yQz\nc+/MPe/XKy/mfu+dOyfjeObm3O8iqooxxpjgaOd1AMYYYxLLEr8xxgSMJX5jjAkYS/zGGBMwlviN\nMSZgLPEbY0zANJn4ReQZESkSkbyItiwRKRSRFe7Pd9z2jiIyU0RWichqEbk94jnDRSRPRNaJyCPx\n+3WMMcY0p7kr/hnABfXaFJiqqqe5P3Pd9h8BqOopwHDgOhE5xt03DRivqoOAQSJS/5zGGGMSpMnE\nr6rvAcUN7JIG2jYDXUQkDegClAGlItIX6KaqS9zjngMubX3Ixhhj2qK1Nf6bRWSliEwXkQwAVX0D\nKMX5AigAclS1BOgHFEY8d5PbZowxxgOtSfzTgIHAMJwkPwVARK4EOgF93f2/FZGBMYrTGGNMjLSP\n9gmqurX6sYg8Dbzmbn4T+I+qVgLbRGQRTq1/IdA/4hT9ca76DyEiNnGQMca0gqo2VIJvUNRX/G7N\nvtplQHWPn7XAt9xjugBnAWtVdQtOrX+EiAhwFfBKY+dXVftRZfLkyZ7H4Jcfey/svbD3oumfaDV5\nxS8iM4GxQC8R2QhMBkIiMgynd08+cJ17+BPAdLfrZzvgGVX92N13A/AsTilojqrmRh2pMcaYmGgy\n8avquAaan2nk2IPAlY3sWwYMjTo6Y4wxMWcjd30qFAp5HYJv2HtRy96LWvZetJ60pj4ULyKiforH\nGGOSgYig8by5a4wxJrlZ4jfGmICxxG+MMQFjid8YYwLGEr8xxgSMJX5jjAkYS/zGGBMwlviNMSZg\nLPEbY0zAWOI3xpiAscRvjDEBY4nfGGMCxhK/McYEjCV+Y0zChQvCXocQaJb4jTEJZ4nfW5b4jTEJ\ntXLLStbvXO91GIHW5NKLxhgTC/vL9/PBpg8IF4QpPVjKC3kvcHyP4wEIZYYIZYa8DdBD4YJwwn9/\nS/zGmLiq0irOeOoM5v90fk2CO/yww8kKZQGwY98O74LzAS8Sv5V6jDEx92XJl3y+83MA2kk7PvzF\nh/Tp2ueQ43bs28G3nvsW5ZXliQ7RFzbu2ujJ69oVvzEm5t78/E06d+jMcT2OA6Bzh8519ldf4fbs\n3JOlv1hKh7QOiQ7RU+GCMOGCMLnrc/lg0wc17Ykqe9li68aYNtu6dytPLnuSu86+q9XnUFV++spP\nuWvMXQzpNSSG0flHlVaxfud6BvccXNOWFc6qKXu1VkwXWxeRZ0SkSETyItqyRKRQRFa4PxdE7DtF\nRP4rIh+LyCoRSXfbh4tInoisE5FHWvOLGWO8M3v2bEpKSuq0FRcX8/rrrwOQ0TGDHp160JYLNxHh\nlm/cUvNXQipat2MdE96c0Kb3KRaaq/HPAC6o16bAVFU9zf3JBRCR9sDzwC9V9WRgLFDhPmcaMF5V\nBwGDIr8sjDH+N2rUKCZNmlST/EtKSjjrx2fR/linWpyels4NZ96ASIsvOht0Zr8zad/OOeeiDYvI\nK8pr5hn+t7dsL3vL9gIwpNcQXv3Rq3XeJy96NDWZ+FX1PaC4gV0N/dc9D1ilqnnuc4tVtUpE+gLd\nVHWJe9xzwKVtiNkYk2AZGRnce9+9/GbibygoKGDSpEn8Y9o/OPekc+P2ml/t/oqivUVxO3+i/P6d\n3/Oftf+p2a7/5ehF4m/tzd2bReSnwFJggqqWAIMAFZFcoDfwkqrmAP2AwojnbnLbjDFJZEHRAipH\nVjJw4EDy8/PJzMyM6+v94KQf1DyurKpk5/6d9O7SO66vGSsHKg7QsX1HAP507p9q/orxi9ZEMw24\n1318HzAFGA90AEYDZwD7gfkisgzYFc3Js7Kyah6HQiFCoVArQjTGtFV5ZTkv5L3Az079GSJCqE+I\neUvnkZ+fT05ODtnZ2WRkZCQkloUbFvLXZX/lyo5XMmrUqDqvW1JSwqJFi7jwwgsTEktz9pfv57Qn\nTmPpL5fSNb1rXJJ+OBwmHA63/gSq2uQPkAnkNbcPuAJ4NmLfXcBvgT7Amoj2ccBfGzmfGmP8oaqq\nSm/NvVV3HdilxcXFesMNN2hxcbGq6iHbiVBRWeGLOF5//fVDXq+4uFhfe+21mu1dB3YlLB5VVTd3\nNpvPq3+iTvxA34jHtwIvuo+7A8uATjh/ScwDvuPu+wAYgXNvYA5wQSOvlYC3yBjTmCeXPqmvrHnl\nkPbGkt3rr7+eqNBqFGwu0AHfHqBr161NeNJXPfTLpri4WM//0fl6/cvXJzSOSNEm/ib78YvITJze\nOb2AImAyEAKG4fTuyQeuU9Ui9/ifAHe4+2ar6u1u+3DgWfdLYY6q3tLI62lT8RhjYq+ssoz0tHQA\nPtryEd07dufYjGM9jqpxqsrsD2dz8YiLE3KvoSElJSVMmjSJiRMnkpOTw+2/v50uh3ehR6ceCY8F\nou/H32TxSVXHNdD8TBPHvwC80ED7MmBoS4MyxsRH/Xlh1m5fy89n/ZzF4xcDMKzPMI8ia7ldu3Yx\n929za+41DPvxMC4ZdglHdjkyYTFkZGSwaeimmhvdA44akLDXjgWbq8eYAAkXhFm+eTlllWUADOk5\nhLk/metxVC1XfaWdnZ1NZmYm2dnZTH9wOjt37kx4HOnvp/PZ+s/Iyck5ZHCb31niNyZgHv3g0Zr5\n8EWEIzoe4XFELbdo0aI6vYkyMjLInZHL56ucCeG279vOgYoDcXntz3d+zoQ3JtR8+Tz58JMMOm4Q\n2dnZdQa3JQObq8eYFJe7Ppd/r/k3R3c7mnvevYfJYycDqTkP/p8W/onOHTpz84ibY37ugxUHeeuL\nt2AdvutSGm2N3xK/MSlu+eblzFgxg8e++1hMJgTzM1VFUdpJbIoZL338EscecSwjB4yMyfniJaaT\ntBljktNXu79iX/k+AE7vezqPffcxjyNKDBGpSfrvF77Pta9e26bz9ejUg67pXWMRmq9Y4jcmBd37\n7r0s2rDokPZUK+00ZXjf4Uz85sSonnOw4iDTPpxWM3vmecedx9CjUq9DoiV+Y1LEzv21PVumXTiN\nc487dAK1ICX+Dmkdaub1r6yq5JH3H+FgxcEmn9O+XXsKSwtr/lpKVZb4jUkBO/btYOyzY2uWMGzr\n9Mip5kDFAUoPltaZNydcEAagsLSQDzd9CEBauzSyv51Nl/QuXoSZMJb4jUlSqlrTH79n554s/+Xy\nwC1h2FJd0rtw99i7SWuXBkBBSQHv5L8DwOptq+ssfxgElviN8bmGVr8qKSnhVw/9ivvfu7+mzZJ+\ny936xq1s2bMFcOr4N33jJo8jSix/TRJtjDlE9epX2dnZHHHEEezatYtJkyZx9+S7ObJn4qYpSAXV\ni5yfcuQp3LvgXvp26wuk5piGplg/fmOSQElJCXfeeSfLvraMIZ8O4dGcRxM2F36qSqUxDTGdpM0Y\n4w8ZGRncdtttDBw4kJe+eMmSvmkTq/Eb42PLvlrGjbNvpKSkhJycHPLz83nwwQeTal4YvwpSaac+\nK/UY42P7y/ez9IulvPToSzWTk0XOUGlX/gZsrh5jkt4Dix7gnMxzOLPfmYDTq8dvk4IZf7HE7wP1\nF7swJhrvFrzL4J6Da3qcGNMcm6TNB6pHBBrTEiUHSshekF0zP8zYzLGW9E1cWeKPsYqqCiqrKr0O\nwySRLh260DW9K5VqnxuTGNadM0aqB4Ys3riYeV/MqxkaHrSBIaZllm9ezsGKg4wcMJIOaR349Vm/\n9jokEyCW+Nuo5EAJD/33IbJCWYQyQ5RVlpG9IDtlBoaY2Im897Nt77a4LRFoTHMs8bdCeWU5ae3S\naCft6Jbejd5delOplbSX9qSnpdvMiOYQqspTy55izDFjSGuXxvnHn+91SCbArMbfCt//x/dZvHEx\n4EzjetM3bqoz3WsoM8Rrn77Gp9s/9SpE40Nrtq+pmRjMGC81ecUvIs8AFwJbVXWo25YFXAtscw+7\nQ1VzI55zDLAamKyqU9y24cCzQEdgjqomVUFzzbY1bN+3nTHHjgFg5uUzm1yOLZQZ4sW8F9lTtidR\nIRof2r5vOy9/8jJFe4sAWLFlBU8tfwqwez/GW82VemYAjwHPRbQpMFVVpzbynKnA7Hpt04DxqrpE\nROaIyAWRXxZea26AzJY9W9i0exNjcBJ/S9bg/PHQH8ctXpMc1m5fy479O+rc77F7P8YPmiz1qOp7\nQHEDuxosYovIpcAXOFf81W19gW6qusRteg64tFXRxkn1tLfV8598VvgZJ37/REZ+cyQA5ww8hytP\nubJV506FAWmm5XLX59Ys7zf6mNHcdfZdHkdkzKFaW+O/WURWish0EckAEJGuwG1AVr1j+wGFEdub\n3DbfyMjI4Ac3/YDf3fE7CgoKeOQPjzDrqVn06N6jTedVVcY+O5bPd34eo0iN3+Wuz2Vj6cYG91lp\nx/hFa3r1TAPudR/fB0wBxuMk/IdUdZ+0oVtLVlZWzeNQKEQoFGrtqaIyc/1MLvv5ZQwcOJD8/Hwy\nMzPbfE4R4cXLX6RfN199z5kY2r5vO59s/YSxmWMBePiChxs91hK/iZVwOEw4HG7185udq0dEMoHX\nqm/uNrZPRBYAA9xdGUAVcDfwb+AdVT3Bfc44YKyqXt/A+RI2V8+rn77K1r1bufb0awFqZjycOHEi\nOTk5NvOhaZG8ojxeXv0y955zb/MHGxMncZ+rx63ZV7sMyANQ1bNVdaCqDgQeBrJV9XFV3QKUisgI\n9y+Bq4BXon3dWIjsSndi7xMZ2d+p4UdOc5uZmUl2dnadmn9blVeWM/+L+TE5l0mcxta6zXoyi9KD\npQAMPWqoJX2TdJpM/CIyE1gMDBGRjSLyc+BPIrJKRFYCY4FbW/A6NwBPA+uA9V706CksLeSiFy+q\nudl6fI/jOenIkwBYtGhRnSv8jIwMsrOzWbRoUUxe+2DlQZ5Y9gTlleUxOZ9JjPo3/asvEA4cfYDC\n0sJmnm2Mf6XctMyRw+IfXPwgVw+7ml6dewFQpVW0ExuzZlqupKSECb+bwOgfjWbpP5daCdD4UuDX\n3H07/+2axN+rcy/KKstq9lnSN9HKyMjgmhuvYcypY8jPz7ekb1JCSmXCRz94lAVfLqjZvnrY1Rzd\n7WgPI6pr+ebl3JrbksqY8drcdXMp2lNESUkJM5+YSX5+Pjk5ObbWrUkJKXHFXz0lclllGe9++S5Z\n4SzAf8PiB/ccbCN6k8Sa7WvoWtW1zlq31Tf9rdxjkl3K1fizwlk2LN5ErXh/MW/nv83lJ15e02Zr\n3ZpkYUsvJoF95fso2lPkdRgmQkVVBR9+9WGdKTYuvPDCQ67sMzIyLOmbpJdyV/zJsND5Q/99CIBb\nR1q930u563MZ0nMIA7sP9DoUY9ok2it+Xyb+VP9zWlVtsRYfmL58OqccdQpn9jvT61CMaZOkL/VU\nD5IZNWqU16HEjSX9xAoXhAHYdWAXf/vobzXt408fb0nfBJLvEn+Qek08/uHjrN62uvkDTZtUJ/60\ndml8uuNTqrTK24CM8ZjvEv/EiRMDkfQBBhw+gE7tO3kdhi81Nk/O7Nn11/ipa+f+nTXz6AD8aeGf\n2LBrA+AsoHP/t++3gXwm8Hz3f0CQBslcPORiu7HYiMbmyfnmN79JZVVlzXEvf/IyizbUzql0T/ge\n5n8xn3BBmKxwFut3rmfGRzOcbr7hrJqrf2OCzHc3d4uLiwNV7gHYX76fTh3syr++kpISxv96PFf9\n6irmPT+P7Oxs/rj0j/Tp2offnPUbAN764i16d+7NqX1ObfQ8NrbDpDrr1ZNkqrSKodOGMv+n8+nT\ntY/X4cRcc91rC0oK2FO2h5OPPBmAJ5c9yc79O7l99O0ATH97Otd++9qaxXHKKsvo0K5DVDfILfGb\nVJf0vXogWINk2kk7lly7JCWTPsAb699g467apQjnrpvLlMVTarZXblnJwg0La7YvP+Fyrht+HeBc\n8S//1/I68+Skp6VH3SvK7+M6jEk0X17xm9Qx7p/j6H94f3LOywGcdRFKD5ZyYu8Tm3xe5OI4GRkZ\nh2wbY2qlRKkniEoPljL7s9mMGzrO61DaLFwQ5p38dxAR7nn3HiaPnQxEN2mezZNjTMsFfj7+ZNVO\n2rFs8zKuOPmKpO9uWFFVwda9W5l20TSAVtXXG0ruQSoBGhNPlvh9omt6Vx4870Gvw4iJs489m+O6\nH+d1GMaYRiT3pWWKqR60FFnuasmgJb9JT0uvGZ9gN1aN8R9L/D4yatQofnbzz7j6pauB5Jq3qLKq\nkotnXszm3ZvrtFviN8Z/LPH7SEZGBn+Z8hfkbaGgoCCperGktUtj8tjJKdst1ZhUYr16fKigoICB\nAwfWDFoyxpimpMQAriArKSkhJyeHRSsX8fv//b3v5y2amTeTp5c/7XUYxpgoNJn4ReQZESkSkbyI\ntiwRKRSRFe7PBW77uSKyVERWuf+eE/Gc4SKSJyLrROSR+P06yS1ykNLS/Us5+5qz60xU5kejjxnN\nmGPGeB2GMSYKTZZ6RGQMsAd4TlWHum2Tgd2qOrXescOALaq6RUROAt5Q1f7uviXATaq6RETmAI+q\nam4DrxfoUo8NWjLGtEZMSz2q+h5Q3NDrNHDsR6q6xd1cDXQSkQ4i0hfopqpL3H3PAZe2NMAgSZbF\nvXcf3M31r1/PwYqDXodijGmF1tb4bxaRlSIyXUQa6nJyObBMVcuBfkBhxL5Nbptpgdvm3ea7Vbo6\nd+jMecedx2HtD/M6FGNMK7Rm5O404F738X3AFGB89U63zPNH4NzWBJSVlVXzOBQKEQqFWnOalHHx\n4Ivp27Wv12HUkdYuje+f8H2vwzAmsMLhMOFwuNXPb7Y7p4hkAq9V1/ib2ici/YH5wNWq+l+3rS/w\ntqqe4G6PA8aq6vUNnC/QNX6/e/ajZ/l6r69zVv+zvA7FGBMh7t053URe7TIgz23PAGYDv6tO+gCq\nuhkoFZER4kykfhXwSrSvG3Rb9mxp/qA4O7rb0fTu3NvrMIwxbdRcr56ZwFigF1AETAZCwDBAgXzg\nOlUtEpG7gNuBdRGnOFdVt4vIcOBZoBMwR1VvaeT17Iq/AeWV5Zz+5OksvGYhR3Q8wutwjDE+Y/Px\np6gqrUrIdM31u5Tu2LeDJxY+wSllp3DRRRfF/fWNMdGzkbspKlFz9I8aNarOoLGdxTt546k3GD16\ndEJe3xgTf5b4k8jm3Zt55P34DnzOyMggOzubSZMmUVBQwMP3P8ysp2YlxURxxpiWsYVYkki3w7pV\n/0kX9YLj0cjIyODr3/t6zURxlvSNSS12xZ9EuqZ35ZYRt8Q16YMzTcTMJ2aycOVCcnJyfD1XkDEm\nepb4k1R5ZXlczls9UdycZ+Yw6pRRNWUfS/7GpA7r1ZOEXv7kZd7Of7tmMfNYKdpTxNP/72luvPxG\nmyjOmCRi3TkDYH/5fkSEju07xvS8S79ayrsF7zLhmxNiel5jTHxZ4jfGmICxfvwB8kHhB6zbsa75\nA5sx7/N5VFZVxiAiY0wysMSfxPK25rFh14Y2naOiqoLnVz3Pzv07YxSVMcbvrNRjjDFJzko9pkUK\nSwspLC1s/kBjTMqxxJ8CrvjnFazZtiaq57yd/zavrLXZsY0JIiv1pIBPtn7CkF5DaN/OZuAwJois\n1BNAJx15UouTvt/W7zXGJJ4l/hShqny05aMmjyneX8wNs2/gYMXBBEVljPEjK/WkiL1le7lo5kXM\n+fEcOnXo1Ohx8Z7Z0xiTeDZy1xxi8+7N9Ozck/S0dK9DMcbEgdX4zSEeWPQAr336mtdhGGN8wq74\nU8yn2z/lxbwXueecewgXhAllhhK2Xq8xxht2xR9w/Q7vx/CjhwMwZ90cIHHr9RpjkoNd8aeo1dtW\n890Xvkv+r/PtZq4xKS7aK34b8ZNiwgVhwgVhAL7c9SX3vHsPAKHMEKHMkHeBGWN8o8nELyLPABcC\nW1V1qNuWBVwLbHMPu1NV57r77gB+DlQCt6jqm277cOBZoCMwR1V/HfPfxACHJvisUJZnsRhj/Km5\n4u8M4IJ6bQpMVdXT3J/qpH8icAVwovucx6W2xjANGK+qg4BBIlL/nMYYYxKkycSvqu8BxQ3saqiW\ndAkwU1XLVbUAWA+MEJG+QDdVXeIe9xxwaetDNi1lpR1jTENa293jZhFZKSLTRaR6Ve6jgch5fguB\nfg20b3LbTZxZ4jfGNKQ1N3enAfe6j+8DpgDjYxVQVlZWzeNQKEQoFIrVqY1JeuEw2P8SJhwOEw6H\nW/38Zrtzikgm8Fr1zd3G9onI7QCq+kd3Xy4wGfgSeEdVT3DbxwFjVfX6Bs5n3TmNaUJWlvNjTKS4\nD+Bya/bVLgPy3MevAj8SkXQRGQgMApao6hagVERGuDd7rwJsBRBjWqGsDEpKvI2hDReaxiea6845\nExgL9BKRjThX8CERGYbTuycfuA5AVVeLyD+A1UAFcEPE5fsNON05O+F058yNw+9iTEp6+21YsMB5\n/Ic/ONsXXOCUfL76CpYsgYcfdvbn5sKaNXDrrc72kiWwYQP8z/84259/Djt3wplnOtvFxc6XyVFH\ntTweKzclvyYTv6qOa6D5mSaOvx+4v4H2ZcAhpSJjTNNKS50kvmQJHHaY0xZZ6jlwAC68sHb75JOh\nX0TXicMOgy5darcLCmD9+trE/+abkJcH//u/zvZTTznbjz7qbP/nP87xEyc62wsWwKpVsfwNjRds\nygZjfKiyEtLSnMdffQVHH+08jneNv7ISysuhY0dne8sW2LMHCgudK/0tW+CJJ2DyZGd/KGRX/35g\nUzYYk+SmTHES8G23OdvVSR/in2TT0mq/cAD69HH+Pf742tfu08duMCc7S/zG+EBFBbR3/2/8+c+h\nc+eGj/PT1fWkSXDllXDCCV5HYqJl8/XWYz0WTKIdOACnngq7dzvb3bvX1vP9qPrLJxSq+9eISR6W\n+OuxxG8Spbzc+bdjR+dz162bp+G0WHXiP/dcOOIIT0MxrWSJP0JVFezd63UUJgj++le4557a7d69\nvYulrYqK4JZbnPsSJjlYjR/naischu3b4S9/ceqrItZjwcRWWRmku+vdX3kldOjgbTyx0rMnjB1b\n96aw8TfrzlnP5Ml1r8SMiYWKCqeO/+670KuX19GYVGNr7rbCvHlQ/X1TvYLAtm3OTTdj2uLgQeff\n9u1h0aLUT/qzZsHMmV5HYZoT+MR/8CA8+2zt/CfVpZ3sbJg926uoTCr4+99r++IDZGQ0fmyqOP54\nOPFEr6MwzbFSTyOqqqBd4L8WTbT274dOnZzHBw44f0H6uWumSQ1W6mmhWbOc3giNiUz6O3bEPx6T\n/FRh5EjYuNHZ7tgxuElf1ZnfZ/NmryOJr2Tt/h3YxL9unTNLYXNKSuCcc2prtcbUt2+f868ILFwI\nAwZ4G48fiMBZZ6V+P39L/Enmt79t2VDzjAxYujS4V26mabNmwU031W537epdLH5z+eWNTz2RCqoH\n4CWjQNX4w2Hn6uy7323d8ysrnb7Y1TVcExyRc9Dv2VOb4CsqnM+FXRg0Lj/f6UCRKt2kq8f9rFwJ\nr7zij5lKbXbOJrQ1Yf/lL7BrF9x9d2ziMckjMvGffz5Mnw5f/7rTTbN9oP4vil6vXjA0RVbjUK1N\n8KrOl1kyzlQaqI/siBFte/7119f28zfBEnmP56237K++aHTrVrsCWLK77jrndznvvOTOBSmf+D/7\nDGbMcJasa6vq4fZQd/i9SU3Vf9J/+aVTqqgu59hUHq33/PPQv7/TYSIZ3X039I1YdTxZPwcpn/j7\n9nXmEYmlykqnx8Ls2XU/BCa1RCb4fv2S8096vzn22OjW9/Vaeblz0XjbbU733Po9tizx+1S3bs7C\n1LGUluZM89CzZ2zPa/zL6vixcfbZXkcQnfbtoUcP5yZ+KknJ7px798LFFzu9L+IlMulXVcXvdYx3\nli93buhD8l7Z+VVZGfzud7WLz/hJVRWsXu08FnG666ZaN92UTPxdusCddybmP9a8efCLX8T/dUzi\n9egBgwc7jy3xx1aHDk6vKD92g123Du66q3bixlSUUv34VRN/p7283JnH32r9xiS3ZJ6fK6Zz9YjI\nMyJSJCJ5DeybICJVItLD3e4oIjNFZJWIrBaR2yOOHS4ieSKyTkQeieYXisYPfgDLlsXr7A3r0KE2\n6fvoO9S0QUmJ82VuEuOTT+DPf/Y2hhkzagdiBUFz328zgENujYrIAOBc4MuI5h8BqOopwHDgOhE5\nxt03DRivqoOAQSIS49utjpwcOO20eJy5ZX7yk8R/8ZjYe+steOABr6MIju7dnV5TXrr8cmdSuaBo\nttQjIpnAa6o6NKLtZeA+YBYwXFV3isj5wI3AZUB3YBEwAugEvK2qJ7jP/REQUtXrG3itqEs95eXO\nn2d+WPbts8/guOP8EYtpGy/Khiax/vxnZxT2oEFeR9J2cZ+WWUQuAQpVdVVku6q+AZQCm4ECIEdV\nS4B+QGHEoZvctph44AF47LFYna1tBg+2pJ8qLOl7Y9o0p/STCP36BXcQZlS9k0WkM3AnTpmnptnd\ndyXO1X1foAfwnojMjzagrIhRMqFQiFAz3SluvdV/N2TWr3dG+L34oiWQZPL8805X4OsP+VvUJErf\nvvHtjbdiRW05+LLL4vc68RYOhwm3YU7oqEo9IjIUeAtwZyCnP84V/AhgMrBYVf/uPm86MBdYCLwT\nUeoZB4xtS6lHFYqLne52flRZ6XzAzjjD60hMNIqKoLQ0Nf70N4c6eBC+9z146SXnvkIqiWupR1Xz\nVPUoVR2oqgNxSjinq2oRsBb4lhtEF+AsYK2qbgFKRWSEiAhwFfBKNK9b38KF8MtftuUM8ZWWZkk/\nGR11lCV9v9izx5kiIxYjZisrnX8POwzeeCP1kn5rNNedcyawGBgsIhtF5Jp6h0Renj8BpLtdP5cA\nz6jqx+7LqCWfAAALg0lEQVS+G4CngXXAelXNbUvQY8Y4ZZRkMGUK5LbptzXxtnatMxGb8Y/DDoMj\nj2z7efLy4Dvfaft5Uk1SDeBavx6OPz6BAcXAihXOTaRYfIhNfPztb869mJ/+1OtITKypOiW8Pn28\njiS+oi31JE3iLylxVs565x1/DvM2xsTHkiVQUAA//GHLjn/zTWcOoMsvj2tYvhL37pxeyciARYuS\nN+mXlsKNN9qi7X7io2se04QuXaLr6XPkkXD00fGLJxX4PvF/8EHtDZ5k7hrZtauzLkBQ+w370R13\nwN//7nUUpjknndT8OtkrV8L+/c7jYcNg5Mj4x5XMfJ34VZ0BHV984XUkbdeunfOnajJ/eaWaO+5w\npu82yePhh+Grr5zHkd3Yn3rKuZFrWsbXy0uIOEvepZr58517FkGqQfrREUd4HYGJ1hFH1A7YnD+/\ndrpsryd5Sza+vLm7ZImz0Mlxx3kdUXx89JHTT3n0aK8jCaa5c+FrX4MhQ7yOxLTWgQNwzDFOT7/D\nD/c6Gu9Fe3PXl1f8H3/sdIFM1cQ/bJjXEQTbli3Jte6rqRUO15Z4tm2DqVOdx5HrI5vm+fKKPyiq\nquChh5y5Ybp08ToaY5JLVpbzY1K4O2cqEnHqleXlXkcSDPFcg9mYZOK7K/7Jk5147E83E2vjxsE1\n18B553kdiYmFcNhyRLWUHbmb6jZuhOXL4ZJLvI4kdR044Iyj8Ns03sa0lZV6ktTu3bBhQ+12G6ba\nNo3o2NGSvjFgid83TjwRbr65dtsSf+xMmQKffup1FMb4hyV+H5o71ylLmNjo3x969/Y6CmP8w5f9\n+IOquo9yOAzvvuuUJsBudLfVFVd4HYEx/mI3d33K+ii3XX4+ZGba/Egm9dnN3RRSUeFMSmX9/Fvn\nxhud1bWMMXVZqcenQiFnZO/evc4XQIcOXkeUfGbPtqt9YxpipR5jjElyVupJUQUFcN99Xkfhf6rw\nk584A+KMMQ2zxJ8kune3aYRbQsSZ9K5fP68jMca/rNRjjDFJzko9AfDPf8Jf/+p1FP4TDteuz2yM\naVyTiV9EnhGRIhE5ZDVLEZkgIlUi0iOi7RQR+a+IfCwiq0Qk3W0fLiJ5IrJORB6J/a8RLN/4hg3o\nqq+iAh5/3JnzyBjTtCZLPSIyBtgDPKeqQyPaBwBPAUOA4aq6U0TaA8uAK1U1T0S6A7tUtUpElgA3\nqeoSEZkDPKqquQ28npV6oqRqXRaNCbqYlnpU9T2guIFdU4Hb6rWdB6xS1Tz3ucVu0u8LdFPVJe5x\nzwGXtjRA07Rf/QoWLPA6Cm/ZtYIx0Ym6xi8ilwCFqrqq3q5BgIpIrogsE5GJbns/oDDiuE1um4mB\n//t/4ayzvI7CO/v3wxlnWInHmGhENXJXRDoDdwLnRja7/3YARgNnAPuB+SKyDNgVzWtkRUxQEwqF\nCFkxu0mDB9c+rqiA9gEbi92pE/z739Ctm9eRGJM44XCYcBvmbm+2O6eIZAKvqepQERkKvAXsc3f3\nx7mCHwGEgO+o6tXu8+4CDgB/B95R1RPc9nHAWFW9voHXshp/K5WVwciR8MYb0KuX19EYYxIprt05\nVTVPVY9S1YGqOhCnhHO6qhYBbwBDRaSTe6N3LPCJqm4BSkVkhIgIcBXwSjSva5qXnu7MTROUpF9V\nBU8+aRPYGdMazXXnnAksBgaLyEYRuabeITWX56pagnPT90NgBbBMVee6u28AngbWAesb6tFj2q5P\nn9rHe/Z4F0ci7NnjTGORluZ1JMYkHxu5m4LWrYOrr4aFC62rpzFBEG2pxxJ/itq7F7p08TqK+Ejl\n382Y1rApGwxQmxgPHoRt27yNJZa2bHG6r1ZWeh2JMcnLEn+Ke/llmDLF6yhip08f+OADq+0b0xZW\n6klxqk4PGEuUxqQuK/WYOkRqk/66dU59PBkVF8OkSTY9gzGxYIk/QJ54AhYt8jqK1jvpJOulZEws\nWKnH+F44bNNQG9MUK/WYFsnNTY6eMVVVMGuW11EYk1os8QdQZaXT2ycZunl+8okzFYUxJnas1GN8\n58UXobTU6bMPcM89MHmy8zgUsrKPMfVFW+oJ2CS+pr59++Af/3CmePDKv/7l9Dy61F2ep2dPOPVU\n52ZutYjZuo0xbWSJP+D274cvvnBq6e0SVPibO9fpWnrLLc72175Wd5zB+ecnJg5jgspKPSbuFi92\n/qp4+GFn+8svnRWzTj65Zc+3Xj3GNM0maTOttmIFrFzZurJP5KLva9fChAm1N2WLi2HrVhgyJGah\nGmMiWHdO02rdukGPHrXbTa3sVlZW+3jHDmcJyOrv7K99DR5/vHZ/9+6W9I3xE0v8psbxx8P3vle7\nHZn4S0qc+wDgdAc95hjnxjA4N2Pff7/2ij89HY49NiEhG2NawRK/adBjjzk3faudcw5s2OA8Tktz\n6vSdO9fu79kzsfEZY1rPevWYOsJh52fTJnj+eadsAzB1KmRm1h532GEeBGeMiQm7uWsalZVl/eeN\nSQZ2c9cYY0yTLPGbRlnfeWNSk5V6jDEmycW01CMiz4hIkYjkNbBvgohUiUiPeu3HiMgeEZkQ0TZc\nRPJEZJ2IPNLS4IwxxsRec6WeGcAF9RtFZABwLvBlA8+ZCtSfSHcaMF5VBwGDROSQc5q6wk2NngoY\ney9q2XtRy96L1msy8avqe0BxA7umArfVbxSRS4EvgNURbX2Bbqq6xG16Dri0tQEHhX2oa9l7Ucve\ni1r2XrRe1Dd3ReQSoFBVV9Vr74rzZZBV7yn9gMKI7U1umzHGGA9ENYBLRDoDd+KUeWqa3X+zgIdU\ndZ+ILYltjDF+1WyvHhHJBF5T1aEiMhR4C3BnaaE/zhX8COBlYIDbngFUAXcD/wbeUdUT3PONA8aq\n6vUNvJZ16THGmFaI2wpcqpoHHFW9LSL5wHBV3QmcHdE+Gditqo+726UiMgJYAlwFPNrWwI0xxrRO\nc905ZwKLgcEislFErql3SEuv0G8AngbWAetVNTfqSI0xxsSErwZwGWOMiT9fTNkgIheIyFp3gNfv\nvI7HSyJSICKrRGSFiCxp/hmpo6EBgyLSQ0TmichnIvKmiGR4GWOiNPJeZIlIofvZWBGU8TAiMkBE\n3hGRT0TkYxG5xW0P3Gejifciqs+G51f8IpIGfAr8H5wbxR8C41R1jaeBeaTefZNAEZExwB7gOVUd\n6rY9AGxX1Qfci4Luqnq7l3EmQiPvRfW9s6meBpdgItIH6KOqH7ndxpfhjAW6hoB9Npp4L35IFJ8N\nP1zxfwOn7l+gquXAS8AlHsfktUDe5G5kwOD3gL+5j/9GQAb/NTF4MnCfDVXdoqofuY/3AGtwxgIF\n7rPRxHsBUXw2/JD4+wEbI7YLCfYALwXeEpGlIvILr4PxgaNUtch9XEREr7KAullEVorI9CCUNupz\nu5efBnxAwD8bEe/F+25Tiz8bfkj8dne5rlGqehrwHeBG909+A7hTtwb58zINGAgMAzYDU7wNJ7Hc\n0sa/gF+r6u7IfUH7bLjvxT9x3os9RPnZ8EPi30TtwC/cx4WNHJvyVHWz++824D84pbAgK3LrmtXz\nPm31OB7PqOpWdeF0jw7MZ0NEOuAk/edV9RW3OZCfjYj34u/V70W0nw0/JP6lODN2ZopIOnAF8KrH\nMXlCRDqLSDf3cRfgPOCQKbED5lXgZ+7jnwGvNHFsSnOTW7XLCMhnw50CZjqwWlUfjtgVuM9GY+9F\ntJ8Nz3v1AIjId4CHgTRguqr+weOQPCEiA3Gu8sEZVf1CkN4Ld8DgWKAXTs3298As4B/AMUAB8ENV\nLfEqxkRp4L2YDIRw/pRXIB+4LqLGnbJEZDSwAFhFbTnnDpyZAAL12WjkvbgTGEcUnw1fJH5jjDGJ\n44dSjzHGmASyxG+MMQFjid8YYwLGEr8xxgSMJX5jjAkYS/zGGBMwlviNMSZgLPEbY0zA/H9Xebhl\nQ0eXHgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], - "metadata": {} + "source": [ + "# Plot first 10 values of slice 0 times against first 10 of slice 0 time series\n", + "# Plot first 10 values of slice 1 times against first 10 of slice 1 time series\n", + "# Plot first 10 values of slice 0 times against first 10 of interpolated slice 1 time series\n", + "plt.plot(slice_0_times[0:10],voxel_0_timecourse[0:10],'b:+')\n", + "plt.plot(slice_1_times[0:10],voxel_1_timecourse[0:10],'g:+')\n", + "plt.plot(slice_0_times[0:10],voxel_1_timecourse_est[0:10],'kx')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make a new data matrix to contain the slice time corrected values for all voxels. We want to keep the values for z slice 0 unchanged, so make the new data matrix by copying the old one (and therefore getting the slice 0 values):" + ] + }, + { + "cell_type": "code", + "execution_count": 86, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Copy old data to a new array\n", + "new_data = good_data.copy()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Loop over every x voxel coordinate, and then loop over every y voxel coordinate.\n", + "\n", + "For each x, y voxel coordinate:\n", + "\n", + "* extract the time series at this x, y coordinate for slice 1;\n", + "* make a linear interpolator object with the slice 1 times and the extracted time series;\n", + "* resample this interpolator at the slice 0 times;\n", + "* put this new resampled time series into the new data at the same position" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# loop over all x coordinate values\n", + "for x in range(good_data.shape[0]):\n", + " # loop over all y coordinate values\n", + " for y in range(good_data.shape[1]):\n", + " # extract the time series at this x, y coordinate for slice 1;\n", + " current_timecourse = good_data[x,y,1,:]\n", + " # make a linear interpolator object with the slice 1 times and the extracted time series;\n", + " current_interp = IUS(slice_1_times, current_timecourse, k=1)\n", + " # resample this interpolator at the slice 0 times;\n", + " current_timecourse_estimate = current_interp(slice_0_times)\n", + " # put this new resampled time series into the new data at the same position\n", + " new_data[x,y,1,:] = current_timecourse_estimate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to do the same thing for all the z slices." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To do this, we want to construct an offset vector (call it `time_offset`) of length (number of z slices) such that adding the `time_offset[z]` to the acquisition time of the the first slice will give us the time of acquisition of slice `z`. The next few steps are to get to that `time_offset` vector." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, make a new vector `acquisition_order` that is length 30, where `acquisition_order[i]` is the order of acquisition of slice index `i`. For example, the first 4 elements of `acqusition_order` should be 0, 15, 1, 16." + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0., 15., 1., 16., 2., 17., 3., 18., 4., 19., 5.,\n", + " 20., 6., 21., 7., 22., 8., 23., 9., 24., 10., 25.,\n", + " 11., 26., 12., 27., 13., 28., 14., 29.])" + ] + }, + "execution_count": 88, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Make acquisition_order vector, length 30, with values:\n", + "# 0, 15, 1, 16 ... 14, 29\n", + "first_half = np.arange(0,15)\n", + "second_half = np.arange(15,30)\n", + "acquisition_order = np.array([])\n", + "for i in first_half:\n", + " acquisition_order = np.append(acquisition_order,first_half[i])\n", + " acquisition_order = np.append(acquisition_order,second_half[i])\n", + "acquisition_order" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Divide the acquisition order vector by number of slices, and multiply by the TR, to get the time offset for each z slice, relative to the start of the scan:" + ] + }, + { + "cell_type": "code", + "execution_count": 89, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0. , 1.25 , 0.08333333, 1.33333333, 0.16666667,\n", + " 1.41666667, 0.25 , 1.5 , 0.33333333, 1.58333333,\n", + " 0.41666667, 1.66666667, 0.5 , 1.75 , 0.58333333,\n", + " 1.83333333, 0.66666667, 1.91666667, 0.75 , 2. ,\n", + " 0.83333333, 2.08333333, 0.91666667, 2.16666667, 1. ,\n", + " 2.25 , 1.08333333, 2.33333333, 1.16666667, 2.41666667])" + ] + }, + "execution_count": 89, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Divide acquisition_order by number of slices, multiply by TR\n", + "time_offset = acquisition_order/acquisition_order.shape * TR\n", + "time_offset" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can do our whole slice time correction, for every slice.\n", + "\n", + "* For each z coordinate (slice index):\n", + " * Make a time vector by adding the slice time offset for this slice, to the slice_0 times. Call this the `slice_z_times` vector;\n", + " * For each x coordinate:\n", + " * For each y coordinate:\n", + " * extract the time series at this x, y, z coordinate;\n", + " * make a linear interpolator object with the `slice_z_times` and the extracted time series;\n", + " * resample this interpolator at the slice 0 times;\n", + " * put this new resampled time series into the new data at the same position" + ] + }, + { + "cell_type": "code", + "execution_count": 90, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# For each z coordinate (slice index):\n", + "for z in range(good_data.shape[2]):\n", + " # Make `slice_z_times` vector for this slice\n", + " slice_z_times = slice_0_times + time_offset[z]\n", + " # For each x coordinate:\n", + " for x in range(good_data.shape[0]):\n", + " # For each y coordinate:\n", + " for y in range(good_data.shape[1]):\n", + " # extract the time series at this x, y, z coordinate;\n", + " current_timecourse = good_data[x,y,z,:]\n", + " # make a linear interpolator object with the `slice_z_times` and the extracted time series;\n", + " current_interp = IUS(slice_z_times, current_timecourse, k=1)\n", + " # resample this interpolator at the slice 0 times;\n", + " current_timecourse_estimate = current_interp(slice_0_times)\n", + " # put this new resampled time series into the new data at the same position\n", + " new_data[x,y,z,:] = current_timecourse_estimate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Congratulations - you have just done slice timing correction on this 4D image." + ] } - ] -} \ No newline at end of file + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/day3/slice_timing_solutions.ipynb b/day3/slice_timing_solutions.ipynb index 0becd97..bbda0a1 100644 --- a/day3/slice_timing_solutions.ipynb +++ b/day3/slice_timing_solutions.ipynb @@ -646,7 +646,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.9" + "version": "2.7.10" } }, "nbformat": 4,