-
Notifications
You must be signed in to change notification settings - Fork 2
/
test_mref_gpu_align.py
1464 lines (1248 loc) · 67.8 KB
/
test_mref_gpu_align.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/local/EMAN2/bin/python
#******************************************************************************
#*
#* GPU based multireference alignment using CuPy
#*
#* Author (C) Szu-Chi Chung 2020 ([email protected])
#* Cheng-Yu Hung 2020 ([email protected])
#* Huei-Lun Siao 2020 ([email protected])
#* Hung-Yi Wu 2020 ([email protected])
#* Markus Stabrin 2019 ([email protected])
#* Fabian Schoenfeld 2019 ([email protected])
#* Thorsten Wagner 2019 ([email protected])
#* Tapu Shaikh 2019 ([email protected])
#* Adnan Ali 2019 ([email protected])
#* Luca Lusnig 2019 ([email protected])
#* Toshio Moriya 2019 ([email protected])
#* Pawel A.Penczek, 09/09/2006 ([email protected])
#*
#* Copyright (C) 2020 SABID Laboratory, Institute of Statistical Science, Academia Sinica
#* Copyright (c) 2019 Max Planck Institute of Molecular Physiology
#* Copyright (c) 2000-2006 The University of Texas - Houston Medical School
#*
#* This program is free software: you can redistribute it and/or modify it
#* under the terms of the GNU General Public License as published by the Free
#* Software Foundation, either version 3 of the License, or (at your option) any
#* later version.
#*
#* This program is distributed in the hope that it will be useful, but WITHOUT
#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License along with
#* this program. If not, please visit: http://www.gnu.org/licenses/
#*
#******************************************************************************/
from __future__ import print_function
import os
import global_def
from global_def import *
from optparse import OptionParser
import sys
import sparx as spx
###############################################################
def kernel_sum_oe(gpu_ptr, nima, nx, gpu_refs, numref):
class CustomGpuArray(object):
def __init__(self, shape, dtype):
if (shape.__class__ is not (1,).__class__):
shape = (shape,)
self.shape = shape
self.ptr = gpu_ptr# <-- a RAW GPU pointer
# def __del__(self):
# cp.cuda.runtime.free(self.ptr) # <-- malloc/free should be a pair
@property
def __cuda_array_interface__(self):
desc = {
'shape': self.shape,
'version': 0,
'data': (self.ptr, False),
'typestr': '<f4', # <-- because of float32, set '<f4',
'descr': [('', '<f4')]
}
return desc
obj = CustomGpuArray((nima,nx,nx), cp.float32)
gpu_batch_img = cp.asarray(obj, cp.float32)
gpu_ref_ave = cp.zeros((numref, 2, nx, nx), cp.float32)
for ref_id in range(numref):
even_indices = (cp.where(gpu_refs[::2]==ref_id)[0])*2
gpu_ref_ave[ref_id][0] = cp.sum(gpu_batch_img[even_indices], axis=0)
odd_indices = (cp.where(gpu_refs[1::2]==ref_id)[0])*2+1
gpu_ref_ave[ref_id][1] = cp.sum(gpu_batch_img[odd_indices], axis=0)
obj = None
return gpu_ref_ave
########################################################### GPU
import cupy as cp
import math
import numpy as np
import ctypes
from EMAN2 import EMNumPy
from cupy.cuda.nvtx import RangePush, RangePushC, RangePop
CUDA_PATH = os.path.join(os.path.dirname(__file__), "cuda", "")
cu_module = ctypes.CDLL( CUDA_PATH + "gpu_aln_pack.so" )
float_ptr = ctypes.POINTER(ctypes.c_float)
cu_module.ref_free_alignment_2D_init.restype = ctypes.c_ulonglong
cu_module.pre_align_init.restype = ctypes.c_ulonglong
cu_module.mref_align_run.restype = ctypes.c_ulonglong
# cu_module.get_gpu_img_ptr.restype = ctypes.c_ulonglong
# cu_module.get_gpu_img_ptr.restype = ctypes.POINTER(ctypes.c_float)
# cu_module.get_gpu_img_ptr.restype = ctypes.POINTER(ctypes.c_void_p)
class Freezeable( object ):
def freeze( self ):
self._frozen = None
def __setattr__( self, name, value ):
if hasattr( self, '_frozen' )and not hasattr( self, name ):
raise AttributeError( "Error! Trying to add new attribute '%s' to frozen class '%s'!" % (name,self.__class__.__name__) )
object.__setattr__( self, name, value )
class AlignConfig( ctypes.Structure, Freezeable ):
_fields_ = [ # data param
("sbj_num", ctypes.c_uint), # number of subject images we want to align
("ref_num", ctypes.c_uint), # number of reference images we want to align the subjects to
("img_dim", ctypes.c_uint), # image dimension (in both x- and y-direction)
# polar sampling parameters
("ring_num", ctypes.c_uint), # number of rings when converting images to polar coordinates
("ring_len", ctypes.c_uint), # number of rings when converting images to polar coordinates
# shift parameters
("shift_step", ctypes.c_float), # step range when applying translational shifts
("shift_rng_x", ctypes.c_float), # translational shift range in x-direction
("shift_rng_y", ctypes.c_float) ] # translational shift range in y-direction
class AlignParam( ctypes.Structure, Freezeable ):
_fields_ = [ ("sbj_id", ctypes.c_int),
("ref_id", ctypes.c_int),
("shift_x", ctypes.c_float),
("shift_y", ctypes.c_float),
("angle", ctypes.c_float),
("mirror", ctypes.c_bool) ]
def __str__(self):
return "s_%d/r_%d::(%d,%d;%.2f)" % (self.sbj_id, self.ref_id, self.shift_x, self.shift_y, self.angle) \
+("[M]" if self.mirror else "")
aln_param_ptr = ctypes.POINTER(AlignParam)
def get_c_ptr_array( emdata_list ):
ptr_list = []
for img in emdata_list:
img_np = EMNumPy.em2numpy( img )
assert img_np.flags['C_CONTIGUOUS'] == True
assert img_np.dtype == np.float32
img_ptr = img_np.ctypes.data_as(float_ptr)
ptr_list.append(img_ptr)
return (float_ptr*len(emdata_list))(*ptr_list)
def print_gpu_info( cuda_device_id ):
cu_module.print_gpu_info( ctypes.c_int(cuda_device_id) )
def send_EMData_w_ID(img, dst, tag, comm=-1):
from mpi import mpi_send, MPI_INT, MPI_FLOAT, MPI_COMM_WORLD
from sp_utilities import get_image_data
if comm==-1:
comm = MPI_COMM_WORLD
img_head = []
img_head.append(img.get_xsize())
img_head.append(img.get_ysize())
img_head.append(img.get_zsize())
img_head.append(img.is_complex())
img_head.append(img.is_ri())
img_head.append(img.get_attr("changecount"))
img_head.append(img.is_complex_x())
img_head.append(img.get_attr("is_complex_ri"))
img_head.append(int(img.get_attr("apix_x")*10000))
img_head.append(int(img.get_attr("apix_y")*10000))
img_head.append(int(img.get_attr("apix_z")*10000))
img_head.append(img.get_attr("ID"))
head_tag = 2*tag
mpi_send(img_head, 12, MPI_INT, dst, head_tag, comm)
img_data = get_image_data(img)
data_tag = 2*tag + 1
ntot = img_head[0]*img_head[1]*img_head[2]
mpi_send(img_data, ntot, MPI_FLOAT, dst, data_tag, comm)
def recv_EMData_w_ID(src, tag, comm=-1):
from mpi import mpi_recv, MPI_INT, MPI_FLOAT, MPI_COMM_WORLD
from numpy import reshape
from EMAN2 import EMNumPy
if comm==-1:
comm = MPI_COMM_WORLD
head_tag = 2*tag
img_head = mpi_recv(12, MPI_INT, src, head_tag, comm)
nx = int(img_head[0])
ny = int(img_head[1])
nz = int(img_head[2])
is_complex = int(img_head[3])
is_ri = int(img_head[4])
data_tag = 2*tag+1
ntot = nx*ny*nz
img_data = mpi_recv(ntot, MPI_FLOAT, src, data_tag, comm)
if nz != 1:
img_data = reshape(img_data, (nz, ny, nx))
elif ny != 1:
img_data = reshape(img_data, (ny, nx))
else:
pass
img = EMNumPy.numpy2em(img_data)
img.set_complex(is_complex)
img.set_ri(is_ri)
img.set_attr_dict({"changecount":int(img_head[5]),
"is_complex_x":int(img_head[6]),
"is_complex_ri":int(img_head[7]),
"apix_x":int(img_head[8])/10000.0,
"apix_y":int(img_head[9])/10000.0,
"apix_z":int(img_head[10])/10000.0,
"ID":int(img_head[11])})
return img
########################################################### GPU
###############################################################
def mref_ali2d_gpu(
filename,stack, refim, outdir, maskfile=None, ir=1, ou=-1, rs=1, xrng = 0, yrng = 0, step = 1,
center=-1, maxit=0, CTF=False, snr=1.0,
user_func_name="ref_ali2d", rand_seed= 1000, number_of_proc = 1, myid = 0, main_node = 0, mpi_comm = None,
mpi_gpu_proc=False, gpu_class_limit=0, cuda_device_occ=0.9):
from sp_utilities import model_circle, combine_params2, inverse_transform2, drop_image, get_image, get_im
from sp_utilities import reduce_EMData_to_root, bcast_EMData_to_all, bcast_number_to_all
from sp_utilities import send_attr_dict
from sp_utilities import center_2D
from sp_statistics import fsc_mask
from sp_alignment import Numrinit, ringwe, search_range
from sp_fundamentals import rot_shift2D, fshift
from sp_utilities import get_params2D, set_params2D
from random import seed, randint
from sp_morphology import ctf_2
from sp_filter import filt_btwl, filt_params
from numpy import reshape, shape
from sp_utilities import print_msg, print_begin_msg, print_end_msg
import os
import alignment
import fundamentals
import sys
import mpi
import time
import utilities as util
from sp_applications import MPI_start_end
from mpi import mpi_comm_size, mpi_comm_rank, MPI_COMM_WORLD
from mpi import mpi_reduce, mpi_bcast, mpi_barrier, mpi_recv, mpi_send
from mpi import MPI_SUM, MPI_FLOAT, MPI_INT
import sp_user_functions
from sp_applications import MPI_start_end
user_func = sp_user_functions.factory[user_func_name]
# sanity check: gpu procs sound off
if not mpi_gpu_proc: return []
#----------------------------------[ setup ]
print( time.strftime("%Y-%m-%d %H:%M:%S :: ", time.localtime()) + "mref_ali2d_gpu() :: MPI proc["+str(myid)+"] run pre-alignment CPU setup" )
# mpi communicator sanity check
assert mpi_comm != None
xrng = [xrng]
yrng = [yrng]
step = [step]
# polar conversion (ring) parameters
first_ring = int(ir)
last_ring = int(ou)
rstep = int(rs)
max_iter = int(maxit)
if max_iter ==0:
max_iter = 10
auto_stop = True
else:
auto_stop = False
# determine global index values of particles handled by this process
data = stack
total_nima = len(data)
total_nima = mpi.mpi_reduce(total_nima, 1, mpi.MPI_INT, mpi.MPI_SUM, 0, mpi_comm)
total_nima = mpi.mpi_bcast (total_nima, 1, mpi.MPI_INT, main_node, mpi_comm)[0]
list_of_particles = list(range(total_nima))
image_start, image_end = MPI_start_end(total_nima, number_of_proc, myid)
list_of_particles = list_of_particles[image_start:image_end] # list of global indices
nima = len(list_of_particles)
assert( nima == len(data) ) # sanity check
# read nx and broadcast to all nodes
# NOTE: nx is image size and images are assumed to be square
if myid == main_node:
nx = data[0].get_xsize()
else:
nx = 0
nx = util.bcast_number_to_all(nx, source_node=main_node, mpi_comm=mpi_comm)
if CTF:
phase_flip = True
else:
phase_flip = False
CTF = False # okay..?
# set default value for the last ring if none given
if last_ring == -1: last_ring = nx/2-2
# sanity check for last_ring value
if last_ring + max([max(xrng),max(yrng)]) > (nx-1) // 2:
ERROR( "Shift or radius is too large - particle crosses image boundary", "ali2d_MPI", 1 )
if maskfile:
import types
if type(maskfile) is bytes: mask = get_image(maskfile)
else: mask = maskfile
else : mask = model_circle(last_ring, nx, nx)
# references, do them on all processors
refi = []
numref = EMUtil.get_image_count(refim)
# image center
cny = cnx = nx/2+1
mode = "F"
RangePush("Preprocess data")
# prepare reference images on all nodes
ima = data[0].copy()
ima.to_zero()
for j in range(numref):
# even, odd, numer of even, number of images. After frc, totav
refi.append([get_im(refim,j), ima.copy(), 0])
refi[j][0].process_inplace("normalize.mask", {"mask":mask, "no_sigma":1}) # normalize reference images to N(0,1)
# create/reset alignment parameters
for im in range(nima):
#data[im].set_attr( 'ID', list_of_particles[im] )
util.set_params2D( data[im], [0.0, 0.0, 0.0, 0, 1.0], 'xform.align2d' )
data[im].process_inplace("normalize.mask", {"mask":mask, "no_sigma":0}) # subtract average under the mask
#data[im] -= st[0]
if phase_flip:
data[im] = filt_ctf(data[im], data[im].get_attr("ctf"), binary = True)
# precalculate rings & weights
numr = alignment.Numrinit( first_ring, last_ring, rstep, mode )
wr = alignment.ringwe( numr, mode )
RangePop()
if myid == main_node:
seed(rand_seed)
a0 = -1.0
ref_data = [mask, center, None, None]
again = True
Iter = 0
#----------------------------------[ gpu setup ]
print( time.strftime("%Y-%m-%d %H:%M:%S :: ", time.localtime()) + "mref_ali2d_gpu() :: MPI proc["+str(myid)+"] run pre-alignment GPU setup" )
# alignment parameters
aln_cfg = AlignConfig(
len(data), numref, # no. of images & averages
data[0].get_xsize(), # image size (images are always assumed to be square)
numr[-3], 256, # polar sampling params (no. of rings, no. of sample on ring)
step[0], xrng[0], xrng[0]) # shift params (step size & shift range in x/y dim)
aln_cfg.freeze()
# find largest batch size we can fit on the given card
gpu_batch_limit = 0
RangePush("Determine batch size")
for split in [ 2**i for i in range(int(math.log(len(data),2))+1) ][::-1]:
aln_cfg.sbj_num = min( gpu_batch_limit + split, len(data) )
if cu_module.pre_align_size_check( ctypes.c_uint(len(data)), ctypes.byref(aln_cfg), ctypes.c_uint(myid), ctypes.c_float(cuda_device_occ), ctypes.c_bool(False) ) == True:
gpu_batch_limit += split
gpu_batch_limit = aln_cfg.sbj_num
RangePop()
print( time.strftime("%Y-%m-%d %H:%M:%S :: ", time.localtime()) + "mref_ali2d_gpu() :: GPU["+str(myid)+"] pre-alignment batch size: %d/%d" % (gpu_batch_limit, len(data)) )
# initialize gpu resources (returns location for our alignment parameters in CUDA unified memory)
gpu_aln_param = cu_module.pre_align_init( ctypes.c_uint(len(data)), ctypes.byref(aln_cfg), ctypes.c_uint(myid) )
gpu_aln_param = ctypes.cast( gpu_aln_param, aln_param_ptr ) # cast to allow Python-side r/w access
gpu_batch_count = len(data)/gpu_batch_limit if len(data)%gpu_batch_limit==0 else len(data)//gpu_batch_limit+1
print( time.strftime("%Y-%m-%d %H:%M:%S :: ", time.localtime()) + "mref_ali2d_gpu() :: GPU["+str(myid)+"] batch count:", gpu_batch_count )
# if the local stack fits on the gpu we only fetch the img data and the reference data once before we loop
if( gpu_batch_count == 1 ):
cu_module.pre_align_fetch(
get_c_ptr_array(data),
ctypes.c_uint(len(data)),
ctypes.c_char_p("sbj_batch") )
#----------------------------------[ alignment ]
print( time.strftime("%Y-%m-%d %H:%M:%S :: ", time.localtime()) + "mref_ali2d_gpu() :: MPI proc["+str(myid)+"] start alignment iterations" )
N_step = 0
# set gpu search parameters
cu_module.reset_shifts( ctypes.c_float(xrng[N_step]), ctypes.c_float(step[N_step]) )
while Iter < max_iter and again:
cu_module.pre_align_fetch(
get_c_ptr_array([tmp[0] for tmp in refi]),
ctypes.c_int(numref),
ctypes.c_char_p("ref_batch"))
# backup last iteration's alignment parameters
RangePush("Iteration's alignment parameters")
old_ali_params = []
for im in range(nima):
alpha, sx, sy, mirror, scale = util.get_params2D(data[im])
old_ali_params.extend([alpha, sx, sy, mirror])
RangePop()
############################################GPU
#
# FOR gpu_batch_i DO ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ >
gpu_ref_sum = cp.zeros((numref, 2, nx, nx), cp.float32)
for gpu_batch_idx in range(gpu_batch_count):
# determine next gpu batch and shove to the gpu (NOTE: if we
# only have a single batch, this already happened earlier)
RangePush("Batch fetch data")
if( gpu_batch_count > 1 ):
gpu_batch_start = gpu_batch_idx * gpu_batch_limit
gpu_batch_end = gpu_batch_start + gpu_batch_limit
if gpu_batch_end > len(data): gpu_batch_end = len(data)
cu_module.pre_align_fetch(
get_c_ptr_array( data[gpu_batch_start:gpu_batch_end] ),
ctypes.c_int( gpu_batch_end-gpu_batch_start ),
ctypes.c_char_p("sbj_batch") )
else:
gpu_batch_start = 0
gpu_batch_end = len(data)
RangePop()
# run the alignment on gpu
# Get ptr of alligned images
RangePush("Batch alignment")
gpu_img_ptr = cu_module.mref_align_run( ctypes.c_int(gpu_batch_start), ctypes.c_int(gpu_batch_end))
RangePop()
gpu_refs = cp.asarray([gpu_aln_param[k].ref_id for k in range(gpu_batch_start,gpu_batch_end)])
RangePush("Batch sum oe fast cp.sum")
gpu_ref_sum += kernel_sum_oe(gpu_img_ptr, gpu_batch_end-gpu_batch_start, nx, gpu_refs, numref)
RangePop()
# print progress bar
gpu_calls_ttl = 1 * max_iter * gpu_batch_count
#gpu_calls_ttl = len(xrng) * max_iter * gpu_batch_count - 1
gpu_calls_cnt = N_step*max_iter*gpu_batch_count + Iter*gpu_batch_count + gpu_batch_idx
gpu_calls_prc = int( float(gpu_calls_cnt+1)/gpu_calls_ttl * 50.0 )
sys.stdout.write( "\r[MREF-ALIGN][GPU"+str(myid)+"][" + "="*gpu_calls_prc + "-"*(50-gpu_calls_prc) + "]~[%d/%d]~[%.2f%%]\n" % (gpu_calls_cnt+1, gpu_calls_ttl, (float(gpu_calls_cnt+1)/gpu_calls_ttl)*100.0) )
sys.stdout.flush()
if gpu_calls_cnt+1 == gpu_calls_ttl: print("")
# < ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ FOR gpu_batch_i DO
#set reference image to zero
RangePush("transfer angle and average")
assign = [[] for i in range(numref)]
# ------------------------------------------------------------
for j in range(numref):
refi[j][0].to_zero()
refi[j][1].to_zero()
refi[j][2] = 0
for k, img in enumerate(data):
iref = int(gpu_aln_param[k].ref_id)
img.set_attr('assign',iref)
assign[iref].append(list_of_particles[k])
refi[iref][2] += 1.0
ref_sum = cp.asnumpy(gpu_ref_sum)
for iref in range(numref):
for it in [0,1]:
Util.add_img(refi[iref][it], EMNumPy.numpy2em(ref_sum[iref][it]))
RangePop()
# ------------------------------------------------------------
# print("="*30)
# print("gpu: 10,5,75:85 ", gpu_batch_img[10,5,75:85])
# print("cpu: 10,5,75:85 ", data[10].numpy()[5,75:85])
# print("="*30)
#===========
RangePush("MPI average")
for j in range(numref):
reduce_EMData_to_root(refi[j][0],myid,main_node,comm = mpi_comm)
reduce_EMData_to_root(refi[j][1],myid,main_node,comm = mpi_comm)
refi[j][2] = mpi_reduce(refi[j][2], 1, MPI_FLOAT, MPI_SUM, main_node, mpi_comm)
if(myid == main_node): refi[j][2] = int(refi[j][2][0])
RangePop()
RangePush("set param")
#gather assignments
for j in range(numref):
if myid == main_node:
for n in range(number_of_proc):
if n != main_node:
import sp_global_def
ln = mpi_recv(1, MPI_INT, n, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, mpi_comm)
lis = mpi_recv(ln[0], MPI_INT, n, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, mpi_comm)
for l in range(ln[0]): assign[j].append(int(lis[l]))
else:
import sp_global_def
mpi_send(len(assign[j]), 1, MPI_INT, main_node, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, mpi_comm)
mpi_send(assign[j], len(assign[j]), MPI_INT, main_node, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, mpi_comm)
if myid == main_node:
# replace the name of the stack with reference with the current one
refim = os.path.join(outdir,"aqm%03d.hdf"%Iter)
a1 = 0.0
ave_fsc = []
for j in range(numref):
if refi[j][2] < 4:
#ERROR("One of the references vanished","mref_ali2d_MPI",1)
# if vanished, put a random image (only from main node!) there
assign[j] = []
assign[j].append(randint(0,nima-1))
refi[j][0] = data[assign[j][0]].copy()
else:
#frsc = fsc_mask(refi[j][0], refi[j][1], mask, 1.0, os.path.join(outdir,"drm%03d%04d"%(Iter, j)))
from sp_statistics import fsc
frsc = fsc(refi[j][0], refi[j][1], 1.0, os.path.join(outdir,"drm%03d%04d.txt"%(Iter,j)))
Util.add_img( refi[j][0], refi[j][1] )
Util.mul_scalar( refi[j][0], 1.0/float(refi[j][2]) )
if ave_fsc == []:
for i in range(len(frsc[1])): ave_fsc.append(frsc[1][i])
c_fsc = 1
else:
for i in range(len(frsc[1])): ave_fsc[i] += frsc[1][i]
c_fsc += 1
#print 'OK', j, len(frsc[1]), frsc[1][0:5], ave_fsc[0:5]
#print 'sum', sum(ave_fsc)
if sum(ave_fsc) != 0:
for i in range(len(ave_fsc)):
ave_fsc[i] /= float(c_fsc)
frsc[1][i] = ave_fsc[i]
for j in range(numref):
ref_data[2] = refi[j][0]
ref_data[3] = frsc
refi[j][0], cs = user_func(ref_data)
# write the current average
TMP = []
for i_tmp in range(len(assign[j])): TMP.append(float(assign[j][i_tmp]))
TMP.sort()
refi[j][0].set_attr_dict({'ave_n': refi[j][2], 'members': TMP })
del TMP
refi[j][0].process_inplace("normalize.mask", {"mask":mask, "no_sigma":1})
refi[j][0].write_image(refim, j)
Iter += 1
msg = "ITERATION #%3d %d\n\n"%(Iter,again)
print_msg(msg)
for j in range(numref):
msg = " group #%3d number of particles = %7d\n"%(j, refi[j][2])
print_msg(msg)
Iter = bcast_number_to_all(Iter, main_node, mpi_comm) # need to tell all
if again:
for j in range(numref):
bcast_EMData_to_all(refi[j][0], myid, main_node, mpi_comm)
RangePop()
# clean up
for k, img in enumerate(data):
# this is usually done in ormq()
angle = gpu_aln_param[k].angle
sx_neg = -gpu_aln_param[k].shift_x
sy_neg = -gpu_aln_param[k].shift_y
c_ang = math.cos( math.radians(angle) )
s_ang = -math.sin( math.radians(angle) )
shift_x = sx_neg*c_ang - sy_neg*s_ang
shift_y = sx_neg*s_ang + sy_neg*c_ang
mn = gpu_aln_param[k].mirror
set_params2D( img, [angle, shift_x, shift_y, int(mn), 1.0], "xform.align2d" )
del assign
RangePush("disk")
# write out headers and STOP, under MPI writing has to be done sequentially (time-consumming)
mpi_barrier(mpi_comm)
# if CTF and data_had_ctf == 0:
# for im in range(nima): data[im].set_attr('ctf_applied', 0)
# par_str = ['xform.align2d', 'assign', 'ID']
# if myid == main_node:
# from sp_utilities import file_type
# if(file_type(filename) == "bdb"):
# from sp_utilities import recv_attr_dict_bdb
# recv_attr_dict_bdb(main_node, filename, data, par_str, image_start, image_end, number_of_proc)
# else:
# from sp_utilities import recv_attr_dict
# recv_attr_dict(main_node, filename, data, par_str, image_start, image_end, number_of_proc)
# else: send_attr_dict(main_node, data, par_str, image_start, image_end)
#
# # free gpu resources
# cu_module.gpu_clear()
# mpi.mpi_barrier(mpi_comm)
RangePop()
if myid == main_node:
print_end_msg("mref_ali2d_GPU")
def mref_ali2d_MPI(stack, refim, outdir, maskfile = None, ir=1, ou=-1, rs=1, xrng=0, yrng=0, step=1, center=1, maxit=10, CTF=False, snr=1.0, user_func_name="ref_ali2d", rand_seed=1000):
# 2D multi-reference alignment using rotational ccf in polar coordinates and quadratic interpolation
from sp_utilities import model_circle, combine_params2, inverse_transform2, drop_image, get_image, get_im
from sp_utilities import reduce_EMData_to_root, bcast_EMData_to_all, bcast_number_to_all
from sp_utilities import send_attr_dict
from sp_utilities import center_2D
from sp_statistics import fsc_mask
from sp_alignment import Numrinit, ringwe, search_range
from sp_fundamentals import rot_shift2D, fshift
from sp_utilities import get_params2D, set_params2D
from random import seed, randint
from sp_morphology import ctf_2
from sp_filter import filt_btwl, filt_params
from numpy import reshape, shape
from sp_utilities import print_msg, print_begin_msg, print_end_msg
import os
import sys
import shutil
from sp_applications import MPI_start_end
from mpi import mpi_comm_size, mpi_comm_rank, MPI_COMM_WORLD
from mpi import mpi_reduce, mpi_bcast, mpi_barrier, mpi_recv, mpi_send
from mpi import MPI_SUM, MPI_FLOAT, MPI_INT
number_of_proc = mpi_comm_size(MPI_COMM_WORLD)
myid = mpi_comm_rank(MPI_COMM_WORLD)
main_node = 0
# create the output directory, if it does not exist
if os.path.exists(outdir): ERROR('Output directory exists, please change the name and restart the program', "mref_ali2d_MPI ", 1, myid)
mpi_barrier(MPI_COMM_WORLD)
import sp_global_def
if myid == main_node:
os.mkdir(outdir)
sp_global_def.LOGFILE = os.path.join(outdir, sp_global_def.LOGFILE)
print_begin_msg("mref_ali2d_MPI")
nima = EMUtil.get_image_count(stack)
image_start, image_end = MPI_start_end(nima, number_of_proc, myid)
nima = EMUtil.get_image_count(stack)
ima = EMData()
ima.read_image(stack, image_start)
first_ring=int(ir); last_ring=int(ou); rstep=int(rs); max_iter=int(maxit)
if max_iter == 0:
max_iter = 10
auto_stop = True
else:
auto_stop = False
if myid == main_node:
print_msg("Input stack : %s\n"%(stack))
print_msg("Reference stack : %s\n"%(refim))
print_msg("Output directory : %s\n"%(outdir))
print_msg("Maskfile : %s\n"%(maskfile))
print_msg("Inner radius : %i\n"%(first_ring))
nx = ima.get_xsize()
# default value for the last ring
if last_ring == -1: last_ring=nx/2-2
if myid == main_node:
print_msg("Outer radius : %i\n"%(last_ring))
print_msg("Ring step : %i\n"%(rstep))
print_msg("X search range : %f\n"%(xrng))
print_msg("Y search range : %f\n"%(yrng))
print_msg("Translational step : %f\n"%(step))
print_msg("Center type : %i\n"%(center))
print_msg("Maximum iteration : %i\n"%(max_iter))
print_msg("CTF correction : %s\n"%(CTF))
print_msg("Signal-to-Noise Ratio : %f\n"%(snr))
print_msg("Random seed : %i\n\n"%(rand_seed))
print_msg("User function : %s\n"%(user_func_name))
import sp_user_functions
user_func = sp_user_functions.factory[user_func_name]
if maskfile:
import types
if type(maskfile) is bytes: mask = get_image(maskfile)
else: mask = maskfile
else : mask = model_circle(last_ring, nx, nx)
# references, do them on all processors...
refi = []
numref = EMUtil.get_image_count(refim)
# IMAGES ARE SQUARES! center is in SPIDER convention
cnx = nx/2+1
cny = cnx
mode = "F"
#precalculate rings
numr = Numrinit(first_ring, last_ring, rstep, mode)
wr = ringwe(numr, mode)
# prepare reference images on all nodes
ima.to_zero()
for j in range(numref):
# even, odd, number of even, number of images. After frc, totav
refi.append([get_im(refim,j), ima.copy(), 0])
# for each node read its share of data
data = EMData.read_images(stack, list(range(image_start, image_end)))
for im in range(image_start, image_end):
data[im-image_start].set_attr('ID', im)
if myid == main_node: seed(rand_seed)
a0 = -1.0
again = True
Iter = 0
ref_data = [mask, center, None, None]
while Iter < max_iter and again:
ringref = []
mashi = cnx-last_ring-2
for j in range(numref):
refi[j][0].process_inplace("normalize.mask", {"mask":mask, "no_sigma":1}) # normalize reference images to N(0,1)
cimage = Util.Polar2Dm(refi[j][0] , cnx, cny, numr, mode)
Util.Frngs(cimage, numr)
Util.Applyws(cimage, numr, wr)
ringref.append(cimage)
# zero refi
refi[j][0].to_zero()
refi[j][1].to_zero()
refi[j][2] = 0
assign = [[] for i in range(numref)]
# begin MPI section
for im in range(image_start, image_end):
alpha, sx, sy, mirror, scale = get_params2D(data[im-image_start])
# Why inverse? 07/11/2015 PAP
alphai, sxi, syi, scalei = inverse_transform2(alpha, sx, sy)
# normalize
data[im-image_start].process_inplace("normalize.mask", {"mask":mask, "no_sigma":0}) # subtract average under the mask
# If shifts are outside of the permissible range, reset them
if(abs(sxi)>mashi or abs(syi)>mashi):
sxi = 0.0
syi = 0.0
set_params2D(data[im-image_start],[0.0,0.0,0.0,0,1.0])
ny = nx
txrng = search_range(nx, last_ring, sxi, xrng, "mref_ali2d_MPI")
txrng = [txrng[1],txrng[0]]
tyrng = search_range(ny, last_ring, syi, yrng, "mref_ali2d_MPI")
tyrng = [tyrng[1],tyrng[0]]
# align current image to the reference
[angt, sxst, syst, mirrort, xiref, peakt] = Util.multiref_polar_ali_2d(data[im-image_start],
ringref, txrng, tyrng, step, mode, numr, cnx+sxi, cny+syi)
iref = int(xiref)
# combine parameters and set them to the header, ignore previous angle and mirror
[alphan, sxn, syn, mn] = combine_params2(0.0, -sxi, -syi, 0, angt, sxst, syst, (int)(mirrort))
set_params2D(data[im-image_start], [alphan, sxn, syn, int(mn), scale])
data[im-image_start].set_attr('assign',iref)
# apply current parameters and add to the average
temp = rot_shift2D(data[im-image_start], alphan, sxn, syn, mn)
it = im%2
Util.add_img( refi[iref][it], temp)
assign[iref].append(im)
#assign[im] = iref
refi[iref][2] += 1.0
del ringref
# end MPI section, bring partial things together, calculate new reference images, broadcast them back
for j in range(numref):
reduce_EMData_to_root(refi[j][0], myid, main_node)
reduce_EMData_to_root(refi[j][1], myid, main_node)
refi[j][2] = mpi_reduce(refi[j][2], 1, MPI_FLOAT, MPI_SUM, main_node, MPI_COMM_WORLD)
if(myid == main_node): refi[j][2] = int(refi[j][2][0])
# gather assignements
for j in range(numref):
if myid == main_node:
for n in range(number_of_proc):
if n != main_node:
import sp_global_def
ln = mpi_recv(1, MPI_INT, n, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, MPI_COMM_WORLD)
lis = mpi_recv(ln[0], MPI_INT, n, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, MPI_COMM_WORLD)
for l in range(ln[0]): assign[j].append(int(lis[l]))
else:
import sp_global_def
mpi_send(len(assign[j]), 1, MPI_INT, main_node, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, MPI_COMM_WORLD)
mpi_send(assign[j], len(assign[j]), MPI_INT, main_node, sp_global_def.SPARX_MPI_TAG_UNIVERSAL, MPI_COMM_WORLD)
if myid == main_node:
# replace the name of the stack with reference with the current one
refim = os.path.join(outdir,"aqm%03d.hdf"%Iter)
a1 = 0.0
ave_fsc = []
for j in range(numref):
if refi[j][2] < 4:
#ERROR("One of the references vanished","mref_ali2d_MPI",1)
# if vanished, put a random image (only from main node!) there
assign[j] = []
assign[j].append( randint(image_start, image_end-1) - image_start )
refi[j][0] = data[assign[j][0]].copy()
#print 'ERROR', j
else:
#frsc = fsc_mask(refi[j][0], refi[j][1], mask, 1.0, os.path.join(outdir,"drm%03d%04d"%(Iter, j)))
from sp_statistics import fsc
frsc = fsc(refi[j][0], refi[j][1], 1.0, os.path.join(outdir,"drm%03d%04d.txt"%(Iter,j)))
Util.add_img( refi[j][0], refi[j][1] )
Util.mul_scalar( refi[j][0], 1.0/float(refi[j][2]) )
if ave_fsc == []:
for i in range(len(frsc[1])): ave_fsc.append(frsc[1][i])
c_fsc = 1
else:
for i in range(len(frsc[1])): ave_fsc[i] += frsc[1][i]
c_fsc += 1
#print 'OK', j, len(frsc[1]), frsc[1][0:5], ave_fsc[0:5]
#print 'sum', sum(ave_fsc)
if sum(ave_fsc) != 0:
for i in range(len(ave_fsc)):
ave_fsc[i] /= float(c_fsc)
frsc[1][i] = ave_fsc[i]
for j in range(numref):
ref_data[2] = refi[j][0]
ref_data[3] = frsc
refi[j][0], cs = user_func(ref_data)
# write the current average
TMP = []
for i_tmp in range(len(assign[j])): TMP.append(float(assign[j][i_tmp]))
TMP.sort()
refi[j][0].set_attr_dict({'ave_n': refi[j][2], 'members': TMP })
del TMP
refi[j][0].process_inplace("normalize.mask", {"mask":mask, "no_sigma":1})
refi[j][0].write_image(refim, j)
Iter += 1
msg = "ITERATION #%3d %d\n\n"%(Iter,again)
print_msg(msg)
for j in range(numref):
msg = " group #%3d number of particles = %7d\n"%(j, refi[j][2])
print_msg(msg)
Iter = bcast_number_to_all(Iter, main_node) # need to tell all
if again:
for j in range(numref):
bcast_EMData_to_all(refi[j][0], myid, main_node)
# clean up
del assign
# write out headers and STOP, under MPI writing has to be done sequentially (time-consumming)
mpi_barrier(MPI_COMM_WORLD)
if CTF and data_had_ctf == 0:
for im in range(len(data)): data[im].set_attr('ctf_applied', 0)
par_str = ['xform.align2d', 'assign', 'ID']
if myid == main_node:
from sp_utilities import file_type
if(file_type(stack) == "bdb"):
from sp_utilities import recv_attr_dict_bdb
recv_attr_dict_bdb(main_node, stack, data, par_str, image_start, image_end, number_of_proc)
else:
from sp_utilities import recv_attr_dict
recv_attr_dict(main_node, stack, data, par_str, image_start, image_end, number_of_proc)
else: send_attr_dict(main_node, data, par_str, image_start, image_end)
if myid == main_node:
print_end_msg("mref_ali2d_GPU")
def mref_ali2d(stack, refim, outdir, maskfile=None, ir=1, ou=-1, rs=1, xrng=0, yrng=0, step=1, center=1, maxit=0, CTF=False, snr=1.0, user_func_name="ref_ali2d", rand_seed=1000, MPI=False):
"""
Name
mref_ali2d - Perform 2-D multi-reference alignment of an image series
Input
stack: set of 2-D images in a stack file, images have to be squares
refim: set of initial reference 2-D images in a stack file
maskfile: optional maskfile to be used in the alignment
inner_radius: inner radius for rotational correlation > 0
outer_radius: outer radius for rotational correlation < nx/2-1
ring_step: step between rings in rotational correlation >0
x_range: range for translation search in x direction, search is +/xr
y_range: range for translation search in y direction, search is +/yr
translation_step: step of translation search in both directions
center: center the average
max_iter: maximum number of iterations the program will perform
CTF: if this flag is set, the program will use CTF information provided in file headers
snr: signal-to-noise ratio of the data
rand_seed: the seed used for generating random numbers
MPI: whether to use MPI version
Output
output_directory: directory name into which the output files will be written.
header: the alignment parameters are stored in the headers of input files as 'xform.align2d'.
"""
# 2D multi-reference alignment using rotational ccf in polar coordinates and quadratic interpolation
if MPI:
mref_ali2d_MPI(stack, refim, outdir, maskfile, ir, ou, rs, xrng, yrng, step, center, maxit, CTF, snr, user_func_name, rand_seed)
return
from sp_utilities import model_circle, combine_params2, inverse_transform2, drop_image, get_image
from sp_utilities import center_2D, get_im, get_params2D, set_params2D
from sp_statistics import fsc
from sp_alignment import Numrinit, ringwe, fine_2D_refinement, search_range
from sp_fundamentals import rot_shift2D, fshift
from random import seed, randint
import os
import sys
from sp_utilities import print_begin_msg, print_end_msg, print_msg
import shutil
# create the output directory, if it does not exist
if os.path.exists(outdir): shutil.rmtree(outdir) #ERROR('Output directory exists, please change the name and restart the program', "mref_ali2d", 1)
os.mkdir(outdir)
import sp_global_def
sp_global_def.LOGFILE = os.path.join(outdir, sp_global_def.LOGFILE)
first_ring=int(ir); last_ring=int(ou); rstep=int(rs); max_iter=int(maxit)
print_begin_msg("mref_ali2d")
print_msg("Input stack : %s\n"%(stack))
print_msg("Reference stack : %s\n"%(refim))
print_msg("Output directory : %s\n"%(outdir))
print_msg("Maskfile : %s\n"%(maskfile))
print_msg("Inner radius : %i\n"%(first_ring))
ima = EMData()
ima.read_image(stack, 0)
nx = ima.get_xsize()
# default value for the last ring
if last_ring == -1: last_ring = nx/2-2
print_msg("Outer radius : %i\n"%(last_ring))
print_msg("Ring step : %i\n"%(rstep))
print_msg("X search range : %i\n"%(xrng))
print_msg("Y search range : %i\n"%(yrng))
print_msg("Translational step : %i\n"%(step))
print_msg("Center type : %i\n"%(center))
print_msg("Maximum iteration : %i\n"%(max_iter))
print_msg("CTF correction : %s\n"%(CTF))
print_msg("Signal-to-Noise Ratio : %f\n"%(snr))
print_msg("Random seed : %i\n\n"%(rand_seed))
print_msg("User function : %s\n"%(user_func_name))
output = sys.stdout
import sp_user_functions
user_func = sp_user_functions.factory[user_func_name]
if maskfile:
import types
if type(maskfile) is bytes: mask = get_image(maskfile)
else: mask = maskfile
else: mask = model_circle(last_ring, nx, nx)
# references
refi = []
numref = EMUtil.get_image_count(refim)
# IMAGES ARE SQUARES! center is in SPIDER convention
cnx = nx/2+1
cny = cnx
mode = "F"
#precalculate rings
numr = Numrinit(first_ring, last_ring, rstep, mode)
wr = ringwe(numr, mode)
# reference images
params = []
#read all data
data = EMData.read_images(stack)
nima = len(data)
# prepare the reference
ima.to_zero()
for j in range(numref):
temp = EMData()
temp.read_image(refim, j)
# eve, odd, numer of even, number of images. After frc, totav
refi.append([temp, ima.copy(), 0])
seed(rand_seed)
again = True