Commit 8c9c6556 authored by Jakob Knollmueller's avatar Jakob Knollmueller

paper plots

parent b6ed499e
......@@ -64,7 +64,7 @@ if __name__ == '__main__':
plt.imsave('point_sex.png', (data_unmod - sex), vmax = lin_max,vmin =lin_min)
alpha = 1.4
Starblade = sb.build_starblade(data, alpha=alpha,cg_iterations=100,newton_iterations=10)
Starblade = sb.build_starblade(data, alpha=alpha, cg_steps=100, newton_steps=10)
for i in range(1000):
Starblade = sb.starblade_iteration(Starblade, samples=5)
......
......@@ -70,7 +70,7 @@ if __name__ == '__main__':
plt.imsave('log_point_sex.png', np.log((data_unmod - sex).clip(min=0.0001)), vmax = vmax,vmin =vmin)
alpha = 1.4
Starblade = sb.build_starblade(data, alpha=alpha,cg_iterations=100,newton_iterations=3)
Starblade = sb.build_starblade(data, alpha=alpha, cg_steps=100, newton_steps=3)
for i in range(10):
Starblade = sb.starblade_iteration(Starblade, samples=5)
......
This diff is collapsed.
from keras.layers import Activation, Dense, Input
from keras.layers import Conv2D, Flatten
from keras.layers import Reshape, Conv2DTranspose, BatchNormalization, merge, Dropout
from keras.models import Model
def DDCAE_model(image_size):
input_shape = (image_size, image_size,1)
batch_size = 32
kernel_size = 3
latent_dim = 8
# Encoder/Decoder number of CNN layers and filters per layer
layer_filters = [32,16]
# Build the Autoencoder Model
# First build the Encoder Model
inputs = Input(shape=input_shape, name='encoder_input')
x = inputs
# Stack of Conv2D blocks
# Notes:
# 1) Use Batch Normalization before ReLU on deep networks
# 2) Use MaxPooling2D as alternative to strides>1
# - faster but not as good as strides>1
for filters in layer_filters:
x = Conv2D(filters=filters,
kernel_size=kernel_size,
strides=1,
activation='relu',
padding='same',
data_format='channels_last')(x)
# x = BatchNormalization()(x)
x = Dropout(0.1)(x)
for filters in layer_filters[::-1]:
x = Conv2DTranspose(filters=filters,
kernel_size=kernel_size,
strides=1,
activation='relu',
padding='same',
data_format='channels_last')(x)
# x = BatchNormalization()(x)
x = Dropout(0.1)(x)
x = Conv2DTranspose(filters=1,
kernel_size=kernel_size,
padding='same')(x)
outputs = Activation('sigmoid', name='decoder_output')(x)
outputs = merge([inputs,outputs],mode='mul')
return Model(inputs, outputs, name='autoencoder')
......@@ -9,32 +9,40 @@ from matplotlib.colors import LogNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1 import AxesGrid
from astropy.io import fits
from scipy.ndimage import zoom
import starblade as sb
from DDCAE import DDCAE_model
np.random.seed(42)
if __name__ == '__main__':
path = 'demos/data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
path = '../demos/data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
data = fits.open(path)[1].data
data = data.clip(min=0.001)
# data = data.clip(min=0.001)
data = np.ndarray.astype(data, float)
alpha = 1.28
#data = zoom(data, 0.25)
data = data.clip(min=0.001)
hdu = fits.PrimaryHDU(data)
hdul = fits.HDUList([hdu])
hdul.writeto('my_hubble_data.fits',overwrite=True)
alpha = 1.1
myEnergy = sb.build_starblade(data, alpha=alpha)
for i in range(10):
myEnergy = sb.starblade_iteration(myEnergy)
A = FFTSmoothingOperator(myEnergy.s.domain, sigma=2.)
for i in range(30):
myEnergy = sb.starblade_iteration(myEnergy, samples=5,newton_steps=30,cg_steps=10, sampling_steps=100)
print i
A = FFTSmoothingOperator(myEnergy.s.domain, sigma=.002)
plt.magma()
size = 15
vmin = data.min()+0.01
vmax = 0.01*data.max()
vmax = data.max()*0.01
plt.figure()
plt.title('diffuse emission', size=size)
plt.axis('off')
......@@ -81,16 +89,16 @@ if __name__ == '__main__':
plt.suptitle('zoomed in section', size=size)
# fig.tight_layout()
vmin = data.min() + 0.0001
vmax = 0.001 * data.max()
im = ax[0].imshow(data[600:700, 650:750], norm=LogNorm(vmin=vmin, vmax=vmax))
vmin_small = data.min() + 0.0001
vmax_small = 0.001 * data.max()
im = ax[0].imshow(data[600:700, 650:750], norm=LogNorm(vmin=vmin_small, vmax=vmax_small))
ax[0].set_title('data', size=15)
ax[0].axis('off')
ax[1].imshow(exp(myEnergy.s).val[600:700, 650:750], norm=LogNorm(vmin=vmin, vmax=vmax))
ax[1].imshow(exp(myEnergy.s).val[600:700, 650:750], norm=LogNorm(vmin=vmin_small, vmax=vmax_small))
ax[1].set_title('diffuse', size=15)
ax[1].axis('off')
ax[2].imshow(exp(myEnergy.u).val[600:700, 650:750], norm=LogNorm(vmin=vmin, vmax=vmax))
ax[2].imshow(exp(myEnergy.u).val[600:700, 650:750], norm=LogNorm(vmin=vmin_small, vmax=vmax_small))
ax[2].set_title('point-like', size=15)
ax[2].axis('off')
......@@ -140,7 +148,7 @@ if __name__ == '__main__':
k_lengths = power.domain[0].k_lengths
plt.plot(k_lengths, power.val, 'k-', label='diffuse')
plt.plot(k_lengths, power_d.val, 'k:', label='data')
plt.plot(k_lengths, power_d.val, 'k+-', label='data')
plt.legend()
plt.yscale('log')
plt.xscale('log')
......@@ -148,3 +156,55 @@ if __name__ == '__main__':
plt.ylabel('power',size=15)
plt.xlabel('harmonic mode',size=15)
plt.savefig('hubble_log_power.pdf')
DDCAE=DDCAE_model(None)
DDCAE.load_weights('DDCAC',by_name=True)
ddd = data.reshape([1,1000,1000,1])
separation = DDCAE.predict(ddd)
plt.figure()
plt.title('DAE diffuse emission', size=size)
plt.axis('off')
ax = plt.gca()
im = ax.imshow(separation[0,:,:,0],norm=LogNorm(vmin=vmin, vmax=vmax))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('flux', size=size)
plt.tight_layout()
plt.savefig('DAE_hubble_diffuse.pdf')
plt.figure()
plt.title('DAE point-like emission', size=size)
plt.axis('off')
ax = plt.gca()
im = ax.imshow(A(ift.Field(myEnergy.position.domain,data-separation[0,:,:,0])).val,norm=LogNorm(vmin=vmin, vmax=vmax))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('flux', size=size)
plt.tight_layout()
plt.savefig('DAE_hubble_point_like.pdf')
sextracted = fits.open('check_hubble64.fits')[0].data
plt.figure()
plt.title('SExtractor diffuse emission', size=size)
plt.axis('off')
ax = plt.gca()
im = ax.imshow(sextracted.clip(0.001),norm=LogNorm(vmin=vmin, vmax=vmax))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('flux', size=size)
plt.tight_layout()
plt.savefig('SExtractor_hubble_diffuse.pdf')
plt.figure()
plt.title('SExtractor point-like emission', size=size)
plt.axis('off')
ax = plt.gca()
im = ax.imshow(A(ift.Field(myEnergy.position.domain,(data-sextracted))).val.clip(0.001),norm=LogNorm(vmin=vmin, vmax=vmax))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('flux', size=size)
plt.tight_layout()
plt.savefig('SExtractor_hubble_point_like.pdf')
......@@ -17,7 +17,7 @@
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
from nifty4 import Energy, Field, log, exp, DiagonalOperator,\
create_power_operator, SandwichOperator, ScalingOperator, InversionEnabler
create_power_operator, SandwichOperator, ScalingOperator, InversionEnabler, SamplingEnabler
from nifty4.library import WienerFilterCurvature
from nifty4.library.nonlinearities import PositiveTanh, Tanh
......@@ -56,6 +56,7 @@ class StarbladeEnergy(Energy):
self.parameters = parameters
self.inverter = parameters['inverter']
self.sampling_inverter = parameters['sampling_inverter']
self.d = parameters['data']
self.FFT = parameters['FFT']
self.power_spectrum = parameters['power_spectrum']
......@@ -115,10 +116,11 @@ class StarbladeEnergy(Energy):
point = self.q * exp(-self.u) * self.u_p ** 2
# R = self.FFT.inverse * self.s_p
# N = self.correlation
N_inv = DiagonalOperator(point + 1/self.var_x )#+ 2*self.a_p))
O_x = DiagonalOperator(Field(self.position.domain,val=1./self.var_x))
N_inv = DiagonalOperator(point )#+ 2*self.a_p))
R = ScalingOperator(1., point.domain)
S_p = DiagonalOperator(self.s_p)
my_S_inv = SandwichOperator.make(self.FFT.adjoint.inverse.adjoint * S_p, self.correlation.inverse)
curv = InversionEnabler(N_inv + my_S_inv, self.inverter)
curv = InversionEnabler(SamplingEnabler(my_S_inv+N_inv, O_x, self.sampling_inverter), self.inverter)
return curv
# return WienerFilterCurvature(R=R, N=N, S=S, inverter=self.inverter)
......@@ -23,8 +23,8 @@ import nifty4 as ift
from .starblade_energy import StarbladeEnergy
from .starblade_kl import StarbladeKL
def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=100, newton_iterations = 3,
manual_power_spectrum = None):
def build_starblade(data, alpha=1.5, q=1e-40, cg_steps=100, newton_steps = 3,
manual_power_spectrum = None, sampling_steps = 100):
""" Setting up the StarbladeEnergy for the given data and parameters
Parameters
----------
......@@ -34,9 +34,9 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=100, newton_iteratio
The slope parameter of the point source prior (default: 1.5).
q : float
The cutoff parameter of the point source prior (default: 1e-40).
cg_iterations : int
cg_steps : int
Maximum number of conjugate gradient iterations for numerical operator inversion (default: 500).
newton_iterations : int
newton_steps : int
Number of consecutive Newton steps within one algorithmic step.(default: 3)
manual_power_spectrum : None, Field or callable
Option to set a manual power spectrum which is kept constant during the separation.
......@@ -63,21 +63,24 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=100, newton_iteratio
initial_x = ift.Field(s_space, val=-1.)
alpha = ift.Field(s_space, val=alpha)
q = ift.Field(s_space, val=q)
ICI = ift.GradientNormController(iteration_limit=cg_iterations,
ICI = ift.GradientNormController(iteration_limit=cg_steps,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(controller=ICI)
IC_samples = ift.GradientNormController(iteration_limit=sampling_steps,
tol_abs_gradnorm=1e-5)
sampling_inverter = ift.ConjugateGradient(controller=IC_samples)
parameters = dict(data=data, power_spectrum=initial_spectrum,
alpha=alpha, q=q,
inverter=inverter, FFT=FFT,
newton_iterations=newton_iterations, update_power=update_power)
newton_iterations=newton_steps, sampling_inverter=sampling_inverter, update_power=update_power)
Starblade = StarbladeEnergy(position=initial_x, parameters=parameters)
return Starblade
def starblade_iteration(starblade, samples=3):
def starblade_iteration(starblade, samples=3, cg_steps=10, newton_steps=3, sampling_steps=100):
""" Performing one Newton minimization step
Parameters
----------
......@@ -88,16 +91,26 @@ def starblade_iteration(starblade, samples=3):
"""
controller = ift.GradientNormController(name="Newton",
tol_abs_gradnorm=1e-8,
iteration_limit=starblade.newton_iterations)
iteration_limit=newton_steps)
minimizer = ift.RelaxedNewton(controller=controller)
ICI = ift.GradientNormController(iteration_limit=cg_steps,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(controller=ICI)
IC_samples = ift.GradientNormController(iteration_limit=sampling_steps,
tol_abs_gradnorm=1e-5)
sampling_inverter = ift.ConjugateGradient(controller=IC_samples)
# minimizer = ift.VL_BFGS(controller=controller)
energy = starblade
para = starblade.parameters
para['inverter'] = inverter
para['sampling_inverter'] = sampling_inverter
energy = StarbladeEnergy(starblade.position,parameters=para)
sample_list = []
for i in range(samples):
sample = energy.curvature.inverse.draw_sample()
sample_list.append(sample)
if len(sample_list)>0:
energy = StarbladeKL(starblade.position, samples=sample_list, parameters=starblade.parameters)
energy = StarbladeKL(starblade.position, samples=sample_list, parameters=energy.parameters)
else:
energy = starblade
energy, convergence = minimizer(energy)
......@@ -140,9 +153,9 @@ def build_multi_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500,
"""
MultiStarblade = []
for i in range(data.shape[-1]):
starblade = build_starblade(data[...,i],alpha=alpha, q=q,
cg_iterations=cg_iterations,
newton_iterations=newton_iterations,
starblade = build_starblade(data[...,i], alpha=alpha, q=q,
cg_steps=cg_iterations,
newton_steps=newton_iterations,
manual_power_spectrum = manual_power_spectrum)
MultiStarblade.append(starblade)
return MultiStarblade
......
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