"""Various visualizations and summary statistics for a one layer representation"""

if column_label is None:

column_label = list(map(str, range(data.shape[1])))

if row_label is None:

row_label = list(map(str, range(corex.n_samples)))

alpha = corex.alpha[:, :, 0]

print('Groups in sorted_groups.txt')

output_groups(corex.tcs, alpha, corex.mis, column_label, prefix=prefix)

output_labels(corex.labels, row_label, prefix=prefix)

output_strong(corex.tcs, alpha, corex.mis, corex.labels, prefix=prefix)

anomalies(corex.log_z, row_label=row_label, prefix=prefix)

plot_convergence(corex.tc_history, prefix=prefix)

if data is not None:

print('Pairwise plots among high TC variables in "relationships"')

data_to_plot = np.where(data == corex.missing_values, np.nan, data)

log_p_y_given_x = calculate_log_latent(corex, data)

cont = cont3(log_p_y_given_x)

output_cont_labels(log_p_y_given_x, row_label, prefix=prefix)

plot_heatmaps(data_to_plot, corex.labels, alpha, corex.mis, column_label, cont, prefix=prefix, focus=focus)

plot_pairplots(data_to_plot, corex.labels, alpha, corex.mis, column_label, prefix=prefix, focus=focus, topk=topk)

plot_top_relationships(data_to_plot, corex.labels, alpha, corex.mis, column_label, cont, prefix=prefix, topk=topk)

"""Visualize a hierarchy of representations."""

if column_label is None:

column_label = list(map(str, range(corexes[0].alpha.shape[1])))

if row_label is None:

row_label = list(map(str, range(corexes[0].labels.shape[0])))

f = safe_open(prefix + '/text_files/higher_layer_group_tcs.txt', 'w+')

params = ['dim_hidden', 'eps', 'marginal_description', 'max_iter', 'max_samples', 'missing_values', 'n_cpu',

'n_hidden', 'n_repeat', 'n_samples', 'n_visible', 'ram', 'smooth_marginals']

parameter_dict = {}

for j, corex in enumerate(corexes):

parameter_dict[j] = {}

for param in params:

parameter_dict[j][param] = getattr(corex, param, None)

f.write('At layer: %d, Total TC: %0.3f\n' % (j, corex.tc))

f.write('Individual TCS:' + str(corex.tcs) + '\n')

plot_convergence(corex.tc_history, prefix=prefix, prefix2=j)

g = safe_open('{}/text_files/mis_layer{}.csv'.format(prefix, j), 'w+')

h = safe_open('{}/text_files/weights_layer{}.csv'.format(prefix, j), 'w+')

if j == 0:

g.write('factor,' + ','.join(column_label) + '\n')

h.write('factor,' + ','.join(column_label) + '\n')

else:

g.write('factor,'+ ','.join(map(str, list(range(len(corex.mis[0,:]))))) + '\n')

h.write('factor,'+ ','.join(map(str, list(range(len(corex.mis[0,:]))))) + '\n')

mis = corex.mis / np.log(2)

alpha = corex.alpha

for ir, r in enumerate(mis):

g.write(str(ir) + ',' + ','.join(map(str, mis[ir])) + '\n')

h.write(str(ir) + ',' + ','.join(map(str, mis[ir] * alpha[ir].ravel())) + '\n')

g.close()

h.close()

f.close()

with open('{}/text_files/parameters.json'.format(prefix), 'w') as fp:

json.dump(parameter_dict, fp)

import textwrap

column_label = list(map(lambda q: '\n'.join(textwrap.wrap(q, width=17, break_long_words=False)), column_label))

# Construct non-tree graph

weights = [corex.alpha[:, :, 0].clip(0, 1) * corex.mis for corex in corexes]

node_weights = [corex.tcs for corex in corexes]

g = make_graph(weights, node_weights, column_label, max_edges=max_edges)

# Display pruned version

h = g.copy() # trim(g.copy(), max_parents=max_parents, max_children=max_children)

edge2pdf(h, prefix + '/graphs/graph_prune_' + str(max_edges), labels='label', directed=True, makepdf=True)

# Display tree version

tree = g.copy()

tree = trim(tree, max_parents=1, max_children=False)

edge2pdf(tree, prefix + '/graphs/tree', labels='label', directed=True, makepdf=True)

cmap = sns.cubehelix_palette(as_cmap=True, light=.9)

m, nv = mis.shape

for j in range(m):

inds = np.where(np.logical_and(alpha[j] > 0, mis[j] > 0.))[0]

inds = inds[np.argsort(- alpha[j, inds] * mis[j, inds])][:topk]

if focus in column_label:

ifocus = column_label.index(focus)

if not ifocus in inds:

inds = np.insert(inds, 0, ifocus)

if len(inds) >= 2:

plt.clf()

order = np.argsort(cont[:,j])

subdata = data[:, inds][order].T

subdata -= np.nanmean(subdata, axis=1, keepdims=True)

subdata /= np.nanstd(subdata, axis=1, keepdims=True)

columns = [column_label[i] for i in inds]

sns.heatmap(subdata, vmin=-3, vmax=3, cmap=cmap, yticklabels=columns, xticklabels=False, mask=np.isnan(subdata))

filename = '{}/heatmaps/group_num={}.png'.format(prefix, j)

if not os.path.exists(os.path.dirname(filename)):

os.makedirs(os.path.dirname(filename))

plt.title("Latent factor {}".format(j))

plt.savefig(filename, bbox_inches='tight')

plt.close('all')

#plot_rels(data[:, inds], list(map(lambda q: column_label[q], inds)), colors=cont[:, j],

cmap = sns.cubehelix_palette(as_cmap=True, light=.9)

plt.rcParams.update({'font.size': 32})

m, nv = mis.shape

for j in range(m):

inds = np.where(np.logical_and(alpha[j] > 0, mis[j] > 0.))[0]

inds = inds[np.argsort(- alpha[j, inds] * mis[j, inds])][:topk]

if focus in column_label:

ifocus = column_label.index(focus)

if not ifocus in inds:

inds = np.insert(inds, 0, ifocus)

if len(inds) >= 2:

plt.clf()

subdata = data[:, inds]

columns = [column_label[i] for i in inds]

subdata = pd.DataFrame(data=subdata, columns=columns)

try:

sns.pairplot(subdata, kind="reg", diag_kind="kde", height=5, dropna=True)

filename = '{}/pairplots_regress/group_num={}.pdf'.format(prefix, j)

if not os.path.exists(os.path.dirname(filename)):

os.makedirs(os.path.dirname(filename))

plt.suptitle("Latent factor {}".format(j), y=1.01)

plt.savefig(filename, bbox_inches='tight')

plt.clf()

except:

pass

subdata['Latent factor'] = labels[:,j]

try:

sns.pairplot(subdata, kind="scatter", dropna=True, vars=subdata.columns.drop('Latent factor'), hue="Latent factor", diag_kind="kde", height=5)

filename = '{}/pairplots/group_num={}.pdf'.format(prefix, j)

if not os.path.exists(os.path.dirname(filename)):

os.makedirs(os.path.dirname(filename))

plt.suptitle("Latent factor {}".format(j), y=1.01)

plt.savefig(filename, bbox_inches='tight')

plt.close('all')

except:

all_edges = np.hstack(list(map(np.ravel, weights)))

max_edges = min(max_edges, len(all_edges))

w_thresh = np.sort(all_edges)[-max_edges]

print('weight threshold is %f for graph with max of %f edges ' % (w_thresh, max_edges))

g = nx.DiGraph()

max_node_weight = max([max(w) for w in node_weights])

for layer, weight in enumerate(weights):

m, n = weight.shape

for j in range(m):

g.add_node((layer + 1, j))

g._node[(layer + 1, j)]['weight'] = 0.3 * node_weights[layer][j] / max_node_weight

for i in range(n):

if weight[j, i] > w_thresh:

if weight[j, i] > w_thresh / 2:

g.add_weighted_edges_from([( (layer, i), (layer + 1, j), 10 * weight[j, i])])

else:

g.add_weighted_edges_from([( (layer, i), (layer + 1, j), 0)])

# Label layer 0

for i, lab in enumerate(column_label):

g.add_node((0, i))

g._node[(0, i)]['label'] = lab

g._node[(0, i)]['name'] = lab # JSON uses this field

g._node[(0, i)]['weight'] = 1

if float(nx.__version__) < 2:

edgedict = g.edge

else:

edgedict = g.adj

for node in g:

if max_parents:

parents = list(g.successors(node))

weights = [edgedict[node][parent]['weight'] for parent in parents]

for weak_parent in np.argsort(weights)[:-max_parents]:

g.remove_edge(node, parents[weak_parent])

if max_children:

children = g.predecessors(node)

weights = [edgedict[child][node]['weight'] for child in children]

for weak_child in np.argsort(weights)[:-max_children]:

g.remove_edge(children[weak_child], node)

f = safe_open(prefix + '/text_files/groups.txt', 'w+')

g = safe_open(prefix + '/text_files/groups_no_overlaps.txt', 'w+')

m, nv = mis.shape

for j in range(m):

f.write('Group num: %d, TC(X;Y_j): %0.3f\n' % (j, tcs[j]))

g.write('Group num: %d, TC(X;Y_j): %0.3f\n' % (j, tcs[j]))

inds = np.where(alpha[j] * mis[j] > thresh)[0]

inds = inds[np.argsort(-alpha[j, inds] * mis[j, inds])]

for ind in inds:

f.write(column_label[ind] + ', %0.3f, %0.3f, %0.3f\n' % (

mis[j, ind], alpha[j, ind], mis[j, ind] * alpha[j, ind]))

inds = np.where(alpha[j] == 1)[0]

inds = inds[np.argsort(- mis[j, inds])]

for ind in inds:

g.write(column_label[ind] + ', %0.3f\n' % mis[j, ind])

f.close()

f = safe_open(prefix + '/text_files/labels.txt', 'w+')

ns, m = labels.shape

for l in range(ns):

f.write(row_label[l] + ',' + ','.join(map(str, labels[l, :])) + '\n')

f = safe_open(prefix + '/text_files/cont_labels.txt', 'w+')

m, ns, k = log_p_y_given_x.shape

# assert k==2, 'More complicated if k>2... use cont3_test to generate if k==3'

labels = cont3(log_p_y_given_x)

for l in range(ns):

f.write(row_label[l] + ',' + ','.join(map(lambda q: '%0.6f' % q, labels[l, :])) + '\n')

f = safe_open(prefix + '/text_files/most_deterministic_groups.txt', 'w+')

m, n = alpha.shape

topk = 5

ixy = np.clip(np.sum(alpha * mis, axis=1) - tcs, 0, np.inf)

hys = np.array([entropy(labels[:, j]) for j in range(m)]).clip(1e-6)

ntcs = [(np.sum(np.sort(alpha[j] * mis[j])[-topk:]) - ixy[j]) / ((topk - 1) * hys[j]) for j in range(m)]

f.write('Group num., NTC\n')

for j, ntc in sorted(enumerate(ntcs), key=lambda q: -q[1]):

f.write('%d, %0.3f\n' % (j, ntc))

m, nv = mis.shape

for j in range(m):

# inds = np.where(alpha[j] * mis[j] > 0)[0]

#inds = inds[np.argsort(-alpha[j, inds] * mis[j, inds])][:topk]

inds = np.where(alpha[j] > athresh)[0]

inds = inds[np.argsort(- alpha[j, inds] * mis[j, inds])][:topk]

if len(inds) >= 2:

plot_rels(data[:, inds], list(map(lambda q: column_label[q], inds)), colors=cont[:, j],

from scipy.special import erf

ns = log_z.shape[1]

if row_label is None:

row_label = list(map(str, range(ns)))

a_score = np.sum(log_z[:, :, 0], axis=0)

mean, std = np.mean(a_score), np.std(a_score)

a_score = (a_score - mean) / std

percentile = 1. / ns

anomalies = np.where(0.5 * (1 - erf(a_score / np.sqrt(2)) ) < percentile)[0]

f = safe_open(prefix + '/text_files/anomalies.txt', 'w+')

for i in anomalies:

f.write(row_label[i] + ', %0.1f\n' % a_score[i])

# Return indices and data matrix only for rows with at least 50% of the columns non-empty

ns, nv = data.shape

inds = np.sum(data == missing_values, axis=1) < fraction * nv

# sample entropy for one discrete var

xsamples = np.asarray(xsamples)

xsamples = xsamples[xsamples >= 0] # by def, -1 means missing value

xs = np.unique(xsamples)

ns = len(xsamples)

ps = np.array([float(np.count_nonzero(xsamples == x)) / ns for x in xs])

if directed:

os.system(

"sfdp " + fname + ".dot -Tpdf -Earrowhead=none -Nfontsize=12 -GK=2 -Gmaxiter=1000 -Goverlap=False -Gpack=True -Gpackmode=clust -Gsep=0.01 -Gsplines=False -o " + fname + "_sfdp.pdf")

os.system(

"sfdp " + fname + ".dot -Tpdf -Earrowhead=none -Nfontsize=12 -GK=2 -Gmaxiter=1000 -Goverlap=False -Gpack=True -Gpackmode=clust -Gsep=0.01 -Gsplines=True -o " + fname + "_sfdp_w_splines.pdf")

return True

if position is None:

os.system("neato " + fname + ".dot -Tpdf -o " + fname + ".pdf")

os.system("fdp " + fname + ".dot -Tpdf -o " + fname + "fdp.pdf")

else:

os.system("neato " + fname + ".dot -Tpdf -n -o " + fname + ".pdf")

colors = 'indigo,gold,hotpink,firebrick,indianred,yellow,mistyrose,darkolivegreen,darkseagreen,pink,tomato,lightcoral,orangered,navajowhite,palegreen,darkslategrey,greenyellow,burlywood,seashell,mediumspringgreen,papayawhip,blanchedalmond,chartreuse,dimgray,black,peachpuff,springgreen,aquamarine,white,orange,lightsalmon,darkslategray,brown,ivory,dodgerblue,peru,lawngreen,chocolate,crimson,forestgreen,slateblue,lightseagreen,cyan,mintcream,antiquewhite,mediumorchid,skyblue,gray,darkturquoise,goldenrod,darkgreen,floralwhite,darkviolet,moccasin,saddlebrown,grey,darkslateblue,lightskyblue,lightpink,mediumvioletred,slategrey,red,deeppink,limegreen,palegoldenrod,plum,turquoise,lightgrey,lightgoldenrodyellow,darkgoldenrod,lavender,maroon,yellowgreen,sandybrown,thistle,violet,navy,magenta,dimgrey,tan,rosybrown,blue,lightblue,ghostwhite,honeydew,cornflowerblue,linen,powderblue,seagreen,darkkhaki,snow,sienna,mediumblue,royalblue,lightcyan,green,mediumpurple,midnightblue,cornsilk,paleturquoise,bisque,slategray,khaki,wheat,darkorchid,deepskyblue,salmon,steelblue,palevioletred,lightslategray,aliceblue,lightslategrey,orchid,gainsboro,mediumseagreen,lightgray,mediumturquoise,lemonchiffon,cadetblue,lightyellow,lavenderblush,coral,purple,whitesmoke,mediumslateblue,darkorange,mediumaquamarine,darksalmon,beige,blueviolet,azure,lightsteelblue,oldlace'.split(',')

parts = label.split('_')

for part in parts:

if part in colors:

parts.remove(part)

return '_'.join(parts), part

#This function will takes list of edges and a filename

#and write a file in .dot format. Readable, eg. by omnigraffle

# OR use "neato file.dot -Tpng -n -o file.png"

# The -n option says whether to use included node positions or to generate new ones

# for a grid, positions = [(i%28,i/28) for i in range(784)]

def cnn(node):

#change node names for dot format

if type(node) is tuple or type(node) is list:

return u'n' + u'_'.join(map(str, node))

else:

return str(node)

if connected:

touching = list(set(sum([[a, b] for a, b in g.edges()], [])))

g = nx.subgraph(g, touching)

print('non-isolated nodes,edges', len(list(g.nodes())), len(list(g.edges())))

f = safe_open(filename + '.dot', 'wb+')

if directed:

f.write("strict digraph {\n".encode('utf-8'))

else:

f.write("strict graph {\n".encode('utf-8'))

#f.write("\tgraph [overlap=scale];\n".encode('utf-8'))

f.write("\tnode [shape=point];\n".encode('utf-8'))

for a, b, d in g.edges(data=True):

if 'weight' in d:

if directed:

f.write(("\t" + cnn(a) + ' -> ' + cnn(b) + ' [penwidth=%.2f' % float(

np.clip(d['weight'], 0, 9)) + '];\n').encode('utf-8'))

else:

if d['weight'] > threshold:

f.write(("\t" + cnn(a) + ' -- ' + cnn(b) + ' [penwidth=' + str(3 * d['weight']) + '];\n').encode(

'utf-8'))

else:

if directed:

f.write(("\t" + cnn(a) + ' -> ' + cnn(b) + ';\n').encode('utf-8'))

else:

f.write(("\t" + cnn(a) + ' -- ' + cnn(b) + ';\n').encode('utf-8'))

for n in g.nodes():

if labels is not None:

if type(labels) == dict or type(labels) == list:

thislabel = labels[n].replace(u'"', u'\\"')

lstring = u'label="' + thislabel + u'",shape=none'

elif type(labels) == str:

if 'label' in g._node[n]:

thislabel = g._node[n][labels].replace(u'"', u'\\"')

# combine dupes

#llist = thislabel.split(',')

#thislabel = ','.join([l for l in set(llist)])

thislabel, thiscolor = extract_color(thislabel)

lstring = u'label="%s",shape=none,fontcolor="%s"' % (thislabel, thiscolor)

else:

weight = g._node[n].get('weight', 0.1)

if n[0] == 1:

lstring = u'shape=circle,margin="0,0",style=filled,fillcolor=black,fontcolor=white,height=%0.2f,label="%d"' % (

2 * weight, n[1])

else:

lstring = u'shape=point,height=%0.2f' % weight

else:

lstring = 'label="' + str(n) + '",shape=none'

#lstring = unicode(lstring)

else:

lstring = False

if position is not None:

if position == 'grid':

position = [(i % 28, 28 - i / 28) for i in range(784)]

posstring = 'pos="' + str(position[n][0]) + ',' + str(position[n][1]) + '"'

else:

posstring = False

finalstring = u' [' + u','.join([ts for ts in [posstring, lstring] if ts]) + u']\n'

#finalstring = u' ['+lstring+u']\n'

f.write((u'\t' + cnn(n) + finalstring).encode('utf-8'))

f.write("}".encode('utf-8'))

f.close()

if makepdf:

neato(filename, position=position, directed=directed)

tvalue=0.1):

alpha, labels, lpygx, mis, lasttc = out[:5]

ns, m = labels.shape

m, nv = mis.shape

hys = [entropy(labels[:, j]) for j in range(m)]

#alpha = np.array([z[2] for z in zs]) # m by nv

nmis = []

ixys = []

for j in range(m):

if hys[j] > 0:

#ixy = np.dot((alpha[j]>0.95).astype(int),mis[j])-lasttc[-1][j]

ixy = max(0., np.dot(alpha[j], mis[j]) - lasttc[-1][j])

ixys.append(ixy)

tcn = (np.sum(np.sort(alpha[j] * mis[j])[-topk:]) - ixy) / ((topk - 1) * hys[j])

nmis.append(tcn) #ixy) #/hys[j])

else:

ixys.append(0)

nmis.append(0)

f = safe_open(prefix + outfile, 'w+')

print(list(enumerate(np.argsort(-np.array(nmis)))))

print(','.join(map(str, list(np.argsort(-np.array(nmis))))))

for i, top in enumerate(np.argsort(-np.array(nmis))):

f.write('Group num: %d, Score: %0.3f\n' % (top, nmis[top]))

inds = np.where(alpha[top] > athresh)[0]

inds = inds[np.argsort(-mis[top, inds])]

for ind in inds:

f.write(wdict[ind] + ', %0.3f\n' % (mis[top, ind] / np.log(2)))

if wdict:

print(','.join(map(lambda q: wdict[q], inds)))

print(','.join(map(str, inds)))

print(top)

print(nmis[top], ixys[top], hys[top], ixys[top] / hys[top])

#,lasttc[-1][top],hys[top],lasttc[-1][top]/hys[top]

if graphs:

print(inds)

if len(inds) >= 2:

plot_rels(data[:, inds[:5]], list(map(lambda q: wdict[q], inds[:5])),

outfile='relationships/' + str(i) + '_group_num=' + str(top), latent=out[1][:, top],

alpha=tvalue)

f.close()

s, _ = extract_color(s)

if len(s) > 2 * n:

return s[:n] + '..' + s[-n:]

pylab.plot(tc_history)

pylab.xlabel('# iterations')

filename = '{}/text_files/convergence{}.pdf'.format(prefix, prefix2)

if not os.path.exists(os.path.dirname(filename)):

os.makedirs(os.path.dirname(filename))

pylab.savefig(filename)

pylab.close('all')

ns, n = data.shape

if labels is None:

labels = list(map(str, range(n)))

ncol = 5

# ncol = 4

nrow = int(np.ceil(float(n * (n - 1) / 2) / ncol))

#nrow=1

#pylab.rcParams.update({'figure.autolayout': True})

fig, axs = pylab.subplots(nrow, ncol)

fig.set_size_inches(5 * ncol, 5 * nrow)

#fig.set_canvas(pylab.gcf().canvas)

pairs = list(combinations(range(n), 2)) #[:4]

pairs = sorted(pairs, key=lambda q: q[0]**2+q[1]**2) # Puts stronger relationships first

if colors is not None:

colors = (colors - np.min(colors)) / (np.max(colors) - np.min(colors)).clip(1e-7)

for ax, pair in zip(axs.flat, pairs):

if latent is None:

ax.scatter(data[:, pair[0]], data[:, pair[1]], marker='.', edgecolors='none', alpha=alpha)

else:

# cs = 'rgbcmykrgbcmyk'

markers = 'x+.o,<>^^<>,+x.'

for j, ind in enumerate(np.unique(latent)):

inds = (latent == ind)

ax.scatter(data[inds, pair[0]], data[inds, pair[1]], c=colors[inds], cmap=pylab.get_cmap("jet"),

marker=markers[j], alpha=0.5, edgecolors='none', vmin=0, vmax=1)

ax.set_xlabel(shorten(labels[pair[0]]))

ax.set_ylabel(shorten(labels[pair[1]]))

for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:

ax.scatter(data[:, 0], data[:, 1], marker='.')

pylab.rcParams['font.size'] = 12 #6

pylab.draw()

#fig.set_tight_layout(True)

fig.tight_layout()

for ax in axs.flat[axs.size - 1:len(pairs) - 1:-1]:

ax.set_visible(False)

filename = outfile + '.png'

if not os.path.exists(os.path.dirname(filename)):

os.makedirs(os.path.dirname(filename))

fig.savefig(outfile + '.png') #df')

pylab.close('all')

"""

nv, ns, k = log_p_y_given_x.shape

p_y_given_x = np.exp(log_p_y_given_x)

allcont = []

if k == 3:

for j in range(nv):

ends = sorted([(np.dot(p_y_given_x[j, :, end[0]], p_y_given_x[j, :, end[1]]), end) for end in

[(0, 1), (0, 2), (1, 2)]])[0][1]

cont = log_p_y_given_x[j, :, ends[0]] - log_p_y_given_x[j, :, ends[1]]

cont = np.where(np.isnan(cont), 0, cont)

allcont.append(cont)

out = np.array(allcont).T # Bounds approximately reflect log of float precision

top = np.max(np.abs(np.where(np.isinf(out), 0, out)))

out = np.clip(out, -top, top)

elif k == 2:

for j in range(nv):

Y = p_y_given_x[j, :, 0] #np.log(p_y_given_x[j,:,0].clip(1e-40))

allcont.append(Y)

out = np.array(allcont).T

else:

print('error, not able to visualize with k > 3')

"""

nv, ns, k = p_y_given_x.shape

allcont = []

assert k==3

for j in range(nv):

ends = sorted([(np.dot(p_y_given_x[j, :, end[0]], p_y_given_x[j, :, end[1]]), end) for end in

[(0, 1), (0, 2), (1, 2)]])[0][1]

if logp:

cont = p_test[j, :, ends[0]] - p_test[j, :, ends[1]]

else:

cont = np.log(p_test[j, :, ends[0]]) - np.log(p_test[j, :, ends[1]])

cont = np.where(np.isnan(cont), 0, cont)

allcont.append(cont)

out = np.array(allcont).T # Bounds approximately reflect log of float precision

top = np.max(np.abs(np.where(np.isinf(out), 0, out)))

out = np.clip(out, -top, top)

""""Calculate the log probabilities for hidden factors for each sample, with high precision."""

Xm = ma.masked_equal(X, corex.missing_values)

n_samples, n_visible = Xm.shape

n_hidden, dim_hidden, ram = corex.n_hidden, corex.dim_hidden, corex.ram

log_p_y_given_x_unnorm = np.empty((n_hidden, n_samples, dim_hidden))

memory_size = float(n_samples * n_visible * n_hidden * dim_hidden * 64) / 1000**3 # GB

batch_size = np.clip(int(ram * n_samples / memory_size), 1, n_samples)

for l in range(0, n_samples, batch_size):

log_marg_x = corex.calculate_marginals_on_samples(corex.theta, Xm[l:l+batch_size]) # LLRs for each sample, for each var.

log_p_y_given_x_unnorm[:, l:l+batch_size, :] = corex.log_p_y + np.einsum('ikl,ijkl->ijl', corex.alpha, log_marg_x, optimize=False)

"""We need to normalize log p(y|x) with high precision... This is a hack to recalculate those since

log_z = logsumexp(log_p_y_given_x_unnorm, axis=2) # Essential to maintain precision.

log_z = log_z.reshape((n_hidden, -1, 1))

if z=='':

return str(options.missing)

else: