Refactor out some of the "op" calls #8

ManifoldEM/CC/runGlobalOptimization.py

@@ -1,4 +1,3 @@
-import logging
 import pickle
 
 import numpy as np
@@ -7,8 +6,6 @@
 from ManifoldEM.CC import MRFGeneratePotentials, MRFBeliefPropagation
 from ManifoldEM.CC.MRFBeliefPropagation import createBPalg
 
-_logger = logging.getLogger(__name__)
-_logger.setLevel(logging.DEBUG)
 '''
 function list = rearrange(seeds,nn)
 % function to return a list of nodes ordered according to the

ManifoldEM/DMembeddingII.py

@@ -3,7 +3,8 @@
 from collections import namedtuple
 from scipy.sparse import csr_matrix
 
-from ManifoldEM import p, fergusonE, sembeddingonFly
+from ManifoldEM import p, sembeddingonFly
+from ManifoldEM.core import fergusonE
 '''
 Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -133,7 +134,7 @@ def op(D, k, tune, prefsigma):  #*arg
     resnorm = np.inf
 
     logEps = np.arange(-150, 150.2, 0.2)
-    popt, logSumWij, resnorm, R_squared = fergusonE.op(np.sqrt(yVal), logEps, a0)
+    popt, logSumWij, resnorm, R_squared = fergusonE(np.sqrt(yVal), logEps, a0)
     nS = D.shape[0]
     nEigs = min(p.num_eigs, nS - 3)  #number of eigenfunctions to compute
     nA = 0  #autotuning parameter

ManifoldEM/NLSA.py

@@ -5,7 +5,8 @@
 from scipy.fftpack import fft2
 from scipy.fftpack import ifft2
 
-from ManifoldEM import DMembeddingII, myio, get_wiener, L2_distance, svdRF
+from ManifoldEM import DMembeddingII, myio, get_wiener
+from ManifoldEM.core import L2_distance, svdRF
 from ManifoldEM.fit_1D_open_manifold_3D import fit_1D_open_manifold_3D
 '''
 Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
@@ -84,7 +85,7 @@ def op(NLSAPar, DD, posPath, posPsi1, imgAll, msk2, CTF, ExtPar):  #pass the msk
         print('A is an imaginary matrix!')
         sys.exit()
 
-    U, S, V = svdRF.op(A)
+    U, S, V = svdRF(A)
     VX = np.matmul(V.T, psiC.T)
 
     sdiag = np.diag(S)
@@ -132,7 +133,7 @@ def op(NLSAPar, DD, posPath, posPsi1, imgAll, msk2, CTF, ExtPar):  #pass the msk
         IMGT[:, i] = ttmp
 
     nSrecon = min(IMGT.shape)
-    Drecon = L2_distance.op(IMGT, IMGT)
+    Drecon = L2_distance(IMGT, IMGT)
     k = nSrecon
 
     lamb, psirec, sigma, mu, logEps, logSumWij, popt, R_squared = DMembeddingII.op((Drecon**2), k, tune, 30)

ManifoldEM/NLSAmovie.py

@@ -8,9 +8,9 @@
 
 from functools import partial
 from contextlib import contextmanager
-from subprocess import Popen
 
-from ManifoldEM import myio, p, set_params, myio, makeMovie
+from ManifoldEM import myio, p, set_params, myio
+from ManifoldEM.core import makeMovie
 '''
 % scriptPsiNLSAmovie
 % Matlab Version V1.2
@@ -65,7 +65,7 @@ def movie(input_data, out_dir, dist_file, psi2_file, fps):
         IMG1All.append(data['IMG1'])
         Topo_mean.append(data['Topo_mean'])
         # make movie
-        makeMovie.op(IMG1All[psinum], prD, psinum, fps)
+        makeMovie(IMG1All[psinum], prD, psinum, fps)
 
         ######################
         # write topos

ManifoldEM/core.py

@@ -0,0 +1,158 @@
+"""Utilities that can be defined with basically one function."""
+import os
+import imageio
+import warnings
+import numpy as np
+
+from scipy.optimize import curve_fit, OptimizeWarning
+
+from ManifoldEM import myio, p
+
+warnings.simplefilter(action='ignore', category=OptimizeWarning)
+
+
+def clusterAvg(clust, PrD):
+    dist_file = p.get_dist_file(PrD)
+    data = myio.fin1(dist_file)
+
+    imgAll = data['imgAll']
+    boxSize = np.shape(imgAll)[1]
+
+    imgAvg = np.zeros(shape=(boxSize, boxSize), dtype=float)
+
+    for i in clust:
+        imgAvg += imgAll[i]
+
+    return imgAvg
+
+
+def L2_distance(a, b):
+    """Computes the Euclidean distance matrix between a and b.
+    Inputs:
+        A: D x M array.
+        B: D x N array.
+    Returns:
+        E: M x N Euclidean distances between vectors in A and B.
+    Author   : Roland Bunschoten
+               University of Amsterdam
+               Intelligent Autonomous Systems (IAS) group
+               Kruislaan 403  1098 SJ Amsterdam
+               tel.(+31)20-5257524
+               [email protected]
+    Last Rev : Wed Oct 20 08:58:08 MET DST 1999
+    Tested   : PC Matlab v5.2 and Solaris Matlab v5.3
+    Copyright notice: You are free to modify, extend and distribute
+       this code granted that the author of the original code is
+       mentioned as the original author of the code.
+    Fixed by JBT (3/18/00) to work for 1-dimensional vectors
+    and to warn for imaginary numbers.  Also ensures that
+    output is all real, and allows the option of forcing diagonals to
+    be zero.
+    Basic functionality ported to Python 2.7 by JCS (9/21/2013).
+
+    Copyright (c) Columbia University Hstau Liao 2019
+    """
+    eps = 1e-8
+
+    if a.shape[0] != b.shape[0]:
+        raise ValueError("A and B should be of same dimensionality")
+
+    aa = np.sum(a**2, axis=0)
+    bb = np.sum(b**2, axis=0)
+    ab = np.matmul(a.T, b)
+    tmp = aa[:, np.newaxis] + bb[np.newaxis, :] - 2 * ab
+    tmp[tmp < eps] = 0
+    return np.sqrt(tmp)
+
+
+def svdRF(A):
+    # Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
+    # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    # Copyright (c) Columbia University Hstau Liao 2018 (python version)
+
+    def tidyUp(D, EV):
+        order = np.argsort(D)[::-1]
+        D = np.sort(D)[::-1]
+        EV = EV[:, order]
+        sqrtD = np.sqrt(D)
+        S = np.diag(sqrtD)
+        invS = np.diag(1. / sqrtD)
+        return (D, EV, S, invS)
+
+    D1, D2 = A.shape
+    if D1 > D2:
+        D, V = np.linalg.eigh(np.matmul(A.T, A))
+        D, V, S, invS = tidyUp(D, V)
+        U = np.matmul(A, np.matmul(V, invS))
+    else:
+        D, U = np.linalg.eigh(np.matmul(A, A.T))
+        D, U, S, invS = tidyUp(D, U)
+        V = np.matmul(A.T, np.matmul(U, invS))
+
+    return (U, S, V)
+
+
+def makeMovie(IMG1, prD, psinum, fps):
+    # Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
+    # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    # Copyright (c) Columbia University Hstau Liao 2018 (python version)
+    # Copyright (c) Columbia University Evan Seitz 2019 (python version)
+
+    dim = int(np.floor(np.sqrt(max(IMG1.shape))))  # window size
+    nframes = IMG1.shape[1]
+    images = -IMG1
+    gif_path = os.path.join(p.out_dir, "topos", f"PrD_{prD + 1}", f'psi_{psinum + 1}.gif')
+    frame_dt = 1.0/fps
+    with imageio.get_writer(gif_path, mode='I', duration=frame_dt) as writer:
+        for i in range(nframes):
+            img = images[:, i].reshape(dim, dim)
+            frame = np.round(255 * (img - np.min(img)) / (np.max(img) - np.min(img))).astype(np.uint8)
+
+            frame_path = p.out_dir + '/topos/PrD_{}/psi_{}/frame{:02d}.png'.format(prD + 1, psinum + 1, i)
+            imageio.imwrite(frame_path, frame)
+            writer.append_data(frame)
+
+
+def fergusonE(D, logEps, a0):
+    # Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
+    # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    # Copyright (c) Columbia University Hstau Liao 2018 (python version)
+    # Copyright (c) Columbia University Evan Seitz 2019 (python version)
+
+    def fun(xx, aa0, aa1, aa2, aa3):
+        return aa3 + aa2 * np.tanh(aa0 * xx + aa1)
+
+    def find_thres(logEps, D2):
+        eps = np.exp(logEps)
+        d = 1. / (2. * np.max(eps)) * D2
+        sg = np.sort(d)
+        ss = np.sum(np.exp(-sg))
+        thr = max(-np.log(0.01 * ss / len(D2)), 10)  # taking 1% of the average (10)
+        return thr
+
+    # range of values to try:
+    logSumWij = np.zeros(len(logEps))
+    D2 = D * D
+    thr = find_thres(logEps, D2)
+    for k in range(len(logEps)):
+        eps = np.exp(logEps[k])
+        d = 1. / (2. * eps) * D2
+        d = -d[d < thr]
+        Wij = np.exp(d)  #see Coifman 2008
+        logSumWij[k] = np.log(np.sum(Wij))
+
+    # curve fitting of a tanh():
+    resnorm = np.inf
+    cc = 0
+    while (resnorm > 100):
+        cc += 1
+        popt, pcov = curve_fit(fun, logEps, logSumWij, p0=a0)
+        resnorm = np.sum(np.sqrt(np.fabs(np.diag(pcov))))
+        a0 = 1 * (np.random.rand(4, 1) - .5)
+
+        residuals = logSumWij - fun(logEps, popt[0], popt[1], popt[2], popt[3])
+        ss_res = np.sum(residuals**2)  #residual sum of squares
+        ss_tot = np.sum((logSumWij - np.mean(logSumWij))**2)  #total sum of squares
+        R_squared = 1 - (ss_res / ss_tot)  #R**2-value
+
+    return (popt, logSumWij, resnorm, R_squared)