import sys
sys.path.append("/home/gabor/github/butools/Python")
import math
import numpy as np
import numpy.matlib as ml
import matplotlib.pyplot as plt
import butools
from butools.utils import *
from butools.ph import *
from butools.dph import *
from butools.map import *
from butools.moments import *
from butools.reptrans import*
from butools.mc import *
from butools.dmap import *
from butools.trace import *
from butools.mam import *
from butools.queues import *
from butools.fitting import *
from contextlib import redirect_stdout
import os


print('---BuTools: DMAP package test file---')
print('Enable the verbose messages with the BuToolsVerbose flag')
butools.verbose = True
print('Enable input parameter checking with the BuToolsCheckInput flag')
butools.checkInput = True
np.set_printoptions(precision=5,linewidth=1024)
print('========================================')
print('Testing BuTools function MarginalDistributionFromDMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('a, A = MarginalDistributionFromDMAP(D0, D1):')
a, A = MarginalDistributionFromDMAP(D0, D1)
print('a = ')
print(a)
print('A = ')
print(A)
assert Length(a)==D0.shape[0] and CheckDPHRepresentation(a, A), "MarginalDistributionFromDMAP returned a wrong DPH representation!"
print('========================================')
print('Testing BuTools function MarginalMomentsFromDMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('moms = MarginalMomentsFromDMAP(D0, D1):')
moms = MarginalMomentsFromDMAP(D0, D1)
print('moms = ')
print(moms)
assert Length(moms)==2*D0.shape[0]-1 and CheckMoments(moms), "MarginalMomentsFromDMAP returned wrong moments!"
print('========================================')
print('Testing BuTools function MarginalDistributionFromDRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1, -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('a, A = MarginalDistributionFromDRAP(H0, H1):')
a, A = MarginalDistributionFromDRAP(H0, H1)
print('a = ')
print(a)
print('A = ')
print(A)
assert Length(a)==H0.shape[0] and CheckMGRepresentation(a, A), "MarginalDistributionFromDRAP returned a wrong MG representation!"
print('========================================')
print('Testing BuTools function MarginalMomentsFromDRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1, -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('moms = MarginalMomentsFromDRAP(H0, H1):')
moms = MarginalMomentsFromDRAP(H0, H1)
print('moms = ')
print(moms)
assert Length(moms)==2*H0.shape[0]-1 and CheckMoments(moms), "MarginalMomentsFromDRAP returned wrong moments!"
print('========================================')
print('Testing BuTools function MarginalDistributionFromDMMAP')
D0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('D1 = ')
print(D1)
D2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('D2 = ')
print(D2)
D3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('D3 = ')
print(D3)
print('Test:')
print('-----')
print('a, A = MarginalDistributionFromDMMAP([D0, D1, D2, D3]):')
a, A = MarginalDistributionFromDMMAP([D0, D1, D2, D3])
print('a = ')
print(a)
print('A = ')
print(A)
assert Length(a)==D0.shape[0] and CheckDPHRepresentation(a, A), "MarginalDistributionFromDMMAP returned a wrong DPH representation!"
print('========================================')
print('Testing BuTools function MarginalMomentsFromDMMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('D1 = ')
print(D1)
D2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('D2 = ')
print(D2)
D3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('D3 = ')
print(D3)
print('Test:')
print('-----')
print('moms = MarginalMomentsFromDMMAP([D0, D1, D2, D3]):')
moms = MarginalMomentsFromDMMAP([D0, D1, D2, D3])
print('moms = ')
print(moms)
assert Length(moms)==2*D0.shape[0]-1 and CheckMoments(moms), "MarginalMomentsFromDMMAP returned wrong moments!"
print('========================================')
print('Testing BuTools function MarginalDistributionFromDMRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0.15, 0.2, 0.18],[-0.23, 0.17, 0.22],[0.19, 0.15, 0.16]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0.01, 0.08, 0.16],[0.02, 0.2, 0.07],[0.02, 0.15, 0.17]])
print('H1 = ')
print(H1)
H2 = ml.matrix([[0.14, 0.07, 0.01],[0.19, 0.02, 0.34],[0.06, 0.1, 0]])
print('H2 = ')
print(H2)
print('Test:')
print('-----')
print('a, A = MarginalDistributionFromDMRAP([H0, H1, H2]):')
a, A = MarginalDistributionFromDMRAP([H0, H1, H2])
print('a = ')
print(a)
print('A = ')
print(A)
assert Length(a)==H0.shape[0] and CheckMGRepresentation(a, A), "MarginalDistributionFromDMRAP returned a wrong MG representation!"
print('========================================')
print('Testing BuTools function MarginalMomentsFromDMRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0.15, 0.2, 0.18],[-0.23, 0.17, 0.22],[0.19, 0.15, 0.16]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0.01, 0.08, 0.16],[0.02, 0.2, 0.07],[0.02, 0.15, 0.17]])
print('H1 = ')
print(H1)
H2 = ml.matrix([[0.14, 0.07, 0.01],[0.19, 0.02, 0.34],[0.06, 0.1, 0]])
print('H2 = ')
print(H2)
print('Test:')
print('-----')
print('moms = MarginalMomentsFromDMRAP([H0, H1, H2]):')
moms = MarginalMomentsFromDMRAP([H0, H1, H2])
print('moms = ')
print(moms)
assert Length(moms)==2*H0.shape[0]-1 and CheckMoments(moms), "MarginalMomentsFromDMRAP returned wrong moments!"
print('========================================')
print('Testing BuTools function LagCorrelationsFromDMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('corr = LagCorrelationsFromDMAP(D0, D1, 3):')
corr = LagCorrelationsFromDMAP(D0, D1, 3)
print('corr = ')
print(corr)
assert Length(corr)==3 and np.all(corr<=1) and np.all(corr>=-1), "LagCorrelationsFromDMAP returned wrong autocorrelation coefficients!"
print('========================================')
print('Testing BuTools function LagCorrelationsFromDRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1, -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('corr = LagCorrelationsFromDRAP(H0, H1, 3):')
corr = LagCorrelationsFromDRAP(H0, H1, 3)
print('corr = ')
print(corr)
assert Length(corr)==3 and np.all(corr<=1) and np.all(corr>=-1), "LagCorrelationsFromDRAP returned wrong autocorrelation coefficients!"
print('========================================')
print('Testing BuTools function LagkJointMomentsFromDMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('Nm = LagkJointMomentsFromDMAP(D0, D1, 4, 1):')
Nm = LagkJointMomentsFromDMAP(D0, D1, 4, 1)
print('Nm = ')
print(Nm)
print('moms = MarginalMomentsFromDMAP(D0, D1, 4):')
moms = MarginalMomentsFromDMAP(D0, D1, 4)
print('moms = ')
print(moms)
cjm = np.zeros(3)
for i in range(1,4,1):
    Nx = LagkJointMomentsFromDMAP(D0, D1, 1, i)
    cjm[i-1] = (Nx[1, 1]-moms[0]**2)/(moms[1]-moms[0]**2)
print("Correlation from joint moments:")
print(cjm)
print('corr = LagCorrelationsFromDMAP(D0, D1, 3):')
corr = LagCorrelationsFromDMAP(D0, D1, 3)
print('corr = ')
print(corr)
mNm1 = Nm[0, 1:].A.flatten()
mNm2 = Nm[1:, 0].A.flatten()
assert np.all(Nm>0) and la.norm(np.array(moms)-mNm1)<10**-12 and la.norm(np.array(moms)-mNm2)<10**-12 and la.norm(corr-cjm)<10**-12, "Joint moment matrix is invalid!"
print('========================================')
print('Testing BuTools function LagkJointMomentsFromDRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1, -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('Nm = LagkJointMomentsFromDRAP(H0, H1, 4, 1):')
Nm = LagkJointMomentsFromDRAP(H0, H1, 4, 1)
print(Length(Nm))
print('moms = MarginalMomentsFromDRAP(H0, H1, 4):')
moms = MarginalMomentsFromDRAP(H0, H1, 4)
print('moms = ')
print(moms)
cjm = np.zeros(3)
for i in range(1,4,1):
    Nx = LagkJointMomentsFromDRAP(H0, H1, 1, i)
    cjm[i-1] = (Nx[1, 1]-moms[0]**2)/(moms[1]-moms[0]**2)
print("Correlation from joint moments:")
print(cjm)
print('corr = LagCorrelationsFromDRAP(H0, H1, 3):')
corr = LagCorrelationsFromDRAP(H0, H1, 3)
print('corr = ')
print(corr)
mNm1 = Nm[0, 1:].A.flatten()
mNm2 = Nm[1:, 0].A.flatten()
assert np.all(Nm>0) and la.norm(np.array(moms)-mNm1)<10**-12 and la.norm(np.array(moms)-mNm2)<10**-12 and la.norm(corr-cjm)<10**-12, "Joint moment matrix is invalid!"
print('========================================')
print('Testing BuTools function LagkJointMomentsFromDMMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('D1 = ')
print(D1)
D2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('D2 = ')
print(D2)
D3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('D3 = ')
print(D3)
print('Test:')
print('-----')
print('Nm = LagkJointMomentsFromDMMAP([D0, D1, D2, D3], 3, 1):')
Nm = LagkJointMomentsFromDMMAP([D0, D1, D2, D3], 3, 1)
print(Nm[0])
print(Nm[1])
print(Nm[2])
assert Length(Nm)==3 and la.norm(Nm[0]+Nm[1]+Nm[2]-LagkJointMomentsFromDMAP(D0, D1+D2+D3, 3, 1))<10**-12, "Joint moment matrix is invalid!"
print('========================================')
print('Testing BuTools function LagkJointMomentsFromDMRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0.15, 0.2, 0.18],[-0.23, 0.17, 0.22],[0.19, 0.15, 0.16]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0.01, 0.08, 0.16],[0.02, 0.2, 0.07],[0.02, 0.15, 0.17]])
print('H1 = ')
print(H1)
H2 = ml.matrix([[0.14, 0.07, 0.01],[0.19, 0.02, 0.34],[0.06, 0.1, 0]])
print('H2 = ')
print(H2)
print('Test:')
print('-----')
print('Nm = LagkJointMomentsFromDMRAP([H0, H1, H2], 3, 2):')
Nm = LagkJointMomentsFromDMRAP([H0, H1, H2], 3, 2)
print(Nm[0])
print(Nm[1])
assert Length(Nm)==2 and la.norm(Nm[0]+Nm[1]-LagkJointMomentsFromDRAP(H0, H1+H2, 3, 2))<10**-12, "Joint moment matrix is invalid!"
print('========================================')
print('Testing BuTools function RandomDMAP')
print('Test:')
print('-----')
print('D0, D1 = RandomDMAP(4, 5.62, 10):')
D0, D1 = RandomDMAP(4, 5.62, 10)
print('D0 = ')
print(D0)
print('D1 = ')
print(D1)
print('m = MarginalMomentsFromDMAP(D0, D1, 1)[0]:')
m = MarginalMomentsFromDMAP(D0, D1, 1)[0]
print('m = ')
print(m)
assert CheckDMAPRepresentation(D0, D1), "RandomDMAP failed to return a valid DMAP representation!"
assert np.abs(m-5.62)<10**-14, "RandomDMAP failed to match the given mean value!"
print('========================================')
print('Testing BuTools function RandomDMMAP')
print('Test:')
print('-----')
print('Dm = RandomDMMAP(4, 3, 5.62, 10):')
Dm = RandomDMMAP(4, 3, 5.62, 10)
print(Dm[0])
print(Dm[1])
print(Dm[2])
print(Dm[3])
print('m = MarginalMomentsFromDMMAP(Dm, 1)[0]:')
m = MarginalMomentsFromDMMAP(Dm, 1)[0]
print('m = ')
print(m)
assert CheckDMMAPRepresentation(Dm), "RandomDMMAP failed to return a valid MMAP representation!"
assert np.abs(m-5.62)<10**-14, "RandomDMMAP failed to match the given mean value!"
print('========================================')
print('Testing BuTools function CheckDMAPRepresentation')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0],[0, 0.17, 0.2],[0.16, 0.17, 0.02]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('flag = CheckDMAPRepresentation(D0, D1):')
flag = CheckDMAPRepresentation(D0, D1)
print('flag = ')
print(flag)
assert flag==False, "CheckDMAPRepresentation failed to detect the incompatible shapes of D0 and D1!"
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0],[0, 0.17, 0.2],[0.16, 0.17, 0.02]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1],[0.18, 0.07, 0.14],[0.13, 0.15, 0.15]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('flag = CheckDMAPRepresentation(D0, D1):')
flag = CheckDMAPRepresentation(D0, D1)
print('flag = ')
print(flag)
assert flag==False, "CheckDMAPRepresentation failed to detect invalid rowsums!"
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('flag = CheckDMAPRepresentation(D0, D1):')
flag = CheckDMAPRepresentation(D0, D1)
print('flag = ')
print(flag)
assert flag==True, "CheckDMAPRepresentation failed to recognize a valid DMAP representation!"
print('========================================')
print('Testing BuTools function CheckDRAPRepresentation')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02],[0.2, 0, 0]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1., -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29],[0, 0.8, 0]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==False, "CheckDRAPRepresentation failed to detect the incompatible shapes of D0 and D1!"
print('Input:')
print('------')
H0 = ml.matrix([[0.2, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1., -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==False, "CheckDRAPRepresentation failed to detect invalid rowsums!"
print('Input:')
print('------')
H0 = ml.matrix([[-1., 0, 0],[0, -2., 2.],[0, 3., -3.]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 0, 1.],[0, -1., 1.],[1., 0, -1.]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==False, "CheckDRAPRepresentation failed to detect invalid eigenvalues!"
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 15.],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1., -15.],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==False, "CheckDRAPRepresentation failed to detect invalid eigenvalues!"
print('Input:')
print('------')
H0 = ml.matrix([[0, 0.5, 0.1],[0, -1.4, 3.1],[0.67, 0, 0.4]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 0.4, 0],[0, -0.2, -0.5],[0.3, -0.7, 0.33]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==False, "CheckDRAPRepresentation failed to detect invalid eigenvalues!"
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1., -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('flag = CheckDRAPRepresentation(H0, H1):')
flag = CheckDRAPRepresentation(H0, H1)
print('flag = ')
print(flag)
assert flag==True, "CheckDRAPRepresentation failed to recognize a valid RAP representation!"
print('========================================')
print('Testing BuTools function CheckDMMAPRepresentation')
print('Input:')
print('------')
D0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('D1 = ')
print(D1)
D2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('D2 = ')
print(D2)
D3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('D3 = ')
print(D3)
print('Test:')
print('-----')
print('flag = CheckDMMAPRepresentation([D0, D1, D2, D3]):')
flag = CheckDMMAPRepresentation([D0, D1, D2, D3])
print('flag = ')
print(flag)
assert flag==True, "CheckDMMAPRepresentation failed to recognize a valid DMMAP representation!"
print('========================================')
print('Testing BuTools function CheckDMRAPRepresentation')
print('Input:')
print('------')
H0 = ml.matrix([[0.15, 0.2, 0.18],[-0.23, 0.17, 0.22],[0.19, 0.15, 0.16]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0.01, 0.08, 0.16],[0.02, 0.2, 0.07],[0.02, 0.15, 0.17]])
print('H1 = ')
print(H1)
H2 = ml.matrix([[0.14, 0.07, 0.01],[0.19, 0.02, 0.34],[0.06, 0.1, 0]])
print('H2 = ')
print(H2)
print('Test:')
print('-----')
print('flag = CheckDMRAPRepresentation([H0, H1, H2]):')
flag = CheckDMRAPRepresentation([H0, H1, H2])
print('flag = ')
print(flag)
assert flag==True, "CheckDMRAPRepresentation failed to recognize a valid DMRAP representation!"
print('========================================')
print('Testing BuTools function DRAPFromMoments')
print('Input:')
print('------')
G0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('G0 = ')
print(G0)
G1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('G1 = ')
print(G1)
print('moms = MarginalMomentsFromDRAP(G0, G1, 5):')
moms = MarginalMomentsFromDRAP(G0, G1, 5)
print('moms = ')
print(moms)
print('Nm = LagkJointMomentsFromDRAP(G0, G1, 2, 1):')
Nm = LagkJointMomentsFromDRAP(G0, G1, 2, 1)
print('Nm = ')
print(Nm)
print('Test:')
print('-----')
print('H0, H1 = DRAPFromMoments(moms, Nm):')
H0, H1 = DRAPFromMoments(moms, Nm)
print('H0 = ')
print(H0)
print('H1 = ')
print(H1)
print('rmoms = MarginalMomentsFromDRAP(H0, H1, 5):')
rmoms = MarginalMomentsFromDRAP(H0, H1, 5)
print('rmoms = ')
print(rmoms)
print('rNm = LagkJointMomentsFromDRAP(H0, H1, 2, 1):')
rNm = LagkJointMomentsFromDRAP(H0, H1, 2, 1)
print('rNm = ')
print(rNm)
assert la.norm(np.array(moms)-np.array(rmoms))<10**-11 and la.norm(Nm-rNm)<10**-12, "The moments and joint moments returned by RAPFromMoments are not the same as given!"
print('========================================')
print('Testing BuTools function DMRAPFromMoments')
print('Input:')
print('------')
G0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('G0 = ')
print(G0)
G1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('G1 = ')
print(G1)
G2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('G2 = ')
print(G2)
G3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('G3 = ')
print(G3)
print('G = [G0, G1, G2, G3]:')
G = [G0, G1, G2, G3]
print('moms = MarginalMomentsFromDMRAP(G, 5):')
moms = MarginalMomentsFromDMRAP(G, 5)
print('moms = ')
print(moms)
print('Nm = LagkJointMomentsFromDMRAP(G, 2, 1):')
Nm = LagkJointMomentsFromDMRAP(G, 2, 1)
print('Nm1, Nm2, Nm3 = Nm:')
Nm1, Nm2, Nm3 = Nm
print('Nm1 = ')
print(Nm1)
print('Nm2 = ')
print(Nm2)
print('Nm3 = ')
print(Nm3)
print('Test:')
print('-----')
print('H = DMRAPFromMoments(moms, Nm):')
H = DMRAPFromMoments(moms, Nm)
print('H[0]:')
print(H[0])
print('H[1]:')
print(H[1])
print('H[2]:')
print(H[2])
print('H[3]:')
print(H[3])
butools.checkPrecision = 10.**-10
print('rmoms = MarginalMomentsFromDMRAP(H, 5):')
rmoms = MarginalMomentsFromDMRAP(H, 5)
print('rmoms = ')
print(rmoms)
print('rNm = LagkJointMomentsFromDMRAP(H, 2, 1):')
rNm = LagkJointMomentsFromDMRAP(H, 2, 1)
print('rNm1, rNm2, rNm3 = rNm:')
rNm1, rNm2, rNm3 = rNm
print('rNm1 = ')
print(rNm1)
print('rNm2 = ')
print(rNm2)
print('rNm3 = ')
print(rNm3)
print('la.norm(np.array(moms)-np.array(rmoms)):')
print(la.norm(np.array(moms)-np.array(rmoms)))
assert la.norm(np.array(moms)-np.array(rmoms))<10**-9 and la.norm(Nm1-rNm1)<10**-9 and la.norm(Nm2-rNm2)<10**-9 and la.norm(Nm3-rNm3)<10**-9, "The moments and joint moments returned by DMRAPFromMoments are not the same as given!"
print('========================================')
print('Testing BuTools function DMAPFromDRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0, 0, 0.13],[0, 0.6, 0.18],[0.31, 0.26, 0.02]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0, 1, -0.13],[0, 0.18, 0.04],[0.03, 0.09, 0.29]])
print('H1 = ')
print(H1)
print('Test:')
print('-----')
print('D0, D1 = DMAPFromDRAP(H0, H1):')
D0, D1 = DMAPFromDRAP(H0, H1)
print('D0 = ')
print(D0)
print('D1 = ')
print(D1)
print('err = la.norm(LagkJointMomentsFromDRAP(D0, D1, 3, 1)-LagkJointMomentsFromDRAP(H0, H1, 3, 1)):')
err = la.norm(LagkJointMomentsFromDRAP(D0, D1, 3, 1)-LagkJointMomentsFromDRAP(H0, H1, 3, 1))
print('err = ')
print(err)
assert err<10**-10, "The DMAP returned by DMAPFromDRAP is not similar to the input!"
assert CheckDMAPRepresentation(D0, D1), "The result of DMAPFromDRAP is not a DMAP, as it should be!"
print('========================================')
print('Testing BuTools function DMMAPFromDMRAP')
print('Input:')
print('------')
H0 = ml.matrix([[0.15, 0.2, 0.18],[-0.20, 0.17, 0.22],[0.19, 0.15, 0.16]])
print('H0 = ')
print(H0)
H1 = ml.matrix([[0.01, 0.08, 0.16],[0.02, 0.2, 0.07],[0.02, 0.15, 0.17]])
print('H1 = ')
print(H1)
H2 = ml.matrix([[0.14, 0.07, 0.01],[0.19, 0.02, 0.31],[0.06, 0.1, 0.]])
print('H2 = ')
print(H2)
H = [H0, H1, H2]
print('moms = MarginalMomentsFromDMRAP(H):')
moms = MarginalMomentsFromDMRAP(H)
print('moms = ')
print(moms)
print('jmom = LagkJointMomentsFromDMRAP(H, 3, 1):')
jmom = LagkJointMomentsFromDMRAP(H, 3, 1)
print('Test:')
print('-----')
print('G = DMMAPFromDMRAP(H):')
G = DMMAPFromDMRAP(H)
print(G[0])
print(G[1])
print(G[2])
print('rmoms = MarginalMomentsFromDMMAP(G):')
rmoms = MarginalMomentsFromDMMAP(G)
print('rmoms = ')
print(rmoms)
print('rjmom = LagkJointMomentsFromDMMAP(G, 3, 1):')
rjmom = LagkJointMomentsFromDMMAP(G, 3, 1)
print('err = la.norm(rjmom[0]-jmom[0])+la.norm(rjmom[1]-jmom[1]):')
err = la.norm(rjmom[0]-jmom[0])+la.norm(rjmom[1]-jmom[1])
print('err = ')
print(err)
assert err<10**-12, "The DMMAP returned by DMMAPFromDMRAP is not similar to the input!"
assert CheckDMMAPRepresentation(G), "The result of DMMAPFromDMRAP is not a DMMAP, as it should be!"
print('========================================')
print('Testing BuTools function CanonicalFromDMAP2')
print('Input:')
print('------')
D0 = ml.matrix([[0.46, 0.28],[0.35, 0.23]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.08, 0.18],[0.14, 0.28]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('H0, H1 = CanonicalFromDMAP2(D0, D1):')
H0, H1 = CanonicalFromDMAP2(D0, D1)
print('H0 = ')
print(H0)
print('H1 = ')
print(H1)
Cm = SimilarityMatrix(H0, D0)
err1 = la.norm(H0*Cm-Cm*D0)
err2 = la.norm(H1*Cm-Cm*D1)
print("Transformation errors:")
print(np.max(err1, err2))
assert CheckDMAPRepresentation(H0, H1), "The result of CanonicalFromDMAP2 is not a valid DMAP representation!"
assert np.max(err1, err2)<10**-12, "The DMAP returned by CanonicalFromDMAP2 is not similar to the input!"
print('Input:')
print('------')
D0 = ml.matrix([[0.26, 0.28],[0.35, 0.23]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.28, 0.18],[0.14, 0.28]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('H0, H1 = CanonicalFromDMAP2(D0, D1):')
H0, H1 = CanonicalFromDMAP2(D0, D1)
print('H0 = ')
print(H0)
print('H1 = ')
print(H1)
Cm = SimilarityMatrix(H0, D0)
err1 = la.norm(H0*Cm-Cm*D0)
err2 = la.norm(H1*Cm-Cm*D1)
print("Transformation errors:")
print(np.max(err1, err2))
assert CheckDMAPRepresentation(H0, H1), "The result of CanonicalFromDMAP2 is not a valid DMAP representation!"
assert np.max(err1, err2)<10**-12, "The DMAP returned by CanonicalFromDMAP2 is not similar to the input!"
print('Input:')
print('------')
D0 = ml.matrix([[0.14, 0.34],[0.35, 0.23]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.22, 0.3],[0.28, 0.14]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('H0, H1 = CanonicalFromDMAP2(D0, D1):')
H0, H1 = CanonicalFromDMAP2(D0, D1)
print('H0 = ')
print(H0)
print('H1 = ')
print(H1)
Cm = SimilarityMatrix(H0, D0)
err1 = la.norm(H0*Cm-Cm*D0)
err2 = la.norm(H1*Cm-Cm*D1)
print("Transformation errors:")
print(np.max(err1, err2))
assert CheckDMAPRepresentation(H0, H1), "The result of CanonicalFromDMAP2 is not a valid DMAP representation!"
assert np.max(err1, err2)<10**-12, "The DMAP returned by CanonicalFromDMAP2 is not similar to the input!"
print('========================================')
print('Testing BuTools function DMAP2FromMoments')
print('Input:')
print('------')
D0 = ml.matrix([[0.2, 0.7],[0.6, 0.1]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.09, 0.01],[0.2, 0.1]])
print('D1 = ')
print(D1)
print('moms = MarginalMomentsFromDMAP(D0, D1, 3):')
moms = MarginalMomentsFromDMAP(D0, D1, 3)
print('moms = ')
print(moms)
print('corr = LagCorrelationsFromDMAP(D0, D1, 1)[0]:')
corr = LagCorrelationsFromDMAP(D0, D1, 1)[0]
print('corr = ')
print(corr)
print('Test:')
print('-----')
print('D0, D1 = DMAP2FromMoments(moms, corr):')
D0, D1 = DMAP2FromMoments(moms, corr)
print('D0 = ')
print(D0)
print('D1 = ')
print(D1)
print('rmoms = MarginalMomentsFromDMAP(D0, D1, 3):')
rmoms = MarginalMomentsFromDMAP(D0, D1, 3)
print('rmoms = ')
print(rmoms)
print('rcorr = LagCorrelationsFromDMAP(D0, D1, 1)[0]:')
rcorr = LagCorrelationsFromDMAP(D0, D1, 1)[0]
print('rcorr = ')
print(rcorr)
assert CheckDMAPRepresentation(D0, D1), "DMAP2FromMoments returned an invalid DMAP representation!"
assert la.norm(np.array(moms)-np.array(rmoms))<10**-11 and np.abs(corr-rcorr)<10**-11, "The moments and the correlation returned by DMAP2FromMoments are not the same as given!"
print('========================================')
print('Testing BuTools function SamplesFromDMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0, 0.02, 0, 0],[0, 0.17, 0.2, 0.14],[0.16, 0.17, 0.02, 0.18],[0, 0, 0, 0.12]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0, 0.88, 0.1, 0],[0.18, 0.07, 0.14, 0.1],[0.13, 0.15, 0.15, 0.04],[0.31, 0.18, 0.12, 0.27]])
print('D1 = ')
print(D1)
print('Test:')
print('-----')
print('x = SamplesFromDMAP(D0, D1, 10000):')
x = SamplesFromDMAP(D0, D1, 10000)
print('mt = MarginalMomentsFromTrace(x, 3):')
mt = MarginalMomentsFromTrace(x, 3)
print('mt = ')
print(mt)
print('mm = MarginalMomentsFromDMAP(D0, D1, 3):')
mm = MarginalMomentsFromDMAP(D0, D1, 3)
print('mm = ')
print(mm)
print('========================================')
print('Testing BuTools function SamplesFromDMMAP')
print('Input:')
print('------')
D0 = ml.matrix([[0.34, 0, 0],[0.06, 0.05, 0.03],[0.11, 0.13, 0]])
print('D0 = ')
print(D0)
D1 = ml.matrix([[0.3, 0, 0],[0.16, 0.18, 0.05],[0.15, 0.04, 0.09]])
print('D1 = ')
print(D1)
D2 = ml.matrix([[0, 0.01, 0],[0.1, 0.07, 0.08],[0.13, 0.12, 0.13]])
print('D2 = ')
print(D2)
D3 = ml.matrix([[0.35, 0, 0],[0, 0.18, 0.04],[0.06, 0.03, 0.01]])
print('D3 = ')
print(D3)
Dm = [D0, D1, D2, D3]
print('Test:')
print('-----')
print('x = SamplesFromDMMAP(Dm, 10000):')
x = SamplesFromDMMAP(Dm, 10000)
print('mt = MarginalMomentsFromTrace(x[:, 0], 3):')
mt = MarginalMomentsFromTrace(x[:, 0], 3)
print('mt = ')
print(mt)
print('mm = MarginalMomentsFromDMMAP(Dm, 3):')
mm = MarginalMomentsFromDMMAP(Dm, 3)
print('mm = ')
print(mm)

