Skip to content

Commit

Permalink
revert safe_division for sparx
Browse files Browse the repository at this point in the history
  • Loading branch information
shadowwalkersb committed Jul 18, 2018
1 parent 036a899 commit 425ce17
Show file tree
Hide file tree
Showing 67 changed files with 1,991 additions and 2,126 deletions.
21 changes: 10 additions & 11 deletions sparx/bin/sx3dvariability.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
#
from past.utils import old_div
from builtins import range
from EMAN2 import *
from sparx import *
Expand Down Expand Up @@ -243,7 +242,7 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):
img = get_image(stack, 0)
nx = img.get_xsize()
ny = img.get_ysize()
if(min(nx, ny) < options.window): keepgoing = 0
if int(options.decimate*min(nx, ny)+0.5)< options.window:keepgoing = 0
keepgoing = bcast_number_to_all(keepgoing, main_node, MPI_COMM_WORLD)
if keepgoing == 0: ERROR("The target window size cannot be larger than the size of decimated image", "sx3dvariability", 1, myid)

Expand Down Expand Up @@ -311,9 +310,9 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):
ERROR("Input stack is not prepared for symmetry, please follow instructions", "sx3dvariability", 1, myid)
from utilities import get_symt
i = len(get_symt(options.sym))
if((old_div(nima,i))*i != nima):
if((nima/i)*i != nima):
ERROR("The length of the input stack is incorrect for symmetry processing", "sx3dvariability", 1, myid)
symbaselen = old_div(nima,i)
symbaselen = nima/i
else: symbaselen = nima
else:
nima = 0
Expand All @@ -336,14 +335,14 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):
nx = int(nx*options.decimate+0.5)
ny = int(ny*options.decimate+0.5)
else:
nx = int(options.window*options.decimate+0.5)
nx = int(options.window*decimate+0.5)
ny = nx
Tracker["nx"] = nx
Tracker["ny"] = ny
Tracker["nz"] = nx

symbaselen = bcast_number_to_all(symbaselen)
if radiuspca == -1: radiuspca = old_div(nx,2)-2
if radiuspca == -1: radiuspca = nx/2-2

if myid == main_node: log_main.add("%-70s: %d\n"%("Number of projection", nima))
img_begin, img_end = MPI_start_end(nima, number_of_proc, myid)
Expand Down Expand Up @@ -455,7 +454,7 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):
log_main.add("%-70s: %.2f\n"%("Finding neighboring projections lasted [s]", time()-t2))
log_main.add("%-70s: %d\n"%("Number of groups processed on the main node", len(proj_list)))
if options.VERBOSE:
log_main.add("Grouping projections took: ", old_div((time()-t2),60) , "[min]")
log_main.add("Grouping projections took: ", (time()-t2)/60 , "[min]")
log_main.add("Number of groups on main node: ", len(proj_list))
mpi_barrier(MPI_COMM_WORLD)

Expand Down Expand Up @@ -529,7 +528,7 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):

if not options.no_norm:
#print grp_imgdata[j].get_xsize()
mask = model_circle(old_div(nx,2)-2, nx, nx)
mask = model_circle(nx/2-2, nx, nx)
for k in range(img_per_grp):
ave, std, minn, maxx = Util.infomask(grp_imgdata[k], mask, False)
grp_imgdata[k] -= ave
Expand Down Expand Up @@ -582,7 +581,7 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):

var = model_blank(nx,ny)
for q in grp_imgdata: Util.add_img2( var, q )
Util.mul_scalar( var, old_div(1.0,(len(grp_imgdata)-1)))
Util.mul_scalar( var, 1.0/(len(grp_imgdata)-1))
# Switch to std dev
var = square_root(threshold(var))
#if options.CTF: ave, var = avgvar_ctf(grp_imgdata, mode="a")
Expand Down Expand Up @@ -749,7 +748,7 @@ def params_3D_2D_NEW(phi, theta, psi, s2x, s2y, mirror):
if myid == main_node:
from fundamentals import fpol
res = fpol(res, int(org_nx*options.decimate+0.5), int(org_ny*options.decimate+0.5), int(org_nz*options.decimate+0.5))
res = resample(res, old_div(1.,options.decimate))
res = resample(res, 1./options.decimate)
set_pixel_size(res, 1.0)
res.write_image(os.path.join(options.output_dir, options.var3D))
log_main.add("%-70s: %.2f\n"%("Reconstructing 3D variability took [s]", time()-t6))
Expand Down
19 changes: 9 additions & 10 deletions sparx/bin/sxcenter_projections.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@


from __future__ import print_function
from past.utils import old_div
from builtins import range
from builtins import object
from EMAN2 import *
Expand Down Expand Up @@ -64,7 +63,7 @@ def fuselowf(vs, fq):
a = vs[0].copy()
for i in range(1,n):
Util.add_img(a, vs[i])
Util.mul_scalar(a, old_div(1.0,float(n)))
Util.mul_scalar(a, 1.0/float(n))
a = filt_tophatl(a, fq)
for i in range(n):
vs[i] = fft(Util.addn_img(a, filt_tophath(vs[i], fq)))
Expand Down Expand Up @@ -187,8 +186,8 @@ def run3Dalignment(paramsdict, partids, partstack, outputdir, procid, myid, main
ali3d_options.ts = paramsdict["ts"]
ali3d_options.xr = paramsdict["xr"]
# low pass filter is applied to shrank data, so it has to be adjusted
ali3d_options.fl = old_div(paramsdict["lowpass"],paramsdict["shrink"])
ali3d_options.initfl = old_div(paramsdict["initialfl"],paramsdict["shrink"])
ali3d_options.fl = paramsdict["lowpass"]/paramsdict["shrink"]
ali3d_options.initfl = paramsdict["initialfl"]/paramsdict["shrink"]
ali3d_options.aa = paramsdict["falloff"]
ali3d_options.maxit = paramsdict["maxit"]
ali3d_options.mask3D = paramsdict["mask3D"]
Expand All @@ -199,7 +198,7 @@ def run3Dalignment(paramsdict, partids, partstack, outputdir, procid, myid, main
projdata = getindexdata(paramsdict["stack"], partids, partstack, myid, nproc)
onx = projdata[0].get_xsize()
last_ring = ali3d_options.ou
if last_ring < 0: last_ring = int(old_div(onx,2)) - 2
if last_ring < 0: last_ring = int(onx/2) - 2
mask2D = model_circle(last_ring,onx,onx) - model_circle(ali3d_options.ir,onx,onx)
if(shrinkage < 1.0):
# get the new size
Expand Down Expand Up @@ -274,8 +273,8 @@ def run3Dalignment(paramsdict, partids, partstack, outputdir, procid, myid, main
# store params
if(myid == main_node):
for im in range(nima):
params[im][0] = old_div(params[im][0],shrinkage) +oldshifts[im][0]
params[im][1] = old_div(params[im][1],shrinkage) +oldshifts[im][1]
params[im][0] = params[im][0]/shrinkage +oldshifts[im][0]
params[im][1] = params[im][1]/shrinkage +oldshifts[im][1]
line = strftime("%Y-%m-%d_%H:%M:%S", localtime()) + " =>"
print(line,"Executed successfully: ","3D alignment"," number of images:%7d"%len(params))
write_text_row(params, os.path.join(outputdir,"params.txt") )
Expand Down Expand Up @@ -381,11 +380,11 @@ def main():
if ali3d_options.CTF:
i = a.get_attr('ctf')
pixel_size = i.apix
fq = old_div(pixel_size,fq)
fq = pixel_size/fq
else:
pixel_size = 1.0
# No pixel size, fusing computed as 5 Fourier pixels
fq = old_div(5.0,nxinit)
fq = 5.0/nxinit
del a
else:
total_stack = 0
Expand Down Expand Up @@ -437,7 +436,7 @@ def main():
else: yr = options.yr

delta = float(options.delta)
if(delta <= 0.0): delta = "%f"%round(degrees(atan(old_div(1.0,float(radi)))), 2)
if(delta <= 0.0): delta = "%f"%round(degrees(atan(1.0/float(radi))), 2)
else: delta = "%f"%delta

paramsdict = { "stack":stack,"delta":delta, "ts":"1.0", "xr":xr, "an":angular_neighborhood, \
Expand Down
9 changes: 4 additions & 5 deletions sparx/bin/sxchains.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@
#
#

from past.utils import old_div
from builtins import range
import global_def
from global_def import *
Expand Down Expand Up @@ -63,7 +62,7 @@ def TotalDistance(city, lccc):

def reverse(city, n):
nct = len(city)
nn = old_div((1+ ((n[1]-n[0]) % nct)),2) # half the lenght of the segment to be reversed
nn = (1+ ((n[1]-n[0]) % nct))/2 # half the lenght of the segment to be reversed
# the segment is reversed in the following way n[0]<->n[1], n[0]+1<->n[1]-1, n[0]+2<->n[1]-2,...
# Start at the ends of the segment and swap pairs of cities, moving towards the center.
for j in range(nn):
Expand Down Expand Up @@ -101,7 +100,7 @@ def tsp(lccc):

# ncity = 100 # Number of cities to visit
from math import sqrt
ncity = int( old_div((1+sqrt(1+8*len(lccc))),2) ) # Number of cities to visit
ncity = int( (1+sqrt(1+8*len(lccc)))/2 ) # Number of cities to visit
# sanity check
if( ncity*(ncity-1)/2 != len(lccc) ): return [-1]

Expand Down Expand Up @@ -153,7 +152,7 @@ def tsp(lccc):
de = Distance(city[n[2]], city[n[1]], lccc) + Distance(city[n[3]], city[n[0]], lccc)\
- Distance(city[n[2]], city[n[0]], lccc) - Distance(city[n[3]] ,city[n[1]], lccc)

if de<0 or exp(old_div(-de,T))>random.rand(): # Metropolis
if de<0 or exp(-de/T)>random.rand(): # Metropolis
accepted += 1
dist += de
reverse(city, n)
Expand All @@ -169,7 +168,7 @@ def tsp(lccc):
de += Distance( city[n[0]], city[n[4]], lccc) + Distance( city[n[1]], city[n[5]], lccc) \
+ Distance( city[n[2]], city[n[3]], lccc)

if de<0 or exp(old_div(-de,T))>random.rand(): # Metropolis
if de<0 or exp(-de/T)>random.rand(): # Metropolis
accepted += 1
dist += de
city = transpt(city, n)
Expand Down
19 changes: 9 additions & 10 deletions sparx/bin/sxcompute_isac_avg.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
# New version of sort3D.
#
from __future__ import print_function
from past.utils import old_div
from builtins import range
import os
import sys
Expand Down Expand Up @@ -101,22 +100,22 @@ def adjust_pw_to_model(image, pixel_size, roo):
c2 = 15.0
c3 = 0.2
c4 = -1.0
c5 = old_div(1.,5.)
c5 = 1./5.
c6 = 0.25 # six params are fitted to Yifan channel model
rot1 = rops_table(image)
fil = [None]*len(rot1)
if roo is None: # adjusted to the analytic model, See Penczek Methods Enzymol 2010
pu = []
for ifreq in range(len(rot1)):
x = float(ifreq)/float(len(rot1))/pixel_size
v = exp(c1+old_div(c2,(old_div(x,c3)+1)**2)) + exp(c4-0.5*((old_div((x-c5),c6**2))**2))
v = exp(c1+c2/(x/c3+1)**2) + exp(c4-0.5*(((x-c5)/c6**2)**2))
pu.append(v)
s =sum(pu)
for ifreq in range(len(rot1)): fil[ifreq] = sqrt(old_div(pu[ifreq],(rot1[ifreq]*s)))
for ifreq in range(len(rot1)): fil[ifreq] = sqrt(pu[ifreq]/(rot1[ifreq]*s))
else: # adjusted to a given 1-d rotational averaged pw2
if roo[0]<0.1 or roo[0]>1.: s =sum(roo)
else: s=1.0
for ifreq in range(len(rot1)):fil[ifreq] = sqrt(old_div(roo[ifreq],(rot1[ifreq]*s)))
for ifreq in range(len(rot1)):fil[ifreq] = sqrt(roo[ifreq]/(rot1[ifreq]*s))
return filt_table(image, fil)

def get_optimistic_res(frc):
Expand All @@ -134,13 +133,13 @@ def get_optimistic_res(frc):

def apply_enhancement(avg, B_start, pixel_size, user_defined_Bfactor):
guinierline = rot_avg_table(power(periodogram(avg),.5))
freq_max = old_div(1.,(2.*pixel_size))
freq_min = old_div(1.,B_start)
freq_max = 1./(2.*pixel_size)
freq_min = 1./B_start
b, junk, ifreqmin, ifreqmax = compute_bfactor(guinierline, freq_min, freq_max, pixel_size)
print(ifreqmin, ifreqmax)
global_b = b*4. #
if user_defined_Bfactor < 0.0: global_b = user_defined_Bfactor
sigma_of_inverse = sqrt(old_div(2.,global_b))
sigma_of_inverse = sqrt(2./global_b)
avg = filt_gaussinv(fft(avg), sigma_of_inverse)
return avg, global_b

Expand Down Expand Up @@ -338,7 +337,7 @@ def main():
Tracker = wrap_mpi_bcast(Tracker, Blockdata["main_node"], communicator = MPI_COMM_WORLD)

#print(Tracker["constants"]["pixel_size"], "pixel_size")
x_range = max(Tracker["constants"]["xrange"], int(old_div(1.,Tracker["ini_shrink"]))+1)
x_range = max(Tracker["constants"]["xrange"], int(1./Tracker["ini_shrink"])+1)
y_range = x_range

if(Blockdata["myid"] == Blockdata["main_node"]): parameters = read_text_row(os.path.join(Tracker["constants"]["isac_dir"], "all_parameters.txt"))
Expand All @@ -363,7 +362,7 @@ def main():
abs_id = members[im]
global_dict[abs_id] = [iavg, im]
P = combine_params2( init_dict[abs_id][0], init_dict[abs_id][1], init_dict[abs_id][2], init_dict[abs_id][3], \
parameters[abs_id][0], old_div(parameters[abs_id][1],Tracker["ini_shrink"]), old_div(parameters[abs_id][2],Tracker["ini_shrink"]), parameters[abs_id][3])
parameters[abs_id][0], parameters[abs_id][1]/Tracker["ini_shrink"], parameters[abs_id][2]/Tracker["ini_shrink"], parameters[abs_id][3])
if parameters[abs_id][3] ==-1:
print("WARNING: Image #{0} is an unaccounted particle with invalid 2D alignment parameters and should not be the member of any classes. Please check the consitency of input dataset.".format(abs_id)) # How to check what is wrong about mirror = -1 (Toshio 2018/01/11)
params_of_this_average.append([P[0], P[1], P[2], P[3], 1.0])
Expand Down
21 changes: 10 additions & 11 deletions sparx/bin/sxconsistency.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@

# VERSION 2 09/25/2014

from past.utils import old_div
from builtins import range
from EMAN2 import *
from sparx import *
Expand Down Expand Up @@ -93,13 +92,13 @@ def average2dtransform(data, return_avg_pixel_error=False):

# Get ave transform params
H = Transform({"type":"2D"})
H.set_matrix([old_div(sum_cosa,sqrtP), old_div(sum_sina,sqrtP), 0.0, sx, old_div(-sum_sina,sqrtP), old_div(sum_cosa,sqrtP), 0.0, sy, 0.0, 0.0, 1.0, 0.0])
H.set_matrix([sum_cosa/sqrtP, sum_sina/sqrtP, 0.0, sx, -sum_sina/sqrtP, sum_cosa/sqrtP, 0.0, sy, 0.0, 0.0, 1.0, 0.0])
dd = H.get_params("2D")

H = Transform({"type":"2D","alpha":dd[ "alpha" ],"tx":dd[ "tx" ],"ty": dd[ "ty" ],"mirror":0,"scale":1.0})
dd = H.get_params("2D")
if return_avg_pixel_error:
return old_div(sum(sqr_pixel_error),N)
return sum(sqr_pixel_error)/N
else:
return [dd[ "alpha" ], dd[ "tx" ], dd[ "ty" ]]

Expand Down Expand Up @@ -300,7 +299,7 @@ def average_trans(params):
nphi = 0.0
ntheta = 0.0
else:
ntheta = degrees(acos(old_div(nas[2],nom)))%360.0
ntheta = degrees(acos(nas[2]/nom))%360.0
if(sym>1 and ntheta>90.0): nphi = (degrees(atan2( nas[1], nas[0] ))-180.0)%qsym + 180.0
else: nphi = degrees(atan2( nas[1], nas[0] ))%qsym

Expand Down Expand Up @@ -397,7 +396,7 @@ def main():
#print ll,rw[0][j], rw[2][j]
break
if isame: thesame += 1
qt = old_div(float(thesame),nn)
qt = float(thesame)/nn
print("Proportion of the same orientations ",qt)
write_text_file([qt], os.path.join(outdir,howmanythesame) )

Expand All @@ -408,7 +407,7 @@ def main():
radius = options.ou
thresherr = options.thresherr #for chc5 1, for ribo 4#4.0
sym = int(options.sym[1:])
qsym = old_div(360.0,sym)
qsym = 360.0/sym

blocks = ['A','B','C','D']
params = {}
Expand Down Expand Up @@ -487,7 +486,7 @@ def main():
nphi = 0.0
ntheta = 0.0
else:
ntheta = degrees(acos(old_div(nas[2],nom)))%360.0
ntheta = degrees(acos(nas[2]/nom))%360.0
if(sym>1 and ntheta>90.0): nphi = (degrees(atan2( nas[1], nas[0] ))-180.0)%qsym + 180.0
else: nphi = degrees(atan2( nas[1], nas[0] ))%qsym

Expand Down Expand Up @@ -619,7 +618,7 @@ def main():
radius = options.ou
thresherr = options.thresherr
sym = int(options.sym[1:])
qsym = old_div(360.0,sym)
qsym = 360.0/sym
#params = [[None for i in xrange(3)] for j in xrange(4)]
ll = 3 # this is hardwired as we have three groups. however, I would like to keep the code general.
for jj in range(4):
Expand Down Expand Up @@ -671,7 +670,7 @@ def main():
nphi = 0.0
ntheta = 0.0
else:
ntheta = degrees(acos(old_div(nas[2],nom)))%360.0
ntheta = degrees(acos(nas[2]/nom))%360.0
if(sym>1 and ntheta>90.0): nphi = (degrees(atan2( nas[1], nas[0] ))-180.0)%qsym + 180.0
else: nphi = degrees(atan2( nas[1], nas[0] ))%qsym

Expand Down Expand Up @@ -709,7 +708,7 @@ def main():
radius = options.ou
thresherr = options.thresherr
sym = int(options.sym[1:])
qsym = old_div(360.0,sym)
qsym = 360.0/sym
#params = [[None for i in xrange(3)] for j in xrange(4)]
ll = 3 # this is hardwired as we have three groups. however, I would like to keep the code general.
for jj in range(1):
Expand Down Expand Up @@ -762,7 +761,7 @@ def main():
nphi = 0.0
ntheta = 0.0
else:
ntheta = degrees(acos(old_div(nas[2],nom)))%360.0
ntheta = degrees(acos(nas[2]/nom))%360.0
if(sym>1 and ntheta>90.0): nphi = (degrees(atan2( nas[1], nas[0] ))-180.0)%qsym + 180.0
else: nphi = degrees(atan2( nas[1], nas[0] ))%qsym

Expand Down
5 changes: 2 additions & 3 deletions sparx/bin/sxcter.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
#
from past.utils import old_div
import os
import sys
from optparse import OptionParser
Expand Down Expand Up @@ -162,11 +161,11 @@ def main():

if options.f_start >0.0:
if options.f_start <=0.5: ERROR("f_start should be in Angstrom","sxcter", 1) # exclude abs frequencies and spatial frequencies
else: freq_start = old_div(1.,options.f_start)
else: freq_start = 1./options.f_start

if options.f_stop >0.0:
if options.f_stop <=0.5: ERROR("f_stop should be in Angstrom","sxcter", 1) # exclude abs frequencies and spatial frequencies
else: freq_stop = old_div(1.,options.f_stop)
else: freq_stop = 1./options.f_stop

while True:
# --------------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 425ce17

Please sign in to comment.