3.区間分割アルゴリズム

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv("E003_FractionalMethylatioLength.txt")

df.head()

R = (df.unconverted - df.converted).values 

print(f"データ数: {len(R)}")

from kerasy.Bio.maxsets import MSS

model = MSS()
model.run(R, limit=5, verbose=1)
np.save(f"{limit:>02}.npy", model.cotinuous_area)

def decorate_axes(ax, stats, xlabel):
    ax.set_xlabel(xlabel, fontsize=16)
    ax.set_title(stats, fontsize=18)
    ax.set_ylabel(stats, fontsize=16)
    return ax

fig, ((ax1,ax2),(ax3,ax4),(ax5,ax6)) = plt.subplots(nrows=3,ncols=2,figsize=(15,18), sharex="all")

for limit in range(5,11):
    detected_area = np.load(f"{limit:>02}.npy")
    detected_length = detected_area[:,1]
    
    median, tile90 = np.percentile(a=detected_length, q=[50,90])
    mean = np.mean(detected_length)
    var = np.var(detected_length)
    score = np.sum(detected_length)
    num = len(detected_length)
    
    ax1.scatter(limit, median, color="red",    s=100)
    ax2.scatter(limit, tile90, color="blue",   s=100)
    ax3.scatter(limit, mean,   color="green",  s=100)
    ax4.scatter(limit, var,    color="orange", s=100)
    ax5.scatter(limit, score,  color="black",  s=100)
    ax6.scatter(limit, num,    color="black",  s=100)

ax1 = decorate_axes(ax1, stats="median", xlabel="limit")
ax2 = decorate_axes(ax2, stats="90%-tile", xlabel="limit")
ax3 = decorate_axes(ax3, stats="mean", xlabel="limit")
ax4 = decorate_axes(ax4, stats="variance", xlabel="limit")
ax5 = decorate_axes(ax5, stats="score", xlabel="limit")
ax6 = decorate_axes(ax6, stats="num detected area", xlabel="limit")

plt.suptitle("The relationship between\nlimit length and Statistics.", fontsize=20)
# plt.tight_layout()
plt.show()

plt.hist(detected_length, bins=30)
plt.title(f"The histgram of detected length (limit={limit})")
plt.xlabel("detected length"), plt.ylabel("frequency")
plt.show()

from kerasy.ML.HMM import BinomialHMM
from kerasy.clib._pyutils import extract_continuous_area

X = df[["converted", "unconverted"]].values

model = BinomialHMM(n_hstates=2)

model.fit(X)

# Not necessary
model.save_params("kadai03_HMM.json")

# Not necessary
model.load_params("kadai03_HMM.json")

predictions = model.predict(X).astype(np.int32)
detected_area = np.asarray(extract_continuous_area(predictions), dtype=np.int32)
detected_length = detected_area[:,1]

plt.hist(detected_length, bins=30)
plt.title(f"The histgram of detected length (HMM)")
plt.xlabel("detected length"), plt.ylabel("frequency")
plt.show()

median, tile90 = np.percentile(a=detected_length, q=[50,90])
print(f"- median  : {median}")
print(f"- 90%-tile: {tile90}")
print(f"- mean    : {np.mean(detected_length)}")
print(f"- variance: {np.var(detected_length)}")
print(f"- score   : {np.sum(detected_length)}")
print(f"- num     : {len(detected_length)}")

raise NotImplementedError("Not Checked.")

# 内訳を整形して記述する関数
def breakdown(combs):
    """ display breakdowns """
    use_coins = sorted(set(combs))
    num_coins = [combs.count(coin) for coin in use_coins]
    total_pay = [n*coin for n,coin in zip(use_coins,num_coins)]
    width_coin  = max([len(str(e)) for e in use_coins]+[len("coins")])
    width_num   = max([len(str(e)) for e in num_coins]+[len("number")])
    width_total = max([len(str(e)) for e in total_pay]+[len("pay"),len(str(sum(total_pay)))])
    width_line  = width_coin+width_num+width_total+2
    
    print_func = lambda c,n,p: print(f"{c:^{width_coin}}|{n:>{width_num}}|{p:>{width_total}}")
    print_func('coins','number','pay')
    print("="*width_line)
    for coin,num,t in zip(use_coins,num_coins,total_pay):
        print_func(coin,num,t)
    print("-"*width_line)
    print_func('total',sum(num_coins),sum(total_pay))

def not_restricted_smart_pay(coins, total):
    """
    Find the minimum number of coin combinations by using Dynamic Programming.
    @params coins: (list) Coin Types.
    @params total: (int) Amount of Payment.
    """
    total += 1 # because 0-origin.
    coins = np.asarray(coins, dtype=int)
    if len(np.unique(coins)) < len(coins):
        raise ValueError("All elements of `coins` must be different integers.")

    # Initialization.
    B = np.zeros(shape=(total), dtype=int) # Memory for Traceback.
    m = np.full(fill_value=np.inf, shape=(total))
    m[0] = 0
    for coin in coins: m[coin] = 1
        
    # Recursion
    for t in range(1,total):
        cands = [m[t-coin] if (t-coin)>=0 else np.inf for coin in coins]
        if np.all(np.array(cands)==np.inf): continue
        k = np.argmin(cands)
        m[t] = m[t-coins[k]]+1
        B[t] = t-coins[k]
    # Trace Back.
    idx = total-1
    combs = []
    while idx:
        last = B[idx]
        combs.append(idx-last)
        idx = last
    breakdown(combs)

not_restricted_smart_pay([11,7,4,2], 100)

def restricted_smart_pay(coins, total, L):
    # Initialization.
    coins = np.asarray(coins, dtype=int)
    if len(np.unique(coins)) < len(coins):
        raise ValueError("All elements of `coins` must be different integers.")
    K = len(coins)
    total += 1 # because 0-origin.

    B = np.zeros(shape=(total*2), dtype=int) # Memory for Traceback.
    m = np.full(fill_value=np.inf, shape=(total,2))
    m[0,0] = 0
    for k in range(1,K): m[coins[k],0] = 1
    
    # Recursion
    for t in range(1,total):
        # m0
        cands = [m[t-coins[k], 0] if (t-coins[k])>=0 else np.inf for k in range(1,K)]
        if np.all(np.array(cands)==np.inf): continue
        k = np.argmin(cands)+1
        m[t,0] = m[t-coins[k], 0] + 1
        B[t] = t-coins[k]

        # mL
        cands = [m[t-coins[k], 1] if (t-coins[k])>=0 else np.inf for k in range(1,K)]
        k = np.argmin(cands)+1
        if t-L*coins[0]>=0 and m[t-L*coins[0],0]+L<m[t-coins[k],1]+1:
            m[t,1] = m[t-L*coins[0],0]+L
            B[t+total] = t-L*coins[0]
        else:
            if np.all(np.array(cands)==np.inf): continue
            m[t,1] = m[t-coins[k],1]+1
            B[t+total] = t-coins[k]+total
    
    num_d1 = np.argmin([m[-1-coins[0]*l,0]+l for l in range(L+1)])
    idx = total-1-num_d1*coins[0] if num_d1!=L else total*2-1
    combs = [coins[0] for _ in range(num_d1)]
    while idx:
        last = B[idx]
        combs.append(idx%total-last%total)
        idx = last
    if num_d1==L: combs.pop(-1)
    breakdown(combs)

restricted_smart_pay([11,7,4,2], 100, 1)

def smart_pay(coins, total, limit=None, verbose=1):
    """
    Find the minimum number of coin combinations by using Dynamic Programming.
    @params coins: (int list) Coins.
    @params total: (int) Amount of Payment.
    @params limit: (int) Maximum number of times a restricted coin can be used.
    """
    total += 1 # because 0-origin.
    if len(set(coins)) < len(coins):
        raise ValueError("All elements of `coins` must be different integers.")
    restricted = coins[0]
    free_coins = coins[1:]
    
    if limit is None:
        limit = total//restricted+1
    elif verbose:
        print(f'{restricted} coin can only be used up to {limit} times at the same time.')

    # Initialization.
    B = [0 for _ in range(total)] # Memory for Traceback.
    m = [0 if t==0 else 1 if t in free_coins else float('inf') for t in range(total)]
    
    # Recursion
    for t in range(1,total):
        cands = [m[t-coin] if (t-coin)>=0 else float('inf') for coin in free_coins]
        if sum([e!=float('inf') for e in cands])==0:
            continue # B[t]=0; m[t]=float('inf') : default.            
        minnum = min(cands)
        m[t],B[t] = [(e+1,t-coin) for e,coin in zip(cands,free_coins) if e==minnum][0]

    ms = [(l,m[-1-restricted*l]+l) for l in range(limit+1) if restricted*l<=total]
    num_restricted, num_total = min(ms, key=lambda x:x[1])    
    idx = total-1-restricted*num_restricted
    combs = [restricted for _ in range(num_restricted)]
    while idx:
        last = B[idx]
        combs.append(idx-last)
        idx = last
    if verbose: breakdown(combs)

for l in range(0,9,2):
    smart_pay([11,7,4,2], 100, limit=l)
    print()

from kerasy.utils import measure_complexity

seed   = 0
coin_cands = np.arange(1,100)

total = 100000
K = 10
coins = np.random.RandomState(seed).choice(coin_cands, K, replace=False)

for l in [10,50,100]:
    kwargs = {'limit':l, 'verbose':0}
    time = measure_complexity(smart_pay, coins, total, **kwargs)
    print(f"L:{l:>3}, processing time: {time:>5.3f}")

kwargs = {'limit':10, 'verbose':0}
Tcands = [1000,10000,100000,1000000]
Kcands = [3,30,300]
coin_cands = np.arange(1,1000)

print_func = lambda K,T,time: print(f"{K:>5}|{T:>7}|{time:>7.3f}[s]")

memory = []
print(f"{'coins':^5}|{'total':^7}|{'time[s]':>10}")
print("="*24)
for K in Kcands:
    coins = np.random.RandomState(seed).choice(coin_cands, K, replace=False)
    for T in Tcands:      
        time = measure_complexity(smart_pay, coins, T, repetitions_=3, **kwargs)
        print_func(K,T,time)
        memory.append([K,T,time])

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
for ax,k in zip(axes,Kcands):
    tmp = [(e[1],e[2]) for e in memory if e[0]==k]
    for t,time in tmp:
        ax.scatter(t,time, color="red")
    ax.grid()
    ax.set_xscale("log"), ax.set_yscale("log")
    ax.set_title(f"$K={k}$"), ax.set_xlabel("total amount"), ax.set_ylabel("Processing time")    
plt.show()

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(24,6))
for ax,T in zip(axes,Tcands):
    tmp = [(e[0],e[2]) for e in memory if e[1]==t]
    for k,time in tmp:
        ax.scatter(k,time, color="blue")
    ax.grid()
    ax.set_xscale("log"), ax.set_yscale("log")
    ax.set_title(f"$T={t}$"), ax.set_xlabel("total amount"), ax.set_ylabel("Processing time")    
plt.show()