Stability of FDH-BOD on hydrogel data on 2016/3/9

201703110 stability of FDH-BOD on hydrogel.py

# -*- coding: utf-8 -*-

import csv
import numpy as np
import matplotlib.pyplot as plt
import math
import os
import operator


#
# File manipulation
#

def ambiguous_path(path_head):
    import glob
    files = glob.glob(path_head)
    assert len(files) == 1, 'Only one file must be found (%d files found).' % len(files)
    return files[0]


colors10 = """#1f77b4
#ff7f0e
#2ca02c
#d62728
#9467bd
#8c564b
#e377c2
#7f7f7f
#bcbd22
#17becf""".split('\n')

os.chdir(os.path.dirname(__file__))


def calc_time_cycle(ts, ys, count=620):
    # count = 0
    # last_t = 0
    # # print(ys.shape[0])
    # for i in range(ys.shape[0]):
    #     if i >= 2 and ys[i] - ys[i - 2] > 400:
    #         # print(ts[i]-last_t)
    #         last_t = ts[i]
    #         if count == 0:
    #             t0 = ts[i]
    #         count += 1
    # t1 = last_t
    # # return (t1 - t0) / count
    return (4.966633E+4 - 2.900019E+1) / count


def scanl(f, base, l):
    for x in l:
        base = f(base, x)
        yield base


def prepare_cycles(num_samples, time_start, time_cycle, num_cycles):
    time_interval = [time_cycle - ((num_samples - 1) * 10) if i == 0 else 10 for i in range(num_samples)]
    time_interval_accum = [0] + list(scanl(operator.add, 0, time_interval))[0:-1]
    print(time_interval_accum)
    # [[]] * num won't work. If you do so, elements would be just references to the same list.
    _rs = [[] for _ in range(num_samples)]
    _ts2 = [[] for _ in range(num_samples)]
    time_points = np.tile(time_interval_accum, num_cycles)
    for i in range(num_cycles):
        i0 = i * num_samples
        i1 = i0 + len(time_interval)
        time_points[i0:i1] += (time_start + i * time_cycle)
    return time_interval, _rs, _ts2, time_points


def read_data():
    with open(ambiguous_path("../02 *.DY20")) as f:
        reader = csv.reader(f, delimiter='\t')
        for _ in range(21):
            reader.next()
        vs = []
        for row in reader:
            if len(row) < 2:
                break
            vs.append(map(float, row))

    vs = np.array(vs)
    ts = vs[:, 0]
    ys = vs[:, 1] * 1000
    return ts, ys


def split_data(ts, ys, num_cycles, time_start):
    num_samples = 7

    time_cycle = calc_time_cycle(ts, ys, num_cycles)
    time_interval, rs, ts2, time_points = prepare_cycles(num_samples, time_start, time_cycle, num_cycles * 2)
    assert len(time_interval) == num_samples
    assert np.abs(np.sum(time_interval) - time_cycle) <= 0.01

    c = 0
    for i in range(len(ts)):
        timepoint = time_points[c]
        if ts[i] >= timepoint:
            ts2[c % num_samples].append(ts[i])
            rs[c % num_samples].append(ys[i])
            c += 1

    l = min([len(r) for r in rs])
    rs = np.array([r[0:l] for r in rs])
    ts = np.array([t[0:l] for t in ts2])

    time_offset_min = np.array([-45, -30, -23, 9, 22, 37, 50])
    offset_matrix = np.tile(time_offset_min * 60, (ts.shape[1], 1)).transpose()
    ts -= offset_matrix
    np.savetxt('time.tsv', ts.transpose(), delimiter='\t', fmt='%.2f')
    np.savetxt('voltage.tsv', rs.transpose(), delimiter='\t', fmt='%.2f')
    return ts, rs


def do_job(num_cycles, time_start, show=False):
    ts_raw, ys_raw = read_data()
    ts, ys = split_data(ts_raw, ys_raw, num_cycles, time_start)

    # print('Voltage max [mV]: ')
    vmax = np.amax(ys, axis=1)
    # print(vmax)
    # print('Voltage after 12 h [mV]: ')
    v12h = ys[:, -1]
    # print(v12h)
    print('Voltage max and after 12 h [mV]')
    print(np.array([vmax, v12h]).transpose())

    plt.subplot(3,1,1)
    for i, v in enumerate(zip(ts, ys)):
        t, r = v
        # if i in [0,6]:
        t_in_min = np.array(t) / 60
        plt.plot(t_in_min, r, color=colors10[i % 10], lw=1, label=str(i + 1))

    # http://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
    # ax = plt.axes()
    # box = ax.get_position()
    # ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
    plt.legend(loc='center left', ncol=1, bbox_to_anchor=(1, 0.5))

    plt.xlabel('Time [min]')
    plt.ylabel('Voltage [mV]')
    plt.xlim([0, 60 * 12])
    plt.ylim([0, 800])
    # plt.savefig('20170309 voltage over 12 h.pdf')
    plt.title('Voltage vs time')

    plt.subplot(3,1,2)
    resistance_list = [100, 100, 100, 1, 1, 1, 10]
    for i, v in enumerate(zip(ts, ys)):
        t, voltage = v
        resistance = resistance_list[i]
        current = voltage / resistance / (0.05 * 1)
        # if i in [0,6]:
        t_in_min = np.array(t) / 60
        plt.plot(t_in_min, current, color=colors10[i % 10], lw=1, label=str(i + 1))

    # http://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
    # ax = plt.axes()
    # box = ax.get_position()
    # ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
    plt.legend(loc='center left', ncol=1, bbox_to_anchor=(1, 0.5))

    plt.xlabel('Time [min]')
    plt.ylabel('Current density [uA/cm^2]')
    plt.xlim([0, 60 * 12])
    plt.ylim([0, 5000])
    # plt.savefig('20170309 voltage over 12 h.pdf')
    plt.title('Current vs time')

    plt.subplot(3,1,3)
    resistance_list = [100, 100, 100, 1, 1, 1, 10]
    for i, v in enumerate(zip(ts, ys)):
        t, voltage = v
        resistance = resistance_list[i]
        current = voltage / resistance / (0.05 * 1)
        # if i in [0,6]:
        t_in_min = np.array(t) / 60
        plt.plot(t_in_min, current, color=colors10[i % 10], lw=1, label=str(i + 1))

    # http://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
    # ax = plt.axes()
    # box = ax.get_position()
    # ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
    plt.legend(loc='center left', ncol=1, bbox_to_anchor=(1, 0.5))

    plt.xlabel('Time [min]')
    plt.ylabel('Current density [uA/cm^2]')
    plt.xlim([0, 60])
    plt.ylim([0, 5000])
    # plt.savefig('20170309 voltage over 12 h.pdf')
    plt.title('Current vs time')
    if show:
        plt.show()


do_job(num_cycles=620, time_start=2.900019E+1 + 5, show=True)