<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 DCGraM Algorithm
This notebook implements the D-Markov with Clustering and Graph Minimization (DCGraM) Algorithm. Its objective is to model a discrete dynamical system using a Probabilistic Finite State Machine (PFSA). 
Given a sequence X over the alphabet $\Sigma$ of length N that is an output of the original dynamical system, DCGraM works by
Step1:  Parameters
The next cell initializes the parameters that are used throughout the code. They are listed as
Step2:  Original Sequence Analysis
Make sure that the original sequence of length N is stored in the correct directory and run the cell to load it to X. After this, run the cells corresponding to the computation of the subsequence probabilities and the conditional probabilites for the value d_max, which is the last value in drange. Additional results can also be computed in the respective cells (autocorrelation and conditional entropy).
Step3:  D-Markov Machines
The next step of DCGraM consists of generating D-Markov Machines for each value of D in drange defined above. The values of p_cond for each of these values is then needed, so it is necessary to compute it above. A D-Markov Machine is a PFSA with $|\Sigma|^D$ states, each one labeled with one of the subsquences of length $D$. Given a state $\omega = \sigma_1\sigma_2\ldots\sigma_D$, for each $\sigma \in \Sigma$, it transitions to the state $\sigma_2\sigma_3\ldots\sigma_D\sigma$ with probability $\Pr(\sigma|\omega)$. This is done for all states in the D-Markov machine.
Step4:  D-Markov Machine Analysis
First of all, sequences should be generated from the D-Markov Machines. The same parameters computed in the analysis of the original sequence should be computed for the D-Markov Machines' sequences. Besides those parameters, the Kullback-Leibler Divergence and Distribution Distance between these sequences and the original sequence.
Step5:  Clustering
Now that we have obtained the D-Markov Machines, the next step of DCGraM is to cluster the states of these machines. For a given D-Markov Machine G$_D$, its states $q$ are considered points in a $\Sigma$-dimensional space, in which each dimension is labeled with a symbol $\sigma$ from the alphabet and the position of the state $q$ in this dimension is its probability of transitioning with this symbol. These point-states are then clustered together in $K$ clusters using a variation of the K-Means clustering algorithm that instead of using an Euclidean distance between points, uses the Kullback-Leibler Divergence between the point-state and the cluster centroids.
Step6:  Graph Minimization
Once that the states of the D-Markov Machines are clustered, these clusterings are then used as initial partitions of the D-Markov Machines' states. To these machines and initial partitions, a graph minimization algorithm (in the current version, only Moore) is applied in order to obtain a final reduced PFSA, the DCGraM PFSA.
Step7:  DCGraM Analysis
Now that the DCGraM machines have been generated, the same analysis done for the D-Markov Machines is used for them. Sequences are generated for each of the DCGraM machines and afterwards all of the analysis is applied to them so the comparison can be made between regular D-Markov and DCGraM.
Step8:  Plots
Once all analysis have been made, the plots of each of those parameters is created to visualize the performance. The plots have the x-axis representing the number of states of each PFSA and the y-axis represents the parameters being observed. There are always two curves
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
import yaml
import sequenceanalyzer as sa
#import dmarkov
name = 'ternary_even_shift'
tag = 'v1'
import os
if not os.path.exists(name):
    os.makedirs(name)
    os.makedirs(name + '/sequences')
    os.makedirs(name + '/pfsa')
    os.makedirs(name + '/results')
    os.makedirs(name + '/results/probabilities')
    os.makedirs(name + '/results/probabilities/conditional')
    os.makedirs(name + '/results/cond_entropies')
    os.makedirs(name + '/results/kldivergences')
    os.makedirs(name + '/results/autocorrelations')
    os.makedirs(name + '/results/prob_distances')
    os.makedirs(name + '/results/plots')
Explanation: DCGraM Algorithm
This notebook implements the D-Markov with Clustering and Graph Minimization (DCGraM) Algorithm. Its objective is to model a discrete dynamical system using a Probabilistic Finite State Machine (PFSA). 
Given a sequence X over the alphabet $\Sigma$ of length N that is an output of the original dynamical system, DCGraM works by:
Creating a D-Markov model for the original system for a given D;
Using a clustering algorithm on the D-Markov model states in order to create an initial partition;
Using a graph minimization algorithm to refine the initial partition until the final reduced PFSA is obtained.
Initialization
First, it is necessary to create the directories that store the working files for the current system. The first cell sets the system's name and the tag to be used in the current run. The following cell only has to be ran when creating modeling a new system.  A directory is then created with this tag and inside it subdirectories that contain the sequence, PFSA and result files.
End of explanation
N = 10000000
drange = range(4,11)
a = 20
Explanation: Parameters
The next cell initializes the parameters that are used throughout the code. They are listed as:
N: The original sequence length N, which is also the length of the sequences that are going to be generated by the PFSA generated by DCGraM;
drange: range of values of D for which D-Markov and DCGraM machines that will be generated;
a: value up to which the autocorrelation is computed.
End of explanation
#Open original sequence from yaml file
with open(name + '/sequences/original_len_' + str(N) + '_' + tag + '.yaml', 'r') as f:
    X = yaml.load(f)
    
#Value up to which results are computed
d_max = drange[-1]
#Initialization of variables:
p = None
p_cond = None
#Compute subsequence probabilities of occurrence up to length d_max
p, alphabet = sa.calc_probs(X, d_max)
with open(name + '/results/probabilities/original_' + tag + '.yaml', 'w') as f:
    yaml.dump(p, f)
with open(name + '/alphabet.yaml', 'w') as f:
    yaml.dump(alphabet, f)
#If p has been previously computed, use this cell to load the values
if not p:
    with open(name + '/results/probabilities/original_' + tag + '.yaml', 'r') as f:
        p = yaml.load(f)
    with open(name + '/alphabet.yaml', 'r') as f:
        alphabet = yaml.load(f)
#Compute conditional probabilities of subsequences occurring after given each symbol of the alphabet
#One of the two previous cells needs to be executed first.
if p:
    p_cond = sa.calc_cond_probs(p, alphabet, d_max - 1) 
    with open(name + '/results/probabilities/conditional/original_' + tag + '.yaml', 'w') as f:
        yaml.dump(p_cond, f)
else:
    print("Run a cell that either computes or opens the probabilities.")
#If p_cond has been previously computed, use this cell to load the values
if not p_cond:
    with open(name + '/results/probabilities/conditional/original_' + tag + '.yaml', 'r') as f:
        p_cond = yaml.load(f)
#Compute conditional entropy
if p and p_cond:
    h = sa.calc_cond_entropy(p, p_cond, d_max)
    h.to_csv(name + '/results/cond_entropies/original_' + tag + '.csv')
else:
    print("Run the conditional probabilities cell first.")
#If p_cond has been previously computed, use this cell to load the values
if not h:
    h = pd.read_csv(name + '/results/cond_entropies/original_' + tag + '.csv')
#Compute autocorrelation
aut = sa.calc_autocorr(X, a)
aut.to_csv(name + '/results/autocorrelations/original_' + tag + '.csv')
#If aut has been previously computed, use this cell to load the values
if not aut:
    aut = pd.read_csv(name + '/results/autocorrelations/original_' + tag + '.csv')
Explanation: Original Sequence Analysis
Make sure that the original sequence of length N is stored in the correct directory and run the cell to load it to X. After this, run the cells corresponding to the computation of the subsequence probabilities and the conditional probabilites for the value d_max, which is the last value in drange. Additional results can also be computed in the respective cells (autocorrelation and conditional entropy).
End of explanation
dmark_machines = []
#If the D-Markov machines have not been previously created, generate them with this cell
for D in list(map(str,drange)):
    dmark_machines.append(dmarkov.create(p_cond, D))
    dmark_machines[-1].to_csv(name + '/pfsa/dmarkov_D' + D + '_' + tag + '.csv')
#On the other hand, if there already are D-Markov machines, load them with this cell
if not dmark_machines:
    for D in drange:
        dmark_machines.append(pd.read_csv(name + '/pfsa/dmarkov_D' + D + '_' + tag + '.csv'))
Explanation: D-Markov Machines
The next step of DCGraM consists of generating D-Markov Machines for each value of D in drange defined above. The values of p_cond for each of these values is then needed, so it is necessary to compute it above. A D-Markov Machine is a PFSA with $|\Sigma|^D$ states, each one labeled with one of the subsquences of length $D$. Given a state $\omega = \sigma_1\sigma_2\ldots\sigma_D$, for each $\sigma \in \Sigma$, it transitions to the state $\sigma_2\sigma_3\ldots\sigma_D\sigma$ with probability $\Pr(\sigma|\omega)$. This is done for all states in the D-Markov machine.
End of explanation
dmark_seqs = []
#Generate sequences:
count = 0
for machine in dmark_machines:
    seq = machine.generate_sequence(N)
    with open(name + '/sequences/dmarkov_D' + str(drange[count]) + '_' + tag + '.yaml', 'w') as f:
        yaml.dump(seq, f)
    dmark_seqs.append(seq)
    count += 1
#If the sequences have been previously generated, load them here:
if not dmark_seqs:
    for D in list(map(str,drange)):
        with open(name + '/sequences/dmarkov_D' + D + '_' + tag + '.yaml', 'w') as f:
            dmark_seqs.append(yaml.load(f))
#Compute subsequence probabilities of occurrence of the D-Markov sequences
count = 0
p_dmark = []
for seq in dmark_seqs:
    p_dm, alphabet = sa.calc_probs(seq, d_max)
    p_dm.to_csv(name + '/results/probabilities/dmarkov_D'+ str(drange[count])  + '_' + tag + '.csv')
    p_dmark.append(p_dm)
    count += 1
#If p_dmark has been previously computed, use this cell to load the values
if not p_dmark:
    for D in list(map(str,drange)):
        p_dm = pd.read_csv(name + '/results/probabilities/dmarkov_D' + D + '_' + tag + '.csv')
        p_dmark.append(p_dm)
    with open(name + '/alphabet.yaml', 'r') as f:
        alphabet = yaml.load(f)
#Compute conditional probabilities of subsequences occurring after given each symbol of the alphabet
#One of the two previous cells needs to be executed first.
p_cond_dmark = []
count = 0
if p_dmark:
    for p_dm in p_dmark:
        p_cond_dm = sa.calc_cond_probs(p_dm, alphabet, d_max) 
        p_cond_dm.to_csv(name + '/results/probabilities/conditional/dmarkov_D' + str(drange[count]) + '_' + tag + '.csv')
        p_cond_dmark.append(p_cond_dm)
        count += 1
else:
    print("Run a cell that either computes or opens the probabilities.")
#If p_cond has been previously computed, use this cell to load the values
if not p_cond_dmark:
    for D in list(map(str,drange)):
        p_cond_dmark.append(pd.read_csv(name + '/results/probabilities/conditional/dmarkov_D' + D + '_' + tag + '.csv'))
#Compute conditional entropy
count = 0
h_dmark = []
if p_dmark and p_cond_dmark:
    for p_dm in p_dmark:
        h_dm = sa.calc_cond_entropy(p_dm, p_cond_dmark[count], d_max)
        h_dm.to_csv(name + '/results/cond_entropies/dmarkov_D' + str(drange[count]) + '_' + tag + '.csv')
        h_dmark.append(h_dm)
        count += 1
else:
    print("Run the conditional probabilities cell first.")
#If h_dmark has been previously computed, use this cell to load the values
if not h_dmark:
    for D in list(map(str,drange)):
        h_dmark.append(pd.read_csv(name + '/results/cond_entropies/dmarkov_D' + D + '_' + tag + '.csv'))
#Compute autocorrelation
aut_dmark = []
count = 0
for dseq in dmark_seqs:
    aut_dm = sa.calc_autocorr(dseq, a)
    aut_dm.to_csv(name + '/results/autocorrelations/dmarkov_D' + str(drange[count]) + '_'  + tag + '.csv')
    aut_dmark.append(aut_dm)
    count += 1
#If aut has been previously computed, use this cell to load the values
if not aut_dmark:
    for D in list(map(str,drange)):
        aut_dmark.append(pd.read_csv(name + '/results/autocorrelations/dmarkov_D' + D + '_' + tag + '.csv'))
#Compute the Kullback-Leibler Divergence between the sequences generated by the D-Markov Machines and the original
#sequence.
kld_dmark = []
for dseq in dmark_seqs:
    kld_dm = sa.calc_kld(dseq, X, d_max)
    kld_dmark.append(kld_dm)
    
kld_dmark.to_csv(name + '/results/kldivergences/dmarkov_' + tag + '.csv')
#If the D-Markov Kullback-Leibler divergence has been previously computed, use this cell to load the values
if not kld_dmark:
    kld_dmark = pd.read_csv(name + '/results/kldivergences/dmarkov_' + tag + '.csv')
#Compute the Probability Distances between the sequences generated by the D-Markov Machines and the original
#sequence.
pdist_dmark = []
for p_dm in p_dmark:
    pdist_dm = sa.calc_pdist(p_dm, p, d_max)
    pdist_dmark.append(pdist_dm)
    
pdist_dmark.to_csv(name + '/results/prob_distances/dmarkov_' + tag + '.csv')
#If the Probability Distances of the D-Markov Machines have been previously computed, load them with this cell.
if not pdist_dmark:
    pdist_dmark = pd.read_csv(name + '/results/prob_distances/dmarkov_' + tag + '.csv')
Explanation: D-Markov Machine Analysis
First of all, sequences should be generated from the D-Markov Machines. The same parameters computed in the analysis of the original sequence should be computed for the D-Markov Machines' sequences. Besides those parameters, the Kullback-Leibler Divergence and Distribution Distance between these sequences and the original sequence.
End of explanation
clustered = []
K = 4
for machine in dmark_machines:
    clustered.append(clustering.kmeans_kld(machine, K))
Explanation: Clustering
Now that we have obtained the D-Markov Machines, the next step of DCGraM is to cluster the states of these machines. For a given D-Markov Machine G$_D$, its states $q$ are considered points in a $\Sigma$-dimensional space, in which each dimension is labeled with a symbol $\sigma$ from the alphabet and the position of the state $q$ in this dimension is its probability of transitioning with this symbol. These point-states are then clustered together in $K$ clusters using a variation of the K-Means clustering algorithm that instead of using an Euclidean distance between points, uses the Kullback-Leibler Divergence between the point-state and the cluster centroids.
End of explanation
dcgram_machines = []
for ini_part in clustered:
    dcgram_machines.append(graphmin.moore(clustered))
Explanation: Graph Minimization
Once that the states of the D-Markov Machines are clustered, these clusterings are then used as initial partitions of the D-Markov Machines' states. To these machines and initial partitions, a graph minimization algorithm (in the current version, only Moore) is applied in order to obtain a final reduced PFSA, the DCGraM PFSA.
End of explanation
dcgram_seqs = []
#Generate sequences:
count = 0
for machine in dcgram_machines:
    seq = machine.generate_sequence(N)
    with open(name + '/sequences/dcgram_D' + str(drange[count]) + '_' + tag + '.yaml', 'w') as f:
        yaml.dump(seq, f)
    dcgram_seqs.append(seq)
    count += 1
#If the sequences have been previously generated, load them here:
if not dcgram_seqs:
    for D in list(map(str,drange)):
        with open(name + '/sequences/dcgram_D' + D + '_' + tag + '.yaml', 'w') as f:
            dcgram_seqs.append(yaml.load(f))
#Compute subsequence probabilities of occurrence of the DCGraM sequences
count = 0
p_dcgram = []
for seq in dcgram_seqs:
    p_dc, alphabet = sa.calc_probs(seq, d_max)
    p_dc.to_csv(name + '/results/probabilities/dcgram_D'+ str(drange[count])  + '_' + tag + '.csv')
    p_dcgram.append(p_dc)
    count += 1
#If p_dcgram has been previously computed, use this cell to load the values
if not p_dcgram:
    for D in list(map(str,drange)):
        p_dc = pd.read_csv(name + '/results/probabilities/dcgram_D' + D + '_' + tag + '.csv')
        p_dcgram.append(p_dm)
    with open(name + '/alphabet.yaml', 'r') as f:
        alphabet = yaml.load(f)
#Compute conditional probabilities of subsequences occurring after given each symbol of the alphabet
#One of the two previous cells needs to be executed first.
p_cond_dcgram = []
count = 0
if p_dcgram:
    for p_dc in p_dcgram:
        p_cond_dc = sa.calc_cond_probs(p_dc, alphabet, d_max) 
        p_cond_dc.to_csv(name + '/results/probabilities/conditional/dcgram_D' + str(drange[count]) + '_' + tag + '.csv')
        p_cond_dcgram.append(p_cond_dc)
        count += 1
else:
    print("Run a cell that either computes or opens the probabilities.")
#If p_cond_dcgram has been previously computed, use this cell to load the values
if not p_cond_dcgram:
    for D in list(map(str,drange)):
        p_cond_dcgram.append(pd.read_csv(name + '/results/probabilities/conditional/dcgram_D' + D + '_' + tag + '.csv'))
#Compute conditional entropy
count = 0
h_dcgram = []
if p_dcgram and p_cond_dcgram:
    for p_dc in p_dcgram:
        h_dc = sa.calc_cond_entropy(p_dc, p_cond_dcgram[count], d_max)
        h_dc.to_csv(name + '/results/cond_entropies/dcgram_D' + str(drange[count]) + '_' + tag + '.csv')
        h_dcgram.append(h_dc)
        count += 1
else:
    print("Run the conditional probabilities cell first.")
#If h_dcgram has been previously computed, use this cell to load the values
if not h_dcgram:
    for D in list(map(str,drange)):
        h_dcgram.append(pd.read_csv(name + '/results/cond_entropies/dcgram_D' + D + '_' + tag + '.csv'))
#Compute autocorrelation
aut_dcgram = []
count = 0
for dcseq in dcgram_seqs:
    aut_dc = sa.calc_autocorr(dcseq, a)
    aut_dc.to_csv(name + '/results/autocorrelations/dcgram_D' + str(drange[count]) + '_'  + tag + '.csv')
    aut_dcgram.append(aut_dc)
    count += 1
#If aut has been previously computed, use this cell to load the values
if not aut_dcgram:
    for D in list(map(str,drange)):
        aut_dmark.append(pd.read_csv(name + '/results/autocorrelations/dcgram_D' + D + '_' + tag + '.csv'))
#Compute the Kullback-Leibler Divergence between the sequences generated by the DCGraM Machines and the original
#sequence.
kld_dcgram = []
for dcseq in dcgram_seqs:
    kld_dc = sa.calc_kld(dcseq, X, d_max)
    kld_dcgram.append(kld_dc)
    
kld_dcgram.to_csv(name + '/results/kldivergences/dcgram_' + tag + '.csv')
#If the DCGraM Kullback-Leibler divergence has been previously computed, use this cell to load the values
if not kld_dcgram:
    kld_dcgram = pd.read_csv(name + '/results/kldivergences/dcgram_' + tag + '.csv')
#Compute the Probability Distances between the sequences generated by the DCGraM Machines and the original
#sequence.
pdist_dcgram = []
for p_dc in p_dcgram:
    pdist_dc = sa.calc_pdist(p_dc, p, d_max)
    pdist_dcgram.append(pdist_dc)
    
pdist_dcgram.to_csv(name + '/results/prob_distances/dcgram_' + tag + '.csv')
#If the Probability Distances of the DCGraM Machines have been previously computed, load them with this cell.
if not pdist_dcgram:
    pdist_dcgram = pd.read_csv(name + '/results/prob_distances/dcgram_' + tag + '.csv')
Explanation: DCGraM Analysis
Now that the DCGraM machines have been generated, the same analysis done for the D-Markov Machines is used for them. Sequences are generated for each of the DCGraM machines and afterwards all of the analysis is applied to them so the comparison can be made between regular D-Markov and DCGraM.
End of explanation
#initialization
import matplotlib.pyplot as plt
#Labels to be used in the plots' legends
labels = ['D-Markov Machines, D from ' + str(drange[0]) + ' to ' + str(d_max),
          'DCGraM Machines, D from ' + str(drange[0]) + ' to ' + str(d_max),
          'Original Sequence Baseline']
#Obtaining number of states of the machines to be used in the x-axis:
states_dmarkov = []
for dm in dmark_machines:
    states_dmarkov.append(dm.shape[0])
    
states_dcgram = []
for dc in dcgram_machines:
    states_dcgram.append(dc.shape[0])
    
states = [states_dmarkov, states_dcgram]
#Conditional Entropy plots
H = 10
h_dmark_curve = []
for h_dm in h_dmarkov:
    h_dmark_curve.append(h_dm[H])
plt.semilogx(states[0], h_dmark_curve, marker='o', label=labels[0])
    
h_dcgram_curve = []
for h_dc in h_dcgram:
    h_dcgram_curve.append(h_dc[H])
plt.semilogx(states[1], h_dcgram_curve, marker='x', label=labels[1])
    
#Opening original sequence baseline:
h_base = h[H]
plt.axhline(y=h_base, color='k', linewidth = 3, label=labels[2])
plt.xlabel('Number of States', fontsize=16)
plt.yalbel('$h_' + str(H) + '$', fontsize=16)
plt.legend(loc='upper right', shadow=False, fontsize='large')
plt.title('Conditional Entropy',fontsize=18,weight='bold')
plt.savefig(name + '/plots/conditional_entropy_' + tag + '.eps' , bbox_inches='tight', format='eps',dpi=1000)
plt.show()
#Kullback-Leibler plots
plt.semilogx(states[0], kld_dmark, marker='o', label=labels[0])
plt.semilogx(states[1], kld_dcgram, marker='x', label=labels[1])
plt.xlabel('Number of States', fontsize=16)
plt.yalbel('$k_' + str(H) + '$', fontsize=16)
plt.legend(loc='upper right', shadow=False, fontsize='large')
plt.title('Kullback-Leibler Divergence',fontsize=18,weight='bold')
plt.savefig(name + '/plots/kldivergence_' + tag + '.eps' , bbox_inches='tight', format='eps',dpi=1000)
plt.show()
#Probability Distance plots
plt.semilogx(states[0], pdist_dmark, marker='o', label=labels[0])
plt.semilogx(states[1], pdist_dcgram, marker='x', label=labels[1])
plt.xlabel('Number of States', fontsize=16)
plt.yalbel('$P_' + str(H) + '$', fontsize=16)
plt.legend(loc='upper right', shadow=False, fontsize='large')
plt.title('Probability Distance',fontsize=18,weight='bold')
plt.savefig(name + '/plots/prob_distance_' + tag + '.eps' , bbox_inches='tight', format='eps',dpi=1000)
plt.show()
#TODO: Think how to have good plots for autocorrelation
Explanation: Plots
Once all analysis have been made, the plots of each of those parameters is created to visualize the performance. The plots have the x-axis representing the number of states of each PFSA and the y-axis represents the parameters being observed. There are always two curves: one for the DCGraM machines and one for the D-Markov Machines. Each point in these curves represents a machine of that type for a certain value of $D$. The further right a point is in the curve, the higher its $D$-value. On the curve for conditional entropy there is also a black representing the original sequence's conditional entropy for the $L$ being used as a baseline.
End of explanation
<END_TASK>