diff --git a/ManifoldEM/CC/runGlobalOptimization.py b/ManifoldEM/CC/runGlobalOptimization.py index 09131df..382856b 100644 --- a/ManifoldEM/CC/runGlobalOptimization.py +++ b/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 diff --git a/ManifoldEM/DMembeddingII.py b/ManifoldEM/DMembeddingII.py index 963b6ab..2eb6d39 100644 --- a/ManifoldEM/DMembeddingII.py +++ b/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 diff --git a/ManifoldEM/L2_distance.py b/ManifoldEM/L2_distance.py deleted file mode 100644 index ff85f90..0000000 --- a/ManifoldEM/L2_distance.py +++ /dev/null @@ -1,42 +0,0 @@ -import numpy as np - -eps = 1e-8 -''' -Copyright (c) Columbia University Hstau Liao 2019 -''' - - -def op(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 - bunschot@wins.uva.nl - 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). - """ - - 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) diff --git a/ManifoldEM/NLSA.py b/ManifoldEM/NLSA.py index 8cbde1b..f42b41a 100644 --- a/ManifoldEM/NLSA.py +++ b/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) diff --git a/ManifoldEM/NLSAmovie.py b/ManifoldEM/NLSAmovie.py index 90270af..8aeb62e 100644 --- a/ManifoldEM/NLSAmovie.py +++ b/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 diff --git a/ManifoldEM/clusterAvg.py b/ManifoldEM/clusterAvg.py deleted file mode 100644 index 0cd4002..0000000 --- a/ManifoldEM/clusterAvg.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np - -from ManifoldEM import myio, p -''' -Copyright (c) Columbia University Evan Seitz 2019 -''' - - -def op(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 diff --git a/ManifoldEM/core.py b/ManifoldEM/core.py new file mode 100644 index 0000000..b45229e --- /dev/null +++ b/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 + bunschot@wins.uva.nl + 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) diff --git a/ManifoldEM/fergusonE.py b/ManifoldEM/fergusonE.py deleted file mode 100644 index 4b7fbe9..0000000 --- a/ManifoldEM/fergusonE.py +++ /dev/null @@ -1,63 +0,0 @@ -import numpy as np - -from scipy.optimize import curve_fit, OptimizeWarning - -import warnings - -warnings.simplefilter(action='ignore', category=OptimizeWarning) -""" -%-------------------------------------------------------------------------- -% function ferguson(D,s) -% D: Distance matrix -% logEps: Range of values to try -% Adapted from Chuck, 2011 -%-------------------------------------------------------------------------- -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): - F = aa3 + aa2 * np.tanh(aa0 * xx + aa1) - return F - - -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 - - -def op(D, logEps, a0): - # 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) diff --git a/ManifoldEM/makeMovie.py b/ManifoldEM/makeMovie.py deleted file mode 100644 index da25c31..0000000 --- a/ManifoldEM/makeMovie.py +++ /dev/null @@ -1,30 +0,0 @@ -import time - -from ManifoldEM import p -''' -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 op(IMG1, prD, psinum, fps): - import os - import imageio - import numpy as np - - 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) diff --git a/ManifoldEM/svdRF.py b/ManifoldEM/svdRF.py deleted file mode 100644 index da76dae..0000000 --- a/ManifoldEM/svdRF.py +++ /dev/null @@ -1,30 +0,0 @@ -import numpy as np -''' -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) - - -def op(A): - 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) diff --git a/notes/TODO.md b/notes/TODO.md index d9501a1..3754811 100644 --- a/notes/TODO.md +++ b/notes/TODO.md @@ -23,3 +23,4 @@ * TODO proper relative pathing * Shouldn't have to cache all these paths either... * TODO Resume after eigenvector selection broken +* TODO Parallelize 3d trajectory calculation diff --git a/scripts/ManifoldEM_GUI b/scripts/ManifoldEM_GUI index 3358446..5c25ee9 100755 --- a/scripts/ManifoldEM_GUI +++ b/scripts/ManifoldEM_GUI @@ -51,8 +51,9 @@ import pickle from ManifoldEM import p p.user_dir = output_dir #default -from ManifoldEM import (p, myio, Data, GetDistancesS2, manifoldAnalysis, psiAnalysis, NLSAmovie, embedd, clusterAvg, +from ManifoldEM import (p, myio, Data, GetDistancesS2, manifoldAnalysis, psiAnalysis, NLSAmovie, embedd, FindConformationalCoord, EL1D, backup, PrepareOutputS2, set_params) +from ManifoldEM.core import clusterAvg import logging from traits.api import HasTraits, Instance, on_trait_change, List, Str, Int, Float, Range, Button, Enum @@ -6762,7 +6763,7 @@ class Manifold2dCanvas(QtGui.QDialog): print('Encircled Points:', index_enc) - Manifold2dCanvas.imgAvg = clusterAvg.op(idx_encircled, P4.user_PrD - 1) + Manifold2dCanvas.imgAvg = clusterAvg(idx_encircled, P4.user_PrD - 1) self.ClusterAvg() def ClusterAvg(self): diff --git a/scripts/mrc2svd.py b/scripts/mrc2svd.py index ce62d91..296f9e9 100755 --- a/scripts/mrc2svd.py +++ b/scripts/mrc2svd.py @@ -40,7 +40,7 @@ def op(svd_dir, proj_name, user_dir): print('Input volume:', (bin + 1)) topoNum = 8 #number of topos considered print('Performing SVD...') - U, S, V = svdRF.op(b) + U, S, V = svdRF(b) print('SVD complete. Preparing volumes...') sdiag = np.diag(S) @@ -103,6 +103,6 @@ def op(svd_dir, proj_name, user_dir): sys.path.append(modDir) import p import mrcfile - import svdRF + from ManifoldEM.core import svdRF import set_params op(os.path.splitext(sys.argv[0])[0], sys.argv[1], sys.argv[2]) #enter the params file name