## Equations
s_{pi, tau=1}=sigma((1/2)(lnD+ln(B^dagger_{pi, tau}s_{pi, tau+1}))+lnA^T*o_tau)
s_{pi, tau>1}=sigma((1/2)(ln(B_{pi, tau-1}s_{pi, tau-1})+ln(B^dagger_{pi, tau}s_{pi, tau+1}))+lnA^T*o_tau)
G_pi=sum_tau(A*s_{pi, tau}*(ln(A*s_{pi, tau})-lnC_tau)-diag(A^TlnA)*s_{pi, tau})
pi=sigma(-G)