Skip to content

Instantly share code, notes, and snippets.

@hirokai

hirokai/201703110 stability of FDH-BOD on hydrogel.py

Created March 10, 2017 02:46
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save hirokai/15a94f4420c72af505117277bf45d39c to your computer and use it in GitHub Desktop.
Save hirokai/15a94f4420c72af505117277bf45d39c to your computer and use it in GitHub Desktop.
Download ZIP
Stability of FDH-BOD on hydrogel data on 2016/3/9
Raw
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment