Commit 814dec9b authored by Maximilian Beckers's avatar Maximilian Beckers
Browse files

Python3 support

parent 0993f51c
#Author: Maximilian Beckers, EMBL Heidelberg, Sachse Group (2017)
#import some stuff
from confidenceMapUtil import mapUtil, locscaleUtil, FDRutil
from confidenceMapUtil import FDRutil, confidenceMapMain
import numpy as np
import argparse, os, sys
import subprocess
......@@ -151,12 +151,13 @@ def main():
else:
mpi = False;
if args.stepSize > args.window_size_locscale:
print("Step Size cannot be bigger than the window_size. Job is killed ...")
return;
if (args.stepSize is not None) & (args.window_size_locscale is not None):
if args.stepSize > args.window_size_locscale:
print("Step Size cannot be bigger than the window_size. Job is killed ...")
return;
#run the actual analysis
confidenceMap, locFiltMap, locScaleMap, binMap, maskedMap = FDRutil.calculateConfidenceMap(mapData, args.apix, args.noiseBox, args.testProc, args.ecdf, args.lowPassFilter, args.method, args.window_size, locResMapData, meanMapData, varMapData, args.fdr, modelMapData, stepSize, windowSizeLocScale, mpi);
confidenceMap, locFiltMap, locScaleMap, binMap, maskedMap = confidenceMapMain.calculateConfidenceMap(mapData, args.apix, args.noiseBox, args.testProc, args.ecdf, args.lowPassFilter, args.method, args.window_size, locResMapData, meanMapData, varMapData, args.fdr, modelMapData, stepSize, windowSizeLocScale, mpi);
if locFiltMap is not None:
locFiltMapMRC = mrcfile.new(splitFilename[0] + '_locFilt.mrc', overwrite=True);
......@@ -196,7 +197,7 @@ def main():
end = time.time();
totalRuntime = end -start;
mapUtil.printSummary(args, totalRuntime);
FDRutil.printSummary(args, totalRuntime);
if (__name__ == "__main__"):
main()
......
......@@ -4,163 +4,9 @@ import math
import gc
import os
import sys
import mapUtil
import locscaleUtil
#Author: Maximilian Beckers, EMBL Heidelberg, Sachse Group
#--------------------------------------------------------------------------
def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter_resolution, method, window_size, locResMap,
meanMap, varMap, fdr, modelMap, stepSize, windowSizeLocScale, mpi):
#*********************************************
#******* this function calc. confMaps ********
#*********************************************
#print pixelSize and give feedback to the user
output = "Pixel size was read as " + "%.2f" %apix + " Angstrom. If this is incorrect, please specify with -p pixelSize";
print(output);
# get boxCoordinates
if noiseBox is None:
boxCoord = 0;
else:
boxCoord = noiseBox;
# set test procdure
if testProc is not None:
testProc = testProc;
else:
testProc = 'rightSided';
# set ECDF
if ecdf:
ECDF = 1;
else:
ECDF = 0;
sizeMap = em_map.shape;
if lowPassFilter_resolution is not None:
frequencyMap = calculate_frequency_map(em_map);
providedRes = apix/float(lowPassFilter_resolution);
em_map = lowPassFilter(np.fft.rfftn(em_map), frequencyMap, providedRes, em_map.shape);
# handle FDR correction procedure
if method is not None:
method = method;
else:
# default is Benjamini-Yekutieli
method = 'BY';
if window_size is not None:
wn = window_size;
wn = int(wn);
if wn < 20:
print("Provided window size is quite small. Please think about potential inaccuracies of your noise estimates!");
else:
wn = max(int(0.05 * sizeMap[0]), 10);
if windowSizeLocScale is not None:
wn_locscale = windowSizeLocScale;
if window_size is not None:
print("Window size for noise estimation is automatically set to the size of LocScale sliding window.");
else:
wn_locscale = None;
if stepSize is None:
stepSize = 5;
# generate a circular Mask
sphere_radius = (np.min(sizeMap) // 2);
circularMaskData = mapUtil.makeCircularMask(np.copy(em_map), sphere_radius);
# plot locations of noise estimation
if modelMap is None:
pp = mapUtil.makeDiagnosticPlot(em_map, wn, False, boxCoord);
pp.savefig("diag_image.pdf");
pp.close();
else:
pp = mapUtil.makeDiagnosticPlot(em_map, wn, True, boxCoord);
pp.savefig("diag_image.pdf");
pp.close();
# estimate noise statistics
if ((locResMap is None) & (modelMap is None)): # if no local Resolution map is given,don't do any filtration
checkNormality(em_map, wn, boxCoord);
mean, var, _ = estimateNoiseFromMap(em_map, wn, boxCoord);
if varMap is not None:
var = varMap;
if meanMap is not None:
mean = meanMap;
if np.isscalar(mean) and np.isscalar(var):
output = "Estimated noise statistics: mean: " + repr(mean) + " and variance: " + repr(var);
else:
output = "Using user provided noise statistics";
print(output);
locFiltMap = None;
locScaleMap = None;
elif (locResMap is not None) & (modelMap is None): # do localFiltration and estimate statistics from this map
checkNormality(em_map, wn, boxCoord);
em_map, mean, var, ECDF = mapUtil.localFiltration(em_map, locResMap, apix, True, wn, boxCoord, ECDF);
locFiltMap = em_map;
locScaleMap = None;
else:
em_map, mean, var, ECDF = locscaleUtil.launch_amplitude_scaling(em_map, modelMap, apix, stepSize, wn_locscale, wn, method, locResMap, boxCoord, mpi, ECDF );
locScaleMap = em_map;
locFiltMap = None;
# calculate the qMap
qMap = calcQMap(em_map, mean, var, ECDF, wn, boxCoord, circularMaskData, method, testProc);
# if a explicit thresholding is wished, do so
if fdr is not None:
fdr = fdr;
# threshold the qMap
binMap = binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, em_map);
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
maskedMap = np.multiply(maskedMap, circularMaskData);
if (locResMap is None) & (modelMap is None): # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
print(output);
else:
# threshold the qMap
fdr = 0.01;
binMap = binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, np.copy(em_map));
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
if (locResMap is None) & (modelMap is None): # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
print(output);
binMap = None;
maskedMap = None;
# invert qMap for visualization tools
confidenceMap = np.subtract(1.0, qMap);
# apply lowpass-filtered mask to maps
confidenceMap = np.multiply(confidenceMap, circularMaskData);
return confidenceMap, locFiltMap, locScaleMap, binMap, maskedMap;
#-------------------------------------------------------------------------------------
def estimateNoiseFromMap(map, windowSize, boxCoord):
......@@ -438,7 +284,7 @@ def normalizeMap(map, mean, var):
#****************************************
if np.isscalar(var):
normMap = np.subtract(map, mean);
normMap = np.subtract(map, mean);
normMap = np.multiply(normMap, (1.0/(math.sqrt(var))));
else: #if local variances are known, use them
var[var==0] = 1000;
......
import numpy as np
import subprocess
import math
import gc
import os
import sys
from confidenceMapUtil import mapUtil, locscaleUtil, FDRutil
#--------------------------------------------------------------------------
def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter_resolution, method, window_size, locResMap,
meanMap, varMap, fdr, modelMap, stepSize, windowSizeLocScale, mpi):
#*********************************************
#******* this function calc. confMaps ********
#*********************************************
#print pixelSize and give feedback to the user
output = "Pixel size was read as " + "%.2f" %apix + " Angstrom. If this is incorrect, please specify with -p pixelSize";
print(output);
# get boxCoordinates
if noiseBox is None:
boxCoord = 0;
else:
boxCoord = noiseBox;
# set test procdure
if testProc is not None:
testProc = testProc;
else:
testProc = 'rightSided';
# set ECDF
if ecdf:
ECDF = 1;
else:
ECDF = 0;
sizeMap = em_map.shape;
if lowPassFilter_resolution is not None:
frequencyMap = FDRutil.calculate_frequency_map(em_map);
providedRes = apix/float(lowPassFilter_resolution);
em_map = FDRutil.lowPassFilter(np.fft.rfftn(em_map), frequencyMap, providedRes, em_map.shape);
# handle FDR correction procedure
if method is not None:
method = method;
else:
# default is Benjamini-Yekutieli
method = 'BY';
if window_size is not None:
wn = window_size;
wn = int(wn);
if wn < 20:
print("Provided window size is quite small. Please think about potential inaccuracies of your noise estimates!");
else:
wn = max(int(0.05 * sizeMap[0]), 10);
if windowSizeLocScale is not None:
wn_locscale = windowSizeLocScale;
if window_size is not None:
print("Window size for noise estimation is automatically set to the size of LocScale sliding window.");
else:
wn_locscale = None;
if stepSize is None:
stepSize = 5;
# generate a circular Mask
sphere_radius = (np.min(sizeMap) // 2);
circularMaskData = mapUtil.makeCircularMask(np.copy(em_map), sphere_radius);
# plot locations of noise estimation
if modelMap is None:
pp = mapUtil.makeDiagnosticPlot(em_map, wn, False, boxCoord);
pp.savefig("diag_image.pdf");
pp.close();
else:
pp = mapUtil.makeDiagnosticPlot(em_map, wn, True, boxCoord);
pp.savefig("diag_image.pdf");
pp.close();
# estimate noise statistics
if ((locResMap is None) & (modelMap is None)): # if no local Resolution map is given,don't do any filtration
FDRutil.checkNormality(em_map, wn, boxCoord);
mean, var, _ = FDRutil.estimateNoiseFromMap(em_map, wn, boxCoord);
if varMap is not None:
var = varMap;
if meanMap is not None:
mean = meanMap;
if np.isscalar(mean) and np.isscalar(var):
output = "Estimated noise statistics: mean: " + repr(mean) + " and variance: " + repr(var);
else:
output = "Using user provided noise statistics";
print(output);
locFiltMap = None;
locScaleMap = None;
elif (locResMap is not None) & (modelMap is None): # do localFiltration and estimate statistics from this map
FDRutil.checkNormality(em_map, wn, boxCoord);
em_map, mean, var, ECDF = mapUtil.localFiltration(em_map, locResMap, apix, True, wn, boxCoord, ECDF);
locFiltMap = em_map;
locScaleMap = None;
else:
em_map, mean, var, ECDF = locscaleUtil.launch_amplitude_scaling(em_map, modelMap, apix, stepSize, wn_locscale, wn, method, locResMap, boxCoord, mpi, ECDF );
locScaleMap = em_map;
locFiltMap = None;
# calculate the qMap
qMap = FDRutil.calcQMap(em_map, mean, var, ECDF, wn, boxCoord, circularMaskData, method, testProc);
# if a explicit thresholding is wished, do so
if fdr is not None:
fdr = fdr;
# threshold the qMap
binMap = FDRutil.binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, em_map);
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
maskedMap = np.multiply(maskedMap, circularMaskData);
if (locResMap is None) & (modelMap is None): # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
print(output);
else:
# threshold the qMap
fdr = 0.01;
binMap = FDRutil.binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, np.copy(em_map));
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
if (locResMap is None) & (modelMap is None): # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
print(output);
binMap = None;
maskedMap = None;
# invert qMap for visualization tools
confidenceMap = np.subtract(1.0, qMap);
# apply lowpass-filtered mask to maps
confidenceMap = np.multiply(confidenceMap, circularMaskData);
return confidenceMap, locFiltMap, locScaleMap, binMap, maskedMap;
\ No newline at end of file
import numpy as np
import FDRutil
from mapUtil import *
from confidenceMapUtil import FDRutil
import argparse, math, os, sys
from argparse import RawTextHelpFormatter
import time
......@@ -23,16 +22,16 @@ def pad_or_crop_volume(vol, dim_pad=None, pad_value = None, crop_volume=False):
crop_volume = True
if crop_volume:
crop_vol = vol[vol.shape[0]/2-dim_pad[0]/2:vol.shape[0]/2+dim_pad[0]/2+dim_pad[0]%2, :, :]
crop_vol = crop_vol[:, vol.shape[1]/2-dim_pad[1]/2:vol.shape[1]/2+dim_pad[1]/2+dim_pad[1]%2, :]
crop_vol = crop_vol[:, :, vol.shape[2]/2-dim_pad[2]/2:vol.shape[2]/2+dim_pad[2]/2+dim_pad[2]%2]
crop_vol = vol[int(vol.shape[0]/2-dim_pad[0]/2):int(vol.shape[0]/2+dim_pad[0]/2+dim_pad[0]%2), :, :]
crop_vol = crop_vol[:, int(vol.shape[1]/2-dim_pad[1]/2):int(vol.shape[1]/2+dim_pad[1]/2+dim_pad[1]%2), :]
crop_vol = crop_vol[:, :, int(vol.shape[2]/2-dim_pad[2]/2):int(vol.shape[2]/2+dim_pad[2]/2+dim_pad[2]%2)]
return crop_vol
else:
pad_vol = np.pad(vol, ((dim_pad[0]/2-vol.shape[0]/2, dim_pad[0]/2-vol.shape[0]/2+dim_pad[0]%2), (0,0), (0,0) ), 'constant', constant_values=(pad_value,))
pad_vol = np.pad(pad_vol, ((0,0), (dim_pad[1]/2-vol.shape[1]/2, dim_pad[1]/2-vol.shape[1]/2+dim_pad[1]%2 ), (0,0)), 'constant', constant_values=(pad_value,))
pad_vol = np.pad(pad_vol, ((0,0), (0,0), (dim_pad[2]/2-vol.shape[2]/2, dim_pad[2]/2-vol.shape[2]/2+dim_pad[2]%2)), 'constant', constant_values=(pad_value,))
pad_vol = np.pad(vol, ((int(dim_pad[0]/2-vol.shape[0]/2), int(dim_pad[0]/2-vol.shape[0]/2+dim_pad[0]%2)), (0,0), (0,0) ), 'constant', constant_values=(pad_value,))
pad_vol = np.pad(pad_vol, ((0,0), (int(dim_pad[1]/2-vol.shape[1]/2), int(dim_pad[1]/2-vol.shape[1]/2+dim_pad[1]%2) ), (0,0)), 'constant', constant_values=(pad_value,))
pad_vol = np.pad(pad_vol, ((0,0), (0,0), (int(dim_pad[2]/2-vol.shape[2]/2), int(dim_pad[2]/2-vol.shape[2]/2+dim_pad[2]%2))), 'constant', constant_values=(pad_value,))
return pad_vol
......@@ -188,11 +187,11 @@ def calculate_scaled_map(emmap, modmap, mask, wn, wn_locscale, apix, locFilt, lo
#prepare windows of particle for scaling
frequencyMap_mapWindow = FDRutil.calculate_frequency_map(np.zeros((wn_locscale, wn_locscale, wn_locscale)));
numSteps = len(xrange(0, sizeMap[0] - int(wn_locscale), stepSize))*len(xrange(0, sizeMap[1] - int(wn_locscale), stepSize))*len(xrange(0, sizeMap[2] - int(wn_locscale), stepSize));
numSteps = len(range(0, sizeMap[0] - int(wn_locscale), stepSize))*len(range(0, sizeMap[1] - int(wn_locscale), stepSize))*len(range(0, sizeMap[2] - int(wn_locscale), stepSize));
counterSteps = 0;
for k in xrange(0, sizeMap[0] - int(wn_locscale), stepSize):
for j in xrange(0, sizeMap[1] - int(wn_locscale), stepSize):
for i in xrange(0, sizeMap[2] - int(wn_locscale), stepSize):
for k in range(0, sizeMap[0] - int(wn_locscale), stepSize):
for j in range(0, sizeMap[1] - int(wn_locscale), stepSize):
for i in range(0, sizeMap[2] - int(wn_locscale), stepSize):
#print progress
counterSteps = counterSteps + 1;
......@@ -267,8 +266,12 @@ def calculate_scaled_map(emmap, modmap, mask, wn, wn_locscale, apix, locFilt, lo
sharpened_map[k + halfStep : k + halfStep + stepSize, j + halfStep : j + halfStep + stepSize, i + halfStep : i + halfStep + stepSize] = np.copy(map_b_sharpened[halfStep:halfStep+stepSize, halfStep:halfStep+stepSize, halfStep:halfStep+stepSize]);
sharpened_mean_vals[k + halfStep : k + halfStep + stepSize, j + halfStep : j + halfStep + stepSize, i + halfStep : i + halfStep + stepSize] = mean;
sharpened_var_vals[k + halfStep : k + halfStep + stepSize, j + halfStep : j + halfStep + stepSize, i + halfStep : i + halfStep + stepSize] = var;
sharpened_ecdf_vals[k + halfStep : k + halfStep + stepSize, j + halfStep : j + halfStep + stepSize, i + halfStep : i + halfStep + stepSize] = ecdf[halfStep:halfStep+stepSize, halfStep:halfStep+stepSize, halfStep:halfStep+stepSize];
if ecdfBool:
sharpened_ecdf_vals[k + halfStep : k + halfStep + stepSize, j + halfStep : j + halfStep + stepSize, i + halfStep : i + halfStep + stepSize] = ecdf[halfStep:halfStep+stepSize, halfStep:halfStep+stepSize, halfStep:halfStep+stepSize];
else:
sharpened_ecdf_vals[k + halfStep: k + halfStep + stepSize, j + halfStep: j + halfStep + stepSize,
i + halfStep: i + halfStep + stepSize] = 0.0;
return sharpened_map, sharpened_mean_vals, sharpened_var_vals, sharpened_ecdf_vals;
def get_central_scaled_pixel_vals_after_scaling(emmap, modmap, masked_xyz_locs, wn, wn_locscale, apix, locFilt, locResMap, boxCoord, ecdfBool):
......@@ -518,7 +521,7 @@ def launch_amplitude_scaling(em_map, model_map, apix, stepSize, wn_locscale, wn,
if not mpi:
stepSize = int(stepSize);
if stepSize == 1:
LocScaleVol, meanVol, varVol, ecdfVol = run_window_function_including_scaling(emmap, modmap, mask, wn, wn_locscale , apix, locFilt, locResMap, boxCoord, ecdf);
LocScaleVol, meanVol, varVol, ecdfVol = run_window_function_including_scaling(emmap, modmap, mask, wn, wn_locscale , apix, locFilt, locResMap, boxCoord, ecdf);
elif stepSize <= 0:
print("Invalid step size parameter. It has to be greater than 0! Quit program ...");
return;
......
from FDRutil import *
from confidenceMapUtil import FDRutil
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
......@@ -57,14 +57,14 @@ def localFiltration(map, locResMap, apix, localVariance, windowSize, boxCoord, E
#get initial noise statistics
initMapData = np.copy(map);
initMean, initVar, _ = estimateNoiseFromMap(initMapData, windowSize, boxCoord);
initMean, initVar, _ = FDRutil.estimateNoiseFromMap(initMapData, windowSize, boxCoord);
noiseMapData = np.random.normal(initMean, math.sqrt(initVar), (100, 100, 100));
#do FFT of the respective map
mapFFT = np.fft.rfftn(map);
#get frequency map
frequencyMap = calculate_frequency_map(map);
frequencyMap = FDRutil.calculate_frequency_map(map);
# Initial call to print 0% progress
#printProgressBar(counter, numRes, prefix = 'Progress:', suffix = 'Complete', bar_length = 50)
......@@ -89,7 +89,7 @@ def localFiltration(map, locResMap, apix, localVariance, windowSize, boxCoord, E
xInd, yInd, zInd = indices[0], indices[1], indices[2];
#do local filtration
tmpFilteredMapData = lowPassFilter(mapFFT, frequencyMap, tmpRes, map.shape);
tmpFilteredMapData = FDRutil.lowPassFilter(mapFFT, frequencyMap, tmpRes, map.shape);
#set the filtered voxels
filteredMapData[xInd, yInd, zInd] = tmpFilteredMapData[xInd, yInd, zInd];
......@@ -97,7 +97,7 @@ def localFiltration(map, locResMap, apix, localVariance, windowSize, boxCoord, E
else:
xInd, yInd, zInd = indices[0], indices[1], indices[2];
#do local filtration
tmpFilteredMapData = lowPassFilter(mapFFT, frequencyMap, tmpRes, map.shape);
tmpFilteredMapData = FDRutil.lowPassFilter(mapFFT, frequencyMap, tmpRes, map.shape);
#set the filtered voxels
filteredMapData[xInd, yInd, zInd] = tmpFilteredMapData[xInd, yInd, zInd];
if localVariance == True:
......@@ -105,13 +105,13 @@ def localFiltration(map, locResMap, apix, localVariance, windowSize, boxCoord, E
if ECDF == 1:
#if ecdf shall be used, use if to p-vals
tmpECDF, sampleSort = estimateECDFFromMap(tmpFilteredMapData, windowSize, boxCoord);
tmpECDF, sampleSort = FDRutil.estimateECDFFromMap(tmpFilteredMapData, windowSize, boxCoord);
vecECDF = np.interp(tmpFilteredMapData[xInd, yInd, zInd], sampleSort, tmpECDF, left=0.0, right=1.0);
ECDFmap[xInd, yInd, zInd] = vecECDF;
else:
ECDFmap = 0;
tmpMean, tmpVar, _ = estimateNoiseFromMap(tmpFilteredMapData, windowSize, boxCoord);
tmpMean, tmpVar, _ = FDRutil.estimateNoiseFromMap(tmpFilteredMapData, windowSize, boxCoord);
mean[xInd, yInd, zInd] = tmpMean;
var[xInd, yInd, zInd] = tmpVar;
......@@ -238,6 +238,7 @@ def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1,
sys.stdout.write('\n');
sys.stdout.flush();
return;
#---------------------------------------------------------------------------------
def makeCircularMask(map, sphereRadius):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment